Advances in Numerical Analysis

Volume 2013 (2013), Article ID 303952, 13 pages

http://dx.doi.org/10.1155/2013/303952

## Parallel Nonoverlapping DDM Combined with the Characteristic Method for Incompressible Miscible Displacements in Porous Media

School of Mathematics, Shandong University, Jinan, Shandong 250100, China

Received 28 May 2012; Revised 27 October 2012; Accepted 1 November 2012

Academic Editor: J. Rappaz

Copyright © 2013 Keying Ma and Tongjun Sun. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Two types of approximation schemes are established for incompressible miscible displacements in porous media. First, standard mixed finite element method is used to approximate the velocity and pressure. And then parallel non-overlapping domain decomposition methods combined with the characteristics method are presented for the concentration. These methods use the characteristic method to handle the material derivative term of the concentration equation in the subdomains and explicit flux calculations on the interdomain boundaries by integral mean method or extrapolation method to predict the inner-boundary conditions. Thus, the velocity and pressure can be approximated simultaneously, and the parallelism can be achieved for the concentration equation. The explicit nature of the flux prediction induces a time step limitation that is necessary to preserve stability. These schemes hold the advantages of nonoverlapping domain decomposition methods and the characteristic method. Optimal error estimates in -norm are derived for these two schemes, respectively.

#### 1. Introduction

The two-phase fluid displacements in porous media is one of the most important basic problems in the oil reservoir numerical simulation. It is governed by a nonlinear coupled system of partial differential equations with initial and boundary values. In this paper, we will consider the following incompressible miscible case: the pressure is governed by an elliptic equation and the concentration is governed by a convection-diffusion equation [1–5].

where is a bounded domain in , , and is nonzero at injection wells only. The variables in (1a)–(1d) are the pressure in the fluid mixture, the Darcy velocity , and the relative concentration of the injected fluid. The is the unit outward normal vector on boundary .

The coefficients and data in (1a)–(1d) are : the permeability of the porous media; , the viscosity of the fluid mixture; : representing flow rates at wells; and , the gravity coefficient and vertical coordinate; : the porosity of the rock; , the injected concentration at injection wells () and the resident concentration at production wells (). Here, is a tensor matrix and generally has the form where matrix satisfies , , is the molecular diffusivity, and , are longitudinal and transverse dispersivities, respectively. Furthermore, a compatibility condition must be imposed to determine the pressure.

The pressure equation is elliptic and easily handled by standard mixed finite element method, which has been proven to be an effective numerical method for solving fluid problems. It has an advantage to approximate the unknown variable and its diffusive flux simultaneously. There are many research articles on this method [6–9]. The concentration equation is parabolic and normally convection dominated. It is well known that standard Galerkin scheme applied to the convection dominated problems does not work well and produces excessive numerical diffusion or nonphysical oscillation. A variety of numerical techniques have been introduced to obtain better approximations for (1a)–(1d), such as characteristic finite difference method [10], characteristic finite element method [11], the modified of characteristic finite element method (MMOC-Galerkin) [12], and the Eulerian-Lagrangian localized adjoint method (ELLAM) [13].

It is well known that parallel algorithms, based upon overlapping or non-overlapping domain decompositions, are effective ways to solve the large scale of PDE systems for most practical problems in engineering (e.g., see [14–17]). We have presented parallel Galerkin domain decomposition procedures for parabolic equation [18–20]. These procedures use implicit Galerkin methods in the subdomains and explicit flux calculations on the interdomain boundaries by integral mean method or extrapolation-integral mean method. Some constraints for time step are still needed for these procedures to preserve stability, but less severe than that for fully explicit methods. With respect to the accuracy order of , -norm error estimates are optimal for higher-order finite element spaces and almost optimal for linear finite element space in twodimensional domain. Compared with Dawson-Dupont's schemes [21], these -norm error estimates avoid the loss of factor.

We also have considered using the procedures in [18] for wave equation [22] and convection-diffusion equation [23]. This paper is one of our sequent research papers. The main purpose is to use parallel Galerkin domain decomposition procedures in [18] combined with the characteristic method for the concentration equation of incompressible miscible displacements in porous media. This paper is organized as follows. In Section 2, we first present mixed finite element method for the velocity and pressure and then formulate parallel non-overlapping domain decomposition methods combined with the characteristics method for the concentration. We establish the combined approximation schemes. Auxiliary lemmas are listed in Section 3, which show some properties of finite element spaces and projections. In Section 4, we derive the optimal-order -norm error estimates. Finally in Section 5, we extend the consideration for another approximate scheme by using extrapolation method. It is worthwhile to specially emphasize that the research of this paper is creative. No former researchers discussed it. These schemes not only hold the advantages of non-overlapping domain decomposition methods, but also hold the advantages of the characteristic method.

#### 2. Formulation of the Methods

In this paper, we adopt notations and norms of usual Sobolev spaces [3]: In particular, and . The inner product on is denoted by .

We also use the following spaces that incorporate time dependence. Let and let be any of the spaces just defined. For suitably smooth on , we let Similarly, and the norm are defined. If , we simplify our notation and write for , and so forth.

If is a vector function, we note that if and . We also use the vector-function spaces and norms:

We need some assumptions. The regularity assumptions on the solution of (1a)–(1d) are noted by [4]:

We also require the following assumptions on the coefficients in (1a)–(1d) [4]. Let , , , , , , , and be positive constants such that Further assumptions will be made in individual theorems as necessary.

For convenience, we assume that (1a)–(1d) is periodic (see [4]); that is, all functions will be assumed to be spatially -periodic throughout the rest of this paper. This is physically reasonable, because no-flow conditions (1c) are generally treated by reflection, and in general, interior flow patterns are much more important than boundary effects in reservoir simulation. Thus, the boundary conditions (1c) can be dropped.

Throughout the analysis, and will denote generic positive constants, independent of , , , and , but possibly depending on constants in (), norms in (). Similarly, will denote a generic small positive constant.

##### 2.1. A Mixed Finite Element Method for the Pressure and Velocity

Let

As in [4], the pressure equation (1a) is equal to the following saddle-point problem of finding a map such that where the bilinear forms

For , we discrete (10) in space on a quasi-uniform mesh of with diameter of element . Let be Raviart-Thomas spaces [6, 7] of index for this mesh.

The mixed method for pressure and velocity, given a concentration approximation at a time , consists of and such that Existence and uniqueness of and are proved in [1], based on ideas of [24].

##### 2.2. A Characteristic DDM for the Concentration

In this section, we assume in (1a)–(1d) is given. Define

Let and let the characteristic direction associated with the operator be denoted by , where

The weak form of (1a)–(1d) is finding a map such that

Part into , with . The analysis is valid for variable time steps, but we drop the superscript from for convenience. For functions on , we write for . Approximate by a backward difference quotient in the -direction,

If we let and , then we get Since the problem is -periodic [4], is always defined and the tangent to the characteristic (i.e., the segment) cannot cross a boundary to an undefined location. The difference quotient relates the concentration at a given at time to the concentration that would flow to from time if the problem were purely hyperbolic.

The time difference (18) will be combined with a characteristic DDM in the space variables. We recall the domain decomposition procedures in [18] here. For simplicity and without losing generality, we only discuss the case of two subdomains. But the algorithms and theories can be extended to the case of many subdomains. Divide into two subdomains () by an interdomain boundary , which is a surface of dimension , see Figure 1. We denote by the part of the boundary of the subdomains which coincides with . Denote the unit vector normal to as , which points from toward .

Let be quasi-uniform partitions of () and . Here, denotes the maximal element diameter of . We construct the finite element space on which satisfies the following condition (I) (1)For , let be a finite element subspace of , and let such that if , then . (2)For , , where is a polynomial space of degree at most . (3)For , , the integer , and , there exists a positive constant independent of such that where .

From the definition above, we note that functions in have a well-defined jump on : where .

To construct parallel algorithm, for a small constant , we introduce an integral mean value of a given function on the interdomain boundary as Furthermore, we define the extrapolation of on as

Generally, near the intersection of boundary and inner-boundary , the value of outside is needed for and . If , let denote the symmetric point of with respect to . For a given function , we define By (23), we know and have the values on a strip domain , see Figure 2.

Now, we present the non-overlapping characteristic DDM for the concentration equation: where The coefficient will be given in Lemma 3.

##### 2.3. The Combined Approximate Schemes

We now present our sequential time-stepping procedure that combines (24) and (12). As in [4], let us part into pressure time steps , with . Each pressure step is also a concentration step; that is, for each there exists such that ; in general, . We may vary , but except for we drop the superscript. For functions on , we write for ; thus, subscripts refer to pressure steps and superscripts to concentration steps.

If concentration step relates to pressure steps by , we require a velocity approximation for (25) based on and earlier values. If , take the linear extrapolation of and defined by If , set We retain the superscript on because is first order correct in time during the first pressure step and second order during later steps.

The combined time-stepping procedure is finding a map and a map such that where is given by (45).

The steps of the above calculation are as follows: Step 1. known solve by (28b) and (28c); Step 2. On each domain, use (28d) to parallelly compute until ; Step 3. Then by (28b) and (28c) for ; Step 4. Calculate the approximations in turn analogically to get the pressure, velocity, and concentration at other time-step, respectively.

The convergence analysis will make use of an analogue of defined for the exact velocity . If is a function on , set The time step will be clear from the context.

*Remark 1. * In the scheme (28a)–(28d), the numerical flux on is computed explicitly from , so that can be computed on and fully parallel once has been got.

#### 3. Auxiliary Lemmas

We adopt some auxiliary lemmas about the finite element spaces, which will be used in the next section.

##### 3.1. For the Pressure and Velocity

The Raviart-Thomas spaces in Section 2.1 possess the following approximation and inverse properties [7]: where is a positive constant independent of and is an element of the mesh .

Define the map by (see [7]) where is the exact solution of (1a)–(1d).

By arguments in [1, 24], the map exists and implies that The positive constant depends on constant in () but independent of , , , and .

The estimate (33) and imply that As in [1], comparison of (32) and (12) implies that The estimates (33) and (35) will handle the coupling of concentration and velocity errors in the convergence analysis.

##### 3.2. For the Concentration

In this section, we adopt some following lemmas [18].

Lemma 2. * For smooth enough function , there holds estimates:
**
where and are the second and fourth normal derivative of on , respectively. *

Lemma 3. * Let . If and is small, one has
**
where
*

Define an elliptic projection for : , , where . Let . From [25–28], we see the following.

Lemma 4. *For defined by (39), there holds the following.**-norm error estimate
**-norm error estimate
*

For functions with restrictions in , we define a norm

Lemma 5. * There exists a positive constant such that for small ,
*

As we have shown, the combined approximation scheme (28d) includes two terms on the interdomain boundary by integral mean method to present explicit flux calculation. These terms are distinct ones different from Dawson-Dupont's schemes such that the standard elliptic projection (39) is insufficient for optimal error estimates. To get optimal error estimates, we need a new elliptic projection including terms on . This new elliptic projection of the solution is defined as It follows from Lemma 5 that the projection problem (44) has a unique solution for small . Choose the initial to be the projections of by (44) and take the initial conditions Let The following lemma gives the bounds of .

Lemma 6. *There holds a priori estimates:
*

#### 4. A Priori Error Estimates in -Norm

Now, we turn to derive an optimal priori error estimate in -norm for the concentration of approximation (28d). Optimal error estimates for velocity in and pressure in follow at once from (33) and (35).

It follows from trace theorem in [29] that which will be used in the following proof. Then, we can derive an -norm error estimate for .

Lemma 7. * For defined by (46), there holds the following error estimate:
**
provided
*

*Proof. *Combining (39) and (44), we have ,
Subtracting (52) from (28d) and taking , we adopted the skill of [4] to obtain
Since
by summing (53) over we have
Furthermore, noting that
where , and taking , we have
We estimate the terms on the right-hand side of (57) one by one. We denote the terms as from the third term on the right-hand side of (57). We turn to analyze them one by one.

For the term , similarly as that of [18], by taking
we can obtain
The analyses of the terms are similar to that of [4]
Combining the above analyses, we have
Taking , applying the Gronwall Lemma and Lemmas 5 and 6, we can derive (50).

Applying the triangle inequality, Lemmas 4 and 7, we have the following error theore.

Theorem 8. * Suppose that the assumptions , , , , and hold. Assume that the discretization parameters obey the relations
**
Then for , , and sufficiently small, the errors of the approximation (28a)–(28d) for (1a)–(1d) satisfy
**
provided
**
where is given by (51). Here is a positive constant dependent on , , , , , , , and , but independent of , , , and . *

By combining Theorem 8 with (33) and (35), we obtain at once the following result.

Corollary 9. * Under the assumption of Theorem 8, the errors for velocity , and pressure are obtained by
*

*Remark 10. * As for the accuracy order of and , we know that (63) and (65) are optimal for the concentration , the velocity and the pressure .

#### 5. Extensions

To get higher-order accuracy with respect to , we apply the definition (22) to an extrapolation value on the interdomain boundary . Similar to the analysis of [18], we can present another approximate scheme for (1a)–(1d):

Since the differences between two schemes (28a)–(28d) and (66) are the second and third term to calculate the flux on the inner-domain boundary , the convergence analysis of the scheme (66) is similar to that of the scheme (28a)–(28d). For the sake of brevity, we describe the processes of analysis for (66) simply.

The proof of Theorem 15 is based on the following four lemmas. The former three lemmas are adopted from [18]. We omit the proof of Lemma 14, which is similar to that of Lemma 7 in Section 4.

Lemma 11. *For sufficiently smooth function , there holds estimates:
**
where is the fourth normal derivative of on . *

For functions with restrictions in , we use the definition of the norm

Lemma 12. * There exists a positive constant such that for small ,
*

Similarly as the elliptic projection (44), in order to get optimal error estimates, we introduce an elliptic projection of the solution as follows: It follows from Lemma 12 that the projection problem (70) has a unique solution for small .

Let The following lemma gives the bounds of .

Lemma 13. * There holds a priori estimates:
*

Now, we turn to derive an -norm error estimate for .

Lemma 14. * For defined by (71), there exists the following error estimate
**
provided
*

Theorem 15. * Suppose that the assumptions , , , , and , hold. Assume that the discretization parameters obey the relations
**
Then for , and sufficiently small, the errors of the approximation (66) for (1a)–(1d) satisfy
**
provided
**
where is given by (74). Here is a positive constant dependent on , , , , , , , , but independent of , , and . *

By combining Theorem 15 with (33) and (35), we obtain at once the following result.

Corollary 16. *Under the assumption of Theorem 15, the errors for velocity and pressure are obtained by
*

*Remark 17. * As for the accuracy order of and , we know that (76) and (78) are optimal for the concentration , the velocity and the pressure .

*Remark 18. * From Theorem 15, we can know that the scheme (66) has the same accuracy order as that of the scheme (28a)–(28d) with respect to , and , except that has an accuracy of higher order for . This shows that the scheme (66) can use larger width of middle strip domain than that of the scheme (28a)–(28d) so that the time step constraint is more weaker.

#### 6. Numerical Experiments

In this section, we present some numerical experiments for the procedures described above. All computer programs below are written by Fortran 90 code and run on a Lenovo PC with Intel(R) Core i5 3.2 GHz CPU and 4 GB memory. The resulting linear systems of algebraic equations are solved by banded Gaussian elimination. Single precision is used for all calculations.

The main purpose of this paper is to analyze the integral mean non-overlapping DDM combined with the characteristic method for the concentration equation. There were many experimental results for the integral mean non-overlapping DDM in [18–20, 22, 23]. For simplicity, we consider the following convection dominated diffusion equation where, , , , , , . We choose suitably so that the exact solution of (79) is , see Figure 3.

We consider two scenarios: (1) the characteristic implicit Galerkin scheme (*Charac-teristic-FEM scheme*) on uniform mesh; that is, no domain decomposition; (2) the characteristic integral mean DDM scheme (*Characteristic-DDM scheme*) on global uniform mesh with two equal sub-domains , , with the interdomain boundary . In these runs, we approximate (79) by using linear finite element with 4-node quadrilateral mesh on ,