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

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 L-norm are derived for these two schemes, respectively.


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][2][3][4][5].
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; c(, ), the injected concentration at injection wells ( > 0) and the resident concentration at production wells ( < 0).Here, (u) is a tensor 2 × 2 matrix and generally has the form where 2 × 2 matrix  = (  ) satisfies   =     /|u| 2 ,  ⊥ =  − ,   is the molecular diffusivity, and   ,   are longitudinal and transverse dispersivities, respectively.Furthermore, a compatibility condition ∫ Ω (, ) = 0 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][7][8][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][15][16][17]).We have presented parallel Galerkin domain decomposition procedures for parabolic equation [18][19][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 ℎ,  2 -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  2 -norm error estimates avoid the loss of  −1/2 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  2 -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.
Similarly,   ∞ (, ; ) and the norm ‖‖   ∞ (,;) are defined.If [, ] = , we simplify our notation and write  ∞ ( 1 ∞ ) for  ∞ (0, ;  1 ∞ (Ω)), and so forth.If  = ( 1 ,  2 ) is a vector function, we note that  ∈  if  1 ∈  and  2 ∈ .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  0 ,  1 ,  0 ,  1 ,  0 ,  1 ,  1 , and  2 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.

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 (u, ) :  →  ×  such that For ℎ  > 0, we discrete (10) in space on a quasiuniform mesh T ℎ  of Ω with diameter of element ≤ ℎ  .Let  ℎ  ×  ℎ  ⊂  ×  be Raviart-Thomas spaces [6,7] of index  = 0 for this mesh.

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 / + u ⋅ ∇ be denoted by , where The weak form of (1a)-(1d) is finding a map  :  →  2 (Ω) such that Part 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  =  − (u()/())Δ  and () = (), then we get Since the problem is Ω-periodic [4],  −1 is always defined and the tangent to the characteristic (i.e., the  segment) Advances in Numerical Analysis 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  −1 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 Ω  ( = 1, 2) by an interdomain boundary Γ, which is a surface of dimension  − 1, 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 Ω 1 toward Ω 2 .
From the definition above, we note that functions  in M ℎ c have a well-defined jump [] on Γ: where (x ± ) := lim  → 0 ± (x +  Γ ).
To construct parallel algorithm, for a small constant 0 <  < min{diam(Ω 1 ), diam(Ω 2 )}, we introduce an integral mean value of a given function  ∈  2 (Ω) on the interdomain boundary Γ as Furthermore, we define the extrapolation of   (x) on Γ as Generally, near the intersection of boundary Ω and inner-boundary Γ, the value of  outside Ω is needed for   and V .If x ∉ Ω, let x ∈ Ω denote the symmetric point of x with respect to Ω.For a given function  ∈  2 (Ω), we define By ( 23), we know   and V have the values on a strip domain  = {y | y = x +  Γ ,  ∈ [−, ], ∀x on Γ}, see Figure 2. Now, we present the non-overlapping characteristic DDM for the concentration equation: where The coefficient  will be given in Lemma 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 0 =  0 <  1 < ⋅ ⋅ ⋅ <   = , with Δ   =   −  −1 .Each pressure step is also a concentration step; that is, for each  there exists  such that   =   ; in general, Δ  > Δ  .We may vary Δ  , but except for Δ 1   we drop the superscript.For functions  on Ω × , we write   () for (,   ); thus, subscripts refer to pressure steps and superscripts to concentration steps.
If  = 1, set We retain the superscript on Δ 1  because   is first order correct in time during the first pressure step and second order during later steps.
The convergence analysis will make use of an analogue of x defined for the exact velocity u  .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  −1 , so that   can be computed on Ω 1 and Ω 2 fully parallel once  −1 has been got.

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  3 is a positive constant independent of ℎ  and   is an element of the mesh T ℎ p .
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.

For the Concentration.
In this section, we adopt some following lemmas [18].
Lemma 2. For smooth enough function , there holds estimates: where   2 Γ and   4 Γ are the second and fourth normal derivative of  on Γ, respectively. where where For functions  with restrictions in  1 (Ω 1 ) ∪ 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 C ∈ M ℎ  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 C0 to be the projections of  0 () by ( 44) and take the initial conditions Let The following lemma gives the bounds of   .Lemma 6.There holds a priori estimates: Proof.Combining (39) and (44), we have ∀ ∈ M, Subtracting ( 52) from (28d) and taking  =   , we adopted the skill of [4] to obtain Since Advances in Numerical Analysis by summing (53) over  we have We estimate the terms on the right-hand side of (57) one by one.We denote the terms as  1 ,  2 ,  3 , . . .,  10 from the third term on the right-hand side of (57).We turn to analyze them one by one.
For the term  1 , similarly as that of [18], by taking we can obtain The analyses of the terms  2 , . . .,  10 are similar to that of [4] Taking  0 = 0, 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  ≥ 1,  ≥ 0 hold.Assume that the discretization parameters obey the relations Remark 10.As for the accuracy order of ℎ c and ℎ p , we know that (63) and (65) are optimal for the concentration , the velocity u and the pressure .

Extensions
To get higher-order accuracy with respect to , we apply the definition (22) to an extrapolation value Ĉ (x) 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   4 Γ is the fourth normal derivative of  on Γ.
For functions  with restrictions in  1 (Ω 1 ) ∪ Similarly as the elliptic projection (44), in order to get optimal error estimates, we introduce an elliptic projection C ∈ M ℎ c 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: provided Theorem 15.Suppose that the assumptions (), (), (), (  ), (  ) and  ≥ 1,  ≥ 0 hold.Assume that the discretization parameters obey the relations Then for ℎ c , ℎ p and Δ  sufficiently small, the errors of the approximation (66) for (1a)-(1d) satisfy Remark 17.As for the accuracy order of ℎ c and ℎ  , we know that (76) and (78) are optimal for the concentration , the velocity u 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.

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 following Figure 4 presents the approximate solution of the characteristic integral mean DDM scheme at time  = 0.5 with space step size ℎ = 1/40.
Figure 5 shows the convergence order of the characteristic integral mean DDM scheme.

Remark 19.
From Table 1, we see that the error estimate of the characteristic integral mean DDM scheme is better than that of the characteristic implicit Galerkin scheme.From Table 1 and Figure 5, we see that the convergence order is near 2 which is consisting with the analytical results.

Figure 2 :
Figure 2: The strip domain  with width 2.
Now, we turn to derive an  2 (Ω)-norm error estimate for   .