A Stabilized Mixed Finite Element Method for Single-Phase Compressible Flow

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 P1 − P1 pair . An optimal error estimate in divergence norm for the velocity and suboptimal error estimates in the L2-norm for both velocity and pressure are obtained. Numerical results are given in support of the developed theory.


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 P 2 − P 1 pair i.e., the Taylor-Hood element 1, 2 and the MINI-element 1, 3 i.e., the P 1 − P 1 pair with the addition of the cubic bubble functions on triangles to the P 1 -velocity space , where P r stands for the space of polynomials of degree r ≥ 0. 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 P r − P r .Furthermore, assume that K k ij d i,j 1 is a bounded, symmetric, and positive definite matrix in Ω; that is, there exist positive constants α 1 and α 2 such that where ξ T is the transpose of ξ.
In order to introduce a mixed formulation on Ω, set with the norms Introduce u −K∇p and transform 2.1 and 2.3 into the standard mixed form: for any t ∈ 0, T , find u t ∈ V and p t ∈ W such that cp t , q div u, q f, q , ∀q ∈ W, 2.8 This is the mixed variational form of 2.1 and 2.3 .The initial condition 2.2 remains valid.A generalized bilinear form on V, W × V, W is defined by Consequently, the system 2.8 and 2.9 is written as follows: for any t ∈ 0, T , find u t ∈ V and p t ∈ W such that For convenience, we state the Gronwall Lemma that will be used later.
Proof.We divide the proof into three parts.
1 Taking q p and v u in 2.8 and 2.9 , respectively, and adding the resulting equations, we obtain 1 2 2 Taking q p t in 2.8 , differentiating 2.9 with respect to t, taking v u, and adding the resulting equations, we see that The proof has been completed.

The Stabilized Mixed Finite Element Method
For h > 0, we introduce finite-dimensional subspaces V h , W h ⊂ V, W , which are associated with Γ h , 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 equalorder finite element spaces:

3.1
As noted earlier, this choice of the approximate spaces V h and W h does not satisfy the inf-sup condition uniformly in h 10 : where the constant β > 0 is independent of h.The pressure projection stabilization method for elliptic problems and the Stokes and Navier-Stokes problems 10, 12, 15, 20 will be adopted.Let Π h : W → W h be the standard L 2 -projection, which satisfies

3.3
where Now, the stabilized mixed finite element method reads as follows: for any t ∈ 0, T , find

3.5
Introduce the bilinear form

3.6
Using this notation, the stabilized mixed method of problem 2.8 and 2.9 reads as follows: where p 0h is a proper approximation of p 0 in W h .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 V h , W h defined as above, there exist positive constants C and β, independent of h, such that 3.10

Error Analysis
To derive error estimates for the mixed finite element solution u h , p h , we define the mixed method elliptic projection operator R h , which is well defined by Theorem 3.1 and satisfies the following approximation properties, where we define u 0h , p 0h R h u 0 , p 0 , Q h u 0 , p 0 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
Proof.Let I h denote the standard P 1 -interplant.Then we see from the definition of R h , Q h , 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 v, q v h , q h , we have and v h , q h 2 h , η h in 4.5 and using 4.1 , we see that

4.6
Then, by integrating 4.6 from 0 to t and noting that η h 0 0, we obtain Proof.Differentiating the term 5 with respect to time t and taking v h , q h h , η ht in 4.5 , we see that Using the Young inequality, we have Then, integrating 4.11 from 0 to t and using Lemma 4.1, we obtain

4.21
Consequently, using Lemmas 4.1 and 4.6 and the triangle inequality, we complete the proof.The result in Theorem 4.8 is optimal for u in the divergence norm and suboptimal for u and p in the L 2 -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 W h always pollutes the global accuracy.

Set e h t
Π h p t − p h t , t ∈ 0, T .

5.1
Then, by the definition of the projection operator Π h and the fact that div V h W h , the error equation 4.5 can be written as We recall the projection operator 4.1 as follows:

5.3
Under the assumption of Lemma 2.2, it follows from 20 that For simplicity of the presentation, we assume that the coefficient c is piecewise constant in the next theorem; for a variable coefficient c the same result holds if it is projected into the space W h using the projection operator Π h .
Theorem 5.1.Under the assumptions 2.4 , 2.5 , and 2.15 , it holds that , and v h , q h 2 h , Π h η h in 5.2 and using 5.3 , we see that Then, by integrating 5.6 from 0 to t and noting that η h 0 0, we obtain Applying Lemmas 2.1 and 2.2 and 5.4 to 5.7 gives

5.8
As a result, 5.5 comes from 5.4 and the triangle inequality.

Numerical Results
Numerical results are presented to check the theory developed in the previous sections.In all the experiments, the triangulations Γ h are based on the partition of the unit square Ω 0, 1 × 0, 1 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.

Accuracy
The purpose of this example is to check the convergence rates for the solution p and u.Here, let c 1 and the coefficient tensor K be the identity tensor.We choose the exact solution as follows: p x, y x 2 x − 1 y y − 1 cos t .6.1 The numerical experiments have been performed by using the equal-order finite element pair P 1 − P 1 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 dt 0.0005; the final time is T 1.
Detailed numerical results are shown in Tables 1 and 2, which clearly present the expected convergence rates.In these tables, the p L 2 error, for example, is defined by p L 2 error p − p h 0 .6.2 A selected numerical pressure and velocity for this example at T 1 are shown in Figure 1.From these tables, better computational results, up to the convergence of order O h 2 for p and u in the L 2 -norm, occur.

A Single-Phase Flow Problem
Here, we consider the unsteady-state single-phase flow of oil taking place in a twodimensional, 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 r r w and r 200 ft after 20 days.The pressure decreases as time progresses, which is reasonable; as the well produces, the reservoir pressure decreases.

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

Figure 2 :c
Figure 2: Conditions adopted in simulation.

Figure 3 :
Figure 3: Numerical results pressure contours with respect to time.

Figure 4 :
Figure 4: Numerical results at r r w and r 200 ft.
In addition to the assumptions 2.4 and 2.5 , if Lemma 2.1.Let g t , l t , and ξ t be three nonnegative functions satisfying, for t ∈ 0, T , where G t is a nonnegative function on 0, T and C is a number.Then ξ t G t ≤ C t 0 l ds exp t 0 g ds , t ∈ 0, T .2.14 Lemma 2.2.
Under the assumptions 2.4 , 2.5 , and 4.14 , the following estimate holds Theorem 4.4.Under the assumptions 2.4 , 2.5 , and 2.15 , it holds that u t − u h t 0 p t − p h t 0 ≤ Ch, t ∈ 0, T .Differentiating 2.8 and 2.9 with respect to t, the following result follows from a similar proof to that of Lemma 2.2.Lemma 4.5.ht , η ht ds ≤ Ch 2 .4.19 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 u − u h H div;Ω ≤ Ch, t ∈ 0, T .4.20 Proof.It follows from the inf-sup condition 3.10 , Theorem 3.1, and 4.5 that Theorem 4.8.Under the assumptions 2.4 , 2.5 , 2.15 , and 4.14 , it holds that, for t ∈ 0, T ,

Table 1 :
Errors and convergence rates for the velocity.

Table 2 :
Errors and convergence rates for the pressure.

Table 3 :
Parameters for a reservoir.