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

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 L-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.


Introduction
Let Ω be a bounded domain in Here, where 0 1 , a a 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 convectiondominated 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 L 2 -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 esti- mates were proved in the framework of L 2 -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 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 L 2 -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.

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 the closed subspace of ( ) H Ω defined by ; : , .

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 ( ) [ ] , 0, x t T ∈ Ω × , the characteristic line through ( ) We adopted some propositions and lemmas from [23].
and it is n C with respect to the x variable.
In order to compute second order approximations of matrices e F and 1 e F − , we need the following equations: , ; , , , 0, .
, then e F satisfies the Taylor expansion and its inverse, By using the Liouville's theorem and the chain rule, we obtain ( )) )

Variational Formulation
From the definition of the characteristic lines and by using the chain rule, it follows e e e e X x t By introducing the flux σ φ = ∇ A and using (12), we rewrite equation (1.1) at point

( )
, ; e X x t τ and time τ as follows Before giving a week formulation of ( 13), we adopted a lemma from [23], which can be considered as a Green's formula.
Ω be an invertible vector valued function.Let F X = ∇ and assume that ( ) a vector valued function and ( ) Now, we can multiply ( 13) by test functions ( ) ( ) integrate in Ω respecetively, and apply the usual Green's formula and ( 14) with and assume that ( ) where n is the outward unit normal vector to Γ .Now, replacing in ( 15) formula ( 16) with

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 n t n t = ∆ for 0,1 2,1,3 2, , n N = .
Taking (20) into (19), we can obtain ( ) ) We propose two explicit numerical schemes to approximate ( ) Euler scheme , : Runge-Kutta scheme . 2 A similar notation to the one in Section 2.2 is used for the Jacobian of .
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.
Corollary 1.Under the assumptions of Lemma 9, x ∀ ∈Ω , we have then there exists a positive constant c such that , for 0, , and , , ) Thus, in the case where the characteristic lines and their gradients are not explicitly known, we propose the following time approximation of ( ) The time difference approximation (27) will be combined with a standard Galerkin finite element and mixed finite element in the space for ( ) We define a bilinear form and a linear form ( ) where , , Then, the fully discrete scheme reads: Given .
Throughout the analysis, K will denote a generic positive constant, independent of , , h h t φ σ ∆ .Similarly, ε will denote a generic small positive constant.

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

A .
As a consequence of Claim 3, there exists a unique positive definite symmetric : Clearly, Claim 3, we have ( ) Claim 4. The reaction function, ( ) in Ω , where γ is a constant.
Similar definitions are considered for functional spaces l L Γ associated to the Robin boundary condition.Moreover, we define Let the above Claims 2 -5 and 1 1 2 be the solution of (29).Then it holds that the constant c is given by 6 .

I I I + +
(32) Lemma 4.1 in [23] implies that ( ) ( ) Then when 3 0 I ≥ and 3 0 I < , we have Similarly, for 4 I we obtain the estimate For the boundary integral term 6 I , we first use some properties of the inner product in the space and the inequality ( ) L ψ ∈ Γ .Thus, we obtain ( ) 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 0 h φ .Then there exist two positive constants c and ( ) , , ,  d d c c  c δ γ = , such that, for t d ∆ < , we have .

Error Estimate Theorem
Now, we turn to derive a priori error estimate in L 2 -norm for the solutions of (29).In order to state error estimates, we need two following Lagrange interpolation operators ([29] [30]) Lemma 13.There exist positive constants 1 K and 2 K , independent of h φ and h σ respectively, such that ( ) ( ) ( ) Corresponding to , we introduce a bilinear form ; Ω and a linear form − as follows ( ) If φ and σ φ = ∇ A are the solutions of (1), we have for 0, , ) 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) Then, for each 0,1, , 1, n N = − there exist two functions ( ) Γ and the following estimates hold: ( where 1 c 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 . Then, for each : , and the following estimates hold: where 1 c denotes a constant independent of t ∆ .Now, we turn to bound the third term on the right-hand side of (44).
From the arguments of Lemma 4.1 in [24], we get similarly the bounds of the above terms 1 5 6 , , I I I as the followings: To estimate 2 I , we divide it into three parts ( ) ( ) where , : .
We now turn to estimate .
= ∅ .Let T be a positive constant.In this paper, we will consider the following linear convection-diffusion-reaction equations: find ( ) -th derivative of ϕ with respect to time.considered for the boundary sets Γ R and Γ D .Denote by( )

,
27][28].We discrete ( )x φin space on a quasi-uniform finite element mesh h φ  of Ω with maximal element diameter h φ .For ( ) element spaces with index k and l, respectively.

3 . 1 ,
this paper 1 c denotes the maximum between the positive constant appearing in Lemma 10 and the norm of the velocity in The diffusion matrix coefficients, ij A , belong to ( ) W ∞ Ω .Moreover, 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 1 − For a given series of functions { } 0 N n ψ = , we define the following norms

3 Z=
s inequality, we can derive the following error estimate: Theorem 17 (Error Estimate) Let Claims 2 -5 be assumed.Let be the solution of(29) subject to the initial values is the convection vector field; ) , . )