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

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

where RTENOTITLE 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).

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

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

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

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

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

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

Similarity transformations

Eq. 5 relates flow rate RTENOTITLE and the pressure gradient term, RTENOTITLE. We can use a similarity transformation to diagonalize the matrix RTENOTITLE while preserving the form of the relationship between RTENOTITLE and RTENOTITLE. Let us first show that a similarity transformation preserves the form of Eq. 5.

We begin by multiplying Eq. 5 on the left by RTENOTITLE to find

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

where RTENOTITLE is a nonsingular, n × n square matrix. Because RTENOTITLE is nonsingular, it is invertible; that is, it satisfies the equality

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

where RTENOTITLE is the n × n identity matrix. Substituting Eq. 8 into Eq. 7 gives

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

Defining the transformed matrices

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

and using the similarity transformation

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

in Eq. 9 yields

RTENOTITLE....................(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 RTENOTITLE. We diagonalize the matrix RTENOTITLE by finding and applying a similarity transformation matrix RTENOTITLE. The procedure for finding a matrix RTENOTITLE that diagonalizes an n × n matrix RTENOTITLE is as follows:[4]

  • Find the eigenvalues {λi: i = 1, ..., n} of RTENOTITLE from the eigenvalue equation detRTENOTITLE = 0.
  • Find n linearly independent eigenvectors RTENOTITLE .
  • Form the similarity transformation matrix RTENOTITLE with the eigenvectors as columns.
  • Calculate the diagonalized matrix RTENOTITLE´. The diagonal entries of RTENOTITLE´ are the eigenvalues corresponding to the eigenvectors RTENOTITLE .

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

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

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

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

as viewed in the two-dimensional (2D) Cartesian coordinate system RTENOTITLE ={ x1,x2} shown in Fig. 1. For this example, we require that the elements of RTENOTITLE satisfy the relations

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

The relation k12 = k21 for off-diagonal elements is necessary to assure that the matrix RTENOTITLE is symmetric. The requirement that RTENOTITLE is symmetric is important when we consider a coordinate transformation. To find the diagonal matrix RTENOTITLE´ corresponding to RTENOTITLE, we must first solve the eigenvalue problem

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

The two characteristic roots or eigenvalues λ+ and λ of Eq. 17 are the diagonal elements of the diagonalized 2 × 2 matrix RTENOTITLE´, thus

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

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

The eigenvalue problem is detRTENOTITLE = 0. Using Eq. 15 gives

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

or

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

We expand the characteristic equation to get

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

The two eigenvalues are found from the quadratic equation to be

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

The sum of the eigenvalues satisfies the relation

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

Eigenvectors

The matrix RTENOTITLE is composed of orthonormal eigenvectors RTENOTITLE found from

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

The basis vector, RTENOTITLE, satisfies

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

with the identity matrix RTENOTITLE. Expanding Eq. 25 gives

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

for the eigenvalue λ+, and

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

for the eigenvalue λ. Rearranging Eq. 26 gives

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

Eq. 28 and the normalization condition,

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

provide the two equations that are necessary for determining the components of RTENOTITLE+; thus,

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

and

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

Similar calculations for RTENOTITLE- yield the results

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

and

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

where the relation

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

from Eq. 33, has been used.

To show that RTENOTITLE+ and RTENOTITLE- are orthogonal, we must show that

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

Substituting in the expressions for RTENOTITLE+ and RTENOTITLE- gives

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

as expected.

Coordinate transformation

We now use the orthonormal eigenvectors RTENOTITLE+ and RTENOTITLE- to construct the transformation matrix RTENOTITLE. According to the algorithm for diagonalizing a square matrix presented previously, we form RTENOTITLE as

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

or

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

A coordinate vector in the transformed coordinate system RTENOTITLE is given by RTENOTITLE. Rewriting the matrix equation for coordinate transformations in algebraic form gives

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

or

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

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

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

The coordinate systems RTENOTITLE and RTENOTITLE are related by the counterclockwise rotation shown in Fig. 1.

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

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

and

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

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

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

or

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

Note that the equality,

RTENOTITLE....................(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 RTENOTITLE where only the diagonal elements of a square matrix RTENOTITLE´ are nonzero to a coordinate system RTENOTITLE in which a 2 × 2 square matrix RTENOTITLE has nonzero off-diagonal elements. We do this by performing a similarity transformation on the matrix RTENOTITLE. The coordinate systems RTENOTITLE and RTENOTITLE are related by the similarity transformation matrix RTENOTITLE such that

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

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

The coordinate systems RTENOTITLE and RTENOTITLE are related by the counterclockwise rotation shown in Fig. 1. We have an aligned coordinate system RTENOTITLE with the principal axes of the permeability tensor. The diagonal tensor in the coordinate system RTENOTITLE 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 RTENOTITLE.

The relationship between the elements of RTENOTITLE and RTENOTITLE is

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

where RTENOTITLE is

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

If we multiply Eq. 50 from the left by RTENOTITLE-1 and from the right by RTENOTITLE, we obtain

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

We find the elements of RTENOTITLE-1 by solving

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

The result is

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

where RTENOTITLE is the transpose of RTENOTITLE. Substituting Eqs. 51 and 54 into 52 gives

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

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

RTENOTITLE....................(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 RTENOTITLE 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 RTENOTITLE′ in the y1 - y2 plane of a channel sand is the matrix

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

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

RTENOTITLE....................(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 RTENOTITLE is more closely aligned with the spatial orientation of Region I than coordinate system RTENOTITLE, while coordinate system RTENOTITLE is more closely aligned with the spatial orientation of Region II than coordinate system RTENOTITLE. The coordinate system RTENOTITLE is obtained by rotating the coordinate system RTENOTITLE through an angle θ as in Fig. 1.

We consider four grid orientations for each case: (1) grid RTENOTITLE in Regions I and II; (2) grid RTENOTITLE in Region I and grid RTENOTITLE in Region II; (3) grid RTENOTITLE in Region I and grid RTENOTITLE in Region II; and (4) grid RTENOTITLE 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

RTENOTITLE = rotation matrix
RTENOTITLE = transpose of matrix A , Eq. 54
A = cross-sectional area, Eq. 1
RTENOTITLE = 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
RTENOTITLE = 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
RTENOTITLE = flow rate, Eqs. 1 and 5
RTENOTITLE = column vectors, Eq. 41
Δt = time interval
Δx = length
Δy = width
Δz = thickness
θ = angle
λ = eigenvalues
μ = fluid viscosity, Eq. 1
RTENOTITLE = 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