Domain Decomposition Modified with Characteristic Mixed Finite Element and Numerical Analysis for Three-Dimensional Slightly Compressible Oil-Water Seepage Displacement

A parallel algorithm is presented to solve three-dimensional slightly compressible seepage displacement where domain decomposition and characteristics-mixed finite element are combined. Decomposing the computational domain into several subdomains, we define a special function to approximate the derivative at interior boundary explicitly and obtain numerical solutions of the saturation implicitly on subdomains in parallel. The method of characteristics can confirm strong stability at the fronts, and can avoid numerical dispersion and nonphysical oscillation. It can adopt large-time step but can obtain small time truncation error. So a characteristic domain decomposition finite element scheme is put forward to compute the saturation. The flow equation is computed by the method of mixed finite element and numerical accuracy of Darcy velocity is improved one order. For a model problem we apply some techniques such as variation form, domain decomposition, the method of characteristics, the principle of energy, negative norm estimates, induction hypothesis, and the theory of priori estimates of differential equations to derive optimal error estimate in l2 norm. Numerical example is given to testify theoretical analysis and numerical data show that this method is effective in solving actual applications. Then it can solve the well-known problem.


Mathematical Model and Physical Background
High-pressure pump injects water into oil storage and displaces crude oil out from production wells.This technique is popular and important in modern oil exploration.The displacement of two phase describes the physical phenomena how injected water displaces crude oil in reservoir and crude oil is produced.In modern oil recovery we try to make remaining crude oil out by adopting a third-recovery technique (chemical displacement).It is necessary to consider the compressibility in numerical simulation to avoid numerical distortion.Douglas and other scholars put forward a miscible mathematical model with slight-compressibility and discuss the methods of characteristic finite element and characteristic mixed element, then they give the outline of modern numerical simulation in oil recovery (Douglas, & Roberts, 1983;Ewing, 1983;Yuan, 1992Yuan, ,1993Yuan, ,2013)).
The mathematical model is stated by the following nonlinear system of partial differential equations with initial-boundary value conditions (Douglas, & Roberts, 1983;Ewing, 1983;Yuan, 1992Yuan, ,1993Yuan, ,2013)): where Ω is a bounded domain in R 3 .d(c) = ϕ(X) The given function c denotes the saturation of injection well.k(X) denotes the permeability and µ(c) is the viscosity, then a(c) = k(X)µ −1 (c).The porosity, the pressure and the production rate are denoted by ϕ = ϕ(X), p(X, t) and q(X, t), respectively.q = max{q, 0}.u = u(X, t) is Darcy velocity and D = ϕ(X)d m I, where d m is diffusion coefficient and I is a unit matrix.
Suppose that the fluid is not permeable at the boundary, that is to say that the following conditions hold where ν is the normal vector of ∂Ω, the boundary surface of Ω.
Combining normal finite difference or finite element with the method of characteristics, they present two different composite schemes.These schemes can reflect the first-order hyperbolic nature of convection-diffusion equation, decrease truncation error, overcome numerical oscillation and dispersion and they can improve the stability and accuracy greatly.The compressibility must be considered in new numerical simulation of modern enhanced oil displacement (Douglas, & Roberts, 1983;Yuan, 1992Yuan, ,2013)).Under periodic assumption Douglas and Yuan firstly discuss characteristic finite element and characteristic mixed finite element, obtain optimal order error estimate in L 2 norm and give a powerful tool to solve the well-known problem (Douglas, & Roberts, 1983;Yuan, 1999Yuan, ,2003)).
In numerical simulation of modern oil field exploration and development (special for enhanced chemical recovery), the computation is large-scaled, and it runs not only on a three-dimensional domain but also on a long time interval.Its nodes amount up to tens of thousands or millions, therefore traditional methods can not solve this problem well.So new modern parallel computation methods should be introduced (Ewing, 1983;Shen, Liu, & Tang, 2002).Dawson, Dupont and Du firstly present Galerkin domain decomposition procedure and give convergence analysis (Dawson, & Du, 1990;Dawson, Du, & Dupont, 1991;Daswon, & Dupont, 1992,1994) for a simple parabolic equation.For heat conductor transient problem we have published many research results (Yuan, Chang, Li, & Sun, 2015) about domain decomposition modified by characteristic finite element and characteristic mixed finite element.We have considered the enhanced oil recovery simulation with incompressible condition and give the primary study (Yuan, Chang, Li, & Sun, 2015).Since computational task of the saturation is dominated and large-scaled in the whole numerical simulation, so parallel computation of the saturation is considered mainly in this paper.Based on the above research we put forward a modified characteristic domain decomposition method to solve three-dimensional compressible seepage displacement in this paper.Decomposing computational domain into several subdomains, we define a special function to approximate the value at interior boundary explicitly and obtain numerical solution implicitly in parallel in subdomains.The flow equation is discretized by the method of mixed finite element and the saturation is approximated by a domain decomposition scheme of modified characteristic finite element.For the model problem we use variation form, domain decomposition, the characteristics, the principle of energy, induction hypothesis, the theory and technique of priori estimates of partial differential equations to get optimal order error estimate in L 2 norm.Numerical experiment is consistent with theoretical analysis and confirms that the method is efficient and feasible in actual computation.It is an important and powerful tool in model analysis, numerical method, principle research and engineering applicable software design of modern oil reservoir exploration and development and it can solve the well-known problem (Douglas, & Roberts, 1983;Ewing, 1983;Shen, Liu, & Tang, 2002).
In the following discussion the symbols K and ε denote a generic positive constant and a generic small positive number, respectively.They have different definitions at different places.

Some Preliminary Notations
For simplicity, decompose To approximate the normal derivative at Γ, we define two special functions Φ 2 and Φ 4 as follows (Dawson, & Du, 1990;Dawson, & Dupont, 1992), (5a) Note that if p(x 1 ) is a polynomial of degree at most one, then and if p(x 1 ) is a polynomial of degree at most three, then Let N h, j denote a finite-dimensional finite element space of H 1 (Ω j )( j = 1, 2), and let N h (Ω) be an l-dimensional subspace of L 2 (Ω).Moreover, for a function v ∈ N h , then v| Definition 2: A bilinear function D(u, v) is defined by where u, v ∈ H 1 (Ω j ), j = 1, 2. D(X) is a positive definite function, D T = 1 and λ s is a positive constant.Definition 3: An integral operator is defined to approximate normal derivative at interior boundary, where Φ(x 1 ) is defined by (7).
Let (•, •) denote inner product in L 2 (Ω 1 ∪ Ω 2 ), and omit the subscript (ψ, ρ) = (ψ, ρ) Ω on Ω = Ω 1 ∪ Ω 2 .For a function ψ in H 1 (Ω 1 ) and H 1 (Ω 2 ), define Noting that then, we have Rewrite the second term of the above expression as follows, ) is the value of Taylor remainder and ξ 1 (X) is a point between 1 2 and x 1 .Then where Therefore, there exists a positive constant M 0 such that Similarly, we have the following estimates for 0 where M 1 , M 2 , M 3 are positive constants, m = 2, 4, and ∂u ∂γ denotes the normal derivative of u across interior boundary Γ.

Modified Characteristic Mixed Finite Element Domain Decomposition Procedure
The variation of ( 2) is defined by where The pressure is solved by the method of mixed finite element.For a vector f , periodic} and L2 (Ω) = {g : g ∈ L 2 (Ω), periodic}.For convenience we omit the symbol " ∼ " and let V = H(div; Ω), W = L 2 (Ω).Considering the equations on Ω 1 , Ω 2 , we introduce the following compatibility condition where then we get the following variation of (1) on subdomains ) Making summation of (17a) and (17b) on i = 1, 2, we get the variation of (1) on the whole domain Ω ( where The elliptic projections of the saturation, Darcy velocity and the pressure are defined as follows. Definition 4: The elliptic projection of c(X, t) is defined by c(X, t) : where λ is a positive constant.
Definition 5: The projections of u(X, t) and p(X, t) are defined by { ũ, p} : where W h × V h is Raviart-Thomas space with the index k and the mesh step h p .The approximations are stated as follows (Ewing, Russell, & Wheeler, 1984;Raviart, & Thomas, 1977) inf Lemma 1.By using Galerkin method we get error estimate of elliptic projection for the saturation (Ewing, & Wheeler, 1981;Wheeler, 1973) in the finite element space N h with the index l and the mesh step h c , Lemma 2. By Brezzi theory (Brezzi, 1974;Ciarlet, 1978) we get error estimates of mixed finite element elliptic projections for Darcy and the pressure Considering that the fluid flows along the characteristics ϕρ ∂c ∂t + u • ∇c, so we introduce the method of characteristics.Define The characteristics depends on c, p and Darcy velocity u.Then ( 2) is expressed by Let N h ⊂ N denote an l-dimensional subspace with mesh step h c , W h × V h denote a k-order Raviart-Thomas mixed finite element space with h p , and let Λ h = {β : β| Γ ∈ P k (Γ)} denote a subspace of Λ.
The parallel procedures of characteristic mixed element are stated as follows.Given numerical solutions where where Ĉn  (10).By the definition of N h = N h,1 ∪ N h,2 , the problem is decomposed into two independent subproblems on subdomains, then the saturation is solved in parallel.The solutions exist and are sole according to (C).

Numerical Analysis of a Model Problem
In this section, we show how to finish the convergence analysis and error estimate of the composite method.For convenience, we only discuss a model problem by simplifying the coefficients and descriptions.Suppose that the problem is incompressible, a(c) = k(X)µ −1 (c) ≈ k(X)µ −1 0 = a(X), i.e. the viscosity µ(c) is approximately taken as a constant such as the displacement of lower seepage (Douglas, & Roberts, 1983;Yuan, 2001).Then (1) and ( 2) are degraded into The variation of ( 27) at saddle point is formulated as follows ( The parallel algorithm of characteristic mixed element is defined as follows.
If the following constraint condition holds for sufficiently small ∆t An induction hypothesis is introduced By using (48) we get estimates of the other terms on the right-hand side of (43) as follows Collecting all the estimates (44)-( 49) for ( 43), we have Applying ( 22) in (50), multiplying the resulting expression by 2∆t, summing them on n (0 ≤ n ≤ L − 1), and noting that ξ 0 = 0, we have Then it follows from ( 40) and (51), Using Gronwall Lemma, we have It remains to testify the hypothesis (48).As n = 0, it is true because of ξ 0 = 0 and σ 0 = 0.If (48) holds for 0 ≤ n ≤ L − 1.
For n = L and k, l ≥ 1, we get the conclusion from (53) and the following constraint From ( 53), ( 22) and ( 33), we obtain the following convergence theorem.
Numerical approximations of normal derivative at interior boundary ∂c ∂x (0.5) = e T sin(π) = 0 is shown in Table 3.The approximation at interior boundary is quite well so that the parallel computations on subdomains have high order accuracy.4 (unit: second).We can see that the computational scale becomes large and domain decomposition is more efficient as h → 0. One reason is that domain decomposition can solve the problem in parallel and saves computational time.The other is the computational scale of each subdomain becomes a half of the computation of whole domain.It is hoped that domain decomposition is used in actual applications efficiently and can solve the complicated problem quickly as the partition becomes refined and the number of subdomains becomes large.In numerical simulation of enhanced oil recovery (chemical recovery), we should consider multicomponent seepage displacement and formulate its mathematical model by the following nonlinear system of partial differential equations with initial-boundary value conditions (Douglas, & Roberts, 1983;Yuan, 2013) The pressure and the saturation of α-th component are denoted by p(X, t) and c α (X, t), c α (X, t) = 1, so only n c − 1 component saturations are independent.In this section we aim to find numerical solutions of the pressure p(X, t) and the saturations c(X The parallel algorithm of characteristic mixed element is defined as follows.If numerical solutions at After a similar analysis, we derive the following convergence statement.

Conclusions and Discussions
A domain decomposition method of modified characteristic mixed finite element is discussed for slightly compressible oilwater seepage displacement in this paper.In §1 Induction, the mathematical model, physical background and some related research studies are stated.In §2 some preliminary work and elementary error estimates are given.In §3 we formulate the parallel procedures of characteristic mixed finite element.Numerical analysis and optimal order l 2 norm estimate are shown in §4.Finally a numerical example is given to show the efficiency and application in §5.In §6 we construct the domain decomposition algorithm of characteristic mixed element for multi-component seepage displacement and give numerical analysis.In this paper some features are concluded as follows.(I) This method is suitable to solve oil reservoir numerical simulation on three-dimensional irregular geometric region especial enhanced oil numerical computation.(II) By using mixed finite element, the method improves an order accuracy to compute Darcy velocity.It is most important in numerical simulation of oil reservoir.(III) The method of characteristics can confirm strong stability of numerical simulation at the fronts, and can avoid numerical dispersion and nonphysical oscillation.It can adopt large-time step but can obtain small time truncation error and can improve the computation accuracy.The parallel scheme is quite efficient in computing the saturation because the whole computation on a large-scaled domain is decomposed into several smallscaled computations on subdomains.(IV) It is easy to compute numerical simulation in parallel with high order accuracy on modern parallel machine.
j , b(c) = ϕ(X)c i {z 1 − 2 ∑ j=1 z j c j }, where c 1 , c 2 denote two different components of the saturation, and z 1 , z 2 denote the compressibility.Let c = c 1 = 1 − c 2 denote the first component of the saturation at production well.
n∆t≤T ||g n || X , ||g|| L2 (J;L 2 (Ω)) on p, c and their derivatives.Note.The method discussed above can be used by decomposing the whole domain Ω into several subdomains such as four subdomains in Fig. 2, Ω = 4 ∪ i=1 Ω i .Therefore, it is important in solving actual problems numerically.

Table 3
Domain decomposition with two subdomains and whole domain scheme have similar discretizations and their computational efficiency are compared.Computational time cost of different schemes under the fixed error level are shown in Table