arXiv:0704.1286v1  [math.NA]  10 Apr 2007
IMA Journal of Numerical Analysis (2007) Page 1 of 26
doi: 10.1093/imanum/dri017
Study of a finite volume- finite element scheme for a nuclear transport
model
CATHERINE CHOQUET
LATP, Université Aix-Marseille 3, 13013 Marseille Cedex 20, France
SÉBASTIEN ZIMMERMANN
Laboratoire de mathématiques, École Centrale de Lyon, 36 avenue Guy de Collongue, 69134
Ecully, France.
[Received on 24 July 2007]
We consider a problem of nuclear waste contamination. It takes into account the thermal effects. The
temperature and the contaminant’s concentration fulfill convection-diffusion-reaction equations. The ve-
locity and the pressure in the flow satisfy the Darcy equation, with a viscosity depending on both con-
centration and temperature. The equations are nonlinear and strongly coupled. Using both finite volume
and nonconforming finite element methods, we introduce a scheme adapted to this problem. We prove
the stability and convergence of this scheme and give some error estimates.
Keywords: porous media, miscible flow, nonconforming finite element, finite volume.
1. Introduction
A part of the high-level nuclear waste is now stored in environmentally safe locations. One has to
consider the eventuality of a leakage through the engineering and geological barriers. It may cause
the contamination of underground water sources far away from the original repository’s location. In the
present paper, we consider such a problem of nuclear waste contamination in the basement. We take into
account the thermal effects. The evolution in time of the temperature and of the contaminants concen-
tration is then governed by convection-diffusion-reaction equations. The velocity and the pressure in the
flow satisfy the Darcy equation, with a viscosity depending on both concentrations and temperature in
a nonlinear way. The velocity satisfies an incompressibility constraint. We introduce a scheme adapted
to this problem. We use both finite volume and nonconforming finite element methods. It ensures that
a maximum principle holds and that the associated linear systems have good-conditioned matrices. We
prove the stability and convergence of the scheme and give some error estimates.
Let us briefly point out some previous works. A complete model coupling concentrations and pres-
sure equations is very seldom studied, since the system is strongly coupled. Instead, each equation
is considered separately. In the general context of convection-diffusion-reaction equations, numerous
schemes are available (see [15] or [19] and the references therein). Finite difference schemes are some-
times used for the convective term (in [20] for instance). But they are not adapted for the complex
geometry of a reservoir. More recently, finite volume methods were developed and analysed. Let us just
cite the book [2] or [13] and the references therein. Finite elements (for the diffusive term) and finite
volumes (for the convective term) are coupled for instance in [3, 8] . In convection dominant problems,
the equations are of degenerate parabolic type. This setting is considered in [9, 18] . The reaction
terms are specifically studied through operator splitting methods in [14] . Now, in the specific context of
porous media flow, we mention [17] who consider only the evolution of the pressure. In [11, 12] a more
complete set of equations is used, and a mixed finite element approximation is developed. We stress that
IMA Journal of Numerical Analysis c© Institute of Mathematics and its Applications 2007; all rights reserved.
http://arxiv.org/abs/0704.1286v1
2 of 26 C. CHOQUET AND S. ZIMMERMANN
in all theses works, as in most, the thermic effects are not taken into account.
The present paper is organized as follows. Section 2 is devoted to the derivation of the model. In
section 3, we introduce the discrete tools used in this paper. It allows us to define the numerical scheme
of section 5. The analysis of the scheme uses the properties of section 4. We then prove the stability and
convergence of the scheme, in sections 6 and 7 respectively. We conclude with some error estimates in
section 8.
2. Model of contamination
The thickness of the medium is significantly smaller than its length and width. Hence it is reasonable to
average the medium properties vertically and to describe the far-field repository by a polyhedral domain
Ω of R2 with a smooth boundary ∂Ω . It is characterized by a porosity φ and a permeability tensor K.
The time interval of interest is [0,T ]. We denote by p the pressure, by (ai)
i=1 the concentrations of the
Nr radionuclides involved in the flow and by θ the temperature. The Darcy velocity is represented by
u. We assume a miscible and incompressible displacement. Due to the mass and energy conservation
principles, the flow is governed by the following system satisfied in Ω × [0,T ], with i = 1, ..,Nr −1 (see
[10]).
φ Ri ∂tai + div(ai u)− div(φ Dc ∇ai) = si − sai −λi Ri φ ai +
j=1, j 6=i
ki j λ j R j φ a j , (2.1)
φ Cp ∂tθ + div(θ u)− div(φ Dθ ∇θ ) =−sθ − s(θ −θ∗) , (2.2)
divu = s, u+
(a j)
j=1 ,θ
) ∇p = f. (2.3)
In (2.1) the retardation factors Ri > 0 are due to the sorption mechanism. The real λ−1i > 0 denotes the
half life time of radionuclide i. The term −λi Ri φ ai describes the radioactive decay of the i–th specy.
Meanwhile, the quantity ∑ j 6=i ki j λ j R j φ a j is created by radioactive filiation. The molecular diffusion
effects are given by the coefficient Dc > 0. The contamination is represented by the source term si and
s=∑Nri=1 si. In (2.2) the coefficient Cp > 0 is the relative specific heat of the porous medium. The thermic
diffusion coefficient is denoted by Dθ > 0. The real θ∗ > 0 is a reference temperature. The constitutive
relation (2.3) is the Darcy law and f is a density of body forces. For a large range of temperatures µ has
the form
µ(a,θ ) = µR(a)exp
where µR is a nonlinear function. For instance, in the Koval model for a two-species mixture [16], we
µR(a) = µ(0)(1+(M1/4 − 1)a1)−4
where M = µ(0)/µ(1) is the mobility ratio.
We notice that the equations (2.1)-(2.3) are strongly coupled. Moreover, every concentration equation
(2.1) involves a different time scale. Therefore, it is difficult to build a numerical scheme that captures
all the physical phenomena. We have to transform these equations. We first assume that only serial or
parallel first-order reactions occur, so that ki j = yi with y1 = 0. Next, following [4], we assume that no
two isotopes have identical decay rates and we set
c1 = a1, ci = ai +
yl+1λl
λl −λi
a j for i = 2, ..,Nr. (2.4)
Study of a finite volume-finite element scheme for a nuclear transport model 3 of 26
Lastly, without losing any mathematical difficulty (see remark 2.2 below), we set Ri = 1 for i= 1, ..,Nr−
1 and φ = 1, Cp = 1. We also set sci = si +∑
∏i−1l= j
yl+1λl
λl−λi
s j and κ(c,θ ) = K/µ
i=1 ,θ
. The
contamination problem is now modelized by the following parabolic-elliptic system
∂tci + div(ci u)−Dc ∆ci = sci − sci −λi ci , (2.5)
∂tθ + div(θ u)−Dθ ∆θ =−sθ − s(θ −θ∗) , (2.6)
divu = s, u+κ(c,θ )∇p = f , (2.7)
with i = 1, ..,Nr − 1. These equations are completed with the boundary and initial conditions
∇ci ·ν = 0 , ∇θ ·ν = 0 , u ·n = 0 , (2.8)
ci|t=0 = ci0, θ |t=0 = θ0. (2.9)
The pressure p is normalized by
Ω pdx = 0. Equations (2.5) are all similar. Thus, for the sake of
simplicity, we will assume that there is only one. We set Nr = 2 and c = c1, c0 = c
0, sc = sc1 , λ = λ1.
The results of this paper readily extend to the general case.
We conclude with some notations and hypothesis. Let D be a bounded open set of Rk with k > 1.
We denote by C ∞0 (D) the set of functions that are continuous on D together with all their derivatives,
and have a compact support in D. For p ∈ {2,∞}, we use the Lebesgue spaces
Lp(D),‖.‖Lp(D)
Lp(D),‖.‖Lp(D)
with Lp = (Lp)2. We also use the Sobolev spaces W p,q(D) for p ∈ [1,∞[ and
q ∈ [1,∞[. In the case D = Ω we use the following conventions. We drop the domain dependancy.
We denote by |.| (resp. ‖.‖∞) the norms associated to L2 = L2(Ω) and L2 = L2(Ω) (resp. L∞ and
L∞). We set L20 = {v ∈ L2;
Ω v(x)dx = 0}. For p ∈ [1,∞[ we define (H p,‖.‖p) and (Hp,‖.‖p) with
H p = W p,2 and Hp = (H p)2. Now let (X , |.|) be a Banach space. The functions g : [0,T ] → X such
that t → ‖g(t)‖X is continuous (resp. bounded and square integrable) form the set C (0,T ;X) (resp.
L∞(0,T ;X) and L2(0,T ;X)). The associated norm for the space L∞(0,T ;X) (resp. L2(0,T ;X)) is
defined by ‖g‖L∞(0,T ;X) = supt∈[0,T ] ‖g(t)‖X (resp. ‖g‖L2(0,T ;X) =
0 ‖ f (t)‖2X dt
). Finally in all
computations we use C > 0 as a generic constant. It depends only on the data of the problem.
We assume the following regularities for the data in (2.5)–(2.7)
κ ∈W 1,∞((0,1)× (0,∞)) , s, sc, sθ ∈ L2 , f ∈ C (0,T ;L2). (2.10)
We also assume that κ > κin f with κin f > 0. For the initial data, we assume that c0 ∈ H1, θ0 ∈ H1, and
that we have a.e. in Ω
0 6 c0(x)6 1 , θ− 6 θ0(x)6 θ+ (2.11)
with θ− > 0. Finally we assume that we have a.e. in Ω
2s(x)+λ > sc(x)> 0 , 2(θ−−θ ∗)s(x)+ sθ (x)6 0 , 2(θ+−θ ∗)s(x)+ sθ (x)> 0. (2.12)
These conditions ensure a maximum principle (proposition 6.1 below).
REMARK 2.1 We have assumed that first-order reactions occur, so that the coefficients ki j in (2.1)
depend only on i. If ki j depends on i and j, one can still uncouple the equations by iterating the transfor-
mation (2.4), provided that k1 j = 0 for j = 2, ..,Nr −1. This assumption means that the first long-lasting
isotope disappears and is not created anymore. It is satisfied by many radionuclides.
4 of 26 C. CHOQUET AND S. ZIMMERMANN
REMARK 2.2 We have assumed that the retardation factors R j are identical. If it is not the case, the
difficulty and the approach remain the same. Indeed, let us consider the Fourier transform of (2.1). For
a Fourier mode â j(k, t) we obtain
(φ R j â j) =−φ (λ j R j + k2 Dc − ik u) â j + yi R j−1 φ â j−1 =−λ ′j R j φ â j + yi R j−1φ â j−1
with λ ′j = λ j +(k
2 Dc − ik u)/R j for j = 1, ..,Nr − 1. A transform analogous to (2.4) uncouples the
problem. By taking the partial differential equation counterpart, we obtain an equation similar to (2.5).
3. Discrete tools
3.1 Mesh and discrete spaces
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 π2 , so that xK ∈ K. The set of the edges of the triangle K ∈ Th is EK . The symbol
nK,σ denotes the unit normal vector 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 is |σ |. For each edge σ ∈ E inth let Kσ
and Lσ be the two triangles having σ in common; we set dσ = d(xKσ ,xLσ ). For all σ ∈ E exth only the
triangle Kσ located inside Ω is defined and we set dσ = d(xKσ ,xσ ). Then for all σ ∈ Eh we set τσ =
We assume the following on the mesh (see [13]): 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. (3.1)
We define on the mesh the following spaces. The usual space for finite volume schemes is
P0 = {q ∈ L2 ; ∀K ∈ Th, q|K is a constant}.
For any function qh ∈ P0 and any K ∈ Th we set qK = qh|K . We also consider
Pd1 = {q ∈ L2 ; ∀K ∈ Th, q|K is affine} ,
Pc1 = {qh ∈ Pd1 ; qh is continuous over Ω} ,
Pnc1 = {qh ∈ Pd1 ; ∀σ ∈ E inth , qh is continuous at the middle of σ}.
We have Pc1 ⊂ H1. On the other hand Pnc1 6⊂ H1, but Pnc1 ⊂ H1d with
H1d = {q ∈ L2 ; ∀K ∈ Th, q|K ∈ H1(K)}.
Thus we define ∇̃h : H1d → L2 by setting
∀qh ∈ H1d , ∀K ∈ Th, ∇̃hqh|K = ∇qh|K (3.2)
and the associated norm ‖.‖1,h is given by
∀qh ∈ H1d , ‖qh‖1,h = (|qh|2 + |∇̃hqh|2)1/2.
We then have the following Poincaré-like inequality for the space Pnc1 ∩L20 (see [1]).
Study of a finite volume-finite element scheme for a nuclear transport model 5 of 26
PROPOSITION 3.1 There exists C > 0 such that |qh|6C |∇̃hqh| for all qh ∈ Pnc1 ∩L20.
We also define discrete analogues of the norms H1 and H−1 for the space P0 by setting
‖qh‖h =
σ∈E inth
τσ (qLσ − qKσ )2
and ‖qh‖−1,h = sup
ψh∈P0
(qh,ψh)
‖ψh‖h
for any function qh ∈ P0. Note that for any ph ∈P0 and qh ∈ P0, (ph,qh)6 ‖qh‖−1,h‖qh‖h. The following
Poincaré-like inequality holds for the space P0 ∩L20 (see [13]).
PROPOSITION 3.2 There exists C > 0 such that |qh|6C‖qh‖h for all qh ∈ P0 ∩L20.
Finally we set P0 = (P0)2, Pd1 = (P
2 and use the Raviart-Thomas spaces [7]
RTd0 = {vh ∈ Pd1 ; ∀K ∈ Th , ∀σ ∈ EK , vh|K ·nK,σ is a constant} ,
RT0 = {vh ∈ RTd0 ; ∀σ ∈ E inth , vh|Kσ ·nKσ ,σ = vh|Lσ ·nKσ ,σ and vh ·n|∂Ω = 0}.
For all vh ∈ RT0, K ∈ Th and σ ∈ EK , we set (vh ·nK,σ )σ = vh|K ·nK,σ .
3.2 Projection operators
We associate with the spaces of section 3.1 some projection operators. First, we define ΠPc1 : H
d → Pc1
by setting
∀q ∈ H1d , ∀φh ∈ Pc1 ,
∇(ΠPc1 q),∇φh
= (∇q,∇φh). (3.3)
Next, we consider the space P0. Let Cd = {qh ∈ L2 ; qh is equal a.e. to a continuous function}. We
define ΠP0 : L
2 → P0 and Π̃P0 : Cd → P0 by setting
(ΠP0 p)K =
p(x)dx , (Π̃P0q)K = q(xK) , (3.4)
for all p ∈ L2, q ∈ Cd and K ∈ Th. We also set ΠP0 = (ΠP0)2. For the space Pnc1 , we define Π̃Pnc1 : L
Pnc1 and ΠPnc1 : H
1 → Pnc1 . For all p ∈ L2 and q ∈ H1, Π̃Pnc1 p and ΠPnc1 q satisfy
∀ψh ∈ Pnc1 , (Π̃Pnc1 p,ψh) = (p,ψh) ; ∀σ ∈ Eh ,
(ΠPnc1 q)dσ =
qdσ . (3.5)
For the space RT0, we define Π̃RT0 : L
2 → RT0 and ΠRT0 : H1 → RT0. For all v ∈ L2 and w ∈ H1,
Π̃RT0v and ΠRT0w satisfy
∀wh ∈ RT0 , (Π̃RT0v,wh) = (v,wh) ; ∀σ ∈ E
(w−ΠRT0w) ·nKσ ,σ dσ = 0. (3.6)
The operators ΠP0 , Π̃Pnc1 (resp. ΠP0 , Π̃RT0) are L
2 (resp. L2) projection operators. They are stable for
the L2 (resp. L2) norms. The operators Π̃P0 , ΠPnc1 and ΠRT0 are interpolation operators. The following
estimates are classical ([6] p.109 and [7]).
PROPOSITION 3.3 There exists C > 0 such that for all q ∈ H1 and v ∈ H1
|q−ΠP0q|6C h‖q‖1, |v−ΠRT0v|6C h‖v‖1.
6 of 26 C. CHOQUET AND S. ZIMMERMANN
For all p ∈ H1 and q ∈ H2 we have
|p−ΠPnc1 p|6C h‖p‖1 , |∇̃h(q−ΠPnc1 q)|6C h‖q‖2.
For all q ∈ Hd1 we have
|q−ΠPc1 q|6C h‖q‖1,h.
Finally, using the Sobolev embedding theorem, one checks that
|ΠP0q− Π̃P0q|6C h‖q‖W1,r (3.7)
for all q ∈W 1,r with r > 1 (see [22]).
3.3 Discrete operators
Equations (2.5)–(2.7) use the differential operators gradient, divergence and laplacian. We have to define
analogous operators in the discrete setting. The discrete gradient operator ∇h : Pnc1 →P0 is the restriction
to Pnc1 of the operator ∇̃h given by (3.2). The discrete divergence operator divh : P0 → Pnc1 is defined by
∀σ ∈ E inth , (divh vh)(xσ ) =
3 |σ |
|Kσ |+ |Lσ |
(vLσ − vKσ ) ·nK,σ ,
∀σ ∈ E exth , (divh vh)(xσ ) =−
3 |σ |
|Kσ |
vKσ ·nK,σ ,
for all vh ∈ P0. It is adjoint to ∇h (proposition 4.1 below). The discrete laplacian operator ∆h : P0 → P0
is the usual one for finite volume schemes (see [13]). For all qh ∈ P0 and K ∈ Th we have
|K| ∑
σ∈EK∩E inth
qLσ − qKσ
. (3.8)
Let us now consider the convection terms in (2.5) and (2.6). We define b̃ : H1 ×H1 → L2 by
b̃(v,q) = div(qv) (3.9)
for all q ∈ H1 and v ∈ H1. In order to define a discrete counterpart to b̃ we use the classical upwind
scheme (see [13]). The discrete operator b̃h : RT0 ×P0 → P0 is such that
b̃h(vh,qh)
|K| ∑
σ∈EK∩E inth
(vh ·nK,σ )+σ qK +(vh ·nK,σ )−σ qLσ
(3.10)
for all vh ∈ RT0, qh ∈ P0 and K ∈ Th. We have set a+ = max(a,0) and a− = min(a,0) for all a ∈ R.
Integrating by parts the convection terms also leads to consider b : L2 ×L2 ×L∞ →R defined by
b(v, p,q) =−
pv ·∇qdx (3.11)
for all v ∈ L2, p ∈ L2 and q ∈ L∞. The discrete counterpart is bh : RT0 ×P0×P0 → R with
bh(vh, ph,qh) = ∑
σ∈EK∩E inth
(vh ·nK,σ )+σ pK +(vh ·nK,σ )−σ pLσ
(3.12)
for all vh ∈ RT0, ph ∈ P0 and qh ∈ P0.
Study of a finite volume-finite element scheme for a nuclear transport model 7 of 26
4. Properties of the discrete operators
The properties of the discrete operators are analogous to the ones satisfied by their continuous counter-
part. The gradient and divergence operators are adjoint. For the operators ∇h and divh we state in [21]
the following.
PROPOSITION 4.1 For all vh ∈ P0 and qh ∈ Pnc1 we have (vh,∇hqh) =−(qh,divhvh).
Let us now consider the convection terms. Let q ∈ L∞ ∩H1, v ∈ L2 with divv ∈ L2 and divv(x) > 0
a.e. in Ω . We obtain b(v,q,q) =
2/2)divvdx > 0 by integration by parts. For bh we state in [21] a
similar result.
PROPOSITION 4.2 Let vh ∈ RT0 with divvh > 0. We have bh(vh,qh,qh)> 0 for all qh ∈ P0 .
The following stability properties are used to prove the error estimates in section 8.
PROPOSITION 4.3 There exists C > 0 such that for all ph ∈ P0, qh ∈ P0 and vh ∈ RT0 with divvh = 0
|bh(vh, ph,qh)|6C |vh|‖ph‖h‖qh‖h. (4.1)
There exists C > 0 such that for all ph ∈ P0, qh ∈ P0 ∩L20, vh ∈ RT0
|bh(vh, ph,qh)|6C (|vh|‖ph‖h + |divvh|‖ph‖∞ )‖qh‖h. (4.2)
PROOF. For all K ∈ Th and σ ∈ EK ∩E inth we write
(vh ·nK,σ )+σ pK +(vh ·nK,σ )−σ pLσ = (vh ·nK,σ )σ pK −|(vh ·nK,σ )σ |(pLσ − pK).
Thus (3.12) reads bh(vh, ph,qh) = S1 + S2 with
S1 =− ∑
σ∈EK∩E inth
|σ | |(vh ·nK,σ )σ |(pLσ − pK) , S2 = ∑
pK qK ∑
σ∈EK∩E inth
|σ |(vh ·nK,σ )σ .
Using the Cauchy-Schwarz inequality we write
|S1|=
∣∣∣ ∑
σ∈E inth
|σ | |vh(xσ ) ·nK,σ |(pLσ − pK)(qLσ − qK)
6 h‖vh‖∞
σ∈E inth
(pLσ − pKσ )2
)1/2( ∑
σ∈E inth
(qLσ − qKσ )2
Since vh ∈ RT0 ⊂ (Pd1 )2 we have h‖vh‖∞ 6C |vh| ([6] p. 112). Moreover (3.1) implies ∑σ∈E inth (pLσ −
pKσ )
6C ∑σ∈E inth τσ (pLσ − pKσ )
2 =C‖ph‖2h and ∑σ∈E inth (qLσ − qKσ )
6C‖qh‖2h. Thus
|S1|6C |vh|‖ph‖h ‖qh‖h. (4.3)
We now consider S2. We have (vh ·nK,σ )σ = 0 for all K ∈ Th and σ ∈ EK ∩E exth . Thus we write
σ∈EK∩E inth
|σ |(vh ·nK,σ )σ = ∑
|σ |(vh ·nK,σ )σ =
divvh dx.
8 of 26 C. CHOQUET AND S. ZIMMERMANN
It gives the following relation.
S2 = ∑
pK qK
divvh dx =
ph qh divvh dx.
Thus if divvh = 0 then S2 = 0 and estimate (4.3) gives (4.1). Let us prove (4.2). Since qh ∈ P0 ∩L20 we
can apply proposition 3.2. Using the Cauchy-Schwarz inequality we get
|S2|6 ‖ph‖∞ |qh| |divvh|6C‖ph‖∞ ‖qh‖h |divvh|.
This latter estimate together with (4.3) gives (4.2).
Lastly, we claim that b̃h is a consistent approximation of b̃ [21].
PROPOSITION 4.4 Let r > 0. There exists C > 0 such that for all functions q ∈ H2 and v ∈ H1+r with
v ·n|∂Ω = 0
‖ΠP0 b̃(v,q)− b̃h(ΠRT0v,Π̃P0q)‖−1,h 6C h‖q‖1‖v‖1+r.
Let us now consider the discrete laplacian. We have a coercivity and stability result.
PROPOSITION 4.5 For all ph ∈ P0 and qh ∈ P0, we have
−(∆h ph, ph) = ‖ph‖2h , |(∆h ph,qh)|6 ‖ph‖h‖qh‖h.
PROOF. Definition (3.8) implies
(∆h ph,qh) = ∑
σ∈EK∩E inth
τσ (pLσ − pKσ ) =− ∑
σ∈E inth
τσ (pLσ − pKσ )(qLσ − qKσ ). (4.4)
Setting qh = ph gives the first part of the result. Using the Cauchy-Schwarz inequality, we get the second
We also deduce from (4.4) the following property.
PROPOSITION 4.6 For all ph ∈ P0 and qh ∈ P0 we have (∆h ph,qh) = (ph,∆hqh).
Lastly, we state that ∆h is a consistent approximation of the laplacian. The proof follows the lines of the
one of proposition 1.14 in [22].
PROPOSITION 4.7 There exists C > 0 such that for all q ∈ H2 with ∇q ·n|∂Ω = 0 we have
‖ΠP0(∆q)−∆h(Π̃P0q)‖−1,h 6C h‖q‖2.
5. The finite element-finite volume scheme
We now introduce the scheme for (2.5)-(2.9). The interval [0,T ] is split with a constant time step
k = T/N. We set [0,T ] =
m=0[tm, tm+1] with tm = mk. The time derivatives are approximated using a
first order Euler scheme. The convection terms are discretized semi-implicitly in time and the other ones
in an implicit way. We set sh = ΠP0s, s
h = ΠP0sc, s
h = ΠP0sθ and f
h = ΠP0f(tm) for all m ∈ {0, . . . ,N}.
Since ΠP0 (resp. ΠP0) is stable for the L
2 (resp. L2) norm we have
|sh|6 |s| , |sch|6 |sc| , |sθh |6 |sθ | , |fmh |6 |f(tm)|6 ‖f‖L∞(0,T ;L2). (5.1)
Study of a finite volume-finite element scheme for a nuclear transport model 9 of 26
The initial values are c0h = ΠP0c0 and θ
h = ΠP0θ0. Then for all n ∈ {0, . . . ,N − 1}, the quantities
cn+1h ∈ P0, θ
h ∈ P0, p
h ∈ P
1 ∩L20, u
h ∈ RT0 are the solutions of the following problem.
cn+1h − c
−Dc ∆hcn+1h = s
h − (sh +λ )cn+1h − b̃h(u
h ), (5.2)
θ n+1h −θ
−Dθ ∆hθ n+1h =−s
h − sh (θ
h −θ∗)− b̃h(u
h ), (5.3)
divh(κn+1h ∇h p
h ) = divh f
h − Π̃Pnc1 sh , (5.4)
un+1h = Π̃RT0(f
h ∇h p
h ), (5.5)
with κn+1h = κ(c
h ) ∈ P0. This term is defined thanks to proposition 6.1 below. Note also that
the boundary conditions are implicitly included in the definition of the discrete operators (section 3.3).
The existence of a unique solution to (5.2) and (5.3) is classical (see [13]). Since κn+1h > κmin > 0 and
pn+1h ∈ L
0 equation (5.4) also has a unique solution (see [6]). We have a discrete equivalent for the
divergence condition (2.7).
PROPOSITION 5.1 For all m ∈ {1, . . . ,N} we have divumh = sh .
PROOF. Let m ∈ {1, . . . ,N} and n = m− 1. We compare the solution of (5.4)–(5.5) with the solution
of the following mixed hybrid problem. Let E0 =
µh : ∪σ∈Eh → R ; ∀σ ∈ Eh , µh|σ is constant
. Then
ũmh ∈ RT
0 , p
h ∈ P0 and λ
h ∈ E0 are the solution of (see [7])
∀vh ∈ RTd0 , (ũmh ,vh)+ ∑
κmK ∑
σ ′∈EK
|σ ′|λ mσ ′ (vh|K ·nK,σ ′)− ∑
|K|κmK pmK divvh|K = (fmh ,vh) , (5.6)
∀µh ∈ E0, ∑
µh (ũmh ·n)dσ = 0 , ∀K ∈ Th ,
div ũmh dx =
sdx , (5.7)
and p̃mh ∈ Pnc1 is defined by
h dσ = λ
σ for all σ ∈ Eh. Let σ ∈ Eh. We define φσ ∈ Pnc1 by setting
φσ (xσ ) = 1 and φσ (xσ ′) = 0 for all σ ′ ∈ Eh\{σ}. We set vh = ∇hφσ ∈ P0 ⊂ RTd0 in (5.6). We have
κmK ∑
σ ′∈EK
|σ ′|λ mσ ′ ∇hφσ |K ·nK,σ ′ = ∑
κmK ∇hφσ |K · ∑
σ ′∈EK
|σ ′|λ mσ ′ nK,σ ′
and according to the gradient formula
σ ′∈EK
|σ ′|λ mσ ′ nK,σ ′ = ∑
σ ′∈EK
p̃mh nK,σ ′ dσ
∇h p̃
h dx.
Thus we get from (5.6)
(ũmh ,∇hφσ )+ (κ
h ∇h p̃
h ,∇hφσ ) = (f
h ,∇hφσ ). (5.8)
The first term in (5.8) is treated as follows. Integrating by parts we get
(ũmh ,∇hφσ ) =−(φσ ,div ũmh )+ ∑
σ ′∈EK
φσ (ũmh |K ·nK,σ ′)dσ ′.
Since (5.7) implies that ũmh ∈ RT0, we have
σ ′∈EK
φσ (ũmh ·nK,σ ′)dσ ′ = ∑
σ∈E inth
|σ |φσ (xσ )(ũmh |Lσ ·nKσ ,σ − ũmh |Kσ ·nKσ ,σ ) = 0.
10 of 26 C. CHOQUET AND S. ZIMMERMANN
Thus (ũmh ,∇hφσ ) =−(φσ ,div ũmh ). Then, using (5.7), we get
(ũmh ,∇hφσ ) =−(φσ ,sh) =−(φσ ,Π̃Pnc1 sh).
Furthermore, according to proposition 4.1, we have (κmh ∇h p̃
h ,∇hφσ ) = −
φσ ,divh(κmh ∇h p̃
(fmh ,∇hφσ ) =−(φσ ,divhfmh ). Hence we deduce from (5.8) that
∀φσ ∈ Pnc1 ,
φσ ,divh(κmh ∇h p̃
h )− divhf
h + Π̃Pnc1 s
Since (φσ )σ∈Eh is a basis of P
1 , we get divh(κ
h ∇h p̃
h ) = divhf
h − Π̃Pnc1 s
h . Thus, by (5.4), there exists
a real C such that p̃mh = p
h +C. We now compare ũ
h with u
h . Since for all vh ∈ RT0 we have
κmK ∑
σ ′∈EK
|σ ′|λ mσ ′ (vh|K ·nK,σ ′) = ∑
σ∈E inth
|σ |φσ (xσ )(ũmh |Lσ ·nKσ ,σ − ũ
h |Kσ ·nKσ ,σ ) = 0,
it follows from (5.6) that (ũmh ,v) = (f
h −κmh ∇h p̃mh ,v) for any v ∈ RT0. It means that
ũmh = Π̃RT0(f
h −κmh ∇h p̃mh ) = Π̃RT0(f
h −κmh ∇h pmh ) = umh .
Thus umh = ũ
h satisfies (5.7) and divu
h = sh.
6. Stability analysis
We first check that a maximum principle holds. It ensures that the computed concentration and temper-
ature are physically relevant.
PROPOSITION 6.1 For any m ∈ {0, . . . ,N} we have 0 6 cmh 6 1 and θ− 6 θ
h 6 θ+.
PROOF. We prove the result by induction. Since c0h = ΠP0c0 and θ
h = ΠP0θ0 the result holds for m = 0
thanks to (2.11) and (3.4). Let us assume that it is true for m = n ∈ {0, . . . ,N−1}. Let K ∈Th. Equation
(5.2) implies
(1+ k sK + k λ )cn+1K = c
K + k s
K + k Dc ∑
σ∈EK∩E inth
τσ (cn+1Lσ − c
K )− k b̃h(u
We consider the last term of this relation. Since for any σ ∈ EK ∩E inth we have
(unh ·nK,σ )+σ cn+1K +(u
h ·nK,σ )−σ cn+1Lσ = (u
h ·nK,σ )σ cn+1K +(−u
h ·nK,σ )+σ (cn+1K − c
We deduce from (3.10)
−b̃h(unh,cn+1h )
−cn+1K ∑
σ∈EK∩E inth
|σ |(unh ·nK,σ )σ + ∑
σ∈EK∩E inth
(−unh ·nK,σ )+σ (cn+1Lσ − c
Since unh ∈ RT0, (unh ·nK,σ )σ = 0 for any σ ∈ EK ∩E exth . It implies that ∑σ∈EK∩E inth |σ |(u
h ·nK,σ )σ =
∑σ∈EK |σ |(u
h ·nK,σ )σ . Thus using the divergence formula and proposition 5.1 we obtain
|K| ∑
σ∈EK∩E inth
|σ |(unh ·nK,σ )σ =
divunh dx = sK .
Study of a finite volume-finite element scheme for a nuclear transport model 11 of 26
Therefore we get
(1+ 2k sK + k λ )cn+1K = c
K + k s
K + k Dc ∑
σ∈EK∩E inth
τσ (cn+1Lσ − c
|K| ∑
σ∈EK∩E inth
(−unh ·nK,σ )+σ (cn+1Lσ − c
K ). (6.1)
We consider Ki ∈ Th such that cn+1Ki = minK∈Th c
K . According to hypothesis (2.11) and definition
(3.4) we have 2sKi +λ > 0 and s
> 0. Thus, using the induction hypothesis, we deduce from (6.1)
cn+1K = c
cnKi + ks
1+ 2k sKi + k λ
k scKi
1+ 2k sKi + k λ
We now consider Ks ∈ Th such that cn+1Ks = maxK∈Th c
K . Using again hypothesis (2.11) we have
2sKs +λ > s
> 0. Thus, using the induction hypothesis, we deduce from (6.1)
cn+1K = c
cnKs + k s
1+ 2k sKs + k λ
1+ k scKs
1+ 2k sKs + k λ
A similar work for equation (5.3) proves that θ− 6 minK∈Th θ
K and maxK∈Th θ
K 6 θ+. Thus the
induction hypothesis still holds for m = n+ 1.
We now state the stability of the scheme (5.2)-(5.5).
PROPOSITION 6.2 For any 1 6 m 6 N we have
‖cnh‖2h + k
‖θ nh ‖2h 6C , (6.2)
|umh |+ |∇h pmh |6C. (6.3)
PROOF. Let 0 6 n 6 N − 1. Multiplying (5.2) by 2k cn+1h we get
(cn+1h − c
h )− 2k Dc (∆hc
h )+ k
(sh +λ )cn+1h ,c
+ k bh(u
h ) = k (s
We have (cn+1h − c
h ) = |c
2 −|cnh|
2 + |cn+1h − c
2. Thanks to propositions 4.2 and 4.5
−2k (∆hcn+1h ,c
h ) = 2k‖c
h , bh(u
h )> 0.
Using the Cauchy-Schwarz and Young inequalities we write
k (sch,c
h )6 k |s
h| |cn+1h |6C k |c
h |6 k
|cn+1h |
2 +C k.
Finally thanks to (2.11) and (3.4) we have sh > 0. Thus we obtain
|cn+1h |
2 −|cnh|
2 + 2k Dc‖cn+1h ‖
h + k
|cn+1h |
6C k.
12 of 26 C. CHOQUET AND S. ZIMMERMANN
Let m ∈ {1, . . . ,N}. Summing up the latter relation from n = 0 to m− 1 we get
|cmh |2 + 2k Dc
‖cnh‖2h 6 |c0h|
k 6C ,
thanks to proposition 6.1. With a similar work on equation (5.3), we get (6.2). We now prove (6.3). Let
n = m− 1 ∈ {0, . . . ,N − 1}. Multiplying equation (5.4) by −pn+1h and using proposition 4.1, we get
(κn+1h ∇h p
h ,∇h p
h ) = (f
h ,∇h p
h )+ (Π̃Pnc1 sh, p
h ). (6.4)
The left-hand side term satisfies (κn+1h ∇h p
h ,∇h p
h ) > κin f |∇h p
2. We now consider the right-
hand side. Using (5.1), the Cauchy-Schwarz and Young inequalities we write
|(fn+1h ,∇h p
h )|6 |f
h | |∇h p
κin f
|∇h pn+1h |
2 +C‖f‖2L∞(0,T ;L2).
Also, the stability of Π̃Pnc1 for the L
2-norm, proposition 3.1 and the Young inequality lead to
|(Π̃Pnc1 sh, p
h )|6 |sh| |p
h |6C |p
h |6C |∇h p
κin f
|∇h pn+1h |
2 +C.
Thus we deduce from (6.4) that |∇h pn+1h |
2 = |∇h pmh |
6C. Then (5.5) imply
|umh |= |un+1h |6 |f
h |+ |κ
h ∇h p
h |6 ‖f‖L∞(0,T ;L2)+ ‖κ‖W1,∞((0,1)×(0,∞)) |∇h p
h |6C.
Estimate (6.3) is proven.
7. Convergence analysis
Let ε = max(h,k). In this section we study the behavior of the scheme (5.2)-(5.5) as ε → 0. We first
define the applications cε : R→ P0, c̃ε : R→ P0, θε : R → P0, pε : R→ Pnc1 , sε : R→ P0, scε : R→ P0
and uε : R→ RT0, fε : R→ P0 by setting for all n ∈ {0, . . . ,N − 1} and t ∈ [tn, tn+1]
cε(t) = c
h , c̃ε(t) = c
(t − tn)(cn+1h − c
h), θε(t) = θ
pε(t) = p
h , sε(t) = sh, s
ε (t) = s
h, uε(t) = u
h, fε(t) = f
and for all t 6∈ (0,T )
cε(t) = c̃ε(t) = θε(t) = pε(t) = sε (t) = scε (t) = 0, uε(t) = fε(t) = 0.
We recall that the Fourier transform f̂ of a function f ∈ L1(R) is defined for any τ ∈ R by
f̂ (τ) =
e−2iπτt f (t)dt. (7.1)
We begin with the following estimate.
PROPOSITION 7.1 Let 0 < γ < 14 . There exists C > 0 such that for all ε > 0
|τ|2γ (|ĉε (τ)|2 + |θ̂ε(τ)|2)dτ 6C.
Study of a finite volume-finite element scheme for a nuclear transport model 13 of 26
PROOF. Since equations (5.2) and (5.3) are similar we only prove the estimate on ĉε . We first define
gε : R→ P0 ∩L20 as the solution of
∆hgε = Dc ∆hcε + scε − (sε +λ )cε − b̃h(uε ,cε ).
Multiplying this equation by −gε we obtain
−(∆hgε ,gε) =−Dc
∆hcε ,gε
scε − (sε +λ )cε ,gε
+ bh(uε ,cε ,gε). (7.2)
Proposition 4.5 allows us to write
−(∆hgε ,gε) = ‖gε‖2h, −(∆hcε ,gε)6 ‖cε‖h ‖gε‖h.
Thanks to the Cauchy-Schwarz inequality, (5.1) and proposition 3.2 we have
∣∣(scε − (sε +λ )cε ,gε
)∣∣6C (|sc|+ |s|+λ ) |gε |6C‖gε‖h.
According to proposition 4.3, then proposition 6.1 and (6.3), we have
uε ,cε ,gε
)∣∣6C‖cε‖∞‖gε‖h |divuε |+C‖cε‖h ‖gε‖h |uε |
6C‖gε‖h |divuε |+C‖cε‖h ‖gε‖h.
Let us plug these estimates into (7.2) and integrate from 0 to T . We get
‖gε‖h dt 6C
|divuε |dt +C
‖cε‖h dt 6C ,
because of proposition 5.1 and (6.2). Definition (7.1) then leads to
∀τ ∈ R , ‖ĝε(τ)‖h 6C. (7.3)
We now use this estimate to prove (7.1). Equation (5.2) reads
c̃ε = ∆hgε +(c0hδ0 − cNh δT )
where δ0 and δT are Dirac distributions respectively localized in 0 and T . Let τ ∈ R. Applying the
Fourier transform to the latter equation we obtain
−2iπτ ̂̃cε (τ) = ∆hĝε(τ)+ (c0h − cNh e−2iπτT ).
Let us take the scalar product of this relation with isign(τ)̂̃cε(τ). Applying propositions 3.2 and 4.5
leads to
2π |τ| |̂c̃ε(τ)|2 6C
‖ĝε(τ)‖h + |c0h|+ |c
‖̂̃cε(τ)‖h.
We assume that τ 6= 0 and multiply this estimate by |τ|2γ−1. Using proposition 6.1 and (7.3) we get
|τ|2γ |̂c̃ε(τ)|2 6C |τ|2γ−1 ‖̂̃cε(τ)‖h.
Using the Young inequality and integrating over {τ ∈ R ; |τ|> 1}, we obtain
|τ|>1
|τ|2γ |̂c̃ε(τ)|2 dτ 6
|τ|>1
|τ|4γ−2 dτ +C
|τ|>1
‖̂̃cε(τ)‖2h dτ.
14 of 26 C. CHOQUET AND S. ZIMMERMANN
For |τ|6 1, we have |τ|2γ |̂c̃ε (τ)|2 6 |̂c̃ε(τ)|2 6C‖̂̃cε(τ)‖2h according to proposition 3.2. Thus
|τ|61
|τ|2γ |̂c̃ε(τ)|2 dτ 6C
|τ|61
‖̂̃cε(τ)‖2h dτ.
By combining the bounds for |τ|> 1 and |τ|6 1 we get
|τ|2γ |̂c̃ε(τ)|2 dτ 6
|τ|>1
|τ|4γ−2 dτ +C
‖̂̃cε(τ)‖2h dτ.
Since 4γ − 2 <−1, we have
|τ|>1 |τ|4γ−2 dτ 6C. Thanks to the Parseval theorem and (6.2)
‖̂̃cε(τ)‖2h dτ 6
‖c̃ε‖2h dt 6C
k‖c0h‖
h + k
‖cnh‖2h
because ‖c0h‖h = ‖ΠP0c0‖h 6C‖c0‖1 (see [13] p. 776). Hence the result.
We can now prove the following convergence result.
PROPOSITION 7.2 There exists a subsequence of (cε ,θε , pε ,uε)ε>0, not relabeled for convenience, such
that the following convergences hold for ε → 0
cε → c in L2(0,T ;L2), θε → θ in L2(0,T ;L2), (7.4)
pε ⇀ p weakly in L
2(0,T ;H1), uε ⇀ u weakly in L
2(0,T ;L2). (7.5)
The limits (c,θ , p,u) satisfy the following properties. We have c ∈ L2(0,T ;H1), θ ∈ L2(0,T ;H1),
p ∈ L∞(0,T ;H1) and u ∈ L∞(0,T ;L2). We also have 0 6 c(x, t) 6 1 and θ− 6 θ (x, t) 6 θ+ a.e. in
Ω × [0,T ]. For all φ ∈ C ∞0 (Ω × (−1,T )), c and θ satisfy
(c,∂tφ)+Dc (∇c,∇φ)− c(u ·∇φ)−
sc − (s+λ )c
c0,φ(·,0)
, (7.6)
(θ ,∂tφ)+Dθ (∇θ ,∇φ)−θ (u ·∇φ)+
sθ + s(θ −θ∗)
θ0,φ(·,0)
. (7.7)
Lastly we have
u = f−κ(c,θ )∇p in L∞(0,T ;L2) , divu = s in L2. (7.8)
PROOF. In what follows, the convergence results hold for extracted subsequences. They are not rela-
beled for convenience. We begin by proving (7.4). According to proposition 6.1, the sequence (cε)ε>0
is uniformly bounded in L∞(0,T ;L2). Thus there exists c ∈ L∞(0,T ;L2) such that
cε ⇀ c weakly in L
2(0,T ;L2).
Using the Fourier transform, we prove that this convergence is strong. Let dε = cε − c and M > 0. We
use the following splitting
|d̂ε(τ)|2 dτ =
|τ|>M
|d̂ε(τ)|2 dτ +
|τ|6M
|d̂ε(τ)|2 dτ = IMε + JMε . (7.9)
Study of a finite volume-finite element scheme for a nuclear transport model 15 of 26
Since |d̂ε(τ)|2 6 2|ĉε(τ)|2 + 2|ĉ(τ)|2 we have
IMε 6 2
|τ|>M
|ĉε(τ)|2 dτ + 2
|τ|>M
|ĉ(τ)|2 dτ.
Using proposition 7.1 we write
|τ|>M
|ĉε(τ)|2 dτ 6
|τ|>M
|τ|2γ |ĉε(τ)|2 dτ 6
Hence
IMε 6
|τ|>M
|ĉ(τ)|2 dτ.
This implies that for all ε > 0, IMε → 0 when M → ∞. We now consider JMε . Let τ ∈ [−M,M]. Since
cε(t) ∈ P0 for all t ∈ R, and cε ⇀ c weakly in L2(0,T ;L2), we deduce from (7.1) that ĉε(τ) ∈ P0 and
ĉε(τ)⇀ ĉ(τ) weakly in L2. Extanding ĉε(τ) by 0 outside Ω , one checks ([13] p.811) that
∀η ∈ R2 , |ĉε(τ)(·+η)− ĉε(τ)|6C‖ĉε(τ)‖h |η |(|η |+ h). (7.9)
Then, using estimate (6.2), we deduce from [13] (p.834) that ĉε(τ)→ ĉ(τ) strongly in L2. Thus d̂ε(τ) =
ĉε(τ)− ĉ(τ) → 0 in L2, so that JεM → 0 when ε → 0. Now, let us report the limits for IεM and JεM into
(7.9). Using the Parseval identity we get
|d̂ε(τ)|2 dτ =
|dε |2 dt =
|cε − c|2 dt → 0.
Thus we have proven that cε → c in L2(0,T ;L2). A similar work proves that θε → θ in L2(0,T ;L2)
with θ ∈ L∞(0,T ;L2). Hence (7.4) is proven. Moreover using proposition 6.1 we obtain 0 6 c(x, t)6 1
and θ− 6 θ (x, t) 6 θ+ a.e. in Ω × [0,T ]. Lastly, using (6.2) and (7.9), we get as in [13] (p.811) that
c ∈ L2(0,T ;H1) and θ ∈ L2(0,T ;H1).
Let us now consider the sequences (pε)ε>0 and (uε )ε>0. According to (3.3) and (6.3) the sequence
(ΠPc1 pε)ε>0 is bounded in L
∞(0,T ;H1). It implies that there exists p ∈ L∞(0,T ;H1) such that ΠPc1 pε ⇀
p weakly in L2(0,T ;H1). Using proposition 3.3 we get pε ⇀ p weakly in L
2(0,T ;H1). Moreover,
according to (6.3), the sequence (uε)ε>0 is bounded in L∞(0,T ;L2) . Thus we have uε ⇀ u weakly
in L2(0,T ;L2) with u ∈ L∞(0,T ;L2). We check the properties of u. Using a Taylor expansion, the
Cauchy-Schwarz inequality, and a density argument, we have
‖κ(c,θ )−κ(cε,θε )‖L2(0,T ;L2) 6 ‖κ‖W1,∞((0,1)×(0,∞)) (‖c− cε‖L2(0,T ;L2)+ ‖θ −θε‖L2(0,T ;L2)).
Thus, using the strong convergence of the sequences (cε)ε>0 and (θε )ε>0, we have κ(cε ,θε )→ κ(c,θ )
in L2(0,T ;L2). Since ∇h pε ⇀ ∇p weakly in L2(0,T ;L2), we deduce from this
κ(cε ,θε )∇h pε ⇀ κ(c,θ )∇p weakly in L2(0,T ;L2). (7.10)
Now let v ∈ L2(0,T ;(C ∞0 )2). According to (5.5) we have
(uε ,ΠRT0v) = (fε −κ(cε ,θε)∇h pε ,v).
Using proposition 3.3 one checks easily that (fε ,ΠRT0v)→ (f,v) and (uε ,ΠRT0v)→ (u,v) in L1(0,T ).
Using moreover convergence (7.10) and a density argument, we deduce from this that u= f−κ(c,θ )∇p.
And since divuε = sε by proposition 5.1, we also have divu = s.
16 of 26 C. CHOQUET AND S. ZIMMERMANN
We finally prove that c satisfies (7.6). For all t ∈ (0,T ) equation (5.2) reads
c̃ε −Dc ∆hcε + b̃h(uε ,cε) = scε − (sε +λ )cε .
Let ψ ∈ C ∞0 (Ω × (−1,T )) and ψh = Π̃P0ψ . Multiplying the latter equation by ψh and integrating over
[0,T ], we obtain
c̃ε ,ψh
dt −Dc
(∆hcε ,ψh)dt +
bh(uε ,cε ,ψh)dt =
(scε − (sε +λ )cε ,ψh) dt. (7.11)
We now pass to the limit ε → 0 in this equation. We begin with the term
0 bh(uε ,cε ,ψh)dt. We use
the splitting b(u,c,ψ)− bh(uε ,cε ,ψh) = Aε1 +Aε2 +Aε3 with
Aε1 = b(u,c,ψ)− b(uε ,c,ψ), Aε2 = b(uε ,c,ψ)−
div(cuε)ψh dx ,
Aε3 =
div(cuε)ψh dx− bh(uε ,cε ,ψh).
According to definition (3.11)
Aε1 = b(u,c,ψ)− b(uε ,c,ψ) =−
c(u−uε) ·∇ψ dx.
We know that c∇ψ ∈ L2(0,T ;L2). Since uε ⇀ u in L2(0,T ;L2) we get
1 dt → 0. We now consider
Aε2. We have
Aε2 =
(ψ −ψh)div(cuε)dx =
(ψ −ψh)(uε ·∇c+ cdivuε)dx.
Using the Cauchy-Schwarz inequality we get
|Aε2|dt 6 ‖ψ −ψh‖L∞(Ω×(0,T)) (‖uε‖L2(0,T ;L2)+ ‖divuε‖L2(0,T ;L2))‖c‖L2(0,T ;H1).
Using a Taylor expansion, one checks that ‖ψ −ψh‖L∞(Ω×(0,T)) 6 h‖∇ψ‖L∞(Ω×(0,T )). Thus
2 dt →
0. Finally we estimate Aε3. For all triangles K ∈ Th and L ∈ Th sharing an edge σ , we set cK,L = cK if
uε ·nK,σ > 0 and cK,L = cL otherwise. Using the divergence formula, we deduce from definition (3.12)
Aε3 = ∑
σ∈EK∩E inth
(c− cK,Lσ )(uε ·nK,σ )dσ
σ∈E inth
(ψKσ −ψLσ )
(c− cKσ ,Lσ )(uε ·nKσ ,σ )dσ .
Using definition (3.5) this reads
Aε3 = ∑
σ∈E inth
(ψKσ −ψLσ )
(ΠPnc1 c− cKσ ,Lσ )(uε ·nKσ ,σ )dσ
σ∈E inth
(ψKσ −ψLσ ) |σ |
(ΠPnc1 c)(xσ )− cKσ ,Lσ
(uε ·nKσ ,σ )σ .
Study of a finite volume-finite element scheme for a nuclear transport model 17 of 26
Using a Taylor expansion, one checks that |ψKσ −ψLσ |6 h‖∇ψ‖L∞(Ω×(0,T )). Moreover |σ |6 h. Thus,
using the Cauchy-Schwarz inequality, we have
|Aε3| 6 C h2 ∑
σ∈E inth
|uε (xσ )|
∣∣∣(ΠPnc1 c)(xσ )− cKσ ,Lσ
6 C h2
σ∈E inth
|uε (xσ )|2
)1/2( ∑
σ∈E inth
|(ΠPnc1 c)(xσ )− cKσ ,Lσ |
Using the assumption on the mesh, one checks that |K| > C h2 for all K ∈ Th. Thus, thanks to a
quadrature formula, we have
|Aε3| 6 C
σ∈EK∩E inth
|uε (xσ )|2
)1/2(
σ∈EK∩E inth
|(ΠPnc1 c)(xσ )− cKσ |
6 C |uε | |ΠPnc1 c− cε |.
We write ΠPnc1 c− cε = (ΠPnc1 c− c)+ (c− cε) and we use proposition 3.3 . We obtain with (6.3)
|Aε3|dt 6C‖uε‖L∞(0,T ;L2) (h‖c‖L2(0,T ;H1)+ ‖c− cε‖L2(0,T ;L2)).
Since cε → c in L2(0,T ;L2) when ε = max(h,k) → 0, we conclude that
3 dt → 0. Gathering the
limits for Aε1, A
3, we obtain
bh(uε ,cε ,Π̃P0ψ)dt →
b(u,c,ψ)dt.
We now consider the other terms in (7.11). Proposition 4.6 leads to
(∆hcε ,Π̃P0ψ) =
cε ,∆h(Π̃P0ψ)
cε ,∆h(Π̃P0 ψ)−∆ψ
+(cε ,∆ψ). (7.12)
According to proposition 4.7
cε ,∆h(Π̃P0ψ)−∆ψ
)∣∣∣ 6 ‖cε‖h ‖∆h(Π̃P0ψ)−∆ψ‖−1,h 6C h‖cε‖h‖ψ‖2.
We then apply the Cauchy-Schwarz inequality and use (6.2). We obtain
cε ,∆h(Π̃P0 ψ)−∆ψ
)∣∣∣ dt 6C h
‖cε‖2h dt
‖cnh‖
6C h.
Moreover, since cε → c in L2(0,T ;L2), we have
0 (cε ,∆ψ)dt →
0 (c,∆ψ)dt. Thus we deduce from
(7.12) ∫ T
(∆hcε ,Π̃P0ψ)dt →
(c,∆ψ)dt.
We are left with two terms. First, using Taylor expansions, one checks that
ψh → ψ , ∂tψh → ∂tψ in L2(Ω × (−1,T)) , ψh(·,0)→ ψ(·,0) in L2. (7.13)
18 of 26 C. CHOQUET AND S. ZIMMERMANN
We know that cε → c in L2(0,T ;L2). Thus
(sc − (s+λ )cε ,ψh) dt →
(sc − (s+λ )c,ψ) dt.
Finally, integrating by parts the first term of (7.11), we get
c̃ε ,ψh
dt = (c̃ε ,ψh)t=T − (c̃ε ,ψh)t=0 −
(c̃ε ,∂tψh)dt.
Since ψ ∈ C ∞0 (Ω × (−1, ,T )) we have (c̃ε ,ψh)t=T = 0. Using proposition 3.3 one checks that c0h =
ΠP0c0 → c0 in L2; using moreover (7.13) we get
(c̃ε ,ψh)t=0 =
c0h,ψh(·,0)
ΠP0c0,ψh(·,0)
c0,ψ(·,0)
For the last term, one easily checks that ‖c̃ε − cε‖L2(0,T ;L2) → 0. Thus, since cε → c in L2(0,T ;L2),
we also have c̃ε → c in L2(0,T ;L2). Using moreover (7.13) we get
0 (c̃ε ,∂t ψh)dt →
0 (c,∂tψ)dt.
Therefore ∫ T
c̃ε ,ψh
dt →−
c0,ψ(·,0)
(c,∂tψ)dt.
By gathering the limits we have obtained in (7.11) we get (7.6). The relation (7.7) for θ is proven in a
similar way.
8. Error estimates
We have proven in section 7 that the problem (2.5)–(2.9) has a weak solution (c,θ , p,u). From now on,
we assume the following regularity for this solution:
c,θ ∈ C (0,T ;H2) , ct ,θt ∈ L2(0,T ;H1+r)∩C (0,T ;L2) ,
ctt ,θtt ∈ L2(0,T ;L2) , p ∈ C (0,T ;H2) , u ∈ C (0,T ;H1+s) ,
with r > 0 and s > 0. We also assume that f ∈ C (0,T ;H1).
8.1 Definitions
For all m ∈ {0, . . . ,N}, we define the following errors
emh,c = c(tm)− cmh , emh,θ = c(tm)−θ mh ,
emh,p = p(tm)− pmh , emh,u = u(tm)−umh .
We have the following splittings
emh,c = ε
h,c +η
h,c, e
h,θ = ε
h,θ +η
h,θ ,
emh,p = ε
h,p +η
h,p, e
h,u = ε
h,u+η
with the discrete errors
εmh,c = Π̃P0c(tm)− c
h , ε
h,θ = Π̃P0θ (tm)−θ
εmh,p = ΠPnc1 p(tm)− p
h , ε
h,u= ΠRT0u(tm)−u
Study of a finite volume-finite element scheme for a nuclear transport model 19 of 26
and the interpolation errors
ηmh,c = c(tm)− Π̃P0c(tm), η
h,θ = θ (tm)− Π̃P0θ (tm),
ηmh,p = p(tm)−ΠPnc1 p(tm), η
h,u= u(tm)−ΠRT0u(tm).
The interpolation errors are estimated as follows. We write |ηmh,c| 6 |c(tm)−ΠP0c(tm)|+ |ΠP0c(tm)−
Π̃P0c(tm)| and the same for ηmh,θ . Using proposition 3.3 and (3.7) we obtain
|ηmh,c|6C h‖c(tm)‖1 6C h‖c‖L∞(0,T ;H1) , |ηmh,θ |6C h‖θ‖L∞(0,T ;H1). (8.1)
According to proposition 3.3 we also have
|ηmh,p|+ |∇̃hηmh,p|6C h‖p(tm)‖2 6C h‖p‖L∞(0,T ;H2), (8.2)
|ηmh,u|6C h‖u(tm)‖1 6C h‖u‖L∞(0,T ;H1). (8.3)
We now have to estimate the discrete errors.
PROPOSITION 8.1 For all n ∈ {0, . . . ,N − 1} and ψh ∈ Pnc1 we have
εn+1h,c − ε
−Dc ∆hεn+1h,c + b̃h
εnh,u,Π̃P0c(tn+1)
+ b̃h(u
h,c )+ (s
h +λ )ε
h,c =C
h,1 +C
h,2 , (8.4)
εn+1h,θ − ε
−Dθ ∆hεn+1h,θ + b̃h
εnh,u,Π̃P0θ (tn+1)
+ b̃h(u
h,θ )+ sh ε
h,θ =Θ
h,1 +Θ
h,2 , (8.5)
κ(cn+1h ,θ
h )∇hε
h,p ,∇hψh
(κn+1h,1 ε
h,c +κ
h,2 ε
h,θ )∇p(tn+1),∇hψh
Pn+1h ,∇hψh
, (8.6)
εn+1h,u =−Π̃RT0
(κn+1h,1 ε
h,c +κ
h,2 ε
h,θ )∇p(tn+1)+κ(c
h )∇hε
−Un+1h . (8.7)
For all m ∈ {0, . . . ,N}, the consistency errors Cmh,1, C
h,2, Θ
h,1, Θ
h,2, P
h and U
h are defined in (8.9),
(8.10), (8.14), (8.15) and the terms κmh,1 and κ
h,2 are given by (8.13) below.
PROOF. Let n ∈ {0, . . . ,N − 1}. Equation (2.5) for t = tn+1 reads
∂tc(tn+1)−Dc ∆c(tn+1)+ b̃
u(tn+1),c(tn+1)
= sc − (s+λ )c(tn+1).
We introduce the time discretization by setting
Rn+1 =
(c(tn+1)− c(tn)
− ct(tn+1)
u(tn)−u(tn+1),c(tn+1)
We get
c(tn+1)− c(tn)
−Dc ∆c(tn+1)+ b̃
u(tn),c(tn+1)
= sc − (s+λ )c(tn+1)+Rn+1.
We apply ΠP0 to this equation. By subtracting the result from (5.2) we get
(c(tn+1)− c(tn)
cn+1h − c
ΠP0∆c(tn+1)−∆hc
+ΠP0 b̃
u(tn),c(tn+1)
− b̃h(unh,cn+1h )+ΠP0
(s+λ )c(tn+1)
− (sh +λ )cn+1h = ΠP0R
n+1. (8.8)
20 of 26 C. CHOQUET AND S. ZIMMERMANN
We now introduce the discrete errors as follows. Since c(tn+1)− c(tn) =
∫ tn+1
tn ct(s)ds one checks that
(c(tn+1)− c(tn)
cn+1h − c
∫ tn+1
ΠP0ct(s)− Π̃P0ct(s)
(εn+1h,c − ε
h,c).
We also have
ΠP0∆c(tn+1)−∆hc
h = ΠP0∆c(tn+1)−∆h
Π̃P0c(tn+1)
+∆hεn+1h,c .
Using the linearity of b̃h, one easily checks that
ΠP0 b̃
u(tn),c(tn+1)
− b̃h(unh,cn+1h ) = b̃h(u
h,c )+ b̃h
εnh,u,Π̃P0c(tn+1)
+ΠP0 b̃
u(tn),c(tn+1)
− b̃h
ΠRT0u(tn),Π̃P0c(tn+1)
Lastly
(s+λ )c(tn+1)
− (sh +λ )cn+1h = ΠP0
(s+λ )ηn+1h,c
+(sh +λ )εn+1h,c .
Using these relations in (8.8) we get (8.4). For any m ∈ {1, . . . ,N}, the consistency errors Cmh,1 ∈ P0 and
Cmh,2 ∈ P0 are given by
Cmh,1 = ΠP0
(c(tm)− c(tm−1)
− ct(tm)+ b̃
u(tm−1)−u(tm),c(tm)
+ ΠP0
(s+λ )ηmh,c
Π̃P0ct(s)−ΠP0ct(s)
ds, (8.9)
Cmh,2 = Dc
ΠP0∆c(tm)−∆h
Π̃P0c(tm)
ΠP0 b̃
u(tm−1),c(tm)
− b̃h(ΠRT0u(tm−1),Π̃P0c(tm)
A similar proof leads to (8.5) where the consistence errors Θ mh,1 ∈ P0 and Θ mh,2 ∈ P0 are defined for any
m ∈ {1, . . . ,N} by
Θ mh,1 = ΠP0
(θ (tm)−θ (tm−1)
−θt(tm)+ b̃
u(tm−1)−u(tm),θ (tm)
+ ΠP0(sη
h,θ )−
Π̃P0θt(s)−ΠP0θt(s)
ds, (8.10)
Θ mh,2 = Dθ
ΠP0∆θ (tm)−∆h
Π̃P0θ (tm)
ΠP0 b̃
u(tm−1),θ (tm)
− b̃h
ΠRT0u(tm−1),Π̃P0θ (tm)
We now consider the problem associated with the pressure. Let n ∈ {0, . . . ,N − 1} and ψh ∈ Pnc1 . Mul-
tiplying equation (2.7) written for t = tn+1 by ψh and integrating by parts, we get
κ(c(tn+1),θ (tn+1))∇p(tn+1),∇hψh
= (f(tn+1),∇hψh)+ (s,ψh). (8.11)
On the other hand, using (5.4) and proposition 4.1, we have
κ(cn+1h ,θ
h )∇h p
h ,∇hψh
= (fn+1h ,∇hψh)+ (Π̃Pnc1 sh,ψh).
Since ∇hψh ∈ P0, one checks that (fn+1h ,∇hψh) = (ΠP0f(tn+1),∇hψh) = (f(tn+1),∇hψh). According to
(3.5) we also have (Π̃Pnc1 sh,ψh) = (sh,ψh). Thus
κ(cn+1h ,θ
h ) ∇h p
h ,∇hψh
= (f(tn+1),∇hψh)− (sh,ψh).
Study of a finite volume-finite element scheme for a nuclear transport model 21 of 26
Substracting (8.11) from the latter relation, we obtain
κ(c(tn+1),θ (tn+1))∇p(tn+1)−κ(cn+1h ,θ
h )∇h p
h ,∇hψh
=−(s− sh,ψh). (8.12)
We split the left-hand side as
κ(cn+1h ,θ
h )(∇p(tn+1)−∇h p
κ(c(tn+1),θ (tn+1))−κ(cn+1h ,θ
∇p(tn+1).
Using a Taylor expansion, one can check that
κ(c(tn+1),θ (tn+1))−κ(cn+1h ,θ
h ) = (ε
h,c +η
h,c )κ
h,1 +(ε
h,θ +η
h,θ )κ
h,2 .
We have set for any m ∈ {0, . . . ,N} and s ∈ [0,1]
cmh (s) = c
h +(c(tm)− c
h )s, θ
h (s) = θ
h +(θ (tm)−θ
h )s,
κmh,1 =
κx(cmh (s),θ
h (s))ds, κ
h,2 =
κy(cmh (s),θ
h (s))ds.
(8.13)
We also have
∇p(tn+1)−∇hpn+1h = ∇hε
h,p + ∇̃hη
h,p .
Plugging these relations into (8.12) we get (8.6). For all m ∈ {0, . . . ,N} we have
Pmh = (κ
h,c +κ
h,θ )∇p(tm)+κ(c
h )∇hη
h,p. (8.14)
We end with the equation associated with u. Let n ∈ {0, . . . ,N − 1}. Applying the operator Π̃RT0 to
(2.7) for t = tn+1 we obtain
Π̃RT0u(tn+1) = Π̃RT0f(tn+1)− Π̃RT0
κ(c(tn+1),θ (tn+1))∇p(tn+1)
Let us substract this equation from (5.5). Since fn+1h = ΠP0f(tn+1) we get
Π̃RT0u(tn+1)−u
h = Π̃RT0
f(tn+1)−ΠP0f(tn+1)
− Π̃RT0
c(tn+1),θ (tn+1)
∇p(tn+1)−κn+1h ∇h p
One easily checks that
Π̃RT0u(tn+1)−u
h = Π̃RT0(u(tn+1)−ΠRT0u(tn+1))+ ε
h,u .
Thus we get (8.7). For all m ∈ {0, . . . ,N}, we have
Umh = Π̃RT0
(f(tm)−ΠP0f(tm))−ηmh,u−P
. (8.15)
This ends the proof of proposition 8.1.
22 of 26 C. CHOQUET AND S. ZIMMERMANN
8.2 Error estimates
We first estimate the consistency errors.
PROPOSITION 8.2 For all m ∈ {1, . . . ,N} the consistency errors satisfy
|Cnh,1|
2 + k
|Θ nh,1|
6C (h2 + k2), (8.16)
‖Cnh,2‖2−1,h + k
‖Θ nh,2‖2−1,h 6C h2 , (8.17)
|Pmh |+ |Umh |6C h. (8.18)
PROOF. Let n ∈ {1, . . . ,N}. Since the operator ΠP0 is stable for the L2-norm we have
|ΠP0R
n|6 |Rn|6
c(tn)− c(tn−1)
− ct(tn)
∣∣∣b̃
u(tn−1)−u(tn),c(tn)
)∣∣∣ .
Using a Taylor expansion and the Cauchy-Schwarz inequality, we get
c(tn)− c(tn−1)
− ct(tn)
∫ tn−1
|tn−1 − s| |ctt(s)|ds 6
(∫ tn−1
|ctt(s)|2 ds
On the other hand, since ∇c(tn)|∂Ω = 0, we deduce from (3.10) by integrating by parts
b̃(u(tn−1)−u(tn),c(tn)) = (u(tn−1)−u(tn)) ·∇c(tn).
Using a Taylor expansion and the Cauchy-Schwarz inequality, we get
|(u(tn−1)−u(tn)) ·∇c(tn)|6
k‖c‖L∞(0,T ;H1)
(∫ tm−1
|ut(s)|2 ds
|ΠP0R
(∫ tm−1
|ctt(s)|2 ds
k‖c‖L∞(0,T ;H1)
(∫ tm−1
|ut(s)|2 ds
Thanks to the stability of ΠP0 for the L
2-norm and to (8.1) we have
∣∣ΠP0
(s+λ )ηnh,c
)∣∣6C h(‖s‖L∞(0,T ;L2)+λ )‖c‖L∞(0,T ;H1).
The Cauchy-Schwarz inequality and (3.7) allow us to write
∣∣ΠP0ct(s)− Π̃P0ct(s)
∣∣ds 6C h
(∫ tm
‖ct(s)‖21+r ds
By plugging these estimates into definition (8.9) we get
k |Cnh,1|2 6 k2 ‖c‖2L∞(0,T ;H1)
|ut(s)|2 ds+ k2
|ctt(s)|2 ds
+ C h2
‖ct(s)‖21+r ds+C k h2‖c‖2L∞(0,T ;H1).
Study of a finite volume-finite element scheme for a nuclear transport model 23 of 26
Summing up the latter relation for n = 1 to m ∈ {1, . . . ,N} and using a similar work on Θ mh,1 we get
(8.16). Now let n ∈ {1, . . . ,N}. Using propositions 4.7 and 4.4 we have
‖ΠP0∆c(tn)−∆h
Π̃P0c(tn)
‖−1,h 6C h‖c‖L∞(0,T ;H2)
‖ΠP0 b̃(u(tn−1),c(tn))− b̃h(ΠRT0u(tn−1),Π̃P0c(tn))‖−1,h 6C h‖c‖L∞(0,T ;H1)‖u‖L∞(0,T ;H1+s).
Plugging these estimates into definition (8.9) and summing up from n = 1 to m, we obtain
‖Cnh,2‖2−1,h 6C h2
‖c‖2L∞(0,T ;H1) ‖u‖
L∞(0,T ;H1+s)+ ‖c‖
L∞(0,T ;H2)
A similar work on Θ mh,2 then leads to (8.17). We finally prove (8.18). Let m ∈ {1, . . . ,N}. On the one
hand, we have by (8.14)
|Pmh |6 (‖κmh,1‖∞ |ηmh,c|+ ‖κmh,2‖∞ |ηmh,θ |) |∇p(tm)|+ ‖κ(cmh ,θ mh )‖∞|∇hηmh,p|.
Using estimates (8.1)–(8.3) we get
|Pmh |6C h‖κ‖W1,∞((0,1)×(0,∞))(‖c‖L∞(0,T ;H1)+ ‖θ‖L∞(0,T ;H1)+ ‖p‖L∞(0,T ;H2)).
On the other hand definition (8.15) leads to
|Umh |6
∣∣Π̃RT0
f(tm)−ΠP0f(tm)
)∣∣+ |ηmh,u|+ |P
Using the stability of Π̃RT0 for the L
2-norm and proposition 3.3 we have
∣∣∣Π̃RT0
f(tm)−ΠP0f(tm)
)∣∣∣6 |f(tm)−ΠP0f(tm)|6C h‖f‖L∞(0,T ;H1).
Using moreover (8.3) we obtain
|Umh |6C h‖κ‖W1,∞((0,1)×(0,∞))(‖c‖L∞(0,T ;H1)+ ‖θ‖L∞(0,T ;H1))
+C (‖p‖L∞(0,T ;H2)+ ‖f‖L∞(0,T ;H1)+ ‖u‖L∞(0,T ;H1)).
We have proven (8.18).
Using the former proposition we are now able to estimate the discrete errors.
PROPOSITION 8.3 There exists some real k0 > 0 such that for any k < k0 and m ∈ {1, . . . ,N}
|εmh,c|2 + |εmh,θ |2 + k
‖εmh,c‖2h + ‖εmh,θ‖2h
6C (h2 + k2), (8.19)
|∇hεmh,p|+ |εmh,u|6C (h+ k). (8.20)
PROOF. Multiplying (8.4) by 2k εn+1h,c , we obtain
(εn+1h,c − ε
, 2k εn+1h,c
− 2Dc k (∆hεn+1h,c ,ε
h,c )+ 2k bh
unh,ε
h,c ,ε
+ 2k λ |εn+1h,c |
= 2k (Cn+1h+1 +C
h,2 ,ε
h,c )− 2k (sh, |ε
h,c |
2)− 2k bh
εnh,u,Π̃P0c(tn+1),ε
. (8.21)
24 of 26 C. CHOQUET AND S. ZIMMERMANN
Using an algebraic identity we have
(εn+1h,c − ε
,2k εn+1h,c
= |εn+1h,c |
2 −|εnh,c|2 + |εn+1h,c − ε
h,c|2.
We know by propositions 4.2 and 4.5 that
−2k (∆hεn+1h,c ,ε
h,c ) = 2k‖ε
h,c ‖
h , bh
unh,ε
h,c ,ε
We have
2k (sh, |εn+1h,c |
2)6 2k‖sh‖∞ |εn+1h,c |
6C k |εn+1h,c |
Using the Young inequality, we also write
∣∣∣2k (Cn+1h,1 ,ε
h,c )
∣∣∣6 2k |Cn+1h,1 | |ε
h,c |6 k λ |ε
h,c |
2 +C k |Cn+1h,1 |
|2k (Cn+1h,2 ,ε
h,c )|6 2k‖C
h,2 ‖−1,h‖ε
h,c ‖h 6 Dc
‖εn+1h,c ‖
h +C k‖Cn+1h,2 ‖
−1,h.
We are left with the term bh
εnh,u,Π̃P0c(tn+1),ε
. We have εnh,u∈ RT0. Using the divergence formula,
proposition 5.1 and (3.5), one easily checks that divεnh,u= 0. Thus we can apply proposition 4.3 to get
∣∣bh(εnh,u,Π̃P0c(tn+1),ε
h,c )
∣∣6C |εnh,u|‖Π̃P0c(tn+1)‖h‖ε
h,c ‖h. (8.22)
Let us first bound ‖Π̃P0c(tn+1)‖h. We have
‖Π̃P0c(tn+1)‖h 6 ‖Π̃P0c(tn+1)−ΠP0c(tn+1)‖h + ‖ΠP0c(tn+1)‖h.
Using an inverse inequality (see proposition 1.2 in [22]) and (3.7)
‖Π̃P0c(tn+1)−ΠP0c(tn+1)‖h 6
|Π̃P0c(tn+1)−ΠP0c(tn+1)|6C.
Moreover, according to [13] (p. 776), we have ‖ΠP0c(tn+1)‖h 6C‖c‖L∞(0,T ;H1). Thus ‖Π̃P0c(tn+1)‖h 6
C. We now estimate |εnh,u|. Using the stability of Π̃RT0 for the L
2-norm and the Cauchy-Schwarz
inequality in (8.7) we get
|εnh,u|6 ‖κ‖W1,∞((0,1)×(0,∞))
‖∇p‖L∞(0,T ;L2) (|ε
h,c |+ |ε
h,θ |)+ |∇hε
h,p |
+ |Un+1h |. (8.23)
We bound εn+1h,p as follows. Setting ψh = ε
h,p in (8.6) and using the Cauchy-Schwarz inequality, we get
κ(cn+1h ,θ
h )∇hε
h,p ,∇hε
6 ‖κ‖W1,∞((0,1)×(0,∞))‖∇p‖L∞(0,T ;L2) (|εn+1h,c |+ |ε
h,θ |) |∇hε
h,p |
+ |Pn+1h | |∇hε
h,p |+ |s− sh| |ε
h,p |.
The left-hand side is such that (κ(cn+1h ,θ
h )∇hε
h,p ,∇hε
h,p )> κin f |∇hε
h,p |
2. As for the right-hand
side, we have sh = ΠP0s and ε
h,p ∈ P
1 ∩L20. Thus, according to propositions 3.1 and 3.3
|s− sh| |εn+1h,p |6C h‖s‖L∞(0,T ;H1) |∇hε
h,p |.
Study of a finite volume-finite element scheme for a nuclear transport model 25 of 26
Finally |Pn+1h |6C h thanks to (8.18). Therefore we obtain
|∇hεn+1h,p |6C (h+ |ε
h,c |+ |ε
h,θ |). (8.24)
Let us plug this estimate into (8.23). Since |Un+1h |6C h thanks to (8.18), we get
|εnh,u|6C (h+ |ε
h,c |+ |ε
h,θ |). (8.25)
Now, plugging this bound into (8.23) and using the Young inequality, we obtain
∣∣bh(εnh,u,Π̃P0c(tn+1),ε
h,c )
∣∣ 6 C k (h+ |εn+1h,c |+ |ε
h,θ |)‖ε
h,c ‖h
‖εn+1h,c ‖
h +C k (h
2 + |εn+1h,c |
2 + |εn+1h,θ |
Now we have treated all the terms in (8.21). This equation implies
|εn+1h,c |
2 −|εnh,c|2 +Dc k‖εn+1h ‖
h 6C k
h2 + |εn+1h,c |
2 + |εn+1h,θ |
2 + |Cn+1h,1 |
2 + ‖Cn+1h,2 ‖
Let m ∈ {1, . . . ,N}. Let us sum up the latter estimate from n = 0 to m− 1. Thanks to (3.3)
|ε0h,c|= |Π̃P0c0 − c
h|= |Π̃P0c0 −ΠP0c0|6C h‖c‖L∞(0,T ;H2).
Using moreover estimates (8.16) and (8.17) we get
|εmh,c|2 +Dc k
‖εnh,c‖2h 6C k
(|εnh,c|2 + |εnh,θ |2)+C (h2 + k2).
Summing this relation with the one obtained by a similar work on (8.5) we obtain
|εmh,c|2 + |εmh,θ |2 + k
(Dc ‖εnh,c‖2h +Dθ ‖εnh,θ‖2h)6C k
(|εnh,c|2 + |εnh,θ |2)+C (h2 + k2).
Using a discrete Gronwall lemma (see lemma 5.2 in [22]) we get (8.19). Then (8.24) and (8.25) imply
(8.20).
By combining proposition 8.3 with estimates (8.1)-(8.3), we obtain finally the following result.
Theorem 8.1 There exists a real k0 > 0 such that for all k < k0 and m ∈ {1, . . . ,N}
|emh,c|2 + |emh,θ |2 + k
‖Π̃P0e
h,c‖2h + ‖Π̃P0e
h,θ‖2h
6C (h2 + k2) ,
|∇̃hemh,p|+ |emh,u|6C (h+ k).
REFERENCES
[1] ACHDOU, Y., BERNARDI, C. & COQUEL, F. (2003) A priori and a posteriori analysis of finite volume discretiza-
tions of Darcy’s equations, Numer. Math., 96(1) 17–42.
[2] AFIF, M.,AMAZIANE, B. (2002) Convergence of finite volume schemes for a degenerate convection-diffusion
equation arising in flow in porous media, Comput. Methods Appl., 191 5265–5286.
26 of 26 C. CHOQUET AND S. ZIMMERMANN
[3] ANGOT, P., DOLJŠI, V., FEISTAUER, M. & FELCMAN, J. (1998) Analysis of a combined barycentric finite
volume-nonconforming finite element method for nonlinear convection-diffusion problems, Appl. Math., 43(4)
263–310.
[4] BATEMAN, H. (1910) The solution of a system of differential equations occuring in the theory of radioactive
transformations, Proc. Camb. Philos. Soc., 15 423–427.
[5] BEAR, J. (1972) Dynamics of Fluids in Porous Media, New York: American Elsevier.
[6] BRENNER, S. C., SCOTT, L. R. & RIDGWAY, L. (2002) The mathematical theory of finite element methods, Texts
in applied mathematics, 16, New-York: Springer-Verlag.
[7] BREZZI, F. & FORTIN, M. (1991) Mixed and hybrid finite element methods, Springer Series in Computational
Mathematics, 15, New-York: Springer-Verlag.
[8] CHAVENT, G., JAFFRÉ, J. & ROBERTS, J. E. (1995) Mixed hybrid finite elements and cell-centered finite volumes
for two-phase flow in porous media. In: A. P. BOURGEAT AND AL, ed. Mathematical modeling of flow
through porous media, London: World Scientific, 100–114.
[9] CHEN, Z. (2000) Formulation and numerical methods of the black oil model in porous media, SIAM J. Numer.
Anal., 38(2) 489–514.
[10] CHOQUET C. (2004) Existence result for a radionuclide transport model with unbounded viscosity , J. Math. Fluid
Mech., 6(4) 365–388.
[11] DOUGLAS J., JR. & SPAGNUOLO, A. M. (2001) The approximation of nuclear contaminant transport in porous
media, J. Korean Math. Soc., 38 723–761.
[12] EWING, R. E., YUAN, Y. & LI, G. (1989) A time-discretization procedure for a mixed finite element approxi-
mation of contamination by incompressible nuclear waste in porous media, In: Mathematics for large scale
computing, Lecture Notes in Pure and Appl. Math., 120, 127–145.
[13] EYMARD, R., GALLOUËT, T. & HERBIN, R. (2000) The finite volume method (chapter III), In P.G. CIARLET
AND J. L. LIONS, eds Handbook of numerical analysis, 7, 762–835.
[14] GEISER, J. (2002) Numerical simulation of a model for transport and reaction of radionuclides, In: Algoritmy 2002,
Vysoké Tatry - Podbanské, Slovakia, Bratislava: Slovak university of technology.
[15] JOHN, V., MAUBACH, J. M. & TOBISKA, L. (1997) Nonconforming streamline-diffusion finite-element methods
for convection-diffusion problems, Numer. Math., 78(2) 165–188.
[16] KOVAL, E. J. (1963) A method for predicting the performance of unstable miscible displacements in heterogeneous
media, SPEJ trans. AIME, 228 145–154.
[17] LIN P. & YANG, D. (1998) An iterative perturbation method for the pressure equation in the simulation of miscible
displacement in porous media, SIAM J. Sci. Comput., 19(3) 893–911.
[18] MICHEL, A. (1999) Convergence of a finite volume scheme for a nonlinear convection-diffusion problem, In:
Proceedings of the Second international symposium on finite volumes for complex applications, Duisburg.
[19] MORTON, K. W. (1996) Numerical solution of convection-diffusion problems, In: Applied Mathematics and math-
ematical Computation, 12, London: Chapman and Hall.
[20] VERWER, J. G., BLOM, J. G. & HUNDSDORFER, W. (1996) An implicit-explicit approach for atmospheric trans-
port chemistry problems, Appl. Numer. Math., 20(1-2) 191–209.
[21] ZIMMERMANN, S. (2006) Stability of a finite volume scheme for incompressible fluids, preprint.
[22] ZIMMERMANN, S. (2006) Étude et implémentation de méthodes de volumes finis pour les fluides incompressibles.
Thesis (PhD). Blaise Pascal university (France).

  We consider a problem of nuclear waste contamination. It takes into account
the thermal effects. The temperature and the contaminant's concentration
fulfill convection-diffusion-reaction equations. The velocity and the pressure
in the flow satisfy the Darcy equation, with a viscosity depending on both
concentration and temperature. The equations are nonlinear and strongly
coupled. Using both finite volume and nonconforming finite element methods, we
introduce a scheme adapted to this problem. We prove the stability and
convergence of this scheme and give some error estimates.

Introduction
A part of the high-level nuclear waste is now stored in environmentally safe locations. One has to
consider the eventuality of a leakage through the engineering and geological barriers. It may cause
the contamination of underground water sources far away from the original repository’s location. In the
present paper, we consider such a problem of nuclear waste contamination in the basement. We take into
account the thermal effects. The evolution in time of the temperature and of the contaminants concen-
tration is then governed by convection-diffusion-reaction equations. The velocity and the pressure in the
flow satisfy the Darcy equation, with a viscosity depending on both concentrations and temperature in
a nonlinear way. The velocity satisfies an incompressibility constraint. We introduce a scheme adapted
to this problem. We use both finite volume and nonconforming finite element methods. It ensures that
a maximum principle holds and that the associated linear systems have good-conditioned matrices. We
prove the stability and convergence of the scheme and give some error estimates.
Let us briefly point out some previous works. A complete model coupling concentrations and pres-
sure equations is very seldom studied, since the system is strongly coupled. Instead, each equation
is considered separately. In the general context of convection-diffusion-reaction equations, numerous
schemes are available (see [15] or [19] and the references therein). Finite difference schemes are some-
times used for the convective term (in [20] for instance). But they are not adapted for the complex
geometry of a reservoir. More recently, finite volume methods were developed and analysed. Let us just
cite the book [2] or [13] and the references therein. Finite elements (for the diffusive term) and finite
volumes (for the convective term) are coupled for instance in [3, 8] . In convection dominant problems,
the equations are of degenerate parabolic type. This setting is considered in [9, 18] . The reaction
terms are specifically studied through operator splitting methods in [14] . Now, in the specific context of
porous media flow, we mention [17] who consider only the evolution of the pressure. In [11, 12] a more
complete set of equations is used, and a mixed finite element approximation is developed. We stress that
IMA Journal of Numerical Analysis c© Institute of Mathematics and its Applications 2007; all rights reserved.
http://arxiv.org/abs/0704.1286v1
2 of 26 C. CHOQUET AND S. ZIMMERMANN
in all theses works, as in most, the thermic effects are not taken into account.
The present paper is organized as follows. Section 2 is devoted to the derivation of the model. In
section 3, we introduce the discrete tools used in this paper. It allows us to define the numerical scheme
of section 5. The analysis of the scheme uses the properties of section 4. We then prove the stability and
convergence of the scheme, in sections 6 and 7 respectively. We conclude with some error estimates in
section 8.
2. Model of contamination
The thickness of the medium is significantly smaller than its length and width. Hence it is reasonable to
average the medium properties vertically and to describe the far-field repository by a polyhedral domain
Ω of R2 with a smooth boundary ∂Ω . It is characterized by a porosity φ and a permeability tensor K.
The time interval of interest is [0,T ]. We denote by p the pressure, by (ai)
i=1 the concentrations of the
Nr radionuclides involved in the flow and by θ the temperature. The Darcy velocity is represented by
u. We assume a miscible and incompressible displacement. Due to the mass and energy conservation
principles, the flow is governed by the following system satisfied in Ω × [0,T ], with i = 1, ..,Nr −1 (see
[10]).
φ Ri ∂tai + div(ai u)− div(φ Dc ∇ai) = si − sai −λi Ri φ ai +
j=1, j 6=i
ki j λ j R j φ a j , (2.1)
φ Cp ∂tθ + div(θ u)− div(φ Dθ ∇θ ) =−sθ − s(θ −θ∗) , (2.2)
divu = s, u+
(a j)
j=1 ,θ
) ∇p = f. (2.3)
In (2.1) the retardation factors Ri > 0 are due to the sorption mechanism. The real λ−1i > 0 denotes the
half life time of radionuclide i. The term −λi Ri φ ai describes the radioactive decay of the i–th specy.
Meanwhile, the quantity ∑ j 6=i ki j λ j R j φ a j is created by radioactive filiation. The molecular diffusion
effects are given by the coefficient Dc > 0. The contamination is represented by the source term si and
s=∑Nri=1 si. In (2.2) the coefficient Cp > 0 is the relative specific heat of the porous medium. The thermic
diffusion coefficient is denoted by Dθ > 0. The real θ∗ > 0 is a reference temperature. The constitutive
relation (2.3) is the Darcy law and f is a density of body forces. For a large range of temperatures µ has
the form
µ(a,θ ) = µR(a)exp
where µR is a nonlinear function. For instance, in the Koval model for a two-species mixture [16], we
µR(a) = µ(0)(1+(M1/4 − 1)a1)−4
where M = µ(0)/µ(1) is the mobility ratio.
We notice that the equations (2.1)-(2.3) are strongly coupled. Moreover, every concentration equation
(2.1) involves a different time scale. Therefore, it is difficult to build a numerical scheme that captures
all the physical phenomena. We have to transform these equations. We first assume that only serial or
parallel first-order reactions occur, so that ki j = yi with y1 = 0. Next, following [4], we assume that no
two isotopes have identical decay rates and we set
c1 = a1, ci = ai +
yl+1λl
λl −λi
a j for i = 2, ..,Nr. (2.4)
Study of a finite volume-finite element scheme for a nuclear transport model 3 of 26
Lastly, without losing any mathematical difficulty (see remark 2.2 below), we set Ri = 1 for i= 1, ..,Nr−
1 and φ = 1, Cp = 1. We also set sci = si +∑
∏i−1l= j
yl+1λl
λl−λi
s j and κ(c,θ ) = K/µ
i=1 ,θ
. The
contamination problem is now modelized by the following parabolic-elliptic system
∂tci + div(ci u)−Dc ∆ci = sci − sci −λi ci , (2.5)
∂tθ + div(θ u)−Dθ ∆θ =−sθ − s(θ −θ∗) , (2.6)
divu = s, u+κ(c,θ )∇p = f , (2.7)
with i = 1, ..,Nr − 1. These equations are completed with the boundary and initial conditions
∇ci ·ν = 0 , ∇θ ·ν = 0 , u ·n = 0 , (2.8)
ci|t=0 = ci0, θ |t=0 = θ0. (2.9)
The pressure p is normalized by
Ω pdx = 0. Equations (2.5) are all similar. Thus, for the sake of
simplicity, we will assume that there is only one. We set Nr = 2 and c = c1, c0 = c
0, sc = sc1 , λ = λ1.
The results of this paper readily extend to the general case.
We conclude with some notations and hypothesis. Let D be a bounded open set of Rk with k > 1.
We denote by C ∞0 (D) the set of functions that are continuous on D together with all their derivatives,
and have a compact support in D. For p ∈ {2,∞}, we use the Lebesgue spaces
Lp(D),‖.‖Lp(D)
Lp(D),‖.‖Lp(D)
with Lp = (Lp)2. We also use the Sobolev spaces W p,q(D) for p ∈ [1,∞[ and
q ∈ [1,∞[. In the case D = Ω we use the following conventions. We drop the domain dependancy.
We denote by |.| (resp. ‖.‖∞) the norms associated to L2 = L2(Ω) and L2 = L2(Ω) (resp. L∞ and
L∞). We set L20 = {v ∈ L2;
Ω v(x)dx = 0}. For p ∈ [1,∞[ we define (H p,‖.‖p) and (Hp,‖.‖p) with
H p = W p,2 and Hp = (H p)2. Now let (X , |.|) be a Banach space. The functions g : [0,T ] → X such
that t → ‖g(t)‖X is continuous (resp. bounded and square integrable) form the set C (0,T ;X) (resp.
L∞(0,T ;X) and L2(0,T ;X)). The associated norm for the space L∞(0,T ;X) (resp. L2(0,T ;X)) is
defined by ‖g‖L∞(0,T ;X) = supt∈[0,T ] ‖g(t)‖X (resp. ‖g‖L2(0,T ;X) =
0 ‖ f (t)‖2X dt
). Finally in all
computations we use C > 0 as a generic constant. It depends only on the data of the problem.
We assume the following regularities for the data in (2.5)–(2.7)
κ ∈W 1,∞((0,1)× (0,∞)) , s, sc, sθ ∈ L2 , f ∈ C (0,T ;L2). (2.10)
We also assume that κ > κin f with κin f > 0. For the initial data, we assume that c0 ∈ H1, θ0 ∈ H1, and
that we have a.e. in Ω
0 6 c0(x)6 1 , θ− 6 θ0(x)6 θ+ (2.11)
with θ− > 0. Finally we assume that we have a.e. in Ω
2s(x)+λ > sc(x)> 0 , 2(θ−−θ ∗)s(x)+ sθ (x)6 0 , 2(θ+−θ ∗)s(x)+ sθ (x)> 0. (2.12)
These conditions ensure a maximum principle (proposition 6.1 below).
REMARK 2.1 We have assumed that first-order reactions occur, so that the coefficients ki j in (2.1)
depend only on i. If ki j depends on i and j, one can still uncouple the equations by iterating the transfor-
mation (2.4), provided that k1 j = 0 for j = 2, ..,Nr −1. This assumption means that the first long-lasting
isotope disappears and is not created anymore. It is satisfied by many radionuclides.
4 of 26 C. CHOQUET AND S. ZIMMERMANN
REMARK 2.2 We have assumed that the retardation factors R j are identical. If it is not the case, the
difficulty and the approach remain the same. Indeed, let us consider the Fourier transform of (2.1). For
a Fourier mode â j(k, t) we obtain
(φ R j â j) =−φ (λ j R j + k2 Dc − ik u) â j + yi R j−1 φ â j−1 =−λ ′j R j φ â j + yi R j−1φ â j−1
with λ ′j = λ j +(k
2 Dc − ik u)/R j for j = 1, ..,Nr − 1. A transform analogous to (2.4) uncouples the
problem. By taking the partial differential equation counterpart, we obtain an equation similar to (2.5).
3. Discrete tools
3.1 Mesh and discrete spaces
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 π2 , so that xK ∈ K. The set of the edges of the triangle K ∈ Th is EK . The symbol
nK,σ denotes the unit normal vector 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 is |σ |. For each edge σ ∈ E inth let Kσ
and Lσ be the two triangles having σ in common; we set dσ = d(xKσ ,xLσ ). For all σ ∈ E exth only the
triangle Kσ located inside Ω is defined and we set dσ = d(xKσ ,xσ ). Then for all σ ∈ Eh we set τσ =
We assume the following on the mesh (see [13]): 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. (3.1)
We define on the mesh the following spaces. The usual space for finite volume schemes is
P0 = {q ∈ L2 ; ∀K ∈ Th, q|K is a constant}.
For any function qh ∈ P0 and any K ∈ Th we set qK = qh|K . We also consider
Pd1 = {q ∈ L2 ; ∀K ∈ Th, q|K is affine} ,
Pc1 = {qh ∈ Pd1 ; qh is continuous over Ω} ,
Pnc1 = {qh ∈ Pd1 ; ∀σ ∈ E inth , qh is continuous at the middle of σ}.
We have Pc1 ⊂ H1. On the other hand Pnc1 6⊂ H1, but Pnc1 ⊂ H1d with
H1d = {q ∈ L2 ; ∀K ∈ Th, q|K ∈ H1(K)}.
Thus we define ∇̃h : H1d → L2 by setting
∀qh ∈ H1d , ∀K ∈ Th, ∇̃hqh|K = ∇qh|K (3.2)
and the associated norm ‖.‖1,h is given by
∀qh ∈ H1d , ‖qh‖1,h = (|qh|2 + |∇̃hqh|2)1/2.
We then have the following Poincaré-like inequality for the space Pnc1 ∩L20 (see [1]).
Study of a finite volume-finite element scheme for a nuclear transport model 5 of 26
PROPOSITION 3.1 There exists C > 0 such that |qh|6C |∇̃hqh| for all qh ∈ Pnc1 ∩L20.
We also define discrete analogues of the norms H1 and H−1 for the space P0 by setting
‖qh‖h =
σ∈E inth
τσ (qLσ − qKσ )2
and ‖qh‖−1,h = sup
ψh∈P0
(qh,ψh)
‖ψh‖h
for any function qh ∈ P0. Note that for any ph ∈P0 and qh ∈ P0, (ph,qh)6 ‖qh‖−1,h‖qh‖h. The following
Poincaré-like inequality holds for the space P0 ∩L20 (see [13]).
PROPOSITION 3.2 There exists C > 0 such that |qh|6C‖qh‖h for all qh ∈ P0 ∩L20.
Finally we set P0 = (P0)2, Pd1 = (P
2 and use the Raviart-Thomas spaces [7]
RTd0 = {vh ∈ Pd1 ; ∀K ∈ Th , ∀σ ∈ EK , vh|K ·nK,σ is a constant} ,
RT0 = {vh ∈ RTd0 ; ∀σ ∈ E inth , vh|Kσ ·nKσ ,σ = vh|Lσ ·nKσ ,σ and vh ·n|∂Ω = 0}.
For all vh ∈ RT0, K ∈ Th and σ ∈ EK , we set (vh ·nK,σ )σ = vh|K ·nK,σ .
3.2 Projection operators
We associate with the spaces of section 3.1 some projection operators. First, we define ΠPc1 : H
d → Pc1
by setting
∀q ∈ H1d , ∀φh ∈ Pc1 ,
∇(ΠPc1 q),∇φh
= (∇q,∇φh). (3.3)
Next, we consider the space P0. Let Cd = {qh ∈ L2 ; qh is equal a.e. to a continuous function}. We
define ΠP0 : L
2 → P0 and Π̃P0 : Cd → P0 by setting
(ΠP0 p)K =
p(x)dx , (Π̃P0q)K = q(xK) , (3.4)
for all p ∈ L2, q ∈ Cd and K ∈ Th. We also set ΠP0 = (ΠP0)2. For the space Pnc1 , we define Π̃Pnc1 : L
Pnc1 and ΠPnc1 : H
1 → Pnc1 . For all p ∈ L2 and q ∈ H1, Π̃Pnc1 p and ΠPnc1 q satisfy
∀ψh ∈ Pnc1 , (Π̃Pnc1 p,ψh) = (p,ψh) ; ∀σ ∈ Eh ,
(ΠPnc1 q)dσ =
qdσ . (3.5)
For the space RT0, we define Π̃RT0 : L
2 → RT0 and ΠRT0 : H1 → RT0. For all v ∈ L2 and w ∈ H1,
Π̃RT0v and ΠRT0w satisfy
∀wh ∈ RT0 , (Π̃RT0v,wh) = (v,wh) ; ∀σ ∈ E
(w−ΠRT0w) ·nKσ ,σ dσ = 0. (3.6)
The operators ΠP0 , Π̃Pnc1 (resp. ΠP0 , Π̃RT0) are L
2 (resp. L2) projection operators. They are stable for
the L2 (resp. L2) norms. The operators Π̃P0 , ΠPnc1 and ΠRT0 are interpolation operators. The following
estimates are classical ([6] p.109 and [7]).
PROPOSITION 3.3 There exists C > 0 such that for all q ∈ H1 and v ∈ H1
|q−ΠP0q|6C h‖q‖1, |v−ΠRT0v|6C h‖v‖1.
6 of 26 C. CHOQUET AND S. ZIMMERMANN
For all p ∈ H1 and q ∈ H2 we have
|p−ΠPnc1 p|6C h‖p‖1 , |∇̃h(q−ΠPnc1 q)|6C h‖q‖2.
For all q ∈ Hd1 we have
|q−ΠPc1 q|6C h‖q‖1,h.
Finally, using the Sobolev embedding theorem, one checks that
|ΠP0q− Π̃P0q|6C h‖q‖W1,r (3.7)
for all q ∈W 1,r with r > 1 (see [22]).
3.3 Discrete operators
Equations (2.5)–(2.7) use the differential operators gradient, divergence and laplacian. We have to define
analogous operators in the discrete setting. The discrete gradient operator ∇h : Pnc1 →P0 is the restriction
to Pnc1 of the operator ∇̃h given by (3.2). The discrete divergence operator divh : P0 → Pnc1 is defined by
∀σ ∈ E inth , (divh vh)(xσ ) =
3 |σ |
|Kσ |+ |Lσ |
(vLσ − vKσ ) ·nK,σ ,
∀σ ∈ E exth , (divh vh)(xσ ) =−
3 |σ |
|Kσ |
vKσ ·nK,σ ,
for all vh ∈ P0. It is adjoint to ∇h (proposition 4.1 below). The discrete laplacian operator ∆h : P0 → P0
is the usual one for finite volume schemes (see [13]). For all qh ∈ P0 and K ∈ Th we have
|K| ∑
σ∈EK∩E inth
qLσ − qKσ
. (3.8)
Let us now consider the convection terms in (2.5) and (2.6). We define b̃ : H1 ×H1 → L2 by
b̃(v,q) = div(qv) (3.9)
for all q ∈ H1 and v ∈ H1. In order to define a discrete counterpart to b̃ we use the classical upwind
scheme (see [13]). The discrete operator b̃h : RT0 ×P0 → P0 is such that
b̃h(vh,qh)
|K| ∑
σ∈EK∩E inth
(vh ·nK,σ )+σ qK +(vh ·nK,σ )−σ qLσ
(3.10)
for all vh ∈ RT0, qh ∈ P0 and K ∈ Th. We have set a+ = max(a,0) and a− = min(a,0) for all a ∈ R.
Integrating by parts the convection terms also leads to consider b : L2 ×L2 ×L∞ →R defined by
b(v, p,q) =−
pv ·∇qdx (3.11)
for all v ∈ L2, p ∈ L2 and q ∈ L∞. The discrete counterpart is bh : RT0 ×P0×P0 → R with
bh(vh, ph,qh) = ∑
σ∈EK∩E inth
(vh ·nK,σ )+σ pK +(vh ·nK,σ )−σ pLσ
(3.12)
for all vh ∈ RT0, ph ∈ P0 and qh ∈ P0.
Study of a finite volume-finite element scheme for a nuclear transport model 7 of 26
4. Properties of the discrete operators
The properties of the discrete operators are analogous to the ones satisfied by their continuous counter-
part. The gradient and divergence operators are adjoint. For the operators ∇h and divh we state in [21]
the following.
PROPOSITION 4.1 For all vh ∈ P0 and qh ∈ Pnc1 we have (vh,∇hqh) =−(qh,divhvh).
Let us now consider the convection terms. Let q ∈ L∞ ∩H1, v ∈ L2 with divv ∈ L2 and divv(x) > 0
a.e. in Ω . We obtain b(v,q,q) =
2/2)divvdx > 0 by integration by parts. For bh we state in [21] a
similar result.
PROPOSITION 4.2 Let vh ∈ RT0 with divvh > 0. We have bh(vh,qh,qh)> 0 for all qh ∈ P0 .
The following stability properties are used to prove the error estimates in section 8.
PROPOSITION 4.3 There exists C > 0 such that for all ph ∈ P0, qh ∈ P0 and vh ∈ RT0 with divvh = 0
|bh(vh, ph,qh)|6C |vh|‖ph‖h‖qh‖h. (4.1)
There exists C > 0 such that for all ph ∈ P0, qh ∈ P0 ∩L20, vh ∈ RT0
|bh(vh, ph,qh)|6C (|vh|‖ph‖h + |divvh|‖ph‖∞ )‖qh‖h. (4.2)
PROOF. For all K ∈ Th and σ ∈ EK ∩E inth we write
(vh ·nK,σ )+σ pK +(vh ·nK,σ )−σ pLσ = (vh ·nK,σ )σ pK −|(vh ·nK,σ )σ |(pLσ − pK).
Thus (3.12) reads bh(vh, ph,qh) = S1 + S2 with
S1 =− ∑
σ∈EK∩E inth
|σ | |(vh ·nK,σ )σ |(pLσ − pK) , S2 = ∑
pK qK ∑
σ∈EK∩E inth
|σ |(vh ·nK,σ )σ .
Using the Cauchy-Schwarz inequality we write
|S1|=
∣∣∣ ∑
σ∈E inth
|σ | |vh(xσ ) ·nK,σ |(pLσ − pK)(qLσ − qK)
6 h‖vh‖∞
σ∈E inth
(pLσ − pKσ )2
)1/2( ∑
σ∈E inth
(qLσ − qKσ )2
Since vh ∈ RT0 ⊂ (Pd1 )2 we have h‖vh‖∞ 6C |vh| ([6] p. 112). Moreover (3.1) implies ∑σ∈E inth (pLσ −
pKσ )
6C ∑σ∈E inth τσ (pLσ − pKσ )
2 =C‖ph‖2h and ∑σ∈E inth (qLσ − qKσ )
6C‖qh‖2h. Thus
|S1|6C |vh|‖ph‖h ‖qh‖h. (4.3)
We now consider S2. We have (vh ·nK,σ )σ = 0 for all K ∈ Th and σ ∈ EK ∩E exth . Thus we write
σ∈EK∩E inth
|σ |(vh ·nK,σ )σ = ∑
|σ |(vh ·nK,σ )σ =
divvh dx.
8 of 26 C. CHOQUET AND S. ZIMMERMANN
It gives the following relation.
S2 = ∑
pK qK
divvh dx =
ph qh divvh dx.
Thus if divvh = 0 then S2 = 0 and estimate (4.3) gives (4.1). Let us prove (4.2). Since qh ∈ P0 ∩L20 we
can apply proposition 3.2. Using the Cauchy-Schwarz inequality we get
|S2|6 ‖ph‖∞ |qh| |divvh|6C‖ph‖∞ ‖qh‖h |divvh|.
This latter estimate together with (4.3) gives (4.2).
Lastly, we claim that b̃h is a consistent approximation of b̃ [21].
PROPOSITION 4.4 Let r > 0. There exists C > 0 such that for all functions q ∈ H2 and v ∈ H1+r with
v ·n|∂Ω = 0
‖ΠP0 b̃(v,q)− b̃h(ΠRT0v,Π̃P0q)‖−1,h 6C h‖q‖1‖v‖1+r.
Let us now consider the discrete laplacian. We have a coercivity and stability result.
PROPOSITION 4.5 For all ph ∈ P0 and qh ∈ P0, we have
−(∆h ph, ph) = ‖ph‖2h , |(∆h ph,qh)|6 ‖ph‖h‖qh‖h.
PROOF. Definition (3.8) implies
(∆h ph,qh) = ∑
σ∈EK∩E inth
τσ (pLσ − pKσ ) =− ∑
σ∈E inth
τσ (pLσ − pKσ )(qLσ − qKσ ). (4.4)
Setting qh = ph gives the first part of the result. Using the Cauchy-Schwarz inequality, we get the second
We also deduce from (4.4) the following property.
PROPOSITION 4.6 For all ph ∈ P0 and qh ∈ P0 we have (∆h ph,qh) = (ph,∆hqh).
Lastly, we state that ∆h is a consistent approximation of the laplacian. The proof follows the lines of the
one of proposition 1.14 in [22].
PROPOSITION 4.7 There exists C > 0 such that for all q ∈ H2 with ∇q ·n|∂Ω = 0 we have
‖ΠP0(∆q)−∆h(Π̃P0q)‖−1,h 6C h‖q‖2.
5. The finite element-finite volume scheme
We now introduce the scheme for (2.5)-(2.9). The interval [0,T ] is split with a constant time step
k = T/N. We set [0,T ] =
m=0[tm, tm+1] with tm = mk. The time derivatives are approximated using a
first order Euler scheme. The convection terms are discretized semi-implicitly in time and the other ones
in an implicit way. We set sh = ΠP0s, s
h = ΠP0sc, s
h = ΠP0sθ and f
h = ΠP0f(tm) for all m ∈ {0, . . . ,N}.
Since ΠP0 (resp. ΠP0) is stable for the L
2 (resp. L2) norm we have
|sh|6 |s| , |sch|6 |sc| , |sθh |6 |sθ | , |fmh |6 |f(tm)|6 ‖f‖L∞(0,T ;L2). (5.1)
Study of a finite volume-finite element scheme for a nuclear transport model 9 of 26
The initial values are c0h = ΠP0c0 and θ
h = ΠP0θ0. Then for all n ∈ {0, . . . ,N − 1}, the quantities
cn+1h ∈ P0, θ
h ∈ P0, p
h ∈ P
1 ∩L20, u
h ∈ RT0 are the solutions of the following problem.
cn+1h − c
−Dc ∆hcn+1h = s
h − (sh +λ )cn+1h − b̃h(u
h ), (5.2)
θ n+1h −θ
−Dθ ∆hθ n+1h =−s
h − sh (θ
h −θ∗)− b̃h(u
h ), (5.3)
divh(κn+1h ∇h p
h ) = divh f
h − Π̃Pnc1 sh , (5.4)
un+1h = Π̃RT0(f
h ∇h p
h ), (5.5)
with κn+1h = κ(c
h ) ∈ P0. This term is defined thanks to proposition 6.1 below. Note also that
the boundary conditions are implicitly included in the definition of the discrete operators (section 3.3).
The existence of a unique solution to (5.2) and (5.3) is classical (see [13]). Since κn+1h > κmin > 0 and
pn+1h ∈ L
0 equation (5.4) also has a unique solution (see [6]). We have a discrete equivalent for the
divergence condition (2.7).
PROPOSITION 5.1 For all m ∈ {1, . . . ,N} we have divumh = sh .
PROOF. Let m ∈ {1, . . . ,N} and n = m− 1. We compare the solution of (5.4)–(5.5) with the solution
of the following mixed hybrid problem. Let E0 =
µh : ∪σ∈Eh → R ; ∀σ ∈ Eh , µh|σ is constant
. Then
ũmh ∈ RT
0 , p
h ∈ P0 and λ
h ∈ E0 are the solution of (see [7])
∀vh ∈ RTd0 , (ũmh ,vh)+ ∑
κmK ∑
σ ′∈EK
|σ ′|λ mσ ′ (vh|K ·nK,σ ′)− ∑
|K|κmK pmK divvh|K = (fmh ,vh) , (5.6)
∀µh ∈ E0, ∑
µh (ũmh ·n)dσ = 0 , ∀K ∈ Th ,
div ũmh dx =
sdx , (5.7)
and p̃mh ∈ Pnc1 is defined by
h dσ = λ
σ for all σ ∈ Eh. Let σ ∈ Eh. We define φσ ∈ Pnc1 by setting
φσ (xσ ) = 1 and φσ (xσ ′) = 0 for all σ ′ ∈ Eh\{σ}. We set vh = ∇hφσ ∈ P0 ⊂ RTd0 in (5.6). We have
κmK ∑
σ ′∈EK
|σ ′|λ mσ ′ ∇hφσ |K ·nK,σ ′ = ∑
κmK ∇hφσ |K · ∑
σ ′∈EK
|σ ′|λ mσ ′ nK,σ ′
and according to the gradient formula
σ ′∈EK
|σ ′|λ mσ ′ nK,σ ′ = ∑
σ ′∈EK
p̃mh nK,σ ′ dσ
∇h p̃
h dx.
Thus we get from (5.6)
(ũmh ,∇hφσ )+ (κ
h ∇h p̃
h ,∇hφσ ) = (f
h ,∇hφσ ). (5.8)
The first term in (5.8) is treated as follows. Integrating by parts we get
(ũmh ,∇hφσ ) =−(φσ ,div ũmh )+ ∑
σ ′∈EK
φσ (ũmh |K ·nK,σ ′)dσ ′.
Since (5.7) implies that ũmh ∈ RT0, we have
σ ′∈EK
φσ (ũmh ·nK,σ ′)dσ ′ = ∑
σ∈E inth
|σ |φσ (xσ )(ũmh |Lσ ·nKσ ,σ − ũmh |Kσ ·nKσ ,σ ) = 0.
10 of 26 C. CHOQUET AND S. ZIMMERMANN
Thus (ũmh ,∇hφσ ) =−(φσ ,div ũmh ). Then, using (5.7), we get
(ũmh ,∇hφσ ) =−(φσ ,sh) =−(φσ ,Π̃Pnc1 sh).
Furthermore, according to proposition 4.1, we have (κmh ∇h p̃
h ,∇hφσ ) = −
φσ ,divh(κmh ∇h p̃
(fmh ,∇hφσ ) =−(φσ ,divhfmh ). Hence we deduce from (5.8) that
∀φσ ∈ Pnc1 ,
φσ ,divh(κmh ∇h p̃
h )− divhf
h + Π̃Pnc1 s
Since (φσ )σ∈Eh is a basis of P
1 , we get divh(κ
h ∇h p̃
h ) = divhf
h − Π̃Pnc1 s
h . Thus, by (5.4), there exists
a real C such that p̃mh = p
h +C. We now compare ũ
h with u
h . Since for all vh ∈ RT0 we have
κmK ∑
σ ′∈EK
|σ ′|λ mσ ′ (vh|K ·nK,σ ′) = ∑
σ∈E inth
|σ |φσ (xσ )(ũmh |Lσ ·nKσ ,σ − ũ
h |Kσ ·nKσ ,σ ) = 0,
it follows from (5.6) that (ũmh ,v) = (f
h −κmh ∇h p̃mh ,v) for any v ∈ RT0. It means that
ũmh = Π̃RT0(f
h −κmh ∇h p̃mh ) = Π̃RT0(f
h −κmh ∇h pmh ) = umh .
Thus umh = ũ
h satisfies (5.7) and divu
h = sh.
6. Stability analysis
We first check that a maximum principle holds. It ensures that the computed concentration and temper-
ature are physically relevant.
PROPOSITION 6.1 For any m ∈ {0, . . . ,N} we have 0 6 cmh 6 1 and θ− 6 θ
h 6 θ+.
PROOF. We prove the result by induction. Since c0h = ΠP0c0 and θ
h = ΠP0θ0 the result holds for m = 0
thanks to (2.11) and (3.4). Let us assume that it is true for m = n ∈ {0, . . . ,N−1}. Let K ∈Th. Equation
(5.2) implies
(1+ k sK + k λ )cn+1K = c
K + k s
K + k Dc ∑
σ∈EK∩E inth
τσ (cn+1Lσ − c
K )− k b̃h(u
We consider the last term of this relation. Since for any σ ∈ EK ∩E inth we have
(unh ·nK,σ )+σ cn+1K +(u
h ·nK,σ )−σ cn+1Lσ = (u
h ·nK,σ )σ cn+1K +(−u
h ·nK,σ )+σ (cn+1K − c
We deduce from (3.10)
−b̃h(unh,cn+1h )
−cn+1K ∑
σ∈EK∩E inth
|σ |(unh ·nK,σ )σ + ∑
σ∈EK∩E inth
(−unh ·nK,σ )+σ (cn+1Lσ − c
Since unh ∈ RT0, (unh ·nK,σ )σ = 0 for any σ ∈ EK ∩E exth . It implies that ∑σ∈EK∩E inth |σ |(u
h ·nK,σ )σ =
∑σ∈EK |σ |(u
h ·nK,σ )σ . Thus using the divergence formula and proposition 5.1 we obtain
|K| ∑
σ∈EK∩E inth
|σ |(unh ·nK,σ )σ =
divunh dx = sK .
Study of a finite volume-finite element scheme for a nuclear transport model 11 of 26
Therefore we get
(1+ 2k sK + k λ )cn+1K = c
K + k s
K + k Dc ∑
σ∈EK∩E inth
τσ (cn+1Lσ − c
|K| ∑
σ∈EK∩E inth
(−unh ·nK,σ )+σ (cn+1Lσ − c
K ). (6.1)
We consider Ki ∈ Th such that cn+1Ki = minK∈Th c
K . According to hypothesis (2.11) and definition
(3.4) we have 2sKi +λ > 0 and s
> 0. Thus, using the induction hypothesis, we deduce from (6.1)
cn+1K = c
cnKi + ks
1+ 2k sKi + k λ
k scKi
1+ 2k sKi + k λ
We now consider Ks ∈ Th such that cn+1Ks = maxK∈Th c
K . Using again hypothesis (2.11) we have
2sKs +λ > s
> 0. Thus, using the induction hypothesis, we deduce from (6.1)
cn+1K = c
cnKs + k s
1+ 2k sKs + k λ
1+ k scKs
1+ 2k sKs + k λ
A similar work for equation (5.3) proves that θ− 6 minK∈Th θ
K and maxK∈Th θ
K 6 θ+. Thus the
induction hypothesis still holds for m = n+ 1.
We now state the stability of the scheme (5.2)-(5.5).
PROPOSITION 6.2 For any 1 6 m 6 N we have
‖cnh‖2h + k
‖θ nh ‖2h 6C , (6.2)
|umh |+ |∇h pmh |6C. (6.3)
PROOF. Let 0 6 n 6 N − 1. Multiplying (5.2) by 2k cn+1h we get
(cn+1h − c
h )− 2k Dc (∆hc
h )+ k
(sh +λ )cn+1h ,c
+ k bh(u
h ) = k (s
We have (cn+1h − c
h ) = |c
2 −|cnh|
2 + |cn+1h − c
2. Thanks to propositions 4.2 and 4.5
−2k (∆hcn+1h ,c
h ) = 2k‖c
h , bh(u
h )> 0.
Using the Cauchy-Schwarz and Young inequalities we write
k (sch,c
h )6 k |s
h| |cn+1h |6C k |c
h |6 k
|cn+1h |
2 +C k.
Finally thanks to (2.11) and (3.4) we have sh > 0. Thus we obtain
|cn+1h |
2 −|cnh|
2 + 2k Dc‖cn+1h ‖
h + k
|cn+1h |
6C k.
12 of 26 C. CHOQUET AND S. ZIMMERMANN
Let m ∈ {1, . . . ,N}. Summing up the latter relation from n = 0 to m− 1 we get
|cmh |2 + 2k Dc
‖cnh‖2h 6 |c0h|
k 6C ,
thanks to proposition 6.1. With a similar work on equation (5.3), we get (6.2). We now prove (6.3). Let
n = m− 1 ∈ {0, . . . ,N − 1}. Multiplying equation (5.4) by −pn+1h and using proposition 4.1, we get
(κn+1h ∇h p
h ,∇h p
h ) = (f
h ,∇h p
h )+ (Π̃Pnc1 sh, p
h ). (6.4)
The left-hand side term satisfies (κn+1h ∇h p
h ,∇h p
h ) > κin f |∇h p
2. We now consider the right-
hand side. Using (5.1), the Cauchy-Schwarz and Young inequalities we write
|(fn+1h ,∇h p
h )|6 |f
h | |∇h p
κin f
|∇h pn+1h |
2 +C‖f‖2L∞(0,T ;L2).
Also, the stability of Π̃Pnc1 for the L
2-norm, proposition 3.1 and the Young inequality lead to
|(Π̃Pnc1 sh, p
h )|6 |sh| |p
h |6C |p
h |6C |∇h p
κin f
|∇h pn+1h |
2 +C.
Thus we deduce from (6.4) that |∇h pn+1h |
2 = |∇h pmh |
6C. Then (5.5) imply
|umh |= |un+1h |6 |f
h |+ |κ
h ∇h p
h |6 ‖f‖L∞(0,T ;L2)+ ‖κ‖W1,∞((0,1)×(0,∞)) |∇h p
h |6C.
Estimate (6.3) is proven.
7. Convergence analysis
Let ε = max(h,k). In this section we study the behavior of the scheme (5.2)-(5.5) as ε → 0. We first
define the applications cε : R→ P0, c̃ε : R→ P0, θε : R → P0, pε : R→ Pnc1 , sε : R→ P0, scε : R→ P0
and uε : R→ RT0, fε : R→ P0 by setting for all n ∈ {0, . . . ,N − 1} and t ∈ [tn, tn+1]
cε(t) = c
h , c̃ε(t) = c
(t − tn)(cn+1h − c
h), θε(t) = θ
pε(t) = p
h , sε(t) = sh, s
ε (t) = s
h, uε(t) = u
h, fε(t) = f
and for all t 6∈ (0,T )
cε(t) = c̃ε(t) = θε(t) = pε(t) = sε (t) = scε (t) = 0, uε(t) = fε(t) = 0.
We recall that the Fourier transform f̂ of a function f ∈ L1(R) is defined for any τ ∈ R by
f̂ (τ) =
e−2iπτt f (t)dt. (7.1)
We begin with the following estimate.
PROPOSITION 7.1 Let 0 < γ < 14 . There exists C > 0 such that for all ε > 0
|τ|2γ (|ĉε (τ)|2 + |θ̂ε(τ)|2)dτ 6C.
Study of a finite volume-finite element scheme for a nuclear transport model 13 of 26
PROOF. Since equations (5.2) and (5.3) are similar we only prove the estimate on ĉε . We first define
gε : R→ P0 ∩L20 as the solution of
∆hgε = Dc ∆hcε + scε − (sε +λ )cε − b̃h(uε ,cε ).
Multiplying this equation by −gε we obtain
−(∆hgε ,gε) =−Dc
∆hcε ,gε
scε − (sε +λ )cε ,gε
+ bh(uε ,cε ,gε). (7.2)
Proposition 4.5 allows us to write
−(∆hgε ,gε) = ‖gε‖2h, −(∆hcε ,gε)6 ‖cε‖h ‖gε‖h.
Thanks to the Cauchy-Schwarz inequality, (5.1) and proposition 3.2 we have
∣∣(scε − (sε +λ )cε ,gε
)∣∣6C (|sc|+ |s|+λ ) |gε |6C‖gε‖h.
According to proposition 4.3, then proposition 6.1 and (6.3), we have
uε ,cε ,gε
)∣∣6C‖cε‖∞‖gε‖h |divuε |+C‖cε‖h ‖gε‖h |uε |
6C‖gε‖h |divuε |+C‖cε‖h ‖gε‖h.
Let us plug these estimates into (7.2) and integrate from 0 to T . We get
‖gε‖h dt 6C
|divuε |dt +C
‖cε‖h dt 6C ,
because of proposition 5.1 and (6.2). Definition (7.1) then leads to
∀τ ∈ R , ‖ĝε(τ)‖h 6C. (7.3)
We now use this estimate to prove (7.1). Equation (5.2) reads
c̃ε = ∆hgε +(c0hδ0 − cNh δT )
where δ0 and δT are Dirac distributions respectively localized in 0 and T . Let τ ∈ R. Applying the
Fourier transform to the latter equation we obtain
−2iπτ ̂̃cε (τ) = ∆hĝε(τ)+ (c0h − cNh e−2iπτT ).
Let us take the scalar product of this relation with isign(τ)̂̃cε(τ). Applying propositions 3.2 and 4.5
leads to
2π |τ| |̂c̃ε(τ)|2 6C
‖ĝε(τ)‖h + |c0h|+ |c
‖̂̃cε(τ)‖h.
We assume that τ 6= 0 and multiply this estimate by |τ|2γ−1. Using proposition 6.1 and (7.3) we get
|τ|2γ |̂c̃ε(τ)|2 6C |τ|2γ−1 ‖̂̃cε(τ)‖h.
Using the Young inequality and integrating over {τ ∈ R ; |τ|> 1}, we obtain
|τ|>1
|τ|2γ |̂c̃ε(τ)|2 dτ 6
|τ|>1
|τ|4γ−2 dτ +C
|τ|>1
‖̂̃cε(τ)‖2h dτ.
14 of 26 C. CHOQUET AND S. ZIMMERMANN
For |τ|6 1, we have |τ|2γ |̂c̃ε (τ)|2 6 |̂c̃ε(τ)|2 6C‖̂̃cε(τ)‖2h according to proposition 3.2. Thus
|τ|61
|τ|2γ |̂c̃ε(τ)|2 dτ 6C
|τ|61
‖̂̃cε(τ)‖2h dτ.
By combining the bounds for |τ|> 1 and |τ|6 1 we get
|τ|2γ |̂c̃ε(τ)|2 dτ 6
|τ|>1
|τ|4γ−2 dτ +C
‖̂̃cε(τ)‖2h dτ.
Since 4γ − 2 <−1, we have
|τ|>1 |τ|4γ−2 dτ 6C. Thanks to the Parseval theorem and (6.2)
‖̂̃cε(τ)‖2h dτ 6
‖c̃ε‖2h dt 6C
k‖c0h‖
h + k
‖cnh‖2h
because ‖c0h‖h = ‖ΠP0c0‖h 6C‖c0‖1 (see [13] p. 776). Hence the result.
We can now prove the following convergence result.
PROPOSITION 7.2 There exists a subsequence of (cε ,θε , pε ,uε)ε>0, not relabeled for convenience, such
that the following convergences hold for ε → 0
cε → c in L2(0,T ;L2), θε → θ in L2(0,T ;L2), (7.4)
pε ⇀ p weakly in L
2(0,T ;H1), uε ⇀ u weakly in L
2(0,T ;L2). (7.5)
The limits (c,θ , p,u) satisfy the following properties. We have c ∈ L2(0,T ;H1), θ ∈ L2(0,T ;H1),
p ∈ L∞(0,T ;H1) and u ∈ L∞(0,T ;L2). We also have 0 6 c(x, t) 6 1 and θ− 6 θ (x, t) 6 θ+ a.e. in
Ω × [0,T ]. For all φ ∈ C ∞0 (Ω × (−1,T )), c and θ satisfy
(c,∂tφ)+Dc (∇c,∇φ)− c(u ·∇φ)−
sc − (s+λ )c
c0,φ(·,0)
, (7.6)
(θ ,∂tφ)+Dθ (∇θ ,∇φ)−θ (u ·∇φ)+
sθ + s(θ −θ∗)
θ0,φ(·,0)
. (7.7)
Lastly we have
u = f−κ(c,θ )∇p in L∞(0,T ;L2) , divu = s in L2. (7.8)
PROOF. In what follows, the convergence results hold for extracted subsequences. They are not rela-
beled for convenience. We begin by proving (7.4). According to proposition 6.1, the sequence (cε)ε>0
is uniformly bounded in L∞(0,T ;L2). Thus there exists c ∈ L∞(0,T ;L2) such that
cε ⇀ c weakly in L
2(0,T ;L2).
Using the Fourier transform, we prove that this convergence is strong. Let dε = cε − c and M > 0. We
use the following splitting
|d̂ε(τ)|2 dτ =
|τ|>M
|d̂ε(τ)|2 dτ +
|τ|6M
|d̂ε(τ)|2 dτ = IMε + JMε . (7.9)
Study of a finite volume-finite element scheme for a nuclear transport model 15 of 26
Since |d̂ε(τ)|2 6 2|ĉε(τ)|2 + 2|ĉ(τ)|2 we have
IMε 6 2
|τ|>M
|ĉε(τ)|2 dτ + 2
|τ|>M
|ĉ(τ)|2 dτ.
Using proposition 7.1 we write
|τ|>M
|ĉε(τ)|2 dτ 6
|τ|>M
|τ|2γ |ĉε(τ)|2 dτ 6
Hence
IMε 6
|τ|>M
|ĉ(τ)|2 dτ.
This implies that for all ε > 0, IMε → 0 when M → ∞. We now consider JMε . Let τ ∈ [−M,M]. Since
cε(t) ∈ P0 for all t ∈ R, and cε ⇀ c weakly in L2(0,T ;L2), we deduce from (7.1) that ĉε(τ) ∈ P0 and
ĉε(τ)⇀ ĉ(τ) weakly in L2. Extanding ĉε(τ) by 0 outside Ω , one checks ([13] p.811) that
∀η ∈ R2 , |ĉε(τ)(·+η)− ĉε(τ)|6C‖ĉε(τ)‖h |η |(|η |+ h). (7.9)
Then, using estimate (6.2), we deduce from [13] (p.834) that ĉε(τ)→ ĉ(τ) strongly in L2. Thus d̂ε(τ) =
ĉε(τ)− ĉ(τ) → 0 in L2, so that JεM → 0 when ε → 0. Now, let us report the limits for IεM and JεM into
(7.9). Using the Parseval identity we get
|d̂ε(τ)|2 dτ =
|dε |2 dt =
|cε − c|2 dt → 0.
Thus we have proven that cε → c in L2(0,T ;L2). A similar work proves that θε → θ in L2(0,T ;L2)
with θ ∈ L∞(0,T ;L2). Hence (7.4) is proven. Moreover using proposition 6.1 we obtain 0 6 c(x, t)6 1
and θ− 6 θ (x, t) 6 θ+ a.e. in Ω × [0,T ]. Lastly, using (6.2) and (7.9), we get as in [13] (p.811) that
c ∈ L2(0,T ;H1) and θ ∈ L2(0,T ;H1).
Let us now consider the sequences (pε)ε>0 and (uε )ε>0. According to (3.3) and (6.3) the sequence
(ΠPc1 pε)ε>0 is bounded in L
∞(0,T ;H1). It implies that there exists p ∈ L∞(0,T ;H1) such that ΠPc1 pε ⇀
p weakly in L2(0,T ;H1). Using proposition 3.3 we get pε ⇀ p weakly in L
2(0,T ;H1). Moreover,
according to (6.3), the sequence (uε)ε>0 is bounded in L∞(0,T ;L2) . Thus we have uε ⇀ u weakly
in L2(0,T ;L2) with u ∈ L∞(0,T ;L2). We check the properties of u. Using a Taylor expansion, the
Cauchy-Schwarz inequality, and a density argument, we have
‖κ(c,θ )−κ(cε,θε )‖L2(0,T ;L2) 6 ‖κ‖W1,∞((0,1)×(0,∞)) (‖c− cε‖L2(0,T ;L2)+ ‖θ −θε‖L2(0,T ;L2)).
Thus, using the strong convergence of the sequences (cε)ε>0 and (θε )ε>0, we have κ(cε ,θε )→ κ(c,θ )
in L2(0,T ;L2). Since ∇h pε ⇀ ∇p weakly in L2(0,T ;L2), we deduce from this
κ(cε ,θε )∇h pε ⇀ κ(c,θ )∇p weakly in L2(0,T ;L2). (7.10)
Now let v ∈ L2(0,T ;(C ∞0 )2). According to (5.5) we have
(uε ,ΠRT0v) = (fε −κ(cε ,θε)∇h pε ,v).
Using proposition 3.3 one checks easily that (fε ,ΠRT0v)→ (f,v) and (uε ,ΠRT0v)→ (u,v) in L1(0,T ).
Using moreover convergence (7.10) and a density argument, we deduce from this that u= f−κ(c,θ )∇p.
And since divuε = sε by proposition 5.1, we also have divu = s.
16 of 26 C. CHOQUET AND S. ZIMMERMANN
We finally prove that c satisfies (7.6). For all t ∈ (0,T ) equation (5.2) reads
c̃ε −Dc ∆hcε + b̃h(uε ,cε) = scε − (sε +λ )cε .
Let ψ ∈ C ∞0 (Ω × (−1,T )) and ψh = Π̃P0ψ . Multiplying the latter equation by ψh and integrating over
[0,T ], we obtain
c̃ε ,ψh
dt −Dc
(∆hcε ,ψh)dt +
bh(uε ,cε ,ψh)dt =
(scε − (sε +λ )cε ,ψh) dt. (7.11)
We now pass to the limit ε → 0 in this equation. We begin with the term
0 bh(uε ,cε ,ψh)dt. We use
the splitting b(u,c,ψ)− bh(uε ,cε ,ψh) = Aε1 +Aε2 +Aε3 with
Aε1 = b(u,c,ψ)− b(uε ,c,ψ), Aε2 = b(uε ,c,ψ)−
div(cuε)ψh dx ,
Aε3 =
div(cuε)ψh dx− bh(uε ,cε ,ψh).
According to definition (3.11)
Aε1 = b(u,c,ψ)− b(uε ,c,ψ) =−
c(u−uε) ·∇ψ dx.
We know that c∇ψ ∈ L2(0,T ;L2). Since uε ⇀ u in L2(0,T ;L2) we get
1 dt → 0. We now consider
Aε2. We have
Aε2 =
(ψ −ψh)div(cuε)dx =
(ψ −ψh)(uε ·∇c+ cdivuε)dx.
Using the Cauchy-Schwarz inequality we get
|Aε2|dt 6 ‖ψ −ψh‖L∞(Ω×(0,T)) (‖uε‖L2(0,T ;L2)+ ‖divuε‖L2(0,T ;L2))‖c‖L2(0,T ;H1).
Using a Taylor expansion, one checks that ‖ψ −ψh‖L∞(Ω×(0,T)) 6 h‖∇ψ‖L∞(Ω×(0,T )). Thus
2 dt →
0. Finally we estimate Aε3. For all triangles K ∈ Th and L ∈ Th sharing an edge σ , we set cK,L = cK if
uε ·nK,σ > 0 and cK,L = cL otherwise. Using the divergence formula, we deduce from definition (3.12)
Aε3 = ∑
σ∈EK∩E inth
(c− cK,Lσ )(uε ·nK,σ )dσ
σ∈E inth
(ψKσ −ψLσ )
(c− cKσ ,Lσ )(uε ·nKσ ,σ )dσ .
Using definition (3.5) this reads
Aε3 = ∑
σ∈E inth
(ψKσ −ψLσ )
(ΠPnc1 c− cKσ ,Lσ )(uε ·nKσ ,σ )dσ
σ∈E inth
(ψKσ −ψLσ ) |σ |
(ΠPnc1 c)(xσ )− cKσ ,Lσ
(uε ·nKσ ,σ )σ .
Study of a finite volume-finite element scheme for a nuclear transport model 17 of 26
Using a Taylor expansion, one checks that |ψKσ −ψLσ |6 h‖∇ψ‖L∞(Ω×(0,T )). Moreover |σ |6 h. Thus,
using the Cauchy-Schwarz inequality, we have
|Aε3| 6 C h2 ∑
σ∈E inth
|uε (xσ )|
∣∣∣(ΠPnc1 c)(xσ )− cKσ ,Lσ
6 C h2
σ∈E inth
|uε (xσ )|2
)1/2( ∑
σ∈E inth
|(ΠPnc1 c)(xσ )− cKσ ,Lσ |
Using the assumption on the mesh, one checks that |K| > C h2 for all K ∈ Th. Thus, thanks to a
quadrature formula, we have
|Aε3| 6 C
σ∈EK∩E inth
|uε (xσ )|2
)1/2(
σ∈EK∩E inth
|(ΠPnc1 c)(xσ )− cKσ |
6 C |uε | |ΠPnc1 c− cε |.
We write ΠPnc1 c− cε = (ΠPnc1 c− c)+ (c− cε) and we use proposition 3.3 . We obtain with (6.3)
|Aε3|dt 6C‖uε‖L∞(0,T ;L2) (h‖c‖L2(0,T ;H1)+ ‖c− cε‖L2(0,T ;L2)).
Since cε → c in L2(0,T ;L2) when ε = max(h,k) → 0, we conclude that
3 dt → 0. Gathering the
limits for Aε1, A
3, we obtain
bh(uε ,cε ,Π̃P0ψ)dt →
b(u,c,ψ)dt.
We now consider the other terms in (7.11). Proposition 4.6 leads to
(∆hcε ,Π̃P0ψ) =
cε ,∆h(Π̃P0ψ)
cε ,∆h(Π̃P0 ψ)−∆ψ
+(cε ,∆ψ). (7.12)
According to proposition 4.7
cε ,∆h(Π̃P0ψ)−∆ψ
)∣∣∣ 6 ‖cε‖h ‖∆h(Π̃P0ψ)−∆ψ‖−1,h 6C h‖cε‖h‖ψ‖2.
We then apply the Cauchy-Schwarz inequality and use (6.2). We obtain
cε ,∆h(Π̃P0 ψ)−∆ψ
)∣∣∣ dt 6C h
‖cε‖2h dt
‖cnh‖
6C h.
Moreover, since cε → c in L2(0,T ;L2), we have
0 (cε ,∆ψ)dt →
0 (c,∆ψ)dt. Thus we deduce from
(7.12) ∫ T
(∆hcε ,Π̃P0ψ)dt →
(c,∆ψ)dt.
We are left with two terms. First, using Taylor expansions, one checks that
ψh → ψ , ∂tψh → ∂tψ in L2(Ω × (−1,T)) , ψh(·,0)→ ψ(·,0) in L2. (7.13)
18 of 26 C. CHOQUET AND S. ZIMMERMANN
We know that cε → c in L2(0,T ;L2). Thus
(sc − (s+λ )cε ,ψh) dt →
(sc − (s+λ )c,ψ) dt.
Finally, integrating by parts the first term of (7.11), we get
c̃ε ,ψh
dt = (c̃ε ,ψh)t=T − (c̃ε ,ψh)t=0 −
(c̃ε ,∂tψh)dt.
Since ψ ∈ C ∞0 (Ω × (−1, ,T )) we have (c̃ε ,ψh)t=T = 0. Using proposition 3.3 one checks that c0h =
ΠP0c0 → c0 in L2; using moreover (7.13) we get
(c̃ε ,ψh)t=0 =
c0h,ψh(·,0)
ΠP0c0,ψh(·,0)
c0,ψ(·,0)
For the last term, one easily checks that ‖c̃ε − cε‖L2(0,T ;L2) → 0. Thus, since cε → c in L2(0,T ;L2),
we also have c̃ε → c in L2(0,T ;L2). Using moreover (7.13) we get
0 (c̃ε ,∂t ψh)dt →
0 (c,∂tψ)dt.
Therefore ∫ T
c̃ε ,ψh
dt →−
c0,ψ(·,0)
(c,∂tψ)dt.
By gathering the limits we have obtained in (7.11) we get (7.6). The relation (7.7) for θ is proven in a
similar way.
8. Error estimates
We have proven in section 7 that the problem (2.5)–(2.9) has a weak solution (c,θ , p,u). From now on,
we assume the following regularity for this solution:
c,θ ∈ C (0,T ;H2) , ct ,θt ∈ L2(0,T ;H1+r)∩C (0,T ;L2) ,
ctt ,θtt ∈ L2(0,T ;L2) , p ∈ C (0,T ;H2) , u ∈ C (0,T ;H1+s) ,
with r > 0 and s > 0. We also assume that f ∈ C (0,T ;H1).
8.1 Definitions
For all m ∈ {0, . . . ,N}, we define the following errors
emh,c = c(tm)− cmh , emh,θ = c(tm)−θ mh ,
emh,p = p(tm)− pmh , emh,u = u(tm)−umh .
We have the following splittings
emh,c = ε
h,c +η
h,c, e
h,θ = ε
h,θ +η
h,θ ,
emh,p = ε
h,p +η
h,p, e
h,u = ε
h,u+η
with the discrete errors
εmh,c = Π̃P0c(tm)− c
h , ε
h,θ = Π̃P0θ (tm)−θ
εmh,p = ΠPnc1 p(tm)− p
h , ε
h,u= ΠRT0u(tm)−u
Study of a finite volume-finite element scheme for a nuclear transport model 19 of 26
and the interpolation errors
ηmh,c = c(tm)− Π̃P0c(tm), η
h,θ = θ (tm)− Π̃P0θ (tm),
ηmh,p = p(tm)−ΠPnc1 p(tm), η
h,u= u(tm)−ΠRT0u(tm).
The interpolation errors are estimated as follows. We write |ηmh,c| 6 |c(tm)−ΠP0c(tm)|+ |ΠP0c(tm)−
Π̃P0c(tm)| and the same for ηmh,θ . Using proposition 3.3 and (3.7) we obtain
|ηmh,c|6C h‖c(tm)‖1 6C h‖c‖L∞(0,T ;H1) , |ηmh,θ |6C h‖θ‖L∞(0,T ;H1). (8.1)
According to proposition 3.3 we also have
|ηmh,p|+ |∇̃hηmh,p|6C h‖p(tm)‖2 6C h‖p‖L∞(0,T ;H2), (8.2)
|ηmh,u|6C h‖u(tm)‖1 6C h‖u‖L∞(0,T ;H1). (8.3)
We now have to estimate the discrete errors.
PROPOSITION 8.1 For all n ∈ {0, . . . ,N − 1} and ψh ∈ Pnc1 we have
εn+1h,c − ε
−Dc ∆hεn+1h,c + b̃h
εnh,u,Π̃P0c(tn+1)
+ b̃h(u
h,c )+ (s
h +λ )ε
h,c =C
h,1 +C
h,2 , (8.4)
εn+1h,θ − ε
−Dθ ∆hεn+1h,θ + b̃h
εnh,u,Π̃P0θ (tn+1)
+ b̃h(u
h,θ )+ sh ε
h,θ =Θ
h,1 +Θ
h,2 , (8.5)
κ(cn+1h ,θ
h )∇hε
h,p ,∇hψh
(κn+1h,1 ε
h,c +κ
h,2 ε
h,θ )∇p(tn+1),∇hψh
Pn+1h ,∇hψh
, (8.6)
εn+1h,u =−Π̃RT0
(κn+1h,1 ε
h,c +κ
h,2 ε
h,θ )∇p(tn+1)+κ(c
h )∇hε
−Un+1h . (8.7)
For all m ∈ {0, . . . ,N}, the consistency errors Cmh,1, C
h,2, Θ
h,1, Θ
h,2, P
h and U
h are defined in (8.9),
(8.10), (8.14), (8.15) and the terms κmh,1 and κ
h,2 are given by (8.13) below.
PROOF. Let n ∈ {0, . . . ,N − 1}. Equation (2.5) for t = tn+1 reads
∂tc(tn+1)−Dc ∆c(tn+1)+ b̃
u(tn+1),c(tn+1)
= sc − (s+λ )c(tn+1).
We introduce the time discretization by setting
Rn+1 =
(c(tn+1)− c(tn)
− ct(tn+1)
u(tn)−u(tn+1),c(tn+1)
We get
c(tn+1)− c(tn)
−Dc ∆c(tn+1)+ b̃
u(tn),c(tn+1)
= sc − (s+λ )c(tn+1)+Rn+1.
We apply ΠP0 to this equation. By subtracting the result from (5.2) we get
(c(tn+1)− c(tn)
cn+1h − c
ΠP0∆c(tn+1)−∆hc
+ΠP0 b̃
u(tn),c(tn+1)
− b̃h(unh,cn+1h )+ΠP0
(s+λ )c(tn+1)
− (sh +λ )cn+1h = ΠP0R
n+1. (8.8)
20 of 26 C. CHOQUET AND S. ZIMMERMANN
We now introduce the discrete errors as follows. Since c(tn+1)− c(tn) =
∫ tn+1
tn ct(s)ds one checks that
(c(tn+1)− c(tn)
cn+1h − c
∫ tn+1
ΠP0ct(s)− Π̃P0ct(s)
(εn+1h,c − ε
h,c).
We also have
ΠP0∆c(tn+1)−∆hc
h = ΠP0∆c(tn+1)−∆h
Π̃P0c(tn+1)
+∆hεn+1h,c .
Using the linearity of b̃h, one easily checks that
ΠP0 b̃
u(tn),c(tn+1)
− b̃h(unh,cn+1h ) = b̃h(u
h,c )+ b̃h
εnh,u,Π̃P0c(tn+1)
+ΠP0 b̃
u(tn),c(tn+1)
− b̃h
ΠRT0u(tn),Π̃P0c(tn+1)
Lastly
(s+λ )c(tn+1)
− (sh +λ )cn+1h = ΠP0
(s+λ )ηn+1h,c
+(sh +λ )εn+1h,c .
Using these relations in (8.8) we get (8.4). For any m ∈ {1, . . . ,N}, the consistency errors Cmh,1 ∈ P0 and
Cmh,2 ∈ P0 are given by
Cmh,1 = ΠP0
(c(tm)− c(tm−1)
− ct(tm)+ b̃
u(tm−1)−u(tm),c(tm)
+ ΠP0
(s+λ )ηmh,c
Π̃P0ct(s)−ΠP0ct(s)
ds, (8.9)
Cmh,2 = Dc
ΠP0∆c(tm)−∆h
Π̃P0c(tm)
ΠP0 b̃
u(tm−1),c(tm)
− b̃h(ΠRT0u(tm−1),Π̃P0c(tm)
A similar proof leads to (8.5) where the consistence errors Θ mh,1 ∈ P0 and Θ mh,2 ∈ P0 are defined for any
m ∈ {1, . . . ,N} by
Θ mh,1 = ΠP0
(θ (tm)−θ (tm−1)
−θt(tm)+ b̃
u(tm−1)−u(tm),θ (tm)
+ ΠP0(sη
h,θ )−
Π̃P0θt(s)−ΠP0θt(s)
ds, (8.10)
Θ mh,2 = Dθ
ΠP0∆θ (tm)−∆h
Π̃P0θ (tm)
ΠP0 b̃
u(tm−1),θ (tm)
− b̃h
ΠRT0u(tm−1),Π̃P0θ (tm)
We now consider the problem associated with the pressure. Let n ∈ {0, . . . ,N − 1} and ψh ∈ Pnc1 . Mul-
tiplying equation (2.7) written for t = tn+1 by ψh and integrating by parts, we get
κ(c(tn+1),θ (tn+1))∇p(tn+1),∇hψh
= (f(tn+1),∇hψh)+ (s,ψh). (8.11)
On the other hand, using (5.4) and proposition 4.1, we have
κ(cn+1h ,θ
h )∇h p
h ,∇hψh
= (fn+1h ,∇hψh)+ (Π̃Pnc1 sh,ψh).
Since ∇hψh ∈ P0, one checks that (fn+1h ,∇hψh) = (ΠP0f(tn+1),∇hψh) = (f(tn+1),∇hψh). According to
(3.5) we also have (Π̃Pnc1 sh,ψh) = (sh,ψh). Thus
κ(cn+1h ,θ
h ) ∇h p
h ,∇hψh
= (f(tn+1),∇hψh)− (sh,ψh).
Study of a finite volume-finite element scheme for a nuclear transport model 21 of 26
Substracting (8.11) from the latter relation, we obtain
κ(c(tn+1),θ (tn+1))∇p(tn+1)−κ(cn+1h ,θ
h )∇h p
h ,∇hψh
=−(s− sh,ψh). (8.12)
We split the left-hand side as
κ(cn+1h ,θ
h )(∇p(tn+1)−∇h p
κ(c(tn+1),θ (tn+1))−κ(cn+1h ,θ
∇p(tn+1).
Using a Taylor expansion, one can check that
κ(c(tn+1),θ (tn+1))−κ(cn+1h ,θ
h ) = (ε
h,c +η
h,c )κ
h,1 +(ε
h,θ +η
h,θ )κ
h,2 .
We have set for any m ∈ {0, . . . ,N} and s ∈ [0,1]
cmh (s) = c
h +(c(tm)− c
h )s, θ
h (s) = θ
h +(θ (tm)−θ
h )s,
κmh,1 =
κx(cmh (s),θ
h (s))ds, κ
h,2 =
κy(cmh (s),θ
h (s))ds.
(8.13)
We also have
∇p(tn+1)−∇hpn+1h = ∇hε
h,p + ∇̃hη
h,p .
Plugging these relations into (8.12) we get (8.6). For all m ∈ {0, . . . ,N} we have
Pmh = (κ
h,c +κ
h,θ )∇p(tm)+κ(c
h )∇hη
h,p. (8.14)
We end with the equation associated with u. Let n ∈ {0, . . . ,N − 1}. Applying the operator Π̃RT0 to
(2.7) for t = tn+1 we obtain
Π̃RT0u(tn+1) = Π̃RT0f(tn+1)− Π̃RT0
κ(c(tn+1),θ (tn+1))∇p(tn+1)
Let us substract this equation from (5.5). Since fn+1h = ΠP0f(tn+1) we get
Π̃RT0u(tn+1)−u
h = Π̃RT0
f(tn+1)−ΠP0f(tn+1)
− Π̃RT0
c(tn+1),θ (tn+1)
∇p(tn+1)−κn+1h ∇h p
One easily checks that
Π̃RT0u(tn+1)−u
h = Π̃RT0(u(tn+1)−ΠRT0u(tn+1))+ ε
h,u .
Thus we get (8.7). For all m ∈ {0, . . . ,N}, we have
Umh = Π̃RT0
(f(tm)−ΠP0f(tm))−ηmh,u−P
. (8.15)
This ends the proof of proposition 8.1.
22 of 26 C. CHOQUET AND S. ZIMMERMANN
8.2 Error estimates
We first estimate the consistency errors.
PROPOSITION 8.2 For all m ∈ {1, . . . ,N} the consistency errors satisfy
|Cnh,1|
2 + k
|Θ nh,1|
6C (h2 + k2), (8.16)
‖Cnh,2‖2−1,h + k
‖Θ nh,2‖2−1,h 6C h2 , (8.17)
|Pmh |+ |Umh |6C h. (8.18)
PROOF. Let n ∈ {1, . . . ,N}. Since the operator ΠP0 is stable for the L2-norm we have
|ΠP0R
n|6 |Rn|6
c(tn)− c(tn−1)
− ct(tn)
∣∣∣b̃
u(tn−1)−u(tn),c(tn)
)∣∣∣ .
Using a Taylor expansion and the Cauchy-Schwarz inequality, we get
c(tn)− c(tn−1)
− ct(tn)
∫ tn−1
|tn−1 − s| |ctt(s)|ds 6
(∫ tn−1
|ctt(s)|2 ds
On the other hand, since ∇c(tn)|∂Ω = 0, we deduce from (3.10) by integrating by parts
b̃(u(tn−1)−u(tn),c(tn)) = (u(tn−1)−u(tn)) ·∇c(tn).
Using a Taylor expansion and the Cauchy-Schwarz inequality, we get
|(u(tn−1)−u(tn)) ·∇c(tn)|6
k‖c‖L∞(0,T ;H1)
(∫ tm−1
|ut(s)|2 ds
|ΠP0R
(∫ tm−1
|ctt(s)|2 ds
k‖c‖L∞(0,T ;H1)
(∫ tm−1
|ut(s)|2 ds
Thanks to the stability of ΠP0 for the L
2-norm and to (8.1) we have
∣∣ΠP0
(s+λ )ηnh,c
)∣∣6C h(‖s‖L∞(0,T ;L2)+λ )‖c‖L∞(0,T ;H1).
The Cauchy-Schwarz inequality and (3.7) allow us to write
∣∣ΠP0ct(s)− Π̃P0ct(s)
∣∣ds 6C h
(∫ tm
‖ct(s)‖21+r ds
By plugging these estimates into definition (8.9) we get
k |Cnh,1|2 6 k2 ‖c‖2L∞(0,T ;H1)
|ut(s)|2 ds+ k2
|ctt(s)|2 ds
+ C h2
‖ct(s)‖21+r ds+C k h2‖c‖2L∞(0,T ;H1).
Study of a finite volume-finite element scheme for a nuclear transport model 23 of 26
Summing up the latter relation for n = 1 to m ∈ {1, . . . ,N} and using a similar work on Θ mh,1 we get
(8.16). Now let n ∈ {1, . . . ,N}. Using propositions 4.7 and 4.4 we have
‖ΠP0∆c(tn)−∆h
Π̃P0c(tn)
‖−1,h 6C h‖c‖L∞(0,T ;H2)
‖ΠP0 b̃(u(tn−1),c(tn))− b̃h(ΠRT0u(tn−1),Π̃P0c(tn))‖−1,h 6C h‖c‖L∞(0,T ;H1)‖u‖L∞(0,T ;H1+s).
Plugging these estimates into definition (8.9) and summing up from n = 1 to m, we obtain
‖Cnh,2‖2−1,h 6C h2
‖c‖2L∞(0,T ;H1) ‖u‖
L∞(0,T ;H1+s)+ ‖c‖
L∞(0,T ;H2)
A similar work on Θ mh,2 then leads to (8.17). We finally prove (8.18). Let m ∈ {1, . . . ,N}. On the one
hand, we have by (8.14)
|Pmh |6 (‖κmh,1‖∞ |ηmh,c|+ ‖κmh,2‖∞ |ηmh,θ |) |∇p(tm)|+ ‖κ(cmh ,θ mh )‖∞|∇hηmh,p|.
Using estimates (8.1)–(8.3) we get
|Pmh |6C h‖κ‖W1,∞((0,1)×(0,∞))(‖c‖L∞(0,T ;H1)+ ‖θ‖L∞(0,T ;H1)+ ‖p‖L∞(0,T ;H2)).
On the other hand definition (8.15) leads to
|Umh |6
∣∣Π̃RT0
f(tm)−ΠP0f(tm)
)∣∣+ |ηmh,u|+ |P
Using the stability of Π̃RT0 for the L
2-norm and proposition 3.3 we have
∣∣∣Π̃RT0
f(tm)−ΠP0f(tm)
)∣∣∣6 |f(tm)−ΠP0f(tm)|6C h‖f‖L∞(0,T ;H1).
Using moreover (8.3) we obtain
|Umh |6C h‖κ‖W1,∞((0,1)×(0,∞))(‖c‖L∞(0,T ;H1)+ ‖θ‖L∞(0,T ;H1))
+C (‖p‖L∞(0,T ;H2)+ ‖f‖L∞(0,T ;H1)+ ‖u‖L∞(0,T ;H1)).
We have proven (8.18).
Using the former proposition we are now able to estimate the discrete errors.
PROPOSITION 8.3 There exists some real k0 > 0 such that for any k < k0 and m ∈ {1, . . . ,N}
|εmh,c|2 + |εmh,θ |2 + k
‖εmh,c‖2h + ‖εmh,θ‖2h
6C (h2 + k2), (8.19)
|∇hεmh,p|+ |εmh,u|6C (h+ k). (8.20)
PROOF. Multiplying (8.4) by 2k εn+1h,c , we obtain
(εn+1h,c − ε
, 2k εn+1h,c
− 2Dc k (∆hεn+1h,c ,ε
h,c )+ 2k bh
unh,ε
h,c ,ε
+ 2k λ |εn+1h,c |
= 2k (Cn+1h+1 +C
h,2 ,ε
h,c )− 2k (sh, |ε
h,c |
2)− 2k bh
εnh,u,Π̃P0c(tn+1),ε
. (8.21)
24 of 26 C. CHOQUET AND S. ZIMMERMANN
Using an algebraic identity we have
(εn+1h,c − ε
,2k εn+1h,c
= |εn+1h,c |
2 −|εnh,c|2 + |εn+1h,c − ε
h,c|2.
We know by propositions 4.2 and 4.5 that
−2k (∆hεn+1h,c ,ε
h,c ) = 2k‖ε
h,c ‖
h , bh
unh,ε
h,c ,ε
We have
2k (sh, |εn+1h,c |
2)6 2k‖sh‖∞ |εn+1h,c |
6C k |εn+1h,c |
Using the Young inequality, we also write
∣∣∣2k (Cn+1h,1 ,ε
h,c )
∣∣∣6 2k |Cn+1h,1 | |ε
h,c |6 k λ |ε
h,c |
2 +C k |Cn+1h,1 |
|2k (Cn+1h,2 ,ε
h,c )|6 2k‖C
h,2 ‖−1,h‖ε
h,c ‖h 6 Dc
‖εn+1h,c ‖
h +C k‖Cn+1h,2 ‖
−1,h.
We are left with the term bh
εnh,u,Π̃P0c(tn+1),ε
. We have εnh,u∈ RT0. Using the divergence formula,
proposition 5.1 and (3.5), one easily checks that divεnh,u= 0. Thus we can apply proposition 4.3 to get
∣∣bh(εnh,u,Π̃P0c(tn+1),ε
h,c )
∣∣6C |εnh,u|‖Π̃P0c(tn+1)‖h‖ε
h,c ‖h. (8.22)
Let us first bound ‖Π̃P0c(tn+1)‖h. We have
‖Π̃P0c(tn+1)‖h 6 ‖Π̃P0c(tn+1)−ΠP0c(tn+1)‖h + ‖ΠP0c(tn+1)‖h.
Using an inverse inequality (see proposition 1.2 in [22]) and (3.7)
‖Π̃P0c(tn+1)−ΠP0c(tn+1)‖h 6
|Π̃P0c(tn+1)−ΠP0c(tn+1)|6C.
Moreover, according to [13] (p. 776), we have ‖ΠP0c(tn+1)‖h 6C‖c‖L∞(0,T ;H1). Thus ‖Π̃P0c(tn+1)‖h 6
C. We now estimate |εnh,u|. Using the stability of Π̃RT0 for the L
2-norm and the Cauchy-Schwarz
inequality in (8.7) we get
|εnh,u|6 ‖κ‖W1,∞((0,1)×(0,∞))
‖∇p‖L∞(0,T ;L2) (|ε
h,c |+ |ε
h,θ |)+ |∇hε
h,p |
+ |Un+1h |. (8.23)
We bound εn+1h,p as follows. Setting ψh = ε
h,p in (8.6) and using the Cauchy-Schwarz inequality, we get
κ(cn+1h ,θ
h )∇hε
h,p ,∇hε
6 ‖κ‖W1,∞((0,1)×(0,∞))‖∇p‖L∞(0,T ;L2) (|εn+1h,c |+ |ε
h,θ |) |∇hε
h,p |
+ |Pn+1h | |∇hε
h,p |+ |s− sh| |ε
h,p |.
The left-hand side is such that (κ(cn+1h ,θ
h )∇hε
h,p ,∇hε
h,p )> κin f |∇hε
h,p |
2. As for the right-hand
side, we have sh = ΠP0s and ε
h,p ∈ P
1 ∩L20. Thus, according to propositions 3.1 and 3.3
|s− sh| |εn+1h,p |6C h‖s‖L∞(0,T ;H1) |∇hε
h,p |.
Study of a finite volume-finite element scheme for a nuclear transport model 25 of 26
Finally |Pn+1h |6C h thanks to (8.18). Therefore we obtain
|∇hεn+1h,p |6C (h+ |ε
h,c |+ |ε
h,θ |). (8.24)
Let us plug this estimate into (8.23). Since |Un+1h |6C h thanks to (8.18), we get
|εnh,u|6C (h+ |ε
h,c |+ |ε
h,θ |). (8.25)
Now, plugging this bound into (8.23) and using the Young inequality, we obtain
∣∣bh(εnh,u,Π̃P0c(tn+1),ε
h,c )
∣∣ 6 C k (h+ |εn+1h,c |+ |ε
h,θ |)‖ε
h,c ‖h
‖εn+1h,c ‖
h +C k (h
2 + |εn+1h,c |
2 + |εn+1h,θ |
Now we have treated all the terms in (8.21). This equation implies
|εn+1h,c |
2 −|εnh,c|2 +Dc k‖εn+1h ‖
h 6C k
h2 + |εn+1h,c |
2 + |εn+1h,θ |
2 + |Cn+1h,1 |
2 + ‖Cn+1h,2 ‖
Let m ∈ {1, . . . ,N}. Let us sum up the latter estimate from n = 0 to m− 1. Thanks to (3.3)
|ε0h,c|= |Π̃P0c0 − c
h|= |Π̃P0c0 −ΠP0c0|6C h‖c‖L∞(0,T ;H2).
Using moreover estimates (8.16) and (8.17) we get
|εmh,c|2 +Dc k
‖εnh,c‖2h 6C k
(|εnh,c|2 + |εnh,θ |2)+C (h2 + k2).
Summing this relation with the one obtained by a similar work on (8.5) we obtain
|εmh,c|2 + |εmh,θ |2 + k
(Dc ‖εnh,c‖2h +Dθ ‖εnh,θ‖2h)6C k
(|εnh,c|2 + |εnh,θ |2)+C (h2 + k2).
Using a discrete Gronwall lemma (see lemma 5.2 in [22]) we get (8.19). Then (8.24) and (8.25) imply
(8.20).
By combining proposition 8.3 with estimates (8.1)-(8.3), we obtain finally the following result.
Theorem 8.1 There exists a real k0 > 0 such that for all k < k0 and m ∈ {1, . . . ,N}
|emh,c|2 + |emh,θ |2 + k
‖Π̃P0e
h,c‖2h + ‖Π̃P0e
h,θ‖2h
6C (h2 + k2) ,
|∇̃hemh,p|+ |emh,u|6C (h+ k).
REFERENCES
[1] ACHDOU, Y., BERNARDI, C. & COQUEL, F. (2003) A priori and a posteriori analysis of finite volume discretiza-
tions of Darcy’s equations, Numer. Math., 96(1) 17–42.
[2] AFIF, M.,AMAZIANE, B. (2002) Convergence of finite volume schemes for a degenerate convection-diffusion
equation arising in flow in porous media, Comput. Methods Appl., 191 5265–5286.
26 of 26 C. CHOQUET AND S. ZIMMERMANN
[3] ANGOT, P., DOLJŠI, V., FEISTAUER, M. & FELCMAN, J. (1998) Analysis of a combined barycentric finite
volume-nonconforming finite element method for nonlinear convection-diffusion problems, Appl. Math., 43(4)
263–310.
[4] BATEMAN, H. (1910) The solution of a system of differential equations occuring in the theory of radioactive
transformations, Proc. Camb. Philos. Soc., 15 423–427.
[5] BEAR, J. (1972) Dynamics of Fluids in Porous Media, New York: American Elsevier.
[6] BRENNER, S. C., SCOTT, L. R. & RIDGWAY, L. (2002) The mathematical theory of finite element methods, Texts
in applied mathematics, 16, New-York: Springer-Verlag.
[7] BREZZI, F. & FORTIN, M. (1991) Mixed and hybrid finite element methods, Springer Series in Computational
Mathematics, 15, New-York: Springer-Verlag.
[8] CHAVENT, G., JAFFRÉ, J. & ROBERTS, J. E. (1995) Mixed hybrid finite elements and cell-centered finite volumes
for two-phase flow in porous media. In: A. P. BOURGEAT AND AL, ed. Mathematical modeling of flow
through porous media, London: World Scientific, 100–114.
[9] CHEN, Z. (2000) Formulation and numerical methods of the black oil model in porous media, SIAM J. Numer.
Anal., 38(2) 489–514.
[10] CHOQUET C. (2004) Existence result for a radionuclide transport model with unbounded viscosity , J. Math. Fluid
Mech., 6(4) 365–388.
[11] DOUGLAS J., JR. & SPAGNUOLO, A. M. (2001) The approximation of nuclear contaminant transport in porous
media, J. Korean Math. Soc., 38 723–761.
[12] EWING, R. E., YUAN, Y. & LI, G. (1989) A time-discretization procedure for a mixed finite element approxi-
mation of contamination by incompressible nuclear waste in porous media, In: Mathematics for large scale
computing, Lecture Notes in Pure and Appl. Math., 120, 127–145.
[13] EYMARD, R., GALLOUËT, T. & HERBIN, R. (2000) The finite volume method (chapter III), In P.G. CIARLET
AND J. L. LIONS, eds Handbook of numerical analysis, 7, 762–835.
[14] GEISER, J. (2002) Numerical simulation of a model for transport and reaction of radionuclides, In: Algoritmy 2002,
Vysoké Tatry - Podbanské, Slovakia, Bratislava: Slovak university of technology.
[15] JOHN, V., MAUBACH, J. M. & TOBISKA, L. (1997) Nonconforming streamline-diffusion finite-element methods
for convection-diffusion problems, Numer. Math., 78(2) 165–188.
[16] KOVAL, E. J. (1963) A method for predicting the performance of unstable miscible displacements in heterogeneous
media, SPEJ trans. AIME, 228 145–154.
[17] LIN P. & YANG, D. (1998) An iterative perturbation method for the pressure equation in the simulation of miscible
displacement in porous media, SIAM J. Sci. Comput., 19(3) 893–911.
[18] MICHEL, A. (1999) Convergence of a finite volume scheme for a nonlinear convection-diffusion problem, In:
Proceedings of the Second international symposium on finite volumes for complex applications, Duisburg.
[19] MORTON, K. W. (1996) Numerical solution of convection-diffusion problems, In: Applied Mathematics and math-
ematical Computation, 12, London: Chapman and Hall.
[20] VERWER, J. G., BLOM, J. G. & HUNDSDORFER, W. (1996) An implicit-explicit approach for atmospheric trans-
port chemistry problems, Appl. Numer. Math., 20(1-2) 191–209.
[21] ZIMMERMANN, S. (2006) Stability of a finite volume scheme for incompressible fluids, preprint.
[22] ZIMMERMANN, S. (2006) Étude et implémentation de méthodes de volumes finis pour les fluides incompressibles.
Thesis (PhD). Blaise Pascal university (France).
