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 the error estimates are obtained for both the velocity and the pressure. Numerical examples are presented to confirm the theoretical analysis.

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 [36]), 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 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 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: where is the permeability, 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 is the Darcy velocity vector, and is the pressure. The source must satisfy the consistency constraint

Let 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 and : 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 be a closed subspace of via For the bilinear forms and , we have the standard result.

Lemma 1. The bilinear form is bounded on and coercive on , and the bilinear form is bound on . Namely,(1) there exist two constants and such that (2) there is a constant such that For the space and , the Ladyzhenskaya-Babuka-Brezzi(L-B-B) condition holds; see [15, 16], for example.

Lemma 2. There is a constant such that
It is clear that there exists a unique solution 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.

Let be a given integer and . We construct the finite-dimensional subspaces of and by introducing three different quadrangulations , , of .

First, we divide into uniform squares where and . The corresponding quadrangulation is denoted by . See Figure 1(a).

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 (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 (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 and , we construct a subspace of . Denote by , , , and the south, right, north, and left sides on the boundary of . Set where 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 , such that

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 , 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 , 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.

Lemma 3. For any , the interpolating is uniquely determined by (18).

Proof. It is easy to see that (18) is equivalent to an equation of , where is a matrix and , are vectors. Direct calculation shows that and the form of submatrix is as follows We can see that the matrix is invertible and the equation is solvable, and therefore 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 and independent of , such that
(ii) There exists a constant independent of 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, , 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 , such that

Proof. From the process above, we obtain that , any . For any , there exists , such that where is a constant independent of ; then we obtain Using Lemma 4, we have that Taking , 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 mixed element spaces [3, 5, 6] based on the partition

Define the interpolation operator by the following equations: Denote by the projection operator and by the vector projection operator. The following properties of the projections hold: Then, we have an important property about the operator .

Lemma 6. For any , there holds the equivalence ; namely,

Proof. As the definition of is based on each element , we focus our discussion on arbitrary element , . Firstly, we consider the -component (see Figure 4). The analysis for -component is similar.
For a function , on an element , it is uniquely given by its node values . As is a continuous bilinear function on each of the two parts as shown in Figure 4. Then, from (32), we know that is given by We deduce that It is clear that we just need to verify (34) for both and .
We first consider . 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, so
By direct computation, we can easily see that has the same value, so When , we have that where are defined in (36). Next, we compare the coefficients of in (40) with the coefficients in , which determine as the coefficient of . With similar computation, we obtain that Comparing with (40), we can find that (34) is true with .
So, we certify the lemma.

Theorem 7. If satisfy (6) and satisfy (17), then there exists a positive constant independent of such that the following error estimates hold:

Proof. First, we focus on the error . From (6), (17), (18), and (32), we derive that Let in (6); then Namely, Here, we used the property . Subtracting from (17), we get that Take Then
Due to (44), we find that
Now, we analyze the error based on the equations above where , are positive constants. Take the value of , 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: . So
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 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 on the normal space, not in the subspace of , and the same issue also occurs in [11]. The problem is that testing 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 . Failing to obtain convergence rate of 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 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 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 superconvergence in the 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 , the degrees of freedom of Arbogast’s scheme are and for our formulation. The convergence rate of 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: where and are the gravity coefficient and vertical coordinate, is the porosity of the rock, and represents a known source. is the molecular diffusion and mechanical dispersion coefficient. For convenience, we denote that and . Let be the solution of the ordinary differential equation Let , , ; then, we derive the entire weak formulation for the model: find , such that Let be the time step for both concentration and pressure; define 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 We can see that the continuous flux is indispensable for . Let in (64), and summing it up from 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: Now, we select , and the following analytical solution of the problem is 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.

Acknowledgment

The work is supported by the National Natural Science Foundation of China Grant no. 11171190.