Rock Flow Simulation by High-Order Quasi-Characteristics Scheme

. A pure second-order scheme of quasi-characteristics based on a pyramidal stencil is applied to the numerical modelling of non-stationary two-phase flows through porous media with the essentially heterogeneous properties. In contrast to well-known other high-resolution schemes with monotone properties, this scheme preserves a second-order approximation in regions, where discontinuities of solutions arise, as well as monotone properties of numerical solutions in those regions despite of well-known Godunov theorem. It is possible because the scheme under consideration is defined on a non-fixed stencil and is a combination of two high-order approximation scheme solutions with different dispersion properties. A special criterion according to which, one or another admissible solution is chosen, plays a key role in this scheme. A simple criterion with local character suitable for parallel computations is proposed. Some numerical results showing the efficiency of present approach in computations of two-phase flows through porous media with strongly discontinuous penetration coefficients are presented.


Introduction
In recent fifteen years many high-resolution numerical schemes modifying Godunov scheme have been proposed (see, for instance, [1][2][3][4][5]).However, the problem of development of high-order schemes with monotone properties in regions near the discontinuities of solutions remains in the focus of activities for many researchers in numerical methods for partial differential equations (PDE) and in computational fluid dynamics (CFD).According to the well-known Godunov's theorem, second-order explicit monotone schemes on the fixed stencils do not exist.Up to now, two different ways to resolve this restriction are known.The first one uses the idea of lowering the approximation order in the narrow regions near the discontinuities of solution.In fact, this approach has been realized in most of the modern high resolution schemes, because they set some restrictions on the recovery functions or limiters to provide the 198 monotone properties of solutions in zones where the discontinuities could arise.An excellent analysis of this approach is presented in [3,5].Therefore, most of highresolution schemes cited above are hybrid schemes, because their approximation order is lowered in the zones near the discontinuities.Various hybrid quasicharacteristics schemes for the solution of supersonic aerodynamics problems and two-phase porous media problems were developed and considered in [6][7][8][9][10][11].The second way consists in the construction of high-resolution schemes on the nonfixed stencils.For instance, one can apply two or more high-order schemes defined on different stencils and choose a final solution in each nodal point among admissible solutions to provide a monotone properties in regions where the discontinuities could arise.Such approach was considered in [12][13][14].In these articles, various quasi-characteristics schemes of the second-order were proposed and considered.All these schemes use a combination of two second-order approximation scheme solutions having the different dispersion properties.A special criterion according to which, one or another admissible solution is chosen between two admissible solutions to provide the monotone properties near the discontinuities, plays a key role in this scheme.In [12][13], a heuristic criterion based on the third-order theoretical estimation of the average value of the governing equation operator with respect to the grid cell is proposed.Unfortunately, it has a non-local and directed character and could not be easily adopted in multi-dimensional case and in parallel computations.In [14], simpler local and non-directed heuristic criterions suitable for parallel computations are proposed.As is shown in [12], the quasi-characteristics schemes are more accurate than fourth-order approximation schemes in computing of the initialvalue problems for the PDE of hyperbolic type, because the quasi-characteristics schemes are generalization of well-known back-ward characteristics schemes which are essentially more accurate in comparison with all other well-known numerical schemes.The reason of this consists in the naturally accurate treatment of the characteristic properties of the governing equations by the quasi-characteristics schemes in comparison with Godunov's type schemes based on the conservation laws treatment.Therefore, in recent years, various numerical schemes based on characteristics were proposed [15][16][17][18][19] for the solution of initial boundary value problems for reaction-diffusion equation and for correct setting of boundary conditions in decomposition of such problems.In this article, we consider the application of a new multidimensional scheme of quasi-characteristics to the numerical simulation of two-phase flows through porous media with strongly discontinuous penetration coefficients.This scheme approximates a transport equation in the system of the porous media equations on the pyramidal stencil without any splitting.A simple criterion suitable for the selection of final solution among two admissible solutions to provide the monotone properties of the final solution without spuriuos oscillations is proposed.Numerical results for various ratio of penetration coefficients are presented.These results show that the technique considered here could be efficiently used for the accurate numerical modelling of flows through the essentially heterogeneous porous media.

Governing equations, initial and boundary conditions
Let us consider the problem of a numerical simulation of two-phase flows through essentially heterogeneous porous medium with piece-wise constant absolute penetration factor.In the two-phase case, the governing equations [20] can be presented in the following form with respect to the water saturation and the pressure as unknown functions Here is a porosity factor, is an absolute penetration factor of porous medium, and are a relative penetration factors of water and oil, and are a viscosity of water and oil.Let us notice that the oil saturation can be evaluated by the water saturation according to the following simple formula .Since the relative penetration factors and are functions only of water saturation , then equation ( 1) can be presented as follows (3) Now we see that the system of governing equations ( 23) is of mixed type.Equation ( 3) is a nonlinear transport equation of hyperbolic type and the equation ( 2) is of elliptic type.Let us consider the transport equation (3) as a main governing equation and the elliptic equation ( 2) as a nonlinear restriction for coefficients of the main governing equation.Then we can apply the quasi-characteristics technique to solve the initial boundary-value problem for the transport equation and also on each time level we need to solve the boundary-value problem for elliptic equation to define the coefficients of the governing equation.In our approach, for the solution of the boundary-value problem for the elliptic equation we use the well-known five points finite difference conservative scheme and bi-conjugate gradient algorithm as in [7,11].

Let us consider rectangle flow regions divided into two subregions and
The absolute penetration factor in each subregion is a constant function, therefore in all region we have (4) Initial conditions for the transport equation ( 3) are (5) and the boundary conditions are ( 6)

Numerical scheme
In this section, we consider a non-splitting quasi-characteristics scheme on the pyramidal stencil applied to the solutions of the transport equation (3).Non-splitting scheme means that we do not use in our scheme any splitting technique for solution of the couple system of finite difference equations approximating the governing partial differential equation.It is very important in application to problems with heterogeneous coefficients, because in such problems sometimes splitting leads to the lowering of exactness of solutions.We develop this scheme with respect to the 3D transport equation written in the generalized form as follows (10) satisfying the following initial conditions (11) Here is a searching function and and are given.
In quasi-characteristics schemes [10], we approximate the governing equation written in the expanded characteristics form along some spatial grid lines (quasicharacteristics) in space as follows (12) Here is a total derivative of the searching function with respect to along line .As quasi-characteristics usually are used some grid lines belonging to the considering stencil.They should lie in close vicinity with respect to the characteristics of governing equation, and sometimes can coincide with them.Now we consider a uniform, for simplicity, in each direction finite-difference grid in space .We denote grid steps , and respectively.Let us consider a pyramidal stencil in the grid space.suppose that its basement belongs to some data layer and vertex belongs to the new layer .Coordinates of the above mentioned vertices are follows: , , , , .Also we take into consideration a center point of the basement of the pyramid stencil and denote the nodal points corresponding to the central points of the pyramid basement ribs as follows: , , , .
We suppose that the characteristics of the transport equation going through the point is lying inside the considering stencil and as a quasi-characteristics we can choose ribs of the pyramidal stencil.Then approximating the expanded characteristic form of the governing equation along these lines, we obtain 202 (13) where .

Fig.2. The pyramidal stencil
According to [10] we take the following approximation of the outward derivatives at the middle layer ( 14) Here we take a center point of the middle section of our stencil as the point (or ) and choose values and (or and ) at the middle layer by the formulas (16) (17) We denote as and the finite difference operators approximating the second order derivatives of the searching function with the second order approximation error on the appropriate stencil with middle point .
By substitution of relations (14-17) into ( 13), we obtain a system of four equations with respect to , , and .Solving it we obtain .In non-linear case, we need to do three iterations on nonlinear coefficients as is usually done in the method of characteristics.We call the scheme considered above scheme I. Now we construct a second scheme of the second order approximation with different dispersive properties in comparison with scheme I.For this purpose we choose point and values , according to the following formulas ( 18) By substitution of relations (14)(15) and ( 18) into ( 13), we obtain a system of four equations with respect to , , and .Solving it we obtain .As was mentioned above in non-linear case, it is necessary to do three iterations in evaluation of .We call this scheme scheme II.In [12][13] for 2D case, it was shown that by choosing one of two non-monotonous admissible solutions of the second order approximation with different dispersive properties, one can construct the final solution with monotone properties.As in the papers cited before, the criterion for the choosing the final solution is based on the analysis of the average value of the governing transport equation operator evaluated on each elementary mesh cell by the high order quadrature formulas.In this criterion, the history of computations in previous grid points is taken into account and therefore is not suitable for the parallel computations and in multi-dimensional case.For 2D case, a simpler heuristic criterion based on the minimization of the rough approximation of the average value of governing operator was proposed in [14].This criterion does not take into account the history of computations.It has a local character and it is suitable for the parallel realization.In this paper, we construct a simple heuristic criterion as a minimal principle for the increment of searching function over the stencil in following form (19) Here are some constants to be chosen.As our numerical tests show, the best result corresponds to the following set .Thus the final solution in each grid point is chosen among two admissible solutions and according to the following simple minimal principle (20) It is easy to see that this principle has a local character and it is very suitable for parallel computations, because it allows in principle to provide computations of searching function in each grid node independently in separate processors in computers with massive parallel processors and in computers with pipe-line processors it allows to provide the maximal loading of pipe-line.

Results of computations
Now let us consider some numerical results obtained by the proposed method.
Computations were carried out for the following values of parameters , , , , , , .Parameter varies in the range from to .Thus the absolute penetration in the subregion is 2 to 100 times less than those in subregion .Presented results correspond to the uniform grid with 61*61 nodal points in (x,y)-space.The first series of results corresponds to .Fig. 3 shows the isolines of water saturation and appropriate 3D chart for time hours.Fig. 4  and 5 show the same results for hours and hours respectively.Fig. 6 shows two functions characterizing the efficiency of the oil recovery process by the water drive.Line 1 corresponds to the ratio of the recovery oil to the total oil volume in initial moment with respect to time describing the water content in the development mixture at the production well corresponding to the boundary .According to the presented results we can see that the solution of the considering problem has a wave type and the front of water wave solution is spreading faster in the upper part of the flow region with high penetration.Subregion with low penetration plays the role of the partial obstacle and the water wave also is spreading in those region but more slowly.The second series of results corresponds to and the appropriate results are presented on Fig. 7-10.In this case, we can see that there are two shock-type water waves in the considering flow.The first wave corresponds to isolines 0.25 and 0.35 and the second corresponds to 0.45 and 0.55.According to the presented results it is easy to see that the first shock wave is spreading through the region with the low penetration, but the second wave stays near the right border of the low penetration subregion.The third series of results corresponds to and the appropriate results are presented on fig.11-14.In this case, the water is not spreading through the region with the low penetration and the water wave front is stopping near the right border of the low penetration subregion, which plays a role of a solid obstacle in the flow region.The analysis of the efficiency of the oil recovery process shows that after the moment hours, when the water wave in the upper part of the flow region is close to the production well (boundary ), the efficiency falls down and oil from the low penetration subregion and even from the high penetration subregion almost can not be developed by the water drive.According to our results, we see that in the case considered in this paper, it is possible to develop only about 35 percents of oil by the usual water drive technology although 70 percents of oil is contained in the high penetration subregion.These results are in good correspondence with well-known practice.

Conclusions
Our high-precision numerical quasi-characteristics technique developed for the transport equation of hyperbolic allows us to obtain solutions of complicated porous media problem with essentially heterogeneous parameters without mesh fitting procedures on rough spatial meshes.This technique can be implemented even on small computers and workstations for fast evaluation and exact modeling of oil and gas development technological processes.

Fig. 1 .
Fig.1.Flow region For the pressure equation (2) of the elliptic type we set a mixed Neumann and Dirichlet boundary conditions as follows