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

and Applied Analysis 3


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][4][5][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 2 error estimates of velocity and pressure are also established.
The outline of the rest of this paper is organized as follows. In Sections 2 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 2 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.

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: where > 0 is the permeability, > 0 is the viscosity, and is the volumetric flow rate source or sink. Γ is the boundary of Ω, and is the unit outward normal vector to Γ. The variable = ( 1 , 2 ) is the Darcy velocity vector, and p is the pressure. The source must satisfy the consistency constraint Let 2 (Ω) be the space of square integrable function in Ω with inner product (⋅, ⋅) and norm ‖ ⋅ ‖. We use the notation of the Hilbert space with norm Define the following subspaces of (div , Ω) and 2 (Ω): The classical weak variational formulation of Problem (1) is as follows: find ( , ) ∈ × , such that Here, The following discussion and discrete analysis are related to the weak form (6). Let 0 be a closed subspace of via 0 = {V ∈ : (V, ) = 0, ∀ ∈ } . (1) there exist two constants 1 > 0 and > 0 such that (2) there is a constant 2 > 0 such that For the space V 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 It is clear that there exists a unique solution ( , ) ∈ × to the Problem (6).

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.
Let be a given integer and ℎ = 1/ . We construct the finite-dimensional subspaces of and by introducing three different quadrangulations ℎ , 1 ℎ , 2 ℎ of Ω.
where = ℎ and = ℎ. The corresponding quadrangulation is denoted by ℎ . See Figure 1 Then, for all , ∈ ℎ , we connect all the neighbor midpoints of the vertical sides of , 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 1 ℎ (see Figure 1(b)). Similarly, for all , ∈ ℎ , we connect all the neighbor midpoints of the horizontal sides of , by straight line segments if the neighbor midpoints have the same horizontal coordinate. Then, we obtained the third quadrangulation of Ω, which is denoted by 2 ℎ (see Figure 1(c)).
Based on the quadrangulation ℎ , we define the piecewise constant functional space used to approximate the pressure ℎ is a subspace of . Furthermore, using the quadrangulations 1 ℎ and 2 ℎ , we construct a subspace of . Denote by Γ 1 , Γ 2 , Γ 3 , and Γ 4 the south, right, north, and left sides on the boundary of Ω. Set where 1,1 denotes the piecewise bilinear polynomial space with respect to the variables and . Let Obviously, ℎ is a subspace of . Using the subspaces ℎ and ℎ instead of and in the variational Problem (6), we obtain the discrete problem: find (17)

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 ℎ , we divided the edges of all squares into two sets. The first one denoted by contains all vertical edges, and the second one denoted by contains all horizontal edges. We define the interpolation operator Π : → ℎ by Π = (Π 1 ℎ 1 , Π 2 ℎ 2 ) ∈ 1 ℎ × 2 ℎ , which satisfy the following: where is a set of edges of elements got by bisecting the most bottom element edges and the most top element edges of and are got by bisecting the most left element edges and the most right element edges of . See Figures 2 and 3.
Proof. It is easy to see that (18) is equivalent to an equation of = , where A is a matrix and X, B are vectors. Direct calculation shows that and the form of submatrix 1 is as follows ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) . (20) We can see that the matrix is invertible and the equation is solvable, and therefore X can be uniquely determined.
Assume that the solution ( , ) of Problem (6) has the following smoothness properties: Then, we should give the following lemma about the properties of the interpolations defined in (18).

Lemma 4. (i) There exist two constants 3 and 4 independent of h, such that
(ii) There exists a constant 5 independent of h such that (iii) For any ∈ , we have that 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 Here, ⃗ = ( 1 , 2 ), 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 Proof. From the process above, we obtain that (V, ℎ ) = (ΠV, ℎ ), any V ∈ , ℎ ∈ ℎ . For any ℎ ∈ ℎ , there exists where 6 is a constant independent of ℎ ; then we obtain Using Lemma 4, we have that Taking = 1/ 5 6 , 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 0 mixed element spaces ℎ × ℎ [3, 5, 6] based on the partition ℎ Define the interpolation operator Π : → ℎ by the following equations: Denote by : → ℎ the 2 projection operator and by : → ℎ the ( 2 (Ω)) 2 vector projection operator. The following properties of the projections hold: Then, we have an important property about the operator Π .
For a function 1 ∈ 1 ℎ , on an element , it is uniquely given by its node values , = 1, . . . , 6. As 1 is a continuous bilinear function on each of the two parts as shown in Figure 4. Then, from (32), we know that Π 1 = + is given by We deduce that It is clear that we just need to verify (34) for both V = 1 and V = .
We first consider V = 1. Denote by the node basis function at the point , which implies that ( ) = , , which has the value 1 if and only if = ; otherwise, it is zero. By direct calculation, we can get the basis, for example, By direct computation, we can easily see that ∫ Π 1 has the same value, so When V = , we have that where , are defined in (36). Next, we compare the coefficients of in (40) with the coefficients in ∫ 1 * , which determine 1 as the coefficient of 1 . With similar computation, we obtain that 5 = 1 , 2 = 6 = 0 ℎ 2 8 + ℎ 3 12 , Comparing with (40), we can find that (34) is true with V = . So, we certify the lemma. (6) and ( ℎ , ℎ ) satisfy (17), then there exists a positive constant independent of ℎ such that the following error estimates hold:

Theorem 7. If ( , ) satisfy
Proof. First, we focus on the error − ℎ . From (6), (17), (18), and (32), we derive that Let V = Π V ℎ in (6); then Namely, Here, we used the property ∇⋅Π V = ∇⋅V. Subtracting from (17), we get that Take Then Due to (44), we find that Abstract and Applied Analysis 7 Now, we analyze the error − ℎ based on the equations above where > 0, = 1, 2, 3 are positive constants. Take the value of 1 = 3 = /4 , 2 = 1, and combining with (22) and (33), we conclude that We also can obtain a higher order error estimate for ‖ − ℎ ‖. Consider the classical duality argument. Let be the solution of the following elliptical problem: By the elliptic regularity, the estimate holds: Now, we estimate the right hand terms of the above inequality. From (33), (22), and (52), we have It is easy to see that Here, we used the fact that (Π − ℎ , V ∇ ) = 0 which is got from the Green formulation and (44).
Combining the above inequalities, we conclude that 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 for ( ℎ , V ℎ ) on the normal 2 space, not in the subspace of V 0ℎ = {V ℎ ∈ ℎ : (V ℎ , ℎ ) = 0, ∀ ℎ ∈ ℎ }, and the same issue also occurs in [11]. The problem is that testing (∇ ⋅ V, ) by ∈ ℎ does not control the full divergence of , 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 ‖ − ℎ ‖ 0 . Failing to obtain convergence rate of ‖ − ‖ (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.

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. We compare our method to the formulation constructed by Arbogast and Wheeler [11]. Its corresponding discrete finite element spaces are 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. Close results of numerical errors for both formulations are shown in Table 2. From Table 3, we can see that converges to ℎ as (ℎ) and − ℎ as (ℎ 2 ) for our formulation, which both agree with the theorem. From the examples, we can observe that ℎ converges to somewhat better than expected, and it appears that on the uniform grid we attain (ℎ 3/2 ) superconvergence in the 2 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 ‖ − ℎ ‖ (div ,Ω) is first order, but here we cannot give the corresponding analysis.

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: Let Δ be the time step for both concentration and pressure; define (62) 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 We obtain the corresponding full-discrete mass-conservative characteristic (MCC) scheme: find ( ℎ , ℎ , ℎ ) ∈ ℎ × ℎ × ℎ , such that where = det ( ) Abstract and Applied Analysis 9   We can see that the continuous flux is indispensable for . Let ℎ = 1 in (64), and summing it up from = 1 to , we get the mass balance 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: compute mass error : ∫ 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. As can be seen from Tables 4 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.