PEH:Mathematics of Fluid Flow

Publication Information

Petroleum Engineering Handbook

Larry W. Lake, Editor-in-Chief

Volume I – General Engineering

John R. Fanchi, Editor

Copyright 2007, Society of Petroleum Engineers

Chapter 2 – Mathematics of Fluid Flow

John R. Fanchi, Colorado School of Mines

Pgs. 45-76

ISBN 978-1-55563-108-6
Get permission for reuse

The purpose of this chapter is to review the mathematics of fluid flow. We limit our review to essential aspects of partial differential equations, vector analysis, numerical methods, matrices, and linear algebra. These topics are discussed in the context of two fluid flow applications: analysis of the convection/dispersion equation and diagonalization of the permeability tensor. For more details about the mathematics presented here, consult the literature.

Partial Differential Equations

Partial differential equations (PDEs) are frequently encountered in petroleum engineering. We review basic concepts of PDEs by considering the relevant mathematical properties of the continuity equation.

Continuity Equation

Fluid flow through a volume can be described mathematically by the continuity equation. The continuity equation has many uses, and its derivation is provided to illustrate the construction of a partial differential equation from physical reasoning. We begin by considering the flow illustrated in Fig. 2.1. The block in Fig. 2.1 has length (Δx), width (Δy), and depth (Δz). Fluid flux (J) is the rate of flow of mass per unit cross-sectional area normal to the direction of flow. The notation (J_x)_x denotes fluid flux in the x direction at location x. The cross-sectional area perpendicular to the flux direction is Δy Δz. Fluid flows into the block at x with fluid flux J_x and out of the block at x + Δx with fluid flux J_x + Δx. Applying the principle of conservation of mass, we have the mass balance, which is written as

....................(2.1)

The mass entering the block in time interval, Δt, for flux across the face of the block at x is

....................(2.2)

The mass leaving the block in time interval, Δt, across the face of the block at x + Δx is

....................(2.3)

where we have added a source term q. Flow out of the block through q is represented by q > 0, and flow into the block is represented by q < 0. The source term, q, can represent a variety of important physical systems, including wells, aquifer support, fluid flow into a fracture system from matrix blocks in a naturally fractured reservoir, and gas flow into a cleat system from the coal in a coalbed.

• 

  Fig. 2.1 – Coordinate system for continuity equation.

Mass accumulation in the block is the change in concentration C of the mass in the block during the time interval Δt, where concentration, C, is defined as the total mass in the block divided by the block volume. The mass accumulation term is

....................(2.4)

where concentration is evaluated at times t and t + Δt.

Substituting Eqs. 2.2 through 2.4 into Eq. 2.1, dividing by Δx Δy Δz Δt, and rearranging gives

....................(2.5)

In the limit as Δx → 0, Δt → 0, the differences in Eq. 2.5 are replaced by partial derivatives. We assume the fluxes and concentrations are sufficiently smooth and continuous to allow the replacement of differences by partial derivatives. Eq. 2.5 becomes the continuity equation in one space dimension.

....................(2.6)

Eq. 2.6 is an example of a partial differential equation.

Partial Differential Equations

PDEs are an extension of the concept of ordinary differential equations (ODEs). Unlike an ODE, which depends on only one independent variable, a PDE depends on two or more independent variables. In the previous example, Eq. 2.6 depends on two independent variables: one space dimension (x) and time (t). The order of Eq. 2.6 depends on the form of concentration and flux. The order of a PDE is the order of the highest derivative that appears in the equation.

....................(2.7)

for a function (ψ) of two or more independent variables {x,y, ...}. A PDE is linear if it is first order in the unknown function and its partial derivatives, and the coefficients of the partial derivatives, are either constant or depend on the independent variables {x,y, ...}. We illustrate these concepts by considering the continuity equation for flow of a fluid with density (ρ), velocity (v_x), and no source or sink terms. The concentration, C, and flux, J_x, for this example are



and

....................(2.8)

Substituting Eq. 2.8 into Eq. 2.6 gives

....................(2.9)

Eq. 2.9 is a linear, first-order PDE if density is the unknown function and velocity is constant. The situation is not so simple in more physically realistic systems.

Consider, for example, a slightly compressible fluid in which density is given by

....................(2.10)

where P is pressure, c_f is fluid compressibility, and the subscript, 0, refers to a reference value of pressure. Assume, as well, that velocity is proportional to pressure gradient so that

....................(2.11)

where α is the proportionality constant. Substituting Eqs. 2.10 and 2.11 into Eq. 2.9 gives

....................(2.12)

Eq. 2.12 is a nonlinear, second-order PDE. It is second order because of the second-order partial derivative of pressure with respect to x, and it is nonlinear because of the square of the pressure gradient term.

Solutions of PDEs depend on the form of the PDEs and their associated boundary conditions. An important class of second-order PDEs has the form

....................(2.13)

where the functions {A, B, C, G} are known functions of two independent variables {x,y} and the first-order partial derivatives ∂ ψ(x, y) / ∂  x and ∂ ψ(x, y) / ∂  y in a region (R) bounded by a surface (S). The mathematical properties of the second-order PDEs depend on the relationship between the functions {A, B, C, G}. A classification scheme for second-order PDEs is given in Table 2.1.

• 

  Table 2.1

Boundary conditions for second-order PDEs may be written as

....................(2.14)

where ψ(x, y) is the unknown function of two independent variables {x, y}, and ∂ ψ(x, y) / ∂  n is the derivative normal to a boundary. The functions {α, β, γ} are known functions of {x, y}. All of the functions and applicable derivatives are defined in a domain (R) bounded by a surface (S). A classification scheme for the boundary conditions of a second-order PDE is given in Table 2.2. The boundary conditions associated with the examples in Table 2.1 are given in Table 2.3. The significance of PDE classification is considered further in the discussion of the convection/dispersion equation presented next.

• 

  Table 2.2

• 

  Table 2.3

One-Dimensional (1D) Convection/Dispersion Equation

The continuity equation is used to describe the mixing of one substance with another by writing flux in the form

....................(2.15)

The concentration, C, is the concentration of the solute in the solvent. The term v is the velocity of the solute, and D is the dispersion of the solute into the solvent. Substituting Eq. 2.15 into Eq. 2.6, the 1D convection/dispersion equation without sources or sinks, gives

....................(2.16)

If we assume that v and D are constant, Eq. 2.16 simplifies to the form

....................(2.17)

Eq. 2.17 is the 1D convection/dispersion (C/D) equation. The dispersion term is D∂²C/∂x², and the convection term is –v∂C/∂x. If the dispersion term is much larger than the convection term, the solution of Eq. 2.17 can be approximated by the solution of the equation

....................(2.18)

Eq. 2.18 is a parabolic PDE and behaves mathematically like a heat conduction equation. If the convection term is much larger than the dispersion term, the solution of Eq. 2.17 can be approximated by the solution of the equation.

....................(2.19)

Eq. 2.19 is a first-order hyperbolic PDE.

A solution of the 1D C/D equation, given in Eq. 2.17, is obtained as follows. We assume that solute is moving in the x-direction with constant speed (v_x). The concentration C(x,t) is a function of space and time. We must specify two boundary conditions and an initial condition for the concentration C(x,t) to solve the C/D equation. We impose the boundary conditions, C(0, t) = 1 and C(∞, t) = 0, for all time, t > 0, and the initial condition, C (x, 0) = 0, for all values of x > 0. The boundary condition, C (0, t) = 1, says that we are injecting 100% solute at x = 0, and the boundary condition, C(∞, t) = 0, says that the solute never reaches the end of the flow path at x = ∞. The initial condition, C(x, 0) = 0, says that there is no solute in the solvent at the initial time, t = 0. The solution of the C/D equation is

....................(2.20)

where the complementary error function erfc(y) is defined as

....................(2.21)

The integral in Eq. 2.21 can be solved using the series expansion on the right side of Eq. 2.21 or a numerical algorithm. Eq. 2.20 is illustrated in Fig. 2.2 for physical parameters v = 1 ft/day and D = 0.01ft²/day. This solution is used in Sec. 2.4 to evaluate a numerical solution of the C/D equation.

• 

  Fig. 2.2 – Analytical solution of the 1D C/D equation.

Vector Analysis

Fluid flow equations in two and three dimensions can be compactly represented using concepts from vector analysis. For example, the continuity equation in three space dimensions for the Cartesian coordinate system, shown in Fig. 2.1, is

....................(2.22)

The flux terms (J_y) and (J_z) have meanings analogous to (J_x) for flux in the y and z directions, respectively. If we write the components of flux as the flux vector = {J_x, J_y, J_z}, Eq. 2.22 can be written in vector notation as

....................(2.23)

where the divergence of vector = {J_x, J_y, J_z}, in Cartesian coordinates, is

....................(2.24)

The divergence operator ∇• is an example of an operator from vector analysis that determines the spatial variation of a vector or scalar field. Following Fanchi, we first review the concepts of scalar and vector fields and then define gradient (grad), divergence (div), and curl operators.

Scalar and Vector Fields

We define scalar and vector fields in a Cartesian coordinate system with position vector

....................(2.25)

where are unit vectors defined along the orthogonal {x,y,z} coordinate axes. If we can associate a scalar function (f) with every point in a region (R), then the scalar field may be written as

....................(2.26)

Examples of scalar fields include pressure, temperature, and saturation.

If, instead of a scalar function (f), we can associate a vector with every point in the region (R), we can construct a vector field of the form

....................(2.27)

The vector field is a function that assigns a vector to every point in the region R. Examples of vector fields include the Darcy velocity field and seismic velocities.

Gradient, Divergence, and Curl

The spatial variation of a scalar or vector field can be determined with the del operator ∇. The del operator, ∇, is defined in Cartesian coordinates as

....................(2.28)

The gradient of a scalar field (f) is obtained by operating on the scalar field with the del operator, thus

....................(2.29)

The direction of the gradient of the scalar field (f) evaluated at a point is oriented in the direction of maximum increase of the scalar field. In addition, the vector field, ∇f, is perpendicular to a surface that corresponds to a constant value of the scalar field (Fig. 2.3).

• 

  Fig. 2.3 - Gradient ∇f normal to surface S.

Two outcomes are possible when the del operator is applied to a vector field. One outcome is to create a scalar, and the other is to create a vector. A scalar is obtained when we take the dot product of the del operator with a vector field . The result is the divergence of the vector field.

....................(2.30)

A vector is obtained when we take the cross product of the del operator with a vector field . The result is the curl of the vector field .

....................(2.31)

The curl of the vector field is called the rotation of the vector field. It is a vector that is normal to the plane containing the vector field . The divergence of the gradient of a scalar field ( f ) is

....................(2.32)

where we introduce the Laplacian operator,  

....................(2.33)

in Cartesian coordinates.

The gradient, divergence, curl, and Laplacian operators arise in many PDEs that affect petroleum engineering. For example, a vector field is said to be irrotational if curl  = 0, and it is said to be solenoidal if div  = 0. These properties of the vector field are useful for analyzing the propagation of seismic waves. Another useful application of vector analysis is to the mathematical representation of fluid flow in two or three spatial dimensions. Two examples are presented next.

Incompressible Flow

Incompressible flow occurs when the density of a fluid is constant. In this case, the continuity equation for flow of a fluid with density (ρ) and velocity has concentration (C) and flux (J→) given by

....................(2.34)

The concentration and density are scalar fields, and the velocity and flux are vector fields. The continuity equation without source or sink terms becomes

....................(2.35)

A more suitable form of the continuity equation for describing incompressible fluid flow is obtained by substituting the differential operator,

....................(2.36)

into Eq. 2.35 to obtain

....................(2.37)

In the case of incompressible fluid flow, density is constant and Eq. 2.37 reduces to

....................(2.38)

Eq. 2.38 shows that the divergence of the velocity of a flowing, incompressible fluid is zero.

Three-Dimensional (3D) Convection/Dispersion Equation

The convection/dispersion equation in three dimensions is obtained by writing flux in the multidimensional form

....................(2.39)

Substituting Eq. 2.39 into the 3D continuity equation gives

....................(2.40)

If we assume that and D are constant, we can simplify Eq. 2.40 to the form of

....................(2.41)

Eq. 2.41 is the 3D convection/dispersion equation. The term D∇²C is the dispersion term, and the term is the convection term.

Numerical Methods

Systems of nonlinear PDEs are needed to describe realistic multiphase, multidimensional flow in a reservoir. As a rule, these equations cannot be solved analytically; they must be solved with numerical methods. To illustrate the mathematics, we discuss the numerical solution of the 1D C/D equation.

....................(2.42)

as introduced in Sec. 2.2. As a reminder, v is velocity, D is dispersion, and C is concentration. Eq. 2.42 is a good example to use because it illustrates many useful numerical methods that can be compared with the analytical solution given by Eq. 2.20. We first introduce the concept of finite differences to convert Eq. 2.42 to an equation that can be solved numerically. We then present a numerical representation of Eq. 2.42 and illustrate its solution. For more details, you should consult the chapter on reservoir simulation in Vol. V, Reservoir Engineering and Petrophysics, as well as several other sources.

Finite Differences

One way to solve a PDE is to convert the PDE to finite-difference form. The finite-difference form is obtained by replacing the derivatives in the PDE with differences that are obtained from Taylor’s series. To illustrate the procedure, let us suppose that we know the function f(x) at two discrete points x = x_i and x = x_i + Δx, where Δx is an increment along the x-axis (Fig. 2.4). We can approximate the derivative df(x) / dx at x = x_i by solving the Taylor’s series,

....................(2.43)

for df(x) / dx. The result is

....................(2.44)

where E_T is the term

....................(2.45)

If we neglect E_T, we obtain the finite-difference approximation of the first derivative.

....................(2.46)

Eq. 2.46 is an approximation because we are neglecting all of the terms in E_T, which is called the truncation error. In the limit as the increment Δx approaches zero, the truncation error approaches zero, and the finite difference approaches the definition of the derivative.

• 

  Fig. 2.4 – Discrete points in the Taylor series expansion.

The finite difference in Eq. 2.46 is called a forward difference. Other differences are possible. Two that we use next are the backward difference,

....................(2.47)

and the centered difference,

....................(2.48)

Eqs. 2.46 through 2.48 are all derived from Taylor’s series.

Illustration: Numerical Solution of the 1D C/D Equation

We illustrate the application of finite differences in a fluid flow problem by considering a specific finite-difference representation of the 1D C/D equation. For a more detailed discussion of the numerical analysis of Eq. 2.42, see Chap. 4 of Peacemen. In our example, we choose a backward difference for the time derivative in Eq. 2.42, a centered difference for the space derivative in the convection term, and a centered-in-time/centered-in-space difference for the dispersion term. Eq. 2.42 is converted from a PDE to the difference equation

....................(2.49)

The subscripts of concentration C denote points in space, and the superscripts denote points in time. For example, the present time, tⁿ, is denoted by superscript n and future time tⁿ⁺¹ is denoted by n+1. The time increment is Δt = tⁿ⁺¹ - tⁿ. Similarly, the space increment is Δx = xⁱ + 1 - xⁱ. The concentration at time tⁿ⁺¹ and spatial location xi is denoted by .

The future concentration distribution is found from the current concentration distribution by rearranging Eq. 2.49. We collect terms in Cⁿ⁺¹ on the left-hand side and terms in Cⁿ on the right-hand side, thus

....................(2.50)

Eq. 2.50 is now written in the form

....................(2.51)

where the coefficients are

....................(2.52)

All values of the variables in the coefficients are known at time tⁿ. If we assume that the spatial subscript is in the range 1 ≤ I ≤NX, the system of finite-difference equations becomes

....................(2.53)

Eq. 2.53 can be written in matrix form as

....................(2.54)

where   is the NX × NX matrix of coefficients, is the column vector of unknown concentrations at time tⁿ⁺¹, and is the column vector of right-hand-side terms that depend on known concentrations at time tn. Both column vectors and have NX elements.

The system of equations in Eq. 2.54 is called a tridiagonal system because it consists of three lines of nonzero diagonal elements centered about the main diagonal. All other elements are zero. Techniques for solving the tridiagonal system of equations, using the Thomas algorithm, are presented in several sources. A solution of the set of equations for physical parameters v = 1 ft/day and D = 0.01 ft²/day and finite-difference parameters Δx = 0.1 ft and Δt = 0.1 day is shown in Fig. 2.5. The difference between the analytical solution and the numerical solution is because of numerical dispersion, which is beyond the scope of this chapter. What interests us here is the appearance of matrices in the mathematics of fluid flow. Matrices are the subject of the next section.

• 

  Fig. 2.5 – Numerical solution of the 1D C/D equation.

Matrices and Linear Algebra

An example of a matrix was introduced in Sec. 2.4 for the 1D C/D equation. It is often easier to work with many fluid flow equations when they are expressed in terms of matrices. Our review follows the presentation in Fanchi. We begin our discussion with an example of a matrix that is used later in this chapter, namely the matrix associated with the rotation of a coordinate system. We then summarize some important properties of matrices and determinants and review the concepts of eigenvalues and eigenvectors from linear algebra.

Rotation of a Cartesian Coordinate System

Fig. 2.6 illustrates a rotation of Cartesian coordinates from one set of orthogonal coordinates {x₁, x₂} to another set {y₁, y₂} by the angle θ. The equations relating the coordinate systems are

....................(2.55)

The set of equations in Eq. 2.55 has the matrix form

....................(2.56)

which can be written as

....................(2.57)

The column vectors and are

....................(2.58)

with two elements each, and the rotation matrix is the 2 × 2 square matrix,

....................(2.59)

• 

  Fig. 2.6 – Rotation of the coordinate system.

Properties of Matrices

In general, a matrix with m rows and n columns has the order m × n and is referred to as a m × n matrix. The entry in the i^th row and j^th column of the matrix is the ij^th element of the matrix. If the number of rows equals the number of columns so that m = n, the matrix is called a square matrix. On the other hand, if m ≠ n, the matrix is a rectangular matrix.

If the matrix has a single column so that n = 1, it is a column vector as in Eq. 2.58. If the matrix has a single row so that m = 1, it is a row vector. A row vector can be created from a column vector by taking the transpose of the column vector. For example, the transpose of the column vector in Eq. 2.58 may be written as  

....................(2.60)

where the superscript T denotes the transpose of the matrix. In general, we can write a m × n matrix with a set of elements {a_ij: i = 1, 2, ... n; j = 1, 2, ... m} as

....................(2.61)

The transpose of matrix is

....................(2.62)

The conjugate transpose of matrix is obtained by finding the complex conjugate of each element in and then taking the transpose of the matrix . This operation can be written as

....................(2.63)

where * denotes complex conjugation. Recall that the conjugate z* of a complex number z is obtained by replacing the imaginary number with wherever it occurs. If all the elements of matrix are real, the conjugate transpose of matrix is equal to the transpose of matrix .

If the matrix is a square matrix and the elements of matrix satisfy the equality a_ij = a_ji, the matrix is called a symmetric matrix. A square matrix is Hermitian, or self-adjoint, if = ⁺ (i.e, the matrix equals its conjugate transpose).

The set of elements {a_ii} of a square matrix is the principal diagonal of the matrix. The elements {a_ji} with i ≠ j are off-diagonal elements. The matrix is a lower triangular matrix if a_ij = 0 for i < j, and is an upper triangular matrix if a_ij = 0 for i > j. The matrix is a diagonal matrix if a_ij =0 for i ≠ j.

Matrix Operations

Suppose the matrices , , and with elements {a_ij}, {b_ij}, and {c_ij} have the same order m × n. We are using double underlines to denote matrices. Other notations are often used, such as boldface. The addition or subtraction of two matrices may be written as

....................(2.64)

The product of a matrix with a number k may be written as

....................(2.65)

The product of matrix with order m × n and matrix with order n × p is

....................(2.66)

where matrix has order m × p. Notice that matrix multiplication is possible only if the number of columns in equals the number of rows in . This requirement is always satisfied for square matrices.

The transpose of the product of two square matrices and is

....................(2.67)

and the adjoint of the product of two square matrices is

....................(2.68)

Notice that the product of two matrices may not be commutative (i.e., ≠ in general).

The identity matrix, , is a square matrix with all off-diagonal elements equaling zero and all diagonal elements equaling one. The identity matrix preserves the identity of a square matrix in matrix multiplication, thus

....................(2.69)

By contrast, a null matrix is a matrix in which all elements are zero. In this case, the product of the null matrix with a matrix is

....................(2.70)

The matrix, , is singular if the product of matrix with a column vector that has at least one nonzero element yields the null matrix; that is, is singular if

....................(2.71)

The concepts of identity matrix and matrix singularity are needed to define the inverse matrix. Suppose we have two square matrices and that satisfy the product

....................(2.72)

Notice that the matrices and commute. The matrix is nonsingular, and the matrix is the inverse of , thus = ^-1, where ^-1 denotes the inverse of . Eq. 2.72 can be written in terms of the inverse as

....................(2.73)

The inverse matrix is useful for solving systems of equations. For example, suppose we have a system of equations that satisfies

....................(2.74)

where the column vector and the matrix are known, and the column vector contains a set of unknowns. Eq. 2.54 is an example for the 1D C/D equation. We can solve for in Eq. 2.74 by premultiplying Eq. 2.74 by ^-1. The result is

....................(2.75)

Of course, we have to know ^-1 to find . This leads us to a discussion of determinants.

Determinants, Eigenvalues, and Eigenvectors

The determinant (det) of a square matrix is denoted by det or | |. Two examples of determinants are the determinants of a 2×2 matrix and a 3×3 matrix. The determinant of a 2×2 matrix is

....................(2.76)

and the determinant of a 3×3 matrix is

....................(2.77)

Determinants are useful for determining if an inverse matrix ^-1 exists. Inverse matrices are needed to solve finite-difference equations representing fluid flow. The condition det   says that an inverse matrix ^-1 exists, even though we may not know the elements of the inverse matrix. Determinants are also useful for determining eigenvalues and eigenvectors.

Eigenvalues and eigenvectors are useful for understanding the behavior of physical quantities that may be represented by a matrix. An example in fluid flow is permeability, which we discuss in more detail later in this chapter. First, we need to define the concepts of eigenvalue and eigenvector.

Eigenvalues are the values of λ in the eigenvalue equation

....................(2.78)

where is an n × n square matrix and is a column vector with n rows. The eigenvalue equation may be written as

....................(2.79)

where is the n × n identity matrix. Eq. 2.79 has nonzero solutions, , if the eigenvalue, λ, is a characteristic root of , that is, λ must be a solution of

....................(2.80)

Eq. 2.80 is the characteristic equation of , and the n values of λ are the characteristic roots of the characteristic equation. The characteristic roots, λ, are obtained by expanding the determinant in Eq. 2.80 into an nth-degree polynomial and then solving for the n values of λ. These concepts are illustrated in the next section.

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 Sec. 2.5. The relationship between the diagonalized-permeability-tensor assumption and grid orientation is discussed in Sec. 2.7. 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

....................(2.81)

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.82)

In Eq. 2.82, Δ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

....................(2.83)

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

....................(2.84)

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

....................(2.85)

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. 2.83 gives the corresponding set of three equations demonstrating this dependence.

....................(2.86)

Similarity Transformations

Eq. 2.85 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. 2.85.

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

....................(2.87)

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

....................(2.88)

where is the n × n identity matrix. Substituting Eq. 2.88 into Eq. 2.87 gives

....................(2.89)

Defining the transformed matrices

....................(2.90)

and using the similarity transformation

....................(2.91)

in Eq. 2.89 yields

....................(2.92)

Eq. 2.92 is the same form as Eq. 2.85.

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

....................(2.93)

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,

....................(2.94)

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

....................(2.95)

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

....................(2.96)

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

....................(2.97)

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

....................(2.98)

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

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

....................(2.99)

or

....................(2.100)

We expand the characteristic equation to get

....................(2.101)

The two eigenvalues are found from the quadratic equation to be

....................(2.102)

The sum of the eigenvalues satisfies the relation

....................(2.103)

Eigenvectors

The matrix is composed of orthonormal eigenvectors found from

....................(2.104)

The basis vector, , satisfies

....................(2.105)

with the identity matrix . Expanding Eq. 2.105 gives

....................(2.106)

for the eigenvalue λ₊, and

....................(2.107)

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

....................(2.108)

Eq. 2.108 and the normalization condition,

....................(2.109)

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

....................(2.110)

and

....................(2.111)

Similar calculations for ^- yield the results

....................(2.112)

and

....................(2.113)

where the relation

....................(2.114)

from Eq. 2.103, has been used.

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

....................(2.115)

Substituting in the expressions for ⁺ and ^- gives

....................(2.116)

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

....................(2.117)

or

....................(2.118)

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

....................(2.119)

or

....................(2.120)

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

....................(2.121)

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

Equating elements of the transformation matrix in Eqs. 2.120 and 2.121 gives

....................(2.122)

and

....................(2.123)

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

....................(2.124)

or

....................(2.125)

Note that the equality,

....................(2.126)

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

....................(2.127)

The two coordinate systems are shown in Fig. 2.6.

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

....................(2.128)

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

....................(2.129)

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

....................(2.130)

where is

....................(2.131)

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

....................(2.132)

We find the elements of ^-1 by solving

....................(2.133)

The result is

....................(2.134)

where is the transpose of . Substituting Eqs. 2.131 and 2.134 into 2.132 gives

....................(2.135)

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

....................(2.136)

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.7 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. 2.96, 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.7 – 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. 2.8.

• 

  Fig. 2.8 – 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

....................(2.137)

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

....................(2.138)

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. 2.9 and 2.10.

• 

  Fig. 2.9 – Case A: homogeneous and anisotropic permeability.

• 

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

Fig. 2.11 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. 2.6.

• 

  Fig. 2.11 – 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 2.4.

An "ok" in the "Permeability Tensor" column in Table 2.4 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. 2.135. An "ok" in the "Formulation" column in Table 2.4 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 2.4, 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 2.4

Nomenclature

ai, bi, ci, di	=	finite-difference coefficients, Eq. 2.51
aij, bij, cij	=	elements of matrices, Eq. 2.64
	=	matrices, Eq. 2.64
	=	rotation matrix, Eq. 2.57
	=	transpose of matrix A , Eq. 2.134
A	=	cross-sectional area, Eq. 2.81
A, B, C, G	=	functions, Eq. 2.13
cf	=	fluid compressibility, Eq. 2.10
	=	column vector of unknown concentrations at tn +1, Eq. 2.54
C	=	concentration, Eq. 2.4
D	=	dispersion of the solute into solvent, Eq. 2.15
	=	column vector of terms that depend on known concentrations at tn +1, Eq. 2.54
ET	=	truncation error, Eq. 2.45
f	=	scalar function, Eq. 2.26
	=	unit vectors in Cartesian coordinates, Eq. 2.25
	=	identity matrix, Eq. 2.69
Jx, Jy, Jz	=	fluid flux in x-, y-, z-directions
	=	fluid flux vector, Eq. 2.23
(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. 2.85
k	=	permeability, Eq. 2.81
kiso	=	isotropic permeability, Eq. 2.136
kmax	=	maximum permeability, Eq. 2.129
kT	=	transverse permeability, Eq. 2.129
m,n	=	number of rows and columns, Sec. 2.5.2
	=	matrix of coefficients, Eq. 2.54
P	=	pressure, Eq. 2.82
	=	flow rate, Eqs. 2.81 and 2.85
q	=	source term, Eq. 2.3
R	=	region
S	=	surface
t	=	time
tn	=	present time
tn+1	=	future time
	=	vector field, Eq. 2.27
v	=	velocity of solute, Eq. 2.15
vx	=	speed in x-direction, Eq. 2.9
	=	position vector, Eq. 2.25
x,y,z	=	space dimensions
xi	=	discrete point in x-direction, Eq. 2.43
	=	column vectors, Eq. 2.121
α	=	proportionality constant, Eq. 2.11
{ α, β, γ }	=	functions, Eq. 2.14
Δt	=	time interval
Δx	=	length
Δy	=	width
Δz	=	thickness
θ	=	angle, Eq. 2.55
λ	=	eigenvalues, Eq. 2.78
μ	=	fluid viscosity, Eq. 2.81
ρ	=	density, Eq. 2.8
	=	pressure gradient, Eq. 2.85
Φ	=	phase potential, Eq. 2.81
ψ	=	function, Eq. 2.7

Subscripts

i	=	discrete x-direction index
i,j	=	matrix indices, Eq. 2.61
i,j,k	=	x-, y-, z-direction indices
NX	=	range of index, Eq. 2.53
t	=	time index, Eq. 2.4
x	=	x-direction index, Eq. 2.2
0	=	reference value of pressure, Eq. 2.10

Superscripts

*	=	complex conjugation
T	=	transpose of matrix

References

SI Metric Conversion Factor

bbl	×	1.589 873	E – 01	=	m3
cp	×	1.0*	E – 03	=	Pa•s
ft	×	3.048*	E – 01	=	m
ft2	×	9.290 304*	E – 02	=	m2
psi	×	6.894 757	E + 00	=	kPa

*Conversion factor is exact.