TM 4112

Karakterisasi dan Permodelan Reservoir

Laporan Kemajuan #2

Nama : Maria Indira Puspita Sari

NIM : 12215091
Tanggal Penyerahan : 6 November 2018
Dosen : Dr. Eng. Ir. Sutopo, M.Eng
Asisten : Wilson Wiranda (12214008)

PROGRAM STUDI TEKNIK PERMINYAKAN

FAKULTAS TEKNIK PERTAMBANGAN DAN PERMINYAKAN
INSTITUT TEKNOLOGI BANDUNG
2018
FLOW EQUATION DERIVATION

CONTINUITY EQUATION

Start by considerting the mass balance as shown in graph below:

Mass Balance Equation:

Rate In Rate Out

[𝑟𝑎𝑡𝑒 𝑖𝑛] − [𝑟𝑎𝑡𝑒 𝑜𝑢𝑡] = [𝑟𝑎𝑡𝑒 𝑠𝑖𝑛𝑘/𝑠𝑜𝑢𝑟𝑐𝑒] + [𝑟𝑎𝑡𝑒 𝑎𝑐𝑐𝑢𝑚𝑢𝑙𝑎𝑡𝑖𝑜𝑛]

𝑚̇𝑖𝑛 − 𝑚̇𝑜𝑢𝑡 = 𝑚̇𝑠𝑖𝑛𝑘 + 𝑚̇𝑎𝑐𝑐

where :

𝑚̇𝑖𝑛 = (𝜌𝑄)𝑥 + (𝜌𝑄)𝑦 + (𝜌𝑄)𝑧 = (𝜌𝐴𝑣)𝑥 + (𝜌𝐴𝑣)𝑦 + (𝜌𝐴𝑣)𝑧

𝑚̇𝑜𝑢𝑡 = (𝜌𝑄)𝑥+∆𝑥 + (𝜌𝑄)𝑦+∆𝑦 + (𝜌𝑄)𝑧+∆𝑧 = (𝜌𝐴𝑣)𝑥+∆𝑥 + (𝜌𝐴𝑣)𝑦+∆𝑦 + (𝜌𝐴𝑣)𝑧+∆𝑧

∆(𝜌𝜙𝑆)
𝑚̇𝑎𝑐𝑐𝑢𝑚𝑢𝑙𝑎𝑡𝑒𝑑 = 𝑉𝑏
∆𝑡

thus, substitute those to the mass balance equation, we get

∆(𝜌𝜙𝑆)
[(𝜌𝐴𝑣)𝑥 + (𝜌𝐴𝑣)𝑦 + (𝜌𝐴𝑣)𝑧 ] − [(𝜌𝐴𝑣)𝑥+∆𝑥 + (𝜌𝐴𝑣)𝑦+∆𝑦 + (𝜌𝐴𝑣)𝑧+∆𝑧 ] − 𝑚𝑠𝑖𝑛𝑘 = 𝑉𝑏
∆𝑡

∆(𝜌𝑣𝑥 ) ∆(𝜌𝑣𝑦 ) ∆(𝜌𝑣𝑧 ) ∆(𝜌𝜙𝑆)

−𝑉𝑏 [ + + ] − 𝑚𝑠𝑖𝑛𝑘 = 𝑉𝑏
∆𝑥 ∆𝑦 ∆𝑧 ∆𝑡

𝜕(𝜌𝑣𝑥 ) 𝜕(𝜌𝑣𝑦 ) 𝜕(𝜌𝑣𝑧 ) 𝜕(𝜌𝜙𝑆)

−𝑉𝑏 [ + + ] − 𝑚𝑠𝑖𝑛𝑘 = 𝑉𝑏
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑡

Then, we obtain the continuity equation for black oil:

 𝜕(𝜌𝜙𝑆)
−𝑉𝑏 ∇ ∙ (𝜌𝑣⃗) = 𝑚𝑠𝑖𝑛𝑘 + 𝑉𝑏
𝜕𝑡

Fluid Model

In gas reservoir, there are only gas and water:

for Gas:

𝜌𝑔𝑠
𝜌𝑔 =
𝐵𝑔

𝑉𝑔
𝐵𝑔 =
𝑉𝑔𝑠

for Water:
𝑠
𝜌𝑤
𝜌𝑤 =
𝐵𝑤

𝑉𝑤
𝐵𝑤 =
𝑉𝑤𝑠
for sink / source :

𝑚𝑠𝑖𝑛𝑘 = 𝜌 𝑠 𝑞

Based on the EOS stated above, the continuity equation for each phase become:

for Gas:

𝜌𝑔𝑠 𝜌𝑜𝑠 𝜕 𝜌𝑔𝑠 𝜌𝑜𝑠

−𝑉𝑏 ∇ ∙ ( 𝑣⃗𝑔 ) + 𝑅𝑠 (−𝑉𝑏 ∇ ∙ ( 𝑣⃗𝑜 )) = 𝜌𝑔𝑠 𝑞𝑔 + 𝑅𝑠 𝜌𝑜𝑠 𝑞𝑜 + 𝑉𝑏 ( 𝜙𝑆𝑔 ) + 𝑅𝑠 (( 𝜙𝑆𝑜 ))
𝐵𝑔 𝐵𝑜 𝜕𝑡 𝐵𝑔 𝐵𝑜

for Water:

𝑠 𝑠
𝜌𝑤 𝑠
𝜕 𝜌𝑤
−𝑉𝑏 ∇ ∙ ( 𝑣⃗𝑤 ) = 𝜌𝑤 𝑞𝑤 + 𝑉𝑏 ( 𝜙𝑆𝑤 )
𝐵𝑤 𝜕𝑡 𝐵𝑤

Then, divide both arms by surface density.

for Gas:

𝑣⃗𝑔 𝜕 𝜙𝑆𝑔
−𝑉𝑏 ∇ ∙ ( ) = 𝑞𝑔 + 𝑉𝑏 ( )
𝐵𝑔 𝜕𝑡 𝐵𝑔

for Water:
 𝑣⃗𝑤 𝜕 𝜙𝑆𝑤
−𝑉𝑏 ∇ ∙ ( ) = 𝑞𝑤 + 𝑉𝑏 ( )
𝐵𝑤 𝜕𝑡 𝐵𝑤

MOMENTUM BALANCE

Here is Darcy Equation for both phases using the relative permeability:

𝐾𝑘𝑟𝑛
𝑣⃗𝑛 = − (∇𝑃𝑛 − 𝜌𝑛 𝑔∇𝑧)
𝜇𝑛
𝑣⃗𝑛
= −𝜆𝑛 (∇𝑃𝑛 − 𝛾𝑛 ∇𝑧)
𝐵𝑛

Where,
𝐾𝑘𝑟𝑛
𝜆𝑛 =
𝜇𝑛 𝐵𝑛
𝛾𝑛 = 𝜌𝑛 𝑔
𝑛 = 𝑜, 𝑔

With assumption that reservoir is modeled with grid in the form of horizontal block:

𝜕/𝜕𝑥 𝜆𝑥 𝑝𝑥
𝑣⃗
−∇ ∙ ( ) = ∇ ∙ [𝜆(∇𝑝 − 𝛾∇𝑧)] = [𝜕/𝜕𝑦] ∙ [ 𝜆𝑦 𝑝𝑦 ]
𝐵
𝜕/𝜕𝑧 𝜆𝑧 (𝑝𝑧 − 𝛾∇𝑧)

𝑣⃗ 𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕
−∇ ∙ ( ) = (𝜆𝑥 ) + (𝜆𝑦 ) + (𝜆𝑧 ) − (𝜆𝑧 𝛾∇𝑧)
𝐵 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑧

MULTIPHASE FLOW EQUATION

Combining the continuity equation for each phase with momentum balance equation from
Darcy’s Equation, The equation for each phase becomes:

for Gas

𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕 𝜕 𝜙𝑆𝑔
𝑉𝑏 [ (𝜆𝑔𝑥 ) + (𝜆𝑔𝑦 ) + (𝜆𝑔𝑧 ) − (𝜆𝑔𝑧 𝛾𝑔 ∇𝑧)] = 𝑞𝑔 + 𝑉𝑏 ( )
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑧 𝜕𝑡 𝐵𝑔

for Water:

𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕 𝜕 𝜙𝑆𝑤
𝑉𝑏 [ (𝜆𝑤𝑥 ) + (𝜆𝑤𝑦 ) + (𝜆𝑤𝑧 ) − (𝜆𝑤𝑧 𝛾𝑤 ∇𝑧)] = 𝑞𝑤 + 𝑉𝑏 ( )
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑧 𝜕𝑡 𝐵𝑤

Unknown Parameters:

1) Saturation (Sw and Sg)

2) Pressure (Pw and Pg).
The additional equation necessary for solving the multiphase flow equation are :

𝑆𝑔 + 𝑆𝑤 = 1

𝑝𝑐𝑔𝑤 = 𝑝𝑔 − 𝑝𝑤

DISCRETIZATION OF DIFUSIVITY EQUATION

 LEFT HAND SIDE

Left term of of gas equation ( the transmissibility term ) is

𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕
𝑉𝑏 [ (𝜆𝑔𝑥 ) + (𝜆𝑔𝑦 ) + (𝜆𝑔𝑧 ) − (𝜆𝑔𝑧 𝛾𝑔 ∇𝑧)]
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑧

For the z-direction :

𝜕 𝜕𝑝
𝑉𝑏 . (𝜆𝑔𝑧 ) ≡ 𝑉𝑏 . ∆𝑧 (𝜆𝑔𝑧 . ∆𝑧 𝑝𝑔 )
𝜕𝑧 𝜕𝑧

Since

𝑉𝑏 = ∆𝑧𝑘 . 𝐴𝑘

Therefore,

𝑉𝑏 . ∆𝑧 (𝜆𝑔𝑧 . ∆𝑧 𝑝𝑔 ) = ∆𝑧 (∆𝑧𝑘 . 𝐴𝑘 . 𝜆𝑔𝑧 . ∆𝑧 𝑝𝑔 )

Using concept of central difference in space :

1
∆𝑧 (∆𝑧𝑘 . 𝐴𝑘 . 𝜆𝑔𝑧 . ∆𝑧 𝑝𝑔 ) = [(∆𝑧𝑘 . 𝐴𝑘 . 𝜆𝑔𝑧 . ∆𝑧 𝑝𝑔 )𝑘+1 − (∆𝑧𝑘 . 𝐴𝑘 . 𝜆𝑔𝑧 . ∆𝑧 𝑝𝑔 )𝑘−1 ]
∆𝑧𝑘 2 2

= (𝐴𝑘 . 𝜆𝑔𝑧 . ∆𝑧 𝑝𝑔 ) 1 −(𝐴𝑘 . 𝜆𝑔𝑧 . ∆𝑧 𝑝𝑔 ) 1

𝑘+ 𝑘−
2 2

Using the concept of difference in the pressure term:

 𝑝𝑔,𝑘+1 − 𝑝𝑔,𝑘 𝑝𝑔,𝑘 − 𝑝𝑔,𝑘−1
= (𝐴𝑘 . 𝜆𝑔𝑧 . ) 1 − (𝐴𝑘 . 𝜆𝑔𝑧 . ) 1
∆𝑧 1 𝑘+
2 ∆𝑧 1 𝑘−
2
𝑘+ 𝑘−
2 2

Transmisibility,

𝐴𝑘
𝑇𝑍 1 = 𝜆 1
𝑘+
2 ∆𝑧 1 𝑔,𝑘+2
𝑘+
2

Therefore,

For x-direction

𝜕 𝜕𝑝
𝑉𝑏 (𝜆 ) ≈ 𝑇𝑋 1 (𝑝𝑖+1 − 𝑝𝑖 ) − 𝑇𝑋 1 (𝑝𝑖 − 𝑝𝑖−1 )
𝜕𝑥 𝑥 𝜕𝑥 𝑖+
2
𝑖−
2

For y-direction

𝜕 𝜕𝑝
𝑉𝑏 (𝜆𝑦 ) ≈ 𝑇𝑌 1 (𝑝𝑗+1 − 𝑝𝑗 ) − 𝑇𝑌 1 (𝑝𝑗 − 𝑝𝑗−1 )
𝜕𝑦 𝜕𝑦 𝑗+
2
𝑗−
2

For z-direction:

𝜕 𝜕𝑝
𝑉𝑏 (𝜆𝑧 ) ≈ 𝑇𝑍 1 (𝑝𝑘+1 − 𝑝𝑘 ) − 𝑇𝑍 1 (𝑝𝑘 − 𝑝𝑘−1 )
𝜕𝑧 𝜕𝑧 𝑘+
2
𝑘−
2

For gravity term,

𝜕 𝜕 𝜕𝑧
(𝜆𝑤𝑧 𝛾𝑤 ∇𝑧) = (𝜆𝑧 𝛾𝑧 )
𝜕𝑧 𝜕𝑧 𝜕𝑧

𝑧𝑘+1 − 𝑧𝑘 𝑧𝑘 − 𝑧𝑘−1
𝜆 1(
∆𝑧 ) − 𝜆𝑧,𝑘−12 ( ∆𝑥 )
𝜕 𝜕𝑧 𝑧,𝑘+
2
(𝜆𝑧 𝛾𝑧 ) ≈
𝜕𝑧 𝜕𝑧 ∆𝑧

𝜕 𝜕𝑧 𝑧𝑖+1 − 𝑧𝑖 𝑧𝑘 − 𝑧𝑘−1
𝑉𝑏 (𝜆𝑧 𝛾𝑧 ) ≈ 𝜆 1 ∆𝑥∆𝑦 ( ) − 𝜆 1 ∆𝑥∆𝑦 ( )
𝜕𝑧 𝜕𝑧 𝑘+
2 ∆𝑧 𝑘−
2 ∆𝑧
Submitting again the term of transmissibility,

∆𝑥∆𝑦
𝑇𝑍 1 = 𝜆 1
𝑘+
2 ∆𝑧 𝑘+2

𝜕 𝜕𝑧
𝑉𝑏 (𝜆𝑧 𝛾𝑧 ) ≈ 𝑇𝑍 1 𝛾 1 [𝑧𝑘+1 − 𝑧𝑘 ] − 𝑇𝑍 1 𝛾 1 [𝑧𝑘 − 𝑧𝑘−1 ]
𝜕𝑧 𝜕𝑧 𝑘+ 𝑘+
2 2
𝑘− 𝑘−
2 2

Thus, the complete left hand side can be discretized as :

𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕
𝑉𝑏 [ (𝜆𝑔𝑥 ) + (𝜆𝑔𝑦 ) + (𝜆𝑔𝑧 ) − (𝜆𝑔𝑧 𝛾𝑔 ∇𝑧)]
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑧
= 𝑇𝑋 1 (𝑝𝑔,𝑖+1 − 𝑝𝑔,𝑖 ) − 𝑇𝑋𝑔,𝑖−1 (𝑝𝑔,𝑖 − 𝑝𝑔,𝑖−1 )
𝑔,𝑖+
2 2

+ 𝑇𝑌 1 (𝑝𝑔,𝑗+1 − 𝑝𝑔,𝑗 ) − 𝑇𝑌𝑔,𝑗−1 (𝑝𝑔,𝑗 − 𝑝𝑔,𝑗−1 )

𝑔,𝑗+
2 2

+ 𝑇𝑍 1 (𝑝𝑔,𝑘+1 − 𝑝𝑔,𝑘 ) − 𝑇𝑍𝑔,𝑘−1 (𝑝𝑔,𝑘 − 𝑝𝑔,𝑘−1 )

𝑔,𝑘+
2 2

+ 𝑇𝑍 1𝛾 1 [𝑧𝑘+1 − 𝑧𝑘 ] − 𝑇𝑍 1𝛾 1 [𝑧𝑘 − 𝑧𝑘−1 ]

𝑔,𝑘+ 𝑔,𝑘+ 𝑔,𝑘− 𝑔,𝑘−
2 2 2 2

𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕 𝜕𝑝 𝜕
𝑉𝑏 [ (𝜆𝑔𝑥 ) + (𝜆𝑔𝑦 ) + (𝜆𝑔𝑧 ) − (𝜆𝑔𝑧 𝛾𝑔 ∇𝑧)]
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑧

= 𝑇𝑋 1 𝑝𝑔,𝑖+1 − [𝑇𝑋 1 + 𝑇𝑋 1 ] 𝑝𝑔,𝑖 + 𝑇𝑋 1 𝑝𝑔,𝑖−1

𝑔,𝑖+ 𝑔,𝑖+ 𝑔,𝑖− 𝑔,𝑖−
2 2 2 2

+ 𝑇𝑌 1 𝑝𝑔,𝑗+1 − [𝑇𝑌 1 + 𝑇𝑌 1 ] 𝑝𝑔,𝑗 + 𝑇𝑌 1 𝑝𝑔,𝑗−1

𝑔,𝑗+ 𝑔,𝑗+ 𝑔,𝑗− 𝑔,𝑗−
2 2 2 2

+ 𝑇𝑍 1 𝑝𝑔,𝑘+1 − [𝑇𝑍 1 + 𝑇𝑍 1 ] 𝑝𝑔,𝑘 + 𝑇𝑍 1 𝑝𝑔,𝑘−1

𝑔,𝑘+ 𝑔,𝑘+ 𝑔,𝑘− 𝑔,𝑘−
2 2 2 2

+ +𝑇𝑍 1𝛾 1 [𝑧𝑘+1 − 𝑧𝑘 ] − 𝑇𝑍 1𝛾 1 [𝑧𝑘 − 𝑧𝑘−1 ]

𝑔,𝑘+ 𝑔,𝑘+ 𝑔,𝑘− 𝑔,𝑘−
2 2 2 2

with :

∆𝑦∆𝑧 ∆𝑦∆𝑧 𝑘. 𝑘𝑟𝑔

𝑇𝑋 1 = 𝜆 1= [ ]
𝑔,𝑖±
2 ∆𝑥 1 𝑖±
2 ∆𝑥 1 𝜇𝑔 . 𝐵𝑔 𝑖±1
𝑖± 𝑖± 2
2 2

∆𝑥∆𝑧 ∆𝑥∆𝑧 𝑘. 𝑘𝑟𝑔

𝑇𝑌 1 = 𝜆 1= [ ]
𝑔,𝑗±
2 ∆𝑦 1 𝑗±2 ∆𝑦 1 𝜇𝑔 . 𝐵𝑔 1
𝑗± 𝑗± 𝑗±
2 2 2

∆𝑥∆𝑦 ∆𝑥∆𝑦 𝑘. 𝑘𝑟𝑔

𝑇𝑍 1 = 𝜆 1= [ ]
𝑔,𝑘±
2 ∆𝑧 1 𝑘±
2 ∆𝑧 1 𝜇𝑔 . 𝐵𝑔 𝑘±1
𝑘± 𝑘± 2
2 2
  RIGHT HAND TERM
The right term for this equation is:

𝜕 𝜙𝑆𝑔
( ) + 𝑄𝑔
𝜕𝑡 𝐵𝑔
Since LHS is multiplied by Vb, then the RHS is also multiplied by Vb. Using
𝜕 ∅𝑆
backward method, the value of 𝑉𝑏 𝜕𝑡 [ 𝐵 ] is obtained:

𝜕𝑝 𝑝𝑡+1 − 𝑝𝑡
=
𝜕𝑡 ∆𝑡

𝜕 ∅𝑆 𝑡 𝑉𝑏 ∅𝑆 𝑡+1 ∅𝑆 𝑡
𝑉𝑏 [ ] = [( ) −( ) ]
𝜕𝑡 𝐵 ∆𝑡 𝐵 𝐵

Therefore, for gas equation, the RHS becomes:

𝑡 𝑡+1 𝑡
𝜕 𝜙𝑆𝑔 𝑉𝑏 𝜙𝑆𝑔 𝜙𝑆𝑔
𝑉𝑏 [ ] + 𝑄𝑔 = [( ) −( ) ] + 𝑄𝑔
𝜕𝑡 𝐵𝑔 ∆𝑡 𝐵𝑔 𝐵𝑔

By procedure of discretization for RHS and LHS, thus the diffusivity equation for gas and
water becomes:

- GAS EQUATION
𝜕 𝜙𝑆𝑔
𝑉𝑏 ∇ ∙ (𝜆𝑔 (∇𝑃𝑔 − 𝛾𝑔 ∇𝑧)) = 𝑉𝑏 ( ) + 𝑄𝑔
𝜕𝑡 𝐵𝑔

(𝑇𝑋𝑔,𝑖+1 (𝑃𝑔,𝑖+1 − 𝑃𝑔,𝑖 ) + 𝑇𝑋𝑔,𝑖−1 (𝑃𝑔,𝑖−1 − 𝑃𝑔,𝑖 ))

2 2

+ (𝑇𝑌𝑔,𝑗+1 (𝑃𝑔,𝑗+1 − 𝑃𝑔,𝑗 ) + 𝑇𝑌𝑔,𝑗−1 (𝑃𝑔,𝑗−1 − 𝑃𝑔,𝑗 ))

2 2

+ (𝑇𝑍𝑔,𝑘+1 (𝑃𝑔,𝑘+1 − 𝑃𝑔,𝑘 ) + 𝑇𝑍𝑔𝑘−1 (𝑃𝑔,𝑘 − 𝑃𝑔,𝑘−1 ))

2 2

− (𝑇𝑍𝑔,𝑘+1 𝛾𝑔,𝑘+1 (𝑧𝑘+1− 𝑧𝑘 ) − 𝑇𝑍𝑔,𝑘−1 𝛾𝑔,𝑘−1 (𝑧𝑘− 𝑧𝑘−1 ))

2 2 2 2
𝑡+1 𝑡
𝑉𝑏 𝜙𝑆𝑔 𝜙𝑆𝑔
= [( ) −( ) ] + 𝑄𝑔
∆𝑡 𝐵𝑔 𝐵𝑔
 - WATER EQUATION
𝜕 𝜙(1 − 𝑆𝑔 )
𝑉𝑏 ∇ ∙ (𝜆𝑤 (∇𝑃𝑤 − 𝛾𝑤 ∇𝑧)) = 𝑉𝑏 ( ) + 𝑄𝑤
𝜕𝑡 𝐵𝑤

(𝑇𝑋𝑤,𝑖+1 (𝑃𝑤,𝑖+1 − 𝑃𝑤,𝑖 ) + 𝑇𝑋𝑤,𝑖−1 (𝑃𝑤,𝑖−1 − 𝑃𝑤,𝑖 ))

2 2

+ (𝑇𝑌𝑤,𝑗+1 (𝑃𝑤,𝑗+1 − 𝑃𝑤,𝑗 ) + 𝑇𝑌𝑤,𝑗−1 (𝑃𝑤,𝑗−1 − 𝑃𝑤,𝑗 ))

2 2

+ (𝑇𝑍𝑤,𝑘+1 (𝑃𝑤,𝑘+1 − 𝑃𝑤,𝑘 ) + 𝑇𝑍𝑜𝑘−1 (𝑃𝑤,𝑘 − 𝑃𝑤,𝑘−1 ))

2 2

− (𝑇𝑍𝑤,𝑘+1 𝛾𝑤,𝑘+1 (𝑧𝑘+1− 𝑧𝑘 ) − 𝑇𝑍𝑤,𝑘−1 𝛾𝑤,𝑘−1 (𝑧𝑘− 𝑧𝑘−1 ))

2 2 2 2
𝑡+1 𝑡
𝑉𝑏 𝜙(1 − 𝑆𝑔 ) 𝜙(1 − 𝑆𝑔 )
= [( ) −( ) ] + 𝑄𝑤
∆𝑡 𝐵𝑤 𝐵𝑤

JACOBIAN MATRIX FOR RESIDUAL FLOW FUNCTION

Jacobian matrix can be used to find solution for the flow equation of gas and water contains 2
primary variable, water pressure ( or gas pressure) and water saturation (or gas saturation).
- Gas
𝜕 𝜙𝑆𝑔
𝑉𝑏 ∇[𝜆𝑔 (∇𝑝𝑔 − 𝛾𝑔 ∇𝑧)] = 𝑞𝑔 + 𝑉𝑏 ( )
𝜕𝑡 𝐵𝑔

𝜕 𝜙(1 − 𝑆𝑤 )
𝑉𝑏 ∇[𝜆𝑔 (∇(𝑝𝑐𝑔𝑤 + 𝑝𝑤 ) − 𝛾𝑔 ∇𝑧)] = 𝑞𝑔 + 𝑉𝑏 ( )
𝜕𝑡 𝐵𝑔
- Water
𝜕 𝜙𝑆𝑤
𝑉𝑏 ∇[𝜆𝑤 (∇𝑝𝑤 − 𝛾𝑤 ∇𝑧)] = 𝑞𝑤 + 𝑉𝑏 ( )
𝜕𝑡 𝐵𝑤

Using Taylor expansion:

- Gas
𝜕𝐹𝑔 𝜕𝐹𝑔
𝐹𝑔(𝑆𝑤 + ∆𝑆𝑤 , 𝑝𝑤 + ∆𝑝𝑤 ) = 𝐹𝑔(𝑆𝑤 , 𝑝𝑤 ) + ∆𝑆𝑤 + ∆𝑝𝑤
𝜕𝑝𝑤 𝜕𝑆𝑤

- Water

𝜕𝐹𝑤 𝜕𝐹𝑤
𝐹𝑤(𝑆𝑤 + ∆𝑆𝑤 , 𝑝𝑤 + ∆𝑝𝑤 ) = 𝐹𝑤(𝑆𝑤 , 𝑝𝑤 ) + ∆𝑆𝑤 + ∆𝑝𝑤
𝜕𝑝𝑤 𝜕𝑆𝑤

The first guess for the New Rhapson solving method gives:
 𝐹𝑔(𝑆𝑤 + ∆𝑆𝑤 , 𝑝𝑤 + ∆𝑝𝑤 ) = 𝐹𝑤(𝑆𝑤 + ∆𝑆𝑤 , 𝑝𝑤 + ∆𝑝𝑤 ) = 0
Then, the equations above become:
- Gas
𝜕𝐹𝑔 𝜕𝐹𝑔
−𝐹𝑔(𝑆𝑤 , 𝑝𝑤 ) = ∆𝑆𝑤 + ∆𝑝𝑤
𝜕𝑝𝑤 𝜕𝑆𝑤
- Water
𝜕𝐹𝑤 𝜕𝐹𝑤
−𝐹𝑤(𝑆𝑤 , 𝑝𝑤 ) = ∆𝑆𝑤 + ∆𝑝𝑤
𝜕𝑝𝑤 𝜕𝑆𝑤

The system of the equations of gas and water then can be represented in the form of Jacobian
Matrix:

𝜕𝐹𝑔 𝜕𝐹𝑔
𝜕𝑝𝑤 𝜕𝑆𝑤 ∆𝑆𝑤 𝐹𝑔(𝑆𝑤 , 𝑝𝑤 )
[ ] = −[ ]
𝜕𝐹𝑤 𝜕𝐹𝑤 ∆𝑝𝑤 𝐹𝑤(𝑆𝑤 , 𝑝𝑤 )
[ 𝜕𝑝𝑤 𝜕𝑆𝑤 ]

𝜕𝐹𝑔 𝜕𝐹𝑔 −1
∆𝑆 𝜕𝑝𝑤 𝜕𝑆𝑤 𝐹𝑔(𝑆𝑤 , 𝑝𝑤 )
[ 𝑤] = − [ ]
∆𝑝𝑤 𝜕𝐹𝑤 𝜕𝐹𝑤 𝐹𝑤(𝑆𝑤 , 𝑝𝑤 )
[ 𝜕𝑝𝑤 𝜕𝑆𝑤 ]

𝜕𝐹𝑤 𝜕𝐹𝑔
−
𝜕𝑆𝑤 𝜕𝑆𝑤 𝐹𝑔(𝑆𝑤 , 𝑝𝑤 )
−[ ][ ]
𝜕𝐹𝑤 𝜕𝐹𝑔 𝐹𝑤(𝑆𝑤 , 𝑝𝑤 )
−
∆𝑆𝑤 𝜕𝑝𝑤 𝜕𝑝𝑤
[ ]=
∆𝑝𝑤 𝜕𝐹𝑔 𝜕𝐹𝑔
𝜕𝑝 𝜕𝑆𝑤
| 𝑤 |
𝜕𝐹𝑤 𝜕𝐹𝑤
𝜕𝑝𝑤 𝜕𝑆𝑤

𝜕𝐹𝑤 𝜕𝐹𝑔
−
𝜕𝑆𝑤 𝜕𝑆𝑤 𝐹𝑔(𝑆𝑤 , 𝑝𝑤 )
[ ][ ]
𝜕𝐹𝑤 𝜕𝐹𝑔 𝐹𝑤(𝑆𝑤 , 𝑝𝑤 )
−
∆𝑆 𝜕𝑝𝑤 𝜕𝑝𝑤
[ 𝑤] =
∆𝑝𝑤 𝜕𝐹𝑔 𝜕𝐹𝑤 𝜕𝐹𝑔 𝜕𝐹𝑤
−
𝜕𝑝𝑤 𝜕𝑆𝑤 𝜕𝑆𝑤 𝜕𝑝𝑤

Solving the problem for each primary variable will give:

 𝜕𝐹𝑔 𝜕𝐹𝑤 𝜕𝐹𝑤 𝜕𝐹𝑔
𝐹𝑤− 𝐹𝑔 𝐹𝑔 − 𝐹𝑤
𝜕𝑆𝑤 𝜕𝑆𝑤 𝜕𝑝𝑤 𝜕𝑝𝑤
∆𝑆𝑤 = 𝑎𝑛𝑑 ∆𝑝𝑤 =
𝜕𝐹𝑔 𝜕𝐹𝑤 𝜕𝐹𝑔 𝜕𝐹𝑤 𝜕𝐹𝑔 𝜕𝐹𝑤 𝜕𝐹𝑔 𝜕𝐹𝑤
− −
𝜕𝑝𝑤 𝜕𝑆𝑤 𝜕𝑆𝑤 𝜕𝑝𝑤 𝜕𝑝𝑤 𝜕𝑆𝑤 𝜕𝑆𝑤 𝜕𝑝𝑤

WORKFLOW FOR SOLVING THE PROBLEM IN THE PROJECT

Print the Time =

results (10) tmax C
Start

Evaluate mat.
Read data (1) Balance error (9)
Inputdata.f90 Mbal_error.f90
End
C

Intialize (2)
Check
Init_con.f90 convergence

Check the flow

direction (3)

Poten.f90 Update primary

variable (7)

Time = time + dt

Swn (I, j, k) = swk (I, j k)

Solve for primary
Pn (I, j, k) = Pk (I, j, k) variable (6)

Solver.f90

Evaluate sink/ Evaluate

source tem Jacobian

Well.f90 Jacobian.f90