PEH:Mathematics of Transient Analysis

Publication Information

Petroleum Engineering Handbook

Larry W. Lake, Editor-in-Chief

Volume I – General Engineering

John R. Fanchi, Editor

Chapter 3 – Mathematics of Transient Analysis

Erdal Ozkan, Colorado School of Mines

Pgs. 77-172

ISBN 978-1-55563-108-6
Get permission for reuse

This chapter explains how fluid flow in porous media can be translated into a mathematical statement and how mathematical analysis can be used to answer transient-flow problems. This broad area is common to many other disciplines, such as heat conduction in solids and groundwater hydrology. The objective of this chapter is to introduce the fundamentals of transient analysis, present examples, and guide the interested reader to relevant references.

Introduction

Most physical phenomena in the domain of transient fluid flow in porous media can be described generally by partial differential equations (PDEs). With appropriate boundary conditions and sometimes with simplifying assumptions, the PDE leads to an initial boundary value problem (IBVP) that is solved to find a mathematical statement of the resulting flow in the porous medium. This section briefly discusses the statement of the IBVP for transient fluid flow in porous media.

Equations of Transient Fluid Flow in Porous Media

In essence, fluid motion in porous media can be specified by the knowledge of the velocity vector, , and the density of the fluid, ρ, as a function of the position (x, y, z) and time, t; that is, = (x, y, z, t) and ρ= ρ (x, y, z, t). Relative to the fixed Cartesian axes, the velocity vector can be written as

....................(3.1)

where v_x, v_y, and v_z are the velocity components, and , , and are the unit vectors in the x, y, and z directions, respectively.

The physical law governing the macroscopic fluid-flow phenomena in porous media is the conservation of mass, which states that mass is neither created nor destroyed. The mathematical formula of this rule is developed by considering the flow through a fixed arbitrary closed surface, Γ, lying entirely within a porous medium of porosity Φ, which is filled with a fluid of viscosity μ. Fig. 3.1 illustrates an arbitrary closed surface in porous medium.

Fig. 3.1 – Arbitrary closed surface Γ in porous medium.

The conservation of mass principle requires that the difference between the rates at which fluid enters and leaves the volume through its surface must equal the rate at which mass accumulates within the volume. The total mass within the volume at any time is given by

....................(3.2)

Then, the time rate of change of mass within Γ is

....................(3.3)

which, by the conservation of mass law, is equal to the rate at which mass enters V through the surface.

Consider the differential surface element, dΓ, shown in Fig. 3.1. The mass entering the volume through dΓ at the normal velocity, , in a time increment, Δt, is , and the total mass of the fluid passing through Γ during Δt is

....................(3.4)

The surface integral in Eq. 3.4 accounts for both influx and outflux through the surface of the volume; that is, ΔM_g is the difference between the masses entering and leaving the control volume during the time increment, Δt. Then, the mass rate entering the volume, V, through its surface, Γ, can be written as

....................(3.5)

By the principle of conservation of mass, equating the right sides of Eqs. 3.3 and 3.5 yields

....................(3.6)

A more useful relation is found with the divergence theorem, which states that the flux of through the closed surface, Γ, is identical to the volume integral of (the divergence of ) taken throughout V; that is,

....................(3.7)

Here, ∇ is the gradient operator, which in 3D Cartesian and cylindrical coordinates is given, respectively, by

....................(3.8)

and

....................(3.9)

With the relation in Eq. 3.7, Eq. 3.6 can be recast into

....................(3.10)

If the functions involved in the argument of the integral in Eq. 3.10 are continuous, then the integral is identically zero if and only if its argument is zero (because the volume integral in Eq. 3.10 is identically zero for any arbitrarily chosen volume). Then, the following continuity equation can be obtained.

....................(3.11)

Eq. 3.11 is a PDE that is equivalent to the statement of the conservation of mass for fluid flow in porous media. For practical purposes, however, Eq. 3.11 is expressed in terms of pressure because density and velocity cannot be measured directly. To express density, ρ, and velocity, , in terms of pressure, we use an equation of state and a flux law, known as Darcy’s law, respectively.

The following definition of isothermal fluid compressibility, c, is a useful equation of state that relates density to pressure.

....................(3.12)

If c is a constant (the compressibility of many reservoir liquids may be considered as constant), then Eq. 3.12 can be integrated to yield

....................(3.13)

where subscript 0 indicates the conditions at the datum. Similarly, the compressibility of the porous rock, c_f, is defined by

....................(3.14)

and the total system compressibility, c_t, is given by

....................(3.15)

These definitions of compressibility help recast Eq. 3.11 in terms of pressure.

Darcy’s law for fluid flow in porous media is a flux law. Neglecting the gravity effect, it is expressed by

....................(3.16)

In Eq. 3.16, μ is the viscosity of the fluid, and k is the permeability tensor of the formation given by

....................(3.17)

where α, β, and γ are the directions, and k_ij is the permeability in the i direction as a result of the pressure gradient in the j direction.

If Eqs. 3.13 through 3.16 are used in Eq. 3.11, an alternative statement of the conservation of mass principle for fluid flow in porous media is obtained:

....................(3.18)

Eq. 3.18 is the PDE that governs transient fluid flow in porous media. In the present form, Eq. 3.18 is not very useful in obtaining practical solutions because of the nonlinearity displayed in the second term of the left side. For liquid flow, the viscosity, μ, is constant and Eq. 3.18 can be linearized by assuming that the pressure gradients, ∇p, are small in the reservoir and the compressibility of the reservoir liquids, c, is on the order of 10⁻⁵ or smaller. Then, the second term of the left side of Eq. 3.18 may be neglected compared with the remaining terms and the following linear expression is obtained:

....................(3.19)

Eq. 3.19 (or Eq. 3.18) is known as the diffusivity equation. As an example in Cartesian coordinates, assuming that the coordinate axes can be chosen in the directions of the principal permeabilities, k, in Eq. 3.19, may be represented by the following diagonal tensor:

....................(3.20)

Then, Eq. 3.19 may be written

....................(3.21)

If each coordinate, j = x, y, or z, is multiplied by , where k may be chosen arbitrarily (to preserve the material balance, k is usually chosen to be ), Eq. 3.21 may be transformed into the diffusion equation for an isotropic domain:

....................(3.22)

where η is the diffusivity constant defined by

....................(3.23)

If the same transformation is also applied to the boundary conditions (see Sec. 3.1.2), the problems in anisotropic reservoirs may be transformed into those in isotropic reservoirs provided that the system is infinite or bounded by planes perpendicular to the principal axes of permeability. In all other cases, this transformation distorts the bounding surfaces.

For the flow of gases, the assumptions of small fluid compressibility and pressure gradient may not be appropriate and the c(∇p)² term in Eq. 3.18 may not be negligible. In these cases, an expression similar to Eq. 3.21 may be obtained from Eq. 3.18 in terms of pseudopressure, m, as

....................(3.24)

Here, the pseudopressure is defined by^[1]

....................(3.25)

where Z is the compressibility factor. To define a complete physical problem, Eq. 3.21 (or 3.24) should be subject to the initial and boundary conditions discussed in Sec. 3.1.2.

Initial and Boundary Conditions

The solution of the diffusivity equation (Eq. 3.19) should satisfy the initial condition in the porous medium. The initial condition is normally expressed in terms of a known pressure distribution at time zero; that is,

....................(3.26)

The most common initial condition is the uniform pressure distribution in the entire porous medium; that is, f (x, y, z) = p_i.

The boundary conditions are specified at the inner (wellbore) and outer boundaries of the reservoir. These are usually in the form of prescribed flux or pressure at the boundaries. The condition of prescribed flux can be formulated as

....................(3.27)

where Γ is the surface of the boundary, and n indicates the outward normal direction of the boundary surface. The prescribed flux condition may be used at the inner and outer boundaries of the reservoir. The most common use of the prescribed flux condition at the inner boundary is for the production at a constant rate. In this case, the function, g(t), is related to a constant production rate, q. At the outer boundary, the prescribed flux condition is usually used to indicate impermeable boundaries [g(t)=0] and leads to a pseudosteady state under the influence of boundaries.

For some applications, pressure may be specified at the inner and outer boundaries. In this case,

....................(3.28)

When used at the inner boundary, this condition represents production at a constant pressure, p_wf; that is, h(t) = p_wf. At the outer boundary, specified pressure, p_e, is usually a result of injection or influx from an adjacent aquifer, which usually leads to steady state in the reservoir.

It is also possible to have boundary conditions of mixed type. These usually correspond to interface conditions in porous media. Raghavan^[2] contains more details about the common boundary conditions for the diffusion equation.

Assumptions and Limits

Some assumptions have been made in the derivation of the diffusivity equation given by Eq. 3.19. These assumptions determine the limits of applicability of the solutions obtained from Eq. 3.19. One of the most important assumptions involved is the continuity of the properties involved in Eq. 3.19. (This was required to obtain Eq. 3.19 from the more general integral form in Eq. 3.10.) Therefore, sharp changes in the properties of the reservoir rock and fluid (such as faults and fluid banks) should be incorporated in the form of boundary or interface conditions in the solution of Eq. 3.19.

The second important assumption is that Darcy’s law describes the flux in porous media. This assumption is valid at relatively low fluid velocities that may be appropriate to describe liquid flow. At high velocities (when Reynolds number based on average sand grain diameter approaches unity) such as those observed in gas reservoirs, Darcy’s law is not valid.^[3] In this case, Forchheimer’s equation,^[4] which accounts for the inertial effects, should be used. In petroleum engineering, it is a common practice to consider the additional pressure drop as a result of non-Darcy flow in the form of a pseudoskin because it is usually effective in a small vicinity of the wellbore. Therefore, in this chapter, we do not consider non-Darcy flow in the reservoir.

Bessel Functions

As Sec. 3.3 illustrates, the Laplace transform of the diffusion equation in radial coordinates yields a modified Bessel’s equation, and its solutions are obtained in terms of modified Bessel functions. This section introduces Bessel functions and discusses some of their properties to the extent that they are encountered in the solutions of more common petroleum engineering problems.

Preliminary Definitions

A differential equation of the type

....................(3.29)

is called a Bessel’s equation of order v. A solution of Bessel’s equation of order v is called a Bessel function of order v. A differential equation of the type

....................(3.30)

is called a modified Bessel’s equation of order v. Eq. 3.30 is obtained by substituting λz for z in Eq. 3.29. Of particular interest is the case in which λ=ki so that Eq. 3.30 becomes

....................(3.31)

Eq. 3.31 is called the modified Bessel’s equation of order v. A solution of the modified Bessel’s equation of order v is called a modified Bessel function of order v.

Solutions of Bessel’s Equations and Bessel Functions

There are many methods of obtaining or constructing Bessel functions.^[5] Only the final form of the Bessel functions that are of interest are presented here.

If v is not a positive integer, then the general solution of Bessel’s equation of order v (Eq. 3.29) is given by

....................(3.32)

where A and B are arbitrary constants, and J_v(z) is the Bessel function of order v of the first kind given by

....................(3.33)

In Eq. 3.33, Γ(x) is the gamma function defined by

....................(3.34)

If v is a positive integer, n, then J_v and J−_v are linearly dependent, and the solution of Eq. 3.29 is written as

....................(3.35)

In Eq. 3.35, Y_n(z) is the Bessel function of order n of the second kind and is defined by

....................(3.36)

Similarly, if v is not a positive integer, the general solution of the modified Bessel’s equation of order v (Eq. 3.31) is given by

....................(3.37)

where I_v(z) is the modified Bessel function of order v of the first kind defined by

....................(3.38)

If v is a positive integer, n, I_v, and I_−v are linearly dependent. The solution for this case is

....................(3.39)

where K_n(z) is the modified Bessel function of order n of the second kind and is defined by

....................(3.40)

The modified Bessel functions of order zero and one are of special interest, and Sec. 3.2.3 discusses some of their special features.

Modified Bessel Functions of Order Zero and One

Modified Bessel functions of order zero and one are related to each other by the following relations:

....................(3.41)

and

....................(3.42)

Fig. 3.2 shows these functions graphically.

Fig. 3.2 – Bessel functions of order zero and one.

For small arguments, the following asymptotic expansions may be used for the modified Bessel functions of order zero and one:^[5]

....................(3.43)

....................(3.44)

....................(3.45)

where γ = 0.5772…, and

....................(3.46)

Also, for large arguments, the following relations may be useful:

....................(3.47)

for |arg z| < π / 2, and

....................(3.48)

for |arg z| < 3π / 2. On the basis of the relations given by Eqs. 3.43 through 3.48, the following limiting forms may be written:

....................(3.49)

....................(3.50)

....................(3.51)

....................(3.52)

....................(3.53)

....................(3.54)

....................(3.55)

....................(3.56)

and

....................(3.57)

These relations are useful in the evaluation of the asymptotic behavior of transient-pressure solutions.

Laplace Transformation

Integral transforms are useful in solving differential equations. A special form of the linear integral transforms, known as the Laplace transformation, is particularly useful in the solution of the diffusion equation. The Laplace transformation of a function, F(t), denoted by L{F(t)}, is defined by

....................(3.58)

where s is a number whose real part is positive and large enough for the integral in Eq. 3.58 to exist. In this chapter, a bar over the function indicates the image or the Laplace transform of the function; that is,

....................(3.59)

Fundamental Properties of the Laplace Transformation

The following fundamental properties of the Laplace transformation are useful in the solution of common transient-flow problems.

Transforms of Derivatives.

....................(3.60)

....................(3.61)

....................(3.62)

Transforms of Integrals.

....................(3.63)

Substitution.

....................(3.64)

....................(3.65)

where .

Translation.

....................(3.66)

where H(t - t₀) is Heaviside’s unit step function defined by

....................(3.67)

Convolution.

....................(3.68)

Inverse Laplace Transformation and Asymptotic Forms

For the Laplace transform to be useful, the inverse Laplace transformation must be uniquely defined. L⁻¹ denotes the inverse Laplace transform operator; that is,

....................(3.69)

In this operation, p(t) represents the inverse (transform) of the Laplace domain function, . A uniqueness theorem of the inversion guarantees that no two functions of the class ε have the same Laplace transform.^[6] The class ε is defined as the set of sectionally continuous functions F(t) that are continuous on each bounded interval over the half line t > 0 except at a finite number of points, t_i, where they are defined by

....................(3.70)

and | F(t) | < Me^αt for any constants M and α.

The most rigorous technique to find the inverse Laplace transform of a Laplace domain function is the use of the inversion integral,^[6] but its discussion is outside the scope of this chapter. For petroleum engineering applications, a simple table look-up procedure is usually the first resort. Table 3.1 shows an example table of Laplace transform pairs that may be used to find the Laplace transforms of real-space functions or the inverse Laplace transforms of the Laplace domain functions. Fairly large tables of Laplace transform pairs can be found in a couple of sources.^[6]^[7] The relations given in the Laplace transform tables may be extended to more complex functions with the fundamental properties of the Laplace transforms noted in Sec. 3.3.1.

Table 3.1

When a simple analytical inversion is not possible, numerical inversion of a Laplace domain function is an alternate procedure. Many numerical inversion algorithms have been proposed in the literature. For the inversion of the transient-flow solutions in Laplace domain, the numerical inversion algorithm suggested by Stehfest^[8] is the most popular algorithm.

The Stehfest algorithm is based on a stochastic process and suggests that an approximate value, p_a(T), of the inverse of the Laplace domain function, , may be obtained at time t = T by

....................(3.71)

where

....................(3.72)

In Eqs. 3.71 and 3.72, N is an even integer. Although, theoretically, the accuracy of the inversion should increase as N tends to infinity [p_a (T) should tend to p(T)], the accuracy may be lost because of round-off errors when N becomes large. Normally, the optimum value of N is determined as a result of a numerical experiment. As a reference, however, the range of 6 ≤ N ≤ 18 covers the most common values of N for transient-flow problems. The Stehfest algorithm is not appropriate for the numerical inversion of oscillatory and discontinuous functions. In these cases, a more complex algorithm proposed by Crump^[9] may be used.

In some cases, obtaining asymptotic solutions for small and large values of time may be of interest. These asymptotic results may be obtained without inverting the full solution into the real-time domain. The limiting forms of the full solution as s → ∞ and s → 0 correspond to the limiting forms in the time domain for short and long time, respectively. The inversion of the limiting forms may be easier than the inversion of the full solution. Examples 3.1 through 3.4 demonstrate the use of Laplace transformation in the solution of transient-flow problems.

Example 3.1

Consider transient flow toward a fully penetrating vertical well in an infinite homogeneous reservoir of uniform thickness, h, and initial pressure, p_i.

Solution. This problem may be formulated most conveniently in the radial coordinates. The diffusivity equation governing fluid flow in porous media is given, in radial coordinates, by

....................(3.73)

where ∆p = p_i – p. Eq. 3.73 is the same in absolute (cgs or SI) or Darcy units. (In field units, some conversion coefficients are involved in Eq. 3.73.) The initial condition is

....................(3.74)

which means that the pressure is uniform and equal to p_i initially throughout the reservoir. The outer boundary condition for an infinite reservoir is

....................(3.75)

which physically means that for any given time, t, there is a large enough distance, r, in the reservoir at which the initial pressure, p_i, has been preserved.

The inner boundary condition depends on the production conditions at the surface of the wellbore (r = r_w). Assuming that the well is produced at a constant rate, q, for all times,

....................(3.76)

The inner boundary condition given in Eq. 3.76 is simply a restatement of the flux law (Darcy’s law given by Eq. 3.16) at the surface of the wellbore.

Eqs. 3.73 through 3.76 define the IBVP to be solved to obtain the transient-pressure distribution for the given system. Application of the Laplace transforms to Eq. 3.73 yields

....................(3.77)

or, rearranging, we obtain

....................(3.78)

In obtaining the right side of Eq. 3.77, the initial condition (Eq. 3.74) has been used. Similarly, Eqs. 3.75 and 3.76 are transformed into the following forms, respectively.

....................(3.79)

and

....................(3.80)

Comparing Eq. 3.78 with Eq. 3.31, we recognize Eq. 3.78 as the modified Bessel’s equation of order zero. The solution of Eq. 3.78 may be written directly from Eq. 3.39 as

....................(3.81)

The constants C₁ and C₂ in Eq. 3.81 are obtained from the boundary conditions. The outer boundary condition (Eq. 3.79) indicates that C₁ = 0 [because , Eq. 3.79 is satisfied only if C₁ = 0]; therefore,

....................(3.82)

From Eqs. 3.80 and 3.82, we obtain

....................(3.83)

which yields

....................(3.84)

Then, the solution for the transient-pressure distribution is given, in the Laplace transform domain, by

....................(3.85)

To complete the solution of the problem, Eq. 3.85 should be inverted into the real-time domain. The real inversion of Eq. 3.85, however, is not available in terms of standard functions. One option is to use Stehfest’s numerical inversion algorithm^[8] as discussed in Sec. 3.3.2. The dashed line in Fig. 3.3 represents the numerical inversion of the solution in Eq. 3.85. Another option is to find an approximate inversion. One of these asymptotic forms is known as the line-source solution and commonly used in transient-pressure analysis.

Fig. 3.3 – Finite wellbore radius (Eq. 3.85) and line-source (Eq. 3.87 or 3.91) solutions for Example 3.1.

To obtain the line-source approximation of the solution given in Eq. 3.85, we assume that the radius of the wellbore is small compared with the other dimensions of the reservoir. Thus, if we assume r_w→0 and use the relation given in Eq. 3.56, we obtain

....................(3.86)

Using this relation in Eq. 3.85, we obtain the line-source solution in Laplace domain as

....................(3.87)

The inversion of Eq. 3.87 can be accomplished by using a Laplace transform table. From Table 3.1 (or from the tables in two sources^[6]^[7]), we have

....................(3.88)

With Eq. 3.88 and the Laplace transform property noted in Eq. 3.63, we obtain the following inversion of Eq. 3.87 in the real-time domain:

....................(3.89)

Making the substitution u = r² / (4ηt′) and noting the definition of the exponential integral function, Ei(x), given by

....................(3.90)

we obtain the line-source solution as

....................(3.91)

Fig. 3.3 shows a comparison of the results computed from Eq. 3.85 (finite-wellbore radius) and Eq. 3.91 (line source) for the data noted in the figure. The two solutions yield different results at early times but become the same at later times. In fact, it can be shown analytically that the long-time approximation of the finite-wellbore radius solution (Eq. 3.85) is the same as the line-source well solution. To show this, we note that the long-time approximation of the solution in the time domain corresponds to the limiting form of the solution in the Laplace domain as s → 0. Then, with the property of the Bessel function given in Eq. 3.56, we can show that

....................(3.92)

Example 3.2

Consider transient flow as a result of constant-rate production from a fully penetrating vertical well in a closed cylindrical reservoir initially at uniform initial pressure, p_i.

Solution. Fluid flow in cylindrical porous media is described by the diffusion equation in radial coordinates given by

....................(3.93)

The initial condition corresponding to the uniform pressure distribution equal to p_i is

....................(3.94)

and the inner boundary condition for a constant production rate, q, for all times is

....................(3.95)

The closed outer boundary condition is represented mathematically by zero flux at the outer boundary (r = r_e) that corresponds to

....................(3.96)

The Laplace transforms of Eqs. 3.93 through 3.96 yield, respectively,

....................(3.97)

....................(3.98)

and

....................(3.99)

(The initial condition given by Eq. 3.94 has been used to obtain Eq. 3.97.) Because Eq. 3.97 is the modified Bessel’s equation of order zero, its general solution is given by

....................(3.100)

With the outer boundary condition given by Eq. 3.99, we obtain

....................(3.101)

which yields

....................(3.102)

and thus

....................(3.103)

Using the inner boundary condition given by Eq. 3.98 yields

....................(3.104)

From Eqs. 3.102 and 3.104, we obtain the coefficients C₁ and C₂ as follows:

....................(3.105)

and

....................(3.106)

Substituting C₁ and C₂ into Eq. 3.100 yields

....................(3.107)

The inverse of the solution given by Eq. 3.107 may not be found in the Laplace transform tables. van Everdingen and Hurst^[10] provided the following analytical inversion of Eq. 3.107 with the inversion integral.

....................(3.108)

In Eq. 3.108, β₁, β₂, etc. are the roots of

....................(3.109)

The solution given in Eq. 3.107 may also be inverted numerically with the Stehfest algorithm.^[8] Fig. 3.4 shows the results of the numerical inversion of Eq. 3.107.

Fig. 3.4 – Bounded reservoir solution (Eq. 3.107) for Example 3.2.

Example 3.3

Consider the flowing wellbore pressure of a fully penetrating vertical well with wellbore storage and skin in an infinite reservoir.

Solution. Revisit the case in Example 3.1 and add the effect of a skin zone around the wellbore. Assume that the constant production rate is specified at the surface so that the storage capacity of the wellbore needs to be taken into account. Before presenting the initial-boundary value problem, skin factor and surface production rate should be defined.

Using van Everdingen and Hurst’s thin-skin concept^[10] (vanishingly small skin-zone radius), the skin factor is defined by

....................(3.110)

where q_sf is the sandface production rate, p(r_w +) denotes the reservoir pressure immediately outside the skin-zone boundary, and p_wf is the flowing wellbore pressure measured inside the wellbore. Rearranging Eq. 3.110, we obtain the following relation for the flowing wellbore pressure.

....................(3.111)

When the production rate is specified at the surface, it is necessary to account for the fact that the wellbore can store and unload fluids. The surface production rate, q, is equal to the sum of the wellbore unloading rate, q_wb, and the sandface production rate, q_sf; that is,

....................(3.112)

where

....................(3.113)

and

....................(3.114)

In Eq. 3.113, C is the wellbore-storage coefficient. Substituting Eqs. 3.113 and 3.114 into Eq. 3.112, we obtain the following expression for the surface production rate.

....................(3.115)

The mathematical statement of the problem under consideration is similar to that in Example 3.1, except that the inner-boundary condition should be replaced by Eq. 3.115, and Eq. 3.111 should be incorporated to account for the skin effect. The IBVP is defined by the following set of equations in the Laplace domain:

....................(3.116)

....................(3.117)

....................(3.118)

and

....................(3.119)

The general solution of Eq. 3.116 is

....................(3.120)

The condition in Eq. 3.117 requires that C₁ = 0; therefore,

....................(3.121)

The use of Eq. 3.121 in Eq. 3.119 yields

....................(3.122)

From Eqs. 3.118, 3.121, and 3.122, we obtain

....................(3.123)

which yields

....................(3.124)

Substituting Eq. 3.124 for C₂ in Eq. 3.122, we obtain the solution for the transient-pressure distribution in the Laplace transform domain as

....................(3.125)

The real inversion of the solution in Eq. 3.125 has been obtained by Agarwal et al.^[11] with the inversion integral. It is also possible to invert Eq. 3.125 numerically. Fig. 3.5 shows the results of the numerical inversion of Eq. 3.125 with the Stehfest’s algorithm.^[8] Also shown in Fig. 3.5 are the logarithmic derivatives of Δp_wf. These derivatives are computed by applying the Laplace transformation property given in Eq. 3.60 to Eq. 3.125 as follows:

....................(3.126)

Here, we have used Δp_{w f}(t = 0) = 0. To obtain the logarithmic derivatives, we simply note that

....................(3.127)

Fig. 3.5 – Wellbore-storage and skin solution (Eq. 3.125) for Example 3.3.

Example 3.4

Consider pressure buildup with wellbore storage and skin following a drawdown period at a constant rate in an infinite reservoir.

Solution. This example is similar to Example 3.3 except, at time t_p, the well is shut in and pressure buildup begins. The system of equations to define this problem is

....................(3.128)

....................(3.129)

....................(3.130)

....................(3.131)

where H(t - t_p) is Heaviside’s unit function (Eq. 3.67), and

....................(3.132)

The right side of the boundary condition in Eq. 3.131 accounts for a constant surface production rate, q, for 0 < t < t_p and for shut in (q = 0) for t > t_p. Taking the Laplace transforms of Eqs. 3.128 through 3.132, we obtain

....................(3.133)

....................(3.134)

....................(3.135)

and

....................(3.136)

The general solution of Eq. 3.133 is

....................(3.137)

The condition in Eq. 3.134 requires that C₁= 0; therefore,

....................(3.138)

From Eqs. 3.138 and 3.136, we obtain

....................(3.139)

Substituting the results of Eqs. 3.138 and 3.139 into Eq. 3.135, we have

....................(3.140)

which yields

....................(3.141)

Substituting Eq. 3.141 into Eq. 3.139, we obtain the following solution in the Laplace transform domain, which covers both the drawdown and buildup periods.

....................(3.142)

The term contributed by the discontinuity at time t = t_p causes difficulties in the numerical inversion of the right side of Eq. 3.142 with the use of the Stehfest algorithm.^[8] As suggested by Chen and Raghavan,^[12] this problem may be solved by noting that

....................(3.143)

and applying the Stehfest algorithm term by term to the right side of Eq. 3.143. Fig. 3.6 shows sample results obtained by the numerical inversion of Eq. 3.142.

Fig. 3.6 – Drawdown and buildup results with wellbore-storage and skin solution (Eq. 3.142) for Example 3.4.

Green’s Functions and Source Functions

Green’s function and source functions are used to solve 2D and 3D transient-flow problems that may result from complex well geometries, such as partially penetrating vertical and inclined wells, hydraulically fractured wells, and horizontal wells. Before introducing these techniques, it is useful to clarify the terminology.

In our terminology, a source is a point, line, surface, or volume at which fluids are withdrawn from the reservoir. Strictly speaking, fluid withdrawal should be associated with a sink, and the injection of fluids should be related to a source. Here, however, the term source is used for both production and injection with the convention that a negative withdrawal rate indicates injection.

Green’s functions and source functions are closely related. A Green’s function is defined for a differential equation with specified boundary conditions (prescribed flux or pressure) and corresponds to an instantaneous point-source solution. A source function, on the other hand, is the solution of the given differential equation with specified boundary conditions and source geometry.

The details of the theory and application of Green’s function and source functions for the solution of transient-flow problems in porous media can be found in many sources.^[2]^[13]^[14]^[15]^[16]^[17]^[18]^[19]^[20] A brief account of the use of these techniques is presented here, as well as an introduction of the fundamental solution and point-source concepts.

Fundamental Solution of the Diffusion Equation

The fundamental solution, γ_f(M, M′, t, τ), of the diffusion equation for fluid flow in porous media satisfies the following differential equation:

....................(3.144)

where δ(M, M′, t, τ) is a generalized (symbolic) function^[15] called the Dirac delta function and is defined on the basis of its following properties:

....................(3.145)

and

....................(3.146)

The delta function is symmetric in M and M′ and also in t and τ. In this formulation, the delta function represents the symbolic density of a unit-strength, concentrated source located at M′ and acting at time τ. In physical terms, this source corresponds to an infinitesimally small well (located at point M′) at which a finite amount of fluid is withdrawn (or injected) instantaneously (at time τ). Therefore, the solution of Eq. 3.144 (the fundamental solution) is also known as the instantaneous point-source solution. Formally, the point-source solution corresponds to the pressure drop, Δp = p_i − p, at a point M and time t in an infinite porous medium (reservoir) because of a point source of unit strength located at point M′ and acting at τ <t.

The Source-Function Solutions of the Diffusion Equation

The point-source solution was first introduced by Lord Kelvin^[16] for the solution of heat conduction problems and was extensively discussed by Carslaw and Jaeger.^[14] The point-source solution is usually obtained by finding the limiting form of the pressure drop resulting from a spherical source as the source volume vanishes. To demonstrate the derivation of the instantaneous point-source solution, consider the transient flow of a slightly compressible fluid of constant compressibility and viscosity toward a spherical source of radius r = a in an infinite, homogeneous, and isotropic porous medium. Because of the spherical symmetry of the physical problem, we can conveniently express the governing equation of fluid flow in porous media in spherical coordinates as

....................(3.147)

Assume that the initial pressure drop satisfies

....................(3.148)

and we have the condition that

....................(3.149)

On substitution of u = rΔp, Eqs. 3.147 through 3.149 become, respectively,

....................(3.150)

....................(3.151)

and

....................(3.152)

The solution of the problem described by Eqs. 3.150 through 3.152 is given by^[14]

....................(3.153)

If we expand the exponential terms in the integrand in Eq. 3.153 in powers of r′ and neglect the terms with powers higher than four, we obtain

....................(3.154)

In Eq. 3.154, 4πα³/3=V where V is the volume of the spherical source. If denotes the volume of the liquid released as a result of the change in the volume of the source, ΔV, which is caused by the change in pressure, Δp_i, then . With the definition of compressibility, c = -(1 / V)(ΔV / Δp_i), we obtain . Then, we can show that

....................(3.155)

Substituting Eq. 3.155 into Eq. 3.154, we obtain

....................(3.156)

If we let the radius of the spherical source, a, tend to zero while remains constant, Eq. 3.156 yields the point-source solution in spherical coordinates given by

....................(3.157)

This solution may be interpreted as the pressure drop at a distance r because of a volume of fluid, , instantaneously withdrawn at r = 0 and t = 0. Consistent with this interpretation, is the strength of the source, which is the pressure drop in a unit volume of the porous medium caused by the instantaneous withdrawal of a volume of fluid, (see Eq. 3.155).

Instantaneous Point Source in an Infinite Reservoir. Nisle^[21] presented a more general solution for an instantaneous point source of strength acting at t = τ in an infinite, homogeneous, but anisotropic reservoir as

....................(3.158)

In Eq. 3.158, M and M′ indicate the locations of the observation point and the source, respectively. For a 3D Cartesian coordinate system, with η_x, η_y, and η_z representing the diffusivity constants (defined in Eq. 3.23) in the x, y, and z directions, respectively.

Continuous Point Source in an Infinite Reservoir. If the fluid withdrawal is at a continuous rate, , from time 0 to t, then the pressure drop as a result of a continuous point source in an infinite reservoir is obtained by distributing the point sources of strength over a time interval 0 ≤ τ ≤ t. This is given by

....................(3.159)

where S(M, M′, t−τ) corresponds to a unit-strength , instantaneous point source in an infinite reservoir; that is,

....................(3.160)

Instantaneous and Continuous Line, Surface, and Volumetric Sources in an Infinite Reservoir. Similarly, the distribution of instantaneous point sources of strength over a line, surface, or volume, Γ_w, in an infinite reservoir leads to the following solution corresponding to the pressure drop because of production from a line, surface, or volumetric source, respectively.

....................(3.161)

In Eq. 3.161, M_w indicates a point on the source (Γ_w) and is the instantaneous withdrawal volume of fluids per unit length, area, or volume of the source, depending on the source geometry. For example, the pressure drop as a result of an infinite line source at x′, y′ and -∞≤ z′ ≤ ∞ may be obtained as follows:

....................(3.162)

If we assume that the flux is uniform along the line source and the source strength is unity , then we can write the instantaneous, infinite line-source solution in an infinite reservoir as

....................(3.163)

As another example, if we consider an instantaneous, infinite plane source at x = x′, -∞ ≤ y′ ≤ ∞, and -∞ ≤ z′ ≤ ∞ in an infinite reservoir, we can write

....................(3.164)

which also leads to the following uniform-flux, unit-strength, instantaneous, infinite plane-source solution in an infinite reservoir:

....................(3.165)

If the fluid withdrawal is at a continuous rate from time 0 to t, then the continuous line-, surface-, or volumetric-source solution for an infinite reservoir is given by

....................(3.166)

Source Functions for Bounded Reservoirs. The source solutions discussed previously can be extended to bounded reservoirs. The method of images provides a convenient means of generating the bounded-reservoir solutions with the use of the infinite reservoir solutions, especially when the reservoir boundaries consist of impermeable and constant-pressure planes. The method of images is an application of the principle of superposition, which states that if f₁ and f₂ are two linearly independent solutions of a linear PDE and c₁ and c₂ are two arbitrary constants, then f₃ = c₁f₁ + c₂f₂ is also a solution of the PDE. Examples of source functions in bounded reservoirs are presented here.

Instantaneous Point Source Near a Single Linear Boundary. To generate the effect of an impermeable planar boundary at a distance d from a unit-strength, instantaneous point source in an infinite reservoir (see Fig. 3.7), we can apply the method of images to the instantaneous point-source solution given in Eq. 3.157 as

....................(3.167)

It can be shown from Eq. 3.167 that (∂S/∂x)_x=d = 0; that is, the bisector of the distance between the two sources is a no-flow boundary. Similarly, to generate the effect of a constant-pressure boundary, we use the method of images and the unit-strength, instantaneous point-source solution (Eq. 3.160) as follows:

....................(3.168)

Instantaneous Point Source in an Infinite-Slab Reservoir. Using the method of images and considering the geometry shown in Col. A of Fig. 3.8, we can generate the solution for a unit-strength, instantaneous point source in an infinite-slab reservoir with impermeable boundaries at z = 0 and h. The result is given by

....................(3.169)

which, with Poisson’s summation formula given by^[14]

....................(3.170)

may be transformed into

....................(3.171)

Fig. 3.7 – Application of the method of images to generate the effect of a linear boundary.
Fig. 3.8 – Application of the method of images to generate the solutions for infinite-slab reservoirs.

Following similar lines, if the slab boundaries at z = 0 and h are at a constant pressure equal to p_i, we obtain

....................(3.172)

Similarly, for the case in which the slab boundary at z = 0 is impermeable while the boundary at z = h is at a constant pressure equal to p_i, the following solution may be derived:

....................(3.173)

Instantaneous Point Source in a Closed Parallelepiped. The ideas used previously for slab reservoirs may be used to develop solutions for reservoirs bounded by planes in all three directions. For example, if the reservoir is bounded in all three directions (i.e., 0 ≤ x ≤ x_e, 0 ≤ y ≤ y_e, and 0 ≤ z ≤ h) and the bounding planes are impermeable, then we can use Eq. 3.157 and the method of images to write

....................(3.174)

which, with Poisson’s summation formula (Eq. 3.170), may be recast into the following form:

....................(3.175)

Instantaneous Infinite-Plane Source in an Infinite-Slab Reservoir With Impermeable Boundaries. The instantaneous point-source solutions of Eqs. 3.171 through 3.173 may be extended to different source geometries with Eq. 3.161. For example, the solution for an instantaneous infinite-plane source at z = z′ in an infinite-slab reservoir with impermeable boundaries is obtained by substituting Eq. 3.171 for S in Eq. 3.161. This yields

....................(3.176)

Assuming a unit-strength, uniform-flux source , we obtain the following instantaneous infinite-plane source solution in an infinite-slab reservoir with impermeable boundaries:

....................(3.177)

Instantaneous Infinite-Slab Source in an Infinite-Slab Reservoir With Impermeable Boundaries. Following similar lines, we can obtain the solution for an instantaneous, infinite-slab source of thickness, h_p, located at z = z_w (z_w is the z-coordinate of the midpoint of the slab source) in an infinite-slab reservoir with impermeable boundaries.

....................(3.178)

If we assume a uniform-flux slab source , then Eq. 3.178 yields

....................(3.179)

Uniform-Flux, Continuous, Infinite-Slab Source in an Infinite-Slab Reservoir With Impermeable Boundaries. Solutions for continuous plane and slab sources can be obtained as indicated by Eq. 3.159 (or Eq. 3.166). For example, the solution for a uniform-flux, continuous, infinite-slab source in an infinite-slab reservoir with impermeable top and bottom boundaries may be obtained by substituting the right side of Eq. 3.179 for S in Eq. 3.159 and is given by

....................(3.180)

The same solution could have been obtained by substituting the unit-strength instantaneous point-source solution given by Eq. 3.171 for S in Eq. 3.166.

Example 3.5

Consider transient flow toward a partially penetrating vertical well of penetration length, h_w, in an infinite, homogeneous, slab reservoir of uniform thickness, h, and initial pressure, p_i, with impermeable top and bottom boundaries.

Solution. Fig. 3.9 shows the geometry of the well and reservoir system of interest. The solution for this problem can be obtained by assuming that the well may be represented by a vertical line source. Then, starting with Eq. 3.166 and substituting the unit-strength, instantaneous point-source solution in an infinite-slab reservoir with impermeable boundaries [Eq. 3.171 with ] for S, we obtain

....................(3.181)

If we assume that the strength of the source is uniformly distributed along its length (this physically corresponds to a uniform-flux distribution) and the production rate is constant over time [i.e., , where q is the constant production rate of the well], then Eq. 3.181 yields

....................(3.182)

Fig. 3.9 – Geometry of the well/reservoir system for a partially penetrating vertical well in an infinite-slab reservoir with impermeable boundaries for Example 3.5.

The Use of Green’s Functions and Source Functions in Solving Unsteady-Flow Problems

As discussed in Sec. 3.4.2, 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^[13] 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.^[22] In this section, we discuss the use of Green’s functions and source functions in solving unsteady-flow problems in reservoirs.

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.^[13]^[14] If we let G(M, M′, t − τ) denote the Green’s function, then it should satisfy the diffusion equation; that is,

....................(3.183)

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

....................(3.184)

Green’s function also has the following properties: ^[13]^[14]

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.
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.185)

where D indicates the domain of the porous medium, and φ(M) is any continuous function.

Because G corresponds to the pressure because of an instantaneous point source of unit strength, it satisfies

....................(3.186)

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:

....................(3.187)

and

....................(3.188)

Then, we can write

....................(3.189)

or

....................(3.190)

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

....................(3.191)

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

....................(3.192)

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

....................(3.193)

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

In Eq. 3.193, 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. 3.193 may be written as

....................(3.194)

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, G(Γ_e) = 0. Thus, Eq. 3.194 becomes

....................(3.195)

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

....................(3.196)

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

....................(3.197)

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

....................(3.198)

The substitution of Eqs. 3.197 and 3.198 into Eq. 3.196 yields

....................(3.199)

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

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

....................(3.200)

where

....................(3.201)

The integration in the right side of Eq. 3.201 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. In Sec. 3.4.2, we discussed 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.^[22]

Newman’s product method^[22] may be stated for transient-flow problems in porous media as follows: ^[13] 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

....................(3.202)

as follows

....................(3.203)

As expected, this solution is the same as Eq. 3.163, which was obtained by integration of an instantaneous point source in an infinite reservoir. For a radially isotropic reservoir (η_x = η_y = η_z), Eq. 3.203 yields

....................(3.204)

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

....................(3.205)

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

....................(3.206)

Instantaneous plane sources in slab reservoirs can be generated with the plane sources in infinite reservoirs and the method of images as discussed in Sec. 3.4.2. Similarly, the instantaneous slab sources can be obtained by integrating plane sources over the thickness of the slab source (see Sec. 3.4.2). Table 3.2, compiled from the work of Gringarten and Ramey,^[13] presents the basic instantaneous source functions in infinite reservoirs, and Table 3.3 shows the corresponding geometries of the source-reservoir systems. The basic instantaneous source functions given in Table 3.3 may be used to construct the source functions that represent the appropriate well geometry by Newman’s product method.

Fig. 3.10 – Geometry of a line source in 3D Cartesian and radial coordinates.
Table 3.2
Table 3.3

Gringarten and Ramey^[13] 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.4 presents the short- and long-time approximating forms for instantaneous linear sources and their time ranges. Examples 3.6 and 3.7 present some applications of the product-solution method and the derivation of the approximate solutions for pressure distributions.

Table 3.4

Example 3.6

Consider transient flow toward a partially penetrating vertical fracture of vertical penetration h_f and horizontal penetration 2x_f in an infinite, homogeneous, slab reservoir of uniform thickness, h, and initial pressure, p_i, with impermeable top and bottom boundaries.

Solution. Fig. 3.11 shows the geometry of the well reservoir system of interest. Approximate the fracture by a vertical plane of height h_f and length 2x_f. 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 3.2 and 3.3), an infinite-slab source of thickness 2x_f at x = x′ in an infinite reservoir (Source IV), and an infinite-slab source of thickness h_p = h_f at z = z_w in an infinite-slab reservoir of thickness h (Source VIII). Then, by Newman’s product method, the appropriate source function is given by

....................(3.207)

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

....................(3.208)

If the fracture penetrates the entire thickness of the reservoir (i.e., h_f = h) as shown in Fig. 3.12, then Eq. 3.208 yields

....................(3.209)

The fully penetrating fracture solution given in Eq. 3.209 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 3.2 and 3.3) and an infinite-slab source of thickness 2x_f at x = x′ in an infinite reservoir (Source IV). This source function then could be used in Eq. 3.200.

$Fig. 3.11 – Geometry of the well/reservoir system for a partially penetrating vertical fracture in an infinite-slab reservoir with impermeable boundaries for Example 3.6.$

Fig. 3.11 – Geometry of the well/reservoir system for a partially penetrating vertical fracture in an infinite-slab reservoir with impermeable boundaries for Example 3.6.
$Fig. 3.12 – A fully penetrating vertical fracture in an infinite-slab reservoir with impermeable boundaries.$

Fig. 3.12 – A fully penetrating vertical fracture in an infinite-slab reservoir with impermeable boundaries.

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

....................(3.210)

where

....................(3.211)

and

....................(3.212)

$Fig. 3.13 – Transient pressure responses of a uniform-flux fracture computed from Eq. 3.209.$

Fig. 3.13 – Transient pressure responses of a uniform-flux fracture computed from Eq. 3.209.

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

....................(3.213)

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

....................(3.214)

Example 3.7

Consider transient flow toward a uniform-flux horizontal well of length L_h located at (x′, y′, z_w) in a closed, homogeneous, rectangular parallelepiped of dimensions 0 ≤ x ≤ x_e, 0 ≤ y ≤ y_e, 0 ≤ z ≤ h and of initial pressure, p_i.

Solution. Fig. 3.14 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 L_h, 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 3.2 and 3.3), an infinite-slab source of thickness L_h at x = x′ in an infinite-slab reservoir of thickness x_e with impermeable boundaries (Source VIII), and an infinite-plane source at z = z_w 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

....................(3.215)

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

....................(3.216)

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

....................(3.217)

where Ei(−x) is the exponential integral function defined by Eq. 3.90. Eq. 3.217 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.

Fig. 3.14 – A horizontal well in a closed rectangular parallelepiped (Example 3.7).
Table 3.5

It is also possible to obtain another approximation for Eq. 3.216 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.4, we obtain

....................(3.218)

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 3.2 and 3.3), an infinite-slab source of thickness L_h at x = x′ in an infinite reservoir (Source IV), and an infinite-plane source at z = z_w in an infinite-slab reservoir of thickness h with impermeable boundaries (Source V).

The Use of Source Functions in the Laplace Domain To Solve Unsteady-Flow Problems

There are many advantages of developing transient-flow solutions in the Laplace transform domain. For example, in the Laplace transform domain, Duhamel’s theorem^[23] provides a convenient means of developing transient-flow solutions for variable-rate production problems using the solutions for the corresponding constant-rate production problem. Duhamel’s theorem states that if Δp and Δp_c denote the pressure drawdown corresponding to the variable production rate, q(t), and the constant production rate, q_c, respectively, then

....................(3.219)

Applying the Laplace transform converts the convolution integral in Eq. 3.219 to an algebraic expression, and Duhamel’s theorem is given in the Laplace transform domain as

....................(3.220)

The simplicity of the expression given in Eq. 3.220 explains our interest in obtaining transient-flow solutions in the Laplace transform domain.

Another example to explain the convenience of the Laplace domain solutions is for the naturally fractured reservoirs. Common transient-flow models of naturally fractured reservoirs lead to the following differential equation in radial coordinates in the Laplace transform domain: ^[2]

....................(3.221)

where the subscript f indicates the fracture property, and t_D and r_D are the dimensionless time and distance (as defined in Eqs. 3.230 and 3.234).

The naturally fractured reservoir function, f (s), is a function of matrix and fracture properties and depends on the model chosen to represent the naturally fractured reservoir.^[2] The corresponding differential equation for a homogeneous reservoir is obtained by setting f (s) = 1 and is given by

....................(3.222)

The general solutions for Eqs. 3.221 and 3.222 are given, respectively, by

....................(3.223)

and

....................(3.224)

To obtain a solution for constant-rate production from an infinite reservoir, for example, the following boundary conditions are imposed:

....................(3.225)

and

....................(3.226)

Then, it may be shown that

....................(3.227)

where the right side of Eq. 3.227 indicates the substitution of sf (s) for s in sΔp(s). This discussion demonstrates that it is possible to derive transient-flow solutions for naturally fractured reservoirs by following the same lines as those for the homogeneous reservoirs. Furthermore, if the solution for the corresponding homogeneous reservoir system is known in the Laplace transform domain, then the solution for the naturally fractured reservoir problem may be directly obtained from Eq. 3.227.

Obtaining the Laplace transforms of the Green’s and source function solutions developed in the time domain with the methods explained in Secs. 3.4.2 and 3.4.3 usually poses a difficult problem. The problems arise mainly because of the use of the product method solution. For a specific class of functions, Chen et al.^[24] presented a technique that may be used to apply the Laplace transform to the product-solution technique. For a more general procedure to develop source function solutions in the Laplace transform domain, however, the product solution technique should be avoided.^[20]

Ozkan and Raghavan^[18]^[19] have shown that it is more convenient to develop source-function solutions in the Laplace transform domain if the point-source solution is used as a building block. Then, other source geometries may be obtained by the superposition (integration) of the point sources along the length, surface, or volume of the source.

Point-Source Solution in the Laplace Domain. Consider the flow of a slightly compressible fluid in an infinite, naturally fractured reservoir. We can use the double-porosity model suggested by Barenblatt et al.^[25] and Warren and Root^[26] to develop the governing flow equations for naturally fractured reservoirs. The results, however, will be applicable to the model suggested by Kazemi^[27] and de Swaan-O^[28] with a simple modification.

Flow around a point source in an infinite porous medium may be expressed conveniently in spherical coordinates. The differential equations governing flow in a naturally fractured reservoir are given in spherical coordinates by

....................(3.228)

and

....................(3.229)

In Eqs. 3.228 and 3.229, subscripts f and m indicate the property of the fracture and matrix systems, respectively. Initial pressure, p_i, is assumed to be uniform in the entire system; that is, p_fi = p_mi = p_i. The dimensionless time, t_D, is defined by

....................(3.230)

where ℓ is a characteristic length in the system, and

....................(3.231)

The definitions of the other variables used in Eqs. 3.228 and 3.229 are

....................(3.232)

....................(3.233)

and

....................(3.234)

where

....................(3.235)

The initial and outer-boundary conditions are given, respectively, by

....................(3.236)

and

....................(3.237)

The inner-boundary condition corresponding to the instantaneous withdrawal of an amount of fluid, , at t = 0 from a point source is obtained by considering the mass balance on a small sphere. If we require that at any time t = T > 0, the sum of the flux through the surface of a small sphere around the source location must equal the volume of the fluid, , instantaneously withdrawn from the sphere at t = 0, we can write^[29]

....................(3.238)

Although the withdrawal of fluids from the sphere is instantaneous, the resulting flow in the porous medium, and the flux across the surface of the sphere, is continuous. Therefore, if q represents the total flux across the surface of the small sphere during the time interval 0 ≤ t ≤ T, then the mass balance requires that the cumulative production (flux across the surface of the small sphere) at time T be equal to the instantaneous withdrawal volume of fluid from the sphere. That is,

....................(3.239)

For the condition expressed in Eq. 3.239 to hold for every T ≥ 0, we must have

....................(3.240)

where δ(t) is the Dirac delta function satisfying the properties expressed by Eqs. 3.185 and 3.186.

Using the results given by Eqs. 3.239 and 3.240 in Eq. 3.238, we obtain

....................(3.241)

The Laplace transform of Eqs. 3.228, 3.229, 3.237, and 3.241 yields

....................(3.242)

where

....................(3.243)

....................(3.244)

and

....................(3.245)

In deriving these results, we have used the initial condition given by Eq. 3.236 and noted that

....................(3.246)

In Eq. 3.245, the term represents the strength of the source for the naturally fractured porous medium.

The solution of Eqs. 3.242, 3.244, and 3.245 yields the following solution for the pressure distribution in the reservoir, except at the source location (the origin), because of an instantaneous point source of strength acting at t = 0:

....................(3.247)

If the source is located at x′_D, y′_D, z′_D, then, by translation, we can write

....................(3.248)

where

....................(3.249)

and

....................(3.250)

The instantaneous point-source solution for the model suggested by Barenblatt et al.^[25] and Warren and Root^[26] can also be used for the model suggested by Kazemi^[27] and de Swaan-O,^[28] provided that the appropriate f(s) function is invoked. To obtain the solution for a homogeneous reservoir, f(s) should be set to unity, V_f = 1, and V_m = 0.

If we consider continuous withdrawal of fluids from the point source, then, by the principle of superposition, we should have

....................(3.251)

The Laplace transform of Eq. 3.251 yields the following continuous point-source solution in an infinite reservoir:

....................(3.252)

where we have substituted Eq. 3.249 for S, dropped the subscript f, and defined

....................(3.253)

Line-, Surface-, and Volumetric-Source Solution in the Laplace Domain. The point-source solution in the Laplace domain may be used to obtain the source solutions for different source geometries. If we define

....................(3.254)

where Δp_p represents the appropriate point-source solution, then, by the application of the superposition principle, the solution for the withdrawal of fluids from a line, surface, or volume, Γ_w, is given by

....................(3.255)

If we assume a uniform-flux distribution in time and over the length, surface, or volume of the source, then

....................(3.256)

The following presentation of the source function approach in the Laplace domain assumes that the flux distribution is uniform, and . Also, the constant production rate from the length, area, or the volume of the source, Γ_w, is denoted by q so that .

Only sources in infinite reservoirs have been considered so far. These solutions may be easily extended to bounded reservoirs. The following sections present some useful solutions for transient-flow problems in bounded porous media. The first group of solutions is for laterally infinite reservoirs bounded by parallel planes in the vertical direction (infinite-slab reservoirs). The second and third groups comprise the solutions for cylindrical and rectangular reservoirs, respectively.

Solutions for Infinite-Slab Reservoirs. In this section, we consider one of the most common reservoir geometries used in pressure-transient analysis of wells in porous media. It is assumed that the lateral boundaries of the reservoir are far enough not to influence the pressure response during the time period of interest. The top and bottom boundaries of the reservoir at z = 0 and z = h are parallel planes and may be of impermeable, constant pressure, or mixed type. Table 3.6 presents the solutions for the most common well geometries (point-source, vertical, fractured, and horizontal wells) in infinite-slab reservoirs. Next, we briefly discuss the derivation of these solutions.

Table 3.6
Table 3.6 (continued)

Consider a point source in an infinite-slab reservoir with impermeable boundaries at the bottom, z = 0, and the top, z = h. To obtain the point-source solution for this case, we use the point-source solution in an infinite reservoir given by Eq. 3.252 with the method of images. The result is given by

....................(3.257)

where

....................(3.258)

....................(3.259)

....................(3.260)

and

....................(3.261)

The solution given in Eq. 3.257 is not very convenient for computational purposes. To obtain a computationally convenient form of the solution, we use the summation formula given by^[17]^[29]

....................(3.262)

and recast Eq. 3.257 as

....................(3.263)

The point-source solutions for infinite-slab reservoirs with constant pressure and mixed boundaries at the top and bottom are obtained in a similar manner^[17] and are given in Table 3.6. The point-source solutions can be used with Eqs. 3.254 and 3.256 to generate the solutions for the other well geometries given in Table 3.6. For example, to generate the solution for a partially penetrating vertical line-source well of length hw in an infinite-slab reservoir with impermeable slab boundaries, we integrate the right side of Eq. 3.263 from z_w − h_w / 2 to z_w + h_w / 2 with respect to z′, where z_w is the vertical coordinate of the midpoint of the open interval. If h_w = h (i.e., the well penetrates the entire thickness of the slab reservoir), then this procedure yields the solution for a fully penetrating vertical line-source well. The solution for a partially penetrating fracture of height h_f and half-length x_f is obtained if the point-source solution is integrated once with respect to z′ from z_w − h_f / 2 to z_w + h_f / 2 and then with respect to x′ from x_w – x_f to x_w + x_f, where x_w and z_w are the coordinates of the midpoint of the fracture. Similarly, the solution for a horizontal-line source well of length L_h is obtained by integrating the point-source solution with respect to x′ from x_w − L_h / 2 to x_w + L_h / 2, where x_w is the x-coordinate of the midpoint of the horizontal well.

Solutions for Cylindrical Reservoir. Solutions for cylindrical reservoirs may also be obtained by starting from the point-source solution in the Laplace transform domain. The Laplace domain solution for a point source located at r′_D, θ′, z′_D should satisfy the following diffusion equation in cylindrical coordinates.^[19]

....................(3.264)

where

....................(3.265)

The point-source solution is also required to satisfy the following flux condition at the source location (r_D →0+, θ = θ′, z_D = z′_D):

....................(3.266)

Assuming that the reservoir is bounded by a cylindrical surface at r_D = r_eD and by the parallel planes at z_D = 0 and h_D, we should also impose the appropriate physical conditions at these boundaries.

We seek a point-source solution for a cylindrical reservoir in the following form:

....................(3.267)

In Eq. 3.267, is a solution of Eq. 3.264 that satisfies Eq. 3.266 and the boundary conditions at z_D = 0 and h_D. may be chosen as one of the point-source solutions in an infinite-slab reservoir given in Table 3.6, depending on the conditions imposed at the boundaries at z_D = 0 and h_D. If is chosen such that it satisfies the boundary conditions at z_D = 0 and h_D, its contribution to the flux vanishes at the source location, and + satisfies the appropriate boundary condition at r_D = r_eD, then Eq. 3.267 should yield the point-source solution for a cylindrical reservoir with appropriate boundary conditions.

Consider the example of a closed cylindrical reservoir in which the boundary conditions are given by

....................(3.268)

and

....................(3.269)

According to the boundary condition given by Eq. 3.268, we should choose as the point-source solution given in Table 3.6 (or by Eq. 3.263). Then, with the addition theorem for the Bessel function K₀(aR_D) given by^[14]

....................(3.270)

where

....................(3.271)

we can write

....................(3.272)

for r_D < r′_D. If r_D > r′_D, we interchange r_D and r′_D in Eq. 3.272. If we choose in Eq. 3.267 as

....................(3.273)

where a_k and b_k are constants, then satisfies the boundary condition given by Eq. 3.268, and the contribution of to the flux at the source location vanishes. If we also choose the constants a_k and b_k in Eq. 3.273 as

....................(3.274)

and

....................(3.275)

then satisfies the impermeable boundary condition at r_D = r_eD given by Eq. 3.269. Thus, the point-source solution for a closed cylindrical reservoir is given by

....................(3.276)

This solution procedure may be extended to the cases in which the boundaries are at constant pressure or of mixed type.^[19] Table 3.7 presents the point-source solutions for cylindrical reservoirs for all possible combinations of boundary conditions. Solutions for other source geometries in cylindrical reservoirs may be obtained by using the point-source solutions in Table 3.7 in Eq. 3.255 (or Eq. 3.256).

Table 3.7

Example 3.8

Consider a partially penetrating, uniform-flux fracture of height h_f and half-length x_f in an isotropic and closed cylindrical reservoir. The center of the fracture is at r′ = 0, θ′ =0, z′ = z_w, and the fracture tips extend from (r′ = x_f, θ = α + π) to (r′ = x_f, θ = α).

Solution. Fig. 3.15 shows the geometry of the fracture/reservoir system considered in this example. The solution for this problem is obtained by first generating a partially penetrating line source and then using this line-source solution to generate the plane source. The solution for a partially penetrating line source at r′_D, θ′, z_w is obtained by integrating the corresponding point-source solution given in Table 3.7 with respect to z′ from z_w – h_f / 2 to z_w + h_f / 2 and is given by

....................(3.277)

$Fig. 3.15 – Geometry of a partially penetrating fracture in a closed cylindrical reservoir (Example 3.8).$

Fig. 3.15 – Geometry of a partially penetrating fracture in a closed cylindrical reservoir (Example 3.8).

To generate the solution for a partially penetrating plane source that represents the fracture, the partially penetrating line-source solution given in Eq. 3.277 is integrated with respect to r′ from 0 to x_f with θ′ = α + π in the third quadrant and with θ′ = α in the first quadrant. This procedure yields

....................(3.278)

It is possible to obtain an alternate representation of the solution given in Eq. 3.278. With the addition theorem of the Bessel function K₀(x) given by Eq. 3.270, the solution in Eq. 3.277 may be written as

....................(3.279)

where

....................(3.280)

and

....................(3.281)

The integration of the partially penetrating vertical well solution given in Eq. 3.279 with respect to r′ from 0 to x_f (with θ′ = α + π in the third quadrant and with θ′ = α in the first quadrant) yields the following alternative form of the partially penetrating fracture solution:

....................(3.282)

where

....................(3.283)

Example 3.9 Consider a uniform-flux, horizontal line-source well of length L_h in an isotropic and closed cylindrical reservoir. The well extends from (r′ = L_h/2, θ = α + π) to (r′ = L_h/2, θ = α), and the center of the well is at r′ = 0, θ′ = 0, z′ = z_w.

Solution. The solution for a horizontal line-source well in a closed cylindrical reservoir is obtained by integrating the corresponding point-source solution in Table 3.7 with respect to r′ from 0 to L_h / 2 with θ′ = α + π in the third quadrant and with θ′ = α in the first quadrant. The final form of the solution is given by

....................(3.284)

Solutions for Rectangular Parallelepiped Reservoir. Solutions for rectangular parallelepiped reservoirs may also be obtained by starting from the point-source solution in the Laplace transform domain in an infinite reservoir and using the method of images to generate the effects of the planar boundaries. Although the formal procedure to obtain the solution is fairly easy, the use of the method of images in three directions (x, y, z) yields triple infinite Fourier series, which may pose computational inconveniences. As an example, the solution for a continuous point source located at x′, y′, z′ in a rectangular porous medium occupying the region 0 < x < x_e, 0 < y < y_e, and 0 < z < h is obtained by applying the method of images to the point-source solution given by Eq. 3.252: ^[19]^[29]

....................(3.285)

where

....................(3.286)

and

....................(3.287)

....................(3.288)

....................(3.289)

Ozkan^[29] shows that the triple infinite sums in Eq. 3.285 may be reduced to double infinite sums with

....................(3.290)

where

....................(3.291)

The resulting continuous point-source solution for a closed rectangular reservoir is given by

....................(3.292)

where

....................(3.293)

....................(3.294)

....................(3.295)

....................(3.296)

and

....................(3.297)

Following a procedure similar to the one explained here, it is possible to derive the point-source solutions in rectangular parallelepiped reservoirs for different combinations of boundary conditions.^[17]^[29] Table 3.8 gives these solutions, which may be used to derive the solutions for the other source geometries with Eq. 3.255 (or Eq. 3.256). Examples 3.10 and 3.11 demonstrate the derivation of the solutions for the other source geometries in rectangular reservoirs.

Table 3.8
Table 3.8
Table 3.8
Table 3.8

Example 3.10

Consider a fully penetrating vertical fracture of half-length x_f located at x′ = x_w and y′ = y_w in a closed rectangular reservoir.

Solution. Assuming uniform-flux distribution along the fracture surface, the solution for this problem is obtained by integrating the corresponding point-source solution in Table 3.8, first with respect to z′ from 0 to h and then with respect to x′ from x_w – x_f to x_w + x_f. The result is

....................(3.298)

where , , and ε_k are given respectively by Eqs. 3.293, 3.294, and 3.296.

Example 3.11

Consider a horizontal well of length L_h in the x-direction located at x′ = x_w, y′ = y_w, and z′ = z_w in a closed rectangular reservoir.

Solution. The solution for a horizontal line-source well is obtained by integrating the corresponding point-source solution in Table 3.8, with respect to x′ from x_w–L_h /2 to x_w+L_h /2, and is given by

....................(3.299)

where

....................(3.300)

and

....................(3.301)

In Eq. 3.301, , , ε_n, ε_k, and ε_{k, n} are given by Eqs. 3.293 through 3.297.

Conversion From 3D to 2D Anisotropy. The solutions previously presented assume that the reservoir is anisotropic in all three principal directions, x, y, and z with k_x, k_y, and k_z denoting the corresponding permeabilities. In these solutions, an equivalent isotropic permeability, k, has been defined by

....................(3.302)

For some applications, it may be more convenient to define an equivalent horizontal permeability by

....................(3.303)

and replace k in the solutions given in this section (Sec. 3.4.4) by k_h. Note that k takes place in the definition of the dimensionless time t_D (Eq. 3.230). Then, if we define a dimensionless time based on k_h, the relation between and t_D is given by

....................(3.304)

Because in the solutions given in this section the Laplace transformation is with respect to t_D, conversion from 3D to 2D anisotropy requires the use of the following property of the Laplace transforms:

....................(3.305)

As an example, consider the solution for a horizontal well in an infinite-slab reservoir. Assuming that the midpoint of the well is the origin (x_wD = 0, y_wD = 0) and choosing the half-length of the horizontal well as the characteristic length (i.e., ℓ = L_h / 2), the horizontal-well solution given in Table 3.6 may be written as

....................(3.306)

In Eq. 3.306, s is the Laplace transform variable with respect to dimensionless time, t_D, based on k and

....................(3.307)

....................(3.308)

....................(3.309)

and

....................(3.310)

If we define the following variables based on k_h,

....................(3.311)

....................(3.312)

....................(3.313)

and also note that

....................(3.314)

then, we may rearrange Eq. 3.306 in terms of the dimensionless variables based on k_h as

....................(3.315)

where

....................(3.316)

and

....................(3.317)

If we compare Eqs. 3.306 and 3.315, we can show that

....................(3.318)

where we have used the relation given by Eq. 3.305. If we now define as the Laplace transform variable with respect to , we may write

....................(3.319)

With the relation given by Eq. 3.319 and Eq. 3.306, we obtain the following horizontal-well solution in terms of dimensionless variables based on k_h:

....................(3.320)

Computational Considerations and Applications. The numerical evaluation of the solutions given previously may be sometimes difficult, inefficient, or even impossible. Alternative computational forms of some of these solutions have been presented in a few sources.^[18]^[19]^[29] Here, we present a summary of the alternative formulas to be used in the computation of the source functions in the Laplace transform domain. Some of these formulas are for computations at early or late times and may be useful to derive asymptotic approximations of the solutions during the corresponding time periods.

As Sec. 3.3.1 notes, the short- and long-time approximations of the solutions correspond to the limiting forms of the solution in the Laplace transform domain as s→∞ and s→0, respectively. In the solutions given in this section, we have defined u = sf(s). From elementary considerations, it is possible to show that the definition of f(s) given in Eq. 3.243 yields the following limiting forms:

....................(3.321)

and

....................(3.322)

These limiting forms are used in the derivation of the short- and long-time asymptotic approximations. In the following expressions, homogeneous reservoir solutions are obtained by substituting ω = 1.

The Integral.

....................(3.323)

This integral arises in the computation of many practical transient-pressure solutions and may not be numerically evaluated, especially as y_D→0; however, the following alternate forms of the integral are numerically computable.^[19]

....................(3.324)

....................(3.325)

and

....................(3.326)

The integrals in Eqs. 3.324 through 3.326 may be evaluated with the standard numerical integration algorithms for y_D ≠ 0. For y_D = 0, the polynomial approximations given by Luke^[30] or the following power series expansion given by Abramowitz and Stegun^[7] may be used in the computation of the integrals in Eqs. 3.324 through 3.326:

....................(3.327)

For numerical computations and asymptotic evaluations, it may also be useful to note the following relations: ^[19]

....................(3.328)

and

....................(3.329)

It can be shown from Eqs. 3.328 and 3.329 that, for practical purposes, when z ≥ 20, the right sides of Eqs. 3.327 and 3.328 may be approximated by π/2 and π exp (−|c|)/2, respectively.^[19]^[27]

As a few sources^[18]^[19]^[29] show, it is possible to derive the following short- and long-time approximations (i.e., the limiting forms as s→∞ and s→0, respectively) for the integral given, respectively, by

....................(3.330)

where

....................(3.331)

and

....................(3.332)

where γ=0.5772… and

....................(3.333)

It is also useful to note the real inversions of Eqs. 3.330 and 3.332 given, respectively, by

....................(3.334)

and

....................(3.335)

The Series.

....................(3.336)

Two alternative expressions for the series S₁ may be convenient for the large and small values of u (i.e., for short and long times).^[29] When u is large,

....................(3.337)

and when u + a² << n²π²/h²_D,

....................(3.338)

The Series .

....................(3.339)

Alternative computational forms for the series S₂ are given next.^[29] When u is large,

....................(3.340)

and when u + a² << n²π²/h²_D,

....................(3.341)

The Series .

....................(3.342)

The following alternative forms for the series may be convenient for the large and small values of u (i.e., for short and long times).^[29] When u is large,

....................(3.343)

and when u + a² << (2n − 1)² π²/(4h²_D),

....................(3.344)

The Series .

....................(3.345)

where

....................(3.346)

The series may be written in the following forms with the use of Eqs. 3.324 through 3.326.

....................(3.347)

....................(3.348)

and

....................(3.349)

The computation of the series in Eqs. 3.347 and 3.348 should not pose numerical difficulties; however, the series in Eq. 3.349 converges slowly. With the relation given in Eq. 3.328, we may write Eq. 3.349 as^[29]

....................(3.350)

where

....................(3.351)

Before discussing the computation of the series given in Eq. 3.351, we first discuss the derivation of the asymptotic approximations for the series . When s is large (small times), may be approximated by^[29]

....................(3.352)

where β is given by Eq. 3.331. If s is sufficiently large, then Eq. 3.352 may be further approximated by

....................(3.353)

The inverse Laplace transform of Eq. 3.353 yields

....................(3.354)

For small s (large times), depending on the value of x_D, may be approximated by one of the following equations: ^[29]

....................(3.355)

....................(3.356)

....................(3.357)

where is given by Eq. 3.364.

The Series .

....................(3.358)

where

....................(3.346)

With the relations given in Eqs. 3.337 and 3.338, the following alternative forms for the series may be obtained, respectively, for the large and small values of s (i.e., for short and long times).^[29] When u is large,

....................(3.359)

and when u << n²π²/h²_D,

....................(3.360)

It is also possible to derive asymptotic approximations for the series . When s is large (small times), may be approximated by^[29]

....................(3.361)

If s is sufficiently large, then Eq. 3.361 may be further approximated by

....................(3.362)

The inverse Laplace transform of Eq. 3.362 yields

....................(3.363)

For small s (large times), may be approximated by^[29]

....................(3.364)

The Ratio .

....................(3.365)

By elementary considerations, the ratio may be written as^[29]

....................(3.366)

The expression given in Eq. 3.366 provides computational advantages when s is small (time is large).

Example 3.12

Consider a fully penetrating, uniform-flux fracture of half-length x_f located at x′=0, y′=0 in an infinite-slab reservoir with closed top and bottom boundaries.

Solution.Table 3.6 gives the solution for this problem. For simplicity, assuming an isotropic reservoir, choosing the characteristic length as ℓ = x_f and noting that , the solution becomes

....................(3.367)

First consider the numerical evaluation of Eq. 3.367. We note from Eqs. 3.324 through 3.326 that Eq. 3.367 may be written in one of the following forms, depending on the value of x_D.

....................(3.368)

....................(3.369)

and

....................(3.370)

The numerical evaluation of the integrals in Eqs. 3.368 through 3.370 for y_D ≠ 0 should be straightforward with the use of the standard numerical integration algorithms. For y_D = 0, the polynomial approximations given by Luke^[30] or the power series expansion given by Eq. 3.327 should be useful.

The short- and long-time asymptotic approximations of the fracture solution are also obtained by applying the relations given by Eqs. 3.330 and 3.332, respectively, to the right side of Eq. 3.367. This procedure yields, for short times,

....................(3.371)

or, in real-time domain,

....................(3.372)

where β is given by Eq. 3.331 with a = -1 and b = +1. At long times, the following asymptotic approximation may be used:

....................(3.373)

or, in real-time domain,

....................(3.374)

where γ = 0.5772… and σ(x_D, y_D, -1, +1) is given by Eq. 3.333.

Example 3.13

Consider a horizontal well of length L_h located at x′ = 0, y′ = 0, and z′ = z_w in an infinite-slab reservoir with closed top and bottom boundaries.

Solution.Table 3.6 gives the horizontal-well solution for an infinite-slab reservoir with impermeable boundaries. Assuming an isotropic reservoir, choosing the characteristic length as ℓ = L_h / 2 and noting that , the solution may be written as

....................(3.375)

where is the fracture solution given by the right side of Eq. 3.367 and is given by

....................(3.376)

with

....................(3.346)

....................(3.377) br>
and

....................(3.378)

The computation of the first term in the right side of Eq. 3.375 is the same as the computation of the fracture solution given by Eq. 3.367 and has been discussed in Example 3.12. The computational form of the second term in the right side of Eq. 3.375 is given by Eqs. 3.347 through 3.350. Of particular interest is the case for −1 ≤ x_D ≤ +1. In this case, from Eqs. 3.350 and 3.351, we have

....................(3.379)

where

....................(3.380)

The computational considerations for the series have been discussed previously.

Next, we consider the short- and long-time approximations of the horizontal-well solution given by Eq. 3.375. To obtain a short-time approximation, we substitute the asymptotic expressions for and as s→∞ given, respectively, by Eqs. 3.371 and 3.353. This yields

....................(3.381)

where β is given by Eq. 3.331. The inverse Laplace transform of Eq. 3.381 is given by

....................(3.382)

To obtain the long-time approximation of Eq. 3.375, we substitute the asymptotic expressions for and as s→∞ given, respectively, by Eq. 3.374 and Eqs. 3.355 through 3.357. Of particular interest is the case for −1 ≤ x_D ≤ +1, where we have

....................(3.383)

where γ=0.5772… and σ(x_D, y_D, -1, +1) is given by Eq. 3.333. The inverse Laplace transform of Eq. 3.383 yields

....................(3.384)

Example 3.14 Consider a fully penetrating, uniform-flux fracture of half-length x_f in an isotropic and closed cylindrical reservoir. The center of the fracture is at r′ = 0, θ′ = 0 and the fracture tips extend from (r′ = x_f, θ = α + π) to (r′ = x_f, θ = α).

Solution. The solution for this problem has been obtained in Eq. 3.278 in Example 3.8 with h_w = h. Choosing the characteristic length as ℓ = x_f and noting that , the solution is given by

....................(3.385)

For the computation of the pressure responses at the center of the fracture (r_D = 0), Eq. 3.385 simplifies to

....................(3.386)

It is also possible to find a very good approximation for Eq. 3.385, especially when r_eD is large. If we assume^[19]

....................(3.387)

and use the following relation^[4]

....................(3.388)

then Eq. 3.385 may be written as

....................(3.389)

Because^[19]

....................(3.390)

where

....................(3.391)

Eq. 3.389 may also be written as

....................(3.392)

Although the assumption given in Eq. 3.387 may not be justified by itself, the solution given in Eq. 3.392 is a very good approximation for Eq. 3.385, especially when r_eD is large. For a fracture at the center of the cylindrical drainage region, Eq. 3.392 simplifies to

....................(3.393)

It is also possible to obtain short- and long-time approximations for the solution given in Eq. 3.393. For short times, u→∞ and the second term in the argument of the integral in Eq. 3.393 becomes negligible compared with the first term. Then, Eq. 3.393 reduces to the solution for an infinite-slab reservoir given by Eq. 3.367, of which the short-time approximation has been discussed in Example 3.12.

To obtain a long-time approximation, we evaluate Eq. 3.393 at the limit as s→0 (u→s). As shown in Sec. 3.2.3, for small arguments we may approximate the Bessel functions in Eq. 3.393 by

....................(3.394)

....................(3.395)

....................(3.396)

and

....................(3.397)

where γ = 0.5772…. With Eqs. 3.394 through 3.397 and by neglecting the terms of the order s^3/2, we may write^[29]

....................(3.398)

If we substitute the right side of Eq. 3.398 into Eq. 3.393, we obtain

....................(3.399)

where σ(x_D, y_D, −1, +1) is given by Eq. 3.333 and

....................(3.400)

The inverse Laplace transform of Eq. 3.399 yields the following long-time approximation for a uniform-flux fracture at the center of a closed square:

....................(3.401)

Example 3.15

Consider a fully penetrating, uniform-flux fracture of half-length x_f in an isotropic and closed parallelepiped reservoir of dimensions x_e × y_e × h. The fracture is parallel to the x axis and centered at x_w, y_w, z_w.

Solution. The solution for this problem has been obtained in Example 3.10 and, by choosing ℓ = x_f, is given by

....................(3.402)

where

....................(3.403)

The computation of the ratios of the hyperbolic functions in Eq. 3.402 may be difficult, especially when their arguments approach zero or infinity. When s is small (long times), Eq. 3.366 should be useful to compute the ratios of the hyperbolic functions. When s is large (small times), with Eq. 3.366 the solution given in Eq. 3.402 may be written as^[29]

....................(3.404)

where

....................(3.405)

....................(3.406)

and

....................(3.407)

The last equality in Eq. 3.405 follows from the relation given by Eq. 3.349. The expression given in Eq. 3.405 may also be written as

....................(3.408)

where

....................(3.409)

and

....................(3.410)

Therefore, the solution given by Eq. 3.402 may be written as follows for computation at early times (for large values of s):

....................(3.411)

where is given by Eq. 3.409 and corresponds to the solution for a fractured well in an infinite-slab reservoir (see Eq. 3.367 in Example 3.12) and represents the contribution of the lateral boundaries and is given by

....................(3.412)

In Eq. 3.412, , , and are given, respectively, by Eqs. 3.406, 3.407, and 3.410. The integrals appearing in Eqs. 3.409 and 3.410 may be evaluated by following the lines outlined by Eqs. 3.324 through 3.326.

It is also possible to derive short- and long-time approximations for the fracture solution in a closed rectangular parallelepiped. The short-time approximation corresponds to the limit of the solution as s→∞. It can be easily shown that the term in Eq. 3.411 becomes negligible compared with the term for which a short-time approximation has been obtained in Example 3.12 (see Eqs. 3.371 and 3.372).

To obtain a long-time approximation (small values of s), the solution given in Eq. 3.402 may be written as^[27]

....................(3.413)

where

....................(3.414)

and

....................(3.415)

The second equality in Eq. 3.414 results from^[31]

....................(3.416)

For small values of s, replacing u by s and s + α by α, and with^[31]

....................(3.417)

the term H given by Eq. 3.414 may be approximated by

....................(3.418)

The long-time approximation of the second term in Eq. 3.413 is obtained by assuming u << k²π²/x²_eD and taking the inverse Laplace transform of the resulting expressions; therefore, we can obtain the following long-time approximation

....................(3.419)

Example 3.16 Consider a uniform-flux horizontal well of length L_h in an isotropic and closed parallelepiped reservoir of dimensions x_e × y_e × h. The center of the well is at x_w, y_w, z_w, and the well is parallel to the x axis.

Solution. The solution for this problem was obtained in Example 3.11 and, by choosing ℓ = L_h / 2, is given by

....................(3.420)

where is the solution discussed in Example 3.15, and is given by

....................(3.421)

In Eq. 3.421, and are given by Eqs. 3.377 and 3.378, respectively,

....................(3.346)

and

....................(3.422)

The computation and the asymptotic approximations of the term have been discussed in Example 3.15. To compute the term for long times (small s), the relation for the ratios of the hyperbolic functions given by Eq. 3.366 should be useful. For computations at short times (large values of s), following the lines similar to those in Example 3.15, the term in Eq. 3.421 may be written as

....................(3.423)

where

....................(3.424)

....................(3.425)

....................(3.426)

....................(3.427)

and

....................(3.428)

The computational form of the term in Eq. 3.424 is obtained by applying the relations given by Eqs. 3.347 through 3.350 and Eq. 3.328. Of particular interest is the case for −1 ≤ x_D ≤ +1 and y_D = y_wD given by

....................(3.429)

where

....................(3.430)

which can be written as follows by using the relation given in Eq. 3.337:

....................(3.431)

Similarly, for −1 ≤ x_D ≤ +1 and y_D = y_wD, the term given in Eq. 3.428 may be written as

....................(3.432)

where

....................(3.433)

Example 3.15 discussed the short- and long-time approximations of the term in Eq. 3.420. A small-time approximation for given by Eq. 3.423 is obtained with u = ωs and by noting that as s→∞, . Then, substituting the short-time approximations for and given, respectively, by Eqs. 3.371 and 3.353 into Eq. 3.420, the following short-time approximation is obtained: ^[27]

....................(3.434)

where β is given by Eq. 3.331. The inverse Laplace transform of Eq. 3.434 yields

....................(3.435)

The long-time approximation of Eq. 3.420 is obtained by substituting the long-time approximations of and . The long time-approximation of is obtained in Example 3.15 (see Eq. 3.413 through 3.419). The long-time approximation of is obtained by evaluating the right side of Eq. 3.421 as s→0 (u→0) and is given by

....................(3.436)

where

....................(3.437)

and

....................(3.438)

Thus, the long-time approximation Eq. 3.420 is given by

....................(3.439)

where p_Df and F₁ are given, respectively, by Eqs. 3.419 and 3.436. For computational purposes, however, F₁ may be replaced by

....................(3.440)

In Eq. 3.440, F, F_b1, F_b2, and F_b3 are given, respectively, by

....................(3.441)

....................(3.442)

....................(3.443)

and

....................(3.444)

When computing the integrals and the trigonometric series, the relations given by Eqs. 3.324 through 3.326 and 3.345 through 3.350 are useful.

Nomenclature

a	=	radius of the spherical source, L
B	=	formation volume factor, res cm³/std cm³
c	=	fluid compressibility, atm⁻¹
c_f	=	formation compressibility, atm⁻¹
c_t	=	total compressibility, atm⁻¹
C	=	wellbore-storage coefficient, cm³/atm
d	=	distance to a linear boundary, cm
D	=	domain
Ei(x)	=	exponential integral function
f(s)	=	naturally fractured reservoir function
	=	naturally fractured reservoir function based on
	=	Laplace transform of a function f (t)
G	=	Green’s function
h	=	formation thickness, cm
	=	dimensionless thickness, Eq. 3.313
h_D	=	dimensionless thickness, Eq. 3.314
h_f	=	fracture height (vertical penetration), cm
h_p	=	slab thickness, cm
h_w	=	well length (penetration), cm
H(x - x′)	=	Heaviside’s unit step function
	=	unit normal vector in the ξ direction, ξ = x, y, z, r, θ
I_v(x)	=	modified Bessel function of the first kind of order v
I′_v(x)	=	derivative of I_v(x)
J_v(x)	=	Bessel function of the first kind of order v
k	=	isotropic permeability, md
k_f	=	fracture permeability, md
k_h	=	equivalent horizontal permeability, md
k_{i j}	=	permeability in i-direction as a result of pressure gradient in j-direction, md
k_ξ	=	permeability in ξ-direction, ξ = x, y, z, md
k_ξf	=	fracture permeability in ξ-direction, ξ = x, y, z, md
Ki₁(x)	=	first integral of K₀(z)
K_n(x)	=	modified Bessel function of the second kind of order n
K′_n(x)	=	derivative of K_n(x)
ℓ	=	characteristic length of the system, cm
L	=	Laplace transform operator
L^-1	=	inverse Laplace transform operator
L_h	=	horizontal-well length, cm
m	=	pseudopressure, atm²/cp
M_g	=	mass, g
M	=	point in space
M′	=	source point in space
M_w	=	point in Γ_w
M′_w	=	source point in Γ_w
n	=	outward normal direction of the boundary surface
	=	normal vector
N	=	even integer in Stehfest’s algorithm
p	=	pressure, atm
p_c	=	pressure for constant production rate, q_c, atm
	=	dimensionless fracture pressure
p_e	=	external boundary pressure, atm
p _f	=	fracture pressure, atm
p_{f i}	=	initial pressure in fracture system, atm
p_i	=	initial pressure, atm
p_j	=	pressure in medium j, j=m, f, atm
p_m	=	matrix pressure, atm
p_mi	=	initial pressure in matrix system, atm
p_{w f}	=	flowing wellbore pressure, atm
	=	Laplace transform of p(t)
p(t)	=	inverse of the Laplace domain function
p_a(T)	=	approximate inverse of at t=T, atm
q	=	production rate, cm³/s
	=	instantaneous production rate for a point source, cm³/s
q_c	=	constant production rate, cm³/s
q_{s f}	=	sandface production rate, cm³/s
q_wb	=	wellbore production rate as a result of storage, cm³/s
r	=	radial coordinate and distance, cm
r′	=	source coordinate in r-direction, cm
r_e	=	external radius of the reservoir, cm
r_w	=	wellbore radius, cm
R	=	distance in 3D coordinates, cm
R_D	=	dimensionless radial distance in cylindrical coordinates
s	=	Laplace transform parameter
	=	Laplace transform paraeter based on
s_m	=	skin factor
S	=	source function
t	=	time, s
	=	dimensionless time based on k_h
t_AD	=	dimensionless time based on area
t_p	=	producing time, s
T	=	Temperature, °C
u	=	s f(s)
	=	velocity vector
v_ξ	=	velocity component in the ξ direction, ξ = x, y, z, r, θ, cm/s
V	=	volume, cm³
V_i	=	constant in Stehfest’s algorithm
V_f	=	fraction of the volume occupied by fractures
V_m	=	fraction of the volume occupied by matrix
x	=	distance in x-direction, cm
x′	=	source coordinate in x-direction, cm
x_e	=	distance to the external boundary in x-direction, cm
x_p	=	half slab thickness, cm
x_f	=	fracture half-length, cm
	=	dimensionless fracture half-length
x_w	=	well coordinate in x-direction, cm
y	=	distance in y-direction, cm
y′	=	source coordinate in y-direction, cm
y_e	=	distance to the external boundary in y-direction, cm
y_w	=	well coordinate in y-direction, cm
Y_n(x)	=	Bessel function of the second kind of order n
z	=	distance in z-direction, cm
z′	=	source coordinate in z-direction, cm
	=	dimensionless distance in z-direction, Eq. 3.377
z_w	=	well coordinate in z-direction, cm
	=	dimensionless well coordinate in z-direction, Eq. 3.378
Z	=	compressibility factor
α	=	permeability direction, Eq. 3.17
β	=	permeability direction, Eq. 3.17
Γ	=	boundary surface, cm²
Γ_e	=	external boundary surface
Γ_w	=	length, surface, or volume of the source
Γ(x)	=	Gamma function
γ	=	Euler’s constant (γ = 0.5772...)
γ	=	permeability direction, Eq. 3.17
γ _f	=	fundamental solution of diffusion equation
Δ	=	difference operator
δ(x)	=	Dirac delta function
η	=	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
λ	=	transfer coefficient for a naturally fractured reservoir
	=	λ based on k_h
μ	=	viscosity, cp
ρ	=	density, g/cm³
τ	=	time, s
Φ	=	porosity, fraction
φ(M)	=	any continuous function
ω	=	storativity ratio for a naturally fractured reservoir

Subscripts and Superscripts

D	=	dimensionless
f	=	fracture
i	=	initial
m	=	matrix
w	=	wellbore
	=	Laplace transform indicator

References

↑ Al-Hussainy, R., Ramey Jr., H.J., and Crawford, P.B. 1966. The Flow of Real Gases Through Porous Media. J Pet Technol 18 (5): 624–636. SPE-1243-A-PA. http://dx.doi.org/10.2118/1243-A-PA
↑ ^2.0 ^2.1 ^2.2 ^2.3 Raghavan, R. 1993. Well Test Analysis, 28–31, 336–435. Englewood Cliffs, New Jersey: Petroleum Engineering Series, Prentice-Hall.
↑ Fancher, G.H., Lewis, J.A., and Barnes, K.B. 1933. Some Physical Characteristics of Oil Sands. Bulletin 12, Mineral Industries Experimental Station, Pennsylvania State University, University Park, Pennsylvania, 65–167.
↑ ^4.0 ^4.1 Forchheimer, P.F. 1901. Wasserbewegung durch Boden. Zeitschrift des Vereines deutscher Ingenieure 45 (5): 1781–1788.
↑ ^5.0 ^5.1 Watson, G. N. 1944. A Treatise on the Theory of Bessel Functions. London: Cambridge University Press.
↑ ^6.0 ^6.1 ^6.2 ^6.3 Churchill, R.V. 1972. Operational Mathematics, Vol. 2. New York: McGraw-Hill Book Co.
↑ ^7.0 ^7.1 ^7.2 Abramowitz, M. and Stegun, I.A. ed. 1972. Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables, ninth edition, 1020–1029. New York: Dover Publications.
↑ ^8.0 ^8.1 ^8.2 ^8.3 ^8.4 Stehfest, H. 1970. Algorithm 368: Numerical inversion of Laplace transforms. Commun. ACM 13 (1): 47–49. http://dx.doi.org/10.1145/361953.361969
↑ Crump, K.S. 1976. Numerical Inversion of Laplace Transforms Using a Fourier Series Approximation. J. ACM 23 (1): 89-96. http://dx.doi.org/10.1145/321921.321931
↑ ^10.0 ^10.1 van Everdingen, A.F. and Hurst, W. 1953. The Application of the Laplace Transformation to Flow Problems in Reservoirs. In Transactions of the American Institute of Mining and Metallurgical Engineers, Vol. 198, 171. Dallas, Texas: American Institute of Mining and Metallurgical Engineers Inc.
↑ Agarwal, R.G., Al-Hussainy, R., and Ramey Jr., H.J. 1970. An Investigation of Wellbore Storage and Skin Effect in Unsteady Liquid Flow: I. Analytical Treatment. SPE J. 10 (3): 279-290. SPE-2466-PA. http://dx.doi.org/10.2118/2466-PA
↑ Chen, C.-C. and Raghavan, R. 1996. An Approach To Handle Discontinuities by the Stehfest Algorithm. SPE J. 1 (4): 363-368. SPE-28419-PA. http://dx.doi.org/10.2118/28419-PA
↑ ^13.0 ^13.1 ^13.2 ^13.3 ^13.4 ^13.5 ^13.6 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
↑ ^14.0 ^14.1 ^14.2 ^14.3 ^14.4 ^14.5 ^14.6 Carslaw, H.S. and Jaeger, J.C. 1959. Conduction of Heat in Solids, second edition, 353–386. Oxford, UK: Oxford University Press.
↑ ^15.0 ^15.1 Stakgold, I. 1979. Green’s Functions and Boundary Value Problems, 86–104. New York: John Wiley & Sons.
↑ ^16.0 ^16.1 Kelvin, W.T. 1884. Mathematical and Physical Papers, Vol. 2, 41. Cambridge, UK: Cambridge University Press.
↑ ^17.0 ^17.1 ^17.2 ^17.3 Ozkan, E. and Raghavan, R. 1991a. New Solutions for Well-Test-Analysis Problems: Part 1—Analytical Considerations. SPE Form Eval 6 (3): 359–368. SPE-18615-PA. http://dx.doi.org/10.2118/18615-PA
↑ ^18.0 ^18.1 ^18.2 ^18.3 Ozkan, E. and Raghavan, R. 1991b. New Solutions for Well-Test-Analysis Problems: Part 2—Computational Considerations and Applications. SPE Form Eval 6 (3): 369–378. SPE-18616-PA. http://dx.doi.org/10.2118/18616-PA
↑ ^19.00 ^19.01 ^19.02 ^19.03 ^19.04 ^19.05 ^19.06 ^19.07 ^19.08 ^19.09 ^19.10 ^19.11 Raghavan, R. and Ozkan, E. 1994. A Method for Computing Unsteady Flows in Porous Media, No. 318. Essex, England: Pitman Research Notes in Mathematics Series, Longman Scientific & Technical.
↑ ^20.0 ^20.1 Raghavan, R. 1993. The Method of Sources and Sinks. In Well Test Analysis, Chap. 10, 336-435. Englewood Cliffs, New Jersey: Petroleum Engineering Series, Prentice-Hall.
↑ Nisle, R.G. 1958. The Effect of Partial Penetration on Pressure Build-Up in Oil Wells. In Transactions of the American Institute of Mining, Metallurgical, and Petroleum Engineers, Vol. 213, Paper 971-G, 85-90. Dallas, Texas: Society of Petroleum Engineers.
↑ ^22.0 ^22.1 ^22.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
↑ Duhamel, J.M.C. 1833. Mémoire sur la méthode générale relative au mouvement de la chaleur dans les corps solides polongé dans les milieux dont la température varie avec le temps. Journal de l’École Polytechnique 14 (22): 20-66.
↑ Chen, H.Y., Poston, S.W., and Raghavan, R. 1991. An Application of the Product Solution Principle for Instantaneous Source and Green's Functions. SPE Form Eval 6 (2): 161-167. SPE-20801-PA. http://dx.doi.org/10.2118/20801-PA
↑ ^25.0 ^25.1 Barenblatt, G.I., Zheltov, I.P., and Kochina, I.N. 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. J. Appl. Math. Mech. 24 (5): 1286–1303. http://dx.doi.org/10.1016/0021-8928(60)90107-6
↑ ^26.0 ^26.1 Warren, J.E. and Root, P.J. 1963. The Behavior of Naturally Fractured Reservoirs. SPE J. 3 (3): 245–255. SPE-426-PA. http://dx.doi.org/10.2118/426-PA
↑ ^27.0 ^27.1 ^27.2 ^27.3 ^27.4 Kazemi, H. 1969. Pressure Transient Analysis of Naturally Fractured Reservoir with Uniform Fracture Distribution. SPE J. 9 (4): 451–462. SPE-2156-PA. http://dx.doi.org/10.2118/2156-PA
↑ ^28.0 ^28.1 de Swaan O., A. 1976. Analytical Solutions for Determining Naturally Fractured Reservoir Properties by Well Testing. SPE J. 16 (3): 117–122. SPE-5346-PA. http://dx.doi.org/10.2118/5346-PA
↑ ^29.00 ^29.01 ^29.02 ^29.03 ^29.04 ^29.05 ^29.06 ^29.07 ^29.08 ^29.09 ^29.10 ^29.11 ^29.12 ^29.13 ^29.14 ^29.15 ^29.16 ^29.17 ^29.18 Ozkan, E. 1988. Performance of Horizontal Wells. PhD dissertation, University of Tulsa, Tulsa, Oklahoma. <ref name="r30">Luke, Y.L. 1962. Integrals of Bessel Functions, 64–66. New York: McGraw-Hill Book Co.
↑ ^30.0 ^30.1 Cite error: Invalid <ref> tag; no text was provided for refs named r30
↑ ^31.0 ^31.1 Gradshteyn, I.S. and Ryzhik, I.M. 1980. Table of Integrals, Series, and Products, 40. Orlando, Florida: Academic Press.

SI Metric Conversion Factors

atm	×	1.013 250*	E + 05	=	Pa
cp	×	1.0*	E – 03	=	Pa•s
in.	×	2.54*	E + 00	=	cm
in.²	×	6.451 6*	E + 00	=	cm²
°F		(°F−32)/1.8		=	°C
ft	×	3.048*	E – 01	=	m

*Conversion factor is exact.

[r1-1] Al-Hussainy, R., Ramey Jr., H.J., and Crawford, P.B. 1966. The Flow of Real Gases Through Porous Media. J Pet Technol 18 (5): 624–636. SPE-1243-A-PA. http://dx.doi.org/10.2118/1243-A-PA

[r2-2] 2.0 ^2.1 ^2.2 ^2.3 Raghavan, R. 1993. Well Test Analysis, 28–31, 336–435. Englewood Cliffs, New Jersey: Petroleum Engineering Series, Prentice-Hall.

[r3-3] Fancher, G.H., Lewis, J.A., and Barnes, K.B. 1933. Some Physical Characteristics of Oil Sands. Bulletin 12, Mineral Industries Experimental Station, Pennsylvania State University, University Park, Pennsylvania, 65–167.

[r4-4] 4.0 ^4.1 Forchheimer, P.F. 1901. Wasserbewegung durch Boden. Zeitschrift des Vereines deutscher Ingenieure 45 (5): 1781–1788.

[r5-5] 5.0 ^5.1 Watson, G. N. 1944. A Treatise on the Theory of Bessel Functions. London: Cambridge University Press.

[r6-6] 6.0 ^6.1 ^6.2 ^6.3 Churchill, R.V. 1972. Operational Mathematics, Vol. 2. New York: McGraw-Hill Book Co.

[r7-7] 7.0 ^7.1 ^7.2 Abramowitz, M. and Stegun, I.A. ed. 1972. Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables, ninth edition, 1020–1029. New York: Dover Publications.

[r8-8] 8.0 ^8.1 ^8.2 ^8.3 ^8.4 Stehfest, H. 1970. Algorithm 368: Numerical inversion of Laplace transforms. Commun. ACM 13 (1): 47–49. http://dx.doi.org/10.1145/361953.361969

[r9-9] Crump, K.S. 1976. Numerical Inversion of Laplace Transforms Using a Fourier Series Approximation. J. ACM 23 (1): 89-96. http://dx.doi.org/10.1145/321921.321931

[r10-10] 10.0 ^10.1 van Everdingen, A.F. and Hurst, W. 1953. The Application of the Laplace Transformation to Flow Problems in Reservoirs. In Transactions of the American Institute of Mining and Metallurgical Engineers, Vol. 198, 171. Dallas, Texas: American Institute of Mining and Metallurgical Engineers Inc.

[r11-11] Agarwal, R.G., Al-Hussainy, R., and Ramey Jr., H.J. 1970. An Investigation of Wellbore Storage and Skin Effect in Unsteady Liquid Flow: I. Analytical Treatment. SPE J. 10 (3): 279-290. SPE-2466-PA. http://dx.doi.org/10.2118/2466-PA

[r12-12] Chen, C.-C. and Raghavan, R. 1996. An Approach To Handle Discontinuities by the Stehfest Algorithm. SPE J. 1 (4): 363-368. SPE-28419-PA. http://dx.doi.org/10.2118/28419-PA

[r13-13] 13.0 ^13.1 ^13.2 ^13.3 ^13.4 ^13.5 ^13.6 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

[r14-14] 14.0 ^14.1 ^14.2 ^14.3 ^14.4 ^14.5 ^14.6 Carslaw, H.S. and Jaeger, J.C. 1959. Conduction of Heat in Solids, second edition, 353–386. Oxford, UK: Oxford University Press.

[r15-15] 15.0 ^15.1 Stakgold, I. 1979. Green’s Functions and Boundary Value Problems, 86–104. New York: John Wiley & Sons.

[r16-16] 16.0 ^16.1 Kelvin, W.T. 1884. Mathematical and Physical Papers, Vol. 2, 41. Cambridge, UK: Cambridge University Press.

[r17-17] 17.0 ^17.1 ^17.2 ^17.3 Ozkan, E. and Raghavan, R. 1991a. New Solutions for Well-Test-Analysis Problems: Part 1—Analytical Considerations. SPE Form Eval 6 (3): 359–368. SPE-18615-PA. http://dx.doi.org/10.2118/18615-PA

[r18-18] 18.0 ^18.1 ^18.2 ^18.3 Ozkan, E. and Raghavan, R. 1991b. New Solutions for Well-Test-Analysis Problems: Part 2—Computational Considerations and Applications. SPE Form Eval 6 (3): 369–378. SPE-18616-PA. http://dx.doi.org/10.2118/18616-PA

[r19-19] 19.00 ^19.01 ^19.02 ^19.03 ^19.04 ^19.05 ^19.06 ^19.07 ^19.08 ^19.09 ^19.10 ^19.11 Raghavan, R. and Ozkan, E. 1994. A Method for Computing Unsteady Flows in Porous Media, No. 318. Essex, England: Pitman Research Notes in Mathematics Series, Longman Scientific & Technical.

[r20-20] 20.0 ^20.1 Raghavan, R. 1993. The Method of Sources and Sinks. In Well Test Analysis, Chap. 10, 336-435. Englewood Cliffs, New Jersey: Petroleum Engineering Series, Prentice-Hall.

[r21-21] Nisle, R.G. 1958. The Effect of Partial Penetration on Pressure Build-Up in Oil Wells. In Transactions of the American Institute of Mining, Metallurgical, and Petroleum Engineers, Vol. 213, Paper 971-G, 85-90. Dallas, Texas: Society of Petroleum Engineers.

[r22-22] 22.0 ^22.1 ^22.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

[r23-23] Duhamel, J.M.C. 1833. Mémoire sur la méthode générale relative au mouvement de la chaleur dans les corps solides polongé dans les milieux dont la température varie avec le temps. Journal de l’École Polytechnique 14 (22): 20-66.

[r24-24] Chen, H.Y., Poston, S.W., and Raghavan, R. 1991. An Application of the Product Solution Principle for Instantaneous Source and Green's Functions. SPE Form Eval 6 (2): 161-167. SPE-20801-PA. http://dx.doi.org/10.2118/20801-PA

[r25-25] 25.0 ^25.1 Barenblatt, G.I., Zheltov, I.P., and Kochina, I.N. 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. J. Appl. Math. Mech. 24 (5): 1286–1303. http://dx.doi.org/10.1016/0021-8928(60)90107-6

[r26-26] 26.0 ^26.1 Warren, J.E. and Root, P.J. 1963. The Behavior of Naturally Fractured Reservoirs. SPE J. 3 (3): 245–255. SPE-426-PA. http://dx.doi.org/10.2118/426-PA

[r27-27] 27.0 ^27.1 ^27.2 ^27.3 ^27.4 Kazemi, H. 1969. Pressure Transient Analysis of Naturally Fractured Reservoir with Uniform Fracture Distribution. SPE J. 9 (4): 451–462. SPE-2156-PA. http://dx.doi.org/10.2118/2156-PA

[r28-28] 28.0 ^28.1 de Swaan O., A. 1976. Analytical Solutions for Determining Naturally Fractured Reservoir Properties by Well Testing. SPE J. 16 (3): 117–122. SPE-5346-PA. http://dx.doi.org/10.2118/5346-PA

[r29-29] 29.00 ^29.01 ^29.02 ^29.03 ^29.04 ^29.05 ^29.06 ^29.07 ^29.08 ^29.09 ^29.10 ^29.11 ^29.12 ^29.13 ^29.14 ^29.15 ^29.16 ^29.17 ^29.18 Ozkan, E. 1988. Performance of Horizontal Wells. PhD dissertation, University of Tulsa, Tulsa, Oklahoma. <ref name="r30">Luke, Y.L. 1962. Integrals of Bessel Functions, 64–66. New York: McGraw-Hill Book Co.

[r30-30] 30.0 ^30.1 Cite error: Invalid <ref> tag; no text was provided for refs named r30

[r31-31] 31.0 ^31.1 Gradshteyn, I.S. and Ryzhik, I.M. 1980. Table of Integrals, Series, and Products, 40. Orlando, Florida: Academic Press.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

PEH:Mathematics of Transient Analysis

Contents

Introduction

Equations of Transient Fluid Flow in Porous Media

Initial and Boundary Conditions

Assumptions and Limits

Bessel Functions

Preliminary Definitions

Solutions of Bessel’s Equations and Bessel Functions

Modified Bessel Functions of Order Zero and One

Laplace Transformation

Fundamental Properties of the Laplace Transformation

Inverse Laplace Transformation and Asymptotic Forms

Green’s Functions and Source Functions

Fundamental Solution of the Diffusion Equation

The Source-Function Solutions of the Diffusion Equation

The Use of Green’s Functions and Source Functions in Solving Unsteady-Flow Problems

The Use of Source Functions in the Laplace Domain To Solve Unsteady-Flow Problems

Nomenclature

Subscripts and Superscripts

References

SI Metric Conversion Factors

Navigation menu

PEH:Mathematics of Transient Analysis

Introduction

Equations of Transient Fluid Flow in Porous Media

Initial and Boundary Conditions

Assumptions and Limits

Bessel Functions

Preliminary Definitions

Solutions of Bessel’s Equations and Bessel Functions

Modified Bessel Functions of Order Zero and One

Laplace Transformation

Fundamental Properties of the Laplace Transformation

Inverse Laplace Transformation and Asymptotic Forms

Green’s Functions and Source Functions

Fundamental Solution of the Diffusion Equation

The Source-Function Solutions of the Diffusion Equation

The Use of Green’s Functions and Source Functions in Solving Unsteady-Flow Problems

The Use of Source Functions in the Laplace Domain To Solve Unsteady-Flow Problems

Nomenclature

Subscripts and Superscripts

References

SI Metric Conversion Factors

Navigation menu

Search