# A Rectangular Mixed Finite Element Method with a Continuous Flux for an Elliptic Equation Modelling Darcy Flow

Authors: Xindong Li; Hongxing Rui
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We introduce a mixed finite element method for an elliptic equation modelling Darcy flow in porous media. We use a staggered mesh where the two components of the velocity and the pressure are defined on three different sets of grid nodes. In the present mixed finite element, the approximate velocity is continuous and the conservation law still holds locally. The LBB consistent condition is established, while theL2 error estimates are obtained for both the velocity and the pressure. Numerical examples are presented to confirm the theoretical analysis.


---

## Body

## 1.  Introduction

We consider the discretization technique for the elliptic problem modelling the flow in saturated porous media, or the classical Darcy flow problem, including a system of mass conservation law and Darcy’s law [1, 2]. The most popular numerical methods for this elliptic equation focus on mixed finite element methods, since by this kind of methods the original scalar variable, called pressure, and its vector flux, named Darcy velocity, can be approximated simultaneously and maintain the local conservation. The classical theory for the mixed finite element, which includes the LBB consistent condition, the existence and uniqueness of the approximate solution, and the error estimate, has been established. Some mixed finite element methods such as RT mixed finite element and BDM mixed finite element are introduced (as in [3–6]), which satisfy the consistent condition and have optimal order error estimate [7, 8]. Give some stabilized mixed finite methods by adding to the classical mixed formulation some least squares residual forms of the governing equations.By using the abovementioned mixed finite element methods, the approximate velocity is continuous in the normal direction and discontinuous in the tangential direction on the edges of the element. This is reasonable for the case of heterogenous permeability, yet it is desirable that the flux be continuous in some applications [9]. In particular, when we track the characteristic segment using the approximate velocity, the discontinuities of the velocity may introduce some difficulties when the characteristic line cross the edges of element. While applying mass-conservative characteristic finite element method to the coupled system of compressible miscible displacement in porous media, the continuous flux is crucial [10]. A brief description of this point will be found at the last part of this paper.To overcome this disadvantage, Arbogast and Wheeler [11] introduced a mixed finite element method with an approximate velocity continuous in both the normal direction and the tangential direction, which was got by adding some freedom to the RT mixed finite element. In this paper, we introduced a mixed finite element method with an approximate velocity continuous in all directions. It is based on rectangular mesh and uses continuous piecewise bilinear functions to approximate the velocity components and uses piecewise constant functions to approximate the pressure. We obtain the element by improving a kind of element for Stokes equation and Navier-Stokes equation given by Han [12], Han and Wu [13], and Han and Yan [14]. By using this mixed finite element, we can get continuous velocity vector and maintain the local conservation. Comparing to the mixed finite element method in [11], we need less degrees of freedom for the same convergence rate. The LBB consistent condition and L2 error estimates of velocity and pressure are also established.The outline of the rest of this paper is organized as follows. In Sections2 and 3, we recall the model problem and weak formulation for the mixed finite element method and then establish the discrete inf-sup consistent condition and L2 error estimates for the velocity and the pressure in Section 4. In Section 5, we present some numerical examples which verify the efficiency of the proposed mixed finite element method. A valuable application of this method to mass-conservative characteristic (MCC) scheme for the coupled compressible miscible displacement in porous media closes the paper in Section 6.


## 2. The Mixed Finite Formulation for Darcy Equation

The mathematical model for viscous flow in porous media includes Darcy’s law and conservation law of mass, written as follows:(1)u=-κμ∇ponΩ(Darcy'slaw)divu=ϕonΩ(massconservation)u·n=0onΓ,
where κ>0 is the permeability, μ>0 is the viscosity, and ϕ is the volumetric flow rate source or sink. Γ is the boundary of Ω, and n is the unit outward normal vector to Γ. The variable u=(u1,u2) is the Darcy velocity vector, and p is the pressure. The source ϕ must satisfy the consistency constraint
(2)∫Ω‍ϕdΩ=0.LetL2(Ω) be the space of square integrable function in Ω with inner product (·,·) and norm ∥·∥. We use the notation of the Hilbert space
(3)H(div,Ω)={u∈[L2(Ω)]2;divu∈L2(Ω)},
with norm
(4)∥u∥H(div,Ω)={∥u∥2+∥divu∥2}1/2.
Define the following subspaces of H(div,Ω) and L2(Ω):
(5)V=H0(div,Ω)={u∈H(div,Ω):u·n=0onΓ},S={q∣q∈L2(Ω):∫Ω‍qdΩ=0}.
The classical weak variational formulation of Problem (1) is as follows: find (u,p)∈V×S, such that
(6)a(u,v)-b(v,p)=0∀v∈V,b(u,q)=(ϕ,q)∀q∈S.
Here,
(7)a(u,v)=∫Ω‍μκu·vdxb(v,q)=∫Ω‍qdivvdx.The following discussion and discrete analysis are related to the weak form (6). Let V0 be a closed subspace of V via
(8)V0={v∈V:b(v,q)=0,∀q∈S}.
For the bilinear forms a(u,v) and b(v,q), we have the standard result.Lemma 1.
The bilinear forma(u,v) is bounded on V×V and coercive on V0, and the bilinear form b(v,q) is bound on V×S. Namely,(1)
there exist two constantsC1>0 and α>0 such that
(9)|a(u,v)|≤C1∥u∥H(div,Ω)∥v∥H(div,Ω)∀u,v∈V,a(u,u)≥α∥u∥H(div,Ω)2∀u∈V0,(2)
there is a constantC2>0 such that
(10)|b(v,q)|≤C2∥q∥0,Ω∥v∥H(div,Ω)∀q∈S,v∈V.For the spaceV and S, the Ladyzhenskaya-Babus˘ka-Brezzi(L-B-B) condition holds; see [15, 16], for example.Lemma 2.
There is a constantβ>0 such that
(11)supv∈Vb(v,q)∥v∥H(div,Ω)≥β∥q∥0,Ω,∀q∈S.
It is clear that there exists a unique solution(u,p)∈V×S to the Problem (6).


## 3. Finite Element Discretization

In this section, we present the mixed finite element based on rectangular mesh for the Darcy flow problem.In [13], Han and Wu introduced a mixed finite element for Stokes problem and then extended to solve the Navier-Stokes problem [14]. Based on this element, we introduced the new mixed finite element with a continuous flux approximation for Darcy flow problem.For simplicity, we suppose that the domainΩ is a unit square, and the mixed finite element discussed here can be easily generalized to the case when the domain Ω is a rectangular.LetN be a given integer and h=1/N. We construct the finite-dimensional subspaces of S and V by introducing three different quadrangulations τh, τh1, τh2 of Ω.First, we divideΩ into uniform squares
(12)Ti,j={(x,y):xi-1≤x≤xi,yj-1≤y≤yj},i,j=1,…,N,
where xi=ih and yj=jh. The corresponding quadrangulation is denoted by τh. See Figure 1(a).
(13)τi,j={Ti,j:i,j=1,…,N}.Quadrangulations: (a)τh, (b)τh1, and (c)τh2.
(a)
(b)
(c)Then, for allTi,j∈τh, we connect all the neighbor midpoints of the vertical sides of Ti,j by straight segments if the neighbor midpoints have the same vertical coordinate. Then, Ω is divided into squares and rectangles. The corresponding quadrangulation is denoted by τh1 (see Figure 1(b)). Similarly, for all Ti,j∈τh, we connect all the neighbor midpoints of the horizontal sides of Ti,j by straight line segments if the neighbor midpoints have the same horizontal coordinate. Then, we obtained the third quadrangulation of Ω, which is denoted by τh2 (see Figure 1(c)).Based on the quadrangulationτh, we define the piecewise constant functional space used to approximate the pressure
(14)Sh:={qh:qh∣T=constant,∀T∈τh;∫Ω‍qhdx=0}.Sh is a subspace of S.Furthermore, using the quadrangulationsτh1 and τh2, we construct a subspace of V. Denote by Γ1, Γ2, Γ3, and Γ4 the south, right, north, and left sides on the boundary of Ω. Set
(15)Vh1={vh∈C(0)(Ω¯):vh∣T1∈Q1,1(T1)∀T1∈τh1,vh=0onΓ2∪Γ4C(0)(Ω¯)},Vh2={vh∈C(0)(Ω¯):vh∣T2∈Q1,1(T2)∀T2∈τh2,vh=0onΓ1∪Γ3C(0)(Ω¯)},
where Q1,1 denotes the piecewise bilinear polynomial space with respect to the variables x and y. Let
(16)Vh=Vh1×Vh2.
Obviously, Vh is a subspace of V.Using the subspacesVh and Sh instead of V and S in the variational Problem (6), we obtain the discrete problem: find (uh,ph)∈Vh×Sh, such that
(17)a(uh,vh)-b(vh,ph)=0∀vh∈Vh,b(uh,qh)=(ϕ,qh)∀qh∈Sh.


## 4. Convergence Analysis and Error Estimate

In this section, we give the corresponding convergence analysis and error estimate. Firstly, we define an interpolating for the following analysis.For the quadrangulationτh, we divided the edges of all squares into two sets. The first one denoted by LV contains all vertical edges, and the second one denoted by LH contains all horizontal edges. We define the interpolation operator Π:V→Vh by Πu=(Πh1u1,Πh2u2)∈Vh1×Vh2, which satisfy the following:
(18)∫l‍Πh1u1ds=∫l‍u1ds∀l∈LV′,∫l‍Πh2u2ds=∫l‍u2ds∀l∈LH′,
where LV′ is a set of edges of elements got by bisecting the most bottom element edges and the most top element edges of LV and LH′ are got by bisecting the most left element edges and the most right element edges of LH. See Figures 2 and 3.Some edges onLV and corresponding edges on LV′.
(a)
(b)Some edges onLH and corresponding edges on LH′.
(a)
(b)Lemma 3.
For anyu∈V, the interpolating Πu∈Vh is uniquely determined by (18).Proof.
It is easy to see that (18) is equivalent to an equation of AX=B, where A is a matrix and X, B are vectors. Direct calculation shows that
(19)A=h*diag{A1,A1,…},
and the form of submatrix A1 is as follows
(20)(1414000⋯00000381800⋯000001834180⋯000000183418⋯0000⋯⋯⋯⋯⋯⋯⋯⋯⋯00000⋯183418000000⋯01838000000⋯001414).
We can see that the matrix is invertible and the equation is solvable, and therefore X can be uniquely determined.Assume that the solution(u,p) of Problem (6) has the following smoothness properties:
(21)u∈V′:=V⋂‍H2((Ω))2,p∈S⋂‍H1(Ω).
Then, we should give the following lemma about the properties of the interpolations defined in (18).Lemma 4.
(i)  There exist two constantsC3 and C4 independent of h, such that
(22)|u-Πu|i,2,Ω≤C3hj-i|u|j,2,Ω,i=0,1,i≤j≤2,(23)infqh∈Sh∥p-qh∥≤C4h|p|1,Ω.
(ii)  There exists a constantC5 independent of h such that
(24)∥Πu∥H(div,Ω)≤C5∥u∥1,Ω∀u∈V.
(iii)  For anyu∈V, we have that
(25)∫Ω‍qhdiv(u-Πu)dx=0,∀qh∈Sh.Proof.
The estimates (22), (23), and (24) follow from Definition (18) and the approximation theory; see [1], for example.
For (25), based on Green formulation, we know that
(26)∫Ω‍qhdiv(u-Πu)dx=∑T∈τh‍∫T‍qhdiv(u-Πu)dx=∑T∈τh‍∫∂T‍qh(u-Πu)·n→ds-∑T∈τh‍∫T‍∇qh·(u-Πu)dx=∑l∈LV‍∫l‍qh(u1-Πh1u1)n1ds+∑l∈LH‍∫l‍qh(u2-Πh2u2)n2ds=∑l∈LV′‍∫l‍qh(u1-Πh1u1)n1ds+∑l∈LH′‍∫l‍qh(u2-Πh2u2)n2ds=0.
Here, n→=(n1,n2), and we use (18) for the last step. The proof is completed.Theorem 5.
The discrete Inf-sup condition is valid; namely, there is a constantβ≥0, such that
(27)supvh∈Vhb(vh,qh)∥vh∥H(div,Ω)≥β∥qh∥,∀qh∈Sh.Proof.
From the process above, we obtain thatb(v,qh)=b(Πv,qh), any v∈V,qh∈Sh. For any ph∈Sh, there exists v∈(H01(Ω))2, such that
(28)∇·v=qh,∥v∥1,Ω≤C6∥qh∥,
where C6 is a constant independent of qh; then we obtain
(29)supvh∈Vhb(vh,qh)∥vh∥H(div,Ω)≥b(Πv,qh)∥Πv∥H(div,Ω)=b(v,qh)∥Πv∥H(div,Ω)=∥qh∥02∥Πv∥H(div,Ω).
Using Lemma 4, we have that
(30)supvh∈Vhb(vh,qh)∥vh∥H(div,Ω)≥1C5∥qh∥02∥v∥1,Ω≥1C5C6∥qh∥.
Taking β=1/C5C6, we complete the proof of (27).With the analysis technique presented by Arbogast and Wheeler [11], we consider the numerical analysis of the mixed finite element presented in this paper. Recall the RT0 mixed element spaces Vh′×Sh′ [3, 5, 6] based on the partition τh(31)Vh′=Q1,0(τh)×Q0,1(τh),Sh′=Sh.Define the interpolation operatorΠ′:V→Vh′ by the following equations:
(32)∫l‍Π′u1ds=∫l‍u1ds∀l∈LV,∫l‍Π′u2ds=∫l‍u2ds∀l∈LH.
Denote by PS:S→Sh the L2 projection operator and by PV′:V→Vh′ the (L2(Ω))2 vector projection operator. The following properties of the projections hold:
(33)∥p-PSp∥0≤Ch|p|1∥u-PV′u∥0≤Ch∥u∥1.
Then, we have an important property about the operator Π′.Lemma 6.
For anyu∈Q1,1(τh1)×Q1,1(τh2), there holds the equivalence Π′u=PV′u; namely,
(34)(u-Π′u,v)=0,∀v∈Vh′.Proof.
As the definition ofVh′ is based on each element T, we focus our discussion on arbitrary element e⊂τh, e=[x0,x0+h]×[y0,y0+h]. Firstly, we consider the x-component (see Figure 4). The analysis for y-component is similar.
For a functionU1∈Vh1, on an element e, it is uniquely given by its node values ui,i=1,…,6. As U1 is a continuous bilinear function on each of the two parts as shown in Figure 4. Then, from (32), we know that Π′u1=a+bx is given by
(35)∫l1‍(a+bx)ds=(a+bx0)*h=(u1+2u3+u5)*h4∫l2‍(a+bx)ds=(a+b(x0+h))*h=(u2+2u4+u6)*h4.
We deduce that
(36)a=u1+2u3+u5-4bx04,b=((u2-u1)+2(u4-u3)+(u6-u5))*14h.
It is clear that we just need to verify (34) for both v=1 and v=x.
We first considerv=1. Denote by φi the node basis function at the point i, which implies that φi(xj)=δi,j, which has the value 1 if and only if i=j; otherwise, it is zero. By direct calculation, we can get the basis, for example,
(37)φ1=14(2-2hx+2hx0)(2-4hy+4hy0),
so
(38)∫e‍U1dxdy=∫x0x0+h‍∫y0y0+h2‍(u1φ1+u2φ2+u3φ3+u4φ4)dxdy+∫x0x0+h‍∫y0+h/2y0+h‍(u3φ3+u4φ4+u5φ5+u6φ6)dxdy=h28(u1+u2+u3+u4)+h28(u3+u4+u5+u6)=h28(u1+u2+2u3+2u4+u5+u6).
By direct computation, we can easily see that∫e‍Π′U1dxdy has the same value, so
(39)∫e‍Π′U1dxdy=∫e‍U1dxdy.
When v=x, we have that
(40)∫e‍Π′U1*xdxdy=∫x0x0+h‍∫y0y0+h‍ax+bx2dxdy=a(x0h2+h32)+b(x02h2+x0h3+13h4),
where a,b are defined in (36). Next, we compare the coefficients of ui in (40) with the coefficients in ∫e‍U1*xdxdy,
(41)∫e‍φ1*xdxdy=∫y0y0+h/2‍∫x0x0+h‍φ*xdxdy=h2∫x0x0+h‍14(2x-2hx2+2hx0x)=124h3+18x0h2=k1,
which determine k1 as the coefficient of u1. With similar computation, we obtain that
(42)k5=k1,k2=k6=x0h28+h312,k3=x0h24+h312,k4=x0h24+h36.
Comparing with (40), we can find that (34) is true with v=x.
So, we certify the lemma.An element onτh  and its corresponding portion on  τh1.
(a)
(b)Theorem 7.
If(u,p) satisfy (6) and (uh,ph) satisfy (17), then there exists a positive constant C independent of h such that the following error estimates hold:
(43)∥u-uh∥0≤Ch∥u∥1,∥p-ph∥0≤Ch(∥u∥1+∥p∥1).Proof.
First, we focus on the erroru-uh. From (6), (17), (18), and (32), we derive that
(44)b(u,qh)=b(uh,qh)=b(Πu,qh)=b(Π′u,qh),∀qh∈Sh.
Let v=Π′vh in (6); then
(45)a(u,Π′vh)-b(Π′vh,p)=0.
Namely,
(46)a(Pv′u,vh)-(Ps∇·vh,p)=a(Pv′u,vh)-b(vh,Psp)=0.
Here, we used the property ∇·Π′v=Ps∇·v. Subtracting from (17), we get that
(47)a(Pv′u-uh,vh)-b(vh,Psp-ph)=0.
Take
(48)vh=Πu-uh,qh=Psp-ph.
Then
(49)a(Pv′u-uh,Πu-uh)-b(Πu-uh,Psp-ph)=0.
Due to (44), we find that
(50)b(Πu-uh,Psp-ph)=0.
Now, we analyze the erroru-uh based on the equations above
(51)a(u-uh,u-uh)=a(u-uh,u-Πu)+a(u-uh,Πu-uh)=a(u-uh,u-Πu)+a(u-Pv′u,Πu-uh)+a(Pv′u-uh,Πu-uh)≤ϵ1∥u-uh∥02+1ϵ1∥u-Πu∥02+ϵ2∥Πu-u∥02+1ϵ2∥u-Pv′u∥02+ϵ3∥u-uh∥02+1ϵ3∥u-Pv′u∥02,
where ϵi>0, i=1,2,3 are positive constants. Take the value of ϵ1=ϵ3=μ/4κ,ϵ2=1, and combining with (22) and (33), we conclude that
(52)∥u-uh∥0≤Ch∥u∥1.
We also can obtain a higher order error estimate for∥Psp-ph∥. Consider the classical duality argument. Let ϕ be the solution of the following elliptical problem:
(53)Δϕ=Psp-ph,∂ϕ∂n=0.
By the elliptic regularity, the estimate holds: |ϕ|H2≤C∥Psp-ph∥0. So
(54)∥Psp-ph∥02=(Psp-ph,∇·∇ϕ)=(Psp-ph,∇·Π∇ϕ)=a(Pv′u-uh,Π∇ϕ)=a(Pv′u-uh,Π∇ϕ-Pv′∇ϕ)+a(Pv′u-uh,Pv′∇ϕ)=a(Pv′u-uh,Π∇ϕ-Pv′∇ϕ)+a(Pv′u-u,Pv′∇ϕ)+a(u-uh,Pv′∇ϕ)=a(Pv′u-uh,Π∇ϕ-Pv′∇ϕ)+a(u-uh,Pv′∇ϕ-∇ϕ)+a(u-uh,∇ϕ).
Now, we estimate the right hand terms of the above inequality. From (33), (22), and (52), we have
(55)a(Pv′u-uh,Π∇ϕ-Pv′∇ϕ)=a(Pv′u-u,Π∇ϕ-∇ϕ)+a(u-uh,Π∇ϕ-∇ϕ)+a(Pv′u-u,∇ϕ-Pv′∇ϕ)+a(u-uh,∇ϕ-Pv′∇ϕ)≤Ch2∥u∥1∣ϕ|H2≤Ch2∥u∥1∣∥Psp-ph∥0.
It is easy to see that(56)a(u-uh,Pv′∇ϕ-∇ϕ)≤Ch2∥u∥1|ϕ|H2≤Ch2∥u∥1∥Psp-ph∥0,a(u-uh,∇ϕ)=a(u-Πu,∇ϕ)+a(Πu-uh,∇ϕ-Pv′∇ϕ)+a(Πu-uh,Pv′∇ϕ)≤C(h2|u|2∥ϕ∥2+h2|u|1∥ϕ∥2)≤Ch2∥u∥2∥Psp-ph∥0.
Here, we used the fact that a(Πu-uh,Pv′∇ϕ)=0 which is got from the Green formulation and (44).
Combining the above inequalities, we conclude that(57)∥p-ph∥0≤∥p-Psp∥0+∥Psp-ph∥0≤Ch(∥u∥1+∥p∥0).
We complete the proof.It is worth mentioning that we analyze this mixed finite element method in a direct way as it is not straightforward to apply the classical inf-sup theory. We just have the coercivity property fora(uh,vh) on the normal L2 space, not in the subspace of v0h={vh∈Vh:b(vh,qh)=0,∀qh∈Sh}, and the same issue also occurs in [11]. The problem is that testing (∇·v,w) by w∈Wh does not control the full divergence of V, and it does not occur when this method is applied to Stokes or Navier-Stokes equations (as in [13, 14]). As a result, we just obtain a convergence rate of ∥u-uh∥0. Failing to obtain convergence rate of ∥u-uu∥H(div,Ω) is a weak point of this proposed mixed formulation compared to the classical Raviart-Thomas mixed method. But the significance of continuous flux applied to mass conservation can be found in Section 6.


## 5. Numerical Examples

In this section, we present some numerical results for the model Problem (1). For simplicity, we assume that the domain is a unit square Ω=[0,1]×[0,1] and the test cases are summarized in Table 1. We can choose the boundary conditions and the right hand terms according to the analytical solutions.Table 1
Three numerical test cases.
Case
Coefficientμ/κ
True solutionu
True solutionp
1
1
(
x
2
y
-
x
4
y
x
y
4
-
x
y
2
)
(
x
-
1
/
2
)
(
y
-
1
/
2
)
2
(
e
2
x
y
2
0
0
1
x
+
y
)
(
x
2
y
-
x
4
y
x
y
4
-
x
y
2
)
(
x
-
1
/
2
)
(
y
-
1
/
2
)
3
(
e
2
x
y
2
0
0
1
x
+
y
)
(
e
-
x
y
x
2
cos
y
)
y
e
xWe compare our method to the formulation constructed by Arbogast and Wheeler [11]. Its corresponding discrete finite element spaces are
(58)V¯h={vh∈(C(0)(Ω¯))2:vh∣T∈Q1,2(T)×Q2,1(T),∀T∈τh(C(0)(Ω¯))2},S¯h={∫Ω‍qhdx=0qh:qh∣T=constant,∀T∈τh;∫Ω‍qhdx=0}.
The results of the error estimate with various norms are listed in Table 2, while the corresponding convergence rates of the presented method are shown in Table 3.Table 2
The numerical error for fm. 1 (our formulation) and fm. 2 (Arbogast's formulation).
Case
Mesh
∥
u
-
u
h
∥
∥
∇
·
(
u
-
u
h
)
∥
∥
p
-
p
h
∥
∥
P
s
p
-
p
h
∥
fm.1
fm.2
fm.1
fm.2
fm.1
fm.2
fm.1
fm.2
1
4
4.90
e
-
2
5.67
e
-
2
3.06
e
-
1
3.24
e
-
1
2.93
e
-
2
2.93
e
-
2
4.53
e
-
3
4.17
e
-
3
8
1.78
e
-
2
2.05
e
-
2
1.53
e
-
1
1.62
e
-
1
1.47
e
-
2
1.47
e
-
2
1.24
e
-
3
1.20
e
-
3
16
6.45
e
-
3
7.37
e
-
3
7.67
e
-
2
8.13
e
-
2
7.37
e
-
3
7.37
e
-
3
3.18
e
-
4
3.15
e
-
4
32
2.31
e
-
3
2.64
e
-
3
3.84
e
-
2
4.07
e
-
2
3.68
e
-
3
3.68
e
-
3
8.01
e
-
5
7.98
e
-
5
64
8.25
e
-
4
9.38
e
-
4
1.92
e
-
2
2.03
e
-
2
1.84
e
-
3
1.84
e
-
3
2.01
e
-
5
2.01
e
-
5
2
4
4.70
e
-
2
5.47
e
-
2
2.99
e
-
1
3.22
e
-
1
2.95
e
-
2
2.94
e
-
2
5.42
e
-
3
4.97
e
-
3
8
1.72
e
-
2
1.99
e
-
2
1.53
e
-
1
1.63
e
-
1
1.47
e
-
2
1.47
e
-
2
1.54
e
-
3
1.48
e
-
3
16
6.25
e
-
3
7.19
e
-
3
7.75
e
-
2
8.27
e
-
2
7.37
e
-
3
7.37
e
-
3
4.04
e
-
4
3.98
e
-
4
32
2.25
e
-
3
2.58
e
-
3
3.89
e
-
2
4.15
e
-
2
3.68
e
-
3
3.68
e
-
3
1.03
e
-
4
1.02
e
-
4
64
8.08
e
-
4
9.21
e
-
4
1.95
e
-
2
2.08
e
-
2
1.84
e
-
3
1.84
e
-
3
2.59
e
-
5
2.58
e
-
5
3
4
9.65
e
-
2
1.09
e
-
1
4.14
e
-
1
4.67
e
-
1
1.49
e
-
1
1.49
e
-
1
7.39
e
-
3
6.21
e
-
3
8
3.79
e
-
2
4.31
e
-
2
2.16
e
-
1
2.46
e
-
1
7.44
e
-
2
7.44
e
-
2
2.14
e
-
3
1.89
e
-
3
16
1.42
e
-
2
1.62
e
-
2
1.11
e
-
1
1.28
e
-
1
3.72
e
-
2
3.72
e
-
2
5.72
e
-
4
5.19
e
-
4
32
5.19
e
-
3
5.91
e
-
3
5.63
e
-
2
6.51
e
-
2
1.86
e
-
2
1.86
e
-
2
1.47
e
-
4
1.35
e
-
4
64
1.87
e
-
3
2.13
e
-
3
2.84
e
-
2
3.28
e
-
2
9.31
e
-
3
9.31
e
-
3
3.72
e
-
5
3.44
e
-
5Table 3
The corresponding convergence rates of fm. 1 and fm. 2.
Case
Mesh
∥
u
-
u
h
∥
∥
∇
·
(
u
-
u
h
)
∥
∥
p
-
p
h
∥
∥
P
s
p
-
p
h
∥
fm.1
fm.2
fm.1
fm.2
fm.1
fm.2
fm.1
fm.2
1
8
1.459
1.468
0.997
1.001
0.995
0.993
1.875
1.795
16
1.468
1.476
0.998
0.995
0.999
0.999
1.961
1.934
32
1.479
1.484
1.000
0.998
1.000
1.000
1.987
1.978
64
1.486
1.489
1.000
1.000
1.000
1.000
1.996
1.993
2
8
1.449
1.457
0.968
0.978
0.999
0.996
1.817
1.742
16
1.462
1.471
0.983
0.984
1.001
1.001
1.930
1.901
32
1.471
1.479
0.993
0.993
1.001
1.000
1.976
1.960
64
1.480
1.485
0.997
0.997
1.000
1.000
1.989
1.984
3
8
1.347
1.340
0.942
0.924
0.998
0.997
1.787
1.708
16
1.416
1.414
0.957
0.945
0.999
0.999
1.906
1.870
32
1.452
1.452
0.979
0.973
1.000
1.000
1.959
1.942
64
1.472
1.473
0.990
0.986
1.000
1.000
1.983
1.975Close results of numerical errors for both formulations are shown in Table2. From Table 3, we can see that p converges to ph as O(h) and Psp-ph as O(h2) for our formulation, which both agree with the theorem. From the examples, we can observe that uh converges to u somewhat better than expected, and it appears that on the uniform grid we attain O(h3/2) superconvergence in the L2 norm which is similar to the tests of Arbogast’s formulation [11]. Yet, the degrees of freedom of our method are less than Arbogast’s scheme. As in the case of 64*64, the degrees of freedom of Arbogast’s scheme are 20866 and 12676 for our formulation. The convergence rate of ∥u-uh∥H(div,Ω) is first order, but here we cannot give the corresponding analysis.


## 6. A Valuable Application

In this section, we briefly show an application of the proposed mixed finite element method to the miscible displacement of one incompressible fluid by another in porous media. The model is as follows:(59)μ(C)K-1u+∇p=γ(C)∇d,(x,t)∈Ω×J,ϕ∂C∂t+∇·(uC)-∇·(D(u)∇C)=C~q,(x,t)∈Ω×J,∇·u=g,(x,t)∈Ω×J,u·n=g1,(x,t)∈∂Ω×J,C(x,0)=C0(x),x∈Ω,
where γ(C) and d are the gravity coefficient and vertical coordinate, ϕ(x) is the porosity of the rock, and C~q represents a known source. D(x,u) is the molecular diffusion and mechanical dispersion coefficient. For convenience, we denote that f=C~q and a(C)=μ(C)K-1. Let χ:(0,T]→R2 be the solution of the ordinary differential equation
(60)dχdτ=u(χ(x,t;τ),τ)ϕ(x),χ(x,t;t)=x.
Let V=H(div,Ω), S=L02(Ω), M=H1(Ω); then, we derive the entire weak formulation for the model: find (u,p,C)∈V×S×M, such that
(61)(a(C)u,v)-(p,∇·v)=(γ(C)∇d,v),∀v∈V,(ϕ(x)dC(χ,τ)dτ+gC,w)+(D∇C,∇w)=(f,w),∀w∈M,(∇·u,φ)=(g,φ),∀φ∈S.
Let Δt be the time step for both concentration and pressure; define
(62)Mh={vh∈C(0)(Ω¯):vh∣T∈Q1,1(T),∀T∈τh}.
Combing with the new characteristic finite element method which preserves the mass balance proposed by Rui and Tabata [10], the approximate characteristic line of χ is defined as
(63)χn(x)=x-uhnϕ(x)Δt.
We obtain the corresponding full-discrete mass-conservative characteristic (MCC) scheme: find (uh,ph,Ch)∈Vh×Sh×Mh, such that
(64)(a(Chn)uhn,vh)-(ph,∇·vh)=(γ(Chn)∇d,vh),∀vh∈Vh(ϕChn-(ϕChn-1)∘χnγnΔt,φh)+(D(uhn)∇Chn,∇φh)=(f,φh),∀φh∈Mh(∇·uhn,qh)=(g,qh),∀qh∈ShCh0=C~0,
where
(65)γn=det(∂χn∂x)=1-∇·uhnϕΔt+uhn∇ϕϕ2Δt+∇(uh,1nϕ)·curl(uh,2nϕ)Δt2.
We can see that the continuous flux is indispensable for γn. Let φh=1 in (64), and summing it up from n=1 to N, we get the mass balance
(66)∫Ω‍ϕChNdx=∫Ω‍ϕCh0dx+Δt∑n=1N‍∫Ω‍fndx.
Here, we just give numerical example to show the feasibility of this application, and the theoretical analysis of stability, mass balance, and convergence of this discrete scheme will be discussed in the future. Firstly, we define compute mass error and relative mass error as follows:
(67)computemasserror:∫Ω‍ϕChNdx-(∫Ω‍ϕCh0dx+Δt∑n=1N‍∫Ω‍fndx),relativemasserror:∫Ω‍ϕChNdx-∫Ω‍ϕCNdx∫Ω‍ϕCNdx.
Now, we select μ(C)=C, and the following analytical solution of the problem is
(68)u(x,y,t)=(ex+t,ey+t),p(x,y,t)=e-t(x2+y2),C(x,y,t)=e-t((x-12)2+(y-12)2).
The error results with different norms of this numerical simulation can be listed in Tables 4 and 5, and at last we give a mass error to check the mass conservation in Table 6.Table 4
Numerical error and convergence rate(Δt=Ch).
Mesh
5
×
5
10
×
10
20
×
20
40
×
40
Norm type
Error
Rate
Error
Rate
Error
Rate
Error
Rate
∥
u
∥
l
2
(
L
2
)
1.83
e
-
4
—
7.10
e
-
5
1.36
3.38
e
-
5
1.07
1.65
e
-
5
1.03
∥
u
∥
l
∞
(
L
2
)
1.29
e
-
2
—
5.19
e
-
3
1.31
2.64
e
-
3
0.97
1.37
e
-
3
0.95
∥
p
∥
l
2
(
L
2
)
1.33
e
-
3
—
6.67
e
-
4
1.00
3.33
e
-
4
1.00
1.67
e
-
4
1.00
∥
p
∥
l
∞
(
L
2
)
9.43
e
-
2
—
4.71
e
-
2
1.00
2.35
e
-
2
1.00
1.18
e
-
2
1.00
∥
C
∥
l
2
(
H
1
)
2.32
e
-
3
—
1.16
e
-
3
1.01
5.78
e
-
4
1.00
2.88
e
-
4
1.00
∥
C
∥
l
∞
(
H
1
)
1.63
e
-
1
—
8.18
e
-
2
1.00
4.11
e
-
2
0.99
2.05
e
-
2
0.99Table 5
Numerical error and convergence rate(Δt=Ch2).
Mesh
5
×
5
10
×
10
20
×
20
40
×
40
Norm type
Error
Rate
Error
Rate
Error
Rate
Error
Rate
∥
C
∥
l
2
(
L
2
)
8.48
e
-
5
—
2.13
e
-
5
1.995
5.37
e
-
6
1.986
1.36
e
-
6
1.971
∥
C
∥
l
∞
(
L
2
)
1.34
e
-
2
—
3.37
e
-
3
1.989
8.56
e
-
4
1.978
2.21
e
-
4
1.952Table 6
Mass error for concentrationC(Δt=Ch).
Mesh
5
×
5
10
×
10
20
×
20
40
×
40
Compute  mass error
1.209
e
-
3
1.243
e
-
3
1.269
e
-
3
1.284
e
-
3
Relative mass error
2.068
e
-
2
5.427
e
-
3
1.487
e
-
3
4.371
e
-
4As can be seen from Tables4 and 5, we conjecture that almost all the convergence rates are true in general. From Table 6 we find that mass balance is right as computational mass error resulting from computer is inevitable and nearly invariable for different meshes, while the relative mass error decreases as was expected. The corresponding theoretical analysis about this system will be considered in the future work.



---

*Source: 580461-2013-06-16.xml*