page 1 of 12

NORGES TEKNISK-NATURVITENSKAPELIGE UNIVERSITET

INSTITUTT FOR PETROLEUMSTEKNOLOGI OG
ANVENDT GEOFYSIKK

Kontakt under eksamen:

Jon Kleppe 73 59 49 25
918 97300 (mobil)

SOLUTION

FINAL EXAM IN COURSE TPG4160 RESERVOIR SIMULATION

Monday, May 18, 2009

Time: 0900-1300

Grades will be available on: June 9, 2009

Allowed material: -Type-approved calculator, with empty memory, according to

NTNU´s list of approved calculators
-No printed or hand-written materials
Final Exam page 2 of 12
TPG4160 Reservoir Simulation, May 18, 2009

Question 1 (10, 2, 4, 2, 3, 5, 2 points)

The simple, one-dimensional, linear, horizontal, one-phase diffusivity equation may be

written as:
∂ 2 P φµc ∂P
=( ) .
∂x 2 k ∂t
a) Show the steps involved in deriving the diffusivity equation
b) Sketch the one-dimensional, horizontal porous system that the equation applies to, in both
continuous and discrete form.
c) Using Taylor series expansions, derive the finite difference approximations needed for the
discretization of the equation (for constant grid block size).
d) What are the error terms associated with these approximations?
e) Write the difference equation on explicit form, and outline a procedure for pressure solution.
f) Write the difference equation on implicit form, and outline a procedure for pressure solution.
g) Why is the explicit form seldom used?

Solution
∂ ∂
a) (1) -continuity equation − (ρu) = (φρ)
∂x ∂t
k ∂P
(1) -Darcy´s equation u = −
µ ∂x
1 ∂V 1 ∂ρ dρ
(1) -fluid compressibility c f = −( )( )T = ( )( )T ⇒ ρc f = ( )
V ∂P ρ ∂P dP
1 ∂φ
(1) -rock compressibility c r = ( )( )T ⇒ dφ = φc r
φ ∂P dP
-assume constant permeability and viscosity
∂ ⎛ k ∂P ⎞ ∂
(1) -substitution: ⎜ρ ⎟ = (φρ)
∂x ⎝ µ ∂ x ⎠ ∂t
∂ dφ ∂P dρ ∂P ∂P ∂P
(1) -right side: (φρ ) = ρ +φ = ( ρφc r + φρc f ) = ρφc
∂t dP ∂t dP ∂t ∂t ∂t
∂ ⎛ k ∂P ⎞ k ∂ ⎛ ∂P ⎞
(1) -left side: ⎜ρ ⎟= ⎜ρ ⎟
∂x ⎝ µ ∂ x ⎠ µ ∂x ⎝ ∂x ⎠
(1) -then:

∂ ⎛ ∂P ⎞ ∂ ⎛ ∂P ⎞ ⎛ ∂P ⎞ dρ ∂P ⎛ ∂ 2 P ⎞ ⎛ ∂P ⎞ 2 dρ ⎛ ∂ 2 P ⎞ ⎛ ∂P ⎞ 2
⎜ρ ⎟ = ρ ⎜ ⎟ + ⎜ ⎟ = ρ⎜ 2 ⎟ + ⎜ ⎟ = ρ⎜ 2 ⎟ + ⎜ ⎟ ρc f
∂x ⎝ ∂ x ⎠ ∂x ⎝ ∂x ⎠ ⎝ ∂x ⎠ dP ∂x ⎝ ∂x ⎠ ⎝ ∂x ⎠ dP ⎝ ∂ x ⎠ ⎝ ∂x ⎠

⎛ ∂ 2 P ⎞ ⎛ ∂P ⎞ 2
(1) -assume that ⎜ 2 ⎟ >> ⎜ ⎟ c f
⎝ ∂x ⎠ ⎝ ∂x ⎠
∂ ⎛ ∂P ⎞ ⎛ ∂ 2P ⎞
-so that ⎜ ρ ⎟ ≈ ρ⎜ 2 ⎟
∂x ⎝ ∂ x ⎠ ⎝ ∂x ⎠
∂ 2 P φµc ∂P
(1) -thus, the final equation becomes: =
∂x 2 k ∂t
Final Exam page 3 of 12
TPG4160 Reservoir Simulation, May 18, 2009

b) (1) -continuous system

(1) -discrete system

1 i-1 i i+1 N

Δx

c) (2) Space derivative

At constant time, t+Δt (alternatively, t may be used), the pressure function may be
expanded in forward and backward directions:
Δx (Δx) 2 (Δx) 3
P(x + Δx,t + Δt) = P(x,t + Δt) + P ′(x,t + Δt) + P ′′(x,t + Δt) + P ′′′(x,t + Δt) + .....
1! 2! 3!
(−Δx) (−Δx) 2 (−Δx) 3
P(x − Δx,t + Δt) = P(x,t + Δt) + P ′(x,t + Δt) + P ′′(x,t + Δt) + P ′′′(x,t + Δt) + .....
1! 2! 3!
By adding these two expressions, and solving for the second derivative, we get the
following approximation:
P(x + Δx,t + Δt) − 2P(x,t + Δt) + P(x + Δx,t + Δt) (Δx) 2
P ′′(x,t + Δt) = + P ′′′′(x,t + Δt) + .....
(Δx) 2 12
or, by employing the grid index system, and using superscript to indicate time level:
∂ 2P P t +Δt − 2Pit +Δt + Pi−1t +Δt
( 2 ) ti +Δt = i+1 + O(Δx 2 ) .
∂x (Δx) 2
(2) Time derivative
At constant position, x, the pressure function may be expanded in backward direction
in regard to time:
Δt (Δt) 2 (Δt) 3
P(x,t) = P(x,t + Δt) + P ′(x,t + Δt) + P ′′(x,t + Δt) + P ′′′(x,t + Δt) + .....
1! 2! 3!
By solving for the first derivative, we get the following approximation:
P(x,t + Δt) − P(x,t) (Δt)
P ′(x,t + Δt) = + P ′′(x,t + Δt) + .....
Δt 2
or, employing the index system:
∂P t +Δt Pit +Δt − Pit
( )i = + O(Δt) .
∂t Δt
d) (2) Error terms are derived above:
∂ 2P
( 2 ) ti +Δt : O(Δx 2 )
∂x
∂P t +Δt
( ) i : O(Δt)
∂t
e) (3) Explicit
Using the approximations above at time level t we get the explicit difference equation
(shown only in general form; end blocks will be different):
t
Pi+1 − 2Pit + Pi−1
t
φµc Pit +Δt − Pit
≈ ( ) , i = 1,...,N
Δx 2 k Δt
Final Exam page 4 of 12
TPG4160 Reservoir Simulation, May 18, 2009

Since the equation contains only one unknown, it may be solved explicitly:
Δt k
Pit +Δt = Pit + ( 2 )( t
)(Pi+1 − 2Pit + Pi−1
t
), i = 1,...,N
Δx φµc
f) (5) Implicit
Using the approximations above at time level t+Δt we get the explicit difference
equation:
t +Δt
Pi+1 − 2Pit +Δt + Pi−1t +Δt
φµc Pit +Δt − Pit
= ( ) , i = 1,...,N
Δx 2 k Δt
which is a set of N equations with N unknowns, which may be solved simultaneously,
using a number of solution methods, for instance Gaussian elimination:
t +Δt
ai Pi−1 + bi Pit +Δt + c i Pi+1
t +Δt
= di , i = 1,...N
g) (2) The explicit formulation is seldom used because it becomes unstable for large
time steps. It has the following stability requirement:
1 φµc
Δt ≤ ( )Δx 2 ,
2 k

Question 2 (2+3+4 points)

Sketch the coefficient matrix for the following systems, indicating non-zero diagonals with
approximate lines. Label the diagonals.

a) One-dimensional (x), one phase flow, with the pressure equation:

ai Pi−1 + bi Pi + c i Pi+1 = di , i = 1,N

applicable to the following grid system:

1 i-1 i i+1 N

b) Two-dimensional (x,y), one phase flow, with the pressure equation:

ei, j Pi, j −1 + ai, j Pi−1, j + bi, j Pi, j + c i, j Pi+1, j + f i, j Pi, j +1 = di, j i = 1,N x , j = 1,N y

applicable to the following grid system:

i

1 2 3 4 5 6
4 4

j 7 8 9 10 11 12
4

13 14 15 16 17 18
4

19 20 21 22 23 24
4

25 26 27 28 29 30
4

31 32 33 34 35 36
4

37 38 39 40 41 42
4

43 44 45 46 47 48
4
Final Exam page 5 of 12
TPG4160 Reservoir Simulation, May 18, 2009

c) Three-dimensional (x,y,z), one phase flow, with the pressure equation:

gi, j,k Pi, j,k −1 + ei, j,k Pi, j −1,k + ai, j,k Pi−1, j,k + bi, j,k Pi, j,k

+ c i, j,k Pi+1, j,k + f i, j,k Pi, j +1,k + hi, j,k Pi, j,k +1 = di, j,k i = 1,N x , j = 1,N y ,k = 1,N z

applicable to the following grid system (grid blocks numbered in the sequence of x,y,z)

i,j,k-1

k i-1,j,k i,j,k i+1,j,k

i,j,k+1

Nz Ny
j
1 i Nx 1

Question 3 (5+5+5 points)

For two-phase flow of oil and gas in a horizontal, one dimensional, linear porous medium, the
flow equations may be written as:

∂ ⎛ kk ro ∂Po ⎞ ∂ ⎛ φSo ⎞
⎜ ⎟ − qo′ = ⎜ ⎟
∂x ⎝ µo Bo ∂x ⎠ ∂t ⎝ Bo ⎠
and
∂ ⎛ kk rg ∂Pg kk ∂P ⎞ ∂ ⎛ φS φS ⎞
⎜⎜ + Rso ro o ⎟⎟ − q′g − Rsoq′o = ⎜⎜ g + Rso o ⎟⎟ ,
∂x ⎝ µg Bg ∂x µo Bo ∂x ⎠ ∂t ⎝ Bg Bo ⎠
where
Pcog = Pg − Po

So + Sg = 1.

a) Write the two flow equations on discretized forms in terms of transmissibilities, storage
coefficients and pressure differences (no derivations).
b) List the assumptions for IMPES solution, and outline briefly how we solve for pressures
and saturations
c) Outline briefly how we can solve for pressures and saturations by Newtonian iteration (ie.
fully implicit solution).
------
a) (5)
Final Exam page 6 of 12
TPG4160 Reservoir Simulation, May 18, 2009

Txoi+1 2 ( Poi+1 − Poi ) + Txoi−1 2 ( Poi−1 − Poi ) − q′oi

= Cpooi ( Poi − Poit ) + Csgoi ( Sg i − Sg it ), i = 1,N

Txg i+1 2 [( Poi+1 − Poi ) + ( Pcog i+1 − Pcog i )] + Txg i−1 2 [( Poi−1 − Poi ) + ( Pcog i−1 − Pcog i )] − q′gi
+( RsoTxo) i+1 2 ( Poi+1 − Poi ) + ( RsoTxo) i−1 2 ( Poi−1 − Poi ) − ( Rsoq′o ) i
= Cpog i ( Poi − Poit ) + Csgg i ( Sg i − Sg it ), i = 1,N
b) (5) IMPES solution
Assumptions:
Txo t ,Txg t
Cpoo t ,Cpog t
.
Csgo t ,Csgg t
Pcog t ,Rso t

Having made these approximations, the discretized flow equations become:

2 ( Poi+1 − Poi ) + Txoi−1/ 2 ( Poi−1 − Poi ) − q′oi

t t
Txoi+1/
= Cpooit ( Poi − Poit ) + Csgoit ( Sg i − Sg it ), i = 1,N

t
Txg i+1/ [
2 ( Poi+1 − Poi ) + ( Pcog i+1 − Pcog i )
t
]
+Txg i−1/
t
[
2 ( Poi−1 − Poi ) + ( Pcog i−1 − Pcog i ) − qgi
′
t
]
+( RsoTxo) i+1/ 2 ( Poi+1 − Poi ) + ( RsoTxo) i−1/ 2 ( Poi−1 − Poi ) − Rso t q′o
t t
( ) i

=C t
pog i (P oi −P t
oi ) + C (S t
sgg i gi − S ),t
gi i = 1,N

IMPES pressure solution

The pressure equation for the saturated oil-gas becomes:

{T t
xoi+1/ 2 [
+ α i Txg i+1/
t
2 + ( RsoTxo ) i+1/ 2
t
]}(P oi+1 − Poi ) +

{T t
xoi−1/ 2 + α [T i
t
xg i−1/ 2 + ( RsoTxo ) i−1/ 2
t
]}(P oi−1 − Poi )

+α iTxg i+1/ 2 ( Pcog i+1 − Pcog i ) + α iTxg i−1/ 2 ( Pcog i−1 − Pcog i )
t t t t

(
−q′oi − α i q′g + Rso q′oi =
t
) i

(C t
pooi + α iC t
pog i )(P oi − Poit ), i = 1,N

where

α i = −Csgoit /Csgg it .

The pressure equation may now be rewritten as:

 Final Exam page 7 of 12
TPG4160 Reservoir Simulation, May 18, 2009

ai Poi−1 + bi Poi + c i Poi+1 = di , i = 1,N

and solved for pressures using a number of solution methods, ie. Gaussian elimination.

IMPES saturation solution

Having obtained the oil pressures above, we need to solve for gas saturations using either the
oil equation or the gas equation. Using the oil equation yields:
Sg i = Sg it +
1
t
Csgoi
[ t
Txoi+1/ t
′ − Cpooit ( Poi − Poit ) ,
2 ( Poi+1 − Poi ) + Txoi−1/ 2 ( Poi−1 − Poi ) − qoi i = 1,N ]
c) (5)

Solution by Newtonian iteration

Let us express the oil equation as Foi and the gas equation as Fg i . Each equation will depend on
pressures and saturations in blocks i-1, i and i+1., as indicated below.

Foi (Poi−1,Poi ,Poi+1,Sg i−1,Sg i ,Sg i+1 ) = 0

Fg i (Poi−1,Poi ,Poi+1,Sg i−1,Sg i ,Sg i+1 ) = 0 , i = 1,N

By first-order Taylor series expansions, we obtain the following expressions, where iteration
level is given by k:

∂Foi ∂Foi ∂Foi

Foik +1 = Foik + k +1
(Poi−1 − Poi−1
k
)+ (Poik +1 − Poik ) + k +1
(Poi+1 − Poi+1
k
)
∂Poi−1 ∂Poi ∂Poi+1
∂Foi k +1 ∂Foi k +1 ∂Foi
+ (Sg i−1 − Sg i−1
k
)+ (Sg i − Sg ik ) + (Sg k +1 − Sg i+1
k
)
∂Sg i−1 ∂Sg i ∂Sg i+1 i+1

∂Fg i ∂Fg i ∂Fg i

Fg ik +1 = Fg ik + k +1
(Poi−1 − Poi−1
k
)+ (Poik +1 − Poik ) + (Po k +1 − Poi+1
k
)
∂Poi−1 ∂Poi ∂Poi+1 i+1
∂Fg i k +1 ∂Fg i k +1 ∂Fg i k +1
+ (Sg i−1 − Sg i−1
k
)+ (Sg i − Sg ik ) + (Sg i+1 − Sg i+1
k
)
∂Sg i−1 ∂Sg i ∂Sg i+1
i = 1, N
Thus, for a one-dimensional system we have 2N equations and 2N unknowns, and we can
easily solve for estimates of oil pressures and gas and water saturations. By applying
Newtonian iteration until we converge on a solution within some tolerance, we may obtain a
solution to the equations. Our linear equations for iteration step k+1 would then take the
form:
k +1
apooi Poi−1 + bpooi Poik +1 + cpooi Poi+1
k +1 k +1
+ asgoi Sg i−1 + bsgoi Sg ik +1 + csgoi Sg i+1
k +1
= doi

k +1
apog i Poi−1 + bpog i Poik +1 + cpog i Poi+1
k +1 k +1
+ asgg i Sg i−1 + bsgg i Sg ik +1 + csgg i Sg i+1
k +1
= dg i
i = 1, N
 k +1 ˆ  k +1  k +1 
or, on a compact form: aˆ i X i−1 + bi X i + cˆ i X i+1 = di , i = 1,N
Final Exam page 8 of 12
TPG4160 Reservoir Simulation, May 18, 2009

where

apooi asgoi bpooi bsgoi  k +1 Poik +1  do

aˆ i = bˆi = X i = k +1 di = i
apog i asgg i bpog i bsgg i Sg i dg i

The equations are solved for pressures and saturations iteratively, updating coefficients after
each iteration.

Question 4 (4+4+6 points)

In the following 2-dimensional cross-section of a reservoir (one fluid only), a well is

producing at a constant rate Q (st. vol. oil/unit time) and perforates the grid blocks 4, 8, 11,
14, 17 and 21 in the x-z grid system shown:

1 5 9 13 17 21

2 6 10 14 18 22
2

j 3 7 11 15 19 23
3

4 8 12 16 20 24
4

1 2 3 4 5 6
i
The (unknown) bottom hole pressure Pbh is specified at a reference depth dref . Assume that
hydrostatic pressure equilibrium exists inside the well tubing.
a) Write the expression for oil rate from each perforated block (in terms of productivity
indices, mobility terms, pressure differences and hydrostatic pressure differences)

qij = PIij λij [ Pij − Pbhij + (dij − dref ) ρg]

ie.
q4 = PI4 λ4 (P4 − Pbh + Δd4 )
q8 = PI8 λ8 (P8 − Pbh + Δd8 )
...
...
b) Write the expression for the total oil flow rate for the well (group constants into parameters
A, B, C, D, E, F, G, H)

Qtot = ∑ q = A + BP4 + CP8 + DP11 + EP14 + FP17 + GP21 + HPbh

perf blocks ij

or
BP4 + CP8 + DP11 + EP14 + FP17 + GP21 + HPbh = d25

c) The standard pressure equation for this grid system, without the well terms, is:

ei, j Pi, j−1 + ai, j Pi−1, j + bi, j Pi, j + c i, j Pi+1, j + f i, j Pi, j +1 = di, j i = 1,...,N1, j = 1,...,N 2
Final Exam page 9 of 12
TPG4160 Reservoir Simulation, May 18, 2009

Sketch the coefficient matrix for this system, including the well. Indicate how the
coefficient matrix is altered by the well (approximately, with x´s and lines labelled with the
appropriate coefficient name).

b f c
e b f c
e b f c
e b c x
a b f c
a e b f c
a e b f c
a e b c x
a b f c
a e b f c
a e b c x
a b f c
a e b f c
a e b f c x
a e b f c
a e b c
a b f c x
a e b f c
a e b f c
a e b c
a b f x
a e b f
a e b f
a e b
B C D E F G H

Question 5 (2, 2, 6 points)

The discretized form of the left hand side of the oil equation may be written in terms of
transmissibilities and pressure differences, as
T xoi +1 / 2 ( Poi +1 − Poi ) + Txoi−1 / 2 ( Poi−1 − Poi )
Using the following transmissibility as example,
2ki −1 / 2 λ oi−1 / 2
T xoi −1 / 2 =
Δ xi (Δ xi + Δ xi −1)
a) What type of averaging method is normally applied to absolute permeability between grid
blocks? Why? Write the expression for average permeability between grid blocks (i-1) and
(i).
b) Write an expression for the selection of the conventional upstream mobility term for use in
the transmissibility term of the oil equation above for flow between the grid blocks (i-1) and
(i).
c) Make a sketch of a typical Buckley-Leverett saturation profile resulting from the
displacement of oil by water (ie. analytical solution). Then, show how the corresponding
profile, if calculated in a numerical simulation model, typically is influenced by the choice
of mobilities between the grid blocks (sketch curves for saturations computed with upstream
or average mobility terms, respectively).

Solution 3+3+6 points

a) Harmonic average is used, based on a derivation of average permeability of series
flow, assuming steady flow and using Darcy´s equation

⎧ λo if Poi−1 ≥ Poi
b) λoi−1/ 2 = ⎨ i−1
⎩ λoi if Poi−1 < Poi
Final Exam page 10 of 12
TPG4160 Reservoir Simulation, May 18, 2009

c)
Qw
Sw
1-Swir
exact
average
upstream
S wir

Question 6 (6, 1, 1, 1, 1, 3 points)

For a one-dimensional, horizontal, 3-phase oil, water, gas system, the general flow equations
are (including well terms):

∂ ⎛ kk ro ∂Po ⎞ ∂ ⎛ φSo ⎞
⎜ ⎟ − q′o = ⎜ ⎟,
∂x ⎝ µo Bo ∂x ⎠ ∂t ⎝ Bo ⎠
∂ ⎛ kk rg ∂Pg kk ∂P ⎞ ∂ ⎛ φS φS ⎞
⎜⎜ + Rso ro o ⎟⎟ − q′g − Rsoq′o = ⎜⎜ g + Rso o ⎟⎟
∂x ⎝ µg Bg ∂x µo Bo ∂x ⎠ ∂t ⎝ Bg Bo ⎠
∂ ⎛ kk rw ∂Pw ⎞ ∂ ⎛ φSw ⎞
⎜ ⎟ − q′w = ⎜ ⎟
∂x ⎝ µw Bw ∂x ⎠ ∂t ⎝ Bw ⎠

a) Explain briefly the physical meaning of each term in all three equations.
b) What are the criteria for saturated flow? What are the functional dependencies of
Rso and Bo ?
c) What are the primary unknowns when solving the saturated equations?
d) What are the criteria for undersaturated flow? What are the functional dependencies of
Rso and Bo ?
e) What are the primary unknowns when solving the undersaturated equations?
f) Rewrite the equations above for undersaturated flow conditions.

Solution
∂ ⎛ kk ro ∂Po ⎞ ∂ ⎛ φS o ⎞
a) ⎜ ⎟− q′o = ⎜ ⎟
∂x ⎝ µo Bo ∂x ⎠ ∂t ⎝ Bo ⎠
transport of oil well potential accumulation of oil
∂ ⎛ kk rg ∂Pg kk ∂P ⎞ ∂ ⎛ φSg φSo ⎞
⎜⎜ + Rso ro o ⎟⎟ − q′g − Rsoq′o = ⎜ + Rso ⎟
∂x ⎝ µg Bg ∂x µo Bo ∂x ⎠ ∂t ⎜⎝ Bg Bo ⎟⎠
transport of transport of gas well oil well pot. accumulation. accumulation
free gas sol. gas potential (solution gas) of free gas of solution gas
∂ ⎛ kk rw ∂Pw ⎞ ∂ ⎛ φSw ⎞
⎜ ⎟− q′w = ⎜ ⎟
∂x ⎝ µw Bw ∂x ⎠ ∂t ⎝ Bw ⎠
transport of water well potential accumulation of water
Final Exam page 11 of 12
TPG4160 Reservoir Simulation, May 18, 2009

b) Saturated flow: criteria Po = Pbp and Sg > 0.

dependencies Bo = f ( Po ) and Rso = f ( Po ) .
c) Po and Sg
d) Undersaturated flow: criteria Po > Pbp and Sg = 0.
dependencies Bo = f ( Po ,Pbp ) and Rso = f ( Pbp ) .
e) Po and Pbp

∂ ⎛ kk ro ∂Po ⎞ ∂ ⎛ φSo ⎞
f) ⎜ ⎟ − q′o = ⎜ ⎟
∂x ⎝ µo Bo ∂x ⎠ ∂t ⎝ Bo ⎠
∂⎛ kkro ∂Po ⎞ ∂ ⎛ φSo ⎞
⎜ Rso ⎟ − q′g − Rsoq′o = ⎜ Rso ⎟,
∂x ⎝ µo Bo ∂x ⎠ ∂t ⎝ Bo ⎠
∂ ⎛ kk rw ∂Pw ⎞ ∂ ⎛ φSw ⎞
⎜ ⎟ − q′w = ⎜ ⎟
∂x ⎝ µw Bw ∂x ⎠ ∂t ⎝ Bw ⎠
Final Exam page 12 of 12
TPG4160 Reservoir Simulation, May 18, 2009

Question 7 (2+5+2 points)

For a one-dimensional, vertical (z), 3 phase oil, water, gas system, outline how initial pressures
and saturations may be computed in a simulation model, assuming that equilibrium conditions
apply:
a) Sketch the reservoir, with a grid superimposed, including gas-oil-contact (GOC) and water-
oil-contact (WOC).
b) Sketch the oil-gas and oil-water capillary pressure curves, and show the how the initial
equilibrium pressures and saturations are determined in the contineous system.
c) Sketch the initial saturations as they are applied to the grid blocks.