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

Authors: Liyun Zhang; Zhangxin Chen
Journal: Journal of Applied Mathematics
 (2011)

---

## 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., theP1−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.


---

## Body

## 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 theP2-P1 pair (i.e., the Taylor-Hood element [1, 2]) and the MINI-element [1, 3] (i.e., the P1-P1 pair with the addition of the cubic bubble functions on triangles to the P1-velocity space), where Pr 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 Pr-Pr.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 [10–16]. Many stabilization techniques have been used to stabilize these element pairs such as penalty, pressure projection, and residual stabilization methods [16–19]. 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 P1-P1 (or Q1-Q1, 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 spacesP1-P1. Then we derive error estimates for this method; optimal error estimates in the divergence norm for velocity and suboptimal error estimates in the L2-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Ω⊂ℛd,d≥2, with a Lipschitz-continuous boundary Γ. The system is governed by(2.1)cpt-∇⋅(K∇p)=f,inΩ×(0,T],(2.2)p=p0,inΩ×{0},(2.3)p=0,onΓ×(0,T],
where p(x,t) is the pressure, c(x,t) the compressibility factor, f(x,t) the source or sink function, p0(x) the initial pressure, K(x) the permeability (or mobility) of the medium, T>0 the final time, and pt=∂p/∂t. When gravity is ignored, the single-phase compressible flow problem reduces to (2.1)–(2.3) [22]. We assume that c(x,t) is a smooth function:(2.4)0<c0≤c(x,t)<∞,ct(x,t)<∞,x∈Ω,0≤t≤T.
Furthermore, assume that K=(kij)i,j=1d is a bounded, symmetric, and positive definite matrix in Ω; that is, there exist positive constants α1 and α2 such that(2.5)α1‖ξ‖2≤ξTKξ≤α2‖ξ‖2,∀ξ∈Rd,d≥2,
where ξT is the transpose of ξ.In order to introduce a mixed formulation onΩ, set(2.6)V=H(div⁡;Ω)={v∈(L2(Ω))d:div⁡v∈L2(Ω)},W=L2(Ω),
with the norms(2.7)‖v‖H(div⁡;Ω)=(‖v‖02+‖div⁡v‖02)1/2,‖q‖W=‖q‖0,v∈V,q∈W.
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
(2.8)(cpt,q)+(div⁡u,q)=(f,q),∀q∈W,(2.9)(K-1u,v)-(p,div⁡v)=0,∀v∈V.
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(2.10)B((u,p),(v,q))=(K-1u,v)-(p,div⁡v)+(div⁡u,q),(u,p),(v,q)∈V×W.
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
(2.11)(cpt,q)+B((u,p),(v,q))=(f,q),∀(v,q)∈V×W,(2.12)p(x,0)=p0(x),∀x∈Ω.For convenience, we state the Gronwall Lemma that will be used later.Lemma 2.1.
Letg(t),l(t), and ξ(t) be three nonnegative functions satisfying, for t∈[0,T],(2.13)ξ(t)+G(t)≤C+∫0tlds+∫0tgξds,
where G(t) is a nonnegative function on [0,T] and C is a number. Then
(2.14)ξ(t)+G(t)≤(C+∫0tlds)exp⁡(∫0tgds),t∈[0,T].Lemma 2.2.
In addition to the assumptions (2.4) and (2.5), if(2.15)‖p0‖0+(∫0t‖f(s)‖02ds)1/2≤C,t∈[0,T],
then the following stability estimate holds
(2.16)sup⁡0≤s≤t(‖p(s)‖02+‖u(s)‖02)+∫0t(‖pt(s)‖02+‖u(s)‖02+‖div⁡u(s)‖02)ds≤C,
for t∈[0,T].Proof.
We divide the proof into three parts.
(1)
Takingq=p and v=u in (2.8) and (2.9), respectively, and adding the resulting equations, we obtain
(2.17)12ddt‖cp‖02+(K-1u,u)=(f,p)+12(ctp,p)≤‖c-1/2f‖0‖cp‖0+C‖p‖02,
which, along with Gronwall's inequality, (2.4), and (2.5), gives
(2.18)sup⁡0≤s≤t‖p(s)‖02+∫0t‖u(s)‖02ds≤C,t∈[0,T].(2)
Takingq=pt in (2.8), differentiating (2.9) with respect to t, taking v=u, and adding the resulting equations, we see that
(2.19)‖cpt‖02+12ddt‖K-1u‖02=(f,pt)≤‖c-1/2f‖0‖cpt‖0,
which, along with Gronwall's inequality, (2.4), and (2.5), implies
(2.20)sup⁡0≤s≤t‖u(s)‖02+∫0t‖pt(s)‖02ds≤C,t∈[0,T].(3)
Takingq=divu in (2.8), we obtain
(2.21)‖div⁡u‖02=(f,div⁡u)-(cpt,div⁡u)≤(‖f‖0+‖cpt‖0)‖div⁡u‖0,which, along with (2.20), yields
(2.22)∫0t‖div⁡u(s)‖02ds≤C.
The proof has been completed.


## 3. The Stabilized Mixed Finite Element Method

Forh>0, we introduce finite-dimensional subspaces (Vh,Wh)⊂(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 equal-order finite element spaces:(3.1)Vh={v∈V:v∣T∈(P1(T))2,∀T∈Γh},Wh={q∈W:q∣T∈P1(T),∀T∈Γh}.As noted earlier, this choice of the approximate spacesVh and Wh does not satisfy the inf-sup condition uniformly in h [10]:(3.2)sup⁡v∈Vh(div⁡v,q)‖v‖H(div⁡;Ω)≥β‖q‖0,∀q∈Wh,
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 L2-projection, which satisfies
(3.3)(p,q)=(Πhp,q),∀p∈W,q∈W¯h,‖Πhp‖0≤‖p‖0,∀p∈W,
where W¯h={q∈W:q∣T∈P0(T),∀T∈Γh}. We introduce the pressure projection stabilization term(3.4)G(p,q)=(p-Πhp,q-Πhq),p,q∈W.
Now, the stabilized mixed finite element method reads as follows: for any t∈[0,T], find uh(t)∈Vh and ph(t)∈Wh such that(3.5)(cpht,q)+(div⁡uh,q)+G(ph,q)=(f,q),∀q∈Wh,(K-1uh,v)-(ph,div⁡v)=0,∀v∈Vh.
Introduce the bilinear form(3.6)B((uh,ph),(v,q))=(K-1uh,v)-(ph,div⁡v)+(div⁡uh,q)+G(ph,q),(uh,ph),(v,q)∈Vh×Wh.
Using this notation, the stabilized mixed method of problem (2.8) and (2.9) reads as follows: find (uh,ph)∈Vh×Wh, t∈[0,T], such that, for all (v,q)∈Vh×Wh,(3.7)(cpht,q)+B((uh,ph),(v,q))=(f,q),(3.8)ph(0)=p0h,
where p0h is a proper approximation of p0 in Wh. 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 (Vh,Wh) defined as above, there exist positive constants C and β, independent of h, such that(3.9)|B((uh,ph),(vh,qh))|≤C(‖uh‖H(div⁡;Ω)+‖ph‖0)×(‖vh‖H(div⁡;Ω)+‖qh‖0),(3.10)sup⁡(vh,qh)∈(Vh,Wh)|B((uh,ph),(vh,qh))|‖vh‖H(div⁡;Ω)+‖qh‖0≥β(‖uh‖H(div⁡;Ω)+‖ph‖0),∀(uh,ph),(vh,qh)∈(Vh,Wh).


## 4. Error Analysis

To derive error estimates for the mixed finite element solution(uh,ph), we define the mixed method elliptic projection operator (Rh,Qh):(V,W)→(Vh,Wh) by(4.1)B((Rh(v,q),Qh(v,q)),(vh,qh))=B((v,q),(vh,qh)),∀(v,q)∈(V,W),(vh,qh)∈(Vh,Wh),
which is well defined by Theorem 3.1 and satisfies the following approximation properties, where we define (u0h,p0h)=(Rh(u0,p0),Qh(u0,p0)) 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 Theorem3.1, the projection operator (Rh,Qh) satisfies(4.2)‖v-Rh(v,q)‖H(div⁡;Ω)+‖q-Qh(v,q)‖0≤Ch(‖v‖1+‖div⁡v‖1+‖q‖1),∀(v,q)∈(V,W).Proof.
LetIh denote the standard P1-interplant. Then we see from the definition of (Rh,Qh), the triangle inequality, (3.9), and (3.10) that(4.3)‖v-Rh(v,q)‖H(div⁡;Ω)+‖q-Qh(v,q)‖0≤‖v-Ihv‖H(div⁡;Ω)+‖q-Πhq‖0+‖Ihv-Rh(v,q)‖H(div⁡;Ω)+∥Πhq-Qh(v,q)∥0≤‖v-Ihv‖H(div⁡;Ω)+‖q-Πhq‖0+β-1sup⁡(xh,yh)∈(Vh,Wh)|B((Ihv-Rh(v,q),Πhq-Qh(v,q)),(xh,yh))|‖xh‖H(div⁡;Ω)+‖yh‖0≤‖v-Ihv‖H(div⁡;Ω)+‖q-Πhq‖0+β-1sup⁡(xh,yh)∈(vh,Wh)|B((v-Ihv,q-Πhq),(xh,yh))|‖xh‖H(div⁡;Ω)+‖yh‖0≤C(‖v-Ihv‖H(div⁡;Ω)+‖q-Πhq‖0)≤Ch(‖v‖1+‖div⁡v‖1+‖q‖1),
which completes the proof.Lemma 4.2.
Under the assumptions of Lemma2.2 and Theorem 3.1,  it holds that(4.4)‖p(t)-ph(t)‖02+∫0t‖u(s)-uh(s)‖02ds≤Ch2,t∈[0,T].Proof.
Subtracting (3.7) from (2.11) with (v,q)=(vh,qh), we have(4.5)(c[pt-pht],qh)+B((u-uh,p-ph),(vh,qh))=G(p,qh),
for all (vh,qh)∈(Vh,Wh). Setting (ϵh,ηh)=(Rh(u,p)-uh,Qh(u,p)-ph), E=p-Qh(u,p), and (vh,qh)=2(ϵh,ηh) in (4.5) and using (4.1), we see that
(4.6)ddt‖cηh‖02+2(K-1ϵh,ϵh)+2G(ηh,ηh)=-2(cEt,ηh)+(ctηh,ηh)≤2‖c-1/2Et‖0‖cηh‖0+(ctηh,ηh).
Then, by integrating (4.6) from 0 to t and noting that ηh(0)=0, we obtain
(4.7)‖ηh‖02+∫0t((K-1ϵh,ϵh)+G(ηh,ηh))ds≤C(∫0t‖Et‖02ds+∫0t‖ηh‖02ds).
Applying Lemmas 2.1, 2.2, and 4.1 to (4.7) gives
(4.8)‖ηh‖02+∫0t((K-1ϵh,ϵh)+G(ηh,ηh))ds≤Ch2.
As a result, (4.4) follows from Lemma 4.1 and the triangle inequality.Lemma 4.3.
Under the assumption of Lemma2.2 and Theorem 3.1, it holds that(4.9)‖u(t)-uh(t)‖02+∫0t‖pt-pht‖2ds≤Ch2,t∈[0,T].Proof.
Differentiating the term(K-1(u-uh),vh)-(p-ph,divvh) in (4.5) with respect to time t and taking (vh,qh)=(ϵh,ηht) in (4.5), we see that(4.10)‖cηht‖02+12ddt((K-1ϵh,ϵh)+G(ηh,ηh))=-(cEt,ηht).
Using the Young inequality, we have
(4.11)‖ηht‖02+ddt((K-1ϵh,ϵh)+G(ηh,ηh))≤C‖Et‖02.
Then, integrating (4.11) from 0 to t and using Lemma 4.1, we obtain
(4.12)(K-1ϵh,ϵh)+G(ηh,ηh)+∫0t‖ηht‖02ds≤C∫0t‖Et‖02ds≤Ch2,
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(4.13)‖u(t)-uh(t)‖0+‖p(t)-ph(t)‖0≤Ch,t∈[0,T].This theorem follows from Lemmas4.2 and 4.3.Next, we will estimate∥u(t)-uh(t)∥H(div;Ω). For this, we assume(4.14)‖p0‖0+(∫0T‖ft(t)‖02dt)1/2≤C.
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.
Under the assumptions (2.4), (2.5), and (4.14), the following estimate holds(4.15)sup⁡0≤t≤T(‖pt(t)‖02+‖ut(t)‖02)+∫0T(‖ptt(t)‖02+‖ut(t)‖02+‖div⁡ut(t)‖02)dt≤C.Lemma 4.6.
Under the assumptions of Lemma4.5, it holds that(4.16)‖pt(t)-pht(t)‖02+∫0t‖ut(s)-uht(s)‖02ds≤Ch2,t∈[0,T].Proof.
Differentiating (4.5) with respect to time t, using (4.1), and taking (vh,qh)=2(ϵht,ηht) give(4.17)ddt‖cηht‖02+2(K-1ϵht,ϵht)+2G(ηht,ηht)=-2(ctEt,ηht)-2(cEtt,ηht)-(ctηht,ηht)≤C(‖Et‖0+‖Ett‖0+‖ηht‖0)‖ηht‖0.
Then, by integrating (4.17) from 0 to t, we see that
(4.18)‖ηht‖02+∫0t((K-1ϵht,ϵht)+G(ηht,ηht))ds≤C∫0t(‖Et‖02+‖Ett‖02+‖ηht‖02)ds.
Applying Lemmas 2.1, 4.1, and 4.5 to (4.18) gives
(4.19)‖ηht‖02+∫0t((K-1ϵht,ϵht)+G(ηht,ηht))ds≤Ch2.
As a result, (4.16) comes from Lemma 4.1 and the triangle inequality.Lemma 4.7.
Under the assumptions of Lemma4.5, it holds that(4.20)‖u-uh‖H(div⁡;Ω)≤Ch,t∈[0,T].Proof.
It follows from the  inf-sup condition (3.10), Theorem 3.1, and (4.5) that(4.21)‖ϵh‖0+‖div⁡ϵh‖0+‖ηh‖0≤β-1sup⁡(vh,wh)∈(Vh,Wh)|B((ϵh,ηh),(vh,wh))|‖vh‖H(div⁡;Ω)+‖wh‖0≤C‖pt-pht‖0.
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 t∈[0,T],(4.22)‖u(t)-uh(t)‖0+‖div⁡(u(t)-uh(t))‖0+‖p(t)-ph(t)‖0≤Ch.This theorem follows from Lemmas4.2 and 4.7.The result in Theorem4.8 is optimal for u in the divergence norm and suboptimal for u and p in the L2-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.


## 5. Superconvergence

Set(5.1)eh(t)=Πh(p(t)-ph(t)),t∈[0,T].
Then, by the definition of the projection operator Πh and the fact that divVh=W¯h, the error equation (4.5) can be written as(5.2)(c[pt-pht],qh)+B((u-uh,eh),(vh,qh))=G(p,qh),(vh,qh)∈(Vh,Wh).
We recall the projection operator (4.1) as follows:(5.3)B((v-Rh(v,q),Πh(q-Qh(v,q))),(vh,qh))=G(q,qh),∀(v,q)∈(V,W),(vh,qh)∈(Vh,Wh).
Under the assumption of Lemma 2.2, it follows from [20] that(5.4)‖Πh(p-Qh(u,p))‖0≤Ch2.
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(5.5)‖Πh(p(t)-ph(t))‖0≤Ch2,t∈[0,T].Proof.
Setting(ϵh,ηh)=(Rh(u,p)-uh,Qh(u,p)-ph), E=p-Qh(u,p), and (vh,qh)=2(ϵh,Πhηh) in (5.2) and using (5.3), we see that(5.6)ddt‖cΠhηh‖02+2(K-1ϵh,ϵh)≤C‖ΠhEt‖0‖Πhηh‖0.
Then, by integrating (5.6) from 0 to t and noting that ηh(0)=0, we obtain
(5.7)‖Πhηh‖02+∫0t(K-1ϵh,ϵh)ds≤C(∫0t‖ΠhEt‖02ds+∫0t‖Πhηh‖02ds).
Applying Lemmas 2.1 and 2.2 and (5.4) to (5.7) gives
(5.8)‖Πhηh‖02+∫0t(K-1ϵh,ϵh)ds≤Ch4.
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Γ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.


### 6.1. Accuracy

The purpose of this example is to check the convergence rates for the solutionp and u. Here, let c=1 and the coefficient tensor K be the identity tensor. We choose the exact solution as follows: (6.1)p(x,y)=x2(x-1)y(y-1)cos⁡(t).The numerical experiments have been performed by using the equal-order finite element pairP1-P1 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 Tables1 and 2, which clearly present the expected convergence rates. In these tables, the pL2 error, for example, is defined by (6.2)pL2error=‖p-ph‖0.Table 1
Errors and convergence rates for the velocity.
1/huL2 erroruL2 rateudiv errorudiv rate80.06266820.0953064160.01752181.838580.05145470.889269240.008196441.873770.03551130.914632320.00476151.887980.0272290.923136400.003118951.895940.02213710.927761Table 2
Errors and convergence rates for the pressure.
1/hpL2 errorpL2rateΠhL2 errorΠhL2rate80.03038170.0200016160.007425622.032620.004646062.10603240.003284282.011980.002031242.04055320.001845482.003630.001137892.01428400.001181711.99770.0007287811.99673A selected numerical pressure and velocity for this example atT=1 are shown in Figure 1. From these tables, better computational results, up to the convergence of order O(h2) for p and u in the L2-norm, occur.Numerical results on a 32 × 32 uniform triangular mesh.
(a)
The pressure and velocity(b)
The pointwise pressure error


### 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 Table3 [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.Table 3
Parameters for a reservoir.
ItemDescriptionUnitValuekPermeabilitymd88.7μOil viscositycp10HThicknessft100cCompressibility1/psi0.0002ϕPorosityfraction0.2p0Initial pressurepsia4,000qOil flow rateSTB/D400rwRadius of wellboreft0.25BOil FVFRB/STB1αcTransmissibility conversion factor—1.127βcVolume conversion factor—5.614583xmaxLength in thex-directionft4,000ymaxLength in they-directionft4,000The basic differential equation describing this reservoir is(6.3)ϕμcβcαck∂p∂t=∂2p∂x2+∂2p∂y2.
In Figure 2, we display the rectangular domain, boundary conditions, and mesh adopted in the simulation.Conditions adopted in simulation.
(a)
Domain and boundary condition(b)
Base trianglesIn Figures3 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=rw and r=200 ft after 20 days. The pressure decreases as time progresses, which is reasonable; as the well produces, the reservoir pressure decreases.Numerical results (pressure contours) with respect to time.
(a)
Pressure at timet=1 day(b)
Pressure at timet=10 days(c)
Pressure at timet=20 daysFigure 4
Numerical results atr=rw and r=200ft.


## 6.1. Accuracy

The purpose of this example is to check the convergence rates for the solutionp and u. Here, let c=1 and the coefficient tensor K be the identity tensor. We choose the exact solution as follows: (6.1)p(x,y)=x2(x-1)y(y-1)cos⁡(t).The numerical experiments have been performed by using the equal-order finite element pairP1-P1 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 Tables1 and 2, which clearly present the expected convergence rates. In these tables, the pL2 error, for example, is defined by (6.2)pL2error=‖p-ph‖0.Table 1
Errors and convergence rates for the velocity.
1/huL2 erroruL2 rateudiv errorudiv rate80.06266820.0953064160.01752181.838580.05145470.889269240.008196441.873770.03551130.914632320.00476151.887980.0272290.923136400.003118951.895940.02213710.927761Table 2
Errors and convergence rates for the pressure.
1/hpL2 errorpL2rateΠhL2 errorΠhL2rate80.03038170.0200016160.007425622.032620.004646062.10603240.003284282.011980.002031242.04055320.001845482.003630.001137892.01428400.001181711.99770.0007287811.99673A selected numerical pressure and velocity for this example atT=1 are shown in Figure 1. From these tables, better computational results, up to the convergence of order O(h2) for p and u in the L2-norm, occur.Numerical results on a 32 × 32 uniform triangular mesh.
(a)
The pressure and velocity(b)
The pointwise pressure error


## 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 Table3 [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.Table 3
Parameters for a reservoir.
ItemDescriptionUnitValuekPermeabilitymd88.7μOil viscositycp10HThicknessft100cCompressibility1/psi0.0002ϕPorosityfraction0.2p0Initial pressurepsia4,000qOil flow rateSTB/D400rwRadius of wellboreft0.25BOil FVFRB/STB1αcTransmissibility conversion factor—1.127βcVolume conversion factor—5.614583xmaxLength in thex-directionft4,000ymaxLength in they-directionft4,000The basic differential equation describing this reservoir is(6.3)ϕμcβcαck∂p∂t=∂2p∂x2+∂2p∂y2.
In Figure 2, we display the rectangular domain, boundary conditions, and mesh adopted in the simulation.Conditions adopted in simulation.
(a)
Domain and boundary condition(b)
Base trianglesIn Figures3 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=rw and r=200 ft after 20 days. The pressure decreases as time progresses, which is reasonable; as the well produces, the reservoir pressure decreases.Numerical results (pressure contours) with respect to time.
(a)
Pressure at timet=1 day(b)
Pressure at timet=10 days(c)
Pressure at timet=20 daysFigure 4
Numerical results atr=rw and r=200ft.


## 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 theL2-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.



---

*Source: 129724-2011-01-23.xml*