Modelling Dasar: Heat Equation 1 Dimensi


Baru-baru ini saya jadi tertarik untuk mempelajari NWP (Numerical Weather Prediction), karena itu kemudian saya membaca beberapa referensinya di internet dan di buku...dan saya benar-benar tidak mengerti. Dari sini saya mulai mencari tahu "kenapa bisa saya tidak mengerti?". Ternyata alasannya adalah karena saya belum melihat "ujung" dari pembelajaran seperti "mau dibawa kemana semua fakta-fakta ini" dan cara terburuk (dan paling ampuh mematikan motivasi belajar) adalah mempelajari fakta kosong (tak bermakna).

Singkatnya NWP ini adalah numerical modelling. Dan numerical modelling bertujuan untuk mensimulasikan kerja suatu sistem yang mana dari hasil simulasi ini diharapkan bisa didapatkan pemahaman baru tentang sistem atau memprediksi apa yang akan terjadi pada sistem kedepannya. Dan tentu saja dalam modelling harus benar-benar tahu apa yang ingin kita "lihat" (walaupun kedengarannya sepele, tapi ini cukup krusial).

Jadi, saya berencana untuk menguasai dasar dari numerical modelling dulu sebelum beranjak lebih jauh, untuk dapat intuisi dari pembelajaran NWP nanti. Dan pada postingan kali ini, yang akan di simulasikan adalah panas pada batang logam (1 dimensi). Di beberapa referensi kebanyakan simulasi inilah yang di gunakan untuk para pemula. Karena ini simulasi 1 dimensi dan sangat sederhana, jadi jangan harap hasilnya akan jadi keren seperti ini

atau ini



Ok, sebelum kita mulai "memprediksi" kita harus menetapkan apa yang hendak diprediksi itu. Jadi disini kita akan memprediksi suhu batang logam di 1 menit kedepannya yang akan di tampilkan dalam bentuk grafik. Panjang batang logam adalah 10 meter.

Selanjutnya, cari tahu hal-hal yang bisa mempengaruhi variabel yang ingin di prediksi. Dalam hal ini akan di cari suatu hukum fisika yang mengatur tentang perubahan panas pada batang logam. Dan kita punya:

1. Energi Panas
dimana:
Q = energi panas
m = massa
c = panas spesifik
Δt = perubahan suhu

rumus ini bisa di ubah menjadi

dimana:
ρ = massa jenis benda
A = luas penampang benda
Δx = panjang benda

2. Laju Transfer Panas

dimana:
q = laju transfer panas / area
k = konduktivitas termal
∇t = gradien suhu

karena perhitungan yang akan dilakukan adalah 1 dimensi dimana gradien yang ada hanya pada sumbu x, maka:

3. Kekekalan energi
Berdasarkan hukum kekekalan energi
Perubahan panas pada waktu Δt = panas masuk dari kiri - panas keluar dari kanan
Jadi
maaf tulisannya kecil karena keterbatasan tool
dimana:
s = waktu

Perhatikan di situ suhu (t) mengambil parameter waktu (s) dan lokasi (x). Ini dikarenakan suhu berubah tergantung pada waktu dan lokasinya.

kemudian persamaan ini di sederhanakan menjadi
ambil limit Δs dan Δx -> 0, maka
dengan
maka


K inilah yang disebut difusivitas termal. Persamaan inilah yang akan menjadi governing equation pada simulasi kali ini.

Selanjutnya menentukan kondisi awal dan kondisi batas. Kondisi awal adalah kondisi pada saat waktu belum berjalan (s = 0) sedangkan kondisi batas adalah kondisi pada ujung-ujung batas yang di amati (pada ujung 0 meter dan 10 meter logam batang yang akan disimulasikan).

Rasanya penjelasan di atas akan kabur jika tidak di barengi dengan grid. Perhatikan pada governing equation yang sudah di temukan tadi, itu merupakan persamaan diferensial parsial dan Jadinya disini suhu di turunkan terhadap waktu dan sumbu x. Oleh karena itu dibuatkan sebuah grid dengan x pada sumbu horizontal dan s pada sumbu vertikal.

α sebagai jarak antar grid waktu dan h sebagai jarak antar grid x. Angka-angka tersebut (t) menunjukkan nilai pada titik grid. Sistem koordinat pada grid ini adalah t(i,j) dimana i menunjukkan posisinya pada sumbu horizontal (x) dan j untuk sumbu vertikal (s). ps: maaf ilustrasinya tidak rapi.
Angka-angka yang ada di bawah grid adalah suhu pada kondisi awal (t(i,0)). Dua angka 0 pada masing-masing ujung grid itu mungkin terlihat janggal. Jadi pada kasus kali ini kedua ujung batang logam tersebut akan selalu 0 derajat celcius suhunya dan inilah yang dinamakan kondisi batas.

Kita tahu turunan adalah laju perubahan instan, yang mana nilai "perubahan instan" ini didapatkan dengan membandingkan 2 data yang sangat sangat sangat berdekatan sekali jarak waktunya, untuk ilustrasinya sudah pernah dibahas sebelumnya. Oleh karena itu bersama dengan grid yang sudah dibuat sebelumnya diatas, turunan bisa didapatkan dengan suatu pendekatan seperti
penurunan pendekatan forward difference
Yang merupakan pendekatan untuk penurunan suhu terhadap waktu, dan untuk turunan kedua suhu terhadap sumbu x adalah
penurunan kedua pendekatan central difference
Dengan begini maka, governing equation yang sebelumnya telah ditemukan dibuat bentuk pendekatannya yaitu
Karena yang ingin di ketahui dalam simulasi ini adalah prediksi suhu kedepannya pada batang logam, maka persamaan ini dibentuk ulang menjadi
Ok, itu adalah rumus terakhir yang dibutuhkan dalam simulasi kali ini. Sekarang mari deklarasikan nilai-nilai konstanta pada persamaan terakhir ini. Panjang batang logam 10 meter dan grid pada sumbu x ada 10 titik, maka h = 1 meter. Pada simulasi ini ingin di prediksi suhu batang logam 1 menit kedepannya dan grid sumbu s akan dibagi 60 bagian sehingga tiap titik grid pada sumbu s akan bernilai 1, ya s = 1 detik. Dan yang terakhir, batang logam yang ingin di simulasikan ini adalah besi yang mempunyai difusivitas termal 2.3 × 10−5 m²/det menurut sumbernya.

All set up done untuk rumusan di atas kertas, sekarang ayo beralih ke koding untuk memprogram perhitungan yang sudah dirumuskan tadi. Kali ini saya memakai bahasa fortran dalam menyusun programmnya dan gnuplot untuk membuat grafik hasil akhirnya.

Source code program perhitungan bisa di akses di sini. Jadi program ini akan mengambil input dari file (datasuhu.dat) dan output ke dalam file lain (prediksi.dat). prediksi.dat isinya sebagai berikut

 0  
 25.3  
 26.0  
 28.0  
 30.5  
 33.7  
 31.2  
 27.2  
 26.4  
 25.0  
 0  

Kemudian compile source code tadi dan jalankan

 ~$ gfortran heat.f95 -o heat  
 ~$ ./heat  

Output dari programnya (prediksi.dat) akan seperti ini


  #format  
  #x waktu suhu  
  0  0   0.0000000  
  0  1   25.2999992  
  0  2   26.0000000  
  0  3   28.0000000  
  0  4   30.5000000  
 .  
 .  
 .  
 ...dan terus berlanjut

Setelah ini di plot dengan gnuplot. Caranya

 ~$ gnuplot  
 gnuplot> set hidden3d  
 gnuplot> set dgrid3d 50,50 qnorm 2  
 gnuplot> splot 'prediksi.dat' with lines  

Lalu tampaklah hasilnya seperti ini
dan ini
Sepertinya...grafik ini tidak memberi kita pencerahan apa-apa. Terlihat seperti tidak ada perubahan sama sekali pada data output. Ya, itu karena range data suhu (0 - 33 derajat celcius) yang diberikan jika dibandingkan dengan perubahan suhu yang terlihat (perubahan suhunya kira-kria bermain pada 10-4 atau 10-3) akan sangat beda jauh dimana perubahannya sangat kecil. Dan terlebih lagi panjang batang logam yang di simulasikan (10 meter) itu terlalu panjang sehingga distribusi panasnya akan sangat lambat. Oleh karena itulah, panjang batang logam yang disimulasikan akan diperpendek menjadi 0,1 m atau 10 cm yang kemudian akan dibagi 10 untuk jarak antar gridnya (h = 1cm). Kemudian data suhu yang menjadi input akan di perkecil mendekati dengan kondisi batas seperti

 0  
 0.3  
 0.5  
 0.7  
 0.8  
 1.0  
 0.8  
 0.7  
 0.5  
 0.3  
 0  

Source code versi 2 bisa di akses di sini. Dan berikut grafik hasil dari perubahan yang sudah dilakukan tadi (script untuk gnuplot bisa di akses di sini).

Dari sini bisa di ambil kesimpulan kalau semakin tinggi suhu maka semakin cepat kehilangan panas, begitu juga sebaliknya, semakin rendah suhu maka akan semakin lambat kehilangan panas dan cepat lambarnya kehilangan panas ini tidak linier. Pada saat s = tak hingga , suhu batang logam seharusnya rata 0 derajat celcius.

yup, masih ada tittle yang mengganggu di samping :p
grafik animasi perkembangan suhu pada batang logam besi.


Dan berakhir lah sudah modelling yang sangat sederhana ini. Walaupun begitu saya sendiri masih ragu dengan governing equation yang digunakan disini karena beberapa sumber punya rumus yang beda-beda tapi saya ambil rumus ini karena make a sense di saya. Selain itu masih banyak kekurangan juga pada simulasi kali ini seperti penggunaan simbol yang tidak biasa (s = waktu? t = suhu?) dan ilustrasi-ilustrasi yang masih belum rapi. Ok, sekian.

referensi:
http://digilib.its.ac.id/public/ITS-paper-29484-1206100701-Paper.pdf
http://www.uio.no/studier/emner/matnat/ifi/INF2340/v05/foiler/sim02.pdf
http://spiff.rit.edu/classes/phys559/lectures/insulated_end/insulated_end.html
http://en.openei.org/wiki/Definition:Numerical_Modeling
https://en.wikipedia.org/wiki/Heat_equation
http://superuser.com/questions/716014/create-a-smooth-surface-using-x-y-z-data-gnuplot

Komentar