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


PEH:Mathematics of Fluid Flow

PetroWiki
Revision as of 16:33, 26 April 2017 by Denise Watts (Denisewatts) (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

Publication Information

Vol1GECover.png

Petroleum Engineering Handbook

Larry W. Lake, Editor-in-Chief

Volume I – General Engineering

John R. Fanchi, Editor

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.[1][2][3][4]

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.[5] 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 (Jx)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 Jx and out of the block at x + Δx with fluid flux Jx + Δx. Applying the principle of conservation of mass, we have the mass balance, which is written as

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

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

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

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

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


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

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

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

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

RTENOTITLE....................(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 (vx), and no source or sink terms. The concentration, C, and flux, Jx, for this example are

RTENOTITLE

and

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

Substituting Eq. 2.8 into Eq. 2.6 gives

RTENOTITLE....................(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[6]

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

where P is pressure, cf 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

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

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

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

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


Boundary conditions for second-order PDEs may be written as

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

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

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

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

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

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

Eq. 2.17 is the 1D convection/dispersion (C/D) equation. The dispersion term is D∂2C/∂x2, 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

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

RTENOTITLE....................(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 (vx). 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&percnt; 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

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

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

RTENOTITLE....................(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.[7] Eq. 2.20 is illustrated in Fig. 2.2 for physical parameters v = 1 ft/day and D = 0.01ft2/day. This solution is used in Sec. 2.4 to evaluate a numerical solution of the 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

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

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

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

where the divergence of vector RTENOTITLE = {Jx, Jy, Jz}, in Cartesian coordinates, is

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

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

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

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

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

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

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

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

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

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


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 RTENOTITLE. The result is the divergence of the vector field.

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

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

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

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

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

where we introduce the Laplacian operator,  

RTENOTITLE....................(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 RTENOTITLE is said to be irrotational if curl RTENOTITLE = 0, and it is said to be solenoidal if div RTENOTITLE = 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 RTENOTITLE has concentration (C) and flux (J→) given by

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

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

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

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

into Eq. 2.35 to obtain

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

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

RTENOTITLE....................(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 RTENOTITLE in the multidimensional form

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

Substituting Eq. 2.39 into the 3D continuity equation gives

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

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

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

Eq. 2.41 is the 3D convection/dispersion equation. The term D2C is the dispersion term, and the term RTENOTITLE 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.

RTENOTITLE....................(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.[8][9][10][11][12][13][14]

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 = xi and x = xi + Δx, where Δx is an increment along the x-axis (Fig. 2.4). We can approximate the derivative df(x) / dx at x = xi by solving the Taylor’s series,

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

for df(x) / dx. The result is

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

where ET is the term

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

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

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

Eq. 2.46 is an approximation because we are neglecting all of the terms in ET, 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.


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,

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

and the centered difference,

RTENOTITLE....................(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.[8] 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

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

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

The future concentration distribution is found from the current concentration distribution by rearranging Eq. 2.49. We collect terms in Cn+1 on the left-hand side and terms in Cn on the right-hand side, thus

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

Eq. 2.50 is now written in the form

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

where the coefficients are

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

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

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

Eq. 2.53 can be written in matrix form as

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

where RTENOTITLE  is the NX × NX matrix of coefficients, RTENOTITLE is the column vector of unknown concentrations at time tn+1, and RTENOTITLE is the column vector of right-hand-side terms that depend on known concentrations at time tn. Both column vectors RTENOTITLE and RTENOTITLE 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.[8][9][10][11][15] A solution of the set of equations for physical parameters v = 1 ft/day and D = 0.01 ft2/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,[8][16][17] 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.

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.[4] 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 {x1, x2} to another set {y1, y2} by the angle θ. The equations relating the coordinate systems are

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

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

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

which can be written as

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

The column vectors RTENOTITLE and RTENOTITLE are

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

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

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



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 ith row and jth column of the matrix is the ijth 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 mn, 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 RTENOTITLE in Eq. 2.58 may be written as  

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

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

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

The transpose of matrix RTENOTITLE is

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

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

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

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

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

The set of elements {aii} of a square matrix RTENOTITLE is the principal diagonal of the matrix. The elements {aji} with ij are off-diagonal elements. The matrix RTENOTITLE is a lower triangular matrix if aij = 0 for i < j, and RTENOTITLE is an upper triangular matrix if aij = 0 for i > j. The matrix RTENOTITLE is a diagonal matrix if aij =0 for ij.

Matrix Operations

Suppose the matrices RTENOTITLE, RTENOTITLE, and RTENOTITLE with elements {aij}, {bij}, and {cij} 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

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

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

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

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

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

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

The transpose of the product of two square matrices RTENOTITLE and RTENOTITLE is

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

and the adjoint of the product of two square matrices is

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

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

The identity matrix, RTENOTITLE, 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 RTENOTITLE in matrix multiplication, thus

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Determinants, Eigenvalues, and Eigenvectors

The determinant (det) of a square matrix RTENOTITLE is denoted by det RTENOTITLE or | RTENOTITLE |. 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

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

and the determinant of a 3×3 matrix is

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

Determinants are useful for determining if an inverse matrix RTENOTITLE-1 exists. Inverse matrices are needed to solve finite-difference equations representing fluid flow. The condition det RTENOTITLE  says that an inverse matrix RTENOTITLE-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

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

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

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

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

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

Eq. 2.80 is the characteristic equation of RTENOTITLE, 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,[18] geomechanics,[19] and upscaling[20] 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

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

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

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

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

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

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

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

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

Similarity Transformations

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

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

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

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

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

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

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

Defining the transformed matrices

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

and using the similarity transformation

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

in Eq. 2.89 yields

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

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

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

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

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

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

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

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

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

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

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

or

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

We expand the characteristic equation to get

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

The two eigenvalues are found from the quadratic equation to be

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

The sum of the eigenvalues satisfies the relation

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

Eigenvectors

The matrix RTENOTITLE is composed of orthonormal eigenvectors RTENOTITLE found from

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

The basis vector, RTENOTITLE, satisfies

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

with the identity matrix RTENOTITLE. Expanding Eq. 2.105 gives

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

for the eigenvalue λ+, and

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

for the eigenvalue λ. Rearranging Eq. 2.106 gives

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

Eq. 2.108 and the normalization condition,

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

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

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

and

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

Similar calculations for RTENOTITLE- yield the results

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

and

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

where the relation

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

from Eq. 2.103, has been used.

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

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

Substituting in the expressions for RTENOTITLE+ and RTENOTITLE- gives

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

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

or

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

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

or

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

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

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

The coordinate systems RTENOTITLE and RTENOTITLE 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

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

and

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

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

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

or

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

Note that the equality,

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

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

The coordinate systems RTENOTITLE and RTENOTITLE are related by the counterclockwise rotation shown in Fig. 2.6. 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

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

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

where RTENOTITLE is

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

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

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

We find the elements of RTENOTITLE-1 by solving

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

The result is

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

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

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

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

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


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

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

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


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,[16] 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.

Nomenclature


ai, bi, ci, di = finite-difference coefficients, Eq. 2.51
aij, bij, cij = elements of matrices, Eq. 2.64
RTENOTITLE = matrices, Eq. 2.64
RTENOTITLE = rotation matrix, Eq. 2.57
RTENOTITLE = 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
RTENOTITLE = 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
RTENOTITLE = 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
RTENOTITLE = unit vectors in Cartesian coordinates, Eq. 2.25
RTENOTITLE = identity matrix, Eq. 2.69
Jx, Jy, Jz = fluid flux in x-, y-, z-directions
RTENOTITLE = 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
RTENOTITLE = 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
RTENOTITLE = matrix of coefficients, Eq. 2.54
P = pressure, Eq. 2.82
RTENOTITLE = 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
RTENOTITLE = vector field, Eq. 2.27
v = velocity of solute, Eq. 2.15
vx = speed in x-direction, Eq. 2.9
RTENOTITLE = position vector, Eq. 2.25
x,y,z = space dimensions
xi = discrete point in x-direction, Eq. 2.43
RTENOTITLE = 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
RTENOTITLE = 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


  1. 1.0 1.1 Kreyszig, E. 1999. Advanced Engineering Mathematics, eighth edition. New York: John Wiley & Sons.
  2. Collins, R.E. 1999. Mathematical Methods for Physicists and Engineers, second revised edition. New York: Dover Publications.
  3. Chow, T.L. 2000. Mathematical Methods for Physicists: A Concise Introduction. Cambridge, UK: Cambridge University Press.
  4. 4.0 4.1 4.2 4.3 Fanchi, J.R. 2006. Math Refresher for Scientists and Engineers, third edition. New York: Wiley Interscience.
  5. Fanchi, J.R. 2002. Shared Earth Modeling. Boston, Massachusetts: Butterworth-Heinemann.
  6. Towler, B.F. 2002. Fundamental Principles of Reservoir Engineering, Vol. 8. Richardson, Texas: Textbook Series, SPE.
  7. Abramowitz, M. and Stegun, I.A. ed. 1972. Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables, ninth edition. New York: Dover Publications.
  8. 8.0 8.1 8.2 8.3 Peaceman, D.W. 1977. Fundamentals of Numerical Reservoir Simulation. Oxford, UK: Elsevier Publishing.
  9. 9.0 9.1 Aziz, K. and Settari, A. 1979. Petroleum Reservoir Simulation. Essex, UK: Elsevier Applied Science Publishers.
  10. 10.0 10.1 Mattax, C.C. and Dalton, R.L. 1990. Reservoir Simulation, Vol. 13. Richardson, Texas: Monograph Series, SPE.
  11. 11.0 11.1 Ertekin, T., Abou-Kassem, J.H., and King, G.R. 2001. Basic Applied Reservoir Simulation, Vol. 7. Richardson, Texas: Textbook Series, SPE.
  12. Munka, M. and Pápay, J. 2001. 4D Numerical Modeling of Petroleum Reservoir Recovery. Budapest, Hungary: Akadémiai Kiadó.
  13. Fanchi, J.R. 2006. Principles of Applied Reservoir Simulation, third edition. Burlington, Massachusetts: Gulf Professional Publishing/Elsevier.
  14. Fanchi, J.R. 2000. Integrated Flow Modeling, No. 49. Amsterdam, The Netherlands: Developments in Petroleum Science, Elsevier Science B.V.
  15. Chapra, S.C. and Canale, R.P. 2002. Numerical Methods for Engineers, fourth edition. Boston, Massachusetts: McGraw-Hill Book Co.
  16. 16.0 16.1 Fanchi, J.R. 1983. Multidimentional Numerical Dispersion. SPE J. 23 (1): 143–151. SPE-9018-PA. http://dx.doi.org/10.2118/9018-PA
  17. Lantz, R.B. 1971. Quantitative Evaluation of Numerical Diffusion. SPE J. 11 (3): 315–320. SPE-2811-PA. http://dx.doi.org/10.2118/2811-PA
  18. 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
  19. 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
  20. 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

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.