SIMULASI NUMERIK PERPINDAHAN PANAS KONVEKSI ALAMI PADA LAPIS BATAS ALIRAN LAMINAR DENGAN METODE BEDA HINGGA
SKRIPSI Diajukan sebagai salah satu syarat untuk memperoleh gelar Sarjana Teknik
Oleh: WENDY DESTYANTO NIM. I0401051
JURUSAN TEKNIK MESIN FAKULTAS TEKNIK UNIVERSITAS SEBELAS MARET SURAKARTA 2007
MOTTO DAN PERSEMBAHAN
MOTTO “You don’t know what you’ve got till it’s gone” (Counting Crows) “Bukan hasil yang menjadikan kita besar, tetapi proses yang membuat kita lebih besar”
PERSEMBAHAN Karya ini kupersembahkan untuk : Ibunda Sri Supraptiwi dan Ayahanda Waluyo Adikku : Dhimas Willy Ferdianto My Soulmate : Eka Wiziyanti Orang-orang yang senantiasa berdiri di belakangku
Wendy Destyanto. Komputasi Konversi Energi. SIMULASI NUMERIK PERPINDAHAN PANAS KONVEKSI ALAMI PADA ALIRAN LAMINAR DENGAN METODE BEDA HINGGA
Abstrak Penelitian ini dilakukan untuk menghitung nilai koefisien perpindahan panas konveksi pada aliran laminar, serta distribusi kecepatan dan distribusi temperatur udara disekitar plat datar vertikal panas. Penulisan program menggunakan bahasa pemrograman Fortran Power Station 4.0 dan divisualisasi dengan perangkat lunak Matlab R13. Temperatur plat yang disimulasikan adalah 50 oC dan temperatur arus bebas (free stream) 10 oC dengan panjang plat 3 cm. Sifat-sifat udara dihitung pada temperatur film 30 oC, yaitu : viskositas kinematik, υ = 16 x 10-6 m2/s dan konduktivitas termal, k = 26,38 x 10-3 W/m oC. Properti lain yang digunakan adalah percepatan gravitasi, g = 9,81 m/s2 dan Pr = 0,7 dengan kondisi batas di y = 0 adalah u = 0, v = 0 dan T = Tw; di y = ∞ adalah u = 0 dan T = T∞; dan di x = 0 adalah u = 0 dan T = T∞. Metode yang digunakan adalah metode beda hingga dengan diskritisasi dari persamaan energi, persamaan momentum dan persamaan kontinuitas . Hasil perancangan program dapat dijalankan untuk menghasilkan simulasi numerik distribusi kecepatan dan distribusi temperatur aliran laminar konveksi alami pada plat datar vertikal panas, dan distribusi nilai koefisien perpindahan panas konveksi alami lokal. Semakin jauh jarak dari ujung plat, nilai koefisien perpindahan panas konveksi alami lokalnya semakin kecil. Kata kunci : konveksi alami, udara, plat datar vertikal, metode beda hingga.
Wendy Destyanto. Computation of Energy Conversion. NUMERICAL SIMULATION FOR NATURAL CONVECTION HEAT TRANSFER IN THE LAMINAR FLOW AREA USING FINITE DIFFERENCE METHODE
Abstract The main objective of this study is to calculate the heat transfer coefficient for natural convection in the laminar flows area as well as velocity and temperature distributions along the plate for vertical, heated flat- plate. A Fortran Power Station 4.0 was written to obtain the heat transfer coefficient and the velocity and temperature distributions, and the results are ploted using Matlab R13. The plate temperature set at 50 oC and the free stream temperature set at o 10 C, with plate length 3 cm. The air properties are evaluated at film temperature 30 oC, which is viscous kinematic, v = 16 x 10-6 m2/s and thermal conductivity, k = 26,38 x 10-3 W/m oC. Other properties are gravitational acceleration, g = 9,81 m/s2 and Pr = 0,7, with boundary condition at y = 0 are u = 0, v = 0, T = Tw; at y = ∞ are u = 0 and T = T∞; and at x = 0 are u = 0 and T = T∞. A finite difference methode was used by discritizing the energy, momentum and continuity equations. The program was running properly and showing a good agreement with the classical literature in simulating the velocity and temperature distributions in laminar flows area for natural convection along a vertical, heated flat-plate, and also for the heat transfer convection coefficient. The more distances x from the leading edge of the plate, the less heat transfer convection coefficient obtained. Key words : natural convection, air, vertical-flat-plate, finite difference methode.
KATA PENGANTAR
Alhamdulillah. Puji syukur penulis panjatkan kepada Allah swt., yang telah memberikan rahmat, hidayah serta kekuatan kepada penulis, sehingga penulis dapat melaksanakan penelitian dan menyelesaikan laporan tugas akhir dengan judul “Simulasi Numerik Perpindahan Panas Konveksi Alami pada Lapis Batas Aliran Laminar dengan Metode Beda Hingga”, sebagai salah satu syarat untuk memperoleh gelar Sarjana Teknik di Jurusan Teknik Mesin Fakultas Teknik Universitas Sebelas Maret Surakarta. Dalam kesempatan ini, penulis ingin menyampaikan ucapan terima kasih dan penghargaan yang sebesar-besarnya kepada semua pihak yang telah memberikan bantuan, doa, dukungan dan semangat, baik moril maupun materiil kepada : 1. Ibunda, Ayahanda, Dhimas dan Eka Wiziyanti, yang tanpa jemu dan dengan sabar memberikan doa, semangat, keyakinan, nasehat dan cinta kasihnya. 2. Bapak Ir. Agustinus Sujono, MT., selaku Ketua Jurusan Teknik Mesin UNS. 3. Bapak Eko Prasetya Budiana, ST., MT., selaku Pembimbing I tugas akhir, atas bimbingan, nasehat, kepercayaan dan ilmu pengetahuan yang diajarkannya. 4. Bapak R. Lulus Lambang GH, ST., MT., selaku Pembimbing II tugas akhir, atas bimbingan, kesabaran dan ilmu pengetahuan yang diajarkannya. 5. Bapak Dwi Aries Himawanto, ST., MT., selaku Pembimbing Akademik, atas saran dan nasehatnya. 6. Bapak-bapak dosen dan staf karyawan di lingkungan Teknik Mesin UNS, atas didikan, nasehat, ilmu yang diajarkan dan kerjasamanya. 7. Keluarga besar H. Abdul Kadir HS, SE.; pakde, bude, Mbak Tika, Abang, Mas Caca, Mbak Ayu, Mbak Illa, Mas Zuchri dan Dinda “Endo’” Virajati, yang telah memberikan nasehat, kepercayaan dan dukungan yang besar kepada penulis. 8. Bude Yam, Pakde Joni, kedua Mbah Putri, om dan tante, pakde dan bude, dan sepupu-sepupu penulis, atas doa, dukungan dan dorongan semangatnya.
9. Teman-teman almamater Mahasiswa Mesin angkatan 2001 atas kerjasama, dukungan, keceriaan dan petualangan dari pantai ke pantai yang tak terlupakan; Adit, Ali, Andy P, Andy “supit”, Aris “jenggot”, Arif “ini-itu”, Bambang “bams”, Bambang “bombot”, Budi “cobrut”, Fendy “fenoy”, Fauzi, Jamal “Lou Han”, Joko, Irawan, Imbar, Imam W, Imam Fahad, Irvan “irfun”, Kurniawan, Rizka “karjo”, Risharyanto, Hadi, Said “swat”, “si gede” Sulistyana, Taufik, Uki “ukri”, Tri Wahyudi “kentang”, kakak-adik angkatan Teknik Mesin UNS, Let’s get Solidarity M Forever!! 10. Teman-teman ex-Asrama “Ceria” UNS; Halim “hamil” si lugu tak tahu malu, Aris “sipit”, Indro, Bayu “Cikibul”, Wahyu “Dabull”, Latief, Imam “gundul”, Andika, Toni, Teguh “san”, Okta “kuro”, Pras “sube” buat tiket kretanya, Timbul, Agus “SH”, Raden “SH”, Bowo “pak RT”, Jaka “kecu” yang sering ngajarin komputer, Arfan, Dedi “oe”, Dhana “grepes” sohib dari SMA, Satir, Wisnu “bose”, Gendenk, Azis, Rojak, Doa, Filmon, Wiwid, Rere, Pak Dal. Salam Ceria !! 11. Semua pihak yang belum sempat disebutkan, yang telah membantu penelitian dan penyusunan laporan tugas akhir ini. Akhirnya, penulis menyadari bahwa karya kecil ini masih memiliki kekurangan dan kelemahan. Sehingga kritik dan saran penulis harapkan demi perbaikan dan pembelajaran untuk penelitian selanjutnya. Terimakasih.
Surakarta, Januari 2007
penulis
DAFTAR ISI ABSTRAK ……………………………………………………………….. KATA PENGANTAR ..………………………………………………….. DAFTAR ISI ……………………………………………………………... DAFTAR GAMBAR …………………………………………………….. BAB I
BAB II
iv vi viii x
PENDAHULUAN 1.1 Latar Belakang …………………………………………..... 1.2 Perumusan Masalah .............................................. 1.3 Batasan Masalah ………………………………………….. 1.4 Tujuan Penelitian …………………………………………. 1.5 Manfaat Penelitian ………………………………………... 1.6 Sistematika Penulisan ……………………………………..
1 2 2 2 2 3
LANDASAN TEORI 2.1 Tinjauan Pustaka …………………………………………. 2.2 Konveksi Alami (Natural Convection) ...………………… 2.3 Lapis Batas (Boundary Layer) ...................…………. 2.4 Metode Beda Hingga .................…………………………. 2.4.1 Pendekatan Beda Maju Orde Pertama ................ 2.4.2 Pendekatan Beda Mundur Orde Pertama ………... 2.4.3 Pendekatan Beda Tengah Orde Pertama ............. 2.4.4 Pendekatan Beda Tengah Orde Kedua ................ 2.5 Persamaan Lapis Batas Pada Plat Datar .…………………. 2.6 Angka Grashof dan Angka Rayleigh .......................... 2.7 Koefisien Perpindahan Panas Konveksi Alami ..............
4 4 4 6 6 7 8 8 9 10 10
BAB III PELAKSANAAN PENELITIAN 3.1 Alat dan Bahan ...........................……………......… 3.1.1 Alat ....................................................... 3.1.2 Bahan ..................................................... 3.2 Garis Besar Penelitian ........................................... 3.3 Diskritisasi Persamaan Lapis Batas Dalam Bentuk Tak Berdimensi (Dimensionless) ............................... 3.3.1 Diskritisasi Persamaan Energi ......................... 3.3.2 Diskritisasi Persamaan Momentum ................... 3.3.3 Diskritisasi Persamaan Kontinuitas ................... 3.4 Kondisi Batas ..................................................... 3.5 Perhitungan Koefisien Perpindahan Panas Konveksi Alami 3.6 Penyusunan Algoritma dan Bagan Alir (Flow Chart) Program ............................. BAB IV DATA DAN ANALISIS 4.1 Validasi Program ............…………………….......……... 4.2 Simulasi Konveksi Alami Plat Datar Vertikal Panas ……...
12 12 12 12 14 16 18 21 23 23 24 27 28
BAB V KESIMPULAN DAN SARAN 5.1 Kesimpulan ……………………………………………….. 5.2 Saran ………………………………………………………
35 36
DAFTAR PUSTAKA …………………………………………………….. LAMPIRAN ……………………………………………………………….
37 38
DAFTAR GAMBAR
Gambar 2.1
Daerah Lapis Batas, (a) Profil Kecepatan, dan (b) profil temperatur pada konveksi alami
Gambar 2.2
Ilustrasi Pendekatan Beda Maju Orde Pertama
Gambar 2.3
Ilustrasi Pendekatan Beda Mundur Orde Pertama
Gambar 2.4
Ilustrasi Pendekatan Beda Tengah Orde Pertama
Gambar 2.5
Ilustrasi Pendekatan Beda Tengah Orde Kedua
Gambar 3.1
Diagram Alir Penelitian
Gambar 3.2
Grid yang digunakan dalam analisa
Gambar 3.3
Grid untuk Derivasi Pendekatan Beda Hingga
Gambar 3.4
Nodal untuk Diskritisasi Persamaan Kontinuitas
Gambar 3.5
Kondisi Batas
Gambar 3.6
Diagram Alir Program
Gambar 4.1
Kondisi Batas Penelitian Rolando A. Chavez
Gambar 4.2
Grid Penelitian Rolando A. Chavez
Gambar 4.3
Profil Kecepatan Hasil Penelitian Rolando A. Chavez
Gambar 4.4
Profil Kecepatan pada Penelitian
Gambar 4.5
Profil Temperatur Hasil Penelitian Rolando A. Chavez
Gambar 4.6
Profil Temperatur pada Penelitian
Gambar 4.7
Distribusi Kecepatan Hasil Penelitian Rolando A. Chavez
Gambar 4.8
Distribusi Kecepatan pada Penelitian
Gambar 4.9
Distribusi Temperatur Hasil Penelitian Rolando A. Chavez
Gambar 4.10 Distribusi Temperatur pada Penelitian Gambar 4.11 Distribusi Koefisien Perpindahan Panas Konveksi Lokal
BAB I PENDAHULUAN
1.1 Latar Belakang Masalah Persoalan perpindahan panas serta metode penyelesaiannya mengalami perkembangan pesat diberbagai bidang kehidupan. Bidang teknologi industri banyak menggunakan prinsip-prinsip dasar proses perpindahan panas. Sehingga pendalaman
di
bidang
ini
perlu
ditingkatkan,
terutama
pada
metode
penyelesaiannya. Metode yang lebih cepat, akurat dengan sedikit kesalahan sangat dibutuhkan untuk mendapatkan hasil yang lebih cepat. Perpindahan panas (heat transfer) adalah ilmu untuk memprediksikan perpindahan energi yang terjadi akibat perbedaan suhu pada benda atau material. Proses perpindahan panas dapat terjadi melalui tiga cara, yaitu perpindahan panas secara konduksi, konveksi dan radiasi. Perpindahan panas konveksi adalah perpindahan panas yang terjadi di antara permukaan benda dengan fluida yang bergerak, karena terdapat gradien suhu diantara keduanya. Konveksi alami terjadi karena adanya perubahan densitas (kerapatan) fluida akibat proses pemanasan, yang menyebabkan fluida bergerak ke atas. Gerakan fluida pada konveksi alami (baik gas maupun zat cair) terjadi karena gaya apung (buoyancy force) yang timbul apabila densitas fluida berkurang akibat proses pemanasan. Konveksi alami memegang peranan penting dalam rekayasa industri, seperti pada perancangan alat penukar kalor, pendinginan transformator, dan komponen elektronika. Penelitian mengenai fenomena konveksi alami telah banyak dilakukan, baik secara eksperimen di laboratorium maupun secara numerik. Penelitian secara eksperimen di laboratorium untuk mengetahui fenomena yang terjadi pada proses konveksi alami membutuhkan biaya yang mahal dan proses yang cukup rumit. Oleh karena itu, dikembangkanlah suatu penelitian mengenai metode penyelesaian dengan biaya yang jauh lebih rendah serta waktu yang lebih cepat, yaitu dengan metode simulasi numerik yang didasarkan pada metode bedahingga (finite-difference methode).
Pada konveksi alami, kecepatan aliran fluidanya sangat rendah. Dan pada aliran dengan kecepatan rendah, aliran laminar akan lebih sering terbentuk dibandingkan dengan aliran Turbulen. Oleh sebab itu penelitian ini difokuskan pada aliran laminar.
1.2 Perumusan Masalah Perumusan masalah dalam penelitian ini adalah bagaimana mensimulasikan secara numerik perpindahan panas konveksi alami plat datar vertikal di daerah lapis batas aliran laminar dengan metode beda hingga.
1.3 Batasan Masalah Masalah pada penelitian ini dibatasi pada persoalan konveksi alami pada plat datar vertikal panas yang diselesaikan dengan menggunakan metode beda hingga untuk memperoleh distribusi kecepatan, distribusi temperatur dan koefisien perpindahan panas konveksi alami dengan udara sebagai fluida penghantar kalor.
1.4 Tujuan Penelitian Penelitian ini bertujuan untuk menampilkan distribusi nilai koefisien perpindahan panas konveksi alami lokal, distribusi kecepatan udara dan distribusi temperatur udara pada lapis batas aliran laminar secara kualitatif, dengan menggunakan metode beda-hingga sebagai alternatif metode penyelesaian persoalan perpindahan panas.
1.5 Manfaat Penelitian Manfaat dari penelitian ini adalah : a. Mengembangkan dan menerapkan ilmu pengetahuan terutama ilmu pengetahuan
Komputasi
Perpindahan
Panas,
Metode
Numerik,
Perpindahan Panas dan Mekanika Fluida yang diperoleh di bangku kuliah b. Mempelajari fenomena konveksi alami yang terjadi pada plat datar vertikal panas. c. Sebagai dasar pengembangan untuk penelitian yang lebih kompleks.
1.6 Sistematika Penulisan Sistematika penulisan yang digunakan adalah : BAB I : PENDAHULUAN Berisi latar belakang masalah, batasan dan perumusan masalah, tujuan dan manfaat penelitian serta sistematika penulisan. BAB II : LANDASAN TEORI Berisi tentang tinjauan pustaka, dasar teori konveksi alami dan penjelasan mengenai metode beda hingga. BAB III : PELAKSANAAN PENELITIAN Berisi tentang alat dan bahan yang digunakan dalam penelitian, tata dan cara penelitian, penurunan persamaan kontinuitas, persamaan momentum dan persamaan energi dengan metode beda hingga, dan diagram alir program. BAB IV : HASIL DAN PEMBAHASAN Berisi data hasil penelitian (simulasi) dan pembahasannya. BAB V : PENUTUP Berisi kesimpulan penelitian dan saran-saran untuk penelitian selanjutnya. DAFTAR PUSTAKA LAMPIRAN
BAB II LANDASAN TEORI
2.1 Tinjauan Pustaka Rolando A Chavez (2004) menggunakan metode numerik untuk menghitung koefisien perpindahan panas konveksi alami dan menyelesaikan persamaan lapis batas plat datar vertikal panas pada fluida superkritis. Chavez menggunakan bahasa pemrograman Fortran untuk menghitung kecepatan dan temperatur di sepanjang plat. Hasil dari penelitian Chavez menunjukkan bahwa profil kecepatan meningkat dengan bertambahnya jarak tegak lurus terhadap plat sampai mencapai titik maksimum, kemudian menurun sampai mencapai batas lapis aliran. Profil temperatur mengalami penurunan dengan semakin bertambahnya jarak tegak lurus terhadap plat, hingga menjadi sama dengan temperatur arus bebas.
2.2 Konveksi Alami (Natural Convection) Sudah umum diketahui bahwa plat logam panas akan menjadi lebih cepat dingin ketika diletakkan di depan kipas angin dibandingkan dengan ketika diletakkan di udara tenang. Dikatakan bahwa kalor dikonveksi atau diili keluar dan proses terjadinya perpindahan panas ini disebut perpindahan kalor secara konveksi atau ilian (Holman, 1997). Perpindahan panas konveksi adalah proses perpindahan energi dari permukaan benda ke fluida yang mengalir di atasnya karena perbedaan suhu diantara keduanya (benda – fluida) (Oosthuizen, 1999). Konveksi alami terjadi karena fluida, yang karena proses pemanasan, berubah densitasnya (kerapatannya) dan bergerak ke atas. Gerakan fluida pada konveksi alami terjadi karena gaya apung (buoyancy force) yang timbul apabila densitas fluida berkurang akibat proses pemanasan (Holman, 1997).
2.3 Lapis Batas (Boundary Layer) Daerah lapis batas (Boundary Layer) adalah daerah atau lapisan tipis yang dekat permukaan benda, yang masih berada dalam pengaruh viskositas fluida dan
perpindahan panas fluida. Tebal lapis batas dibagi menjadi dua, yaitu lapis batas kecepatan dan lapis batas termal. Tebal lapis batas kecepatan (δ) adalah jarak yang diukur dari permukaan benda sampai suatu titik dimana efek viskositas sudah tidak berpengaruh lagi. Tebal lapis batas termal (δT) adalah jarak yang diukur dari permukaan benda sampai suatu titik dimana efek perpindahan panas sudah tidak berpengaruh. Fluida disekitar permukaan plat panas menjadi lebih ringan dibandingkan dengan fluida yang lebih jauh dari permukaan plat. Sifat ringan fluida ini yang menyebabkan terjadinya pergerakan ke atas, bergesekan dengan dinding dan memindahkan panas dari dinding (panas). Fluida yang jauh dari dinding tidak terpengaruh oleh efek panas yang ditimbulkan plat (Bejan, 1993). Pada plat datar vertikal panas (gambar 2.1a) akan terbentuk lapis-batas (boundary layer) konveksi alami. Pada dinding, kecepatannya adalah nol karena terdapat kondisi tanpa gelincir (no-slip). Kecepatan bertambah sampai mencapai suatu nilai maksimum, kemudian turun secara bertahap sampai nol pada tepi lapisan batas. Perkembangan awal lapisan batas adalah laminar, tetapi pada suatu jarak tertentu dari tepi depan, bergantung pada sifat-sifat fluida dan beda suhu antara dinding dan lingkungan, terbentuk pusaran-pusaran dan transisi ke lapisan batas turbulen akan terbentuk (Ozisik, 1988). Lapis batas
Lapis batas Termal
Profil kecepatan Plat vertikal panas
Plat vertikal panas
X
X
u
u
V
.y
(a)
Profil Temperatur
V
.y
(b)
Gambar 2.1 Daerah Lapis Batas, (a) Profil Kecepatan, (b) profil temperatur pada konveksi alami
2.4 Metode Beda Hingga Salah satu metode penyelesaian persamaan lapis batas adalah dengan metode beda hingga. Metode beda hingga merupakan suatu cara, selain metode elemen hingga, untuk menentukan penyelesaian numerik dari persamaanpersamaan diferensial parsial. Metode beda hingga didasarkan pada ekspansi deret Taylor, yaitu metode pendekatan agar sebuah persamaan diferensial parsial dapat diubah menjadi operasi aritmatika dan operasi logika yang dapat dibaca oleh komputer (Hoffmann, 1989). Ekspansi deret Taylor menghasilkan pendekatan beda maju orde pertama, beda mundur orde pertama, beda tengah orde pertama dan beda tengah orde kedua. 2.4.1 Pendekatan Beda Maju Orde Pertama Ekspansi deret Taylor untuk f(x + ∆x) pada x : f (x + ∆x ) = f ( x ) + (∆x ) ∞
= f (x ) + ∑ n =1
penyelesaian untuk
∂f (∆x ) ∂ 2 f (∆x ) ∂ 3 f + ... + + 2! ∂x 2 3! ∂x 3 ∂x 3
2
(∆x )n
∂n f n! ∂x n
(2.1)
∂f diperoleh : ∂x
∂f f ( x + ∆x ) − f (x ) ∆x ∂ 2 f (∆x ) ∂ 3 f + ... = − − 2! ∂x 2 3! ∂x 3 ∂x ∆x 2
∂f f ( x + ∆x ) − f ( x ) = + O(∆x ) ∂x ∆x
atau bisa ditulis
f − fi ∂f + O(∆x ) = i +1 ∆x ∂x i
Persamaan di atas disebut dengan pendekatan beda maju orde pertama.
(2.2)
f i −1
fi
f i +1
(−)
(+)
∆x
Gambar 2.2 Ilustrasi Pendekatan Beda Maju Orde Pertama
2.4.2 Pendekatan Beda Mundur Orde Pertama Ekspansi deret Taylor untuk f(x - ∆x) pada x :
∂f (∆x ) ∂ 2 f (∆x ) ∂ 3 f + − + ... ∂x 2! ∂x 2 3! ∂x 3 + untuk n genap n n ∞ ⎡ ⎤ (∆x ) ∂ f = f ( x ) + ∑ ⎢± ⎥ n! ⎦ ∂x n n =1 ⎣ - untuk n ganjil 2
f ( x − ∆x ) = f ( x ) − (∆x )
penyelesaian untuk
3
(2.3)
∂f diperoleh : ∂x
∂f f ( x ) − f ( x − ∆x ) = + O(∆x ) ∂x ∆x atau bisa ditulis
f − f i −1 ∂f = i + O(∆x ) ∂x i ∆x
(2.4)
Persamaan di atas disebut pendekatan beda mundur orde pertama.
f i −1
fi
(− )
(+ )
f i +1
∆x
Gambar 2.3 Ilustrasi Pendekatan Beda Mundur Orde Pertama
2.4.3 Pendekatan Beda Tengah Orde Pertama
Dengan mengurangkan ekspansi deret Taylor untuk f(x + ∆x) (persamaan (2.1)) dengan ekspansi deret Taylor untuk f(x - ∆x) (persamaan 2.3)), diperoleh :
f ( x + ∆x ) − f ( x − ∆x ) = 2∆x penyelesaian untuk
(∆x ) ∂ 3 f + ... ∂f +2 3! ∂x 3 ∂x 3
(2.5)
∂f diperoleh : ∂x
f ( x + ∆x ) − f ( x − ∆x ) ∂f 2 + O(∆x ) = 2∆x ∂x atau bisa ditulis
f − f i −1 ∂f 2 = x+i + O(∆x ) ∂x i 2 ∆x
(2.6)
Persamaan di atas disebut pendekatan beda tengah orde pertama.
f i −1
fi
(− )
f i +1 (+ )
∆x
Gambar 2.4 Ilustrasi Pendekatan Beda Tengah Orde Pertama
2.4.4 Pendekatan Beda Tengah Orde Kedua
Pendekatan beda hingga untuk persamaan turunan orde yang lebih tinggi dapat ditentukan dengan menambahkan ekspansi deret Taylor pada f(x). Dengan menambahkan persamaan (2.1) dan persamaan (2.3) didapat : 2 ( ∆x ) f ( x + ∆x ) + f ( x − ∆x ) = 2∆x + 2
( ∂2 f ∆x ) ∂ 4 f +2 . + .. . 2! ∂x 2 4! ∂x 4
∂2 f diperoleh : penyelesaian untuk ∂x 2 ∂2 f f ( x + ∆x ) − 2 f (x ) + f ( x − ∆x ) 2 = + O(∆x ) 2 2 ∂x ∆x
4
(2.7)
atau bisa ditulis
∂2 f ∂x 2
= i
f i +1 − 2 f i + f i −1 2 + O(∆x ) 2 ∆x
(2.8)
Persamaan di atas disebut Pendekatan Beda Tengah Orde Kedua.
f i −1
fi
(− )
f i +1 (+ )
∆x
Gambar 2.5 Ilustrasi Pendekatan Beda Tengah Orde Kedua
2.5 Persamaan Lapis Batas Pada Plat Datar
Persamaan lapis batas yang berlaku pada perpindahan panas konveksi alami untuk plat datar pada kondisi tunak (steady), tak mampu mampat (incompressible) adalah sebagai berikut : a) Persaman Kontinuitas :
∂u ∂v + =0 ∂x ∂y
(2.9)
Pada konveksi alami berlaku pendekatan Boussinesq, yaitu dalam analisa mengenai aliran pada konveksi alami, sifat-sifat fluida diasumsikan konstan kecuali perubahan densitas terhadap temperatur yang menyebabkan munculnya gaya apung (buoyancy force) (Oosthuizen, 1999). Meskipun gerakan fluida akibat perbedaan densitas, tetapi angka perbedaannya sangat kecil. Sehingga dapat diperoleh penyelesaian dengan mengandaikan aliran tak mampu mampat (Incompressible) (Holman, 1997). b)
Persamaan Momentum : u
∂u ∂u ∂ 2u +v = β g (T − T∞ ) + υ 2 ∂x ∂y ∂y
(2.10)
Dengan asumsi dasar teori lapis batas, apabila (δ/L) kecil maka (v / u ∞ ) bernilai kecil. Untuk konveksi alami pada lapis batas, komponen kecepatan arah y, (v), memiliki besaran yang sangat kecil dibandingkan komponen kecepatan u. Sehingga momentum pada arah y dapat diabaikan (Oosthuizen, 1999).
c)
Persamaan Energi
: u
∂T ∂T ∂ 2T +v =α 2 ∂x ∂y ∂y
(2.11)
2.6 Angka Grashof dan Angka Rayleigh
Angka Grashof adalah satuan rasio perbesaran gaya apung (buoyancy force) terhadap viskositas pada aliran konveksi alami. Secara matematis dituliskan sebagai : Gr =
βg (Tw − T∞ ) L3 v2
(2.12)
Angka Rayleigh didefinisikan sebagai satuan tak berdimensi hasil kali antara angka Grashof dengan angka Prandtl (Pr), yang dirumuskan sebagai : Ra = Gr. Pr = Dimana,
β .g (Tw − T∞ )L3 . Pr υ2
Ra
= Angka Rayleigh
β
= Koefisien ekspansivitas termal ( 1 / oC )
g
= Percepatan gravitasi ( m / s2 )
(2.13)
Tw - Tf = Selisih temperatur plat dengan fluida ( oC ) L
= Panjang plat ( m )
Pr
= Angka Prandtl
υ
= Viskositas kinematik ( m2 / s )
Angka Rayleigh digunakan sebagai salah satu acuan untuk menentukan jenis aliran dalam konveksi alami, yaitu : Ra < 109 : Aliran Laminar Ra = 109 : Aliran Transisi Ra > 109 : Aliran Turbulen
2.7 Koefisien Perpindahan Panas Konveksi Alami
Koefisien perpindahan panas (h) berpengaruh terhadap laju perpindahan panas pada suatu sistem konveksi. Besarnya nilai h dipengaruhi oleh jenis fluida yang digunakan, bentuk permukaan yang dilewati fluida dan kecepatan fluidanya
(laminar, turbulen atau transien). Viskositas mempengaruhi profil kecepatan yang akan berpengaruh terhadap laju perpindahan energi pada daerah dinding (Holman, 1997) Persamaan koefisien perpindahan panas konveksi : hx = −
⎛ ∂T ⎞ k ⎜ ⎟ (Tw − T∞ ) ⎜⎝ ∂y ⎟⎠ y →0
(2.14)
Persamaan di atas disebut sebagai persamaan koefisien perpindahan panas konveksi lokal. Tanda negatif menyatakan gradien temperatur terhadap y mengalami penurunan. Semakin besar y, semakin kecil gradien temperaturnya.
BAB III PELAKSANAAN PENELITIAN
3.1 Alat dan Bahan
3.1.1 Alat a. Komputer pribadi dengan spesifikasi : - Prosesor AMD Sempron 2400+ - Memori DDR RAM 256 MB b. Perangkat lunak Microsoft Fortran PowerStation 4.0 c. Perangkat lunak Matlab 6.5.1 (R13) d. Printer 3.1.2 Bahan Hasil diskritisasi persamaan kontinuitas, persamaan momentum dan persamaan energi dengan metode Beda Hingga dan persamaan koefisien perpindahan panas konveksi.
3.2 Garis Besar Penelitian
Penelitian yang dilakukan menggunakan metode studi pustaka dengan langkah pelaksanaan secara garis besar sebagai berikut : a. Mengumpulkan literatur berupa hasil-hasil penelitian terdahulu dan buku penunjang. b. Mempelajari literatur 1. Mempelajari penelitian-penelitian yang pernah dilakukan. 2. Mempelajari persamaan lapis batas yang digunakan. c. Merencanakan algoritma program 1. Membuat diskritisasi persamaan lapis batas dalam bentuk tak berdimensi (Dimensionless Form). 2. Menyusun bagan alir program. d. Menulis bagan alir dalam bahasa program (Fortran). e. Menjalankan program. f. Memperbaiki kesalahan dalam pemrograman 1. Kesalahan penulisan
2. Kesalahan algoritma g. Membuat visualisasi hasil program dengan perangkat lunak Matlab. h. Menyusun laporan. Diagram alir penelitian yang dilakukan adalah : Mulai Mengumpulkan literatur Mempelajari literatur Membuat diskritisasi persamaan lapis batas tak berdimensi : • • •
Persamaan Energi Persamaan Momentum Persamaan Kontinuitas
Membuat algoritma program
Menulis bagan atur dalam bahasa Fortran
• • •
Menjalankan program : Distribusi Kecepatan Distribusi Temperatur Menghitung Koefisien Perpindahan Panas
• • •
Membuat visualisasi : Distribusi Kecepatan Distribusi Suhu Distribusi Koefisien Perpan
Program benar ya A Gambar 3.1 Diagram alir penelitian
tidak
A Menyusun Laporan
Selesai Gambar 3.1 (lanjutan)
3.3 Diskritisasi Persamaan Lapis Batas Dalam Bentuk Tak Berdimensi (Dimensionless)
Untuk menyederhanakan penyelesaian derivasi persamaan lapis batas, digunakan variabel referensi yang mengubah persamaan lapis batas menjadi persamaan lapis batas tak berdimensi (dimensionless). .x Y
.i = nx X
.y
.i = 1 j=1
j = ny
Gambar 3.2 Grid yang digunakan dalam analisa Variabel referensi tak berdimensi yang digunakan adalah sebagai berikut : • X =
x L.G
x = X .L.G
∂x = L.G.∂X
• Y=
y W
y = W .Y ∂y = W∂Y ∂y 2 = W 2 ∂ 2Y
⎛ uW ⎞⎛ W ⎞ • U =⎜ ⎟⎜ ⎟ ⎝ υ ⎠⎝ LG ⎠
u=
υ .L.G.U
∂u =
W2
υ .L.G
∂ 2u = • V =
vW
υ
v=
T − T∞ Tw − T∞
∂U
υ.L.G W
2
∂ 2U
υ .V
∂v = • θ=
W2
W
υ W
∂V
T − T∞ = θ (Tw − T∞ ) T = θ (Tw − T∞ ) + T∞ ∂T = (Tw − T∞ )∂θ
∂ 2T = (Tw − T∞ )∂ 2θ
β ref .g .(Tw − T∞ )W 4 υ 2 .L
• G=
• β ref =
1 Tf
• Tf =
1 (Tw + T∞ ) 2
• β* =
β β ref
• α=
υ Pr
g=
G.υ 2 .L β ref (Tw − T∞ )W 4
∆y
i, j-1
∆y
i, j
i, j+1 ∆x
i-1, j Gambar 3.3 Grid untuk Derivasi Pendekatan Beda Hingga 3.3.1 Diskritisasi Persamaan Energi
Persamaan dasar energi :
u
∂T ∂T ∂ 2T +v =α 2 ∂x ∂y ∂y
(3.1)
Dengan mensubstitusikan variabel tak berdimensi, tiap suku dari persamaan di atas dapat diubah ke bentuk persamaan tak berdimensi sebagai berikut :
•
u
∂T ∂x
⎛ υ (Tw − T∞ ) ⎞ ∂θ =⎜ ⎟U 2 ⎝ W ⎠ ∂X
(3.2)
•
v
∂T ∂y
⎛ υ (Tw − T∞ ) ⎞ ∂θ =⎜ ⎟V 2 ⎝ W ⎠ ∂Y
(3.3)
•
α
∂ 2T ⎛ υ (Tw − T∞ ) ⎞ 1 ∂ 2θ = ⎜ ⎟ 2 2 ∂y 2 ⎝ W ⎠ Pr ∂Y
(3.4)
Substitusi persamaan (3.2), (3.3) dan (3.4) ke persamaan (3.1), diperoleh :
⎛ υ (Tw − T∞ ) ⎞ ∂θ ⎛ υ (Tw − T∞ ) ⎞ 1 ∂ 2θ ⎛ υ (Tw − T∞ ) ⎞ ∂θ + = U V ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ 2 2 2 2 ⎝ W ⎠ ∂Y ⎝ W ⎠ Pr ∂Y ⎝ W ⎠ ∂X ⎛ υ (Tw − T∞ ) ⎞ Persamaan di atas disederhanakan dengan mengeliminasi ⎜ ⎟ , menjadi : 2 ⎝ W ⎠ U
1 ∂ 2θ ∂θ ∂θ +V = ∂X ∂Y Pr ∂Y 2
(3.5)
Persamaan (3.5) disebut persamaan energi dalam bentuk tak berdimensi. Diskritisasi tiap suku persamaan di atas adalah sebagai berikut :
•
U
∂θ ∂X
⎛ θ ik, j − θ ik−1, j = U ik, −j 1 ⎜ ⎜ ∆X i, j ⎝
⎞ ⎟ ⎟ ⎠
(3.6)
•
•
∂θ V ∂Y ∂ 2θ ∂Y 2
k −1 i, j
=V i, j
i, j
⎛ θ ik, j +1 − θ ik, j −1 ⎞ ⎜ ⎟ ⎜ ⎟ ∆ 2 Y ⎝ ⎠
(3.7)
⎛ θ ik, j +1 − 2θ ik, j + θ ik, j −1 ⎞ ⎟ =⎜ 2 ⎜ ⎟ ∆ ( Y ) ⎝ ⎠
(3.8)
Dengan menyusun ulang persamaan (3.6), (3.7) dan (3.8) identik dengan persamaan (3.5), diperoleh :
U
k −1 i, j
⎛ θ ik, j − θ ik−1, j ⎜ ⎜ ∆X ⎝
⎛ θ k − θ ik, j −1 ⎞ ⎛ θ k − 2θ ik, j + θ ik, j −1 ⎞ ⎞ ⎟ + Vi ,kj−1 ⎜ i , j +1 ⎟ = 1 ⎜ i , j +1 ⎟ 2 ⎜ ⎟ ⎟ Pr ⎜ ⎟ ∆ 2 Y ∆ ( Y ) ⎝ ⎝ ⎠ ⎠ ⎠
Variabel yang sudah diketahui disusun disebelah kanan tanda "=", dan variabel data yang belum diketahui diletakkan di sebelah kiri tanda "=". Diperoleh persamaan baru :
⎛ Vi ,kj−1 1 ⎜− − ⎜ 2∆Y Pr(∆Y ) 2 ⎝
⎛ U k −1 ⎞ k 2 ⎟θ i , j −1 + ⎜ i , j + 2 ⎜ ∆X ⎟ Pr (∆Y ) ⎝ ⎠ ⎛ Vi ,kj−1 1 ⎜ − ⎜ 2∆Y Pr (∆Y )2 ⎝
⎞ k ⎟θ i , j + ⎟ ⎠
⎛ U k −1 ⎞ ⎞ k ⎟θ i , j +1 = ⎜ i , j ⎟θ ik−1, j ⎜ ∆X ⎟ ⎟ ⎝ ⎠ ⎠
(3.9)
Koefisien matriks untuk persamaan di atas adalah :
⎛ Vi ,kj−1 1 aj = ⎜− − ⎜ 2∆Y Pr(∆Y ) 2 ⎝
⎞ ⎟ ⎟ ⎠
(3.10)
⎛ U ik, −j 1 2 + bj = ⎜ 2 ⎜ ∆X Pr (∆Y ) ⎝
⎞ ⎟ ⎟ ⎠
(3.11)
⎛ Vi ,kj−1 1 − cj = ⎜ ⎜ 2∆Y Pr (∆Y )2 ⎝
⎞ ⎟ ⎟ ⎠
(3.12)
⎛ U ik, −j 1 ⎞ k ⎟θ dj =⎜ ⎜ ∆X ⎟ i −1, j ⎝ ⎠
(3.13)
Persamaan (3.9) berubah menjadi : a jθ ik, j −1 + b jθ ik, j + c jθ ik, j +1 = d j
(3.14)
Persamaan (3.14) disebut persamaan diskritisasi energi. Dari persamaan (3.14) dapat dibuat Matriks Tridiagonal pada arah i, untuk j = 1,2, 3, 4, ..., ny ( j = 1 dan j = ny adalah kondisi batas),
θ i ,1 = θ w a 2θ i ,1 + b2θ i , 2 + c 2θ i ,3 = d 2 a3θ i , 2 + b3θ i ,3 + c3θ i , 4 = d 3 a 4θ i ,3 + b4θ i , 4 + c 4θ i ,5 = d 4 . . .
a ny −1θ i ,ny − 2 + bny −1θ i ,ny −1 + c ny −1θ i ,ny = d ny −1
θ i ,ny = 0 Persamaan-persamaan di atas dapat dituliskan dalam bentuk matriks berikut ini : ⎡1 ⎢a ⎢ 2 ⎢0 ⎢ ⎢0 ⎢. ⎢. ⎢ ⎢. ⎢0 ⎢ ⎣0
0
0
0
0
...
0
0
b2
c2
0
0
...
0
0
a3 0 .
b3 a4 .
c3 b4 .
0 c4 .
... ... .
0 0 .
0 0 .
. . 0 0
. . 0 0
. . 0 0
. . 0 0
. . 0 0
. .
. .
a ny −1 0
bny −1 0
0 ⎤ ⎡ θ i ,1 ⎤ ⎡ θ w ⎤ ⎢ ⎥ 0 ⎥⎥ ⎢ θ i , 2 ⎥ ⎢⎢ d 2 ⎥⎥ 0 ⎥ ⎢ θ i ,3 ⎥ ⎢ d 3 ⎥ ⎥ ⎢ ⎥⎢ ⎥ 0 ⎥ ⎢ θ i,4 ⎥ ⎢ d 4 ⎥ . ⎥⎢ . ⎥ = ⎢ . ⎥ ⎢ ⎥ . ⎥⎢ . ⎥ ⎢ . ⎥ ⎥ ⎢ ⎥ . ⎥⎢ . ⎥ ⎢ . ⎥ c ny −1 ⎥ ⎢θ i ,ny −1 ⎥ ⎢d ny −1 ⎥ ⎥ ⎢ ⎥⎢ ⎥ 1 ⎦ ⎢⎣ θ i ,ny ⎥⎦ ⎣ 0 ⎦
Matriks di atas disebut Matriks Tridiagonal untuk persamaan energi. 3.3.2 Diskritisasi Persamaan Momentum
Persamaan dasar momentum : u
∂ 2u ∂u ∂u = β g (T − T∞ ) + υ 2 +v ∂y ∂x ∂y
(3.15)
Dengan mensubstitusikan variabel tak berdimensi, tiap suku dari persamaan di atas dapat diubah kedalam bentuk persamaan tak berdimensi sebagai berikut : •
∂u ⎛ υ 2 .L.G ⎞ ∂U ⎟U =⎜ u ∂x ⎜⎝ W 4 ⎟⎠ ∂X
(3.16)
•
v
∂u ⎛ υ 2 .L.G ⎞ ∂U ⎟⎟V = ⎜⎜ 4 ∂y ⎝ W ⎠ ∂Y
(3.17)
•
β .g =
•
T − T∞ = θ (Tw − T∞ )
β β ref
⎛ υ 2 .L.G ⎜ ⎜ (T − T )W 4 ∞ ⎝ w
⎞ ⎟ ⎟ ⎠
(3.18) (3.19)
υ
•
∂ 2u ⎛ υ 2 .L.G ⎞ ∂ 2U ⎟ =⎜ ∂y 2 ⎜⎝ W 4 ⎟⎠ ∂Y 2
(3.20)
Substitusi persamaan (3.16), (3.17), (3.18), (3.19) dan (3.20) ke persamaan (3.15), diperoleh : ⎛ υ 2 .L.G ⎞ ∂U ⎛ υ 2 .L.G ⎞ ∂U β ⎟ ⎟⎟V ⎜⎜ ⎜⎜ + = U 4 4 ⎟ β ref ⎝ W ⎠ ∂X ⎝ W ⎠ ∂Y
⎛ υ 2 .L.G ⎜ ⎜ (T − T )W 4 ∞ ⎝ w
⎞ ⎟ θ (Tw − T∞ ) ⎟ ⎠
⎛ υ 2 .L.G ⎞ ∂ 2U ⎟⎟ 2 + ⎜⎜ 4 ⎝ W ⎠ ∂Y ⎛ υ 2 .L.G ⎞ ⎟⎟ menjadi : Persamaan di atas disederhanakan dengan mengeliminasi ⎜⎜ 4 W ⎠ ⎝
U
∂ 2U ∂U ∂U +V = θβ * + ∂X ∂Y ∂Y 2
(3.21)
Persamaan (3.21) disebut persamaan momentum dalam bentuk tak berdimensi. Dari persamaan (3.21) dapat dibuat koefisien Matriks tridiagonal berikut ini : ∂U U ∂X
•
∂U V ∂Y
•
∂ 2U ∂Y 2
•
⎛ U ik, j − U ik−1, j ⎜ ⎜ ∆X ⎝
k −1 i, j
⎛ U ik, j +1 − U ik, j −1 ⎞ ⎟ ⎜ ⎟ ⎜ 2 Y ∆ ⎠ ⎝
=U i, j
=V i, j
i, j
⎞ ⎟ ⎟ ⎠
k −1 i, j
(3.22)
(3.23)
⎛ U ik, j +1 − 2U ik, j + U ik, j −1 ⎞ ⎟ =⎜ 2 ⎟ ⎜ ( Y ) ∆ ⎠ ⎝
(3.24)
Dengan menyusun ulang persamaan (3.22), (3.23) dan (3.24) identik dengan persamaan (3.21), diperoleh : U
k −1 i, j
⎛ U ik, j − U ik−1, j ⎜ ⎜ ∆X ⎝
⎞ ⎛ U k − U ik, j −1 ⎞ ⎟ + Vi ,kj−1 ⎜ i , j +1 ⎟ = θ i, j β * + ⎟ ⎟ ⎜ 2 Y ∆ ⎠ ⎠ ⎝
⎛ U ik, j +1 − 2U ik, j + U ik, j −1 ⎞ ⎟ ⎜ 2 ⎟ ⎜ ( Y ) ∆ ⎠ ⎝
Dengan cara yang sama dengan persamaan energi, didapat : ⎛ Vi ,kj−1 1 ⎜− − ⎜ 2∆Y (∆Y )2 ⎝
⎞ k ⎛ U k −1 ⎟U i , j −1 + ⎜ i , j + 2 ⎟ ⎜ ∆X (∆Y ) 2 ⎠ ⎝
⎞ k ⎛ Vi ,kj−1 1 ⎞⎟ k ⎟U i , j + ⎜ − U = ⎜ 2∆Y (∆Y )2 ⎟ i , j +1 ⎟ ⎝ ⎠ ⎠
⎛ U ik, −j 1 ⎞ k ⎟U θ i, j β + ⎜ (3.25) ⎜ ∆X ⎟ i −1, j ⎠ ⎝ *
Koefisien matriksnya adalah : ⎛ Vi ,kj−1 1 aj = ⎜− − ⎜ 2∆Y (∆Y )2 ⎝
⎞ ⎟ ⎟ ⎠
(3.26)
⎛ U ik, −j 1 2 bj = ⎜ + ⎜ ∆X (∆Y ) 2 ⎝
⎞ ⎟ ⎟ ⎠
(3.27)
⎛ Vi ,kj−1 1 cj = ⎜ − ⎜ 2∆Y (∆Y )2 ⎝
⎞ ⎟ ⎟ ⎠
(3.28)
⎛ U ik, −j 1 ⎞ k ⎟U d j = θ i, j β + ⎜ ⎜ ∆X ⎟ i −1, j ⎠ ⎝ *
(3.29)
Sehingga persamaan (3.25) menjadi : a jU ik, j −1 + b jU ik, j + c jU ik, j +1 = d j
(3.30)
Persamaan (3.30) disebut persamaan diskritisasi momentum. Dari persamaan (3.30) dapat dibuat Matriks Tridiagonal pada arah i, untuk j = 1,2, 3, 4, ..., ny ( j = 1 dan j = ny merupakan kondisi batas), U i ,1 = 0 a 2U i ,1 + b2U i , 2 + c 2U i ,3 = d 2 a3U i , 2 + b3U i ,3 + c3U i , 4 = d 3 a 4U i ,3 + b4U i , 4 + c 4U i ,5 = d 4 . . .
a ny −1U i ,ny − 2 + bny −1U i ,ny −1 + c ny −1U i ,ny = d ny −1 U i,ny = U ∞ Persamaan-persamaan di atas dapat dituliskan dalam bentuk matriks sebagai berikut :
⎡1 ⎢a ⎢ 2 ⎢0 ⎢ ⎢0 ⎢. ⎢. ⎢ ⎢. ⎢0 ⎢ ⎣0
0
0
0
0
...
0
0
b2
c2
0
0
...
0
0
a3
b3
c3
0
...
0
0
0 .
a4 .
b4 .
c4 .
... .
0 .
0 .
.
.
.
.
.
.
.
. 0
. 0
. 0
. 0
. 0
.
.
a ny −1
bny −1
0
0
0
0
0
0
0
0 ⎤ ⎡ U i ,1 ⎤ ⎡ 0 ⎤ ⎥ ⎢ 0 ⎥⎥ ⎢ U i , 2 ⎥ ⎢⎢ d 2 ⎥⎥ 0 ⎥ ⎢ U i ,3 ⎥ ⎢ d 3 ⎥ ⎥ ⎢ ⎥ ⎥⎢ 0 ⎥ ⎢ U i,4 ⎥ ⎢ d 4 ⎥ . ⎥⎢ . ⎥ = ⎢ . ⎥ ⎥ ⎢ . ⎥⎢ . ⎥ ⎢ . ⎥ ⎥ ⎢ ⎥ . ⎥⎢ . ⎥ ⎢ . ⎥ c ny −1 ⎥ ⎢U i ,ny −1 ⎥ ⎢d ny −1 ⎥ ⎥ ⎢ ⎥ ⎥⎢ 1 ⎦ ⎣⎢ U i ,ny ⎦⎥ ⎣⎢ U ∞ ⎦⎥
Matriks Tridiagonal untuk persamaan momentum. 3.3.3 Diskritisasi Persamaan Kontinuitas
Untuk mendiskritisasi persamaan kontinuitas digunakan titik-titik nodal pada gambar 3.4. Persamaan kontinuitas di diskritisasi pada midpoint (bukan pada titik nodal), dinotasikan dengan (i,j-1/2), yang terletak pada baris ke-i dan berjarak setengah dari jarak antara j - 1 dengan j.
i, j-1
i, j-1/2
i, j
∆y ∆x
i-1, j-1
i-1, j
Gambar 3.4 Nodal untuk Diskritisasi Persamaan Kontinuitas Persamaan dasar kontinuitas :
∂u ∂v + =0 ∂x ∂y
(3.31)
Dengan mensubstitusikan variabel tak berdimensi, tiap suku dari persamaan di atas dapat diubah kedalam bentuk persamaan tak berdimensi sebagai berikut :
•
∂u = ∂x
⎛υ ⎜⎜ 2 ⎝W
⎞ ∂U ⎟⎟ ⎠ ∂X
(3.32)
•
∂v = ∂y
⎛υ ⎜⎜ 2 ⎝W
⎞ ∂V ⎟⎟ ⎠ ∂Y
(3.33)
Substitusi persamaan (3.32) dan (3.33) ke persamaan (3.31), diperoleh :
⎛υ ⎜⎜ 2 ⎝W
⎛υ ⎞ ∂U ⎟⎟ + ⎜⎜ 2 ⎝W ⎠ ∂X
⎞ ∂V ⎟⎟ =0 ⎠ ∂Y
⎛υ ⎞ Persamaan di atas disederhanakan dengan mengeliminasi ⎜⎜ 2 ⎟⎟ menjadi : ⎝W ⎠
∂U ∂V + =0 ∂X ∂Y
(3.34)
Persamaan (3.34) disebut persamaan kontinuitas dalam bentuk tak berdimensi. Dari persamaan (3.34) dapat dibuat koefisien Matriks tridiagonal berikut ini :
•
Diasumsikan derivatif x dititik (i,j-1/2) sama dengan rata-rata dari derivatif dititik (i,j) dan (i,j-1), yaitu :
∂U ∂X
= i , j −1 / 2
1 ⎡ ∂U ⎢ 2 ⎢⎣ ∂X
+ i, j
∂U ∂X
⎤ ⎥ i , j −1 ⎥ ⎦
dimana, ∂U ∂X
i, j
⎛ U i , j − U i −1, j = ⎜⎜ ∆X ⎝
⎞ ⎟⎟ ; ⎠
∂U ∂X
i , j −1
⎛ U i , j −1 − U i −1, j −1 ⎞ ⎟⎟ = ⎜⎜ ∆X ⎝ ⎠
Sehingga, ∂U ∂X
•
= i , j −1 / 2
1 ⎡U i , j − U i −1, j U i , j −1 − U i −1, j −1 ⎤ + ⎢ ⎥ 2⎣ ∆X ∆X ⎦
(3.35)
Derivatif y menggunakan pendekatan beda tengah, ∂V ∂Y
i , j −1 / 2
⎛ Vi , j − Vi , j −1 ⎞ ⎟⎟ = ⎜⎜ ∆Y ⎠ ⎝
(3.36)
Dengan menyusun ulang persamaan (3.35) dan (3.36) identik dengan persamaan (3.34), diperoleh : ⎛ Vi , j − Vi , j −1 ⎞ 1 ⎡U − U i −1, j U i , j −1 − U i −1, j −1 ⎤ ⎜⎜ ⎟⎟ = − ⎢ i , j + ⎥ ∆Y 2⎣ ∆X ∆X ⎝ ⎠ ⎦
⎛ ∆Y Vi , j = Vi , j −1 − ⎜ ⎝ 2∆X
⎞ ⎟(U i , j − U i −1, j + U i , j −1 − U i −1, j −1 ) ⎠
(3.37)
Persamaan (3.37) adalah persamaan kontinuitas hasil diskritisasi yang digunakan dalam penulisan program.
3.4 Kondisi Batas
Kondisi batas yang digunakan pada penelitian ini adalah daerah pada lapis batas aliran laminar konveksi alami pada plat datar vertikal seperti gambar (3.2) : .x
V = 0 U = 0 0 = 1
U 0
0 0
.y
.x = 0 y=0
U = 0 0 = 0
y
∞
Gambar 3.5 Kondisi Batas U =0
•
Untuk Y = 0,
V =0
θ =1 •
Untuk Y → ∞ ,
•
Untuk X = 0,
U →0
θ →0 U =0
θ =0
Kondisi batas pada Metode Beda Hingga digunakan untuk modifikasi koefisien matriks persamaan lapis batas, yaitu persamaan energi, momentum dan persamaan kontinuitas.
3.5 Perhitungan Koefisien Perpindahan Panas Konveksi Alami
Untuk menghitung koefisien perpindahan panas konveksi, digunakan persamaan (2.14) :
hx = −
⎛ ∂T ⎞ k ⎟ ⎜ (Tw − T∞ ) ⎜⎝ ∂y ⎟⎠ y →0
Dengan menggunakan beda mundur, derivatif dari persamaan diatas adalah :
Ti , 2 − Tw ⎛ ∂T ⎞ ⎟⎟ ⎜⎜ ≅ ∆y ⎝ ∂y ⎠ y →0
3.6 Penyusunan Algoritma dan Bagan Alir (Flow Chart) Program
Algoritma yang digunakan dalam penulisan program adalah sebagai berikut : a. Membaca data-data masukan (input); properti udara, temperatur, tebal lapis batas, syarat awal dan kondisi batas yang digunakan. b. Menghitung parameter tak berdimensi dan grid. c. Mengontrol angka Raleigh. Apabila Raleigh kurang dari 109, lanjutkan pengerjaan. Jika tidak, tulis ”Aliran Turbulen” dan sesuaikan masukan. d. Mengerjakan persamaan energi untuk menghitung temperatur, dengan menggunakan persamaan diskritisasi koefisien matriks (3.10) sampai dengan (3.13) e. Mengerjakan persamaan momentum untuk menghitung kecepatan U pada arah
i dengan menggunakan persamaan diskritisasi koefisien matriks (3.26) sampai (3.29). f. Mengerjakan persamaan kontinuitas menghitung kecepatan V arah j dengan menggunakan persamaan diskritisasi (3.37). g. Melakukan looping iterasi dengan memeriksa konvergensi, jika belum konvergen ulangi langkah c sampai d, jika sudah tulis data. h. Melakukan looping untuk i = 2 sampai dengan i = nx. i. Menulis data. j. Selesai. Diagram alir (Flow Chart) program : Mulai
Baca : Grav, Pr, w, vis, Tw, Tf, PL
A
Gambar 3.6 Diagram Alir Program
A
Menghitung : TA, Betar, Ra, Gr, G,Xmax
Ra < 109
TIDAK
YA Menghitung Grid :
∆ x, ∆ y Baca : Syarat Awal & Kondisi Batas U(1,1)=0, V(1,1)=0, T(1,1)=0 U(1,J)=0, V(1,J)=0, T(1,J)=0 U(i,1)=0, V(i,1)=0, T(i,1)=1, T(i,N)=0 VCHX=0
Mengerjakan Persamaan Energi
Panggil Subroutine Tridag
Mengerjakan Persamaan Momentum
Panggil Subroutine Tridag VCHX = V (i,NY) Mengerjakan Persamaan kontinuitas
Menghitung : VDIFF = ABS ( V(i,NY) – VCHX )
B
Gambar 3.6 (lanjutan)
B
VDIFF < 0.01
YA Menghitung : velu (i,j), velv (i,j), temp (i,j), hx (i)
Tulis : velu, velv, temp, hx
Selesai
Gambar 3.6 (lanjutan)
TIDAK
BAB IV DATA DAN ANALISIS
4.1 Validasi Program
Sebagai validasi program pada penelitian ini, digunakan penelitian yang dilakukan oleh Rolando A. Chavez untuk fluida superkritis. Kondisi batas yang digunakan adalah :
u=0 v=0 T = Tw
u=0 T = T∞
x y
u=0 T = T∞
Gambar 4.1 Kondisi batas penelitian Rolando A. Chavez Grid yang digunakan adalah grid dengan ∆x tidak seragam yang rapat di bagian bawah dan lebih renggang di bagian atas. Sedangkan grid ∆y konstan, seperti pada gambar 4.2. .x Y
.i = nx X
.y
.i = 1 j=1
j = ny
Gambar 4.2 Grid penelitian Rolando A. Chavez
4.2 Simulasi Konveksi Alami Plat Datar Vertikal Panas
Simulasi kasus konveksi alami pada plat datar vertikal ditampilkan dengan kondisi : a. Data ditentukan : - Angka Prandtl, Pr = 0.7 - Temperatur Plat, Tw = 50 oC - Temperatur udara, T∞ = 10 oC - Panjang plat, PL = 3 cm b. Data perhitungan : - Temperatur film, T f =
(Tw + T∞ ) 2
= 30 oC
- Viskositas kinematik pada temperatur film, υ = 16 x 10-6 (m2/s) - Konduktivitas termal pada Tf, k = 26.38 x 10-3 W/m oC - Koefisien ekspansivitas termal, β = - Angka Rayleigh, Ra =
1 Tf
β .g .(Tw − T∞ )(PL )3 Pr υ2
Kondisi batas yang digunakan adalah sebagai berikut : Untuk y = 0
: u = 0, v = 0 dan T = Tw
Untuk y = ∞ : u = 0 dan T = T∞ Untuk x = 0
: u = 0 dan T = T∞
Gambar 4.3 menunjukkan profil kecepatan dari penelitian Rolando A. chavez pada kasus konveksi alami plat datar vertikal dengan fluida air untuk ekspansivitas termal konstan dan ekspansivitas bervariasi. Nilai u diukur pada x = 0.02176 (m) dengan Pr = 1.05.
Gambar 4.3 Profil kecepatan hasil penelitian Rolando A. Chavez 0,40 0,35
(0,0012; 0,347)
u (m/s)
0,30 0,25 0,20 0,15 0,10 0,05 0,00 0,000
0,002
0,004
0,006
0,008
0,010
y (m)
Gambar 4.4 Profil kecepatan udara hasil penelitian Gambar 4.4 menunjukkan profil kecepatan udara pada plat datar vertikal konveksi alami. Nilai u di ukur pada x = 0.03 m dengan Pr = 0.7 dan Tf = 30 oC. Temperatur plat, Tw = 50 oC, temperatur udara T∞ = 10 oC dan Ra = 9.66 x 105 (angka Rayleigh total) yang dihitung dengan menggunakan persamaan 2.13 dengan panjang plat, PL = 0.03 m. Angka Rayleigh yang dihasilkan lebih kecil dari 109, menandakan bahwa aliran tersebut adalah laminar ( Ra < 109 ). Kecepatan fluida meningkat dari u = 0, di y = 0, hingga mencapai u maksimum pada u = 0.347 m/s di y = 1.2 x 10-3 m, kemudian secara bertahap turun hingga mencapai 0 pada lapis batas kecepatan. Secara kualitatif hasil yang diperoleh menunjukkan kesesuaian dengan hasil penelitian Rolando A. Chavez pada gambar 4.3.
Gambar 4.5 Profil temperatur hasil penelitian Rolando A. Chavez Gambar 4.5 menunjukkan profil temperatur dari penelitian Rolando A.
Temperatur (C)
Chavez yang dihitung dengan kondisi yang sama dengan gambar 4.3. 50 45 40 35 30 25 20 15 10 5 0 0,000
0,002
0,004
0,006
0,008
0,010
y (m)
Gambar 4.6 Profil temperatur udara hasil penelitian Gambar 4.6 menunjukkan profil temperatur udara pada kondisi yang sama dengan gambar 4.4. Semakin besar y, temperatur udara semakin mengecil dari 50 o
C pada plat sampai mencapai temperatur konstan 10 oC pada daerah lapis batas.
Hal ini terjadi karena pada daerah lapis batas sudah tidak terjadi transfer panas, dimana efek panas yang ditimbulkan plat sudah tidak ada. Sehingga temperaturnya menjadi sama dengan temperatur arus bebas (free stream), yaitu 10 o
C. Secara kualitatif, grafik gambar 4.6 menunjukkan kesesuian dengan hasil
penelitian Rolando A. Chavez pada gambar 4.5.
Gambar 4.7 Distribusi kecepatan hasil penelitian Rolando A. Chavez Gambar 4.7 menunjukkan plot kontur kecepatan dari penelitian Rolando A. Chavez untuk fluida superkritis air dengan Pr = 1.05. VELU 0.03
0.3 0.025
0.25 0.02
x (m)
0.2 0.015 0.15
0.01 0.1
0.005
0
0.05
0
1
2
3
4
5 y (m)
6
7
8
9 -3
x 10
Gambar 4.8 Distribusi kecepatan udara hasil penelitian Gambar 4.8 menunjukkan kontur kecepatan udara yang di visualisasikan dengan perangkat lunak Matlab R13. Kecepatan diplot dengan domain 0.01 m sepanjang plat 0.03 m. Warna merah tua pada kontur menunjukkan nilai tertinggi dari kecepatan dan warna biru tua menunjukkan besar kecepatan sama dengan nol. Di daerah sekitar plat, kecepatan fluida sama dengan nol. Kemudian naik hingga mencapai suatu titik maksimum dan turun kembali mencapai nol pada batas domain.
Hasil yang didapat menunjukkan kesesuian dengan teori yang ada. Dimana pada dinding, kecepatannya adalah nol karena terdapat kondisi tanpa gelincir (no-
slip condition). Kecepatan bertambah sampai mencapai suatu nilai maksimum, kemudian turun secara bertahap mencapai nol pada tepi lapis batas (Ozisik, 1988). Secara kualitatif kontur gambar 4.8 menunjukkan kesesuaian dengan hasil penelitian Rolando A. Chavez pada gambar 4.7.
Gambar 4.9 Distribusi temperatur hasil penelitian Rolando A. Chavez Gambar 4.9 menunjukkan plot kontur temperatur dari penelitian Rolando A. Chavez untuk fluida superkritis air dengan Pr = 1.05. TEMPERATUR 0.03
45 0.025 40 0.02
x (m)
35
30
0.015
25 0.01 20 0.005 15
0
0
1
2
3
4
5 y (m)
6
7
8
10
9 -3
x 10
Gambar 4.10 Distribusi temperatur udara hasil penelitian
Gambar 4.10 adalah kontur temperatur udara yang diplot dengan kondisi yang sama dengan gambar 4.8. Warna merah tua menunjukkan temperatur udara tertinggi dan warna biru tua menunjukkan temperatur udara terendah. Temperatur udara yang bersentuhan dengan plat adalah 50 oC. Semakin jauh dari plat temperatur udara turun sampai mencapai 10 oC. Temperatur ini sama dengan temperatur arus bebas (free stream). Hasil ini menunjukkan kesesuaian dengan teori yang ada bahwa pada aliran konveksi alami, temperatur udara dekat plat adalah yang tertinggi dan secara bertahap turun hingga memiliki temperatur sama dengan temperatur arus bebas dimana efek panas plat sudah tidak berpengaruh. Secara kualitatif kontur temperatur udara gambar 4.10 menunjukkan kesesuaian dengan hasil penelitian Rolando A. Chavez gambar 4.9. 140
hx (W/m2.C)
120 100 80 60 40 20 0 0,000
0,005
0,010
0,015
0,020
0,025
0,030
0,035
x (m)
Gambar 4.11 Distribusi koefisien perpindahan panas konveksi alami lokal pada Pr = 0.7 dan Tf = 30 oC Gambar 4.11 menunjukkan distribusi koefisien perpindahan panas konveksi alami lokal pada plat. Grafik koefisien perpindahan panas konveksi alami lokal menurun atau berbanding terbalik dengan jarak titik (x) pada plat. Semakin jauh jarak dari ujung plat, semakin kecil harga koefisien perpindahan panas konveksi lokalnya. Hal ini menunjukkan kesesuaian dengan teori yang ada dan dapat dikoreksi dengan persamaan (2.14) :
hx = −
⎛ ∂T ⎞ k ⎟ ⎜ (Tw − T∞ ) ⎜⎝ ∂y ⎟⎠ y →0
Kondisi pada lapis batas termal sangat dipengaruhi oleh gradien temperatur, ∂T / ∂y y →0 . Karena (Tw – T∞) konstan tidak terpengaruh oleh x dan δt
meningkat dengan bertambahnya x, maka gradien temperatur mengecil jika x bertambah. Dengan mengecilnya gradien temperatur
∂T / ∂y
y →0
ketika x
bertambah, mengakibatkan nilai hx mengecil. Atau dengan kata lain koefisien perpindahan panas konveksi alami lokal berbanding terbalik dengan jarak x dari ujung plat.
BAB V PENUTUP
5.1 Kesimpulan
Dari penelitian dan pembahasan yang telah dilakukan dapat ditarik beberapa kesimpulan, yaitu : a. Kode program yang dirancang dalam penelitian ini dapat bekerja dengan baik sesuai dengan tujuan awal penelitian untuk membuat simulasi perpindahan panas konveksi alami pada plat datar vertikal panas. b. Kecepatan udara yang diukur pada x = 0.03 m dengan Pr = 0.7 dan Tf = 30 oC, mengalami peningkatan hingga mencapai titik maksimum di y = 1.2 x 10-3 m dengan nilai u = 0.347 m/s kemudian turun secara bertahap sampai u = 0 pada daerah lapis batas. c. Temperatur udara turun dari 50 oC, yang bersinggungan dengan plat, hingga mencapai 10 oC, yaitu sama dengan temperatur arus bebas, di daerah lapis batas. d. Koefisien perpindahan panas konveksi alami lokal berbanding terbalik dengan jarak x dari ujung plat. Semakin jauh jarak pada plat, semakin kecil nilai koefisien perpindahan panas konveksi lokalnya. e. Secara kualitatif, profil kecepatan udara dan profil temperatur udara yang dihasilkan dari penelitian ini menunjukkan hasil yang sama dengan penelitian pada fluida superkritis yang dilakukan oleh Rolando A. Chavez. f. Secara kualitatif, plot distribusi kecepatan udara dan plot distribusi temperatur udara yang dihasilkan dari penelitian ini menunjukkan hasil yang sama dengan hasil penelitian pada fluida superkritis yang dilakukan oleh Rolando A. Chavez g. Secara kualitatif, hasil perhitungan koefisien perpindahan panas konveksi
alami lokal menunjukkan kesesuaian dengan teori perpindahan panas.
5.2 Saran
Untuk lebih mengembangkan ilmu komputasi perpindahan panas dan simulasi numerik, penulis memberikan saran untuk : a. Melakukan pengembangan penelitian lebih lanjut dengan perhitungan secara kuantitatif. b. Melakukan penelitian kasus perpindahan panas konveksi alami pada plat datar vertikal panas dengan metode yang berbeda.
DAFTAR PUSTAKA
Anderson, J.D. 1995. Computational Fluid Dynamics The Basics With Applications. Singapore: McGraw-Hill, Inc. Bejan, Adrian. 1993. Heat transfer. Singapore: John Wiley & Sons, Inc. Chavez, R.A.C. 2004. Natural – Convection Heat Transfer in Supercritical Fluids. Perto Rico: Mechanical Engineering Dept. University of Puerto Rico. Fox, R. and McDonald, A. 1991. Introduction to Fluid Mechanics. El Hadidi, B.M. 1998. A Computational Study of Flow In Mechanically Ventilated Space. Egypt: Cairo University. Hoffmann, K.A. 1989. Computational Fluid Dynamics for Engineers. Austin, Texas: A Publication of Engineerng Education System. Holman, J.P. 1988. Perpindahan Kalor. Jakarta: Erlangga. Incropera, F.M. 1996. Introduction to Heat Transfer. USA: John Wiley & Sons. Lemos, C.M. 1993. FDFlow : A Fortran-77 Solver for 2-D Incompressible Fluid Flow. Computers & Geosciences, Vol. 20, No.3, pp. 265-291. Oosthuizen, PH. 1999. An Introduction to Convective Heat Transfer Analysis. Queen's University. USA: WCB/McGraw-Hill Book Company. Ozisik, M Necati. 1988. Elements of Heat Transfer. McGraw-Hill Book Company.
Lampiran 1. Program dengan FORTRAN PS 4.0 ************************************************************ ****** PROGRAM ****** ****** SIMULASI NUMERIK PERPINDAHAN PANAS ****** ****** KONVEKSI ALAMI PLAT DATAR VERTIKAL ****** ************************************************************ ************************************************************ ************ KETERANGAN SIMBOL ****************** ************************************************************ *** a = koefisien matriks *** b = koefisien matriks *** BETAR = koefisien ekspansivitas termal *** c = koefisien matriks *** d = koefisien matriks *** dx = jarak grid arah x *** dxmax = jarak grid maksimum arah x *** dy = jarak grid arah y *** G = modifikasi angka Grashof *** Gr = angka grashof *** grav = percepatan gravitasi *** hx = koefisien perpindahan panas lokal *** iter = iterasi *** k = konduktivitas termal udara *** m = indeks baris *** n = indeks kolom *** nump = 100 *** nx = jumlah grid arah x *** ny = jumlah grid arah y *** PL = panjang plat *** Pr = angka Prandtl *** Ra = angka Rayleigh *** rex = under relaxation *** sum = penjumlahan *** T = temperatur non-dimensional *** TA = temperatur rata-rata *** temp = temperatur dimensional *** TF = temperatur fluida *** TW = temperatur dinding (wall) *** U = kecepatan non-dimensional arah x *** V = kecepatan non-dimensional arah y *** VCHX = nilai V kontrol *** vdiff = selisih harga V *** velu = kecepatan dimensional arah x *** velv = kecepatan dimensional arah y *** vis = viskositas kinematik udara *** W = tebal domain *** x = koordinat non-dimensional arah x *** Xmax = panjang domain non-dimensional *** xx = koordinat berdimensi arah x sejajar plat *** y = koordinat non-dimensional arah y *** yy = koordinat berdimensi arah y normal terhadap plat ************************************************************ ************************************************************
Lampiran 1. (sambungan) ************************************************************ parameter(m=500,n=500) dimension U(m,n),V(m,n),T(m,n),a(m),b(m) dimension c(m),d(m),X(m),Y(m) dimension xx(m),yy(n) dimension velu(m,n),velv(m,n),temp(m,n) dimension hx(m) real k ************************************************************ open(2,file='c:\matlab6p5p1\work\temp') open(3,file='c:\matlab6p5p1\work\velu') open(4,file='c:\matlab6p5p1\work\num') open(7,file='c:\matlab6p5p1\work\hx') ************************************************************ grav = 9.81 W = 0.02 Pr = 0.7 VIS = 16.E-6 ******************
DATA SUHU
******************
TW = 50. TF = 10. TA = (1./2)*(TW+TF) BETAR = 1./TA ******************
PANJANG PLAT
******************
PL = 0.03 write(*,*)' Panjang Plat= ',PL *************
PARAMETER TAK BERDIMENSI
Ra = (BETAR*grav*(tw-tf)*pr*pl**3)/vis**2 write(*,*)' Raleigh = ',Ra if(Ra.GT.10E9)then write(*,*)’ Aliran Turbulen’ stop endif Gr=Ra/Pr G=BETAR*grav*(TW-TF)*(W**4)/((VIS**2)*PL) Xmax = 1./G write(*,*)' Grashof = ',Gr write(*,*)' Xmax = ',xmax write(*,*)' G = ',G x(1)=0.0 sum=0 nump=100 do it=1,nump sum=sum+1.05**it enddo
*************
Lampiran 1. (sambungan) dx=0.75*Xmax/sum dxmax=dx*(1.05**nump) do i=2,500 if(dx.lt.dxmax)then dx=1.05*dx else dx=dxmax endif x(i)=x(i-1)+dx nx=i if(x(i).gt.xmax)goto 200 enddo 200
ny=301 y(1)=0.0 dy=1./(ny-1) rex=0.5 do j=2,ny y(j)=y(j-1)+dy enddo write(*,*)' x(nx) =',x(i) write(*,*)' nx =',nx
***********
SYARAT AWAL & KONDISI BATAS
U(1,1) = 0.0 V(1,1) = 0.0 T(1,1) = 1.0 do j= 2,ny U(1,J) = 0.0 T(1,J) = 0.0 V(1,J) = 0.0 enddo do i=2,nx V(i,1) = 0.0 U(i,1) = 0.0 U(i,ny) = 0.0 T(i,ny) = 0.0 T(i,1) = 1.0 enddo do i=2,nx dx=x(i)-x(i-1) ITER=0 VCHX=0.0 500
ITER=ITER+1
**********
Lampiran 1. (sambungan) ******
MENGERJAKAN PERSAMAAN ENERGI MENGHITUNG "T"
*****
a(1)=0.0 b(1)=1.0 c(1)=0.0 d(1)=t(i,1) do j=2,ny-1 a(j)=(-v(i,j)/2/dy) - (1./pr/dy/dy) b(j)=(u(i,j)/dx) + (2./pr/dy/dy) c(j)=(v(i,j)/2/dy) - (1./pr/dy/dy) d(j)=u(i,j)*t(i-1,j)/dx enddo a(ny)=0.0 b(ny)=1.0 c(ny)=0.0 d(ny)=t(i,ny) call tridag(a,b,c,d,1,ny) do j=1,ny t(i,j)=t(i,j)+rex*(d(j)-t(i,j)) enddo
*******
MENGERJAKAN PERSAMAAN MOMENTUM MENGHITUNG "U"
*****
a(1)=0.0 b(1)=1.0 c(1)=0.0 d(1)=u(1,j) do j=2,ny-1 a(j)=(-v(i,j)/2/dy) - (1./dy/dy) b(j)=(u(i,j)/dx) + (2./dy/dy) c(j)=(v(i,j)/2/dy) - (1./dy/dy) d(j)=t(i,j) + u(i,j)*u(i-1,j)/dx enddo a(ny)=0.0 b(ny)=1.0 c(ny)=0.0 d(ny)=u(ny,j) call tridag(a,b,c,d,1,ny) do j=1,ny u(i,j)=u(i,j)+rex*(d(j)-u(i,j)) enddo
*****
MENGERJAKAN PERSAMAAN KONTINUITAS MENGHITUNG "V"
do j=2,ny v(i,j) = v(i,j-1)-(dy/(2.0*dx))*(u(i,j)-u(i-1,j)+ c u(i,j-1)-u(i-1,j-1)) enddo
*****
Lampiran 1. (sambungan) vdiff=v(i,ny)-vchx if(iter.lt.50) goto 500 if(iter.gt.100) goto 101
300
if(vdiff.lt.0.01) goto 300 vchx=v(i,ny) goto 500 CONTINUE write(*,*)' i= ',i,' x = ',X(i),' iter= ',iter IF(X(i).GE.XMAX) goto 102 enddo
101
102
write(6,*)' ny>100' goto 103 write(6,1000) write(*,*)' ny = ',ny write(*,*)' nx = ',i write(4,*)ny write(4,*)nx
*************
PARAMETER BERDIMENSI
***************
do i=1,nx xx(i)=x(i)*PL*G enddo do j=1,ny yy(j)=y(j)*w enddo do j=1,ny do i=1,nx velu(i,j)=u(i,j)*vis*PL*G/(w**2) velv(i,j)=v(i,j)*vis/w temp(i,j)=t(i,j)*(Tw-Tf)+Tf enddo enddo
************
KOEFISIEN PERPAN KONVEKSI
k=26.38E-3 do i=2,nx hx(i)=-k*(temp(i,2)-Tw)/(dy*w*(Tw-Tf)) write(7,*)xx(i),hx(i) enddo do j=1,ny write(3,*)yy(j),velu(100,j) enddo
*************
Lampiran 1. (sambungan) do i=1,nx do j=1,ny write(2,*)yy(j),xx(i),velu(i,j) enddo enddo 103 1000
CONTINUE FORMAT (/,1X,' X > Xmax') STOP END
1
2
subroutine tridag(a,b,c,d,l1,l2) Parameter(m=500) dimension a(m),b(m),c(m),d(m) do 1 i=l1+1,l2 r=-a(i)/b(i-1) b(i)=b(i)+r*c(i-1) d(i)=d(i)+r*d(i-1) d(l2)=d(l2)/b(l2) do 2 j=l2-1,l1,-1 d(j)=(d(j)-c(j)*d(j+1))/b(j) return end
Lampiran 2. Program Visualisasi dengan Matlab R13 Program Matlab untuk Plot 1 Dimensi
load VELU -ascii; im=301; jm=1; imax=im*jm; for j=1:jm for i=1:im is=i+(j-1)*im; x(i,j)=VELU(is); y(i,j)=VELU(is+imax); end end hold on; box on; plot(x,y); clear; hold off;
Lampiran 2. (sambungan)
Program Matlab untuk Plot 2 Dimensi
load temp -ascii; load num -ascii; im=num(1); jm=num(2); imax=im*jm; for j=1:jm for i=1:im is=i+(j-1)*im; x(i,j)=temp(is); y(i,j)=temp(is+imax); o(i,j)=temp(is+2*imax); end end hold on; box on; contourf(x,y,o,100); contour(x,y,o,100); colorbar; clear hold off;
Lampiran 3. Simulasi dengan Fluida Air Air
PL
= 0.02 m
Tw
= 40 oC
T∞
= 20 oC
Tf
= 30 oC
Pr
= 5.42
υ
= 0.8012 x 10-6 m2/s
k
= 0.6150 W/m.C
W
= 0.004 m
48
VELU 0.02
0.1
0.018
0.09
0.016
0.08
x (m)
0.014
0.07
0.012
0.06
0.01
0.05
0.008
0.04
0.006
0.03
0.004
0.02
0.002
0.01
0
0
0.5
1
1.5
2 y (m)
2.5
3
0
3.5 -3
x 10
Gambar L3.1 Plot distribusi kecepatan air pada x = 1,14 x 10-3 m
0,030 0,025
(0,000107; 0,0245)
u (m/s)
0,020 0,015 0,010 0,005 0,000 0,0000 0,0005 0,0010 0,0015 0,0020 0,0025 0,0030 0,0035 0,0040 y (m)
Gambar L3.2 Profil kecepatan air Kecepatan maksimum terjadi pada jarak y = 1.07 x 10-4 dengan u = 0.0245 m/s
Lampiran 3. (sambungan)
49
TEMPERATUR
0.02
38
0.016
36
0.014
34
0.012
32
0.01
30
0.008
28
0.006
26
0.004
24
0.002
22
x (m)
0.018
0
0
0.5
1
1.5
2 y (m)
2.5
3
20
3.5 -3
x 10
Gambar L3.3 Plot distribusi temperatur air pada x = 1,14 x 10-3 m 40
Temperatur (C)
35 30 25 20 15 10 5 0 0,0000 0,0005 0,0010 0,0015 0,0020 0,0025 0,0030 0,0035 0,0040 y (m)
Gambar L3.4 Profil temperatur air 25000
hx (W/m2.C)
20000 15000 10000 5000 0 0,000
0,005
0,010
0,015
0,020
x (m)
Gambar L3.5 Distribusi koefisien perpindahan panas lokal
50
Lampiran 4. Simulasi dengan Fluida Udara Udara
PL
= 0.1 m
Tw
= 60 oC
T∞
= 20 oC
Tf
= 40 oC
Pr
= 0.71
υ
= 16.96 x 10-6 m2/s
k
= 27.10 x 10-3 W/m.oC
W
= 0.02 m
51
VELU
0.1
0.5
0.09
0.45 0.08 0.4 0.07 0.35 0.06 x (m)
0.3 0.05 0.25 0.04 0.2 0.03
0.15
0.02
0.1
0.01
0
0.05
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
y (m)
u (m/s)
Gambar L4.1 Plot distribusi kecepatan udara pada x = 0.071 m
0,50 0,45 0,40 0,35 0,30 0,25 0,20 0,15 0,10 0,05 0,00 0,00
(0,00167; 0,461)
0,01
0,01
0,02
0,02
y (m)
Gambar L4.2 Profil kecepatan udara Kecepatan maksimum terjadi pada jarak y = 0.00167 m dengan u = 0.461 m/s. Lampiran 4. (sambungan)
52
TEMPERATUR 0.1
0.09
55
0.08 50 0.07 45
x (m)
0.06 40
0.05
0.04
35
0.03 30 0.02 25 0.01
0
0
0.002
0.004
0.006
0.008
0.01 y (m)
0.012
0.014
0.016
0.018
20
Gambar L4.3 Plot distribusi temperatur udara pada x = 0.071 m 60
Temperatur (C)
50 40 30 20 10 0 0,00
0,01
0,01
0,02
y (m)
Gambar L4.4 Profil temperatur udara
Lampiran 4. (sambungan)
0,02
53
90 80 hx (W/m2.C)
70 60 50 40 30 20 10 0 0,00
0,02
0,04
0,06
0,08
0,10
x (m)
Gambar L4.5 Distribusi koefisien perpindahan panas lokal
54 Lampiran 5. Simulasi dengan fluida Karbon Dioksida (CO2) Karbon Dioksida (CO2)
PL
= 0.0161 m
Tw
= 67.7 oC
T∞
= 25.12 oC
Tf
= 46.41 oC
Pr
= 2.02
υ
= 2.716 x 10-7 m2/s
k
= 6.596 x 10-5 W/m.oC
W
= 0.001 m
55
VELU 0.016
0.15
0.014
0.012 0.1
x (m)
0.01
0.008
0.006 0.05 0.004
0.002
0
0
1
2
3
4
5 y (m)
6
7
8
0
9 -4
x 10
Gambar L5.1 Plot distribusi kecepatan CO2 pada x = 1.413 x 10-4 m 0,016 0,014
(4,0E-05; 0,0139)
0,012 x (m)
0,010 0,008 0,006 0,004 0,002 0,000 0,0000
0,0002
0,0004
0,0006
0,0008
0,0010
y (m)
Gambar L5.2 Profil kecepatan CO2 Kecepatan maksimum terjadi pada jarak y = 4 x 10-5 dengan u = 0.0139 m/s
Lampiran 5. (sambungan)
56
TEMPERATUR 0.016 65 0.014 60 0.012 55 0.01
x (m)
50 0.008
45
0.006
40
0.004
35
0.002
30
0
0
1
2
3
4
5 y (m)
6
7
8
9 -4
x 10
Gambar L5.3 Plot distribusi temperatur CO2 pada x = 1.413 x 10-4 m 70
Temperatur (C)
60 50 40 30 20 10 0 0,0000
0,0002
0,0004
0,0006
0,0008
0,0010
y (m)
Gambar L5.4 Profil temperatur CO2 4 3,5 hx (W/m2.C)
3 2,5 2 1,5 1 0,5 0 0,000
0,005
0,010
0,015
0,020
x (m)
Gambar L5.5 Distribusi koefisien perpindahan panas lokal
