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

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.

## Contents

- 1 Laplace transform
- 2 Fundamental properties
- 3 Inverse Laplace transformation and asymptotic forms
- 4 Example 1 - Transient flow in a homogeneous reservoir
- 5 Example 2 - Constant rate production in a closed cylindrical reservoir
- 6 Example 3 - Homogeneous reservoir with wellbore storage and skin
- 7 Example 4 - Pressure buildup with wellbore storage and skin following a drawdown period
- 8 Nomenclature
- 9 References
- 10 Noteworthy papers in OnePetro
- 11 External links
- 12 See also
- 13 Category

## Laplace transform

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

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,

## Fundamental properties

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

### Transforms of derivatives

### Transforms of integrals

### Substitution

### Translation

where *H*(*t* - *t*_{0}) is Heaviside’s unit step function defined by

### Convolution

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

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.^{[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, *t*_{i}, where they are defined by

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, *p*_{a}(*T*), of the inverse of the Laplace domain function, , may be obtained at time *t* = *T* by

where

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^{[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, *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^{[5]}

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

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

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,

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.

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

or, rearranging, we obtain

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.

and

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

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

From **Eqs. 24** and **28**, we obtain

which yields

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

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 *r*_{w}→0 and use the relation given in **Eq. 32**, we obtain

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

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

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:

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

we obtain the line-source solution as

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

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

The initial condition corresponding to the uniform pressure distribution equal to *p*_{i} is

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

The closed outer boundary condition is represented mathematically by zero flux at the outer boundary (*r* = *r*_{e}) that corresponds to

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

and

(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

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

which yields

and thus

Using the inner boundary condition given by **Eq. 45** yields

From **Eqs. 49** and **51**, we obtain the coefficients *C*_{1} and *C*_{2} as follows:

and

Substituting *C*_{1} and *C*_{2} into **Eq. 47** yields

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.

In **Eq. 55**, *β*_{1}, *β*_{2}, etc. are the roots of

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

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.

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,

where

and

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.

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:

and

The general solution of **Eq. 63** is

The condition in **Eq. 64** requires that *C*_{1} = 0; therefore,

The use of **Eq. 68** in **Eq. 66** yields

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

which yields

Substituting **Eq. 71** for *C*_{2} in **Eq. 69**, we obtain the solution for the transient pressure distribution in the Laplace transform domain as

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 Δ*p*_{wf}. These derivatives are computed by applying the Laplace transformation property given in **Eq. 3** to **Eq. 72** as follows:

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

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

where *H*(*t* - *t*_{p}) is Heaviside’s unit function (**Eq. 10**), and

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

and

The general solution of **Eq. 80** is

The condition in **Eq. 81** requires that *C*_{1}= 0; therefore,

From **Eqs. 85** and **83**, we obtain

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

which yields

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

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.^{[3]} As suggested by Chen and Raghavan,^{[8]} this problem may be solved by noting that

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

## References

- ↑
^{1.0}^{1.1}^{1.2}^{1.3}Churchill, R.V. 1972. Operational Mathematics, Vol. 2. New York: McGraw-Hill Book Co. - ↑
^{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.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 - ↑ 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
- ↑ 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.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. - ↑ 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

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

Differential calculus refresher

PEH:Mathematics_of_Transient_Analysis