A Second Order Characteristic Mixed Finite Element Method for Convection Diffusion Reaction Equations

Abstract

A combined approximate scheme is defined for convection-diffusion-reaction equations. This scheme is constructed by two methods. Standard mixed finite element method is used for diffusion term. A second order characteristic finite element method is presented to handle the material derivative term, that is, the time derivative term plus the convection term. The stability is proved and the L2-norm error estimates are derived for both the scalar unknown variable and its flux. The scheme is of second order accuracy in time increment, symmetric, and unconditionally stable.

Share and Cite:

Sun, T. (2017) A Second Order Characteristic Mixed Finite Element Method for Convection Diffusion Reaction Equations. Journal of Applied Mathematics and Physics, 5, 1301-1319. doi: 10.4236/jamp.2017.56109.

1. Introduction

Let Ω be a bounded domain in R d ( d = 2,3 ) with Lipschitz boundary Γ = Γ D Γ R , where Γ D Γ R = . Let T be a positive constant. In this paper, we will consider the following linear convection-diffusion-reaction equations: find ϕ : Ω × ( 0, T ) R such that

{ ϕ ( x , t ) t + v ( x , t ) ϕ ( x , t ) d i v ( A ( x ) ϕ ( x , t ) ) + r ( x ) ϕ ( x , t ) = f ( x , t ) , in Ω × ( 0 , T ) , ϕ ( x , t ) = 0 , on Γ D × ( 0 , T ) , α ϕ ( x , t ) + A ( x ) ϕ ( x , t ) n ( x ) = g ( x , t ) , on Γ R × ( 0 , T ) , ϕ ( x , 0 ) = ϕ 0 ( x ) , in Ω . (1)

Here, v : Ω ¯ × [ 0, T ] R d is the convection vector field; r : Ω ¯ × [ 0, T ] R is the reaction function; f : Ω ¯ × [ 0, T ] R and g : Γ R × [ 0, T ] R are given scalar functions; α > 0 is a constant and n is the outward unit normal vector to Γ ; A : Ω ¯ S d denotes the diffusion matrix function, where S d is the space of symmetric d × d matrices, such that

a 0 ξ 2 i , j = 1 d A i j ξ i ξ j a 1 ξ 2 , ξ R d , (2)

where a 0 , a 1 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 Δ t . 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 Δ t , 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 l ( L 2 ) error estimates of second order in time increment Δ t 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 C m ( [ 0, T ] , X ) and H m ( ( 0, T ) , X ) are abbreviated as C m ( X ) and H m ( X ) , respectively, and endowed with norms

φ C m ( X ) : = max t [ 0 , T ] { max j = 0 , , m φ ( j ) ( t ) X } , φ H m ( X ) : = ( 0 T j = 0 m φ ( j ) ( t ) X 2 d t ) 1 2 ,

where φ ( j ) denotes the j-th derivative of φ with respect to time. The Banach space Z m is defined by

Z m = { f C j ( H m j ( Ω ) ) ; j = 0 , , m } ,

equipped with the norm φ Z m : = max { φ C j ( H m j ) ; 0 j m } . Similar spaces

and norms are considered for the boundary sets ΓR and ΓD. Denote by H Γ D 1 ( Ω ) the closed subspace of H 1 ( Ω ) defined by

H Γ D 1 ( Ω ) : = { φ H 1 ( Ω ) , φ | Γ D = 0 } ,

and [25]

H ( d i v ; Ω ) : = { v ( L 2 ( Ω ) ) d , d i v v L 2 ( Ω ) } .

2.1. The Characteristic Lines

Now, we define the characteristic lines associated with vector field v and recall some classical properties satisfied by them. Thus, for given ( x , t ) Ω ¯ × [ 0, T ] , the characteristic line through ( x , t ) is the vector function X e ( x , t ; ) solving the initial value problem

X e τ ( x , t ; τ ) = v ( X e ( x , t ; τ ) , τ ) , X e ( x , t ; τ ) = x . (3)

Next, assuming they exist, we denote by F e (respectively, by L) the gradient of X e (respectively, of v ) with respect to the space variable x, i.e.,

( F e ) r s ( x , t ; τ ) : = ( X e ) r x s ( x , t ; τ ) , L r s ( x , t ) : = v r x s ( x , t ) . (4)

We adopted some propositions and lemmas from [23] .

Proposition 1. If v C 0 ( C n ( Ω ¯ ) ) for n 1 an integer, then X e C 0 ( Ω ¯ × [ 0, T ] × [ 0, T ] ) and it is C n with respect to the x variable.

In order to compute second order approximations of matrices F e and F e 1 , we need the following equations:

{ F e τ ( x , t ; τ ) = L ( X e ( x , t ; τ ) , τ ) F e ( x , t ; τ ) , 2 F e τ 2 ( x , t ; τ ) = ( v t + L v ) ( X e ( x , t ; τ ) , τ ) F e ( x , t ; τ ) . (5)

Proposition 2. If v C 0 ( C 1 ( Ω ¯ ) ) , then

F e ( x , t ; τ ) e v C 0 ( C 1 ( Ω ¯ ) ) | τ t | , x Ω , t , τ [ 0, T ] . (6)

Proposition 3. If v C 0 ( C 2 ( Ω ¯ ) ) C 1 ( C 1 ( Ω ¯ ) ) , then F e satisfies the Taylor expansion

F e ( x , t ; s ) = I + ( s t ) L ( x , t ) + s t ( τ s ) ( v t + L v ) ( X e ( x , t ; τ ) , τ ) F e ( x , t ; τ ) d τ , (7)

and its inverse, F e 1 , satisfies the Liouville’s theorem

F e 1 ( x , t ; s ) = I + ( t s ) L ( X e ( x , t ; s ) , s ) s t ( τ t ) ( v t + L v ) ( X e ( x , t ; τ ) , τ ) F e ( X e ( x , t ; s ) , s ; τ ) d τ . (8)

By using the Liouville’s theorem and the chain rule, we obtain

{ τ d e t F e 1 ( x , t ; τ ) = d e t F e 1 ( x , t ; τ ) d i v v ( X e ( x , t ; τ ) , τ ) , 2 τ 2 d e t F e 1 ( x , t ; τ ) = d e t F e 1 ( x , t ; τ ) ( ( d i v v ) 2 ( X e ( x , t ; τ ) , τ ) d i v ( v t + L v ) ( X e ( x , t ; τ ) , τ ) ( L L T ) ( X e ( x , t ; τ ) , τ ) ) . (9)

Proposition 4. If v C 0 ( C 1 ( Ω ¯ ) ) , then

d e t F e 1 ( x , t ; τ ) e v C 0 ( C 1 ( Ω ¯ ) ) | τ t | , x Ω , t , τ [ 0, T ] . (10)

Proposition 5. If v C 0 ( C 2 ( Ω ¯ ) ) C 1 ( C 1 ( Ω ¯ ) ) , then d e t F e 1 satisfies

d e t F e 1 ( x , t ; s ) = 1 ( s t ) d i v v ( x , t ) + s t ( τ s ) 2 τ 2 ( d e t F e 1 ) ( x , t ; τ ) d τ . (11)

2.2. Variational Formulation

From the definition of the characteristic lines and by using the chain rule, it follows that

d ϕ d τ ( X e ( x , t ; τ ) , τ ) = ϕ t ( X e ( x , t ; τ ) , τ ) + v ( X e ( x , t ; τ ) , τ ) ϕ ( X e ( x , t ; τ ) , τ ) . (12)

By introducing the flux σ = A ϕ and using (12), we rewrite equation (1.1) at point X e ( x , t ; τ ) and time τ as follows

{ d ϕ d τ ( X e ( x , t ; τ ) , τ ) d i v σ ( X e ( x , t ; τ ) , τ ) + r ( X e ( x , t ; τ ) ) ϕ ( X e ( x , t ; τ ) , τ ) = f ( X e ( x , t ; τ ) , τ ) , σ ( X e ( x , t ; τ ) , τ ) = A ( X e ( x , t ; τ ) ) ϕ ( X e ( x , t ; τ ) , τ ) . (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 X : Ω ¯ X ( Ω ) ¯ , X C 2 ( Ω ¯ ) be an invertible vector valued function. Let F = X and assume that F 1 C 1 ( Ω ¯ ) . Then

Ω d i v w ( X ( x ) ) ψ ( x ) d x = Γ F T ( x ) n ( x ) w ( X ( x ) ) ψ ( x ) d Γ Ω F 1 ( x ) w ( X ( x ) ) ψ ( x ) d x Ω d i v F T ( x ) w ( X ( x ) ) ψ ( x ) d x , (14)

with w H 1 ( X ( Ω ) ) a vector valued function and ψ H 1 ( Ω ) a scalar function.

Now, we can multiply (13) by test functions ψ H Γ D 1 ( Ω ) and χ H ( d i v ; Ω ) , integrate in Ω respecetively, and apply the usual Green’s formula and (14) with X ( x ) = X e ( x , t ; τ ) , obtaining

{ Ω d ϕ d τ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x + Ω F e 1 ( x , t ; τ ) σ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x + Ω d i v F e T ( x , t ; τ ) σ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x + Ω r ( X e ( x , t ; τ ) ) ϕ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x Γ R F e T ( x , t ; τ ) n ( x ) σ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d Γ = Ω f ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x , ( σ ( X e ( x , t ; τ ) , τ ) , χ ) = ( A ( X e ( x , t ; τ ) ) ϕ ( X e ( x , t ; τ ) , τ ) , χ ) . (15)

Lemma 7. [23] Let X : Ω ¯ X ( Ω ) ¯ , X C 2 ( Ω ¯ ) be an invertible vector valued function satifying X ( x ) = x , x Γ . Let F = X and assume that F 1 C 1 ( Ω ¯ ) . Then

Γ F T n ( x ) w ( x ) ψ ( x ) d Γ = Γ n ( x ) w ( x ) ψ ( x ) d e t F 1 d Γ , (16)

with w H 1 ( X ( Ω ) ) and ψ H 1 ( Ω ) , where n is the outward unit normal vector to Γ .

Now, replacing in (15) formula (16) with X ( x ) = X e ( x , t ; τ ) , and replacing the Robin condition, we have

{ Ω d ϕ d τ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x + Ω F e 1 ( x , t ; τ ) σ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x + Ω d i v F e T ( x , t ; τ ) σ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x + Ω r ( X e ( x , t ; τ ) ) ϕ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x + Γ R α ϕ ( X e ( x , t ; τ ) , τ ) ψ ( x ) d e t F e 1 ( x , t ; τ ) d Γ = Ω f ( X e ( x , t ; τ ) , τ ) ψ ( x ) d x + Γ R g ( X e ( x , t ; τ ) , τ ) ψ ( x ) d e t F e 1 ( x , t ; τ ) d Γ , ( σ ( X e ( x , t ; τ ) , τ ) , χ ) = ( A ( X e ( x , t ; τ ) ) ϕ ( X e ( x , t ; τ ) , τ ) , χ ) . (17)

2.3. The Combined Approximate Scheme

We now present our time-stepping procedure for (17). Let N be the number of time steps, Δ t = T / N be the time step, and t n = n Δ t for n = 0 , 1 / 2 , 1 , 3 / 2 , , N . We will use the notation ψ n ( x ) : = ψ ( x , t n ) for a function. Moreover, for n = 0 , 1 , , we define

X e n ( x ) : = X e ( x , t n + 1 ; t n ) , F e n ( x ) : = X e ( x , t n + 1 ; t n ) , X e n + 1 2 ( x ) : = X e ( x , t n + 1 ; t n + 1 2 ) , F e n + 1 2 ( x ) : = X e ( x , t n + 1 ; t n + 1 2 ) . (18)

By fixing t = t n + 1 , n = 0 , 1 , , N 1 in (17) and using a Crank-Nicolson method [26] with respect to τ . Thus, from (18) we have

{ Ω ϕ n + 1 ( x ) ϕ n ( X e n ( x ) ) Δ t ψ ( x ) d x + 1 2 Ω σ n + 1 ( x ) ψ ( x ) d x + 1 2 Ω ( F e n ) 1 ( x ) σ n ( X e n ( x ) ) ψ ( x ) d x + 1 2 Ω d i v ( F e n ) T ( x ) σ n ( X e n ( x ) ) ψ ( x ) d x + 1 2 Ω r ( x ) ϕ n + 1 ( x ) ψ ( x ) d x + 1 2 Ω r ( X e n ( x ) ) ϕ n ( X e n ( x ) ) ψ ( x ) d x + 1 2 Γ R α ( ϕ n + 1 ( x ) + ϕ n ( x ) d e t ( F e n ) 1 ( x ) ) ψ ( x ) d Γ = 1 2 Ω ( f n + 1 ( x ) + f n ( X e n ( x ) ) ) ψ ( x ) d x + 1 2 Γ R ( g n + 1 ( x ) + g n ( x ) d e t ( F e n ) 1 ( x ) ) ψ ( x ) d Γ , ( A 1 ( x ) σ n + 1 ( x ) , χ ) = ( ϕ n + 1 ( x ) , χ ) . (19)

By using (8) and (11), we see that

{ ( F e n ) 1 ( x ) = I ( x ) + Δ t L n ( X e n ( x ) ) + O ( ( Δ t ) 2 ) , d e t ( F e n ) 1 ( x ) = 1 + Δ t d i v v n + 1 ( x ) + O ( ( Δ t ) 2 ) , d i v ( F e n ) T ( x ) = Δ t d i v v n ( X e n ( x ) ) + O ( ( Δ t ) 2 ) . (20)

Taking (20) into (19), we can obtain

{ Ω ϕ n + 1 ϕ n ( X e n ( x ) ) Δ t ψ d x + Ω σ n + 1 + σ n ( X e n ( x ) ) 2 ψ d x + Δ t 2 Ω L n ( X e n ) σ n ( X e n ( x ) ) ψ d x + Δ t 2 Ω d i v v n ( X e n ( x ) ) σ n ( X e n ( x ) ) ψ d x + Ω r ϕ n + 1 + r ( X e n ( x ) ) ϕ n ( X e n ( x ) ) 2 ψ d x + Γ R α ϕ n + 1 + ϕ n ( 1 + Δ t d i v v n + 1 ) 2 ψ d Γ = Ω f n + 1 + f n ( X e n ( x ) ) 2 ψ d x + Γ R g n + 1 + g n ( 1 + Δ t d i v v n + 1 ) 2 ψ d Γ , ( A 1 σ n + 1 , χ ) = ( ϕ n + 1 , χ ) . (21)

We propose two explicit numerical schemes to approximate X e n ( x ) :

X E n ( x ) : = x Δ t v n + 1 ( x ) ( Euler scheme ) , X R K n ( x ) : = x Δ t v n + 1 2 ( x Δ t 2 v n + 1 ( x ) ) ( Runge-Kutta scheme ) . (22)

A similar notation to the one in Section 2.2 is used for the Jacobian of X E n , namely,

F E n ( x ) : = X E n ( x ) = I ( x ) Δ t L n + 1 ( x ) .

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 v C 0 ( W 1, ( Ω ) ) and satisfies v = 0 on Γ .

Lemma 8. [23] Under Claim 1, if v C 0 ( W 1, ( Ω ) ) Δ t < 1 / 2 , we can see that

X E n ( Ω ¯ ) = X R K n ( Ω ¯ ) = Ω ¯ . (23)

Lemma 9. [23] Under Claim 1, if v C 0 ( W 1, ( Ω ) ) Δ t < 1 / 2 , we have

( F E n ) 1 ( x ) = I + Δ t L n + 1 ( x ) + ( Δ t ) 2 ( L n + 1 ( x ) ) 2 + O ( ( Δ t ) 3 ) . (24)

Corollary 1. Under the assumptions of Lemma 9, x Ω , we have

{ d e t ( F E n ) 1 ( x ) = 1 + Δ t d i v v n + 1 ( x ) + O ( ( Δ t ) 2 ) , | d e t ( F E n ) 1 ( x ) | 1 + Δ t v C 0 ( W 1, ( Ω ) ) + O ( ( Δ t ) 2 ) . (25)

Lemma 10. [22] Under Claim 1, if ψ L 2 ( Ω ) and v C 0 ( W 1, ( Ω ) ) Δ t < 1 / 2 , then there exists a positive constant c such that

ψ X i n 2 ( 1 + c Δ t ) ψ 2 , for n = 0, , N and i = E , R K , (26)

where ψ X i n = ψ ( X i n ) .

Thus, in the case where the characteristic lines and their gradients are not explicitly known, we propose the following time approximation of (21)

{ Ω ϕ n + 1 ϕ n X R K n Δ t ψ d x + Ω σ n + 1 + σ n X E n 2 ψ d x + Δ t 2 Ω ( L n σ n ) X E n ψ d x + Δ t 2 Ω ( d i v v n σ n ) X E n ψ d x + Ω r ϕ n + 1 + ( r ϕ n ) X E n 2 ψ d x + Γ R α ϕ n + 1 + ϕ n ( 1 + Δ t d i v v n + 1 ) 2 ψ d Γ = Ω f n + 1 + f n X E n 2 ψ d x + Γ R g n + 1 + g n ( 1 + Δ t d i v v n + 1 ) 2 ψ d Γ , ( A 1 σ n + 1 , χ ) = ( ϕ n + 1 , χ ) . (27)

The time difference approximation (27) will be combined with a standard Galerkin finite element and mixed finite element in the space for ϕ ( x ) and σ ( x ) , respectively [27] [28] . We discrete ϕ ( x ) in space on a quasi-uniform finite element mesh T h ϕ of Ω with maximal element diameter h ϕ . For σ ( x ) , we denote as h σ > 0 and T h σ similarly. Let V h ϕ k × W h σ l H Γ D 1 ( Ω ) × H ( d i v ; Ω ) be finite element spaces with index k and l, respectively.

We define a bilinear form A h n + 1 / 2 on V h ϕ k × W h σ l × V h ϕ k and a linear form F h n + 1 / 2 on V h ϕ k for n = 0 , , N 1 by

{ A h n + 1 / 2 ( ϕ h , σ h , ψ h ) ( ϕ h n + 1 ϕ h n X R K n Δ t , ψ h ) + ( σ h n + 1 + σ h n X E n 2 , ψ h ) + Δ t 2 ( ( L n σ h n ) X E n , ψ h ) + Δ t 2 ( ( d i v v n σ h n ) X E n , ψ h ) + ( r ϕ h n + 1 + ( r ϕ h n ) X E n 2 , ψ h ) + α ( ϕ h n + 1 + ϕ h n ( 1 + Δ t d i v v n + 1 ) 2 , ψ h ) Γ R , F h n + 1 / 2 ( ψ h ) ( f n + 1 + f n X E n 2 , ψ h ) + ( g n + 1 + g n ( 1 + Δ t d i v v n + 1 ) 2 , ψ h ) Γ R , (28)

where ϕ h V h ϕ k , σ h W h σ l , ψ h V h ϕ k .

Then, the fully discrete scheme reads: Given ϕ h 0 V h k , find { ϕ h n , σ h n } n = 1 N V h ϕ k × W h σ l such that

{ A h n + 1 / 2 ( ϕ h , σ h , ψ h ) = F h n + 1 / 2 ( ψ h ) , ψ h V h ϕ k , ( A 1 σ h n + 1 , χ h ) = ( ϕ h n + 1 , χ h ) , χ h W h σ l . (29)

Throughout the analysis, K will denote a generic positive constant, independent of h ϕ , h σ , Δ t . 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 v C 0 ( W 2 , ( Ω ) ) and satisfies v = 0 on Γ .

Remark. Throughout this paper c 1 denotes the maximum between the positive constant appearing in Lemma 10 and the norm of the velocity in C 0 ( W 2, ( Ω ) ) .

Claim 3. The diffusion matrix coefficients, A i j , belong to W 1, ( Ω ) . More- over, A is a positive definite symmetric d × d matrix and there exists a strictly positive constant δ which is a uniform lower bound for the eigenvalues of A 1 .

As a consequence of Claim 3, there exists a unique positive definite symmetric d × d matrix function C , such that A 1 = C 2 and C i j W 1, ( Ω ) . Let us in-

troduce the constant c 2 : = m a x i , j { C i j W 1, ( Ω ) 2 } . Clearly, Claim 3, we have

δ w 0 2 ( A 1 w , w ) = C w 0 2 c 2 w 0 2 , w R d . (30)

Claim 4. The reaction function, r W 1, ( Ω ) , satisfies 0 < γ r ( x ) in Ω , where γ is a constant.

Under the previous claims, let c 3 : = r W 1 , ( Ω ) 2 .

Claim 5. The source function f C 0 ( L 2 ( Ω ) ) . In Robin boundary condition, g C 0 ( L 2 ( Γ R ) ) and α > 0 .

For a given series of functions { ψ } n = 0 N , we define the following norms

ψ l ( L 2 ) = m a x 0 n N { ψ n } , ψ l 2 ( L 2 ) = { Δ t n = 0 N ψ n 2 } 1 / 2 , t ψ n = ( ψ n + 1 ψ n ) / Δ t , ψ A 1 2 = ( A 1 ψ , ψ ) .

Similar definitions are considered for functional spaces l ( L 2 ( Γ R ) ) and l 2 ( L 2 ( Γ R ) ) associated to the Robin boundary condition. Moreover, we define

t ψ l 2 ( L 2 ( Γ R ) ) = { Δ t n = 0 N 1 t ψ n 0 , Γ R 2 } 1 / 2 .

Lemma 11. Let the above Claims 2 - 5 and c 1 Δ t 1 / 2 be assumed. If { ϕ h n , σ h n } n = 1 N be the solution of (29). Then it holds that

A h n + 1 / 2 ( ϕ h , σ h , ϕ h n + 1 ) t ( 1 2 ϕ h n 2 + δ Δ t 4 σ h n 2 + Δ t 4 r ϕ h n 2 + α Δ t 4 ϕ h n 0 , Γ R 2 ) + 1 2 Δ t ϕ h n + 1 ϕ h n X R K n ϕ 2 + δ 4 σ h n + 1 + σ h n X E n 2 + 1 4 r ϕ h n + 1 + ( r ϕ h n ) X E n 2 + α 4 ϕ h n + 1 + ϕ h n ( 1 + Δ t d i v v n + 1 ) 0 , Γ R 2 { c 2 ( ϕ h n 2 + ϕ h n + 1 2 + σ h n 2 ) + c Δ t σ h n 2 + c Δ t ϕ h n + 1 2 + c Δ t ( r ϕ h n 2 + r ϕ h n + 1 2 + α ϕ h n 0 , Γ R 2 ) } , (31)

the constant c is given by c = max { 1 , c 1 , δ 4 2 c 2 , c 1 c 3 / γ } .

Proof. Substituting ψ h = ϕ h n + 1 into (28), we have

A h n + 1 / 2 ( ϕ h , σ h , ϕ h n + 1 ) = ( ϕ h n + 1 ϕ h n X R K n Δ t , ϕ h n + 1 ) + ( σ h n + 1 + σ h n X E n 2 , A 1 σ h n + 1 ) + Δ t 2 ( ( L n σ h n ) X E n , ϕ h n + 1 ) + Δ t 2 ( ( d i v v n σ h n ) X E n , ϕ h n + 1 ) + ( r ϕ h n + 1 + ( r ϕ h n ) X E n 2 , ϕ h n + 1 ) + α ( ϕ h n + 1 + ϕ h n ( 1 + Δ t d i v v n + 1 ) 2 , ϕ h n + 1 ) Γ R = I 1 + I 2 + I 3 + I 4 + I 5 + I 6 . (32)

Lemma 4.1 in [23] implies that

I 1 = 1 2 Δ t ( ϕ h n + 1 2 ϕ h n X R K n 2 ) + 1 2 Δ t ϕ h n + 1 ϕ h n X R K n 2 t ( 1 2 ϕ h n 2 ) c 1 2 ϕ h n 2 + 1 2 Δ t ϕ h n + 1 ϕ h n X R K n 2 , (33)

and

I 2 = 1 4 ( σ h n + 1 A 1 2 σ h n X E n A 1 2 ) + 1 4 σ h n + 1 + σ h n X E n A 1 2 1 4 ( δ σ h n + 1 2 ( 1 + c 1 Δ t ) c 2 σ h n 2 ) + 1 4 σ h n + 1 + σ h n X E n A 1 2 t ( δ Δ t 4 σ h n 2 ) + ( δ 4 2 c 2 ) σ h n 2 + 1 4 σ h n + 1 + σ h n X E n A 1 2 . (34)

Next, by using c 1 Δ t < 1 , we obtain

I 3 Δ t 2 1 + c 1 Δ t L n σ h n ϕ h n + 1 c 1 Δ t 2 { σ h n 2 + ϕ h n + 1 2 } .

Then when I 3 0 and I 3 < 0 , we have

I 3 c 1 Δ t 2 { σ h n 2 + ϕ h n + 1 2 } . (35)

Similarly, for I 4 we obtain the estimate

I 4 c 1 Δ t 2 σ h n 2 1 4 ϕ h n + 1 2 . (36)

Analogous computations to term I 2 give

I 5 t ( Δ t 4 r ϕ h n 2 ) + 1 4 r ϕ h n + 1 + ( r ϕ h n ) X E n 2 max { c 1 , c 1 c 3 / γ } Δ t 2 ( r ϕ h n 2 + r ϕ h n + 1 2 ) . (37)

For the boundary integral term I 6 , we first use some properties of the inner product in the space L 2 ( Γ R ) and the inequality ( 1 + c 1 Δ t ) 2 1 + 3 c 1 Δ t to get the estimate

ψ ( 1 + Δ t d i v v n + 1 ) 0, Γ R 2 ( 1 + c 1 Δ t ) 2 ψ 0, Γ R 2 ( 1 + 3 c 1 Δ t ) ψ 0, Γ R 2

for ψ L 2 ( Γ R ) . Thus, we obtain

I 6 t ( α Δ t 4 ϕ h n 0, Γ R 2 ) + α 4 ϕ h n + 1 + ϕ h n ( 1 + Δ t d i v v n + 1 ) 0, Γ R 2 3 4 c 1 α Δ t ϕ h n 0, Γ R 2 . (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 { ϕ h n , σ h n } n = 1 N be the solution of (29) subject to the initial value ϕ h 0 . Then there exist two positive constants c and d = d ( c 1 , c 2 , δ , c 3 , γ ) , such that, for Δ t < d , we have

1 2 ϕ h l ( L 2 ) + δ Δ t 4 σ h l ( L 2 ) + Δ t 4 r ϕ h l ( L 2 ) + α Δ t 16 ϕ h l ( L 2 ( Γ R ) ) c { 1 2 ϕ h 0 + δ Δ t 4 σ h 0 + Δ t 4 r ϕ h 0 + α Δ t 16 ϕ h 0 0 , Γ R + f l 2 ( L 2 ) + g l 2 ( L 2 ( Γ R ) ) + Δ t t g l 2 ( L 2 ( Γ R ) ) } . (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] ) Π h : C 0 ( Ω ¯ ) V h ϕ k , and P h : C 0 ( Ω ¯ ) W h c l .

Lemma 13. There exist positive constants K 1 and K 2 , independent of h ϕ and h σ respectively, such that

Π h ψ ψ s K 1 h ϕ k + 1 s ψ k + 1 , s = 0 , 1 , ψ H k + 1 ( Ω ) C 0 ( Ω ¯ ) , (40)

and

P h σ σ m K 2 h σ l + 1 m σ l + 1 , m = 0 , 1 , σ H l + 1 ( Ω ) C 0 ( Ω ¯ ) . (41)

Let e h n = ϕ h n Π h ϕ n , η h n = ϕ n Π h ϕ n , θ h n = σ h n P h σ n , ρ h n = σ n P h σ n .

Corresponding to A h n + 1 / 2 and F h n + 1 / 2 , we introduce a bilinear form A n + 1 / 2 on H Γ D 1 ( Ω ) × H ( d i v ; Ω ) × H Γ D 1 ( Ω ) and a linear form F n + 1 / 2 on H Γ D 1 ( Ω ) for n = 0 , , N 1 as follows

{ A n + 1 / 2 ( ϕ , σ , ψ ) ( ( d ϕ d t ) n + 1 2 X e n + 1 2 , ψ ) + ( ( F e n + 1 2 ) 1 σ n + 1 2 X e n + 1 2 , ψ ) + ( d i v ( F e n + 1 2 ) T σ n + 1 2 X e n + 1 2 , ψ ) + ( ( r ϕ n + 1 2 ) X e n + 1 2 , ψ ) + α ( d e t ( F e n + 1 2 ) 1 ϕ n + 1 2 , ψ ) Γ R , F n + 1 / 2 ( ψ ) ( f n + 1 2 X e n + 1 2 , ψ ) + ( d e t ( F e n + 1 2 ) 1 g n + 1 2 , ψ ) Γ R . (42)

If ϕ and σ = A ϕ are the solutions of (1), we have for n = 0 , , N 1

A n + 1 / 2 ( ϕ , σ , ψ ) = F n + 1 / 2 ( ψ ) , ψ H Γ D 1 ( Ω ) . (43)

We decompose e h , θ h as

A h n + 1 / 2 ( e h , θ h , ψ h ) = ( F h n + 1 / 2 F n + 1 / 2 ) ( ψ h ) + A h n + 1 / 2 ( η h , ρ h , ψ h ) + ( A n + 1 / 2 A h n + 1 / 2 ) ( ϕ , σ , ψ h ) , (44)

where ψ h V h ϕ k .

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 v C 0 ( W 3, ( Ω ) ) C 1 ( W 2, ( Ω ) ) C 2 ( L ( Ω ) ) ,

v | Γ = 0, A W 3, ( Ω ) , r W 2, ( Ω ) and that v C 0 ( W 1, ( Ω ) ) Δ t < 1 / 2 . Let the solution of (43) satisfy ϕ Z 3 , ϕ Z 3 , ϕ | Γ R C 2 ( L 2 ( Γ R ) ) . Then, for each n = 0 , 1 , , N 1 , there exist two functions ξ 1 n + 1 2 : Ω R and ξ 2 n + 1 2 : Γ R R ,

such that

( A n + 1 / 2 A h n + 1 / 2 ) ( ϕ , σ , ψ ) = ( ξ 1 n + 1 2 , ψ ) + ( ξ 2 n + 1 2 , ψ ) Γ R , (45)

ψ H Γ D 1 ( Ω ) . Moreover, ξ 1 n + 1 2 L 2 ( Ω ) , ξ 2 n + 1 2 L 2 ( Γ R ) and the following estimates hold:

{ ξ 1 n + 1 2 c ˜ 1 ( Δ t ) 2 ( ϕ Z 3 + σ Z 3 + r ϕ Z 2 ) , ξ 2 n + 1 2 0, Γ R c ˜ 1 ( Δ t ) 2 ( σ n Z 2 , Γ R + α ϕ C 2 ( L 2 ( Γ R ) ) ) , (46)

where c ˜ 1 denotes a constant independent of Δ t .

Lemma 15. [23] Assume the above Claims 2 - 5 hold, and that the coefficients of the problem (1) satisfy v C 0 ( W 2, ( Ω ) ) C 1 ( W 1, ( Ω ) ) , f Z 2 , g C 2 ( L 2 ( Γ R ) ) and v C 0 ( W 1, ( Ω ) ) Δ t < 1 / 2 . Then, for each n = 0 , 1 , , N 1 , there

exist function ξ f n + 1 2 : Ω R and ξ g n + 1 2 : Γ R R , such that

( F n + 1 / 2 F h n + 1 / 2 ) ( ψ ) = ( ξ f n + 1 2 , ψ ) + ( ξ g n + 1 2 , ψ ) Γ R , ψ H 1 ( Ω ) . (47)

Moreover, ξ f n + 1 2 L 2 ( Ω ) , ξ g n + 1 2 L 2 ( Γ R ) and the following estimates hold:

{ ξ f n + 1 2 c ˜ 1 ( Δ t ) 2 f Z 2 , ξ g n + 1 2 0, Γ R c ˜ 1 ( Δ t ) 2 g C 2 ( L 2 ( Γ R ) ) , (48)

where c ˜ 1 denotes a constant independent of Δ t .

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 ϕ C 0 ( C 0 ( Ω ¯ ) ) C 0 ( H k + 1 ( Ω ) ) H 1 ( H k ( Ω ) ) and c 1 Δ t < 1 / 2 . There exists

A h n + 1 / 2 ( η h , ρ h , e h n + 1 ) 1 8 C θ h n + 1 + ( C θ h n ) X E n 2 + t ( Δ t 2 ( C ρ h n , C θ h n ) ) + 1 8 r e h n + 1 + ( r e h n ) X E n 2 + t ( Δ t 2 ( r η h n , r e h n ) ) + α 8 e h n + 1 + e h n ( 1 + Δ t d i v v n + 1 ) 0, Γ R 2 + t ( α Δ t 2 ( η h n , e h n ) Γ R ) + c e h n + 1 2 + c Δ t ( θ n 2 + θ n + 1 2 + r e h n 2 + r e h n + 1 2 + α e h n 0, Γ R 2 ) + c ˜ K 1 2 h ϕ 2 k ( 1 Δ t ϕ t L 2 ( t n , t n + 1 ; H k ) 2 + 1 Δ t ϕ L 2 ( t n , t n + 1 ; H k + 1 ) 2 + Δ t ϕ n k + 1 2 ) + c ˜ K 2 2 h σ 2 l σ n l + 1 2 + c ˜ Δ t e h n + 1 2 , (49)

with c = max { 1 , ( c 1 c 2 + 1 ) / 4 δ , c 1 c 3 / 4 , c 1 } , c ˜ a positive constant.

Proof. From the definition of A h n + 1 / 2 , we have

A h n + 1 / 2 ( η h , ρ h , e h n + 1 ) = ( η h n + 1 η h n X R K n Δ t , e h n + 1 ) + ( ρ h n + 1 + ρ h n X E n 2 , A 1 θ h n + 1 ) + Δ t 2 ( ( L n ρ h n ) X E n , e h n + 1 ) + Δ t 2 ( ( d i v v n ρ h n ) X E n , e h n + 1 ) + ( r η h n + 1 + ( r η h n ) X E n 2 , e h n + 1 ) + α ( η h n + 1 + η h n ( 1 + Δ t d i v v n + 1 ) 2 , e h n + 1 ) Γ R = I 1 + I 2 + I 3 + I 4 + I 5 + I 6 . (50)

From the arguments of Lemma 4.1 in [24] , we get similarly the bounds of the above terms I 1 , I 5 , I 6 as the followings:

I 1 ( 1 + c 1 Δ t ) ( 1 + 2 c 1 ) 2 K 1 2 h ϕ 2 k 2 Δ t ( ϕ t L 2 ( t n , t n + 1 ; H k ) 2 + ϕ L 2 ( t n , t n + 1 ; H k + 1 ) 2 ) + 1 2 e h n + 1 2 , (51)

I 5 t ( Δ t 2 ( r η h n , r e h n ) ) + ( c 1 c 3 Δ t 4 + ( 1 + c 1 Δ t ) c 3 ) K 1 2 h ϕ 2 k ϕ n k + 1 2 + ( 2 c 1 + c 1 2 Δ t ) c 3 Δ t 8 ( r e h n 2 + r e h n + 1 2 ) + 1 8 r e h n + 1 + ( r e h n ) X E n 2 , (52)

I 6 t ( α Δ t 2 ( η h n , e h n ) Γ R ) + α 8 e h n + 1 + e h n ( 1 + Δ t d i v v n + 1 ) 0, Γ R 2 + α ( c 1 Δ t ( 2 + c 1 Δ t ) 4 + 1 + 2 c 1 Δ t + c 1 2 ( Δ t ) 2 2 ) c Ω K 1 2 h ϕ 2 k ϕ n k + 1 2 + α c 1 Δ t ( 2 + c 1 Δ t ) 4 e h n 0, Γ R 2 . (53)

To estimate I 2 , we divide it into three parts

I 2 = 1 2 [ ( C ρ h n + 1 , C θ h n + 1 ) ( C ρ h n , C θ h n ) ] + I 21 + I 22 , (54)

where

I 21 = 1 2 [ ( C ρ h n , C θ h n ) 1 2 ( ( C ρ h n ) X E n , ( C θ h n ) X E n ) ] , I 22 = 1 2 [ ( ( C ρ h n ) X E n , C X E n ) ( θ h n + 1 + θ h n X E n ) ] .

Using the transformation y = X E n ( x ) , we have

I 21 = 1 2 Ω ( 1 d e t ( ( F E n ) 1 ( x ) ) ) ( C ρ h n ) ( x ) ( C θ h n ) ( x ) d x c 1 c 2 Δ t 4 ( K 2 2 h σ 2 l σ n l + 1 2 + θ h n 2 ) . (55)

Now, we replace in I 22 equality C ( X E n ( x ) ) = C ( x ) D n ( x ) , where

D i j n ( x ) : = t n t n + 1 C i j ( Y E n ( x , s ) ) v n + 1 ( x ) d s , a . e . x Ω , (56)

with | D i j n ( x ) | c 1 c 2 Δ t , and the function Y E n ( x , ) : [ t n , t n + 1 ] Ω ¯ is defined by Y E n ( x , s ) : = x ( t n + 1 s ) v n + 1 ( x ) . Thus, we get

I 22 = 1 2 ( ( C ρ h n ) X E n , C θ h n + 1 D θ h n + 1 + ( C θ h n ) X E n ) ( C ρ h n ) X E n 2 + 1 8 C θ h n + 1 + ( C θ h n ) X E n 2 + 1 8 D θ h n + 1 2 . (57)

Moreover, we have

{ ( C ρ h n ) X E n 2 ( 1 + c 1 Δ t ) c 2 K 2 2 h σ 2 l σ n l + 1 2 , D θ h n + 1 2 c 1 2 c 2 ( Δ t ) 2 θ h n + 1 2 . (58)

By considering together from (54) to (58), we can state

I 2 t ( Δ t 2 ( C ρ h n , C θ h n ) ) + ( c 1 c 2 Δ t 4 + ( 1 + c 1 Δ t ) c 2 ) K 2 2 h σ 2 l σ n l + 1 2 + ( 2 c 1 + c 1 2 Δ t ) c 2 Δ t 8 δ δ ( θ h n + 1 2 + θ h n 2 ) + 1 8 C θ h n + 1 + ( C θ h n ) X E n 2 . (59)

Similar arguments, i.e. Lemma 5.4 in [23] lead to

{ ( L n ρ h n ) X E n 2 ( 1 + c 1 Δ t ) c 1 2 c 2 2 K 2 2 h σ 2 l σ n l + 1 2 , ( d i v v n ρ h n ) X E n 2 ( 1 + c 1 Δ t ) c 1 2 c 2 2 K 2 2 h σ 2 l σ n l + 1 2 . (60)

Using these inequalities, I 3 and I 4 can be bounded as follows

{ I 3 ( 1 + c 1 Δ t ) c 1 2 c 2 2 K 2 2 h σ 2 l 4 σ n l + 1 2 + Δ t 4 δ δ θ h n + 1 2 , I 4 ( 1 + c 1 Δ t ) c 1 2 c 2 2 K 2 2 h σ 2 l 4 σ n l + 1 2 + 1 4 θ h n + 1 2 . (61)

Gathering (50), (51), (52), (53), (59), and (61) together, we complete the proof.

We now turn to estimate e h n + 1 and σ h n + 1 . From the definition of A h n + 1 / 2 and A n + 1 / 2 , we see

A h n + 1 / 2 ( e h , σ h , e h n + 1 ) = A h n + 1 / 2 ( η h , ρ h , e h n + 1 ) + ( F h n + 1 / 2 F n + 1 / 2 ) ( e h n + 1 ) + ( A n + 1 / 2 A h n + 1 / 2 ) ( ϕ , σ , e h n + 1 ) . (62)

From Lemma 11, we obtain the lower bound for (62)

A h n + 1 / 2 ( e h , σ h , e h n + 1 ) t ( 1 2 e h n 2 + δ Δ t 4 θ h n 2 + Δ t 4 r e h n 2 + α Δ t 4 e h n 0, Γ R 2 ) + 1 2 Δ t e h n + 1 e h n X R K n e 2 + δ 4 θ h n + 1 + θ h n X E n 2 + 1 4 r e h n + 1 + ( r e h n ) X E n 2 + α 4 e h n + 1 + e h n ( 1 + Δ t d i v v n + 1 ) 0, Γ R 2 { c 2 ( e h n 2 + e h n + 1 2 + θ h n 2 ) + c Δ t θ h n 2 + c Δ t e h n + 1 2 + c Δ t ( r e h n 2 + r e h n + 1 2 + α e h n 0, Γ R 2 ) } . (63)

By jointly considering the lower bound (63), the upper bounds given in Lemmas 14 - 16, we deduce

t ( 1 2 e h n 2 + δ Δ t 4 θ h n 2 + Δ t 4 r e h n 2 + α Δ t 4 e h n 0, Γ R 2 ) { c 2 ( e h n 2 + e h n + 1 2 + θ h n 2 ) + c Δ t θ h n 2 + c Δ t e h n + 1 2 + c Δ t ( r e h n 2 + r e h n + 1 2 + α e h n 0, Γ R 2 ) } + t ( Δ t 2 ( C ρ h n , C θ h n ) ) + t ( Δ t 2 ( r η h n , r e h n ) ) + t ( α Δ t 2 ( η h n , e h n ) Γ R ) + c Δ t ( θ n 2 + θ n + 1 2 + r e h n 2 + r e h n + 1 2 + α e h n 0, Γ R 2 ) + c ˜ K 1 2 h ϕ 2 k ( 1 Δ t ϕ t L 2 ( t n , t n + 1 ; H k ) 2 + 1 Δ t ϕ L 2 ( t n , t n + 1 ; H k + 1 ) 2 + Δ t ϕ n k + 1 2 ) + c ˜ K 2 2 h σ 2 l σ n l + 1 2 + c ˜ Δ t e h n + 1 2 + c e h n + 1 2 + c ˜ ( ξ 1 n + 1 2 2 + ξ 2 n + 1 2 0, Γ R 2 + ξ f n + 1 2 2 + ξ g n + 1 2 0, Γ R 2 ) . (64)

Analogous computations to those developed in Lemma 16 give

{ Δ t 2 ( C ρ h j , C θ h j ) c 2 K 2 2 h σ 2 l Δ t 2 σ j l + 1 2 + Δ t 8 C θ h j 2 , Δ t 2 ( r η h j , r e h j ) c 3 K 1 2 h ϕ 2 k Δ t 2 ϕ j k + 1 2 + Δ t 8 r e h j 2 , α Δ t 2 ( η h j , e h j ) Γ R α c Ω K 1 2 h ϕ 2 k Δ t ϕ j k + 1 2 + α Δ t 16 e h j 0 , Γ R 2 , j = 0 , N . (65)

Putting (65) into (64), multiplying (64) by Δ t , summing it about time form t = 0 and t = t N , taking the initial values ϕ h 0 = Π h ϕ 0 , σ h 0 = P h σ 0 and using discrete Gronwall’s inequality, we can derive the following error estimate:

Theorem 17 (Error Estimate) Let Claims 2 - 5 be assumed. Let ϕ Z 3 C 0 ( H k + 1 ( Ω ) ) H 1 ( H k ( Ω ) ) be the solution of (1) with ϕ Z 3 ,

ϕ | Γ R Z 3 ( Γ R ) , { ϕ h n , σ h n } n = 1 N be the solution of (29) subject to the initial values

ϕ h 0 = Π h ϕ 0 , σ h 0 = P h σ 0 . Then there exists a positive constant c independent of h ϕ , h σ and Δ t such that

1 2 ϕ ϕ h l ( L 2 ) + δ Δ t 8 σ σ h l ( L 2 ) + Δ t 8 r ϕ r ϕ h l ( L 2 ) + α Δ t 16 ϕ ϕ h l ( L 2 ( Γ R ) ) c h ϕ k ( ϕ H 1 ( H k ) + ϕ C 0 ( H k + 1 ) ) + c h σ l σ C 0 ( H l + 1 ) + c ( Δ t ) 2 ( ϕ Z 3 + σ Z 3 + r ϕ Z 2 + ϕ Z 2 , Γ R + f Z 2 + g Z 2 , Γ R ) + c ( Δ t ) 5 2 ( ϕ Z 3 , Γ R + g Z 3 , Γ R ) . (66)

Acknowledgement

This research was supported by the NSF of China (Grant Nos. 11271231 and 11301300).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1] Frolkovic, P. and De Schepper, H. (2000) Numerical Modelling of Convection Dominated Transport Coupled with Density Driven Flow in Porous Media. Advances in Water Resources, 24, 63-72.
https://doi.org/10.1016/S0309-1708(00)00025-7
[2] Markowich, P.A. and Szmolyan, P. (1989) A System of Convection-Diffusion Equations with Small Diffusion Coefficient Arising in Semiconductor Physics. Journal of Differential Equations, 81, 234-254.
https://doi.org/10.1016/0022-0396(89)90122-8
[3] Morton, K.W. (1996) Numerical Solution of Convection-Diffusion Problems. Chapman & Hall, London.
[4] Ewing, R.E. and Wang, H. (2001) A Summary of Numerical Methods for Time-Dependent Advection-Fominated Partial Differential Equations. Journal of Computational and Applied Mathematics, 128, 423-445.
https://doi.org/10.1016/S0377-0427(00)00522-7
[5] Wilmott, P., Dewynne, J. and Howison, S. (1993) Option Pricing, Mathematical Models and Computation. Oxford Financial Press, Oxford.
[6] Ewing R. (1983) The Mathematics of Reservoir Simulation. SIAM, Philadelphia.
[7] Johnson, C. (1987) Numerical Solution of Partial Differential Equations by the Finite Element Method. Cambridge University Press, Cambridge.
[8] Douglas, Jr.J. and Russell, T.F. (1982) Numerical Methods for Convection-Dominated Diffusion Problems Based on Combining the Method of Characteristics with Finite Element or Finite Difference Procedures. SIAM Journal on Numerical Analysis, 19, 871-885.
https://doi.org/10.1137/0719063
[9] Russell, T.F. (1985) Time-Stepping along Characteristics with Incomplete Iterationfora Galerkin Approximation of Miscible Displacement in Porous Media. SIAM Journal on Numerical Analysis, 22, 970-1013.
https://doi.org/10.1137/0722059
[10] Pironneau, O. (1982) On the Transport-Diffusion Algorithm and Its Application to the Navier-Stokes Equations. Numerische Mathematik, 38, 309-332.
https://doi.org/10.1007/BF01396435
[11] Süli, E. (1986) Convergence Analysis of the Lagrange-Galerkin Method for the Navier-Stokes Equations, Report 86/3. Oxford University Computing Laboratory, Oxford.
[12] Morton, K.W., Priestley, A. and Süli, E. (1986) Stability Analysis of the Lagrange-Galerkin Method with Non-Exact Integration, Report 86/14. Oxford University Computing Laboratory, Oxford.
[13] Johnson, C. and Thomée, V. (1981) Error Estimates for Some Mixed Finite Element Methods for Parabolic Type Problems. RAIRO Analyse Numérique, 15, 41-78.
https://doi.org/10.1051/m2an/1981150100411
[14] Raviart, P.A. and Thomas, J.M. (1977) A Mixed Finite Element Method for Second Order Elliptic Problems, Mathematical Aspects of the Finite Element Method. Lecture Notes in Mathematics 606, Springer, Berlin.
[15] Nedelec, J.C. (1980) Mixed Finite Elements in R3. Numerische Mathematik, 35, 315-341.
https://doi.org/10.1007/BF01396415
[16] Douglas, Jr.J. and Roberts, J.E. (1985) Global Estimates for Mixed Methods for Second Order Elliptic Equations. Mathematics of Computation, 44, 39-52.
https://doi.org/10.1090/S0025-5718-1985-0771029-9
[17] Ewing, R.E., Russell, T.F. and Wheeler, M.F. (1984) Convergence Analysis of an Approximation of Miscible Displacement in Porous Media by Mixed Finite Elements and a Modified Method of Characteristics. Computer Methods in Applied Mechanics and Engineering, 47, 73-92.
https://doi.org/10.1016/0045-7825(84)90048-3
[18] Wheeler, M.F. and Dawson, C.N. (1984) An Operator-splitting Method for Advection Diffusion-Reaction Problems. The Mathematics of Finite Elements and Applications, V, 122-144.
[19] Arbogast, T. and Wheeler, M.F. (1995) A Charcteristics-Mixed Finite Element Methods for Advection-Dominated Transport Problems. SIAM Journal on Numerical Analysis, 32, 404-424.
https://doi.org/10.1137/0732017
[20] Sun, T.J. and Yuan, Y.R. (2009) An Approximation of Incompressible Miscible Displacement in Porous Media by Mixed Finite Element Method and Characteristics-mixed Finite Element Method. Journal of Computational and Applied Mathematics, 228, 391-411.
https://doi.org/10.1016/j.cam.2008.09.029
[21] Sun, T.J. and Yuan, Y.R. (2015) Mixed Finite Element Method and the Characteristics-mixed Finite Element Method for a Slightly Compressible Miscible Displacement Problem in Porous Media. Mathematics and Computers in Simulation, 107, 24-45. https://doi.org/10.1016/j.matcom.2014.07.005
[22] Rui, H.X. and Tabata, M. (2002) A Second Order Charcteristic Finite Element Scheme for Convection-Diffusion Problems. Numerische Mathematik, 92, 161-177.
https://doi.org/10.1007/s002110100364
[23] Bermúdez, A., Nogueiras, M.R. and Vázquez, C. (2006) Numerical Analysis of Convection-Diffusion-Reaction Problems with Higher Order Characteristic/Finite Elements, Part I: Time Discretization. SIAM Journal on Numerical Analysis, 44, 1829-1853.
https://doi.org/10.1137/040612014
[24] Bermúdez, A., Nogueiras, M.R. and Vázquez, C. (2006) Numerical Analysis of Convection-Diffusion-Reaction Problems with Higher Order Characteristic/Finite Elements, Part II: Fully Discretized Scheme and Quadrature Formulas. SIAM Journal on Numerical Analysis, 44, 1854-1876.
https://doi.org/10.1137/040615109
[25] Brezzi, F. (1974) On the Existence, Uniqueness and Approximation of Saddle-Point Problems Arising from Lagrangian Multipliers. RAIRO Analyse Numérique, 2, 129-151.
https://doi.org/10.1051/m2an/197408R201291
[26] Wheeler, M.F. (1973) A Priori L2 Error Estimates for Galerkin Approximations to Parabolic Partial Differential Equations, SIAM Journal on Numerical Analysis, 10, 723-759.
https://doi.org/10.1137/0710062
[27] Arnolds, D.N., Scott, L.R. and Vogelus, M. (1998) Regular Inversion of the Divergence Operator with Dirichlet Boundary Conditions on a Polygonal. Annali della Scuola Normale Superiore di Pisa-Classe di Scienze, 15, 169-192.
[28] Ciarlet, P.G. (1991) Basic Error Estimates for Ellptic Problems, Handbook of Numerical Analysis, Volume II. North-Holland Publishing Company, Amsterdam, 19-351.
[29] Russell, T.F. (1986) Finite Elements with Characteristics for Two-component Incompressible Miscible Displacement, SPE 10500. Proceedings of 6th SPE Symposium on Reservoir Simulation, New Orleans, 123-135.
[30] Ciarlet, P.G. (1978) The Finite Element Method for Ellptic Problems. North-Holland Publishing Company, Amsterdam.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.