## Perpindahan Panas pada Media Berpori Menggunakan  Metode Elemen Hingga
Surya Ningsih *1 , Wahyu Srigutomo 2
1,2 KK Fisika Bumi dan Sistem Kompleks Program Studi Magister Fisika FMIPA Institut Teknologi Bandung e-mail: suryaningsih885@yahoo.co.id, wahyu@fi.itb.ac.id

## Abstrak
Masuknya  aliran  panas  dari  bawah  reservoir  ke  dalam  medium  berpori  di  dalamnya mempengaruhi  kondisi  tekanan,  entalpi,  dan  temperatur  di  dalam  reservoir.  Hal  ini terjadi karena panas mengalir dari daerah bertemperatur tinggi ke daerah bertemperatur lebih  rendah.  Panas  mengalir  melalui  daerah  patahan  (fracture  zone)  yang  memiliki permeabilitas dan porositas lebih tinggi dari daerah sekitarnya di dalam reservoir untuk kemudian terjadi perpindahan panas secara konduksi menuju daerah sekitarnya. Tulisan ini  mensimulasikan  proses  perpindahan  panas  dengan  mengamati  distribusi  temperatur dan aliran fluida pada medium berpori. Metode numerik yang digunakan adalah metode elemen  hingga.  Hasil  yang  ditampilkan  adalah  distribusi  tekanan,  kecepatan  aliran fluida, dan distribusi temperatur dan kemudian diinterpretasikan.
Kata  kunci :  permeabilitas,  porositas,  kecepatan,  distribusi  temperatur,  metode  elemen hingga.

## 1. PENDAHULUAN
Perpindahan panas pada sistem panas bumi adalah perpindahan energi panas dari daerah yang bertemperatur tinggi ke daerah yang bertemperatur lebih rendah. Perpindahan panas dapat terjadi dalam tiga mekanisme, yaitu: konduksi, konveksi, dan radiasi. Pada sistem panas  bumi,  hanya    konduksi  dan  konveksi  yang  berperan.  Perpindahan  panas  secara konduksi  adalah  perpindahan  energi  karena  adanya  getaran  molekul  pada  batuan. Perpindahan  panas  secara  konveksi  adalah  perpindahan  energi  melalui  pergerakan  dari fluida.
Perpindahan    panas  yang  terjadi  di  dalam  bumi  merupakan  persoalan  kompleks karena  melibatkan  banyak  parameter  dan  bersifat  inhomogen  dan  anisotropi.  Sehingga penyelesaian  persoalan  perpindahan  panas  dalam  bumi  memerlukan  asumsi  -asumsi untuk  menyederhanakan  permasalahan.  Untuk  masalah  yang  kompleks  dapat  dikerjakan dengan menggunakan metode numerik. Salah satu metode numerik yang dapat digunakan adalah Metode Elemen Hingga ( Finite Elemen Method ). Metode elemen hingga adalah salah  satu  metode  numerik  yang  dapat  digunakan  untuk  memecahkan  persoalan  syarat batas  dalam  bentuk  persamaan  differensial  dengan  mencari  solusi  pada  domain-domain kecil dengan menerapkan syarat batas.
Ide  utama  dalam  metode  elemen  hingga  adalah  membagi  geomerti  domain  dari permasalahan  menjadi  sub-sub  domain  yang  lebih  kecil  atau  sering  disebut  sebagai elemen  hingga.  Bentuk  elemen  yang  sering  digunakan  adalah  elemen  segitiga  dan quadrilateral. Pada tulisan ini Penulis melakukan perhitungan solusi permasalahan syarat
batas dengan menggunakan metode elemen hingga pendekatan Galerkin. Metode elemen hingga  diterapkan  pada  persamaan  perpindahan  panas  konduksi-konveksi  dengan  solusi berupa distribusi temperatur pada domain dua dimensi.

## 2. METODE PENELITIAN

## 2.1 Metode Elemen Hingga
Dalam  penelitian  ini  digunakan  elemen  hingga  ( Finite  Elemen )  untuk memecahkan  persamaan  persamaan  perpindahan  panas  dan  kecepatan  aliran  fluida. Metode  ini  digunakan  untuk  menyelesaikan  persamaan  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  mempresentasikan  ciri  solusi  dari  setiap  elemen.  Dalam  penelitian  ini  digunakan elemen segitiga (2-D) dengan fungsi bentuk linear.

## 2.1.1 Elemen Segitiga Linear
Zienkiewick  dkk  (2005)  menjelaskan  bahwa  pada  elemen  segitiga,  interpolasi linear untuk aproksimasi fungsi tertentu dapat dinyatakan sebagai berikut:
y x y x u 3 2 1 ) , (      
(1)
Pada gambar 1, nilai T pada masing-masing  titik adalah:
k k k j j j i i i y x T y x T y x T 3 2 1 3 2 1 3 2 1                   (2)
Yang menghasilkan koefisien-koefisien
                        k u i x j x j u k x i x i u j x k x A k u j y i y j u i y k y i u k y j y A k u i y j x j y i x j u k y i x i y k x i u j y k x k y j x A                   2 1 3 2 1 2 2 1 1    (3) dengan A adalah luas elemen segitiga.
Gambar 1 Elemen Segitiga
Dengan  mensubstitusikan 1  , , 2  dan 3  akan  diperoleh  fungsi  bentuk  untuk elemen segitiga.
                k j i k j i k k j j i i T T T N N N T N T N T N T (4) dengan       y c x b a A N y c x b a A N y c x b a A N k k k k j j j j i i i i          2 1 2 1 2 1 (5) dan i j j i k k i i k j j k k j i y x y x a y x y x a y x y x a       j i k i k j k j i y y b y y b y y b       i j k k i j j k i x x c x x c x x c      

## 2.1.2 Teorema residual dan formulasi Galerkin
Metode  residual  adalah  metode  aproksimasi  untuk  menyelesaikan  persamaan differensial  dengan  memberikan  nilai  bobot  pada  residual  atau  sisa  dari  persamaan aproksimasi dalam Zienkiewicz dkk  (2005). Formulasi Galerkin digunakan agar residual menjadi minimal dengan mengalikan integrasi  residual dengan suatu fungsi bobot w i (x) ;
   0 ) ( ) ( dx x R x w i (6)
Dimana fungsi bobot diganti dengan fungsi bantu atau shape function dari metode elemen hingga  yang  digunakan  sebagai  aproksimasi.  Dengan  meninjau  persamaan  kecepatan aliran  fluida  dan  persamaan  perpindahan  panas  konduksi-konveksi  pada  keadaan  tunak serta menerapkan metode residual dan syarat batas Neumann maka diperoleh:
                                         dxdy y N k x N k dxdy y N u x N u C N K e j y i x j y j x p i e ij  (7)    e Qdxdy N f i e i
(8)
       
      e dl n y T k n x T k N p y y x x i e i (9)
2.1.3 Langkah Pemodelan
Untuk mengetahui distribusi temperatur terlebih dahulu mencari kecepatan aliran fluida dengan menggunakan persamaan (4), kemudian persamaan (4) diterapkan ke dalam metode  elemen  hingga.  Solusi  dari  persamaan  (4)  kemudian  di  masukkan  ke  dalam persamaan (6) dan diterapkan kedalam metode elemen hingga.
Adapun model domain yang digunakan adalah sebagai berikut:
a. Model reservoir bidang datar
Gambar 2 Mesh dan syarat batas pemodelan non topografi
b. Model reservoir topografi lembah
in q n u    0    n u 0    n u out q n u   
Gambar 3. Mesh dan syarat batas pemodelan topografi lembah c. Model reservoir topografi gunung
0    n u in q n u    out q n u    0    n u
Gambar 4 Mesh dan syarat batas pemodelan topografi gunung
Parameter fisis yang digunakan  adalah sebagai berikut:
Tabel 1 Paramter batuan dan fluida pada reservoir panas bumi
Parameter fisis Nilai Satuan
1  1.5x10 -12 m 2
2  9.85x10 -12 m 2 3  1x10 -12
m 2  0.1102x10 3 Pa.s k 1 2.34 W/m. o C k 2 2.0 W/m. o C k 3 0.2 W/m. o C
w  998.2
Kg/m 3
C pw 4200 J/kgK C ps 2250 J/kgK

## 2.2 Tinjauan Pustaka

## 2.2.1 Sistem Panas Bumi
Panas  bumi  merupakan  energi  panas  yang  terbentuk  secara  alami  dan  tersimpan  dalam bentuk air panas atau uap panas pada kondisi geologi tertentu pada kedalaman beberapa kilometer  di  dalam  kerak  bumi.  Sistem  panas  bumi  terdiri  dari  5  bagian,  yaitu:  Sumber panas, reservoir, fluida reservoir, lapisan penudung, dan siklus hidrologi.

## 2.2.2 Sifat Fluida Satu Fasa
Pada  paper  ini  yang  ditinjau  adalah  fluida  satu  fasa  yaitu  fasa  air.  Adapun  sifat dari fluida tersebut yaitu, volume spesifiknya adalah perbandingan antara volume fluida dengan  massanya,  densitas  fasa  fluida  merupakan  perbandingan  antara  massa  dengan volume  fasa  fluida,  dan  viskositas  yang  merupakan  ukuran  kelembaman  suatu  fluida untuk mengalir.

## 2.2.3 Sifat Batuan
Media  berpori  pada  sistem  panas  bumi  merupakan  batuan  yang  memiliki  pori sebagai tempat penyimpanan fluida panas. Reservoir pada sistem panas bumi merupakan batuan  yang  permeabel  yang  memiliki  sifat  porositas  yang  tinggi.  Permeabilitas  suatu batuan merupakan ukuran kemampuan batuan untuk mengalirkan fluida. Parameter yang menyatakan  besarnya  kemampuan  suatu  batuan  untuk  menghantarkan  panas  secara konduksi  disebut  konduktivitas  panas.  Secara  matematis  untuk  satu  dimensi  dinyatakan sebagai berikut:
) / ( dz dT Q K  (10)
dengan Q adalah laju aliran panas per satuan luas dan dT/dz adalah gradien temperatur. Satuan  dari  konduktivitas  panas  adalah  W/m.K.  Panas  spesifik  batuan  adalah  suatu parameter  yang  menyatakan  banyaknya  panas  yang  diperlukan  untuk  menaikkan  suhu persatuan massa batuan tersebut sebesar 1 o C.

## 2.2.4 Aliran Fluida dalam Media Berpori
Fluida  akan  mengalir  dalam  media  berpori  apabila  media  tersebut  empunyai permeabilitas  yang  searh  dengan  tenaga  pendorong.  Persamaan  yang  menggambarkan mengenai  aliran  fluida  dalama  media  berpori  pertama  kali  dikembangkan  oleh  Henry Darcy  (1856).  Persamaan  tersebut  merupakan  persamaan  yang  menunjukkan  kecepatan
aliran  fluida  dengan  permeabilitas  batuan,  viskositas  fluida  serta  gradien  tekanan  antar jarak tempuh aliran. Berikut merupakan persamaan Darcy:
  D      g p u   
(11)
Dengan u,  ,  , p,  , dan D berturut turut adalah  kecepatan aliran fluida, viskositas dinamik fluida, permeabilitas batuan, tekanan, densitas, dan ketinggian yang bergantung pada garavitasi. Dengan mengabaikan efek dari percepatan gravitasi, dan mengasumsikan fluida tak termampatkan, maka persamaan kontinuitas dapat dituliskan sebagai berikut:
0 .   u  (12)
Dengan mensubstitusikan persamaan (2) ke persamaan (3), maka didapatkan:
0 .          p  
(13)
Jika permeabilitas dan viskositas konstan maka persamaan (4), menjadi:
0 2   p (14)

## 2.2.5 Persamaan Konduksi dan Konveksi
Untuk  mengetahui  distribusi  temepartur  pada  sistem  panas  bumi  digunakan persamaan perpindahan panas konduksi-konveksi. Berikut adalah persamaannya
G Q T u p C T k      . ) .(  
(15)
Pada kasus sistem panas bumi Q G menggambarkan sumber panas bumi. Dalam tulisan ini Q G diasumsikan  sama  dengan  nol.  Sedangkan  k  adalah  konduktivitas  panas  yang menggambarkan seberapa mudah panas dapat ditransmisikan. Kemudian ρ densitas fluida dan C p adalah kapasitas panas.

## 3. HASIL DAN PEMBAHASAN

## 3.1 Hasil Pemodelan Aliran Fluida pada Reservoir Panas Bumi
Kecepatan  aliran  fluida  ditentukan  dengan  menggunakan  persamaan  (2).
Hasilnya  berupa  nilai  kecepatan.  Dalam  hal  ini  diandaikan  pada  model  reservoir  non topografi aliran fluida yang masuk pada sisi batas kiri sebesar 10 -4 m/s dan aliran fluida yang keluar sebesar -10 -4 . Pada model ini dibuat model homogen, dua lapis, dan 3 lapis.
Dimana lapisan bawah diberi nilai permeabilitas lebih besar dari pada lapisan atas. Pada reservoir topografi lembah aliran yang masuk pada batas sisi atas dan keluar pada batas sisi  bawah.  Pada  topografi  gunung  aliran  fluida  yang  masuk  pada  batas  sisi  kiri  dan keluar  pada  sisi  kanan  dengan  nilai  q in dan  q out sama  dengan  model  reservoir  non topografi dan topografi lembah. Berikut hasil pemodelannya:
Gambar 5 Distribusi tekanan pada reservoir homogen.
Gambar 6 Arah aliran fluida pada reservoir homogen
Gambar 7 Distribusi tekanan pada reservoir 2 lapis Horizontal Distance (m) D e pt h (m ) Pressure Distribution 0 100 200 300
400 500 600 700 800 -200 -150 -100 -50 0 -4 -3 -2 -1 0 1 2 3 4 x 10 6
Horizontal Distance (m) D e pt h (m ) Velocity 0 100 200
300 400 500
600 700 800 -200 -150 -100 -50 0 1 2 3 4 5 x 10 -4
Horizontal Distance (m) D e pt h (m ) Pressure Distribution 0 100 200
300 400 500
600 700 800
-200 -150 -100 -50 0 -1 -0.5 0 0.5 1 x 10 6
Gambar 8 Arah aliran fluida pada resesvoir 2 lapis
Gambar 9 Distribusi tekanan pada reservoir 3 lapis
Gambar 10 Arah aliran fluida pada reservoir 3 lapis
Dari  gambar  diatas  dapat  dilihat  bahwa  semakin  besar  nilai  permeabilitas  yang  diberi, maka nilai tekanan akan semakin kecil dan dapat juga dilihat bahwa semakin besar nilai permeabilitas yang diberi maka kecepatan aliran fluida akan semakin cepat.
Horizontal Distance (m) D e pt h (m ) Velocity 0 100 200 300
400
500 600 700 800 -200 -150 -100 -50 0 0.5 1 1.5 2 2.5 3 x 10 -3
Horizontal Distance (m) D e pt h (m ) Pressure Distribution 0 100 200 300
400 500
600 700 800
-200 -150 -100 -50 0 -6 -4 -2 0 2 4 6 x 10 5
Horizontal Distance (m) D e pt h (m ) Velocity 0 100 200 300
400 500
600 700 800
-200 -150 -100 -50 0 0.5 1 1.5 2 2.5 3 3.5 x 10 -3
Gambar 11 Distribusi tekanan pada reservoir topografi lembah
Gambar 12 Arah aliran fluida pada reservoir topografi lembah
Gambar 13 Distribusi tekanan pada reservoir topografi gunung
Horizontal Distance (m) D e pt h (m ) Pressure Distribution 0 100 200 300
400
500 600 700 800 -80 -60 -40 -20 0 20 40 60 80 100
-8 -6 -4 -2 0 2 4 6 8 10 x 10 5 Horizontal Distance (m) D e pt h (m ) Velocity 0 100 200 300
400
500 600 700 800 -50 0 50 100 1 2 3 4 5 6 x 10 -3 Horizontal Distance (m) D e pt h (m ) Pressure Distribution 0 100 200 300
400 500
600 700 800
-120 -100 -80 -60 -40 -20 0 20 40 60 80
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x 10 6
Gambar 14 Arah aliran fluida pada reservoir topografi gunung
Pada reservoir topografi lembah dan gunung distribusinya mengikuti model topografinya dan aliran fluidanya semakin cepat pada daerah yang permeabilitasnya besar.

## 3.2 Hasil Pemodelan Perpindahan Panas Konduksi-Konveksi
Distribusi temperatur ditentukan dengan menggunakan persamaan (6). Dalam hal ini  diasumsikan  batas  sisi  bawah  sebagai  bagian  dasar  reservoir  berupa  batuan  panas dengan  temperatur  sebesar  500 o C  yang  mentransfer  panasnya  secara  konduksi  dan konveksi  ke  dalam  reservoar  yang  berisi  fluida  air,  batas  sisi  atas  diasumsikan  sebagai temperatur awal reservoir dengan temperatur 100 o C, sedangkan pada batas sisi kiri dan kanan  diasumsikan  tidak  ada  panas  yang  keluar  pada  reservoir.  Berikut  hasil pemodelannya:
Gambar 15 Distribusi temperatur pada reservoir homogen
Horizontal Distance (m) D e pt h (m ) Velocity 0 100 200 300
400
500 600 700 800 -100 -50 0 50 1 2 3 4 5 6 7 8 x 10 -3
Horizontal Distance (m) D e pt h (m ) Temperature Distribution 0 100 200
300 400 500
600 700 800 -200 -150 -100 -50 0 100 150 200 250 300 350 400 450 Horizontal Distance (m) D e pt h (m ) Temperature Distribution 0 100 200 300
400 500
600 700 800
-200 -150 -100 -50 0 100 150 200 250 300 350 400 450
Gambar 16 Distribusi temperatur pada reservoir 2 lapis
Gambar 17 Distribusi temperatur pada reservoir 3 lapis
Pada  pemodelan  perpindahan  panas  konduksi-konveksi  pada  gambar  diatas    dapat disimpulkan bahwa pada lapisan yang nilai konduktivitas besar, perubahan temperaturnya lebih kecil dibandingkan dengan lapisan yang nilai konduktivitas yang kecil.
Gambar 18 Distribusi temperatur pada topografi lembah
Pada model diatas  terjadi perpindahan panas secara konduksi yaitu antara sumber panas dan lingkungan sekitarnya, sedangkan antara sumber panas dan fluida terjadi perpindahan panas secara konveksi.
Gambar 19 Distribusi temperatur pada topografi gunung
Horizontal Distance (m) D e pt h (m ) Temperature Distribution 0 100 200 300
400 500 600 700 800 -200 -150 -100 -50 0 100 150 200 250 300 350 400 450
Horizontal Distance (m) D e pt h (m ) Temperature Distribution 0 100 200 300
400
500 600 700 800 -80 -60 -40 -20 0 20 40
-100 0 100 200 300 400 500
Horizontal Distance (m) D e pt h (m ) Temperature Distribution 0 100 200 300
400 500
600 700 800
-120 -100 -80 -60 -40 -20 0 20 40 60 80
50 100 150 200 250 300 350 400 450 500
Pada model diatas terjadi perpindahan panas secara konduksi yaitu antara sumber panas dan lingkungan sekitarnya. Aliran konveksi tidak begitu terbentuk. Ini diakibatkan karena penempatan  nilai  permeabilitas  suatu  batuan  pada  model  tidak  diatur  sedemikian  rupa sehingga tidak terbentuk seperti yang real.

## 4. KESIMPULAN
Dari  hasil  pemodelan  kecepatan  aliran  fluida  dan  perpindahan  panas  konduksi- konveksi  disimpulkan  bahwa  nilai  tekanan  bergantung  pada  nilai  permeabilitas  batuan. Semakin kecil nilai permeabiltas maka nilai tekanan semakin besar dan aliran fluida juga ditentukan  pada  nilai  permeabilitas.  Semakin  besar  nilai  permeabilitas  maka  kecepatan aliran  fluida  akan  semakin  cepat.  Distribusi  temperatur  bergantung  pada  nilai konduktivitas batuan. Semakin besar nilai konduktivitas maka distribusi temperatur akan semakin kecil.

## DAFTAR PUSTAKA
Anonim.2007. http://infoenergi.wordpress.com/2007/03/13egs-dan-masa-depan-energi- panasbumi-di-indonesia/
C.  Polycarpu.  Anastasis.  2006. An  Introduction  to  the  Finite  Element  Method  in Electromagnetism . Morgan & Claypool Publisher : USA. Gold, T.L.:” Verttical  Two-Phasa  Steap-Water  flow  in  Geothermal  Well”, J.  Pet.  Tech. (Aug.1974) 833-842.
Kumar  C  P  2009  Groundwater  Flow  Models.  Scientist  Institute  of  hydrology,  Rookie Uttaranchal.
Lewis, R. W., K. Morgan, H.R. Thomas, and K.N. seetharamu, 1996, The Finite Element Method  in  Heat  Transfer  Analysis, John  Wiley  and  sons,  Chichester,  West  Sussex,
England
Lewis,  R.  W.,  P.  Nithiarasu,  and  K.  N.  Seetharamu,  2004, Fundamental  of  the  Finite Element  Method  for  Heat  and  fluid  flow ,  John  Wiley  and  sons,  Chichester,  West Sussex, England.
R. Baigoria, L.O. Ottone, J.L. Pousa, F. Di Leo, J. Maranon, AOSTRA J.R., 8(1992)1-10.
R,  Salvius.  2014. Mathematical  Modelling  of  Groundwater  Flow.  University  of  Dar  es Saalam.
Singarimbun, A. A Numerical Model of Magmatic Hydrothermal System A Case Srudy of Kuju Volcano, Central Kyushu, Japan-Kyushu University, Fukuoka Japan. 1997