Diagonalizing the permeability tensor

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, geomechanics, and upscaling 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

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

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

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

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

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

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

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

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

The diagonal permeability elements {k_xx, k_yy, k_zz} represent the dependence of flow rate in one direction on pressure differences in the same direction. The off-diagonal permeability elements {k_xy, k_xz, k_yx, k_yz, k_zx, k_zy} 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.

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

Similarity transformations

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

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

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

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

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

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

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

Defining the transformed matrices

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

and using the similarity transformation

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

in Eq. 9 yields

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

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

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

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

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

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

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

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

• 

  Fig. 1 – Rotation of the coordinate system.

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

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

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

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

where λ₊ and λ_– are calculated as the solutions to the eigenvalue problem.

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

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

or

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

We expand the characteristic equation to get

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

The two eigenvalues are found from the quadratic equation to be

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

The sum of the eigenvalues satisfies the relation

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

Eigenvectors

The matrix is composed of orthonormal eigenvectors found from

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

The basis vector, , satisfies

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

with the identity matrix . Expanding Eq. 25 gives

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

for the eigenvalue λ₊, and

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

for the eigenvalue λ_–. Rearranging Eq. 26 gives

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

Eq. 28 and the normalization condition,

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

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

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

and

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

Similar calculations for ^- yield the results

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

and

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

where the relation

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

from Eq. 33, has been used.

To show that ⁺ and ^- are orthogonal, we must show that

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

Substituting in the expressions for ⁺ and ^- gives

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

as expected.

Coordinate transformation

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

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

or

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

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

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

or

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

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

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

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

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

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

and

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

The equalities a₁⁺ = a₂^– and a₂⁺ = –a₁^– are in agreement with Eqs. 30 through 33. The angle θ may be found from either

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

or

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

Note that the equality,

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

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

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

The coordinate systems and are related by the counterclockwise rotation shown in Fig. 1. We have an aligned coordinate system with the principal axes of the permeability tensor. The diagonal tensor in the coordinate system has the form

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

where k_max is the maximum permeability in the direction y₁ and k_T is the permeability that is transverse to k max in the direction y₂. We want to know how the elements of the permeability tensor change if we transform to the different coordinate system .

The relationship between the elements of and is

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

where is

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

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

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

We find the elements of ^-1 by solving

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

The result is

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

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

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

We can use Eq. 55 to calculate the elements of for any rotation angle θ. If the permeability is isotropic, we have k_max = k_T = k_iso and Eq. 55 simplifies to the form

....................(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 k_max = 200 md and k_T = 50 md. The values of elements k₁₁, k₁₂ and k₂₂ of are presented in the figure. The off-diagonal terms satisfy the equality k₁₂ = k₂₁ for a symmetric matrix given in Eq. 16, so it is sufficient to show only k₁₂. 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 x₂, and k_T is aligned along x₁.

• 

  Fig. 2 – Coordinate rotation in 2D.

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.

• 

  Fig. 3 – Plan view of channel sand with two permeability regions.

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 k_T that is transverse to the direction of k_max. The diagonalized, anisotropic permeability tensor ′ in the y₁ - y₂ plane of a channel sand is the matrix

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

and Darcy’s law for flow in the y₁ - y₂ plane is

....................(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. 4 – Case A: homogeneous and anisotropic permeability.

• 

  Fig. 5 – Case B: heterogeneous and anisotropic permeability.

Fig. 6 shows two possible coordinate systems for orienting the grid. Coordinate system is more closely aligned with the spatial orientation of Region I than coordinate system , while coordinate system is more closely aligned with the spatial orientation of Region II than coordinate system . The coordinate system is obtained by rotating the coordinate system through an angle θ as in Fig. 1.

• 

  Fig. 6 – Possible coordinate systems for grid orientation.

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

• 

  Table 1

Nomenclature

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

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