ELEMENTS OF RESERVOIR SIMULATION )

BY

PROF. M.O. ONYEONWU

&

DR. WACHUKU PRINCE O.

 INTRODUCTION TO RESERVIOR SIMULATION

• The American Heritage dictionary defines “simulate” as “to give

appearance of”.

• In reservoir engineering, to simulate is to utilize a model to obtain

some insights into the behaviour of a reservoir.

• Two types of models used for reservoir simulation are physical and
mathematical models.

• In this course, interest will be only on mathematical models.

 INTRODUCTION TO RESERVIOR SIMULATION

• In reservoir engineering, many mathematical models are used.

• These include simple mathematical models such as, material

balance model, Buckley-Leverett model, pressure analysis
models, Dykstra-Parsons model, etc.

• For this course, the development of the complex mathematical

models are discussed.

• The development of such models require very large digital

computer.

• Hence, the models are known as computer or simulators.

 INTRODUCTION TO RESERVIOR SIMULATION
 The growth in computer systems has been a precursor to the development of
computer models known as simulators.

 The simulators are used for different types of studies.

 Information that may be obtained form simulation studies in the petroleum

industry include;
 development strategy for maximum hydrocarbon recovery,
 ultimate economic recovery
 Best completion practices
 Suitable enhanced recovery scheme for the reservoir
 and design parameters

 In addition, the computer models are used to complement the physical

models and also for studies on the effect of a wide range of parameters
which would involve a substantial effort to investigate using physical
 INTRODUCTION TO RESERVIOR SIMULATION
• Although there is increasing reliance on simulators in the petroleum
industry, it is necessary to be aware of the limitations.

LIMITATIONS OF SIMULATORS

• The simulators cannot be used to determine the physics of the

problem.

• Also, the quality of simulation results depend on the quality of the

input.

• In addition, the simulators is just a tool and cannot replace good

engineering judgment.

• Other ways in which computer models can be misused in the

1.1 Development of Reservoir Simulator
• Before the advent of computers, and hence computer models,
petroleum engineers have relied on other models such as the analog
models and physical models.

• The most common analog models are the electrical models in which
the similarities between electrical flow and fluid flow in petroleum
reservoir are utilized to develop electrical analogs of petroleum
reservoirs.

• These similarities are discussed by Chrichlow (1977). Due to the

uniqueness of petroleum reservoirs, and the size of the analog models,
their use has been discontinued.
 Development of Reservoir Simulator
• The physical models are laboratory models or scaled-down models,
which operate the same set of physical law as the reservoir.

• Such models are used in the petroleum industry to aid understanding

of the physical processes in the reservoir

and

 also produce results that can be applied in the field.

• One shortcoming of the physical models is the inability to use

such models for investigation of complex phenomenon.
 DEVELOPMENT OF RESERVOIR SIMULATOR
• With the advances in computer hardware and software technology,
petroleum engineers have sophisticated computer programs to
simulate some of the complex processes that take place in the
reservoir during the implementation of recovery schemes.
 DEVELOPMENT OF RESERVOIR SIMULATOR
• The processes/steps involved in the development of the computer
models are
 DEVELOPMENT OF RESERVOIR SIMULATOR
 The starting point in developing a reservoir simulator is to
identify the reservoir system and physical processes
occurring in the reservoir.

 Development of the mathematical model follows.

 This involves derivation of the mathematical equations that

describe the physical processes in the reservoir.

 Some assumptions may be made at this stage to make the

problem tractable.
 DEVELOPMENT OF RESERVOIR SIMULATOR
 The derived equations are usually non-linear partial
differential equations with appropriate initial and boundary
conditions.

 The equations constitute a mathematical model and are

usually too complex to solve using analytical methods.

 Hence, the equations are solved numerically.

 The procedure involved in numerical modeling consists of

discretizing the reservoir into blocks and performing mass
and energy balances on all the blocks simultaneously.
 DEVELOPMENT OF RESERVOIR SIMULATOR

• This gridding of reservoir allows for a more realistic

representation of rock and fluid properties which can vary
manner.

• In addition, dynamic effects of fluid movement within the

system are accounted for

• Discussion on different numerical techniques is given by

Aziz and Settari (1979)
 DEVELOPMENT OF RESERVOIR SIMULATOR
A typical grid is shown in below

Reservoir Boundary

Gridding the reservoir

 DEVELOPMENT OF RESERVOIR SIMULATOR
 A computer program written to solve the equations of the numerical
models constitutes a computer model of the reservoir.

 This is also generally known as a reservoir simulator.

 Hence the use of a computer model to solve reservoir engineering

problems is generally known as reservoir simulation.

 The outlined procedure will be discussed further with examples on

the development of gas reservoir simulator.
 Workflow for Conducting a Simulation Study
Seismic Data Acquisition

Data Processing

Seismic Interpretation
Petrophysical Evaluation
Geologic Model

Production Performance
Review BHP Analysis
Reservoir Engineering
Material Balance Pressure Data Review

Simulation Grid Model PVT Model

Simulation Model
Production Model Relative Permeability Model

Initialization

History Matching

Prediction Runs &

Economic Analysis

Development Plan
 Sources of Errors in Reservoir Simulators
 Simulators are interesting to use, but it is important to know that
results from the simulator may not be correct because of the following
reasons:

 The mathematical equations may not adequately represent the

processes taking place in the reservoir

 Blunders due to arithmetic errors or programming errors.

 Computer round-off errors resulting from the fact that a computer can
store a limited number of digits.

 Truncation errors as the numerical forms of the mathematical

equations are approximations of the mathematical equations.
 Classification of Reservoir Simulators
• Reservoir simulators can be classified on the basis of dimensionality,
geometry, phase, reservoir system

• Dimensionality: A model’s dimensionality is simply the number of

dimensions in which changes can occur.

• For most simulators, the number of dimensions can range from one to
three.

• Figures below show the three possible dimensions when Cartesian

coordinates are used.

• The conditions when the different dimensions will be used are

summarized below the figures
 Reservoir Simulator Selection
SIMULATOR SELECTOR

SYSTEM

Gas Black oil Chemical Compositional Thermal

Dimensionality

0-D 1-D 2-D 3-D

Geometry

X X-Y X-Z R-Z-θ X-Y-Z

PHASES

1-PHASE 2 - PHASE 3 - PHASE N-COMPONENT

18
Increasing Complexity and Expense
 Uses of 1-D Horizontal Model

Gas

Oil

Water

USED FOR
 Reservoir sections
 Specialized studies
 Line Drive behavior
 Miscible flooding
 Pilot Flood
 Sensitivity
 Effect of Heterogeneity in flow direction 19
 Uses of 1-D Radial Model

1-D Radial

 Single well study

20
 Uses of 1-D Dipping Model
gas

oil

Water

USED FOR
 Reservoir sections
 Specialized studies
 Line Drive behavior
 Miscible flooding
 Pilot Flood
 Sensitivity
 Effect of Heterogeneity in flow direction
 Laboratory studies
 Significant gravity override
 Up dip gas injection
 Flank Injection of water
 Uses of 2-D Horizontal Reservoir Dimension

 Most general dimension.

USED FOR
• Large multi-well structure
• Heterogeneous rock properties
• Slight vertical variation in rock and fluid properties
• Analysis of migration across lease lines
• Selection of optimum operational schemes/pattern in secondary recovery and
pressure maintenance
• 3-D cases (with pseudo-functions)
22
 Uses of 2-D Radial Reservoir Dimension

• Single well optimization studies

2 – D Radial
• Location of completion interval
• Maximum efficient rate
• Deliverability studies
• Well test analysis
23
Uses of a 2-D Vertical Reservoir Dimension

Gas

Oil

Water

USED FOR:
 Single or multiple well reservoirs
 Cross-section analysis of reservoirs
Investigates: - gravity segregation
- effect of anisotropy
24
 Uses of 3 Dimensional Layered Reservoir Model

Three dimensional (3D) reservoir layers

USED FOR
 Simulation of large reservoirs consisting of several producing horizons.
 Commingled or non-commingled production.
 Multiple completions.
 This model is in effect several 2-D models stacked together with special well bore hydraulics
25
routines.
Uses of 3-D Continuous Model
Used for:
 Simulation of large multi well
systems
 Thick reservoir pay sections
 Significant variation of rock
properties
 Significant variation of fluid
properties vertically
 Layered systems with common
aquifer or significant vertical cross
flow
 Areal and vertical detail needed

26
 Uses of 2-D Coning Model

USED FOR
 Single Well Optimization Study
 Location of Completion Interval
 Maximum Efficient Rates
 Deliverability Studies
 Well Test Analysis
 Reservoir Simulator Selection Implications

 With increasing dimensions, the cost of simulation

increases.

 This implies that lower dimensions should be used

whenever possible.

 Recently, methods have been proposed whereby 2-D models

can be used to simulate 3-D flow by use of pseudo-relative
permeability curves.
 Reservoir Geometry
 Simulators can be developed to model flow in system with different
coordinate system.

 The commonly used coordinate systems are the cartesian and cylindrical
coordinates.

 With the Cartesian coordinates, flow can be modeled in the x-y-z directions
while with cylindrical coordinate system, flow can occur in r-z-θ
direction.

 A typical cylindrical coordinates system where flow occurs in directions is

shown in the Figure on slide 27 for 2-D coning model

 Simulators that model flow in cylindrical systems are best suited for single
well studies.
 Phases
 Reservoir simulators can be classified on the basis of whether they can
model single, two or three phases.

 The phases modeled are oil, water, and gas. With the development of
thermal simulator, the fourth phase (solid phase) is now being considered.

 The cost of Simulation increase with the number of phases considered.

 Hence, it is important to use the minimum phases that adequately represent

the situation in the reservoirs.

 For example, with little or immobile gas phase, there is no need to include
the gas phase.

 Also, condensate reservoirs with immobile liquid phase can simply be

modeled as a single phase gas reservoir.
 Reservoir System
 Most commercial Simulators can be used to model fluid (oil, gas and water) in any
chosen dimension and geometry.

 Rather than classify these Simulators based on dimensionality, geometry and fluid
phase, these Simulators are classified on the basis of the reservoir system they can
be used to model.

 The different classifications are:

 gas reservoir Simulator,
 back oil reservoir Simulator,
 condensate reservoir Simulator,
 thermal Simulator and
 chemical flood Simulator.

 No single reservoir simulator can be used to effectively simulate all the reservoir
systems because of the variety of the processes going on in the reservoirs.
 Gas Reservoir Simulator

This is one of the simplest types of reservoir

Simulators.

The development of gas reservoir Simulators will be

discussed in details in this course.

In gas reservoirs with no mobile water phase, a gas

reservoir Simulators simply models single-phases
gas flow.
 Black Oil Simulator
 Black oil Simulators are widely used in the petroleum industry.

 The Simulators are used to model fluid flow in black oil reservoirs.

 The characteristics of black oil reservoirs are:

 Three distinct phases, oil, water and gas, are present.

 Water is the wetting phase while gas is the non-wetting phase. Oil has
intermediate wettability.

 Water and oil are immiscible and there is no mass transfer between the phases.

 Gas is soluble in oil and water.

 The reservoir is at uniform temperature and thermodynamic equilibrium is

attained throughout the reservoir.
 Black Oil Simulator Cont’d
Ideally, black oil simulators cannot be used to
model reservoir where there is mass transfer
between phases.

In situations where the mass transfer is minimal,

black oil simulators will be adequate.

 Black oil Simulators are less expensive to run than

Simulators where mass transfer is considered.
 Condensate Simulators

 These are generally known as compositional models.

 The major difference between black oil simulators and

compositional model is that the latter are used when the
mass transfer between hydrocarbon phases is significant.

 This is the case in very volatile hydrocarbon systems

OR
 In situations where gas-cycling operations are undertaken.
 Thermal Simulators

 Thermal Simulators are more complex than compositional models.

 In addition to modeling the mass transfer between phases, thermal

simulators also model the generation of heat and the effects of non-
uniform temperature in reservoirs undergoing thermal flooding
 such as:
 steam flooding,
 hot water flooding
and
 in-situ combustion.
 In thermal Simulators used for modeling in-situ combustion processes,
the chemical reactions are considered.

 For some of them, a fourth phase, the solid (coke) phase is also
 Chemical Flood Simulators:
 These Simulators model processes that occur in reservoirs
undergoing chemical flooding
 such as
 surfactant flooding,
 carbon dioxide flooding, etc.
 The major process modeled is the miscibility between the
injected fluid and in-situ reservoir fluids.

• Others: Reservoir Simulators can be classified on the basis of

formulation (e.g. IMPES, sequential, implicit, variational, finite-
difference, five-point, nine-point, etc) and solution technique
(e.g. direct solution, SOR, ADIP, etc).

• Costs (1982) present an overview of the various Simulators types,

 Uses of Reservoir Simulators
The use of reservoir simulators in the petroleum industry
has gained wide acceptance.

There is a cliché that states, “When all things fail,

simulate”.

Actually, it is not in all cases that one should simulate but

when judiciously applied, simulators are invaluable tools.

One of the advantages of simulators is that whereas the

field can be produced only once and at considerable
expense over a short period of time
 The uses of reservoir Simulators in reservoir engineering
(a) Estimation of Oil-In-Place and Reservoir Efficiency
 Simulators can be used for determination of the oil-in-place and
recovery efficiency for different recovery schemes.

 The oil-in-place and recovery efficiency are important economic

parameters for calculations of the commercial value of reserves.

 In a multi-zone reservoir, the oil-in-place in each zone is required by

the engineer for effective production scheduling and completion
operations.

 In cases where there may be used to determine the oil-in-place in each

lease.

 This information is essential in planning an equitable unitization

 The uses of reservoir Simulators Cont’d
(b) For Production Forecast
 A computer model can be used to make forecast on the oil
and gas production rate with time.

 This is required for planning delivery schedules and

development of cash flow for the project.

 In cases where filed data are available, simulated production

rate data are made to match the field data before a forecast of
future production is made.

Assignment
 Write a Term Paper on the uniqueness of the match between
 The uses of reservoir Simulators Cont’d
(c) Optimization of Petroleum Systems

 Simulation studies can provide information on bottom-hole pressures

of the wells.

 These pressures are used to plan the installation of downhole or

surface equipment.

 In secondary recovery projects such as water injection or gas

injection projects, simulation studies can provide information on the
required optimum water or gas injection rate and well locations.

 These parameters are required to design the size of the injection units,
water supply, water treatment, or gas processing plant size.
 The uses of reservoir Simulators Cont’d
d) To monitor Fluid Movement in Reservoirs
 In large reservoirs or reservoirs which are common to several
operators, it is possible that during the life of the project significant
quantities will move large distances from one lease to another.

 The movement is based on potential gradient regardless of what

subjective boundaries exist on the surface.

 With simulation studies on the reservoir, the migration of fluids can

be monitored and the location of the wells and the required production
rates chosen to control migration.

 In addition to monitoring migration pattern of the fluids in the

reservoir, the simulator enables the determination of the sweep
efficiencies of different recovery schemes.
 Once the flood fronts are located, the movable oil in the unswept
areas can be calculated and the location of new production wells
determined to maximize recovery.

 Hence, this provides information on optimal drilling sequence, the

conversion sequence from production to injection, and the optimal
water-oil ratio at which wells are shut in or converted to injectors.
 The uses of reservoir Simulators Cont’d
e) For Single-Well Studies

 The ability to design an optimum completion program is essential for

the exploitation of a reservoir.

 In some operations, it is not feasible to carry out a full-blown

reservoir simulation study, and a single-well study can be used to
obtain parameters that allow the engineer to determine the following:
1. Critical flow rates required to prevent coning of gas or water.

2. Effects of perforation interval and fracture penetration on well

productivity.

3. Maximum efficient rates to ensure optimum well response.

 The uses of reservoir Simulators Cont’d
(f) Well Testing
 Another area where simulators are now extensively used in the oil industry is
in well testing.
 In well testing, the petroleum engineer uses well pressure information to
obtain vital information about the reservoir.

 Such information includes the reservoir permeability, average pressure,

drainage volume, etc.

 Methods for analyzing well pressure data are discussed by Earlougher (1977).

 Simulators are used for well testing in complex systems.

 Examples of systems are fractured reservoirs, stratified reservoirs,
discontinuous reservoirs, reservoirs containing non-Newtonian fluids, and
reservoirs with peculiar conditions around the wellbore.
 The uses of reservoir Simulators Cont’d
(g) Other Applications

 In addition to the uses discussed, simulators are used for sensitivity

studies.

 Such studies reveal the important parameters that greatly influence

recovery.

 Hence, substantial effort is then made to determine the parameters more

accurately.

 Simulators are also used as educational tools with which the mechanics of
fluid flow in porous media can be investigate.

 For example, the complex interactions between gravity, viscous and

capillary forces can be investigated with reservoir simulators.
 DERIVATION OF FLUID FLOW EQUATIONS

 To develop a reservoir simulator, we need to derive equations

governing fluid flow in the reservoir.
 These equations are derived by combining the following basic
equations:

 Law of conservation equation (Darcy law, Ohm’s law, Fourier’s law,

Fick’s law)

 Transport rate equation

 Equation of state

 Equation of equilibrium.
These equations are discussed in details.
 Conservation Equation:
 The law of conservation simply states that mass, energy
and momentum can neither be created nor destroyed.

 Hence, for a chosen control volume, the mass, energy and

momentum transferred into the control volume must be
equal to the mass, energy and momentum that leaves the
control volume plus the mass, energy, and momentum that
remains in the control volume.

 The conservation equation can also be stated in terms of rate

of transfer.
 Let us use the conservation of mass equation to derive the equation
for single-phase flow in a linear system.

 The control volume is shown below

 The conservation equation is:

 Rate of mass in – Rate of mass out = Rate of mass storage

Rate of mass in = Rate of mass transport in + Rate
of mass injection into control volume
 q x0  qin Ax

Rate of mass out = Rate of mass transport out + Rate

of mass production out of control volume

 q x0  x  qout Ax
Mass of fluid stored in the control volume = Ax

Rate of mass storage in the control volume

=  
Ax   Ax  
t t
 Substituting into the conservation equation

   
q x0  qin Ax  q x0  x  qout Ax  Ax  
t

 To simplify the above Equation, Taylor’s Theorem. This theorem can
be used

 This Theorem gives a relatively simple method for approximating a

function f(x) that has continuous derivatives on a given interval by a
polynomial.

 Neglecting higher order terms ( x 2 , x 3 , etc) the Taylor’s series

expansion of f(x) about x is:
0

 f x 
f x   f x0   x  x0  x0
x
 By mere changing notation x  x0  x the above
equation becomes
 f x0 
f x0  x   f x0   x
x

 Using Taylor’s theorem,

 q 
q x0  x  q x0  x
x

 On substitution, we have that:

 q    
  qin  qout A  A
x t
 But
q  uA

 Therefore, the continuity equation becomes:

 u    
   qout  qint 
x t
 The equation holds for gases and liquids in a linear system.

 Following the same procedure, the continuity equation for

flow in the radial equation can be obtained.

 The equation is: 1  ur    

   qout  qint 
r r t
 Transport Rate Equation
This equation describes how mass, energy, etc,
moves when potential is applied.

For fluid flow in porous media, the transport rate

equation is the Darcy’s law which is expressed as:
k
u   p  z  ;

g
 
gc
The term z is positive in the vertical direction. The
Del operator, , is defined as:       
x y z
 This implies that: p p p
p   
x y z
k p
 For flow in the x-direction only, Darcy’s law is: u  
 x
 If oil, water and gas are present, the rate at which the oil
phase will move is:
kk ro p
uo  
 o x
 The term k is the relative permeability to oil and it will be
ro
discussed in detail later.

 Similar equation can be written for gas and water phases.

 Equation of State:
 The equation of state gives the relationship between fluid density,
pressure, temperature and other fluid properties.

 This equation makes it possible to cast the fluid flow equation in

different forms. pM
 
 For gases, the equation of state is zRT

 For fluid constant compressibility,

c( pthe
 pequation of state is
  be b)

 For incompressible fluid, c = 0.

 Hence the density is constant and mass balance will be equivalent

 Statement of Equation:
 If in a system, a component undergoes some phase changes, this must
be accounted for in the conservation equation.

 For example, in a reservoir that contains volatile oil, some oil will be
in the gas phase while some will remain in the liquid phase.

 The equilibrium statement is required to calculate the amounts in both

phases.
 The equilibrium equation is:
y i  K i x i

 The term yi is the amount of component i in the gaseous phase while

xi is the amount in the liquid phase.

 Ki is the equilibrium constant which depends on temperature,

pressure, and composition.
 Derivation of Gas Flow Equation in Linear System
 The continuity for flow in the linear system is:
 u    
   qout  qint 
x t
 In equation above, qout and qin are the mass rate of
production and injection per unit volume.

 Let us define q  qout  qin

 Neglecting gravity effect, which is usually small for gas

flow and substitute transport equation into the continuity
equation,
   g k g p    
 The resulting equation is:   q
x   g x  t
 But
B g  VRC /VSTC

 Hence

 STC
g  and q   gSTC qSTC
Bg

 Substituting for gas density and q , we have that:

  p     
  qSTC
 g    
x  x  t  Bg 
 Where kg
g 
 g Bg
 Porosity as a function of pressure is given as:
   0 1  c r ( p  p 0 ) 
  p    0 1  cr ( p  p0 )
 g      qSTC
x  x  t 
 Bg 


  p  p  d  1  0 c r


       qSTC
x  B 
g 0
x   t   dp  g  Bg 

  p  p
  g   c ( p )  qSTC
x  x  t

where
d  1  0 c r
C ( p )  0  
dp B  Bg
 g 
Generally,
p
   (p  z )  c ( p )  qSTC
t
 Other forms of the gas equation can be derived by substituting

pM
 into the continuity equation
zRT

  PM kp    PM 
    
x  zRT x  t  zRT 
 MULTIPHASE FLOW EQUATIONS

 Previously, single phase fluid in the porous medium was

considered.

 Here, cases with more than one fluid phase flowing will be
considered.

 That is, when there is water, oil and gas in the porous media.

 Equations for these cases will be derived.

 Before that, there is need to briefly define some relevant fluid and
rock properties.
 MULTIPHASE FLOW EQUATIONS cont’d
 Definition of Rock and Fluid Properties

 Definitions of rock and fluid properties such as porosity and

permeability with their types are discussed here.

 Porosity: This is a measure of the void spaces within a rock.

Porosity is expressed as a fraction (or percentage) of the bulk
volume of the rock.

 Absolute Porosity: This is the total porosity regardless of whether

the void spaces are connected or not.

 Effective Porosity: Porosity due to void spaces that are connected.

We are actually interested in effective porosity because fluid can
 MULTIPHASE FLOW EQUATIONS cont’d
 Porosity may be obtained from core samples or from logs.
 Porosity is dependent on pressure due to rock
compressibility .
 Mathematically, porosity at any a given pressure is:

  0 1  cr ( p  p0 )
0 = is porosity at the reference pressure p0
Viscosity of Fluid: Viscosity is a measure of fluid resistance
to flow.
 Mathematically, Viscosity is defined as:    yx
 yx = shearing stress d v x /d y

d vx /d y = rate of shearing strain

 MULTIPHASE FLOW EQUATIONS cont’d
Shearing stress has units of force per unit area and
1 2
dimensions of ML T

Rate of shearing strain1has units of velocity per length

T
and dimensions of

ML1T 1
Hence, viscosity has dimension of

In Oilfield units, viscosity is measured in centipoises

(cp)
and
In SI units, it is measured in kilopascal-second (Kpa-s)
 MULTIPHASE FLOW EQUATIONS cont’d
 For Newtonian fluid, viscosity is constant.
 For other fluids, viscosity is a function of shear stress and
rate of shear as shown in the Figure below

 Viscosity can be measured with Marsh funnel, Stormer viscosimeter

and multispeed rotational viscosimeter.

 The Marsh funnel gives a comparative value of viscosity based on

the time of flow of the fluid.
 MULTIPHASE FLOW EQUATIONS cont’d
 Absolute Permeability k: The absolute permeability of a rock is defined as the
permeability to 100% saturating fluid that does not react with mineral
components of the rock.

 It is a function of pore size distribution only and it is commonly measured in

Darcy units.

 Effective Permeability ko, kw, kg : Effective permeability is the permeability

to the fluid flowing under the saturation condition existing in the reservoir.

 In this case, there is always more than one fluid existing in this reservoir.

 The effective permeability depends on:

• Pore size distribution
• Wettability
• Saturation
• Saturation history
 MULTIPHASE FLOW EQUATIONS cont’d
 Relative Permeability kro, krw, krg: Relative permeability characteristics
are a measure of the ability of the porous system to conduct one fluid
when one or more fluids are present.
effective permeabili ty to phase i
 Mathematically, it is defined as: k ri 
base permeabili ty

 Either of the following may be used as the base permeability

• Absolute permeability, k

• The effective permeability of one of the hydrocarbon phases at

irreducible water saturation, Swi.

• The dry air absolute permeability kg measured at near atmospheric

pressure
 MULTIPHASE FLOW EQUATIONS cont’d

 The dry gas absolute permeability kg is greater than the

liquid absolute permeability kL due to gas slippage (a
phenomenon known as Klinkenberg or slippage effect).

 b
 The relationship between kg and kL is k g  k L 1  
 p

 b is a constant that depends on gas properties.

 The relationship between kg and kL is shown graphically in

Figure below
 MULTIPHASE FLOW EQUATIONS cont’d

 Depending on the base permeability used, the calculated relative

permeability values are different.
 For example, if a core has effective permeability to oil of 50mD and
absolute permeability as base permeability is
50
k ro   0.5
100
 If the same core has an air permeability of 115mD the relative
permeability to oil with air permeability as base is 50
k ro   0.435
115
 MULTIPHASE FLOW EQUATIONS cont’d

 If the irreducible water saturation of this core is 30%, and

the oil permeability at irreducible water saturation and zero
gas saturation is 0.7 of the absolute permeability,

 then the relative permeability based on effective

permeability to oil at irreducible water saturation is

50
k ro   0.713
100  0.7
 MULTIPHASE FLOW EQUATIONS cont’d
Drainage and Imbibition Relative Permeability

 Relative permeability is known to depend on saturation history

 That is, the direction of wetting phase saturation change processes

that result to decreasing wetting phase saturation are known as
drainage processes

 while processes that result in increasing wetting phase saturation are

know as imbibition processes.

 For example, flooding a water-wet reservoir with water is an

imbibition process while injecting gas into a water-wet reservoir is
drainage process.
 MULTIPHASE FLOW EQUATIONS cont’d
 Usually relative permeability data is obtained from core samples in the
laboratory by saturating or de-saturating the core with different fluids.

 Data obtained while the wetting phase is decreasing is known as draining

relative permeability while data obtained during increasing wetting phase
saturation is known as imbibition relative permeability

Some deductions that may be made about the effect of saturation history on
relative permeability include:

 Wetting phase relative permeability is not really affected by saturation

history

 At any given wetting phases saturation, the drainage relative permeability to

the non-wetting phase is greater than imbibition relative permeability.
• This is probably due to trapped saturation because during imbibition, the
 MULTIPHASE FLOW EQUATIONS cont’d
• A sketch of imbibition and drainage relative permeability are shown
in Figure below

• The imbibition curves are indicated by using arrows pointing in the

direction of increasing wetting phase saturation.

• The reverse is true for the drainage relative permeability curves

 MULTIPHASE FLOW EQUATIONS cont’d

Fluid Distribution in Porous Media

 We can understand the manner in which different fluids are distributed in the porous
medium by considering a model of porous medium in Figure below published by Fat
1958.

 In this model, an individual pore is represented by a line segment between

connecting points.
 MULTIPHASE FLOW EQUATIONS cont’d
 Each pore has its own size, shape and surface characteristics.

 The average flow path from one point to another is greater than the
straight-line distance between the points.

 This is due to tortusoity.

 If only one fluid is present in the porous medium, all the flow areas
are open to flow for the fluid.

 But if two different immiscible fluids are present, the flow area to
each is reduced and flow path length is increased.

 Hence the relative permeability to each fluid is affected.

 MULTIPHASE FLOW EQUATIONS cont’d
 Mathematically if one of the fluid occupies just the big pores, then

 Ai  where
k ro  f
 l  Ai = area open to flow for phase i
 i  li = path length for phase i

 Relative permeability to the other fluid will be greatly reduced.

 The wetting phase preferentially occupies the smaller pores.

 Based on these, some relative permeability characteristics may be

explained.

 Discussion on these follows.

 MULTIPHASE FLOW EQUATIONS cont’d
Water/Oil Relative Characteristics

 Consider a drainage process where oil is displacing water in a water-

wet core that was initially 100% saturated with water.

 The fluid distribution in this system is shown in the Figure below

 Water, the wetting phases, occupies the smaller pores while oil occupies
the bigger pores
 MULTIPHASE FLOW EQUATIONS cont’d
 This implies that immediately oil that occupies the bigger pores enters, the relative
permeability to water, drops drastically.

 Relative permeability to water becomes zero at irreducible water saturation, as no

water is mobile.

 Although no water is mobile at the irreducible water saturation, the whole flow path
is not open to oil flow because water blocks some paths.
 A sketch of the relative permeability curves is shown in the Figure below
 MULTIPHASE FLOW EQUATIONS cont’d
 The saturation of the wetting phase depends on the wetting phase
saturation.

 The inferences made for the oil/water system may also be made for
the gas/liquid system with the liquid as the wetting phases.

 Usually, relative permeability data are obtained experimentally.

 There are also theoretical methods which show that porosity, capillary
pressure, and relative permeability are related.

 An approximate relationship between them is represented by the

Leverett J-function given as follows
pc k
J S  
 cos  
 MULTIPHASE FLOW EQUATIONS cont’d
Three Phase Relative Permeability Characteristics
 In a water-wet porous media occupied by oil, gas and water, the
wetting phase (water in this case) occupied the smallest pores and the
relative permeability of the wetting phase is a function of the wetting
phase saturation only.
 The non-wetting phase (gas in this case) occupies only the biggest
pores and the relative permeability to the non-wetting phase is
dependent on the non-wetting phase saturation only.
 Oil which is of intermediate wettability to oil is dependent both on the
gas and water saturations.
 Figure below is a schematic of the fluid distribution in oil/gas/water
system
 MULTIPHASE FLOW EQUATIONS cont’d

 From experiments, it has been found that the relative permeability to

water obtained from an oil/water system is same as the relative
permeability to water in an oil/gas/water system.

 Similarly, the relative permeability of gas in gas/liquid system can be

used to represent the relative permeability to gas in oil/gas/water system.

 These follow from the fact that non-wetting phase relative permeability
is dependent on non-wetting phase saturation only.

 The problem in an oil/gas/water system is how to predict the relative

permeability to oil as a function of water and gas saturations.

 Some experiments with three-phase system have been performed to

determine the oil relative permeability
 MULTIPHASE FLOW EQUATIONS cont’d
 However, there are theoretical models that utilize two phase relative
permeability data (oil/water and gas/liquid) to correctly predict the
relative permeability to oil in an oil/gas/water system.

 One of these models is the Stone’s (1973) model that is derived from
channel-flow considerations.
 In this, the relative permeability to oil in an oil/gas/water system is

k ro  k row  k rw k rog  k rg  k rw  k rg 

 With the restriction that k ro  0
 Negative values of k ro implies immobile oil.
 This model will reduce exactly to two-phase data if the relative
permeability at the end points is equal to one.
 That is
k ro S wc   k rog S g  0  1
 MULTIPHASE FLOW EQUATIONS cont’d
 The equation above can also be satisfied if the gas-liquid relative
.
permeability data are measured in the presence of irreducible water
and
 effective permeability to oil at irreducible water saturation used as
based permeability.

 Alternatively if we denote k row S wc   k rog S L  1  k rocw

Where S L  1  S g  S o  S wc

 Stone’s model can be modified as

 
k ro  k rocw k row / k rocw  k rw k rog / k rocw  k rg  k rw  k rg 
 In the above equation, any of the base-permeability may be used to
define relative permeability.
 MULTIPHASE FLOW EQUATIONS cont’d
Capillary Pressure
 This is the pressure difference existing across the interface separating
two immiscible fluids, one of which wets the surface of the rock in
preference to the other.
 Mathematically;
pc  p nw  p w
 In oil/water system with water as the wetting phase;
pcow  po  p w
 For oil/gas: system with oil as the wetting phase;

pcog  pg  po
 MULTIPHASE FLOW EQUATIONS cont’d
 Thus, capillary pressure can have either a positive or a negative value
depending upon the wettability preference of the rock.

 Like relative permeability, capillary pressure is affected by saturation

history.

 The difference in multiphase rock properties that depend upon the

direction of saturation change is known as hysteresis effect.
 A typical capillary pressure curve for an oil/water system is shown in
the figure below for drainage and imbibition processes.
 MULTIPHASE FLOW EQUATIONS cont’d
Some deductions that may be made from the Figure above are:
1.The pressure in oil phase, , must be greater than the pressure in water
phase, p o , before oil enters water-wet rock initially saturated with water.

2. The entrance pressure is called threshold or displacement pressure, pT ,

and it depends on rock wettability, oil-water interfacial tension, and pore
sizes on exterior of the rock sample.

3. The slope of the drainage capillary pressure curve is a measure of the

range of pore size distribution. The more nearly horizontal or the flatter
the capillary pressure curve, the more uniform the pore sizes within the
rock.

4. At irreducible saturation, the wetting phase can no longer be

displaced by applying pressure.
 MULTIPHASE FLOW EQUATIONS cont’d
Formation Volume Factor
 Formation volume factors appear in fluid equations

 They are used to relate the volume of a given phase at

reservoir conditions to the volume at stock tank conditions.

 The volume of a given phase at reservoir conditions include

the volume of all materials (or other components) dissolved
in it.

 In our discussion, we shall concentrate on defining the

formation volume factors for black oil model.
 MULTIPHASE FLOW EQUATIONS cont’d
The characteristics of the black oil model are:
1. Three distinct phases present: oil, water and gas.

2. Water is the wetting phase while gas is the non-wetting phase.

3. Oil has intermediate wettability.

4. Water and oil are immiscible and there is no mass transfer between the
phases.

5. Gas is soluble in oil but usually not in water.

6.The reservoir is at uniform temperature and thermodynamic

equilibrium attained throughout the reservoir.
 MULTIPHASE FLOW EQUATIONS cont’d
 For the black oil model, the formation volume factors are
then defined as:
 MULTIPHASE FLOW EQUATIONS cont’d
 The below Figures are schematics of variation of formation volume
factors and fluid viscosities with pressure.
 MULTIPHASE FLOW EQUATIONS cont’d
 From Equations on slide 90, the fluid densities at reservoir
condition is related to fluid densities at stock tank
conditions by the following equations.

 For oil, o 
1
Bo

 oSTC  Rs  gSTC 

 For water, 1
w   wSTC
Bw

 For gas 1
g   gSTC
Bb
 Continuity Equation for Multiphase system

 When there are more than one fluid in a system, the continuity
equation must be written for each of the fluids.

 That is, each fluid in the system must be conserved irrespective of the
form or phase in which the fluid exist.

 Let us consider flow in a linear reservoir in which multiphase exist.

 The mass conservation equation for fluid in the shaded control

volume shown in the Figure below is
 Continuity Equation for Multiphase system

Rate of mass in of fluid i 

plus() 

Rate of mass injection of fluid i 

minus(-)   Rate of mass storage of fluid i
Rate of mass out of fluid i 

minus (-) 
Rate of mass of production of fluid i 

 Continuity Equation for Multiphase system
 Let us assume that fluid i is being injected and produced at volumetric
rate per unit volume of qinSC and qoutSC respectively.

 The volumetric rates are measured at standard conditions.

 If the fluid has a saturation of S i in the system, the conservation

equation becomes

 q
i i x   inSC qinSC V   i qi x  x   intSCq outSCV    iVS i 
t

 Expanding the above with Taylor’s series and evaluating gives:

 
 x  i qi   V  iS i    inSC qoutSC  qinSC V
x t
 Continuity Equation for Multiphase system
 But
qi  Aui

 And noting that V  Ax and cancelling the common

terms,


  i u i     iS i    inSC qoutSC  qinSC 
x t

 Therefore, the general continuity equation is written as:


  i ui    iSi   inSC qoutSC  qinSC 
t
 Continuity Equation for Multiphase system

  i ui    iSi   inSC qoutSC  qinSC 
t
 The above equation applies to any types of reservoir as long as the
terms  i ui and  iS i are defined appropriately.

 If the fluid in the system is oil, water, and gas, the continuity equation
will be

 FOR OIL   o u o    oS o    oSC q oiloutSC  q oilinSC 
t

 FOR WATER   w u w    wS w    wSC q wateroutSC  q waterinSC 
t

 FOR GAS   g u g  


t
 
 gS g   gSC qgasoutSC  qgasinSC 
 Continuity Equation for Multiphase system
 The general form of Darcy velocity of fluid i is
 kk ri
ui  pi   i z 
i
g
Where i  i
gc

 It implies that:
 Continuity Equation for Multiphase system
 The pressure gradient terms p in the above five equations refer to
the pressure gradient in each phase.

 The pressures in each phase are not the same because of capillary
pressure.

 The interpretation of the term  i ui for the black oil simulator is

explained in the next section.
 Multiphase Flow Equation; Black Oil Reservoir
 The characteristics of black oil have been given.

 The term is simply the mass flow rate of fluid i per unit area open to
flow.

 This is defined for black oil as follows;

Multiphase Flow Equation; Black Oil Reservoir
 Multiphase Flow Equation; Black Oil Reservoir
 Substituting the above three equations into the general
continuity equation gives the gas flow equation as:

  gSTC uo      gSTC S o  
  ug   gSTC Rso     S g   gSTC Rso  
 Bg Bo  t   Bg Bo  
  gSC qgasoutSC  qgasinSC 
 The above equation simplifies to:

 ug   
uo    Sg So  
   Rso      Rso    qgSC

 Bg  t 
Bo    Bg Bo  

 Where qgSC  qgasoutSC  qgasinSC
 Multiphase Flow Equation; Black Oil Reservoir
 Substituting the general equations of Darcy Velocity for
fluid i (SLIDE 98) into the above equation, we have that:

 kk rg Rso kk rg     S g S o  
 pg   gz  po   oz      Rso    qgSC
 Bg  g Bo  o  t   Bg Bo  

If we consider that gas dissolves in water, the

equation on slide 102 will be modified to account
for this.

This is left as an exercise.

 Multiphase Flow Equation; Black Oil Reservoir
 For oil

 Substituting the above three equations into the continuity

equation for oil, we have     S o

  o      q oSC
 Bo  t  Bo 
 Where q
oSC  q oiloutSC  q oilinSC
 Multiphase Flow Equation; Black Oil Reservoir
 Substituting for fluid velocity, we have
 kk ro    S o 
 po   o z      qoSC
 Bo  o  t  Bo 
 For water
 Following the same procedure used for oil, the flow
equation for water is   w    S w 
       q wSC
 Bw  t  Bw 
 where q wSC  q wateroutSC  q waterinSC

 Substituting the fluid velocity, we have

 kk rw    S w 
 p w   w z      q wSC
 Bw  w  t  Bw 
 Summary of Multiphase Flow Equation
• Neglecting the injection and production terms, the
fluid flow equation for black oil can be summarized
as
  S o 
  o po   o z    
t  Bo 
  S w 
  w p w   w z    
t  Bw 
• And
   S g S o  
  g pg   g z  Rso o po   oz     Rso  
t   Bg Bo  
 ,
Summary of Multiphase Flow Equation
,
,
,

 On the previous slide there are six unknowns po , p , pg , S , S w , and

w o
Sg
 But we have only three equations.

 Hence, additional three independent equations are required for a

unique solution to be obtained.
 Summary of Multiphase Flow Equation
 The additional three equations are:

 The first equation is the saturation constraint while

 The last two equations are capillary pressure equations.

 They relate the wetting and non-wetting phase pressures.

 Summary of Multiphase Flow Equation
 The equations are solved for each grid block for every time
step.

 This is quite expensive as it requires a lot of computer time.

 Hence, different forms of the three phase flow equations are

usually used to minimize solution cost.

 This is discussed further in the next section.

 ,
Useable Forms of Multiphase Flow Equations
,
, and

 The cost of solving several equations simultaneously to obtain

required reservoir properties made it imperative to seek methods of
reducing the number of equations to be solved.

 An example of the reduced form of the equations will be given with

the two-phase formulation given as:
 Water is the wetting phase in this
formulation.
 The terms q w and q o are the
source/sink terms measured in units of
volume at standard conditions.

 The unknowns in the Equations are po

pw , S , and S w
o

 Different forms of the equations with

reduced number of unknowns and
equations are given.
 Formulation in po and pc

 Suppose that a unique inverse function to exist. That is,

S w  Fw  p c 
 Then,

S o  1  S w  1  Fw  pc 
 Function exists if is monotonically increasing or decreasing with
saturation

 From the fact that:

p w  po  pc
 Formulation in po and Sw

 On substitution we have that:

   Fw 
  w  po  pc    w z      q w
t  Bw 
 and
  1  Fw 
  o po   o z       qo
t  Bo 
 Therefore the first equation is written as:

   Fw 
  w po  pc   w z      q w
t  Bw 
 The new form of the formulation has two equations (the last two
equations) and two unknowns p and
o Fw
 Formulation in po and Sw

 The finite difference form of the above equations can be decoupled

under suitable conditions.

 This is discussed further.

 Multiplying the equation for oil with Bo and the one for water with
Bw gives:

 1 Sw 
Bo   o po   o z   Bo   
  Bo q o
t  Bo 

       S w 
Bw   w  po  pc S w   w z   Bw    Bw q w
   t  Bw 
 Formulation in po and S w
 The above equation can be written
=
as:

       S w 
Bw   w  po   oz  pc S w   w z   oz   Bw    Bw qw
   t  Bw 

 By addition,

Bo  o po   oz  Bw   w po   oz  B w   w  pc S w  z 
  
   1  Sw      S w  
 Bo      qo   Bw     qw 
 t  Bo    t  Bw  
 Where
   w   o
 ,

Formulation in po and S w
 In a finite difference, if the saturations in the last equation are taken
explicitly, the equation becomes:
Bo   o po   o z  Bw   w po   o z 
           
 Bo 1  S w  
B   q o  B S
w w

B   q w
 t  o   t  w  

o , w , pc , S w and  / Bw 
 With saturation explicit, t are known and
po
may be taken as function of

 This is known as the “implicit pressure explicit saturation” or IMPES

method.
 This method is widely used in reservoir simulation.
 When the explicit treatment of saturations is not justifiable, as in the
case of coning simulation, the equations remain coupled
 CLASSIFICATION OF DIFFERENTIAL EQUATIONS
AND BOUNDARY CONDITIONS
 The procedures for obtaining solution to the equations depend on the
nature of equation and associated boundary conditions.

 A differential equation is called an ordinary differential equation if

there is only one independent variable in the equation.

 For example,

2
d y dy 4.1
2
 0
dx dx

 In Equation 4.1 the independent variable is x while the dependent

variable is y
 CLASSIFICATION OF DIFFERENTIAL EQUATIONS AND
BOUNDARY CONDITIONS Cont’d
 A differential equation is a partial differential equation if the
equation contains more than one independent variable.

 For example,  2u  2u 4.2

 0
x 2
y 2

 The independent variable in Equation 4.2 are x and y while the

dependent variable is u.

 The order of a differential equation is the highest order of the

derivative in the equation.

 The derivation n y or d n y is an n-order differential

x n d xn
 CLASSIFICATION OF DIFFERENTIAL EQUATIONS AND
BOUNDARY CONDITIONS Cont’d

 Equations 4.1 and 4.2 are second-order differential equations while a

third-order differential equation is
d3 y
3
y0 4.3
dx

 Linearity of a differential equation depends on the degree of the

dependent variable.

 For example, the following equations are linear differential equations:

 2 p p
  f x  4.4
x 2
t

dn y d y
n
 60 4.5
dx dx
 CLASSIFICATION OF DIFFERENTIAL EQUATIONS
AND BOUNDARY CONDITIONS Cont’d
 Examples of non-linear equations are

 2 p p 2
 4.6
x 2
t

d2 y d y
y 2 6 0
dx dx 4.7
 Homogeneity of a differential equation depends on the form of the
independent variable in the equation.

d2 y d y
 For example 2
  y0 4.8
dx dx
is homogenous while
2
4.9
d y dy
2
  y  f x  is non- homogenous
dx dx
 Discussion on Classification of Differential Equation
 Elliptic, Parabolic and hyperbolic Differential Equations

 The general linear homogenous equation of the second-

order is given as:
 2u  2u  2u u u
a 2  2h  b 2  2 f  2 g  cu  0
x xy y x y 4.10

 Where a, h, b, f, g, and c may be constants or functions of x and y.

 The form of Equation 4.10 resembles that of a general conic section

which is
4.11
ax  2hxy  by  2 fx  2 gy  c  0
2 2
 Examples of PDE’S of Different Types
Linear, homogenous
 Examples of PDE’S of Different Types
 Non Linear,
Examples of PDE’S of Different Types
Examples of PDE’S of Different Types
 Examples of PDE’S of Different Types
 With the general equation (Equation 4.10), a = 1, h = 0, and b = 0.

 Hence, Equation 4.12 is parabolic as ab  h 2  0

 Note that the independent variable t in Equation 4.12 is equivalent to

the independent variable y in equation 4.10.

 Similarly, the equation is ellipse because

4.13
 p
2
 p 2
 0
x 2
y 2

which is greater than zero

ab  h  1
2
 Examples of PDE’S of Different Types
 Equation with variable coefficient may be elliptic, parabolic or
hyperbolic depending on the region.

 A similar but complicated classification can be carried out for linear

equation in three or more independent variables,

 the terms elliptic, parabolic and hyperbolic should be replaced by

their three dimension analogues (ellipsoidal, etc).

 However, it is customary to continue to use the two -dimensional

terms.
 For example,
 2u  0 4.14
is still considered elliptic even when the third dimension is considered.
 Examples of PDE’S of Different Types
 Carnaham (1969) presented an easy method of classification of
equations, which applicable to equation with more than two
independent variables.

 In this the general linear second-order partial differential equation in

the coefficients Ai,
n
 2u n
u

i 1
Ai
xi
2
 
i 1
Bi
xi
2
 cu  D  0 4.15

 Evaluated at the point may be 1, -1 or zero.

 Note the absence of mixed derivatives  2 u / xi x j i, j  in this

form
 Examples of PDE’S of Different Types
 Based on Equation 4.15, the following classification may be
made.
Examples of PDE’S of Different Types
 Classification of Boundary Conditions
 In deriving the differential equation that governs the physical processes of
interest, a control volume is chosen within the medium.

 The conditions at the boundaries are not considered.

 The solution to the resulting differential equation will therefore be a

general solution.

 The general solution will contain arbitrary constants, which are evaluated
by imposing the boundary conditions to produce a unique solution.

 Hence, although the differential equation is the same for phenomenon

controlled by the same physical process, the solutions will be different
under different conditions at the boundary.

 For example the general solution to the differential equation

 Classification of Boundary Conditions

 For example the general solution to the differential equation

2 y
 6 4.19
x 2

is
y  3 x  Ax  B
2

4.20

where A and B are arbitrary constants.

 The values of A and B for a particular problem for which
Equation 4.19 holds are determined by the boundary
conditions for that problem.
Classification of Boundary Conditions
 Classification of Boundary Conditions

4. Mixed boundary condition: In this, a combination of the above boundary conditions

may be used.

• For example, the Dirichlet boundary condition may be specified at one end while the
Neumann boundary condition is specified at another part of the boundary.

• A differential equation with associated boundary conditions is said to be well posed if

the solution is unique and stable.

• For a stable solution, small changes in the boundary conditions produce small changes
in the solution.

• Unique and stable solutions of both elliptical and parabolic differential equations may
be obtained by using either dirichlet or Neumann boundary conditions.

• For hyperbolic equations, the Cauchy-Type boundary condition may be used.

• Understanding the physical processes modeled by the equations and the significance of
 Boundary and Initial Conditions used with Fluid Flow
Equations in Petroleum Engineering
 One of the fluid flow equations is the diffusivity equation for flow of
slight but constant compressibility fluid.

 This is given in cylindrical coordinates as:

 2 p 1 p ct p
 
r 2
r r k t
 In Equation 4.21, flow in both angular and Z directions are neglected.

 The boundaries in system where Equation 4.21 holds are shown in

Figure 4.1.

 The well location is the inner boundary while the outer limit of the
reservoir constitutes the outer boundary.
Boundary and Initial Conditions used with Fluid Flow Equations in
Petroleum Engineering
 The possible conditions that may be imposed on these boundaries are
discussed.
 The mathematical representations of these conditions are also given
Boundary and Initial Conditions used with Fluid Flow Equations in
Petroleum Engineering
Boundary and Initial Conditions used with Fluid Flow Equations in Petroleum
Engineering
 Boundary and Initial Conditions used with Fluid Flow Equations in
Petroleum Engineering
• Equation 4.31 is simply a mathematical expression used to represent
systems for which the outer boundary has not noticed the effect of
fluid withdrawal at the well.

• Hence, every reservoir goes through this stage which is known as

infinite acting stage.

• For equations with time as one of the independent variables, it is

required that the initial condition be given in addition to the boundary
conditions.
• The commonly used initial condition states that the pressure in the
reservoir is the initial pressure at time zero.

• This is p r ,0   pi
4.32
Boundary and Initial Conditions used with Fluid Flow Equations in
Petroleum Engineering

• With the appropriate boundary conditions the solution to the problem

may be obtained.

• Further discussion is given by Crichlow (1978).

• Although different types of boundary conditions have been given but

for “n” order differential equation, only “n” boundary conditions will
be required.

• For time-dependent problems, the initial condition will be required.

 FINITE DIFFERENCE APPROXIMATIONS (FDA)

 We have discussed methods of deriving the governing differential

equation using the conservation equations, transport rate equations
and equations of state.

 Some of the equations we derived can be solved analytically while

others can be solved numerically.

 The difference between an analytical solution and the numerical

solution is discussed later.

 One of the numerical methods used extensively in reservoir

simulation is the finite difference method.

 This will be the subject of further discussion.

 FINITE DIFFERENCE APPROXIMATIONS (FDA)
 Here, the development of the FDA will be illustrated as we try to
obtain solution to the equation of the following forms:

 2u
 q x  5.1
x 2

 2u u 5.2
  q x , t 
x 2
t

  u  u
  x, u    Cx,u   qx, t  5.3
x  x  t

1   u  u
 r  r , u    C r , u   q r , t 
5.4
r r  r  t
 FINITE DIFFERENCE APPROXIMATIONS (FDA)

 Equation 5.1 is the steady state equation

 Equation 5.2 is the transient linear flow equation.

 Equation 5.3 is the transient non-linear flow equation with variable

coefficients and

 Equation 5.4 is the transient radial flow equation with variable

coefficient.
 Analytical and Numerical Solutions

 The analytical solution of a differential equation is one that expresses

the dependent variable as a continuous function of the independent
variables.

 But with the numerical solution the value of the dependent variable is
known only at chosen discrete points within the medium.

Figure 5.1 shows the forms of an analytical solution and numerical

solution.
Analytical and Numerical Solutions
 Analytical and Numerical Solutions
 For a one dimensional fluid flow problem
2 p c t p

x 2 Ak t
 The initial and boundary conditions are

p x,0   pi

p 0, t   p0
 The numerical solution is an approximation to the correct solution of
the problem, which is usually given by the analytical solution.

 Hence, at a given time and location, the numerical solution may not
be exactly the same as the analytical solution.
 Analytical and Numerical Solutions
 The objective of every problem solving technique is to simplify the
problem so that a solution can be obtained.

 With numerical methods, differential equations are reduced to

algebraic equation before the solution to the differential equation can
be obtained.

 With numerical method, the solution is obtained at chosen points

within the domain of interest.

 This is known as discretization process.

Analytical and Numerical Solutions
 Analytical and Numerical Solutions
 Discretization in space is done by imposing a grid system on the domain of interest
and then obtaining solutions at discrete points within the domain of interest.

 For example, to obtain the solution to the fluid flow equation in the reservoir
shown in Figure 5.2a, a grid system may be imposed on the reservoir as shown in
Figure 5.2b.

 The imposed grid system may not necessarily be regular

 Discretization in time is done by selecting time steps at which the solution is

sought.

 This is done for problems with time as one of the independent variables.

 The chosen time steps may not be of the same size

 The smaller the time steps and grid sizes, the more accurate the numerically
obtained solutions
 FINITE DIFFERENCES

 With finite-difference method, the differential equation is replaced with

its finite-difference equivalent, which is always in algebraic form.

 This is so because the derivatives are expressed as functions of

independent parameter at discrete points in the grid system.

 The finite-difference equations can be derived using

• Taylor’s series method

• The integral method

• Variational method
 In this course, we shall use the Taylor’s series method to derive the
finite-difference approximations of derivatives
 Discretization in Space

 In this section, discretization in space is considered.

 The different approximations for discretization in space are given

 Forward Difference Approximation

Discretization in Space
 Discretization in Space
 The implications of Equation 5.8 is that using finite-difference
method, the derivative at point i can be expressed as a function of the
independent parameters at points i and i+1.

 Because the point i+1 which is ahead of point i is used,

 the formulation given by Equation 5.8 is called a forward difference

approximation.

 The truncation error is usually not known.

 Hence the approximation is simply written as

du u i 1  u i 5.9
i 
dx x
 Discretization in Space

 The forward difference approximation is shown graphically in Figure 5.4.

 In this figure, the actual derivative du/dx, at point i is the slope of the tangent to the
curve FEBDG, at point i.

 Using forward difference approximation, the derivative at point is taken as the slope
of the BDH
 Backward Difference Approximation
 Taylor’s series expansion

u x  x   u x   x u ' x  
x 2
u '' x  
x  u ''' x    5.10
2! 3!

 From Equation 5.10

u x   u x  x 
u x  
'
 b 5.11
x

Where  
x ''
u x  
x 2
u ''' x    5.12
b
2! 3!

 Ox  5.13
 Backward Difference Approximation

 Using the notation shown in figure 5.3, Equation 5.11 is

du u i  u i 1
i   b
dx x
5.14
 In Equation 5.14, the derivative at point i is approximated using
points i and i-1.

 Hence, the approximation is known as backward difference formula

 The order of the error in this case is unity as shown by Equation 5.13.

 A graphical interpretation of the backward difference approximation

is also shown in figure 5.4,

 the slope of the line IEB is used to approximate the derivative given
 Central Difference Approximation
 Subtracting Equation 5.10 from Equation 5.5 gives

u x  x   u x  x   2x u ' x   2
x 3
u ''' x   2
x 5
u v x    5.15
3! 5!

 Solving u’(x) for from Equation 5.15 gives

u x  x   u x  x 
u ' x    c 5.16
2x
where
c  
x  ''
2
u x  
x  v
4
u x   
5.17
3! 5!

 0x 
2 5.18
 Central Difference Approximation

 Using the notations in figure 5.3, gives

du u i 1  u i 1
i   c
dx 2x 5.19

 Equation 5.19 is the central difference approximation of the

derivative at point i.

 It has an error of order two as shown in Equation 5.18.

 As the central difference approximation has an error of order two, it is

more accurate than the forward or backward difference
approximations with error of order one.

 This is so because the higher the order of the error the smaller the
 Central Difference Approximation
 Graphical interpretation of the central difference approximation is
also shown in Figure 5.4.

 The slope of the line JEDK is used to approximate the derivative at

point i given by the slope of the line ABC.

5.20
 FINITE-DIFFERENCE APPROXIMATION OF SECOND DERIVATIVE

 Adding Equations 5.5 and 5.10 gives

2x  '' 2x  iv

2 4
u x  x   u x  x   2u x   u x   u x    5.20
2! 4!
 Solving Equation 5.20 for u11(x) gives

u x  x   2u x   u x  x 
u x  
''
 5.21
x 2

where  
x  iv
4
u x  
x  vi
6
u x  5.22
12 360

 O x 
2
5.23
 FINITE-DIFFERENCE APPROXIMATION OF SECOND DERIVATIVE

 Using the notation in Figure 5.3, Equation 5.21 is:

d2 u u i 1  2u i  u i 1
 
dx 2 i
x 2
5.24

 Equation 5.24 is the finite difference approximation of the second derivative at

point i.

 The formulation has an error of order two as shown in Equation 5.23.

 FINITE-DIFFERENCE APPROXIMATION OF SECOND DERIVATIVE

 Example: Write the finite-difference approximation of the

differential equation

d 2u du
 5  0 E5.1
dx 2 dx

 The differential equation holds in the interval

 Using the more accurate central difference approximation to represent

the first derivative, the finite-difference equation is:
ui i  2ui  ui 1
3
ui 1  U i 1 
5  0
x 2 2x E5.2
 The above equation is:

 3   3 
 1  x  i 1
u  2u i   1  x u i 1  5x 2
E5.3
 2   2 
 FINITE-DIFFERENCE APPROXIMATION OF SECOND DERIVATIVE

 If we define
3
A  1 x
2

3
B  1 x
2

and
C  x 
2

 Equation E5.3 becomes

E5.4
Au i 1  2u i  Bu i 1  5C
 FINITE-DIFFERENCE APPROXIMATION OF SECOND DERIVATIVE

 The given differential equation holds in an interval which we might

divide into three parts to obtain a finite- difference solution.

 Hence i will take values 1, 2, and 3.

 The resulting equations for the values of i are

Au0  2u1  Bu2  5C

Au1  2u 2  Bu3  5C

Au 2  2u3  Bu 4  5C
 FINITE-DIFFERENCE APPROXIMATION OF SECOND DERIVATIVE

,
 Hence, using the finite difference method, the resulting system of
algebraic equations will be solved to obtain a numerical solution to
the differential equation at discrete points corresponding to i =1, 2 and
3.

 Note that uthe system

u1 u 2 of uequation
u 4 has three equations and five
0 3
unknowns , , , ,
u0 u4
 The unknowns and are eliminated using specified
boundary conditions.

 This is discussed in details later.

 Discretization in Time
 In the discussed example, the only independent variable is x.

 This is the case with steady state problems.

 But many petroleum engineering equations are not only space

dependent but also time dependent.

 Hence, we also need to discretize the time by choosing time steps and
obtaining solutions at various times.

 Before discussing discretization in time we need to introduce some

notations.
 Discretization in Time

 We shall use the notation u in

in to mean the value of u at
point i (in space) and time level n.
n
 u 
 Also will  
 t i mean the derivative with respect to
time at point i and time level n

 With this notation, the forward and backward difference

approximations in time at point i are:

 Forward Difference in Time

n n 1
 u  ui  u in
    Ot  5.25
 t  i t
 Discretization in Time

 Backward Difference in Time

n 1 n 1
 u  ui  u in
    Ot  5.26
 t  t

 Equation 5.25 is called the forward difference approximation in time

because to obtain the derivative of u with respect to time at time level
n, the value of u at time level n+1 is used

 Similarly Equation 5.26 is the backward difference approximation

because to obtain the derivative at time level n+1 , use is made of the
value of u at time level n
 Discretization in Time
 Both the forward and backward differences have an error of
order one.

 This can easily be shown by Taylor’s series expansion.

 The implication of using either the forward difference of

backward difference in time is discussed next.
 Explicit and Implicit Formulations
 Assume that we need to solve the equation
5.27
 u
2
u

x 2 t

using a finite difference method

 Using the forward difference in time, Equation 5.27 can be

written as
n n
 u
2
 u  5.28
 2   
 x i  t i
 Discretizing Equation 5.28 becomes
 Explicit and Implicit Formulations

n 1in
u n
 2u  u
n n
u  u in
i 1 i i 1
 i

x  2
t 5.29
 If we define
t
  5.30
x 2
 Equation 5.29 then becomes
5.31
u n 1
i   u  1  2 u   u
n
i 1
n
i
n
i 1

• Equation 5.31 implies that the value of u at point i at time level i+1
can be calculated from values of u at old time level n.

• Hence, this formulation is known as an explicit formulation

 Explicit and Implicit Formulations

 Note that in the explicit formulation, the values of u at

different points at any time level have no relationship with
each other.

 Also in the explicit formulation there is just one equation

and one unknown

 Using the backward difference in time, Equation 5.27 can

be written as:
n 1 n 1
 u
2
 u 
 2    5.32
 x  i  t  i
 Explicit and Implicit Formulations

 Discretizing, Eq. 5.32 becomes

u in11  2u in 1  u in11 u in 1  u in

x 2 t

 Therefore, equation 5.32 can be written as:

u n 1
i 1  1  2 u n 1
i 1  u n 1
i 1  u n
i 3.33

 The values of u at the new time level n+1 are unknown.

 Hence, Equation 5.33 contain three unknowns.

 Solution can only be obtained if the equations for all the points are written
and solved simultaneously.
 Explicit and Implicit Formulations

 If we decide to break the domain in which Equation 5.27 applies into

three, then i in Equation 5.33 will take values of 1, 2, and 3.

 The resulting equations to be solved simultaneously are:

u n 1
0  1  2 u n 1
1  u n 1
2  u n
1

 u1n 1  1  2 u2n1  u3n1  u2n

u n 1
2  1  2 u n 1
3  u n 1
4  u n
3

n 1 n 1
u0 u4
 and are specified boundary conditions.

 Hence, resulting system of equations has three equations and three

 Explicit and Implicit Formulations
 The resulting system of equations can be written in matrix form as
follows:

 1  2 
n 1
 0   u1   u1n  u 0 
      
  1  2   u
 2     u2  n

 0   1  2  u3   u3n  u 4 

 
 u 0n 1 and u 4n 1 are in the right-hand side because they are known.

 Also u n is known.
i

 At the start, n is zero and n corresponds to the initial

ui
condition.
 Explicit and Implicit Formulations
 The matrix is a tridiagonal matrix because it has only one main
diagonal element, one upper diagonal element and one lower diagonal
element.

 This is further discussed in the section on matrix algebra.

 Explicit methods are attractive because they do not require solving a

system of linear equations.

 However, explicit methods are not generally used in reservoir

simulation because of the usually severe restriction on the time step
size.

 The reason for this is given in section on stability.

 Crank-Nicolson Method

 This is a second order method in time obtained by forming

an approximation to  2 u / x 2 atn the
1/ 2 time
u n1  u n / t is central difference ofu / t
level for which
i i

 Using Crank-Nicolson, Equation 5.27 is

n 1 / 2 n 1 / 2
 u
2
 u 
 
2 
  5.35
 x  i  t  i

 Using central difference in time

n 1 / 2
 u  u in 1  u in
 Ot 
2
   5.36
 t  t
 Crank-Nicolson Method
 The left hand side of Equation 5.35 is approximated to

1  u  
n 1 / 2 n 1 n
 u    2 u 

 2 
  
2  x 2  i
  

  O x
2
  5.37
 x  i 
2
 x i  

1  u in11  2u in 1  u in11 u in11  2u in  u n


   i 1

2 x 2 x 2 
5.38

 Combining Equations 5.36 and 5.37 and dropping the truncation error term
5.39
 ui 1  2uin1  uin11  uin1  2uin  uin1   uin1  uin
 n1
2

t
where 
x 2
 Equation 5.39 is the Crank-Nicolson method and it is an implicit
formation.
 Crank-Nicolson Method

 The implicit, explicit and Crank-Nicolson implicit methods can be

written in general form as
n n 1
 u
2
 u
2
u in 1  u in
  2   1    2   5.40
 x  i  x  i t

 where  is the weighting factor that lies within 0 and 1

 For different values of  , we have implicit and Crank-Nicolson

formulation as shown below:
Implicit
 0
Other explicit and
Crank-Nicolson implicit methods are
  1/ 2 discussed by Aziz and
Settari:
 1 Explicit
 Crank-Nicolson Method
Examples
 So far we have not considered the source term or discretization of
equation in two dimensions.

 In this section, some problems and solutions to some of the cases are
given.
 Problem 1: Discretize
 2u  2u
  q
x 2
t
 Solution 1
 Using the explicit formulation
u in1  2u in  u in1 u in 1  u in
  qin
x 2 t
 Problem 11 discretize

 2u  2u u
 
x 2 y 2 t
 Solution 11
 The problem is a two dimensional problem and the grid system to be
used is shown in Figure E.1

 Figure E.1: Two Dimensional Grid

 Using the implicit formulation, the discretized equation is

u in11, j  2u in, j 1  u in11, j u in, j 11  2u in, j 1  u in, j 11 u in, j 1  u in, j
 
x  2
y  2
t

 For simplicity, let us assume that

x  y

and define
t
 
x 2
 The discretized equation becomes

uin11, j  1  4 uin, j 1  uin11, j  uin, j 11  uin, j 11  uin, j

 To obtain the unknowns, we need to write the equation for every
node.
 For a three by three system (that is, I = J = 3), the nodes are shown as
figure E.2

Figure E.2: Nodes for Two Dimensional Problem

• In Figure E.2, the dots represent the nodes where u’s are known while
u’s represent the boundary nodes.

• The values of the independent variable u

 CONSISTENCY, CONVERGENCE AND STABILITY

 Some methods for discretizing a given differential equation have been

discussed.

 The next thing we need to do is to ascertain that the finite difference

solution is close to the actual solution.

 This then brings us to the question of consistency, convergence and

stability of finite different schemes.

 Consistency: A finite-difference approximation is consistent if as the

grid spacing, say gets smaller the truncation error (local
discretization) tends to diminish

 and the finite difference approximation tends to be the true solution

 . This implies that for consistency, the order of the error of the finite-
difference approximation of the derivatives must be greater than zero.

 For example, if  is the order of the error, the following deductions

can be made about finite-difference formulations with the given order
of the error.

  0x  2
consistent
  00 x  consistent
  0x  0
inconsistent
 Convergence:
 Convergence: Let e i be the error in the finite-difference
approximate solution at i, then
e i  U i  ui 6 .1

Where U i  exact solution at point i

u i  finite - difference solution at point i

ei
 The error is called the error of the solution or the global
discretization error.

 A finite difference scheme is said to converge if the error tends

to zero as tends to zero.
 Convergence cont’d
 It has been shown by Aziz and Sattari (1979) that the errors of
solution e i satisfy the difference equation for with the source term
replaced by the truncation error  R i

 For example, the differential equation

d 2u
 q x   0 6.2
dx 2

 can be discretized as

u  2u  u
i 1 i i 1
 q 0 6.3
i
x 2
 Convergence cont’d
 As the global error esatisfies
i Equation 6.3 with source term
q replaced by negative of the truncation error  R , that implies,
i i

e i 1  2 e i  e i 1
 Ri  0 6.4
x 2

 In shorthand form, Equation 6.4 is

1
2 e i  Ri  0 6.5
x  2

 Similarly, for the differential equation

 2u u
  q  0 6.6
x 2 t
 Convergence cont’d
 the explicitly finite difference formulation is
i u in1  2u in  u in1 u in 1  u in
  qin  0 6.7
x  2
t

 Hence, the global error  i satisfies

 in1  2 in   in1   in 1   in
  Ri  0 6.8
x 2 t

 Note that the truncation error Ri in Eqn. (6.8) is the truncation error
due to both space and time discretization

 Since the error equations also contain the truncation error,

investigation on the convergence of a finite-difference scheme is
always difficult
 Convergence cont’d
 This is so because, the truncation errors contain derivatives that
cannot be evaluated or bounded.

 Hence, convergence of a finite-difference scheme can be determined

by investigating consistency and stability of the scheme.

 Both consistency and stability are easier to investigate and some

deductions that may be made about convergence of the scheme are:

 For a consistent approximation, Richtmyer and Morton (1967) stated

that stability is a necessary and sufficient condition for convergence.

 This is quite true for linear differential equations and not true for
some non-linear problems
 Stability
 Stability: A finite-difference scheme is stable if any error introduced
at some stage of computation do not amplify during subsequent
computations.

 The error is made up of discretization errors and round-off errors.

 Quite unlike discretization error, round-off is arbitrary and cannot be

minimized by choosing x and t to be small

 Rather the opposite is the case because smaller mesh sizes and time
steps imply an increase in the number of arithmetic operations and
thus an increase in the round-off error.
 Stability cont’d
 If we define the followings;

U i
n
 true solution
u in*  finite difference approximation that includes both
the discretization and round - off errors
u i  finite difference solution that includes both discretization error ,
the total error at po int i will be given as

 Tin  U in  u in*


 U u i
n n
i   u n
i u n*
i  6.9
Discretization Round-off
 Numerical studies show that the rounding errors are smaller than
 Stability cont’d
 Eqn. (6.9) gives the total error at point i at time level n.

 The total error at point i and time level n+1 is

 n 1
Ti u n 1
i u n 1
i 6.10
 The mathematical expression of the stability criterion is

 n 1
Ti
 1 6.11
 n
Ti
 Eqn. (6.12) implies that the error cannot grow with time. For an
unstable scheme,
 Tin 1
 1 6.12
 Tin
 Stability cont’d
 Fig. 6.1 shows the forms of a stable and an unstable
scheme.

Fig. 6.1: Stable and Unstable Scheme

 Methods for investigating the stability of a scheme is

 Stability cont’d

Von Neumann (Fourier Series) Analysis

 The Von Neumann or Fourier series analysis method is one of the

methods for investigating the stability of a finite-difference scheme.

 The Fourier series method is based on the fact that:

i. The error term satisfies the finite-difference equation with the

truncation error replacing the source terms.

ii. The initial error in the finite-difference approximation can be

expressed as a finite Fourier series. The implication of this fact is
further investigated.
 Stability cont’d
 Using the finite Fourier series, a continuous function f(x) can be
expressed as
N
 n x 
f x    a n cos  6.13
n 0  L 
 Or
N
 n x 
f x    bn sin   6.14
n 0  L 

 depending on the nature of the function.

 But
e i  cos   i sin  6.15
 Stability cont’d
 Hence, a general finite Fourier series representation of the function
is
N
f x    A e inx / L 6.16
n
n0

 where

i  1
 Using Eqn. (6.16), the error at point p (we are using p to avoid
conflict of nomenclature) can be written as

N inpx / L
E   A e 6.17
p n
n0
 Stability cont’d
 The x in Eqn. 6.17 is now replaced by pΔx and p takes values 0, 1, 2,3
……..

 Let us define
n
  6.18
n L

 Then, Eqn. (6.18) becomes

N iB px
E   A e n 6.19
p n
n0
 As the error is additive, we can forget the summation sign and simply
express the error as
iB px
E  A e n 6.20
p n
 Stability cont’d
 The error term we have developed so far is
independent on time.

 Actually the error is also dependent on time except

at time zero.

 The error term that is time dependent will satisfiy

the criterion iB px t
E A e n e 6.21
p n
 And iB px nt
E A e n e 6.22
p n
 Stability cont’d
 In Equation 6.23, n is the time level and takes vaues 0, 1, 2….

 We can drop the constant An in Eqn. 6.22 as it will eventually drop out
in our analysis.

 Hence, the error term for investigation of the stability of a finite

difference scheme in simply written as:

 N
P  c inpx
 n
6.23
 Where  is the amplification factor

 If the scheme is to be stable, the error must not grow with time. Form
Eq 6.24, this implies that
  1 6. 24
 Stability cont’d
 Examples on stability of some finite-difference scheme now follow

 Example 1: Discretize the differential equation:

 2u u
 6.25
x 2
t
 Using the implicit formulation establish the stability criterion for this
Scheme

 The finite-difference equation is

u in11  2u in 1  u in11 u in 1  u in
 6.26
x 2 t
 Stability cont’d
 Let us define
t
  6.27
x 2

 Hence Eq. (6.26) becomes

 u in11  2u in 1  u in11   u in 1  u in 6.28

 The error term satisfies Eq. (6.29)

 Hence

  p 1 n  1  2 np1   p 1 n  1   np1   np 6..29

 Stability cont’d

 Substituting Eq.. 6.23 into Eq. 6.29 gives

 iBn  p  1x  n 1  2iBnpx  n 1iBn  p 1x  n 1   iBn px  n 1  iBnpx  n

 Dividing both sides by iBn p  n gives


 iBnx   2  iBnx      1
 But
iBn x  iBn x
 
Cos Bn x 
2
 On substitution yields

2 Cos  n x  1    1
 Stability cont’d
 Also,
Z
1  Cos Z  2 sin 2

2
 n x
 On substitution yields  4   Sin   1
2

1
 
 n x
 Thus, 1  4 Sin 2
2

 The term on the right hand side always less that 1,

Thus:   1 always

 This implies that finite-difference formulation is unconditionally

 Stability cont’d
 Example 2: Discretize Eq.(6.25) using the explicit formulation and
establish the stability criterion.

 The finite-difference equation is

 u n
i 1  2u
n
i u n
i 1  u n 1
i u n
i

 The error term satisfies the above eqn., thus

  np 1  2 np  np 1    np1   np
 substituting Eqn. (6.23)

 iBn p 1x  n  2iBpx  n  iBn p 1x  n   iBnpx  n 1  iBnpx  n

 Stability cont’d
 Dividing both sides by iBnx  n gives

 iBnx  2  iBnx     1

 But, iBnx   iBnx

Cos Bn x 
2
 On substitution yields  n x
  1  4  Sin 2

2
 For the scheme to be stable

  1
 This implies that
 n x
 1  1  4  Sin 2
 1
2
 Stability cont’d
 The above Equation gives
 n x
0  4  Sin 2
 2
2

1
 And hence  
 n x
2 Sin 2
2

 n x
 The maximum value Sin 2 can take is unity. Hence the criterion
2

1
 
2
 Will guarantee the stability of the scheme.

 This implies that the explicit formulation is conditionally stable.

 Stability cont’d
 Example 3: Establish the stability criterion for the implicit
formulation to be used in solving the equation
uc c

x t

 Discretizing, we have

cin11  cin 1 cin 1  cin

u 
x t
 Let
ut

x
 We have that

 cin11    1cin 1  cin

 Stability cont’d
 The error satisfies the above equation, hence
  np11    1 np1    np

 Therefore,

 iBn  p 1x  n 1    1 iBnpx  n 1    iBnpx  n

 Dividing both sides by  iBnpx n

 C iBnx     1   1

 Simplifying, we have that:

1
 
1      iBn x
 Stability cont’d
 For stability
  1

 Therefore,
1
 1
1     iBnx

 Thus, 1
 1
1     c iBnx

 Using the identity,

1      Cos   i  Sin  1
Stability cont’d