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


Diagonalizing the permeability tensor

PetroWiki
Revision as of 09:11, 10 April 2014 by Petrowiki User (Petrowikistaff) (talk | contribs) (Added links to books)
Jump to navigation Jump to search

The form of Darcy’s law that is most widely used in formulating fluid flow equations in reservoir simulators assumes that the coordinate system is aligned with the principal axes of the permeability tensor. The resulting diagonalized permeability greatly simplifies the fluid flow equations. The simplified equations are easier to code and can be solved with less computation time than fluid flow equations that include the full permeability tensor.

Research in naturally-fractured-reservoir modeling,[1] geomechanics,[2] and upscaling[3] has demonstrated that the full permeability tensor is needed to correctly solve fluid flow problems in a variety of realistic settings. The mathematical procedure for diagonalizing the permeability tensor is presented here as an illustration of the mathematics discussed in matrices and linear algebra. The relationship between the diagonalized-permeability-tensor assumption and grid orientation is discussed later on this page. An understanding of the relationship between grid orientation and the permeability tensor can help us decide how to orient a fluid flow grid to most accurately represent the permeability distribution in a reservoir. The directional dependence of permeability and the permeability tensor are first introduced. The procedure for diagonalizing the permeability tensor is then presented.

Darcy’s law and the permeability tensor

In one dimension, Darcy’s law says that flow rate is proportional to pressure gradient. This can be expressed in oilfield units for single-phase flow as

Vol1 page 0063 eq 001.png....................(1)

where Vol1 page 0063 inline 001.png is flow rate (B/D), x is length (ft), A is cross-sectional area (ft2), μ is fluid viscosity (cp), k is permeability (md), and Φ is the phase potential (psia).

Vol1 page 0063 eq 002.png....................(2)

In Eq. 2, Δz is depth from a datum (ft), P is fluid pressure (psia), and γ is the pressure gradient associated with the gravity term (psia/ft). The form of Darcy’s law with full permeability tensor in Cartesian coordinate system {x, y, z} is

Vol1 page 0063 eq 003.png....................(3)

where we have treated the cross-sectional area, A, as a constant with respect to direction. Eq. 3 can be rewritten as either a dyadic equation,

Vol1 page 0063 eq 004.png....................(4)

by treating permeability as a dyadic, or as a matrix equation,

Vol1 page 0063 eq 005.png....................(5)

by treating permeability as a matrix. We are interested here in the matrix representation.

The diagonal permeability elements {kxx, kyy, kzz} represent the dependence of flow rate in one direction on pressure differences in the same direction. The off-diagonal permeability elements {kxy, kxz, kyx, kyz, kzx, kzy} account for the dependence of flow rate in one direction on pressure differences in orthogonal directions. Expanding Eq. 3 gives the corresponding set of three equations demonstrating this dependence.

Vol1 page 0063 eq 006.png....................(6)

Similarity transformations

Eq. 5 relates flow rate Vol1 page 0064 inline 001.png and the pressure gradient term, Vol1 page 0064 inline 002.png. We can use a similarity transformation to diagonalize the matrix Vol1 page 0064 inline 003.png while preserving the form of the relationship between Vol1 page 0064 inline 004.png and Vol1 page 0064 inline 002.png. Let us first show that a similarity transformation preserves the form of Eq. 5.

We begin by multiplying Eq. 5 on the left by Vol1 page 0059 inline 006.png to find

Vol1 page 0064 eq 001.png....................(7)

where Vol1 page 0059 inline 002.png is a nonsingular, n × n square matrix. Because Vol1 page 0059 inline 002.png is nonsingular, it is invertible; that is, it satisfies the equality

Vol1 page 0064 eq 002.png....................(8)

where Vol1 page 0060 inline 006.png is the n × n identity matrix. Substituting Eq. 8 into Eq. 7 gives

Vol1 page 0064 eq 003.png....................(9)

Defining the transformed matrices

Vol1 page 0064 eq 004.png....................(10)

and using the similarity transformation

Vol1 page 0064 eq 005.png....................(11)

in Eq. 9 yields

Vol1 page 0064 eq 006.png....................(12)

Eq. 12 is the same form as Eq. 5.

Matrix diagonalization procedure

It is mathematically possible to find a coordinate system {x′, y′, z′} in which the permeability tensor has the diagonal form Vol1 page 0064 inline 005.png. We diagonalize the matrix Vol1 page 0064 inline 006.png by finding and applying a similarity transformation matrix Vol1 page 0058 inline 003.png. The procedure for finding a matrix Vol1 page 0058 inline 003.png that diagonalizes an n × n matrix Vol1 page 0064 inline 007.png is as follows:[4]

  • Find the eigenvalues {λi: i = 1, ..., n} of Vol1 page 0064 inline 007.png from the eigenvalue equation detVol1 page 0064 inline 008.png = 0.
  • Find n linearly independent eigenvectors Vol1 page 0064 inline 009.png .
  • Form the similarity transformation matrix Vol1 page 0059 inline 002.png with the eigenvectors as columns.
  • Calculate the diagonalized matrix Vol1 page 0064 inline 006.png´. The diagonal entries of Vol1 page 0064 inline 006.png´ are the eigenvalues corresponding to the eigenvectors Vol1 page 0065 inline 001.png .

The coordinate axes {x′, y′, z′} are the principal axes of the diagonalized tensor, and the diagonal form of the permeability tensor is obtained by a principal-axis transformation. The flow equations along the principal axes are

Vol1 page 0065 eq 001.png....................(13)

The form of the permeability tensor depends on the properties of the porous medium. If the medium is anisotropic, at least two elements of the diagonalized permeability tensor are not equal. If permeability does not depend on direction, then permeability is isotropic, and the elements of the diagonalized permeability tensor are equal, that is,

Vol1 page 0065 eq 002.png....................(14)

If the magnitude of the elements of the permeability tensor varies from one point in the medium to another, the permeability tensor is heterogeneous; otherwise, permeability is homogeneous. The principal axes of the permeability tensor may also vary from point to point in the medium if permeability is heterogeneous.

Diagonalizing a symmetric 2 x 2 matrix

The ideas previously presented are implemented by applying the matrix diagonalization algorithm to the 2 × 2 symmetric matrix

Vol1 page 0065 eq 003.png....................(15)

as viewed in the two-dimensional (2D) Cartesian coordinate system Vol1 page 0059 inline 001.png ={ x1,x2} shown in Fig. 1. For this example, we require that the elements of Vol1 page 0064 inline 003.png satisfy the relations

Vol1 page 0065 eq 004.png....................(16)

The relation k12 = k21 for off-diagonal elements is necessary to assure that the matrix Vol1 page 0064 inline 007.png is symmetric. The requirement that Vol1 page 0064 inline 006.png is symmetric is important when we consider a coordinate transformation. To find the diagonal matrix Vol1 page 0065 inline 002.png´ corresponding to Vol1 page 0065 inline 002.png, we must first solve the eigenvalue problem

Vol1 page 0065 eq 005.png....................(17)

The two characteristic roots or eigenvalues λ+ and λ of Eq. 17 are the diagonal elements of the diagonalized 2 × 2 matrix Vol1 page 0065 inline 002.png´, thus

Vol1 page 0066 eq 001.png....................(18)

where λ+ and λ are calculated as the solutions to the eigenvalue problem.

The eigenvalue problem is detVol1 page 0066 inline 001.png = 0. Using Eq. 15 gives

Vol1 page 0066 eq 002.png....................(19)

or

Vol1 page 0066 eq 003.png....................(20)

We expand the characteristic equation to get

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

The two eigenvalues are found from the quadratic equation to be

Vol1 page 0066 eq 005.png....................(22)

The sum of the eigenvalues satisfies the relation

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

Eigenvectors

The matrix Vol1 page 0059 inline 006.png is composed of orthonormal eigenvectors Vol1 page 0066 inline 002.png found from

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

The basis vector, Vol1 page 0067 inline 001.png, satisfies

Vol1 page 0066 eq 008.png....................(25)

with the identity matrix Vol1 page 0066 inline 003.png. Expanding Eq. 25 gives

Vol1 page 0066 eq 009.png....................(26)

for the eigenvalue λ+, and

Vol1 page 0066 eq 010.png Vol1 page 0067 eq 001.png....................(27)

for the eigenvalue λ. Rearranging Eq. 26 gives

Vol1 page 0067 eq 002.png....................(28)

Eq. 28 and the normalization condition,

Vol1 page 0067 eq 003.png....................(29)

provide the two equations that are necessary for determining the components of Vol1 page 0067 inline 001.png+; thus,

Vol1 page 0067 eq 004.png....................(30)

and

Vol1 page 0067 eq 005.png....................(31)

Similar calculations for Vol1 page 0067 inline 001.png- yield the results

Vol1 page 0067 eq 006.png....................(32)

and

Vol1 page 0067 eq 007.png....................(33)

where the relation

Vol1 page 0067 eq 008.png....................(34)

from Eq. 33, has been used.

To show that Vol1 page 0067 inline 001.png+ and Vol1 page 0067 inline 001.png- are orthogonal, we must show that

Vol1 page 0067 eq 009.png....................(35)

Substituting in the expressions for Vol1 page 0067 inline 002.png+ and Vol1 page 0067 inline 002.png- gives

Vol1 page 0068 eq 001.png....................(36)

as expected.

Coordinate transformation

We now use the orthonormal eigenvectors Vol1 page 0067 inline 002.png+ and Vol1 page 0067 inline 002.png- to construct the transformation matrix Vol1 page 0059 inline 002.png. According to the algorithm for diagonalizing a square matrix presented previously, we form Vol1 page 0059 inline 002.png as

Vol1 page 0068 eq 002.png....................(37)

or

Vol1 page 0068 eq 003.png....................(38)

A coordinate vector in the transformed coordinate system Vol1 page 0068 inline 002.png is given by Vol1 page 0068 inline 003.png. Rewriting the matrix equation for coordinate transformations in algebraic form gives

Vol1 page 0068 eq 004.png....................(39)

or

Vol1 page 0068 eq 005.png....................(40)

An angle (θ) can be associated with the linear transformation by writing the 2D coordinate transformation as

Vol1 page 0068 eq 006.png....................(41)

The coordinate systems Vol1 page 0068 inline 004.png and Vol1 page 0068 inline 005.png are related by the counterclockwise rotation shown in Fig. 1.

Equating elements of the transformation matrix in Eqs. 40 and 41 gives

Vol1 page 0068 eq 007.png....................(42)

and

Vol1 page 0069 eq 001.png....................(43)

The equalities a1+ = a2 and a2+ = –a1 are in agreement with Eqs. 30 through 33. The angle θ may be found from either

Vol1 page 0069 eq 002.png....................(44)

or

Vol1 page 0069 eq 003.png....................(45)

Note that the equality,

Vol1 page 0069 eq 004.png....................(46)

demonstrates that the eigenvectors are orthonormal.

Rotational transformation of a 2 x 2 permeability tensor

We want to calculate changes to the permeability tensor when we transform from a coordinate system Vol1 page 0068 inline 002.png where only the diagonal elements of a square matrix Vol1 page 0064 inline 006.png´ are nonzero to a coordinate system Vol1 page 0068 inline 004.png in which a 2 × 2 square matrix Vol1 page 0064 inline 006.png has nonzero off-diagonal elements. We do this by performing a similarity transformation on the matrix Vol1 page 0064 inline 006.png. The coordinate systems Vol1 page 0068 inline 004.png and Vol1 page 0068 inline 005.png are related by the similarity transformation matrix Vol1 page 0059 inline 002.png such that

Vol1 page 0069 eq 005.png....................(47)

The two coordinate systems are shown in Fig. 1.

An angle (θ) is associated with the transformation in Eq. 47 by writing the 2D coordinate transformation as

Vol1 page 0069 eq 006.png....................(48)

The coordinate systems Vol1 page 0068 inline 004.png and Vol1 page 0068 inline 005.png are related by the counterclockwise rotation shown in Fig. 1. We have an aligned coordinate system Vol1 page 0068 inline 002.png with the principal axes of the permeability tensor. The diagonal tensor in the coordinate system Vol1 page 0068 inline 005.png has the form

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

where kmax is the maximum permeability in the direction y1 and kT is the permeability that is transverse to k max in the direction y2. We want to know how the elements of the permeability tensor change if we transform to the different coordinate system Vol1 page 0070 inline 001.png.

The relationship between the elements of Vol1 page 0070 inline 002.png and Vol1 page 0064 inline 003.png is

Vol1 page 0070 eq 001.png....................(50)

where Vol1 page 0059 inline 002.png is

Vol1 page 0070 eq 002.png....................(51)

If we multiply Eq. 50 from the left by Vol1 page 0058 inline 003.png-1 and from the right by Vol1 page 0058 inline 003.png, we obtain

Vol1 page 0070 eq 003.png....................(52)

We find the elements of Vol1 page 0059 inline 002.png-1 by solving

Vol1 page 0070 eq 004.png....................(53)

The result is

Vol1 page 0070 eq 005.png....................(54)

where Vol1 page 0070 inline 003.png is the transpose of Vol1 page 0059 inline 002.png. Substituting Eqs. 51 and 54 into 52 gives

Vol1 page 0070 eq 006.png....................(55)

We can use Eq. 55 to calculate the elements of Vol1 page 0065 inline 002.png for any rotation angle θ. If the permeability is isotropic, we have kmax = kT = kiso and Eq. 55 simplifies to the form

Vol1 page 0070 eq 007.png....................(56)

In the special case of isotropic permeability, the orientation of the coordinate system does not affect the values of the elements of the permeability tensor.

Fig. 2 shows the results for an anisotropic case in which kmax = 200 md and kT = 50 md. The values of elements k11, k12 and k22 of Vol1 page 0064 inline 006.png are presented in the figure. The off-diagonal terms satisfy the equality k12 = k21 for a symmetric matrix given in Eq. 16, so it is sufficient to show only k12. The values of the diagonal elements change most as θ approaches 45°, and the values of the off-diagonal elements are greatest at θ = 45°. A rotation of 90° recovers a diagonal permeability tensor, but k max is now aligned along x2, and kT is aligned along x1.

Gridding a channel sand

The ideas discussed are now considered in the context of a realistic application. Our problem is to find a coordinate system that lets us accurately model fluid flow in a channel sand with the two permeability regions in Fig. 3.

We can highlight important features of the relationship between grid orientation and the assumption of diagonalized permeability by assuming the permeability in each region is anisotropic with a maximum permeability k max and a permeability kT that is transverse to the direction of kmax. The diagonalized, anisotropic permeability tensor Vol1 page 0064 inline 006.png′ in the y1 - y2 plane of a channel sand is the matrix

Vol1 page 0072 eq 001.png....................(57)

and Darcy’s law for flow in the y1 - y2 plane is

Vol1 page 0072 eq 002.png....................(58)

Consider two cases. In Case A, permeability is homogeneous and anisotropic, and in Case B, permeability is heterogeneous and anisotropic. The two cases are illustrated in Figs. 4 and 5.

Fig. 6 shows two possible coordinate systems for orienting the grid. Coordinate system Vol1 page 0068 inline 005.png is more closely aligned with the spatial orientation of Region I than coordinate system Vol1 page 0070 inline 001.png, while coordinate system Vol1 page 0070 inline 001.png is more closely aligned with the spatial orientation of Region II than coordinate system Vol1 page 0068 inline 002.png. The coordinate system Vol1 page 0068 inline 005.png is obtained by rotating the coordinate system Vol1 page 0068 inline 004.png through an angle θ as in Fig. 1.

We consider four grid orientations for each case: (1) grid Vol1 page 0068 inline 002.png in Regions I and II; (2) grid Vol1 page 0068 inline 005.png in Region I and grid Vol1 page 0068 inline 004.png in Region II; (3) grid Vol1 page 0068 inline 004.png in Region I and grid Vol1 page 0068 inline 002.png in Region II; and (4) grid Vol1 page 0068 inline 002.png in Regions I and II. The grid orientation cases allow us to consider the effect of different coordinate systems on the permeability tensor in each region. We assume in our analysis that the reservoir simulator is a typical simulator with a formulation of fluid flow equations that uses Darcy’s law with a diagonal permeability tensor. We also assume the simulator allows different grid orientations in different regions of the model; otherwise, grid orientation cases 2 and 3 are not feasible. Our analysis does not include multidimensional numerical dispersion,[5] which can also affect the accuracy of flow calculations. Results of the analysis are summarized in Table 1.

An "ok" in the "Permeability Tensor" column in Table 1 indicates that the diagonal permeability tensor is aligned with the grid. An "X" indicates that the magnitudes of the diagonal terms in the permeability tensor must be corrected with Eq. 55. An "ok" in the "Formulation" column in Table 1 indicates that the formulation of the fluid flow equations is correct. An "X" indicates that the formulation of the fluid flow equations is incorrect because the formulation does not include off-diagonal terms in the permeability tensor. Based on the results in Table 1, we observe that the grid orientation in Case A.1 provides the most faithful representation of the permeability tensor in Case A, and the grid orientation in Case B.2 provides the most faithful representation of the permeability tensor in Case B.

Nomenclature

Vol1 page 0059 inline 002.png = rotation matrix
Vol1 page 0074 inline 002.png = transpose of matrix A , Eq. 54
A = cross-sectional area, Eq. 1
Vol1 page 0074 inline 006.png = identity matrix
Jx, Jy, Jz = fluid flux in x-, y-, z-directions
(Jx)x = fluid flux in x-direction at location x
(Jy)y = fluid flux in y-direction at location y
(Jz)z = fluid flux in z-direction at location z
Vol1 page 0064 inline 006.png = permeability matrix, Eq. 5
k = permeability, Eq. 1
kiso = isotropic permeability, Eq. 56
kmax = maximum permeability, Eq. 49
kT = transverse permeability, Eq. 49
P = pressure, Eq. 2
Vol1 page 0074 inline 007.png = flow rate, Eqs. 1 and 5
Vol1 page 0075 inline 002.png = column vectors, Eq. 41
Δt = time interval
Δx = length
Δy = width
Δz = thickness
θ = angle
λ = eigenvalues
μ = fluid viscosity, Eq. 1
Vol1 page 0064 inline 002.png = pressure gradient, Eq. 5
Φ = phase potential, Eq. 1

References

  1. Gupta, A., Penuela, G., and Avila, R. 2001. An Integrated Approach to the Determination of Permeability Tensors for Naturally Fractured Reservoirs. J Can Pet Technol 40 (12): 43. PETSOC-01-12-02. http://dx.doi.org/10.2118/01-12-02
  2. Settari, A., Walters, D.A., and Behie, G.A. 2001. Use of Coupled Reservoir and Geomechanical Modelling for Integrated Reservoir Analysis and Management. J Can Pet Technol 40 (12): 55. PETSOC-01-12-04. http://dx.doi.org/10.2118/01-12-04
  3. Young, L.C. 1999. Rigorous Treatment of Distorted Grids in 3D. Presented at the SPE Reservoir Simulation Symposium, Houston, 14-17 February. SPE-51899-MS. http://dx.doi.org/10.2118/51899-MS
  4. Fanchi, J.R. 2006. Math Refresher for Scientists and Engineers, third edition. New York: Wiley Interscience.
  5. Fanchi, J.R. 1983. Multidimentional Numerical Dispersion. SPE J. 23 (1): 143–151. SPE-9018-PA. http://dx.doi.org/10.2118/9018-PA

Noteworthy papers in OnePetro

Ramey, H.J., Jr. Interference Analysis for Anisotropic Formations - A Case History (includes associated paper 6406 ). http://dx.doi.org/10.2118/5319-PA.

Noteworthy books

Carslaw, H., & Jaeger, J. (1986). Conduction of heat in solids (2nd ed.). Oxford Oxforshire, New York: Clarendon Press. 85026963 Worldcat

Earlougher, R. C. 1977. Advances in Well Test Analysis, SPE Monograph Series Vol. 5, Society of Petroleum Engineers, Richardson, TX, 264 pp. SPE Bookstore or Worldcat

Hahn, D., & Özisik, M. (2012). Heat Conduction. (3rd). Hoboken, New Jersey: John Wiley & Sons, Inc. Worldcat

Kreyszig, E. O. (2011). Advanced Engineering Mathematics. (10th). New York, New York: John Wiley & Sons Inc. Online Resource or Worldcat

External links

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

See also

Mathematics of fluid flow

Numerical methods analysis of fluid flow

Vector analysis of fluid flow

PEH:Mathematics of Fluid Flow