SIMULASI NUMERIK KONVEKSI ALAMI PADA SINGLE FIN DAN MULTIPLE FINS DALAM KOTAK 2D DENGAN METODE BEDA HINGGA

SIMULASI NUMERIK KONVEKSI ALAMI PADA SINGLE FIN DAN MULTIPLE FINS DALAM KOTAK 2D DENGAN METODE BEDA HINGGA

SKRIPSI

Diajukan sebagai salah satu syarat untuk memperoleh gelar Sarjana Teknik

Oleh : FERRY ENDHARTA NIM. I1407514

JURUSAN TEKNIK MESIN FAKULTAS TEKNIK UNIVERSITAS SEBELAS MARET SURAKARTA 2010

SIMULASI NUMERIK KONVEKSI ALAMI PADA SINGLE FIN DAN MULTIPLE FINS DALAM KOTAK 2D DENGAN METODE BEDA HINGGA Di susun oleh : Ferry Endharta NIM. I1407514

Dosen Pembimbing I

Dosen Pembimbing II

Eko Prasetya Budiana, S.T.M.T. NIP. 197109261999031002

Budi Kristiawan, S.T.M.T. NIP. 197104251999031001

Telah dipertahankan di hadapan dosen tim penguji pada hari Selasa, tanggal 13 Juli 2010.

1. Syamsul Hadi, S.T.M.T. NIP. 197106151998021002

…………………………..

2. Purwadi Joko Widodo, S.T.M.Kom NIP. 197301261997021001

…………………………..

3. Rendy Adhi Rachmanto, S.T.M.T. NIP. 197101192000121006

…………………………..

Mengetahui,

Ketua Jurusan Teknik Mesin

Koordinator Tugas Akhir

Dody Ariawan, S.T.M.T. NIP. 197308041999031003

Wahyu Purwo Raharjo, S.T.M.T. NIP. 197202292000121001

ii

MOTTO DAN PERSEMBAHAN Motto

Capailah cita–citamu dengan usaha maksimal niscaya akan berhasil. Percayalah Allah akan selalu ada untuk kita, sesulit apapun masalah yang kita hadapi pasti ada jaln untuk kita Berbuatlah yang terbaik untuk setiap detiknya Dan tak perlu kau pikirkan esok akan jadi apa Bila keadaan sudah berada dititik buruk yang paling buruk, rendah yang paling rendah dan sakit yang paling sakit, maka tidak lain yang akan terjadi adalah menjadi baik.

Persembahan karya ini kupersembahkan untuk : Bapak, Ibu dan adikku tercinta

iii

ABSTRAK

FERRY ENDHARTA, Simulasi Numerik Konveksi Alami Pada Single Fin Dan Multiple Fins Dalam Kotak 2D Dengan Metode Beda Hingga

Simulasi numerik untuk konveksi alami dalam kotak 2D pada single fin dan multiple fins dilakukan untuk mengetahui fenomena yang terjadi pada pola aliran dan distribusi temperatur di sekitar fin pada kotak 2D. Tulisan ini menguraikan metode untuk penyelesaiaan konveksi alami steady dalam kotak 2D. Metode ini didasarkan pada skema Runge – Kutta untuk diskritisasi waktu dan skema kompak beda hingga orde-4 untuk diskritisasi ruang. Kesulitan pada penyelesaian tekanan diselesaikan dengan metode kompresibilitas tiruan. Metode beda hingga dituliskan dengan bahasa Fortran dan divisualisasikan dengan perangkat lunak Matlab. Hasil dari visualisasi menunjukkan bahwa separasi muncul karena adanya pusaran kecepatan disekitar sirip. Semakin besar angka Rayleigh dengan domain dan panjang sirip yang sama, kerapatan pepindahan panas disekitar dinding akan meningkat. Hasil menunjukkan pada angka Rayleigh 106 dan 107 telah terjadi separasi disekitar sirip. Kata kunci : konveksi alami, skema kompak, kompresibilitas tiruan, separasi.

iv

ABSTRACT

FERRY ENDHARTA, Numerical Simulation of natural convection with single fin and multiple fins in 2D cavity by finite different method

Numerical simulation for natural convection in 2D cavity with single fin and multiple fins done to know the phenomenon that happened at stream function and distribution of temperature around fin in 2D cavity. This paper present numerical method for solving steady natural convection in 2D cavity. The method is based on Runge-Kutta schemes for temporal discretization and fourh-order compact finite difference schemes for spatial discretization. Difficulty related to the pressure can be overcome by using artificial compressibility method. Finite difference written by Fortran language and visualized by Matlab. Result shown that separation appear caused by vortex in the velocity vector around fin. Increasing of Rayleigh number with same fin length and domain, density of heat transfer around wall will be increase. Result show for Rayleigh number 106 and 107 have happened separation around fin.

Key words : natural convection, compact schemes, artificial compressibility, separation

v

KATA PENGANTAR 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 Konveksi Alami Pada Single Fin Dan Multiple Fins Dalam Kotak 2D 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 yang sebesar – besarnya kepada semua pihak yang telah memberikan bantuan, doa, dukungan dan semangat, baik moral maupun spiritual kepada : 1. Bapak Dody Ariawan, ST., MT., selaku Ketua Jurusan Teknik Mesin UNS. 2. Bapak Eko Prasetya Budiana, ST.,MT., selaku Pembimbing I tugas akhir, atas bimbingan, nasehat, kesabaran, dan ilmu pengetahuan yang diajarkannya. 3. Bapak Budi Kristiawan, ST.,MT., selaku Pembimbing II tugas akhir , atas bimbingan , nasehat, kesabaran dan ilmu pengetahuan yang diajarkannya. 4. Bapak

Purwadi

Joko

Widodo,

S.T.M.Kom,

selaku

Pembimbing

Akademik. 5. Bapak – bapak dosen dan staf karyawan di lingkungan Teknik Mesin UNS, atas didikan, nasehat, ilmu yang diajarkan dan kerjasamanya. 6. Teman – teman Teknik Mesin transfer angkatan 2007 Teman – teman Teknik Mesin UNS 7. Teman – teman kos Oblong. 8. Dan semua pihak yang telah membantu Penulis dalam menyelesaikan skripsi ini.

Penulis menyadari bahwa skripsi ini masih banyak terdapat kekurangan, untuk itu masukan dan saran yang membangun akan penulis terima dengan ikhlas

vi

dan penulis mengucapkan terima kasih. Penulis berharap semoga skripsi ini dapat memberikan manfaat bagi penulis khususnya dan bagi pembaca pada umumnya.

Surakarta, juli 2010

Penulis

vii

DAFTAR ISI

HALAMAN JUDUL........................................................................................ i HALAMAN PENGESAHAN ……………………………………………..... ii MOTTO DAN PERSEMBAHAN……………………………………….. ..... iii ABSTRAK…………………………………………………………………. .. iv ABSTRACT………………………………………………………………… ... v KATA PENGANTAR…………………………………………………… ..... vi DAFTAR ISI………………………………………………………………. ... viii DAFTAR GAMBAR……………………………………………………… ... x DAFTAR TABEL ............................................................................................ xii DAFTAR NOTASI.................................………………………………….. ... xiii DAFTAR LAMPIRAN……. ........................................................................... xv BAB I

PENDAHULUAN ............................................................................ 1 1.1 Latar Belakang Masalah .............................................................. 1 1.2 Perumusan Masalah .................................................................... 2 1.3 Batasan Masalah.......................................................................... 2 1.4 Tujuan Penelitian ........................................................................ 2 1.5 Manfaat Penelitian ...................................................................... 2 1.6 Sistematika Penulisan ................................................................. 3

BAB II LANDASAN TEORI ....................................................................... 4 2.1 Tinjauan Pustaka ......................................................................... 4 2.2 Dasar Teori…………………………………………………. ..... 6 2.2.1 Sirip…………………… .................................................... 6 2.2.2 Konveksi Alami………………………………………… . 7 2.2.2.1 Persamaan Atur Konveksi Alami ........................... 8 2.2.2.2 Diskritisasi Waktu .................................................. 8 2.2.2.3 Diskritisasi Ruang .................................................. 9 2.2.2.4 Metode Kompresibilitas Tiruan ............................. 13 BAB III PELAKSANAAN PENELITIAN ................................................ 15

viii

3.1 Alat dan Bahan ....................................................................... 15 3.2 Garis Besar Penelitian ............................................................ 15 3.3 Diskritisasi Persamaan Atur ................................................... 17 3.3.1 Diskritisasi Persamaan Momentum ............................... 17 3.3.2 Diskritisasi Persamaan Energi ....................................... 18 3.3.3 Diskritisasi Metode Kompresibilitas Tiruan ................. 19 3.4 Diskritisasi Syarat Batas ......................................................... 19 3.5 Algoritma Pemrograman ........................................................ 24 BAB IV HASIL DAN PEMBAHASAN 4.1 Simulasi Konveksi Alami Sirip Tunggal dan Multiple Sirip . 24 4.2 Validasi Program .................................................................... 32 BAB V PENUTUP ..................................................................................... 37 5.1 Kesimpulan ............................................................................................ 37 5.2 Saran ...................................................................................................... 37 DAFTAR PUSTAKA.................................................................................. 38 LAMPIRAN ................................................................................................ 39

ix

DAFTAR GAMBAR

Gambar 2.1

stream function dan isotherm pada Ra = 106 ............................ 4

Gambar 2.2

separasi di sekitar sirip ............................................................. 5

Gambar 2.3

Beberapa contoh jenis extended surface: (a) sirip longitudinal (memanjang) dengan profil segiempat (b)pipa silindris dengan sirip berprofil segiempat (c) sirip longitudinal dengan profil trapezioda (d) sirip longitudinal dengan profil parabola (e) pipa silindris dengan sirip radial berprofil segiempat (f) pipa silindris dengan sirip radial berprofil trapezoida (g)cylindrical spine (h)truncated conical spine (i) truncated parabolic spine ......................................................................................... 6

Gambar 2.4

Galat dispersi untuk pendekatan numerik dari turunan Pertama ..................................................................................... 10

Gambar 2.5

Galat disipasi untuk pendekatan numerik dari turunan kedua ........................................................................................ 12

Gambar 3.1

Domain dan syarat batas ........................................................... 18

Gambar 4.1

Isothermal pada ra = 106 pada sirip tunggal ............................. 24

Gambar 4.2

Isothermal pada Ra = 106 pada 2 sirip ...................................... 24

Gambar 4.3

Isothermal pada Ra = 106 pada 3 sirip ...................................... 25

Gambar 4.4

Isothermal pada Ra = 106 pada 3 sirip ...................................... 25

Gambar 4.5

Isothermal pada Ra = 107 pada sirip tunggal ............................ 25

Gambar 4.6

Isothermal pada Ra = 107 pada 2 sirip ...................................... 26

Gambar 4.7

Isothermal pada Ra = 107 pada 3 sirip ...................................... 26

Gambar 4.8

Isothermal pada Ra = 107 pada 4 sirip ...................................... 26

Gambar 4.9

Stream function pada Ra = 106 pada sirip tunggal ................... 27

Gambar 4.10 Stream function pada Ra = 106 pada 2 sirip.............................. 28 Gambar 4.11 Stream function pada Ra = 106 pada 3 sirip.............................. 28 Gambar 4.12 Stream function pada Ra = 106 pada 4 sirip ............................ 28 Gambar 4.13 Stream function pada Ra = 107 pada sirip tunggal................... 29 Gambar 4.14 Stream function pada Ra = 107 pada 1 sirip ............................. 29

x

Gambar 4.15 Stream function pada Ra = 107 pada 2 sirip ............................ 29 Gambar 4.16 Stream function pada Ra = 107 pada 4 sirip ............................. 30 Gambar 4.17 Vektor kecepatan pada Ra = 106 pada 1 sirip ........................... 30 Gambar 4.18 Kurva konvergensi untuk Ra=106 sirip tunggal ...................... 31 Gambar 4.19. Kurva konvergensi untuk Ra=107 sirip tunggal ........................ 31 Gambar 4.20 Domain dan Syarat Batas penelitian Pranowo dan Tri Iswanto 32 Gambar 4.21 Distribusi temperatur (a) dan stream function (b) yang dibuat penulis pada Ra = 106.................................................... 33 Gambar 4.22 Distribusi temperatur (a) dan stream function (b) yang dibuat penulis pada Ra = 107.................................................... 33 Gambar 4.23 Distribusi temperatur (a) dan stream function (b) yang dibuat Pranowo dan Tri Iswanto pada Ra = 106....................... 34 Gambar 4.24 Distribusi temperatur (a) dan stream function (b) yang dibuat Pranowo dan Tri Iswanto pada Ra = 107....................... 34 Gambar 4.25 Separasi di sekitar sirip yang dibuat penulis ............................ 36 Gambar 4.26 Separasi di sekitar sirip oleh F. Xu, J.C. Patterson dan C.Lei (2007) ...................................................................... 36

xi

DAFTAR TABEL Tabel 2.1 Koefisien Runge-Kutta orde-4 dari Carpenter dan Kennedy ......... 9 Tabel 2.2 Perbandingan skema beda hingga dan skema kompak turunan pertama.............................................................................. 10 Tabel 2.3 Perbandingan skema beda hingga dan skema kompak turunan kedua ................................................................................. 12 Tabel 4.1 Hasil Perhitungan dan Perbandingan untuk Ra=106 ....................... 35 Tabel 4.2 Hasil Perhitungan dan Perbandingan untuk Ra=107 ....................... 35

xii

DAFTAR NOTASI

a

: koefisien skema kompak

aM

: koefisien skema Runge-Kutta

b

: koefisien skema kompak

bM

: koefisien skema Runge-Kutta

c

: konstanta persamaan konveksi 1-D

g

: percepatan gravitasi (m/s2)

H

: tinggi kotak

HM

: variabel untuk skema Runge-Kutta

i,j

: indek nodal

k

: numerical wave number

Lr

: variabel referensi untuk panjang kotak

nx

: jumlah index arah x

ny

: jumlah indek arah y

Nu

: bilangan Nusselt

p

: tekanan

u

: kecepatan arah x

v

: kecepatan arah y

Vr

: variabel referensi untuk kecepatan

x,y

: koordinat

Pr

: bilangan Prandtl

Ra

: bilangan Rayleigh

t

: variabel waktu

tr

: variabel reverensi untuk waktu

Huruf Yunani 

: koefisien skema kompak



: koefisien ekspansi volumetri



: operator diferensial



: operator diferensial parsial

xiii



: konstanta metode kompresibilitas tiruan



: variabel generik

’

: variabel turunan pertama

”

: variabel turunan kedua



: densitas



: variabel temperatur



: vortisitas



: stream function



: jumlah

xiv

DAFTAR LAMPIRAN

Lampiran 1. Non-Dimensional Persamaan Atur .............................................. 39 Lampiran 2. Skema kompak beda-hingga........................................................ 42 Lampiran 3. Kurva Konvergensi ...................................................................... 44 Lampiran 4. Numerical Wavenumber .............................................................. 47 Lampiran 5. Program Konveksi Alami Pada Sirip Tunggal ............................ 49 Lampiran 6. Program Konveksi Alami Pada 2 Sirip........................................ 62 Lampiran 7. Program Konveksi Alami Pada 3 Sirip........................................ 79 Lampiran 8. Program Konveksi Alami Pada 4 Sirip........................................ 98 Lampiran 9. Program Tambahan Untuk Sirip Tunggal ................................... 120 Lampiran 10. Program Tambahan Untuk 2 Sirip ............................................ 127 Lampiran 11. Program Tambahan Untuk 3 Sirip ............................................ 136 Lampiran 12. Program Tambahan Untuk 4 Sirip ............................................ 146 Lampiran 13. Kode Program MATLAB ......................................................... 157

xv

BAB I PENDAHULUAN

1.1 Latar Belakang Masalah Dalam kehidupan sehari – hari sering sekali kita jumpai aplikasi mengenai perpindahan panas,salah satunya adalah perpindahan panas secara konveksi. Perpindahan panas secara konveksi berdasarkan jenis penyebab aliran fluida yang terjadi dikategorikan menjadi dua kategori, yaitu konveksi paksa dan konveksi alami. Konveksi paksa (forced convection) adalah konveksi dimana aliran fluida yang terjadi disebabkan adanya alat – alat eksternal,seperti fan,pompa,aliran udara atmosfer (angin). Sedangkan konveksi alami (natural convection) adalah konveksi yang terjadi karena fluida yang berubah densitasnya (kerapatannya) disebabkan proses pemanasan dan fluida ini bergerak naik karena adanya gaya apung (bouyancy force). Untuk meningkatkan perpindahan panas antara permukaan utama dan fluida di sekitarnya biasanya kita menggunakan sirip (fin). Sirip biasa digunakan dalam berbagai macam aplikasi,misalnya pada sistem pendingin ruangan,peralatan elektronik,tubin gas dan sebagainya,dengan udara sebagai media perpindahan panasnya. Penelitian mengenai fenomena perpindahan panas konveksi alami dengan menggunakan sirip tunggal maupun multiple sirip telah banyak dilakukan baik secara eksperimental maupun secara numerik. Penelitian secara eksperimen untuk mengetahui fenomena yang terjadi pada proses perpindahan panas konveksi alami dengan menggunakan sirip tunggal maupun multiple sirip membutuhkan biaya yang cukup mahal dan proses yang cukup rumit. Oleh karena itu, di zaman komputerisasi ini percobaan-percobaan dengan program komputer atau simulasi sangat diperlukan, hal ini bertujuan untuk mendukung penelitian secara eksperimen, menghemat waktu dan biaya serta keakuratannya.

Dan

dikembangkanlah

membutuhkan biaya yang jauh lebih murah.

1

penelitian

secara

numerik

yang

2

1.2 Perumusan Masalah Perumusan masalah dalam penelitian ini adalah bagaimana mensimulasikan secara numerik konveksi alami pada sirip tunggal dan multiple sirip dalam kotak 2D dengan menggunakan metode beda hingga orde-4.

1.3 Batasan Masalah Dalam penelitian ini masalah dibatasi sebagai berikut : a. Masalah pada penelitian ini dibatasi pada persoalan konveksi alami pada sirip tunggal dan multiple sirip dalam kotak 2D dengan menyelesaikannya menggunakan metode beda hingga orde-4 untuk memperoleh distribusi temperatur dan pola aliran (stream function). b. Aliran fluida diasumsikan sebagai aliran fluida tak mampat (incompressible flow). c. Penelitian dibatasi pada ruang 2 dimensi.

1.4 Tujuan Penelitian Adapun tujuan dari penelitian ini adalah untuk mengetahui fenomena yang terjadi pada aliran dan perpindahan panas konveksi alami pada sirip tunggal dan multiple sirip dalam kotak 2D,meliputi profil aliran fluida dan distribusi temperatur.

1.5 Manfaat Penelitian Manfaat dari penelitian ini : a. Untuk mengembangkan ilmu pengetahuan, terutama dalam bidang komputasi numerik dan perpindahan panas. b. Untuk mempelajari fenomena yang terjadi pada konveksi alami dengan menggunakan sirip tunggal maupun multiple sirip di sekitar kotak 2D.

3

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 perpindahan panas konveksi dan penjelasan mengenai metode beda hingga orde 4.

BAB III

: PELAKSANAAN PENELITIAN Berisi tentang alat dan bahan yang digunakan dalam penelitian, cara penelitian, diskritisasi persamaan atur.

BAB IV

: HASIL DAN PEMBAHASAN Berisi 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 N. Kasayapanand (2008) menyelidiki efek medan listrik konveksi alami pada sirip tunggal dan multiple sirip dalam sebuah kotak bujur sangkar dengan pemodelan numerik,dimana dipengaruhi oleh medan listrik, aliran dan temperaturnya. Parameter yang digunakan adalah tegangan, Rayleigh number, ukuran kotak, penyusunan elektroda, jumlah sirip, dan panjang sirip. Dihasilkan bahwa dengan angka Rayleigh yang sama maka koefisien perpindahan panas akan bertambah besar dengan menambah jumlah sirip dan panjang sirip.

Gambar 2.1. stream function dan isotherm pada Ra = 106

F. Xu, J.C. Patterson dan C. Lei (2007) mempelajari tentang konveksi alami pada sebuah kotak yang dipanaskan pada dinding samping dengan menggunakan tiga sirip yang diselesaikan secara numerik. Berkaitan dengan adanya sirip, maka terjadi separasi (pemisahan) aliran panas di sekitar sirip tersebut.

4

5

Gambar2. 2. separasi di sekitar sirip

Wilson dan Demuren(1998) menggunakan skema kompak beda hingga untuk diskritasi ruang dan skema Runge-Kutta untuk diskritasi waktu pada simulasi aliran fluida tak mampat. Pada penelitian ini skema kompak beda hingga digunakan untuk diskritisasi turunan ruang dan skema Runge-Kutta orde-empat untuk diskritasi turunan waktu. Le Querre (1990) menggunakan algoritma pseudo – spectral Chebsyev untuk meneliti konveksi alami pada kotak 2D dengan dinding bawah di panasi dan dinding atas adiabatic. Dengan metode ini dapat menghilangkan osilasi numeric dan mencapai hasil yang akurat hingga nilai Ra 108. Pranowo dan Priyo tri Iswanto (1999) menyelesaikan persamaan Navier – Stoke 2 dimensi dengan menggunakan primitive variable pada non staggered grid dengan diskritisasi beda hingga. Hasil simulasi dengan metode ini menunjukkan hasil yang akurat untuk Ra = 106 dan Ra = 107. Aris Sulistyono (2006) melakukan penelitian untuk mengetahui fenomena yang terjadi pada konveksi alami kotak 2D dengan berbagai variasi kemiringan.

6

2.2 Dasar Teori 2.2.1 Sirip Sirip digunakan dalam banyak alat penukar kalor untuk meningkatkan luasan perpindahan panas. Aplikasi sirip banyak dijumpai dalam sistem pendingin ruangan,peralatan elektronik,turbin gas dan sebagainya, dengan udara merupakan media perpindahan panasnya. Dalam alat penukar panas sirip terbagi dari berbagai macam tipe, mulai dari bentuk yang sederhana, seperti sirip segiempat (rectangular), silindris, annular, tirus (tapered) atau pin, sampai kombinasi dari berbagai geometri yang berbeda telah digunakan. Tipe sirip yang digunakan tergantung dari proses permesinan dan ruang yang tersedia dalam peralatan pembangkit panas yang terlibat dalam proses pendinginan. Dalam desain dan konstruksi dari berbagai macam peralatan perpindahan panas, bentuk-bentuk sederhana seperti;

silinder,

batang

dan

plat biasa

diterapkan pada aliran panas antara sumber panas dan penyerap panas (heat source

and

heat sink). Permukaan-permukaan penyerap panas maupun

pembuang panas masing-masing dikenal sebagai

permukaan

utama

(prime

surface). Apabila permukaan utama diperluas dengan permukaan tambahan seperti dalam gambar 2.3, maka gabungan antara kedua permukaan tersebut dinamakan

permukaan

yang diperluas

(extended

surface). Elemen

yang

digunakan untuk memperluas permukaan utama dikenal sebagai sirip (sirip).

Gambar 2.3. Beberapa contoh jenis extended surface: (a) sirip longitudinal (memanjang) dengan profil segiempat (b)pipa silindris dengan sirip berprofil segiempat (c) sirip longitudinal dengan profil trapezioda (d) sirip longitudinal dengan profil parabola (e) pipa silindris dengan sirip radial berprofil segiempat (f) pipa silindris dengan sirip radial berprofil trapezoida (g)cylindrical spine (h)truncated conical spine (i) truncated parabolic spine

7

2.2.2 Konveksi Alami Konveksi alami adalah perpindahan panas antara suatu permukaan dan fluida yang mengalir diatasnya. Aliran fluida disebabkan oleh adanya perbedaan densitas fluida yang ditimbulkan oleh pemanasan dan pendinginan. Densitas fluida akan berkurang jika fluida mendapat pemanasan sehingga fluida akan mengapung dan daerah yang ditinggalakan akan diisi oleh fluida yang relatif dingin. Fluida yang relatif panas jika mendekati dinding yang relatif dingin densitasnya akan meningkat sehingga akan mengalir turun akibat tarikan gaya grafitasi. Dengan demikian densitas merupakan driving force sirkulasi fluida. Konveksi alami memegang peranan penting dalam rekayasa industri seperi: perancangan alat penukar kalor, perancangan ventilasi, pendinginan transformator, pendinginan kabel bawah tanah dan pendinginan komponen elektronika. Pada penelitian konveksi alami model matematika yang dipakai adalah persamaan kontinyuitas, persamaan Navier-Stokes dan persamaan energi.

2.2.2.1 Persamaan Atur Konveksi Alami Untuk permasalahan 2-D persamaan atur konveksi alami dalam bentuk variabel tak berdimensi adalah sebagai berikut (Le Quere,1990): u v  0 x y

(1)

u u u p Pr   2u  2u   u  v    0.5  2  2  t x y x Ra  x y 

(2)

v v v p Pr   2 v  2 v   u  v    0.5  2  2   Pr t x y y Ra  x y 

(3)

θ θ θ 1   2θ  2θ   u  v  0.5  2  2  t x y Ra  x y 

(4)

Persamaan di atas diperoleh dengan membagi variabel berdimensi dengan variabel referensi, untuk panjang adalah Lr=H, untuk kecepatan Vr=(/H)Ra-0.5, dimana Ra=(gTH3)/(), untuk variabel waktu tr=(H2/) Ra-0.5, untuk temperatur didefinisikan sebagai berikut : =(T-Tr)/(Th-Tc), Tr=(Th+Tc)/2 dan Pr=(/).

8

2.2.2.2 Diskritisasi Waktu Diskritisasi waktu untuk persamaan momentum menggunakan skema RungeKutta orde-4 dari Williamson(Wilson dan Demuren,1998) yang didefinisikan sebagai berikut :

u M 1  u M  b M 1tH iM

(5)

dimana : t

= langkah waktu

bM = koefisien skema Runge-Kutta aM = koefisien skema Runge-Kutta uiM = komponen kecepatan arah xi pada sub tingkat ke-M Pi M = Tekanan M

Hi

 u j  x ui  Pi  M

M

Pr Ra

0.5

M 1

 xxui  a H i M

M

Tabel 2.1 Koefisien Runge-Kutta orde-4 dari Carpenter dan Kennedy M 1 2 3 4 5

aM 0 -0.41789047 -1.19215169 -1.69778469 -1.51418344

bM 0.14965902 0.37921031 0.82295502 0.69945045 0.15305724

2.2.2.3 Diskritisasi Ruang Skema beda-hingga orde-2 untuk turunan pertama memiliki galat dispersi yang besar, sedangkan skema kompak beda hingga memiliki kelebihan yaitu akurasi tinggi, fleksibel dan pengoperasiannya lebih mudah. a. Turunan pertama. Bentuk diskritisasi turunan pertama dengan pendekatan skema kompak beda hingga orde-4 dirumuskan oleh Lele(Wilson dan Demuren, 1998). Bentuk persamaannya adalah seperti berikut :

i' 1  i'  i' 1 

a b (i 1  i 1 )  ( i  2   i  2 ) 2x 4x

(6)

9

dengan : x

= Lx/Nx

Nx

= jumlah grid point

 i'

= turunan pertama dari variabel  i terhadap x

, a, b = koefisien skema kompak Turunan terhadap y dan z dapat dilakukan dengan cara yang sama. Untuk skema orde-empat maka ; =1/4, a=3/2 dan b=0. Perbandingan skema ekplisit beda-hingga dan skema kompak beda hingga dari turunan pertama ditunjukkan dalam tabel 2.2. Di sini terlihat bahwa skema kompak beda hingga memiliki grid stensil yang lebih sedikit, koefisien galat pemenggalan berkurang menjadi ¼ untuk orde-4 dari koefisien beda tengah ekplisit untuk orde yang sama. Tabel 2.2 Perbandingan skema beda hingga dan skema kompak turunan pertama Skema

Kesalahan pemenggalan

Jumlah stensil

Beda tengah orde-4

(-4/5!)( x) 4  5 

5

Kompak orde-4

(-1/5!)( x) 4  5 

3

Menurut Hu dkk(1996) resolusi dari diskritisasi turunan pertama

dapat

dianalisa dengan mentransformasi persamaan konveksi 1-D sebagai berikut:

  c 0 t x 1 N  u   al  j l     x  j x l  N

(7)

(8)

~ Dalam mode Fourier   (t)e ikx maka : ~   ikx  e t t ~  j l  eik ( x  x )

(10)

Sehingga persamaan konveksi 1-D menjadi : ~  ikx c N ~ ik ( x lx ) e  0  al e t x l  N

(11)

(9)

10

~  c N ~ iklx   al e  0 t x l  N

(12)

~  ~  ick *  0 t

(13)

dimana : k* 

i N  al eiklx x l  N

(14)

k* adalah numerical wave number. Numerical wave number untuk skema kompak beda hingga dari turunan pertama dalah :

b   a sin(kx)  sin(2kx)   1 2 k*    x  2 cos(kx)  1   

(15)

Simpangan dari kurva real(k*) terhadap k menunjukkan galat dispersi. 3.5

eksak 2nd-order central

3

4th-order central

real(k*dx)

2.5

4th-order compact 6th-order compact

2 1.5 1 0.5 0 0

0.5

1

1.5

2

2.5

3

3.5

kdx

Gambar 2.4. Galat dispersi untuk pendekatan numeric dari turunan pertama

Syarat batas diselesaikan dengan skema kompak orde-3 dengan persamaan sebagai berikut : 1'   bs '2 

1 3  absii x i 1

(16)

11

 bs  2 dan abs1  5 / 2 , abs2  2 , abs1  1 / 2 adalah koefisien orde-3 dari syarat batas pada i=1. Persamaan yang sama juga digunakan untuk syarat batas pada i=N. b. Turunan kedua. Persamaan skema kompak beda hingga untuk turunan kedua adalah sebagai berikut : a i 1  2i  i 1   b 2 i 2  2i  i 2 (17) "i 1  "i  "i 1  2 x  4x  dimana :

 "i

= turunan kedua dari variabel  i terhadap x

 , a, b = koefisien skema kompak beda hingga turunan kedua Untuk orde-empat, =1/10, a  6 / 5 , b=0 Perbandingan antara skema beda-hingga ekplisit dan skema kompak beda hingga ditunjukkan dalam tabel 2.3. Di sini terlihat bahwa skema kompak beda hingga memiliki stensil lebih sedikit, koefisien galat pemenggalan berkurang menjadi ½ untuk orde-4 dari koefisien beda tengah ekplisit untuk orde yang sama. Tabel 2.3 Perbandingan skema beda hingga dan skema kompak turunan kedua Skema Beda tengah orde-4 Kompak orde-4

Kesalahan pemenggalan

Jumlah stensil

(-8/6!)(x)4(6)

5

(-3.6/6!)(x)4(6)

3

Analisa resolusi untuk turunan kedua dari pendekatan numerik skema kompak beda hingga dilakukan dengan cara yang sama dengan analisa turunan pertama. Numerical wave number untuk skema kompak beda hingga dari turunan kedua adalah :

b   a1  cos(kx)  1  cos(2kx)  1 2 (k*)2   2  (x)  1  2 cos(kx)   

(18)

Deviasi dari (k*x)2 terhadap (kx)2 ditampilkan dalam bentuk grafik gambar 2.5 untuk beberapa skema beda-hingga.

12

10

Eksak "2nd-order central"

9

"4th-order central"

8

"4th-order compact" "6th-order compact"

(k*dx)**2

7 6 5 4 3 2 1 0 0

0.5

1

1.5

2

2.5

3

3.5

kdx

Gambar 2.5. Galat disipasi untuk pendekatan numerik dari turunan kedua

Galat disipasi dari berbagai skema beda-hingga tampak pada gambar 2.5 dapat diketahui bahwa nilai numerical wavenumber untuk skema kompak lebih mendekati nilai exact wavenumber. Kondisi batas pada i=1 dan i=N diselesaikan dengan skema kompak orde-3 sebagai berikut :

"i   bs"2 

1 x 2

4

a i 1

i

(19)

bsi

dimana, bs =11 dan abs1=13, abs2=-27, abs3=15 dan abs4 =-1 adalah koefisien skema kompak orde-3.

2.2.2.4 Metode Kompresibilitas Tiruan (Artificial Compressibility) Konsep metode kompresibilitas tiruan adalah menambahkan turunan terhadap waktu pada persamaan kontinyuitas. Bentuk modifikasi persamaan adalah : p  V  0 t

(20)

13

dimana  adalah konstanta positip. Persamaan ini tidak mempunyai arti fisik jika kondisi tunak belum tercapai.

BAB III PELAKSANAAN PENELITIAN

3.1 Alat dan Bahan a. Komputer pribadi dengan spesifikasi : 

Intel Pentium Dual CPU E2140 @ 1.6GHz



Memori 512 MB

b. Perangkat lunak Mikrosoft Fortran Power Station 4.0 c. Perangkat lunak Matlab 6.1 Realease 12 d. Printer Canon iP 1880

3.2 Garis Besar penelitian Penelitian dilakukan dengan cara membuat implementasi program untuk menyelesaikan persamaan momentum, persamaan energi dan persamaan kontinyuitas dengan pendekatan skema kompak orde-empat dan skema RungeKutta orde-empat. Program dibuat dalam Bahasa Fortran dan untuk visualisasi hasil program menggunakan perangkat lunak Matlab 6. Langkah-langkah penelitian yang dilakukan adalah seperti berikut : 1. Mengumpulkan literatur 2. Mempelajari literatur a. Mempelajari penelitan-penelitian yang pernah dilakukan b. Mempelajari persamaan atur yang berhubungan dengan permasalahan 3. Merencanakan algoritma program a. Membuat diskritisasi persamaan atur b. Menyusun bagan alir program 4. Menulis bagan alir dalam bahasa program (Fortran) 5. Menjalankan program 6. Memperbaiki kesalahan pemrograman a. kesalahan penulisan b. kesalahan algoritma 7. Membuat visualisasi hasil program dengan perangkat lunak Matlab

14

15

8. Menyusun laporan penelitian Garis besar penelitian tersebut dapat dibuat diagram alir sebagai berikut :

Mulai

Mengumpulkan dan Mempelajari literatur - literatur

Membuat diskritisasi persamaan atur

Membuat algoritma program

Menulis bagan alir dalam bahasa fortran

Menjalankan program tidak Program benar ya Membuat visualisasi dengan Matlab

Analisa hasil

Kesimpulan

Selesai

16

3.3 Diskritisasi Persamaan Atur Persamaan atur konveksi alami terdiri dari persamaan kontinuitas, persamaan momentum dan persamaan energi. Model matematika dari persamaan atur konveksi alami terdiri dari persamaan diferensial parsial orde-1 dan orde-2. Agar persamaan atur konveksi alami dapat diaplikasikan dalam bahasa program maka terlebih dahulu dibuat diskritisasi persamaan atur. Diskritisasi waktu dilakukan dengan skema Runge-Kutta orde-4 dan diskritisasi ruang dengan skema kompak beda hingga orde-4. Matrik yang terbentuk dari diskritisasi turunan ruang adalah matrik tridiagonal yang bisa diselesaikan dengan algoritma Thomas.

3.3.1 Diskritisasi persamaan momentum Diskritisasi persamaan momentum dengan skema Runge-Kutta adalah seperti berikut : Kecepatan arah x (u) (21)

uiM, j  uiM, j  b M 1tH iM, j H i,Mj  uiM, j uxiM, j  viM, j uyiM, j  pxiM, j 

Pr uxxiM, j  uyyiM, j   a M H iM, j 0.5 Ra

(22)

Kecepatan arah y (v) (23)

viM, j 1  viM, j  b M 1tH iM, j

H iM, j  ui , j vxi , j  vi , j vyi , j  pyi , j 

Pr vxxi, j  vyyi, j  Pri, j  a M HiM, j 0.5 Ra

(24)

Diskritasi turunan ruang dengan skema kompak orde-4 adalah seperti berikut : Diskritisasi turunan pertama 1 M 1 3 uxi 1, j  uxiM, j  uxiM1, j  uiM1, j  uiM1, j 4 4 4x





(25)

1 M 1 3  uyi , j 1  uyiM, j  uxiM, j 1  uiM, j 1  uiM, j 1  4 4 4y

(26)

17





(27)

1 M 1 3 M vxi 1, j  vxiM, j  vxiM1, j  vi 1, j  viM1, j 4 4 4x 1 M 1 3 viM, j 1  viM, j 1  vyi , j 1  vyiM, j  vxiM, j 1  4 4 4y



1 M 1 3 pxi 1, j  pxiM, j  pxiM1, j  piM1, j  piM1, j 4 4 4x

(28)



(29)

1 M 1 3  pyi , j 1  pyiM, j  pxiM, j 1  piM, j 1  piM, j 1  4 4 4y

(30)

Diskritisasi turunan kedua





(31)





(32)

1 1 6 uxxiM1, j  uxxiM, j  uxxiM1, j  uiM1, j  2uiM, j  uiM1, j 2 10 10 5x  1 1 6 uyyiM, j 1  uyyiM, j  uyyiM, j 1  uiM, j 1  2uiM, j  uiM, j 1 2 10 10 5y 





(33)





(34)

1 1 6 vxxiM1, j  vxxiM, j  vxxiM1, j  viM1, j  2viM, j  viM1, j 2 10 10 5x  1 1 6 vyyiM, j 1  vyyiM, j  vyyiM, j 1  viM, j 1  2viM, j  viM, j 1 2 10 10 5y 

3.3.2. Diskritisasi persamaan energi

iM, j  iM, j  b M 1tH iM, j H iM, j  uiM, jxiM, j  viM, jyiM, j 

(35)





1 xxiM, j  yyiM, j  a M H iM, j Ra 0.5

(36)

Diskritasi turunan ruang dengan skema kompak orde-4 adalah seperti berikut : Diskritisasi turunan pertama



1 M 1 3 xi 1, j  xiM, j  xiM1, j  iM1, j  iM1, j 4 4 4x



1 M 1 3 iM, j 1  iM, j 1  yi , j 1  yiM, j  xiM, j 1  4 4 4y

(37)

(38)

18

Diskritisasi turunan kedua



1 1 6 xxiM1, j  xxiM, j  xxiM1, j  iM1, j  2iM, j  iM1, j 2 10 10 5x 



1 1 6 yyiM, j 1  yyiM, j  yyiM, j 1  iM, j 1  2iM, j  iM, j 1 2 10 10 5y 





(39)

(40)

3.3.3. Diskritisasi metode kompresibilitas tiruan.

piM, j 1  piM, j  b M 1tH iM, j

(41)

H iM, j  uxiM, j  uyiM, j  a M H iM, j

(42)

3.4 Diskritisasi Syarat Batas Dalam penelitian ini kasus yang dibahas adalah konveksi alami dalam kotak 2D dengan dinding bawah dan atas merupakan dinding adiabatis, dinding kiri mendapat pemanasan dan dinding kanan mendapat pendinginan,dengan grid 201 x 201. Pada seluruh dinding kecepatan bernialai nol sedangkan syarat batas tekanan dan temperatur adalah seperti berikut :  0 y p 0 y

p 0 x

p 0 x

 = 0.5

 = 0.5 p 0 y

 = -0.5

 0 y Gambar 3.1 Domain dan syarat batas

19

a. Syarat batas kecepatan Turunan pertama. Untuk i=1 dan i=nx

ux1M, j 



1  3u5M, j  16u4M, j  36u3M, j  48u2M, j  25u1M, j 12x

(43)



1 3unxM4, j  16unxM3, j  36unxM2, j  48unxM1, j  25unxM, j 12x 1   3v5M, j  16v4M, j  36v3M, j  48v2M, j  25v1M, j 12x

uxnxM, j 

vx1M, j





(44)



(45)



(46)

uyiM,1 

1   3uiM,5  16uiM, 4  36uiM,3  48uiM, 2  25uiM,1  12y

(47)

uyiM,98 

1  3uiM,94  16uiM,95  36uiM,96  48uiM,97  25uiM,98  12y

(48)

vxnxM, j 





1 3vnxM4, j  16vnxM3, j  36vnxM2, j  48vnxM1, j  25vnxM, j 12x

Untuk j=1 dan j=ny

uyiM,102 

1  3uiM,106  16uiM,105  36uiM,104  48uiM,103  25uiM,102  12y

(49)

uyiM,ny 

1 3uiM,ny4 16uiM,ny3  36uiM,ny2  48uiM,ny1  25uiM,ny  12y

(50)

vyiM,1 

1   3viM,5  16viM, 4  36viM,3  48viM, 2  25viM,1  12y

(51)

vyiM,98 



1 3viM,94  16viM,95  36viM,96  48viM,97  25viM,98 12y

vyiM,102  vyiM,ny 



1   3viM,106  16viM,105  36viM,104  48viM,103  25viM,102  12y

1  3viM,ny4  16viM,ny3  36viM,ny2  48viM,ny1  25viM,ny  12y

(52)

(53)

(54)

Turunan kedua. Untuk i=1 dan i=nx uxx1M, j  11uxx2M, j 



1 13u1M, j  27u2M, j  15u3M, j  u4M, j 2 x 



(55)

20

uxxnxM, j  11uxxnxM1, j 



vxx1M, j  11vxx2M, j



1 13unxM, j  27unxM1, j  15unxM2, j  unxM3, j 2 x 1  13v1M, j  27v2M, j  15v3M, j  v4M, j 2 x 

vxxnxM, j  11vxxnxM1, j 



(56)



(57)



1 13vnxM, j  27vnxM1, j  15vnxM2, j  vnxM3, j 2 x



(58)

Untuk j=1 dan j=ny





(59)





(60)

uyyiM,1  11uyyiM, 2 

1 13uiM,1  27uiM, 2  15uiM,3  uiM, 4 2 y 

uyyiM,98  11uyyiM,97 

1 13uiM,98  27uiM,97  15uiM,96  uiM,95 2 y 



1 13uiM,102  27uiM,103  15uiM,104  uiM,105 2 y 

uyyiM,102  11uyyiM,103  uyyiM,ny  11uyyiM,ny1 



1 13uiM,ny  27uiM,ny1  15uiM,ny2  uiM,ny3 2 y 



1 13viM,1  27viM, 2  15viM,3  viM, 4 2 y 

vyyiM,98  11vyyiM,97 

1 13v M  27viM,97  15viM,96  viM,95 y 2 i,98

vyyiM,102  11vyyiM,103  vyyiM,ny  11vyyiM,ny1 

(63)



(64)



1 13viM,102  27viM,103  15viM,104  viM,105 2 y 



(61)

(62)



vyyiM,1  11vyyiM, 2 







1 13viM,ny  27viM,ny1  15viM,ny2  viM,ny3 2 y 





(65)

(66)

b. Syarat batas tekanan Untuk i=1 dan i=nx

px1M, j  0

(67)

pxnxM, j  0

(68)

21

Untuk j=1 dan j=ny

pyiM,1  0

(69)

pyiM,ny  0

(70)

c. Syarat batas temperatur Turunan pertama Untuk i=1 dan i=nx



x1M, j 

1  35M, j  16 4M, j  363M, j  48 2M, j  251M, j 12x

xnxM, j 

1 3 nxM4, j  16 nxM3, j  36 nxM2, j  48 nxM1, j  25 nxM, j 12x





(71)



(72)

Untuk j=1 dan j=ny

yiM,1  0 xiM,98 

(73)



1 3iM,94  16iM,95  36iM,96  48iM,97  25iM,98 12x

xiM,102 



(74)



1  3iM,106  16iM,105  36iM,104  48iM,103  25iM,102 12x

yiM,ny  0



(75) (76)

Turunan kedua Untuk i=1 dan i=nx

xx1M, j  11xx2M, j 



1 131M, j  27 2M, j  153M, j   4M, j 2 x 

xxnxM, j  11xxnxM1, j 





1 13 nxM, j  27 nxM1, j  15 nxM2, j   nxM3, j 2 x

(77)



(78)

Untuk j=1 dan j=ny

yyiM,1  11yyiM, 2 



1 13iM,1  27iM, 2  15iM,3  iM, 4 2 y 



(79)

22





1 13iM,98  27iM,97  15iM,96  iM,95 2 y  1 yyiM,102  11yyiM,103  13iM,102  27iM,103  15iM,104  iM,105 2 y 

yyiM,98  11yyiM,97 



yyiM,ny  11yyiM,ny1 



1 13iM,ny  27iM,ny1  15iM,ny2  iM,ny3 2 y 

(80)



(81)



(82)

3.5 Algoritma Pemrograman Algoritma pemrograman dari sistem persamaan diatas adalah sebagai berikut : 1. Tentukan kondisi awal (t=0), dan kondisi batas untuk semua variabel (u,v,,p). 2. Hitung

turunan

pertama

(ux,uy,vx,vy,x,y,px,py)

dari dan

kecepatan, turunan

temperatur

kedua

dari

dan

tekanan

kecepatan

dan

temperatur(uxx,uyy,vxx,vyy,xx,yy) dengan skema kompak orde-empat. 3. Hitung kecepatan(u,v) dengan skema Runge-Kutta orde-empat. 4. Hitung tekanan dengan metode artificial compressibility. 5. Hitung temperatur()dengan skema Runge-Kutta orde-empat. 6. Periksa apakah sudah mencapai batas perhitungan atau belum, jika belum kembali ke langkah 2, jika sudah ke langkah 7. 7. Tulis hasil. 8. Selesai.

23

Bagan alir program yang akan dibuat adalah sebagai berikut: MULAI

DATA AWAL

SYARAT BATAS

TENTUKAN TURUNAN PERTAMA UNTUK u,v,p, DAN TURUNAN KEDUA UNTUK u,v,

SELESAIKAN PERSAMAAN MOMENTUM UNTUK MEMPEROLEH Um+1 DAN vm+1 HITUNG TEKANAN pm+1 DENGAN METODE ARTIFICIAL COMPRESSIBILITY

SELESAIKAN PERSAMAAN ENERGI UNTUK MEMPEROLEH m+1

PERIKSA KONVERGENSI ? Y TULIS HASIL

SELESAI

T

BAB IV HASIL DAN PEMBAHASAN

4.1 Simulasi Konveksi Alami Sirip Tunggal Dan Multiple Sirip Simulasi secara konveksi alami pada sirip tunggal dan multiple sirip ditampilkan dengan susunan grid sebesar 201 x 201, bilangan Prandtl (Pr) = 0.71 dan langkah waktu dt = 0.0025 serta angka Rayleigh yang digunakan adalah Ra = 106 dan 107. Hasil simulasi disrtibusi temperatur selengkapnya dapat dilihat pada gambar berikut :

Gambar 4.1 isothermal pada Ra = 106 pada sirip tunggal

Gambar 4.2 isothermal pada Ra = 106 pada 2 sirip

24

25

Gambar 4.3 isothermal pada Ra = 106 pada 3 sirip

Gambar 4.4 isothermal pada Ra = 106 pada 3 sirip

Gambar 4.5 isothermal pada Ra = 107 pada sirip tunggal

26

Gambar 4.6 isothermal pada Ra = 107 pada 2 sirip

Gambar 4.7 isothermal pada Ra = 107 pada 3 sirip

Gambar 4.8 isothermal pada Ra = 107 pada 4 sirip

27

Dari hasil di atas maka dapat kita tinjau secara visual mengenai gambar distribusi temperatur. Dari gambar 4.1 sampai 4.8 dapat dilihat bahwa pergerakan fluida panas akan bergerak keatas karena adanya gaya apung (buoyancy force),hal ini disebabkan karena density yang turun karena temperatur, sedangkan fluida dingin bergerak ke bawah karena density lebih besar dan karena adanya gaya gravitasi. Semakin banyak sirip maka dapat dilihat arah perpindahan panasnya. Distribusi temperatur yang relatif panas dibagian kiri atas akan semakin condong ke kanan dan distribusi temperatur yang relatif dingin pada bagian kiri bawah semakin condong ke kiri. Dengan adanya peningkatan Ra maka akan membuat lapis batas termal di dinding menipis sehingga gradien temperatur di dinding meningkat,hal ini terjadi karena kecepatan fluida yang membawa panas juga meningkat seiring peningkatan Ra.

Pola aliran fuida dapat dilihat pada gambar berikut :

Gambar 4.9 stream function pada Ra = 106 pada sirip tunggal

28

Gambar 4.10 stream function pada Ra = 106 pada 2 sirip

Gambar 4.11 stream function pada Ra = 106 pada 3 sirip

Gambar 4.12 stream function pada Ra = 106 pada 4 sirip

29

Gambar 4.13 stream function pada Ra = 107 pada sirip tunggal

Gambar 4.14 stream function pada Ra = 107 pada 1 sirip

Gambar 4.15 stream function pada Ra = 107 pada 2 sirip

30

Gambar 4.16 stream function pada Ra = 107 pada 4 sirip

Dari gambar 4.9 sampai gambar 4.16 dapat kita lihat pola stream fungtionnya, dari gambar diatas dapat dilihat bahwa pada Ra = 106 telah muncul separasi, separasi ini timbul akibat adanya suatu gerakan berputar pada vektor kecepatan,untuk lebih jelasnya dapat ditampilkan pada gambar 4.17.

Gambar 4.17 vektor kecepatan pada Ra = 106 pada 1 sirip

31

Kurva Konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 4.18. Kurva konvergensi untuk Ra=106 sirip tunggal Kurva Konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 4.19. Kurva konvergensi untuk Ra=107 sirip tunggal Gambar 4.18 dan 4.19 menunjukkan peningkatan Ra menyebabkan kecepatan konvergensi berkurang. Hal ini disebabkan peningkatan Ra membuat suku difusi berkurang pengaruhnya terhadap perhitungan. Suku difusi mempunyai sifat sebagai peredam, sehingga perhitungan stabil. Dengan melemahnya suku difusi maka suku adveksi menguat dan menggantikan dominasi suku difusi. Penguatan suku adveksi membuat sistem persamaan atur cenderung bersifat hiperbolik.

32

4.2 Validasi Program Untuk meguji validitas dari program yang telah dibuat, hasil dari proses simulasi dibandingkan secara visual maupun perhitungan dengan hasil yang telah dilakukan oleh Pranowo dan Tri Iswanto pada domain kotak bujur sangkar 2D dengan aspek rasio 1 :1, dengan kondisi dinding bawah dan atas merupakan dinding adiabatis, dinding kiri mendapat pemanasan dan dinding kanan mendapat pendinginan.  0 y

p 0 y  = 0.5

p 0 x

p 0 x

 = -0.5

p 0 y  0 y

Gambar 4.20. Domain dan Syarat Batas penelitian Pranowo dan Tri Iswanto

Disrtibusi temperatur dan pola aliran hasil penelitian yang dibuat Pranowo dan Tri Iswanto ditunjukkan pada gambar 4.23 dan gambar 4.24 sedangkan distribusi temperatur dan pola aliran hasil penelitian yang dibuat penulis ditunjukkan oleh gambar 4.21 dan 4.22.

33

(a)

(b)

Gambar 4.21. Distribusi temperatur (a) dan stream function (b) yang dibuat penulis pada Ra = 106

(a)

(b)

Gambar 4.22. Distribusi temperatur (a) dan stream function (b) yang dibuat penulis pada Ra = 107

34

(a)

(b)

Gambar 4.23. Distribusi temperatur (a) dan stream function (b) yang dibuat Pranowo dan Tri Iswanto pada Ra = 106

Gambar 4.24. Distribusi temperatur (a) dan stream function (b) yang dibuat Pranowo dan Tri Iswanto pada Ra = 107

Hasil perhitungan pada penelitian kali ini dibandingkan dengan penelitian Le Querre pada Ra = 106 dan Ra = 107 yang diperliatkan pada tabel 4.1 dan tabel 4.2.

35

Tabel.4.1Hasil Perhitungan dan Perbandingan untuk Ra=106 middle max X Y umax(1/2,y) Y vmax(x,1/2) X Nuwall Numiddle Numax Y Numin Y

Sekarang 0.0164552 0.0169117 0.150000 0.550000 0.0649050 0.850000 0.220236 0.040 8.73394 8.82299 17.1575 0.030 0.98427 1.0

Le Quere 0.016384 0.016811 0.1500 0.5470 0.064834 0.850 0.2206 0.038 8.8252 8.8244 17.5343 0.039 0.97948 1.0

Beda (%) 0.43 0.60 0.00 0.55 0.11 0.00 0.17 5.26 1.03 0.02 2.15 23.08 0.49 0.00

Tabel 4.2 Hasil Perhitungan dan Perbandingan untuk Ra=107 middle max X Y umax(1/2,y) Y vmax(x,1/2) X Nuwall Numiddle Numax Y Numax Y

sekarang 0.00936819 0.00964076 0.0866667 0.553333 0.0473129 0.88 0.221048 0.02 16.2068 16.5638 40.3192 0.00667 1.37516 1.0

Le Quere 0.00928496 0.00953872 0.086 0.556 0.046986 0.879 0.21118 0.021 16.523 16.523 39.3947 0.018 1.36635 1.0

Beda (%) 0.90 1.07 0.78 0.48 0.70 0.11 4.67 4.76 1.91 0.25 2.35 62.94 0.64 0.00

Dari gambar 4.21 dan 4.22 dapat dilihat bahwa secara visual hasil penelitian ini menunjukkan kesamaan yang baik dengan penelitian dari Pranowo dan Tri Iswanto dan dari tabel 4.1 dan 4.2 hasil perhitungan pada penelitian ini menunjukkan kedekatan yang baik dengan hasil penelitian Le Querre sehingga dengan metode ini dapat diterima sebagai validasi. Untuk domain dengan sirip, akan dinbandingkan dengan hasil percobaan dibuat oleh F. Xu, J.C. Patterson dan C.Lei (2007) dengan mengacu timbulnya

36

sutatu separasi disekitar sirip, ditunjukkan pada gambar 4.26 sedangkan pola aliran hasil penelitian yang dibuat penulis ditunjukkan oleh gambar 4.25.

Gambar 4.25 Separasi di sekitar sirip yang dibuat penulis

Gambar 4.26 separasi di sekitar sirip oleh F. Xu, J.C. Patterson dan C.Lei (2007)

BAB V PENUTUP

5.1 Kesimpulan Dari penelitian dan pembahasan yang telah dilakukan, dapat ditarik beberapa kesimpulan yaitu : a. Hasil penelitian pada domain busur sangkar tanpa sirip mempunyai kedekatan visual yang baik dengan hasil penelitian dari Pranowo dan Tri Iswanto baik pada Ra = 106 maupun pada Ra = 107. b. Metode yang digunakan untuk kasus konveksi alami pada single fin dan multiple fins dalam kotak 2D mampu mensimulasikan pola aliran dan distribusi temperatur sampai Ra = 107. c. Separasi timbul karena adanya pusaran pada pola aliran disekitar sirip. d. Separasi muncul pada daerah sekitar sirip saat Ra 106 dan 107. e. Semakin besar angka Rayleigh maka kerapatan perpindahan panas disekitar dinding akan meningkat. f. Semakin banyak jumlah sirip maka perpindahan panas akan meningkat. 5.2 Saran Skema kompak orde-tinggi selain memiliki akurasi yang baik mempunyai bentuk yang sederhana dan mudah diaplikasikan. Bagi para pembaca yang berminat penelitian ini masih terbuka kemungkinan untuk penyelesaian kasus 3-D atau penyelesaian kasus 2-D dengan menggunaan grid yang tidak seragam.

37

DAFTAR PUSTAKA

Hoffmann, K.A. 1989. Computational Fluid Dynamic for Engineers. Austin, Texas: A Publication of Engineering Education System. Holman, J.P. 1988. Perpindahan Kalor. Jakarta: Erlangga. Kasayapanand, N. 2008. A Computational Fluid Dynamics Modeling of Natural Convection in Siripned Enclosure Under Electric Field. Applied Thermal Engineering 29 (2009) 131-141. Lemos, C.M. FDFlow: a Fortran-77 Solver for 2-D Incompressible Fluid Flow. Computers & Geosciences, Vol. 20 (1994): pp.265-261. Pranowo dan Priyo Tri Iswanto. 1999. Analisis Numerik Konveksi Alami Dalam Kotak dengan Primitive Variable pada Grid Kolokasi. Makalah Seminar Regional Universitas Sanata Dharma Yogyakarta. Prijono, Arko. 1985. Prinsip-prinsip Perpindahan Panas, Terjemahan dari Principles Of Heat Transfer Third Edition (Frank Kreith, 1958). Quere, P.L. 1990. Accurate Solutions to The Square Thermally Driven Cavity at High Rayleigh Number. International Journal of Computers & Fluids, Vol.20, No. 1, hal. 29-41. Sulistiyono, Aris. 2006. Simulasi Numeric Konveksi Alami Dalam Kotak 3 Dimnsi Dengan Variasi Kemiringan Dengan Metode Beda Hingga. Wilson, Robert V., and Demuren, Ayodeji O., 1998, Higher-Order Compact Schemes for Numerical Simulation of Incompressible Flows, ICASE Report No. 98-13, NASA Langley Research Center, Hampton. Xu, F., J.C. Patterson and C. Ley. 2007. Natural convection adjacent to sidewall with three fins in a differentially heated cavity. Anziam J. 48 pp. C806C819.

38

38

Lampiran 1. Non-Dimensionalisasi Persamaan Atur

NON-DIMENSIONALISASI PERSAMAAN ATUR Persamaan atur dibuat tanpa dimensi dengan membagi variabel yang ada dengan parameter referensi yang mempunyai dimensi samam (Le Quere,1990), seperti dibawah ini : x  *

P  P x y u v t * * * * * , y  , p  ,u  ,v  ,t  Lr Lr Vr Vr tr Lr

(L.1.1)

dimana : Lr = H

tr

= (H2/)Ra-0.5

Vr = (/H)Ra0.5

Ra

= (gTH3)/()

 = (T-Tr)/(T2-T1)

Tr

= (T1+T2)/2

1. Non-dimensionalisasi Persamaan Kontinyuitas

u v  0 x y

(L1.2)

Persamaan (L1.1) dimasukkan ke dalam persamaan (L1.2) : Vr u Vr v  0 * Lr x Lr y * *

u

*

x

*



*

v

*

y

*

0

(L1.3)

(L1.4)

Tanda (*) pada persamaan (L1.4) dihilangkan sehingga diperoleh :

u v  0 x y

(L1.5)

39

2. Non-dimensionalisasi Persamaan Navier-Stokes Persamaan momentum arah x

  2u  2u  u u u 1 p u v    2  2  t x y  x  x y 

(L1.6)

Persamaan (L1.1) dimasukkan kepersamaan (L1.6) : * 2 * * 2 *   2u *  2u *  Vr u Vr  * u 1  Vr  p Vr * u    u v   2  *2  *2  * tr t * Lr  x *  Lr x * y  Lr  x y 

(L1.7)

2 * 2 * * 2 * 0.5 2 * 2 *  Ra u  Ra  * u  Ra p Ra   u  u  * u    3 u *  v *    3   *2  (L1.8) 3 * * 3 *2 H t H  x y  H x H  x y 

u t

*

*

u

*

u

*

x

*

v

*

u

*

y

*



p

*

x

*



 Ra

0.5

  2u *  2u *   *2  * 2   x y  

(L1.9)

Tanda (*) pada persamaan (L1.9) dihilangkan sehingga diperoleh : 2 2 u u u p Pr   u  u  u v    0.5  2  2  t x y x Ra  x y 

(L1.10)

Persamaan momentum arah y

  2v  2v  v v v 1 p u v     2  2   g T  Tr   x t x y  y y  

(L1.11)

Persamaan (L1.1) dimasukkan kepersamaan (L1.11) : * 2 * * 2 *   2v*  2v*  Vr v Vr  * v 1  Vr  p Vr * v    v v   2  *2  *2  * tr t * Lr  x *  Lr x * y  Lr  x y 

+ gT

(L1.12)

2 * 2 * * 2 * 0.5   2v*  2v*   Ra v  Ra  * v  Ra  p  Ra * v   * 2  *2   3  u *  v *    3  3 * * 3 H t H  x y  H y H  x y 

40

+ gT

(L1.13)

2 * 2 *    v  v  gTH 3 1  u v  *   *2  +  (L1.14) * * * 0.5 *2 v Ra  t x y y Ra  x y 

v

*

*

v

*

*

v

*

p

*

Tanda (*) pada persamaan (L1.14) dihilangkan sehingga diperoleh : 2 2 v v v p Pr   v  v   u  v    0.5  2  2   Pr t x y y Ra  x y 

(L1.15)

3. Non-dimensionalisasi Persamaan Energi 2  2  T T T  T  T  u v   2  2  t x y  x y   

(L1.16)

Persamaan (L1.1) dimasukkan kepersamaan (L1.16) :

T

2

2  2     T 1  Vr T 2  T 1  *   T 2  T 1       *    v  * 2  *2  *2  (L1.17)  u * *  tr Lr  t Lr x y y   x  

Ra H

0.5



2

t

*



Ra H

0.5

2

2  2          *    *  v  2  *2  *2   u * *   y  y  H  x  x 

2   2 1     u v   * * 0.5  *2 *2  *  t  y Ra  x y  x 



*



*



Tanda (*) pada persamaan (L1.19) dihilangkan sehingga diperoleh : 2   2    1      u v   t x y Ra 0.5  x 2 y 2   

(L1.18)

(L1.19)

41

Lampiran 2. Skema kompak beda-hingga SKEMA KOMPAK BEDA-HINGGA Skema kompak beda-hingga dapat diturunkan dengan menggunakan formula PADE (Hirsch.1961). Beberapa operator diferensial yang digunakan untuk penurunan skema kompak beda hingga adalah seperti berikut :  adalah operator beda tengah :

ui  ui 1 / 2  ui 1 / 2

(L2.1)

 adalah operator beda tengah :

 ui 

1 ui 1  ui 1  2

(L2.2)

E adalah operator pergeseran :

Eu i  ui 1

(L2.3)

Hubungan antara operator-operator tersebut adalah seperti berikut :





1 1 EE 2



  E 2 E 2

(L2.4)

1

(L2.5)

  

(L2.6)

Diskritisasi turunan pertama uxi 

1 u i 4  Ox  2 x  1 6 i

(L2.7)

 2 1  uxi   u i  Ox 4  6  x 







(L2.8)



1 1 1 1 4 E  4  E uxi  E  E ui  Ox  6 2x

(L2.9)

42

1 1 3 ui 1  ui 1   Ox 4 uxi 1  uxi  uxi 1  4 4 4x

(L2.10)

Diskritisasi turunan kedua 1  ui 4 uxx i   Ox  2 x  1 12 2

2  1    12 



(L2.11)

2  uxxi   u i  Ox 4  x  





E 2 E 1 1 E  10  E uxxi  12 x 2

1

u

i

(L2.12)

 Ox 

4

1 1 6 ui 1  2ui  ui 1   Ox 4 uxxi 1  uxxi  uxxi 1  2 10 10 5x 

(L2.13)

(L2.14)

43

Lampiran 3. Kurva Konvergensi Kurva Konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 1. Kurva konvergensi untuk Ra=106 sirip tunggal Kurva Konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 2. Kurva konvergensi untuk Ra=107 sirip tunggal

Kurva Konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 3. Kurva konvergensi untuk Ra=106 pada 2 sirip

44

Kurva Konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 4. Kurva konvergensi untuk Ra=107 pada 2 sirip Kurva Konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 5. Kurva konvergensi untuk Ra=106 pada 3 sirip Kurva Konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 6. Kurva konvergensi untuk Ra=107 pada 3 sirip

45

kurva konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 7. Kurva konvergensi untuk Ra=106 pada 4 sirip

Kurva Konvergensi

Log10(du/dx+dv/dy)

0 -0,5 -1 -1,5 -2 -2,5 -3 0

2000

4000

6000

8000

10000

Jumlah Iterasi

Gambar 8. Kurva konvergensi untuk Ra=107 pada 4 sirip

46

Lampiran 4. Numerical Wavenumber

NUMERICAL WAVENUMBER Untuk mengetahui akurasi berbagai skema beda hingga dapat dilakukan dengan analisa numerical wavwnumber. Perhitungan numerical wavenumber beberapa skema beda hingga dijelaskan di bawah ini. Sebelumnya kita definisikan variabel  dalam mode Fourier adalah :

  t e

ikx

(L4.1)

(t ) adalah koefisien Fourier darai , i=  1 dan k adalah wavenumber.

 i 1  e

ik ( x  x )

 i 1  e

ik ( x  x )

(L4.2) (L4.3)

ik ( x  x )       ik e  x  i 1

(L4.4)

ik ( x  x )       ik e  x  i 1

(L4.5)

Skema beda tengah orde-dua   i 1   i 1  x 2x

ik e

ikx

kx 



e

e

ikx

(L4.6)

ik ( x  x )

e 2i

 e 2x

ik ( x  x )

(L4.7)

 ikx

kx  sinkx 

(L4.8) (L4.9)

47

Skema beda tengah orde-empat   i 2  8 i 1  8 i 1   i  2  x 12x

ik e

ik 

ikx

e



e

 i 2 kx

ik ( x  2 x )

 8e

 8e

 ikx

 8e 12x

(L4.10)

ik ( x  x )

ikx

 8e 12x

e

ik ( x  x )

ik ( x  2 x )

(L4.11)

i 2 kx

(L4.12)

i 2 kx  i 2 kx ikx  ikx 1 e e e e kx    8 6  2i 2i

kx 

 e

   

(L4.13)

1  sin2kx   8 sinkx  6

(L4.14)

Skema kompak beda hingga

 i 1   i   i 1  '

'

ik e

ik ( x  x )

'

 ik e

ikx

a  i 1     b  i 2   i 2  2x 4x  ik e

ik ( x  x )





(L4.15)



a ik ( x  x ) ik ( x  x ) + e  e 2x



b ik ( x  2 x ) ik ( x  2 x ) e  e 4x









a b ikx ikx ) i 2 kx i 2 kx e e  e e 2x 4x ikx  ikx ikx  ikx e e   bi  e i 2 kx  e i 2 kx  ai  e  e      2ik  ik      2x   2  x 2 i 2 i       a b k 2 coskx   1  sin kx   sin 2kx  x 2x b a sin kx   sin 2kx  2 kx  2 cos kx   1

ike

ikx

 ik  ike

ikx





(L4.16)

(L4.17) (L4.18) (L4.19)

(L4.20)

Kurva numerical wavenumber untuk berbagai skema beda hingga ditampilkan dalam gambar 2.4.

48

Lampiran 5. Program Konveksi Alami Pada Sirip Tunggal C------------------------------------------------------------------------C Program penyelesaian konveksi alami pada sirip tunggal C------------------------------------------------------------------------PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 OPEN(8,FILE='C:\matlab6p1\work\temp1') OPEN(7,FILE='C:\matlab6p1\work\pv1') OPEN(4,FILE='C:\matlab6p1\work\temp') OPEN(3,FILE='C:\matlab6p1\work\pv') OPEN(2,FILE='C:\matlab6p1\work\num') OPEN(1,FILE='C:\matlab6p1\work\div') call awal do k=1,nt ck=0. call rkv do i=1,nx do j=1,ny ck=ck+abs(ux(i,j)+vy(i,j)) end do end do if(ck.gt.0) then WRITE(*,*)k,LOG10(ck/nx/ny) WRITE(1,*)k,LOG10(ck/nx/ny) endif end do call hasil stop end C---------------------------------------C Syarat awal dan syarat batas C---------------------------------------subroutine awal PARAMETER(m=500,n=500) CHARACTER mul*2 COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 WRITE(*,*)' t= ' READ(*,*)tt

49

nx=201 ny=201 i1=1 i2=nx j1=1 j2=ny dx=1./200. dy=1./200. dt=0.0025 nt=tt/dt+1 pr=0.71 ra=1000000. WRITE(*,*)' Dari awal ?' READ(*,'(a)')mul IF(mul.eq.'y')then do i=1,nx do j=1,ny x(i,j)=(i-1)*dx y(i,j)=(j-1)*dy u(i,j)=0. v(i,j)=0. p(i,j)=0. o(i,j)=0. enddo enddo do j=1,ny o(1,j)=0.5 o(nx,j)=-0.5 enddo else do j=1,ny do i=1,nx READ(7,*)x(i,j),y(i,j),u(i,j),v(i,j),p(i,j) READ(8,*)x(i,j),y(i,j),o(i,j) end do end do end if return end C-------------------------------------------------------C Penyelesaian persamaan Navier - Stokes C-------------------------------------------------------subroutine rkv PARAMETER(m=500,n=500)

50

COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), vx(m,n),vy(m,n),px(m,n),py(m,n), uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 DIMENSION h1(m,n,6),h2(m,n,6),s1(m,n,6),s2(m,n,6),a(5),b(5) DATA a(1),a(2),a(3),a(4),a(5)/0.,-0.41789047,-1.19215169, 1 -1.69778469,-1.51418344/ DATA b(1),b(2),b(3),b(4),b(5)/0.14965902,0.37921031,0.82295502, 1 0.69945045,0.15305724/ 1 1 1

do i=2,nx-1 do j=2,ny-1 s1(i,j,1)=u(i,j) s2(i,j,1)=v(i,j) end do end do do k=1,5 call derv1(u,ux,uy) call derv1(v,vx,vy) call derv2(u,uxx,uyy) call derv2(v,vxx,vyy) call derp1(p,px,py) do i=2,nx-1 do j=2,ny-1 h1(i,j,k)=-u(i,j)*ux(i,j)-v(i,j)*uy(i,j)-px(i,j) 1 +(uxx(i,j)+uyy(i,j))*pr/SQRT(ra) h2(i,j,k)=-u(i,j)*vx(i,j)-v(i,j)*vy(i,j)-py(i,j) 1 +(vxx(i,j)+vyy(i,j))*pr/SQRT(ra)+pr*o(i,j) end do end do do i=2,nx-1 do j=2,ny-1 IF(i.le.101.and.j.ge.98.and.j.le.102)then u(i,j)=0. v(i,j)=0 else h1(i,j,k)=h1(i,j,k)+a(k)*h1(i,j,k-1) h2(i,j,k)=h2(i,j,k)+a(k)*h2(i,j,k-1) s1(i,j,k+1)=s1(i,j,k)+b(k)*dt*h1(i,j,k) s2(i,j,k+1)=s2(i,j,k)+b(k)*dt*h2(i,j,k) u(i,j)=s1(i,j,k+1) v(i,j)=s2(i,j,k+1) endif end do end do

51

call derv1(u,ux,uy) call derv1(v,vx,vy) c------------------------------------------------------------------c Tekananan dihitung dengan artifisial kompresibiliti c------------------------------------------------------------------do i=2,nx-1 do j=2,ny-1 IF(i.lt.101.and.j.gt.98.and.j.lt.102)then p(i,j)=0 else p(i,j)=p(i,j)-(ux(i,j)+vy(i,j))*dt*0.5 endif end do end do do i=1,nx IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny else j1=1 j2=ny endif p(i,j1)=(48*p(i,j1+1)-36*p(i,j1+2)+16*p(i,j1+3)-3*p(i,j1+4))/25 p(i,j2)=(48*p(i,j2-1)-36*p(i,j2-2)+16*p(i,j2-3)-3*p(i,j2-4))/25 p(i,j3)=(48*p(i,j3+1)-36*p(i,j3+2)+16*p(i,j3+3)-3*p(i,j3+4))/25 p(i,j4)=(48*p(i,j4-1)-36*p(i,j4-2)+16*p(i,j4-3)-3*p(i,j4-4))/25 end do i1=1 i2=nx j1=1 j2=ny do j=1,ny IF(j.ge.98.and.j.le.102)then i1=101 else i1=1 i2=nx endif p(i1,j)=(48*p(i1+1,j)-36*p(i1+2,j)+16*p(i1+3,j)-3*p(i1+4,j))/25 p(i2,j)=(48*p(i2-1,j)-36*p(i2-2,j)+16*p(i2-3,j)-3*p(i2-4,j))/25 enddo call rkt c-----------------------------------------------------------

52

end do return end C-------------------------------------------C Penyelesaian persamaan energi C------------------------------------------subroutine rkt PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 DIMENSION h(m,n,6),s(m,n,6),a(5),b(5) DATA a(1),a(2),a(3),a(4),a(5)/0.,-0.41789047,-1.19215169, 1 -1.69778469,-1.51418344/ DATA b(1),b(2),b(3),b(4),b(5)/0.14965902,0.37921031,0.82295502, 1 0.69945045,0.15305724/ do i=2,nx-1 do j=2,ny-1 s(i,j,1)=o(i,j) end do end do do k=1,5 call dero1(o,ox,oy) call derv2(o,oxx,oyy) do i=2,nx-1 do j=2,ny-1 IF(i.lt.101.and.j.gt.98.and.j.lt.102)then o(i,j)=0. else h(i,j,k)=-u(i,j)*ox(i,j)-v(i,j)*oy(i,j)+ 1 (oxx(i,j)+oyy(i,j))/SQRT(ra)+a(k)*h(i,j,k-1) s(i,j,k+1)=s(i,j,k)+b(k)*dt*h(i,j,k) o(i,j)=s(i,j,k+1) endif end do end do do i=1,nx IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny

53

o(i,j1)=(48*o(i,j1+1)-36*o(i,j1+2)+16*o(i,j1+3)-3*o(i,j1+4))/25 o(i,j2)=0.5 o(i,j3)=0.5 o(i,j4)=(48*o(i,j4-1)-36*o(i,j4-2)+16*o(i,j4-3)-3*o(i,j4-4))/25 else j1=1 j2=ny o(i,j1)=(48*o(i,j1+1)-36*o(i,j1+2)+16*o(i,j1+3)-3*o(i,j1+4))/25 o(i,j2)=(48*o(i,j2-1)-36*o(i,j2-2)+16*o(i,j2-3)-3*o(i,j2-4))/25 endif end do do j=1,ny IF(j.ge.98.and.j.le.102)then i1=101 else i1=1 i2=nx endif o(i1,j)=0.5 end do i1=1 i2=nx j1=1 j2=ny do i=1,nx do j=1,ny IF(i.lt.101.and.j.gt.98.and.j.lt.102)then o(i,j)=0. endif end do end do end do return end C-----------------------------------------C Turunan pertama variabel kecepatan C-----------------------------------------subroutine derv1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.98.and.j.le.102)then

54

i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1

55

210 qy(i,j)=d(j) else j1=1 j2=ny endif qy(i,j1)=(-3*q(i,j1+4)+16*q(i,j1+3)-36*q(i,j1+2)+48*q(i,j1+1) 1 -25*q(i,j1))/12/dy qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end C-----------------------------------------C turunan pertama variabel temperatur C-----------------------------------------subroutine dero1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.98.and.j.le.102)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx

56

qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny qy(i,j1)=0. qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=0. do 110 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 210j=j1+1,j2-1 210 qy(i,j)=d(j) do 111 j=j3+1,j4-1 a(j)=alp b(j)=1.

57

c(j)=alp 111 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 211 j=j3+1,j4-1 211 qy(i,j)=d(j) else j1=1 j2=ny qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) endif enddo return end C---------------------------------------C Turunan pertama variabel tekanan C---------------------------------------subroutine derp1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.98.and.j.le.102)then

58

i1=101 else i1=1 i2=nx endif qx(i1,j)=0. qx(i2,j)=0. do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny qy(i,j3)=0. qy(i,j4)=0. do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) else j1=1 j2=ny

59

endif qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end C-------------------------------------------------------C Turunan kedua variabel kecepatan dan temperatur C-------------------------------------------------------subroutine derv2(q,qxx,qyy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qxx(500,500),qyy(500,500) alp=0.1 ak=1.2 do 10 j=2,ny-1 IF(j.ge.98.and.j.le.102)then i1=101 else i1=1 i2=nx endif s1=(13*q(i1,j)-27*q(i1+1,j)+15*q(i1+2,j)-1*q(i1+3,j))/dx/dx s2=(-1*q(i2-3,j)+15*q(i2-2,j)-27*q(i2-1,j)+13*q(i2,j))/dx/dx do 20 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 20 d(i)=ak*(q(i+1,j)-2*q(i,j)+q(i-1,j))/dx/dx b(i1+1)=b(i1+1)-11*alp d(i1+1)=d(i1+1)-s1*alp

60

b(i2-1)=b(i2-1)-11*alp d(i2-1)=d(i2-1)-s2*alp l1=i1+2 l2=i2-1 call tridi do 30 i=i1+1,i2-1 30 qxx(i,j)=d(i) qxx(i1,j)=s1-11*qxx(i1+1,j) qxx(i2,j)=s2-11*qxx(i2-1,j) 10 continue do 40 i=1,nx-1 IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny s3=(13*q(i,j3)-27*q(i,j3+1)+15*q(i,j3+2)-q(i,j3+3))/dy/dy s4=(-q(i,j4-3)+15*q(i,j4-2)-27*q(i,j4-1)+13*q(i,j4))/dy/dy do 51 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 51 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j3+1)=b(j3+1)-11*alp d(j3+1)=d(j3+1)-s3*alp b(j4-1)=b(j4-1)-11*alp d(j4-1)=d(j4-1)-s4*alp l1=j3+2 l2=j4-1 call tridi do 61 j=j3+1,j4-1 61 qyy(i,j)=d(j) qyy(i,j3)=s1-11*qyy(i,j3+1) qyy(i,j4)=s2-11*qyy(i,j4-1) else j1=1 j2=ny endif s1=(13*q(i,j1)-27*q(i,j1+1)+15*q(i,j1+2)-q(i,j1+3))/dy/dy s2=(-q(i,j2-3)+15*q(i,j2-2)-27*q(i,j2-1)+13*q(i,j2))/dy/dy do 50 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 50 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j1+1)=b(j1+1)-11*alp d(j1+1)=d(j1+1)-s1*alp

61

b(j2-1)=b(j2-1)-11*alp d(j2-1)=d(j2-1)-s2*alp l1=j1+2 l2=j2-1 call tridi do 60 j=j1+1,j2-1 60 qyy(i,j)=d(j) qyy(i,j1)=s1-11*qyy(i,j1+1) qyy(i,j2)=s2-11*qyy(i,j2-1) 40 continue return end C------------------------C Algoritma Thomas C------------------------SUBROUTINE TRIDI parameter (m=500) COMMON /AA4/A(m),B(m),C(m),D(m),L1,l2 DO 1 I=L1,L2 RT=-A(I)/B(I-1) B(I)=B(I)+RT*C(I-1) 1 D(I)=D(I)+RT*D(I-1) D(L2)=D(L2)/B(L2) DO 2 I=L2-1,L1-1,-1 2 D(I)=(D(I)-C(I)*D(I+1))/B(I) RETURN END C-------------------------C Hasil perhitungan C-------------------------subroutine hasil PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 WRITE(2,*)NX WRITE(2,*)NY do 10 j=1,ny do 10 i=1,nx WRITE(3,*)x(i,j),y(i,j),u(i,j),v(i,j)!,p(i,j) 10 WRITE(4,*)x(i,j),y(i,j),o(i,j) return end

62

Lampiran 6. Program Konveksi Alami Pada 2 Sirip C------------------------------------------------------------------C Program penyelesaian konveksi alami pada 2 sirip C------------------------------------------------------------------PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 OPEN(8,FILE='C:\matlab6p1\work\temp1') OPEN(7,FILE='C:\matlab6p1\work\pv1') OPEN(4,FILE='C:\matlab6p1\work\temp') OPEN(3,FILE='C:\matlab6p1\work\pv') OPEN(2,FILE='C:\matlab6p1\work\num') OPEN(1,FILE='C:\matlab6p1\work\div') call awal do k=1,nt ck=0. call rkv do i=1,nx do j=1,ny ck=ck+abs(ux(i,j)+vy(i,j)) end do end do if(ck.gt.0) then WRITE(*,*)k,LOG10(ck/nx/ny) WRITE(1,*)k,LOG10(ck/nx/ny) endif end do call hasil stop end C------------------------------------C Syarat awal dan syarat batas C------------------------------------subroutine awal PARAMETER(m=500,n=500) CHARACTER mul*2 COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 WRITE(*,*)' t= ' READ(*,*)tt nx=201 ny=201

63

i1=1 i2=nx j1=1 j2=ny dx=1./200. dy=1./200. dt=0.0025 nt=tt/dt+1 pr=0.71 ra=10000000. WRITE(*,*)' Dari awal ?' READ(*,'(a)')mul IF(mul.eq.'y')then do i=1,nx do j=1,ny x(i,j)=(i-1)*dx y(i,j)=(j-1)*dy u(i,j)=0. v(i,j)=0. p(i,j)=0. o(i,j)=0. enddo enddo do j=1,ny o(1,j)=0.5 o(nx,j)=-0.5 enddo else do j=1,ny do i=1,nx READ(7,*)x(i,j),y(i,j),u(i,j),v(i,j),p(i,j) READ(8,*)x(i,j),y(i,j),o(i,j) end do end do end if return end C----------------------------------------------C Penyelesaian persamaan Navier - Stokes C----------------------------------------------subroutine rkv PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n),

64

1 1

uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 DIMENSION h1(m,n,6),h2(m,n,6),s1(m,n,6),s2(m,n,6),a(5),b(5) DATA a(1),a(2),a(3),a(4),a(5)/0.,-0.41789047,-1.19215169, 1 -1.69778469,-1.51418344/ DATA b(1),b(2),b(3),b(4),b(5)/0.14965902,0.37921031,0.82295502, 1 0.69945045,0.15305724/ do i=2,nx-1 do j=2,ny-1 s1(i,j,1)=u(i,j) s2(i,j,1)=v(i,j) end do end do do k=1,5 call derv1(u,ux,uy) call derv1(v,vx,vy) call derv2(u,uxx,uyy) call derv2(v,vxx,vyy) call derp1(p,px,py) do i=2,nx-1 do j=2,ny-1 h1(i,j,k)=-u(i,j)*ux(i,j)-v(i,j)*uy(i,j)-px(i,j) 1 +(uxx(i,j)+uyy(i,j))*pr/SQRT(ra) h2(i,j,k)=-u(i,j)*vx(i,j)-v(i,j)*vy(i,j)-py(i,j) 1 +(vxx(i,j)+vyy(i,j))*pr/SQRT(ra)+pr*o(i,j) end do end do do i=2,nx-1 do j=2,ny-1 IF(i.le.101.and.j.ge.65.and.j.le.69)then u(i,j)=0. v(i,j)=0 ELSEIF(i.le.101.and.j.ge.131.and.j.le.135)then u(i,j)=0. v(i,j)=0 else h1(i,j,k)=h1(i,j,k)+a(k)*h1(i,j,k-1) h2(i,j,k)=h2(i,j,k)+a(k)*h2(i,j,k-1) s1(i,j,k+1)=s1(i,j,k)+b(k)*dt*h1(i,j,k) s2(i,j,k+1)=s2(i,j,k)+b(k)*dt*h2(i,j,k) u(i,j)=s1(i,j,k+1) v(i,j)=s2(i,j,k+1) endif end do

65

end do call derv1(u,ux,uy) call derv1(v,vx,vy) c----------------------------------------------------------c Tekananan dihitung dengan artifisial kompresibiliti c----------------------------------------------------------do i=2,nx-1 do j=2,ny-1 IF(i.lt.101.and.j.gt.65.and.j.lt.69)then p(i,j)=0 ELSEIF(i.lt.101.and.j.gt.131.and.j.lt.135)then p(i,j)=0 else p(i,j)=p(i,j)-(ux(i,j)+vy(i,j))*dt*0.5 endif end do end do do i=1,nx IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135 j6=ny else j1=1 j2=ny endif p(i,j1)=(48*p(i,j1+1)-36*p(i,j1+2)+16*p(i,j1+3)-3*p(i,j1+4))/25 p(i,j2)=(48*p(i,j2-1)-36*p(i,j2-2)+16*p(i,j2-3)-3*p(i,j2-4))/25 p(i,j3)=(48*p(i,j3+1)-36*p(i,j3+2)+16*p(i,j3+3)-3*p(i,j3+4))/25 p(i,j4)=(48*p(i,j4-1)-36*p(i,j4-2)+16*p(i,j4-3)-3*p(i,j4-4))/25 p(i,j5)=(48*p(i,j5+1)-36*p(i,j5+2)+16*p(i,j5+3)-3*p(i,j5+4))/25 p(i,j6)=(48*p(i,j6-1)-36*p(i,j6-2)+16*p(i,j6-3)-3*p(i,j6-4))/25 end do i1=1 i2=nx j1=1 j2=ny do j=1,ny IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else i1=1

66

i2=nx endif p(i1,j)=(48*p(i1+1,j)-36*p(i1+2,j)+16*p(i1+3,j)-3*p(i1+4,j))/25 p(i2,j)=(48*p(i2-1,j)-36*p(i2-2,j)+16*p(i2-3,j)-3*p(i2-4,j))/25 enddo call rkt c----------------------------------------------------------end do return end C---------------------------------------C Penyelesaian persamaan energi C---------------------------------------subroutine rkt PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 DIMENSION h(m,n,6),s(m,n,6),a(5),b(5) DATA a(1),a(2),a(3),a(4),a(5)/0.,-0.41789047,-1.19215169, 1 -1.69778469,-1.51418344/ DATA b(1),b(2),b(3),b(4),b(5)/0.14965902,0.37921031,0.82295502, 1 0.69945045,0.15305724/ do i=2,nx-1 do j=2,ny-1 s(i,j,1)=o(i,j) end do end do do k=1,5 call dero1(o,ox,oy) call derv2(o,oxx,oyy) do i=2,nx-1 do j=2,ny-1 IF(i.lt.101.and.j.gt.65.and.j.lt.69)then o(i,j)=0. ELSEIF(i.lt.101.and.j.gt.131.and.j.lt.135)then o(i,j)=0. else h(i,j,k)=-u(i,j)*ox(i,j)-v(i,j)*oy(i,j)+ 1 (oxx(i,j)+oyy(i,j))/SQRT(ra)+a(k)*h(i,j,k-1)

67

s(i,j,k+1)=s(i,j,k)+b(k)*dt*h(i,j,k) o(i,j)=s(i,j,k+1) endif end do end do do i=1,nx IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135 j6=ny o(i,j1)=(48*o(i,j1+1)-36*o(i,j1+2)+16*o(i,j1+3)-3*o(i,j1+4))/25 o(i,j2)=0.5 o(i,j3)=0.5 o(i,j4)=0.5 o(i,j5)=0.5 o(i,j6)=(48*o(i,j6-1)-36*o(i,j6-2)+16*o(i,j6-3)-3*o(i,j6-4))/25 else j1=1 j2=ny o(i,j1)=(48*o(i,j1+1)-36*o(i,j1+2)+16*o(i,j1+3)-3*o(i,j1+4))/25 o(i,j2)=(48*o(i,j2-1)-36*o(i,j2-2)+16*o(i,j2-3)-3*o(i,j2-4))/25 endif end do do j=1,ny IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else i1=1 i2=nx endif o(i1,j)=0.5 end do i1=1 i2=nx j1=1 j2=ny do i=1,nx do j=1,ny IF(i.lt.101.and.j.gt.65.and.j.lt.69)then o(i,j)=0. ELSEIF(i.lt.101.and.j.gt.131.and.j.lt.135)then

68

o(i,j)=0. endif end do end do end do return end C-----------------------------------------C Turunan pertama variabel kecepatan C-----------------------------------------subroutine derv1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo

69

do i=1,nx IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135 j6=ny qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1) 1 +25*q(i,j6))/12/dy do 222 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 222 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 322 j=j5+1,j6-1 322 qy(i,j)=d(j) else j1=1

70

j2=ny endif qy(i,j1)=(-3*q(i,j1+4)+16*q(i,j1+3)-36*q(i,j1+2)+48*q(i,j1+1) 1 -25*q(i,j1))/12/dy qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end C-----------------------------------------C turunan pertama variabel temperatur C-----------------------------------------subroutine dero1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx

71

qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135 j6=ny qy(i,j1)=0. qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=0. do 110 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1

72

call tridi do 210j=j1+1,j2-1 210 qy(i,j)=d(j) do 111 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 111 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 211 j=j3+1,j4-1 211 qy(i,j)=d(j) do 345 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 345 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 789 j=j5+1,j6-1 789 qy(i,j)=d(j) else j1=1 j2=ny qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2

73

l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) endif enddo return end C---------------------------------------C Turunan pertama variabel tekanan C---------------------------------------subroutine derp1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else i1=1 i2=nx endif qx(i1,j)=0. qx(i2,j)=0. do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx

74

IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135 j6=ny qy(i,j3)=0. qy(i,j4)=0. do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) qy(i,j5)=0. qy(i,j6)=0. do 222 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 222 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 322 j=j5+1,j6-1 322 qy(i,j)=d(j) else j1=1 j2=ny endif qy(i,j1)=0.

75

qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end C-------------------------------------------------------C Turunan kedua variabel kecepatan dan temperatur C-------------------------------------------------------subroutine derv2(q,qxx,qyy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qxx(500,500),qyy(500,500) alp=0.1 ak=1.2 do 10 j=2,ny-1 IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else i1=1 i2=nx endif s1=(13*q(i1,j)-27*q(i1+1,j)+15*q(i1+2,j)-1*q(i1+3,j))/dx/dx s2=(-1*q(i2-3,j)+15*q(i2-2,j)-27*q(i2-1,j)+13*q(i2,j))/dx/dx do 20 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 20 d(i)=ak*(q(i+1,j)-2*q(i,j)+q(i-1,j))/dx/dx b(i1+1)=b(i1+1)-11*alp d(i1+1)=d(i1+1)-s1*alp

76

b(i2-1)=b(i2-1)-11*alp d(i2-1)=d(i2-1)-s2*alp l1=i1+2 l2=i2-1 call tridi do 30 i=i1+1,i2-1 30 qxx(i,j)=d(i) qxx(i1,j)=s1-11*qxx(i1+1,j) qxx(i2,j)=s2-11*qxx(i2-1,j) 10 continue do 40 i=1,nx-1 IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135 j6=ny s3=(13*q(i,j3)-27*q(i,j3+1)+15*q(i,j3+2)-q(i,j3+3))/dy/dy s4=(-q(i,j4-3)+15*q(i,j4-2)-27*q(i,j4-1)+13*q(i,j4))/dy/dy do 51 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 51 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j3+1)=b(j3+1)-11*alp d(j3+1)=d(j3+1)-s3*alp b(j4-1)=b(j4-1)-11*alp d(j4-1)=d(j4-1)-s4*alp l1=j3+2 l2=j4-1 call tridi do 61 j=j3+1,j4-1 61 qyy(i,j)=d(j) qyy(i,j3)=s1-11*qyy(i,j3+1) qyy(i,j4)=s2-11*qyy(i,j4-1) s5=(13*q(i,j5)-27*q(i,j5+1)+15*q(i,j5+2)-q(i,j5+3))/dy/dy s6=(-q(i,j6-3)+15*q(i,j6-2)-27*q(i,j6-1)+13*q(i,j6))/dy/dy do 71 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 71 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j5+1)=b(j5+1)-11*alp d(j5+1)=d(j5+1)-s5*alp b(j6-1)=b(j6-1)-11*alp

77

d(j6-1)=d(j6-1)-s6*alp l1=j5+2 l2=j6-1 call tridi do 81 j=j5+1,j6-1 81 qyy(i,j)=d(j) qyy(i,j5)=s1-11*qyy(i,j5+1) qyy(i,j6)=s2-11*qyy(i,j6-1) else j1=1 j2=ny endif s1=(13*q(i,j1)-27*q(i,j1+1)+15*q(i,j1+2)-q(i,j1+3))/dy/dy s2=(-q(i,j2-3)+15*q(i,j2-2)-27*q(i,j2-1)+13*q(i,j2))/dy/dy do 50 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 50 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j1+1)=b(j1+1)-11*alp d(j1+1)=d(j1+1)-s1*alp b(j2-1)=b(j2-1)-11*alp d(j2-1)=d(j2-1)-s2*alp l1=j1+2 l2=j2-1 call tridi do 60 j=j1+1,j2-1 60 qyy(i,j)=d(j) qyy(i,j1)=s1-11*qyy(i,j1+1) qyy(i,j2)=s2-11*qyy(i,j2-1) 40 continue return end C------------------------C Algoritma Thomas C-------------------------

1

SUBROUTINE TRIDI parameter (m=500) COMMON /AA4/A(m),B(m),C(m),D(m),L1,l2 DO 1 I=L1,L2 RT=-A(I)/B(I-1) B(I)=B(I)+RT*C(I-1) D(I)=D(I)+RT*D(I-1) D(L2)=D(L2)/B(L2) DO 2 I=L2-1,L1-1,-1

78

2

D(I)=(D(I)-C(I)*D(I+1))/B(I) RETURN END C-------------------------C Hasil perhitungan C-------------------------subroutine hasil PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 WRITE(2,*)NX WRITE(2,*)NY do 10 j=1,ny do 10 i=1,nx WRITE(3,*)x(i,j),y(i,j),u(i,j),v(i,j)!,p(i,j) 10 WRITE(4,*)x(i,j),y(i,j),o(i,j) return end

79

Lampiran 7. Program Konveksi Alami Pada 3 Sirip C------------------------------------------------------------------C Program penyelesaian konveksi alami pada 3 sirip C------------------------------------------------------------------PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 OPEN(8,FILE='C:\matlab6p1\work\temp1') OPEN(7,FILE='C:\matlab6p1\work\pv1') OPEN(4,FILE='C:\matlab6p1\work\temp') OPEN(3,FILE='C:\matlab6p1\work\pv') OPEN(2,FILE='C:\matlab6p1\work\num') OPEN(1,FILE='C:\matlab6p1\work\div') call awal do k=1,nt ck=0. call rkv do i=1,nx do j=1,ny ck=ck+abs(ux(i,j)+vy(i,j)) end do end do if(ck.gt.0) then WRITE(*,*)k,LOG10(ck/nx/ny) WRITE(1,*)k,LOG10(ck/nx/ny) endif end do call hasil stop end C---------------------------------------C Syarat awal dan syarat batas C---------------------------------------subroutine awal PARAMETER(m=500,n=500) CHARACTER mul*2 COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 WRITE(*,*)' t= '

80

READ(*,*)tt nx=201 ny=201 i1=1 i2=nx j1=1 j2=ny dx=1./200. dy=1./200. dt=0.0025 nt=tt/dt+1 pr=0.71 ra=10000000. WRITE(*,*)' Dari awal ?' READ(*,'(a)')mul IF(mul.eq.'y')then do i=1,nx do j=1,ny x(i,j)=(i-1)*dx y(i,j)=(j-1)*dy u(i,j)=0. v(i,j)=0. p(i,j)=0. o(i,j)=0. enddo enddo do j=1,ny o(1,j)=0.5 o(nx,j)=-0.5 enddo else do j=1,ny do i=1,nx READ(7,*)x(i,j),y(i,j),u(i,j),v(i,j),p(i,j) READ(8,*)x(i,j),y(i,j),o(i,j) end do end do end if return end C-------------------------------------------------------C Penyelesaian persamaan Navier - Stokes C-------------------------------------------------------subroutine rkv PARAMETER(m=500,n=500)

81

COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), vx(m,n),vy(m,n),px(m,n),py(m,n), uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 DIMENSION h1(m,n,6),h2(m,n,6),s1(m,n,6),s2(m,n,6),a(5),b(5) DATA a(1),a(2),a(3),a(4),a(5)/0.,-0.75789047,-1.19215169, 1 -1.69778469,-1.51418344/ DATA b(1),b(2),b(3),b(4),b(5)/0.14965902,0.37921031,0.82295502, 1 0.69945045,0.15305724/ 1 1 1

do i=2,nx-1 do j=2,ny-1 s1(i,j,1)=u(i,j) s2(i,j,1)=v(i,j) end do end do do k=1,5 call derv1(u,ux,uy) call derv1(v,vx,vy) call derv2(u,uxx,uyy) call derv2(v,vxx,vyy) call derp1(p,px,py) do i=2,nx-1 do j=2,ny-1 h1(i,j,k)=-u(i,j)*ux(i,j)-v(i,j)*uy(i,j)-px(i,j) 1 +(uxx(i,j)+uyy(i,j))*pr/SQRT(ra) h2(i,j,k)=-u(i,j)*vx(i,j)-v(i,j)*vy(i,j)-py(i,j) 1 +(vxx(i,j)+vyy(i,j))*pr/SQRT(ra)+pr*o(i,j) end do end do do i=2,nx-1 do j=2,ny-1 IF(i.le.101.and.j.ge.48.and.j.le.52)then u(i,j)=0. v(i,j)=0 ELSEIF (i.le.101.and.j.ge.98.and.j.le.102)then u(i,j)=0. v(i,j)=0 ELSEIF(i.le.101.and.j.ge.148.and.j.le.152)then u(i,j)=0. v(i,j)=0 else h1(i,j,k)=h1(i,j,k)+a(k)*h1(i,j,k-1) h2(i,j,k)=h2(i,j,k)+a(k)*h2(i,j,k-1) s1(i,j,k+1)=s1(i,j,k)+b(k)*dt*h1(i,j,k)

82

s2(i,j,k+1)=s2(i,j,k)+b(k)*dt*h2(i,j,k) u(i,j)=s1(i,j,k+1) v(i,j)=s2(i,j,k+1) endif end do end do call derv1(u,ux,uy) call derv1(v,vx,vy) c----------------------------------------------------------c Tekananan dihitung dengan artifisial kompresibiliti c----------------------------------------------------------do i=2,nx-1 do j=2,ny-1 IF(i.lt.101.and.j.gt.48.and.j.lt.52)then p(i,j)=0 ELSEIF (i.lt.101.and.j.gt.98.and.j.lt.102)then p(i,j)=0 ELSEIF(i.lt.101.and.j.gt.148.and.j.lt.152)then p(i,j)=0 else p(i,j)=p(i,j)-(ux(i,j)+vy(i,j))*dt*0.5 endif end do end do do i=1,nx IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny else j1=1 j2=ny endif p(i,j1)=(48*p(i,j1+1)-36*p(i,j1+2)+16*p(i,j1+3)-3*p(i,j1+4))/25 p(i,j2)=(48*p(i,j2-1)-36*p(i,j2-2)+16*p(i,j2-3)-3*p(i,j2-4))/25 p(i,j3)=(48*p(i,j3+1)-36*p(i,j3+2)+16*p(i,j3+3)-3*p(i,j3+4))/25 p(i,j4)=(48*p(i,j4-1)-36*p(i,j4-2)+16*p(i,j4-3)-3*p(i,j4-4))/25 p(i,j5)=(48*p(i,j5+1)-36*p(i,j5+2)+16*p(i,j5+3)-3*p(i,j5+4))/25 p(i,j6)=(48*p(i,j6-1)-36*p(i,j6-2)+16*p(i,j6-3)-3*p(i,j6-4))/25 p(i,j7)=(48*p(i,j7+1)-36*p(i,j7+2)+16*p(i,j7+3)-3*p(i,j7+4))/25 p(i,j8)=(48*p(i,j8-1)-36*p(i,j8-2)+16*p(i,j8-3)-3*p(i,j8-4))/25

83

end do i1=1 i2=nx j1=1 j2=ny do j=1,ny IF(j.ge.48.and.j.le.52)then i1=101 elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then i1=101 else i1=1 i2=nx endif p(i1,j)=(48*p(i1+1,j)-36*p(i1+2,j)+16*p(i1+3,j)-3*p(i1+4,j))/25 p(i2,j)=(48*p(i2-1,j)-36*p(i2-2,j)+16*p(i2-3,j)-3*p(i2-4,j))/25 enddo call rkt c----------------------------------------------------------end do return end C-------------------------------------------C Penyelesaian persamaan energi C------------------------------------------subroutine rkt PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 DIMENSION h(m,n,6),s(m,n,6),a(5),b(5) DATA a(1),a(2),a(3),a(4),a(5)/0.,-0.75789047,-1.19215169, 1 -1.69778469,-1.51418344/ DATA b(1),b(2),b(3),b(4),b(5)/0.14965902,0.37921031,0.82295502, 1 0.69945045,0.15305724/ do i=2,nx-1 do j=2,ny-1 s(i,j,1)=o(i,j) end do end do

84

do k=1,5 call dero1(o,ox,oy) call derv2(o,oxx,oyy) do i=2,nx-1 do j=2,ny-1 IF(i.lt.101.and.j.gt.48.and.j.lt.52)then o(i,j)=0. elseif(i.lt.101.and.j.gt.98.and.j.lt.102)then o(i,j)=0. ELSEIF(i.lt.101.and.j.gt.148.and.j.lt.152)then o(i,j)=0. else h(i,j,k)=-u(i,j)*ox(i,j)-v(i,j)*oy(i,j)+ 1 (oxx(i,j)+oyy(i,j))/SQRT(ra)+a(k)*h(i,j,k-1) s(i,j,k+1)=s(i,j,k)+b(k)*dt*h(i,j,k) o(i,j)=s(i,j,k+1) endif end do end do do i=1,nx IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny o(i,j1)=(48*o(i,j1+1)-36*o(i,j1+2)+16*o(i,j1+3)-3*o(i,j1+4))/25 o(i,j2)=0.5 o(i,j3)=0.5 o(i,j4)=0.5 o(i,j5)=0.5 o(i,j6)=0.5 o(i,j7)=0.5 o(i,j8)=(48*o(i,j8-1)-36*o(i,j8-2)+16*o(i,j8-3)-3*o(i,j8-4))/25 else j1=1 j2=ny o(i,j1)=(48*o(i,j1+1)-36*o(i,j1+2)+16*o(i,j1+3)-3*o(i,j1+4))/25 o(i,j2)=(48*o(i,j2-1)-36*o(i,j2-2)+16*o(i,j2-3)-3*o(i,j2-4))/25 endif end do do j=1,ny

85

IF(j.ge.48.and.j.le.52)then i1=101 elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then i1=101 else i1=1 i2=nx endif o(i1,j)=0.5 end do i1=1 i2=nx j1=1 j2=ny do i=1,nx do j=1,ny IF(i.lt.101.and.j.gt.48.and.j.lt.52)then o(i,j)=0. elseif(i.lt.101.and.j.gt.98.and.j.lt.102)then o(i,j)=0. ELSEIF(i.lt.101.and.j.gt.148.and.j.lt.152)then o(i,j)=0. endif end do end do end do return end C-----------------------------------------C Turunan pertama variabel kecepatan C-----------------------------------------subroutine derv1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.48.and.j.le.52)then i1=101 elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then

86

i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0.

87

l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1) 1 +25*q(i,j6))/12/dy do 222 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 222 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 322 j=j5+1,j6-1 322 qy(i,j)=d(j) qy(i,j7)=(-3*q(i,j7+4)+16*q(i,j7+3)-36*q(i,j7+2)+48*q(i,j7+1) 1 -25*q(i,j7))/12/dy qy(i,j8)=(3*q(i,j8-4)-16*q(i,j8-3)+36*q(i,j8-2)-48*q(i,j8-1) 1 +25*q(i,j8))/12/dy do 333 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp 333 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8) c(j8-1)=0. l1=j7+2 l2=j8-1 call tridi do 433 j=j7+1,j8-1 433 qy(i,j)=d(j) else j1=1 j2=ny

88

endif qy(i,j1)=(-3*q(i,j1+4)+16*q(i,j1+3)-36*q(i,j1+2)+48*q(i,j1+1) 1 -25*q(i,j1))/12/dy qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end C-----------------------------------------C turunan pertama variabel temperatur C-----------------------------------------subroutine dero1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.48.and.j.le.52)then i1=101 elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx

89

qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny qy(i,j1)=0. qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1) 1 +25*q(i,j6))/12/dy qy(i,j7)=(-3*q(i,j7+4)+16*q(i,j7+3)-36*q(i,j7+2)+48*q(i,j7+1) 1 -25*q(i,j7))/12/dy qy(i,j8)=0. do 110 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1)

90

a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 210j=j1+1,j2-1 210 qy(i,j)=d(j) do 111 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 111 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 211 j=j3+1,j4-1 211 qy(i,j)=d(j) do 345 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 345 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 789 j=j5+1,j6-1 789 qy(i,j)=d(j) do 444 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp 444 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8)

91

c(j8-1)=0. l1=j7+2 l2=j8-1 call tridi do 544 j=j7+1,j8-1 544 qy(i,j)=d(j) else j1=1 j2=ny qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) endif enddo return end C---------------------------------------C Turunan pertama variabel tekanan C---------------------------------------subroutine derp1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.48.and.j.le.52)then i1=101 elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then

92

i1=101 else i1=1 i2=nx endif qx(i1,j)=0. qx(i2,j)=0. do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny qy(i,j3)=0. qy(i,j4)=0. do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1

93

210 qy(i,j)=d(j) qy(i,j5)=0. qy(i,j6)=0. do 222 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 222 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 322 j=j5+1,j6-1 322 qy(i,j)=d(j) qy(i,j7)=0. qy(i,j8)=0. do 555 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp 555 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8) c(j8-1)=0. l1=j7+2 l2=j8-1 call tridi do 655 j=j7+1,j8-1 655 qy(i,j)=d(j) else j1=1 j2=ny endif qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy

94

d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end C-------------------------------------------------------C Turunan kedua variabel kecepatan dan temperatur C-------------------------------------------------------subroutine derv2(q,qxx,qyy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qxx(500,500),qyy(500,500) alp=0.1 ak=1.2 do 10 j=2,ny-1 IF(j.ge.48.and.j.le.52)then i1=101 elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then i1=101 else i1=1 i2=nx endif s1=(13*q(i1,j)-27*q(i1+1,j)+15*q(i1+2,j)-1*q(i1+3,j))/dx/dx s2=(-1*q(i2-3,j)+15*q(i2-2,j)-27*q(i2-1,j)+13*q(i2,j))/dx/dx do 20 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 20 d(i)=ak*(q(i+1,j)-2*q(i,j)+q(i-1,j))/dx/dx b(i1+1)=b(i1+1)-11*alp d(i1+1)=d(i1+1)-s1*alp b(i2-1)=b(i2-1)-11*alp d(i2-1)=d(i2-1)-s2*alp l1=i1+2 l2=i2-1

95

call tridi do 30 i=i1+1,i2-1 30 qxx(i,j)=d(i) qxx(i1,j)=s1-11*qxx(i1+1,j) qxx(i2,j)=s2-11*qxx(i2-1,j) 10 continue do 40 i=1,nx-1 IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny s3=(13*q(i,j3)-27*q(i,j3+1)+15*q(i,j3+2)-q(i,j3+3))/dy/dy s4=(-q(i,j4-3)+15*q(i,j4-2)-27*q(i,j4-1)+13*q(i,j4))/dy/dy do 51 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 51 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j3+1)=b(j3+1)-11*alp d(j3+1)=d(j3+1)-s3*alp b(j4-1)=b(j4-1)-11*alp d(j4-1)=d(j4-1)-s4*alp l1=j3+2 l2=j4-1 call tridi do 61 j=j3+1,j4-1 61 qyy(i,j)=d(j) qyy(i,j3)=s1-11*qyy(i,j3+1) qyy(i,j4)=s2-11*qyy(i,j4-1) s5=(13*q(i,j5)-27*q(i,j5+1)+15*q(i,j5+2)-q(i,j5+3))/dy/dy s6=(-q(i,j6-3)+15*q(i,j6-2)-27*q(i,j6-1)+13*q(i,j6))/dy/dy do 71 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 71 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j5+1)=b(j5+1)-11*alp d(j5+1)=d(j5+1)-s5*alp b(j6-1)=b(j6-1)-11*alp d(j6-1)=d(j6-1)-s6*alp l1=j5+2

96

l2=j6-1 call tridi do 81 j=j5+1,j6-1 81 qyy(i,j)=d(j) qyy(i,j5)=s1-11*qyy(i,j5+1) qyy(i,j6)=s2-11*qyy(i,j6-1)

91

82

50

60 40

s7=(13*q(i,j7)-27*q(i,j7+1)+15*q(i,j7+2)-q(i,j7+3))/dy/dy s8=(-q(i,j8-3)+15*q(i,j8-2)-27*q(i,j8-1)+13*q(i,j8))/dy/dy do 91 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j7+1)=b(j7+1)-11*alp d(j7+1)=d(j7+1)-s7*alp b(j8-1)=b(j8-1)-11*alp d(j8-1)=d(j8-1)-s8*alp l1=j7+2 l2=j8-1 call tridi do 82 j=j7+1,j8-1 qyy(i,j)=d(j) qyy(i,j7)=s1-11*qyy(i,j7+1) qyy(i,j8)=s2-11*qyy(i,j8-1) else j1=1 j2=ny endif s1=(13*q(i,j1)-27*q(i,j1+1)+15*q(i,j1+2)-q(i,j1+3))/dy/dy s2=(-q(i,j2-3)+15*q(i,j2-2)-27*q(i,j2-1)+13*q(i,j2))/dy/dy do 50 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j1+1)=b(j1+1)-11*alp d(j1+1)=d(j1+1)-s1*alp b(j2-1)=b(j2-1)-11*alp d(j2-1)=d(j2-1)-s2*alp l1=j1+2 l2=j2-1 call tridi do 60 j=j1+1,j2-1 qyy(i,j)=d(j) qyy(i,j1)=s1-11*qyy(i,j1+1) qyy(i,j2)=s2-11*qyy(i,j2-1) continue

97

return end C------------------------C Algoritma Thomas C------------------------SUBROUTINE TRIDI parameter (m=500) COMMON /AA4/A(m),B(m),C(m),D(m),L1,l2 DO 1 I=L1,L2 RT=-A(I)/B(I-1) B(I)=B(I)+RT*C(I-1) 1 D(I)=D(I)+RT*D(I-1) D(L2)=D(L2)/B(L2) DO 2 I=L2-1,L1-1,-1 2 D(I)=(D(I)-C(I)*D(I+1))/B(I) RETURN END C-------------------------C Hasil perhitungan C-------------------------subroutine hasil PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 WRITE(2,*)NX WRITE(2,*)NY do 10 j=1,ny do 10 i=1,nx WRITE(3,*)x(i,j),y(i,j),u(i,j),v(i,j)!,p(i,j) 10 WRITE(4,*)x(i,j),y(i,j),o(i,j) return end

98

Lampiran 8. Program Konveksi Alami Pada 4 Sirip C------------------------------------------------------------------------C Program penyelesaian konveksi alami pada 4 sirip C------------------------------------------------------------------------PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 OPEN(8,FILE='C:\matlab6p1\work\temp1') OPEN(7,FILE='C:\matlab6p1\work\pv1') OPEN(4,FILE='C:\matlab6p1\work\temp') OPEN(3,FILE='C:\matlab6p1\work\pv') OPEN(2,FILE='C:\matlab6p1\work\num') OPEN(1,FILE='C:\matlab6p1\work\div') call awal do k=1,nt ck=0. call rkv do i=1,nx do j=1,ny ck=ck+abs(ux(i,j)+vy(i,j)) end do end do if(ck.gt.0) then WRITE(*,*)k,LOG10(ck/nx/ny) WRITE(1,*)k,LOG10(ck/nx/ny) endif end do call hasil stop end C---------------------------------------C Syarat awal dan syarat batas C---------------------------------------subroutine awal PARAMETER(m=500,n=500) CHARACTER mul*2 COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 WRITE(*,*)' t= ' READ(*,*)tt

99

nx=201 ny=201 i1=1 i2=nx j1=1 j2=ny dx=1./200. dy=1./200. dt=0.0025 nt=tt/dt+1 pr=0.71 ra=1000000. WRITE(*,*)' Dari awal ?' READ(*,'(a)')mul IF(mul.eq.'y')then do i=1,nx do j=1,ny x(i,j)=(i-1)*dx y(i,j)=(j-1)*dy u(i,j)=0. v(i,j)=0. p(i,j)=0. o(i,j)=0. enddo enddo do j=1,ny o(1,j)=0.5 o(nx,j)=-0.5 enddo else do j=1,ny do i=1,nx READ(7,*)x(i,j),y(i,j),u(i,j),v(i,j),p(i,j) READ(8,*)x(i,j),y(i,j),o(i,j) end do end do end if return end C-------------------------------------------------------C Penyelesaian persamaan Navier - Stokes C-------------------------------------------------------subroutine rkv PARAMETER(m=500,n=500)

100

COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), vx(m,n),vy(m,n),px(m,n),py(m,n), uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 DIMENSION h1(m,n,6),h2(m,n,6),s1(m,n,6),s2(m,n,6),a(5),b(5) DATA a(1),a(2),a(3),a(4),a(5)/0.,-0.41789047,-1.19215169, 1 -1.69778469,-1.51418344/ DATA b(1),b(2),b(3),b(4),b(5)/0.14965902,0.37921031,0.82295502, 1 0.69945045,0.15305724/ 1 1 1

do i=2,nx-1 do j=2,ny-1 s1(i,j,1)=u(i,j) s2(i,j,1)=v(i,j) end do end do do k=1,5 call derv1(u,ux,uy) call derv1(v,vx,vy) call derv2(u,uxx,uyy) call derv2(v,vxx,vyy) call derp1(p,px,py) do i=2,nx-1 do j=2,ny-1 h1(i,j,k)=-u(i,j)*ux(i,j)-v(i,j)*uy(i,j)-px(i,j) 1 +(uxx(i,j)+uyy(i,j))*pr/SQRT(ra) h2(i,j,k)=-u(i,j)*vx(i,j)-v(i,j)*vy(i,j)-py(i,j) 1 +(vxx(i,j)+vyy(i,j))*pr/SQRT(ra)+pr*o(i,j) end do end do do i=2,nx-1 do j=2,ny-1 IF(i.le.101.and.j.ge.38.and.j.le.42)then u(i,j)=0. v(i,j)=0 ELSEIF (i.le.101.and.j.ge.78.and.j.le.82)then u(i,j)=0. v(i,j)=0 ELSEIF(i.le.101.and.j.ge.118.and.j.le.122)then u(i,j)=0. v(i,j)=0 elseif(i.le.101.and.j.ge.158.and.j.le.162)then u(i,j)=0. v(i,j)=0 else

101

h1(i,j,k)=h1(i,j,k)+a(k)*h1(i,j,k-1) h2(i,j,k)=h2(i,j,k)+a(k)*h2(i,j,k-1) s1(i,j,k+1)=s1(i,j,k)+b(k)*dt*h1(i,j,k) s2(i,j,k+1)=s2(i,j,k)+b(k)*dt*h2(i,j,k) u(i,j)=s1(i,j,k+1) v(i,j)=s2(i,j,k+1) endif end do end do call derv1(u,ux,uy) call derv1(v,vx,vy) c----------------------------------------------------------c Tekananan dihitung dengan artifisial kompresibiliti c----------------------------------------------------------do i=2,nx-1 do j=2,ny-1 IF(i.lt.101.and.j.gt.38.and.j.lt.42)then p(i,j)=0 ELSEIF (i.lt.101.and.j.gt.78.and.j.lt.82)then p(i,j)=0 ELSEIF(i.lt.101.and.j.gt.118.and.j.lt.122)then p(i,j)=0 elseif(i.lt.101.and.j.gt.158.and.j.lt.162)then p(i,j)=0 else p(i,j)=p(i,j)-(ux(i,j)+vy(i,j))*dt*0.5 endif end do end do do i=1,nx IF(i.le.101)then j1=1 j2=38 j3=42 j4=78 j5=82 j6=118 j7=122 j8=158 j9=162 j10=ny else j1=1 j2=ny endif p(i,j1)=(48*p(i,j1+1)-36*p(i,j1+2)+16*p(i,j1+3)-3*p(i,j1+4))/25 p(i,j2)=(48*p(i,j2-1)-36*p(i,j2-2)+16*p(i,j2-3)-3*p(i,j2-4))/25

102

p(i,j3)=(48*p(i,j3+1)-36*p(i,j3+2)+16*p(i,j3+3)-3*p(i,j3+4))/25 p(i,j4)=(48*p(i,j4-1)-36*p(i,j4-2)+16*p(i,j4-3)-3*p(i,j4-4))/25 p(i,j5)=(48*p(i,j5+1)-36*p(i,j5+2)+16*p(i,j5+3)-3*p(i,j5+4))/25 p(i,j6)=(48*p(i,j6-1)-36*p(i,j6-2)+16*p(i,j6-3)-3*p(i,j6-4))/25 p(i,j7)=(48*p(i,j7+1)-36*p(i,j7+2)+16*p(i,j7+3)-3*p(i,j7+4))/25 p(i,j8)=(48*p(i,j8-1)-36*p(i,j8-2)+16*p(i,j8-3)-3*p(i,j8-4))/25 p(i,j9)=(48*p(i,j9+1)-36*p(i,j9+2)+16*p(i,j9+3)-3*p(i,j9+4))/25 p(i,j10)=(48*p(i,j10-1)-36*p(i,j10-2)+16*p(i,j10-3)-3*p(i,j10-4)) 1 /25 end do i1=1 i2=nx j1=1 j2=ny do j=1,ny IF(j.ge.38.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101 elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif p(i1,j)=(48*p(i1+1,j)-36*p(i1+2,j)+16*p(i1+3,j)-3*p(i1+4,j))/25 p(i2,j)=(48*p(i2-1,j)-36*p(i2-2,j)+16*p(i2-3,j)-3*p(i2-4,j))/25 enddo call rkt c----------------------------------------------------------end do return end C-------------------------------------------C Penyelesaian persamaan energi C------------------------------------------subroutine rkt PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 DIMENSION h(m,n,6),s(m,n,6),a(5),b(5)

103

DATA a(1),a(2),a(3),a(4),a(5)/0.,-0.41789047,-1.19215169, 1 -1.69778469,-1.51418344/ DATA b(1),b(2),b(3),b(4),b(5)/0.14965902,0.37921031,0.82295502, 1 0.69945045,0.15305724/ do i=2,nx-1 do j=2,ny-1 s(i,j,1)=o(i,j) end do end do do k=1,5 call dero1(o,ox,oy) call derv2(o,oxx,oyy) do i=2,nx-1 do j=2,ny-1 IF(i.lt.101.and.j.gt.38.and.j.lt.42)then o(i,j)=0. elseif(i.lt.101.and.j.gt.78.and.j.lt.82)then o(i,j)=0. ELSEIF(i.lt.101.and.j.gt.118.and.j.lt.122)then o(i,j)=0. ELSEIF(i.lt.101.and.j.gt.158.and.j.lt.162)then o(i,j)=0. else h(i,j,k)=-u(i,j)*ox(i,j)-v(i,j)*oy(i,j)+ 1 (oxx(i,j)+oyy(i,j))/SQRT(ra)+a(k)*h(i,j,k-1) s(i,j,k+1)=s(i,j,k)+b(k)*dt*h(i,j,k) o(i,j)=s(i,j,k+1) endif end do end do do i=1,nx IF(i.le.101)then j1=1 j2=38 j3=42 j4=82 j5=78 j6=118 j7=122 j8=158 j9=162 j10=ny o(i,j1)=(48*o(i,j1+1)-36*o(i,j1+2)+16*o(i,j1+3)-3*o(i,j1+4))/25 o(i,j2)=0.5 o(i,j3)=0.5

104

o(i,j4)=0.5 o(i,j5)=0.5 o(i,j6)=0.5 o(i,j7)=0.5 o(i,j8)=0.5 o(i,j9)=0.5 o(i,j10)=(48*o(i,j10-1)-36*o(i,j10-2)+16*o(i,j10-3)-3*o(i,j10-4)) 1 /25 else j1=1 j2=ny o(i,j1)=(48*o(i,j1+1)-36*o(i,j1+2)+16*o(i,j1+3)-3*o(i,j1+4))/25 o(i,j2)=(48*o(i,j2-1)-36*o(i,j2-2)+16*o(i,j2-3)-3*o(i,j2-4))/25 endif end do do j=1,ny IF(j.ge.38.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101 elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif o(i1,j)=0.5 end do i1=1 i2=nx j1=1 j2=ny do i=1,nx do j=1,ny IF(i.lt.101.and.j.gt.38.and.j.lt.42)then o(i,j)=0. elseif(i.lt.101.and.j.gt.78.and.j.lt.82)then o(i,j)=0. ELSEIF(i.lt.101.and.j.gt.118.and.j.lt.122)then o(i,j)=0. ELSEIF(i.lt.101.and.j.gt.158.and.j.lt.162)then o(i,j)=0. endif end do

105

end do end do return end C-----------------------------------------C Turunan pertama variabel kecepatan C-----------------------------------------subroutine derv1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.38.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101 elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo

106

do i=1,nx IF(i.le.101)then j1=1 j2=38 j3=42 j4=78 j5=82 j6=118 j7=122 j8=158 j9=162 j10=ny qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1) 1 +25*q(i,j6))/12/dy do 222 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 222 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 322 j=j5+1,j6-1

107

322 qy(i,j)=d(j) qy(i,j7)=(-3*q(i,j7+4)+16*q(i,j7+3)-36*q(i,j7+2)+48*q(i,j7+1) 1 -25*q(i,j7))/12/dy qy(i,j8)=(3*q(i,j8-4)-16*q(i,j8-3)+36*q(i,j8-2)-48*q(i,j8-1) 1 +25*q(i,j8))/12/dy do 333 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp 333 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8) c(j8-1)=0. l1=j7+2 l2=j8-1 call tridi do 433 j=j7+1,j8-1 433 qy(i,j)=d(j) qy(i,j9)=(-3*q(i,j9+4)+16*q(i,j9+3)-36*q(i,j9+2)+48*q(i,j9+1) 1 -25*q(i,j9))/12/dy qy(i,j10)=(3*q(i,j10-4)-16*q(i,j10-3)+36*q(i,j10-2)-48*q(i,j10-1) 1 +25*q(i,j10))/12/dy do 666 j=j9+1,j10-1 a(j)=alp b(j)=1. c(j)=alp 666 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j9+1)=d(j9+1)-a(j9+1)*qy(i,j9) a(j9+1)=0. d(j10-1)=d(j10-1)-c(j10-1)*qy(i,j10) c(j10-1)=0. l1=j9+2 l2=j10-1 call tridi do 766 j=j9+1,j10-1 766 qy(i,j)=d(j) else j1=1 j2=ny endif qy(i,j1)=(-3*q(i,j1+4)+16*q(i,j1+3)-36*q(i,j1+2)+48*q(i,j1+1) 1 -25*q(i,j1))/12/dy qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1)

108

1 +25*q(i,j2))/12/dy do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end C-----------------------------------------C turunan pertama variabel temperatur C-----------------------------------------subroutine dero1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.38.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101 elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx

109

do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=38 j3=42 j4=78 j5=82 j6=118 j7=122 j8=158 j9=162 j10=ny qy(i,j1)=0. qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1) 1 +25*q(i,j6))/12/dy qy(i,j7)=(-3*q(i,j7+4)+16*q(i,j7+3)-36*q(i,j7+2)+48*q(i,j7+1) 1 -25*q(i,j7))/12/dy qy(i,j8)=(3*q(i,j8-4)-16*q(i,j8-3)+36*q(i,j8-2)-48*q(i,j8-1) 1 +25*q(i,j8))/12/dy qy(i,j9)=(-3*q(i,j9+4)+16*q(i,j9+3)-36*q(i,j9+2)+48*q(i,j9+1) 1 -25*q(i,j9))/12/dy qy(i,j10)=0. do 110 j=j1+1,j2-1 a(j)=alp

110

b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 210j=j1+1,j2-1 210 qy(i,j)=d(j) do 111 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 111 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 211 j=j3+1,j4-1 211 qy(i,j)=d(j) do 345 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 345 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 789 j=j5+1,j6-1 789 qy(i,j)=d(j) do 444 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp

111

444

d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8) c(j8-1)=0. l1=j7+2 l2=j8-1 call tridi do 544 j=j7+1,j8-1 544 qy(i,j)=d(j) do 255 j=j9+1,j10-1 a(j)=alp b(j)=1. c(j)=alp 255 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j9+1)=d(j9+1)-a(j9+1)*qy(i,j9) a(j9+1)=0. d(j10-1)=d(j10-1)-c(j10-1)*qy(i,j10) c(j10-1)=0. l1=j9+2 l2=j10-1 call tridi do 665 j=j9+1,j10-1 665 qy(i,j)=d(j) else j1=1 j2=ny qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) endif enddo

112

return end C---------------------------------------C Turunan pertama variabel tekanan C---------------------------------------subroutine derp1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.48.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101 elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif qx(i1,j)=0. qx(i2,j)=0. do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then

113

j1=1 j2=38 j3=42 j4=78 j5=82 j6=118 j7=122 j8=158 j9=162 j10=ny qy(i,j3)=0. qy(i,j4)=0. do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) qy(i,j5)=0. qy(i,j6)=0. do 222 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 222 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 322 j=j5+1,j6-1 322 qy(i,j)=d(j) qy(i,j7)=0. qy(i,j8)=0. do 555 j=j7+1,j8-1

114

a(j)=alp b(j)=1. c(j)=alp 555 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8) c(j8-1)=0. l1=j7+2 l2=j8-1 call tridi do 655 j=j7+1,j8-1 655 qy(i,j)=d(j) qy(i,j9)=0. qy(i,j10)=0. do 777 j=j9+1,j10-1 a(j)=alp b(j)=1. c(j)=alp 777 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j9+1)=d(j9+1)-a(j9+1)*qy(i,j9) a(j9+1)=0. d(j10-1)=d(j10-1)-c(j10-1)*qy(i,j10) c(j10-1)=0. l1=j9+2 l2=j10-1 call tridi do 877 j=j9+1,j10-1 877 qy(i,j)=d(j) else j1=1 j2=ny endif qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2

115

l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end C-------------------------------------------------------C Turunan kedua variabel kecepatan dan temperatur C-------------------------------------------------------subroutine derv2(q,qxx,qyy) COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qxx(500,500),qyy(500,500) alp=0.1 ak=1.2 do 10 j=2,ny-1 IF(j.ge.38.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101 elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif s1=(13*q(i1,j)-27*q(i1+1,j)+15*q(i1+2,j)-1*q(i1+3,j))/dx/dx s2=(-1*q(i2-3,j)+15*q(i2-2,j)-27*q(i2-1,j)+13*q(i2,j))/dx/dx do 20 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 20 d(i)=ak*(q(i+1,j)-2*q(i,j)+q(i-1,j))/dx/dx b(i1+1)=b(i1+1)-11*alp d(i1+1)=d(i1+1)-s1*alp b(i2-1)=b(i2-1)-11*alp d(i2-1)=d(i2-1)-s2*alp l1=i1+2 l2=i2-1 call tridi do 30 i=i1+1,i2-1 30 qxx(i,j)=d(i)

116

qxx(i1,j)=s1-11*qxx(i1+1,j) qxx(i2,j)=s2-11*qxx(i2-1,j) 10 continue do 40 i=1,nx-1 IF(i.le.101)then j1=1 j2=38 j3=42 j4=78 j5=82 j6=118 j7=122 j8=158 j9=162 j10=ny s3=(13*q(i,j3)-27*q(i,j3+1)+15*q(i,j3+2)-q(i,j3+3))/dy/dy s4=(-q(i,j4-3)+15*q(i,j4-2)-27*q(i,j4-1)+13*q(i,j4))/dy/dy do 51 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 51 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j3+1)=b(j3+1)-11*alp d(j3+1)=d(j3+1)-s3*alp b(j4-1)=b(j4-1)-11*alp d(j4-1)=d(j4-1)-s4*alp l1=j3+2 l2=j4-1 call tridi do 61 j=j3+1,j4-1 61 qyy(i,j)=d(j) qyy(i,j3)=s1-11*qyy(i,j3+1) qyy(i,j4)=s2-11*qyy(i,j4-1) s5=(13*q(i,j5)-27*q(i,j5+1)+15*q(i,j5+2)-q(i,j5+3))/dy/dy s6=(-q(i,j6-3)+15*q(i,j6-2)-27*q(i,j6-1)+13*q(i,j6))/dy/dy do 71 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 71 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j5+1)=b(j5+1)-11*alp d(j5+1)=d(j5+1)-s5*alp b(j6-1)=b(j6-1)-11*alp d(j6-1)=d(j6-1)-s6*alp l1=j5+2 l2=j6-1

117

call tridi do 81 j=j5+1,j6-1 81 qyy(i,j)=d(j) qyy(i,j5)=s1-11*qyy(i,j5+1) qyy(i,j6)=s2-11*qyy(i,j6-1) s7=(13*q(i,j7)-27*q(i,j7+1)+15*q(i,j7+2)-q(i,j7+3))/dy/dy s8=(-q(i,j8-3)+15*q(i,j8-2)-27*q(i,j8-1)+13*q(i,j8))/dy/dy do 91 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp 91 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j7+1)=b(j7+1)-11*alp d(j7+1)=d(j7+1)-s7*alp b(j8-1)=b(j8-1)-11*alp d(j8-1)=d(j8-1)-s8*alp l1=j7+2 l2=j8-1 call tridi do 82 j=j7+1,j8-1 82 qyy(i,j)=d(j) qyy(i,j7)=s1-11*qyy(i,j7+1) qyy(i,j8)=s2-11*qyy(i,j8-1) s9=(13*q(i,j9)-27*q(i,j9+1)+15*q(i,j9+2)-q(i,j9+3))/dy/dy s10=(-q(i,j10-3)+15*q(i,j10-2)-27*q(i,j10-1)+13*q(i,j10))/dy/dy do 92 j=j9+1,j10-1 a(j)=alp b(j)=1. c(j)=alp 92 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j9+1)=b(j9+1)-11*alp d(j9+1)=d(j9+1)-s9*alp b(j10-1)=b(j10-1)-11*alp d(j10-1)=d(j10-1)-s10*alp l1=j9+2 l2=j10-1 call tridi do 83 j=j9+1,j10-1 83 qyy(i,j)=d(j) qyy(i,j9)=s1-11*qyy(i,j9+1) qyy(i,j10)=s2-11*qyy(i,j10-1) else j1=1 j2=ny endif

118

s1=(13*q(i,j1)-27*q(i,j1+1)+15*q(i,j1+2)-q(i,j1+3))/dy/dy s2=(-q(i,j2-3)+15*q(i,j2-2)-27*q(i,j2-1)+13*q(i,j2))/dy/dy do 50 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 50 d(j)=ak*(q(i,j+1)-2*q(i,j)+q(i,j-1))/dy/dy b(j1+1)=b(j1+1)-11*alp d(j1+1)=d(j1+1)-s1*alp b(j2-1)=b(j2-1)-11*alp d(j2-1)=d(j2-1)-s2*alp l1=j1+2 l2=j2-1 call tridi do 60 j=j1+1,j2-1 60 qyy(i,j)=d(j) qyy(i,j1)=s1-11*qyy(i,j1+1) qyy(i,j2)=s2-11*qyy(i,j2-1) 40 continue return end C------------------------C Algoritma Thomas C------------------------SUBROUTINE TRIDI parameter (m=500) COMMON /AA4/A(m),B(m),C(m),D(m),L1,l2 DO 1 I=L1,L2 RT=-A(I)/B(I-1) B(I)=B(I)+RT*C(I-1) 1 D(I)=D(I)+RT*D(I-1) D(L2)=D(L2)/B(L2) DO 2 I=L2-1,L1-1,-1 2 D(I)=(D(I)-C(I)*D(I+1))/B(I) RETURN END C-------------------------C Hasil perhitungan C-------------------------subroutine hasil PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),p(m,n),ux(m,n),uy(m,n), 1 vx(m,n),vy(m,n),px(m,n),py(m,n), 1 uxx(m,n),uyy(m,n),vxx(m,n),vyy(m,n),x(m,n),y(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n)

119

COMMON/aa2/nx,ny,nt,dx,dy,dt,pr,ra,i1,i2,j1,j2 WRITE(2,*)NX WRITE(2,*)NY do 10 j=1,ny do 10 i=1,nx WRITE(3,*)x(i,j),y(i,j),u(i,j),v(i,j)!,p(i,j) 10 WRITE(4,*)x(i,j),y(i,j),o(i,j) return end

120

Lampiran 9. Program Tambahan Untuk Sirip Tunggal c -------------------------------------------------------------------c PROGRAM TAMBAHAN UNTUK SIRIP TUNGGAL c -------------------------------------------------------------------PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),ux(m,n),uy(m,n),vx(m,n),vy(m,n) 1 ,p(m,n),x(m,n),y(m,n),sf(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt OPEN(1,FILE='C:\matlab6p1\work\pv') OPEN(2,FILE='C:\matlab6p1\work\temp') OPEN(3,FILE='C:\matlab6p1\work\streamf') OPEN(4,FILE='C:\matlab6p1\work\num') nx=201 ny=201 dx=1./200. dy=1./200. do j=1,ny do i=1,nx READ(1,*)x(i,j),y(i,j),u(i,j),v(i,j)!,p(i,j) READ(2,*)x(i,j),y(i,j),o(i,j) end do end do call derv1(u,ux,uy) call derv1(v,vx,vy) call stream WRITE(4,*)nx WRITE(4,*)ny do 10 j=1,ny do 10 i=1,nx 10 WRITE(3,*)x(i,j),y(i,j),-1*sf(i,j) stop end subroutine dero1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.98.and.j.le.102)then i1=101 else i1=1 i2=nx endif

121

qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny qy(i,j1)=0. qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=0. do 110 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 210j=j1+1,j2-1 210 qy(i,j)=d(j)

122

do 111 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 111 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 211 j=j3+1,j4-1 211 qy(i,j)=d(j) else j1=1 j2=ny qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) endif enddo return end cxx subroutine derv1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny

123

IF(j.ge.98.and.j.le.102)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1

124

call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) else j1=1 j2=ny endif qy(i,j1)=(-3*q(i,j1+4)+16*q(i,j1+3)-36*q(i,j1+2)+48*q(i,j1+1) 1 -25*q(i,j1))/12/dy qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end

1 2

SUBROUTINE TRIDI parameter (m=500) COMMON /AA4/A(m),B(m),C(m),D(m),L1,l2 DO 1 I=L1,L2 RT=-A(I)/B(I-1) B(I)=B(I)+RT*C(I-1) D(I)=D(I)+RT*D(I-1) D(L2)=D(L2)/B(L2) DO 2 I=L2-1,L1-1,-1 D(I)=(D(I)-C(I)*D(I+1))/B(I) RETURN END subroutine stream parameter (m=500,n=500) COMMON/aa1/u(m,n),v(m,n),ux(m,n),uy(m,n),vx(m,n),vy(m,n) 1 ,p(m,n),x(m,n),y(m,n),sf(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n)

125

COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/AA4/A(m),B(m),C(m),D(m),L1,l2 DIMENSION f(m,n) ermx=0.000001 k=0 c----------------------------------------c initial condition c----------------------------------------do i=2,nx-1 IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny do j=j3+1,j4-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do else j1=1 j2=ny endif do j=j1+1,j2-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do end do c----------------------------------------c boundary condition c----------------------------------------do i=1,nx IF(i.le.101)then j1=1 j2=98 j3=102 j4=ny sf(i,j3)=0. sf(i,j4)=0. else j1=1 j2=ny endif sf(i,j1)=0. sf(i,j2)=0. end do do j=1,ny

126

IF(j.ge.98.and.j.le.102)then i1=101 else i1=1 i2=nx endif sf(i1,j)=0. sf(i2,j)=0. end do c----------------------------------------15 ERR=0. k=k+1 IF(k.gt.20000)stop c----------------------------------------c matric coefficient c----------------------------------------do j=2,ny-1 IF(j.ge.98.and.j.le.102)then i1=101 else i1=1 i2=nx endif do i=i1+1,i2-1 a(i)=1. b(i)=-4. c(i)=1. d(i)=-sf(i,j+1)-sf(i,j-1)+f(i,j)*dx*dx enddo c----------------------------------------c modification matric coefficien c----------------------------------------d(i1+1)=d(i1+1)-a(i1+1)*sf(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*sf(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do i=i1+1,i2-1 ERR=ERR+ABS(d(i)-sf(i,j)) sf(i,j)=d(i) end do enddo WRITE(*,*)k,err IF(err.gt.ermx)GOTO 15 return end

127

Lampiran 10. Program Tambahan Untuk 2 Sirip c c c

-------------------------------------------------------PROGRAM TAMBAHAN UNTUK 2 SIRIP -------------------------------------------------------PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),ux(m,n),uy(m,n),vx(m,n),vy(m,n) 1 ,p(m,n),x(m,n),y(m,n),sf(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt OPEN(1,FILE='C:\matlab6p1\work\pv') OPEN(2,FILE='C:\matlab6p1\work\temp') OPEN(3,FILE='C:\matlab6p1\work\streamf') OPEN(4,FILE='C:\matlab6p1\work\num') nx=201 ny=201 dx=1./200. dy=1./200. do j=1,ny do i=1,nx READ(1,*)x(i,j),y(i,j),u(i,j),v(i,j)!,p(i,j) READ(2,*)x(i,j),y(i,j),o(i,j) end do end do call derv1(u,ux,uy) call derv1(v,vx,vy) call stream WRITE(4,*)nx WRITE(4,*)ny do 10 j=1,ny do 10 i=1,nx 10 WRITE(3,*)x(i,j),y(i,j),-1*sf(i,j) stop end subroutine dero1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else

128

i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135 j6=ny qy(i,j1)=0. qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=0. do 110 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp

129

110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 210j=j1+1,j2-1 210 qy(i,j)=d(j) do 111 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 111 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 211 j=j3+1,j4-1 211 qy(i,j)=d(j) do 333 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 333 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 433 j=j5+1,j6-1 433 qy(i,j)=d(j) else j1=1 j2=ny qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1.

130

c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) endif enddo return end cxx

subroutine derv1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5

do j=1,ny IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1

131

call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135 j6=ny qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1) 1 +25*q(i,j6))/12/dy do 330 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 330 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2

132

l2=j6-1 call tridi do 340 j=j5+1,j6-1 340 qy(i,j)=d(j) else j1=1 j2=ny endif qy(i,j1)=(-3*q(i,j1+4)+16*q(i,j1+3)-36*q(i,j1+2)+48*q(i,j1+1) 1 -25*q(i,j1))/12/dy qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end

1 2

SUBROUTINE TRIDI parameter (m=500) COMMON /AA4/A(m),B(m),C(m),D(m),L1,l2 DO 1 I=L1,L2 RT=-A(I)/B(I-1) B(I)=B(I)+RT*C(I-1) D(I)=D(I)+RT*D(I-1) D(L2)=D(L2)/B(L2) DO 2 I=L2-1,L1-1,-1 D(I)=(D(I)-C(I)*D(I+1))/B(I) RETURN END subroutine stream parameter (m=500,n=500)

133

COMMON/aa1/u(m,n),v(m,n),ux(m,n),uy(m,n),vx(m,n),vy(m,n) ,p(m,n),x(m,n),y(m,n),sf(m,n) ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/AA4/A(m),B(m),C(m),D(m),L1,l2 DIMENSION f(m,n) ermx=0.000001 k=0 c----------------------------------------c initial condition c----------------------------------------do i=2,nx-1 IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135 j6=ny 1 1

do j=j3+1,j4-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do do j=j5+1,j6-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do else j1=1 j2=ny endif do j=j1+1,j2-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do end do c----------------------------------------c boundary condition c----------------------------------------do i=1,nx IF(i.le.101)then j1=1 j2=65 j3=69 j4=131 j5=135

134

j6=ny sf(i,j3)=0. sf(i,j4)=0. sf(i,j5)=0. sf(i,j6)=0. else j1=1 j2=ny endif sf(i,j1)=0. sf(i,j2)=0. end do do j=1,ny IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else i1=1 i2=nx endif sf(i1,j)=0. sf(i2,j)=0. end do c----------------------------------------15 ERR=0. k=k+1 IF(k.gt.20000)stop c----------------------------------------c matric coefficient c----------------------------------------do j=2,ny-1 IF(j.ge.65.and.j.le.69)then i1=101 elseif(j.ge.131.and.j.le.135)then i1=101 else i1=1 i2=nx endif do i=i1+1,i2-1 a(i)=1. b(i)=-4. c(i)=1. d(i)=-sf(i,j+1)-sf(i,j-1)+f(i,j)*dx*dx enddo c-----------------------------------------

135

c modification matric coefficien c----------------------------------------d(i1+1)=d(i1+1)-a(i1+1)*sf(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*sf(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do i=i1+1,i2-1 ERR=ERR+ABS(d(i)-sf(i,j)) sf(i,j)=d(i) end do enddo WRITE(*,*)k,err IF(err.gt.ermx)GOTO 15 return end

136

Lampiran11. Program Tambahan Untuk 3 Sirip c -------------------------------------------------------c PROGRAM TAMBAHAN UNTUK 3 SIRIP c -------------------------------------------------------PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),ux(m,n),uy(m,n),vx(m,n),vy(m,n) 1 ,p(m,n),x(m,n),y(m,n),sf(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt OPEN(1,FILE='C:\matlab6p1\work\pv') OPEN(2,FILE='C:\matlab6p1\work\temp') OPEN(3,FILE='C:\matlab6p1\work\streamf') OPEN(4,FILE='C:\matlab6p1\work\num') nx=201 ny=201 dx=1./200. dy=1./200. do j=1,ny do i=1,nx READ(1,*)x(i,j),y(i,j),u(i,j),v(i,j)!,p(i,j) READ(2,*)x(i,j),y(i,j),o(i,j) end do end do call derv1(u,ux,uy) call derv1(v,vx,vy) call stream WRITE(4,*)nx WRITE(4,*)ny do 10 j=1,ny do 10 i=1,nx 10 WRITE(3,*)x(i,j),y(i,j),-1*sf(i,j) stop end subroutine dero1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.48.and.j.le.52)then i1=101 elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then

137

i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny qy(i,j1)=0. qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1)

138

1 +25*q(i,j6))/12/dy qy(i,j7)=(-3*q(i,j7+4)+16*q(i,j7+3)-36*q(i,j7+2)+48*q(i,j7+1) 1 -25*q(i,j7))/12/dy qy(i,j8)=0. do 110 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 210j=j1+1,j2-1 210 qy(i,j)=d(j) do 111 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 111 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 211 j=j3+1,j4-1 211 qy(i,j)=d(j) do 333 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 333 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 433 j=j5+1,j6-1

139

433 qy(i,j)=d(j) do 555 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp 555 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8) c(j8-1)=0. l1=j7+2 l2=j8-1 call tridi do 655 j=j7+1,j8-1 655 qy(i,j)=d(j) else j1=1 j2=ny qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) endif enddo return end cxx subroutine derv1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5

140

do j=1,ny IF(j.ge.48.and.j.le.52)then i1=101 elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy do 110 j=j3+1,j4-1

141

a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1) 1 +25*q(i,j6))/12/dy do 330 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 330 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 340 j=j5+1,j6-1 340 qy(i,j)=d(j) qy(i,j7)=(-3*q(i,j7+4)+16*q(i,j7+3)-36*q(i,j7+2)+48*q(i,j7+1) 1 -25*q(i,j7))/12/dy qy(i,j8)=(3*q(i,j8-4)-16*q(i,j8-3)+36*q(i,j8-2)-48*q(i,j8-1) 1 +25*q(i,j8))/12/dy do 440 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp 440 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8) c(j8-1)=0. l1=j7+2

142

l2=j8-1 call tridi do 450 j=j7+1,j8-1 450 qy(i,j)=d(j) else j1=1 j2=ny endif qy(i,j1)=(-3*q(i,j1+4)+16*q(i,j1+3)-36*q(i,j1+2)+48*q(i,j1+1) 1 -25*q(i,j1))/12/dy qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end

1 2

SUBROUTINE TRIDI parameter (m=500) COMMON /AA4/A(m),B(m),C(m),D(m),L1,l2 DO 1 I=L1,L2 RT=-A(I)/B(I-1) B(I)=B(I)+RT*C(I-1) D(I)=D(I)+RT*D(I-1) D(L2)=D(L2)/B(L2) DO 2 I=L2-1,L1-1,-1 D(I)=(D(I)-C(I)*D(I+1))/B(I) RETURN END subroutine stream

143

parameter (m=500,n=500) COMMON/aa1/u(m,n),v(m,n),ux(m,n),uy(m,n),vx(m,n),vy(m,n) 1 ,p(m,n),x(m,n),y(m,n),sf(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/AA4/A(m),B(m),C(m),D(m),L1,l2 DIMENSION f(m,n) ermx=0.000001 k=0 c----------------------------------------c initial condition c----------------------------------------do i=2,nx-1 IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny do j=j3+1,j4-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do do j=j5+1,j6-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do do j=j7+1,j8-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do else j1=1 j2=ny endif do j=j1+1,j2-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do end do c----------------------------------------c boundary condition c-----------------------------------------

144

do i=1,nx IF(i.le.101)then j1=1 j2=48 j3=52 j4=98 j5=102 j6=148 j7=152 j8=ny sf(i,j3)=0. sf(i,j4)=0. sf(i,j5)=0. sf(i,j6)=0. sf(i,j7)=0. sf(i,j8)=0. else j1=1 j2=ny endif sf(i,j1)=0. sf(i,j2)=0. end do do j=1,ny IF(j.ge.48.and.j.le.52)then i1=101 elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then i1=101 else i1=1 i2=nx endif sf(i1,j)=0. sf(i2,j)=0. end do c----------------------------------------15 ERR=0. k=k+1 IF(k.gt.20000)stop c----------------------------------------c matric coefficient c----------------------------------------do j=2,ny-1 IF(j.ge.48.and.j.le.52)then i1=101

145

elseif(j.ge.98.and.j.le.102)then i1=101 elseif(j.ge.148.and.j.le.152)then i1=101 else i1=1 i2=nx endif do i=i1+1,i2-1 a(i)=1. b(i)=-4. c(i)=1. d(i)=-sf(i,j+1)-sf(i,j-1)+f(i,j)*dx*dx enddo c----------------------------------------c modification matric coefficien c----------------------------------------d(i1+1)=d(i1+1)-a(i1+1)*sf(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*sf(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do i=i1+1,i2-1 ERR=ERR+ABS(d(i)-sf(i,j)) sf(i,j)=d(i) end do enddo WRITE(*,*)k,err IF(err.gt.ermx)GOTO 15 return end

146

Lampiran 12. Program Tambahan Untuk 4 Sirip c -------------------------------------------------------c PROGRAM TAMBAHAN UNTUK 4 SIRIP c -------------------------------------------------------PARAMETER(m=500,n=500) COMMON/aa1/u(m,n),v(m,n),ux(m,n),uy(m,n),vx(m,n),vy(m,n) 1 ,p(m,n),x(m,n),y(m,n),sf(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n) COMMON/aa2/nx,ny,nt,dx,dy,dt OPEN(1,FILE='C:\matlab6p1\work\pv') OPEN(2,FILE='C:\matlab6p1\work\temp') OPEN(3,FILE='C:\matlab6p1\work\streamf') OPEN(4,FILE='C:\matlab6p1\work\num') nx=201 ny=201 dx=1./200. dy=1./200. do j=1,ny do i=1,nx READ(1,*)x(i,j),y(i,j),u(i,j),v(i,j)!,p(i,j) READ(2,*)x(i,j),y(i,j),o(i,j) end do end do call derv1(u,ux,uy) call derv1(v,vx,vy) call stream WRITE(4,*)nx WRITE(4,*)ny do 10 j=1,ny do 10 i=1,nx 10 WRITE(3,*)x(i,j),y(i,j),-1*sf(i,j) stop end subroutine dero1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.38.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101

147

elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=38 j3=42 j4=78 j5=82 j6=118 j7=122 j8=158 j9=162 j10=ny qy(i,j1)=0. qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1)

148

1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1) 1 +25*q(i,j6))/12/dy qy(i,j7)=(-3*q(i,j7+4)+16*q(i,j7+3)-36*q(i,j7+2)+48*q(i,j7+1) 1 -25*q(i,j7))/12/dy qy(i,j8)=(3*q(i,j8-4)-16*q(i,j8-3)+36*q(i,j8-2)-48*q(i,j8-1) 1 +25*q(i,j8))/12/dy qy(i,j9)=(-3*q(i,j9+4)+16*q(i,j9+3)-36*q(i,j9+2)+48*q(i,j9+1) 1 -25*q(i,j9))/12/dy qy(i,j10)=0. do 110 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 210j=j1+1,j2-1 210 qy(i,j)=d(j) do 111 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 111 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 211 j=j3+1,j4-1 211 qy(i,j)=d(j) do 333 j=j5+1,j6-1 a(j)=alp b(j)=1. c(j)=alp 333 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0.

149

d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 433 j=j5+1,j6-1 433 qy(i,j)=d(j) do 555 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp 555 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8) c(j8-1)=0. l1=j7+2 l2=j8-1 call tridi do 655 j=j7+1,j8-1 655 qy(i,j)=d(j) do 666 j=j9+1,j10-1 a(j)=alp b(j)=1. c(j)=alp 666 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j9+1)=d(j9+1)-a(j9+1)*qy(i,j9) a(j9+1)=0. d(j10-1)=d(j10-1)-c(j10-1)*qy(i,j10) c(j10-1)=0. l1=j9+2 l2=j10-1 call tridi do 755 j=j9+1,j10-1 755 qy(i,j)=d(j) else j1=1 j2=ny qy(i,j1)=0. qy(i,j2)=0. do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1)

150

a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) endif enddo return end cxx subroutine derv1(q,qx,qy) COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/aa4/a(500),b(500),c(500),d(500),l1,l2 DIMENSION q(500,500),qx(500,500),qy(500,500) alp=0.25 ak=1.5 do j=1,ny IF(j.ge.38.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101 elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif qx(i1,j)=(-3*q(i1+4,j)+16*q(i1+3,j)-36*q(i1+2,j)+48*q(i1+1,j) 1 -25*q(i1,j))/12/dx qx(i2,j)=(3*q(i2-4,j)-16*q(i2-3,j)+36*q(i2-2,j)-48*q(i2-1,j) 1 +25*q(i2,j))/12/dx do 10 i=i1+1,i2-1 a(i)=alp b(i)=1. c(i)=alp 10 d(i)=ak*(q(i+1,j)-q(i-1,j))/2/dx d(i1+1)=d(i1+1)-a(i1+1)*qx(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*qx(i2,j)

151

c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do 20 i=i1+1,i2-1 20 qx(i,j)=d(i) enddo do i=1,nx IF(i.le.101)then j1=1 j2=38 j3=42 j4=78 j5=82 j6=118 j7=122 j8=158 j9=162 j10=ny qy(i,j3)=(-3*q(i,j3+4)+16*q(i,j3+3)-36*q(i,j3+2)+48*q(i,j3+1) 1 -25*q(i,j3))/12/dy qy(i,j4)=(3*q(i,j4-4)-16*q(i,j4-3)+36*q(i,j4-2)-48*q(i,j4-1) 1 +25*q(i,j4))/12/dy do 110 j=j3+1,j4-1 a(j)=alp b(j)=1. c(j)=alp 110 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j3+1)=d(j3+1)-a(j3+1)*qy(i,j3) a(j3+1)=0. d(j4-1)=d(j4-1)-c(j4-1)*qy(i,j4) c(j4-1)=0. l1=j3+2 l2=j4-1 call tridi do 210 j=j3+1,j4-1 210 qy(i,j)=d(j) qy(i,j5)=(-3*q(i,j5+4)+16*q(i,j5+3)-36*q(i,j5+2)+48*q(i,j5+1) 1 -25*q(i,j5))/12/dy qy(i,j6)=(3*q(i,j6-4)-16*q(i,j6-3)+36*q(i,j6-2)-48*q(i,j6-1) 1 +25*q(i,j6))/12/dy do 330 j=j5+1,j6-1 a(j)=alp b(j)=1.

152

c(j)=alp 330 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j5+1)=d(j5+1)-a(j5+1)*qy(i,j5) a(j5+1)=0. d(j6-1)=d(j6-1)-c(j6-1)*qy(i,j6) c(j6-1)=0. l1=j5+2 l2=j6-1 call tridi do 340 j=j5+1,j6-1 340 qy(i,j)=d(j) qy(i,j7)=(-3*q(i,j7+4)+16*q(i,j7+3)-36*q(i,j7+2)+48*q(i,j7+1) 1 -25*q(i,j7))/12/dy qy(i,j8)=(3*q(i,j8-4)-16*q(i,j8-3)+36*q(i,j8-2)-48*q(i,j8-1) 1 +25*q(i,j8))/12/dy do 440 j=j7+1,j8-1 a(j)=alp b(j)=1. c(j)=alp 440 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j7+1)=d(j7+1)-a(j7+1)*qy(i,j7) a(j7+1)=0. d(j8-1)=d(j8-1)-c(j8-1)*qy(i,j8) c(j8-1)=0. l1=j7+2 l2=j8-1 call tridi do 450 j=j7+1,j8-1 450 qy(i,j)=d(j) qy(i,j9)=(-3*q(i,j9+4)+16*q(i,j9+3)-36*q(i,j9+2)+48*q(i,j9+1) 1 -25*q(i,j9))/12/dy qy(i,j10)=(3*q(i,j10-4)-16*q(i,j10-3)+36*q(i,j10-2)-48*q(i,j10-1) 1 +25*q(i,j10))/12/dy do 550 j=j9+1,j10-1 a(j)=alp b(j)=1. c(j)=alp 550 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j9+1)=d(j9+1)-a(j9+1)*qy(i,j9) a(j9+1)=0. d(j10-1)=d(j10-1)-c(j10-1)*qy(i,j10) c(j10-1)=0. l1=j9+2 l2=j10-1 call tridi do 560 j=j9+1,j10-1

153

560 qy(i,j)=d(j) else j1=1 j2=ny endif qy(i,j1)=(-3*q(i,j1+4)+16*q(i,j1+3)-36*q(i,j1+2)+48*q(i,j1+1) 1 -25*q(i,j1))/12/dy qy(i,j2)=(3*q(i,j2-4)-16*q(i,j2-3)+36*q(i,j2-2)-48*q(i,j2-1) 1 +25*q(i,j2))/12/dy do 11 j=j1+1,j2-1 a(j)=alp b(j)=1. c(j)=alp 11 d(j)=ak*(q(i,j+1)-q(i,j-1))/2/dy d(j1+1)=d(j1+1)-a(j1+1)*qy(i,j1) a(j1+1)=0. d(j2-1)=d(j2-1)-c(j2-1)*qy(i,j2) c(j2-1)=0. l1=j1+2 l2=j2-1 call tridi do 21 j=j1+1,j2-1 21 qy(i,j)=d(j) enddo return end

1 2

SUBROUTINE TRIDI parameter (m=500) COMMON /AA4/A(m),B(m),C(m),D(m),L1,l2 DO 1 I=L1,L2 RT=-A(I)/B(I-1) B(I)=B(I)+RT*C(I-1) D(I)=D(I)+RT*D(I-1) D(L2)=D(L2)/B(L2) DO 2 I=L2-1,L1-1,-1 D(I)=(D(I)-C(I)*D(I+1))/B(I) RETURN END subroutine stream parameter (m=500,n=500) COMMON/aa1/u(m,n),v(m,n),ux(m,n),uy(m,n),vx(m,n),vy(m,n) 1 ,p(m,n),x(m,n),y(m,n),sf(m,n) 1 ,o(m,n),ox(m,n),oy(m,n),oxx(m,n),oyy(m,n)

154

COMMON/aa2/nx,ny,nt,dx,dy,dt COMMON/AA4/A(m),B(m),C(m),D(m),L1,l2 DIMENSION f(m,n) ermx=0.000001 k=0 c----------------------------------------c initial condition c----------------------------------------do i=2,nx-1 IF(i.le.101)then j1=1 j2=38 j3=42 j4=78 j5=82 j6=118 j7=122 j8=158 j9=162 j10=ny do j=j3+1,j4-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do do j=j5+1,j6-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do do j=j7+1,j8-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do do j=j9+1,j10-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do else j1=1 j2=ny endif do j=j1+1,j2-1 f(i,j)=-(vx(i,j)-uy(i,j)) sf(i,j)=0. end do end do c-----------------------------------------

155

c boundary condition c----------------------------------------do i=1,nx IF(i.le.101)then j1=1 j2=38 j3=42 j4=78 j5=82 j6=118 j7=122 j8=158 j9=162 j10=ny sf(i,j3)=0. sf(i,j4)=0. sf(i,j5)=0. sf(i,j6)=0. sf(i,j7)=0. sf(i,j8)=0. sf(i,j9)=0. sf(i,j10)=0. else j1=1 j2=ny endif sf(i,j1)=0. sf(i,j2)=0. end do do j=1,ny IF(j.ge.38.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101 elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif sf(i1,j)=0. sf(i2,j)=0. end do c-----------------------------------------

156

15

ERR=0. k=k+1 IF(k.gt.20000)stop c----------------------------------------c matric coefficient c----------------------------------------do j=2,ny-1 IF(j.ge.38.and.j.le.42)then i1=101 elseif(j.ge.78.and.j.le.82)then i1=101 elseif(j.ge.118.and.j.le.122)then i1=101 elseif(j.ge.158.and.j.le.162)then i1=101 else i1=1 i2=nx endif do i=i1+1,i2-1 a(i)=1. b(i)=-4. c(i)=1. d(i)=-sf(i,j+1)-sf(i,j-1)+f(i,j)*dx*dx enddo c----------------------------------------c modification matric coefficien c----------------------------------------d(i1+1)=d(i1+1)-a(i1+1)*sf(i1,j) a(i1+1)=0. d(i2-1)=d(i2-1)-c(i2-1)*sf(i2,j) c(i2-1)=0. l1=i1+2 l2=i2-1 call tridi do i=i1+1,i2-1 ERR=ERR+ABS(d(i)-sf(i,j)) sf(i,j)=d(i) end do enddo WRITE(*,*)k,err IF(err.gt.ermx)GOTO 15 return end

157

Lampiran 13. Kode Program MATLAB  Program Visualisasi Isotermal load eko -ascii; load num -ascii; im=num(1); jm=num(2); imax=im*jm; for i=1:im for j=1:jm is=i+(j-1)*im; x(i,j)=eko(is); y(i,j)=eko(is+imax); u(i,j)=eko(is+2*imax); end end hold on; box on; contour(x,y,u,30); colorbar; hold off; clear  Program Visualisasi Steam function load streamf -ascii; load num -ascii; im=num(1); jm=num(2); imax=im*jm; for i=1:im for j=1:jm is=i+(j-1)*im; x(i,j)=streamf(is); y(i,j)=streamf(is+imax); u(i,j)=streamf(is+2*imax); end end hold on; box on; contour(x,y,u,40); colorbar; hold off; clear