`International Journal of Mathematics and Mathematical SciencesVolumeΒ 2012, Article IDΒ 870402, 21 pageshttp://dx.doi.org/10.1155/2012/870402`
Research Article

## A Second Order Characteristic Method for Approximating Incompressible Miscible Displacement in Porous Media

School of Mathematics, Shandong University, Jinan 250100, China

Received 25 March 2012; Revised 5 July 2012; Accepted 28 July 2012

Copyright Β© 2012 Tongjun Sun and Keying Ma. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

An approximation scheme is defined for incompressible miscible displacement in porous media. This scheme is constructed by two methods. Under the regularity assumption for the pressure, cubic Hermite finite element method is used for the pressure equation, which ensures the approximation of the velocity smooth enough. A second order characteristic finite element method is presented to handle the material derivative term of the concentration equation. It is of second order accuracy in time increment, symmetric, and unconditionally stable. The optimal -norm error estimates are derived for the scalar concentration.

#### 1. Introduction

In this paper, we will consider the following incompressible miscible displacement in porous media, which is governed by a coupled system of partial differential equations with initial and boundary values. The pressure is governed by an elliptic equation and the concentration is governed by a convection-diffusion equation [16] as follows:where is a bounded domain in , , and is nonzero at injection wells only. The variables in (1.1a)–(1.1d) are the pressure in the fluid mixture, the Darcy velocity , and the relative concentration of the injected fluid. The is the unit outward normal vector on boundary .

The coefficients and data in (1.1a)–(1.1d) are , the permeability of the porous media; , the viscosity of the fluid mixture; , representing flow rates at wells; , the porosity of the rock; , the injected concentration at injection wells () and the resident concentration at production wells (); . Here, for the diffusion coefficient, we consider a dispersion-free case [5, 6] as follows: where is the molecular diffusivity and is a identity matrix. Furthermore, a compatibility condition must be imposed to determine the pressure.

The pressure equation is elliptic and easily handled, but the concentration equation is parabolic and normally convection dominated. It is well known that standard Galerkin scheme applied to the convection-dominated problems does not work well and produces excessive numerical diffusion or nonphysical oscillation. The characteristic method has been introduced to obtain better approximations for (1.1a)–(1.1d), such as characteristic finite element method [37], characteristic finite difference method [8], the modified of characteristic finite element method (MMOC-Galerkin) [9], and the Eulerian-Lagrangian localized adjoint method (ELLAM) [10].

We had considered a combined numerical approximation for (1.1) in [11]. 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 -norm error estimates were derived.

It should be pointed out that the works mentioned above 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 [12] 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 -theory. Numerical analysis of convection-diffusion-reaction problems with higher order characteristic/finite elements were analyzed in [13, 14], which extended the work [12]. The error estimates of second order in time increment were obtained.

The goal of this paper is to present a second order characteristic finite element method in time increment to handle the material derivative term of the concentration equation of (1.1a)–(1.1d). It is organized as follows. In Section 2, we formulate the second order characteristic finite element method for the concentration and cubic Hermite finite element method for the pressure, respectively. Then, we present a combined approximation scheme. In Section 3, we analyze the stability of the approximation scheme for the concentration equation. In Section 4, we derive the optimal-order -norm error estimates for the scalar concentration. They are of second order accuracy in time increment, symmetric, and unconditionally stable. We conclude our results in Section 5.

#### 2. Formulation of the Method

##### 2.1. Statements and Assumptions

In this paper, we adopt notations and norms of usual Sobolev spaces. For periodic functions, we use the notation as in [4] as follows. Let: be the periodic Sobolev space on with the usual norm. If , we write with norm

Moreover, we adopt some notations for the functional spaces involved, which were used in [1214]. For a Banach space and a positive integer , spaces and are abbreviated as and , respectively, and endowed with norms where denotes the th derivative of with respect to time. The Banach space is defined by equipped with the norm

We also require the following assumptions on the coefficients in (1.1) [3]. Let , and be positive constants such that Other assumptions will be made in individual theorems as necessary.

##### 2.2. A Cubic Hermite Element Method for the Pressure

The variational form for the pressure equation (1.1a) is equal to the following:

For , we discretize (2.7) in space on a quasiuniform triangle mesh of with diameter of element . Let be a cubic Hermite finite element space for this mesh. The finite element method for the pressure, given at a time , consists of such that

Since the left-hand side of (2.7) ((2.8), resp.) is a continuous and coercive bilinear form and the right-hand side of (2.7) ((2.8), resp.) is a continuous linear functional, existence and uniqueness of (, resp.) are ensured obviously [15].

By the theory of the cubic Hermit element, the finite element space possess the following approximation property [15, 16]: where is a positive constant independent of .

Define Π: is the interplant operator. Then, we have the following theorem.

Theorem 2.1. Let and be the solutions of (1.1a)–(1.1d) and (2.8), respectively, and assume . Then, there exists a positive constant independent of such that

Proof. By Céa lemma, we have Then by (2.9), we can derive

##### 2.3. A Second Order Characteristic Method for the Concentration

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 following initial value problem: where .

Next, assuming they exist, we denote by (by , resp.) the gradient of (of , resp.) with respect to the space variable , that is, We adopted some propositions and lemmas from [13].

Proposition 2.2. If for an integer, then and it is with respect to the variable.

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

Proposition 2.3. If , then

Proposition 2.4. If , then satisfies the Taylor expansions as and its inverse, , satisfies the Liouville theorem as

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

Proposition 2.5. If , then

Proposition 2.6. If , then satisfies

Variational Formulation
From the definition of the characteristic curves and by using the chain rule, it follows that By writing (1.1b) at point and time and using (2.22), we have
Before giving a week formulation of (2.23), we adopted a lemma from [13], which can be considered as Green’s formula.

Lemma 2.7. Let , be an invertible vector valued function. Let and assume that . Then with a vector valued function and a scalar function.

Now, we can multiply (2.23) by a test function , integrate in , and apply the usual Green’s formula and (2.24) with , obtaining

##### 2.4. The Combined Approximation Scheme

We now present our sequential time-stepping procedure that combines (2.8) and (2.25). Part into , with . The analysis is valid for variable time steps, but we drop the superscript from for convenience. For functions on , we write for . As in [3], let us part into pressure time steps , with . Each pressure step is also a concentration step, that is, for each there exists such that , in general, . We may vary , but except for we drop the superscript. For functions on , we write for ; thus, subscripts refer to pressure steps and superscripts to concentration steps.

If concentration step relates to pressure steps by , we require a velocity approximation for (2.25) based on and earlier values. If , take the linear extrapolation of and defined by if , set We retain the superscript on because is first-order correct in time during the first pressure step and second-order during later steps.

Let be a solution of the ordinary differential equation Then, we can write Subject to an initial condition , we get approximate values of at by the Euler method and the second order Runga-Kutta method, respectively,

Next, assuming they exist, we denote by (resp., by ) the gradient of (resp., of ) with respect to the space variable , that is, [13, 14]

Hypothesis 2.8. For convenience, we assume that (1.1a)–(1.1d) is -periodic ([3]), that is, all functions will be assumed to be spatially -periodic throughout the rest of this paper.

This is physically reasonable, because no-flow conditions (1.1c) are generally treated by reflection, and, in general, interior flow patterns are much more important than boundary effects in reservoir simulation. Thus, the boundary conditions (1.1c) can be dropped.

Lemma 2.9. Under Hypothesis 2.8, if , we can see that

Lemma 2.10. Under Hypothesis 2.8, if , we have

Corollary 2.11. Under the assumptions of Lemma 2.10, , we have

The time difference for (2.25) will be combined with a standard finite element in the space variables. For , we discrete (2.25) in space on a quasi-uniform mesh of with diameter of element . Let be a finite element space.

A characteristic discretization of the weak form (2.25) is given by such thatwhere and are compositions is an initial approximation of exact solution into , which will be defined in Section 5.

At each pressure time step , we define is the truncation of to . Then at , (2.8) is the following: The steps of calculation are as follows.

Step 1. known solve by (2.37), (2.39a) and (2.39b);

Step 2. by (2.35a) and (2.35b) to solve then by (2.35a) and (2.35b) to solve ;

Step 3. analogously, known such that ;

Step 4. then by (2.37), (2.39a) and (2.39b) for ;

Step 5. calculate the approximations in turn analogously to get the pressure, velocity, and concentration at other time step, respectively.

Throughout the analysis, will denote a generic positive constant, independent of , , , and , but possibly depending on constants in . Similarly, will denote a generic small positive constant.

#### 3. Stability for the Concentration Equation

In this section, we derive the stability for the concentration equation. For a given series of functions , we define the following norms and seminorm: In the following sections, we use positive constants as In our analysis, we need some lemmas.

Lemma 3.1. Under the definitions (2.26) and (2.30), for and , it holds that

Proof. We only need to show a proof in the case . Let be the Jacobian matrix of the transformation as According to the proof of (3.23) in [4], we have Since , for sufficiently small, we see that
Following (3.6), we have We complete the proof.

Suppose that and be given. We define linear forms and on for by where .

Lemma 3.2. Let be a solution of (2.35a) and (2.35b). Then it holds that where is the forward difference operator defined by

Proof. Substituting into (3.8), we have Lemma 3.1 implies that Next, by using , we obtain Then when and , we have Similarly, for , we obtain the estimate Analogous computations to term give which completes the proof.

From Lemma 3.1, we have Combining (3.18) with (3.8), we get which completes the proof of the stability by virtue of Gronwall’s inequality.

Theorem 3.3 3.3 (stability). Let be the solution of (2.35a) and (2.35b) subject to the initial value . Then there exists a positive constant independent of and such that

#### 4. Error Estimate Theorem

Now, we turn to derive an optimal priori error estimate in -norm for the concentration of approximation (2.35a) and (2.35b). In order to state error estimates, we need the following Lagrange interpolation operator [15] .

Lemma 4.1. There exists a positive integer such that Let and . By Lemma 4.1, it holds that

Let , , and be given. Corresponding to and , we introduce linear forms and on for by (2.25) as follows: If is the solution of (1.1a)–(1.1d), we have for

We decompose as In order to estimate the terms of the right-hand side of (4.5), we need the following lemmas.

Lemma 4.2 (see [13]). Assume the above Hypotheses hold, and that the coefficients of the problem (1.1a)–(1.1d) satisfy , , and . Let the solution of (4.4) satisfy , . Then, for each , there exists function , such that Moreover, and the following estimate holds: where denotes a constant independent of .

Lemma 4.3 (see [13]). Assume the above Hypotheses hold, and that the coefficients of the problem (1.1a)–(1.1d) satisfy , , and . Let the solution of (4.4) satisfy , . Then, for each , there exists function , such that Moreover, and the following estimate holds: where denotes a constant independent of .

Lemma 4.4. Let . There exists

Proof. From the definition of , we have Since it holds that we have To estimate , we divide it into three parts as where The term is estimated as Using the transformation , we have It is clear that
Term can be divided it into three parts like term as where Thus, the same kind of computations used for leads to the following estimate: Gathering (4.11)–(4.19) together, we complete the proof.

We now turn to estimate . From the definition of and , we see From Lemmas 4.2, 4.3 and 4.4, we obtain Summing up the above equation about time from to , we get where Using the estimates for and , and Gronwall’s inequality, we derive the following error estimate.

Theorem 4.5 (error estimate). Let be the solution of (2.35a) and (2.35b) subject to the initial value . Then there exists a positive constant independent of and such that