Published: 31 August 2016

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
Corresponding Author:
V. E. Lyalin
Views 30
Reads 12
Downloads 1167

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.

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
hrλrP-γZr=-hrl=1Lqltδr-rl,

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Φx2+2Φy2=-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=DGr,ρQρdσρ,

where:

4
Gr,ρ=12πln1r-ρ.

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=kPμ=12πl=1Lqllnr-ρ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
Pija=Pka=l=1LαlδllGllλl+βl1-δlkGklλkqlΔxΔy,

where:

7
Gkl=12π-ln(ΔxΔy, l=k,ln1rk-ρl, lk, δlk=1, l=k,0, lk.

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:

ρl1=a-xl,yl, ρl2=xl,2d-yl, ρl3=2b-xl,yl, ρl4=xl,c-yl,

thus we can rewrite Eq. (6) as:

Pka=l=1LαlδllGllλl+βl1-δlkGklrk-ρlλkqlΔxΔy
+m=14l=1Lβl1-δlkGklrk-ρlmλkqlΔxΔy.

We use the following notation:

A=δllGllλll, B=1-δlkGklrk-ρlλk, l=1,L;¯ k=1,L¯, d=P1PN, X=α1αLβ1βL,
C=1-δlkGklrk-ρlλk+m=141-δlkGklrk-ρlmλk, l=L,2L¯, 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=AB0C=gkl, k=1,N¯, l=1,2L¯.

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
GX=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=GTG-1GTd.

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
Pa=GX.

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. 1Pressure distribution obtained by finite difference method

Pressure distribution obtained  by finite difference method

Fig. 2Approximate pressure distribution

Approximate pressure distribution

Fig. 3Comparison 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. 4Comparison 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. 5Comparison 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
hrλrP-γZr=hrβrPt-hrl=1Lqltδr-rl,

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
P0,r=P0r.

For parabolic model equation:

13
2Φ=α2Pt+ft,r, Φ0,r=Φ0r,

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

14
Φt,r=0tVGt,r;τ,ξfτ,ξdτdVξ-α2VGt,r;0,ξΦ0ξdVξ,

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, rV, Gt,r;τ,ξ=-1α2ke-λkα2t-τψ-kξψkr.

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:

Gt,r;τ,ξ=-1a2e-λα2t-τΨξ,r,

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

Time integration of Eq. (14) gives:

16
Φt,r=1-e-λα2tλVΨξ,rfξdVξ+e-λα2tVΨξ,rΦ0ξdVξ.

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

ΦSr=VΨξ,rfξdVξ.

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:

Gt,r;τ,ξ=-1a2ke-λkα2t-τΨξ,r,

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

17
Φt,r=k1-e-λkα2tλkVGsξ,rfξdVξ+ke-λkα2tΦ0r.

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:

Δqij=-l=1Lhrlqlb-ad-cΔxiΔyj, i=1,nx¯, j=1,ny¯,

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=k1-e-λkα2tλkVGsξ,rfξdVξ+ke-λkα2tΦ0r+δqα2t.

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
Patn+1,r
=k1-e-λkα2tn+1λkVGsξ,rftn+1ξ)dV(ξ+ke-λkα2tn+1Patn,r+δqα2tn+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:

Patn,rij=θ0tn,rij+θ1tn,rije-λαij2tn.

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:

λ=-1Nnxnyijnaij2tnlnPijn-θ0ijnθ1ijn

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. 6The 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. 7Comparison 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

  • 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.
  • Basniev K. S., Kochina I. N., Maksimov V. M. Underground Hydrodynamics. Nedra, Moscow, 1993, p. 416, (in Russian).
  • Kanevskaya R. D. Mathematical modeling of hydrodynamic processes of hydrocarbon deposit development. Izhevsk, 2003, p. 128, (in Russian).
  • Osowski S. Neural Networks for Information Processing. Finance and Statistics, Moscow, 2002, p. 344, (in Russian).
  • Chanson H. Applied Hydrodynamics: An Introduction to Ideal and Real Fluid Flows. CRC Press, Taylor & Francis Group, Leiden, The Netherlands, 2009, p. 478.
  • Economides Michael J., Nolte Kenneth G. Reservoir Stimulation. Wiley, USA, 2000, p. 856.

About this article

Received
05 August 2016
Accepted
08 August 2016
Published
31 August 2016
SUBJECTS
Flow induced structural vibrations
Keywords
approximation
Green’s function
two-dimensional fluid flow
porous medium