You must log in to edit PetroWiki. Help with editing

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


Laplace transformation for solving transient flow problems: Difference between revisions

PetroWiki
Jump to navigation Jump to search
 
Line 523: Line 523:
==Noteworthy papers in OnePetro==
==Noteworthy papers in OnePetro==
Use this section to list papers in OnePetro that a reader who wants to learn more should definitely read
Use this section to list papers in OnePetro that a reader who wants to learn more should definitely read
Carslaw, H.S. and Jaeger, J.C. 1986. Conduction of heat in solids, 2nd. Oxford Oxfordshire New York: Clarendon Press ;
Oxford University Press. 85026963


==External links==
==External links==

Revision as of 14:43, 7 March 2014

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

Vol1 page 0086 eq 010.png....................(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,

Vol1 page 0086 eq 011.png....................(2)

Fundamental properties

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

Transforms of derivatives

Vol1 page 0087 eq 001.png....................(3)

Vol1 page 0087 eq 002.png....................(4)

Vol1 page 0087 eq 003.png....................(5)

Transforms of integrals

Vol1 page 0087 eq 004.png....................(6)

Substitution

Vol1 page 0087 eq 005.png....................(7)

Vol1 page 0087 eq 006.png....................(8)

where Vol1 page 0087 inline 001.png

Translation

Vol1 page 0087 eq 007.png....................(9)

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

Vol1 page 0087 eq 008.png....................(10)

Convolution

Vol1 page 0087 eq 009.png....................(11)

Inverse Laplace transformation and asymptotic forms

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

Vol1 page 0088 eq 001.png....................(12)

In this operation, p(t) represents the inverse (transform) of the Laplace domain function, Vol1 page 0088 inline 001.png. A uniqueness theorem of the inversion guarantees that no two functions of the class ε have the same Laplace transform.[1] 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, ti, where they are defined by

Vol1 page 0088 eq 002.png....................(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,[1] 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.[1][2] 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.

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[3] is the most popular algorithm.

The Stehfest algorithm is based on a stochastic process and suggests that an approximate value, pa(T), of the inverse of the Laplace domain function, Vol1 page 0088 inline 002.png, may be obtained at time t = T by

Vol1 page 0088 eq 003.png....................(14)

where

Vol1 page 0088 eq 004.png....................(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 [pa (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[4] 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, pi.

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

Vol1 page 0089 eq 001.png....................(16)

where ∆p = pip. 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

Vol1 page 0090 eq 001.png....................(17)

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

Vol1 page 0090 eq 002.png....................(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, pi, has been preserved.

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

Vol1 page 0090 eq 003.png....................(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.

Vol1 page 0080 eq 003.png....................(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

Vol1 page 0090 eq 004.png....................(21)

or, rearranging, we obtain

Vol1 page 0090 eq 005.png....................(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.

Vol1 page 0090 eq 006.png....................(23)

and

Vol1 page 0090 eq 007.png....................(24)

Vol1 page 0083 eq 003.png....................(25)

Vol1 page 0084 eq 004.png....................(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

Vol1 page 0090 eq 008.png....................(27)

The constants C1 and C2 in Eq. 27 are obtained from the boundary conditions. The outer boundary condition (Eq. 23) indicates that C1 = 0 [because Vol1 page 0091 inline 001.png, Eq. 23 is satisfied only if C1 = 0]; therefore,

Vol1 page 0091 eq 001.png....................(28)

From Eqs. 24 and 28, we obtain

Vol1 page 0091 eq 002.png....................(29)

which yields

Vol1 page 0091 eq 003.png....................(30)

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

Vol1 page 0091 eq 004.png....................(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[3] 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.

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 rw→0 and use the relation given in Eq. 32, we obtain

Vol1 page 0086 eq 008.png....................(32)

Vol1 page 0092 eq 001.png....................(33)

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

Vol1 page 0092 eq 002.png....................(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[1][2]), we have

Vol1 page 0092 eq 003.png....................(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:

Vol1 page 0092 eq 004.png....................(36)

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

Vol1 page 0092 eq 005.png....................(37)

we obtain the line-source solution as

Vol1 page 0092 eq 006.png....................(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

Vol1 page 0092 eq 007.png....................(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, pi.

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

Vol1 page 0093 eq 001.png....................(40)

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

Vol1 page 0093 eq 002.png....................(41)

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

Vol1 page 0093 eq 003.png....................(42)

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

Vol1 page 0093 eq 004.png....................(43)

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

Vol1 page 0093 eq 005.png....................(44)

Vol1 page 0093 eq 006.png....................(45)

and

Vol1 page 0093 eq 007.png....................(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

Vol1 page 0093 eq 008.png....................(47)

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

Vol1 page 0093 eq 009.png....................(48)

which yields

Vol1 page 0094 eq 001.png....................(49)

and thus

Vol1 page 0094 eq 002.png....................(50)

Using the inner boundary condition given by Eq. 45 yields

Vol1 page 0094 eq 003.png....................(51)

From Eqs. 49 and 51, we obtain the coefficients C1 and C2 as follows:

Vol1 page 0094 eq 004.png....................(52)

and

Vol1 page 0094 eq 005.png....................(53)

Substituting C1 and C2 into Eq. 47 yields

Vol1 page 0094 eq 006.png....................(54)

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

Vol1 page 0094 eq 007.png Vol1 page 0095 eq 001.png....................(55)

In Eq. 55, β1, β2, etc. are the roots of

Vol1 page 0095 eq 002.png....................(56)

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

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

Vol1 page 0095 eq 003.png....................(57)

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

Vol1 page 0096 eq 001.png....................(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, qwb, and the sandface production rate, qsf; that is,

Vol1 page 0096 eq 002.png....................(59)

where

Vol1 page 0096 eq 003.png....................(60)

and

Vol1 page 0096 eq 004.png....................(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.

Vol1 page 0096 eq 005.png....................(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:

Vol1 page 0096 eq 006.png....................(63)

Vol1 page 0096 eq 007.png....................(64)

Vol1 page 0096 eq 008.png....................(65)

and

Vol1 page 0096 eq 009.png....................(66)

The general solution of Eq. 63 is

Vol1 page 0096 eq 010.png....................(67)

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

Vol1 page 0097 eq 001.png....................(68)

The use of Eq. 68 in Eq. 66 yields

Vol1 page 0097 eq 002.png....................(69)

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

Vol1 page 0097 eq 003.png....................(70)

which yields

Vol1 page 0097 eq 004.png....................(71)

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

Vol1 page 0097 eq 005.png....................(72)

The real inversion of the solution in Eq. 72 has been obtained by Agarwal et al.[6] 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.[3] Also shown in Fig. 3 are the logarithmic derivatives of Δpwf. These derivatives are computed by applying the Laplace transformation property given in Eq. 3 to Eq. 72 as follows:

Vol1 page 0097 eq 006.png....................(73)

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

Vol1 page 0097 eq 007.png....................(74)

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

Vol1 page 0098 eq 001.png....................(75)

Vol1 page 0098 eq 002.png....................(76)

Vol1 page 0098 eq 003.png....................(77)

Vol1 page 0098 eq 004.png....................(78)

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

Vol1 page 0098 eq 005.png....................(79)

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

Vol1 page 0098 eq 006.png....................(80)

Vol1 page 0098 eq 007.png....................(81)

Vol1 page 0099 eq 001.png....................(82)

and

Vol1 page 0099 eq 002.png....................(83)

The general solution of Eq. 80 is

Vol1 page 0099 eq 003.png....................(84)

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

Vol1 page 0099 eq 004.png....................(85)

From Eqs. 85 and 83, we obtain

Vol1 page 0099 eq 005.png....................(86)

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

Vol1 page 0099 eq 006.png....................(87)

which yields

Vol1 page 0099 eq 007.png....................(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.

Vol1 page 0099 eq 008.png....................(89)

The Vol1 page 0100 inline 001.png term contributed by the discontinuity at time t = tp causes difficulties in the numerical inversion of the right side of Eq. 89 with the use of the Stehfest algorithm.[3] As suggested by Chen and Raghavan,[7] this problem may be solved by noting that

Vol1 page 0100 eq 001.png....................(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.

Nomenclature

B = formation volume factor, res cm3/std cm3
C = wellbore-storage coefficient, cm3/atm
Vol1 page 0168 inline 002.png = 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
Vol1 page 0169 inline 003.png = Laplace transform of p(t)
p(t) = inverse of the Laplace domain function
pa(T) = approximate inverse of Vol1 page 0169 inline 004.png 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
Vol1 page 0145 inline 001.png = Laplace transform paraeter based on Vol1 page 0170 inline 002.png
t = time, s
Vol1 page 0077 inline 001.png = velocity vector
η = diffusivity constant

References

  1. 1.0 1.1 1.2 1.3 Churchill, R.V. 1972. Operational Mathematics, Vol. 2. New York: McGraw-Hill Book Co.
  2. 2.0 2.1 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.
  3. 3.0 3.1 3.2 3.3 3.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
  4. 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
  5. 5.0 5.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.
  6. 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
  7. 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

Noteworthy papers in OnePetro

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

Carslaw, H.S. and Jaeger, J.C. 1986. Conduction of heat in solids, 2nd. Oxford Oxfordshire New York: Clarendon Press ; Oxford University Press. 85026963

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