1 Introduction

For the nonstationary flow computation of multilayer dynamics of fluids in porous media,we consider a kind of three-layer problems. The flows move horizontally at the first and third layers and move vertically at the middle layer (the weak percolation layer). Here,we need to obtain the solutions to the following three-dimensional convection-dominated diffusion nonlinear section coupled systems with initial-boundary values^{[1, 2, 3, 4, 5]}:

where

Ω ₀ and Ω₀ are the base lines of the section and the base region of the three-dimensional problem (see Fig. 1),respectively. ∂Ω and ∂Ω _i (i = 1,2,3) are the external boundaries of Ωand Ω_i, respectively.

The initial conditions are

The boundary conditions (the Dirichlet boundary conditions) are

and the interior boundary conditions are

During the simulation of multilayer dynamics of fluids in porous media,these unknown functions u,w,and v and the potential functions should be computed. Some other functions Δu,Δv,and ∂w / ∂z denote Darcy’s velocities,Φ_α (α= 1,2,3) represent the porosity,K₁(x,z,u,t), K₂(x,z,w,t),and K₃(x,z,v,t) are the stratigraphical permeabilities which are defined in the model,and a(x,z,u,t) = (a₁(x,z,u,t),a₂(x,z,u,t))T and b(x,z,v,t) = (b₁(x,z,v,t),b₂(x,z, v,t))^T,which are the convective coefficients,are given by the following conditions:

Q₁(x,z,t,u) and Q₃(x,z,t,v) are the external volumetric flow rates.

For convection-diffusion problems,especially those that are convection-dominated,Douglas and Russell^[6] have published some famous papers that discuss the characteristic finite difference method and the characteristic finite element method,which are used to overcome oscillation and faults likely arising in the traditional method^{[1, 7, 8, 9, 10]}. Later,Douglas^[7] considered the application of the characteristic finite difference method to solve incompressible two-phase dis- placement problems and obtained some convergence results and theoretical foundations by the way of energy mathematics. However,he did not obtain an optimal-order theorem in l2 norm and made it necessary that the problem is -periodic. The authors of this paper have extended Douglas’s results,obtained optimal-order error estimates in l2 norm,and used this method in the numerical simulation of oil deposit and the numerical simulation of enhanced oil recovery simulation^{[11, 12]} with moving boundary value conditions. For large-scale scientific and engi- neering computation (the mode number may amount up to tens of thousands or even hundreds of thousands),fractional step techniques are needed to convert a multi-dimensional problem into a series of successive one-dimensional problems^{[1, 10, 11, 13]}. Though Douglas and Peaceman introduced the alternating-direction method into the two-phase displacement problem,they did not obtain theoretical analysis results. By Fourier analysis,they succeeded in proving the results of stability and convergence of the constant coefficient case,but this method cannot be used for the common equations with variable coefficients^{[1, 14, 15]}. The authors have used the characteristic factional steps method in the compressible two-phase displacement problems and obtained optimal order error estimates in l² norm. However,the problems should be supposed to be -periodic in the numerical analysis^[16].

This paper,in the point of view of the practical exploration and development of oil-gas resources and groundwater percolation computation,discusses the nonstationary flow compu- tation of the multilayer ground percolation coupled system displacement problems. For the nonlinear section coupled system of multilayer fluid dynamics in porous media,a parallel algo- rithm derived by characteristic finite difference fractional steps is stated in this paper. Some techniques,such as calculus of variations,energy methods,twofold-quadratic interpolation of piecewise product type,multiplicative commutation rule of difference operators,decomposition of high order difference operators,and prior estimates are adopted to accomplish the theoretical analysis of optimal order in l² norm. This problem can be simplified into a two-dimensional expression. Thus,we suppose that the initial problem is linear. We have got some prelimi- nary results^{[17, 18]}. These methods have been used in numerical simulations of multilayer oil resources migration-accumulation problems. Thus,we have thoroughly and completely solved the well-known problem^{[1, 4, 5, 19, 20, 21, 22]}.

Generally speaking,this is a positive definite problem as follows:

where Φ_∗,Φ^∗,K_∗,K^∗,and A^∗ are constants,and the equality K^′ ₂(w) ≡ 0 holds nearby the planes z = H₁ and z = H₂.

Our assumption on the regularity of the solutions to (1)-(4) is denoted collectively by

and the functions Q₁(x,z,t,u) and Q₃(x,z,t,v) are Lipschitz continuous within a distance "0 from the solutions u and v. That is,when |ε_i| 6 ε₀ (1 6 i 6 4),we have

In this paper,M and " express the general positive constant and the general positive small constant,respectively,and they have different meanings at different places. 2 Second-order characteristic finite difference fractional step scheme

In order to get the solutions by the finite difference method,let the mesh region Ω 0,h replace Ω₀ . On the two-dimensional plane Oxy,let h₁ be the step length,x_i = ih₁,y_j = jh₁,and let Ω ₀,h replaceΩ ₀.

Let h₂ = (H₃ − H₂)/N₁ (N₁ is an integer) be the step length of the partition in the _z-axis,_zk = H₂ + kh₂,and tⁿ = nΔt. UseΩ ₁,_h to replace Ω ₁,

Let ∂Ω ₁,_h and ∂Ω ₀,_h stand for the boundaries ofΩ ₁,_h and Ω ₀,_h,respectively. Similarly, the length of the partition of Ω ₂ in the _z-direction is denoted by h₃ = (H₂ − H₁)/N₂,and zk = H₁ + kh₃. Use Ω ₂,_h to replace Ω ₂,

In Ω ₃,the length of the _z-direction is taken as h₄ = H₁/N₃,and ^z_k = kh₄. Use Ω ₃,_h to replace Ω ₃,

where N₂ and N₃ are integers. ∂Ω ₂,_h and ∂Ω ₃,_h stand for the external boundaries of Ω ₂,_h andΩ ₃,_h,respectively.

Let δ_x,δ_z and δ_x,δz be the forward and backward difference quotients,respectively. dtUn is the forward quotient of mesh function Uⁿ _ik.

and

As the flow is essentially in the characteristic direction,we use the modified method for the characteristic procedure to the first-order hyperbolic part of (1a),thus ensuring the high-order accuracy of numerical results^{[6, 7, 8, 9]}. Let

We write E_q. (1a) in the following form:

Approximate δ_uⁿ⁺¹/ δ_Τ1 = δ_u/ δ_Τ1 (X,tⁿ⁺¹) by a backward difference quotient in the _Τ1-direction,

For E_q. (1a),the characteristic finite difference fractional step scheme is given by

where we interpret Uⁿ(X) as the piecewise product twofold-quadratic interpolation^{[8, 9]}.

For E_q. (1b),the finite difference scheme is

Similarly,let

Then,we rewrite E_q. (1c) as follows:

Approximate δ_vⁿ⁺¹/ δ_Τ1 = δ_v/ δ_Τ1 (X,tⁿ⁺¹) by a backward difference quotient in the _Τ3-direction,

For E_q. (1c),the characteristic finite difference fractional step scheme is given by

where we interpret V ⁿ(X) as the piecewise product twofold-quadratic interpolation,

The algorithm runs a time interval as follows. Assume that the approximate solution {Uⁿ _ik,Wⁿ _ik,V ⁿ _ik} at the time level t = tⁿ is known. It is needed to find out the approximate solution {Uⁿ⁺¹ _ik ,Wⁿ⁺¹ _ik ,V ⁿ⁺¹ _ik } at the next time level t = tⁿ⁺¹. First,by (7a) and (7b),we can get the solution of transition sheaf {U ^n+1/2 _ik } by the method of speedup,take Uⁿ⁺¹ _i,₀ = Wⁿ _i,N2 , and obtain the solution {Uⁿ⁺¹ _ik } by (7c) and (7d). Next,by (9a),(9b),(9c),and (9d),we get the solution of transition sheaf {V ^n+1/2 _ik} and {V ⁿ⁺¹ _ik} in a similar way. By (8) and the interior boundary conditions (7e) and (9e),we can obtain {Wⁿ⁺¹ _ik}. Only in this way can we calculate continuously,so that a complete time step can be carried out. Finally,the positive definite condition guarantees that this difference solution exists,and it is unique. Computation runs for all j in the increasing order from j ₁( j ) to j ₂( j ) to produce the solutions in the section. Then, the whole approximate solution to the three-dimensional problem can be obtained.

Here,let

where Δt is suitably small. Using the condition of a(X,Uⁿ,tⁿ)|δΩ 1 = 0 and b(X,V ⁿ,tⁿ)|_{δΩ 3} = 0,we know that Y = _g(X) and Z = _f(X) are homeomorphism maps of Ω₁ into Ω₁ and Ω₃ into Ω₃ ,respectively. This can pledge that and still belong to Ω1 and Ω3,and the Ω-periodic condition can be ignored. 3 Convergence analysis

Let

Define the inner products on the net function space H_h ^{[19, 20]}. Considering a two-dimensional mesh region,

The norms on the discrete space are defined by

In general,for brevity,we omit the subscript symbol i of inner products and norms

Consider (7),(8),and (9) firstly. Let , ,and ,where u,v,and w are exact solutions to the problem (1)-(4),and U,V ,and W are the difference solutions to (7),(8) and (9),respectively.

Dispelling U^{n+1 2} of (7a)-(7e),we get the following equivalent forms:

Next to dispel V ^{n+1 2} of (9a)-(9e),we get the following equivalent forms: Then,we get the error equations of (1)-(4), where where where

Suppose that the space and time steps satisfy

First,consider the error equations (12),and interpret (X) as the piecewise product twofold-quadratic interpolation of value {}. Then,

where I is an identity operator,I₂ is a twofold-quadratic interpolation operator,and ˇ

By (12),(15),and (16),we can obtain

Test both sides of (17) against ,sum it up by parts,and note the deformation of the second term on the left-hand side,

We can obtain

Rewriting (18),we can obtain Estimate the first term on the right-hand side of (19). Let Supposing that (15) holds for the spacial and temporal discretizations,we have where

Estimate the third term on the right-hand side of (19),we can obtain

An induction hypothesis is given as follows: For the sixth term on the right hand side of (19),the operators − and − are self-conjugate,positive definite,and bounded,and the computational region is a regular cube,but generally their products vary in a commutative order. Note that , ,and the flow (the intermediate layer) only moves vertically. Then,we have

Owing to the induction hypothesis (23),it is sure that K₁(Uⁿ)_ik,K₁(Uⁿ)_i+1/2 ,_k,δxK₁(Un)ik, δz(Φ⁻¹ ₁,_ikK₁(Uⁿ)_{i+1/ 2} ,_k),· · · are bounded. For the first and second terms on the right hand side of (24),owing to the fact that K₁(Uⁿ) and Φ^-1 ₁ are positive definite and Cauchy inequalities, we can obtain

For the third term,we have Similarly,for the other terms,we can obtain

For the error equation (19),together with (20)-(26),we can obtain

Next,consider the error equation (13). Similarly,we have Testing both sides of (28) against and summing it up by parts, we have Note that Together with (30),we can rewrite (29) as follows: At last,consider the error equations (14). Testing both sides of (14) against δ_tωⁿ _ik = ωⁿ⁺¹ _ik −ωⁿ _ik = d_tωⁿ _ikΔt and summing it up by parts,we have Then,we have Adding (27),(31),and (33),summing over 0 ≤ n ≤ L,noting the initial conditions ω⁰ = 0 and the inequalities ,· · · ,we have where

By the discrete Gronwall inequality,we have

At last,it remains to check the induction hypothesis (23). As n = 0,we have = 0. Thus,the inductive assumption (23) is evidently correct. Assume (23) holds for n = 1,2,· · · ,L. Let n equal L + 1 and the spacial and temporal parameters satisfy (15). By (35),we have Then,the induction hypothesis (23) holds as n = L + 1. The proof of (23) is completed.

Theorem 1 Suppose that the exact solutions to the problems (1)-(4) satisfy the following smooth conditions:

Adopt the upwind finite difference fractional steps schemes (7),(8),and (9). Then,the following error estimate holds:

where

and the constant M is dependent of the solutions u,v and their derivatives,M^∗ = . 4 Actual applications

Recently,we have introduced these methods into the software system of multilayer oil re- sources migration-accumulation and oil resources estimation of Shengli oil field. The mathe- matical model is

where Ψ _o and Ψ _w are the potential functions,and k_ro and k_rw are the relative permeabilities for oil and water phases,respectively. s′ = ds/dpc,where s is the water concentration,and pc is the capillary pressure. In the simulations,the data of the related parameters can be referred in Table 1 4.1 Numerical simulation and analysis of Dongying hollow section of Shengli oil field

The flow of oil and water migrates and accumulates in the three-dimensional space which proceeds in different directions due to special environments such as geological structures. If the spacial region is divided into lots of sections in suitable direction,the secondary migration- accumulation problems can be similarly interpreted by the sections. Numerical simulation has been done on several sections in the middle Shasan segment of Dongying Hollow.

The industrial district of the middle Shasan region is illustrated by two kinds of vertical sections. One is measured in the east-west direction (named Section I),whose geodetic coor- dinate is from (20 621 246,4 151 110) to (20 641 246,4 151 110). The other is measured from(20 627 246,4 135 110) to (20 627 246,4 153 110) in the south-north direction (named Section II). Their geological structures are displayed in Fig. 3(a) and Fig. 3(b),respectively.

Then,the numerical simulation results are as follows. The computational parameters are 2 000 meters for the horizontal mesh step,30 for the number of mesh dots in the vertical partitions,and 1 000 years for the temporal step. The computation begins to simulate the migration-accumulation occurring twenty-five million years ago. The initial condition of water flow is assumed to be stationary,and the flow moves outwards on both sides. Numerical results on the 25 million years period include the potentials of water and oil and the concentration of water of different kinds of sections,whose isograms of the present concentration and potential of water and the potential of oil of the East-West sections and the South-North sections are shown in Fig. 4 and Fig. 5,respectively.

The analysis of numerical results is as follows.

(i) The numerical results are reasonable and give expression for the physical-mechanical properties of secondary migration-accumulation. That is to say,oil moves from the higher potential zone to the lower zone,accumulates and shapes a store in the lower potential zones of upper and local tuberositate regions.

(ii) The numerical results depend on the direction choice of sections,and incorrect data might be illustrated due to inappropriate plane cuts of the three-dimensional space.