Considering two-dimensional compressible miscible displacement flow in porous media, finite difference schemes on grids with local refinement in time are constructed and studied. The construction utilizes a modified upwind approximation and linear interpolation at the slave nodes. Error analysis is presented in the maximum norm and numerical examples illustrating the theory are given.

1. Introduction

Numerical models of percolation flow are almost built up on a basis of the finite difference method to solve the system of partial differential equations. Usually, grids that we used are thinner, then the truncation error is smaller and the computing accuracy is higher. In order to assure certain computation accuracy, the grid number cannot be too little. But on the other hand, along with the increment of the grid number, the computation cost is greatly increased and the algebraic system which is formed finally cannot be resolved, even with the largest of today’s supercomputers. Actually, we only need to refine grids around wells, cracks, obstacles, domain boundaries, and so forth, where the pressure changes radically. But because the finite difference grid is composed of straight lines and the grid density cannot be varied with space, it limits the simulating scale and the simulating accuracy. For the local grid refinement technique, we still make use of the finite difference grid system and divide partial grids which are needed to be refined into fine grids. In this way, we can resolve problems, such as small well spacing, fault, and boundary, and we can improve the simulating accuracy and extend the simulating scale [1].

Ewing et al. construct some finite difference approximations on grids with local refinement in space for the ellipse equation and obtain error estimates in the -norm [2]. Cai et al. analyze stationary local grid refinement for the diffusion equation [3, 4]. Ewing et al. derive implicit schemes on the basis of a finite volume approach by approximation of the balance equation. This approach leads to schemes that are locally conservative and are absolutely stable [5]. Ewing et al. construct and study finite difference schemes for transient convection-diffusion problems on grids with local refinement in time and space. The proposed schemes are unconditionally stable and use linear interpolation along the interface [6]. Respectively for incompressible miscible displacement flow in porous media and the semiconductor device problem, authors discuss discrete schemes, error estimates, and numerical examples on composite triangular grids [7, 8].

In this paper, we study a finite difference scheme on grids with local refinement in time for two-dimensional compressible miscible displacement flow in porous media. The pressure equation is approximated by a five-point difference scheme, and the saturation equation is discretized by a modified upwind scheme. At the slave nodes, the construction utilizes linear interpolation. Finally, error analysis in the maximum norm is derived and numerical examples are given to support the numerical method and its convergence.

The paper is organized as follows. In Sections 2 and 3, we formulate the problem and introduce the necessary notations. In Section 4 the construction of the finite difference scheme is presented. The error analysis is addressed in Section 5. Finally, in Section 6 we present numerical experiments that conform our theoretical results.

2. Problem Formulation

We will consider a system of three nonlinear partial differential equations in a bounded domain , which forms a basic model of compressible miscible displacement flow in porous media [911]: where is the saturation of the th component in mixed liquid, is the th component of compression constant factor, is the porosity of the rock, is the permeability of the rock, is the viscosity of the fluid, which is the diffusion matrix, is the diffusion coefficient, and is the unit matrix. The unknowns are the pressure function and the saturation function .

In addition, we have boundary conditions and initial conditions where is a plane bounded domain and is the boundary of .

Usually this question is positive. Suppose the coefficients of (1) satisfy where , , , , , , are positive constants and , , and are Lipschitz continuous in the neighborhood of the solution.

We suppose that the exact solutions of (1) are distributed smoothly; and satisfy Throughout this paper, the notations are used to denote generic constants.

3. Grids, Grid Functions, and Associated Notations

First, is discretized using a regular grid with a parameter . The spatial nodes of the grid on are then defined by , where , . Next, we introduce closed domains , which are subsets of with boundaries aligned with the spatial discretization already defined. Further, it is required that , and we set . In order to avoid unnecessary complications, for ,  , we assume that , where is an integer.

With each subdomain , we associate corresponding sets of nodal points: is defined to be the set of all nodes of the discretization of that are in . We require , for , . And assume that there is no spatial refinement. In each , , we define a subset of boundary nodes as the nodes which have at least one neighbor not in . Then set .

A discrete time-step is associated with each such that, for integers , Consequently, discrete time levels for are defined by , . Finally, we define the grid points by setting We continue by specifying the nodes in between time levels and as Correspondingly, the boundary nodes of are defined by

The grid function is a function defined at the grid points of . we denote the nodal values of a grid function between time levels and as for , , . For we define , and , are the divided forward and backward difference operators, respectively, in and direction. Also define the divided backward time difference by

4. Construction of the Finite Difference Schemes

Let , , and be the numerical approximations to the pressure , the velocity , and the saturation , respectively. The approximation for the pressure and the concentration approximation are done on composite grids in time.

First for the pressure equation, we let Similarly we define , then let

For regular coarse grids, the five-point difference scheme is where the difference operator . The Darcy velocity is computed as follows: corresponds to another direction, and the computational formula is similar to .

Next we consider the saturation equation (1)(c). The positive and negative of the function are defined as and . For regular coarse grids, the upwind difference scheme of the saturation equation is where

In the region that is refined in time, we can construct finite difference schemes similar to (16)–(18). It is obvious that at time , , when the difference operators defined above are applied to the points of , not all space-time positions required correspond to actual nodes in . For such cases, we define In Figure 1, the slave nodes represent the missing space-time positions in the stencil of nodes in . The values there are computed by the interpolation formula (20).

The discretization schemes of (1)–(4) on composite grids are given by , correspond to another direction, and computational formulae are similar to (22).

5. Error Analysis

The discrete inner product and -norm of grid functions are defined, respectively, by We also use the standard notation for the discrete -norm of a grid function in the Sobolev space :

Define the error of the above scheme by

Firstly consider the pressure equation. Using (1)(a) and (21), we get the error equation: Using (27)(b), we can get Then using the maximum principle With the method of recursion and noticing that , we obtain the error estimate:

Similarly using (27)(c) and (27)(d),

Then we consider the saturation equation. Suppose that where is a positive constant. In the end of error analysis, we will prove the supposition (32). Under the supposition (32), we can prove the discretizations (23) satisfy the following property.

Theorem 1. The finite difference schemes (23) comply with the requirements of the maximum principle and the difference operator is coercive in , that is, such that

Considering and using (1)(c), (23), we can get the error equation of the saturation equation

We need an induction hypothesis. Let , we assume that for , . When , we can obtain (36) by using (4) and (23). At the end of error analysis, we will prove (36) for , .

From error estimates of the pressure equation (30) and (31) and the induction hypothesis (36), we have

In order to get an estimate for the error , we need two types of auxiliary functions and . The grid functions and , respectively, satisfy where is the characteristic function of , is the characteristic function of , and is the characteristic function of . On condition that is coercive in , Ewing et al. have proved and satisfy the following lemma [6].

Lemma 2. and exist and are nonnegative, and the following estimates hold:

Theorem 3. Let the exact solutions of (1) satisfy the condition (6), then the discretization scheme (23) is stable, and if the following estimate for the error holds:

Proof. Define where , . Using induction over , it easy to observe that Moreover, Using the maximum principle, it follows that Similarly, Therefore, Then in view of (39), we conclude (40).
It remains to check (32) and the induction hypothesis (36) for . From and the error estimate (40), it is easy to obtain (36). Then using (30) and (31), (34), and (40), we can obtain the supposition (32). The proof of Theorem 3 is complete.

6. Numerical Results

We consider a system of coupled partial differential equations: where , . The following functions are used as exact solutions of (47): When and , exact solutions of (47) are shown in Figures 2 and 3. Specific forms of , , , , and are derived by exact solutions.

From Figures 2 and 3, we can see exact solutions of (47) possess highly localized properties in . First, is discretized using a regular grid. Let the space-step and the time-step . Then choose the subregion , which is refined in time. Let the discretization parameters in the refined region , where is a positive integer. We denote by the numerical approximation to obtained by (23). And let the error estimate in maximum norm .

Example 4. Let . Choosing , computational results obtained by (23) are shown in Figures 4 and 5 and Table 1.

Example 5. Let . Choosing , computational results obtained by (23) are shown in Table 2.

From Figures 4 and 5 and Tables 1 and 2, we can see that numerical results produced by using the local refinement technique are more accurate than those produced without refinement. From Tables 1 and 2, we observe a monotonic improvement of the accuracy in maximum norm when using difference refinement factors . These results are of great importance for the research on numerical simulation of the fluid flow problem and also indicate that the method proposed in this paper can be widely applied to some application fields, such as energy numerical simulation and environmental science.


The author thanks Professor Yirang Yuan for his valuable constructive suggestions which lead to a significant improvement of this paper. This work is supported in part by the National Natural Science Foundation of China (Grant no. 71071088) and the Natural Science Foundation of Shandong Province of China (Grant no. ZR2011AQ021).