Numerical Methods in Reservoir Simulation

6
CONTENTS

1. INTRODUCTION 7. NUMERICAL DISPERSION - A

MATHEMATICAL APPROACH
2. REVIEW OF FINITE DIFFERENCES 7.1. Introduction to the Problem
7.2. Mathematical Derivation of Numerical
3. APPLICATION OF FINITE DIFFERENCES TO Dispersion
PARTIAL DIFFERENTIAL EQUATIONS
(PDEs) 8. CLOSING REMARKS
3.1. Explicit Finite Difference Approximation
of the Linear Pressure Equation APPENDIX A: Some Useful Matrix Theorems.
3.2. Implicit Finite Difference Approximation
of the Linear Pressure Equation
3.3. Implicit Finite Difference Approximation
of the 2D Pressure Equation
3.3.1 Discretisation of the 2D Pressure Equation
3.3.2 Numbering Schemes in Solving the 2D
Pressure Equation
3.4. Implicit Finite Difference Approximation
of Non-linear Pressure Equations

4. APPLICATION OF FINITE DIFFERENCES TO

TWO-PHASE FLOW
4.1 Discretisation of the Two-Phase Pressure
and Saturation Equations
4.2 IMPES Strategy for Solving the Two-Phase
Pressure and Saturation Equations

5. THE NUMERICAL SOLUTION OF LINEAR

EQUATIONS
5.1. Introduction to Linear Equations
5.2. General Methods for Solving Linear
Equations
5.3. Direct Methods for Solving Linear
Equations
5.4. Iterative Methods for Solving Linear
Equations
5.5. A Comparison of Iterative and Direct
Methods for Solving Linear Equations

6. DIRECT SOLUTION OF THE NON-LINEAR

EQUATIONS OF MULTI-PHASE FLOW
6.1. Introduction to Sets of Non-linear
Equations
6.2. Newton’s Method for Solving Sets of Non-
linear Equations
6.3. Newton’s Method Applied to the Non-linear
Equations of Two-Phase Flow
LEARNING OBJECTIVES:

Having worked through this chapter the student should be able to:

• write down from memory simple finite difference expressions for derivatives,
(∂P/∂x), (∂P/∂t) and (∂2P/∂x2) explaining your spatial (space) and temporal (time)
n n +1 n +1 n +1
notation ( Pi , Pi , Pi +1 , Pi −1 etc. ); for (∂P/∂x), the student should know the
meaning of the forward difference, the backward difference and the central difference
and the order of the error associated with each, O(Δx) or O(Δx2).

• apply finite difference approximations to a simple partial differential equation

(PDE) such as the diffusion equation and explain what is meant by an explicit and
an implicit numerical scheme.

• write a simple spreadsheet to solve the explicit numerical scheme for a simple
PDE for given boundary and initial conditions and be able to describe the effect of
time step size, Δt.

• show how the implicit finite difference scheme applied to a simple linear PDE
leads to a set of linear equations which are tridiagonal in 1D and pentadiagonal in
2D.

• derive the structure of the pentadiagonal A-matrix in 2D for a given numbering

scheme going from (i, j) notation to m-notation where m is an ordered numbering
scheme e.g. for the natural numbering scheme, m = (j - 1).NX + i

• describe a solution strategy for the non-linear single phase 2D pressure equation
where the fluid and rock compressibility (and density, ρ, and viscosity, μ) are functions
of the dependent variable, pressure, P(x,y,t).

• write down the discretised form of both the pressure and saturation equation for
two-phase flow given the governing equations (in simplified form in 1D), and be able
to explain why these lead to sets of non-linear algebraic equations.

• outline with an explanation and a simple flow chart the main idea behind the IMPES
solution strategy for the discretised two-phase flow equations.

• write down the expanded expressions for a set of linear equations which, in compact
form are written A.x = b, where the matrix A is an nxn matrix of (known) coefficients
(aij; i = rows, j = columns), b is a vector of n (known) values and x is the vector of n
unknowns which we are solving for.

• explain clearly the main differences between a direct and an iterative solution
method for the set of linear equations, A.x = b.

• write down the algorithm for a very simple iterative scheme for solving A.x = b,
and be able to describe the significance of the initial guess, x(0) , what is meant by
iteration (and iteration counter, ν), the idea of convergence of x(ν) as ν → ∞; and be
able to comment on the number of iterations required for convergence, Niter.

2
Numerical Methods in Reservoir Simulation
6
• explain how to apply (without derivation) the Newton-Raphson method for
solving a single non-linear algebraic equation, f ( x ) = 0 ; Newton-Raphson scheme
f ( x (ν ) )
=> x (ν +1) = x (ν ) − .
f ' ( x (ν ) )
• extend the application of the Newton-Raphson to sets of non-linear algebraic
equations such as those arising from the discretisation of the two-phase pressure and
S
saturation equations; F( X ) = 0, where X =   and S and P are the (unknown)
 P
vectors of saturation and pressure; (given the Newton-Raphson expression,
[ ]
X (ν +1) = X (ν ) − J ( X (ν ) ) .F( X (ν ) )).
−1

• understand, but not be able to reproduce the detailed derivation of, the more
mathematical explanation of numerical dispersion.

NUMERICAL METHODS IN RESERVOIR SIMULATION

As noted in Chapter 5, the multi-phase flow equations for real systems are so complex
that it is not remotely possible to solve them analytically. In practice these equations
can only be solved numerically. The most commonly applied numerical methods
are based on finite difference approximations of the flow equations and this approach
will be followed in this chapter.

After a brief review of finite differences, we go on to apply them to very simple

systems such as for the simplified 1D pressure equation (derived in Chapter 5). This
equation does not need to be solved numerically but it demonstrates how explicit
and implicit finite difference solution can be developed. We then show how sets of
linear equations arise in solving implicit equations and we consider solution of the
linear equations in an elementary manner. We then consider how the 2D pressure
equation is solved.

In the numerical solution of the multi-phase flow equations, we need to solve for
pressure and flow. An outline of how this is done is presented but we do not go into
great detail.

1. INTRODUCTION

At the very start of this course, we considered a very simple “simulation model” for a
growing colony of bacteria. The number of bacteria, N, grew with a rate proportional
to N itself i.e. (dN/dt) = α.N, where α is a constant. We saw that this simple equation
α.t
had a well-known analytical solution, N( t ) = N o .e , where No is the number of
bacteria at time, t =0. This exponential growth law provides the solution to our
model. However, we also introduced the idea of a numerical model where, even
although we could solve the problem analytically, we formulated it this way, in any
case. The numerical version of the model came up with the algorithm (or recipe)
N n +1 = (1 + α.∆t ).N n . In the exercise in Chapter 1, you should have compared the
results of the analytical and numerical models and found out that, they get closer as

Institute of Petroleum Engineering, Heriot-Watt University 3

we take successively smaller time steps, Δt. In this case, we say that the numerical
model converges to the analytical model. In fact, in many areas of science and
engineering, we often apply a numerical technique to a problem we can already
solve analytically i.e. where we “know the answer”. Why would we do this? The
answer is that we might be testing the numerical method to see how closely it gives
the right answer. More commonly, we might test several - maybe 3 or 4 - numerical
techniques to determine which one works best. The phrase “works best” in the context
of a numerical method usually means gives the most accurate numerical agreement
with the analytical answer for the least amount of computational work. Note the
importance of this balance between accuracy and work for a numerical method.
There may be no point in having a numerical method that is “twice” as accurate (in
some sense) for ten times the amount of computational work.

In Chapter 5, we already met the equations for single-phase and two-phase flow of
compressible fluids through porous media. These turned out to be non-linear partial
differential equations (PDEs). Recall that a non-linear PDE is one where certain
coefficients in the equation depend on the answer we are trying to find e.g. Sw(x,t),
P(x,y,z,t) etc. For example, for single phase compressible flow, the equation for
pressure, P(x,t), in 1D is given by:

 ∂P  ∂  k.ρ  ∂P  
c( P)  =  
 ∂t  ∂x  µ  ∂x  
(1)

In this equation, both the generalised compressibility term, c(P) (of both the rock and
the fluid), and the density, ρ(P), terms depend on the unknown pressure, P, which
we are trying to find. As noted previously, such non-linear PDEs are very difficult
to solve analytically and we must usually resort to numerical methods. The main
topic of this module is on how we solve the reservoir flow equations numerically.
This process involves the following steps:

(i) Firstly, we must take the PDE describing the flow process and “chop it up” into
grid blocks in space. This is known as spatial discretisation and, in this course,
we will exclusively use finite difference methods for this purpose.

(ii) When we apply finite differences, we usually end up with sets of non-linear
algebraic equations which are still quite difficult to solve. In some cases, we do
solve these non-linear equations. However, we usually linearise these equations
in order to obtain a set of linear equations.

(iii) We then solve the resulting sets of linear equations. Many numerical options
are available for solving sets of linear equations and these will be discussed below.
This is often done iteratively by repeatedly solving them until the numerical solution
converges.

This module will deal successively with each of the parts of the numerical solution
process, (i) - (iii) above.

4
Numerical Methods in Reservoir Simulation
6

2. REVIEW OF FINITE DIFFERENCES

Definition: A finite difference scheme is simply a way of approximating

derivatives of a function (e.g.  dP  ,  ∂P  ,  ∂ P  etc.)

2

 dx   ∂t   ∂x 
2

numerically from either point or block values of the function.

The main concept of finite difference approximation is best illustrated by the

following simple example where we refer to Figure 1. Study the Notation in this
figure, since it is the basis of that used throughout this chapter. The main task of the
finite difference approach is to represent the derivatives of the function, P(x), in an
approximate manner
 dP  ,  ∂P  ,  ∂ P  etc.
2
   2
i.e.  dx   ∂t   ∂x 
First consider how we might approximate (dP/dx) at xi using the quantities in Figure
1. In fact, it is easily seen that there are three ways we may do this as follows:

Approximation 1 -Forward Differences (fd): we may take the slope between Pi and
Pi+1 as being approximately  dP  at xi to obtain:
 dx  i

 dP  = Pi +1 − Pi
 dx  i fd ∆x (2)

Approximation 2 -Backward Differences (bd): we may take the slope between Pi-1
and Pi as being approximately 
dP  at x to obtain:
 dx  i
i

 dP  = Pi − Pi −1
 dx  i bd ∆x (3)

Approximation 3 -Central Differences (cd): thirdly, we could take the average of

the forward difference (fd) and backward difference (bd) approximations to give us
 dP  at x . This is known as central differences (cd) and is given by:
 dx  i i

 dP  = 1  dP  +  dP   = 1  Pi +1 − Pi + Pi − Pi −1 
 dx  i cd 2  dx  i fd  dx  i bd  2  ∆x ∆x 
 dP  P −P
=  i +1 i −1 
 dx  i cd  2.∆x  (4)

Institute of Petroleum Engineering, Heriot-Watt University 5

Pi+1

P(x)
Pi

Pi-1

∆x = constant
∆x ∆x
xi-1 xi xi-1 Figure 1
x
Notation for the application
Notation:
∆x = the x- grid spacing or distance between grid points of finite difference methods
(i - 1), i, (i + 1) = the x - label (subscript) for the grid point or block numbers
for approximating
Pi-1, Pi, Pi+1 = the corresponding function values at grid point i-1, i, i+1 etc.
derivatives.

 dP 
Each of the above approximations to is shown graphically in Figure 2.
 dx  i
You may wonder why we bother with three different numerical approximations to
 dP 
 dx  i and ask: which is “best”? There is not an unqualified answer to this question
just yet, but let us take a simple numerical example where we know the right answer
and examine each of the forward, backward and central difference approximations.
The values given in Figure 3 will illustrate how the methods perform. Firstly, simply
calculate the values given by each methods for the data in Figure 3:

 dP  ≈  3.0042 − 2.7183  = 2.859 [err ≈ +0.14]

 dx  i fd  0.1 

 dP  ≈  2.7183 − 2.4596  = 2.587 [err ≈ −0.13]

 dx  i bd  0.1 

 dP  ≈  3.0042 − 2.7183  = 2.723 [err ≈ −0.005]

 dx  i cd  2 x 0.1 
(5)

The quantity in square brackets after each of the finite difference approximations
above is the error i.e. the difference between that method and the right answer which
is 2.7183. As expected from the figure for this case, the forward difference answer
is a little too high (by ≈ +0.14) and the backward difference answer is a little too
low (by ≈ -0.13). The central difference approximation is rather better than either
of the previous ones. In fact, we note that the fd and bd methods give an error of
order Δx, the grid spacing (where, as engineers, we are saying 0.14 ≈ 0.1 !). The
cd approximation, on the other hand, has an error of order Δx2 (where again we are

6
Numerical Methods in Reservoir Simulation
6
saying 0.005 (0.1x0.1 = 0.01). More formally, we say that the error in the fd and bd
approximations are "of order Δx" which we denote, O(Δx), and the cd approximation
has an error "of order Δx2" which we denote, O(Δx2).
 ∂2 P 
 2
Now consider the finite difference approximation of the second derivative,  ∂x 
 ∂P 
 
Going back to Figure 1, the definition of second derivative at  ∂x  xi is the rate of
change of slope (dP/dx) at xi. Therefore, we can evaluate this derivative between xi-1
and xi (i.e. the bd approximation) and then do the same between xi and xi+1 (i.e. the
fd approximation) and simply take the rate of change of these two quantities with
respect to x, as follows:

 dP  −  dP 
 ∂ P
2  dx  i fd  dx  i bd
 2 ≈
 ∂x  ∆x

 Pi +1 − Pi  −  Pi − Pi −1 
 ∆x  i fd  ∆x  i bd
≈
∆x

 ∂2 P  ( Pi +1 + Pi −1 − 2 Pi )
 2 ≈
 ∂x  ∆x 2 (6)

 ∂2 P 
 2
Calculating the numerical value of  ∂x  for the example in Figure 3 (where again
 ∂2 P 
 2  = 2.7183 since the example is the exponential function), gives:
 ∂x 

 ∂ 2 P  3.0042 + 2.4596 − 2 x 2.7183

 2 ≈ = 2.7200 [err ≈ 0.0017]
 ∂x  0.12
(7)

Thus, we can see that the error in this case (err ≈ 0.0017 ≈ 0.12) is actually O(Δx2).

In the introductory section of Chapter 1, we already applied the idea of finite differences
(although we didn’t call it that at the time) to the simple ordinary differential equation
dN 
(ODE);  = α .N .
 dt 
Here, N (number of bacteria in the colony) as a function of time, N(t), is the unknown
we want to find. Denoting Nn and Nn+1 as the size of the colony at times t (labelled
n) and t+Δt (labelled n+1), we applied finite differences to obtain:

Institute of Petroleum Engineering, Heriot-Watt University 7

n +1
 dN  ≈ N − N ≈ α . N n
n

 dt  ∆t

N n +1 ≈ (1 + α .∆t ) N n (8)

This gave us our very simple algorithm to explicitly calculate Nn+1. We then took
this as the current value of N and applied the algorithm repeatedly.

dP Pi+1 - Pi-1
=
dx i cd 2.∆x
Pi+1

P(x)
dP Pi+1 - Pi
=
Pi dx i fd ∆x Figure 2
Pi-1 Graphical illustration of the
dP Pi - Pi-1 finite difference derivatives
=
dx i bd ∆x calculated by backward
∆x ∆x differences (bd), forward
xi-1 xi xi+1 differences (cd) and central
x ∆x = constant
differences.

dP = 2.7183 d2P = 2.7183

dx x=1.0 dx2 x=1.0

Pi+1 3.0042

P(x)

Pi 2.7183

Pi-1 2.4596
Figure 3
A numerical example where
0.1 0.1
the function values and
0.9 1.0 1.1 derivative values are known
x
(P(x) = ex)

8
Numerical Methods in Reservoir Simulation
6

EXERCISE 1.

Apply finite differences to the solution of the equation:

 dy  = 2. y 2 + 4
 dt 
where, at t = 0, y(t = 0) = 1. Take time steps of Δt = 0.001 (arbitrary time units) and
step the solution forward to t = 0.25. Use the notation yn+1 for the (unknown) y at
n+1 time level and yn for the (known) y value at the current, n, time level.

Plot the numerically calculated y as a function of t between t = 0 and t = 0.25 and

plot it against the analytical value (do the integral to find this).

Answer: is given below where the working is shown in spreadsheet CHAP6Ex1.xls.

This gives the finite difference formula, a spreadsheet implementing it and the analytical
solution for comparison.

3. APPLICATION OF FINITE DIFFERENCES TO PARTIAL

DIFFERENTIAL EQUATIONS (PDEs)

3.1 Explicit Finite Difference Approximation of the Linear Pressure

Equation
We have seen in Chapter 5 that the flow equations are actually partial differential
equations (PDEs) since the unknowns, P(x,t) and Sw(x,t) say, depend on both space and
time. As an example of a linear PDE, we will take the simplified pressure equation
(equation 26; Chapter 5) as follows:

 ∂P  k  ∂2 P 
 =
 ∂t  cφµ  ∂x 2 
(9)

k
where the constant cφµ is the hydraulic diffusivity, which we denoted previously
by Dh. As we noted in Chapter 5, this PDE is linear which has known analytical
solutions for various boundary conditions. However, we will again neglect these
and apply numerical methods as an example of how to use finite differences to solve
PDEs numerically. To make things even simpler, we will take Dh = 1, giving the
equation:

 ∂P   ∂ P 
2
  = 2
 ∂t   ∂x 
(10)

This is the pressure equation for a 1D system where 0 ≤ x ≤ L, where L is the length
of the system. We can visualise this physically - much as we did in Chapter 5, Section
2.1 - using Figure 4. After the system is held constant at P = Po, the inlet pressure is
raised (at x = 0) instantly to P = Pin while the outlet pressure is held at Pout = Po.

Institute of Petroleum Engineering, Heriot-Watt University 9

Pin

t = t1 t = t2
P Figure 4
Physical picture of pressure
Po Pout = Po propagation in a 1D (com-
t=0 pressible) system described
0 x L
by.  ∂P   ∂ 2 P 
  = 2
 ∂t   ∂x 
These pressures, Pin and Pout, represent the constant pressure boundary conditions
(Mathematically these are sometimes called Dirichlet boundary conditions). In other
words, these are set, as if by experiment, and the system between 0 < x < L must “sort
itself out” or "respond" by simply obeying equation 10 above. We now approach
this problem using finite differences as follows:

• discretise the x-direction by dividing it into a numerical grid of size Δx;

• choose a time step, Δt;
• use the following notation

→ time level n = 0, 1, 2 ...

Pin
→ x-grid block label, i = 1, 2, 3 ... NX (at x = L)
Pinn +1
P
Piin current (known) P at time level
n +1
P
Pin +1 new (unknown) P at time level
i

• fix the boundary conditions which, in this case, are as follows (see Figure
4):

P1 = P in and PNX = Pout = Po which are fixed for all t.

• apply finite differences to equation 10 using the above notation to obtain:

 ∂P  P n +1 − Pi n
  ≈ i
 ∂t  i ∆t (11)

and

 ∂2 P  Pi ??+1 + Pi ??−1 − 2 Pi ??
 2 ≈
 ∂x  i ∆x 2 (12)

However, an issue arises in equation 12 above as shown by the question marks on

the spatial derivative time levels. It is simply: which time level should we take for
the spatial derivative terms in equation 12? This is important and we will return to
this matter soon. However, for the moment, let us take these spatial derivatives at
time level n (the “known” level) since this will turn out to be the simplest thing we
can do. Thus, we obtain:

10
Numerical Methods in Reservoir Simulation
6
 ∂2 P  Pi n+1 + Pi n−1 − 2 Pi n
 2 ≈
 ∂x  i ∆x 2 (13)

Equating the numerical finite difference approximations of each of the above derivatives
as required by the original PDE, equation 10 (i.e. equating the expressions in equations
11 and 13) gives:

Pi n +1 − Pi n Pi n+1 + Pi n−1 − 2 Pi n
≈
∆t ∆x (14)

which easily rearranges to obtain an explicit expression for, Pi n +1, the only
unknown in the above equation:

∆t
Pi n +1 = Pi n +
∆x 2
( Pi n+1 + Pi n−1 − 2 Pi n )
(15)

In words, we can interpret this above algorithm as saying:

New (n+1 level) value of Pi = Old (n level) value of Pi + a “correction term”

Equation 15 gives the algorithm for propagating the solution of the PDE forward in
time from the given set of initial conditions.

We now consider how to set the initial conditions. The initial conditions are the values
of all the P0i (i = 1, 2, 3 ... , NX) at t = 0. From Figure 4, these are clearly:
P10
P10
P1 = 1 for all time (boundary condition)
Pi 00
Pi 0 = 0 at time t = 0 for 2 < i < NX-1
Pi
0
PNX 0
= Po for all time (boundary condition)
PNX 0
PNX
Let us now apply the above algorithm to the solution of equation 10. For this
problem, suppose we take the following data:

0 ≤ x ≤ 1.0

Δx = 0.1 => implies 11 grid points, P1 , P2 , P3, .... , P11 (NX = 11).

Δt = 0.001

n +1
The problem is then to calculate the solution, P(x,t), - that is Pi for all i (at all
grid points) and all future times up to some final time (possibly when the equation
comes to a steady-state as will happen in this example). This can be done by filling
in the “solution chart” Table 1 - see Exercise 2 below.

Institute of Petroleum Engineering, Heriot-Watt University 11

Grid Blocks

Time
x= 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
t n (BC↓) (BC↓)
i= 1 2 3 4 5 6 7 8 9 10 11
0 0 P0
i
→ 1 0 0 0 0 0 0 0 0 0 0
(IC)
0.1 100 P100
i 1 0

0.2 200 P200

i
1 0

0.3 300 P300

i
1 0

0.4 400 P400

i
1 0

0.5 500 P500

i
1 0

0.6 600 P600 1 0

i

0.7 700 P700

i
1 0

0.8 800 P800

i
1 0

0.9 900 P900

i
1 0
Table 1.
1 1000 P1000 1 0 Solution chart for the solu-
i
Note: BC = Boundary Conditions - these points are fixed; tion of the simple pressure
IC = Initial Conditions, i.e. values of P(x, t = 0) for all i at t= 0 equation

EXERCISE 2.

Fill in the above table using the algorithm (where (Δt/Δx2)=0.1):

∆t
2 ( i +1
Pi n +1 = Pi n + P n + Pi n−1 − 2 Pi n )
∆x

Hint: make up a spread sheet as above and set the first unknown block (shown grey
shaded in table above) with the above formula. Copy this and paste it into all of the
cells in the entire unknown area (surrounded by bold border above).

Answer: If you get stuck, look at spreadsheet CHAP6Ex2.xls on the disk.

12
Numerical Methods in Reservoir Simulation
6
The assumption we made in equation 12 was that the spatial derivative was taken at
the n (known) time level. This allowed us to develop an explicit formula for Pi n +1
(for all i). This method is therefore known as an explicit finite difference method.
We can learn about some interesting and useful features of the explicit method
which is encapsulated in equation 15 by simply experimenting with the spreadsheet
CHAP6Ex2.xls. Three features of this method can be illustrated by “numerical
experiment” as follows:

(i) the effect of time step size, Δt;

(ii) the effect of refining the spatial grid size, Δx (or number of grid cells, NX);

(iii) the effect of running the calculation to steady-state as t → ∞ (in practice,

until the numerically computed solution stops changing).

It is best if you do these yourself by modifying the spreadsheet (CHAP6Ex2.xls).

EXERCISE 3.

Experiment with the spreadsheet in CHAP6Ex2.xls to examine the effects of the

three quantities above - Δt, Δx (or NX) and the solution as t → ∞.

HINTS for Exercise 6.3:

• Sensitivity to Δt: When you make Δt too big, the predicted numerical solution of
the PDE goes badly wrong - indeed negative pressures can occur which is physically
impossible. In fact, the solution has become unstable for the larger time steps.
This means that this explicit numerical method does have some time step limitations
which we must be careful of.

• Sensitivity to Δx: the effects of grid refinement are that, as the grid blocks get
smaller (i.e. Δx → 0 or NX → ∞). The answer should get more accurate although,
to make this happen, you need to reduce the time step as well.

• Behaviour as t → ∞ : Finally, considering the long-time behavior of the

solution of the PDE, you should find that P(x, t → ∞) tends to a straight line. Some
other intermediate pressure profiles can also be plotted using CHAP6Ex2.xls. In
fact, a little bit of analysis shows us that this is quite expected. As t → ∞, then if
steady-state is reached, then:

 ∂P   ∂2 P 
  = 0, which implies =>  2  = 0
 ∂t   ∂x 

But the “curve” with a zero second derivative (i.e. a first derivative which is constant)
is a straight line. Therefore, this result is physically reasonable and our numerical
model appears to be behaving correctly.

Institute of Petroleum Engineering, Heriot-Watt University 13

Another way to represent this explicit finite difference solution to the PDE is shown
in Figure 5, where we indicate the time levels for the solution of the equations and
we also show the dependency of the unknown pressure ( Pi n +1).

n n n
Time P i-1 Pi P i+1
∆x ∆x
Time level n
i-1 i i+1 Figure 5.
∆t Schematic of the explicit
finite difference algorithm
Time level n + 1
n+1 n+1 n+1 for solving the simple
P i-1 Pi P i+1
pressure equation (a PDE).

3.2 Implicit Finite Difference Approximation of the Linear Pressure

Equation
We now return to the original finite difference equation 12, where we had to
 ∂2 P 
 2
make a choice of time level for the spatial derivative,  ∂x  . We now examine
the consequences of taking this derivative at the (n+1) time level - this is the
“unknown” time level. The finite difference equation for this case is as follows:
The time derivative is the same, i.e.

 ∂P  P n +1 − Pi n
  ≈ i
 ∂t  i ∆t (16)

but the spatial derivative now becomes:

 ∂2 P  Pi n+1+1 + Pi n−1+1 − 2 Pi n +1
 2 ≈
 ∂x  i ∆x 2 (17)

As before, we now equate the numerical finite difference approximations of each of the
above derivatives (as required by the original PDE, equation 10) to obtain:

Pi n +1 − Pi n Pi n+1+1 + Pi n−1+1 − 2 Pi n +1
≈
∆t ∆x (18)

n +1
This equation shouldPibe compared with equation 14. In the previous case, this could
easily be rearranged intoPi +1 an explicit equation for Pi n +1 (equation 15). However,
n

equation 18 above cannot be rearranged to give a simple expression for the pressure
n +1 n +1
at the new time step P i −1 , Pi
(n+1), .and Pi n+1+1 there now appear to be three unknowns at
Indeed,
Pi n−1+1 , Ptoi n +be
time level (n+1), viz. Pi n−1+1 , Pi n +1 and Pi n+1+1 . This appears
1 n +1
a bitPiof
and +1 a paradox:
how do we find three unknowns from a single equation (equation 18 above)? The
answer is not really too difficult: Basically, we have an equation - like equation 18
- at every grid point. We will show below that this leads to a set of linear equations
where we have exactly the same number of unknowns as we have linear equations.

14
Numerical Methods in Reservoir Simulation
6
Therefore, if these equations are all linearly independent, then we can solve them
numerically. This is a bit more trouble than our earlier simple explicit finite difference
method. Because, we do not get an explicit expression for our unknowns - instead
we get an implicit set of equations - this method is known as an implicit finite
difference method.

To see where these linear equations come from, rearrange equation 18 above
to obtain:

 ∆x 2  n +1  ∆x 2  n
Pi n−1+1 −  2 + P
 i
n + 1 + P n +1
= −  P
∆t Pi
i +1
  ∆t  i (19)

where all the unknowns ( Pi n−1+1 , Pi n +1 and Pi n+1+1 ) are on the LHS of the equation
and the term on the RHS is “known”, since it is at the old time level, n. We
can write this equation as follows:

ai −1 Pi n−1+1 + ai Pi n +1 + ai +1 Pi n+1+1 = bi (20)

 ∆x 2   ∆x 2  n
and where ai −1 = 1; ai = − 2 + ; a
 i +1 = 1 and bi = −  P
 ∆t   ∆t  i
are all constants. The ai do not change throughout the calculation but the quantity bi is
updated at each time step as the newly calculated Pn+1 is set to the Pn for the next time
step. We can see how this works for the 5 grid block system in Figure 6 below:

1 2 3 4 5
n n n n n
P1 P2 P3 P4 P5
Time level n ∆x

∆t
Time level n=1
Figure 6.
P n+1
1
P n+1
2
P n+1
3
P n+1
4
P n+1
5
Simple example of a 5 grid
block system showing how
Boundary Boundary
the implicit finite difference condition condition
scheme is applied (fixed) (fixed)

At each grid point, then (a) we know the value from the boundary condition (i =
1 and i = 5), (b) it uses a boundary condition (i=2 and i=4) or (c) it is specified
completely by equation 20 (only i = 3, in this case - but it would be most points
for a large number of grid points).

Consider each grid point in turn as follows:

i = 1 a boundary; therefore P1 is fixed, say as P1

Institute of Petroleum Engineering, Heriot-Watt University 15

 ∆x 2  n
P1 + a2 P2n +1 + P3n +1 = b2 = − P
 ∆t  2
i=2
 ∆x 2  n
=> a2 P2n +1 + P3n +1 = −  P − P1
 ∆t  2

 ∆x 2  n
i=3 P2n +1 + a3 P3n +1 + P4n +1 = b3 = − P
 ∆t  3

 ∆x 2  n
P3n +1 + a4 P4n +1 + P5 = b4 = − P
 ∆t  4
i=4
 ∆x 2  n
=> P3n +1 + a4 P4n +1 = −  P − P5
 ∆t  4

iP=5 5 a boundary; therefore P5 is fixed, say as P5

Therefore there are only three unknowns in the above set of linear equations
( P2n +as
( P2n +1 , P3n +1 and P4n +1 ) which can be summarised
1
,P n +1
and P4n +1 )
follows:
3

 ∆x 2  n
i = 2: a2 P2n +1 + P3n +1 = − P − P1
 ∆t  2

n +1 n +1 n +1  ∆x 2  n
i = 3: P2 + aP 3 3 + P 4 = − P
 ∆t  3

n +1 n +1  ∆x 2  n
i = 4: P 3 + aP 4 4 = − P − P5
 ∆t  4
(21)

Note that the set of linear equations above can be represented as a simple
matrix equation as follows:

  ∆x 2  n 
−  P2 − P1 
 P2n +1   ∆t  
 a2 1 0
  ∆x 2  
   n +1 
a3  P3 = −  P3
n

1 1 
  ∆t  
 0 1 a4   n +1 
 P4   
− ∆x  P n −
2
P5 
  ∆t  4 
(22)

16
Numerical Methods in Reservoir Simulation
6
The structure of this matrix equation
P1 and Pis clearer when there are more equations
12
involved. For example, it is quite easy to show that, if we take 12 grid points instead
of the 5 above, we obtain 10 equations (using the two fixed boundary conditions,
P1 and P12 ) for the quantities, P2n +1 , P3n +1 , P4n +1 ....P11n +1, of the form:

 − ( ∆x 2 ∆t ) P2n − P1 
n +1 n +1 n +1 n +1  P2n +1   
P2a2 , P
13 ,0P 4 0 ....P110 0 0 0 0 0   n +1   − ( ∆x 2 ∆t ) P3n 
1 a3 1 0 0 0 0 0 0 0   P3   
   n +1   − ( ∆x 2 ∆t ) P4n 
0 1 a4 0 0 0 0 0 0 0   P4   
   P n +1   − ( ∆x 2 ∆t ) P5n

0 0 1 a5 1 0 0 0 0 0   5   
0 0 0 1 a6 1 0 0 0 0   P6n +1   − ( ∆x 2 ∆t ) P6n 
   n +1  =  
0 0 0 0 1 a7 1 0 0 0   P7   − ( ∆x 2 ∆t ) P7n 
0  
0 0 0 0 1 a8 1 0 0   n +1 
 − ( ∆x 2
   P8  ∆t ) P8n 
0 0 0 0 0 0 1 a9 1 0   P n +1   
   9   − ( ∆x 2 ∆t ) P9n 
0 0 0 0 0 0 0 1 a10 1   P10n +1   
0 0 0 0 0 0 0 0 1 a11   n +1   − ( ∆x 2 ∆t ) P10n 
  
 P11 
 − ( ∆x 2 ∆t ) P11n − P12 

(23)

And 20 grid points in 1D would lead to the following set of 18 equations:

 − ( ∆x 2 ∆t ) P2n − P1 
 
 P2n +1   − ( ∆x 2 ∆t ) P3n 
 n +1   
a2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 P3   − ( ∆x 2 ∆t ) P4n

1 a3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   n +1   
   P4   − ( ∆x 2 ∆t ) P5n 
0 1 a4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0   P n +1   

0 1 a5 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
  5   − ( ∆x 2 ∆t ) P6n

0  P6n +1   
0 0 0 1 a6 1 0 0 0 0 0 0 0 0 0 0 0 0   n +1   − ( ∆x 2 ∆t ) P7n 
   P7   
0 0 0 0 1 a7 1 0 0 0 0 0 0 0 0 0 0 0   − ( ∆x 2
 n +1  ∆t ) P8n 
0 0 0 0 0 1 a8 1 0 0 0 0 0 0 0 0 0 0   P8   
   P n +1   − ( ∆x 2 ∆t ) P9n 
0 0 0 0 0 0 1 a9 1 0 0 0 0 0 0 0 0 0   9   
   P10n +1   − ( ∆x 2 ∆t ) P10n 
0 0 0 0 0 0 0 1 a10 1 0 0 0 0 0 0 0 0 
0
 n +1  =  
0 0 0 0 0 0 0 1 a11 1 0 0 0 0 0 0 0   P11  −( ∆x 2 ∆t ) P11n 
   n +1   
0 0 0 0 0 0 0 0 0 1 a12 1 0 0 0 0 0 0   − ( ∆x 2
0
 P12  ∆t ) P12n 
0 0 0 0 0 0 0 0 0 1 a13 1 0 0 0 0 0   P n +1   
   − ( ∆x 2
0 0 0 0 0 0 0 0 0 0 0 1 a14 1 0 0 0 0 
 13  ∆t ) P13n 
 P n +1   
   − ( ∆x 2
0 0 0 0 0 0 0 0 0 0 0 0 1 a15 1 0 0 0   14  ∆t ) P14n 
 P15n +1   
0 0 0 0 0 0 0 0 0 0 0 0 0 1 a16 1 0 0   − ( ∆x 2
 
 n +1  ∆t ) P15n 
 P16   
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 a17 1 0   − ( ∆x 2
0
 n +1  ∆t ) P16n 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 a18 1   P17   
   − ( ∆x 2
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 a19 
 P n +1  ∆t ) P17n 
 18   
 P19n +1   − ( ∆x 2 ∆t ) P18n 
 
 ( ∆t ) P19 − P20 
 − ∆x 2 n

(24)

Institute of Petroleum Engineering, Heriot-Watt University 17

For this simple 1D PDE, it is clear that the matrix arising from our implicit finite
difference method has the following properties:

(i) It is tridiagonal - that is, it has a maximum of three non-zero elements in any
row and these are symmetric around the central diagonal;

(ii) It is very sparse - that is, most of the elements are zero. In an MxM matrix,
there are only 3M non-zero terms but M2 actual elements. If M = 100, then the
matrix is only (300/1002)x100% = 3% filled with non-zero terms.

As it happens, a very simple computational procedure, called the Thomas algorithm,

can be used to solve tridiagonal systems very quickly. The FORTRAN code
for this is shown for interest in Figure 7. Note that it is very compact and
quite simple in structure. You are not expected to know this or to necessarily
understand how this algorithm works.

18
Numerical Methods in Reservoir Simulation
6
THE THOMAS ALGORITHM

subroutine thomas
c Thomas algorithm for tridiagonal systems
+ (N ! matrix size
+ ,a ! diagonal
+ ,b ! super-diagonal
+ ,c ! sub-diagonal
+ ,s ! rhs
+ ,x ! solution
+ )
c----------------------------------------
c This program accompanies the book:
c C. Pozrikidis; Numerical Computation in Science and Engineering
c Oxford University Press, 1998
c------------------------------------------------
c Coefficient matrix:
c | a1 b1 0 0 ... 0 0 0 |
c | c2 a2 b2 0 ... 0 0 0 |
c | 0 c3 a3 b3 ... 0 0 0 |
c | .............................. |
c | 0 0 0 0 ... cn-1 an-1 bn-1 |
c | 0 0 0 0 ... 0 cn an |
c------------------------------------------
Implicit Double Precision (a-h,o-z)
Dimension a(200),b(200),c(200),s(200),x(200)
Dimension d(200),y(200)
Parameter (tol=0.000000001) this is a measure of how close to coverage we are
c prepare
Na = N-1
c reduction to upper bidiagonal
d(1) = b(1)/a(1)
y(1) = s(1)/a(1)
DO i=1,Na
i1 = i+1
Den = a(i1)-c(i1)*d(i)
d(i1) = b(i1)/Den
y(i1) =(s(i1)-c(i1)*y(i))/Den
End DO
c Back substitution
x(N) = y(N)
DO i=Na,1,-1
x(i)= y(i)-d(i)*x(i+1)
End DO
c Verification and alarm
Res = s(1)-a(1)*x(1)-b(1)*x(2)
If(abs(Res).gt.tol) write (6,*) “ thomas: alarm”
DO i=2,Na
Res = s(i)-c(i)*x(i-1)-a(i)*x(i)-b(i)*x(i+1)
If(abs(Res).gt.tol) write (6,*) “ thomas: alarm”
End DO
Figure 7. Res = s(N)-c(N)*x(N-1)-a(N)*x(N)
The Thomas algorithm for If(abs(Res).gt.tol) write (6,*) “ thomas: alarm”
the solution of tridiagonal c Done
matrix systems. 100 Format (1x,f15.10)
Return
End
Institute of Petroleum Engineering, Heriot-Watt University 19
Finally, we note that the implicit finite difference method can be viewed as shown
in Figure 8 below. The choice of the (n+1) time level for the spatial discretisation
leads to a set of linear equations which are somewhat more difficult to solve than the
simple algebraic equation that arises in the explicit method (equation 15). However,
there are exactly the same number of linear equations as there are unknowns and so
it is possible to solve these.

n n n
Time P i-1 Pi P i+1
∆x ∆x
Time level n
i-1 i i+1 Figure 8.
∆t
Schematic of the implicit
finite difference algorithm
Time level n + 1
n+1 n+1 n+1 for solving the simple pres-
P i-1 Pi P i+1
sure PDE

3.3 Implicit Finite Difference Approximation of the 2D Pressure

Equation
Before we go on to discuss how to solve the sets of linear equations that arise in the
finite difference approximation of PDEs arising in reservoir simulation, we will first
consider the discretisation of the 2D single-phase pressure equation. This raises
some additional important issues which occur when we try to solve more complicated
systems, as follows:

(i) the complication of the more “connected” grid block system in a 2D

domain;

(ii) the possibility of heterogeneity in the permeability field which leads us to

the matter of how to take average properties in grid-to-grid flows between
blocks of different permeability (dealt with in Chapter 4);

(iii) the issue of non-linearity for a compressible system e.g. c(P) is clearly a
function of P(x,t), which is the “unknown”.

3.3.1 Discretisation of the 2D Pressure Equation

We take as the basic pressure equation for single phase slightly compressible flow,
the following (see the solution for exercise 2 at the end of Chapter 5):

cµφ  ∂P  ∂  ˜  ∂P   ∂   ∂P  
 = kx   +  k̃ y   
k  ∂t  ∂x   ∂x   ∂y   ∂y   (25)
k˜ x and k˜ y  cφµ 
 k 
 cφµ  k˜ x and k˜ y
(a  
k simplified form of equation 39, Chapter 5) where the term k is a constant
which we will denote as β below, k is an average permeability of the entire
reservoir and k˜ x and k˜ y are the local permeabilities normalised by k ; i.e.
k non-linearities for
k˜ x = k x / k and k˜ y = ky / k. Note that we have avoided the
the moment (the quantities c(P), μ and ρ usually depend on pressure). However, the
k k˜ x = k x / k and k˜ y = ky / k.

20
k˜ x = k x / k and k˜ y = ky / k.
Numerical Methods in Reservoir Simulation
6
system may be anisotropic (kx ≠ ky) and heterogeneous (the permeability may vary
from grid block to grid block).

Equation 25 above can be discretised in a similar way to that applied to the 1D

linear pressure PDE discussed above. However, we will need to be quite clear
about our notation in 2D and, for this purpose, we refer to Figure 9. This shows the
discretisation grid for the above PDE - note that this is essentially the opposite of
what we did when we derived the equation in the first place! In Chapter 5, we used
a control volume (or grid block) to express the mass conservation and then inserted
Darcy’s law for the block to block flows; we then took limits as Δx, Δy and Δt →
0. Here, we are starting with the PDE and going back to the local conservation of
flows and introducing finite size Δx and Δy.

i, j + 1

(j + 1/2)

∆y
i - 1, j i, j i + 1, j

(j - 1/2)

Figure 9. i, j - 1
y (j)
Discretisation and notation
∆x
for the 2D pressure 1/ )
(i - 2 (i + 1/2)
equation. x (i)

Discretising the following equation using the above notation

 ∂P  ∂   ∂P   ∂   ∂P  
β   = k˜x    + k˜y   
 ∂t  ∂x   ∂x   ∂y   ∂y  
(26)

we obtain:
 ˜  ∂P     ∂P    ˜  ∂P     ∂P  
−  k˜ x    k y    −  k˜ y   
 k x  ∂x  
n +1
P −P  
n
 i +1/ 2   ∂x   i −1/ 2   ∂y   j +1/ 2   ∂y   j −1/ 2
β  ≈ +
 ∆t  ∆x ∆y
(27)

where the (i ± 1/2) and (j ± 1/2) subscripts refer to quantities at the boundaries as
shown in Figure 9. We can now expand these boundary flows as follows:

Institute of Petroleum Engineering, Heriot-Watt University 21

 P n +1 − Pi ,nj   ˜  Pi n+1+,1j − Pi ,nj+1    ˜  Pi ,nj+1 − Pi n−1+,1j  
β  i, j
 ∆t  =  kx
 
( ) 
i +1 / 2  ∆x 2   
( )
  −  kx 
i −1 / 2  ∆x 2 
 

  Pi ,nj++11 − Pi ,nj+1    ˜  Pi ,nj+1 − Pi ,nj+−11  


( )
+  k˜ y 
j +1 / 2  ∆y 2
( )
  −  ky
  

j −1 / 2  ∆y 2 
 
(28)

(˜ ) (˜ ) (˜ )
where the quantities k x i +1/ 2 , k x i −1/ 2 , k y j +1/ 2 and k y j −1/ 2 are some type of (˜ )
average permeabilities between the two neighbouring grid blocks - as discussed in
Chapter 4. Note also we have chosen the spatial discretisation terms at the new (n+1)
time level making the this an implicit finite difference scheme.

We can rearrange equation 28 above by taking all the unknown terms (at n+1) to the
LHS and the known terms (at time level n) to the RHS. This gives the following:

−
( )
 k˜ x
i −1 / 2

 P n +1
− 
( )
 k˜ x
i +1 / 2

 P n +1
+ 
( )
 k˜
x
i −1 / 2
+
( )
k˜ x
i +1 / 2
+
( )
k˜ y
j −1 / 2
+
( )
k˜ y
j +1 / 2
−

β n +1
 Pi , j
i −1, j i +1, j
 ∆x 2   ∆x 2   ∆x
2
∆x 2 ∆y 2 ∆y 2 ∆t 
     

−
( )
 k˜ x
j −1 / 2

 P n +1
− 
( )
 k˜ y
j +1 / 2

 Pi ,nj++11 =
β n
Pi , j
i , j −1
 ∆y   ∆y  ∆t
2 2
   
(29)

Since the coefficients in equation 29 above are constants, then this defines a set of
linear equations similar to those found in 1D. However, here we have up to five
non-zero terms per grid block to deal with rather than the three we found for the
1D system. The matrix which arises in this 2D case is known as a pentadiagonal
matrix. This set of linear equations can be written as follows:

ai −1, j .Pi n−1+,1j + ai +1, j .Pi n+1+,1j + ai , j .Pi ,nj+1 + ai , j −1 .Pi ,nj+−11 + ai , j +1 .Pi ,nj++11 = bi , j
(30)

where the constant coefficients, ai −1, j , ai +1, j , ai , j , ai , j −1 and ai , j +1 , are given by the
coefficients in equation 29; bi, j is also a (know) constant.

3.3.2 Numbering Schemes in Solving the 2D Pressure Equation

It is quite convenient to label the pressures as Pi, j when we are working out the
discretisation of the equations, but this is not helpful when we are arranging the
linear equations. Here, it is useful first to consider the numbering scheme for the
2D system that allows us to dispense with the (i,j) subscripting in equations 29 or 30
above. The structure of the A - matrix in equation 30 can be made clearer by working
out a specific 2D example as shown in Figure 10.

22
Numerical Methods in Reservoir Simulation
6
j=5 m = 17 m = 18 m = 19 m = 20
j=4 m = 13 m = 14 m = 15 m = 16
j=3 m=9 m = 10 m = 11 m = 12
j=2 m=5 m=6 m=7 m=8
j=1 m=1 m=2 m=3 m=4
i=1 i=2 i=3 i=4
Notation:
NX = maximum number of grid blocks in x - direction, i = NX;
Figure 10. NY = maximum number of grid blocks in y - direction, j = NY
Numbering system for 2D m = grid block number in the natural ordering scheme shown
m = (j - 1).NX + i
grid conversion from (i,j) →
e.g. for i = 3, j = 4 and NX = 4, m = (4 - 1).4 + 3 = 15 (as above)
m counter.

In the m-notation shown in Figure 10 (m = (j-1)NX +i), the reordered equations 30

become the following:

a˜ m−1 .Pmn−+11 + a˜ m+1 .Pmn++11 + a˜ m .Pmn +1 + a˜ m− NX .Pmn−+1NX + a˜ m+ NX .Pmn++1NX = bm

(31)

where the ãm are the reordered coefficients where the subscript is calculated from
the m-formula in Figure 10. For example:

ai , j −1 → a˜ ( j −1−1) NX +i =( j −1) NX +i − NX → a˜ m− NX (32)

ai , j +1 → a˜ ( j +1−1) NX +i = ( j −1) NX +i + NX → a˜ m+ NX (33)

Note that, when we apply the above equation numbering scheme to the example in
Figure 10 (NX = 4, NY = 5 and therefore, 1 ≤ m ≤ 20), certain “neighbours” are
“missing” since a block is at the boundary (or in a corner where two neighbours are
missing). For example, for block (i = 3; j = 3), that is block m = 9, the “j-1 block” is
not there. Therefore, the coefficient Am-1 =0 in this case. This is best seen by writing
out the structure of the 20 x 20 A-matrix by referring to Figure 10; the A-matrix
structure is shown in Figure 11.

Note that the A-matrix structure in Figure 11 is sparse and has a maximum of five
non-zero coefficients in a given row - it is a pentadiagonal matrix.

All implicit methods for discretising the pressure equation lead to sets of linear
equations. These have the general matrix form:

A.x = b (34)

where A is a matrix like the examples shown above, x is the column vector of unknowns
(like the pressures) and b is a column vector of the RHSs. This is just like equation
31 but it is in shorthand form. We will discuss methods for solving these equations
later in this chapter. For the meantime, we will just assume that it can be solved. We
next consider when the PDEs describing a phenomenon are non-linear PDEs.

Institute of Petroleum Engineering, Heriot-Watt University 23

m 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1 x x x
2 x x x x
3 x x x x
4 x x x
5 x x x x
6 x x x x x
7 x x x x x
8 x x x x
9 x x x x
10 x x x x x
11 x x x x x
12 x x x x
13 x x x x
14 x x x x x
15 x x x x x
16 x x x x
17 x x x
18 x x x x
19 x x x x
20 x x x

m - ordering scheme used in figure 10 - shown here for reference.

Figure 11.
j=5 m = 17 m = 18 m = 19 m = 20
j=4 m = 13 m = 14 m = 15 m = 16 Structure of the A-matrix
j=3 m=9 m = 10 m = 11 m = 12 for the m-ordering scheme
j=2 m=5 m=6 m=7 m=8
j=1 m=1 m=2 m=3 m=4
in the table shown below for
i=1 i=2 i=3 i=4 reference.

3.4 Implicit Finite Difference Approximation of Non-linear Pressure

Equations
We have noted previously that using numerical methods are usually the only way
which we can solve non-linear PDEs of the type that arise in reservoir simulation.
We will use the single-phase 2D equation for the flow of compressible fluids (and
rocks) as the example for discussing the numerical solution of non-linear PDEs.
Equation 35 below was derived in Chapter 5 and has the form:

 ∂P  ∂  kx ρ  ∂P   ∂  ky ρ  ∂P  
c( P)  =   +
 ∂t  ∂x  µ  ∂x   ∂y  µ  ∂y  
(35)

where the non-linearities arise in this equation due to the dependence of the generalised
compressibility, c(P), the density, ρ(P) and the viscosity, μ(P), on pressure, P(x,t). It
is the pressure that is the main unknown in this equation and, hence, if these quantities
depend on it, then they introduce difficulties into the equations. In fact, we will see
shortly that the equations that arise are not linear equations - they are non-linear
algebraic equations.

24
Numerical Methods in Reservoir Simulation
6
Proceeding as before, we can discretise the above equation as follows:

 kx ρ  ∂P    kx ρ  ∂P    ky ρ  ∂P    k ρ  ∂P  
    −         −  y  
 P n +1 − P n   µ  ∂x   i +1/ 2  µ  ∂x   i −1/ 2  µ  ∂y   j +1/ 2  µ  ∂y   j −1/ 2
c( P) = +
 ∆t  ∆x ∆y
(36)

where we have not yet decided on the time level of the non-linearities (i.e. the time
level of the pressure, Pn or Pn+1, at which we evaluate c, ρ and μ). The most difficult
case would be if these were set at the “unknown” time level Pn+1 and this is what we
will do as follows:

 k ρ  n +1  ∂P    k ρ  n +1  ∂P  
 x     −  x    
 P n +1
− P n
 
 µ   ∂x    µ   ∂x   i −1/ 2
c( P n +1 ) =
i +1 / 2

 ∆t  ∆x
 ky ρ  n +1  ∂P    ky ρ  n +1  ∂P  
    −    
 µ   ∂y   j +1/ 2  µ   ∂y   j −1/ 2
+
∆y
(37)
 ∂P   ∂P 
Now expanding up the   and   terms gives:
 ∂x   ∂y 

 k ρ  n +1  Pi n+1+,1j − Pi ,nj+1    k ρ  n +1  Pi ,nj+1 − Pi n−1+,1j  

 x    −  µ  
x
 µ   ∆ x  ∆x
n +1  P n +1
i, j − P n
i, j  
 i +1 / 2   
 
 i −1 / 2   
c( P )  =
 ∆t  ∆x
 ky ρ  n +1  Pi ,nj++11 − Pi ,nj+1    ky ρ  n +1  Pi ,n1+1 − Pi ,nj+−11  
     −  µ   
 µ  j +1/ 2  ∆y    j −1 / 2 
∆y  
+
∆y
(38)
which can be simplified to the following:
n +1 n +1
 P n +1 − Pi ,nj   kx ρ   kρ 
c( P n +1
) i , j
 ∆t  = 2
  µ.∆x  i +1/ 2
( )
Pi n+1+,1j − Pi ,nj+1 −  x 2 
 µ.∆x  i −1/ 2
(
Pi ,nj+1 − Pi n−1+,1j )

n +1 n +1
 kρ   kρ   Pi ,nj+1 − Pi ,nj+−11 
+ y 2
 µ.∆y  j +1/ 2
(
Pi ,nj++11 − Pi ,nj+1 ) −  y 2
 µ.∆y  j −1/ 2  ∆y 


(39)

Institute of Petroleum Engineering, Heriot-Watt University 25

This can be treated as before to gather together similar unknown P-terms as
follows:
n +1 n +1
 kρ   kρ 
−  x 2 .Pi n−1+,1j −  y 2  .Pi ,nj+−11
 µ.∆x  i −1/ 2  µ.∆y  j −1/ 2
 c( P n +1 )   k ρ  n +1  k ρ  n +1  ky ρ  n +1  ky ρ 
n +1
 n +1
+  +
 
x
2
+ 
x
2
+  2
+  2 .Pi , j
 ∆t   µ.∆x  i +1/ 2  µ.∆x  i −1/ 2  µ.∆y  j +1/ 2  µ.∆y  j −1/ 2 
n +1 n +1
 kρ   kρ   c( P n +1 ).Pi ,nj 
− x 2 .Pi n+1+,1j −  y 2  .Pi ,nj++11 =  
 µ.∆x  i +1/ 2  µ.∆y  j +1/ 2  ∆t 
(40)

As before, we can write this in a compact form as follows:

α in−+11, j .Pi n−1+,1j + α in++11, j .Pi n+1+,1j + α in, +j 1 .Pi,nj+1 + α in, +j −11 .Pi,nj+−11 + α in, +j +11 .Pi,nj++11 = βin, +j 1
(41)

where the elements of theA-matrix (now denoted α in−+11, j , α in, +j −11 , α in, +j 1 , α in++11, j and α in, +j +11 )
are not constants since they depend on the (unknown) value of Pn+1.

How do we go about solving the non-linear set of algebraic equations in 39 or 40

above? It turns out that we have two choices in tackling this more difficult problem
as follows:

(i) We can actually use a numerical equation solver which is specifically designed
to solve more difficult non-linear problems. An example of this type of approach
is in using the Newton-Raphson method. We will return to this method later (once
we have seen how to solve linear equations).

(i) We can choose to apply a more pragmatic algorithms such as the following:
(a) Although our α-terms in equation 41 are strictly at time level (n+1), simply
take them at time level n, as a first guess. So we approximate equation 41 by the
following first guess:

α in−1, j .Pi n−1+,1j + α in+1, j .Pi n+1+,1j + α in, j .Pi,nj+1 + α in, j −1 .Pi,nj+−11 + α in, j +1 .Pi,nj++11 = βin, j
Pi ,nj+1 α in−1, j , α in+1, j , α in, j , α i(42)
n
, j −1 , α i , j +1 and β i , j
n n

n +1
Pi , at
where the quantities α in−1, j , α in+1, j , α in, j , α in, j −1 , α in, j +1 and βin, j are evaluated j the

(known) time level

α ν n - i.e. at values of pressure = Pi, j .
i, j
n
Pi ,nj+1
Pi , 42 n +1 αν
Pi ,nj equations
(b) Solve the now linear j above to obtain a first estimate - or ia, j first
( ) α iν, j
ν
iteration - of the Pi ,nj+1at each grid point. We will use the following notation ⇒
( )
ν ν
where ν is the iteration counter and α i , j is the value of the a-coefficient at the Pi ,nj+1
value after ν iterations; that is:
(P ) n +1 ν
i, j

[(
α iν, j = α i , j Pi,nj+1 )]
ν
(P )
n +1 ν
i, j
(43)

26
Numerical Methods in Reservoir Simulation
6
Pi ,nj+1

ν
(c) Use the latest iterated values of the α i , j to solve the (now linear) equations to
(P ) n +1 ν
( )
ν +1
go from i, j to Pi ,nj+1 .

(d) Keep iterating the above scheme untilPit ( )

n +1 ν
i , j converges; i.e. the difference between
the sum of two successive iterated values of the pressure (Err.) is sufficiently small
(< Tol, which is an acceptably small value)

∑ (P ) n +1 ν +1
( )
ν
Err = i, j − Pi ,nj+1
all i , j (44)

Stop if Err < Tol. Otherwise continue through steps above.

The algorithm outlined in (ii) above for solving the non-linear pressure equation
is represented in Figure 12.

Set the iteration counter ν = 0

Calculate the α iν−1, j , α iν+1, j , α iν, j , α iν, j −1 , α iν, j +1 and β iν, j at

( )
n ν
values of pressure = Pi , j

[(
α iν, j = α i , j Pin, j+1 )]
ν
where ν = current iteration number.

Solve the set of linear equations:

ν = ν +1
α iν−1, j .Pin−+11, j + α iν+1, j .Pin++11, j + α iν, j .Pin, j+1 + α iν, j −1 .Pin, j+−11 + α iν, j +1 .Pin, j++11 = β iν, j

( )
ν +1
to obtain Pin, j+1 = Pin, j+1

Calculate
∑ (P ) n +1 ν +1
( )
ν
Err = i, j − Pi ,nj+1
all i , j

No - continue iterations
Err < Tol ?

Figure 12.
Algorithm for the numerical Yes (the method has converged)
solution of the non-linear
Stop
2D pressure equations:

Institute of Petroleum Engineering, Heriot-Watt University 27

4 APPLICATION OF FINITE DIFFERENCES TO TWO-PHASE
FLOW

4.1 Discretisation of the Two-Phase Pressure and Saturation

Equations
We now consider how finite difference methods are applied to the two-phase flow
equations. We have seen how these equations were derived in Chapter 5. Recall that,
in two-phase flow, we have two coupled equations to solve - a pressure equation and
a saturation equation. For example, we may solve for the oil pressure, Po(x,t) and the
oil saturation, So (x,t); we would then find the water saturation and water pressure by
using the constraint equations, So+Sw = 1 and the capillary pressure relation, Pc(Sw)
= Po - Pw, respectively.

From Chapter 5, we use the highly simplified form of the 1D pressure and saturation
equations, where we choose to solve for ⇒ P(x,t) and So (x,t) as shown in Chapter
5, equations 81 and 82. Note that we have taken zero capillary pressure (therefore
P = Po = Pw) and zero gravity which gives:

∂   ∂P  
PRESSURE EQUATION  λ T ( So )   = 0 (45)
∂x   ∂x  

 ∂S  ∂   ∂P  
SATURATION EQUATION φ  o  = λ o ( So )   (46)
 ∂t  ∂x   ∂x  

(Equations 45 and 46 are the same as equations 81 and 82 from Chapter 5,

λ (S )
respectively).
T o

The quantity λ T ( So ) is the total mobility and is the sum of the oil and water mobilities;
λ T ( So ) = λ o ( So ) + λ w ( So ) . The above two equations are clearly coupled together since
the oil saturation appears in the non-linear coefficient of the pressure equation. Likewise,
λ T ( So ) in
the pressure appears (So ) + λterm
= λthe
o flow w
(Soon) the RHS of the saturation equation.
We can now apply finite differences to each of these equations in the same
way as discussed above to obtain, for the pressure equation (see notation in
Figure 8):

  ∂P     ∂P  
λ T ( So ) ∂x   − λ T ( So )  
 ∂x   i −1/ 2
  i +1 / 2 
=0
∆x (47)

which can be expanded further to give:

  Pi n+1+1 − Pi n +1     Pi n +1 − Pi n−1+1  
λ ( S
 T o  )
∆x
 − λ ( S
  i +1/ 2  T o 
)
∆x

  i −1 / 2

=0
∆x (48)

28
Numerical Methods in Reservoir Simulation
(λ (S )) T o i +1 / 2
and λ T (So ) ( ) i −1 / 2
6
(λ (S )) ( )
n +1
T o and λSTo(So )
i +1 / 2 i −1 / 2

(So )of the and

λ Tlevel
In equation 48 above, we must now specify the time λ T (Somobility
non-linear ) ( ) ( )
(
Son +1 λ T (So ) ) ( ( ) )
λ T (Sλon +T1()So )
i +1 / 2 i −1 / 2
terms, and
; we must choose whether we take these terms
i +1 / 2 i +1 / 2i −1 / 2
at time n or at time level (n+1)? Again, we go straight to the most difficult case by
n +1
defining these terms at the later time level i.e. at So . These terms are denoted as,
(λ (S S )) and (λ (S )) , in equation 48 above to obtain:
n +1
n +1 n +1
o
T o T o
i +1 / 2 i −1 / 2

  (λ ((λS (S)) ))
n +1
 (λ ( S )) 
( )
n +1 nT+1 o

( ) n +1 o i +1 / 2
λ
(λ (S )) ∆x ( P − P ) −  ∆x ( P − P ) = 0
S
n +1 T o
T   i +1 / 2 n +1 n +1 T o
i −1 / 2 n +1 n +1
i +21 / 2 i +1 i 2 i i −1
T o
i −1 / 2
  
(λ (S ))  n +1

(λ (S ))
o
n +1 (49) T i −1 / 2
T o
i −1 / 2
The above equation can now be arranged into the usual order as follows:

(λ (S ))
T
n +1
o
i −1 / 2
P n +1
−
T (
 λ ( Son +1 ) ) i −1 / 2
+
(λ (S )) T
n +1
o
i +1 / 2

 P n +1 + T o
(
λ ( S n +1 ) ) i +1 / 2
Pi n+1+1 = 0
i −1
∆x 2
 ∆x 2 ∆x 2
 i ∆x 2
 
(50)
n +1 n +1 n +1
P i −1 , Pi and P i +1
This is a non-linear set of algebraic equations for the unknown pressures,
Pi n−1+1 , Pi n +1 and Pi n+1+1 , but the coefficients depend on the - also unknown - saturations,
Son +1. Note that equation 50 represents one of a set of equations since there is one
at each grid point (and at the ends of the 1D system the values may be set by the
Son +1
boundary conditions).

We now consider the discretisation of the saturation equation 46. Finite differences
may be applied to this equation as follows:

  ∂P     ∂P  
λ o ( So ) ∂x  
n +1
− λ o ( So )  
 ∂x   i −1/ 2
S −S   n
 i +1 / 2 
φ =
oi oi
 ∆t  ∆x (51)

Again, we can expand the derivative terms at the (i+1/2) and (i-1/2) boundaries (Figure
9) and we can take the mobility terms at the (n+1) time level to obtain:

 Soin+1 − Soin   λ o ( So )
 n +1
( ) (
  λ ( Son +1 ) ) 
φ
 ∆t  
=
∆x 2
i +1 / 2
(P n +1
i +1 − Pi n +1
) −  o ∆x 2 i −1 / 2
(P
i
n +1
−P n +1
i −1 )
   
(52)

The above non-linear algebraic set of equations can be written in various ways. Two
particularly useful ways to write the above equations are as follows:

Institute of Petroleum Engineering, Heriot-Watt University 29

Form A:

∆t  λ o ( So )
 n +1
( ) (
  λ ( Son +1 ) ) 
Sn +1
oi =S + 
n

φ 
oi
∆x 2
i +1 / 2
(P n +1
i +1 − Pi n +1
) −  o ∆x 2 i −1 / 2
(P
i
n +1
−P n +1
i −1 )
    
(53)
or Form B:

∆t  λ o ( So )
 n +1
( ) (
  λ ( Son +1 ) ) 
n +1  
Sn +1
oi −S − 
n

φ 
oi
∆x 2
i +1 / 2
(P n +1
i +1 − Pi n +1
) −  o ∆x 2 i −1 / 2
(i
P n +1
− Pi −1 )  = 0

    
(54)

The reason for writing the above two forms of the saturation equation is that each is
useful, depending on how we intend to approach the solution of the coupled pressure
(equation 45) and saturation (equation 46) equations.

As with the case of the solution of the non-linear single phase pressure equation for
a compressible system, we have two strategies which we can use to solve the two-
phase equations above, as follows:

(i) We can actually use a numerical equation solver which is specifically designed to
solve more difficult non-linear problems; e.g. the Newton-Raphson method. We will
return to this method later (once we have seen how to solve linear equations).

(i) We can again choose to apply a more pragmatic algorithm and this is the subject
of section 4.2 below.

4.2 IMPES Strategy for Solving the Two-Phase Pressure and

Saturation Equations
The form of the discretised pressure and saturation equations which we will take are
as follows (equations 50 and 53):

Pressure:

(λ (S ))
T
n +1
o
i −1 / 2
P n +1
−
(
 λ ( Son +1 )
T ) i −1 / 2
+
(λ (S ))
T
n +1
o
i +1 / 2

 P n +1 + T o
(
λ ( S n +1 ) ) i +1 / 2
Pi n+1+1 = 0
i −1
∆x 2
 ∆x 2 ∆x 2
 i ∆x 2
 
(50)

Saturation:

 (
∆t  λ o ( So )
n +1
)   λ ( Son +1 ) ( ) 
n +1  
Soin+1 + Soin − 
φ  ∆x 2
i +1 / 2
( i +1 i ) −  o ∆x 2
P n +1
− P n +1  i −1 / 2
(i
P n +1
− Pi −1 ) 

   
(53)

30
Numerical Methods in Reservoir Simulation
6
Pressure equation 50 would be a set of linear equations if the coefficients were
known at the current time step (n), rather than being specified at the (n+1) - i.e. the
unknown - time level. However, we could solve equation 50 above as if it were
Pi n−1+1 , Pi n +1 and Pi n+1+1
a linear system of equations by time-lagging Pi n−1+1as
the coefficients n +1
, Pbefore and
P 1 n +nthe
n +(see
P +i1 1+1 and Pi n+1+1
i i −1 , iP
algorithm in Figure 13). This would give us a “first guess” (or first iteration) to find
Pi n−n1++11, Pi nn++11and Pi n+n1++11
the unknowns i.e. the quantities Pin−+11 , Pi and Pi +1 . The same problem exists
P
forPthen +1saturation
, P n +1 equation
and P n +1 53, if wei had the first guess at the Pi n +1 (by thePprocedure i
n +1
i −1 i i +1
just mentioned), then we could use Pi nthese
+1 latest values of pressure and still time lag
n +1
the coefficients (the oil mobility Sterms)Pni +1 and use the saturation equation 53 as if it
+1
werePi n an
o
+1 explicit expression. This would give us an updated Sonvalue of Son +1 which
could then be used back in the pressure Sonn++11 equation and the whole process could be
So n +1 νpreviously.
( ) ( )
iterated until convergence, as discussed
( () ) ( ) A flow chart of this procedure
ν ν
Pi −1 , Pi n +1 and Pi n+1+1 .n +1 ν n +1n +ν1 ν n +1n +ν1 ν
hasSalready been outlined in Chapter 5 (Figure 9) but is elaborated
o
n + 1
( ) (( ))
P i −1 ,herePi Pin
i −1
Figure
and, Pi Pi +1 and.

is linearlised and can then be solved

(( )) (( )) (( ))
13. Note that by taking time-laggednvalues +1 ν ofn +the
1 ν saturations,
Pi −n1 +1 ν, Pi n +1 ν and Pi +n1 +1 ν.
Pin−+implicitly
n +1 νthe pressure equation
for the pressure for that iteration, ν
11 ν , Pin +1 ν and Pin++11 ν .
n +1 ν n +1 ν ( P ) , ( P ) and ( P )P obtained
⇒ ( P ) , ( P ) and ( P ) . The saturations can then (be
n +1 ν i −1 i
) , ( P( Pexplicitly
) )and, ( P( P ) )and
i +1 n +1 ν ν ν
n +1n + 1 n +1n +ν1
ν
i −1 i i +1 i −1 i i −1 i i +1
using the latest pressures (i.e. the ( P ) , ( P ) and ( P ) ) and the most recent
n +1 ν n +1 ν n +1 ν

( P ) , ( P ) and ( P )
i −n1 +1 ν i n +1 ν i +n1 +1 ν

iteration of the saturation (i.e. the ( S ) , ( S ) etc.). Therefore, this approach is

ni −+1 ν in +1 ν i +1

( P as) the
known i −1 , ( PIMPES
n +1 ν
) and
i ( P ) which stands for IMplicit(Sin Pressure,
n +1 ν
approximation,
n +1 ν
i +1
i i +1
) , (S(S Explicit
) )etc.
, ( S ) etc.
i
n +1 ν n +1n +ν1 ν
i +1i
n +1 ν
i +1
in Saturation. This can be iterated(until
S )it,converges
((SS )) etc.
n +1 ν n +1 ν
although there are limitations in
( S ) , i n +1 ν
etc. i +n1+1 ν
the size of time step, Δt, which can be taken. If Δt is too large, the IMPES method
i i +1

may(Sbecome
i ) , (Sunstable
n +1 ν
) etc.and give unphysical results like the example in Exercise 3
n +1 ν
i +1
above.

Institute of Petroleum Engineering, Heriot-Watt University 31

Set the iteration counter to zero, ν = 0

Take current values of Son and P n as the ν level

iteration values: ⇒ Soν and P ν

PRESSURE EQUN. Set the total mobilities ⇒ λ T (Soν )i +1/ 2 and λ T (Soν )i −1/ 2
obtain the linearised pressure equation:
(λ (S ))
T
ν
o
 λ (Soν )
Pin−+11 − 
T
i −1 / 2 i −1 / 2
+
( )
λ T (Soν )  (
i +1 / 2  n +1
Pi +
λ T (Soν ) )
i +1 / 2
Pin++11 = 0
( )
∆x 2
 ∆x 2
∆x 2
 ∆x 2
 
Solve this linear implicit pressure equation for pressure at iteration, ν, to
obtain:

⇒ ( Pin−+11 ) , ( Pin +1 ) and ( Pin++11 )

ν ν ν

SATURATION EQUN. Now update the saturation equation as if it were explicit by taking
ν ν
(i) the oil mobility terms at the latest iteration, ν ⇒ λ 0 (So )i λ 0 (So )i +1 , etc.
(ii) the latest values of the pressures just calculated implicity ⇒ ( Pin−+11 ) , ( Pin +1 ) and ( Pin++11 )
ν ν ν

Update saturation explicity as follows:

(
∆t  λ o (So ) i +1/ 2
 ν
) (
  λ (Soν ) ) 
Soin+1 = Soin + 
φ  ∆x 2
(
(Pin++11 ) − (Pin +1 )
ν ν
) − o
  ∆x 2
i −1 / 2
(
(Pin +1 ) − (Pin−+11 )
ν ν
)  

   

⇒ (Sin +1 ) , (Sin++11 ) etc.

ν ν
to obtain latest iteration of saturations

ν = ν +1

Calculate the error, Err., by comparing the change in pressure or saturations

(or both as here) at the current (ν) and last (ν-1) iterations:

NX NX
Err. = ∑ ( Pin +1 ) − ( Pin +1 ) + ∑ (Sin +1 ) − (Sin +1 )
ν ν −1 ν ν −1

i =1 i =1
Figure 13.
Note that the Err. term above would have some scaling factors weighing the IMPES Algorithm for the
pressure and saturation terms because of the units of pressure.
numerical solution of the
two-phase pressure and sat-
No - continue iteration Yes
uration equations (IMPES
Is Err. < Tol ? Stop
⇒ IMplicit in Pressure,
Explicit in Saturation)

32
Numerical Methods in Reservoir Simulation
6
5 THE NUMERICAL SOLUTION OF LINEAR EQUATIONS

5.1 Introduction to Linear Equations

In all implicit finite difference methods, we end up with a set of linear equations to
solve. In fact, we have shown above that the equations involved - e.g. the discretised
pressure equation - may be a set of non-linear algebraic equations. However, these
can usually be linearised (say, by time lagging the coefficients as discussed above)
in order to obtain a set of linear equations. The solution to the original non-linear
equations may then be found by repeating or iterating through the cycle of solving
the linear equations. In the following section, we will address the issue of solving
the non-linear equations which arise in the discretised two-phase flow equations
directly.

A generalised set of linear equations may be written in full as follows:

a11 . x1 + a12 . x2 + a13 . x3 + a14 . x4 + a15 . x5 ... + a1n . xn = b1

a21 . x1 + a22 . x2 + a23 . x3 + a24 . x4 + a25 . x5 ... + a2 n . xn = b2

a31 . x1 + a32 . x2 + a33 . x3 + a34 . x4 + a35 . x5 ... + a3n . xn = b3

..... ...
..... ...
..... ...
an1 . x1 + an 2 . x2 + an 3 . x3 + an 4 . x4 + an 5 . x5 ... + ann . xn = bn (55)

where the aij and bi are known constants and the xi (i = 1, 2, ... , n) are the unknown
quantities which we are trying to find. The xi in our pressure equation, for example,
would be the unknown pressures at the next time step.

The set of linear equations can be written in matrix form as follows:

a11 a12 a13 a14 a15 ... a1n   x1  b1 

    
    
a21 a22 a23 a24 a25 ... a2 n   x2  b2 
    
    
a31 a32 a33 a34 a35 ... a3n   x3  = b3 
    
..... ...  .  . 
..... ...  .  . 
    
..... ...  .  . 
 
an1 an 2 an 3 an 4 an 5 ... ann   xn  bn  (56)

Institute of Petroleum Engineering, Heriot-Watt University 33

We can write the matrix A and the column vectors x and b as follows:

a11 a12 a13 a14 a15 ... a1n   x1  b1 

     
     
a21 a22 a23 a24 a25 ... a2 n   x2  b2 
     
     
A = a31

a32 a33 a34 a35 ... a3n  ;

x =  x3 ; and b = b3 
 
..... ...  .  . 
..... ...  .  . 
     
..... ...  .  . 
  x  b 
an1 an 2 an 3 an 4 an 5 ... ann   n  n
(57)

and the set of linear equations can be written in very compact form as follows:

A.x = b (58)

where A in an n x n square matrix and x and b are n x 1 column vectors.

A very simple example of a set of linear equations is given below:

4. x1 + 2. x2 + 1. x3 = 13

1. x1 + 3. x2 + 3. x3 = 14

2. x1 + 1. x2 + 2. x3 = 11 (59)

How do we solve these? In the above case, you can do a simple trial and error
solution to find that x1 = 2, x2 = 1 and x3 = 3, although this would be virtually
impossible if there were hundreds of equations. Remember, there is one equation for
every grid block in a linearised pressure equation in reservoir simulation (although
there are lots of zeros in the A matrix). This implies that we need a clear
numerical algorithm that a computer can work through and solve the linear
equations. The overview of approaches to solving these equations is discussed
in the next section.

5.2 General Methods for Solving Linear Equations

Firstly, let us note that there is a vast literature on solving sets of linear equations
and many books on the underlying theory and on the numerical techniques have
appeared. Indeed, petroleum reservoir simulation has led to the development of
some of these numerical techniques. We will not cover much of this huge subject
but we will cover enough that the student can appreciate - rather than understand in
detail - how the linear equation are solved in a simulator.

34
Numerical Methods in Reservoir Simulation
6
Basically, there are two main approaches for solving sets of linear equations involving
direct methods and iterative methods as follows:

(i) Direct Methods: this group of methods involves following a specific algorithm
and taking a fixed number of steps to get to the answer (i.e. the numerical values
of x1, x2 .. xn). Usually, this involves a set of forward elimination steps in order to
get the equations into a particularly suitable form for solution (see below), followed
by a back substitution set of steps which give us the answer. An example of such
a direct method is Gaussian Elimination.

(ii) Iterative Methods: in this type of method, we usually start off with a first guess
(or estimate) to the solution vector, say x(o) . Where we take this as iteration zero,
v = 0. We then have a procedure - an algorithm - for successively improving this
guess by iteration to obtain, x(1) → x(2) → x(3) .... → x(ν) etc. If it is successful, then
this iterative method should converge to the correct x as ν increases. The solution
should get closer and closer to the “correct” answer but, in many cases, we cannot
say exactly how quickly it will get there. Therefore, iterative methods do not
have a fixed number of steps in them as do direct methods. Examples of iterative
methods are the Jacobi iteration, the LSOR (Line Successive Over Relaxation)
method, etc.

The following two sections will discuss each of these approaches for solving
sets of linear equations in turn.

5.3 Direct Methods for Solving Linear Equations

In fact, we will not present the details of any direct method for solving linear equations.
Instead, we will discuss a general outline of how these methods work. Starting with
the basic linear equation in matrix form:

A.x = b (60)

where A in an n x n square matrix and x and b are n x 1 column vectors. We can

write the n x (n+1) augmented matrix of A and b coefficients as follows:

a11 a12 a13 a14 a15 ... a1n b1 

 
 
a21 a22 a23 a24 a25 ... a2 n b2 
 
 
a31 a32 a33 a34 a35 ... a3n b3 
 
..... ... ...
..... ... ...
 
..... ... ...
 
an1 an 2 an 3 an 4 an 5 ... ann bn  (61)

Remember that any mathematical operation we perform on one of the linear equations
(e.g. multiplying through by x2) must be performed on both sides of the equation i.e.

Institute of Petroleum Engineering, Heriot-Watt University 35

on the aij and the bi. Therefore, suppose we can transform the above augmented matrix
into the following form (we do not say how this is done, just that it can be done):

c11 c12 c13 c14 c15 ... c1n e1 

 
 
0 c22 c23 c24 c25 ... c2 n e2 
 
 
0 0 c33 c34 c35 ... c3n e3 
 
..... ... ...
..... ... ...
 
..... ... ...
0 0 0 0 ... 0 cnn en 
 (62)

The C-matrix above is in upper triangular form i.e. only the diagonals and above
are non-zero (cij ≠ 0, if j ≥ i) and all elements below the diagonal are zero (cij = 0,
if i > j), as shown above. Now consider why this particular form is of interest in
solving the original equations. The reason is that, if we have formed the augmented
matrix correctly (i.e. doing the same operations to the A and b coefficients), then the
following matrix equation is equivalent to the original matrix equation:

i.e. A.x = b <=> C.x = e (63)

Therefore,

c11 c12 c13 c14 c15 ... c1n   x1  e1 

    
    
0 c22 c23 c24 c25 ... c2 n   x2  e2 
    
    
0 0 c33 c34 c35 ... c3n   x3  = e3 
    
..... ..  .  . 
..... ..  .  . 
    
..... ..  .  . 
0 0 0 0 ... 0 cnn   xn  en 
 (64)

This matrix equation is very easy to solve, as can be seen by writing it out
in full as follows:

36
Numerical Methods in Reservoir Simulation
6
c11 . x1 + c12 . x2 + c13 . x3 + c14 . x4 + c15 . x5 ... + c1n . xn = e1

c22 . x2 + c23 . x3 + c24 . x4 + c25 . x5 ... + c2 n . xn = e2

c33 . x3 + c34 . x4 + c35 . x5 ... + c3n . xn = e3

c44 . x4 + c45 . x5 ... + c4 n . xn = e4

...... ...
...... ...

cn −1, n −1 . xn −1 + cn −1, n . xn = en −1
cnn . xn = en
(65)

We can easily solve this equation by back-substitution, starting from the equation n
which only has one term to find xn, this xn is then used in the (n-1) equation (which
has 2 terms) to find xn-1 etc. as follows:

cnn . xn = en => xn = en / cnn now use xn in ..

 e − c .x 
cn −1, n −1 . xn −1 + cn −1, n . xn = en −1 => xn −1 =  n −1 n −1, n n  use xn −1 , xn in ..
 cn −1, n −1 

cn − 2, n − 2 . xn − 2 + cn − 2, n −1 . xn −1 + cn − 2, n . xn = en − 2
e - cn − 2, n −1 . xn −1 − cn − 2, n . xn 
=> xn − 2 =  n − 2  etc.
 cn − 2, n − 2 

(66)

Hence, by working back through these equations in upper triangular form, we can
calculate xn, xn-1, , xn-2.... back to x1.

Institute of Petroleum Engineering, Heriot-Watt University 37

EXERCISE 4.

Solve
4 x -the1following
x − 2 x simple
+ 1example
x + 1ofx an
= upper
5 triangular matrix:
1 2 3 4 5
4 x1 - 1x2 − 2 x3 + 1 x 4 + 1x 5 = 5
2 x 2 − 2 x3 + 1x4 + 2 x5 = 12
2 x2 − 2 x3 + 1x4 + 2 x5 = 12
3 x3 + 4 x4 + 1x5 = 30
3 x3 + 4 x4 + 1x5 = 30
2 x4 − 1x 5 = 3
2 x 4 − 1x 5 = 3
3 x5 = 15
3 x5 = 15

Answer
Answer:
Answer
x1 = ..........; x2 = ..........; x3 = ..........; x4 = ..........; x5 = ..........
x1 = ..........; x2 = ..........; x3 = ..........; x4 = ..........; x5 = ..........

5.4 Iterative Methods for Solving Linear Equations

As noted above, the idea in an iterative method for solving a set of linear equations is
to make a first guess at the solution and then to refine it in a stepwise manner using a
suitable algorithm. This procedure should gradually converge to the correct answer.
We illustrate the idea of such a method using a very simple iterative scheme. Note
that this simple point iterative scheme will work for some of the examples we
try here but it is not one that is recommended for use in reservoir simulation.
However, it adequately illustrates the main ideas which you need to know for the
purposes of this part of the course.

Our simple scheme starts with the longhand version of a set of linear equations and,
for our purposes, we will just take a set of four linear equations as follows:

a11 . x1 + a12 . x2 + a13 . x3 + a14 . x4 = b1

a21 . x1 + a22 . x2 + a23 . x3 + a24 . x4 = b2

a31 . x1 + a32 . x2 + a33 . x3 + a34 . x4 = b3

a41 . x1 + a42 . x2 + a43 . x3 + a44 . x4 = b4 (67)

We can rearrange the above set of equations as follows:

38
Numerical Methods in Reservoir Simulation
6
1
x1 =
a11
[
b1 − ( a12 x2 + a13 x3 + a14 x4 )]

1
x2 =
a22
[
b2 − ( a21 x1 + a23 x3 + a24 x4 ) ]

1
x3 =
a33
[
b3 − ( a31 x1 + a32 x2 + a34 x4 ) ]

1
x4 =
a44
[
b4 − ( a41 x1 + a42 x2 + a43 x3 ) ]
(68)

The resulting equations are precisely equivalent to the original set although, if
anything, they look a bit more complicated. Why would we deliberately complicate
the situation? In fact, the above reorganised equations forms the basis for an
iterative scheme, which we can use to solve the equations numerically. Firstly, we
need to introduce some notation as follows:

Notation:
(ν )
xi
xi(v) - the solution for xi at iteration ν; where ν is the iteration
counter.

x ( 0 ) , x (ν ) - the first guess and νth iteration of the solution vector, x.

Err. - an estimate of the error in the iteration scheme from the ν

to the (ν+1) iteration

Tol - some “small” quantity which determines the acceptable

error in the iteration scheme; if Err. < Tol, then the
scheme has converged.

Using the above notation in equations 68 above, gives the following simple iteration
scheme:

Institute of Petroleum Engineering, Heriot-Watt University 39

x1(ν +1) =
1
a11
[ ]
b1 − ( a12 . x2(ν ) + a13 . x3(ν ) + a14 . x4(ν ) )

x2(ν +1) =
1
a22
[
b2 − ( a21 . x2(ν ) + a23 . x3(ν ) + a24 . x4(ν ) ) ]

x3(ν +1) =
1
a33
[
b3 − ( a31 . x1(ν ) + a32 . x2(ν ) + a34 . x4(ν ) ) ]

x4(ν +1) =
1
a44
[
b3 − ( a41 . x1(ν ) + a42 . x2(ν ) + a43 . x3(ν ) ) ] (69)

This iteration scheme may now be applied as follows:

(i) Make an initial guess at the solution, iteration

ν = 0: x ( 0 ) = x1( 0 ) , x2( 0 ) , x3( 0 ) , x4( 0 )

(ii) Update the solution to the next iteration ν+1 using equations 69
4
Err. = ∑ xik +1 − x ki
(iii) Estimate an Error term (Err.) by comparing the latest with the previous
i =1
iteration of the unknowns, e.g.:
4
Err. = ∑ xiν +1 − xiν
i =1

(iv) Is Err. < Tol.

If yes - the scheme has converged; if no - go back to step (ii) and continue
the iterations.

This scheme is now illustrated with a practical example.

Example: Solve the following set of equations using the iterative scheme above.

3.1x1 - 0.32 x2 + 0.5 x3 + 0.1x4 = 13.92

0.2 x1 + 2.1x2 + 0.33 x3 + 0.21x4 = 5.63

0.23 x1 - 0.32 x2 + 4.0 x3 + 0.3 x4 = 15.19

0.42 x1 + 0.22 x2 + 0. 5 x 3 + 5.2 x4 = 16.76

40
Numerical Methods in Reservoir Simulation
6
(0) (0) (0) (0)
Take the first guess: x1 = 2, x2 = 2, x3 = 2, x4 = 2

Hint: reorganise the above equations as follows:

(
x1(ν +1) = (1 / 3.1) * 13.92 - (- 0.32 x2(ν ) + 0.5 x3(ν ) + 0.1x4(ν ) ) )

(
x2(ν +1) = (1 / 2.1) * 5.63 - (0.20 x ( ν + 1)
1 + 0.33 x3(ν ) + 0.21x4(ν ) ))

(
x3(ν +1) = (1 / 4.0) * 15.19 - (0.23x (ν +1)
1 - 0.32 x2(ν +1) + 0.30 x4(ν ) ))

(
x4(ν +1) = (1 / 5.2) * 16.76 - ( 0.42 x (ν +1)
1 + 0.22 x2(ν +1) + 0.5 x3(ν +1) ))

EXERCISE 5.

Fill in the table below for 10 iterations using a calculator or a spreadsheet:

POINT ITERATIVE SOLUTION OF LINEAR EQUATIONS

Iteration
counter x at iter. ν
ν x1 x2 x3 x4
First guess = 0 2.0 2.0 2.0 2.0
1
2
3
4
5
6
7
8
9
10

Iteration
counter x at iter. ν
ν x1 x2 x3 x4
In some First
cases, it may
guess = 0be possible
2.0 to develop
2.0 improved
2.0 iteration 2.0 schemes by using
the latest information 1that is available. For
4.309677 1.97619 example, in
3.6925the above iteration scheme,
2.784615
we could use the very latest value of x1
2 4.008926 (ν+1)
when we are calculating x2(ν+1) since we
1.411795 3.498943 2.436331
3 3.993119 1.505683
already have this quantity. Likewise, we can use both 3.497206
x1(ν+1) 2.503112
and x2(ν+1) when we
4 4.000937 1.500783 3.500617 2.500584
calculate x3 (ν+1)
etc as shown below:
5 3.999963 1.499755 3.499965 2.499832
6 3.999986 1.500026 3.499995 2.500017
7 4.000003 1.5 3.500002 2.500001
8 4 1.499999 3.5 2.5
Converged ⇒ 9 4 1.5 3.5 2.5
10 4 1.5 3.5 2.5

Note - converges fully at k = 9 iterations.

Institute of Petroleum Engineering, Heriot-Watt University 41
(
x1(ν +1) = (1 / 3.1) * 13.92 - (- 0.32 x2(ν ) + 0.5 x3(ν ) + 0.1x4(ν ) ) )

(
x2(ν +1) = (1 / 2.1) * 5.63 - (0.20 x ( ν + 1)
1 + 0.33 x3(ν ) + 0.21x4(ν ) ))

(
x3(ν +1) = (1 / 4.0) * 15.19 - (0.23x (ν +1)
1 - 0.32 x2(ν +1) + 0.30 x4(ν ) ))

(
x4(ν +1) = (1 / 5.2) * 16.76 - ( 0.42 x (ν +1)
1 + 0.22 x2(ν +1) + 0.5 x3(ν +1) ))
x1(ν +1) , x2(ν +1) , x3(ν +1) => underlined terms, imply they are at the latest time
available).

When this is done, we find the following results:

IMPROVED POINT ITERATIVE SOLUTION OF LINEAR EQUATIONS (Using
latest information)

- see spreadsheet CHAP6Ex5.xls - Sheet 2

Iteration
counter x at iter. ν
ν x1 x2 x3 x4
First guess = 0 2 2 2 2
1 4.309677 1.756221 3.540191 2.460283
2 4.021247 1.495632 3.501408 2.498333
3 3.999376 1.500005 3.500161 2.500035
4 3.999973 1.499974 3.499997 2.500004
5 3.999998 1.5 3.5 2.5
Converged ⇒ 6 4 1.5 3.5 2.5
7 4 1.5 3.5 2.5
8 4 1.5 3.5 2.5
9 4 1.5 3.5 2.5
10 4 1.5 3.5 2.5

Note - converges fully at ν = 6 iterations.

Clearly, comparing this table with the previous one using the simple point iterative
method, we see that this method does indeed converge quicker. Although this
is just one example, it is generally true that using the latest information in an
iterative scheme improves convergence.

Notes on iterative schemes: there are several point to note about iterative
schemes, which we have not really demonstrated or explained here. However,
you can confirm some of these points by using the supplied spreadsheets.
These points are:

(i) Iterative solution schemes for linear equations are often relatively simply to apply
- and to program on a computer (usually in FORTRAN). If you study the supplied
spreadsheets (CHAP6Ex5.xls), you can see that they are quite simple in structure
for this set of equations.

42
Numerical Methods in Reservoir Simulation
6
(ii) The convergence rate of an iterative scheme may depend on how good the initial
guess (x(0)) is. If you want to demonstrate this for yourself, run the spreadsheet
(CHAP6Ex5.xls) with a more remote initial guess. Here is the same example as
that presented above with a completely absurd initial guess:

POINT ITERATIVE SOLUTION OF LINEAR EQUATIONS - Bad First Guess

(you can use CHAP6Ex5.xls to confirm these results)

Iteration
counter x at iter. ν
ν x1 x2 x3 x4
First guess = 0 900 -2000 456 23000
1 -1017.45 -2454.69 -1932.95 -28.7
2 63.79526 406.2002 -131.922 375.1144
3 55.59796 -20.1756 4.491703 -6.43019
4 1.890636 -2.67691 -0.53117 0.84584
5 4.326953 2.668945 3.538073 3.234699
6 4.090824 1.389409 3.519613 2.420476
7 3.987986 1.49622 3.491895 2.495457
8 4.001064 1.502872 3.500729 2.50191
9 4.000117 1.499593 3.500025 2.499722
10 3.999963 1.500013 3.499982 2.500005
11 4.000004 1.500006 3.500003 2.500004
12 4 1.499999 3.5 2.499999
Converged ⇒ 13 4 1.5 3.5 2.5
14 4 1.5 3.5 2.5
15 4 1.5 3.5 2.5

(iii) We cannot tell in advance how many iterations may be required in order to
converge a given iterative scheme. Clearly, from the results above, a good initial
guess helps. In reservoir simulation, when we solve the pressure equation, a good
enough guess of the new pressure is the values of the old pressures at the last
time step. If nothing radical has changed in the reservoir, then this may be
fine. However, if new wells have started up in the model or existing wells have
changed rate very significantly, then the “pressure at last time step” guess may
not be very good. However, with a very robust numerical method, convergence
should still be achieved.

(iv) It helps in an iterative method to use the latest information available. This was
demonstrated in the iterative schemes in spreadsheet CHAP6Ex5.xls.

5.5 A Comparison of Iterative and Direct Methods for Solving Linear

Equations
In this Chapter, we have described two general methods - direct and iterative - for
solving the linear equations which arise when we discretise the flow equations of
reservoir simulation. We have not indicated which of these methods is usually “best”
for reservoir simulation problems. This is because, it depends on the problem. The
final numerical problem which is solved in a reservoir simulator could be quite small
and simple, or it could be very large and have some intrinsically difficult numerical
problems within it.

Institute of Petroleum Engineering, Heriot-Watt University 43

In any numerical computational method, such as those for solving a set of linear
equations, we need to have some way of calculating the amount of “work” involved.
To some extent this depends on the computer architecture since new generation parallel
processing machines are becoming available as discussed in Chapter 1. However, we
will take quite a simple view based on a conventional serial processing machine and
will define computational work as simply the number of multiplications and divisions
(addition and subtraction is cheap!). In this way, we can define the amount of work
required for any given direct and iterative scheme for solving linear equations. We
present results without any proof for two such schemes called the BAND (a direct
method) and LSOR (Line Successive Over-Relaxation; an iterative method) schemes.
For a 2D problem with NX and NY grid blocks in the x and y directions, respectively
(where we choose NY < NX)

BAND (direct), Work, WB ≈ NX.NY3 (70)

LSOR (iterative), Work, WLSOR ≈ NX.NY. Niter (71)

where Niter is the number of iterations until convergence is reached. It is already

quite clear that the amount of work for the direct method - which only has to be done
once - is larger than that for a single iteration of the iterative method. However, we
do not know in advance the number of iterations, Niter. We illustrate some typical
results for these two methods with a simple exercise below.

EXERCISE 6.

Which is best method (i.e. that requiring the lowest amount of computational work),
BAND or LSOR, for the following problems?

(i) NX = 5, NY = 3, Niter = 50 (a small problem)

(ii) NX = 20, NY = 5, Niter = 50

(iii) NX = 100, NY = 20, Niter = 70

(iv) NX = 400, NY = 100, Niter = 150

NX NY Niter BAND LSOR, Comment

WB WLSOR
5 3 50
20 5 50
100 20 70
400 100 150

44
Numerical Methods in Reservoir Simulation
6
We can summarise our comparison between direct and iterative methods for solving
the sets of linear equations that arise in reservoir simulation as follows:

(i) A direct method for solving a set of linear equations has an algorithm that
involves a fixed number of steps for a given size of problem. Given that the
equations are properly behaved (i.e. the problem has a stable solution), the direct
method is guaranteed to get to the solution in this fixed number of steps. No first
guess is required for a direct method.

(ii) An iterative method on the other hand, starts from a first guess at the solution (x(o))
and then applied a (usually simpler) algorithm to get better and better approximations
to the true solution of the linear equations. If successful, the method will converge
in a certain number of iterations, Niter, which we hope will be as small as possible.
However, we cannot usually tell what this number will be in advance. Also, in some
cases the iterative method may not converge for certain types of “difficult” problem.
We may need to have good first guess to make our iterative method fast. Also, it
often helps to use the latest computed information that is available (see example
above).

(iii) Usually the amount of “work” required for a direct method is smaller for
smaller problems but iterative methods usually win out for larger problems. For
an iterative method, the amount of work per iteration is usually relatively small
but the number of iterations (Niter) required to reach convergence may be large and
is usually unknown in advance.

Institute of Petroleum Engineering, Heriot-Watt University 45

EXERCISE 7.

The linear equations which arise in reservoir simulation may be solved by a direct
solution method or an iterative solution method. Fill in the table below:

Direct solution method Iterative solution method

Give a very
brief description
of each method

Main advantages 1. 1.

2. 2.

Main 1. 1.
disadvantages

2. 2.

46
Numerical Methods in Reservoir Simulation
6

6. DIRECT SOLUTION OF THE NON-LINEAR EQUATIONS OF

MULTI-PHASE FLOW

6.1 Introduction to Sets of Non-linear Equations

We noted in Section 4 above that the equations which we obtain when we discretise
the two-phase flow equations are actually non-linear in nature. However, the strategy
we discussed above (the IMPES method) involved tackling the problem almost as
if it were a linear set of equations since we time-lagged the coefficients to reduce
the problem to a linear set of equations. Then we repeated this process - we applied
repeated iterations - until it (hopefully) converged. Therefore, we took a non-linear
problem and solved it as if it were a series of linear problems.

In this section, we will introduce the general idea of solving sets of non-linear equations
numerically and indicate how this can be applied to the two-phase flow equations.
We will do this in a simplified manner that shows the basic principles without going
into too much detail. Firstly, compare the difference between solving the following
two sets of two equations, a set of 2 linear equations:

2 x1 + x2 = 10
3 x1 − x2 = 5 (72)

and a set of 2 non-linear equations:

x1 + 2 x2 .e − x1 = 2
x2 .( x1 )2 − x1 .( x2 ) = −5
2
(73)

It is immediately obvious that the second set of non-linear equations is more difficult.
The first set of linear equations can be rearranged easily to show that x1 = 3 and x2 =
4. However, it is not as straightforward to do the same thing for the non-linear set
of equations. As a first attempt to solve these, we might try to develop an iterative
scheme by rearranging the equations as follows
(ν )
x1(ν +1) = 2 − 2 x2(ν ) .e − x1

(ν +1)
x2(ν ) .( x1(ν ) )2 +5
x =
x1(ν )
2
(74)

where, as before, ν denotes the iteration counter.

The actual solution in this case is, x1 =2.51068, x2 = 3.144089. Applying the above
iterative scheme with a first guess, x1(0) = x2(0) = 1.0, gives the results shown in Table
2 below (where we have removed some of the intermediate iterations).

Table 2: Non-linear scheme applied to solution of equations 73 using the point iterative
scheme of equation 74 (spreadsheet, CHAP6Ex5.xls - Sheet 3)

Institute of Petroleum Engineering, Heriot-Watt University 47

Iteration
counter x at iter. ν
ν x1 x2
First guess = 0 1 1
1 2.735759 2.44949
2 2.317673 2.920421
3 2.575338 2.987627
4 2.454885 3.104133
.... .... ....
9 2.513337 3.141814
10 2.508958 3.144173
.... .... ....
25 2.510682 3.144089
26 2.510681 3.144089

You can check the results in Table 2 or experiment with other first guesses using
Sheet 3 of spreadsheet CHAP6Ex5.xls. We note that the convergence rate in Table
2 is not very fast but, in this case, it does reach a solution. In general, it is usually
very difficult to establish for certain whether a given scheme will converge for non-
linear equations although there is a large body of mathematics associated with the
solution of such systems (which is beyond the scope of this course).

6.2 Newton’s Method for Solving Sets of Non-linear Equations

We take a step back to the solution of a single non-linear equation such as:

x 2 .e − x = 0.30 (75)

which has the solution x = 0.829069 (you can verify this by calculating x 2 .e − x
and making sure it is 0.30 - to an accuracy of ~ 4.6x10-8). This equation can be
represented as:

f ( x ) = 0 where f ( x ) = x 2 .e − x − 0.30 (76)

and the solution we require is the value of x for which the function f(x) is zero.
Expand f(x) as a Taylor series as follows:

δx 2 ' '
f ( x ) ≈ f ( xo ) + δx. f ' ( x0 ) + . f ( x0 ) + ...
2 (77)

Thinking in terms of an iteration scheme (x(ν) → x(ν+1)) as a basis for calculating

better and better guesses of the solution of f(x) = 0, we can rewrite the Taylor series
above as:

f ( x ) ≈ f ( x (ν ) ) + ( x (ν +1) − x (ν ) ). f ' ( x (ν ) ) (78)

where we have neglected all the second order and higher terms. Since we require
the solution for f(x) = 0, then we obtain:

f ( x (ν ) ) + ( x (ν +1) − x (ν ) ). f ' ( x (ν ) ) = 0 (79)

48
Numerical Methods in Reservoir Simulation
6
which can be rearranged to the following algorithm to estimate our updated guess,
x(ν+1) as follows:

f ( x (ν ) )
x (ν +1) = x (ν ) −
f ' ( x (ν ) ) (80)

This equation is the basis of the Newton-Raphson algorithm for obtaining better and
better estimates of the solution of the equation, f(x) = 0. Note that we need both a first
guess, x(0), and also an expression for the derivative f ' (x(ν)) at iteration ν. In the simple
example in equation 76, we can obtain the derivative analytically as follows:

df ( x )
= f ' ( x ) = 2 x.e − x − x 2 .e − x
dx (81)

Therefore the Newton-Raphson algorithm for solving equation 75 above is as

follows:

 ( x (ν ) )2 .e − x (ν ) − 0.30 
x ( ν + 1)
=x (ν )
-  (ν ) − x ( ν ) 
− ( x (ν ) ) .e − x 
2 (ν )
 2 x .e
(82)

The point iterative method of the previous section and the Newton-Raphson method
of equation 82 have both been applied to the solution of equation 75 (see Sheet 4
of spreadsheet CHAP6Ex5.xls for details). The results are shown in Table 3 for a
first guess of x(0) = 1.

Table 3: Comparison of the point iterative and Newton-Raphson methods for solving
non-linear equation 75; (see spreadsheet CHAP6Ex5.xls - Sheet 4 for details)

Iteration Value of x at iteration ν

Number

ν Point Newton-
Iteration Raphson
method method
x(ν) x(ν)

0 1 1
1 0.903042 0.815485
2 0.860307 0.825272
3 0.84212 0.82806
4 0.834497 0.828805
5 0.831322 0.829
6 0.830003 0.829051
7 0.829456 0.829064
8 0.82923 0.829068
9 0.829136 0.829069
10 0.829097 0.829069
11 0.82908 0.829069
12 0.829074 0.829069
13 0.829071 0.829069
14 0.82907 0.829069
15 0.829069 0.829069
16 0.829069 0.829069

Institute of Petroleum Engineering, Heriot-Watt University 49

The correct solution to equation 75 is x = 0.829069. The results in Table 3 show
that the Newton-Raphson method converges to this solution (to an accuracy better
than 10-6) in 9 iterations, whereas it takes 15 iteration for the point iterative method
to converge at the same level accuracy. This performance can vary quite a bit from
problem to problem but the Newton-Raphson method is generally better if the derivative
term, f '(x(ν)), is not too close to zero. When this derivative gets too small, it can be
seen that the second term in equation 80 would start to get very large or “blow up”,
as it is sometimes described. We will not go into details about the convergence
properties of the Newton-Raphson but it is meant to be “quadratic” meaning that the
error should decrease quite rapidly.

We now go on to see how the Newton-Raphson method can be applied to sets of

non-linear equations. Returning to the set of two non-linear equations in the previous
section, we note that another way of writing this set of equations in a general way
is as follows:

F1 ( x1 , x2 ) = 0

F2 ( x1 , x2 ) = 0 (83)

where, in the case of equation 73 above, these functions would be given by:

F1 ( x1 , x2 ) = x1 + 2 x2 .e − x1 − 2

F2 ( x1 , x2 ) = x2 .( x1 )2 − x1 .( x2 ) + 5
2
(84)

The problem is then to find the values of x1 and x2 that make F1 = F2 = 0. In this
section, we present Newton’s method for the solution of sets of non-linear equations.
Basically, we will simply state the Newton-Raphson algorithm without proof (although
it will resemble the simple form of the Newton-Raphson method above) for sets of
non-linear equations. We will then illustrate what this means for the example of the
set of two non-linear equations above (equation 73).

Definitions:

As before,
(ν )
x (νx) and andx (νx+(ν1)+1=) the solution vectors at iteration levels ν and (ν+1)

New terms are:

xx((ν ) and
F(F x) ) ) (ν +1)
(ν )
(ν
x
N = the number of non-linear equations (and hence
unknowns, x1, x2, .... xN)
νx)(ν )
x (F ( x (ν ) ) = the vector of function values, F1, F2 .. FN at x(ν)

F(Fx((ν)x) =) =0 0forforthe
the"true"
"true"solution
solutionofofthe
thesetsetofofnon
non - linearequations
- linear equations
x

50
F( x ) = 0 for the "true" solution of the set of non - linear equations
Numerical Methods in Reservoir Simulation
6
That is:
 F1 ( x1ν , x2ν , x3ν ,....., x νN ) 
 
 F2 ( x1ν , x2ν , x3ν ,....., x νN ) 
 
 F3 ( x1ν , x2ν , x3ν ,....., x νN ) 
F( x ) = 
(ν )

 ..... 
 ..... 
 
 F ( x ν , x ν , x ν ,....., x ν )
 N 1 2 3 N 
(85)

J ( x (ν ) ) = the NxN Jacobian matrix defined as follows:

 ∂F1 ( x (ν ) ) ∂F1 ( x (ν ) ) ∂F1 ( x (ν ) ) ∂F1 ( x (ν ) ) 

 ...... 
 ∂x1 ∂x2 ∂x3 ∂x N 
 
 
 ∂F2 ( x (ν ) ) ∂F2 ( x (ν ) ) ∂F2 ( x (ν ) ) ∂F2 ( x (ν ) ) 
 ...... 
 ∂x1 ∂x2 ∂x3 ∂x N 
 
J ( x (ν ) ) =  
 ∂F3 ( x (ν ) ) ∂F3 ( x (ν ) ) ∂F3 ( x (ν ) ) ∂F3 ( x (ν ) ) 
 ...... 
 ∂ x 1 ∂x 2 ∂ x 3 ∂x N 
 ..... 
 
 ..... 
 (ν ) 
 ∂FN ( x ) ∂FN ( x ) ∂FN ( x ) ...... ∂FN ( x ) 
(ν ) (ν ) (ν )

 ∂x1 ∂x2 ∂x3 ∂x N  (86)

[ J ( x )]
−1
(ν )
= the inverse of the Jacobian matrix. (Recall that the inverse of a
matrix A is denoted by A-1 and by definition, A-1 .A = I, where I is the identity matrix
- all diagonals 1 and all other elements zero - see below). The above matrix is also
an NxN matrix with the property:
[ J ( x )] .[ J ( x )] = I, the identity matrix
−1
(ν ) (ν )

[ J ( x )] .[ J ( x )] = I, the identity matrix

−1
(ν ) (ν )

1 0 0 0 0 ...... 0 
10 0 1 0 0 0 0 0 0 ...... 0 0 
......
0 1 0 0 0 ...... 
 0 0 1 0 0 ...... 0 0 
00 0 0 1 0 0 1 0 0 ...... 0 0 
......
[ ] [
That is → J ( x (ν ) ) . J ( x (ν ) ) = I = ]
−1
0 0 0 1 0  (87)
0 0 0 0 1 ...... 0 0 
......
[ ] [
J(x ) . J(x ) = I =   ]
−1
(ν ) (ν )
0 .......
0 0 0 1 ...... 0  
 
.......
....... 

....... 
0 0 0 0 0 ...... 1 
 
0 0 0 0 0 ...... 1 
Institute of Petroleum Engineering, Heriot-Watt University 51
Using all of the above definitions, the Newton-Raphson algorithm for a set of non-
linear equations, F( x ) = 0, is given by the following expression:

[ ]
x (ν+1) = x (ν) − J ( x (ν) ) .F( x (ν) )
−1

(88)

We first demonstrate how this is applied to a simple example (equations 73) before
going on to show how it is applied to the more complicated non-linear sets of
equations which arise in the fully implicit discretisation of the reservoir simulation
equations.

Example: Returning to the simple example, where N = 2 (rearranged equations

73):

F1 ( x1 , x2 ) = x1 + 2 x2 .e − x1 − 2

F2 ( x1 , x2 ) = x2 .( x1 )2 − x1 .( x2 ) + 5
2
(89)

In the notation developed above, these equations become:

 F1 ( x1 , x2 )   x1 + 2 x2 .e − x1 − 2 
F( x ) =   =  
 F2 ( x1 , x2 )  x2 .( x1 )2 − x1 .( x2 ) + 5
2
(90)

We note that the 2x2 Jacobian matrix is given by:

 ∂F1 ( x (ν ) ) ∂F1 ( x (ν ) ) 
 
 ∂x1 ∂x2 
J(x ) = 
(ν ) 
 
 ∂F2 ( x (ν ) ) ∂F2 ( x (ν ) ) 
 
 ∂x1 ∂x2 
(91)

where each of the elements can be evaluated analytically to obtain:

(
 1 − 2 x (ν ) .e − x1(ν )
2 ) (2.e )
− x1( ν ) 

J ( x (ν ) ) =  
 
(
 2 1 (2 )
 2 x (ν ) . x (ν ) − x (ν ) 2
) (( x ) − 2x
(ν ) 2
1
(ν )
1 . x2(ν ) ) 
 (92)

Hence, the Newton-Raphson method for this simple system becomes:

52
Numerical Methods in Reservoir Simulation
6
( ) (2.e )
−1
 1 − 2 x (ν ) .e − x1(ν ) − x1( ν )   x (ν ) + 2 x (ν ) .e − x1(ν ) − 2 
 2
  1 2 
x ( ν + 1) = x ( ν ) −    
   
(
 2 1
(2 )
 2 x (ν ) . x (ν ) − x (ν ) 2
) (( x ) − 2x
(ν ) 2
1
(ν )
1 . x2(ν ) ) 

 x2(ν ) .( x1(ν ) )2 − x1(ν ) .( x2(ν ) ) + 5
2

(93)

[ ( )]
−1
(ν )
which we could solve if we knew how to invert the Jacobian matrix to find J x
This is a detail which we don’t need to know for this course but it can be done.
Indeed, for a 2x2 matrix - such as in the above case - it is quite easy to work out the
inverse. The inverse of a simple 2 x 2 matrix A given by:

A =  a b
 
 c d

is well known to be

A-1 = 1  d − b
 
( ad − bc)  − c a 

This can easily be proven by multiplying out there matrices and showing that
A-1 .A = I. The factor (ad-bc) is known as the determinant of the matrix and it must be
non-zero for A-1 to exist. Rather than using equation above to write out the analytical
[ ]
form of the J ( x (ν ) ) matrix we first simply numerically evaluate the matrix J at a
[]
−1
given iteration and apply equation above no get J . This is done in the spreadsheet
Chap6Exxx.xls where results are shown.

[ J ( x )]
−1
(k )
= *** (94)

which allows us to apply the Newton-Raphson method directly to our simple example.
For details see Sheet 5 spreadsheet ***.xls. The results are shown in Table 6.xx.
(NOTES & SPEADSHEET STILL UNDER CONSTRUCTION)

Institute of Petroleum Engineering, Heriot-Watt University 53

6.3 Newton’s Method Applied to the Non-linear Equations of Two-
Phase Flow
We now return to the discretised equations of two phase flow which we noted at the
time, formed a set of non-linear equations. The simplified form of these equations
is as follows (see Section 4.1):

Pressure (equation 50):

(λ (S ))
T
n +1
o
i −1 / 2
Pn +1
−
(
 λ (Son +1 )
T ) i −1 / 2
+
(λ (S ))
T
n +1
o
i +1 / 2

 P n +1 + T o
(
λ (S n + 1 ) ) i +1 / 2
Pin++11 = 0
i −1
∆x 2
 ∆x 2 ∆x 2
 i ∆x 2
 

Saturation (Form B - equation 54):


(
∆t  λ o ( So )
n +1
)   λ ( Son +1 )( ) 
n +1  
Soin+1 − Soin − 
φ  ∆x 2
i +1 / 2
( i +1 i ) −  o ∆x 2
P n +1
− P n +1  i −1 / 2
(i
P n +1
− Pi −1 )  = 0

    

n +1
Given that the unknowns that we are trying to find are Pi , Sin +1 i = 1,2, NX. ,
then we can write the above equations in general form as:

FP, i ( S n +1 , P n +1 ) = 0
FS , i ( S n +1 , P n +1 ) = 0
(95)
FP, i ( S n +1 , P n +1 ) andn +F1S , i (nS+n1+1 , P n +1 ) n +1 n +1
where the two non-linear equations, FP, i ( S , P ) and FS , i ( S , P ) , arise from
the pressure (P) and saturation (S) equations, respectively, as given above. The
vectors of unknowns, S n +1 and P n +1, nare
+1 given nby
+1 :
S and P
n +1
S  1  P1n +1 
 n +1   n +1 
S2   P2 
 n +1   n +1 
S3   P3 
...  ... 
   
...  ... 
 S n +1   P n +1 
S n +1 =  i −1  and P n +1 =  i −1 
Sin +1   Pi n +1 
 n +1   n +1 
Si +1   Pi +1 
   
...  ... 
...  ... 
   
...  ... 
 S n +1   P n +1 
 NX   NX  (96)

54
Numerical Methods in Reservoir Simulation
6
FP, i ( S n +1 , P n +1 ) and FS , i ( S n +1 , P n +1 )

However, the given FP, i S , P ( n +1 n +1

)
and FS , i S , P ( n +1 n +1
)
at grid block i only depend
on the quantities Sin−+11 , Sin +1 , Sin++11 and Pi n−1+1 , Pi n +1 , Pi n+1+1 and, rather than on all of the
other saturations and pressures in the system, since these are the nearest neighbours
coupled together in the Sin−+11equations
, Sin +1 , Sinabove.
+1 n +1
+1 and Pi −1 , Pi
n +1
, Pi n+1+1

We may also write one total unknows vector Xn+1 using a combination of the saturation
and pressure vectors as follows:

 S1 n +1 
 n +1 
 P1 
 S n +1 
 2 
 P2 n +1 
 
X n +1 =  
 S n +1 
 i 
 Pi n +1 
 
 
 S n +1 
 NX 
 n +1 
 PNX 

and the equation to be solved is then F(Xn+1)=0

The saturations and pressures are coupled together to their nearest neighbours through
the discretisation equations (50 and 54 above). The general form of the solution using
the Newton iteration is then to take a starting guess (iteration, ν=0) Xn+1(0) and then
apply the formulation as above to obtain:

[ ]
−1
X n +1 (ν +1) = X n +1(ν ) + J ( X )(ν ) . F ( X (ν ) )

It is rather more complex to constuct the Jacobian matrix J ( X ) of derivatives but

[]
−1
it can be done. We also need to invert this matrice to obtain J in order to apply
the above algorithm. However, in practice, there are various methods that try to
[]
−1
construct the J matrix more directly, often in an approximate manner. Note that
the Jacobian matrix is very sparse since there are just nearest neighbour interaction
(as was the case for the matrices associated with single phase flow discussed earlier
in this Chapter). A consequence of this is that the inverse matrix is also quite sparse
and there are many numerical techniques available to solve such problems.

In this course, we will not give any more detail on the solution of the non-linear
equations which arise in reservoir simulation.

Institute of Petroleum Engineering, Heriot-Watt University 55

7 NUMERICAL DISPERSION - A MATHEMATICAL APPROACH

7.1 Introduction to the Problem

In Chapter 4, we discussed the physical idea of numerical dispersion, which is also
sometimes referred to as numerical diffusion. Numerical dispersion is essentially
the artificial spreading or “diffusion” of a front - e.g. a waterfront in a water/oil
displacement - due to the coarse grid used in the simulation. It is a numerical effect
and it can be reduced by taking a finer grid. We note that in real reservoirs there are
some real dispersive physical mechanisms which result in the spreading of fronts.
These arise due to the effects of capillary pressure and also from the interaction of the
fluid flow with small scale (cm - m) permeability heterogeneity of the reservoir rock.
Indeed, there are also some quite complex interactions between the capillary forces
and the small scale heterogeneity that also lead to types of physical spreading in the
reservoir. In an ideal calculation, we would take a sufficiently fine grid that the level
of physical (i.e. real) spreading or diffusion was correctly represented by our grid;
i.e. the numerical diffusion would be less that the physical diffusion. This is almost
never possible in a field scale simulation, although it can be achieved in modelling
laboratory scale experiments on flows through small rock samples or bead packs.
For such a fine scale simulation, the level of numerical diffusion can realistically be
made much less than that from physical sources.

In this section, we return to the issue of numerical dispersion - a term which we will
now use interchangeably with “numerical diffusion”. Indeed, “diffusion” is often
referred to as quite a specific physical effect which is described by certain well-known
equations which are, not surprisingly, known as diffusion equations. For example,
the simplified pressure equation derived in Chapter 5, equation 27, is an example of
a diffusion equation. This equation has the form:

 ∂P   ∂2 P 
  = Dh  2 
 ∂t   ∂x  (97)

where Dh is the hydraulic diffusivity (Dh = k/( cf .φ.μ)) and this is a standard form
of the classical diffusion equation. We now consider how the effects of a grid can
lead to “diffusion-like” terms when we try to solve the flow equations numerically.
 ∂2 P 
In the above equations, it is specifically the  2  term that is the “dispersive” or
“diffusive” part.  ∂x 

7.2 Mathematical Derivation of Numerical Dispersion

In order to show mathematically how a “diffusive” term arises when we solve certain
transport equations numerically, we use a simple form of the two-phase saturation
equation. In 1D, the transport equation for a simple waterfront (with Pc = 0 and zero
gravity) is the well-known fractional flow equation (see Chapter 2) as follows:

 ∂S   ∂f 
φ  w  = − v w 
 ∂t   ∂x  (98)

where Sw is the water saturation and fw ( Sw ) is the fractional flow of water

 qw 
56 fw ( Sw ) =  
 qw + qo 
Numerical Methods in Reservoir Simulation
6
fw ( Sw )

 qw 
( fw ( Sw ) =   ) which is a function only of water saturation, Sw. We can
 qw + qo 
differentiate fw ( Sw ) by the chain rule to obtain:

fw ( Sw )  ∂f   ∂S 
 ∂S   ∂f 
φ w  = − v w  = − v w   w 
 ∂ft S =  ∂xqw  ∂Sw   ∂x 
w( w)   (99)
 qw + qo 
 ∂f 
− v w 
where the termfw ( Sw)∂Sw  is a non-linear term giving the water velocity, v w (Sw ) ,
which is also a function of water saturation: that is:

 ∂f 
v w (S w ) = − v  w  fw (Sw )
 ∂Sw  (100)
To simplify the problem even further for our purposes here, we take a straight line
 ∂fw 
 
fractional flow, fw (Sw ) , which means that  ∂Sw  is a constant and, hence, so is
the water velocity, Vw. Hence, the simple 1D transport equation for a convected
waterfront is:
 ∂fw 
 w  ∂S   ∂Sw 
∂S
  = −wv w  
 ∂t   ∂x  (101)

where, as noted above, vw is now constant. This equation describes the physical
situation illustrated schematically in Figure 14.

t1 t2 Governing equation
1 ∂Sw ∂Sw
= - Vw
∂t ∂x
Figure 14
Sw Vw = constant
The advance of a sharp
Water Vw Oil
saturation front governed
0 x2 - x1
by the transport equation x1 x2 Vw =
x t2 - t1
with a constant velocity vw.

Starting from the transport equation 101 above, we can easily apply finite differences
using our familiar notation (Chapter 5) to obtain:

 Sn +1 − Swi
n
  Swi − Swi −1 
φ wi  = − vw   + Error terms
 ∆t   ∆x  (102)

where we have used the backward difference (sometimes referred to as the upstream
 ∂Sw 
 
difference) to discretise the spatial term,  ∂x  . Note that we have not yet specified
the time level of the spatial terms. If we take these at the n time level (known), then it

Institute of Petroleum Engineering, Heriot-Watt University 57

would be an explicit scheme and if we take them at the n+1 time level (unknown), it
would be an implicit scheme. It does not matter for our current purposes here, since
we are principally interested in the “Error terms” which are indicated in equation
102. To determine what these terms look like, we need to go back to the original
finite differences based on Taylor expansion of the underlying function, Sw(x,t), in
this case.

We can expand Sw (x,t) either in space (x) or in time (t) as follows:

SPACE
 ∂S  δx  ∂ Sw 
2 2
Sw (x, t ) = Sw (x 0 ) + δx. w  +   + higher order terms
(t fixed)  ∂x  2  ∂x 2  (103)

 ∂S  δt 2  ∂ 2S 
TIME Sw ( x, t ) = Sw ( t 0 ) + δt. +  2  + higher order terms
w w
 ∂t  2  ∂t 
(x fixed) (104)

We can rearrange each of the above equations 103 and 104 to obtain the finite
 ∂S   ∂S 
difference approximations for the  ∂x 
  and  
 ∂t  terms with the leading error
terms as follows (now ignoring the higher order terms):

∂S  S(x, t ) − S(x 0 )  δx  ∂ 2S 
SPACE   ≈   −  2 (105)
 ∂x   δx  2  ∂x 

 ∂S   S(x, t ) − S( t 0 )  δt  ∂ S 
2
TIME   ≈  − 2  ∂t 2  (106)
 ∂t   δt 
To return to our usual notation, we make the identities:

S( x, t ) ⇔ Sin +1
S( x 0 ) ⇔ Si −1
S( t 0 ) ⇔ Sin
δx ⇔ ∆x
δt ⇔ ∆t (107)

Equation 105 now becomes:

 ∂S   Si − Si −1  ∆x  ∂ S 
2
  ≈ −
 ∂x   ∆x  2  ∂x 2 
 
(108)

and equation 106 becomes:

n +1
 ∂S   Si − Si  ∆t  ∂ S 
n 2
  ≈ −
 ∂t   ∆t  2  ∂t 2 
(109)

58
Numerical Methods in Reservoir Simulation
6
We now substitute the above finite difference approximations into the governing
equation 102 with their leading error terms as follows:

 Sn +1 − Sin  ∆t  ∂ 2S   Si − Si −1  v w .∆x  ∂ 2S 
φ i  −   ≈ − v w +  
 ∆t  2  ∂t 2   ∆x  2  ∂x 2  (110)

Collecting the error terms together on the RHS gives:

 Sn +1 − Sin   Si − Si −1  v w .∆x  ∂ 2S  ∆t  ∂ 2S 
φ i  ≈ − v w +  +  
 ∆t   ∆x  2  ∂x 2  2  ∂t 2 
Error term (111)

The error term in the finite difference scheme is now clear and is shown in equation
111. However, it still needs some simplifying since it is a strange mixed term with
 ∂ 2S   ∂ 2S   ∂ 2S 
 2  and  2 
both  ∂x   ∂t  terms in it. We now want to eliminate the  ∂t 2  term
and this is done by using the original governing equation.

 ∂ 2S   ∂Sw   ∂S 
To obtain  
2  , first note from the original equation 101 that 
 = − vw  w 
 ∂t  ∂t   ∂x 
 ∂Sw 
and we can then differentiate   with respect to t as follows:
 ∂t 
∂  ∂Sw   ∂ 2Sw  ∂  ∂S 
  =  2  = − vw  w 
∂t  ∂t   ∂t  ∂t  ∂x  (112)

Thus, we can rearrange the RHS of equation 112 as follows:

∂  ∂Sw  ∂  ∂S 
− vw   = − vw  w 
∂t  ∂x  ∂x  ∂t  (113)

 ∂S 
We now return again to the governing equation (equation 101) and use it for  w 
in equation 113 to obtain:  ∂t 

∂  ∂Sw  ∂  ∂Sw  2  ∂ Sw 
2
− vw   = − v −
 w v  = v w 
∂x  ∂t  ∂x  ∂x 
w
 ∂x 2  (114)

and hence:

 ∂ 2 Sw  2  ∂ Sw 
2

 2  = v w 
 ∂t   ∂x 2  (115)

We now substitute this expression into the error term in equation 111 above to obtain
the following:

Institute of Petroleum Engineering, Heriot-Watt University 59

v w .∆x  ∂ 2S  ∆t  ∂ 2S  v w .∆x  ∂ 2S  v 2w ∆t  ∂ 2S 
 +  =  +  
2  ∂x 2  2  ∂t 2  2  ∂x 2  2  ∂x 2 

 v .∆x v 2w ∆t   ∂ 2S 
= w +  
 2 2   ∂x 2  (116)

The finite difference equation with its error term in equation 111 now becomes:

 Sin +1 − Sin   Si − Si −1   v w .∆x v 2w ∆t   ∂ 2S 

  ≈ − v w  + +  
 ∆t   ∆x   2 2   ∂x 2  (117)

In this equation, we now see that the form of the error term is exactly like a “diffusive”
term i.e. it multiplies a (∂2Sw/∂x2) term. Hence we identify the level of numerical
dispersion or diffusion, Dnum, arising from our simple finite difference scheme as:

 v .∆x v 2w ∆t   ∆x v w ∆t   ∆x v w ∆t 
D num =  w +  = v w . +  D num = v w . + 
 2 2   2 2  (118)  2 2 

 ∆x v w ∆t 
If D num = v w . +  , we can take such a small time step that v w ∆t << ∆x
For this case:  2 2 

v w .∆x
∆x≈
D num
v w ∆t <<
2 (119)

We can now solve two problems as follows

(i) Firstly, we may apply an explicit finite difference method to obtain an accurate
solution of the following convection-dispersion equation.

 ∂Sw   ∂Sw   ∂ 2 Sw 
  = − v w  + D. 2 
 ∂t   ∂x   ∂x  (120)

That is, where the numerical dispersion is much less than the physical dispersion due
to the explicit D term. If this solution is converged (i.e. does not change on refining
Δx and Δt) then we can plot this at a suitable time, t1, as shown in Figure 15.

Governing equation:
2
1 t = t1 Sw Sw Sw
= - Vw +D
t x x2
Sw Figure 15.
Vw = water velocity
Dispersed frontal displace-
0 ment with a known and con-
x
verged level of dispersion.

60
Numerical Methods in Reservoir Simulation
6
n +1
 Swi − Swi
n
  Swi
n
− Swi
n
−1 
(ii) Now take only the explicit finite difference φapproximation
 =−
tov wthe
 convection 
 ∆t   ∆x 
equation only, that is:

 Sn +1 − Swi
n
  Swi
n
− Swi
n
−1   v .∆t  n
φ wi  = − vw   which gives ⇒ Swi
n +1
= Swi
n
−  w  (Swi − Swi
n
−1 )
 ∆t   ∆x   ∆x.φ 
Solve this with a relatively coarser grid, Δx, but with a fine time step thus predicting
 vwv.∆ x n
w .∆t 
2x.φ.( Choose 1)
n +1
S = S n
Dwinum to be −
wi   Swi − Swi
n
−this level of dispersion to deliberately match that in
∆
part (i) above i.e. choose Δx such that Dnum = D. Plot out a saturation profile like
that in Figure 15 at the same time t1 and compare these.

8 CLOSING REMARKS

In this Chapter, we have introduced the student to finite difference approximations

of the partial differential equations (PDEs) that describe both single- and two-phase
flow through porous media. These discretised equations have led to systems of either
linear or non-linear equations which are then solved numerically. Methods for solving
these equations have been discussed in principle and some idea has been given of
how these are applied in practical reservoir simulation. At the end of the unit, some
discussion was presented on grid-to-grid flows and how these inter-gridblock properties
are averaged. A more mathematical of numerical dispersion was also given.

Institute of Petroleum Engineering, Heriot-Watt University 61

APPENDIX A: Some useful Matrix Theorems

This lists (without proof) some useful theorems from matrix algebra which underpin
much of the practical application work which we have described in this section.
A.x = b
The central problem which we are interested in is:
A.x = b
A.x = b where A is an NxN square matrix

1. If A.x = b actually has a solution, then A is said to be non-singular or invertible

Aand the following five conditions apply and are equivalent (i.e. if one of them hold,
they all hold): I
A
(i) A −1 exists where A −1 . A = I ; I is the identity matrix which has 1s on the
diagonal and zeros elsewhere. I is like the number “1” in that I.A = A.I = A .

(ii) det( A ) ≠ 0 ; the determinant of the matrix is non-zero. This is a number found
by “multiplying out” the matrixI.in A a=specific
A.I = Away.

A.x = 0vector x such that, A.x = 0 . In other words, if A.x = 0

(iii) There is no non-zero
, then the vector x must be zero - have 0s for all its elements.
A.x = 0
(iv) The rows of A are linearly independent.AAnd .. A

(v) the columns of A are linearly independent. Where (iv) and (v) say that we
cannot get one
Aof
.x the
= 0rows or columns by making some combination of the existing
ones.

2. A matrix A which is square and symmetric

( A = A T , whereAA.x
T
=is0the transpose of A; i.e. a ij = a ji ) is said to be positive
definite if: x .A.x > 0
T

3. If the matrixAA=isGpositive
.G T definite then it can be decomposed in exactly one
way into a product: A = G.G T

such that matrix G is lower triangular and has positive entries on the main diagonal.
G decomposition.
This is known as Cholesky

62
Numerical Methods in Reservoir Simulation
6
SOLUTIONS TO EXERCISES

EXERCISE 1.

Apply finite differences to the solution of the equation:

 dy  = 2. y 2 + 4
 dt 

where, at t = 0, y(t = 0) = 1. Take time steps of Δt = 0.001 (arbitrary time units) and
step the solution forward to t = 0.25. Use the notation yn+1 for the (unknown) y at
n+1 time level and yn for the (known) y value at the current, n, time level.

Plot the numerically calculated y as a function of t between t = 0 and t = 0.25 and

plot it against the analytical value (do the integral to find this).

Answer: is given below where the working is shown in spreadsheet CHAP6Ex1.xls.

This gives the finite difference formula, a spreadsheet implementing it and the analytical
solution for comparison.

SOLUTION 1.

Discretise the equation using the suggested notation as follows:

 y n +1 − y n 
 = 2.( y ) + 4
n 2

 n +1∆t n 
y −y 
 = 2.( y ) + 4
n 2

 ∆t 
which is rearranged to the explicit formula:
y n +1 = y n + ∆t (2.( y n )2 + 4)
 ny+n1+1 − yn n 
y = y + ∆=t (22.(.(yy )) ++ 44)
nn 22

 n +du ∆t   y n +1 − y n 
 y 1 − y n = 1 tan −n1 2 u 
∫ u 2 + α 2  α  = 2.( y ) + 4
n 2
= 2.( y ) +4 
which can be set  1 easily α
∆t up quite on a spreadsheet.  ∆ t 
du −1n 2u 
∫y 2 += αy 2 + ∆αt(tan
n +1 n =
2.( y ) + 4)
α might need the standard form of the following
To find the uanalytical answer, you
integral: y =2 y + ∆t2(2=.( y ) +tan
1n +1 dy n1 2
)
−1  y  y n++1C= y n + ∆t (2.( y n )2 + 4)
2 ∫ ydu+ ( 2 1)  2 ∫
n
4 = dt = t
2 2
1 dy= tan −11 u  −1  y 
∫2 u∫2y+2 α+ 2( 2α)2 = 2 2α tan   2  = ∫ dt = t + Cdu 1 u
du 1 −1  u  ∫ = tan −1  
α∫ u=2 +2α 2 α = tan
α u +α
2 2
α α
Using this1standard form gives:
α = 2 2 dy 2 = 1 tan −1  y  = dt = t + C
2 1∫ y + (−1 2 )y  2 2  2 ∫ 1 dy 1  y 
1 tandy 1+C −1  y 
= =
22 ∫ 2y 2 + ( 2 )22  2 2
t tan
 
= ∫ dt = t
2 +∫ Cy 2
+ ( 2 ) 2
=
2 2
tan −1
 2 
=
1 −1  y 
2
α2 =2 tan 2  2  = t +C
where C isy(the [
α t=) = 2 2 tan 2 2 (t + C )
this becomes:
]
constant of integration and we identify α = 2 above. Therefore,

2 2
−1
[ = t +C ]
 y2 2 (t + C )
y(1t ) =tan2 tan
 2
1  y−1  1  1  y 
C 1= tan −1tan = t+=C0.2176049 tan −1 = t +C
2  2 2Heriot-Watt
2 22 1Engineering,
Institute of Petroleum  University
2 2  2 63
−1  1 
C =
2 2
tan[
y(t ) = 2 tan 2 2 (t + C )
 2
= 0 .]2176049
∫ u + α = α tan  α 
2 2
du 1  u
∫ u + α = α tan  α 
−1
2 2
du 1  u
∫1 u + αdy= α tan
−1
1   y 
21 ∫ y +dy
α tan −1
= ∫ dt = t + C
2 2
=
2
( 2) 2
2 12  y2 
−1 

21 ∫ y +dy  ∫
= tan = dt = t + C
2
( 2) 2
2 12  2
−1  y 
α2 ∫= y 2+ ( 2 ) = 2 2 tan
2 2
 2 ∫
= dt = t + C

α= 2
α 1= tan 2 −1  y  = t +C
2 12  y2 
 
tan −1 = t +C
2 12  2 
−1  yto:
which easily rearranges 
[
y2(t )2=tan2 tan = t +C
 2  2 (t + C )
2 ]
[
y(t ) = 2 tan 2 2 (t + C ) ]
tan[2 12 (t += C0).2176049
y(t=) = 1 2 tan
C −1 ]
 conditions since y(t = 0) = 1; it is easy to see that:
We now find C from the initial

2 12 12 
C= tan −1   = 0.2176049
2 12  2
1 
C= tan −1  = 0.2176049
2 2  2 

which is then used in the analytic solution for y(t) above.

This is implemented in the spreadsheet CHAP6Ex1.xls. Note that at t > 0.3, both
the analytical and numerical solutions go a bit strange (very large and they start to
disagree) since the tan function has a singularity y(t) → ∞.

EXERCISE 2.

Fill in the above table using the algorithm:

∆t
Pi n +1 = Pi n +
∆x 2
( Pi n+1 + Pi n−1 − 2 Pi n )
Hint: make up a spread sheet as above and set the first unknown block (shown grey
shaded in table above) with the above formula. Copy this and paste it into all of the
cells in the entire unknown area (surrounded by red border above).

SOLUTION 2.

If you get stuck, look at spreadsheet CHAP6Ex2.xls on the disk.

EXERCISE 3.

Experiment with the spreadsheet in CHAP6Ex2.xls to examine the effects of the

three quantities above - Δt, Δx (or NX) and the solution as t → ∞.

64
Numerical Methods in Reservoir Simulation
6
EXERCISE 4.

Solve the following simple example of an upper triangular matrix:

4 x1 - 14xx21 −- 21xx32 −+ 21xx43 +
+ 11xx4 =+ 5 1x5 = 5
5

2 x − 22xx2 −+ 21xx3
2 3 4
++ 12xx4 = +122 x5 = 12
5

3 x3 + 34xx34 ++ 41xx54 =+301x5 = 30

2 x − 21xx4 =− 3 1x5 = 3
4 5

3 x5 = 153 x5 = 15

SOLUTION 4.
Answer Answer
x1 = ..........;
x1 = ..........; x2 = ..........;
x2 = ..........; x3 = ..........;
x3 = ..........; x4 = ..........;
x4 = ..........; x5 = ..........
x5 = ..........

EXERCISE 5.

fill in the table below for 10 iterations using a calculator or a spreadsheet:

POINT ITERATIVE SOLUTION OF LINEAR EQUATIONS

Iteration
counter x at iter. ν
ν x1 x2 x3 x4
First guess = 0 2.0 2.0 2.0 2.0
1
2
3
4
5
6
7
8
9
10

SOLUTION 5. Iteration
counter x at iter. ν
ν x1 x2 x3 x4
Answer First
- Exercise
guess =5: 0see spreadsheet
2.0 CHAP6Ex5.xls
2.0 2.0 - Sheet2.0
1
1 4.309677 1.97619 3.6925 2.784615
POINT ITERATIVE SOLUTION
2 OF LINEAR
4.008926 1.411795EQUATIONS
3.498943 2.436331
3 3.993119 1.505683 3.497206 2.503112
4 4.000937 1.500783 3.500617 2.500584
5 3.999963 1.499755 3.499965 2.499832
6 3.999986 1.500026 3.499995 2.500017
7 4.000003 1.5 3.500002 2.500001
8 4 1.499999 3.5 2.5
Converged ⇒ 9 4 1.5 3.5 2.5
10 4 1.5 3.5 2.5

Note - converges fully at k = 9 iterations.

Institute of Petroleum Engineering, Heriot-Watt University 65

1
2
3
4
5
6
7
8
9
10

Iteration
counter x at iter. ν
ν x1 x2 x3 x4
First guess = 0 2.0 2.0 2.0 2.0
1 4.309677 1.97619 3.6925 2.784615
2 4.008926 1.411795 3.498943 2.436331
3 3.993119 1.505683 3.497206 2.503112
4 4.000937 1.500783 3.500617 2.500584
5 3.999963 1.499755 3.499965 2.499832
6 3.999986 1.500026 3.499995 2.500017
7 4.000003 1.5 3.500002 2.500001
8 4 1.499999 3.5 2.5
Converged ⇒ 9 4 1.5 3.5 2.5
10 4 1.5 3.5 2.5

Note - converges fully at k = 9 iterations.

EXERCISE 6.

Which is best method (i.e. that requiring the lowest amount of computational work),
BAND or LSOR, for the following problems?

(i) NX = 5, NY = 3, Niter = 50 (a small problem)

(ii) NX = 20, NY = 5, Niter = 50

(iii) NX = 100, NY = 20, Niter = 70

(iv) NX = 400, NY = 100, Niter = 150

NX NY Niter BAND LSOR, Comment

WB WLSOR
5 3 50
20 5 50
100 20 70
400 100 150

SOLUTION 6.

NX NY Niter BAND LSOR, Comment

WB WLSOR
5 3 50 135 750 For a very small problem like this, a
direct method is usually better,
although both would take very little
time even on a PC.
20 5 50 2500 5000 Again a direct method is somewhat
better for a small problem.
100 20 70 8x105 1.4x105 As the problem grows in size,
iterative methods start to overtake
direct method.
400 100 150 4x108 6x106 For very large problems, iterative
methods virtually always win over
direct methods by more than an
order of magnitude.

66
Numerical Methods in Reservoir Simulation
6
EXERCISE 7.

The linear equations which arise in reservoir simulation may be solved by a direct
solution method or an iterative solution method. Fill in the table below:

Direct solution method Iterative solution method

Give a very
brief description
of each method

Main advantages 1. 1.

2. 2.

Main 1. 1.
disadvantages

2. 2.

Institute of Petroleum Engineering, Heriot-Watt University 67

68
Numerical Methods in Reservoir Simulation
6

Institute of Petroleum Engineering, Heriot-Watt University 69