A Broken P 1-Nonconforming Finite Element Method for Incompressible Miscible Displacement Problem in Porous Media

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 broken P 1 -nonconforming piecewise linear polynomials on interface triangular elements and utilizing finite element method for the concentration equation. Error estimates for pressure in brokenH norm and for concentration in L2 norm are presented.


Introduction
Miscible displacement of one incompressible fluid in a porous medium Ω over time interval  = [0, ] is modeled by the system ( 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 Ω ⊂  2 .() = (, ) = ()/(), () is the permeability tensor of the medium and () is the viscosity of the fluid mixture that may be discontinuous across some interfaces; () is the porosity of the medium; (, ) is the Darcy velocity of the mixture; (, ) represents flow rates at wells, commonly a linear combination of Dirac measures; () = (, , ∇) = (){   + ||(  () +     ())} is the diffusion-dispersion tensor, with   ,   , and   being the molecular diffusion, the transverse, and longitudinal dispersivities, respectively,  is the identity tensor; c(, ) is specified at sources and c(, ) = (, ) at sinks. 0 (, ) is the initial concentration.The dependent variable is (, ) is the pressure in the fluid mixture, and (, ) is the concentration of a solvent injected into resident reservoir.
For (1) we assume that (, 1) = ∫ Ω  (, )  = 0,  ∈ . ( Let  = (, , ∇) = ( 1 (, ,   ),  2 (, ,   )); for some ε > 0 restrict the variable  1 to lie between −ε ⩽  1 ⩽ 1 + ε and  2 ∈ , so we can suppose that the following regularity for , , , and  holds: There exist much literature concerning numerical method and numerical analysis of miscible placement problem in porous media, for example, [1][2][3][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 () 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 , 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  1+ , 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  1 -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 apply  1 -nonconforming finite element method ( 1 FEM) 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  1 FEM for pressure, we will discuss the concentration equation in the single case.
The rest of the paper is organized as follows: in Section 2, we first introduce some preliminaries.In Section 3, we briefly describe the  1 FEM-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.

Preliminaries
In this paper, without loss of generality, we consider the case in which Ω ⊂  2 is a rectangular domain and the interface Γ is a smooth curve separating Ω into two subdomains Ω + and Ω − such that Ω = Ω + ⋃ Ω − ⋃ Γ; see Figure 1 for an illustration.
By the Sobolev embedding theorem, for any  ∈ H2 (Ω), we can get  ∈  1  (Ω) for any  > 2. Then we have the following regularity lemma for the weak solution  of the variational formulation (8).

Formation of the Method
In this section, firstly, we define a broken  1 -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 broken  1 -nonconforming finite element procedure for pressure equation (1)(a), we assume the following situation.
Let  ℎ be the usual regular triangulation of the domain Ω such that the elements have diameters bounded by ℎ  .Without loss of generality, we assume that the triangles in the partition used have the following features.
( 1 ) If Γ meets one edge at more than two points, then this edge is one part of Γ. ( 2 ) If Γ meets a triangle at two points, then the two points must be on the different edges of this triangle.
( 3 ) The interface curve Γ is defined by a piecewise  2 function, and the mesh  ℎ is formed such that the subset of Γ in any interface element is  2 .
Then we can separate the triangles of partition  ℎ into two classes.For an element  ∈  ℎ , if the interface Γ passes through the interior of , we call it an interface element and denote it by   ; otherwise, we call it a noninterface element and denote it by   , respectively.
As usual, we will define linear finite element spaces on each element of the partition  ℎ , so we have to construct the local basis functions.For a noninterface element   , we can simply use the standard linear shape functions on   with degrees of freedom at average values along edges   of   and construct the linear finite element space  ℎ (  ) as follows: For this space, we define the interpolation operator  ℎ :  2 (  ) →  ℎ (  ), where the following well-known approximation property is satisfied: and we use  ℎ (Ω  ) to denote the space of piecewise  1 -nonconforming finite element space on the domain Ω  = ⋃   .For interface element   , 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   whose geometric configuration is given in Figure 2 in which the three vertices are given by  1 = (0, 0),  2 = (1, 0), and  3 = (0, 1), and the curve between points  and  is part of the interface with  = (0, ) for 0 <  ⩽ 1 and  = (, 0) for 0 <  ⩽ 1.Let  denote the line segment connecting points  and , which divides   into two parts  +  and  −  with   =  +  ⋃  −  ⋃ .Let   ,  = 1, 2, 3 be the edges of   , and let    denote the average of  along   , that is, , and then we would like to construct a new function which is linear on  +  and  −  , respectively, and satisfies the jump condition (5) on .For this purpose, we write the modified basis function φ on an interface element   as follows: with the following constraints: where  − ,  + are averages along  and  =  − / + ,   ,  = 1, 2, 3, are given values, and   is the unit normal vector on the line , that is,   = (, )/ √  2 +  2 .By similar calculation to that given in Theorem 2.2 of [8], we can rewrite (12) and (13) in the matrix form: where the coefficient matrix is defined by and the determinant of the matrix  is det Therefore, the coefficients of ( 12) are uniquely determined by conditions φ  =   ,  = 1, 2, 3, respectively, and when   ,  = 1, 2, 3, have the same value , we can get the piecewise linear function φ =  by uniqueness.Therefore we can define the following finite element space on an interface element   : = span { φ : φ|   is well defined by the above construction}. ( Also we use Ŝℎ (Ω  ) to denote the finite element space defined on the domain Ω  = ⋃   .
Remark 2. Although for functions in Ŝℎ (  ), 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 which is proved by an application of divergence theorem as in [11].
The question at hand is to discretize the concentration equation.
Multiplying the concentration equation (1)(b) by a test function V ∈  1 (Ω) and integrating over Ω which leads to the following weak formulation for concentration : Let  ℎ denote the standard piecewise linear finite element space associated with  ℎ , and ℎ  denote the bound of elements diameters.Then the discrete procedure of ( 22) is to find  :  →  ℎ such that where c0 () is the elliptic projection of  0 ().

Projection and Some Lemmas
For convergence analysis conducted in the following section, we need an interpolation operator  ℎ : H2 () → Sℎ () using the average of  on each edge by and when  ∈ H2 (Ω), we can define  ℎ  by ( ℎ )|  =  ℎ (|  ).
Then the following interpolation estimate hold for the  1 FEM approximation [8]      −  ℎ Next, let C :  →  ℎ be the elliptic projection of  defined by The function  will be chosen to assure coercivity of the form.Then, the following estimates hold [12]: ( For convergence analysis conducted in Section 5, we need to apply quote some lemmas from [9,13] and references therein. Lemma 3 (see [9] (the second Strang lemma)).Let  and  be the solution of (8) Lemma 4 (see [13]).Let  be an edge of .Then there exists a constant  > 0 such that for all , V ∈  1 (): where

Convergence Analysis
In this section, we will present convergence analysis for the pressure and concentration.
Theorem 5. Let  ∈ H2 (Ω),  ∈ Sℎ (Ω) be the solution of ( 8) and (21); then there exists a constant  > 0 independent of ℎ  and the location of the interface, such that the following result holds: 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: inf For the second part, it is a consistency error estimate.By the definition of  ℎ (⋅, ⋅) and Green's formula, we can get where  is the unit outward normal vector on each , φ ∈ Sℎ (Ω), and by the construction of the space Sℎ (Ω) we have a well-defined property on the interior edges, that is, 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 (, ) of (1) satisfies  ∈  ∞ ( 2 ) ∩  2 ( 2 ),  ∈  2 ( H2 ).Let C be the projection of  satisfying (26) and let  be the solution of (23).Suppose that the mesh parameters satisfy ℎ (36) The terms on the left side can be easily bounded by In order to prove the forth term, we need the following induction hypothesis: As the statement in [1], we assume that for some  > 0,  ⩽ ε, Combining ( 27) with ( 30) and (45), we deduce the error estimates for the pressure and the concentration.(47) Remark 8.In this paper, we just get a priori error estimates for the coefficient  = (, , ∇).However, some mathematical models for miscible displacements, which are currently being used by oil companies, make the assumption that the coefficient  = (, ) and it does not depend on ∇; with this assumption, we can prove the optimal order convergence similarly.
and (21).Then there exists a constant  > 0 is the common edge of adjacent element  1 and  2 .By boundary condition, we note that ()(/) vanishing average on the boundary.So we rewrite ∑ ∈ ℎ < ()(/), φ>  and use Lemma 4 to derive that For the fifth term, we can derive the estimate by using the boundness of C and (27), as well as (30), provided that  is Lipschitz continuous. 2 ) .