Content of PetroWiki is intended for personal use only and to supplement, not replace, engineering judgment. SPE disclaims any and all liability for your use of such content. More information

# Solving unsteady flow problems with Green's and source functions

As discussed in Source function solutions of the diffusion equation, the conventional development of the source function solutions uses the instantaneous point-source solution as the building block with the appropriate integration (superposition) in space and time. In 1973, Gringarten and Ramey[1] introduced the use of the source and Green’s function method to the petroleum engineering literature with a more efficient method of developing the source solutions. Specifically, they suggested the use of infinite-plane sources as the building block with Newman’s product method.[2] In this page we discuss the use of Green’s functions and source functions in solving unsteady-flow problems in reservoirs.

## Green's functions and source functions in solving unsteady flow problems

Green’s function for transient flow in a porous medium is defined as the pressure at M (x, y, z) at time t because of an instantaneous point source of unit strength generated at point M′(x′, y′, z′) at time τ < t with the porous medium initially at zero pressure and the boundary of the medium kept at zero pressure or impermeable to flow.[1][3] If we let G(M, M′, tτ) denote the Green’s function, then it should satisfy the diffusion equation; that is,

....................(1)

Because G is a function of tτ, it should also satisfy the adjoint diffusion equation,

....................(2)

Green’s function also has the following properties: [1][3]

1. G is symmetrical in the two points M and M′; that is, Green’s function is invariant as the source and the observation points are interchanged. 2. As tτ, G vanishes at all points in the porous medium; that is, , except at the source location, M = M′, where it becomes infinite, so that G satisfies the delta function property,

....................(3)

where D indicates the domain of the porous medium, and φ(M) is any continuous function. 3. Because G corresponds to the pressure because of an instantaneous point source of unit strength, it satisfies

....................(4)

4. G or its normal derivative, ∂G/∂n, vanishes at the boundary, Γ, of the porous medium. If the porous medium is infinite, then G vanishes as M or M′→∞.

Let p(M′ , τ) be the pressure in the porous medium and G(M, M′, t - τ) be the Green’s function. Let D denote the domain of the porous medium. Then, p and G satisfy the following differential equations:

....................(5)

and

....................(6)

Then, we can write

....................(7)

or

....................(8)

where ε is a small positive number. Changing the order of integration and applying the Green’s theorem,

....................(9)

where D and Γ indicate the volume and boundary surface of the domain, respectively; S denotes the points on the boundary; and /∂n indicates differentiation in the normal direction of the surface Γ; we obtain

....................(10)

Taking the limit as ε→0 and noting the delta-function property of the Green’s function (Eq. 3), Eq. 10 yields

....................(11)

where pi(M) = p(M, t = 0) is the initial pressure at point M.

In Eq. 11, the boundary of the porous medium consists of two surfaces: the inner boundary that corresponds to the surface of the wellbore, Γw, and the outer boundary of the reservoir, Γe. Eq. 11 may be written as

....................(12)

As the fourth property of Green’s function noted previously requires, if the outer boundary of the reservoir is impermeable, or at infinity, then G vanishes at the outer boundary; that is, Ge) = 0. Thus, Eq. 12 becomes

....................(13)

Similarly, if the flux, , is specified at the inner boundary, then the normal derivative of Green’s function, , vanishes at that boundary. This yields

....................(14)

If the initial pressure, pi, is uniform over the entire domain (porous medium), D, then, by the third property of Green’s function (Eq. 4), we should have

....................(15)

Also, the flux law for porous medium (Darcy’s law) requires that the volume of fluid passing through the point, M′w, on the inner-boundary surface, Γw, at time t be equal to

....................(16)

The substitution of Eqs. 15 and 16 into Eq. 14 yields

....................(17)

Not surprisingly, Eq. 17 is the same as Eq. 18 because G in Eq. 17 is the instantaneous point-source solution of unit strength denoted by S in Eq. 18.

....................(18)

The expression given in Eq. 17 may be simplified further by assuming that the flux, , is uniformly distributed on the inner-boundary surface (wellbore), Γw. This yields

....................(19)

where,

....................(20)

The integration in the right side of Eq. 20 represents the distribution of instantaneous point sources over the length, area, or volume of the source (well), and S denotes the corresponding instantaneous source function. Source function solutions of the diffusion equation discusses the conventional derivation of the source functions starting from the basic instantaneous point-source solution. Here, we discuss the use of infinite-plane sources as the building block with Newman’s product method.[2]

Newman’s product method[2] may be stated for transient-flow problems in porous media as follows: [1] if a well/reservoir system can be visualized as the intersection of 1D or 2D systems, then the instantaneous source or Green’s function for this well/reservoir system can be constructed by the product of the source or Green’s functions for the 1D and/or 2D systems. For example, an infinite line-source well at x = x′, y = y′, and −∞ ≤ z′ ≤ +∞ in an infinite reservoir may be visualized as the intersection of two infinite, 1D plane sources; one at x = x′, −∞ ≤ y′ ≤ +∞, and −∞ ≤ z′ ≤ +∞, and the other at −∞ ≤ x′ ≤ +∞, y = y′, and −∞ ≤ z′ ≤ +∞. Then, the instantaneous source function for this infinite line-source well, S(x, x′, y, y′, tτ), may be obtained as the product of two infinite, 1D plane sources, given by

....................(21)

as follows

....................(22)

As expected, this solution is the same as Eq. 23, which was obtained by integration of an instantaneous point source in an infinite reservoir. For a radially isotropic reservoir (ηx = ηy = ηz), Eq. 22 yields

....................(23)

....................(24)

where d is the distance between the line source and the observation line in the x-y plane (see Fig. 1) and is given by

....................(25)

Similarly, intersecting three infinite instantaneous plane sources (or a line source and a plane source), we can obtain the instantaneous point-source solution in an infinite reservoir as

....................(26)

Instantaneous plane sources in slab reservoirs can be generated with the plane sources in infinite reservoirs and the method of images as discussed in Source function solutions of the diffusion equation. Similarly, the instantaneous slab sources can be obtained by integrating plane sources over the thickness of the slab source. Table 1, compiled from the work of Gringarten and Ramey,[1] presents the basic instantaneous source functions in infinite reservoirs, and Table 2 shows the corresponding geometries of the source-reservoir systems. The basic instantaneous source functions given in Table 2 may be used to construct the source functions that represent the appropriate well geometry by Newman’s product method.

Gringarten and Ramey[1] have also presented the approximating forms of the instantaneous linear sources and the time ranges for these approximations to be valid. The approximate solutions are very useful in obtaining expressions for pressure distributions at early and late times and identifying the flow regimes during the corresponding time periods. Table 3 presents the short- and long-time approximating forms for instantaneous linear sources and their time ranges. Examples 1 and 2 present some applications of the product-solution method and the derivation of the approximate solutions for pressure distributions.

## Example 1 - Ppartially penetrating vertical fracture in an infinite homogeneous slab reservoir

Considering transient flow toward a partially penetrating vertical fracture of vertical penetration hf and horizontal penetration 2xf in an infinite, homogeneous, slab reservoir of uniform thickness, h, and initial pressure, pi, with impermeable top and bottom boundaries.

Solution. Fig. 2 shows the geometry of the well reservoir system of interest. Approximate the fracture by a vertical plane of height hf and length 2xf. The corresponding source geometry may be visualized as the intersection of an infinite plane source at y = y′ in an infinite reservoir (Source I in Tables 1 and 2), an infinite-slab source of thickness 2xf at x = x′ in an infinite reservoir (Source IV), and an infinite-slab source of thickness hp = hf at z = zw in an infinite-slab reservoir of thickness h (Source VIII). Then, by Newman’s product method, the appropriate source function is given by

....................(27)

Assuming that the production is at a constant rate, and using Eq. 27 for the source function, S, in Eq. 19, we obtain

....................(28)

If the fracture penetrates the entire thickness of the reservoir (i.e., hf = h) as shown in Fig. 3, then Eq. 28 yields

....................(29)

The fully penetrating fracture solution given in Eq. 29 also could be obtained by constructing the source function as the product of an infinite plane source at y = y′ in an infinite reservoir (Source I in Tables 1 and 2) and an infinite-slab source of thickness 2xf at x = x′ in an infinite reservoir (Source IV). This source function then could be used in Eq. 19.

Fig. 4 presents an example of transient-pressure responses computed from Eq. 29. To obtain the results shown in Fig. 4, numerical integration has been used to evaluate the right side of Eq. 29. It is also of interest to obtain an early-time approximation for the solution given in Eq. 29. Substituting the early-time approximating forms for the slab sources in an infinite reservoir (approximations given in Table 3 for Source Functions IV and VIII), we obtain

....................(30)

where

....................(31)

and

....................(32)

Assuming a constant production rate, , and substituting the source function given by Eq. 30 in Eq. 19, we obtain

....................(33)

where erfc (z) is the complementary error function defined by

....................(34)

## Example 2 - Uniform-flux horizontal well in a closed, homogeneous, rectangular parallelepiped reservoir

Consider transient flow toward a uniform-flux horizontal well of length Lh located at (x′, y′, zw) in a closed, homogeneous, rectangular parallelepiped of dimensions 0 ≤ xxe, 0 ≤ yye, 0 ≤ zh and of initial pressure, pi.

Solution. Fig. 5 shows the sketch of the horizontal well/reservoir system considered in this example. If we approximate the horizontal well by a horizontal line source of length Lh, then the resulting source/reservoir system may be visualized as the intersection of three sources: an infinite plane source at y = y′ in an infinite-slab reservoir of thickness ye with impermeable boundaries (Source V in Tables 1 and 2), an infinite-slab source of thickness Lh at x = x′ in an infinite-slab reservoir of thickness xe with impermeable boundaries (Source VIII), and an infinite-plane source at z = zw in an infinite-slab reservoir of thickness h with impermeable boundaries (Source V). Then, by Newman’s product method, the appropriate source function can be obtained as

....................(35)

Assuming that the production is at a constant rate, , and using Eq. 35 for the source function, S, in Eq. 19, we obtain

....................(36)

Table 4 presents the pressure responses for a uniform-flux horizontal well in a closed square computed from Eq. 36. We may obtain a short-time approximation for Eq. 3.216 with the early-time approximations given in Table 3 for Source Functions V and VIII. This yields

....................(37)

where Ei(−x) is the exponential integral function defined by Eq. 38. Eq. 37 indicates that the early-time flow is radial in the y-z plane around the axis of the horizontal well. This solution corresponds to the time period during which none of the reservoir boundaries influence the pressure response.

....................(38)

It is also possible to obtain another approximation for Eq. 36 that covers the intermediate time-flow behavior. If we approximate the source function in the x direction (Source Function VIII) by its early and intermediate time approximation and the source function in the y direction (Source Function V) by its early time approximation given in Table 3, we obtain

....................(39)

This approximation should correspond to the time period during which the influence of the top and/or bottom boundaries may be evident but the lateral boundaries in the x and y directions do not have an influence on the pressure response. This solution also could have been obtained by assuming a laterally infinite reservoir. In this case, the source function would have been constructed as the product of three source functions: an infinite-plane source at y = y′ in an infinite reservoir (Source I in Tables 1 and 2), an infinite-slab source of thickness Lh at x = x′ in an infinite reservoir (Source IV), and an infinite-plane source at z = zw in an infinite-slab reservoir of thickness h with impermeable boundaries (Source V).

## Nomenclature

 c = fluid compressibility, atm−1 d = distance to a linear boundary, cm D = domain Ei(x) = exponential integral function G = Green’s function h = formation thickness, cm hf = fracture height (vertical penetration), cm k = isotropic permeability, md Lh = horizontal-well length, cm M = point in space M′ = source point in space Mw = point in Γw M′w = source point in Γw p = pressure, atm pi = initial pressure, atm q = production rate, cm3/s = instantaneous production rate for a point source, cm3/s r = radial coordinate and distance, cm r′ = source coordinate in r-direction, cm S = source function t = time, s x = distance in x-direction, cm x′ = source coordinate in x-direction, cm xe = distance to the external boundary in x-direction, cm xf = fracture half-length, cm y = distance in y-direction, cm y′ = source coordinate in y-direction, cm ye = distance to the external boundary in y-direction, cm z = distance in z-direction, cm z′ = source coordinate in z-direction, cm zw = well coordinate in z-direction, cm α = permeability direction, β = permeability direction, Γ = boundary surface, cm2 Γe = external boundary surface Γw = length, surface, or volume of the source Δ = difference operator η = diffusivity constant ηi = diffusivity constant in i direction, i = x, y, z, or r θ = angle from positive x-direction, degrees radian θ′ = source coordinate in θ-direction, degrees radian μ = viscosity, cp τ = time, s Φ = porosity, fraction φ(M) = any continuous function

## References

1. Gringarten, A.C. and Ramey Jr., H.J. 1973. The Use of Source and Green's Functions in Solving Unsteady-Flow Problems in Reservoirs. SPE J. 13 (5): 285-296. SPE-3818-PA. http://dx.doi.org/10.2118/3818-PA
2. Newman, A.B. 1936. Heating and Cooling Rectangular and Cylindrical Solids. Ind. Eng. Chem. 28 (5): 545–548. http://dx.doi.org/10.1021/ie50317a010
3. Carslaw, H.S. and Jaeger, J.C. 1959. Conduction of Heat in Solids, second edition, 353–386. Oxford, UK: Oxford University Press.