Published: 07 October 2016

Solutions of seepage equations in curvilinear coordinates

K. A. Sidelnikov1
A. M. Gubanov2
V. E. Lyalin3
M. A. Sharonov4
1JSC Izhevsk Petroleum Research Center, Izhevsk, Russia
2Institute of Oil and Gas, Udmurt State University, Izhevsk, Russia
3, 4Kalashnikov Izhevsk State Technical University, Izhevsk, Russia
Corresponding Author:
V. E. Lyalin
Views 18
Reads 7
Downloads 1244


The article provides a method for solving equations of filtration in curvilinear coordinate systems through the application of the algorithm for constructing orthogonal difference grids on curved surfaces, allowing cost-effective solution of the quasi-three-dimensional non-stationary filtration problem with the possibility of application of these solvers to optimization problems.

1. Introduction

In the case of complex shapes study area is desirable to use of computational grids borders, adapted to current conditions. The use of grids with rectangular border approximation leads to considerable complication of calculation algorithms to implement the boundary conditions.

2. Mathematical model

To build a grid on a curved surface, consider the geometry of a regular piece of the surface S defined by the equation in three-dimensional Euclidean space:

x=xu,v, y=yu,v, z=zu,v,

where u, v some arbitrary coordinates. The orthogonal grid on the surface given by the equations ξ=ξ(u,v)=const, η=η(u,v)=const. For orthogonal lines on the surface of the relation:

ηn=ξτ ,

where n, τ – normal and tangent to the line.

Functions ξ(u,v), η(u,v) must satisfy the system of equations [6]:

Gηv-FηuW=ξu, -Eηu-FηvW=ξv,

which is a generalization of the Cauchy-Riemann conditions.

To solve the system (3) differentiate the first equation by v, second by u, and add on them. As a result, we obtain a second differential operator:


where E(u,v)=xu2+yu2+zu2, G(u,v)=xv2+yv2+zv2, F(u,v)=xuxv+yuyv+zuzv – factors of fundamental quadratic form of the surface.

3. Algorithm and computational experiments

An algorithm for constructing orthogonal grid consists of consecutive solutions of the equation (4) and equation 2ξ=0 with the appropriate boundary conditions. For variable η at one boundary set values η=0, on the other η=1, on the other borders ηn=0, where n – normal to the boundary. From equation (2) are determined by the boundary conditions for the variable ξ and the equation 2ξ=0 is solved. Equations (1) and ξ=ξ(u,v)=const, η=η(u,v)=const define the position of the mesh nodes on the curved surface.

Equation (4) approximated the system of difference equations in the nine-pattern. The system of differential equations was solved by an iterative method of conjugate gradient regularization [1].

In the area shown in Fig. 1, a joint motion of the two phases: wetting and wetting liquids is considered. In this case, there are four unknown variables and hence four equations: two differential and two algebraic ones [2-4].

Po-Pw=Pc, Sw+So=1,

where f=w,o – the index for the wetting and non-wetting phase; Sf – saturation of the pore volume phase f; Pc – capillary pressure. It is believed that there is between Sf and Pc the functional dependence λfr=kfr/μf, γf=ρfg, rΩΩ-3; Pfr – pressure at the point rΩ; Ω- – filtering domain with boundary Γ; n – outward normal; kf – diagonal tensor absolute permeability; h – seam thickness; Z – depth; μf – liquid viscosity; ρ f – liquid density; g – acceleration of gravity; Qf – power sources; φr=φ°r1+cφP-P° ; φ(r) – porosity; Bf - volumetric coefficient; P – characteristic pressure.

Fig. 1Difference grid on curved surface

Difference grid on curved surface

In curvilinear orthogonal coordinates (ξ,η), the divergence of the gradient operator is as follows:


where e=xη2+yη2+zη2, g=xξ2+yξ2+zξ2, w=eg.

Then the difference approximation of two-dimensional filtering equation difference is a five-point pattern as in the Cartesian coordinates. The difference in rates is the difference in accounting factors 1w, ew, gw. Difference boundary conditions maintain their appearance. Therefore, the entire numerical algorithm, which is implemented into the Cartesian coordinates, is fully maintained.

Fig. 2Isolines S=const

Isolines S=const

In the following figures, the results of calculations for the reservoir, which has a curved shape, as in Figure 1 are shown. There is a problem of two-phase filtration. When t=3.2 contour lines S=const are shown (Figure 2), as well as the movement of water (Figure 3) and oil (Figure 4).

Fig. 3Water flow

Water flow

Fig. 6Oil flow

Oil flow

4. Conclusions

The developed approach to the construction of the difference of the orthogonal curvilinear grid allows the calculation of filtration processes in formations with complex geometry. Accounting reservoir capacity with a relatively small thickness allows solving problems in a quasi-three-dimensional setting.


  • Golub J., Van Loan Charles Matrix Computations (3rd Edition). Johns Hopkins University Press, Baltimore, MD, USA, 1996, p. 728.
  • Aziz H., Settarov E. Mathematical Modeling of Reservoir Systems. Institute of Computer Science, Moscow-Izhevsk, 2004, p. 416. Original Publication: “Nedra”, Moscow, 1982, (in Russian).
  • Chanson, H. Applied Hydrodynamics: An Introduction to Ideal and Real Fluid Flows. CRC Press, Taylor and Francis Group, Leiden, The Netherlands, 2009, p. 478.
  • Economides Michael J., Nolte Kenneth G. Reservoir Stimulation. Wiley, USA, 2000, p. 856.

About this article

18 August 2016
25 August 2016
07 October 2016
Mathematical models in engineering
two-phase filtration
curvilinear coordinates
orthogonal grid
oil reservoir