Pemodelan Temperatur Keluaran Sistem Downhole Heat Exchanger dengan Metoda Elemen Hingga

Jurnal Matematika & Sains, Agustus 2012, Vol. 17 Nomor 2

Pemodelan Temperatur Keluaran Sistem Downhole Heat Exchanger dengan Metoda Elemen Hingga Alamta Singarimbun dan Gilang Satria Prayoga KK Fisika Sistem Kompleks, Program Studi Fisika Fakultas Matematika dan Ilmu Pengetahuan Alam, Institut Teknologi Bandung e-mail: [email protected] Diterima 3 Oktober 2011, disetujui untuk dipublikasikan 2 Maret 2012 Abstrak Downhole Heat Exchanger (DHE) merupakan salah satu teknologi untuk memanfaatkan energi panas bumi. Secara umum DHE terdiri dari pipa berbentuk U yang ditanam secara vertikal di dalam tanah. Ke dalam pipa DHE dialirkan air dan jika temperatur di bawah tanah cukup tinggi, maka air tersebut akan terpanaskan. Air yang telah terpanaskan bertemperatur tinggi dapat dimanfaatkan misalnya untuk kebutuhan rumah tangga. Dalam penelitian ini, kami mengembangkan sebuah model dalam bentuk simulasi numerik untuk memperkirakan temperatur sebagai keluaran DHE berdasarkan parameter fisik optimal yang mempengaruhi sistem. Sebagai hasilnya, diperoleh kaitan antara debit air yang dimasukkan dengan temperatur keluaran untuk temperatur masukan tertentu. Juga diperoleh hasil berupa kaitan antara temperatur masukan dengan temperatur keluaran untuk debit air tertentu. Parameter ini diharapkan dapat dimanfaatkan sebagai kontrol untuk parameter keluaran yang diinginkan. Kata kunci : DHE, Pipa berbentuk U, Simulasi numerik, Output termal, Debit air, Temperatur keluaran, Temperatur masukan.

A Model for Temperature Estimation as Output of Downhole Heat Exchanger by Using Finite Element Method Abstract Downhole Heat Exchanger (DHE) is one of the technologies for harnessing geothermal energy. DHE generally consists of a U-shaped pipe planted in the ground vertically. Water is flowed into the DHE pipe and if the temperature below ground is high enough, then the water will be heated. Water that has been heated by high temperature can be used for example for household needs. In this research, we develop a model in the form of numerical simulations to estimate temperature as output of DHE based on the optimal physical parameters that affect the system. As a result, the relationship between water debit discharge and output is obtained as function of temperature. The relationship between temperature insert with the water discharge temperature are obtained also for a specific output. The parameter input can be used as a control for the desired output parameter. Keywords: DHE, U-shaped, Numerical simulation, Thermal output, Water discharge, Temperature output, Input temperature. disirkulasikan didalamnya. Air menyerap panas dari reservoir (Lienau, 1984). Pada sistem DHE tidak ada fluida atau substansi lain dari reservoir yang digunakan secara langsung di permukaan, hanya kalor atau energi panasnya yang diambil. Tahap perpindahan panas pada sistem DHE adalah sebagai berikut: mula-mula air dari permukaan dialirkan melalui pipa ke bawah permukaan bumi. Air yang mengalir di dalam pipa akan menyerap panas reservoir. Air di dalam pipa akan mengalir ke luar melalui pipa keluar. Dalam proses ini temperatur air akan mengalami penurunan karena bersinggungan dengan dinding-dinding pipa yang temperaturnya lebih rendah. Temperatur yang keluar air dipengaruhi oleh laju atau debit aliran fluida Q (m3 sec-1) dan temperatur masukan fluida.

1. Pendahuluan Downhole Heat Exchanger (DHE) merupakan suatu sistem penukar panas dengan menggunakan energi panas bumi. Sistem DHE dapat diaplikasikan dalam pengembangan energi panas bumi sektor nonlistrik, misalnya untuk beberapa keperluan praktis sehari-hari dalam pengolahan hasil pertanian dan peternakan. Untuk beberapa negara tertentu, mereka menggunakannya sebagai sumber air panas yang dialirkan ke rumah-rumah penduduk. DHE biasa digunakan sebagai sumber panas dalam penghangat ruangan seperti rumah, kantor, sekolah, atau ruangan yang terletak di area dimana di bawahnya terdapat reservoir dengan kedalaman rendah dan bertemperatur sedang yaitu antara 60 hingga 150 oC (Tago, 2006). Sistem DHE terdiri dari pipa berbentuk "U" yang ditanam secara vertikal (lihat Gambar 1). Ke dalam pipa dialirkan fluida, biasanya air, yang 55

56

Jurnal Matematika & Sains, Agustus 2012, Vol. 17 Nomor 2

u 2  x 2

u1

u 2  ... x 2

y

u1  x1

u1  ... x1

x

u 2

Gambar 2. Blok kontrol volume fluida. Gambar 1. Konfigurasi pipa U sistem DHE. Perpindahan panas merupakan dasar dalam sistem DHE. Perpindahan panas adalah perpindahan energi antar material yang disebabkan oleh perbedaan temperatur (Bejan, 2003). Perpindahan panas dibagi menjadi tiga jenis, yakni konduksi, konveksi, dan radiasi. Pada sistem penukar panas DHE, hanya konduksi dan konveksi yang berperan. Konduksi merupakan proses perpindahan energi antar molekul, tanpa adanya pergerakan aktual dari molekul tersebut. Perpindahan panas secara konduksi bergantung kepada karakteristik medium penghantar. Dalam konveksi, perpindahan energi yang diakibatkan pergerakan makroskopik molekul, biasanya terjadi pada zat cair dan gas. Untuk menganalisis secara kuantitatif besaran energi dalam proses perpindahan panas, diperlukan persamaan laju perubahan. Untuk proses konduksi, laju perubahan energi dapat dinyatakan melalui Hukum Fourier (Gregory, 2009)  q  kT (1) dimana q adalah laju fluks kalor, k adalah konduktivitas kalor dan T adalah temperatur. Proses konveksi laju perubahan dijelaskan melalui hukum pendinginan Newton yang dituliskan sebagai : q  h(Tw  Ta )

(2)

dengan h adalah koefisien konveksi dan Tw – Ta adalah perubahan temperatur fluida dengan zat padat. Model matematis yang menggambarkan dinamika fluida kontinu adalah persamaan Navier-Stokes yang terdiri dari tiga persamaan dasar: persamaan kontinuitas, persamaan momentum, dan persamaan energi. Tinjaulah kontrol volume pada Gambar 2. Laju fluks massa yang masuk ke kontrol volume pada arah x1 adalah ρu1, dan pada arah x2 adalah ρu2. Diasumsikan tidak ada reaksi atau produksi massa di dalam kontrol volume.

Dengan menggunakan deret Taylor dapat ditentukan perubahan fluks massa di dalam kontrol volume untuk masing-masing arah kecepatan i=1,2 : xi [( ui ) xi  ( ui ) xi  xi ]    ( u i ) x12  2 ( ui )    xi xi  ... 2 xi 2 xi  

(3)

Laju perubahan massa di dalam kontrol volume didefinisikan sebagai  x1x2 . Dengan menerapkan t

hukum kekekalan massa, maka jumlah fluks massa di dalam kontrol volume sama dengan laju perubahan massa di dalam kontrol volume, sehingga didapat persamaan kontinuitas Navier-Stokes :   ( ui )  0 t xi

(4)

dimana ρ adalah kerapatan massa fluida dan ta adalah waktu. Selanjutnya dapat diturunkan perubahan momentum setiap waktu, yaitu :  ( ui )  ( ui2  ( ui u j ) p  ij (5)     t xi xi xi x j

Dari Persamaan (5)  i dapat dituliskan sebagai fungsi dari viskositas dinamik dan kecepatan sebagai berikut:  ui ui 2 ui   ij     xi x j 3 x j   

 ij   

(6)

Laju perubahan energi persatuan waktu di dalam kontrol volume didefinisikan sebagai : E  Q in  Q out  Q gen (7) dengan dT E  c dV dt Q  Q  QdV

(8b)

E gen  Qgen dV

(8c)

in

out

(8a)

Dalam Persamaan (7), (8a), (8b), dan (8c) E adalah netto perubahan energi, sedangkan subskrib in, out dan gen menunjukkan perubahan energi yang masuk, ke luar dan yang terakumulasi dalam blok yang

Singarimbun dan Prayoga, Pemodelan Temperatur Keluaran Sistem Downhole Heat Exchanger ................. 57 ditinjau. Dengan mensubstitusikan hukum konduksi Fourier, maka didapat persamaan energi untuk kontrol volume 2 dimensi sebagai :

c

  2T  2T  dT  k  2  2   Qgen dt  x1 x21 

dimana c adalah kapasitas konduktifitas panas fluida.

panas,

(9)

k

adalah

2. Metodologi dan Pemodelan Dalam penelitian ini digunakan Elemen Hingga (Finite Element) untuk memecahkan persamaan transfer energi dan aliran fluida. Metode ini digunakan untuk menyelesaikan persamaan diferensial parsial dengan membagi setiap sistem menjadi elemen-elemen dengan geometri tertentu. Dari elemen-elemen tersebut akan tersusun matriks persamaan linear. Dalam proses penyelesaian elemen hingga diperlukan elemen dengan geometri tertentu beserta fungsi bentuk atau shape function, yang merepresentasikan ciri solusi dari setiap elemen. Dalam penelitian ini digunakan elemen garis (1-D) dan elemen segitiga (2-D) dengan fungsi bentuk linear.

Gambar 3. Elemen segitiga. Dengan mensubstitusikan  1 ,  2 , dan  3 , ajan diperoleh fungsi bentuk untuk elemen segitiga.



T  N iTi  N j T j  N k Tk  N i

T ( x, y )   1   2 x   3 y

(10)

Pada Gambar 3, nilai T pada masing-masing titik adalah:

 Ti  N k T j  (13) Tk 



dimana

2.1 Elemen segitiga linear Pada elemen segitiga, interpolasi linear untuk aproksimasi fungsi tertentu dapat dinyatakan sebagai berikut (Zienkiewicz dkk., 2005):

Nj

Ni 

1 ai  bi x  ci x  2A

Nj 

1 a j  bjx  cjx 2A

Nk 

1 a k  bk x  c k x  2A



(14a)



(14b) (14c)

dan ai  x j y k  x k y j ; bi  y j  y k ; ci  x k  x j a j  x k y i  xi y k ; b j  y k  y i ; c j  xi  x k a k  x i y j  x j y i ; bk  y i  y j ; c k  x j  xi

Ti   1   2 xi   3 y i

(11a)

2.2 Teorema residual dan formulasi Galerkin

T j  1   2 x j   3 y j

(11b)

Tk   1   2 x k   3 y k

(11c)

Metode residual adalah metode aproksimasi untuk menyelesaikan persamaan diferensial dengan memberikan nilai bobot pada residual atau sisa dari persamaan aproksimasi (Zienkiewicz dkk., 2005). Formulasi Galerkin digunakan agar residual menjadi minimal dengan mengalikan integrasi residual dengan suatu fungsi bobot wi(x);

yang menghasilkan koefisien-koefisien

1 





1 x j y k  xk y j Ti  2A  xk yi  xk yi T j  x j yi  xi y j Tk



1 

 w ( x)R( x)dx  0

(12b)



  

(12c)

dengan A adalah luas elemen segitiga.

(15)

i

  

1 xk  x j Ti  xi  xk T j  2A  x j  xi Tk



(12a)



1 1  y j  y k Ti   y k  yi T j  2A  yi  y j Tk



 



dimana fungsi bobot diganti dengan fungsi bantu atau shape function dari metode elemen hingga yang digunakan sebagai aproksimasi. Dengan meninjau persamaan energi Navier-Stokes pada keadaan tunak serta menerapkan metode residual dan formulasi Galerkin dan syarat batas Neumann maka diperoleh:  N k N m N k N m   k x  ky x x y y 



 kz

N k N m   z z 

Tm d   GN k d   N k qdS  

S

58

Jurnal Matematika & Sains, Agustus 2012, Vol. 17 Nomor 2

  hN k n m {Tm }ds   hT0 ds  0 S

S

Dengan menyederhanakan Persamaan berdasarkan koefisien Tm didapatkan :  N k N m N k N m [ K ]     k x  ky x x y y   kz

(16) (16)

(17)

N k N m  d   hN k N m dS z z  S

 f    GN k d   N k qdS   hTo dS 

S

(18)

S

2.3 Metode karakteristik base split Skema karakteristik Galerkin atau biasa disebut CBS (Characteristics Base Split) merupakan skema pemisahan bagian tekanan p dari persamaan momentum pada persamaan Navier-Stokes. Selanjutnya Persamaan (17) dan (18) didiskritisasi dengan menggunakan deret Taylor, sehingga didapatkan bentuk semi-diskrit untuk persamaan momentum u’ dan v’ yang merupakan komponenkomponen kecepatan pertengahan. Setelah itu dilakukan perhitungan tekanan, dari perhitungan ini maka komponen kecepatan un+1, dan vn+1 didapatkan dengan melibatkan komponen tekanan yang telah dihitung pada tahap sebelumnya.

dua langkah pemodelan, yaitu pemodelan temperatur latar atau sumber dan pemodelan temperatur air sistem DHE. Pemodelan temperatur latar dilakukan karena tidak ada data temperatur lapangan sesungguhnya yang diperlukan sebagai sumber energi panas dari sistem DHE. Pemodelan keseluruhan dilakukan secara terpisah, yakni pemodelan temperatur sumber dilakukan lebih dahulu kemudian dilakukan pemodelan temperatur sistem fluida DHE. Pemisahan ini dilakukan dengan asumsi bahwa waktu “operasional” DHE dilakukan setelah distribusi temperatur sumber pada sumur relatif stabil. 2.6.1 Penentuan distribusi temperatur sumber Untuk menentukan distribusi temperatur sumber digunakan mesh seperti tertera pada Gambar 4. Fluida yang terdapat pada daerah sumber diasumsikan bersifat tak termampatkan dan memiliki viskositas konstan.

2.4 Diskritisasi temporal Untuk kasus aliran dua dimensi, persamaan Navier-Stokes didiskritisasi menurut langkah-langkah berikut (Zienkiewicz dkk., 2005) : langkah 1: intermediasi momentum u’ dan v’, langkah 2: perhitungan tekanan, langkah 3: koreksi kecepatan untuk un+1 dan vn+1, langkah 4: perhitungan temperatur. 2.5 Diskritisasi spasial Sebelum proses diskritisasi dilakukan, persamaan-persamaan semi diskrit Navier-Stokes diubah terlebih dahulu ke dalam bentuk nondimensional, dan diasumsikan fluida bersifat inkompresibel (Zienkiewicz dkk., 2005) tu x y ; y ;t a L L L u v p ;v ;p u ua ua u a2

x

T

T  Ta Tw  Ta

Dari persamaan-persamaan diskiritisasi temporal didapat aproksimasi Galerkin dalam bentuk momentum, tekanan, koreksi kecepatan dan temperatur. 2.6 Langkah pemodelan Untuk mengetahui distribusi temperatur air dan temperatur keluaran dari sistem DHE, dilakukan

Gambar 4. Mesh dan syarat batas pemodelan (Massarotti, 2004).

Sumber panas dimodelkan sebagai sebuah sumur persegi empat, dimana terdapat lubang masukan fluida dari reservoir atau sistem hidrotermal yang terletak di sekitar sumber. Syarat batas untuk pemodelan distribusi temperatur sumber secara lebih jelas dapat dilihat pada Gambar 4. Pada lubang input temperatur dan laju aliran para arah horizontal diset selalu konstan, begitu juga kecepatan horizontal pada lubang output. Pada dinding-dinding sumur diasumsikan kondisi tidak slip, begitu juga pada daerah kontak fluida pada sumur dengan pipa. Pada permukaan atas sumur, temperatur dan kedua komponen kecepatan diatur sama dengan nol. Nilainilai tekanan diatur konstan (p=1) pada titik-titik ujung pipa dan titik-titik pada input dan output. Langkah pemodelan yang dilakukan menggunakan skema CBS sehingga digunakan besaran-besaran fisis non-dimensional. 2.6.2 Penentuan distribusi temperatur fluida DHE Temperatur air yang mengalir pada pipa bergantung pada jumlah fluks kalor yang masuk dan mengalir dari pipa ke air, demikian juga dengan temperatur pipa akan mengalami perubahan temperatur ketika air mengalir pada dinding pipa. Dengan demikian dapat dikatakan perpindahan panas antara pipa dengan air yang mengalir di dalamnya merupakan sistem yang terkopel. Untuk memecahkan

Singarimbun dan Prayoga, Pemodelan Temperatur Keluaran Sistem Downhole Heat Exchanger ................. 59

dT (19) dx Ruas kanan persamaan di atas dapat disubstitusikan menggunakan persamaan pendinginan Newton:

memberikan perbedaan temperatur antara masukan dengan keluaran yang semakin tinggi.

Kedalaman (m)

sistem yang terkopel seperti ini diasumsikan temperatur awal pipa sama dengan temperatur perbatasan antara pipa dengan sumber yang telah ditentukan dari proses sebelumnya. Proses perpindahan panas diasumsikan hanya terjadi antara dinding bagian dalam pipa dengan air yang mengalir didalamnya. Perpindahan panas antara bagian luar pipa dengan fluida sumber diabaikan. Untuk menentukan besar perubahan energi antara dinding pipa dengan air, diberikan syarat batas Neumann pada dinding bagian dalam pipa: q  k

q  h(T f  TP )

(20)

dimana Tf adalah temperatur fluida dan Tp adalah temperatur pipa. 2.6.3 Penentuan selang waktu dan parameter kompresibilitas Selang waktu untuk menentukan seluruh parameter yang dapat ditentukan dari skema CBS adalah selang waktu lokal, yakni selang waktu minimum antara selang waktu perubahan kecepatan air dan selang waktu perubahan temperatur (Massarotti, 2004), yang didefinisikan sebagai

Distribusi temperatur dengan variasi temperatur masukan dapat dilihat pada Gambar 6. Temperatur masukan yang besar memberikan temperatur yang besar. Semakin besar temperatur masukan maka perbedaan antara temperatur masukan dengan temperatur keluaran semakin besar. Dengan kata lain, sistem DHE yang menggunakan model sumur seperti pada Gambar 6 menghasilkan temperatur output yang besar jika temperatur inputnya diperbesar.

h 2 Re h , t d  u 2

dengan h adalah panjang minimum elemen. Parameter kompresibilitas β ditentukan menggunakan hubungan (Nithiarasu, 2003) :

  max( , vd , ve )

(21)

Dengan ε adalah konstanta (0.1 – 0.5), dan vc  u  v vd 

Gambar 5. Distribusi temperatur air pada pipa U dengan variasi laju aliran.

Kedalaman (m)

t  min(t c , t d ) , t c 

Temperatur (oC)

(22)

2 hRe

(23) Temperatur (oC)

3. Hasil dan Diskusi

Distribusi temperatur sumber ditentukan dengan menggunakan algoritma CBS. Dalam hal ini diandaikan temperatur pada lubang masukan adalah 120o C dan temperatur permukaan adalah 20o C dengan parameter fisis tertera pada Tabel 1. Tabel 1. Konstanta Termal Air (Tago, 2006).

Μ (Darcy)

ρ (kg/m3)

Cp (J/KgK)

k (W/mK)

8 × 10-4

999.8

4,207.5

0.58

Distribusi temperatur pipa selama satu jam ditunjukkan pada Gambar 5. Distribusi air di dalam pipa meningkat terus seiring bertambahnya kedalaman pipa, mengikuti distribusi temperatur pipa yang juga meningkat dengan bertambahnya kedalaman. Temperatur air di dalam pipa bergantung dengan laju alirannya. Laju aliran yang lebih lambat

Gambar 6. Distribusi temperatur air pada pipa U dengan variasi temperatur masukan.

Sebagai perbandingan, dilakukan simulasi distribusi temperatur dengan model sumur yang memiliki temperatur linear terhadap kedalaman, sehingga temperatur pipa juga linear, tidak bergantung pada dekat atau tidaknya pipa dengan sumber panas. Hasilnya diperoleh berupa distribusi temperatur yang memiliki karakteristik sama dengan model sumur sebelumnya. Perbedaan temperatur keluaran dengan temperatur masukan hampir sama untuk seluruh temperatur masukan. 6. Kesimpulan

Dari hasil pemodelan dan simulasi temperatur latar dan temperatur keluaran pipa dari sistem DHE dapat disimpulkan bahwa distribusi temperatur sistem

60

Jurnal Matematika & Sains, Agustus 2012, Vol. 17 Nomor 2

DHE bergantung kepada laju aliran air. Semakin kecil lajunya maka distribusi temperatur air pada pipa DHE akan semakin tinggi. Hasil pemodelan menunjukkan bahwa dengan temperatur masukan 20 o C; untuk debit 4 m3/s menghasilkan temperatur keluaran 29,5 oC, debit 8 m3/s memberi temperatur keluaran 24,5 oC, debit 10 m3/s memberi temperatur keluaran 23,5 oC dan debit 12 m3/s memberikan temperatur keluaran 23 oC. Diperoleh juga untuk debit air 5 m3/s temperatur keluaran hasil simulasi dalam sistem DHE untuk temperatur masukan 15 oC memberi temperatur keluaran 22 oC, temperatur masukan 20 oC memberi temperatur keluaran 27,5 oC dan temperatur masukan 25 oC memberi temperatur keluaran 33,5 oC. Dari model yang dibuat terlihat perbedaan temperatur antara input dan output besar ketika input diperbesar. Perubahan temperatur yang signifikan terjadi jika temperatur air masukan diperkecil, dengan kata lain efisiensi sistem DHE akan semakin besar jika input diperkecil. Hasil-hasil ini diharapkan dapat dipakai sebagai kontrol untuk mengatur parameter masukan agar diperoleh parameter keluaran yang optimal sesuai dengan temperatur keluaran yang diinginkan. Sebagai aplikasi, sistem DHE dapat diaplikasikan untuk pengembangan energi panas bumi sektor nonlistrik, misalnya untuk mengolah hasil pertanian dan peternakan. Beberapa negara tertentu menggunakannya sebagai sumber air panas yang dialirkan ke rumah-rumah.

Daftar Pustaka

Bejan, A. and A. D. Kraus, 2003, Heat Transfer Handbook, John Wiley & Sons, New Jersey. Gregory, N. S. K., 2009, Heat Transfer, Cambridge University Press. Lewis, R. W., P. Nithiarasu, and K. N. Seetharamu, 2004, Fundamentals of the Finite Element Method for Heat and Fluid Flow, John Wiley and Sons, Chichester, West Sussex, England. Lienau, P. J., 1984, Geothermal district heating institutional factors—the Klamath falls experience, Geotherm Resources Councill Bull, 6-11. Massarotti, N. F. A., 2004, Fully Explicit and SemiImplicit CBS Procedures For Incompressible Flows, European Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS 2004, 6. Nithiarasu, P., 2003, An Efficient Artifcial Compressibility (AC) Scheme Based on Split (CBS) Method For Incompressible Flows, Int. J. Numer. Meth. Eng., 56:13, 1815–1845. Tago, M. K. M., 2006, Heat Extraction Characteristics of a Single U-Tube Downhole Heat Exchanger With Square Cross Section, Heat Mass Transfer, 42, 608616. Zienkiewicz, R. T., R. L. Taylor, and P. Nithiarasu, 2005, The Finite Element Method for Fluid Dynamics, 6th ed., Elsevier ButterworhHeinemann, Jordan Hill, Oxford.