# A BrokenP1-Nonconforming Finite Element Method for Incompressible Miscible Displacement Problem in Porous Media

Authors: Fengxin Chen; Huanzhen Chen
Journal: ISRN Applied Mathematics
 (2013)

---

## Abstract

An approximate scheme is defined for incompressible miscible displacement in porous media. This scheme is constructed by using immersed interface finite element method for the pressure equation which is based on the brokenP1-nonconforming piecewise linear polynomials on interface triangular elements and utilizing finite element method for the concentration equation. Error estimates for pressure in broken H1 norm and for concentration in L2 norm are presented.


---

## Body

## 1. Introduction

Miscible displacement of one incompressible fluid in a porous mediumΩ over time interval J=[0,T] is modeled by the system
(1)(a)-∇·(α(c)∇p)=∇·u=q,x∈Ω,t∈J,(b)ϕ∂c∂t+u·∇c-∇·(D(u)∇c)=(c~-c)q=g(x,t,c),x∈Ω,t∈J,(c)u·n=D(u)∇c·n=0,x∈∂Ω,t∈J,(d)c(x,0)=c0(x),x∈Ω.In the formulation, we have utilized the fact that in applications the vertical dimension of a subsurface geological formation is often much smaller than the horizontal dimension. We thus simplify the problem via a vertical average to rewrite the problem as a two-dimensional problem in the horizontal dimension, so the physical domainΩ⊂R2. α(c)=α(x,c)=k(x)/μ(c), k(x) is the permeability tensor of the medium and μ(c) is the viscosity of the fluid mixture that may be discontinuous across some interfaces; ϕ(x) is the porosity of the medium; u(x,t) is the Darcy velocity of the mixture; q(x,t) represents flow rates at wells, commonly a linear combination of Dirac measures; D(u)=D(x,c,∇p)=ϕ(x){dmI+|u|(dlE(u)+dtET(u))} is the diffusion-dispersion tensor, with dm, dl, and dt being the molecular diffusion, the transverse, and longitudinal dispersivities, respectively, I is the identity tensor; c~(x,t) is specified at sources and c~(x,t)=c(x,t) at sinks. c0(x,t) is the initial concentration. The dependent variable is p(x,t) is the pressure in the fluid mixture, and c(x,t) is the concentration of a solvent injected into resident reservoir.For (1) we assume that
(2)(q,1)=∫Ω‍q(x,t)dx=0,t∈J.Letu=u(x,c,∇p)=(u1(x,c,px),u2(x,c,py)); for some ε^>0 restrict the variable q1 to lie between -ε^⩽q1⩽1+ε^ and q2∈R, so we can suppose that the following regularity for α, D, ϕ, and u holds:
(3)0<α*⩽α(x,q1)⩽α*,D*⩽D(u)⩽D*,ϕ*⩽ϕ(x)⩽ϕ*;(4)|ui(x,q1,q2)|⩽K(1+|q2|),i=1,2,q2∈R.There exist much literature concerning numerical method and numerical analysis of miscible placement problem in porous media, for example, [1–4], and so on. As we know the porous media equations used to model the interface between oil and an injected fluid in simulations of secondary recovery in oil reservoirs. For the pressure equation (1)(a), the coefficient α(c) often changes rapidly across fluid interfaces, and this sharp change is accompanied by large changes in the pressure gradient, as a compensatory, yielding a fairly smooth Darcy velocity u, so we can trade (1)(a) as an interface problem. When the interface is smooth enough, the solution of the interface problem is also very smooth in individual regions where the coefficient is smooth, but due to the jump of the coefficient across the interface, the global regularity is usually very low and has order of H1+γ, 0<γ<1. Due to the low global regularity and the irregular geometry of the interface, it seems to be difficult for the standard finite element method to achieve high accuracy.For better approximation, the fitted finite element methods whose mesh depends on the smooth interface are developed [5]. However, this method using fitted grids is costly for more complicated time dependent problems in which the interface moves with time and repeated grid generation is called for. Compared with the fitted finite element methods, the immersed interface method proposed by LeVeque and Li [6] allows the mesh to be independent of the interface, such as a Cartesian mesh. In recent years, Li et al. [7] studied an immersed finite element method using uniform grid, and they proved the approximation property of the finite element space of this scheme. On the other hand, Kwak et al. [8] introduced an immersed finite element method based on the broken P1-nonconforming piecewise linear polynomials on interface triangular elements; this method uses edge averages as degrees of freedom, and the basis functions are C-R type [9]. Theory analysis and numerical experiments also show the optimal-order convergence of the method.In this paper, we applyP1-nonconforming finite element method (P1FEM) to pressure equation, while a FEM was used to approximate the concentration equation. Other methods could also be employed for the discretization of the concentration equation, for example, characteristic Galerkin method and so on. However, since the main point of this paper is to show the feasibility of the use of P1FEM for pressure, we will discuss the concentration equation in the single case.The rest of the paper is organized as follows: in Section2, we first introduce some preliminaries. In Section 3, we briefly describe the P1FEM-FEM scheme. In Section 4, we present some projections and lemmas used in the following error analysis. In Section 5, we prove the main error estimate.


## 2. Preliminaries

In this paper, without loss of generality, we consider the case in whichΩ⊂R2 is a rectangular domain and the interface Γ is a smooth curve separating Ω into two subdomains Ω+­ and Ω- such that Ω¯=Ω+¯⋃Ω-¯⋃Γ; see Figure 1 for an illustration.Figure 1The proper jump conditions for elliptic (1)(a) are given by
(5)[p]=0,[α(c)∇p·n]=0acrossΓ,
where p∈H1(Ω). We assume α(c) is a positive function bounded below and above by two positive constants.LetWqm(Ω), 1⩽q⩽+∞, be the Sobolev spaces consisting of functions whose derivatives up to order-m are qth integrable on Ω, and Hm(Ω):=W2m(Ω) and H0(Ω)=L2(Ω).For the analysis, we introduce the space(6)H~2(Ω)={p∈H1(Ω);p∈H2(Ωs),s=+,-},
equipped with the norm
(7)∥p∥H~2(Ω)2=∥p∥H2(Ω+)2+∥p∥H2(Ω-)2,|p|H~2(Ω)2=|p|H2(Ω+)2+|p|H2(Ω-)2,111111111111111∀p∈H~2(Ω).We integrate (1)(a) multiplied by any test function v∈H1(Ω) over the domain Ω and apply the divergence theorem; then the variational formulation of the pressure equation is as follows Find p∈H1(Ω) such that
(8)a(p,v)=(q,v),v∈H1(Ω),
where the bilinear formation a(p,v)=(α(c)∇p,∇v).By the Sobolev embedding theorem, for anyp∈H~2(Ω), we can get p∈Ws1(Ω) for any s>2. Then we have the following regularity lemma for the weak solution p of the variational formulation (8).Lemma 1 (see [10]).
Assume thatq∈L2(Ω). Then there exists a unique solution p∈H~2(Ω) to the variational formulation (8) such that
(9)∥p∥H~2(Ω)⩽∥q∥L2(Ω).


## 3. Formation of the Method

In this section, firstly, we define a brokenP1-nonconforming finite element method for pressure equation, and then, we use a finite element method to approximate the concentration equation.In order to construct a brokenP1-nonconforming finite element procedure for pressure equation (1)(a), we assume the following situation.LetTh be the usual regular triangulation of the domain Ω such that the elements have diameters bounded by hp. Without loss of generality, we assume that the triangles in the partition used have the following features.(H1)
IfΓ meets one edge at more than two points, then this edge is one part of Γ.(H2)
IfΓ meets a triangle at two points, then the two points must be on the different edges of this triangle.(H3)
The interface curveΓ is defined by a piecewise C2 function, and the mesh Th is formed such that the subset of Γ in any interface element is C2.Then we can separate the triangles of partitionTh into two classes. For an element T∈Th, if the interface Γ passes through the interior of T, we call it an interface element and denote it by Tm; otherwise, we call it a noninterface element and denote it by Tn, respectively.As usual, we will define linear finite element spaces on each element of the partitionTh, so we have to construct the local basis functions. For a noninterface element Tn, we can simply use the standard linear shape functions on Tn with degrees of freedom at average values along edges ej of Tn and construct the linear finite element space S¯h(Tn) as follows:
(10)S¯h(Tn)=span{1|ej|∫ej‍φi:φi∈P1(Tn),hhhhhh1|ej|∫ej‍φids=δij,i,j=1,2,3}.For this space, we define the interpolation operatorIh:H2(Tn)→S¯h(Tn), where the following well-known approximation property is satisfied:
(11)∥p-Ihp∥L2(Tn)+hp∥p-Ihp∥H1(Tn)⩽Chp2∥p∥H2(Tn),
and we use S¯h(Ωn) to denote the space of piecewise P1-nonconforming finite element space on the domain Ωn=⋃‍Tn.For interface elementTm, since the finite element space on a general element can be obtained from the counterpart of a reference element through an affine mapping, we consider a typical interface element Tm whose geometric configuration is given in Figure 2 in which the three vertices are given by A1=(0,0), A2=(1,0), and A3=(0,1), and the curve between points D and E is part of the interface with D=(0,a) for 0<a⩽1 and E=(b,0) for 0<b⩽1. Let DE¯ denote the line segment connecting points D and E, which divides Tm into two parts Tm+ and Tm- with Tm=Tm+⋃Tm-⋃DE¯.Figure 2Letei, i=1,2,3 be the edges of Tm, and let φ¯ei denote the average of φ along ei, that is, φ¯ei=(1/|ei|)∫ei‍φds, φ∈H1(Tm), and then we would like to construct a new function which is linear on Tm+ and Tm-, respectively, and satisfies the jump condition (5) on DE¯. For this purpose, we write the modified basis function φ^ on an interface element Tm as follows:
(12)φ^={φ^-=a0+b0x+c0y,inTm-,φ^+=a1+b1x+c1y,inTm+,
with the following constraints:
(13)φ^¯ei=Vi,i=1,2,3,φ^-(D)=φ^+(D),φ^-(E)=φ^+(E),α-¯∂φ^-∂nDE¯=α+¯∂φ^+∂nDE¯,
where α-¯, α+¯ are averages along DE¯ and ρ=α-¯/α+¯, Vi, i=1,2,3, are given values, and nDE¯ is the unit normal vector on the line DE¯, that is, nDE¯=(a,b)/a2+b2.By similar calculation to that given in Theorem 2.2 of [8], we can rewrite (12) and (13) in the matrix form:
(14)A(a0b0c0a1b1c1)=(V1V2V3000),
where the coefficient matrix is defined by
(15)A=(00011212a012a21-a012(1-a2)b12b201-b12(1-b2)010a-10-ah21b0-1-b00ρaρb0-a-b),
and the determinant of the matrix A is
(16)det(A)=-14a3b+14a3bρ-14a2ρ-14ab3+14ab3ρ-14b2ρ=14(a2+b2){ρ(ab-1)-ab}<0.Therefore, the coefficients of (12) are uniquely determined by conditions φ^¯ei=Vi, i=1,2,3, respectively, and when Vi, i=1,2,3, have the same value V, we can get the piecewise linear function φ^=V by uniqueness.Therefore we can define the following finite element space on an interface elementTm:
(17)S^h(Tm)=span{φ^:φ^|Tmiswelldefinedbytheaboveconstruction}.Also we useS^h(Ωm) to denote the finite element space defined on the domain Ωm=⋃‍Tm.Remark 2.
Although for functions inS^h(Tm), the flux jump condition is enforced on line segments, they actually satisfy a weak flux jump condition along the interface Γ when α is a piecewise constant such that
(18)∫Γ∩T‍(α-∇φ--α+∇φ+)·nΓds=0,∀φ∈S^h(Tm),
which is proved by an application of divergence theorem as in [11].Combining the definitions ofS¯h(Ωn) and S^h(Ωm),we can describe the immersed finite element space S~h(Ω) on the whole domain Ω,
(19)S~h(Ω)={φ~:φ~|Ωn∈S~h(Ωn),φ~|Ωm∈S^h(Ωm)},
endowed with the broken norm,
(20)∥φ~∥1,h2=∑T‍∫T‍|∇φ~|2dx,φ~∈S~h(Ω).Now, assume that the approximation of concentrationC is known. Then, the pressure P∈S~h(Ω) can be determined by the following system:
(21)ah(P,φ~)=(q,φ~),φ~∈S~h(Ω),
where ah(P,φ~)=∑T∈Th∫T‍α(C)∇P∇φ~dx is the bilinear formulation defined on Hh(Ω)×Hh(Ω) with Hh(Ω)=H1(Ω)⊕S~h(Ω).The question at hand is to discretize the concentration equation.Multiplying the concentration equation (1)(b) by a test function v∈H1(Ω) and integrating over Ω which leads to the following weak formulation for concentration c:
(22)(ϕ∂c∂t,v)+(u(c,∇p)·∇c,v)+(D(c,∇p)∇c,∇v)=(g(c),v),v∈H1(Ω).LetMh denote the standard piecewise linear finite element space associated with Th, and hc denote the bound of elements diameters. Then the discrete procedure of (22) is to find C:J→Mh such that
(23)(ϕ∂C∂t,vh)+(u(C,∇P)·∇C,vh)+(D(C,∇P)∇C,∇vh)=(g(C),vh),vh∈Mh,C(x,0)=c~0(x),x∈Ω,
where c~0(x) is the elliptic projection of c0(x).


## 4. Projection and Some Lemmas

For convergence analysis conducted in the following section, we need an interpolation operatorIh:H~2(T)→S~h(T) using the average of p on each edge by
(24)(Ihp¯)ei=p¯ei,i=1,2,3,
and when p∈H~2(Ω), we can define Ihp by (Ihp)|T=Ih(p|T). Then the following interpolation estimate hold for the P1FEM approximation [8]
(25)∥p-Ihp∥L2(Ω)+hp∥p-Ihp∥1,h⩽Chp2∥p∥H~2(Ω),∀p∈H~2(Ω).Next, letC~:J→Mh be the elliptic projection of c defined by
(26)(u(c,∇p)·∇(C~-c),vh)+(D(c,∇p)∇(C~-c),∇vh)+(λ(C~-c),vh)=0,∀vh∈Mh.The functionλ will be chosen to assure coercivity of the form. Then, the following estimates hold [12]:
(27)∥c-C~∥L2(Ω)+hc∥c-C~∥H1(Ω)⩽Khc2∥c∥H2(Ω),∥∂(c-C~)∂t∥L2(Ω)⩽Khc2{∥c∥H2(Ω)+∥∂c∂t∥H2(Ω)}.For convergence analysis conducted in Section5, we need to apply quote some lemmas from [9, 13] and references therein.Lemma 3 (see [9] (the second Strang lemma)).
Letp and P be the solution of (8) and (21). Then there exists a constant C>0 such that
(28)∥p-P∥1,h⩽C{infwh∈S~h(Ω)∥p-wh∥1,h+supφ~∈S~h(Ω)|ah(p,φ~)-(q,φ~)|∥φ~∥1,h}.Lemma 4 (see [13]).
Lete be an edge of  T. Then there exists a constant C>0 such that for all φ,v∈H1(T):
(29)|∫e‍φ(v-v¯e)ds|⩽Chp|φ|1,T|v|1,T,
where v¯e:=(1/|e|)∫e‍vds.


## 5. Convergence Analysis

In this section, we will present convergence analysis for the pressure and concentration.Theorem 5.
Letp∈H~2(Ω), P∈S~h(Ω) be the solution of (8) and (21); then there exists a constant K>0 independent of hp and the location of the interface, such that the following result holds:
(30)∥p-P∥1,h⩽Khp∥p∥H~2(Ω)+K∥c-C∥L2(Ω).Proof.
Since the immersed finite element formulation (21) is nonconforming, we can use Lemma 3 to prove the error bound. The first part in (28) is an approximation error, which can be estimated simply:
(31)infwh∈S~h(Ω)∥p-wh∥1,h⩽Khp∥p∥H~2(Ω).
For the second part, it is a consistency error estimate. By the definition ofah(·,·) and Green’s formula, we can get
(32)|ah(p,φ~)-(q,φ~)|=|∑T∈Th‍∫T‍α(C)∇p∇φ~dx-∫Ω‍qφ~dx|=|∑T∈Th‍∫T‍α(C)∇p∇φ~dxffffffsssffff-∑T∈Th‍∫T‍α(c)∇p∇φ~dxffffffffffsss+∑T∈Th‍〈α(c)∂p∂n,φ~〉∂T|⩽|∑T∈Th‍∫T‍(α(C)-α(c))∇p∇φ~dx|ffffffffff+|∑T∈Th‍〈α(c)∂p∂n,φ~〉∂T|⩽∥∇p∥L∞(Ω)∥c-C∥L2(Ω)∥φ~∥1,hffffffffff+|∑T∈Th‍〈α(c)∂p∂n,φ~〉∂T|⩽K∥c-C∥L2(Ω)∥φ~∥1,hffffffffff+|∑T∈Th‍〈α(c)∂p∂n,φ~〉∂T|,
where n is the unit outward normal vector on each ∂T, φ~∈S~h(Ω), and by the construction of the space S~h(Ω) we have a well-defined property on the interior edges, that is, ∫e‍φ~|T1=∫e‍φ~|T2 for e is the common edge of adjacent element T1 and T2. By boundary condition, we note that α(c)(∂p/∂n) vanishing average on the boundary. So we rewrite ∑T∈Th<α(c)(∂p/∂n),φ~>∂T and use Lemma 4 to derive that
(33)∑T∈Th‍〈α(c)∂p∂n,φ~〉∂T=∑T∈Th‍∑e∈∂T‍〈α(c)∂p∂n-α(c)∂p∂n¯,φ~〉e⩽∑T∈Th‍Chp|α(c)∂p∂n|1,T|φ~|1,T⩽Chp∥p∥H~2(Ω)∥φ~∥1,h.
Combining (28) with (31)–(33) yields the theorem result.We now are in the position to prove the error estimate for concentration.Theorem 6.
Assume that the true solution(c,p) of (1) satisfies c∈L∞(H2)∩L2(H2), p∈L2(H~2). Let C~ be the projection of c satisfying (26) and let C be the solution of (23). Suppose that the mesh parameters satisfy hc2=O(hp). Then the following result holds:
(34)∥C-C~∥L∞(L2)+∥∇(C-C~)∥L2(L2)⩽K(hc2+hp).Proof.
We decompose the errorc-C to ξ=C-C~, and η=c-C~. Then, combining (22), (23), and (26) results in the following error equation for the concentration:
(35)(ϕ∂ξ∂t,vh)+(D(C,∇P)∇ξ,∇vh)=(ϕ∂η∂t,vh)-(λη,vh)+(g(C)-g(c),vh)-(u(C,∇P)∇ξ,vh)+((D(c,∇p)-D(C,∇P))∇C~,∇vh)+((u(c,∇p)-u(C,∇P))∇C~,vh).
We choose the test function asvh=ξ in (35) and integrate from 0 to t on both sides to derive
(36)∫0t‍(ϕ∂ξ∂t,ξ)dτ+∫0t‍(D(C,∇P)∇ξ,∇ξ)dτ=∫0t‍(ϕ∂η∂t,ξ)dτ-∫0t‍(λη,ξ)dτ+∫0t‍(g(C)-g(c),ξ)dτ-∫0t‍(u(C,∇P)∇ξ,ξ)dτ+∫0t‍((D(c,∇p)-D(C,∇P))∇C~,∇ξ)dτ+∫0t‍((u(c,∇p)-u(C,∇P))∇C~,ξ)dτ.
The terms on the left side can be easily bounded by(37)∫0t‍(ϕ∂ξ∂t,ξ)dτ+∫0t‍(D(C,∇P)∇ξ,∇ξ)dτ⩾12ϕ*∥ξ∥2+D*∫0t‍∥∇ξ∥2dτ.
Next, we should estimate terms on the right side one by one.
By (27), it is trivial that
(38)|∫0t‍(ϕ∂η∂t,ξ)dτ|⩽K(hc4+∫0t‍∥ξ∥2dτ).
Similarly, we can bound the second and third terms as follows:(39)|∫0t‍(λη,ξ)dτ|⩽K(hc4+∫0t‍∥ξ∥2dτ),∫0t‍(g(C)-g(c),ξ)dτ⩽K(hc4+∫0t‍∥ξ∥2dτ),
provided that g is Lipschitz continuous.
For the fifth term, we can derive the estimate by using the boundness ofC~ and (27), as well as (30),
(40)|∫0t‍((D(c,∇p)-D(C,∇P))∇C~,∇ξ)dτ|⩽K∥∇C~∥L∞(L∞)∫0t‍(∥ξ∥+∥η∥+∥p-P∥1,h)∥∇ξ∥dτ⩽K∫0t‍(∥ξ∥+hc2+hp)∥∇ξ∥dτ⩽K(∫0t‍∥ξ∥2dτ+hc4+hp2)+D*4∫0t‍∥∇ξ∥2dτ.
Similarly, we can derive(41)|∫0t‍((u(c,∇p)-u(C,∇P))∇C~,ξ)dτ|⩽K(hc4+hp2+∫0t‍∥ξ∥2dτ).
In order to prove the forth term, we need the following induction hypothesis:(42)∥∇P∥L∞(L∞)⩽K*.
As the statement in [1], we assume that for some ε>0,ε⩽ε^,
(43)-ε≤C⩽1+ε
and with this hypothesis and (4), we can get
(44)|∫0t‍(u(C,∇P)∇ξ,ξ)dτ|⩽K(∥∇P∥L∞(L∞)+1)2fffff×∫0t‍∥ξ∥2dτ+D*4∫0t‍∥∇ξ∥2dτ⩽K∫0t‍∥ξ∥2dτ+D*4∫0t‍∥∇ξ∥2dτ.
Therefore, collecting all the bounds derived above and using Gronwall inequality leads to(45)∥ξ∥L∞(L2)2+∥∇ξ∥L2(L2)2⩽K(hc4+hp2).
To complete the argument, we have to check hypothesis (42). We use (25), (27), (30), (45), inverse inequality, and boundness of Ihp to get that
(46)∥∇P∥L∞(L∞)⩽∥∇(P-Ihp)∥L∞(L∞(T))∥∇Ihp∥L∞(L∞(T))⩽K+Khp-1∑T∈Th‍∥∇(P-Ihp)∥L∞(L2(T))⩽K+Khp-1{∥P-p∥L∞(1,h)+∥p-Ihp∥L∞(1,h)}⩽K+Khp-1{hp+hc2}⩽K*.Combining (27) with (30) and (45), we deduce the error estimates for the pressure and the concentration.Theorem 7.
Assume that the true solution(c,p) of (1) satisfies c∈L∞(H2)∩L2(H2), p∈L2(H~2). Let C, P be the solution of (23) and (21). Suppose that the mesh parameters satisfy hc2=O(hp). Then the following result holds:
(47)∥p-P∥L∞(1,h)+∥c-C∥L∞(L2)⩽K(hc2+hp),∥∇(c-C)∥L2(L2)⩽K(hc+hp).Remark 8.
In this paper, we just get a priori error estimates for the coefficientD=D(x,c,∇p). However, some mathematical models for miscible displacements, which are currently being used by oil companies, make the assumption that the coefficient D=D(x,c) and it does not depend on ∇p; with this assumption, we can prove the optimal order convergence similarly.



---

*Source: 498383-2013-12-11.xml*