International Journal of Mathematics and Mathematical Sciences

Volume 2012 (2012), Article ID 870402, 21 pages

doi:10.1155/2012/870402

Research Article

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

Tongjun Sun and Keying Ma

School of Mathematics, Shandong University, Jinan 250100, China

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

Academic Editor: Nistor Victor

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 [1–6] 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 [3–7], 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 [12–14]. 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 that

where 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

5. Conclusions

We have presented an approximation scheme for incompressible miscible displacement in porous media. This scheme was constructed by two methods. This paper is our sequential research work. Cubic Hermite finite element method for the pressure equation was used to ensure the higher regularity of the approximation velocity . A second order characteristic finite element method was presented to handle the material derivative term of the concentration equation. We analyzed the stability of the approximation scheme and derived the optimal-order -norm error estimates for the scalar concentration. They are of second order accuracy in time increment, symmetric, and unconditionally stable.

In the paper, the matrix of the gradient of the characteristic line and the inverse matrix were used to approximate the diffusion term. The properties of these matrices were very important and derived by the complicated theoretical analysis. This paper is the first step of our sequential research work. At current stage, we consider a simple case only for the theoretical aim, which was not related with actual petroleum applications. So, we consider the diffusion coefficient independent of the velocity only in this paper. We will consider the more actual model in petroleum applications later.

Acknowledgments

This paper supported by the Major State Basic Research Development Program of China (Grant no. G19990328), the Scientific Research Award Fund for Excellent Middle-Aged and Young Scientists of Shandong Province (nos. BS2009NJ003 and BS2009HZ015), and the NSF of Shandong Province (no. ZR2010AQ019), China.

References

J. Douglas, Jr., R. E. Ewing, and M. F. Wheeler, “The approximation of the pressure by a mixed method in the simulation of miscible displacement,” RAIRO Analyse Numérique, vol. 17, no. 1, pp. 17–33, 1983.

J. Douglas, Jr., R. E. Ewing, and M. F. Wheeler, “A time-discretization procedure for a mixed finite element approximation of miscible displacement in porous media,” RAIRO Analyse Numérique, vol. 17, no. 3, pp. 249–265, 1983.

R. E. Ewing, T. F. Russell, and M. F. Wheeler, “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, vol. 47, no. 1-2, pp. 73–92, 1984.

T. F. Russell, “Time stepping along characteristics with incomplete iteration for a Galerkin approximation of miscible displacement in porous media,” SIAM Journal on Numerical Analysis, vol. 22, no. 5, pp. 970–1013, 1985.

G. Chavent, “About the identification and modelling of miscible or immiscible displacements in porous media,” in Distributed Parameter Systems: Modelling and Identification, vol. 1 of Lecture Notes in Control and Information Sciences, pp. 196–220, Springer, Berlin, Germany, 1978.

I. Aganovic and A. Mikelic, “On miscible flows in a porous medium,” in Continuum Mechanics and Its Applications, G. A. C. Graham and S. K. Malik, Eds., pp. 569–576, Hemisphere, New York, NY, USA, 1989.

Y. R. Yuan, “Characteristic finite element methods for positive semidefinite problem of two phase miscible flow in three dimensions,” Science Bulletin, vol. 22, pp. 2027–2032, 1996.

Y. R. Yuan, “Characteristic finite difference methods for positive semidefinite problem of two phase miscible flow in porous media,” Systems Science and Mathematical Sciences, vol. 4, pp. 299–306, 1999.

C. N. Dawson, T. F. Russell, and M. F. Wheeler, “Some improved error estimates for the modified method of characteristics,” SIAM Journal on Numerical Analysis, vol. 26, no. 6, pp. 1487–1512, 1989.

M. A. Celia, T. F. Russell, I. Herrera, and R. E. Ewing, “An Eulerian-Lagrangian localized adjoint method for the advection-diffusion equation,” Advances in Water Resources, vol. 13, no. 4, pp. 187–206, 1990.

T. Sun and Y. Yuan, “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, vol. 228, no. 1, pp. 391–411, 2009.

H. Rui and M. Tabata, “A second order characteristic finite element scheme for convection-diffusion problems,” Numerische Mathematik, vol. 92, no. 1, pp. 161–177, 2002.

A. Bermúdez, M. R. Nogueiras, and C. Vázquez, “Numerical analysis of convection-diffusion-reaction problems with higher order characteristics/finite elements. I. Time discretization,” SIAM Journal on Numerical Analysis, vol. 44, no. 5, pp. 1829–1853, 2006.

A. Bermúdez, M. R. Nogueiras, and C. Vázquez, “Numerical analysis of convection-diffusion-reaction problems with higher order characteristics/finite elements. II. Fully discretized scheme and quadrature formulas,” SIAM Journal on Numerical Analysis, vol. 44, no. 5, pp. 1854–1876, 2006.

P. G. Ciarlet, The Finite Element Method for Ellptic Problems, North-Holland Publishing, Amsterdam, The Netherland, 1978.

D. Y. Shi, P. L. Xie, and Z. Y. Yu, “Superclose analysis of the bicubic Hermite element on anisotropic meshes,” Mathematica Numerica Sinica, vol. 30, no. 4, pp. 337–348, 2008.