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


Numerical methods analysis of fluid flow: Difference between revisions

PetroWiki
Jump to navigation Jump to search
No edit summary
 
No edit summary
 
Line 1: Line 1:
Systems of nonlinear partial differential equations (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. This article provides an overview of these methods.
Systems of nonlinear partial differential equations (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. This article provides an overview of these methods.
 
== One-dimensional convection/dispersion equation ==


==One-dimensional convection/dispersion equation==
To illustrate the mathematics, we discuss the numerical solution of the 1D convection/dispersion (C/D) equation
To illustrate the mathematics, we discuss the numerical solution of the 1D convection/dispersion (C/D) equation


[[File:Vol1 page 0054 eq 004.png]]....................(1)
[[File:Vol1 page 0054 eq 004.png|RTENOTITLE]]....................(1)


as introduced in [[Vector analysis of fluid flow|vector analysis of fluid flow]]. As a reminder, ''v'' is velocity, ''D'' is dispersion, and ''C'' is concentration. '''Eq. 1''' 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'''.  
as introduced in [[Vector_analysis_of_fluid_flow|vector analysis of fluid flow]]. As a reminder, ''v'' is velocity, ''D'' is dispersion, and ''C'' is concentration. '''Eq. 1''' 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'''.


[[File:Vol1 page 0050 eq 001.png]]....................(2)
[[File:Vol1 page 0050 eq 001.png|RTENOTITLE]]....................(2)


We first introduce the concept of finite differences to convert '''Eq. 1''' to an equation that can be solved numerically. We then present a numerical representation of '''Eq. 1''' and illustrate its solution. For more details, you should consult the [[Reservoir simulation|reservoir simulation]] page, as well as sources in the literature.<ref name="r2" /><ref name="r3" /><ref name="r4" /><ref name="r5" /><ref name="r6" /><ref name="r7" /><ref name="r8" />
We first introduce the concept of finite differences to convert '''Eq. 1''' to an equation that can be solved numerically. We then present a numerical representation of '''Eq. 1''' and illustrate its solution. For more details, you should consult the [[Reservoir_simulation|reservoir simulation]] page, as well as sources in the literature.<ref name="r1">Peaceman, D.W. 1977. Fundamentals of Numerical Reservoir Simulation. Oxford, UK: Elsevier Publishing.</ref><ref name="r2">Aziz, K. and Settari, A. 1979. Petroleum Reservoir Simulation. Essex, UK: Elsevier Applied Science Publishers.</ref><ref name="r3">Mattax, C.C. and Dalton, R.L. 1990. Reservoir Simulation, Vol. 13. Richardson, Texas: Monograph Series, SPE.</ref><ref name="r4">Ertekin, T., Abou-Kassem, J.H., and King, G.R. 2001. Basic Applied Reservoir Simulation, Vol. 7. Richardson, Texas: Textbook Series, SPE.</ref><ref name="r5">Munka, M. and Pápay, J. 2001. 4D Numerical Modeling of Petroleum Reservoir Recovery. Budapest, Hungary: Akadémiai Kiadó.</ref><ref name="r6">Fanchi, J.R. 2006. Principles of Applied Reservoir Simulation, third edition. Burlington, Massachusetts: Gulf Professional Publishing/Elsevier.</ref><ref name="r7">Fanchi, J.R. 2000. Integrated Flow Modeling, No. 49. Amsterdam, The Netherlands: Developments in Petroleum Science, Elsevier Science B.V.</ref>
 
== Finite differences ==


==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''<sub>''i''</sub> and ''x'' = ''x''<sub>''i''</sub> + Δ''x'', where Δ''x'' is an increment along the ''x''-axis ('''Fig. 1'''). We can approximate the derivative d''f''(''x'')/d''x'' at ''x'' = ''x''<sub>''i''</sub> by solving the Taylor’s series,
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''<sub>''i''</sub> and ''x'' = ''x''<sub>''i''</sub> + Δ''x'', where Δ''x'' is an increment along the ''x''-axis ('''Fig. 1'''). We can approximate the derivative d''f''(''x'')/d''x'' at ''x'' = ''x''<sub>''i''</sub> by solving the Taylor’s series,


[[File:Vol1 page 0054 eq 005.png]]....................(3)
[[File:Vol1 page 0054 eq 005.png|RTENOTITLE]]....................(3)


for d''f''(''x'')/d''x''. The result is
for d''f''(''x'')/d''x''. The result is


[[File:Vol1 page 0055 eq 001.png]]....................(4)
[[File:Vol1 page 0055 eq 001.png|RTENOTITLE]]....................(4)


where ''E''<sub>''T''</sub> is the term
where ''E''<sub>''T''</sub> is the term


[[File:Vol1 page 0055 eq 002.png]]....................(5)
[[File:Vol1 page 0055 eq 002.png|RTENOTITLE]]....................(5)


If we neglect ''E''<sub>''T''</sub>, we obtain the finite-difference approximation of the first derivative.
If we neglect ''E''<sub>''T''</sub>, we obtain the finite-difference approximation of the first derivative.


[[File:Vol1 page 0055 eq 003.png]]....................(6)
[[File:Vol1 page 0055 eq 003.png|RTENOTITLE]]....................(6)


'''Eq. 6''' is an approximation because we are neglecting all of the terms in ''E''<sub>''T''</sub>, 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.  
'''Eq. 6''' is an approximation because we are neglecting all of the terms in ''E''<sub>''T''</sub>, 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.


<gallery widths="300px" heights="200px">
<gallery widths="300px" heights="200px">
Line 37: Line 39:
The finite difference in '''Eq. 6''' is called a forward difference. Other differences are possible. Two that we use next are the backward difference,
The finite difference in '''Eq. 6''' is called a forward difference. Other differences are possible. Two that we use next are the backward difference,


[[File:Vol1 page 0055 eq 004.png]]....................(7)
[[File:Vol1 page 0055 eq 004.png|RTENOTITLE]]....................(7)


and the centered difference,
and the centered difference,


[[File:Vol1 page 0055 eq 005.png]]....................(8)
[[File:Vol1 page 0055 eq 005.png|RTENOTITLE]]....................(8)


'''Eqs. 6''' through '''8''' are all derived from Taylor’s series.
'''Eqs. 6''' through '''8''' are all derived from Taylor’s series.


==Numerical solution of the 1D C/D equation==
== 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. 1''', see Chap. 4 of Peacemen.<ref name="r2" /> In our example, we choose a backward difference for the time derivative in '''Eq. 1''', 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. 1''' is converted from a PDE to the difference equation


[[File:Vol1 page 0056 eq 001.png]]....................(9)
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. 1''', see Chap. 4 of Peacemen.<ref name="r1">Peaceman, D.W. 1977. Fundamentals of Numerical Reservoir Simulation. Oxford, UK: Elsevier Publishing.</ref> In our example, we choose a backward difference for the time derivative in '''Eq. 1''', 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. 1''' is converted from a PDE to the difference equation


The subscripts of concentration ''C'' denote points in space, and the superscripts denote points in time. For example, the present time, ''t''<sup>''n''</sup>, is denoted by superscript ''n'' and future time ''t''<sup>''n''+1</sup> is denoted by ''n''+1. The time increment is Δ''t'' = ''t''<sup>''n''+1</sup> - ''t''<sup>''n''</sup>. Similarly, the space increment is Δ''x'' = ''x''<sup>''i''</sup> + 1 - ''x''<sup>''i''</sup>. The concentration at time ''t''<sup>''n''+1</sup> and spatial location xi is denoted by [[File:Vol1 page 0056 inline 001.png]].  
[[File:Vol1 page 0056 eq 001.png|RTENOTITLE]]....................(9)
 
The subscripts of concentration ''C'' denote points in space, and the superscripts denote points in time. For example, the present time, ''t''<sup>''n''</sup>, is denoted by superscript ''n'' and future time ''t''<sup>''n''+1</sup> is denoted by ''n''+1. The time increment is Δ''t'' = ''t''<sup>''n''+1</sup> - ''t''<sup>''n''</sup>. Similarly, the space increment is Δ''x'' = ''x''<sup>''i''</sup> + 1 - ''x''<sup>''i''</sup>. The concentration at time ''t''<sup>''n''+1</sup> and spatial location xi is denoted by [[File:Vol1 page 0056 inline 001.png|RTENOTITLE]].


The future concentration distribution is found from the current concentration distribution by rearranging '''Eq. 9'''. We collect terms in ''C''<sup>''n''+1</sup> on the left-hand side and terms in ''C''<sup>''n''</sup> on the right-hand side, thus
The future concentration distribution is found from the current concentration distribution by rearranging '''Eq. 9'''. We collect terms in ''C''<sup>''n''+1</sup> on the left-hand side and terms in ''C''<sup>''n''</sup> on the right-hand side, thus


[[File:Vol1 page 0056 eq 002.png]]....................(10)
[[File:Vol1 page 0056 eq 002.png|RTENOTITLE]]....................(10)


'''Eq. 10''' is now written in the form
'''Eq. 10''' is now written in the form


[[File:Vol1 page 0056 eq 003.png]]....................(11)
[[File:Vol1 page 0056 eq 003.png|RTENOTITLE]]....................(11)


where the coefficients are
where the coefficients are


[[File:Vol1 page 0056 eq 004.png]]....................(12)
[[File:Vol1 page 0056 eq 004.png|RTENOTITLE]]....................(12)


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 ≤ ''I'' ≤''NX'', the system of finite-difference equations becomes
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 ≤ ''I'' ≤''NX'', the system of finite-difference equations becomes


[[File:Vol1 page 0057 eq 001.png]]....................(13)
[[File:Vol1 page 0057 eq 001.png|RTENOTITLE]]....................(13)


'''Eq. 13''' can be written in matrix form as
'''Eq. 13''' can be written in matrix form as


[[File:Vol1 page 0057 eq 002.png]]....................(14)
[[File:Vol1 page 0057 eq 002.png|RTENOTITLE]]....................(14)


where [[File:Vol1 page 0057 inline 001.png]]  is the ''NX'' × ''NX'' matrix of coefficients, [[File:Vol1 page 0057 inline 002.png]] is the column vector of unknown concentrations at time ''t''<sup>''n''+1</sup>, and [[File:Vol1 page 0057 inline 003.png]] is the column vector of right-hand-side terms that depend on known concentrations at time tn. Both column vectors [[File:Vol1 page 0057 inline 004.png]] and [[File:Vol1 page 0057 inline 003.png]] have ''NX'' elements.  
where [[File:Vol1 page 0057 inline 001.png|RTENOTITLE]]  is the ''NX'' × ''NX'' matrix of coefficients, [[File:Vol1 page 0057 inline 002.png|RTENOTITLE]] is the column vector of unknown concentrations at time ''t''<sup>''n''+1</sup>, and [[File:Vol1 page 0057 inline 003.png|RTENOTITLE]] is the column vector of right-hand-side terms that depend on known concentrations at time tn. Both column vectors [[File:Vol1 page 0057 inline 004.png|RTENOTITLE]] and [[File:Vol1 page 0057 inline 003.png|RTENOTITLE]] have ''NX'' elements.


The system of equations in '''Eq. 14''' 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.<ref name="r2" /><ref name="r3" /><ref name="r4" /><ref name="r5" /><ref name="r9" /> A solution of the set of equations for physical parameters ''v'' = 1 ft/day and ''D'' = 0.01 ft<sup>2</sup>/day and finite-difference parameters Δ''x'' = 0.1 ft and Δ''t'' = 0.1 day is shown in '''Fig. 2'''. The difference between the analytical solution and the numerical solution is because of numerical dispersion,<ref name="r2" /><ref name="r11" /><ref name="r10" /> 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.  
The system of equations in '''Eq. 14''' 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.<ref name="r1">Peaceman, D.W. 1977. Fundamentals of Numerical Reservoir Simulation. Oxford, UK: Elsevier Publishing.</ref><ref name="r2">Aziz, K. and Settari, A. 1979. Petroleum Reservoir Simulation. Essex, UK: Elsevier Applied Science Publishers.</ref><ref name="r3">Mattax, C.C. and Dalton, R.L. 1990. Reservoir Simulation, Vol. 13. Richardson, Texas: Monograph Series, SPE.</ref><ref name="r4">Ertekin, T., Abou-Kassem, J.H., and King, G.R. 2001. Basic Applied Reservoir Simulation, Vol. 7. Richardson, Texas: Textbook Series, SPE.</ref><ref name="r8">Chapra, S.C. and Canale, R.P. 2002. Numerical Methods for Engineers, fourth edition. Boston, Massachusetts: McGraw-Hill Book Co.</ref> A solution of the set of equations for physical parameters ''v'' = 1 ft/day and ''D'' = 0.01 ft<sup>2</sup>/day and finite-difference parameters Δ''x'' = 0.1 ft and Δ''t'' = 0.1 day is shown in '''Fig. 2'''. The difference between the analytical solution and the numerical solution is because of numerical dispersion,<ref name="r1">Peaceman, D.W. 1977. Fundamentals of Numerical Reservoir Simulation. Oxford, UK: Elsevier Publishing.</ref><ref name="r11">Fanchi, J.R. 2006. Math Refresher for Scientists and Engineers, third edition. New York: Wiley Interscience.</ref><ref name="r10">Lantz, R.B. 1971. Quantitative Evaluation of Numerical Diffusion (Truncation Error). SPE J. 11 (3): 315–320. SPE-2811-PA. http://dx.doi.org/10.2118/2811-PA</ref> 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.


<gallery widths="300px" heights="200px">
<gallery widths="300px" heights="200px">
Line 81: Line 84:


== Matrices and linear algebra ==
== Matrices and linear algebra ==
An example of a matrix was introduced earlier 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.<ref name="r1" /> 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.  
 
An example of a matrix was introduced earlier 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.<ref name="r11">Fanchi, J.R. 1983. Multidimentional Numerical Dispersion. SPE J. 23 (1): 143–151. SPE-9018-PA. http://dx.doi.org/10.2118/9018-PA</ref> 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 ===
=== Rotation of a Cartesian coordinate system ===
'''Fig. 3''' illustrates a rotation of Cartesian coordinates from one set of orthogonal coordinates {''x''<sub>1</sub>, ''x''<sub>2</sub>} to another set {''y''<sub>1</sub>, ''y''<sub>2</sub>} by the angle ''θ''. The equations relating the coordinate systems are
'''Fig. 3''' illustrates a rotation of Cartesian coordinates from one set of orthogonal coordinates {''x''<sub>1</sub>, ''x''<sub>2</sub>} to another set {''y''<sub>1</sub>, ''y''<sub>2</sub>} by the angle ''θ''. The equations relating the coordinate systems are


[[File:Vol1 page 0057 eq 003.png]]....................(15)
[[File:Vol1 page 0057 eq 003.png|RTENOTITLE]]....................(15)


The set of equations in '''Eq. 15''' has the matrix form
The set of equations in '''Eq. 15''' has the matrix form


[[File:Vol1 page 0058 eq 001.png]]....................(16)
[[File:Vol1 page 0058 eq 001.png|RTENOTITLE]]....................(16)


which can be written as
which can be written as


[[File:Vol1 page 0058 eq 002.png]]....................(17)
[[File:Vol1 page 0058 eq 002.png|RTENOTITLE]]....................(17)


The column vectors [[File:Vol1 page 0058 inline 001.png]] and [[File:Vol1 page 0058 inline 002.png]] are
The column vectors [[File:Vol1 page 0058 inline 001.png|RTENOTITLE]] and [[File:Vol1 page 0058 inline 002.png|RTENOTITLE]] are


[[File:Vol1 page 0058 eq 003.png]]....................(18)
[[File:Vol1 page 0058 eq 003.png|RTENOTITLE]]....................(18)


with two elements each, and the rotation matrix [[File:Vol1 page 0058 inline 003.png]] is the 2 × 2 square matrix,
with two elements each, and the rotation matrix [[File:Vol1 page 0058 inline 003.png|RTENOTITLE]] is the 2 × 2 square matrix,


[[File:Vol1 page 0058 eq 004.png]]....................(19)
[[File:Vol1 page 0058 eq 004.png|RTENOTITLE]]....................(19)


<gallery widths="300px" heights="200px">
<gallery widths="300px" heights="200px">
Line 108: Line 113:
</gallery>
</gallery>


===Properties of matrices===
=== 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''<sup>th</sup> row and ''j''<sup>th</sup> column of the matrix is the ''ij''<sup>th</sup> 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. 18'''. 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 [[File:Vol1 page 0059 inline 001.png]] in '''Eq. 18''' may be written as  
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''<sup>th</sup> row and ''j''<sup>th</sup> column of the matrix is the ''ij''<sup>th</sup> 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.


[[File:Vol1 page 0059 eq 001.png]]....................(20)
If the matrix has a single column so that ''n'' = 1, it is a column vector as in '''Eq. 18'''. 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 [[File:Vol1 page 0059 inline 001.png|RTENOTITLE]] in '''Eq. 18''' may be written as  


where the superscript ''T'' denotes the transpose of the matrix. In general, we can write a ''m'' × ''n'' matrix [[File:Vol1 page 0059 inline 002.png]] with a set of elements {''a''<sub>''ij''</sub>: ''i'' = 1, 2, ... ''n''; ''j'' = 1, 2, ... ''m''} as
[[File:Vol1 page 0059 eq 001.png|RTENOTITLE]]....................(20)


[[File:Vol1 page 0059 eq 002.png]]....................(21)
where the superscript ''T'' denotes the transpose of the matrix. In general, we can write a ''m'' × ''n'' matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] with a set of elements {''a''<sub>''ij''</sub>: ''i'' = 1, 2, ... ''n''; ''j'' = 1, 2, ... ''m''} as


The transpose of matrix [[File:Vol1 page 0059 inline 002.png]] is
[[File:Vol1 page 0059 eq 002.png|RTENOTITLE]]....................(21)


[[File:Vol1 page 0059 eq 003.png]]....................(22)
The transpose of matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] is


The conjugate transpose of matrix [[File:Vol1 page 0059 inline 002.png]] is obtained by finding the complex conjugate of each element in [[File:Vol1 page 0059 inline 006.png]] and then taking the transpose of the matrix [[File:Vol1 page 0059 inline 006.png]]. This operation can be written as
[[File:Vol1 page 0059 eq 003.png|RTENOTITLE]]....................(22)


[[File:Vol1 page 0059 eq 004.png]]....................(23)
The conjugate transpose of matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] is obtained by finding the complex conjugate of each element in [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]] and then taking the transpose of the matrix [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]]. This operation can be written as


where * denotes complex conjugation. Recall that the conjugate ''z''* of a complex number ''z'' is obtained by replacing the imaginary number [[File:Vol1 page 0059 inline 003.png]] with [[File:Vol1 page 0059 inline 004.png]] wherever it occurs. If all the elements of matrix [[File:Vol1 page 0059 inline 006.png]] are real, the conjugate transpose of matrix [[File:Vol1 page 0059 inline 006.png]] is equal to the transpose of matrix [[File:Vol1 page 0059 inline 005.png]].  
[[File:Vol1 page 0059 eq 004.png|RTENOTITLE]]....................(23)


If the matrix [[File:Vol1 page 0059 inline 002.png]] is a square matrix and the elements of matrix [[File:Vol1 page 0059 inline 002.png]] satisfy the equality ''a''<sub>''ij''</sub> = ''a''<sub>''ji''</sub>, the matrix [[File:Vol1 page 0059 inline 002.png]] is called a symmetric matrix. A square matrix A¯¯ is Hermitian, or self-adjoint, if [[File:Vol1 page 0059 inline 002.png]]= [[File:Vol1 page 0059 inline 002.png]]<sub>+</sub> (i.e, the matrix equals its conjugate transpose).  
where * denotes complex conjugation. Recall that the conjugate ''z''* of a complex number ''z'' is obtained by replacing the imaginary number [[File:Vol1 page 0059 inline 003.png|RTENOTITLE]] with [[File:Vol1 page 0059 inline 004.png|RTENOTITLE]] wherever it occurs. If all the elements of matrix [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]] are real, the conjugate transpose of matrix [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]] is equal to the transpose of matrix [[File:Vol1 page 0059 inline 005.png|RTENOTITLE]].


The set of elements {''a''<sub>''ii''</sub>} of a square matrix [[File:Vol1 page 0059 inline 005.png]] is the principal diagonal of the matrix. The elements {''a''<sub>''ji''</sub>} with ''i'' ≠ ''j'' are off-diagonal elements. The matrix [[File:Vol1 page 0059 inline 002.png]] is a lower triangular matrix if ''a''<sub>''ij''</sub> = 0 for ''i'' < ''j'', and [[File:Vol1 page 0059 inline 006.png]] is an upper triangular matrix if ''a''<sub>''ij''</sub> = 0 for ''i'' > ''j''. The matrix [[File:Vol1 page 0059 inline 006.png]] is a diagonal matrix if ''a''<sub>''ij''</sub> =0 for ''i'' ≠ ''j''.  
If the matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] is a square matrix and the elements of matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] satisfy the equality ''a''<sub>''ij''</sub> = ''a''<sub>''ji''</sub>, the matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] is called a symmetric matrix. A square matrix A¯¯ is Hermitian, or self-adjoint, if [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]]= [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]]<sub>+</sub> (i.e, the matrix equals its conjugate transpose).
 
The set of elements {''a''<sub>''ii''</sub>} of a square matrix [[File:Vol1 page 0059 inline 005.png|RTENOTITLE]] is the principal diagonal of the matrix. The elements {''a''<sub>''ji''</sub>} with ''i'' ≠ ''j'' are off-diagonal elements. The matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] is a lower triangular matrix if ''a''<sub>''ij''</sub> = 0 for ''i'' < ''j'', and [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]] is an upper triangular matrix if ''a''<sub>''ij''</sub> = 0 for ''i'' > ''j''. The matrix [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]] is a diagonal matrix if ''a''<sub>''ij''</sub> =0 for ''i'' ≠ ''j''.


=== Matrix operations ===
=== Matrix operations ===
Suppose the matrices [[File:Vol1 page 0059 inline 005.png]], [[File:Vol1 page 0060 inline 001.png]], and [[File:Vol1 page 0060 inline 002.png]] with elements {''a''<sub>''ij''</sub>}, {''b''<sub>''ij''</sub>}, and {''c''<sub>''ij''</sub>} 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


[[File:Vol1 page 0060 eq 001.png]]....................(24)
Suppose the matrices [[File:Vol1 page 0059 inline 005.png|RTENOTITLE]], [[File:Vol1 page 0060 inline 001.png|RTENOTITLE]], and [[File:Vol1 page 0060 inline 002.png|RTENOTITLE]] with elements {''a''<sub>''ij''</sub>}, {''b''<sub>''ij''</sub>}, and {''c''<sub>''ij''</sub>} 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
 
[[File:Vol1 page 0060 eq 001.png|RTENOTITLE]]....................(24)


The product of a matrix [[File:Vol1 page 0058 inline 003.png]] with a number ''k'' may be written as
The product of a matrix [[File:Vol1 page 0058 inline 003.png|RTENOTITLE]] with a number ''k'' may be written as


[[File:Vol1 page 0060 eq 002.png]]....................(25)
[[File:Vol1 page 0060 eq 002.png|RTENOTITLE]]....................(25)


The product of matrix [[File:Vol1 page 0058 inline 003.png]] with order ''m'' × ''n'' and matrix [[File:Vol1 page 0060 inline 003.png]] with order ''n'' × ''p'' is
The product of matrix [[File:Vol1 page 0058 inline 003.png|RTENOTITLE]] with order ''m'' × ''n'' and matrix [[File:Vol1 page 0060 inline 003.png|RTENOTITLE]] with order ''n'' × ''p'' is


[[File:Vol1 page 0060 eq 003.png]]....................(26)
[[File:Vol1 page 0060 eq 003.png|RTENOTITLE]]....................(26)


where matrix [[File:Vol1 page 0060 inline 002.png]] has order ''m'' × ''p''. Notice that matrix multiplication is possible only if the number of columns in [[File:Vol1 page 0058 inline 003.png]] equals the number of rows in [[File:Vol1 page 0060 inline 003.png]]. This requirement is always satisfied for square matrices.  
where matrix [[File:Vol1 page 0060 inline 002.png|RTENOTITLE]] has order ''m'' × ''p''. Notice that matrix multiplication is possible only if the number of columns in [[File:Vol1 page 0058 inline 003.png|RTENOTITLE]] equals the number of rows in [[File:Vol1 page 0060 inline 003.png|RTENOTITLE]]. This requirement is always satisfied for square matrices.


The transpose of the product of two square matrices [[File:Vol1 page 0059 inline 002.png]] and [[File:Vol1 page 0060 inline 004.png]] is
The transpose of the product of two square matrices [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] and [[File:Vol1 page 0060 inline 004.png|RTENOTITLE]] is


[[File:Vol1 page 0060 eq 004.png]]....................(27)
[[File:Vol1 page 0060 eq 004.png|RTENOTITLE]]....................(27)


and the adjoint of the product of two square matrices is
and the adjoint of the product of two square matrices is


[[File:Vol1 page 0060 eq 005.png]]....................(28)
[[File:Vol1 page 0060 eq 005.png|RTENOTITLE]]....................(28)


Notice that the product of two matrices may not be commutative (i.e., [[File:Vol1 page 0058 inline 003.png]] [[File:Vol1 page 0060 inline 003.png]] ≠ [[File:Vol1 page 0060 inline 003.png]] [[File:Vol1 page 0058 inline 003.png]] in general).  
Notice that the product of two matrices may not be commutative (i.e., [[File:Vol1 page 0058 inline 003.png|RTENOTITLE]] [[File:Vol1 page 0060 inline 003.png|RTENOTITLE]] ≠ [[File:Vol1 page 0060 inline 003.png|RTENOTITLE]] [[File:Vol1 page 0058 inline 003.png|RTENOTITLE]] in general).


The identity matrix, [[File:Vol1 page 0060 inline 006.png]], 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 [[File:Vol1 page 0059 inline 006.png]] in matrix multiplication, thus
The identity matrix, [[File:Vol1 page 0060 inline 006.png|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 [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]] in matrix multiplication, thus


[[File:Vol1 page 0060 eq 006.png]]....................(29)
[[File:Vol1 page 0060 eq 006.png|RTENOTITLE]]....................(29)


By contrast, a null matrix [[File:Vol1 page 0060 inline 007.png]] is a matrix in which all elements are zero. In this case, the product of the null matrix with a matrix [[File:Vol1 page 0059 inline 002.png]] is
By contrast, a null matrix [[File:Vol1 page 0060 inline 007.png|RTENOTITLE]] is a matrix in which all elements are zero. In this case, the product of the null matrix with a matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] is


[[File:Vol1 page 0061 eq 001.png]]....................(30)
[[File:Vol1 page 0061 eq 001.png|RTENOTITLE]]....................(30)


The matrix, [[File:Vol1 page 0059 inline 002.png]], is singular if the product of matrix [[File:Vol1 page 0059 inline 002.png]] with a column vector [[File:Vol1 page 0061 inline 001.png]] that has at least one nonzero element yields the null matrix; that is, [[File:Vol1 page 0059 inline 002.png]] is singular if
The matrix, [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]], is singular if the product of matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] with a column vector [[File:Vol1 page 0061 inline 001.png|RTENOTITLE]] that has at least one nonzero element yields the null matrix; that is, [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] is singular if


[[File:Vol1 page 0061 eq 002.png]]....................(31)
[[File:Vol1 page 0061 eq 002.png|RTENOTITLE]]....................(31)


The concepts of identity matrix and matrix singularity are needed to define the inverse matrix. Suppose we have two square matrices [[File:Vol1 page 0059 inline 006.png]] and [[File:Vol1 page 0061 inline 003.png]] that satisfy the product
The concepts of identity matrix and matrix singularity are needed to define the inverse matrix. Suppose we have two square matrices [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]] and [[File:Vol1 page 0061 inline 003.png|RTENOTITLE]] that satisfy the product


[[File:Vol1 page 0061 eq 003.png]]....................(32)
[[File:Vol1 page 0061 eq 003.png|RTENOTITLE]]....................(32)


Notice that the matrices [[File:Vol1 page 0059 inline 002.png]] and [[File:Vol1 page 0060 inline 004.png]] commute. The matrix [[File:Vol1 page 0059 inline 002.png]] is nonsingular, and the matrix [[File:Vol1 page 0060 inline 004.png]] is the inverse of [[File:Vol1 page 0059 inline 002.png]], thus [[File:Vol1 page 0060 inline 004.png]] = [[File:Vol1 page 0059 inline 002.png]]<sup>-1</sup>, where [[File:Vol1 page 0059 inline 002.png]]<sup>-1</sup> denotes the inverse of [[File:Vol1 page 0059 inline 002.png]]. '''Eq. 32''' can be written in terms of the inverse as
Notice that the matrices [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] and [[File:Vol1 page 0060 inline 004.png|RTENOTITLE]] commute. The matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] is nonsingular, and the matrix [[File:Vol1 page 0060 inline 004.png|RTENOTITLE]] is the inverse of [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]], thus [[File:Vol1 page 0060 inline 004.png|RTENOTITLE]] = [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]]<sup>-1</sup>, where [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]]<sup>-1</sup> denotes the inverse of [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]]. '''Eq. 32''' can be written in terms of the inverse as


[[File:Vol1 page 0061 eq 004.png]]....................(33)
[[File:Vol1 page 0061 eq 004.png|RTENOTITLE]]....................(33)


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


[[File:Vol1 page 0061 eq 005.png]]....................(34)
[[File:Vol1 page 0061 eq 005.png|RTENOTITLE]]....................(34)


where the column vector [[File:Vol1 page 0061 inline 004.png]] and the matrix [[File:Vol1 page 0058 inline 003.png]] are known, and the column vector [[File:Vol1 page 0058 inline 001.png]] contains a set of unknowns. '''Eq. 13''' is an example for the 1D C/D equation. We can solve for [[File:Vol1 page 0061 inline 001.png]] in '''Eq. 34''' by premultiplying '''Eq. 34''' by [[File:Vol1 page 0059 inline 002.png]]<sup>-1</sup>. The result is
where the column vector [[File:Vol1 page 0061 inline 004.png|RTENOTITLE]] and the matrix [[File:Vol1 page 0058 inline 003.png|RTENOTITLE]] are known, and the column vector [[File:Vol1 page 0058 inline 001.png|RTENOTITLE]] contains a set of unknowns. '''Eq. 13''' is an example for the 1D C/D equation. We can solve for [[File:Vol1 page 0061 inline 001.png|RTENOTITLE]] in '''Eq. 34''' by premultiplying '''Eq. 34''' by [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]]<sup>-1</sup>. The result is


[[File:Vol1 page 0061 eq 006.png]]....................(35)
[[File:Vol1 page 0061 eq 006.png|RTENOTITLE]]....................(35)


Of course, we have to know [[File:Vol1 page 0059 inline 006.png]]<sup>-1</sup> to find [[File:Vol1 page 0061 inline 001.png]]. This leads us to a discussion of determinants.  
Of course, we have to know [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]]<sup>-1</sup> to find [[File:Vol1 page 0061 inline 001.png|RTENOTITLE]]. This leads us to a discussion of determinants.


=== Determinants, eigenvalues, and eigenvectors ===
=== Determinants, eigenvalues, and eigenvectors ===
The determinant (det) of a square matrix [[File:Vol1 page 0059 inline 002.png]] is denoted by det [[File:Vol1 page 0061 inline 002.png]] or | [[File:Vol1 page 0059 inline 002.png]] |. 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


[[File:Vol1 page 0061 eq 007.png]]....................(36)
The determinant (det) of a square matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]] is denoted by det [[File:Vol1 page 0061 inline 002.png|RTENOTITLE]] or | [[File:Vol1 page 0059 inline 002.png|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
 
[[File:Vol1 page 0061 eq 007.png|RTENOTITLE]]....................(36)


and the determinant of a 3×3 matrix is
and the determinant of a 3×3 matrix is


[[File:Vol1 page 0061 eq 008.png]]
[[File:Vol1 page 0061 eq 008.png|RTENOTITLE]] [[File:Vol1 page 0062 eq 001.png|RTENOTITLE]]....................(37)
[[File:Vol1 page 0062 eq 001.png]]....................(37)


Determinants are useful for determining if an inverse matrix [[File:Vol1 page 0059 inline 002.png]]<sup>-1</sup> exists. Inverse matrices are needed to solve finite-difference equations representing fluid flow. The condition det [[File:Vol1 page 0062 inline 001.png]]  says that an inverse matrix [[File:Vol1 page 0059 inline 006.png]]<sup>-1</sup> exists, even though we may not know the elements of the inverse matrix. Determinants are also useful for determining eigenvalues and eigenvectors.  
Determinants are useful for determining if an inverse matrix [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]]<sup>-1</sup> exists. Inverse matrices are needed to solve finite-difference equations representing fluid flow. The condition det [[File:Vol1 page 0062 inline 001.png|RTENOTITLE]]  says that an inverse matrix [[File:Vol1 page 0059 inline 006.png|RTENOTITLE]]<sup>-1</sup> 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 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
Eigenvalues are the values of ''λ'' in the eigenvalue equation


[[File:Vol1 page 0062 eq 002.png]]....................(38)
[[File:Vol1 page 0062 eq 002.png|RTENOTITLE]]....................(38)
 
where [[File:Vol1 page 0059 inline 005.png|RTENOTITLE]] is an ''n'' × ''n'' square matrix and [[File:Vol1 page 0058 inline 001.png|RTENOTITLE]] is a column vector with n rows. The eigenvalue equation may be written as


where [[File:Vol1 page 0059 inline 005.png]] is an ''n'' × ''n'' square matrix and [[File:Vol1 page 0058 inline 001.png]] is a column vector with n rows. The eigenvalue equation may be written as
[[File:Vol1 page 0062 eq 003.png|RTENOTITLE]]....................(39)


[[File:Vol1 page 0062 eq 003.png]]....................(39)
where [[File:Vol1 page 0062 inline 002.png|RTENOTITLE]] is the ''n'' × ''n'' identity matrix. '''Eq. 39''' has nonzero solutions, [[File:Vol1 page 0058 inline 001.png|RTENOTITLE]], if the eigenvalue, ''λ'', is a characteristic root of [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]], that is, ''λ'' must be a solution of


where [[File:Vol1 page 0062 inline 002.png]] is the ''n'' × ''n'' identity matrix. '''Eq. 39''' has nonzero solutions, [[File:Vol1 page 0058 inline 001.png]], if the eigenvalue, ''λ'', is a characteristic root of [[File:Vol1 page 0059 inline 002.png]], that is, ''λ'' must be a solution of
[[File:Vol1 page 0062 eq 004.png|RTENOTITLE]]....................(40)<br/><br/>'''Eq. 40''' is the characteristic equation of [[File:Vol1 page 0059 inline 005.png|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. 40''' into an ''n''th-degree polynomial and then solving for the ''n'' values of ''λ''. These concepts are illustrated in the next section.


[[File:Vol1 page 0062 eq 004.png]]....................(40)
== Nomenclature ==
<br>
<br>
'''Eq. 40''' is the characteristic equation of [[File:Vol1 page 0059 inline 005.png]], and the ''n'' values of ''λ'' are the characteristic roots of the characteristic equation. The characteristic roots, ''λ'', are obtained by expanding the determinant in '''Eq. 40''' into an ''n''th-degree polynomial and then solving for the ''n'' values of ''λ''. These concepts are illustrated in the next section.


==Nomenclature==
{|
{|
|''a''<sub>''i''</sub>, ''b''<sub>''i''</sub>, ''c''<sub>''i''</sub>, ''d''<sub>''i''</sub>
|=
|finite-difference coefficients, '''Eq. 11'''
|-
|-
|''a''<sub>''ij''</sub>, ''b''<sub>''ij''</sub>, ''c''<sub>''ij''</sub>  
| ''a''<sub>''i''</sub>, ''b''<sub>''i''</sub>, ''c''<sub>''i''</sub>, ''d''<sub>''i''</sub>
|=  
| =
|elements of matrices, '''Eq. 24'''  
| finite-difference coefficients, '''Eq. 11'''
|-
|-
|[[File:Vol1 page 0074 inline 001.png]]
| ''a''<sub>''ij''</sub>, ''b''<sub>''ij''</sub>, ''c''<sub>''ij''</sub>
|=  
| =
|matrices, '''Eq. 24'''  
| elements of matrices, '''Eq. 24'''
|-
|-
|[[File:Vol1 page 0059 inline 002.png]]  
| [[File:Vol1 page 0074 inline 001.png|RTENOTITLE]]
|=  
| =
|rotation matrix, '''Eq. 17'''  
| matrices, '''Eq. 24'''
|-
|-
|[[File:Vol1 page 0074 inline 003.png]]  
| [[File:Vol1 page 0059 inline 002.png|RTENOTITLE]]
|=  
| =
|column vector of unknown concentrations at ''t''<sup>''n'' +1</sup>, '''Eq. 14'''  
| rotation matrix, '''Eq. 17'''
|-
|-
|[[File:Vol1 page 0074 inline 004.png]]  
| [[File:Vol1 page 0074 inline 003.png|RTENOTITLE]]
|=  
| =
|column vector of terms that depend on known concentrations at ''t''<sup>''n'' +1</sup>, '''Eq. 14'''  
| column vector of unknown concentrations at ''t''<sup>''n'' +1</sup>, '''Eq. 14'''
|-
|-
|''E''<sub>''T''</sub>  
| [[File:Vol1 page 0074 inline 004.png|RTENOTITLE]]
|=
| =
|truncation error, '''Eq. 5'''  
| column vector of terms that depend on known concentrations at ''t''<sup>''n'' +1</sup>, '''Eq. 14'''
|-
|-
|''i''  
| ''E''<sub>''T''</sub>
|=  
| =
|imaginary operator
| truncation error, '''Eq. 5'''
|-
|-
|[[File:Vol1 page 0074 inline 006.png]]
| ''i''
|=  
| =
|identity matrix, '''Eq. 29'''
| imaginary operator
|-
|-
|''J''<sub>''x''</sub>, ''J''<sub>''y''</sub>, ''J''<sub>''z''</sub>
| [[File:Vol1 page 0074 inline 006.png|RTENOTITLE]]
|=  
| =
|fluid flux in ''x''-, ''y''-, ''z''-directions
| identity matrix, '''Eq. 29'''
|-
|-
|(''J''<sub>''x''</sub>)<sub>''x''</sub>  
| ''J''<sub>''x''</sub>, ''J''<sub>''y''</sub>, ''J''<sub>''z''</sub>
|=  
| =
|fluid flux in ''x''-direction at location ''x''  
| fluid flux in ''x''-, ''y''-, ''z''-directions
|-
|-
|(''J''<sub>''y''</sub>)<sub>''y''</sub>  
| (''J''<sub>''x''</sub>)<sub>''x''</sub>
|=  
| =
|fluid flux in ''y''-direction at location ''y''  
| fluid flux in ''x''-direction at location ''x''
|-
|-
|(''J''<sub>''z''</sub>)<sub>''z''</sub>  
| (''J''<sub>''y''</sub>)<sub>''y''</sub>
|=  
| =
|fluid flux in ''z''-direction at location ''z''  
| fluid flux in ''y''-direction at location ''y''
|-
|-
|[[File:Vol1 page 0074 inline 005.png]]
| (''J''<sub>''z''</sub>)<sub>''z''</sub>
|=  
| =
|matrix of coefficients, '''Eq. 14'''  
| fluid flux in ''z''-direction at location ''z''
|-
|-
|''n''
| [[File:Vol1 page 0074 inline 005.png|RTENOTITLE]]
|=  
| =
|exponent
| matrix of coefficients, '''Eq. 14'''
|-
|-
|''R''  
| ''n''
|=  
| =
|region
| exponent
|-
|-
|''S''  
| ''R''
|=  
| =
|surface
| region
|-
|-
|''t''  
| ''S''
|=  
| =
|time
| surface
|-
|-
|''t''<sup>''n''</sup>
| ''t''
|=  
| =
|present time  
| time
|-
|-
|''t''<sup>''n''+1</sup>  
| ''t''<sup>''n''</sup>
|=  
| =
|future time  
| present time
|-
|-
|''v''  
| ''t''<sup>''n''+1</sup>
|=  
| =
|velocity, L/t, ft/sec
| future time
|-
|-
|''x,y'',''z''  
| ''v''
|=  
| =
|space dimensions
| velocity, L/t, ft/sec
|-
|-
|''x''<sub>''i''</sub>
| ''x,y'',''z''
|=  
| =
|discrete point in ''x''-direction, '''Eq. 3'''
| space dimensions
|-
|-
|Δ''t''  
| ''x''<sub>''i''</sub>
|=  
| =
|time interval
| discrete point in ''x''-direction, '''Eq. 3'''
|-
|-
|Δ''x''  
| Δ''t''
|=  
| =
|length
| time interval
|-
|-
|Δ''y''  
| Δ''x''
|=  
| =
|width
| length
|-
|-
|Δ''z''  
| Δ''y''
|=  
| =
|thickness
| width
|-
|-
|''θ''  
| Δ''z''
|=  
| =
|angle, '''Eq. 15'''
| thickness
|-
|-
|''λ''  
| ''θ''
|=  
| =
|eigenvalues, '''Eq. 38'''  
| angle, '''Eq. 15'''
|-
|-
| ''λ''
| =
| eigenvalues, '''Eq. 38'''
|}
|}


==References==
== References ==
<references>
<ref name="r1">Fanchi, J.R. 2006. ''Math Refresher for Scientists and Engineers'', third edition. New York: Wiley Interscience.</ref>


<ref name="r2">Peaceman, D.W. 1977. ''Fundamentals of Numerical Reservoir Simulation''. Oxford, UK: Elsevier Publishing.</ref>
<references />


<ref name="r3">Aziz, K. and Settari, A. 1979. ''Petroleum Reservoir Simulation''. Essex, UK: Elsevier Applied Science Publishers.</ref>
== Noteworthy papers in OnePetro ==


<ref name="r4">Mattax, C.C. and Dalton, R.L. 1990. ''Reservoir Simulation'', Vol. 13. Richardson, Texas: Monograph Series, SPE. </ref>
Use this section to list papers in OnePetro that a reader who wants to learn more should definitely read


<ref name="r5">Ertekin, T., Abou-Kassem, J.H., and  King, G.R. 2001. ''Basic Applied Reservoir Simulation'', Vol. 7. Richardson, Texas: Textbook Series, SPE.</ref>
== External links ==


<ref name="r6">Munka, M. and Pápay, J. 2001. ''4D Numerical Modeling of Petroleum Reservoir Recovery''. Budapest, Hungary: Akadémiai Kiadó.</ref>
Use this section to provide links to relevant material on websites other than PetroWiki and OnePetro


<ref name="r7">Fanchi, J.R. 2006. ''Principles of Applied Reservoir Simulation'', third edition. Burlington, Massachusetts: Gulf Professional Publishing/Elsevier.</ref>
== See also ==


<ref name="r8">Fanchi, J.R. 2000. ''Integrated Flow Modeling'', No. 49. Amsterdam, The Netherlands: Developments in Petroleum Science, Elsevier Science B.V.</ref>
[[Mathematics_of_fluid_flow|Mathematics of fluid flow]]
 
<ref name="r9">Chapra, S.C. and Canale, R.P. 2002. ''Numerical Methods for Engineers'', fourth edition. Boston, Massachusetts: McGraw-Hill Book Co.</ref>
 
<ref name="r10">Lantz, R.B. 1971. Quantitative Evaluation of Numerical Diffusion (Truncation Error). ''SPE J.'' '''11''' (3): 315–320. SPE-2811-PA. http://dx.doi.org/10.2118/2811-PA</ref>
 
<ref name="r11">Fanchi, J.R. 1983. Multidimentional Numerical Dispersion. ''SPE J.'' '''23''' (1): 143–151. SPE-9018-PA. http://dx.doi.org/10.2118/9018-PA</ref>
</references>
 
==Noteworthy papers in OnePetro==
Use this section to list papers in OnePetro that a reader who wants to learn more should definitely read
 
==External links==
Use this section to provide links to relevant material on websites other than PetroWiki and OnePetro


==See also==
[[Vector_analysis_of_fluid_flow|Vector analysis of fluid flow]]
[[Mathematics of fluid flow]]


[[Vector analysis of fluid flow]]
[[Diagonalizing_the_permeability_tensor|Diagonalizing the permeability tensor]]


[[Diagonalizing the permeability tensor]]
[[PEH:Mathematics_of_Fluid_Flow]]


[[PEH:Mathematics of Fluid Flow]]
[[Category:5.3.1 Flow in porous media]]

Latest revision as of 09:54, 8 June 2015

Systems of nonlinear partial differential equations (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. This article provides an overview of these methods.

One-dimensional convection/dispersion equation

To illustrate the mathematics, we discuss the numerical solution of the 1D convection/dispersion (C/D) equation

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

as introduced in vector analysis of fluid flow. As a reminder, v is velocity, D is dispersion, and C is concentration. Eq. 1 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.

RTENOTITLE....................(2)

We first introduce the concept of finite differences to convert Eq. 1 to an equation that can be solved numerically. We then present a numerical representation of Eq. 1 and illustrate its solution. For more details, you should consult the reservoir simulation page, as well as sources in the literature.[1][2][3][4][5][6][7]

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. 1). We can approximate the derivative df(x)/dx at x = xi by solving the Taylor’s series,

RTENOTITLE....................(3)

for df(x)/dx. The result is

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

where ET is the term

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

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

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

Eq. 6 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. 6 is called a forward difference. Other differences are possible. Two that we use next are the backward difference,

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

and the centered difference,

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

Eqs. 6 through 8 are all derived from Taylor’s series.

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. 1, see Chap. 4 of Peacemen.[1] In our example, we choose a backward difference for the time derivative in Eq. 1, 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. 1 is converted from a PDE to the difference equation

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

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. 9. We collect terms in Cn+1 on the left-hand side and terms in Cn on the right-hand side, thus

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

Eq. 10 is now written in the form

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

where the coefficients are

RTENOTITLE....................(12)

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

Eq. 13 can be written in matrix form as

RTENOTITLE....................(14)

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. 14 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.[1][2][3][4][8] 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. The difference between the analytical solution and the numerical solution is because of numerical dispersion,[1][9][10] 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 earlier 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.[9] 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. 3 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....................(15)

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

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

which can be written as

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

The column vectors RTENOTITLE and RTENOTITLE are

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

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

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

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. 18. 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. 18 may be written as  

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

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

The transpose of matrix RTENOTITLE is

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

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

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 A¯¯ 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....................(24)

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

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

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

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

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

and the adjoint of the product of two square matrices is

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

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

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

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

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

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. 32 can be written in terms of the inverse as

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

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

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

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

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

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

and the determinant of a 3×3 matrix is

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

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

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

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

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

Eq. 40 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. 40 into an nth-degree polynomial and then solving for the n values of λ. These concepts are illustrated in the next section.

Nomenclature

ai, bi, ci, di = finite-difference coefficients, Eq. 11
aij, bij, cij = elements of matrices, Eq. 24
RTENOTITLE = matrices, Eq. 24
RTENOTITLE = rotation matrix, Eq. 17
RTENOTITLE = column vector of unknown concentrations at tn +1, Eq. 14
RTENOTITLE = column vector of terms that depend on known concentrations at tn +1, Eq. 14
ET = truncation error, Eq. 5
i = imaginary operator
RTENOTITLE = identity matrix, Eq. 29
Jx, Jy, Jz = fluid flux in x-, y-, z-directions
(Jx)x = fluid flux in x-direction at location x
(Jy)y = fluid flux in y-direction at location y
(Jz)z = fluid flux in z-direction at location z
RTENOTITLE = matrix of coefficients, Eq. 14
n = exponent
R = region
S = surface
t = time
tn = present time
tn+1 = future time
v = velocity, L/t, ft/sec
x,y,z = space dimensions
xi = discrete point in x-direction, Eq. 3
Δt = time interval
Δx = length
Δy = width
Δz = thickness
θ = angle, Eq. 15
λ = eigenvalues, Eq. 38

References

  1. 1.0 1.1 1.2 1.3 Peaceman, D.W. 1977. Fundamentals of Numerical Reservoir Simulation. Oxford, UK: Elsevier Publishing.
  2. 2.0 2.1 Aziz, K. and Settari, A. 1979. Petroleum Reservoir Simulation. Essex, UK: Elsevier Applied Science Publishers.
  3. 3.0 3.1 Mattax, C.C. and Dalton, R.L. 1990. Reservoir Simulation, Vol. 13. Richardson, Texas: Monograph Series, SPE.
  4. 4.0 4.1 Ertekin, T., Abou-Kassem, J.H., and King, G.R. 2001. Basic Applied Reservoir Simulation, Vol. 7. Richardson, Texas: Textbook Series, SPE.
  5. Munka, M. and Pápay, J. 2001. 4D Numerical Modeling of Petroleum Reservoir Recovery. Budapest, Hungary: Akadémiai Kiadó.
  6. Fanchi, J.R. 2006. Principles of Applied Reservoir Simulation, third edition. Burlington, Massachusetts: Gulf Professional Publishing/Elsevier.
  7. Fanchi, J.R. 2000. Integrated Flow Modeling, No. 49. Amsterdam, The Netherlands: Developments in Petroleum Science, Elsevier Science B.V.
  8. Chapra, S.C. and Canale, R.P. 2002. Numerical Methods for Engineers, fourth edition. Boston, Massachusetts: McGraw-Hill Book Co.
  9. 9.0 9.1 Fanchi, J.R. 2006. Math Refresher for Scientists and Engineers, third edition. New York: Wiley Interscience. Cite error: Invalid <ref> tag; name "r11" defined multiple times with different content
  10. Lantz, R.B. 1971. Quantitative Evaluation of Numerical Diffusion (Truncation Error). SPE J. 11 (3): 315–320. SPE-2811-PA. http://dx.doi.org/10.2118/2811-PA

Noteworthy papers in OnePetro

Use this section to list papers in OnePetro that a reader who wants to learn more should definitely read

External links

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

See also

Mathematics of fluid flow

Vector analysis of fluid flow

Diagonalizing the permeability tensor

PEH:Mathematics_of_Fluid_Flow