Laplace transformation for solving transient flow problems

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 in transient flow.

Laplace transform

The Laplace transformation of a function, F(t), denoted by L{F(t)}, is defined by

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

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

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

Fundamental properties

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

Transforms of derivatives

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

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

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

Transforms of integrals

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

Substitution

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

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

where

Translation

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

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

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

Convolution

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

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,

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

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

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

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, 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 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. The relations given in the Laplace transform tables may be extended to more complex functions with the fundamental properties of the Laplace transforms noted above.

• 

  Table 1 - Laplace transform pairs

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 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

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

where

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

In Eqs. 14 and 15, 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 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 below demonstrate the use of Laplace transformation in the solution of transient flow problems.

Example 1 - Transient flow in a homogeneous reservoir

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

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

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

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

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

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

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,

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

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

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

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

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

or, rearranging, we obtain

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

In obtaining the right side of Eq. 21, the initial condition (Eq. 17) has been used. Similarly, Eqs. 18 and 19 are transformed into the following forms, respectively.

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

and

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

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

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

Comparing Eq. 21 with Eq. 25, we recognize Eq. 22 as the modified Bessel’s equation of order zero. The solution of Eq. 22 may be written directly from Eq. 26 as

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

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

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

From Eqs. 24 and 28, we obtain

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

which yields

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

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

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

To complete the solution of the problem, Eq. 31 should be inverted into the real-time domain. The real inversion of Eq. 31, however, is not available in terms of standard functions. One option is to use Stehfest’s numerical inversion algorithm as discussed above. The dashed line in Fig. 1 represents the numerical inversion of the solution in Eq. 31. 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. 1 – Finite wellbore radius (Eq. 31) and line-source solutions for Example 1.

To obtain the line-source approximation of the solution given in Eq. 31, 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. 32, we obtain

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

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

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

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

The inversion of Eq. 34 can be accomplished by using a Laplace transform table. From Table 1 (or from the tables in two sources), we have

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

With Eq. 35 and the Laplace transform property noted in Eq. 6, we obtain the following inversion of Eq. 34 in the real-time domain:

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

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

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

we obtain the line-source solution as

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

Fig. 1 shows a comparison of the results computed from Eq. 31 (finite-wellbore radius) and Eq. 38 (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. 31) 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. 32, we can show that

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

Example 2 - Constant rate production in a closed cylindrical reservoir

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

....................(40)

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

....................(41)

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

....................(42)

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

....................(43)

The Laplace transforms of Eqs. 40 through 43 yield, respectively,

....................(44)

....................(45)

and

....................(46)

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

....................(47)

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

....................(48)

which yields

....................(49)

and thus

....................(50)

Using the inner boundary condition given by Eq. 45 yields

....................(51)

From Eqs. 49 and 51, we obtain the coefficients C₁ and C₂ as follows:

....................(52)

and

....................(53)

Substituting C₁ and C₂ into Eq. 47 yields

....................(54)

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

....................(55)

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

....................(56)

The solution given in Eq. 54 may also be inverted numerically with the Stehfest algorithm. Fig. 2 shows the results of the numerical inversion of Eq. 54.

• 

  Fig. 2 – Bounded reservoir solution for Example 2

Example 3 - Homogeneous reservoir with wellbore storage and skin

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 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 (vanishingly small skin-zone radius), the skin factor is defined by

....................(57)

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. 57, we obtain the following relation for the flowing wellbore pressure.

....................(58)

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,

....................(59)

where

....................(60)

and

....................(61)

In Eq. 60, C is the wellbore-storage coefficient. Substituting Eqs. 60 and 61 into Eq. 59, we obtain the following expression for the surface production rate.

....................(62)

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

....................(63)

....................(64)

....................(65)

and

....................(66)

The general solution of Eq. 63 is

....................(67)

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

....................(68)

The use of Eq. 68 in Eq. 66 yields

....................(69)

From Eqs. 65, 68, and 69, we obtain

....................(70)

which yields

....................(71)

Substituting Eq. 71 for C₂ in Eq. 69, we obtain the solution for the transient pressure distribution in the Laplace transform domain as

....................(72)

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

....................(73)

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

....................(74)

• 

  Fig. 3 – Wellbore-storage and skin solution (Eq. 72) for Example 3.

Example 4 - Pressure buildup with wellbore storage and skin following a drawdown period

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 except, at time t_p, the well is shut in and pressure buildup begins. The system of equations to define this problem is

....................(75)

....................(76)

....................(77)

....................(78)

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

....................(79)

The right side of the boundary condition in Eq. 78 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. 75 through 79, we obtain

....................(80)

....................(81)

....................(82)

and

....................(83)

The general solution of Eq. 80 is

....................(84)

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

....................(85)

From Eqs. 85 and 83, we obtain

....................(86)

Substituting the results of Eqs. 85 and 86 into Eq. 82, we have

....................(87)

which yields

....................(88)

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

....................(89)

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

....................(90)

and applying the Stehfest algorithm term by term to the right side of Eq. 90. Fig. 4 shows sample results obtained by the numerical inversion of Eq. 89.

• 

  Fig. 4 – Drawdown and buildup results with wellbore-storage and skin solution for Example 4.

Nomenclature

B	=	formation volume factor, res cm3/std cm3
C	=	wellbore-storage coefficient, cm3/atm
	=	Laplace transform of a function f (t)
h	=	formation thickness, cm
H(x - x′)	=	Heaviside’s unit step function
k	=	isotropic permeability, md
Ki1(x)	=	first integral of K0(z)
Kn(x)	=	modified Bessel function of the second kind of order n
K′n(x)	=	derivative of Kn(x)
L	=	Laplace transform operator
L-1	=	inverse Laplace transform operator
p	=	pressure, atm
pw f	=	flowing wellbore pressure, atm
	=	Laplace transform of p(t)
p(t)	=	inverse of the Laplace domain function
pa(T)	=	approximate inverse of at t=T, atm
q	=	production rate, cm3/s
qs f	=	sandface production rate, cm3/s
qwb	=	wellbore production rate as a result of storage, cm3/s
r	=	radial coordinate and distance, cm
re	=	external radius of the reservoir, cm
rw	=	wellbore radius, cm
s	=	Laplace transform parameter
	=	Laplace transform paraeter based on
t	=	time, s
	=	velocity vector
η	=	diffusivity constant

Noteworthy papers in OnePetro

Use this section to list papers in OnePetro that a reader who wants to learn more should definitely read

Chen, H.Y., Poston, S.W., and Raghavan, R. An Application of the Product Solution Principle for Instantaneous Source and Green's Functions. http://dx.doi.org/10.2118/20801-PA. Gringarten, A.C. and Ramey, H.J., Jr. The Use of Source and Green's Functions in Solving Unsteady-Flow Problems in Reservoirs. http://dx.doi.org/10.2118/3818-PA Ozkan, E. and Raghavan, R. New Solutions for Well-Test-Analysis Problems: Part 1-Analytical Considerations(includes associated papers 28666 and 29213 ). http://dx.doi.org/10.2118/18615-PA. Van Everdingen, A.F. and Hurst, W. The Application of the Laplace Transformation to Flow Problems in Reservoirs. http://dx.doi.org/10.2118/949305-G.

External links

Use this section to provide links to relevant material on websites other than PetroWiki and OnePetro

See also

Transient analysis mathematics

Source function solutions of the diffusion equation

Solving unsteady flow problems with Laplace transform and source functions

Mathematics of fluid flow

Differential calculus refresher

PEH:Mathematics_of_Transient_Analysis