Radial basis approximation of single-phase flow in porous media based on the Green’s functions

K. A. Sidelnikov1 , A. M. Gubanov2 , V. E. Lyalin3

1JSC “Izhevsk Petroleum Research Center”, Izhevsk, Russia

2Institute of Oil and Gas, Udmurt State University, Izhevsk, Russia

3Kalashnikov Izhevsk State Technical University, Izhevsk, Russia

3Corresponding author

Vibroengineering PROCEDIA, Vol. 7, 2016, p. 164-170.
Received 5 August 2016; accepted 8 August 2016; published 31 August 2016

Copyright © 2016 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Abstract.

The article discusses the problem of approximating solutions of differential equations describing the process of a two-dimensional fluid flow in porous media. The approximation is presented as a combination of radial basis functions on the basis of the Green’s function is used to solve the Poisson equation with variable coefficients in the case of steady state filtration and parabolic equations in the transient regime. To illustrate the effectiveness of the proposed approximation obtained by the field pressure distribution in the reservoir with a network of injection and production wells. Compare approximated pressure and design points to a satisfactory accuracy of the results.

Keywords: approximation, Green’s function, two-dimensional fluid flow, porous medium.

1. Approximation of the steady-state fluid flow in porous media

Let us consider the problem of approximating solutions of differential equations describing two-dimensional fluid flow through a porous medium. To describe steady-state two-dimensional fluid flow we study the following boundary value problem:

(1)
h r λ r P - γ Z r = - h r l = 1 L q l t δ r - r l ,

with Neumann boundary condition Pr/nΓ=0, where λr=kr/μ,γ=ρg,rΩΩ-2; Pr – pressure at point r=x,yΩ; Ω- – computational domain with boundary Γ; n – outer-pointing normal; k=diagk1,k2 – diagonal tensor of absolute permeability; h – reservoir thickness; Z – reservoir depth; μ – fluid viscosity; ρ – fluid density; g – gravity acceleration; xl – coordinates of lth well locations with flow rate ql, l=1,L¯; L – number of wells; δ – generalized two-dimensional Dirac delta function.

Eq. (1) is Poisson’s equation with variable coefficients. For ordinary Poisson’s equation:

(2)
2 Φ = 2 Φ x 2 + 2 Φ y 2 = - Q .

r = ( x , y ) D , where D denotes the entire plane, we can obtain the solution with Green’s function G(r,ρ)[1]:

(3)
Φ r = D G r , ρ Q ρ d σ ρ ,

where:

(4)
G r , ρ = 1 2 π l n 1 r - ρ .

Such representation of solution to Eq. (2) in the presence of point sources and sinks is examined in [2], and it is known as well interference problem. Under constant coefficients and isotropy in Eq. (1) the solution Eq. (3) is given by:

(5)
Φ r = k P μ = 1 2 π l = 1 L q l l n r - ρ l ,

where ql – well flow rate.

In case of variable coefficients Eq. (5) leads to significant errors in solution. We therefore propose the use of Green’s function Eq. (4) as radial basis functions Gr-ρ to approximate Eq. (1) based on their linear combination. Coefficients of linear combination can be calculated through the numerical solution obtained by finite difference method. Let us consider rectangular domain D=a,b;c,d which is divided into finite blocks Δx=b-a/nx, Δy=d-c/ny, where nx×ny – grid size.

Let us denote the numerical solution to Eq. (1) obtained by finite difference method as Pij, i=1,nx¯; j=1,ny¯, and its approximation as Pija, i=1,nx¯; j=1,ny¯. The approximation is constructed using linear combination of radial functions Gr-ρ. The radial function is undefined at a block contained source when r=ρ. Therefore, its value at this point is set equal to a constant determined by the size of the block [3], e.g. Gll=-lnΔxΔy/2π.

Then the approximate expression can be written as:

(6)
P i j a = P k a = l = 1 L α l δ l l G l l λ l + β l 1 - δ l k G k l λ k q l Δ x Δ y ,

where:

(7)
G k l = 1 2 π - l n ( Δ x Δ y ,           l = k , l n 1 r k - ρ l ,               l k ,           δ l k = 1 ,         l = k , 0 ,         l k .

The sources can be with different signs. Here it is assumed that injection well is a source with positive sign and production well is a source with negative sign. We introduce also index k which corresponds to the point rij, i=1,nx¯; j=1,ny¯; αl, βl – unknown coefficients that is determined through the condition k=1NPk-Pka2minα,β, where N=nxny.

To take into account Neumann boundary conditions Eq. (6) should be completed with additional image sources located outside of the computational domain. Each source ql is associated with four image sources qlm, l=1,L¯; m=1,4¯. Their coordinates are given by:

ρ l 1 = a - x l , y l ,         ρ l 2 = x l , 2 d - y l ,           ρ l 3 = 2 b - x l , y l ,         ρ l 4 = x l , c - y l ,

thus we can rewrite Eq. (6) as:

P k a = l = 1 L α l δ l l G l l λ l + β l 1 - δ l k G k l r k - ρ l λ k q l Δ x Δ y
            + m = 1 4 l = 1 L β l 1 - δ l k G k l r k - ρ l m λ k q l Δ x Δ y .

We use the following notation:

A = δ l l G l l λ l l ,           B = 1 - δ l k G k l r k - ρ l λ k ,         l = 1 , L ; ¯         k = 1 , L ¯ ,           d = P 1 P N ,           X = α 1 α L β 1 β L ,
C = 1 - δ l k G k l r k - ρ l λ k + m = 1 4 1 - δ l k G k l r k - ρ l m λ k ,           l = L , 2 L ¯ ,           k = L , N ¯ ,

where d – vector of known values; X – vector of unknown coefficients; C – calculated matrices. From A, B, C we construct block matrix:

G = A B 0 C = g k l ,         k = 1 , N ¯ ,           l = 1,2 L ¯ .

Here A, B are square matrices, C, G are rectangular matrices (N>2L).

According to the minimum condition of standard deviation we obtain matrix equation:

(8)
G X = d .

We multiply Eq. (8) from left by the transposed matrix GT: GTGX=GTd.

Then we calculate inverse of the matrix GTG and obtain expression to find unknown coefficients:

(9)
X = G T G - 1 G T d .

The operation GTG-1GT=G+ is called pseudoinversion [4], and the approximation based on radial basis function is known as radial basis neural network.

Then radial approximation appears as follows:

(10)
P a = G X .

Coefficients X are calculated for some preselected source locations and their flow rates. Then the approximation Eq. (10) is used for arbitrary source locations and their flow rates. The number of sources should be less than or equal to L. For a smaller number of sources, the flow rates of extra sources are set to zero. According to the Neumann boundary condition the material balance that accounts for source terms should be zero [5, 6].

For illustration, let us consider the following case study. The domain D=a,b;c,d, a= 0; b= 1000 m; c= 0; d= 1000 m was divided into 256×256 blocks of finite difference grid. Reservoir thickness was changed according to h(r)=10(1+0.1sin(4π(x+y))), fluid mobility was changed according to λ(r)=10-91+0.2sin(2πx)+0.15cos(2πy). There were 10 wells (5 production wells and 5 injection wells) in the domain. The absolute values of flow rates of all sources were the same and equal to Q= 50 m3/s. The sources were randomly located in the domain. To calculate pressure in blocks multigrid method with 7 nested grids was used.

Fig. 1 demonstrates distribution of the normalized pressure p=(P-Pmin)/(Pmax-Pmin) in the domain, where Pmax, Pmin – maximum and minimum pressure values. The distribution of the approximate pressure is shown in Fig. 2.

As shown in Fig. 3, comparison of the pressure values in midsection of the domain demonstrates relatively sufficient approximation.

Then we changed locations of the sources and their flow rates in random manner. Comparison of the numerical pressure values with their approximate values in midsection of the domain for different source distribution is shown in Fig. 4.

As can be seen, the radial network is well trained and provides sufficient approximation of the data that is independent of the training set. Fig. 5 demonstrates comparison of the numerical pressure values with their approximate values for all blocks of the domain.

The standard deviation of the approximate pressure values from their numerical values is 0.95 %. The plus signs in Fig. 5 depict pressure values in blocks contained sources. The maximum deviations are observed in these particular blocks where pressure changes more rapidly. The results show that the Neumann boundary conditions are satisfied. It should be noted that Neumann boundary conditions are the most difficult case study especially when they are applied for the entire boundary.

Fig. 1. Pressure distribution obtained by finite difference method

 Pressure distribution obtained  by finite difference method

Fig. 2. Approximate pressure distribution

 Approximate pressure distribution

Fig. 3. Comparison of the numerical pressure values with their approximate values in midsection of the domain

 Comparison of the numerical pressure  values with their approximate values  in midsection of the domain

Fig. 4. Comparison of the numerical pressure values with their approximate values in midsection of the domain in case of different source locations

 Comparison of the numerical pressure values with their approximate values in midsection of the domain in case of different source locations

Fig. 5. Comparison of the numerical pressure values with their approximate values for all blocks of the domain

 Comparison of the numerical pressure values  with their approximate values for all blocks of the domain

2. Approximation of the nonstationary fluid flow in porous media

For the case of compressible fluid flow and rock porosity depending on the pressure Eq. (1) becomes parabolic differential equation:

(11)
h r λ r P - γ Z r = h r β r P t - h r l = 1 L q l t δ r - r l ,

where λr=kr/μB;βx=φrcl/B°+φ°rcφ/B;ρ=ρ°1+clP-P°,B=B°/1+clP-P°,μ=μ°/1-cμP-P°,φr=φ°r1+cφP-P°, φ(r) – rock porosity; B – fluid formation volume factor; P° – reference pressure. Eq. (11) is complemented by the initial condition:

(12)
P 0 , r = P 0 r .

For parabolic model equation:

(13)
2 Φ = α 2 P t + f t , r ,         Φ 0 , r = Φ 0 r ,

there is Green’s function solution as well [1]:

(14)
Φ t , r = 0 t V G t , r ; τ , ξ f τ , ξ d τ d V ξ - α 2 V G t , r ; 0 , ξ Φ 0 ξ d V ξ ,

where G(t,r;τ,ξ) – Green’s function, which can be constructed from eigenfunctions ψk and eigenvalues λk of Helmholtz equation:

(15)
2 ψ r + λ ψ r = 0 ,         r V ,           G t , r ; τ , ξ = - 1 α 2 k e - λ k α 2 t - τ ψ - k ξ ψ k r .

Eq. (14) and Eq. (15) can be used as a basis for approximate solution to Eq. (11). It is important to note that Green’s function (15) can be represented as product of two functions (separation of variables): the first one depends exponentially only on time variable, the second one depends only on spatial variables.

Assume first that flow rate of sources is time-independent. Let us restrict ourselves to the one summand in Eq. (15). Then Eq. (15) can be rewritten as:

G t , r ; τ , ξ = - 1 a 2 e - λ α 2 t - τ Ψ ξ , r ,

where Ψ(ξ,r) – some function sought for approximation.

Time integration of Eq. (14) gives:

(16)
Φ t , r = 1 - e - λ α 2 t λ V Ψ ξ , r f ξ d V ξ + e - λ α 2 t V Ψ ξ , r Φ 0 ξ d V ξ .

It follows from Eq. (16) that solution to nonstationary equation Eq. (13) contains steady-state component:

Φ S r = V Ψ ξ , r f ξ d V ξ .

To which solution Eq. (16) converges asymptotically as t. The second term e-λt/α2VΨ(ξ,r)Φ0(ξ)dV(ξ) corresponds to the initial condition as t0. Therefore, approximation problem is solved in two steps. At first we seek an approximant to a steady-state solution as t. For this purpose the radial function Eq. (7) denoted as GS(ξ,r) is used. For the second term in Eq. (16) we take unit integral transform to satisfy the given initial condition. Then the parameter nonstationarity λ is adjusted. For best approximation on time several terms of series can be used in Eq. (15), so Green’s function can be presented as:

G t , r ; τ , ξ = - 1 a 2 k e - λ k α 2 t - τ Ψ ξ , r ,

and solution Eq. (16) can be written as:

(17)
Φ t , r = k 1 - e - λ k α 2 t λ k V G s ξ , r f ξ d V ξ + k e - λ k α 2 t Φ 0 r .

It should be noted that the underlying considerations are not strictly mathematically sound and needed only to choose basis functions when constructing the approximant.

Let us consider the case where l=1Lh(rl)ql=δq0. The steady-state solution to Eq. (11) as P/t0 doesn’t exist subject to the Neumann boundary condition. To construct function GS(ξ,r) we introduce background sources:

Δ q i j = - l = 1 L h r l q l b - a d - c Δ x i Δ y j ,         i = 1 , n x ¯ ,         j = 1 , n y ¯ ,

for each grid blocks of the computational domain. Now we can define the corresponding steady state and fit the function GS(ξ,r) which describes this state. Eq. (17) is now as follows:

Φ t , r = k 1 - e - λ k α 2 t λ k V G s ξ , r f ξ d V ξ + k e - λ k α 2 t Φ 0 r + δ q α 2 t .

If the flow rate of a source is a function of time f(t,r) then the nonstationary process is divided into time steps of length Δtn by analogy with finite-difference representation, so the approximate expression for pressure is written as:

(18)
P a t n + 1 , r
            = k 1 - e - λ k α 2 t n + 1 λ k V G s ξ , r f t n + 1 ξ ) d V ( ξ + k e - λ k α 2 t n + 1 P a t n , r + δ q α 2 t n + 1 ,

where α2=hrβr; f(t,rl)=hrl=1Lqltδr-rl.

For one summand in sum Eq. (18) can be written in the form:

P a t n , r i j = θ 0 t n , r i j + θ 1 t n , r i j e - λ α i j 2 t n .

According to the condition i,j,nPa(tn,rij)-P(tn,rij)2minλ we obtain the expression for coefficient λ based on the calculation of Pijn by finite-difference method:

λ = - 1 N n x n y i j n a i j 2 t n l n P i j n - θ 0 i j n θ 1 i j n

where Pijn – calculated pressure values; N – number of time steps. The results showed that in most practical cases λ0.25.

Let us consider the case study under the same conditions as in the steady-state case. In addition, we assumed that porosity φ°=0.15, and coefficients cφ= 10-4; cl=2∙10-4; cμ= 0. The stationary function GS(ρ,r) was taken as a result of previous training. The material unbalance was Q= 39 m3/s. The calculated value of λ was λ=0.27. Fig. 6 compares numerical and approximate pressures for different time steps in midsection of the domain.

The solid lines correspond to pressure obtained numerically by finite-difference method, and the circles correspond to approximate pressure. The comparison of the normalized numerical and approximate pressure values for all the grid blocks is shown in Fig. 7.

The standard deviation of the approximate pressure values from their numerical values is 1.84 %, that is approximately twice as high as for steady-state flow.

Fig. 6. The comparison of the numerical and approximate pressures in midsection of the domain in case of nonstationary flow 1 – t-=0.2; 2 – t-=0.4; 3 – t-=0.6; 4 – t-=0.8

The comparison of the numerical and approximate pressures in midsection of the domain  in case of nonstationary flow 1 – t-=0.2;  2 – t-=0.4; 3 – t-=0.6; 4 – t-=0.8

Fig. 7. Comparison of the normalized numerical and approximate pressures

 Comparison of the normalized  numerical and approximate pressures

3. Conclusions

1) The developed model of radial approximation of the steady-state equation of the fluid flow in porous media based on basis Green’s functions provides that standard deviation of the approximate pressure from its numerical values will be within the range 0.95-1.05 %. Such basis functions are used to approximate solution of the nonstationary multiphase flow in porous media as well.

2) To approximate solution of the nonstationary fluid flow in porous media Green’s function is used. It is represented as a product of two terms (separation of variables): the first one exponentially depends only on time; the second one depends only on spatial coordinates. The use of background sources allows the construction of approximants which are applicable for time-varying flow rates of the sources and unbalanced flow rates.

References

  1. Korn G. A., Korn T. M. Mathematical Handbook for Scientists and Engineers; Definitions, Theorems, and Formulas for Reference and Review. McGraw-Hill, New York, 1968, p. 1151. [Search CrossRef]
  2. Basniev K. S., Kochina I. N., Maksimov V. M. Underground Hydrodynamics. Nedra, Moscow, 1993, p. 416, (in Russian). [Search CrossRef]
  3. Kanevskaya R. D. Mathematical modeling of hydrodynamic processes of hydrocarbon deposit development. Izhevsk, 2003, p. 128, (in Russian). [Search CrossRef]
  4. Osowski S. Neural Networks for Information Processing. Finance and Statistics, Moscow, 2002, p. 344, (in Russian). [Search CrossRef]
  5. Chanson H. Applied Hydrodynamics: An Introduction to Ideal and Real Fluid Flows. CRC Press, Taylor & Francis Group, Leiden, The Netherlands, 2009, p. 478. [Search CrossRef]
  6. Economides Michael J., Nolte Kenneth G. Reservoir Stimulation. Wiley, USA, 2000, p. 856. [Search CrossRef]