Abstract

We present and study a stabilized mixed finite element method for single-phase compressible flow through porous media. This method is based on a pressure projection stabilization method for multiple-dimensional incompressible flow problems by using the lowest equal-order pair for velocity and pressure (i.e., the pair). An optimal error estimate in divergence norm for the velocity and suboptimal error estimates in the -norm for both velocity and pressure are obtained. Numerical results are given in support of the developed theory.

1. Introduction

The mixed finite element method is frequently used to obtain approximate solutions to more than one unknown. For example, the Stokes equations are often solved to obtain both pressure and velocity simultaneously. Accordingly, we need a finite element space for each unknown. These two spaces must be chosen carefully so that they satisfy an inf-sup stability condition for the mixed method to be stable. Examples of an appropriate choice for the mixed spaces for the Stokes equations include the pair (i.e., the Taylor-Hood element [1, 2]) and the MINI-element [1, 3] (i.e., the pair with the addition of the cubic bubble functions on triangles to the -velocity space), where stands for the space of polynomials of degree . For second-order partial differential problems of the Darcy flow type, many special mixed finite element spaces have been established, such as those of Nedelec [4], Raviart and Thomas [5], Brezzi et al. [6], Brezzi et al. [7], Chen and Douglas Jr. [8], and Brezzi et al. [9]. These mixed spaces do not include the equal-order pairs such as .

Although the equal-order pairs of mixed finite element spaces do not satisfy the inf-sup stability condition [10], they offer some computational advances such as they are simple and have practical uniform data structure and adequate accuracy. Much attention has recently been attracted to using the equal-order finite element pairs for the fluid mechanics equations, particularly for the Stokes and Navier-Stokes equations [1016]. Many stabilization techniques have been used to stabilize these element pairs such as penalty, pressure projection, and residual stabilization methods [1619]. Among these methods, the pressure projection stabilization method is a preferable choice in that it is free of stabilization parameters, does not require any calculation of high-order derivatives or edge-based data structures, and can be implemented at the element level [10, 12, 14, 15]. These recent studies have been focused on stabilization of the lowest equal-order finite element pair (or , the bilinear function pair) using the pressure projection stabilization method for second-order elliptic problems of the Darcy flow type and the Stokes and Navier-Stokes problems [10, 12, 15, 20]. This pressure projection stabilization method involves the introduction of a local projection of the pressure into a scalar finite element space.

In this paper, we extend the pressure projection stabilization method to solving the single-phase compressible flow problem in porous media that is described by a second-order parabolic equation. We first define this stabilized method using the equal-order pair of mixed finite element spaces . Then we derive error estimates for this method; optimal error estimates in the divergence norm for velocity and suboptimal error estimates in the -norm for pressure and velocity are obtained. For more information on the better features of the present pressure projection stabilization method over other stabilized mixed methods, please refer to [21]. This is the first time that this method is used for the numerical solution of the single-phase compressible flow problem in porous media.

The rest of this paper is organized as follows. In the next section, the basic notation, the differential equation, and its mixed formulation are stated. Then, in the third section, the stabilized mixed finite element method is shown. Error estimates for this stabilized method are derived in the fourth section, and a superconvergence result is proved in the fifth section. In the final section, numerical experiments are given to illustrate the theoretical results.

2. Function Setting

The dynamical problem we consider is single-phase compressible flow in a porous medium domain , with a Lipschitz-continuous boundary . The system is governed by where is the pressure, the compressibility factor, the source or sink function, the initial pressure, the permeability (or mobility) of the medium, the final time, and . When gravity is ignored, the single-phase compressible flow problem reduces to (2.1)–(2.3) [22]. We assume that is a smooth function: Furthermore, assume that is a bounded, symmetric, and positive definite matrix in ; that is, there exist positive constants and such that where is the transpose of .

In order to introduce a mixed formulation on , set with the norms Introduce and transform (2.1) and (2.3) into the standard mixed form: for any , find and such that This is the mixed variational form of (2.1) and (2.3). The initial condition (2.2) remains valid. A generalized bilinear form on is defined by Consequently, the system (2.8) and (2.9) is written as follows: for any , find and such that

For convenience, we state the Gronwall Lemma that will be used later.

Lemma 2.1. Let , and be three nonnegative functions satisfying, for , where is a nonnegative function on and is a number. Then

Lemma 2.2. In addition to the assumptions (2.4) and (2.5), if then the following stability estimate holds for .

Proof. We divide the proof into three parts.
(1)Taking and in (2.8) and (2.9), respectively, and adding the resulting equations, we obtain which, along with Gronwall's inequality, (2.4), and (2.5), gives (2)Taking in (2.8), differentiating (2.9) with respect to , taking , and adding the resulting equations, we see that which, along with Gronwall's inequality, (2.4), and (2.5), implies (3)Taking in (2.8), we obtain which, along with (2.20), yields The proof has been completed.

3. The Stabilized Mixed Finite Element Method

For , we introduce finite-dimensional subspaces , which are associated with , a triangulation of into triangles, assumed to be regular in the usual sense [1, 23]. This paper focuses on the unstable velocity-pressure pair of the lowest equal-order finite element spaces:

As noted earlier, this choice of the approximate spaces and does not satisfy the condition uniformly in [10]: where the constant is independent of . The pressure projection stabilization method for elliptic problems and the Stokes and Navier-Stokes problems [10, 12, 15, 20] will be adopted.

Let be the standard -projection, which satisfies where . We introduce the pressure projection stabilization term Now, the stabilized mixed finite element method reads as follows: for any , find and such that Introduce the bilinear form Using this notation, the stabilized mixed method of problem (2.8) and (2.9) reads as follows: find , , such that, for all , where is a proper approximation of in . The next theorem can be found in [20], which shows the well-posedness of problem (3.7).

Theorem 3.1. Under the assumption (2.5) and with defined as above, there exist positive constants and β, independent of , such that

4. Error Analysis

To derive error estimates for the mixed finite element solution , we define the mixed method elliptic projection operator by which is well defined by Theorem 3.1 and satisfies the following approximation properties, where we define for simplicity of the presentation. The elliptic projection corresponds to the elliptic problem (without the time differentiation term) of the system (2.8) and (2.9).

Lemma 4.1. Under the assumption of Theorem 3.1, the projection operator satisfies

Proof. Let denote the standard -interplant. Then we see from the definition of , the triangle inequality, (3.9), and (3.10) that which completes the proof.

Lemma 4.2. Under the assumptions of Lemma 2.2 and Theorem 3.1, it holds that

Proof. Subtracting (3.7) from (2.11) with , we have for all . Setting , , and in (4.5) and using (4.1), we see that Then, by integrating (4.6) from 0 to and noting that , we obtain Applying Lemmas 2.1, 2.2, and 4.1 to (4.7) gives As a result, (4.4) follows from Lemma 4.1 and the triangle inequality.

Lemma 4.3. Under the assumption of Lemma 2.2 and Theorem 3.1, it holds that

Proof. Differentiating the term in (4.5) with respect to time and taking in (4.5), we see that Using the Young inequality, we have Then, integrating (4.11) from 0 to and using Lemma 4.1, we obtain which, together with Lemma 4.1 and the triangle inequality, implies (4.9).

Theorem 4.4. Under the assumptions (2.4), (2.5), and (2.15), it holds that

This theorem follows from Lemmas 4.2 and 4.3.

Next, we will estimate . For this, we assume Differentiating (2.8) and (2.9) with respect to , the following result follows from a similar proof to that of Lemma 2.2.

Lemma 4.5. Under the assumptions (2.4), (2.5), and (4.14), the following estimate holds

Lemma 4.6. Under the assumptions of Lemma 4.5, it holds that

Proof. Differentiating (4.5) with respect to time , using (4.1), and taking give Then, by integrating (4.17) from 0 to , we see that Applying Lemmas 2.1, 4.1, and 4.5 to (4.18) gives As a result, (4.16) comes from Lemma 4.1 and the triangle inequality.

Lemma 4.7. Under the assumptions of Lemma 4.5, it holds that

Proof. It follows from the inf-sup condition (3.10), Theorem 3.1, and (4.5) that Consequently, using Lemmas 4.1 and 4.6 and the triangle inequality, we complete the proof.

Theorem 4.8. Under the assumptions (2.4), (2.5), (2.15), and (4.14), it holds that, for ,

This theorem follows from Lemmas 4.2 and 4.7.

The result in Theorem 4.8 is optimal for in the divergence norm and suboptimal for and in the -norm in terms of the convergence rates. These are the best estimates one can obtain with method (3.7). The reason is that the approximation property of the lower-order space always pollutes the global accuracy.

5. Superconvergence

Set Then, by the definition of the projection operator and the fact that , the error equation (4.5) can be written as We recall the projection operator (4.1) as follows: Under the assumption of Lemma 2.2, it follows from [20] that For simplicity of the presentation, we assume that the coefficient is piecewise constant in the next theorem; for a variable coefficient the same result holds if it is projected into the space using the projection operator .

Theorem 5.1. Under the assumptions (2.4), (2.5), and (2.15), it holds that

Proof. Setting , , and in (5.2) and using (5.3), we see that Then, by integrating (5.6) from 0 to and noting that , we obtain Applying Lemmas 2.1 and 2.2 and (5.4) to (5.7) gives As a result, (5.5) comes from (5.4) and the triangle inequality.

6. Numerical Results

Numerical results are presented to check the theory developed in the previous sections. In all the experiments, the triangulations are based on the partition of the unit square into triangles. In the first example, the accuracy of the stabilized mixed finite element method is checked; in the second example, a single-phase flow problem is calculated using this method.

6.1. Accuracy

The purpose of this example is to check the convergence rates for the solution and . Here, let and the coefficient tensor be the identity tensor. We choose the exact solution as follows:

The numerical experiments have been performed by using the equal-order finite element pair for velocity and pressure, as defined in the third section. The backward Euler scheme is used for the time discretization. To check the numerical accuracy in space, we use a small time step ; the final time is .

Detailed numerical results are shown in Tables 1 and 2, which clearly present the expected convergence rates. In these tables, the error, for example, is defined by

A selected numerical pressure and velocity for this example at are shown in Figure 1. From these tables, better computational results, up to the convergence of order for and in the -norm, occur.

6.2. A Single-Phase Flow Problem

Here, we consider the unsteady-state single-phase flow of oil taking place in a two-dimensional, homogeneous, isotropic, horizontal reservoir, with its property parameters given in Table 3 [22]. At the internal boundary, there is a well producing at a constant flow rate. On the other hand, at the external boundary, the pressure remains constant.

The basic differential equation describing this reservoir is In Figure 2, we display the rectangular domain, boundary conditions, and mesh adopted in the simulation.

In Figures 3 and 4, the numerical results are shown. Figures 3(a)3(c) show the variation of the pressure with respect to time; Figure 4 gives the pressure distribution in the reservoir at and  ft after 20 days. The pressure decreases as time progresses, which is reasonable; as the well produces, the reservoir pressure decreases.

7. Conclusions

A stabilized mixed finite element method for an unsteady-state single-phase flow problem in a porous medium has been developed and analyzed. An optimal error estimate in divergence norm for the velocity and suboptimal error estimates in the -norm for both velocity and pressure have been proven. A superconvergence result for the pressure has been obtained as well. The numerical results to check the accuracy of this method and calculate single-phase flow have been presented.

Acknowledgment

This work is supported in part by the NSERC/AERI/iCORE/Foundation CMG Chair Funds in Reservoir Simulation.