# Modified Method of Characteristics Combined with Finite Volume Element Methods for  Incompressible Miscible Displacement Problems in Porous Media

Authors: Sarvesh Kumar; Sangita Yadav
Journal: International Journal of Partial Differential Equations
 (2014)

---

## Abstract

The incompressible miscible displacement problem in porous media is modeled by a coupled system of two nonlinear
partial differential equations, the pressure-velocity equation and the concentration equation. In this paper, we present a mixed finite volume element
method (FVEM) for the approximation of the
pressure-velocity equation. Since modified method of characteristics (MMOC) minimizes the grid orientation effect, for the approximation of the concentration
equation, we apply a standard FVEM combined with MMOC.A priori error estimates in L
∞
(
L
2
) norm are derived for velocity, pressure  and concentration.
Numerical results are presented to substantiate the validity of the theoretical results.


---

## Body

## 1. Introduction

A mathematical model describing miscible displacement of one incompressible fluid by another in a horizontal porous medium reservoirΩ
⊂
ℝ
2 with boundary ∂
Ω of unit thickness over a time period of J
=
(
0
,
T
] is given by
(1)
u
=
-
κ
(
x
)
μ
(
c
)
∇
p
∀
(
x
,
t
)
∈
Ω
×
J
,
(2)
∇
·
u
=
q
∀
(
x
,
t
)
∈
Ω
×
J
,
(3)
ϕ
(
x
)
∂
c
∂
t
+
u
·
∇
c
-
∇
·
(
D
(
u
)
∇
c
)
=
g
(
x
,
t
,
c
)
=
(
c
~
-
c
)
q
∀
(
x
,
t
)
∈
Ω
×
J
,
with boundary conditions
(4)
u
·
n
=
0
∀
(
x
,
t
)
∈
∂
Ω
×
J
,
(5)
D
(
u
)
∇
c
·
n
=
0
∀
(
x
,
t
)
∈
∂
Ω
×
J
,
and initial condition
(6)
c
(
x
,
0
)
=
c
0
(
x
)
∀
x
∈
Ω
,
where x
=
(
x
1
,
x
2
)
∈
Ω, u
(
x
,
t
)
=
(
u
1
(
x
,
t
)
,
u
2
(
x
,
t
)
) and p
(
x
,
t
) are, respectively, the Darcy velocity and the pressure of the fluid mixture, c
(
x
,
t
) is the concentration of the fluid, c
~ is the concentration of the injected fluid, μ
(
c
) is the concentration dependent viscosity of the mixture, κ
(
x
) is the 2
×
2 permeability tensor of the medium, q
(
x
,
t
) is the external source/sink term that accounts for the effect of injection and production wells, and ϕ
(
x
) is the porosity of the medium. Further, D
(
u
)
=
D
(
x
,
u
) is the diffusion-dispersion tensor
(7)
D
(
u
)
=
ϕ
(
x
)
[
d
m
I
+
|
u
|
(
d
l
E
(
u
)
+
d
t
(
I
-
E
(
u
)
)
)
]
,
where d
m is the molecular diffusion, d
l and d
t are, respectively, the longitudinal and transverse dispersion coefficients, E
(
u
) is the tensor that projects onto u direction, whose i
jth component is given by
(8)
E
(
u
)
=
u
i
u
j
|
u
|
2
;
1
≤
i
,
j
≤
2
,
|
u
|
2
=
u
1
2
+
u
2
2
,
and I is the identity matrix of order 2. g
(
x
,
t
,
c
)
=
(
c
~
-
c
)
q is a known function representing sources denoted as g
(
c
) for convenience and c
0
(
x
) represents the initial concentration. For physical relevance 0
≤
c
0
(
x
)
≤
1 and n denotes the unit exterior normal to ∂
Ω. The following compatibility condition is imposed on the data:
(9)
∫
Ω
q
(
x
,
t
)
d
x
=
0
‍
∀
t
∈
J
,
which can be easily derived from (1)-(2) and (4). Equation (9) indicates that, for an incompressible flow with an impermeable boundary, the amount of injected fluid and the amount of fluid produced are equal. We also assume that the functions ϕ, μ, κ, and q are bounded; that is, there exist positive constants ϕ
*, ϕ
*, μ
*, μ
*, κ
*, κ
*, q
*, D
* such that
(10)
0
<
ϕ
*
≤
ϕ
(
x
)
≤
ϕ
*
,
0
<
μ
*
≤
μ
(
x
,
c
)
≤
μ
*
,
hhhhhhhh
0
<
κ
*
≤
κ
(
x
)
≤
κ
*
,
(11)
|
q
(
x
)
|
≤
q
*
,
D
(
x
,
u
)
≤
D
*
.The authors in [1–3] have discussed mathematical theory and existence of a unique weak solution of the above system (1)–(6) under suitable assumptions on the data. The pressure-velocity equation is elliptic type while the concentration equation is convection dominated diffusion type. Since, in the concentration equation only velocity is present, one would like to find a good approximation of the velocity. Therefore, for approximating velocity, it is natural to think of some mixed methods, which provide more accurate solution for the velocity compared to the standard finite element methods. Earlier, Douglas et al. [4, 5], Ewing and Wheeler [6], and Darlow et al. [7] have discussed the mixed finite element method for approximating the velocity as well as pressure and a standard Galerkin method for the concentration equation. They have also derived optimal error estimates in L
∞
(
L
2
) norm for the velocity and concentration. Recently, Kumar [8] has discussed a mixed and discontinuous FVEM for incompressible miscible displacement problems in porous media.Since the concentration equation is a convection dominated diffusion equation, the standard numerical methods fail to provide an accurate solution for the concentration and, therefore, suitable numerical methods have been proposed in the past for the approximation of the concentration equation. The standard numerical schemes fail to provide a physically relevant solution because most of these methods suffer from grid orientation effects. The other way to minimize the grid orientation effect is to use modified methods of characteristics (MMOC). Douglas and Russell [9] introduced and analyzed MMOC for the approximation of convection dominated diffusion equations. The authors in [10–12] studied MMOC combined with Galerkin finite element methods for incompressible miscible displacement problems.The basic idea behind the modified method of characteristics for approximating the concentration equation (3) is to set the hyperbolic part, that is, ϕ
(
∂
c
/
∂
t
)
+
u
·
∇
c, as a directional derivative.Set (see Figure1)
(12)
ψ
(
x
,
t
)
=
(
|
u
(
x
,
t
)
|
2
+
ϕ
(
x
)
2
)
1
/
2
=
(
u
1
(
x
,
t
)
2
+
u
2
(
x
,
t
)
2
+
ϕ
(
x
)
2
)
1
/
2
.
The characteristic direction with respect to the operator ϕ
(
∂
/
∂
t
)
+
u
·
∇ is the unit vector
(13)
s
(
x
,
t
)
=
(
u
1
(
x
,
t
)
,
u
2
(
x
,
t
)
,
ϕ
(
x
)
)
ψ
(
x
,
t
)
.
The directional derivative of the concentration c
(
x
,
t
) in the direction of s is given by
(14)
∂
c
∂
s
=
∂
c
∂
t
ϕ
(
x
)
ψ
(
x
,
t
)
+
u
·
∇
c
ψ
(
x
,
t
)
.
This implies that ψ
(
x
,
t
)
(
∂
c
/
∂
s
)
=
ϕ
(
x
)
(
∂
c
/
∂
t
)
+
u
·
∇
c.Figure 1
Direction ofψ
(
x
,
t
).Hence, (3) can be rewritten as
(15)
ψ
(
x
,
t
)
∂
c
∂
s
-
∇
·
(
D
(
u
)
∇
c
)
=
(
c
~
-
c
)
q
∀
(
x
,
t
)
∈
Ω
×
J
.Since (15) is in the form of heat equation, the behavior of the numerical solution of (15) should be better than (3) if the derivative term ∂
c
/
∂
s is approximated properly.We choose the same time steps for pressure and concentration for simplicity. However, the analysis can be extended to the case when different time steps are chosen for velocity and concentration through minor modifications.Let0
=
t
0
<
t
1
<
⋯
t
N
=
T be a given partition of the time interval [
0
,
T
] with the time step size Δ
t. For very small values of Δ
t, the characteristic direction starting from (
x
,
t
n
+
1
) crosses t
=
t
n at (see Figure 2)
(16)
x
ˇ
=
x
-
u
n
+
1
ϕ
(
x
)
Δ
t
,
where u
n
+
1
=
u
(
x
,
t
n
+
1
).Figure 2
An illustration of the definitionx
ˇ.This suggests us to approximate the characteristic directional derivative att
=
t
n
+
1 as
(17)
∂
c
∂
s
|
t
=
t
n
+
1
≈
c
n
+
1
-
c
(
x
ˇ
,
t
n
)
Δ
s
=
c
n
+
1
-
c
(
x
ˇ
,
t
n
)
(
(
x
-
x
ˇ
)
2
+
(
t
n
+
1
-
t
n
)
2
)
1
/
2
,
where c
n
+
1
=
c
(
x
,
t
n
+
1
).Using (16), we obtain
(18)
ψ
(
x
,
t
)
∂
c
∂
s
|
t
=
t
n
+
1
≈
ϕ
(
x
)
c
n
+
1
-
c
ˇ
n
Δ
t
,
where c
˘
n
=
c
(
x
ˇ
,
t
n
).Compared to the conforming finite element methods (FEM), the finite volume methods are conservative in nature, and, hence, they preserve the physical conservative properties. In a mixed FVEM, one uses two different kinds of grids: a primal grid and a dual grid (see Figures3 and 4). Mixed FVEM can also be thought of as a Petrov-Galerkin method. The analysis of these methods is based on the tools borrowed from the mixed FEM. Using a transfer operator which maps the trial space to the test space, the mixed covolume methods can be put in the framework of mixed FEM. This transfer operator plays a vital role in deriving the optimal error estimates. Earlier, Chou et al. [13, 14] have discussed and analyzed mixed covolume or FVEM for the second-order linear elliptic problems in two-dimensional domains. The standard FVEM can also be considered as a Petrov-Galerkin finite element method in which the trial space is chosen as C
0 piecewise linear polynomials on the finite element partition of the domain and the test space, as piecewise constants over the control volumes are to be defined in Section 2. For more details on finite volume methods, we refer to [15–18] and references there in. In this paper, we present a mixed FVEM for approximating the pressure-velocity equations (1)-(2) and (4) and a standard FVEM combined with MMOC for the approximation of the concentration equations (15) and (5)-(6). This paper is organized as follows. In Section 2, FVE approximation procedure is discussed. A priori error estimates for velocity, pressure, and concentration are presented in Section 3. Finally in Section 4, the numerical procedure is discussed and some numerical experiments are presented.Figure 3
Primal grid𝒯
h and dual grid 𝒯
h
*.Figure 4
Primal grid𝒯
h and dual grid 𝒱
h
*.


## 2. Finite Volume Element Approximation

LetU
=
{
v
∈
H
(
div
⁡
;
Ω
)
:
v
·
n
=
0
on
⁡
∂
Ω
}. Note that (1)-(2) with boundary condition (4) has a solution for pressure, which is only unique up to an additive constant. The nonuniqueness of (1)-(2) may be avoided by considering the following quotient space:
(19)
W
=
L
2
(
Ω
)
ℝ
.
Multiply (1) and (2) by v
∈
U and w
∈
W, respectively, and integrate over Ω. A use of Green’s formula and v
·
n
=
0
on
⁡
∂
Ω yields the following weak formulation: Find (
u
,
p
)
:
J
¯
→
U
×
W satisfying
(20)
(
κ
-
1
μ
(
c
)
u
,
v
)
-
(
∇
·
v
,
p
)
=
0
∀
v
∈
U
,
(
∇
·
u
,
w
)
=
(
q
,
w
)
∀
w
∈
W
.
We use a mixed FVEM for the simultaneous approximation of velocity and pressure in (1)-(2) and a standard FVEM for the approximation of the concentration in (15). For this purpose, we introduce three kinds of grids: one primal grid and two dual grids.Let𝒯
h
=
{
T
} be a regular, quasiuniform partition of the domain Ω
- into closed triangles T; that is, Ω
¯
=
∪
T
∈
𝒯
h
T
¯. Let h
T
=
diam
⁡
(
T
) and h
=
max
⁡
T
∈
𝒯
h
h
T. Let the trial function spaces U
h and W
h associated with the approximation of velocity and pressure, respectively, be the lowest order Raviart-Thomas space for triangles defined by
(21)
U
h
=
{
v
h
∈
U
:
v
h
∣
T
=
(
a
+
b
x
,
c
+
b
y
)
∀
T
∈
𝒯
h
}
,
W
h
=
{
w
h
∈
W
:
w
h
∣
T
is
a
constant
∀
T
∈
𝒯
h
}
.
Let us define the discrete norm for v
h
=
(
v
h
1
,
v
h
2
)
∈
U
h as
(22)
∥
v
h
∥
1
,
h
2
=
∥
v
h
∥
(
L
2
(
Ω
)
)
2
2
+
|
v
h
|
1
,
h
2
,
where |
v
h
|
1
,
h
2
=
∑
T
∈
𝒯
h
‍
∥
∇
υ
h
1
∥
0
,
T
2
+
∥
∇
υ
h
2
∥
0
,
T
2. For υ
h
∈
U
h, it is straight forward to check that
(23)
∥
v
h
∥
1
,
h
≤
C
∥
v
h
∥
H
(
div
⁡
;
Ω
)
,
where C is a constant independent of h. For v
h
∈
U
h the inequality
(24)
∥
v
h
∥
(
L
∞
(
Ω
)
)
2
≤
C
(
log
⁡
1
h
)
1
/
2
∥
v
h
∥
1
,
h
also holds true when Ω is in ℝ
2 and the triangulation 𝒯
h is quasiuniform and can be proved using the same arguments as in the proof of Lemma 4 in [19, pp. 67].The dual grid𝒯
h
* consists of interior quadrilaterals and boundary triangles, which are constructed by joining barycenter to the vertices. For the construction of the dual grid 𝒯
h
* we refer to [14]. In general, let T
M
* denote the dual element corresponding to the midside node M. The union of all the dual elements/control volume elements forms a partition 𝒯
h
* of Ω
-. The test space V
h is defined by
(25)
V
h
=
{
v
h
∈
(
L
2
(
Ω
)
)
2
:
v
h
∣
T
M
*
is
a
constant
vector
j
j
j
j
(
L
2
(
Ω
)
)
2
∀
T
M
*
∈
𝒯
h
*
and
v
h
·
n
=
0
on
⁡
∂
Ω
}
,
where T
M
* denotes the dual element corresponding to the midside node M. For connecting our trial and test spaces, we define a transfer operator γ
h
:
U
h
→
V
h by
(26)
γ
h
v
h
(
x
)
=
∑
i
=
1
N
m
v
h
(
M
i
)
χ
i
*
(
x
)
∀
x
∈
Ω
,
where N
m is the total number of midside nodes and χ
i
*’s are the scalar characteristic functions corresponding to the control volume T
M
j
* defined by
(27)
χ
i
*
(
x
)
=
{
1
,
if
x
∈
T
M
i
*
,
0
,
elsewhere
.
Multiplying (1) by γ
h
v
h
∈
V
h, integrating over the control volumes T
M
*
∈
𝒯
h
*, applying the Gauss’s divergence theorem, and summing up over all the control volumes, we obtain
(28)
(
κ
-
1
μ
(
c
)
u
,
γ
h
v
h
)
-
∑
i
=
1
N
m
v
h
(
M
i
)
‍
·
∫
T
M
i
*
p
n
T
M
i
*
d
s
=
0
∀
v
h
∈
U
h
,
where n
T
M
i
* denotes the outward normal vector to the boundary of T
M
i
*. Set
(29)
b
(
γ
h
v
h
,
w
h
)
=
-
∑
i
=
1
N
m
v
h
(
M
i
)
‍
·
∫
∂
T
M
i
*
w
h
n
T
M
i
*
d
s
ggggggggg
∀
v
h
∈
U
h
,
∀
w
h
∈
W
h
.
Then, the mixed FVE approximation corresponding to (1)-(2) can be written as follows: find (
u
h
,
p
h
)
:
J
¯
→
U
h
×
W
h such that, for t
∈
(
0
,
T
],
(30)
(
κ
-
1
μ
(
c
h
)
u
h
,
γ
h
v
h
)
+
b
(
γ
h
v
h
,
p
h
)
=
0
∀
v
h
∈
U
h
,
(
∇
·
u
h
,
w
h
)
=
(
q
,
w
h
)
∀
w
h
∈
W
h
,
where c
h is an approximation to c obtained from (34).Now, we introduce a dual mesh𝒱
h
* based on 𝒯
h which will be used for the approximation of the concentration equation. For construction we refer to [18].For applying the standard FVEM to approximate the concentration, we define the trial spaceM
h on 𝒯
h and the test space L
h on 𝒱
h
* as follows:
(31)
M
h
=
{
z
h
∈
C
0
(
Ω
-
)
:
z
h
∣
T
∈
P
1
(
T
)
∀
T
∈
𝒯
h
}
,
L
h
=
{
w
h
∈
L
2
(
Ω
)
:
w
h
∣
V
P
*
is
a
constant
∀
V
P
*
∈
𝒱
h
*
}
l
,
where V
P
* is the control volume associated with node P. Again, we define a transfer function Π
h
*
:
M
h
→
L
h by
(32)
Π
h
*
z
h
(
x
)
=
∑
j
=
1
N
h
z
h
(
P
j
)
χ
j
(
x
)
∀
x
∈
Ω
,
where N
h is the total number of vertices and χ
j’s are the characteristic functions corresponding to the control volume V
P
j
* given by
(33)
χ
j
(
x
)
=
{
1
,
if
x
∈
V
P
j
*
,
0
,
elsewhere
.
For obtaining a finite volume approximation c
h to the concentration c, we multiply (15) by Π
h
*
z
h
∈
M
h, integrate over the control volumes V
P
j
*
∈
𝒱
h
*, and apply the Gauss’s divergence theorem. Then we sum up over the control volumes to obtain the FVE approximation c
h corresponding to the concentration equation (15) as a solution c
h
:
J
¯
→
M
h such that for t
∈
(
0
,
T
],
(34)
(
ψ
∂
c
h
∂
s
,
Π
h
*
z
h
)
+
a
h
(
u
h
;
c
h
,
z
h
)
+
(
c
h
q
,
Π
h
*
z
h
)
=
(
c
~
q
,
Π
h
*
z
h
)
∀
z
h
∈
M
h
,
c
h
(
0
)
=
c
0
,
h
,
where c
0
,
h is an approximation to c
0 to be defined later and the bilinear form a
h
(
v
;
·
,
·
) is defined by
(35)
a
h
(
v
;
χ
,
ϕ
h
)
=
-
∑
j
=
1
N
h
∫
∂
V
P
j
*
(
D
(
v
)
∇
χ
·
n
P
j
)
Π
h
*
ϕ
h
d
s
,
n
P
j being the unit outward normal to the boundary of V
P
j
* with v
∈
U, χ
∈
H
1
(
Ω
), and ϕ
h
∈
M
h.Remark 1.
The three grids are introduced each for the pressure, velocity, and concentration variables. This is to balance the number of unknowns and the equations in the coupled systems (30) and (34).To approximate the concentration at any time, sayt
n
+
1, we use the approximation to the velocity at the previous time step. The fully discrete scheme corresponding to (30) and (34) is defined as follows. For n
=
0,1
…
N, find (
c
h
n
,
p
h
n
,
u
h
n
)
∈
M
h
×
W
h
×
U
h such that
(36)
c
h
0
=
R
h
c
(
0
)
,
(37)
(
κ
-
1
μ
(
c
h
n
)
u
h
n
,
γ
h
v
h
)
+
b
(
γ
h
v
h
,
p
h
n
)
=
0
∀
v
h
∈
U
h
,
(38)
(
∇
·
u
h
n
,
w
h
)
=
(
q
n
,
w
h
)
∀
w
h
∈
W
h
,
(39)
(
ϕ
c
h
n
+
1
-
c
^
h
n
Δ
t
,
Π
h
*
χ
h
)
+
a
h
(
u
h
n
;
c
h
n
+
1
,
χ
h
)
+
(
q
n
+
1
c
h
n
+
1
,
Π
h
*
χ
h
)
=
(
q
n
+
1
c
~
n
+
1
,
Π
h
*
χ
h
)
∀
χ
h
∈
M
h
,
where c
^
h
n
=
c
h
(
x
^
,
t
n
)
=
c
h
(
x
-
(
u
h
n
/
ϕ
)
Δ
t
,
t
n
) and R
h
c is a projection of c onto M
h which will be defined in (45).Note that in (17). we use the following notation for the exact velocity:
(40)
c
ˇ
n
=
c
(
x
ˇ
,
t
n
)
=
c
(
x
-
u
n
+
1
ϕ
Δ
t
,
t
n
)
.


## 3. Error Estimates

### 3.1. Error Estimates for Velocity

For a givenc, the auxiliary functions (
u
~
h
,
p
~
h
)
:
[
0
,
T
]
→
U
h
×
W
h are defined as follows:
(41)
(
κ
-
1
μ
(
c
)
u
~
h
,
v
h
)
-
(
∇
·
v
h
,
p
~
h
)
=
0
∀
v
h
∈
U
h
,
(
∇
·
u
~
h
,
w
h
)
=
(
q
,
w
h
)
∀
w
h
∈
W
h
.
For a proof of existence and uniqueness of the discrete solution of (41), we refer to [20, pp. 52]. For u
h and u
~
h, the following error estimates can be obtained (see [14]):
(42)
∥
u
h
-
u
~
h
∥
(
L
2
(
Ω
)
)
2
≤
C
(
h
∥
u
~
h
∥
(
L
2
(
Ω
)
)
2
+
∥
c
-
c
h
∥
∥
u
~
h
∥
(
L
∞
(
Ω
)
)
2
)
.
Now, since concentration depends on the velocity and vice versa, to derive the error estimates for the concentration, we also need error estimates for the velocity. For elliptic problems, the authors in [14] have derived error estimates for mixed covolume method by using Raviart-Thomas projection and L
2 projection. In the similar way, for a given c, the following theorem can be shown but the proof is long so we omit it here.Theorem 2.
Assume that the triangulation𝒯
h is quasiuniform. Let (
u
,
p
) and (
u
h
,
p
h
), respectively, be the solutions of (1)-(2) and (30). Then, there exists a positive constant C independent of h but dependent on the bounds of κ
-
1 and μ such that
(43)
∥
u
-
u
h
∥
(
L
2
(
Ω
)
)
2
+
∥
p
-
p
h
∥
≤
C
[
∥
c
-
c
h
∥
+
h
(
∥
u
∥
(
H
1
(
Ω
)
)
2
+
∥
p
∥
1
)
]
,
(44)
∥
∇
·
(
u
-
u
h
)
∥
≤
C
h
∥
∇
·
u
∥
1
.


### 3.2. Error Estimates for Concentration

We writec
-
c
h
=
(
c
-
R
h
c
)
+
(
R
h
c
-
c
h
), where R
h
:
H
1
(
Ω
)
→
M
h is the projection of c defined by
(45)
A
(
u
;
c
-
R
h
c
,
χ
h
)
=
0
∀
χ
h
∈
M
h
,
where
(46)
A
(
u
;
ϕ
,
χ
h
)
=
a
h
(
u
;
ϕ
,
χ
h
)
+
(
q
ϕ
,
χ
h
)
+
(
λ
ϕ
,
χ
h
)
∀
χ
h
∈
M
h
.
The function λ will be chosen such that the coercivity of A
(
u
;
·
,
·
) is assured.We use frequently the following trace inequality [21, pp. 417]: for w
∈
H
1
(
T
),
(47)
∥
w
∥
∂
T
2
≤
C
(
h
T
-
1
∥
w
∥
T
2
+
h
T
|
w
|
1
,
T
2
)
,
where ∥
w
∥
∂
T
2
=
∫
∂
T
‍
w
2
d
s. Further, we need the following inverse inequalities (see [22, pp. 141]): ∀
χ
h
∈
M
h
(48)
∥
χ
h
∥
j
,
∞
≤
C
h
-
1
∥
χ
h
∥
1
,
j
=
0,1
,
g
g
g
g
g
∥
χ
h
∥
1
≤
C
h
-
1
∥
χ
h
∥
.
Using the properties of Π
h
* operator, it is easy to see that, for T
∈
𝒯
h and ϕ
h
∈
M
h the following holds true:
(49)
∫
T
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
‍
=
0
,
∫
∂
T
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
s
=
0
.
Now, using (49), we have the following lemma [17, pp. 317].Lemma 3.
There exists a positive constantC independent of h such that
(50)
j
j
j
j
j
j
j
|
ϵ
h
(
ϕ
,
χ
)
|
≤
C
h
i
+
j
|
ϕ
|
W
p
i
|
χ
|
W
q
j
∀
χ
∈
M
h
w
i
t
h
i
,
j
=
0,1
,
1
p
+
1
q
=
1
,
where
(51)
ϵ
h
(
ϕ
,
χ
h
)
=
(
ϕ
,
χ
h
)
-
(
ϕ
,
Π
h
*
χ
h
)
∀
χ
h
∈
M
h
.Also note that, by the usual interpolation theory, the operatorΠ
h
* has the following approximation property [23, pp. 466]:
(52)
jj
∥
χ
-
Π
h
*
χ
∥
0
,
k
≤
C
h
β
|
χ
|
s
,
k
,
0
≤
β
≤
s
≤
1
,
1
≤
k
≤
∞
.
We also recall two well-known results (Lemmas 4 and 5), which will be used in the the proof of Lemmas 7–9.Lemma 4 (see [24, pp. 240]).
The operatorΠ
h
* has the following properties.(i)
With∥
|
ϕ
h
|
∥
=
(
ϕ
h
,
Π
h
*
ϕ
h
)
1
/
2, the norms ∥
|
·
|
∥ and ∥
·
∥ are equivalent on U
h; that is, there exist positive constants C
7 and C
8, independent of h, such that
(53)
C
7
∥
ϕ
h
∥
≤
∥
|
ϕ
h
|
∥
≤
C
8
∥
ϕ
h
∥
∀
ϕ
h
∈
M
h
.
(ii)
Π
h
* is stable with respect to the L
2 norm; that is, there exists a positive constant C independent of h such that
(54)
∥
Π
h
*
χ
h
∥
≤
C
∥
χ
h
∥
∀
χ
h
∈
M
h
.Lemma 5 (see [25, pp. 1871]).
Assume thatχ
h
,
ϕ
h
∈
M
h. Then one has
(55)
a
h
(
u
;
χ
h
,
ϕ
h
)
=
a
(
u
;
χ
h
,
ϕ
h
)
+
∑
T
∈
𝒯
h
∫
∂
T
(
D
(
u
)
∇
χ
h
·
n
)
(
Π
h
*
ϕ
h
-
ϕ
h
)
d
s
‍
‍
+
∑
T
∈
𝒯
h
∫
T
∇
·
(
D
(
u
)
∇
χ
h
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
.
Moreover, for χ
h
∈
M
h,
(56)
a
h
(
u
;
χ
h
,
χ
h
)
≥
a
(
u
;
χ
h
,
χ
h
)
-
C
h
|
χ
h
|
1
2
.Remark 6.
Note that (55) also holds true for χ
∈
H
1
(
Ω
).SinceD is uniformly positive definite, we obtain following from (56):
(57)
a
h
(
u
;
χ
h
,
χ
h
)
≥
(
α
-
C
h
)
|
χ
h
|
1
2
.
Choose h
0
>
0 such that, for 0
<
h
<
h
0, (
α
-
C
h
)
=
α
0
>
0 and hence
(58)
a
h
(
u
h
;
χ
h
,
χ
h
)
≥
α
0
|
χ
h
|
1
2
∀
χ
h
∈
M
h
.
Now, we derive the error bound in H
1 and L
2 norms for c
-
R
h
c. Let I
h be the continuous interpolant onto M
h satisfying the following approximation properties. For ϕ
∈
H
k
+
1
(
Ω
) with k
=
0,1, we have [22]
(59)
∥
ϕ
-
I
h
ϕ
∥
j
≤
h
k
+
1
-
j
∥
ϕ
∥
k
+
1
j
=
0,1
.
Moreover, if ϕ
∈
W
2
,
∞
(
Ω
), then
(60)
∥
ϕ
-
I
h
ϕ
∥
1
,
∞
≤
C
h
log
⁡
(
1
h
)
∥
ϕ
∥
2
,
∞
.Lemma 7.
There exists a positive constantC independent of h such that
(61)
∥
c
-
R
h
c
∥
1
≤
C
h
∥
c
∥
2
,
provided c
∈
H
2
(
Ω
), for a.e. t
∈
(
0
,
T
].Proof.
The coercivity and boundedness ofA
(
u
;
·
,
·
) with (45) yield
(62)
∥
I
h
c
-
R
h
c
∥
1
2
≤
C
A
(
u
;
I
h
c
-
R
h
c
,
I
h
c
-
R
h
c
)
≤
C
A
(
u
;
I
h
c
-
c
,
I
h
c
-
R
h
c
)
≤
C
∥
c
-
I
h
c
∥
1
∥
I
h
c
-
R
h
c
∥
1
,
and hence
(63)
∥
I
h
c
-
R
h
c
∥
1
≤
C
∥
c
-
I
h
c
∥
1
,
where C depends on the bound of D
(
u
) given in (11). Combine the estimates (63) and (59) and use the triangle inequality to complete the proof.
For deriving theL
2 error bounds for c
-
R
h
c, we need the following lemma.
Lemma 8.Let
ϵ
a
(
u
;
ϕ
,
χ
h
)
=
a
(
u
;
ϕ
,
χ
h
)
-
a
h
(
u
;
ϕ
,
χ
h
)
. There exists a positive constant
C
such that, for
ϕ
∈
H
1
(
Ω
)
and
χ
h
∈
M
h,
(64)
|
ϵ
a
(
u
;
c
-
R
h
c
,
ϕ
h
)
|
≤
C
h
2
(
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
+
∥
c
∥
2
)
|
ϕ
h
|
1
∀
ϕ
h
∈
M
h
.
Using (55) (see also Remark 6), we find that
(65)
|
ϵ
a
(
u
;
c
-
R
h
c
,
ϕ
h
)
|
≤
|
∑
T
∈
𝒯
h
∫
T
∇
·
(
D
(
u
)
∇
(
c
-
R
h
c
)
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
‍
|
+
|
∑
T
∈
𝒯
h
∫
∂
T
(
D
(
u
)
∇
(
c
-
R
h
c
)
·
n
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
s
‍
|
=
J
1
+
J
2
.
To bound J
1, first we use the fact that R
h
c is linear on each triangle T to obtain
(66)
J
1
=
|
∑
T
∈
𝒯
h
∫
T
∇
·
(
D
(
u
)
∇
(
c
-
R
h
c
)
)
g
g
g
g
g
g
g
g
∑
T
∈
𝒯
h
×
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
|
=
|
∑
T
∈
𝒯
h
∫
T
(
∇
·
(
D
(
u
)
∇
c
)
-
(
∇
·
D
(
u
)
)
·
∇
R
h
c
)
g
g
g
g
g
g
g
g
∑
T
∈
𝒯
h
×
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
|
.
Now use (3), (49), and (50) to obtain
(67)
J
1
≤
|
∑
T
∈
𝒯
h
∫
T
(
-
g
+
ψ
∂
c
∂
s
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
|
+
|
∑
T
∈
𝒯
h
∫
T
[
(
∇
·
D
(
u
)
-
(
∇
·
D
(
u
)
)
T
)
·
∇
R
h
c
]
g
g
g
g
g
g
g
g
∑
T
∈
𝒯
h
×
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
|
≤
C
h
2
(
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
+
∥
c
∥
2
)
|
ϕ
h
|
1
,
where (
∇
·
D
(
u
)
)
T denotes the average value of ∇
·
D
(
u
) on triangle T.
Based on the analysis in [25, pp. 1873], we estimate J
2 as follows. Note that an appeal to the continuity of ∇
c
·
n with (49) yields
(68)
J
2
=
|
∑
T
∈
𝒯
h
∫
∂
T
(
(
D
-
D
¯
T
)
∇
(
c
-
R
h
c
)
·
n
)
‍
g
g
g
g
g
g
g
g
∑
T
∈
𝒯
h
×
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
s
|
,
where D
=
D
(
u
) and D
¯
T is a function such that, for any edge of a triangle T
∈
𝒯
h,
(69)
D
¯
T
(
x
)
=
D
(
x
c
)
,
x
∈
E
,
and x
c is the midpoint of E. Since |
D
(
x
)
-
D
¯
T
|
≤
C
h
T
∥
D
∥
1
,
∞, we use trace inequalities (47) and (61) to arrive at
(70)
J
2
≤
C
h
|
∑
T
∈
𝒯
h
∫
∂
T
(
∇
(
c
-
R
h
c
)
·
n
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
s
‍
|
≤
C
h
(
∑
T
∈
𝒯
h
∫
∂
T
|
(
∇
(
c
-
R
h
c
)
·
n
)
|
2
‍
)
1
/
2
×
(
∑
T
∈
𝒯
h
∫
∂
T
|
ϕ
h
-
Π
h
*
ϕ
h
|
2
‍
)
1
/
2
≤
C
h
(
h
T
-
1
/
2
∥
c
-
R
h
c
∥
1
+
h
T
1
/
2
∥
c
∥
2
)
×
(
h
T
-
1
/
2
∥
ϕ
h
-
Π
h
*
ϕ
h
∥
+
h
T
1
/
2
|
ϕ
h
|
1
)
≤
C
h
2
∥
c
∥
2
|
ϕ
h
|
1
.
Substitute the estimates of J
1 and J
2 in (65) to complete  the rest of the proof.Lemma 9.
There exists a positive constantC independent of h such that
(71)
∥
c
-
R
h
c
∥
≤
C
h
2
(
∥
c
∥
2
+
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
)
,
provided c
∈
H
2
(
Ω
)
,
u
·
∇
c
∈
H
1
(
Ω
), and ∂
c
/
∂
t
∈
H
1
(
Ω
) for t
∈
(
0
,
T
] a.e.Proof.
To obtain optimalL
2 error estimates for c
-
R
h
c, we now appeal to Aubin-Nitsche duality argument. Let ϕ
∈
H
2
(
Ω
) be a solution of the following adjoint problem:
(72)
-
∇
·
(
D
(
u
)
∇
ϕ
+
u
ϕ
)
+
λ
ϕ
=
c
-
R
h
c
in
Ω
,
(
D
(
u
)
∇
ϕ
+
u
ϕ
)
·
n
=
0
on
⁡
∂
Ω
,
which satisfies the elliptic regularity condition
(73)
∥
ϕ
∥
2
≤
C
∥
c
-
R
h
c
∥
.
Multiply the above equation by c
-
R
h
c and integrate over Ω. An integration by parts and a use of (45) yield
(74)
∥
c
-
R
h
c
∥
2
=
a
(
u
;
ϕ
,
c
-
R
h
c
)
-
(
u
·
∇
ϕ
,
c
-
R
h
c
)
-
(
∇
·
u
ϕ
,
c
-
R
h
c
)
+
λ
(
ϕ
,
c
-
R
h
c
)
=
[
a
(
u
;
c
-
R
h
c
,
ϕ
-
ϕ
h
)
.
+
(
u
·
∇
(
c
-
R
h
c
)
,
ϕ
-
ϕ
h
)
.
+
λ
(
c
-
R
h
c
,
ϕ
-
ϕ
h
)
]
.
+
ϵ
a
(
u
;
c
-
R
h
c
,
ϕ
h
)
∀
ϕ
h
∈
M
h
=
I
1
+
I
2
.
For I
1, use (61) to find that
(75)
|
I
1
|
=
|
a
(
u
;
c
-
R
h
c
,
ϕ
-
ϕ
h
)
ggl
+
(
u
·
∇
(
c
-
R
h
c
)
,
ϕ
-
ϕ
h
)
ggl
+
λ
(
c
-
R
h
c
,
ϕ
-
ϕ
h
)
|
≤
C
∥
c
-
R
h
c
∥
1
∥
ϕ
-
ϕ
h
∥
1
≤
C
h
∥
c
∥
2
∥
ϕ
-
ϕ
h
∥
1
.
The bound for I
2 follows from Lemma 8 and hence
(76)
|
I
2
|
≤
C
h
2
(
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
+
∥
c
∥
2
)
|
ϕ
h
|
1
.
Substitute (75) and (76) in (74) to find that
(77)
∥
c
-
R
h
c
∥
2
≤
C
[
h
∥
c
∥
2
∥
ϕ
-
ϕ
h
∥
1
|
ψ
∂
c
∂
s
|
1
h
h
h
h
+
h
2
(
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
+
∥
c
∥
2
)
|
ϕ
h
|
1
]
.
Now choose ϕ
h
=
I
h
ϕ in (77). Then use elliptic regularity condition (73) with (59) to obtain
(78)
∥
c
-
R
h
c
∥
≤
C
h
2

(
∥
c
∥
2
+
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
)

,
and this completes the proof.A use of inverse inequalities (48), (59), (60), and (63) yields
(79)
∥
R
h
c
∥
1
,
∞
≤
C
∥
c
∥
2
,
∞
.Lemma 10.
There exists a positive constantC such that ∀
θ
∈
M
h,
(80)
|
a
h
(
u
;
R
h
c
,
θ
)
-
a
h
(
u
h
;
R
h
c
,
θ
)
|
≤
C
(
∥
u
-
u
h
∥
(
L
2
(
Ω
)
)
2
+
h
∥
∇
·
(
u
-
u
h
)
∥
)
|
θ
|
1
.Proof.
Note that(81)
|
a
h
(
u
;
R
h
c
,
θ
)
-
a
h
(
u
h
;
R
h
c
,
θ
)
|
=
|
∑
i
=
1
N
h
∫
∂
V
i
*
(
D
(
u
)
-
D
(
u
h
)
)
∇
R
h
c
·
n
i
Π
h
*
θ
d
s
‍
|
=
|
∑
T
∈
𝒯
h
K
T
‍
|
,
where K
T
=
∑
l
=
1
3
‍
∫
∂
V
l
*
∩
T
‍
(
D
(
u
)
-
D
(
u
h
)
)
∇
R
h
c
·
n
l
θ
l
d
s and θ
l
=
θ
(
P
l
); see Figure 5. For each triangle T, K
T can be written as
(82)
K
T
=
∑
l
=
1
3
∫
M
l
B
(
D
(
u
)
-
D
(
u
h
)
)
∇
R
h
c
·
n
l
‍
‍
hhhhhh
×
(
θ
l
+
1
-
θ
l
)
d
s
(
θ
4
=
θ
1
)
.
Using the Cauchy-Schwarz inequality and (79), we obtain
(83)
K
T
≤
∑
l
=
1
3
|
θ
l
+
1
-
θ
l
|
‍
×
∫
M
l
B
|
(
D
(
u
)
-
D
(
u
h
)
)
∇
R
h
c
·
n
l
|
d
s
‍
≤
C
∑
l
=
1
3
|
θ
l
+
1
-
θ
l
|
∥
D
(
u
)
-
D
(
u
h
)
∥
(
L
2
(
M
l
B
)
)
2
×
2
‍
×
(
meas
⁡
(
M
l
B
)
)
1
/
2
.
It has been proved in [26, pp. 332] that the matrix D is uniformly Lipschitz; that is, there exists a constant C such that for u and v
∈
(
L
2
(
Ω
)
)
2,
(84)
∥
D
(
u
)
-
D
(
v
)
∥
(
L
2
(
Ω
)
)
2
×
2
≤
C
∥
u
-
v
∥
(
L
2
(
Ω
)
)
2
.
A use of the trace inequalities (47) and (84) yields
(85)
K
T
≤
C
∑
l
=
1
3
|
θ
l
+
1
-
θ
l
|
∥
u
-
u
h
∥
(
L
2
(
M
l
B
)
)
2
h
T
1
/
2
‍
≤
C
h
T
1
/
2
∑
l
=
1
3
‍
|
θ
l
+
1
-
θ
l
|
[
h
T
-
1
/
2
∥
u
-
u
h
∥
T
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
l
+
h
T
1
/
2
∥
∇
·
(
u
-
u
h
)
∥
T
]
.
Now, using Taylor series expansion, we find that
(86)
|
θ
l
+
1
-
θ
l
|
≤
h
T
[
|
∂
θ
∂
x
|
+
|
∂
θ
∂
y
|
]
≤
[
(
|
∂
θ
∂
x
|
2
+
|
∂
θ
∂
y
|
2
)
h
T
2
]
1
/
2
=
|
θ
|
1
,
h
,
T
,
l
=
1,2
,
3
.
Let |
ϕ
h
|
1
,
h
=
(
∑
T
∈
𝒯
h
‍
|
ϕ
h
|
1
,
h
,
T
2
)
1
/
2. Since θ is linear on each triangle, seminorms |
·
| and |
·
| are identical.
Substitute (86) in (85) and use the estimates for K
T to arrive at
(87)
|
a
h
(
u
;
R
h
c
,
θ
)
-
a
h
(
u
h
;
R
h
c
,
θ
)
|
≤
C
(
∥
u
-
u
h
∥
(
L
2
(
Ω
)
)
2
+
h
∥
∇
·
(
u
-
u
h
)
∥
)
|
θ
|
1
,
and this completes the rest of the proof.Figure 5
A triangular partition.Now, we prove our main theorem.Theorem 11.
Letc
n and c
h
n be the solutions of (3) and (39) at t
=
t
n, respectively, and let c
h
(
0
)
=
c
0
,
h
=
R
h
c
(
0
). Further assume that Δ
t
=
O
(
h
). Then, for sufficiently small h, there exists a positive constant C
(
T
) independent of h but dependent on the bounds of κ
-
1 and μ such that
(88)
max
⁡
0
≤
n
≤
N
∥
c
n
-
c
h
n
∥
(
L
2
(
Ω
)
)
2
2
≤
C
[
h
4
(
∥
(
c
-
c
~
)
q
∥
L
∞
(
0
,
T
;
H
1
)
2
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
ffffffffffffhhl
+
∥
ψ
∂
c
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
∥
L
∞
(
0
,
T
;
H
2
)
2
ffffffffffffhhl
+
∥
∇
·
u
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
t
∥
L
∞
(
0
,
T
;
H
2
)
2
ffffffffffffhhl
+
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
)
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
0
,
T
;
(
L
2
(
Ω
)
)
2
)
2
+
∥
(
∇
·
u
)
t
∥
L
2
(
0
,
T
;
L
2
)
2
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
ffhhhhffhh
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
)
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
h
2
(
∥
u
∥
L
∞
(
0
,
T
;
(
H
1
(
Ω
)
)
2
)
2
+
∥
p
∥
L
∞
(
0
,
T
;
H
1
)
2
)
]
.Proof.
Writec
n
-
c
h
n
=
(
c
n
-
R
h
c
n
)
-
(
c
h
n
-
R
h
c
n
)
=
ρ
n
-
θ
n. Since the estimates of ρ
n are known, it is enough to estimate θ
n.
Multiplying (3) by Π
h
*
χ
h and subtracting the resulting equation from (39) at t
=
t
n
+
1, we obtain
(89)
(
ϕ
c
h
n
+
1
-
c
^
h
n
Δ
t
,
Π
h
*
χ
h
)
+
a
h
(
u
h
n
;
c
h
n
+
1
,
χ
h
)
-
a
h
(
u
n
+
1
;
c
n
+
1
,
χ
h
)
+
(
q
n
+
1
c
h
n
+
1
,
Π
h
*
χ
h
)
-
(
q
n
+
1
c
n
+
1
,
Π
h
*
χ
)
=
(
u
n
+
1
·
∇
c
n
+
1
+
ϕ
∂
c
n
+
1
∂
t
,
Π
h
*
χ
h
)
∀
χ
h
∈
M
h
.
Choose χ
h
=
θ
n
+
1 in (89) and use the definition of R
h to obtain
(90)
(
ϕ
θ
n
+
1
-
θ
n
Δ
t
,
Π
h
*
θ
n
+
1
)
+
a
h
(
u
h
n
;
θ
n
+
1
,
θ
n
+
1
)
=
(
θ
n
+
1
q
n
+
1
,
Π
h
*
θ
n
+
1
)
+
(
ρ
n
+
1
,
θ
n
+
1
)
+
(
ρ
n
+
1
q
n
+
1
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
+
[
a
h
(
u
n
+
1
;
R
h
c
n
+
1
,
θ
n
+
1
)
-
a
h
(
u
h
n
+
1
;
R
h
c
n
+
1
,
θ
n
+
1
)
]
+
(
u
n
+
1
·
∇
c
n
+
1
+
ϕ
∂
c
n
+
1
∂
t
-
ϕ
(
c
n
+
1
-
c
ˇ
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
+
(
ϕ
(
ρ
n
+
1
-
ρ
ˇ
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
-
(
ϕ
(
θ
n
-
θ
^
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
+
(
ϕ
(
R
^
h
c
n
-
R
ˇ
h
c
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
=
T
1
+
T
2
+
T
3
+
T
4
+
T
5
+
T
6
+
T
7
+
T
8
.
To estimate T
1, T
2, and T
3, we use the Cauchy-Schwarz inequality, boundedness of q, and (54) to obtain
(91)
|
T
1
|
=
|
(
θ
n
+
1
q
n
+
1
,
Π
h
*
θ
n
+
1
)
|
≤
C
∥
θ
n
+
1
∥
2
,
|
T
2
|
=
|
(
ρ
n
+
1
,
θ
n
+
1
)
|
≤
C
∥
ρ
n
+
1
∥
∥
θ
n
+
1
∥
,
|
T
3
|
=
|
(
ρ
n
+
1
q
n
+
1
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
|
≤
C
∥
ρ
n
+
1
∥
∥
θ
n
+
1
∥
.
To bound T
4, we use Lemma 10 to obtain
(92)
|
T
4
|
=
|
a
h
(
u
n
+
1
;
R
h
c
n
+
1
,
θ
n
+
1
)
h
h
-
a
h
(
u
h
n
+
1
;
R
h
c
n
+
1
,
θ
n
+
1
)
|
≤
C
(
∥
u
n
+
1
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
h
h
h
h
∥
u
n
+
1
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
+
h
∥
∇
·
(
u
n
+
1
-
u
h
n
+
1
)
∥
)
|
θ
n
+
1
|
1
≤
C
(
∥
u
n
+
1
-
u
n
∥
(
L
2
(
Ω
)
)
2
h
h
h
h
+
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
h
h
h
h
+
h
∥
∇
·
(
u
n
+
1
-
u
n
)
∥
h
h
h
h
∥
u
n
+
1
-
u
n
∥
(
L
2
(
Ω
)
)
2
+
h
∥
∇
·
(
u
n
-
u
h
n
)
∥
)
|
θ
n
+
1
|
1
.
Since
(93)
|
u
n
+
1
-
u
n
|
2
=
|
∫
t
n
t
n
+
1
u
t
d
s
‍
|
2
≤
Δ
t
∫
t
n
t
n
+
1
‍
|
u
t
|
2
d
s
,
hence,
(94)
∥
u
n
+
1
-
u
n
∥
(
L
2
(
Ω
)
)
2
2
≤
Δ
t
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
(
L
2
(
Ω
)
)
2
)
2
,
and, similarly,
(95)
∥
∇
·
(
u
n
+
1
-
u
n
)
∥
L
2
(
Ω
)
2
≤
Δ
t
∥
∇
·
u

t

∥
L
2
(
t
n
,
t
n
+
1
;
L
2
(
Ω
)
)
2
.
Hence,
(96)
|
T
4
|
≤
C
[
(
Δ
t
)
1
/
2
(
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
,
(
L
2
(
Ω
)
)
2
)
hhhhhhhhhh
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
,
(
L
2
(
Ω
)
)
2
)
+
h
∥
(
∇
·
u
)
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
(
Ω
)
)
)
gggg
+
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
gggg
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
,
(
L
2
(
Ω
)
)
2
)
+
h
∥
∇
·
(
u
n
-
u
h
n
)
∥
]
|
θ
n
+
1
|
.
Using the Cauchy-Schwarz inequality and (54), we obtain
(97)
|
T
5
|
=
|
(
u
n
+
1
·
∇
c
n
+
1
+
ϕ
∂
c
n
+
1
∂
t
-
ϕ
(
c
n
+
1
-
c
^
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
|
≤
C
∥
u
n
+
1
·
∇
c
n
+
1
+
ϕ
∂
c
n
+
1
∂
t
-
ϕ
(
c
n
+
1
-
c
^
n
)
Δ
t
∥
∥
θ
n
+
1
∥
.
Let σ
(
x
)
=
[
ϕ
(
x
)
2
+
u
n
+
1
(
x
)
2
]
1
/
2, so that
(98)
ϕ
∂
c
n
+
1
∂
t
+
u
n
+
1
·
∇
c
n
+
1
=
σ
∂
c
n
+
1
∂
τ
,
where τ approximates the characteristic unit vector s. Now using the same arguments given in [12], we have the following bound:
(99)
∥
σ
∂
c
n
+
1
∂
τ
-
ϕ
c
n
+
1
-
c
ˇ
n
Δ
t
∥
2
≤
C
Δ
t
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2
.
Hence, T
5 is bounded as follows:
(100)
|
T
5
|
≤
C
(
Δ
t
)
1
/
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
∥
θ
n
+
1
∥
.
To bound T
6, we proceed as follows:
(101)
|
T
6
|
=
|
(
ϕ
(
ρ
n
+
1
-
ρ
ˇ
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
|
≤
|
(
ϕ
(
ρ
n
+
1
-
ρ
ˇ
n
)
Δ
t
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
|
+
|
(
ϕ
(
ρ
n
+
1
-
ρ
ˇ
n
)
Δ
t
,
θ
n
+
1
)
|
=
I
1
+
I
2
.
Now I
1 can be written as
(102)
I
1
≤
|
(
ϕ
(
ρ
n
+
1
-
ρ
n
)
Δ
t
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
|
+
|
(
ϕ
(
ρ
n
-
ρ
ˇ
n
)
Δ
t
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
|
.
Following the proof lines of Theorem 4.1 in [12], it can be shown that, for f
∈
H
1
(
Ω
), there exists a positive constant C independent of h and Δ
t such that
(103)
∥
f
-
f
ˇ
∥
≤
C
Δ
t
∥
∇
f
n
∥
.
A use of (103) yields
(104)
∥
(
ρ
n
-
ρ
ˇ
n
)
Δ
t
∥
≤
C
∥
∇
ρ
n
∥
.
It is easy to show that
(105)
∥
(
ρ
n
+
1
-
ρ
n
)
Δ
t
∥
≤
C
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
.
Using (52), (104), (105), and the Cauchy-Schwarz inequality, we obtain
(106)
I
1
≤
C
[
∥
θ
n
+
1
-
Π
h
*
θ
n
+
1
∥
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
h
h
h
h
h
×
(
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
+
∥
∇
ρ
n
∥
)
]
≤
C
h
[
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
g
g
g
g
l
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
+
∥
∇
ρ
n
∥
]
∥
∇
θ
n
+
1
∥
.
I
2 can be bounded as follows:
(107)
I
2
≤
|
(
ϕ
(
ρ
n
+
1
-
ρ
n
)
Δ
t
,
θ
n
+
1
)
|
+
|
(
ϕ
(
ρ
n
-
ρ
ˇ
n
)
Δ
t
,
θ
n
+
1
)
|
=
J
1
+
J
2
.
A use of Cauchy-Schwarz inequality yields
(108)
J
1
≤
C
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
∥
θ
n
+
1
∥
.
In order to bound J
2, we use the following result.
Ifη
∈
L
2
(
Ω
) and η
ˇ
(
x
)
=
η
(
x
ˇ
) with x
ˇ
=
x
-
r
(
x
)
Δ
t, for a nonzero function r
(
x
) such that r and ∇
·
r are bounded, then
(109)
∥
η
-
η
ˇ
∥
-
1
≤
C
∥
η
∥
Δ
t
,
where C is a positive constant independent of h and Δ
t; for a proof, we refer to [9, pp. 875]. Using (109), we obtain
(110)
J
2
≤
C
∥
(
ρ
n
-
ρ
^
n
)
Δ
t
∥
-
1
∥
θ
n
+
1
∥
1
≤
C
∥
ρ
n
∥
∥
θ
n
+
1
∥
1
.
This implies that
(111)
I
2
≤
C
(
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
+
∥
ρ
n
∥
)
∥
θ
n
+
1
∥
.
Now, using (106), (111), and (101), we obtain the following bound for T
6:
(112)
|
T
6
|
≤
C
(
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
h
h
h
h
h
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
+
h
∥
∇
ρ
n
∥
+
∥
ρ
n
∥
)
∥
θ
n
+
1
∥
.
Using the Cauchy-Schwarz inequality, we obtain
(113)
|
T
7
|
≤
|
(
(
θ
n
-
θ
^
n
)
Δ
t
,
Π
h
*
θ
n
+
1
-
θ
n
+
1
)
|
+
|
(
(
θ
n
-
θ
^
n
)
Δ
t
,
θ
n
+
1
)
|
≤
∥
(
θ
n
-
θ
^
n
)
Δ
t
∥
∥
Π
h
*
θ
n
+
1
-
θ
n
+
1
∥
+
∥
(
θ
n
-
θ
^
n
)
Δ
t
∥
-
1
∥
θ
n
+
1
∥
1
.
We use (103) and (109) to bound ∥
(
θ
n
-
θ
^
n
)
/
Δ
t
∥ and ∥
(
θ
n
-
θ
^
n
)
/
Δ
t
∥
-
1, respectively. For this, we need that u
h
n and its first derivative are bounded.
First let us make an induction hypothesis. Assume that there is a constant sayK
*
≥
2
K with ∥
u
~
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
K such that
(114)
∥
u
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
K
*
,
where u
~
h
n is the projection of u
h
n at t
=
t
n defined in (41). To bound ∥
∇
·
u
h
n
∥
∞, we use inverse inequality and (114):
(115)
∥
∇
·
u
h
n
∥
∞
Δ
t
≤
C
h
-
1
Δ
t
∥
u
h
n
∥
∞
≤
C
,
where we have used the assumption that Δ
t
=
O
(
h
).
Using (103), we have
(116)
∥
θ
n
-
θ
^
n
Δ
t
∥
≤
C
(
K
*
)
∥
∇
θ
n
∥
.
Using (109), we have
(117)
∥
θ
n
-
θ
^
n
Δ
t
∥
-
1
≤
C
(
K
*
)
∥
θ
n
∥
.
Now, using (52), (116), inverse inequality (48), and (117), we obtain the following bound for T
7:
(118)
|
T
7
|
≤
C
(
K
*
)
∥
θ
n
∥
∥
θ
n
+
1
∥
1
.
To bound T
8, we use maximum norm estimate of ∇
R
h
c (see (79)):
(119)
∥
R
^
h
c
n
-
R
ˇ
h
c
n
Δ
t
∥
≤
∥
R
h
c
∥
1
,
∞
∥
u
n
-
u
h
n
∥
≤
C
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
,
And, hence, using (54), T
8 is bounded as follows:
(120)
T
8
≤
C
∥
θ
n
+
1
∥
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
.
Substitute the estimates of T
1
,
T
2
⋯
T
8 in (90) and use nonsingular property of ϕ and a kick-back argument with the Young’s inequality to obtain
(121)
1
Δ
t
[
(
θ
n
+
1
,
Π
h
*
θ
n
+
1
)
-
(
θ
n
,
Π
h
*
θ
n
+
1
)
]
≤
C
(
K
*
)
[
Δ
t
-
1
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
hhhhhhhhh
+
∥
θ
n
+
1
∥
2
+
∥
θ
n
+
1
∥
1
2
hhhhhhhhh
+
∥
θ
n
∥
2
+
∥
ρ
n
+
1
∥
2
hhhhhhhhh
+
Δ
t
(
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
hhhhhggggghhhhhl
+

∥
(
∇
·
u
)
t
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2

hhhhhggggghhhhlh
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
)
hhhhhhhhh
+
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
2
hhhhhhhhh
+
h
2
∥
∇
·
(
u
n
-
u
h
n
)
∥
2
hhhhhhhhh
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
+
∥
ρ
n
∥
2
+
h
2
∥
∇
ρ
n
∥
2
]
.
Using (43) and (44), we obtain
(122)
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
≤
C
[
∥
c
n
-
c
h
n
∥
+
h
(
∥
u
n
∥
(
H
1
(
Ω
)
)
2
+
∥
p
n
∥
1
2
)
]
≤
C
[
∥
θ
n
∥
+
∥
ρ
n
∥
+
h
(
∥
u
n
∥
(
H
1
(
Ω
)
)
2
+
∥
p
n
∥
1
2
)
]
,
∥
∇
·
(
u
n
-
u
h
n
)
∥
≤
C
h
∥
∇
·
u
n
∥
1
.
Substitute (122) in (121) to obtain
(123)

∥
|
θ
n
+
1
|
∥
2
-
∥
|
θ
n
|
∥
2

≤
C
(
K
*
)
[
Δ
t
(
∥
θ
n
+
1
∥
2
+
∥
θ
n
∥
2
+
∥
ρ
n
∥
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
ffffffffffffffffffffhfh
∥
θ
n
+
1
∥
2
+
∥
ρ
n
+
1
∥
2
+
h
2
∥
∇
ρ
n
∥
2
)
fffffffffffffffhff
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
fffffffffffffffffffhhhhhhhl
+
∥
∇
·
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2
fffffffffffffffffffhhhhhhhl
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
)
fffffffffffffffhff
+
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
ffffffffffffffhfff
+
h
2
Δ
t
(
h
2
∥
∇
·
u
n
∥
1
2
+
∥
u
n
∥
(
H
1
(
Ω
)
)
2
2
fffffffffffffffffhffhhhhl
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
∥
u
n
∥
(
H
1
(
Ω
)
)
2
2
+
∥
p
n
∥
1
2
)
]
.
Taking summation over n
=
0
⋯
m
-
1, we obtain
(124)

∥
|
θ
m
|
∥
2
-
∥
|
θ
0
|
∥
2

≤
C
(
K
*
)
[
∑
n
=
0
m
-
1
‍
{
Δ
t
(
∥
θ
n
+
1
∥
2
+
∥
θ
n
∥
2
+
∥
ρ
n
∥
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
hhhhhhhhhhhhhhhhhl
+
∥
ρ
n
+
1
∥
2
+
h
2
∥
∇
ρ
n
∥
2
)
hhhhhhhhhhhhh
.
hh
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
hhhhhhhhhhhhhhhhhhhhhh
.
+
∥
∇
·
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2
hhhhhhhhhhhhhhhhhhhhhh
.
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
)
hhhhhhhhhhhhhhh
.
+
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2
hhhhhhhhhhhhhhh
.
+
h
2
Δ
t
(
h
2
∥
∇
·
u
n
∥
1
2
+
∥
u
n
∥
(
H
1
(
Ω
)
)
2
2
hhhhhhhhhhhhhhhhhhhhl
∑
n
=
0
m
-
1
‍
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
∥
∇
·
u
n
∥
1
2
∥
u
n
∥
(
H
1
(
Ω
)
)
2
2
+
∥
p
n
∥
1
2
)
}
]
.
Now, using discrete Gronwall’s, equivalence of the norms (53), and the estimates of ρ, we obtain
(125)
∥
θ
m
∥
2
≤
C
(
K
*
)
×
[
∥
θ
0
∥
2
+
h
4
(
∥
(
c
-
c
~
)
q
∥
L
∞
(
0
,
T
;
H
1
)
2
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
jjjjjjjjjjjjjjjjjjjjjjjhj
+
∥
ψ
∂
c
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
∥
L
∞
(
0
,
T
;
H
2
)
2
jjjjjjjjjjjjjjjjjjjjjjhjj
+
∥
∇
·
u
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
t
∥
L
∞
(
0
,
T
;
H
2
)
2
jjjjjjjjjjjjjjjjjjjhjjjjj
+
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
)
hh
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
0
,
T
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
jjjjjjjjjjhjhjjjjjjk
+
∥
∇
·
u
t
∥
L
2
(
0
,
T
;
L
2
)
2
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
)
jhjjhh
+
h
2
(
∥
u
∥
L
∞
(
0
,
T
;
(
H
1
(
Ω
)
)
2
)
2
jhjjhhhh
.
hh
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
∥
u
∥
L
∞
(
0
,
T
;
(
H
1
(
Ω
)
)
2
)
2
+
∥
p
∥
L
∞
(
0
,
T
;
H
1
)
2
)
]
.
Now, it remains to show the induction hypothesis (114). Using (23) and (24), we have
(126)
∥
u
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
∥
u
h
n
-
u
~
h
n
∥
(
L
∞
(
Ω
)
)
2
+
∥
u
~
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
C
(
log
⁡
1
h
)
1
/
2
∥
u
h
n
-
u
~
h
n
∥
H
(
div
⁡
;
Ω
)
+
K
.
Using ∥
∇
·
(
u
h
n
-
u
~
h
n
)
∥
=
0, we have
(127)
∥
u
h
n
-
u
~
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
C
(
log
⁡
1
h
)
1
/
2
∥
u
h
n
-
u
~
h
n
∥
(
L
2
(
Ω
)
)
2
.
Now using (42) and (125), we obtain for small h
(128)
∥
u
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
C
(
K
*
)
log
⁡
1
h
(
h
+
Δ
t
)
+
K
≤
2
K
.
Here, we have used Δ
t
=
O
(
h
) and h
log
⁡
(
1
/
h
)
→
0 as h
→
0 and this proves our induction hypothesis (114).
Now combine the estimates ofρ and θ and use triangle inequality to complete the rest of the proof.Using (43) and Theorem 11, we obtain the following error estimates for velocity as well as pressure.Theorem 12.
Assume that the triangulation𝒯
h is quasiuniform. Let (
u
,
p
) and (
u
h
,
p
h
) be, respectively, the solutions of (1)-(2) and (30) and let c
h
(
0
)
=
c
0
,
h
=
R
h
c
(
0
). Further assume that Δ
t
=
O
(
h
). Then for sufficiently small h there exists a positive constant C
(
T
) independent of h but dependent on the bounds of κ
-
1 and μ such that
(129)
max
⁡
0
≤
n
≤
N
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
2
≤
C
[
h
4
(
∥
(
c
-
c
~
)
q
∥
L
∞
(
0
,
T
;
H
1
)
2
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
hhhhhhhhhh
.
+
∥
ψ
∂
c
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
∥
L
∞
(
0
,
T
;
H
2
)
2
hhhhhhhhhh
.
+
∥
∇
·
u
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
t
∥
L
∞
(
0
,
T
;
H
2
)
2
hhhhhhhhhh
.
+
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
)
hhhhhhhh
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
0
,
T
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
hhhhhhjjjjjjjjjhhhhh
+
∥
∇
·
u
t
∥
L
2
(
0
,
T
;
L
2
)
2
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
)
hhhhhhhh
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
h
2
(
∥
u
∥
L
∞
(
0
,
T
;
(
H
1
(
Ω
)
)
2
)
2
+
∥
p
∥
L
∞
(
0
,
T
;
H
1
)
2
)
]
.


## 3.1. Error Estimates for Velocity

For a givenc, the auxiliary functions (
u
~
h
,
p
~
h
)
:
[
0
,
T
]
→
U
h
×
W
h are defined as follows:
(41)
(
κ
-
1
μ
(
c
)
u
~
h
,
v
h
)
-
(
∇
·
v
h
,
p
~
h
)
=
0
∀
v
h
∈
U
h
,
(
∇
·
u
~
h
,
w
h
)
=
(
q
,
w
h
)
∀
w
h
∈
W
h
.
For a proof of existence and uniqueness of the discrete solution of (41), we refer to [20, pp. 52]. For u
h and u
~
h, the following error estimates can be obtained (see [14]):
(42)
∥
u
h
-
u
~
h
∥
(
L
2
(
Ω
)
)
2
≤
C
(
h
∥
u
~
h
∥
(
L
2
(
Ω
)
)
2
+
∥
c
-
c
h
∥
∥
u
~
h
∥
(
L
∞
(
Ω
)
)
2
)
.
Now, since concentration depends on the velocity and vice versa, to derive the error estimates for the concentration, we also need error estimates for the velocity. For elliptic problems, the authors in [14] have derived error estimates for mixed covolume method by using Raviart-Thomas projection and L
2 projection. In the similar way, for a given c, the following theorem can be shown but the proof is long so we omit it here.Theorem 2.
Assume that the triangulation𝒯
h is quasiuniform. Let (
u
,
p
) and (
u
h
,
p
h
), respectively, be the solutions of (1)-(2) and (30). Then, there exists a positive constant C independent of h but dependent on the bounds of κ
-
1 and μ such that
(43)
∥
u
-
u
h
∥
(
L
2
(
Ω
)
)
2
+
∥
p
-
p
h
∥
≤
C
[
∥
c
-
c
h
∥
+
h
(
∥
u
∥
(
H
1
(
Ω
)
)
2
+
∥
p
∥
1
)
]
,
(44)
∥
∇
·
(
u
-
u
h
)
∥
≤
C
h
∥
∇
·
u
∥
1
.


## 3.2. Error Estimates for Concentration

We writec
-
c
h
=
(
c
-
R
h
c
)
+
(
R
h
c
-
c
h
), where R
h
:
H
1
(
Ω
)
→
M
h is the projection of c defined by
(45)
A
(
u
;
c
-
R
h
c
,
χ
h
)
=
0
∀
χ
h
∈
M
h
,
where
(46)
A
(
u
;
ϕ
,
χ
h
)
=
a
h
(
u
;
ϕ
,
χ
h
)
+
(
q
ϕ
,
χ
h
)
+
(
λ
ϕ
,
χ
h
)
∀
χ
h
∈
M
h
.
The function λ will be chosen such that the coercivity of A
(
u
;
·
,
·
) is assured.We use frequently the following trace inequality [21, pp. 417]: for w
∈
H
1
(
T
),
(47)
∥
w
∥
∂
T
2
≤
C
(
h
T
-
1
∥
w
∥
T
2
+
h
T
|
w
|
1
,
T
2
)
,
where ∥
w
∥
∂
T
2
=
∫
∂
T
‍
w
2
d
s. Further, we need the following inverse inequalities (see [22, pp. 141]): ∀
χ
h
∈
M
h
(48)
∥
χ
h
∥
j
,
∞
≤
C
h
-
1
∥
χ
h
∥
1
,
j
=
0,1
,
g
g
g
g
g
∥
χ
h
∥
1
≤
C
h
-
1
∥
χ
h
∥
.
Using the properties of Π
h
* operator, it is easy to see that, for T
∈
𝒯
h and ϕ
h
∈
M
h the following holds true:
(49)
∫
T
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
‍
=
0
,
∫
∂
T
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
s
=
0
.
Now, using (49), we have the following lemma [17, pp. 317].Lemma 3.
There exists a positive constantC independent of h such that
(50)
j
j
j
j
j
j
j
|
ϵ
h
(
ϕ
,
χ
)
|
≤
C
h
i
+
j
|
ϕ
|
W
p
i
|
χ
|
W
q
j
∀
χ
∈
M
h
w
i
t
h
i
,
j
=
0,1
,
1
p
+
1
q
=
1
,
where
(51)
ϵ
h
(
ϕ
,
χ
h
)
=
(
ϕ
,
χ
h
)
-
(
ϕ
,
Π
h
*
χ
h
)
∀
χ
h
∈
M
h
.Also note that, by the usual interpolation theory, the operatorΠ
h
* has the following approximation property [23, pp. 466]:
(52)
jj
∥
χ
-
Π
h
*
χ
∥
0
,
k
≤
C
h
β
|
χ
|
s
,
k
,
0
≤
β
≤
s
≤
1
,
1
≤
k
≤
∞
.
We also recall two well-known results (Lemmas 4 and 5), which will be used in the the proof of Lemmas 7–9.Lemma 4 (see [24, pp. 240]).
The operatorΠ
h
* has the following properties.(i)
With∥
|
ϕ
h
|
∥
=
(
ϕ
h
,
Π
h
*
ϕ
h
)
1
/
2, the norms ∥
|
·
|
∥ and ∥
·
∥ are equivalent on U
h; that is, there exist positive constants C
7 and C
8, independent of h, such that
(53)
C
7
∥
ϕ
h
∥
≤
∥
|
ϕ
h
|
∥
≤
C
8
∥
ϕ
h
∥
∀
ϕ
h
∈
M
h
.
(ii)
Π
h
* is stable with respect to the L
2 norm; that is, there exists a positive constant C independent of h such that
(54)
∥
Π
h
*
χ
h
∥
≤
C
∥
χ
h
∥
∀
χ
h
∈
M
h
.Lemma 5 (see [25, pp. 1871]).
Assume thatχ
h
,
ϕ
h
∈
M
h. Then one has
(55)
a
h
(
u
;
χ
h
,
ϕ
h
)
=
a
(
u
;
χ
h
,
ϕ
h
)
+
∑
T
∈
𝒯
h
∫
∂
T
(
D
(
u
)
∇
χ
h
·
n
)
(
Π
h
*
ϕ
h
-
ϕ
h
)
d
s
‍
‍
+
∑
T
∈
𝒯
h
∫
T
∇
·
(
D
(
u
)
∇
χ
h
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
.
Moreover, for χ
h
∈
M
h,
(56)
a
h
(
u
;
χ
h
,
χ
h
)
≥
a
(
u
;
χ
h
,
χ
h
)
-
C
h
|
χ
h
|
1
2
.Remark 6.
Note that (55) also holds true for χ
∈
H
1
(
Ω
).SinceD is uniformly positive definite, we obtain following from (56):
(57)
a
h
(
u
;
χ
h
,
χ
h
)
≥
(
α
-
C
h
)
|
χ
h
|
1
2
.
Choose h
0
>
0 such that, for 0
<
h
<
h
0, (
α
-
C
h
)
=
α
0
>
0 and hence
(58)
a
h
(
u
h
;
χ
h
,
χ
h
)
≥
α
0
|
χ
h
|
1
2
∀
χ
h
∈
M
h
.
Now, we derive the error bound in H
1 and L
2 norms for c
-
R
h
c. Let I
h be the continuous interpolant onto M
h satisfying the following approximation properties. For ϕ
∈
H
k
+
1
(
Ω
) with k
=
0,1, we have [22]
(59)
∥
ϕ
-
I
h
ϕ
∥
j
≤
h
k
+
1
-
j
∥
ϕ
∥
k
+
1
j
=
0,1
.
Moreover, if ϕ
∈
W
2
,
∞
(
Ω
), then
(60)
∥
ϕ
-
I
h
ϕ
∥
1
,
∞
≤
C
h
log
⁡
(
1
h
)
∥
ϕ
∥
2
,
∞
.Lemma 7.
There exists a positive constantC independent of h such that
(61)
∥
c
-
R
h
c
∥
1
≤
C
h
∥
c
∥
2
,
provided c
∈
H
2
(
Ω
), for a.e. t
∈
(
0
,
T
].Proof.
The coercivity and boundedness ofA
(
u
;
·
,
·
) with (45) yield
(62)
∥
I
h
c
-
R
h
c
∥
1
2
≤
C
A
(
u
;
I
h
c
-
R
h
c
,
I
h
c
-
R
h
c
)
≤
C
A
(
u
;
I
h
c
-
c
,
I
h
c
-
R
h
c
)
≤
C
∥
c
-
I
h
c
∥
1
∥
I
h
c
-
R
h
c
∥
1
,
and hence
(63)
∥
I
h
c
-
R
h
c
∥
1
≤
C
∥
c
-
I
h
c
∥
1
,
where C depends on the bound of D
(
u
) given in (11). Combine the estimates (63) and (59) and use the triangle inequality to complete the proof.
For deriving theL
2 error bounds for c
-
R
h
c, we need the following lemma.
Lemma 8.Let
ϵ
a
(
u
;
ϕ
,
χ
h
)
=
a
(
u
;
ϕ
,
χ
h
)
-
a
h
(
u
;
ϕ
,
χ
h
)
. There exists a positive constant
C
such that, for
ϕ
∈
H
1
(
Ω
)
and
χ
h
∈
M
h,
(64)
|
ϵ
a
(
u
;
c
-
R
h
c
,
ϕ
h
)
|
≤
C
h
2
(
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
+
∥
c
∥
2
)
|
ϕ
h
|
1
∀
ϕ
h
∈
M
h
.
Using (55) (see also Remark 6), we find that
(65)
|
ϵ
a
(
u
;
c
-
R
h
c
,
ϕ
h
)
|
≤
|
∑
T
∈
𝒯
h
∫
T
∇
·
(
D
(
u
)
∇
(
c
-
R
h
c
)
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
‍
|
+
|
∑
T
∈
𝒯
h
∫
∂
T
(
D
(
u
)
∇
(
c
-
R
h
c
)
·
n
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
s
‍
|
=
J
1
+
J
2
.
To bound J
1, first we use the fact that R
h
c is linear on each triangle T to obtain
(66)
J
1
=
|
∑
T
∈
𝒯
h
∫
T
∇
·
(
D
(
u
)
∇
(
c
-
R
h
c
)
)
g
g
g
g
g
g
g
g
∑
T
∈
𝒯
h
×
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
|
=
|
∑
T
∈
𝒯
h
∫
T
(
∇
·
(
D
(
u
)
∇
c
)
-
(
∇
·
D
(
u
)
)
·
∇
R
h
c
)
g
g
g
g
g
g
g
g
∑
T
∈
𝒯
h
×
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
|
.
Now use (3), (49), and (50) to obtain
(67)
J
1
≤
|
∑
T
∈
𝒯
h
∫
T
(
-
g
+
ψ
∂
c
∂
s
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
|
+
|
∑
T
∈
𝒯
h
∫
T
[
(
∇
·
D
(
u
)
-
(
∇
·
D
(
u
)
)
T
)
·
∇
R
h
c
]
g
g
g
g
g
g
g
g
∑
T
∈
𝒯
h
×
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
x
|
≤
C
h
2
(
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
+
∥
c
∥
2
)
|
ϕ
h
|
1
,
where (
∇
·
D
(
u
)
)
T denotes the average value of ∇
·
D
(
u
) on triangle T.
Based on the analysis in [25, pp. 1873], we estimate J
2 as follows. Note that an appeal to the continuity of ∇
c
·
n with (49) yields
(68)
J
2
=
|
∑
T
∈
𝒯
h
∫
∂
T
(
(
D
-
D
¯
T
)
∇
(
c
-
R
h
c
)
·
n
)
‍
g
g
g
g
g
g
g
g
∑
T
∈
𝒯
h
×
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
s
|
,
where D
=
D
(
u
) and D
¯
T is a function such that, for any edge of a triangle T
∈
𝒯
h,
(69)
D
¯
T
(
x
)
=
D
(
x
c
)
,
x
∈
E
,
and x
c is the midpoint of E. Since |
D
(
x
)
-
D
¯
T
|
≤
C
h
T
∥
D
∥
1
,
∞, we use trace inequalities (47) and (61) to arrive at
(70)
J
2
≤
C
h
|
∑
T
∈
𝒯
h
∫
∂
T
(
∇
(
c
-
R
h
c
)
·
n
)
(
ϕ
h
-
Π
h
*
ϕ
h
)
d
s
‍
|
≤
C
h
(
∑
T
∈
𝒯
h
∫
∂
T
|
(
∇
(
c
-
R
h
c
)
·
n
)
|
2
‍
)
1
/
2
×
(
∑
T
∈
𝒯
h
∫
∂
T
|
ϕ
h
-
Π
h
*
ϕ
h
|
2
‍
)
1
/
2
≤
C
h
(
h
T
-
1
/
2
∥
c
-
R
h
c
∥
1
+
h
T
1
/
2
∥
c
∥
2
)
×
(
h
T
-
1
/
2
∥
ϕ
h
-
Π
h
*
ϕ
h
∥
+
h
T
1
/
2
|
ϕ
h
|
1
)
≤
C
h
2
∥
c
∥
2
|
ϕ
h
|
1
.
Substitute the estimates of J
1 and J
2 in (65) to complete  the rest of the proof.Lemma 9.
There exists a positive constantC independent of h such that
(71)
∥
c
-
R
h
c
∥
≤
C
h
2
(
∥
c
∥
2
+
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
)
,
provided c
∈
H
2
(
Ω
)
,
u
·
∇
c
∈
H
1
(
Ω
), and ∂
c
/
∂
t
∈
H
1
(
Ω
) for t
∈
(
0
,
T
] a.e.Proof.
To obtain optimalL
2 error estimates for c
-
R
h
c, we now appeal to Aubin-Nitsche duality argument. Let ϕ
∈
H
2
(
Ω
) be a solution of the following adjoint problem:
(72)
-
∇
·
(
D
(
u
)
∇
ϕ
+
u
ϕ
)
+
λ
ϕ
=
c
-
R
h
c
in
Ω
,
(
D
(
u
)
∇
ϕ
+
u
ϕ
)
·
n
=
0
on
⁡
∂
Ω
,
which satisfies the elliptic regularity condition
(73)
∥
ϕ
∥
2
≤
C
∥
c
-
R
h
c
∥
.
Multiply the above equation by c
-
R
h
c and integrate over Ω. An integration by parts and a use of (45) yield
(74)
∥
c
-
R
h
c
∥
2
=
a
(
u
;
ϕ
,
c
-
R
h
c
)
-
(
u
·
∇
ϕ
,
c
-
R
h
c
)
-
(
∇
·
u
ϕ
,
c
-
R
h
c
)
+
λ
(
ϕ
,
c
-
R
h
c
)
=
[
a
(
u
;
c
-
R
h
c
,
ϕ
-
ϕ
h
)
.
+
(
u
·
∇
(
c
-
R
h
c
)
,
ϕ
-
ϕ
h
)
.
+
λ
(
c
-
R
h
c
,
ϕ
-
ϕ
h
)
]
.
+
ϵ
a
(
u
;
c
-
R
h
c
,
ϕ
h
)
∀
ϕ
h
∈
M
h
=
I
1
+
I
2
.
For I
1, use (61) to find that
(75)
|
I
1
|
=
|
a
(
u
;
c
-
R
h
c
,
ϕ
-
ϕ
h
)
ggl
+
(
u
·
∇
(
c
-
R
h
c
)
,
ϕ
-
ϕ
h
)
ggl
+
λ
(
c
-
R
h
c
,
ϕ
-
ϕ
h
)
|
≤
C
∥
c
-
R
h
c
∥
1
∥
ϕ
-
ϕ
h
∥
1
≤
C
h
∥
c
∥
2
∥
ϕ
-
ϕ
h
∥
1
.
The bound for I
2 follows from Lemma 8 and hence
(76)
|
I
2
|
≤
C
h
2
(
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
+
∥
c
∥
2
)
|
ϕ
h
|
1
.
Substitute (75) and (76) in (74) to find that
(77)
∥
c
-
R
h
c
∥
2
≤
C
[
h
∥
c
∥
2
∥
ϕ
-
ϕ
h
∥
1
|
ψ
∂
c
∂
s
|
1
h
h
h
h
+
h
2
(
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
+
∥
c
∥
2
)
|
ϕ
h
|
1
]
.
Now choose ϕ
h
=
I
h
ϕ in (77). Then use elliptic regularity condition (73) with (59) to obtain
(78)
∥
c
-
R
h
c
∥
≤
C
h
2

(
∥
c
∥
2
+
|
g
|
1
+
|
ψ
∂
c
∂
s
|
1
)

,
and this completes the proof.A use of inverse inequalities (48), (59), (60), and (63) yields
(79)
∥
R
h
c
∥
1
,
∞
≤
C
∥
c
∥
2
,
∞
.Lemma 10.
There exists a positive constantC such that ∀
θ
∈
M
h,
(80)
|
a
h
(
u
;
R
h
c
,
θ
)
-
a
h
(
u
h
;
R
h
c
,
θ
)
|
≤
C
(
∥
u
-
u
h
∥
(
L
2
(
Ω
)
)
2
+
h
∥
∇
·
(
u
-
u
h
)
∥
)
|
θ
|
1
.Proof.
Note that(81)
|
a
h
(
u
;
R
h
c
,
θ
)
-
a
h
(
u
h
;
R
h
c
,
θ
)
|
=
|
∑
i
=
1
N
h
∫
∂
V
i
*
(
D
(
u
)
-
D
(
u
h
)
)
∇
R
h
c
·
n
i
Π
h
*
θ
d
s
‍
|
=
|
∑
T
∈
𝒯
h
K
T
‍
|
,
where K
T
=
∑
l
=
1
3
‍
∫
∂
V
l
*
∩
T
‍
(
D
(
u
)
-
D
(
u
h
)
)
∇
R
h
c
·
n
l
θ
l
d
s and θ
l
=
θ
(
P
l
); see Figure 5. For each triangle T, K
T can be written as
(82)
K
T
=
∑
l
=
1
3
∫
M
l
B
(
D
(
u
)
-
D
(
u
h
)
)
∇
R
h
c
·
n
l
‍
‍
hhhhhh
×
(
θ
l
+
1
-
θ
l
)
d
s
(
θ
4
=
θ
1
)
.
Using the Cauchy-Schwarz inequality and (79), we obtain
(83)
K
T
≤
∑
l
=
1
3
|
θ
l
+
1
-
θ
l
|
‍
×
∫
M
l
B
|
(
D
(
u
)
-
D
(
u
h
)
)
∇
R
h
c
·
n
l
|
d
s
‍
≤
C
∑
l
=
1
3
|
θ
l
+
1
-
θ
l
|
∥
D
(
u
)
-
D
(
u
h
)
∥
(
L
2
(
M
l
B
)
)
2
×
2
‍
×
(
meas
⁡
(
M
l
B
)
)
1
/
2
.
It has been proved in [26, pp. 332] that the matrix D is uniformly Lipschitz; that is, there exists a constant C such that for u and v
∈
(
L
2
(
Ω
)
)
2,
(84)
∥
D
(
u
)
-
D
(
v
)
∥
(
L
2
(
Ω
)
)
2
×
2
≤
C
∥
u
-
v
∥
(
L
2
(
Ω
)
)
2
.
A use of the trace inequalities (47) and (84) yields
(85)
K
T
≤
C
∑
l
=
1
3
|
θ
l
+
1
-
θ
l
|
∥
u
-
u
h
∥
(
L
2
(
M
l
B
)
)
2
h
T
1
/
2
‍
≤
C
h
T
1
/
2
∑
l
=
1
3
‍
|
θ
l
+
1
-
θ
l
|
[
h
T
-
1
/
2
∥
u
-
u
h
∥
T
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
l
+
h
T
1
/
2
∥
∇
·
(
u
-
u
h
)
∥
T
]
.
Now, using Taylor series expansion, we find that
(86)
|
θ
l
+
1
-
θ
l
|
≤
h
T
[
|
∂
θ
∂
x
|
+
|
∂
θ
∂
y
|
]
≤
[
(
|
∂
θ
∂
x
|
2
+
|
∂
θ
∂
y
|
2
)
h
T
2
]
1
/
2
=
|
θ
|
1
,
h
,
T
,
l
=
1,2
,
3
.
Let |
ϕ
h
|
1
,
h
=
(
∑
T
∈
𝒯
h
‍
|
ϕ
h
|
1
,
h
,
T
2
)
1
/
2. Since θ is linear on each triangle, seminorms |
·
| and |
·
| are identical.
Substitute (86) in (85) and use the estimates for K
T to arrive at
(87)
|
a
h
(
u
;
R
h
c
,
θ
)
-
a
h
(
u
h
;
R
h
c
,
θ
)
|
≤
C
(
∥
u
-
u
h
∥
(
L
2
(
Ω
)
)
2
+
h
∥
∇
·
(
u
-
u
h
)
∥
)
|
θ
|
1
,
and this completes the rest of the proof.Figure 5
A triangular partition.Now, we prove our main theorem.Theorem 11.
Letc
n and c
h
n be the solutions of (3) and (39) at t
=
t
n, respectively, and let c
h
(
0
)
=
c
0
,
h
=
R
h
c
(
0
). Further assume that Δ
t
=
O
(
h
). Then, for sufficiently small h, there exists a positive constant C
(
T
) independent of h but dependent on the bounds of κ
-
1 and μ such that
(88)
max
⁡
0
≤
n
≤
N
∥
c
n
-
c
h
n
∥
(
L
2
(
Ω
)
)
2
2
≤
C
[
h
4
(
∥
(
c
-
c
~
)
q
∥
L
∞
(
0
,
T
;
H
1
)
2
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
ffffffffffffhhl
+
∥
ψ
∂
c
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
∥
L
∞
(
0
,
T
;
H
2
)
2
ffffffffffffhhl
+
∥
∇
·
u
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
t
∥
L
∞
(
0
,
T
;
H
2
)
2
ffffffffffffhhl
+
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
)
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
0
,
T
;
(
L
2
(
Ω
)
)
2
)
2
+
∥
(
∇
·
u
)
t
∥
L
2
(
0
,
T
;
L
2
)
2
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
ffhhhhffhh
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
)
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
h
2
(
∥
u
∥
L
∞
(
0
,
T
;
(
H
1
(
Ω
)
)
2
)
2
+
∥
p
∥
L
∞
(
0
,
T
;
H
1
)
2
)
]
.Proof.
Writec
n
-
c
h
n
=
(
c
n
-
R
h
c
n
)
-
(
c
h
n
-
R
h
c
n
)
=
ρ
n
-
θ
n. Since the estimates of ρ
n are known, it is enough to estimate θ
n.
Multiplying (3) by Π
h
*
χ
h and subtracting the resulting equation from (39) at t
=
t
n
+
1, we obtain
(89)
(
ϕ
c
h
n
+
1
-
c
^
h
n
Δ
t
,
Π
h
*
χ
h
)
+
a
h
(
u
h
n
;
c
h
n
+
1
,
χ
h
)
-
a
h
(
u
n
+
1
;
c
n
+
1
,
χ
h
)
+
(
q
n
+
1
c
h
n
+
1
,
Π
h
*
χ
h
)
-
(
q
n
+
1
c
n
+
1
,
Π
h
*
χ
)
=
(
u
n
+
1
·
∇
c
n
+
1
+
ϕ
∂
c
n
+
1
∂
t
,
Π
h
*
χ
h
)
∀
χ
h
∈
M
h
.
Choose χ
h
=
θ
n
+
1 in (89) and use the definition of R
h to obtain
(90)
(
ϕ
θ
n
+
1
-
θ
n
Δ
t
,
Π
h
*
θ
n
+
1
)
+
a
h
(
u
h
n
;
θ
n
+
1
,
θ
n
+
1
)
=
(
θ
n
+
1
q
n
+
1
,
Π
h
*
θ
n
+
1
)
+
(
ρ
n
+
1
,
θ
n
+
1
)
+
(
ρ
n
+
1
q
n
+
1
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
+
[
a
h
(
u
n
+
1
;
R
h
c
n
+
1
,
θ
n
+
1
)
-
a
h
(
u
h
n
+
1
;
R
h
c
n
+
1
,
θ
n
+
1
)
]
+
(
u
n
+
1
·
∇
c
n
+
1
+
ϕ
∂
c
n
+
1
∂
t
-
ϕ
(
c
n
+
1
-
c
ˇ
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
+
(
ϕ
(
ρ
n
+
1
-
ρ
ˇ
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
-
(
ϕ
(
θ
n
-
θ
^
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
+
(
ϕ
(
R
^
h
c
n
-
R
ˇ
h
c
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
=
T
1
+
T
2
+
T
3
+
T
4
+
T
5
+
T
6
+
T
7
+
T
8
.
To estimate T
1, T
2, and T
3, we use the Cauchy-Schwarz inequality, boundedness of q, and (54) to obtain
(91)
|
T
1
|
=
|
(
θ
n
+
1
q
n
+
1
,
Π
h
*
θ
n
+
1
)
|
≤
C
∥
θ
n
+
1
∥
2
,
|
T
2
|
=
|
(
ρ
n
+
1
,
θ
n
+
1
)
|
≤
C
∥
ρ
n
+
1
∥
∥
θ
n
+
1
∥
,
|
T
3
|
=
|
(
ρ
n
+
1
q
n
+
1
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
|
≤
C
∥
ρ
n
+
1
∥
∥
θ
n
+
1
∥
.
To bound T
4, we use Lemma 10 to obtain
(92)
|
T
4
|
=
|
a
h
(
u
n
+
1
;
R
h
c
n
+
1
,
θ
n
+
1
)
h
h
-
a
h
(
u
h
n
+
1
;
R
h
c
n
+
1
,
θ
n
+
1
)
|
≤
C
(
∥
u
n
+
1
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
h
h
h
h
∥
u
n
+
1
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
+
h
∥
∇
·
(
u
n
+
1
-
u
h
n
+
1
)
∥
)
|
θ
n
+
1
|
1
≤
C
(
∥
u
n
+
1
-
u
n
∥
(
L
2
(
Ω
)
)
2
h
h
h
h
+
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
h
h
h
h
+
h
∥
∇
·
(
u
n
+
1
-
u
n
)
∥
h
h
h
h
∥
u
n
+
1
-
u
n
∥
(
L
2
(
Ω
)
)
2
+
h
∥
∇
·
(
u
n
-
u
h
n
)
∥
)
|
θ
n
+
1
|
1
.
Since
(93)
|
u
n
+
1
-
u
n
|
2
=
|
∫
t
n
t
n
+
1
u
t
d
s
‍
|
2
≤
Δ
t
∫
t
n
t
n
+
1
‍
|
u
t
|
2
d
s
,
hence,
(94)
∥
u
n
+
1
-
u
n
∥
(
L
2
(
Ω
)
)
2
2
≤
Δ
t
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
(
L
2
(
Ω
)
)
2
)
2
,
and, similarly,
(95)
∥
∇
·
(
u
n
+
1
-
u
n
)
∥
L
2
(
Ω
)
2
≤
Δ
t
∥
∇
·
u

t

∥
L
2
(
t
n
,
t
n
+
1
;
L
2
(
Ω
)
)
2
.
Hence,
(96)
|
T
4
|
≤
C
[
(
Δ
t
)
1
/
2
(
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
,
(
L
2
(
Ω
)
)
2
)
hhhhhhhhhh
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
,
(
L
2
(
Ω
)
)
2
)
+
h
∥
(
∇
·
u
)
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
(
Ω
)
)
)
gggg
+
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
gggg
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
,
(
L
2
(
Ω
)
)
2
)
+
h
∥
∇
·
(
u
n
-
u
h
n
)
∥
]
|
θ
n
+
1
|
.
Using the Cauchy-Schwarz inequality and (54), we obtain
(97)
|
T
5
|
=
|
(
u
n
+
1
·
∇
c
n
+
1
+
ϕ
∂
c
n
+
1
∂
t
-
ϕ
(
c
n
+
1
-
c
^
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
|
≤
C
∥
u
n
+
1
·
∇
c
n
+
1
+
ϕ
∂
c
n
+
1
∂
t
-
ϕ
(
c
n
+
1
-
c
^
n
)
Δ
t
∥
∥
θ
n
+
1
∥
.
Let σ
(
x
)
=
[
ϕ
(
x
)
2
+
u
n
+
1
(
x
)
2
]
1
/
2, so that
(98)
ϕ
∂
c
n
+
1
∂
t
+
u
n
+
1
·
∇
c
n
+
1
=
σ
∂
c
n
+
1
∂
τ
,
where τ approximates the characteristic unit vector s. Now using the same arguments given in [12], we have the following bound:
(99)
∥
σ
∂
c
n
+
1
∂
τ
-
ϕ
c
n
+
1
-
c
ˇ
n
Δ
t
∥
2
≤
C
Δ
t
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2
.
Hence, T
5 is bounded as follows:
(100)
|
T
5
|
≤
C
(
Δ
t
)
1
/
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
∥
θ
n
+
1
∥
.
To bound T
6, we proceed as follows:
(101)
|
T
6
|
=
|
(
ϕ
(
ρ
n
+
1
-
ρ
ˇ
n
)
Δ
t
,
Π
h
*
θ
n
+
1
)
|
≤
|
(
ϕ
(
ρ
n
+
1
-
ρ
ˇ
n
)
Δ
t
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
|
+
|
(
ϕ
(
ρ
n
+
1
-
ρ
ˇ
n
)
Δ
t
,
θ
n
+
1
)
|
=
I
1
+
I
2
.
Now I
1 can be written as
(102)
I
1
≤
|
(
ϕ
(
ρ
n
+
1
-
ρ
n
)
Δ
t
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
|
+
|
(
ϕ
(
ρ
n
-
ρ
ˇ
n
)
Δ
t
,
θ
n
+
1
-
Π
h
*
θ
n
+
1
)
|
.
Following the proof lines of Theorem 4.1 in [12], it can be shown that, for f
∈
H
1
(
Ω
), there exists a positive constant C independent of h and Δ
t such that
(103)
∥
f
-
f
ˇ
∥
≤
C
Δ
t
∥
∇
f
n
∥
.
A use of (103) yields
(104)
∥
(
ρ
n
-
ρ
ˇ
n
)
Δ
t
∥
≤
C
∥
∇
ρ
n
∥
.
It is easy to show that
(105)
∥
(
ρ
n
+
1
-
ρ
n
)
Δ
t
∥
≤
C
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
.
Using (52), (104), (105), and the Cauchy-Schwarz inequality, we obtain
(106)
I
1
≤
C
[
∥
θ
n
+
1
-
Π
h
*
θ
n
+
1
∥
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
h
h
h
h
h
×
(
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
+
∥
∇
ρ
n
∥
)
]
≤
C
h
[
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
g
g
g
g
l
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
+
∥
∇
ρ
n
∥
]
∥
∇
θ
n
+
1
∥
.
I
2 can be bounded as follows:
(107)
I
2
≤
|
(
ϕ
(
ρ
n
+
1
-
ρ
n
)
Δ
t
,
θ
n
+
1
)
|
+
|
(
ϕ
(
ρ
n
-
ρ
ˇ
n
)
Δ
t
,
θ
n
+
1
)
|
=
J
1
+
J
2
.
A use of Cauchy-Schwarz inequality yields
(108)
J
1
≤
C
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
∥
θ
n
+
1
∥
.
In order to bound J
2, we use the following result.
Ifη
∈
L
2
(
Ω
) and η
ˇ
(
x
)
=
η
(
x
ˇ
) with x
ˇ
=
x
-
r
(
x
)
Δ
t, for a nonzero function r
(
x
) such that r and ∇
·
r are bounded, then
(109)
∥
η
-
η
ˇ
∥
-
1
≤
C
∥
η
∥
Δ
t
,
where C is a positive constant independent of h and Δ
t; for a proof, we refer to [9, pp. 875]. Using (109), we obtain
(110)
J
2
≤
C
∥
(
ρ
n
-
ρ
^
n
)
Δ
t
∥
-
1
∥
θ
n
+
1
∥
1
≤
C
∥
ρ
n
∥
∥
θ
n
+
1
∥
1
.
This implies that
(111)
I
2
≤
C
(
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
+
∥
ρ
n
∥
)
∥
θ
n
+
1
∥
.
Now, using (106), (111), and (101), we obtain the following bound for T
6:
(112)
|
T
6
|
≤
C
(
(
Δ
t
)
-
1
/
2
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
h
h
h
h
h
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
+
h
∥
∇
ρ
n
∥
+
∥
ρ
n
∥
)
∥
θ
n
+
1
∥
.
Using the Cauchy-Schwarz inequality, we obtain
(113)
|
T
7
|
≤
|
(
(
θ
n
-
θ
^
n
)
Δ
t
,
Π
h
*
θ
n
+
1
-
θ
n
+
1
)
|
+
|
(
(
θ
n
-
θ
^
n
)
Δ
t
,
θ
n
+
1
)
|
≤
∥
(
θ
n
-
θ
^
n
)
Δ
t
∥
∥
Π
h
*
θ
n
+
1
-
θ
n
+
1
∥
+
∥
(
θ
n
-
θ
^
n
)
Δ
t
∥
-
1
∥
θ
n
+
1
∥
1
.
We use (103) and (109) to bound ∥
(
θ
n
-
θ
^
n
)
/
Δ
t
∥ and ∥
(
θ
n
-
θ
^
n
)
/
Δ
t
∥
-
1, respectively. For this, we need that u
h
n and its first derivative are bounded.
First let us make an induction hypothesis. Assume that there is a constant sayK
*
≥
2
K with ∥
u
~
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
K such that
(114)
∥
u
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
K
*
,
where u
~
h
n is the projection of u
h
n at t
=
t
n defined in (41). To bound ∥
∇
·
u
h
n
∥
∞, we use inverse inequality and (114):
(115)
∥
∇
·
u
h
n
∥
∞
Δ
t
≤
C
h
-
1
Δ
t
∥
u
h
n
∥
∞
≤
C
,
where we have used the assumption that Δ
t
=
O
(
h
).
Using (103), we have
(116)
∥
θ
n
-
θ
^
n
Δ
t
∥
≤
C
(
K
*
)
∥
∇
θ
n
∥
.
Using (109), we have
(117)
∥
θ
n
-
θ
^
n
Δ
t
∥
-
1
≤
C
(
K
*
)
∥
θ
n
∥
.
Now, using (52), (116), inverse inequality (48), and (117), we obtain the following bound for T
7:
(118)
|
T
7
|
≤
C
(
K
*
)
∥
θ
n
∥
∥
θ
n
+
1
∥
1
.
To bound T
8, we use maximum norm estimate of ∇
R
h
c (see (79)):
(119)
∥
R
^
h
c
n
-
R
ˇ
h
c
n
Δ
t
∥
≤
∥
R
h
c
∥
1
,
∞
∥
u
n
-
u
h
n
∥
≤
C
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
,
And, hence, using (54), T
8 is bounded as follows:
(120)
T
8
≤
C
∥
θ
n
+
1
∥
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
.
Substitute the estimates of T
1
,
T
2
⋯
T
8 in (90) and use nonsingular property of ϕ and a kick-back argument with the Young’s inequality to obtain
(121)
1
Δ
t
[
(
θ
n
+
1
,
Π
h
*
θ
n
+
1
)
-
(
θ
n
,
Π
h
*
θ
n
+
1
)
]
≤
C
(
K
*
)
[
Δ
t
-
1
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
hhhhhhhhh
+
∥
θ
n
+
1
∥
2
+
∥
θ
n
+
1
∥
1
2
hhhhhhhhh
+
∥
θ
n
∥
2
+
∥
ρ
n
+
1
∥
2
hhhhhhhhh
+
Δ
t
(
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
hhhhhggggghhhhhl
+

∥
(
∇
·
u
)
t
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2

hhhhhggggghhhhlh
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
)
hhhhhhhhh
+
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
2
hhhhhhhhh
+
h
2
∥
∇
·
(
u
n
-
u
h
n
)
∥
2
hhhhhhhhh
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
+
∥
ρ
n
∥
2
+
h
2
∥
∇
ρ
n
∥
2
]
.
Using (43) and (44), we obtain
(122)
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
≤
C
[
∥
c
n
-
c
h
n
∥
+
h
(
∥
u
n
∥
(
H
1
(
Ω
)
)
2
+
∥
p
n
∥
1
2
)
]
≤
C
[
∥
θ
n
∥
+
∥
ρ
n
∥
+
h
(
∥
u
n
∥
(
H
1
(
Ω
)
)
2
+
∥
p
n
∥
1
2
)
]
,
∥
∇
·
(
u
n
-
u
h
n
)
∥
≤
C
h
∥
∇
·
u
n
∥
1
.
Substitute (122) in (121) to obtain
(123)

∥
|
θ
n
+
1
|
∥
2
-
∥
|
θ
n
|
∥
2

≤
C
(
K
*
)
[
Δ
t
(
∥
θ
n
+
1
∥
2
+
∥
θ
n
∥
2
+
∥
ρ
n
∥
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
ffffffffffffffffffffhfh
∥
θ
n
+
1
∥
2
+
∥
ρ
n
+
1
∥
2
+
h
2
∥
∇
ρ
n
∥
2
)
fffffffffffffffhff
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
fffffffffffffffffffhhhhhhhl
+
∥
∇
·
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2
fffffffffffffffffffhhhhhhhl
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
)
fffffffffffffffhff
+
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
ffffffffffffffhfff
+
h
2
Δ
t
(
h
2
∥
∇
·
u
n
∥
1
2
+
∥
u
n
∥
(
H
1
(
Ω
)
)
2
2
fffffffffffffffffhffhhhhl
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
∥
u
n
∥
(
H
1
(
Ω
)
)
2
2
+
∥
p
n
∥
1
2
)
]
.
Taking summation over n
=
0
⋯
m
-
1, we obtain
(124)

∥
|
θ
m
|
∥
2
-
∥
|
θ
0
|
∥
2

≤
C
(
K
*
)
[
∑
n
=
0
m
-
1
‍
{
Δ
t
(
∥
θ
n
+
1
∥
2
+
∥
θ
n
∥
2
+
∥
ρ
n
∥
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
hhhhhhhhhhhhhhhhhl
+
∥
ρ
n
+
1
∥
2
+
h
2
∥
∇
ρ
n
∥
2
)
hhhhhhhhhhhhh
.
hh
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
hhhhhhhhhhhhhhhhhhhhhh
.
+
∥
∇
·
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2
hhhhhhhhhhhhhhhhhhhhhh
.
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
)
hhhhhhhhhhhhhhh
.
+
∥
∂
ρ
∂
t
∥
L
2
(
t
n
,
t
n
+
1
;
L
2
)
2
hhhhhhhhhhhhhhh
.
+
h
2
Δ
t
(
h
2
∥
∇
·
u
n
∥
1
2
+
∥
u
n
∥
(
H
1
(
Ω
)
)
2
2
hhhhhhhhhhhhhhhhhhhhl
∑
n
=
0
m
-
1
‍
∥
∂
2
c
∂
τ
2
∥
L
2
(
t
n
,
t
n
+
1
,
L
2
)
2
∥
∇
·
u
n
∥
1
2
∥
u
n
∥
(
H
1
(
Ω
)
)
2
2
+
∥
p
n
∥
1
2
)
}
]
.
Now, using discrete Gronwall’s, equivalence of the norms (53), and the estimates of ρ, we obtain
(125)
∥
θ
m
∥
2
≤
C
(
K
*
)
×
[
∥
θ
0
∥
2
+
h
4
(
∥
(
c
-
c
~
)
q
∥
L
∞
(
0
,
T
;
H
1
)
2
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
jjjjjjjjjjjjjjjjjjjjjjjhj
+
∥
ψ
∂
c
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
∥
L
∞
(
0
,
T
;
H
2
)
2
jjjjjjjjjjjjjjjjjjjjjjhjj
+
∥
∇
·
u
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
t
∥
L
∞
(
0
,
T
;
H
2
)
2
jjjjjjjjjjjjjjjjjjjhjjjjj
+
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
)
hh
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
0
,
T
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
jjjjjjjjjjhjhjjjjjjk
+
∥
∇
·
u
t
∥
L
2
(
0
,
T
;
L
2
)
2
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
)
jhjjhh
+
h
2
(
∥
u
∥
L
∞
(
0
,
T
;
(
H
1
(
Ω
)
)
2
)
2
jhjjhhhh
.
hh
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
∥
u
∥
L
∞
(
0
,
T
;
(
H
1
(
Ω
)
)
2
)
2
+
∥
p
∥
L
∞
(
0
,
T
;
H
1
)
2
)
]
.
Now, it remains to show the induction hypothesis (114). Using (23) and (24), we have
(126)
∥
u
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
∥
u
h
n
-
u
~
h
n
∥
(
L
∞
(
Ω
)
)
2
+
∥
u
~
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
C
(
log
⁡
1
h
)
1
/
2
∥
u
h
n
-
u
~
h
n
∥
H
(
div
⁡
;
Ω
)
+
K
.
Using ∥
∇
·
(
u
h
n
-
u
~
h
n
)
∥
=
0, we have
(127)
∥
u
h
n
-
u
~
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
C
(
log
⁡
1
h
)
1
/
2
∥
u
h
n
-
u
~
h
n
∥
(
L
2
(
Ω
)
)
2
.
Now using (42) and (125), we obtain for small h
(128)
∥
u
h
n
∥
(
L
∞
(
Ω
)
)
2
≤
C
(
K
*
)
log
⁡
1
h
(
h
+
Δ
t
)
+
K
≤
2
K
.
Here, we have used Δ
t
=
O
(
h
) and h
log
⁡
(
1
/
h
)
→
0 as h
→
0 and this proves our induction hypothesis (114).
Now combine the estimates ofρ and θ and use triangle inequality to complete the rest of the proof.Using (43) and Theorem 11, we obtain the following error estimates for velocity as well as pressure.Theorem 12.
Assume that the triangulation𝒯
h is quasiuniform. Let (
u
,
p
) and (
u
h
,
p
h
) be, respectively, the solutions of (1)-(2) and (30) and let c
h
(
0
)
=
c
0
,
h
=
R
h
c
(
0
). Further assume that Δ
t
=
O
(
h
). Then for sufficiently small h there exists a positive constant C
(
T
) independent of h but dependent on the bounds of κ
-
1 and μ such that
(129)
max
⁡
0
≤
n
≤
N
∥
u
n
-
u
h
n
∥
(
L
2
(
Ω
)
)
2
2
≤
C
[
h
4
(
∥
(
c
-
c
~
)
q
∥
L
∞
(
0
,
T
;
H
1
)
2
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
hhhhhhhhhh
.
+
∥
ψ
∂
c
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
∥
L
∞
(
0
,
T
;
H
2
)
2
hhhhhhhhhh
.
+
∥
∇
·
u
∥
L
∞
(
0
,
T
;
H
1
)
2
+
∥
c
t
∥
L
∞
(
0
,
T
;
H
2
)
2
hhhhhhhhhh
.
+
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
)
hhhhhhhh
+
(
Δ
t
)
2
(
∥
u
t
∥
L
2
(
0
,
T
;
(
L
2
(
Ω
)
)
2
)
2
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
hhhhhhjjjjjjjjjhhhhh
+
∥
∇
·
u
t
∥
L
2
(
0
,
T
;
L
2
)
2
+
∥
∂
2
c
∂
τ
2
∥
L
2
(
0
,
T
,
L
2
)
2
)
hhhhhhhh
∥
ψ
∂
2
c
∂
t
∂
s
∥
L
∞
(
0
,
T
;
H
1
)
2
+
h
2
(
∥
u
∥
L
∞
(
0
,
T
;
(
H
1
(
Ω
)
)
2
)
2
+
∥
p
∥
L
∞
(
0
,
T
;
H
1
)
2
)
]
.


## 4. Numerical Experiments

For our numerical experiments, we consider (1)–(6), with q
=
q
+
-
q
- and g
(
x
,
t
,
c
)
=
c
-
q
+
-
c
q
-, where c
- is the injection concentration and q
+ and q
- are the production and injection rates, respectively.Experimentally, it has been observed that the velocity is much smoother in time compared to the concentration. It was suggested in [27] that, for a good approximation to the concentration, one should take larger time step for the pressure equation than the concentration equation. Let 0
=
t
0
<
t
1
<
⋯
t
M
=
T a given partition of the time interval (
0
,
T
] with step length Δ
t
m
=
t
m
+
1
-
t
m for the pressure equation and 0
=
t
0
<
t
1
<
⋯
t
N
=
T a given partition of the time interval (
0
,
T
] with step length Δ
t
n
=
t
n
+
1
-
t
n for the concentration equation. We denote C
n
≈
c
h
(
t
n
), C
m
≈
c
h
(
t
m
), U
m
≈
u
h
(
t
m
), and P
m
≈
p
h
(
t
m
).If concentration stept
n relates to pressure steps by t
m
-
1
<
t
n
≤
t
m, we require a velocity approximation at t
=
t
n, which will be used in the concentration equation, based on U
m
-
1 and earlier values. We define a velocity approximation [12, pp. 81] at t
=
t
n by
(130)
E
U
n
=
(
1
+
t
n
-
t
m
-
1
t
m
-
1
-
t
m
-
2
)
U
m
-
1
-
t
n
-
t
m
-
1
t
m
-
1
-
t
m
-
2
U
m
-
2
,
for
m
≥
2
,
E
U
n
=
U
0
,
for
m
=
1
.
Then the combined time stepping procedure is defined as follows: find C
:
{
t
0
,
t
1
,
…
,
t
N
}
→
M
h and (
U
,
P
)
:
{
t
0
,
t
1
,
…
,
t
M
}
→
U
h
×
W
h such that
(131)
(
κ
-
1
μ
(
C
m
)
U
m
,
γ
h
v
h
)
+
b
(
γ
h
v
h
,
P
m
)
=
0
∀
v
h
∈
U
h
,
(132)
(
∇
·
U
m
,
w
h
)
=
(
q
+
-
q
-
,
w
h
)
∀
w
h
∈
W
h
,
m
≥
0
,
(133)
(
ϕ
C
n
+
1
-
C
^
n
Δ
t
,
Π
h
*
χ
h
)
+
a
h
(
E
U
n
+
1
;
C
n
+
1
,
χ
h
)
+
(
q
-
C
n
+
1
,
Π
h
*
χ
h
)
=
(
q
+
c
-
,
Π
h
*
χ
h
)
∀
χ
h
∈
M
h
,
where C
^
n
=
C
n
(
x
-
(
E
U
n
/
ϕ
)
Δ
t
).To solve the pressure equations, that is, (131) and (132), we use the mixed FVEM, and for the concentration equation (133), we use the standard FVEM.For the test problems, we have taken data from [28]. The spatial domain is  Ω
=
(
0,1000
)
×
(
0,1000
)
ft
2, the time period is [
0,3600
] days, and viscosity of oil is   μ
(
0
)
=
1.0 cp. The injection well is located at the upper right corner (
1000,1000
) with the injection rate q
+
=
30
ft
2
/
day and injection concentration c
-
=
1.0. The production well is located at the lower left corner with the production rate q
-
=
30
ft
2
/
day and the initial concentration is c
(
x
,
0
)
=
0. For time discretization, we take Δ
t
p
=
360 days and Δ
t
c
=
120 days; that is, we divide each pressure time interval into three subintervals.Test  1. We assume that the porous medium is homogeneous and isotropic. The permeability is κ
=
80. The porosity of the medium is ϕ
=
.
1 and the mobility ratio between the resident and injected fluid is M
=
1. Furthermore, we assume that the molecular diffusion is d
m
=
1 and the dispersion coefficients are zero. In our numerical simulation, we divide the spatial domain into 20 equal subdivisions along both x and y axis. For time discretization, we take Δ
t
p
=
360 days and Δ
t
c
=
120 days; that is, we divide each pressure time interval into three subintervals.The surface and contour plots for the concentration att
=
3 and t
=
10 years are presented in Figures 6 and 7, respectively. Since only molecular diffusion is present and viscosity is also independent of the velocity, Figure 6 shows that the velocity is radial and the contour plots for the concentration are circular until the invading fluid reaches the production well. Figure 7 shows that when these plots are reached at production well, the invading fluid continues to fill the whole domain until c
=
1.Contour (a) and surface plot (b) in Test1 at t
=
3 years.
(a)
(b)Contour (a) and surface plot (b) in Test1 at t
=
10 years.
(a)
(b)Test  2. In this test we consider the numerical simulation of a miscible displacement problem with discontinuous permeability. Here, the data is the same as given in Test 1 except the permeability of the medium κ
(
x
). We take κ
=
80 on the subdomain Ω
L
:
=
(
0,1000
)
×
(
0,500
) and κ
=
20 on the subdomain Ω
U
:
=
(
0,1000
)
×
(
500,1000
). The contour and surface plot at t
=
3 and t
=
10 years are given in Figures 8 and 9, respectively. Figures 8 and 9 show that when the injecting fluid reaches the lower half domain, it starts moving much faster in the horizontal direction on this domain compared to the low permeability domain, that is, upper half domain. We observe that one should put the production well in a low permeability zone to increase the area swept by the injected fluid.Contour (a) and surface plot (b) in Test2 at t
=
3 years.
(a)
(b)Contour (a) and surface plot (b) in Test2 at t
=
10 years.
(a)
(b)Order of Convergence. In order to verify our theoretical results we also compute the order of convergence for the concentration for this particular test problem. We compute the order of convergence in L
2 norm. To discretize the time interval [
0
,
T
], we take uniform time step Δ
t
=
360 days for pressure and concentration equations. The computed order of convergence is given in Figure 10. Note that the computed order of convergence matches with the theoretical order of convergence derived in Theorem 11.Figure 10
Order of convergence inL
2 norm.Note. This paper has been presented in the International Conference on Numerical Analysis and Applications which was held at Bulgaria during June 15–20, 2012. Moreover, some of the results without proof have been published in the proceeding of that conference, kindly see [29].



---

*Source: 245086-2014-01-19.xml*