5/13/2024

RESERVOIR SIMULATION
FINITE DIFFERENCE SOLUTIONS
TO PDES

GOVERNING EQUATIONS OF TRANSFER PHENOMENA

𝑘 𝑑𝑝
Flow in Porous Media Darcy’s Law 𝑣=−
𝜇 𝑑𝑥

𝑑𝑣
Fluid Mechanics Newton’s Law 𝜏 = −𝜇
𝑑𝑦
Conservation Laws Governing PDE
𝑑𝑇
Heat Transfer Fourier’s Law 𝑞 = −𝑘
𝑑𝑥

𝑑𝜑
Mass Transfer Fick’s Law 𝐽 = −𝐷 How to Solve?
𝑑𝑥

Dr. Khoozan 2

1
 5/13/2024

SOLUTION OF PDES
Analytical Methods
Limited Applicability
Simple Cases
Simple Geometries
Simple Boundary Conditions
Mostly Linear Problems

Numerical Methods
Finite Difference Method (FDM)
Finite Volume Method (FVM)
Finite Element Method (FEM)

Dr. Khoozan 3

FINITE DIFFERENCE METHOD

Finite Difference Method

Structured Grid

Purely Numerical

Do not take into account physical limitations

Simple Geometries

Weak Mathematical Background

Dr. Khoozan 4

2
 5/13/2024

FINITE VOLUME METHOD

Finite Volume Method

Takes into Account Conservation Laws

Uses Integration Over Volumes

Uses Divergence Theorem to Transfer Volume Integrals to Surface (Control Volume

Edges) Integral

Can Handle Complex Geometries

Multiple Physics is a Challenging Task

Requires Flux Approximations that may Introduce Errors

Dr. Khoozan 5

FINITE ELEMENT METHOD

Finite Element Method

Based on Variational Calculus

Strong Mathematical Background

Can Handle Multiple Physics

Can Handle Complex Geometries/Boundary Conditions

Classically Developed for Solid Mechanics

Various Methods for Different Physics: Mixed Finite Element Method, Discontinuous
Galerkin Method, …

Dr. Khoozan 6

3
 5/13/2024

STRUCTURED VS UNSTRUCTURED MESH

Structured Block Structured Unstructured

Local Grid Refinement

Dr. Khoozan 7

STRUCTURED VS UNSTRUCTURED MESH

Domain and Boundary Finite Difference Approach Finite Element Approach

(Structured Grid) (Unstructured Mesh)

Dr. Khoozan 8

4
 5/13/2024

CLASSIFICATION OF PARTIAL DIFFERENTIAL EQUATIONS

The general second-order nonhomogeneous partial differential equation in two

independent variables is:

Dr. Khoozan 9

CLASSIFICATION OF PARTIAL DIFFERENTIAL EQUATIONS

The terminology elliptic, parabolic, and hyperbolic chosen to classify PDEs reflects the analogy
between the form of the discriminant,𝐵2 − 4𝐴𝐶, for PDEs and the form of the discriminant, 𝐵2
− 4𝐴𝐶, which classifies conic sections.

Conic sections are described by the general second-order algebraic equation:

Dr. Khoozan 10

5
 5/13/2024

CLASSIFICATION OF PARTIAL DIFFERENTIAL EQUATIONS

Ellipse Parabola Hyperbola

Dr. Khoozan 11

CHARACTERISTICS OF PDE

What is the significance of the this classification?

What impact, if any, does the classification of a PDE have on the allowable and/or required initial
and boundary conditions?

Does the classification of a PDE have any effect on the choice of numerical method employed to
solve the equation?

Dr. Khoozan 12

6
 5/13/2024

CHARACTERISTICS OF PDE

In two dimensional space, characteristics are paths (curved, in general) in the solution domain
along which information propagates.

Information propagates throughout the solution domain along the characteristics paths.

Discontinuities in the derivatives of the dependent variable (if they exist) also propagate along the
characteristics paths.

Dr. Khoozan 13

CHARACTERISTICS OF PDE

If a PDE possesses real characteristics, then information propagates along these characteristics.

If no real characteristics exist, then there are no preferred paths of information propagation.

Consequently, the presence or absence of characteristics has a significant impact on the solution of
a PDE (by both analytical and numerical methods).

Dr. Khoozan 14

7
 5/13/2024

CHARACTERISTICS OF PDE

A simple physical example can be used to illustrate the physical significance of characteristic paths.

Convection is the process in which a physical property is propagated (convected) through space by
the motion of the medium occupying the space.

Fluid flow is a common example of convection. The convection of a property 𝑓 of a fluid particle in
one dimension is governed by the convection equation:

Dr. Khoozan 15

CHARACTERISTICS OF PDE
The location 𝑥(𝑡) of the fluid particle is related to its velocity 𝑢(𝑡) by the relationship:

Characteristic equation: The differential equation of the

characteristic path

Hence, the path of the fluid particle, called its path-line, is given by:

Along the pathline, the convection equation can be written as:

The fluid property 𝑓 is convected along the path-line, which is the characteristic path associated
with the convection equation.

Dr. Khoozan 16

8
 5/13/2024

CHARACTERISTICS OF PDE

The physical significance of the path-line (i.e., the characteristic path) as the path of propagation of
the fluid property 𝑓 is quite apparent for fluid convection.

Dr. Khoozan 17

CHARACTERISTICS OF PDE

Several procedures exist for determining the characteristics.

Discontinuities in the derivatives of the solution, if they exist, must propagate along the
characteristics.

Are there any paths in the solution domain 𝐷(𝑥, 𝑦) passing through a general point 𝑃 along which
the second derivatives of 𝑓(𝑥, 𝑦) , that is, 𝑓𝑥𝑥 , 𝑓𝑥𝑦 , and 𝑓𝑦𝑦 , are multivalued or discontinuous?

Such paths, if they exist, are the paths of information propagation, that is, the characteristics.

Dr. Khoozan 18

9
 5/13/2024

CHARACTERISTICS OF PDE

Characteristics of the PDE can be obtained by solving the equation:

Dr. Khoozan 19

CHARACTERISTICS OF PDE

Elliptic PDEs have no real characteristic paths, parabolic PDEs have one real repeated characteristic
path, and hyperbolic PDEs have two real distinct characteristic paths.

The presence of characteristic paths in the solution domain leads to the concepts of domain of
dependence and range of influence.

The domain of dependence of point 𝑷 is defined as the region of the solution domain upon which
the solution at point 𝑃, 𝑓 𝑥𝑃 , 𝑦𝑃 , depends. In other words, 𝑓 𝑥𝑃 , 𝑦𝑃 depends on everything that
has happened in the domain of dependence.

The range of influence of point 𝑷 is defined as the region of the solution domain in which the
solution 𝑓(𝑥, 𝑦) is influenced by the solution at point𝑃. In other words, 𝑓 𝑥𝑃 , 𝑦𝑃 influences the
solution at all points in the range of influence.

Dr. Khoozan 20

10
 5/13/2024

CHARACTERISTICS OF PDE

Domain of dependence (horizontal hatching) and range of influence (vertical hatching) of PDEs:

Elliptic Parabolic Hyperbolic

Dr. Khoozan 21

CLASSIFICATION OF PHYSICAL PROBLEMS

Physical problems fall into one of the following three general classifications:

Equilibrium problems

Propagation problems

Eigenproblems

Each of these three types of physical problems has its own special features,
its own particular type of governing partial differential equation, and it’s
own special numerical solution method.

Dr. Khoozan 22

11
 5/13/2024

EQUILIBRIUM PROBLEMS

Equilibrium problems are steady-state problems in closed domains 𝐷(𝑥, 𝑦) in which the solution
𝑓(𝑥, 𝑦) is governed by an elliptic PDE subject to boundary conditions specified at each point on the
boundary 𝐵 of the domain.

Equilibrium problems are jury problems in which the entire solution is passed on by a jury
requiring satisfaction of all internal requirements (i.e., the PDE) and all the boundary conditions
simultaneously.

Elliptic PDEs have no real characteristics. Thus, the solution at every point in the solution domain is
influenced by the solution at all the other points, and the solution at each point influences the
solution at all the other points

Dr. Khoozan 23

EQUILIBRIUM PROBLEMS

A classical example of an equilibrium problem

governed by an elliptic PDE is steady heat
diffusion (i.e., conduction) in a solid.

or

Dr. Khoozan 24

12
 5/13/2024

PROPAGATION PROBLEMS

Propagation problems are initial-value problems in open domains(open with respect to one of the
independent variables) in which the solution 𝑓(𝑥, 𝑡) in the domain of interest 𝐷(𝑥, 𝑡) is marched
forward from the initial state, guided and modified by boundary conditions.

Propagation problems are governed by parabolic or hyperbolic PDEs

The majority of propagation problems are unsteady problems.

The solution of a propagation problem is subject to initial conditions specified at a particular value
of the time-like coordinate and boundary conditions specified at each point on the space-like
boundary.

Dr. Khoozan 25

PROPAGATION PROBLEMS

The domain of interest 𝐷(𝑥, 𝑡) is open in the

direction of the time-like coordinate.

A classical example of a propagation problem

governed by a parabolic PDE is unsteady heat
diffusion in a solid, e.g. in one space dimension:

Dr. Khoozan 26

13
 5/13/2024

CLASSIFICATION OF PHYSICAL PROBLEM

Dr. Khoozan 27

TAYLOR SERIES AND FINITE DIFFERENCES

Any smooth function can be represented by a Taylor series expansion about a single point
using an infinite sum of terms containing the function’s derivatives, where the function, 𝑓,
can be evaluated at any point 𝑥 by expanding around another point 𝑥𝑖 in the infinite series.

It is of course impossible to utilize the infinite number of terms in the series to obtain an
exact function evaluation, but approximate solutions can be found by truncating the series;
that is, using a finite number of terms. This approximation introduces truncation error in
the evaluation of 𝑓(𝑥).

Dr. Khoozan 28

14
 5/13/2024

TAYLOR SERIES AND FINITE DIFFERENCES

Dr. Khoozan 29

TAYLOR SERIES AND FINITE DIFFERENCES

Dr. Khoozan 30

15
 5/13/2024

TAYLOR SERIES AND FINITE DIFFERENCES

Dr. Khoozan 31

FIRST-ORDER FORWARD DIFFERENCE APPROXIMATION

The Taylor series can be rearranged to solve for certain derivatives. Consider, for example,
that the reservoir pressure, 𝑝, is a function of position, 𝑥. Using a Taylor series, the
pressure at position 𝑥𝑖+1 using a reference position of 𝑥𝑖 can be evaluated:

Where ∆𝑥 = 𝑥𝑖+1 − 𝑥𝑖 . Rearranging and using some algebra to solve for the first derivative
of pressure yields:

Dr. Khoozan 32

16
 5/13/2024

FIRST-ORDER FORWARD DIFFERENCE APPROXIMATION

where the remainder term includes the truncated terms. If ∆𝑥 is taken to be sufficiently
small, each term in the remainder is smaller than the preceding term. Neglecting all terms
in the remainder, the derivative is approximated as a first-order finite forward difference.

Since the first (and largest) term in the neglected remainder is proportional to ∆𝑥 , the
forward difference approximation is of first order accuracy, written as 𝑂(∆𝑥). This means
that if ∆𝑥 is decreased by a factor of two, the error in the derivative would decrease by a
factor two, in the limit that ∆𝑥 approaches 0.

Dr. Khoozan 33

FIRST-ORDER FORWARD DIFFERENCE APPROXIMATION

This figure illustrates an intuitive schematic of a forward

difference approximation to a first-order derivative. Recall
that by definition, a function’s derivative is the slope of
the line tangent to the curve at a specific point, 𝑥𝑖 .

The forward difference approximation is the slope of a

“secant” line that connects two points (𝑥𝑖 and 𝑥𝑖+1 ).

In the limit that ∆𝑥 approaches zero, the secant line

becomes the tangent and, mathematically, all remainder
terms become zero, thus resulting in no truncation error.

Dr. Khoozan 34

17
 5/13/2024

FIRST-ORDER BACKWARD DIFFERENCE APPROXIMATION

The Taylor series can be written “backward” in space, i.e., 𝑥𝑖−1 = 𝑥𝑖 − ∆𝑥, where the
terms in the series now alternate between positive and negative as a result of the
subtraction of ∆𝑥 from 𝑥𝑖 (negative terms that are raised to an even power are positive).

Neglecting all terms higher than the first-order derivative and then rearranging to solve for
the derivative gives the first-order backward finite difference approximation,

Dr. Khoozan 35

FIRST-ORDER BACKWARD DIFFERENCE APPROXIMATION

The backward difference, like the forward difference, approximation contains

𝑂(∆𝑥) accuracy meaning that its associated truncation error is proportional to ∆𝑥.

Therefore, the forward and backward differences have the same order of accuracy.

This does not mean that the forward and backward approximations will give the
same answer, nor will they have the same error.

It does mean that both will converge to the true derivative at the same rate.

Dr. Khoozan 36

18
 5/13/2024

FIRST-ORDER BACKWARD DIFFERENCE APPROXIMATION

This figure illustrates the backward difference approximation, which is the
slope of the secant line between 𝑥𝑖 and 𝑥𝑖−1 .

Dr. Khoozan 37

SECOND-ORDER, CENTERED DIFFERENCE APPROXIMATION

Approximations to the first derivative can also be determined using a second order,
centered finite difference, which can be derived by subtracting a backward Taylor series
from a forward Taylor series:

Dr. Khoozan 38

19
 5/13/2024

SECOND-ORDER, CENTERED DIFFERENCE APPROXIMATION

Rearranging for the first-order derivative gives the centered difference approximation:

where the leading term in the remainder, 𝑅𝑐 , is proportional to ∆𝑥 2 . Thus, the centered
difference approximation is 𝑂(∆𝑥 2 ) accurate, which means it converges faster to the true
solution than either a forward or backward difference approximation.

Recall that if ∆𝑥 were decreased by a factor of two, the error in the derivative
approximation would also decrease by a factor of two for either a forward or backward
difference; however, the error would decrease by a factor of four for a centered difference.

Dr. Khoozan 39

SECOND-ORDER, CENTERED DIFFERENCE APPROXIMATION

This figure illustrates a centered difference approximation, which is a secant

line between the points 𝑥𝑖−1 and 𝑥𝑖+1 . Note that a centered difference does
require evaluation of the function at two points other than 𝑥𝑖 (𝑥𝑖−1 and 𝑥𝑖+1 )
which may require additional computations.

Dr. Khoozan 40

20
 5/13/2024

APPROXIMATIONS TO THE SECOND DERIVATIVE

Like the first-order derivative, the second-order approximations can be a
forward, backward, or centered difference. In this course, the centered
difference for the second-order derivative will usually be used.

Dr. Khoozan 41

APPROXIMATIONS TO THE SECOND DERIVATIVE

The leading term in the remainder in the approximation is proportional to ∆𝑥 2 , so the
centered difference approximation for the second order derivative has error of 𝑂(∆𝑥 2 ).

This equation can alternatively be derived by recognizing that the second-order derivative
is simply the derivative of the first-order derivative. Substitution of the centered difference
for the first-order derivatives gives:

Dr. Khoozan 42

21
 5/13/2024

APPROXIMATIONS TO THE SECOND DERIVATIVE

where 𝑥𝑖+1 refers to the midpoint of 𝑥𝑖 and 𝑥𝑖+1 , and 𝑥𝑖−1 refers
2 2
to the midpoint of 𝑥𝑖 and 𝑥𝑖−1 .

The approximation is 𝑂(∆𝑥 2 ) because both the forward and

backward approximations are 𝑂(∆𝑥) and are divided by ∆𝑥.

This approach is useful for deriving reservoir simulation equations

that involve spatial heterogeneities in permeability or grid size.

Dr. Khoozan 43

EXAMPLE: FINITE DIFFERENCE APPROXIMATIONS

Dr. Khoozan 44

22
 5/13/2024

EXAMPLE: FINITE DIFFERENCE APPROXIMATIONS

Dr. Khoozan 45

EXAMPLE: FINITE DIFFERENCE APPROXIMATIONS

Dr. Khoozan 46

23
 5/13/2024

EXAMPLE: FINITE DIFFERENCE APPROXIMATIONS

Dr. Khoozan 47

EXAMPLE: FINITE DIFFERENCE APPROXIMATIONS

Dr. Khoozan 48

24
 5/13/2024

EXAMPLE: FINITE DIFFERENCE APPROXIMATIONS

Dr. Khoozan 49

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

It is possible to derive formulas with higher-order accuracy

by including more terms in the Taylor series.

The additional accuracy, however, generally comes at the

cost of more computational effort.

Here we present a generalized approach for derivation of

these finite differences

Dr. Khoozan 50

25
 5/13/2024

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

The derivative can be approximated as a weighted sum of points where 𝑘 is the derivative
order and 𝑗 is the stencil. For example, 𝑘 = 2 indicates a second-order derivative and the
centered difference would include 𝑗 = −1, 0, 1.

Each of the stencil pressures is then defined using a Taylor series:

Dr. Khoozan 51

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

In the Taylor series, the number of terms used will determine the accuracy. Substitution
and grouping the derivatives yields:

Finally, set the coefficients of 𝑝𝑘𝑖 = 1 (where 𝑘 is the order of the derivative) and all other
coefficients to zero in order to minimize the truncation error. This will lead to a system of
equations that can be solved to obtain 𝛽s. There should be the same number of linearly
independent equations as unknowns (if more equations than unknowns are used they will
not be linearly independent)

Dr. Khoozan 52

26
 5/13/2024

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

Dr. Khoozan 53

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

Dr. Khoozan 54

27
 5/13/2024

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

Dr. Khoozan 55

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

Dr. Khoozan 56

28
 5/13/2024

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

Dr. Khoozan 57

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

Dr. Khoozan 58

29
 5/13/2024

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

Dr. Khoozan 59

GENERALIZATION TO HIGHER-ORDER APPROXIMATIONS

Dr. Khoozan 60

30
 5/13/2024

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

The PDE describing 1D single-phase, slightly compressible flow in a

homogenous porous medium without sources or sinks was derived as:

The PDE is well posed if one initial condition and two boundary conditions
are included, provided at least one of those initial or boundary conditions
defines pressure.

Finite differencing is one numerical technique to find approximate solutions

to these PDEs.

Dr. Khoozan 61

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

The goal is to transform the PDE into algebraic equations using

approximations for the partial derivatives.

This can be accomplished by spatially discretizing the reservoir into

discrete grid blocks (also referred to as cells or elements) each with
length ∆𝑥𝑖 and discretizing the time into intervals of ∆𝑡.

The slide shows a 1D porous medium (e.g., rock core sample) with
length 𝐿 discretized into 𝑁 equally spaced grid blocks.

Dr. Khoozan 62

31
 5/13/2024

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

Dr. Khoozan 63

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

For the numerical solution, it is assumed the pressure is constant in each

block 𝑖 and over each time interval, ∆𝑡; then an algebraic equation is written
for each block.

In 1D flow, the grid blocks are, in general, connected to two neighbors: 𝑖 − 1

to the left of the block and 𝑖 + 1 to the right.

In this course, it is generally assumed that grid points are placed exactly at
the center of the grid block and are therefore referred to as cell-centered or
block-centered grids.

Reservoir properties and state variables (porosity, pressure, saturation, etc.)

are defined at the block centers.
Dr. Khoozan 64

32
 5/13/2024

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

Using the finite difference approximations, the spatial derivatives can be discretized using
a centered difference approximation which was shown to be of 𝑂(∆𝑥 2) in accuracy based
on the truncated terms in the Taylor series:

where 𝑃 is used to denote the approximate solution for pressure. Note the important
switch from lower-case 𝑝 (true, exact pressure) to upper-case 𝑃 (approximate, finite
difference pressure).

Dr. Khoozan 65

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

Likewise, the PDE is solved in discrete increments of time, each separated by an interval
∆𝑡.

A time interval with a superscript “𝑛” is denoted and the next time level as “𝑛 + 1”.

A finite difference approximation in time can also be employed. Note that the superscript
indicates the time level and the subscript the spatial location. Discretization of the time
derivative gives:

Dr. Khoozan 66

33
 5/13/2024

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

This is a forward difference approximation if the current time is envisioned to be 𝑛 and

the future time to be 𝑛 + 1.

It is a backward difference if the current time is 𝑛 + 1 and the previous time is 𝑛.

1
It is a centered difference if the current time is 𝑛 + 2.

The backward (explicit method), forward (implicit method), and centered (Crank-
Nicolson method) time difference approximations are the starting points for the
numerical solutions described in this course.
Dr. Khoozan 67

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

The approximations can be substituted directly into the 1D diffusivity

Equation to give:

where the dimensionless diffusivity constant (sometimes called the Fourier

number), 𝜂, is defined as

Dr. Khoozan 68

34
 5/13/2024

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

Each term in the equation has units of pressure. It is often more convenient to write the
equations in units of flow rate (e.g., ft3/day or RB/day); after all, a primary interest of the
reservoir engineer is to determine the rate of oil, gas, and water produced. Multiplying this
equation by 𝐴 and 𝑐𝑡 :

where 𝑐𝑡 is the total compressibility, and 𝑎 = Δ𝑦Δ𝑧 is the cross-sectional area

perpendicular to the flow. The accumulation, 𝐴, and transmissibility, 𝑇, are defined by:

Dr. Khoozan 69

DISCRETIZATION OF THE PARABOLIC DIFFUSIVITY (HEAT) EQUATION

Given an initial condition (pressures of all blocks at time 𝑛), this equation can be solved for
𝑃𝑖𝑛+1 in all blocks with 𝑖 = 1 to 𝑁.

However, it must be determined at which time level (𝑛 or 𝑛 + 1) the right-hand side

pressures should be evaluated.

Additionally, in blocks 𝑖 = 1 and 𝑖 = 𝑁, a pressure in a nonexistent ghost block, 𝑃𝑖−1 = 𝑃0

and 𝑃𝑖+1 = 𝑃𝑁 , respectively, must be evaluated using boundary conditions.

Dr. Khoozan 70

35
 5/13/2024

BOUNDARY AND INITIAL CONDITIONS

Boundary conditions must be applied to solve any PDE.

For the 1D parabolic diffusivity equation, two boundary conditions (one at 𝑥 = 0 and the
other at 𝑥 = 𝐿) in addition to one initial condition are needed to define the problem.

There are many different possibilities for boundary conditions, but the most common that
are applied to flow in reservoir simulation are Dirichlet (constant pressure) and Neumann
(constant rate).

Dr. Khoozan 71

DIRICHLET BOUNDARY CONDITION

A Dirichlet boundary condition is one of
constant pressure.

At 𝑥 = 0, the pressure is specified as a

constant, 𝑝 = 𝑃𝐵1 , and is an example of a
Dirichlet boundary condition.

For cell-centered grids, boundary pressure is

not defined for any grid point in the reservoir
because the grid point is at the center of the
block while the boundary conditions are
imposed at the edges (i.e., 𝑥 = 0 or 𝑥 = 𝐿).
Dr. Khoozan 72

36
 5/13/2024

DIRICHLET BOUNDARY CONDITION

A crude approximation is to specify the pressure
of the left ghost block as being equal to the
specified boundary condition (in this example, 𝑃0
= 𝑃𝐵1 ).

However, this results in first-order error, 𝑂(∆𝑥),

since the center of the ghost block is actually
located at 𝑥 = − ∆𝑥Τ2 instead of 𝑥 = 0.

A better approximation to the boundary

condition is to specify the pressure at the edge as
an arithmetic average of the pressure in the
boundary block (#1) and the pressure in a ghost
block (#0) outside the reservoir
Dr. Khoozan 73

DIRICHLET BOUNDARY CONDITION

Therefore, for a Dirichlet boundary condition at 𝑥 = 0, one can write:

and an equivalent expression can be written for a Dirichlet condition at 𝑥 = 𝐿. This

approximation is of error 𝑂(∆𝑥 2 ), which is more accurate than using the crude
approximation of making 𝑃0 = 𝑃𝐵1.

The discretized heat equation includes the ghost block pressure (𝑃0 ) if 𝑖 = 1; so, the above
boundary condition equation can be substituted into it if the boundary is subject to a
Dirichlet boundary condition.

Dr. Khoozan 74

37
 5/13/2024

NEUMANN BOUNDARY CONDITION

A Neumann boundary condition is one of constant rate (or flux), e.g., 𝑞 = 𝑞𝐵2 at 𝑥 = 𝐿.
Darcy’s law states that flux is proportional to pressure gradient:

where a centered finite difference approximation is used for the first derivative.

This may “look like” a forward or backward difference, but because the flux is specified at
the boundary (i.e., the interface between block “𝑁” and “𝑁 + 1”) it is a centered
difference approximation with accuracy 𝑂(∆𝑥 2).

Dr. Khoozan 75

NEUMANN BOUNDARY CONDITION

If the Neumann boundary condition is
simply a “no flux” condition, i.e., 𝑞𝐵2 = 0, a
mirroring technique (𝑃𝑁 = 𝑃𝑁+1) can be
used as shown in thee figure (𝑥 = 𝐿).

This occurs for a sealed or insulated

boundary. Physically, this may represent the
beginning of an impermeable zone, like
shale, that does not allow for a significant
amount of flow.

Dr. Khoozan 76

38
 5/13/2024

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

The porous medium is discretized into 𝑁

grid blocks.

Therefore, there are 𝑁 equations, each

associated with a corresponding block in
the discretized reservoir simulator.

To solve for pressures in each grid block 𝑖

(i.e., block 𝑖 = 1, 2, … , 𝑁) using the explicit
method, the discretized heat equation is
written for each grid block
Dr. Khoozan 77

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

Dr. Khoozan 78

39
 5/13/2024

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

The spatial derivatives on the right hand side of the equation are evaluated at the time step 𝑛
in the explicit method.

Since initial conditions (𝑛 = 0) for pressure as a function of position are known, the solution
)
of grid block pressures (𝑃10, 𝑃20 , ⋯ , 𝑃𝑁 ) is known.

These initial conditions can be substituted into the above equation and the righthand side can
be evaluated.

Therefore, new pressures can be computed explicitly at the new time level, 𝑛 + 1.

The equations are not “coupled”, i.e., each equation is independent and does not depend on
the others.

Dr. Khoozan 79

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

Boundary conditions are used to evaluate pressures outside the

boundary (ghost blocks), 𝑃0 and 𝑃𝑁+1 .

For brevity, it is possible to write the algebraic Equations, with

boundary conditions, in matrix form.

The calculation of the new pressure vector, 𝑃𝑛+1 , requires

matrix/vector multiplication and addition/ subtraction, but NO
“system of equations” must be solved.
Dr. Khoozan 80

40
 5/13/2024

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

𝐽𝑖 = 2𝑇 and 𝑄𝑖 = 2𝑇𝑃𝐵𝑖 for a

Dirichlet BC

𝐽𝑖 = 0 and 𝑄𝑖 = 𝑞𝐵𝑖 , for a

Neumann BC

Dr. Khoozan 81

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

In matrix notation:

where 𝐀 is the diagonal accumulation matrix, 𝐜𝐭 is the diagonal

compressibility matrix, 𝐓 is a tridiagonal transmissibility matrix, 𝐉 is a diagonal
productivity matrix, 𝐈 is the identity matrix, and Q is a source vector.

Dr. Khoozan 82

41
 5/13/2024

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

Dr. Khoozan 83

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

Dr. Khoozan 84

42
 5/13/2024

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

Dr. Khoozan 85

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

Dr. Khoozan 86

43
 5/13/2024

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

Dr. Khoozan 87

EXPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

Dr. Khoozan 88

44
 5/13/2024

IMPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

To solve for pressures in each grid block 𝑖 (i.e., block 𝑖 = 1,2, … , 𝑁) using the
implicit method, discretized equation is written for each grid block:

Dr. Khoozan 89

IMPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

New pressures can be computed implicitly at the new time level, 𝑛

+ 1.

The equations are “coupled” (i.e., each equation is dependent on

the others).

A system of linear equations must be solved.

Dr. Khoozan 90

45
 5/13/2024

IMPLICIT SOLUTION TO THE DIFFUSIVITY EQUATION

In expanded matrix form these equations can be written as,

Dr. Khoozan 91

MIXED METHODS AND CRANK-NICOLSON

The time derivative in the diffusivity equation was approximated

as a forward difference in the explicit method and a backward
difference in the implicit method.

As a result, both methods are only first-order accurate in time but

different answers will be obtained when using each of the
methods.

To improve the accuracy of the solution in time without using

smaller timesteps, a weighted average of the two methods can be
used.
Dr. Khoozan 92

46
 5/13/2024

MIXED METHODS AND CRANK-NICOLSON

Multiplying the explicit method by 𝜃 (a weighting factor between 0 and 1)
and the implicit method by 1 − 𝜃 gives:

Dr. Khoozan 93

MIXED METHODS AND CRANK-NICOLSON

This reduces to the implicit method if 𝜃 = 0 and the explicit method if 𝜃

= 1; it is a mixed method for 0 < 𝜃 < 1.

If 𝜃 = 12, the method is known as the Crank-Nicolson (C-N) scheme and

has accuracy 𝑂(Δ𝑡 2 ).

For all mixed methods, a system of equations must still be solved and has
slightly more computational work compared to the implicit method, but
larger timesteps can be taken in C-N and maintain accuracy.

Dr. Khoozan 94

47
 5/13/2024

MIXED METHODS AND CRANK-NICOLSON

The additional accuracy obtained by using C-N might be more easily understood if the
method is viewed as using a centered difference approximation in time, where the time
derivative is centered around the 𝑛 + 1Τ2 time level.

Since this is a centered difference, the time derivative is 𝑂(Δ𝑡 2 ) in accuracy. The spatial
derivative is also at the 𝑛 + 1Τ2 time level; intermediate time levels do not exist but can
be approximated as an average of the “𝑛” and “𝑛 + 1” time intervals.

Dr. Khoozan 95

MIXED METHODS AND CRANK-NICOLSON

After some algebra, one obtains the C-N Eq which is 𝑂(Δ𝑡 2 ) accuracy in time
and 𝑂(Δ𝑥 2 ) accurate in space.

Dr. Khoozan 96

48
 5/13/2024

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 97

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 98

49
 5/13/2024

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 99

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 100

50
 5/13/2024

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 101

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 102

51
 5/13/2024

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 103

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 104

52
 5/13/2024

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 105

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 106

53
 5/13/2024

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 107

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 108

54
 5/13/2024

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 109

EXAMPLE: IMPLICIT AND C-N METHODS

Dr. Khoozan 110

55
 5/13/2024

LINEAR SYSTEMS OF EQUATIONS

Unlike the explicit equations, the implicit (and mixed methods) equations are coupled; a
system of 𝑁 simultaneous linear equations are formed.

The system of linear equations can be solved using direct or indirect methods.

Direct methods attempt to solve the system exactly and are not subject to truncation error.

Examples include Gauss Elimination and LU decomposition, which are theoretically both
𝑂(𝑁 3 )methods, meaning that increasing number of grids by 10-fold results in a computation
time that is 1000-fold larger.
Dr. Khoozan 111

LINEAR SYSTEMS OF EQUATIONS

Direct methods are often not feasible for very large (e.g., tens of thousands or millions
of blocks) systems of equations.

The Thomas algorithm is a fast, 𝑂(𝑁), direct method, but only applicable to 1D,
tridiagonal problems.

Indirect methods begin with a guess value for the solution and iterate until
convergence to a predetermined tolerance and are thus subject to truncation error.

Iterative methods are usually classified as stationary or Krylov subspace methods.

Dr. Khoozan 112

56
 5/13/2024

LINEAR SYSTEMS OF EQUATIONS

Stationary methods such as Jacobi, weighted Jacobi, Gauss-Seidel, or successive

overrelaxation (SOR) are faster than most direct methods for moderately large
systems (thousands of grids).

However, they are only guaranteed to converge for certain matrices (e.g.,
symmetric, positive definite) and are usually slow compared to Krylov subspace
methods for the very large systems often employed in reservoir simulation.

Krylov subspace methods formulate the system of equations as a minimization

problem and find the direction of steepest descent to find the optimum.

Dr. Khoozan 113

LINEAR SYSTEMS OF EQUATIONS

Examples include the conjugate gradient and GMRES methods.

Iterative methods often employ a preconditioner on the matrix which allows

for convergence in less iterations.

Solution to the linear system is often the computational bottleneck of the

simulator and the development of new solvers and preconditioners are an
ongoing area of research.

Dr. Khoozan 114

57
 5/13/2024

STABILITY AND CONVERGENCE

The implicit method requires significantly more computational work per timestep
than the explicit method because a system of linear equations must be solved.

Implicit methods are not, in general, more accurate than explicit methods.

The order of accuracy is the same for both the implicit and explicit methods since
the spatial derivatives in both methods were approximated to second-order
accuracy, 𝑂(∆𝑥 2 ) and the temporal derivatives were first-order accurate, 𝑂(∆𝑡).

Dr. Khoozan 115

STABILITY AND CONVERGENCE

The advantage of the implicit method is that it is more stable and convergent.

A method is said to converge if the solution becomes more accurate by

using smaller step sizes (∆𝑥 or ∆𝑡).

A method is said to be stable if errors (e.g., roundoff error) do not

propagate in the solution.

Stability and convergence are two different numerical issues

Dr. Khoozan 116

58
 5/13/2024

STABILITY AND CONVERGENCE

One can show using von Neumann stability analysis that an explicit solution
to the 1D homogeneous, diffusivity equation is stable and convergent if and
only if the following condition holds.

The solution becomes unstable if the chosen timestep is too large or the grid
block size is too small. This is unfortunate because one may want to improve
accuracy by reducing ∆𝑥 or decrease computation time by using a larger ∆𝑡.

Dr. Khoozan 117

STABILITY AND CONVERGENCE

Dr. Khoozan 118

59
 5/13/2024

STABILITY AND CONVERGENCE

The implicit and Crank-Nicolson solutions, on the other

hand, is unconditionally stable and convergent.

Meaning that any size grid block or any timestep can

be used and the solution is guaranteed to be stable.

However, small timesteps may still be required to obtain

an accurate solution.

Dr. Khoozan 119

60