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
No edit summary
No edit summary
 
(One intermediate revision by one other user not shown)
Line 53: Line 53:
[[File:Vol1 page 0088 eq 001.png|RTENOTITLE]]....................(12)
[[File:Vol1 page 0088 eq 001.png|RTENOTITLE]]....................(12)


In this operation, ''p''(''t'') represents the inverse (transform) of the Laplace domain function, [[File:Vol1 page 0088 inline 001.png|RTENOTITLE]]. A uniqueness theorem of the inversion guarantees that no two functions of the class ''ε'' have the same Laplace transform.<ref name="r1">_</ref> 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''<sub>''i''</sub>, where they are defined by
In this operation, ''p''(''t'') represents the inverse (transform) of the Laplace domain function, [[File:Vol1 page 0088 inline 001.png|RTENOTITLE]]. A uniqueness theorem of the inversion guarantees that no two functions of the class ''ε'' have the same Laplace transform.<ref name="r1">Churchill, R.V. 1972. Operational Mathematics, Vol. 2. New York: McGraw-Hill Book Co.</ref> 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''<sub>''i''</sub>, where they are defined by


[[File:Vol1 page 0088 eq 002.png|RTENOTITLE]]....................(13)
[[File:Vol1 page 0088 eq 002.png|RTENOTITLE]]....................(13)
Line 59: Line 59:
and | ''F''(''t'') | < ''Me''<sup>''αt''</sup> for any constants ''M'' and ''α''.
and | ''F''(''t'') | < ''Me''<sup>''αt''</sup> 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,<ref name="r1">_</ref> 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.<ref name="r1">_</ref><ref name="r2">_</ref> 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.
The most rigorous technique to find the inverse Laplace transform of a Laplace domain function is the use of the inversion integral,<ref name="r1">Churchill, R.V. 1972. Operational Mathematics, Vol. 2. New York: McGraw-Hill Book Co.</ref> 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.<ref name="r1">Churchill, R.V. 1972. Operational Mathematics, Vol. 2. New York: McGraw-Hill Book Co.</ref><ref name="r2">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.</ref> 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.


<gallery widths="300px" heights="200px">
<gallery widths="300px" heights="200px">
Line 65: Line 65:
</gallery>
</gallery>


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<ref name="r3">_</ref> is the most popular algorithm.
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<ref name="r3">Stehfest, H. 1970. Algorithm 368: Numerical inversion of Laplace transforms. Commun. ACM 13 (1): 47–49. http://dx.doi.org/10.1145/361953.361969</ref> is the most popular algorithm.


The Stehfest algorithm is based on a stochastic process and suggests that an approximate value, ''p''<sub>''a''</sub>(''T''), of the inverse of the Laplace domain function, [[File:Vol1 page 0088 inline 002.png|RTENOTITLE]], may be obtained at time ''t'' = ''T'' by
The Stehfest algorithm is based on a stochastic process and suggests that an approximate value, ''p''<sub>''a''</sub>(''T''), of the inverse of the Laplace domain function, [[File:Vol1 page 0088 inline 002.png|RTENOTITLE]], may be obtained at time ''t'' = ''T'' by
Line 75: Line 75:
[[File:Vol1 page 0088 eq 004.png|RTENOTITLE]]....................(15)
[[File:Vol1 page 0088 eq 004.png|RTENOTITLE]]....................(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''<sub>''a''</sub> (''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<ref name="r4">_</ref> may be used.
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''<sub>''a''</sub> (''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<ref name="r4">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</ref> 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.
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.
Line 145: Line 145:
[[File:Vol1 page 0091 eq 004.png|RTENOTITLE]]....................(31)
[[File:Vol1 page 0091 eq 004.png|RTENOTITLE]]....................(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<ref name="r3">_</ref> 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 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<ref name="r3">Stehfest, H. 1970. Algorithm 368: Numerical inversion of Laplace transforms. Commun. ACM 13 (1): 47–49. http://dx.doi.org/10.1145/361953.361969</ref> 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.


<gallery widths="300px" heights="200px">
<gallery widths="300px" heights="200px">
Line 161: Line 161:
[[File:Vol1 page 0092 eq 002.png|RTENOTITLE]]....................(34)
[[File:Vol1 page 0092 eq 002.png|RTENOTITLE]]....................(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<ref name="r1">_</ref><ref name="r2">_</ref>), we have
The inversion of '''Eq. 34''' can be accomplished by using a Laplace transform table. From '''Table 1''' (or from the tables in two sources<ref name="r1">Churchill, R.V. 1972. Operational Mathematics, Vol. 2. New York: McGraw-Hill Book Co.</ref><ref name="r2">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.</ref>), we have


[[File:Vol1 page 0092 eq 003.png|RTENOTITLE]]....................(35)
[[File:Vol1 page 0092 eq 003.png|RTENOTITLE]]....................(35)
Line 243: Line 243:
[[File:Vol1 page 0094 eq 006.png|RTENOTITLE]]....................(54)
[[File:Vol1 page 0094 eq 006.png|RTENOTITLE]]....................(54)


The inverse of the solution given by '''Eq. 54''' may not be found in the Laplace transform tables. van Everdingen and Hurst<ref name="r5">_</ref> provided the following analytical inversion of '''Eq. 54''' with the inversion integral.
The inverse of the solution given by '''Eq. 54''' may not be found in the Laplace transform tables. van Everdingen and Hurst<ref name="r5">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.</ref> provided the following analytical inversion of '''Eq. 54''' with the inversion integral.


[[File:Vol1 page 0094 eq 007.png|RTENOTITLE]] [[File:Vol1 page 0095 eq 001.png|RTENOTITLE]]....................(55)
[[File:Vol1 page 0094 eq 007.png|RTENOTITLE]] [[File:Vol1 page 0095 eq 001.png|RTENOTITLE]]....................(55)
Line 251: Line 251:
[[File:Vol1 page 0095 eq 002.png|RTENOTITLE]]....................(56)
[[File:Vol1 page 0095 eq 002.png|RTENOTITLE]]....................(56)


The solution given in '''Eq. 54''' may also be inverted numerically with the Stehfest algorithm.<ref name="r3">_</ref> '''Fig. 2''' shows the results of the numerical inversion of '''Eq. 54'''.
The solution given in '''Eq. 54''' may also be inverted numerically with the Stehfest algorithm.<ref name="r3">Stehfest, H. 1970. Algorithm 368: Numerical inversion of Laplace transforms. Commun. ACM 13 (1): 47–49. http://dx.doi.org/10.1145/361953.361969</ref> '''Fig. 2''' shows the results of the numerical inversion of '''Eq. 54'''.


<gallery widths="300px" heights="200px">
<gallery widths="300px" heights="200px">
Line 263: Line 263:
''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.
''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<ref name="r5">_</ref> (vanishingly small skin-zone radius), the skin factor is defined by
Using van Everdingen and Hurst’s thin-skin concept<ref name="r5">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.</ref> (vanishingly small skin-zone radius), the skin factor is defined by


[[File:Vol1 page 0095 eq 003.png|RTENOTITLE]]....................(57)
[[File:Vol1 page 0095 eq 003.png|RTENOTITLE]]....................(57)
Line 323: Line 323:
[[File:Vol1 page 0097 eq 005.png|RTENOTITLE]]....................(72)
[[File:Vol1 page 0097 eq 005.png|RTENOTITLE]]....................(72)


The real inversion of the solution in '''Eq. 72''' has been obtained by Agarwal ''et al''.<ref name="r6">_</ref> 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.<ref name="r3">_</ref> Also shown in '''Fig. 3''' are the logarithmic derivatives of Δ''p''<sub>''wf''</sub>. These derivatives are computed by applying the Laplace transformation property given in '''Eq. 3''' to '''Eq. 72''' as follows:
The real inversion of the solution in '''Eq. 72''' has been obtained by Agarwal ''et al''.<ref name="r6">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</ref> 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.<ref name="r3">Stehfest, H. 1970. Algorithm 368: Numerical inversion of Laplace transforms. Commun. ACM 13 (1): 47–49. http://dx.doi.org/10.1145/361953.361969</ref> Also shown in '''Fig. 3''' are the logarithmic derivatives of Δ''p''<sub>''wf''</sub>. These derivatives are computed by applying the Laplace transformation property given in '''Eq. 3''' to '''Eq. 72''' as follows:


[[File:Vol1 page 0097 eq 006.png|RTENOTITLE]]....................(73)
[[File:Vol1 page 0097 eq 006.png|RTENOTITLE]]....................(73)
Line 389: Line 389:
[[File:Vol1 page 0099 eq 008.png|RTENOTITLE]]....................(89)
[[File:Vol1 page 0099 eq 008.png|RTENOTITLE]]....................(89)


The [[File:Vol1 page 0100 inline 001.png|RTENOTITLE]] term contributed by the discontinuity at time ''t'' = ''t''<sub>''p''</sub> causes difficulties in the numerical inversion of the right side of '''Eq. 89''' with the use of the Stehfest algorithm.<ref name="r3">_</ref> As suggested by Chen and Raghavan,<ref name="r7">_</ref> this problem may be solved by noting that
The [[File:Vol1 page 0100 inline 001.png|RTENOTITLE]] term contributed by the discontinuity at time ''t'' = ''t''<sub>''p''</sub> causes difficulties in the numerical inversion of the right side of '''Eq. 89''' with the use of the Stehfest algorithm.<ref name="r3">Stehfest, H. 1970. Algorithm 368: Numerical inversion of Laplace transforms. Commun. ACM 13 (1): 47–49. http://dx.doi.org/10.1145/361953.361969</ref> As suggested by Chen and Raghavan,<ref name="r7">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</ref> this problem may be solved by noting that


[[File:Vol1 page 0100 eq 001.png|RTENOTITLE]]....................(90)
[[File:Vol1 page 0100 eq 001.png|RTENOTITLE]]....................(90)
Line 540: Line 540:
[[PEH:Mathematics_of_Transient_Analysis]]
[[PEH:Mathematics_of_Transient_Analysis]]


[[Category:5.6.3 Pressure Transient analysis]]
==Category==
 
[[Category:5.6.3 Pressure transient analysis]] [[Category:YR]]

Latest revision as of 12:51, 6 July 2015

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

RTENOTITLE....................(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,

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

Fundamental properties

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

Transforms of derivatives

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

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

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

Transforms of integrals

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

Substitution

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

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

where RTENOTITLE

Translation

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

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

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

Convolution

RTENOTITLE....................(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,

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

In this operation, p(t) represents the inverse (transform) of the Laplace domain function, RTENOTITLE. 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

RTENOTITLE....................(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, RTENOTITLE, may be obtained at time t = T by

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

where

RTENOTITLE....................(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[5]

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

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

RTENOTITLE....................(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,

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

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

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

or, rearranging, we obtain

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

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

and

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

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

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

RTENOTITLE....................(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 RTENOTITLE, Eq. 23 is satisfied only if C1 = 0]; therefore,

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

From Eqs. 24 and 28, we obtain

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

which yields

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

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

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

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

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

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

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

RTENOTITLE....................(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:

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

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

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

we obtain the line-source solution as

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

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

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

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

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

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

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

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

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

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

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

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

and

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

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

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

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

which yields

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

and thus

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

Using the inner boundary condition given by Eq. 45 yields

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

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

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

and

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

Substituting C1 and C2 into Eq. 47 yields

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

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

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

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

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

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

RTENOTITLE....................(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,

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

where

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

and

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

RTENOTITLE....................(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:

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

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

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

and

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

The general solution of Eq. 63 is

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

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

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

The use of Eq. 68 in Eq. 66 yields

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

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

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

which yields

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

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

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

The real inversion of the solution in Eq. 72 has been obtained by Agarwal et al.[7] 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:

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

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

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

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

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

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

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

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

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

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

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

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

and

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

The general solution of Eq. 80 is

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

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

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

From Eqs. 85 and 83, we obtain

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

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

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

which yields

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

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

The RTENOTITLE 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,[8] this problem may be solved by noting that

RTENOTITLE....................(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
RTENOTITLE = 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
RTENOTITLE = Laplace transform of p(t)
p(t) = inverse of the Laplace domain function
pa(T) = approximate inverse of RTENOTITLE 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
RTENOTITLE = Laplace transform paraeter based on RTENOTITLE
t = time, s
RTENOTITLE = 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. Carslaw, H.S. and Jaeger, J.C. 1986. Conduction of heat in solids, 2nd. Oxford Oxfordshire New York: Clarendon Press ; Oxford University Press. 85026963
  6. 6.0 6.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.
  7. 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
  8. 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

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

Category