Ikatan Ahli Teknik Perminyakan Indonesia Ikatan Ahli Teknik Perminyakan Indonesia Simposium Nasional IATMI 2009 Bandung, 2-5 Desember 2009
Makalah Profesional IATMI 09 – 021 Prediksi Pola Aliran Multifasa untuk Tiap Segmen pada Jaringan Pipa Kompleks di Lapangan X 2,3
Ucok W.R. Siagian,2,3Leksono Mucharam, Lala S. Riza, 1,3Lira Adiyanti, 2,3Ardian PP
3,4
1)
2)
Program Studi Matematika, Institut Teknologi Bandung Program Studi Teknik Perminyakan, Institut Teknologi Bandung 3) Research Consortium OPPINET, Institut Teknologi Bandung 4) Program Ilmu Komputer, Universitas Pendidikan Indonesia
Abstrak Distribusi tekanan merupakan salah satu pertimbangan dalam perancangan suatu sistem jaringan pipa gas-kondensat. Prediksi mengenai distribusi tekanan dapat dilakukan dengan simulasi komputer dengan menerapkan sejumlah model aliran dalam suatu algoritma numerik. Pada penelitian ini dikembangkan suatu model simulasi untuk prediksi distribusi tekanan pada suatu sistem jaringan pipa gas-kondesat. Model simulasi yang dikembangkan pada penelitian ini terdiri dari dua bagian. Bagian pertama adalah linearisasi model aliran dua fasa Beggs and Brill’s. Bagian kedua adalah pengembangan sistem persamaan berdasarkan susunan jaringan pipa dan metode Kirchoff. Persamaan matematik dari pemodelan ini diselesaikan secara numerik dengan menggunakan metode iterasi Newton. Hasil penelitian ini diaplikasikan untuk prediksi distribusi tekanan jaringan pipa pada kasus Lapangan X. Distribusi tekanan yang diprediksi dengan model hasil penelitian ini sesuai dengan hasil prediksi menggunakan commercial software. Perbandingan hasil tersebut telah menunjukkan kesesuaian. Kata kunci : Flow Pattern, Beggs and Brill’s, Metode Iterasi Newton.
Pendahuluan Pada saat ini, konsumsi minyak bumi di dunia meningkat dengan pesat, seiring dengan meningkatnya permintaan energi secara global. Hal ini juga mempengaruhi harga minyak bumi yang
cenderung tidak stabil jika dibandingkan dengan harga gas bumi yang bergantung pada kontrak. Beberapa negara berencana untuk mensubsitusi minyak bumi dengan energi alternatif yang memiliki cadangan yang lebih banyak, salah satunya adalah gas bumi. Indonesia merupakan salah satu negara yang ikut mengambil langkah ini. Tren produksi minyak bumi di Indonesia cenderung menurun, tetapi produksi gas buminya justru meningkat. Selama beberapa tahun, produksi gas bumi di Indonesia meningkat secara kontinu. Pada tahun 2004, rata-rata produksi gas bumi Indonesia adalah 8.35 BSCFD dan konsumsi gas bumi dalam negeri adalah 4.88 BSCFD sedangkan sisanya, 3.47 BSCFD diekspor ke berbagai negara. Cetak biru Kebijakan Energi Nasional 2005 – 2025, mengharapkan penggunaan energi nasional mencapai titik optimum, yaitu 26.2% untuk minyak bumi, 30.6% untuk gas alam, dan 3.8% untuk energi geothermal. Berdasarkan kasus ini, Research Consortium OPPINET (Optimization on Gas and Oil Transmission and Distribution Pipeline Network) di bawah Pusat Pemodelan Matematika dan Simulasi, sebagai salah satu pusat penelitian di ITB, kolaborasi berbagai bidang ilmu, ikut mendukung pencapaian nasional dalam bidang minyak dan gas bumi. Optimasi produksi gas menjadi topik penelitian yang diharapkan dapat meningkatkan produksi gas. Metode optimasi produksi dikembangkan dari pemodelan matematika dengan membangun simulasi sesuai dengan permasalahan di lapangan. Pendekatan ini dipilih karena biaya yang terjangkau dan dapat direpresentasikan dalam kasus lapangan.
___________________________________________________________________________________________________
1
Ikatan Ahli Teknik Perminyakan Indonesia Sifat fasa
Pengembangan Model Model yang digunakan dalam penelitian ini terdiri dari 3 komponen utama, yaitu model aliran duafasa (metode linier Beggs-Brill), model analisa jaringan Hukum Kirchoff, dan prosedur iterasi. Dalam mengembangan model, digunakan beberapa asumsi, yaitu : 1. Kondisi aliran steady state. 2. Komposisi fluida dan temperatur konstan sepanjang sistem jaringan. 3. Minor lossess diabaikan.
Model aliran dua fasa Model aliran dua fasa yang digunakan adalah metode Beggs-Brill. Dalam penelitian ini, metode Beggs-Brill digunakan untuk memprediksi performa aliran gas kondensat dalam pipa. Model ini akan digunakan dalam suatu sistem jaringan pipa kompleks. Metode tersebut dapat ditulis sebagai berikut :
Wmn +1 = η i ΔP n +1 dimana
ηi
Perhitungan sifat fasa dihampiri oleh perhitungan yang melibatkan berbagai sumber. Sifatsifat fluida yang perlu dikalkulasi adalah densitas dan viskositas fluida (gas, kondensat, dan air). Pertama, akan dihitung faktor kompresibilitas gas. Sifat pseudokritikal gas didefinisikan sebagai (R.Sutton):
Ppc = 756.8 − 131γ g − 3.6γ g2
(7)
Tpc = 169.2 + 349.5γ g − 74γ g2
(8)
Jika terdapat ketidakmurnian dalam gas, maka nilai pseudokritical perlu dikoreksi dengan menggunakan korelasi Wichert-Azis
(
∈ = 120 A − A 0.9
1.6
⎧ C1 ⎩ Wm
P =
⎫
g gc
ρ m sin θ
(3)
C2 =
(2)
Tpr =
(3)
Ppr =
gf m (5)
(
)
2
(10)
T Tpc
(11)
P
(12)
Ppc
z
= 1 + A11 ρ pr + A22 ρ pr − A33 ρ pr 2
5
(
+ 0.6134 1 + 0.721ρ pr
4v sg
π d 2 gc P
Untuk jaringan pipa yang tidak memuat loop, kecepatan massa ( Wm ) pada setiap pipa adalah konstan untuk setiap waktu. Jadi, tidak dibutuhkan prosedur iteratif untuk menentukan Wm sehingga
2
⎛ ρ pr2 ⎞ ) ⎜ T 2 ⎟ exp ( −0.721ρ pr2 ) ⎝ pr ⎠
(13)
dimana A1 1
=
⎛ 1 .0 7 ⎞ ⎛ 0 .5 3 3 9 ⎞ ⎟ − ⎜ ⎟ 3 ⎝ T pr ⎠ ⎝ T pr ⎠
0 .3 2 6 5 − ⎜
⎛ 0 .0 1 5 6 9 ⎞ ⎛ 0 .0 5 1 6 5 ⎞ ⎟ − ⎜ ⎟ T p4r T p5r ⎝ ⎠ ⎝ ⎠
+ ⎜
dapat ditulis sebagai :
(Wm )i = ηi ΔPi
(9)
Faktor kompresibilitas gas dapat ditulis sebagai:
(4)
C3 =
)
Sifat pseudoreduced dapat dihitung dengan
⎭
(4)
g cπ 2 d 5 ρ m
4
Tpc + y H S 1- yH S ∈ 2
dan
C1 =
PpcTpc'
' pc
+ C2 Wm ⎬
−B
Tpc' = Tpc − ∈
(1)
1 − C3 Wm ΔZ ⎨
0.5
dimana A = Jumlah fraksi mol CO2 dan H2S B = Fraksi mol N2
didefinisikan sebagai:
ηi =
) + 15 ( B
⎛ 0 .7 3 6 1 ⎞ ⎛ 0 .1 8 4 4 ⎞ ⎟ + ⎜ ⎟ 2 ⎝ T pr ⎠ ⎝ T pr ⎠
(6)
A 22
=
0 .5 4 7 5 − ⎜
Pressure drop pada setiap pipa dalam jaringan dapat diselesaikan dengan menggunakan teknik solusi matriks. Dalam penelitian ini, kami menggunakan eleminasi Gauss. Matriks persamaan dapat diperoleh dari persamaan (6)
A33
=
0 .1 0 5 6 ⎜
ρ
=
0 .2 7 ⎜
pr
⎛ − 0 .7 3 6 1 ⎝
⎛
Ppr
T pr
+
0 .1 8 4 4 ⎞ T p2r
⎟ ⎠
⎞ ⎟
⎝ z T pr ⎠
___________________________________________________________________________________________________
2
Ikatan Ahli Teknik Perminyakan Indonesia ρo = γ oρw
Densitas gas Berikut adalah persamaan yang digunakan untuk menghitung kepadatan gas: M a = 29γ g (14)
Viskositas kondensat Untuk menentukan viskositas dari kondensat, digunakan korelasi Glasso:
Ma
ρg =
( )
= ⎡⎣3.141 1010 ⎤⎦ (T − 460)−3.444 [ log( API )]
μo
(15)
10.732 z T
μ g 1 = μ g 1uncorrected + ( Δμ ) H S + ( Δμ ) N + ( Δμ )CO 2
2
Densitas air dihitung persamaan berikut:
= yN
( Δμ )CO
= yCO
2
−3
S = Kadar garam dalam persen berat.
Viskositas air
−3
g
Viskositas air diperoleh dari persamaan berikut:
−3
μ w1 = A ( T − 460 )
g
−3
−3
g
2
Viskositas gas pada menggunakan korelasi Dempsey,
μ g = μ g1
tekanan
tinggi
exp ( C11 )
(17)
Tpr
B
(23)
dimana A = A0 + A1 S + A2 S 2 + A3 S 3 A0 = 109.574
A1 = −8.40564
A2 = 0.313314
A3 = 8.72213 10 −3
(
)
dan
dimana
B
C11
(22)
g
−3
2
2
menggunakan
dimana
( ) log γ ⎡⎣8.48(10 ) log γ + 9.59(10 ) ⎤⎦ ⎡⎣9.08(10 ) log γ + 6.24(10 ) ⎤⎦ ⎡⎣8.49(10 ) log γ + 3.73(10 ) ⎤⎦
( Δμ ) N
2
dengan
ρ w = 62.368 + 0.438603 S + 1.60074 (10−3 ) S 2
−3
6.15 10
(21)
Densitas air
(16)
2
(20)
=10.313[ log(T − 460)] − 36.447
a
dimana μg1uncorrected = ⎡⎣1.79(10−5 ) − 2.062(10−6 ) γ g ⎤⎦ (T - 460) + 8.188(10−3 ) − = yH S
(19)
⎛ 141.5 ⎞ ⎟ −131.5 ⎝ γo ⎠
Perhitungan viskositas gas berdasarkan korelasi Standing, dengan mempertimbangkan efek dari viskositas molekul nonhydrokarbon pada tekanan atmosfer.
2
a
API = ⎜
Viskositas gas
( Δμ )H S
(18)
(
B2 = −6.79461 10−4
B11 = −2.46211820 + 2.97054714Ppr − 2.86264054x10−1 Ppr2 + 8.05420522x10−3 Ppr3
B33 = −7.93385684x10−1 + 1.39643306Ppr − 1.49144925x10−1 Ppr2 + 4.41015512x10−3 Ppr3 B44 = 8.39387178x10 − 1.86408848x10 Ppr −4
+ 2.03367881x10 P − 6.09579263x10 P
3 pr
viskositas
(24)
Setelah densitas dan viskositas kondensat dan air diperoleh, maka densitas dan viskositas fluida dapat diperoleh dari persamaan berikut:
CGR ⎛ ⎞ ⎛ WGR ⎞ ρl = ⎜ ⎟ ρo + ⎜ ⎟ ρw CGR WGR CGR WGR + + ⎝ ⎠ ⎝ ⎠
Densitas kondensat
)
Densitas dan viskositas fluida
−1
Perhitungan untuk menentukan kondensat, adalah sebagai berikut:
(
B3 = −5.47119 10−5
μw = μw1 ⎡⎣0.9994 + 4.0295(10−5 ) P + 3.1062 (10−9 ) P2 ⎤⎦
+ 3.60373020x10−1 Ppr2 − 1.044324x10−2 Ppr3
2 pr
)
)
Lalu, viskositas air pada tekanan tertentu dapat diperoleh dengan persamaan berikut:
B22 = 2.80860949 − 3.49803305Ppr
−2
(
B1 = 2.63951 10−2
B0 = −1.12166
= B11 + B22Tpr + B33Tpr2 + B44Tpr3
−2
= B0 + B1S + B2 S 2 + B3 S 3
CGR ⎛ ⎞ ⎛ WGR ⎞ μl = ⎜ ⎟ μo + ⎜ ⎟ μw ⎝ CGR + WGR ⎠ ⎝ CGR + WGR ⎠
___________________________________________________________________________________________________
(25) (26)
3
Ikatan Ahli Teknik Perminyakan Indonesia Prosedur iterasi Metode komputasi yang digunakan di sini adalah untuk jaringan pipa kompleks. Kita tidak perlu melakukan iterasi untuk menentukan kecepatan massa ( Wm ). Berikut adalah prosedur iterasi untuk
ditentukan, maka proses ini telah konvergen. Namun, jika ξ lebih besar daripada τ , ulangi langkah (3) dengan menggunakan nilai-nilai baru dimana tekanan didefinisikan sebagai berikut:
Pi =
Pi k + Pi k +1
jaringan pipa kompleks: 1. Masukan nilai kecepatan massa Wm untuk tiap pipa. 2. Masukan nilai hampiran tekanan di setiap titik node. 3. Tentukan nilai hampiran untuk nilai tekanan rata-rata P pada tiap pipa dengan menggunakan persamaan dibawah ini:
Pj =
P ( i, j ) + P ( i, j + 1)
(27)
2
4. Dengan nilai hampiran P dan T , prediksi dengan nilai-nilai ρ g , ρ L , μ g , μ L dan λL mengunakan model sifat fasa. 5. Dengan menggunakan kecepatan massa Wm dan sifat fluida pada langkah (4), akan ditentukan nilai hampiran densitas campuran ρ m , dan kecepatannya, vm . 6. Tentukan nilai hampiran liquid holdup, HL, untuk setiap konfigurasi horizontal pada tiap pipa, dengan menggunakan persamaan korelasi empiric Beggs-Brill. 7. Tentukan nilai hampiran dari kecepatan superfisial gas, vsg. 8. Tentukan nilai hampiran bilangan Reynolds NRe . 9. Tentukan nilai faktor friksi campuran, fm, dengan menggunakan persamaan Chen 10. Hitung nilai konstanta C1, C2, dan C3 yang telah didefinisikan pada persamaan (3), (4), dan (5), 11. Untuk setiap pipa, hitung nilai koefisien faktor, η , yang telah didefinisikan pada persamaan (2). 12. Susun matriks persamaan untuk menentukan nilai ΔP dari persamaan (6) dengan menggunakan eliminasi Gauss. 13. Tentukan nilai tekanan pada setiap titik node pada sistem jaringan. 14. Hitung nilai Euclidean dari tekanan yang telah dihitung, yaitu
ξ=
∑(P N
i
i =1
m
− Pi m +1
)
2
(28)
15. Periksa apakah tekanan pada iterasi ini sudah konvergen. Jika nilai ξ dari langkah (14) lebih kecil dari nilai toleransi τ yang telah
(29)
2
Studi Kasus Model yang digunakan dalam studi kasus ini adalah jaringan pipa multifasa kompleks di lapangan X (Gambar 1) .Tujuan studi kasus ini adalah untuk memprediksi tekanan pada setiap well head dan stasiun pengumpulan. Dimulai dari nilai tekanan tetap di NMP, perhitungan distribusi tekanan dilakukan seperti langkah-langkah di atas. Data masukan meliputi data dari well head dan geometri pipa tiap segmen
Analisis dan Diskusi Hasil simulasi dengan menggunakan Software OPPINET yang dibandingkan dengan software komersial menunjukkan hasil seperti pada tabel 1.Dari hasil yang diperoleh software OPPINET dan software komersial, terdapat perbedaan rata – rata dalam prediksi tekanan sebesar 2.75 % dari Pipephase dan 2 % dari Pipesim. Software OPPINET mampu memprediksikan pola aliran pada masing – masing segmen (tabel 2) yang dapat dikembangkan lebih lanjut untuk memprediksikan hold up pada masing – masing segmen pipa.
Kesimpulan Berdasarkan penjabaran dalam makalah ini, dapat disimpulkan bahwa: 1. Model dua fasa telah dikembangkan untuk memprediksi distribusi tekanan dalam jaringan pipa aliran multifasa kompleks. 2. Model ini terdiri dari 3 konsep yaitu: model dua-fasa, model sifat fasa, dan prosedur iterasi. 3. Perbedaan prediksi tekanan antara software OPPINET dan software komersial rata – rata sebesar 2.75% (Pipephase) dan 2 % (Pipesim) 4. Perbedaan antara Software OPPINET dan software komersil lainnya, disebabkan oleh perbedaan metode perhitungan.
___________________________________________________________________________________________________
4
Ikatan Ahli Teknik Perminyakan Indonesia Nomenclature
Reference
ft 2
Ap
:
Cross sectional area of pipe,
D d
:
Demand,
:
Inside pipe diameter, ft
f
:
Pipe friction factor
g
:
Gravity acceleration, ft sec
gc
:
Conversion factor, lbm . ft lbf . sec
HL k L M N N FE N RE R
:
Liquid holdup, fraction
: : : :
Time level in the iteration calculation Total number of legs Total number of loop in a network Total number of nodes
:
Froude number
:
Reynolds number
:
lbm / D
2 2
P Qg
:
Right hand side of the system of equations Pressure, psia
:
Volumetric gas flow rate, cuft / D
S T v sg
: : :
Supply, lbm / D Temperature, ° R Superficial gas velocity,
Wm Z
:
Mass flow rate of mixture
:
Length of distance, ft
:
Pipe roughness, ft
: : :
Volume fraction of liquid, fraction Viscosity, cp Two-phase flow equation coefficient
: : :
Density, lbm ft Tolerance Euclidean
ε λL μ η ρ τ ξ
1. Ahmed, T. 1989. Hydrocarbon Phase Behavior. Houston, Texas: Gulf Publishing Company. 2. Mucharam, L. dan Adewumi M.A. 1990. A Compositional Two Phase Flow Model for Analyzing and Designing Complex Pipeline Network System. Pennsylvania: Pennsylvania State U. 3. Mucharam, L. et al. 2007. Two-Phase Flow Model for Predicting Pressure Distribution in Complex Pipeline Network. RC-OPPINET 6th Annual Report. 4. McCain, W.D. 1990. The Properties of Petroleum Fluids. Tulsa, Oklahoma: PennWell Publishing Company.
ft hr
3
Subscript
g
L m T
Relating to the gas phase Relating to the liquid phase Relating to the mixture Relating to the Total
Superscript
k →
Relating to the time level Vector Average
___________________________________________________________________________________________________
5
Ikatan Ahli Teknik Perminyakan Indonesia
Tabel 1. Hasil Pengolahan Data Node 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
Pipephase (psia) 293.6 293.5 293.7 293.7 293.8 293.6 293.5 293.7 293.7 293.5 293.7 293.5 293.5 290.1 290 290.1 290 290 290.5 290 290.3 290 287.9 287.9 287.9 287.9 287.9 287.9 287.9 287.9 287.9 287.9 287.8 287.8 287.8 287.8 287.6 287.6 286.7
Pipesim (psia) 286.76 286.754 286.814 286.813 286.87 286.794 286.75 286.84 286.845 286.755 286.855 286.755 286.748 288.41 288.34 288.39 288.364 288.35 288.57 288.33 288.5 288.328 326.233 326.23 326.23 326.31 326.235 326.23 326.24 326.24 326.231 287.414 287.36 287.36 287.346 287.344 287.249 287.714 286.7
OPPINET (psia) 286.95 287.00 286.99 287.04 287.04 287.10 287.52 287.21 287.25 287.26 287.30 287.32 286.84 290.99 290.99 291.09 291.15 291.19 291.24 291.27 291.29 290.89 306.66 306.71 306.74 308.55 306.82 306.87 307.18 306.98 306.57 295.41 295.45 295.51 295.56 295.27 294.43 286.78 286.70
Differences with Pipephase (%) 2.27 2.22 2.28 2.27 2.30 2.22 2.04 2.21 2.20 2.13 2.18 2.11 2.27 0.31 0.34 0.34 0.40 0.41 0.25 0.44 0.34 0.31 6.52 6.53 6.54 7.17 6.57 6.59 6.70 6.63 6.49 2.61 2.66 2.68 2.70 2.59 2.37 0.28 0.00
Differences with Pipesim (%) 0.07 0.08 0.06 0.08 0.06 0.11 0.27 0.13 0.14 0.18 0.16 0.20 0.03 0.89 0.92 0.94 0.96 0.99 0.93 1.02 0.97 0.89 6.00 5.98 5.97 5.44 5.95 5.93 5.84 5.90 6.03 2.78 2.82 2.83 2.86 2.76 2.50 0.32 0.00
___________________________________________________________________________________________________
6
Ikatan Ahli Teknik Perminyakan Indonesia
Tabel 2. Hasil Pengolahan Data OPPINET Node 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
(psia)
286.95 287.00 286.99 287.04 287.04 287.10 287.52 287.21 287.25 287.26 287.30 287.32 286.84 290.99 290.99 291.09 291.15 291.19 291.24 291.27
Pola
Node
Aliran
segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated
OPPINET (psia)
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
291.29 290.89 306.66 306.71 306.74 308.55 306.82 306.87 307.18 306.98 306.57 295.41 295.45 295.51 295.56 295.27 294.43 286.78 286.70
Pola Aliran
segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated segregated
___________________________________________________________________________________________________
7
Ikatan Ahli Teknik Perminyakan Indonesia
4 5
6
7
Gathering Station
3
Well head 2
8
NMP
1
9 13
10
286,7 Psia 38
11
14
39
37 12
15 22
16
36 23
31 34
24
17
32 18
35
33
25 19 20
21
26
30 27
28
29
Gambar 1. Model Jaringan Pipa Kompleks
Gambar 2. Hasil Perbandingan Distribusi Tekanan Software OPPINET dan Software Komersial ___________________________________________________________________________________________________
8
Ikatan Ahli Teknik Perminyakan Indonesia
___________________________________________________________________________________________________
9