A Second Order Characteristic Mixed Finite Element Method for Convection Diffusion Reaction Equations ()
1. Introduction
Let
be a bounded domain in
with Lipschitz boundary
, where
. Let T be a positive constant. In this paper, we will consider the following linear convection-diffusion-reaction equations: find
such that
(1)
Here,
is the convection vector field;
is the reaction function;
and
are given scalar functions;
is a constant and
is the outward unit normal vector to
;
denotes the diffusion matrix function, where
is the space of symmetric
matrices, such that
(2)
where
are positive constants.
Typically, these Equations (1) express the general mathematical model incorporating different types of transport phenomena in engineering and applied sciences, such as the dispersal of a pollutant through a moving viscous medium (e.g., a river or the atmosphere, [1] ), currents in semiconductor devices [2] , and airflow past an aerofoil (see [3] , for example). When the diffusive term is smaller than the convective one, these equations are the so-called convection dominated problems (see [4] ).
In many practical convection-diffusion processes, convection essentially dominates diffusion (e.g., in some financial models [5] ), and although the governing differential equation is parabolic, it displays several characteristics of hyperbolic problems. When applied to these problems, standard finite element and finite difference methods usually exhibit some combination of nonphysical oscillation and excessive numerical dispersion [6] [7] . It is therefore logical to design numerical procedures that incorporate the parabolic/hyperbolic nature of these problems. One such method is the modified method of characteristics (MMOC) which was first formulated for a scalar parabolic equation by J. Douglas and T. F. Russell in [8] and then extended by Russell [9] to nonlinear coupled systems in two and three spatial dimensions. Similar schemes had been defined by Pironneau [10] for the incompressible Navier-Stokes equations and by Süli [11] and Morton, Priestley, and Süli [12] for first-order hyperbolic equations, with the latter technique being referred to as the Euler characteristic Galerkin method. The intent of the method is to obtain accurate approximations to convection- dominated problems. Basically, in the modified method of characteristics, the time derivative and the convection term are combined as a directional derivative. In other words, the procedure involves time stepping along the characteristics, allowing us to use large, accurate time steps.
Mixed finite element method has been proven to be an effective numerical method for solving fluid problems. It has an advantage to approximate the unknown variable and its diffusive flux simultaneously. There are many research articles on this method ( [13] [14] [15] [16] ). An algorithm combining the mixed finite-element method and the modified method of characteristics was first applied to the miscible displacement problem in porous media by Ewing, Russell, and Wheeler [17] . Then, the scheme had been extended by Wheeler and Dawson [18] to advection-diffusion-reaction problems. Numerical results have verified that large, accurate time steps are possible, and sharp fronts have been resolved (without oscillations or numerical diffusion) by coarser grids that standard procedures can use.
Arbogast and Wheeler [19] defined a characteristics-mixed method to approximate the solution of an advection-dominated transport problem. It used a characteristic approximation that is similar to that of MMOC-Galerkin method to handle advection in time and a lowest order mixed finite element spatial approximation for diffusion term. Piecewise constants were in the space of test function, so mass is conserved element-by-element. It was proved finally that the method was optimally convergent to order 1 in time and at least suboptimally convergent to order 3/2 in space. In [20] , we have considered a combined numerical approximation for incompressible miscible displacement in porous media. Standard mixed finite element was used for Darcy velocity equation and a characteristics-mixed finite element method was presented for approximating the concentration equation. Characteristic approximation was applied to handle the convection term, and a lowest order mixed finite element spatial approximation was adopted to deal with the diffusion term. Thus, the scalar unknown concentration and the diffusive flux can be approximated simultaneously. This approximation conserves mass globally. The optimal L2-norm error estimates were derived. Then, we extended this method to the slightly compressible miscible displacement problem in [21] .
It should be pointed out that the works mentioned above which involved the characteristic method all gave one order accuracy in time increment
. That is to say, the first order characteristic method in time was analyzed. As for higher order characteristic method in time, Rui and Tabata [22] used the second order Runge-Kutta method to approximate the material derivative term for convection-diffusion problems. The scheme presented was of second order accuracy in time increment
, symmetric and unconditionally stable. Optimal error estimates were proved in the framework of L2-theory. Numerical analyses of convection-diffusion-reaction problems with higher order characteristic/finite elements were analyzed in [23] [24] , which extended the work [22] . The
error estimates of second order in time increment
were obtained.
The goal of this paper is to present a second order characteristic mixed finite element method in time increment to handle the material derivative term of (1). It is organized as follows. In Section 2, we formulate an approximate scheme that combines the second order characteristic finite element method for the material derivative term and mixed finite element method for the diffusion term. In Section 3, we prove the stability of the combined approximate scheme. In Section 4, we derive the L2-norm error estimates for both the unknown scalar variable and its flux. The scheme is of second order accuracy in time increment, symmetric, and unconditionally stable.
2. Formulation of the Method
In this paper, we adopt notations and norms of usual Sobolev spaces. Moreover, we adopt some notations for the functional spaces involved, which were used in [22] [23] [24] . For a Banach space X and a positive integer m, spaces
and
are abbreviated as
and
, respectively, and endowed with norms
where
denotes the j-th derivative of
with respect to time. The Banach space
is defined by
equipped with the norm
Similar spaces
and norms are considered for the boundary sets ΓR and ΓD. Denote by
the closed subspace of
defined by
and [25]
2.1. The Characteristic Lines
Now, we define the characteristic lines associated with vector field
and recall some classical properties satisfied by them. Thus, for given
, the characteristic line through
is the vector function
solving the initial value problem
(3)
Next, assuming they exist, we denote by
(respectively, by L) the gradient of
(respectively, of
) with respect to the space variable x, i.e.,
(4)
We adopted some propositions and lemmas from [23] .
Proposition 1. If
for
an integer, then
and it is
with respect to the x variable.
In order to compute second order approximations of matrices
and
, we need the following equations:
(5)
Proposition 2. If
, then
(6)
Proposition 3. If
, then
satisfies the Taylor expansion
(7)
and its inverse,
, satisfies the Liouville’s theorem
(8)
By using the Liouville’s theorem and the chain rule, we obtain
(9)
Proposition 4. If
, then
(10)
Proposition 5. If
, then
satisfies
(11)
2.2. Variational Formulation
From the definition of the characteristic lines and by using the chain rule, it follows that
(12)
By introducing the flux
and using (12), we rewrite equation (1.1) at point
and time
as follows
(13)
Before giving a week formulation of (13), we adopted a lemma from [23] , which can be considered as a Green’s formula.
Lemma 6. Let
be an invertible vector valued function. Let
and assume that
Then
(14)
with
a vector valued function and
a scalar function.
Now, we can multiply (13) by test functions
and
, integrate in
respecetively, and apply the usual Green’s formula and (14) with
, obtaining
(15)
Lemma 7. [23] Let
be an invertible vector valued function satifying
. Let
and assume that
Then
(16)
with
and
, where
is the outward unit normal vector to
.
Now, replacing in (15) formula (16) with
, and replacing the Robin condition, we have
(17)
2.3. The Combined Approximate Scheme
We now present our time-stepping procedure for (17). Let N be the number of time steps,
be the time step, and
for
. We will use the notation
for a function. Moreover, for
, we define
(18)
By fixing
in (17) and using a Crank-Nicolson method [26] with respect to
. Thus, from (18) we have
(19)
By using (8) and (11), we see that
(20)
Taking (20) into (19), we can obtain
(21)
We propose two explicit numerical schemes to approximate
:
(22)
A similar notation to the one in Section 2.2 is used for the Jacobian of
, namely,
Now, we restate three lemmas concerning properties of the characteristic line approximations. For this, we require the time step to be bounded and the velocity to satisfy the following assumption.
Claim 1. The velocity field
and satisfies
on
.
Lemma 8. [23] Under Claim 1, if
, we can see that
(23)
Lemma 9. [23] Under Claim 1, if
, we have
(24)
Corollary 1. Under the assumptions of Lemma 9,
, we have
(25)
Lemma 10. [22] Under Claim 1, if
and
, then there exists a positive constant
such that
(26)
where
.
Thus, in the case where the characteristic lines and their gradients are not explicitly known, we propose the following time approximation of (21)
(27)
The time difference approximation (27) will be combined with a standard Galerkin finite element and mixed finite element in the space for
and
, respectively [27] [28] . We discrete
in space on a quasi-uniform finite element mesh
of
with maximal element diameter
. For
, we denote as
and
similarly. Let
be finite element spaces with index k and l, respectively.
We define a bilinear form
on
and a linear form
on
for
by
(28)
where
.
Then, the fully discrete scheme reads: Given
, find
such that
(29)
Throughout the analysis, K will denote a generic positive constant, independent of
. Similarly,
will denote a generic small positive constant.
3. Stability of the Approximate Scheme
In this section, we derive the stability of the approximate scheme (29). In order to develop the stability, some assumptions on the different terms of (1) are required.
Claim 2. The velocity field
and satisfies
on
.
Remark. Throughout this paper
denotes the maximum between the positive constant appearing in Lemma 10 and the norm of the velocity in
.
Claim 3. The diffusion matrix coefficients,
, belong to
. More- over, A is a positive definite symmetric
matrix and there exists a strictly positive constant
which is a uniform lower bound for the eigenvalues of
.
As a consequence of Claim 3, there exists a unique positive definite symmetric
matrix function
, such that
and
. Let us in-
troduce the constant
. Clearly, Claim 3, we have
(30)
Claim 4. The reaction function,
, satisfies
in
, where
is a constant.
Under the previous claims, let
Claim 5. The source function
. In Robin boundary condition,
and
.
For a given series of functions
, we define the following norms
Similar definitions are considered for functional spaces
and
associated to the Robin boundary condition. Moreover, we define
Lemma 11. Let the above Claims 2 - 5 and
be assumed. If
be the solution of (29). Then it holds that
(31)
the constant c is given by
Proof. Substituting
into (28), we have
(32)
Lemma 4.1 in [23] implies that
(33)
and
(34)
Next, by using
, we obtain
Then when
and
, we have
(35)
Similarly, for
we obtain the estimate
(36)
Analogous computations to term
give
(37)
For the boundary integral term
, we first use some properties of the inner product in the space
and the inequality
to get the estimate
for
. Thus, we obtain
(38)
Then, by summing up from (33) to (38), inequality (31) follows.
By using Lemma 11 and following the arguments to Lemmas 5.6, 5.7 and Theorem 5.8 in [23] , we can get the following stability theorem:
Theorem 12 (Stability Theorem) Let the above Claims 2 - 5 assumed, and
be the solution of (29) subject to the initial value
. Then there exist two positive constants c and
, such that, for
, we have
(39)
4. Error Estimate Theorem
Now, we turn to derive a priori error estimate in L2-norm for the solutions of (29). In order to state error estimates, we need two following Lagrange interpolation operators ( [29] [30] )
, and
.
Lemma 13. There exist positive constants
and
, independent of
and
respectively, such that
(40)
and
(41)
Let
Corresponding to
and
, we introduce a bilinear form
on
and a linear form
on
for
as follows
(42)
If
and
are the solutions of (1), we have for
(43)
We decompose
as
(44)
where
.
In order to estimate the terms on the right-hand side of (44), we adopt the following lemmas.
Lemma 14. [23] Assume the above Claims 2 - 5 hold, and that the coefficients of the problem (1) satisfy
,
and that
. Let the solution of (43) satisfy
,
,
. Then, for each
there exist two functions
and
,
such that
(45)
. Moreover,
,
and the following estimates hold:
(46)
where
denotes a constant independent of
.
Lemma 15. [23] Assume the above Claims 2 - 5 hold, and that the coefficients of the problem (1) satisfy
,
and
. Then, for each
there
exist function
and
, such that
(47)
Moreover,
and the following estimates hold:
(48)
where
denotes a constant independent of
.
Now, we turn to bound the third term on the right-hand side of (44).
Lemma 16. Assume the above Claims 2 - 5 hold, and that the coefficients of the problem (1) satisfy
and
. There exists
(49)
with
,
a positive constant.
Proof. From the definition of
, we have
(50)
From the arguments of Lemma 4.1 in [24] , we get similarly the bounds of the above terms
as the followings:
(51)
(52)
(53)
To estimate
, we divide it into three parts
(54)
where
Using the transformation
, we have
(55)
Now, we replace in
equality
, where
(56)
with
, and the function
is defined by
Thus, we get
(57)
Moreover, we have
(58)
By considering together from (54) to (58), we can state
(59)
Similar arguments, i.e. Lemma 5.4 in [23] lead to
(60)
Using these inequalities,
and
can be bounded as follows
(61)
Gathering (50), (51), (52), (53), (59), and (61) together, we complete the proof.
We now turn to estimate
and
. From the definition of
and
, we see
(62)
From Lemma 11, we obtain the lower bound for (62)
(63)
By jointly considering the lower bound (63), the upper bounds given in Lemmas 14 - 16, we deduce
(64)
Analogous computations to those developed in Lemma 16 give
(65)
Putting (65) into (64), multiplying (64) by
, summing it about time form
and
, taking the initial values
and using discrete Gronwall’s inequality, we can derive the following error estimate:
Theorem 17 (Error Estimate) Let Claims 2 - 5 be assumed. Let
be the solution of (1) with
,
be the solution of (29) subject to the initial values
. Then there exists a positive constant
independent of
and
such that
(66)
Acknowledgement
This research was supported by the NSF of China (Grant Nos. 11271231 and 11301300).