BAB II TINJAUAN PUSTAKA

BAB II TINJAUAN PUSTAKA

2.1

Pendahuluan Minyak mentah yang dihasilkan oleh ladang minyak memiliki campuran

hidrokarbon yang kompleks mulai dari metana dan aspal. Fungsi dari destilasi minyak mentah adalah membuat fraksi minyak mentah menjadi hidrokarbon rendah, neftan/gasoline, kerosin, disel, dan minyak residu. Beberapa hasil tersebut dapat dijual ke pasar secara langsung, sedangkan yang lainnya membutuhkan proses lebih lanjut pada unit-unit kilang minyak agar dapat dijual. Kilang minyak (Oil refinery) merupakan proses yang dilakukan untuk mendapatkan bahan bakar yang memiliki nilai yang lebih tinggi. Ada berbagai macam unit-unit yang terdapat di oil refinery yang saling mendukung proses satu sama lainnya, beberapa diantaranya adalah delayed coker, desalter, penyaringan vakum, penyaringan atmosfir, hydrocracker unit,dll. Unit-unit tersebut saling terintegrasi satu sama lain sehingga apabila terjadi kerusakan tentu saja akan menghambat proses produksi karna semua unit saling berkaitan satu sama lain dan tentunya mengakibatkan kerugian perusahaan. Didalam bidang pertambangan dan perminyakan kerusakan sekecil apapun akan menghasilkan kerugian yang besar, oleh karena itu meramalkan kerusakan alat-alat pada oil refinery sangatlah penting untuk mencegah hal-hal tersebut. Delayed coker adalah sebuah tipe coker (alat pemasak) yang prosesnya memanaskan minyak residu hingga ke temperatur thermal cracking yang menghasilkan berbagai fasa. Delayed coking adalah satu dari unit proses yang digunakan pada banyak kilang minyak. Salah satu komponen alat pada delayed coker adalah coke drum, yaitu tempat dimana terjadinya thermal cracking minyak residu yang menghasilkan fraksi berat dan fraksi ringan. Coke drum merupakan drum berdinding tipis yang digunakan untuk mengolah minyak residu dari hasil destilasi sebelumnya sehingga menghasilkan minyak solar dan fuel oil lainnya sehingga minyak residu tersebut tidak terbuang. Coke drum sering mengalami kerusakan karena drum ini dedesain sesuai dengan ASME Boiler and Pressure vessel Codes Section VIII Division I. Sedangkan coke

6

Universitas Sumatera Utara

drum sendiri dalam pemakaiannya akan mengalami beban berulang dikarenakan proses pemanasan dan pendinginan yang membuat dinding drum mengalami pegembangan dan penyusutan akibat pengaruh temperatur. Berdasarkan survey API, kerusakan pada coke drum muncul setelah mengalami 3000 – 5000 siklus (16 - 27 tahun). Ilustrasi pengembangan dan penyusutan pada dinding coke drum dapat dilihat pada Gambar 2.1

Gambar 2.1. Proses cyclic load pada coke drum[2]

Beberapa istilah yang umum digunakan dalam operasi delayed coker sebagai berikut : 1. Steaming. Coke drum yang penuh di uapkan untuk menghilangkan sisa minyak residu cair. Campuran antara uap dan hidrokarbon dikirmkan ke pemisah (fractinator) 2. Cooling. Coke drum diisi dengan air, mengakibatkan temperatur dibawah 93 oC. 3. Draining. Air pendingin dikeringkan dari drum dan dipulihkan untuk digunakan kembali 4. Unheading. Persiapan untuk mengeluarkan coke dari drum 5. Decoking. Air bertekanan tinggi (2900 psig) digunakan untuk “cut” coke dari coke drum dikeluarkan menuju drain sump.

7


6. Heading and testing. Setelah kepala drum diganti, drum dibersihkan dan dilakukan uji temperatur 7. Heating. Uap dimasukkan kedalam coke drum sebagai pemanas awal dari coke drum yang dingin. Dalam satu siklus memerlukan waktu ± 48 jam. Proses sebagai berikut : minyak residu yang berasal dari unit destilasi vakum (terkadang juga termasuk minyak dengan titik didih yang tinggi dari berbagai sumber di kilang) dipompa kedalam fractinator yang berfungsi memisahkan antara fraksi berat dan fraksi ringan yang mana fraksi ringan (naphtha, light gas oil, heavy gas oil) akan disalurkan ke unit lainnya sedangkan fraksi berat masuk ke coke drum. Sebelum fluida dimasukkan ke coke drum, fluida tersebut dipanaskan di furnace (dapur) untuk mendapatkan thermal cracking dengan temperatur sekitar 482-485 oC. Thermal cracking adalah temperatur yang digunakan untuk mengubah fraksi yang sangat berat atau memproduksi fraksi ringan, bahan bakar minyak dan petroleum coke. Kondisi thermal cracking dibagi atas dua produk yaitu proses temperatur tinggi yang disebut “steam cracking” atau pyrolysis (temperatur proses 750 oC hingga 900 oC) yang mana memproduksi ethilen yang bernilai tinggi, serta pada proses temperatur menengah disebut delayed coking (temperatur 500 oC) yang memproduksi coke (batu bara), yaitu sebuah batubara yang digunakan dalam memproduksi elektroda untuk industri baja dan aluminium. Thermal crack umumnya terjadi di pipa antara furnace dan coke drum.Untuk membuat fluida mengalami thermal cracking ditambahkan uap air kedalam coke drum, sehingga thermal cracking baru akan terjadi ketika fluida tersebut berada didalam coke drum. Ketika thermal cracking terjadi didalam drum, gas oil dan komponen fraksi ringan dalam bentuk uap terpisah dari cairan dan solid, dan fraksi ringan tersebut dialirakan langsung menuju fractinator bercampur dengan minyak residu dari unit vakum destilasi yang kemudian kembali untuk dipisahkan berdasarkan titik didihnya fraksi tersebut. Setelah drum penuh oleh coke yang padat, coke yang masih bercampur dengan hidrokarbon lainnya dipindahkan ke drum selajutnya. Ketika proses pengisian drum selanjutnya, drum yang penuh mengalami pendinginan dengan air.

8


Material yang solid dikeluarkan dengan air bertekanan 2900 psig, sehinga coke drum jatuh ke tempat penampungan coke petroleum. Coke petroleum dapat diklasifikasi kedalam dua katagori, sponge coke ,needle coke, dan Shot coke bergantung pada sifat fisiknya seperti tekstur, densitas, porositas, ketahanan elektrik dan koeffisien konduktifitas termal. Berbagai produk coke petroleum dapat dilihat pada Gambar 2.2 1. Sponge coke Merupakan kualitas coke petroleum standar. Merupakan pencampuran antara needle. Industri aluminum merupakan penggunaan terbesar untuk jenis coke ini, dimana digunakan untuk membuat anode. Secara kasar 0,40,5 ton anode dikonsumsi per ton dari produksi aluminum. Karakteristik dari sponge coke adalah kandungan sulfur 0,5-3 %. Walaupun sulfur dalam

peleburan

petroleum

coke

meningkatkan

performa

anode

dikarenakan mencegah dampak negatif reaksi udara dan karbon dioksida, kandungan sulfur yang diperbolehkan dalam banyak kasus ditentukan berdasarkan peraturan lingkungan pada lokasi dimana peleburan berada. 2. Needle coke Needle coke digunakan dalam produksi batang grafit diameter besar (2428 in) untuk elektroda dapur berkekuatan tinggi. Material yang digunakan memiliki densitas yang tinggi, ketahanan yang rendah, kekuatan tinggi, dan koeffisien ekspansi termal yang rendah (kurang dari 2 x 10-7/OC). Needle coke tidak mengandung aspal dengan ciri-ciri warna silver abu-abu. 3. Shot coke Shot coke diproduksi akibat konsentrasi aspal yang tinggi, kecepatan dan turbulensi dalam coke drum, serta coke drum dengan temperatur yang tinggi. Shot coke berbentuk kecil bundar dengan ukuran diameter 1,5 hingga 5 mm, Shot coke tidak bisa digunakan dalam membuat aluminum anoda karena memiliki ekspansi thermal yang rendah, bersifat isotropik, dan ekspansi termal yang tinggi. Shot coke memiliki kandungan sulfur min. 4% dengan kandungan karbon min. 89%

9


(a)

(b)

(c)

Gambar 2.2 (a) needle coke (b) sponge coke (c) shot coke

Pada tahap awal proses kerja coke drum disebut preheating stage yaitu dengan menginjeksikan uap panas 3500C, kemudian filling stage dengan menginjeksikan feed material (minyak residu/fraksi berat) dengan temperatur 4400C hingga 5000C selama 460-2000 min dimana tekanan didalam coke drum berkisar 300 hingga 350 kPa. Pada tahap ketiga adalah proses pendinginan dengan menginjeksikan uap air dan air selama 2000 – 2400 min, proses pengangkatan hidrokarbon

yang

terperangkap

dalam

fraksi

berat

dilakukan

dengan

menginjeksikan uap air, sedangkan proses pendinginan air untuk mendinginkan fraksi berat dan dinding coke drum dengan aliran massa air hi ngga 63 kg/s. Dan tahap akhir yaitu 2400-2800 min yaitu tahap cutting dan deheading, jika coke drum telah mengalami beberapa proses tersebut maka dapat dikatakan coke drum telah mengalami satu siklus. Dalam penelitian kali ini akan dilakukan analisa coke drum pada bagian sambungan shell skirt yang sebelumnya telah dilakukan perhitungan temperatur dan regangan pada beberapa bagian di coke drum secara manual dengan menggunakan termokopel dan stress gauge. Alur delayed coker dapat dilihat pada Gambar 2.3 sedangkan ilustrasi pengangkatan hidrokarbon yang terperangkap dapat dilihat pada Gambar 2.4

10


Gambar 2.3 Delayed cooker unit[12]

Flange

Fraksi ringan menuju fraksional tower gelembung

3way valve menuj u fraksio nal

Proses cutting dengan menggunakan air bertekanan 2900 psig

Jalur fluida menuju fraksional tower

Coke petroleum

Coke petroleum keluar drum

Flange Fluida masuk

Gambar 2.4 Pengangkatan sisa hidrokarbon ke fraksional tower dan pengeluaran coke dari drum[1]

11


Material yang umum dipakai pada coke drum ada beberapa macam, material tersebut dapat dilihat pada tabel 2.1 Tabel 2.1. Sifat material coke drum [6] Base metal

Jenis

Tarikan

Yield min.

(ksi)

(ksi)

A516-70

Carbon steel

70-90

38

A204 C

Carbon – ½ Moly

75-95

43

A387-11, CL2

1 1/4 Cr – ½ Moly

75-100

45

A387-12,CL2

1 Cr – ½ Moly

65-85

40

A387-22, CL2

2¼ - 1 Moly

75-100

45

405 (clad

13 Cr

60-90

25

410S (clad)

12 Cr

60-90

30

Jenis material coke drum yang digunakan adalah SA-387 Gr.11 CL.2 dengan komposisi 1¼ Cr – ½ Mo – Si. Poison rasio sebesar 0,3, dengan densitas material coke drum 7850 kg/m3.

2.2

Kelelahan Bersiklus Rendah Suatu kelelahan dikatakan bersiklus rendah apabila suatu benda memiliki

ketahanan diberi tegangan bolak balik sampai kira-kira 104 siklus. Pengatahuan tahanan kelelahan pada daerah bersiklus rendah ini diperlukan untuk perencanaan alat-alat berumur pendek atau kemungkinan menerima beberapa beban lebih yang besar selama umurnya. Kegagalan rendah hampir selalu dimulai pada suatu ketidak mulusan setempat seperti sebuah takikan, retak, bidang pemusatan tegangan yang tinggi lainnya. Bila tegangan pada ketidak mulusan tersebut melampaui batas elastis,

12


regangan plastis terjadi. Kalau kepatahan karena lelah terjadi, berarti disana terjadi regangan plastis berulang-ulang Pada tahun 1910 melalui percobaan Bairstow membuktikan teori Bauschinger bahwa batas elastis dari besi dan baja bisa diubah, naik atau turun, dengan pengulangan variasi tegangan. Pada umunya, batas elastis dari baja anil (annealed steel) akan naik apabila diberi pengulangan tegangan yang bolak-balik, sementara baja tarik dingin menujukkan penurunan batas elastis. Benda percobaan untuk beban lentur yang bolak-balik tidak cocok untuk regangan berulang karena sulitnya mengukur regangan plastis. Akibatnya, penelitian yang pernah dilakukan kebanyakan menggunakan benda percobaan untuk beban aksial. Dengan menggunakan transducer listrik adalah mungkin untuk membangkitkan sinyal yang berbanding lurus, masing-masing, terhadap tegangan dan regangan. Sinyal-sinyal ini kemudian dapat diperlihatkan pada osiloskop atau digambarkan pada plotter XY. R.W. Landgraf telah meneliti perilaku dari kelelahan bersiklus rendah pada baja berkukatan tinggi dalam jumlah yang banyak dan selama penelitian tersebut telah banyak membuat gambar teganan regagangan berulang. Gambar 2.5 dibuat untuk memperlihatkan penampilan secara umum dari gambar-gambar untuk beberapa siklus pendahuluan dari regangan berulang terkendali. Dalam kasus ini kekuatan menurun terhadap pengulangan tegangan, seperti dibuktikan oleh kenyataan bahwa pembalikan juga terjadi pada tingkat tegangan yang lebih rendah. Sebagaimana yang sebelumnya dicatat, bahan-bahan lain mungkin bisa menjadi bertambah kuat karena pengulangan tegangan pembalikan tersebut. Hasil yang sedikit berbeda mungkin didapat apabila pembalikan awal terjadi pada daerah tegangan tekan; hal ini mungkin karena pengaruh penguatan kelelahan oleh tegangan tekan. Laporan landgrave mengandung sejumlah gambar yang membandingkan hubungan tegangan regangan yang monoton baik untuk tarikan maupun untuk tekanan dengan kurva tegangan regangan berulang. Pentingnya gambar ini adalah bahwa ini mempertegas betapa sulitnya mencoba menaksir kekuatan lelah dari bahan yang harga kekuatan mengalah (yield stress) atau kekuatan akhirnya (ultimate stress) yang monoton diketahui dalam daerah bersiklus rendah

13


σ

A

Δσ

ε

Ke 4 B

Ke 2

Δε e

Δε p

Δε

Gambar 2.5 Lingkaran histeris tegangan dan regangan[9]

Fatigue Design dan Evaluaton Steering Committee dari SAE tahun 1975 mengeluarkan laporan dimana umur dalam beberapa kali perubahan sampai gagal dikaitkan dengan amplitudo regangan. Laporan tersebut berisi sebuah gambar dari hubungan ini untuk baja SAE 1020 yang dirol panas; grafik ini dicetak kembali pada Gambar 2.6. . Untuk menjelaskan Gambar 2.6 pertama-tama kita menetapkan beberapa istilah berikut 100

ε'F Amplitudo regangan, Δε/2

10

-1

c 1.0

10-2

σ 'F E

Regangan total Regangan plastis b

10

1.0

-3

10-4 100

Regangan elastis

101

102

103

104

105

106

Pengulangan sampai gagal, 2N

Gambar 2.6 log-log menunjukkan bagaimana umur kelelahan dikaitkan dengan amplitudo regangan [9]

14


Koeffisien daktilitas kelelahan (fatigue ductility coefficient) ε’p adalah regangan yang berkaitan dengan kepatahan pada satu pembalikan (titik A pada Gambar 2.5). Koeffisien kekuatan kelelahan (Fatigue strength coefficient) σ’F adalah tegangan sebenarnya yang berkaitan dengan kepatahan pada satu pembalikan (titik A pada Gambar 2.5) perhatikan pada Gambar 2.6 bahwa garis regangan elastis mulai pada σ’F/E Eksponen daktilitas kelelahan (Fatigue ductility exponent) c adalah kemiringan garis regangan plastis pada Gambar 2.6 dan adalah pangkat untuk menaikkan umur 2N agar berbanding lurus dengan amplitudo regangan plastis sebenarnya. Eksponen kekuatan kelelahan (Fatigue strength exponent) b adalah kemiringan garis regangan elastis, dan adalah pangkat yang diberikan untuk menaikkan umur 2N agar berbanding lurus dengan amplitudo tegangan sebenarnya. Pada Gambar 2.6, dapat dilihat bahwa regangan total adalah jumlah dari komponen elastis dan plastis. Maka amplitudo regangan total adalah : Δε Δε e Δε p = + 2 2 2

(2.1)

Persamaan garis regangan plastis dalam Gambar 2.7 adalah Δε p 2

= ε' F (2 N ) c

(2.2)

Persamaan untuk garis regangan elastis adalah Δε e ε ' F = (2 N ) b 2 E

(2.3)

Karena itu, dari persamaan (2.1), untuk amplitudo regangan total didapat : Δε σ 'F = ( 2 N )b + ε 'F ( 2 N ) c 2 E

(2.4)

15


yaitu persamaan Manson-Coffin antara umur kelelahan dan regagan total. Beberapa harga koeffisien dan eksponen tersebut didaftarkan pada laporan SAE J1099. Walaupun persamaan (2.4) adalah suatu persamaan sempurna yang sah untuk umur kelelahan dari suatu bagian bila regangan dan perilaku pengulangan diberikan, ternyata hanya sedikit pemakaiannya dalam perencanaan. Ada kemungkinan faktor pemusatan regangan akan bisa didapat pada buku-buku hasil penelitian menggunakan analisa elemen hingga. Analisa elemen hingga tersendiri dapat menaksir regangan yang akan terjadi pada semua titik pada struktur yang dimaksud. Siklus rentang regangan yang terjadi akibat beban berulang dapat dilihat pada Gambar 2.7, dimana pada setiap siklus dapat dihitung menggunakan persamaan : Δε = Δεmax – Δεmin

(2.5) εmax

regangan

Δεmax Waktu pembebanan

Δεmin εmin

Gambar 2.7 Siklus nilai regangan terhadap waktu pada pembebanan berulang[9]

Coke drum beroperasi dalam kategori produksi dengan temperatur tinggi dengan siklus kelelahan rendah. Diagram kelelahan desain atau rentan regangan dibanding jumlah siklus hingga rusak digunakan untuk menghitung umur kelelahan. Material coke drum adalah SCMV3 (1 1/4Cr – 1/2Mo steel (JIS G 4109)). Persamaan Coffin-Manson digunakan

Δε =

Cp N

α

+

Ce Nβ

(2.6)

16


Pada persamaan (2.6) Δε adalah rentan regangan dan N adalah jumlah siklus hingga rusak. Dimana konstanta yang digunakan berdasarkan material dan temperatur operasional 500 0C adalah Cp= 1,4 , α= 0,85 , Ce= 0,0064, β= 0,1 Untuk hubungan antara switching temperatur dan jumlah siklus hingga rusak, persamaan umur pakai pada sambungan shell skirt adalah sebuah fungsi switching temperatur yang telah dikembangkan. Umur pakai dapat dihitung menggunakan persamaan N= AeBT

(2.7)

Persamaan (2.7) diperoleh dari hasil eksperimen dengan menentukan persamaan garis dari beberapa selisih regangan ΔS pada proses cooking terhadap umur pakai N. Keofisien A dan B adalah 0,6564 dan 0,0249. T adalah switching temperatur dalam derajat Celsius. Persamaan diatas hanya valid untuk subjek dengan dimensi dan karakteristik operasi yang digunakan pada analisa ini.

2.3

Tegangan dan Regangan Termal Beban luar bukanlah satu-satunya sumber tegangan dan regangan di suatu

struktur. Perubahan temperatur menyebabkan ekspansi atau kontraksi bahan, sehingga terjadi regangan termal dan tegangan termal. Pada kebanyakan bahan, regangan termal εt sebanding dengan perubahan temperatur ΔT; jadi , ε x = ε y = ε z = αΔT

(2.8)

Dimana α adalah koeffisien perpanjangan termal dan ΔT adalah perubahan suhu, dalam derajat (o). Dalam hal ini benda mengalami sedikit perubahan volume dimana semua komponen regangan geser sama dengan nol. Jika sebuah batang yang ditahan untuk menahan pertambahan panjang dan karena adanya kenaikan suhu yang merata, akan menyebabkan adanya tegangan tekan akibat desakan arah aksial. Tegangan ini adalah :

σ = εE = α.ΔT .E

(2.9)

17


Dengan cara yang sama, untuk sebuah pelat jika ditahan sisi-sinya dan diberi kenaikan suhu secara merata, tegangan tekan yang timbul dinyatakan σ=

α.ΔT .E 1 μ

(2.10)

Walaupun dikarenakan oleh kenaikan suhu, bukanlah tegangan termal karena tegangan-tegangan itu terjadi akibat penahan pada ujung-ujungnya. Tegangan-tegangan dalam tegangan yang timbul karena adanya perbedaan suhu dalam sebuah benda. Jika material yang mengalami efek termal diasumsikan bahwa bahannya isotropis dan homogen, dan bahwa peningkatan temperatur ΔT adalah seragam di seluruh blok, seperti yang terlihat pada Gambar 2.8

Δl

l

Gambar 2.8 Ekspansi yang terjadi pada elemen akibat beban termal[9]

Kita dapat menghitung bertambahnya dimensi manapun dari blok ini dengan mengalikan dimensi semula dengan regangan termal. Jika salah satu dimensi adalah L, maka dimensi tersebut akan bertambah sebesar ∆L = α(ΔT ) L ∆L = Lα(ΔT )

εT =

(2.11)

Selama pendinginan, tegangan maksimum adalah di permukaan yang berupa tegangan tarik. Pada waktu yang sama, keseimbangan gaya memerlukan tegangan tekan dibagian tengah pelat. Selama pendinginan, permukaan luar panas dan cenderung memuai tetapi ditahan oleh bagian dalam yang dingin. Hal ini menyebabkan tegangan tekan di permukaan dan tegangan tarik disebelah dalam.

18


Bagian-bagian struktur sering diperlukan untuk bisa tahan terhadap suhu, dimana biasanya sifat mekanis bahan berubah. Pada Gambar 2.9 dan 2.10 dapat dilihat perubahan sifat bahan yang timbul akibat suhu.

Gambar 2.9. Pengaruh suhu terhadap nilai tumbukan[9]

Gambar 2.10. Pengaruh nilai perpanjangan pada bahan tarik[9]

Banyak pengujian telah dilakukan pada logam ferro yang diberi beban yang tetap untuk waktu yang lama pada suhu yang dinaikkan. Didapat bahwa benda percobaan ternyata mengalami perubahan yang permanen selama percobaan, walaupun pada saat itu tegangan sebenarnya adalah lebih kecil dari kekuatan yield bahan yang didapat dari percobaan dalam waktu yang tingkat pada suhu yang sama. Deformasi yang terus menerus berkelanjutan karena beban disebut perubahan perlahan-lahan (creep).

19


2.4

Konveksi Paksa Konveksi paksa merupakan konveksi yang diakibatkan oleh fluida yang

terdapat pada permukaan pelat. Aliran fluida pada pelat dengan panjang L pada suatu arah aliran seperti yang ditunjukkan pada Gambar 2.11 Turbulen

T∞ V

Laminar

Ts Xcr L

Gambar 2.11 Daerah batas laminar dan turbulen suatu aliran pada pelat[3]

Koordinat x dihitung sepanjang permukaan pelat dari sisi terdepan pada arah aliran. Fluida mengenai permukaan pelat dalam arah x dengan kecepatan v dan temperatur T∞ yang seragam. Awalnya kecepatan bermula dengan batas aliran laminar , tetapi jika pelat cukup panjang, aliran menjadi turbulen pada jarak xcr dari permukaan depan dimana bilangan Reynold memperoleh nilai kritis untuk daerah transisi. Transisi dari aliran laminar ke turbulen bergantung pada geometri permukaan, kekasaran permukaan, kecepatan, temperatur permukaan, jenis fluida dan lainnya yang menjadi karakter penentu bilangan Reynold. Bilangan Reynold pada jarak x dari sisi terdepan pelat datar dinyatakan sebagai Re x =

ρVx Vx = µ v

(2.12)

Dicatat bahwa nilai bilangan Reynold bervariasi pada sebuah pelat datar sepanjang aliran. Untuk aliran transisi dari laminar ke turbulen diperoleh dengan persamaan Re cr =

ρVx cr = 5 x10 5 µ

(2.13)

20


Nilai bilangan Reynold dari sebuah pelat datar bisanya bervariasi antara 105 hingga 3x106, bergantung pada kekasaran permukaan. Bilangan Nusselt lokal pada lokasi x untuk aliran laminar sepanjang pelat datar ditentukan dengan turunan persamaan energi yaitu Nu x =

hx x = 0.332 Re 0x.5 Pr 1 / 3 k

(2.14)

Sedangkan untuk aliran turbulen adalah Nu x =

hx x = 0.0296 Re 0x.8 Pr 1 / 3 k

(2.15)

Persamaan laminar digunakan apabila bilangan Re dibawah 5x105 dengan nilai Pr>0.6. Hubungan antara koeffisien rata-rata perpindahan panas terhadap jenis aliran dapat dilihat pada Gambar 2.12. hx, turbulen h rata-rata Turbulen hx, laminar

Laminar

Ts Xcr

Gambar 2.12 Grafik yang menunjukkan koeffisien perpindahan panas rata-rata [3]

untuk pelat datar dengan campuran antara aliran laminar dan turbulen

2.5

Konveksi Bebas Konveksi alamiah (natural convection), atau konveksi bebas (free

convection), terjadi karena fluida yang, karena proses pemanasan, berubah densitasnya (kerapatannya), dan bergerak naik. Gerakan fluida dalam konveksi bebas, baik fluida itu gas maupun zat cair, terjadi karena gaya apung (buoyancy force) yang dialaminya apabila densitas fluida didekat permukaan perpindahankalor berkurang sebagai akibat proses pemanasan. Gaya apung itu tidak akan

21


terjadi apabila fluida itu tidak mengalami suatu gaya dari luar seperti gravitasi (gaya berat), walaupun gravitasi bukanlah satu-satunya medan gaya luar yang dapat menghasilkan arus konveksi-bebas; fluida yang terkurung dalam mesin rotasi mengalami medan gaya sentrifugal, dan karena itu mengalami arus konveksi bebas bila salah satu atau beberapa permukannya yang dalam kontak dengan fluida itu dipanaskan. Gaya apung yang menyebabkan arus konveksi bebas disebut gaya badan (body force). Untuk menentukan laju perpindahan panas konveksi, dapat dilakukan berdasarkan teori-teori perpindahan panas

h.L = CRa n k

Nu L =

(2.16)

Disini, NuD, h, L dan k adalah Nusselt number, koeffisien perpindahan panas (W/m. oC), dimensi dari struktur dan konduktifitas termal (W/m. oC). nilai C dan n berbeda beda dari tiap kasus, nilai C kurang dari 1 dan n memiliki nilai ¼ untuk aliran laminar dan 1/3 untuk aliran turbulen. Ra adalah bilangan Rayleigh yang ditentukan dari hasil perhitungan Grashof dan Prandtl. (2.17)

Ra = GrL . Pr

Nilai rata-rata koeffisien natural konveksi perpindahan panasnya didapat dari pengembangan persamaan umum untuk konveksi bebas : GrL =

Pr =

gβ(Ts - T∞ ) L3 v2

(2.18)

c p .μ

(2.19)

k

Dimana g

= percepatan gravitasi, m/s2

β

= koeffisien ekspansi volume, 1/K (β = 1/T untuk gas ideal)

cp

= panas spesifik, kJ/kg. oC

Ts

= temperatur permukaan, oC

T∞

= temperatur fluida saat jauh dari permukaan (ruangan), oC

L

= dimensi geometri (dimeter), m

v

= viskositas kinematik fluida, m2/s

22


Bilangan Grashof merupakan bilangan tak berdimensi yang menggambarkan rasio gaya apung (buoyancy force) terhadap gaya kekentalan (viscous force) yang berkerja pada fluida. Kemudian, koeffisien perpindahan panas konveksi bebas dapat diperoleh dari hubungan antara persamaan (2.16) dan (2.17) : C.k h= L

  ρ 2 .g .β  2 . Pr .ΔT    μ

n

(2.20)

Maka besarnya konveksi perpindahan panas pada permukaan pelat adalah : .

Q conv = hAs (Ts - T∞ )

(2.21)

Dimana As adalah luas permukaan perpindahan panas (m2) dan Qconv dengan satuan Watt. Pada persamaan (2.16) dapat dikembangkan khusus untuk pelat vertikal. Hubungan yang dapat dipakai untuk seluruh nilai Ra direkomendasikan oleh Churchill dan Chu dengan bentuk persamaan 1/ 6  0.387 Ra L Nu L = 0.825 + 9 / 16  1 + (0.492 / Pr )

[

 8 / 27  

2

]

(2.22)

Persamaan (2.22) lebih kompleks tetapi menghasilkan nilai yang lebih akurat. Untuk pelat vertikal miring maka nilai g dalam persamaan (2.18) adalah g cos θ, dimana nilai θ adalah sudut kemiringan dari sumbu y. 2.6

Perpindahan Panas Konduksi pada Dinding Perpindahan panas konduksi pada dinding bergerak dengan arah normal ke

permukaan dinding dan tidak ada perubahan perpindahan panas yang signifikan dari arah lainnya. Perpindahan panas dengan arah tersebut digerakkan oleh gradien temperatur. Tidak akan ada perpindahan panas dalam suatu arah jika tidak ada perubahan temperatur didalamnya. Menghitung temperatur pada beebrapa lokasi di bagian dalam dan luar dinding dilakukan untuk memastikan permukaan dinding

23


isotermal, begitu juga pada bagian atas dan bawah permukaan dinding. Sehingga tidak akan ada perpindahan panas yang melewati dinding baik pada dari atas ke bawah ataupun dari kiri ke kaan dinding., tetapi jika diperhitungkan perbedaan temperatur antara bagian dalam dan luar dinding , akan muncul perpindahan panas yang signifikan dalam arah dari bagian dalam ke luar. Untuk dinding tipis akan mengakibatkan gradien temperatur akan menjadi besar. Jika temperatur udara di dalam dan luar dinding dijaga konstan. Maka perpindahan panas yang melewati dinding dapat dimodelkan secara steady dan satu dimensi. Jika tidak adanya panas yang ditimbulkan paad dinding , maka keseimbangan energi pada dinding dapat dinyatakan sebagai

Atau .

.

Q in − Q out =

dE dinding dt

(2.23)

Untuk kondisi steady dEdinding/dt=0, dimana tidak ada perubahan temperatur pada setiap titik terhadap waktu. Sehingga laju perpindahan panas kedalam dinding dama dengan laju perpindahan panas ke luar dinding. Dalam kata lain, laju perpindahan panas melewati dinding harus konstan. Distribusi temperatur pada dinding dapat dilihat pada Gambar 2.13

24


Q

dT T2

dx 0

L

x

Gambar 2.13 Distribusi temperatur pada dinding dengan arah garis lurus[3]

Dengan

mempertimbangkan

ketebalan

dinding

L

dan

rata-rata

konduktifitas termal K. Kedua permukaan dinding dijaga konstan temperatur nya T1 dan T2. Hukum Fourier untuk perpindahan panas konduksi untuk dinding dapat dinyatakan sebagai .

Q konduksi , dinding = −kA

dT dx

(2.24)

Dimana laju perpindahan panas konduksi dan luas dinding A konstan. Jika persamaan (2.24) dibaut dalam bentuk integral untuk x=0 dimana T(0)=T1, untuk x=L, dimana T(L) = T2, diperoleh . L

∫

x =0

dx =

Q konduksi , dinding

T2

∫

− kAdT

T = T1

Dengan melakukan integral diperoleh hasilnya .

Q konduksi ,dinding = − kA

T2 − T1 T −T = kA 1 2 L L

(2.25)

25


2.6.1

Konsep Ketahanan Termal Persamaan (2.25) untuk perpindahan panas konduksi melewati sebuah

dinding dapat ditulis ulang menjadi .

Q konduksi , dinding =

T1 − T2 R wall

(2.26)

(oC/W)

(2.27)

Dimana R wall =

L kA

Yang merupakan ketahanan termal pada dinding terhadap perpindahan panas konduksi atau disebut ketahanan konduksi pada dinding. Dimana ketahanan termal suatu medium dipengaruhi oleh gemometri dan sifat termal medium. Dengan mempertimbangkan perpindahan panas konveksi dari permukaan solid area As, dan temperatur Ts serta koeffisien perpindahan panas konveksi h. Persamaan

Newton

untuk

laju

perpindahan

panas

konveksi

.

Q konveksi = hAs (Ts − T∞ ) dapat ditulis ulang sebagai .

Q konveksi =

(Ts − T∞ ) R konveksi

(2.28)

1 hAs

(2.29)

Dimana R konveksi =

Adalah ketahanan termal pada permukaan terhadap panas konveksi, atau disebut ketahanan konveksi pada permukaan. Dengan catatan jika koeffisien perpindahan panas konveksi sangat besar, ketahanan konveksi menjadi nol dan Ts~T ∞. Sehingga tidak akan ada ketahanan konveksi dan tidak akan memperlambat proses perpindahan panas.

26


2.6.2 Ketahanan termal dinding berlapis Dalam keseharian sering dijumpai dinding dengan beberapa lapis dengan material yang berbeda. Konsep ketahanan termal masih tetap digunakan untuk menentukan laju perpindahan panas.

T∞1

Dinding 1

Dinding 2

T1 T2 T∞2

Rkonfeksi,1

R1

R2

Rkonfeksi,2

Gambar 2.14 Ketahanan termal dengan dinding dua lapis[3]

Jika dilihat dari Gambar 2.14, merupakan dinding yang memiliki dua lapis dinding. Laju perpindahan panas yang melewati dua lapis dinding dapat dinyatakan .

Q=

(T∞1 − T∞ 2 ) Rtotal

(2.30)

Dimana Rtotal adalah ketahanan termal total, dinyatakan sebagai

Rtotal = Rconv ,1 + R wall ,1 + R wall , 2 + Rconv ,1 =

L L 1 1 + 1 + 2 + h1 A k1 A k2 A h2 A

(2.31)

Subskrip 1 dan 2 dalam Rwall mengindikasikan lapis pertama dan kedua. Dari hubungan total ketahanan termal dilihat bahwa ketahanan berbentuk seri, sehingga ketahanan termal total secara arimatik sederhana langsung dijumlahkan. Hasil dari

27


kasus dua lapis dapat dianalogikan sebagai kasus satu kasus. Hasil ini juga dapat digunakan untuk dinding dengan lebih banyak lapisan dengan menambahkan tambahan tahanan untuk setiap penambahan lapisan. Sedangkan untuk tahanan termal radiasi dapat ditambahkan nilai 1/hradiasi.A, dimana nilai h

radiasi

adalah

1/εσ(Tsur3- Tsur3)(Tsur-Tsur). Nilai ε adalah emisivitas benda dan σ adalah konstanta Stefan Boltzmann yang nilainya 5,67x10-8.

2.7

Perhitungan Dinamika Fluida Cumputational Fluid Dynamic (CFD) merupakan cabang dari mekanika

fluida yang menggunakan metode numerik dan algoritma untuk memecahkan dan menganalisis masalah yang melibatkan aliran fluida. Komputer digunakan untuk melakukan perhitungan yang diperlukan untuk mensimulasikan interaksi antara cairan-cairan dan gas-gas terhadap permukaan yang didefinisikan sebagai kondisi batas. Beberapa aplikasi dibidang industri dan non industri yang berhubungan dengan CFD adalah aerodinamis pesawat dan kendaraan, motor bakar dan turbin gas, meteorologi hingga biomedical engineering Kode-kode CFD tersusun pada algoritma numerik untuk menyelesaikan masalah aliran fluida. Karenanya semua kode-kode mengandung tiga elemen: (i) pre-processor, (ii) solver, (iii) post processor. Dimana penjabaran dari setiap fungsi elemen-elemen ini dijabarkan sebagai berikut :

1. Pre-processor Merupakan kumpulan data-data yang diketahui dari masalah aliran fluida ke program CFD. Aktivitas pemakai pada tahap pre processing seperti mendifinisikan

geometri

yang

dipakai

(computational

domain),

menentukan grid (mesh), memilih fenomena fisik dan kimia yang dibutuhkan untuk dimodelkan serta sifat fluida Solusi dari sebuah masalah aliran fluida (kecepatan, tekanan, temperatur, dsb) didefiniskan pada node didalam setiap sel sehingga disebut Volume Element Method (VEM). Akurasi solusi CFD ditentukan jumlah sel pada grid. Secara umum, semakin banyak jumlah sell semakin baik akurasi penyelesaiannya.

28


2. Solver Finite element method merupakan formulasi diferensial terbatas yang stabil pada penyelesaian berbagai masalah CFD, empat dari lima kode CFD komersial antara lain : PHOENICS, FLUENT,

FLOW3D dan

STAR-CD. 3. Post processor Dalam tahap post processor merupakan elemen untuk menampilkan grafik dan menunjukkan hasil-hasil yang dapat dilihat secara visual. Termasuk didalamnya adalah : tampilan geometri dan grid, tampilan vektor, jejak partikel, fasa fluida,dsb.

Fasilitas ini juga termasuk animasi untuk

tampilan hasil masalah dynamik dan fasilitas mengekspor data untuk diproses ketahap lainnya.

2.8

Persamaan Umum Untuk Aliran Fluida dan Perpindahan Panas Persamaan umum untuk aliran fluida menggambarkan persamaan

matematika dari hukum konservasi fisika yaitu : a. Konservasi massa b. laju perubahan momentum sama dengan total gaya pada partikel fluida c. laju perubahan energi sama dengan total laju penambahan panas dan laju kinerja yang dilakukan partikel fluida

T

N

S (x,y,z)

δz

B

δy

W

E

z

δx

y x

Gambar 2.15 Element fluida[11]

Di keenam permukaan di sebut sebagai N,S,E,W,T,B yang mana merupakan simbol dari North (utara), South (selatan), East (timur), West (barat). Arah positif sepanjang sumbu co-ordinat juga digunakan. Dari Gambar 2.15 dapat

29


dilihat bahwa pusat dari elemen terletak pada posisi (x,y,z). Sebuah perhitungan sistematik dilakukan berupa perubahan massa, momentum, dan energi dari elemen fluida seiring dengan aliran fluida melewati batas akan membuat pergerakan pada bagian dalam elemen, yang mengacu pada persamaan airan fluida. Semua sifat fluida merupakan fungsi dari jarak dan waktu sehingga kita dapat menulisnya ρ(x,y,z,t), p(x,y,z,t), T(x,y,z,t) dan u (x,y,z,t) untuk vektor densitas, tekanan, temperatur dan kecepatan. Sifat-sifat pada salah satu permukaan dapat dinyatakan sebagai persamaan Taylor. Misal tekanan pada permukaan E dan W , yang mana keduanya berjarak 1/2δx dari pusat elemen, dapat dinyatkaan sebagai : p-

2.9

∂p 1 ∂p 1 δx δx dan p + ∂x 2 ∂x 2

Konservasi Massa Langkah pertama dalam derivasi persamaan konservasi massa adalah

menulis keseimbangan massa pada elemen fluida

Laju peningkatan masa dalam elemen fluida = jumlah laju aliran massa kedalam elemen fluida

Laju peningkatan massa kedalam elemen fluida adalah : ∂ ∂ ρ (ρ∂ x∂ y∂ z ) = ∂ x∂ y∂ z ∂t ∂t

(2.32)

Selanjutnya dapat ditentukan jumlah aliran massa pada elemen terhadap kecepatan dan densitas. Pada Gambar 2.16 dapat dilihat laju aliran massa kedalam elemen fluida yang melewati batas dinyatakan sebagai

 ∂ (ρu ) 1  ∂ (ρu ) 1  ∂ (ρv) 1    ∂x ∂y∂z -  ρu + ∂x ∂y∂z +  ρv ∂y ∂x∂z  ρu ∂x 2  ∂x 2  ∂y 2      ∂ (ρv) 1  ∂ (ρw) 1  ∂ (ρw) 1    -  ρv + ∂y ∂x∂z +  ρw ∂z ∂x∂y -  ρw + ∂z ∂x∂y ∂y 2  ∂z 2  ∂z 2    

30

(2.33)


Aliran-aliran yang mana menuju elemen menghasilkan peningkatan massa kedalam elemen dan mendapatkan tanda positif dan aliran-aliran yang meninggalkan elemen mendapatkan tanda negatif.

ρv +

ρw +

∂ (ρv) 1 δy ∂y 2

∂ (ρw) 1 δz ∂z 2

ρu +

(x,y,z)

∂ (ρu ) 1 ρu δx ∂x 2

∂ (ρu ) 1 δx ∂x 2

ρv -

z

y ρw -

x

∂ (ρv) 1 δy ∂y 2

∂ (ρw) 1 δz ∂z 2

Gambar 2.16 Aliran massa masuk dan keluar elemen fluida[11]

Laju peningkatan massa didalam elemen (pers 2.32) disamakan dengan jumlah aliran massa kedalam elemen yang melewati permukaan (pers 2.33). Semua kondisi menghasilkan keseimbangan massa yang diatur pada sisi kiri persamaan dari tanda keseimbangan dan dibagi oleh elemen volume δxδyδz. Ini menghasilkan persamaan

∂ρ ∂ (ρu ) ∂ (ρv) ∂ (ρw) + + + =0 ∂t ∂x ∂y ∂z

(2.34)

Atau dalam bentuk vektor dituliskan :

∂ρ +div(ρu ) = 0 ∂t

(2.35)

Untuk persamaan (2.35) untuk kondisi transien, tiga dimensi konservasi massa atau pesamaan kontinitas pada suatu titik dalam sebuah fluida mampumampat (incompresible). Kondisi kedua dideskribsikan sebagai jumlah aliran massa yang keluar dari elemen melewati batas dan disebut kondisi konvektif Untuk aliran tidak mampumampat (incompressible) (contoh sebuah cairan) densitas ρ bernilai konstan dan persamaan (2.35) menjadi

div.u = 0

(2.36)

31


Atau dalam notasi yang lebih panjang

∂u ∂v ∂w + + =0 ∂x ∂y ∂z

2.10

(2.37)

Laju Perubahan pada Partikel Fluida dan Elemen Fluida Hukum konservatif momentum dan energi membuat persetujuan akan

perubahan sifat sebuah partikel fluida. Setiap sifat pada sebuah partikel adalah sebuah fungsi dari posisi (x,y,z) dari partikel dan waktu t. Nilai suatu sifat per unit massa dinotasikan sebagai ϕ. Turunan total dari ϕ terhadap waktu pada partikel fluida, dituliskan sebagai Dϕ/Dt adalah

D ∂φ ∂ ∂ x ∂ ∂y ∂ ∂z = + + + Dt ∂ t ∂ x ∂ t ∂ y ∂ t ∂ z ∂ t

(2.38)

Sebuah partikel fluida akan mengalir, sehingga dx/dt=u, dy/dt=v dan dz/dt=w. oleh karena itu turunana subtantif dari ϕ dinyatakan sebagai :

Dφ ∂ φ ∂ φ ∂ φ ∂ φ = + u + v + w = u.grad φ Dt ∂t ∂x ∂y ∂z

(2.39)

Dϕ/Dt didefinisikan laju perubahan sifat ϕ per unit massa. Dalam masalah persamaan konservatif massa, persamaan dikembangkan untuk laju perubahan per unit volume. Laju perubahan sifat ϕ per unit volume untuk suatu partikel fluida dinyatakan sebagai hasil dari Dϕ/Dt dan densitas ρ, oleh karena itu

ρ

∂ φ Dφ = ρ +u.grad φ ∂t Dt

(2.40)

Persamaan konservatif massa mengandung massa per unit volume (contoh: densitas ρ) sebagai jumlah yang dikonservasi. Total laju perubahan densitas dan konfektif dalam persamaan konservatif massa untuk sebuha elemen fluida adalah ∂ ρ +div(ρu ) ∂t Dengan mengkondisikan terhadap berbagai sifat adalah ∂ (ρφ) +div(ρφu ) ∂t

(2.41)

32


Persamaan (2.41) menunjukkan laju perubahan dari ϕ per unit volume dijumlahkan dengan jumlah aliran ϕ keluar elemen fluida per unit volume. Sekarang kembali ditulis ulang untuk mengilustrasikan hubungan dengan turunan derifatif dari ϕ:

∂ ( ρφ ) Dφ  ∂φ   ∂ρ  + div ( ρφu ) = ρ  + u.grad φ  + φ  + div( ρu ) = ρ ∂t Dt  ∂t   ∂t 

(2.42)

 ∂ρ  + div( ρu ) = 0 disebabkan konservasi massa. Dalam kata lain. Kondisi φ   ∂t  Hubungan persamaan (2.42) dinyatakan sebagai Laju peningkatan ϕ dari elemen fluida + jumlah laju aliran ϕ keluar elemen fluida = laju peningkatan ϕ untuk sebuah partikel fluida Untuk mengkonstruksi komponen-komponen persamaan momentum dan energi berkaitan dengan ϕ dan laju perubahan per unit volume seperti yang didefiniskan oleh persamaan (2.42) dan (2.40) dapat dilihat pada table 2.2

Tabel 2.2 Persamaan momentum dan energi

 ∂ ( ρu )  + div( ρuu )  ∂t 

Momentum-x

y

ρ

Du Dt

φ

Momentum-y

v

ρ

Dv Dt

φ

Momentum-z

w

ρ

Dw Dt

φ

Energi

E

ρ

DE Dt

φ

 ∂ ( ρv)  + div( ρvu )  ∂t 

 ∂ ( ρw)  + div( ρwu )  ∂t 

 ∂ ( ρE )  + div( ρEu )  ∂t 

33


2.11

Persamaan Momentum Hukum Newton kedua mengatakan laju perubahan momentum sebuah

partikel fluida sama dengan total penjumlahan gaya pada partikel.

Laju peningkatan momentum partikel fluida = total gaya pada partikel fluida

Laju peningkatan momentum pada x-, y- dan z- per unit volume sebuah partikel fluida dinyatakan sebagai

ρ

Du Dt

ρ

Dv Dt

ρ

Dw Dt

(2.43)

Disini dibedakan dua jenis gaya yang bekerja pada partikel fluida: a. gaya permukaan yang meliputi : gaya tekanan, gaya kekentalan b. gaya bodi yang meliputi : gaya gravitasi, gaya sentrifugal, gaya Coriolis, gaya elektromagnetik tegangan pada elemen fluida dinyatakan dalam tekanan dan sembilan viscous stress (tegangan kekentalan) yang komponennya telihat pada Gambar 2.17. Tekanan, tegangan normal, dinyatakan dalam p. viscous stresses dinyatakan dalam τ. Umumnya notasi akhir τij digunakan untuk mengindikasikan arah dari viscous stress. Akhiran i dan j pada τij mengindikasikan bahwa tegangan bergerak kearah j pada permukaan normal i.

τ zz

τ zy τ yz τ yy

τ zx

τ yx

τ xy

τ xx

τ xy

τ xz

z

τ zz

y

τ xx

τ yx τ zx

τ xz

τ yy

τ yz

τ zy

x

Gambar 2.17 Komponen viscous stress[11]

34


Jika dipertimbangkan gaya pada komponen x akibat tekanan p dan tegangan komponen τxx , τyx , τzx yang dapat dilihat pada Gambar 2.18. Gaya yang sejajar dengan arah sebuah sumbu co-ordinat mempunyai tanda positif dan yang gaya berlawanan arah memperoleh tanda negatif. Total gaya pada arah x adalah penjumlahan dari gaya di komponen-komponen elemen fluida.

τ yx +

p−

τ xx −

∂τ yx 1 δy ∂y 2

τ zx +

∂τ zx 1 δz ∂z 2

τ yx −

∂p 1 δx ∂x 2

∂τ yx 1 δy ∂y 2

p+

∂τ xx 1 δx ∂x 2 z

τ xx +

∂p 1 δx ∂x 2

∂τ xx 1 δx ∂x 2

y x

τ zx −

∂τ zx 1 δz ∂z 2

Gambar 2.18 Tegangan pada komponen-komponen pada arah X [11]

Pada permukaan yang berpasangan (E,W) kita peroleh  ∂p 1   δx  − τ xx  p − ∂x 2     ∂p ∂τ xx = − + ∂x  ∂x

  ∂p 1   −  p + ∂x 2 δx   ∂τ 1     − xx δx  ∂y∂z +  ∂y∂z   ∂x 2  ∂τ xx 1  δx  + τ xx + ∂x 2   

 ∂x∂y∂z 

(2.44a)

Total gaya pada arah x pada permukaan yang berpasangan (N,S) adalah ∂τ ∂τ 1  ∂τ yx 1    δy δxδz + τ yx + yx δy δxδz = yx δxδyδz − τ yx − ∂y ∂y 2  ∂y 2   

(2.44b)

Dan total gaya pada arah x pada permukaan T dan B adalah

∂τ 1  ∂τ 1  ∂τ   − τ zx − zx δz δxδy + τ zx + zx δz δxδy = zx δxδyδz ∂z 2  ∂z 2  ∂z  

35

(2.44c)


Total gaya per unit volum pada fluida disebabkan tegangan-tegangan permukaan ini sama dengan penjumlahan persamaan (2.44a,b,c) dibagi oleh volume δxδyδz :

∂ (− p + τ xx ) ∂τ yx ∂τ zx + + ∂x ∂y ∂z

(2.45)

tanpa mempertimbangkan body force pada detail perhitungan effek secara keseluruhan dapat dimasukkan dengan menyatakan sebuah sumber (source) SMx dari momentum x per unit volume per unit waktu. Persamaan momentum pada komponen x diperoleh dengan mengatur laju perubahan pada momentum x pada partkel fluida yang jumlahnya sama dengan total gaya pada arah x di elemen akibat dari tegangan permukaan ditambah dengan laju peningkatan pada momentum akibat sumber (source):

ρ

Du ∂ (− p + τ xx ) ∂τ yx ∂τ zx = + + + S Mx ∂x ∂y Dt ∂z

(2.46a)

Maka persamaan momentum untuk komponen y dan z dapat dituliskan:

ρ ρ

Dv ∂τ xy ∂ (− p + τ yy ) ∂τ zy = + + + S My ∂x ∂y Dt ∂z

(2.46b)

Dw ∂τ xz ∂τ yz ∂ (− p + τ zz ) = + + + S Mz ∂x ∂y Dt ∂z

(2.46c)

Tanda pada tekanan menandakan tekanan bekerja berlawanan dengan norma viscous stress, dikarenakan pada umumnya tanda umum yang digunakan untuk beban tarik adalah positif tegangan normal sehingga pada tekanan yang mana bekerja sebagai tekanan beban normal maka memiliki tanda negatif. Efek pada tegangan permukaan dihitung secara terpisah; kondisi SMx , SMy dan SMz pada persamaan (2.46a,b,c) dimasukkan pada persamaan hanya diakibatkan oleh gaya bodi. Sebagai contoh , gaya bodi diakibatkan oleh gravitasi akan membuat nilai SMx =0, SMy =0 dan SMz = -ρg

36


2.12

Persamaan-Persamaan Energi Persamaan energi diturunkan dari hukum pertama termodinamika yang

mana menyebutkan bahwa laju perubahan energi sebuha partikel fluida sama dengan laju pertambahan panas ke partikel fluida ditambah dengan laju kerja yang dilakukan pada partikel.

Laju pertambahan energi pada partikel fluida = jumlah total panas yang ditambahkan ke partikel fluida + jumlah total kerja yang dilakukan pada partikel fluida

Seperti sebelumnya dimana untuk menurunkan sebuah persamaan dari laju pertambahan energi dari sebuah partikel fluida per unit volume dinyatakan sebagai

ρ

DE Dt

(2.47)

2.12.1 Kerja Yang Dilakukan Oleh Gaya Permukaan Laju kerja yang dilakukan pada partikel fluida di elemen disebabkan sebuah tegangan permukaan adalah sama dengan produk gaya dan kecepatan komponen dalam arah pada gaya yang bekerja. Sebagai contoh, gaya yang diberikan pada persamaan (2.44a,b,c) semuanya bergerak pada arah x. Maka kerja yang dilakukan oleh gaya-gaya tersebut dinyatakan sebagai

  ∂ (τ yx u ) 1    ∂ (τ xx u ) 1   ∂ ( pu ) 1    u τ δy   − − pu x u x δ τ δ − − − −      yx xx    y ∂ 2 x x ∂ ∂ 2 2          ∂y∂z +  ∂x∂z   ∂ (τ xx u ) 1  ∂ (τ yx u ) 1  ∂ ( pu ) 1     + τ u + δx  + τ xx u + δx    pu + δy  ∂x 2   ∂x 2     yx ∂y 2     ∂ (τ zx u ) 1   δz   − τ zx u − ∂z 2     + ∂x∂y   ∂ (τ zx u ) 1  δz  + τ zx u + ∂z 2    Jumlah kerja yang dilakukan oleh gaya permukaan ini yang bekerja pada arah dinyatakan sebagai

37


 ∂[u (− p + τ xx )] ∂ (uτ yx ) ∂ (uτ zx )  + +  δxδyδz ∂x ∂y ∂z  

(2.48a)

Tegangan permukaan pada komponen-komponen pada arah y dan z juga melakukan kerja pada partikel fluida. Pengulangan dari proses diatas memberikan tambahan laju kerja yang dilakukan pada partikel fluida yang diakibatkan kerja yang dilakukan oleh gaya permukaan tersebut :

[

]

 ∂ (vτ xy ) ∂ v(− p + τ yy ) ∂ (vτ zy )  + +  δxδyδz ∂y ∂z   ∂x

(2.48b)

Dan  ∂ ( wτ xz ) ∂ ( wτ yz ) ∂[w(− p + τ zz )] + +  δxδyδz ∂z ∂y  ∂x 

(2.48c)

Total laju kerja yang dilakukan per unit volume pada partikel fluida oleh semua gaya permukaan dnyatakan oleh penjumlahan persamaan (2.48a,b,c) yang dibagi oleh volume δxδyδz. Dalam kondisi terdapat tekanan di setiap arah maka tekanan dapat di satukan bersaman dan ditulis lebih ringkas dalam bentuk vektor :

−

∂ (up ) ∂ (vp) ∂ ( wp ) − − = −div( pu ) ∂x ∂y ∂z

Sehingga total laju kerja yang dilakukan pada partikel fluida disebabkan tegangan-tegangan permukaan:  ∂ (uτ xx ) ∂ (uτ yx ) ∂ (uτ zx ) ∂ (vτ xy ) ∂ (vτ yy ) ∂ (vτ zy )  + + + + +   ∂x ∂y ∂z ∂x ∂y ∂z   [−div( pu )] +   ∂ ( wτ xz ) ∂ ( wτ yz ) ∂ ( wτ ) zz + +  + ∂x ∂y ∂z  (2.49) 

38


2.12.2 Fluks Energi Akibat Panas Flux panas (heat flux) vekor q memiliki tiga komponen qx, qy dan qz

qy +

qz +

∂q y 1 δy ∂y 2

∂q z 1 δz ∂z 2

qx + qx −

∂q x 1 δx ∂x 2

qy − z

y x

qz −

∂q x 1 δx ∂x 2

∂q y 1 δy ∂y 2

∂qz 1 δz ∂z 2

Gambar 2.19 Komponen dari vektor heat flux[11]

Jumlah perpindahan panas ke partikel fluida disebabkan aliran panas pada arah x dinyatakan oleh perbedaan antara laju panas masuk melewati permukaan W dan laju panas yang hilang melewati permukaan E:  ∂q x 1   ∂q ∂q 1  δx  −  q x + x δx δyδz = − x δxδyδz  q x − ∂x 2   ∂x 2  ∂x 

(2.50a)

Maka jumlah laju perpindahan panas ke fluida akibat aliran panas dalam arah y dan z adalah :

−

∂q y ∂y

δxδyδz dan −

∂q z δxδyδz ∂z

(2.50b-c)

Total laju panas yang ditambahkan ke partikel fluida per unit volume akibat alirna panas melewati batas (boundaries) merupakan penjumlahan dari persamaan (2.50a,b,c) yang dibagi oleh volume δxδyδz.

−

∂q x ∂q y ∂q z − − = −div q ∂y ∂x ∂z

(2.51)

39


Hukum Fourier untuk panas konduksi behubungan dengan heat flux pada gradient temperatur lokal. Sehingga

q x = −k

∂T ∂x

q y = −k

∂T ∂y

q z = −k

∂T ∂z

Ini dapat ditulis dalam bentuk vektor sebagai berikut : Q = -k grad T

(2.52)

Menggabungkan persamaan (2.51) dan (2.52) memperoleh bentuk akhir dari laju panas yang ditambahkan ke partikel fluida akibat panas konduksi yang melewati batas elemen: -div q = div (k grad T)

(2.53)

2.12.3 Persamaan Energi Sumber (source) energi per unit volume per unit waktu pada momentum dinyatakan dalam SE. Dimana energi pada partikel fluida dipastikan dengan menyamakan laju perubahan energi partikel fluida dengan total laju kerja yang dilakukan pada partikel fluida dan total laju panas yang ditambahkan pada fluida dan laju pertambahan energi. Persamaan energi dinyatakan sebagai :

 ∂ (uτ xx ) ∂ (uτ yx ) ∂ (uτ zx ) ∂ (vτ xy ) ∂ (vτ yy )  + + + +   ∂x ∂y ∂z ∂x ∂y  DE  = −div( pu ) + ρ  (2.54)  ∂ (vτ zy ) ∂ ( wτ xz ) ∂ ( wτ yz ) ∂ ( wτ ) Dt zz + + +  + ∂z ∂x ∂y ∂z   + div(kgrad

T) + SE

Pada persamaan (2.54) nilai E = I + ½ (u2 + v2 + w2) Walaupun persamaan (2.54) merupakan persamaan energi yang sangat layak, tapi source dapat dirubah menjadi (mechanical) energi kinetik untuk memperoleh sebuah persamaan untuk energi dalam i atau temperaut T. bagian dari persamaan energi dapat dijadikan acuan terhadap energi kinetik dapat diperoleh dengan mengalikan persamaan x momentum dengan komponen kecepatan u, persamaan y momentum dengan v dan persamaan momentum z dengan w dan

40


menambahkan hasilnya bersamaan. Hal ini menghasilkan persamaan konservasi untuk energi kinetic:

1  D  (u 2 + v 2 + w 2 ) ∂τ yx ∂τ zx  ∂τ 2  = −u.grad p + u  xx + + ρ  ∂ ∂ ∂z Dt x y   ∂τ xy ∂τ yy ∂τ zy   + v + +  ∂ ∂ ∂ x y z   ∂τ yz ∂τ zz   ∂τ  + u.S M + w xz + +  ∂ ∂ x y ∂ z  

    (2.55)

Dengan melakukan subtraksi antara persamaan (2.54) dan (2.55) dan mendefiniskan

sebuah sumber baru sebagai Si = SE – u.SM menghasilkan

persamaan energi dalam

ρ

Di = − p.div u Dt + τ yy

∂u ∂u ∂u ∂v + τ yx + τ zx + τ xy ∂x ∂y ∂z ∂x ∂w ∂w ∂w + τ yz + τ zz + Si ∂x ∂y ∂z (2.56)

+ div(k .grad T ) + τ xx ∂v ∂v + τ zy + τ xx ∂y ∂z

Dalam kondisi spesial sebuah fluida inkompresibel memiliki nilai i = cT, dimana c adalah panas spesifik, dan div u =0. Ini mengakibatkan pengulangan kembali persamaan (2.56) menjadi persamaan temperatur

ρc

DT ∂u ∂u ∂u ∂v ∂v = div(k .grad T ) + τ xx + τ yx + τ zx + τ xy + τ yy Dt ∂x ∂y ∂z ∂x ∂y ∂v ∂w ∂w ∂w + τ zy + τ xx + τ yz + τ zz + Si ∂z ∂x ∂y ∂z

(2.57)

Untuk aliran kompresibel persamaan (2.54) sering diatur ulang untuk memberikan sebuah persamaan entalpi. Panas spesifik h dan total entalpi spesifik ho sebuah fluida didefiniskan sebagai : h = I + p/ρ dan ho = h + ½ (u2 + v2 + w2 ) mengkombinasikan dua definisi ini menjadi satu untuk energi spesifik E diperoleh ho = I + p/ρ + ½ (u2 + v2 + w2 ) = E + p/ρ

41

(2.58)


dengan mensubtitusikan persamaan (2.52) dan persamaan (2.48) dan beberapa pengaturan ulang dihasilkan persamaan (total) entalpi:

 ∂ (uτ xx ) ∂ (uτ yx ) ∂ (uτ zx ) ∂ (vτ xy ) ∂ (vτ yy )  + + + +   ∂ ∂ ∂ ∂ ∂y  x y z x ∂ ( ρho )  + div( ρ .ho .u ) =   ∂ (vτ zy ) ∂ ( wτ xz ) ∂ ( wτ yz ) ∂ ( wτ ) ∂t zz + + +  + ∂z ∂x ∂y ∂z   + div( kgrad

2.13

T) +

∂p + Sh ∂t

(2.59)

Konduksi Transien Untuk konduksi transien dibentuk oleh persamaan

∂T ∂  ∂T  = k +S ∂t ∂x  ∂x 

ρc

(2.60)

Sebagai tambahan biasanya ditambahkan variabel panas spesifik c (J/kg/K) dari material.

δxwP w

W

δxPe e

P

E

Gambar 2.20 Kontrol volume untuk kondisi satu dimensi[11]

Untuk kondisi satu dimensi seperti pada Gambar 2.20, integral dari persamaan (2.60) terhadap batas control volume dan interval waktu t hingga t+Δt adalah t + Δt

∫ t

∫ ρc

CV

δT dV .dt = δt

t + Δt

∫ t

∂T ∂ ∫CV ∂x =(k ∂x )dV .dt +

t + Δt

∫ ∫ sdV .dt t CV

Atau dapat ditulis e

∫ w

t + Δt t + Δt δT  c dt dV ρ =  ∫  ∫t δt   t

 ∂T    kA  ∂x  e

t + Δt  ∂T    kA   dt + ∫ SΔV .dt  ∂x  w  t

42

(2.61)


Pada persamaan (2.61), A adalah luas permukaan dari control volume, ΔV adalah volume, yang mana besarnya sama dengan AΔx dimana Δx adalah lebar dari control volume, dan S adalah kekuatan rata-rata source. Jika node pada temperatur diasumsikan merata disemua kontrol volume, pada sisi kiri dapat dituliskan

∫

CV

t + ∆t δT  o dt  dV = ρc(T p − T p )∆V  ∫ ρc t δ  t 

(2.62)

Pada persamaan (2.62) superskrip ‘o’ mengacu pada temperatur pada waktu t; temperatur pada waktu t + Δt tidak di subskrip. Jika kita menerapkan pusat perbedaan untuk kondisi difusi pada sisi kanan persamaan (2.61) dapat dituliskan t + ∆t

∫ t

 TE − TP  k e A δ XPE 

  T − TW  −  k w A P δ XWP  

t + ∆t   dt + ∫ S ∆V .dt = ρc(T p − T p o )∆V  t

(2.63)

Untuk mengevaluasi sisi kanan dari persamaan kita perlu membuat suatu asumsi mengenai variasi dari TP, TE dan TW dengan waktu. Kita dpat menggunakan temperatur pada waktu t atau pada waktu t + Δt untuk menghitung penigkatan waktu atau , mengkombinasikan temperatur pada waktu t dan t + Δt. Kita dapat menyamaratakan pendekatan dengan cara menggunakan parameter θ yang nilainya antara 0 dan 1 dan menulis integral IT dari temperatur TP terhadap waktu t + ∆t

It =

∫T

P

[

]

dt = θTP + (1 − θ )TP0 ∆t

(2.64)

t

Oleh karena itu dengan mensubtitusikan nilai θ terhadap persamaan (2.64), maka dapat ditulis beberapa persamaan baru yang dapat dilihat pada tabel 2.3 Tabel 2.3 Hasil dari pendekatan persamaan parameter θ θ

0

0,5

1

IT

TPo ∆t

1 / 2(TP + TPo )∆t

TP ∆t

43


Dapat menggaris bawahi nilai dari integral IT ; jika θ=0 temperatur pada waktu t digunakan; jika θ=1 temepratur pada waktu yang baru t+Δt digunakan; dan jika θ=1/2, temperatur pada t dan t+Δt bernilai sama. Menggunakan persamaan (2.64) untuk TW dan TE pada persamaan (2.63) dan membaginya dengan AΔt, diperoleh:

 TP − TPo  ∆t

ρc

  k (T − TP ) kW (TP − TW )   ∆x = θ  e E − +  δ XWP  δ XPE  

(2.65)

 k (T o − TP0 ) kW (TPo − TW0 )  − (1 − θ )  e E  + S ∆x δ XPE δ XWP  

Kita dapat mengatur ulangnya mejadi bentuk  ∆x  k k + θ  e + W  ρc  δ XPE δ XEP  ∆t

 k k TP = e θTE + (1 − θ )TE0 + W θT\W + (1 − θ )TW0 δ XPE δ XEP   ∆x k k  +  ρc − (1 − θ ) e − (1 − θ ) W TP0 + S ∆x δ XPE δ XEP   ∆t

[

]

[

]

(2.66)

Sekarang identifikasikan koeffisien TW dan TE sebagai aW dan aE dan menulis persamaan (2.66) menjadi bentuk yang lebih mudah:

[ + [a

]

[

a P TP = aW θTW + (1 − θ )TWo + a E θTE + (1 − θ )TEo o P

]

]

dimana a P = θ (aW + a E ) + a Po , b = S∆x , aE = a Po = ρc

(2.67)

− (1 − θ )aW − (1 − θ )a E T + b o P

ke

δ XPE

,aW =

kW

δ XWP

, dan

∆x ∆t

Bentuk akhir dari persamaan diskretisasi bergantung pada nilai θ. Dimana θ bernilai 0, kita hanya menggunakan temperatur TPo, Two dan TEo pada waktu t persamaan (2.67) untuk menentukan TP pada waktu yang waktu yang baru. Hasil dari skema ini disebut eksplisit. Dimana 0<θ≤1 waktu yang baru digunakan pada persamaan; hasil dari skema ini disebut implisit. Untuk kasus dimana θ=1 dinamakan implisit penuh dan kasus untuk θ=1/2 disebut skema Crank-Nicolson.

44


2.13.1 Skema Eksplisit Pada

skema

eksplisit

sumber

dilinirkan

sebagai

b=Su+SPTPo .

mensubtitusikan θ=0 kedalam persamaan (2.67) memberikan diskretisasi eksplisit dari persamaan perpindahan panas konduski transien:

[

]

a P TP = aW TWo + a E TEo + a Po − (aW + a E − S P ) TPo + S u

Dimana aP = aPo , aE =

ke

δ XPE

, aW =

kW

δ XWP

dan a Po = ρc

(2.68)

∆x ∆t

2.13.2 Skema Crank Nicolson Metode Crank-Nicolson diperoleh dari mengatur nilai θ=1/2 pada persamaan (2.67). Sekarang persamaan diskretisasi konduksi transien adalah :  TW + TW0   o a E aW  o  TE + TE0  − a P TP = a E  ) TP + b  + a P −  + aW  2 2 2   2    

Dimana a P = dan a Po = ρc

(2.69)

k k 1 1 1 (aW + a E ) + a Po − S P ,b = S u + S P T p0 , aE = e , aW = W δ XPE 2 2 2 δ XWP

∆x ∆t

2.13.3 Skema Implisit Penuh Ketika nilai θ diatur sama dengan 1 , diperoleh skema implisit penuh. Persamaan menjadi

a P TP = aW TW + a E TE + a Po TPo + S u

(2.70)

Dimana a P = a Po + aW + a E − S P , a Po = ρc

k k ∆x , aW= W dan aE = e ∆t δ XWP δ XPE

Metode implisit penuh direkomendasikan untuk perhitungan CFD dengan tujuan yang umum dikarenakan hasil perhitungan yang stabil.

45


Untuk masalah 2 dimensi (2D) persamaan (2.60) dapat dikembangkan menjadi sebuah persamaan baru menjadi

ρc

∂φ ∂  ∂φ  ∂  ∂φ  S = k  + k ∂t ∂x  ∂x  ∂y  ∂y 

(2.71)

Jika persamaan (2.71) dideskritisasikan terhadap kontrol volume, persamaan yang dihasilkan adalah :

a P φ P = aW φW + a E φ E + a S φ S + a N φ N + a Po TPo + S u

(2.72)

∆V o Dimana a P = aW + a E + a S + a N + a p − S P dan a Po = ρc ∆t N

δxPe e

δysP

P

δysn

δy Pn w

W

n

δxwP

E

s

δxwe S

Gambar 2.21 Kontrol volume untuk situasi dua dimensi[11]

Dimana koeffisien aW, aE untuk masalah 1D, dan aW, aE, aS, aN untuk 2D seperti dapat dilihat pada Gambar 2.20, serta aW, aE, aS, aN, aB, aT untuk 3D; b=(Su+Spϕp) adalah sumber. Hasil dari koeffisien-koeffisien yang diberikan adalah :

46


Tabel 2.4 implisit untuk 1D, 2D, dan 3D aW

aE

aS

aN

aB

aT

1D

ΓW AW δxWP

Γe Ae δx PE

-

-

-

-

2D

ΓW AW δxWP

Γe Ae δx PE

Γs As δy SP

Γn An δy PN

-

-

3D

ΓW AW δxWP

Γe Ae δx PE

Γs As δy SP

Γn An δy PN

Γb Ab δz BP

Γt At δz PT

Nilai untuk volume dan luas permukaan cell untuk ketiga kasus:

Tabel 2.5 Nilai volume dan luas permukaan cell 1D

2D

3D

ΔV

Δx

ΔxΔy

ΔxΔyΔz

Aw = Ae

1

Δy

ΔyΔz

An = As

-

Δx

ΔxΔz

Ab = At

-

-

ΔxΔy

2.14

Konveksi Transien Untuk pendekatan turunan implisit penuh untuk masalah multi dimensi

maka dibutuhkan tambahan apo lada koeffisien pusat ap dan kontribusi apoϕpo sebagai sumber tambahan pada sisi kanan persamaan. Sedangkan koeffisienkoefisin lainnya sama dengan persamaan turunan untuk masalah steady state. Persamaan perpindahan sifat ϕ transien dinyatakan sebagai : ∂ ( ρφ ) + div( ρφu ) = div(Γ grad φ ) + Sφ ∂t

(2.73)

Atau dapat difinisikan sebagai

47


Laju peningkatan ϕ pada elemen fluida + jumlah laju aliran ϕ keluar elemen fluida = laju peningkatan ϕ disebabkan difusi + laju peningkatan ϕ disebabkan sumber Persamaan (2.73) disebut persamaan angkut (transport) untuk sifat ϕ. Untuk tiga dimensi konveksi-difusi untuk sifat ϕ dalam suatu kecepatan bidang u dinyatakan sebagai:

∂ (ρφ) ∂ (ρuφ) ∂ (ρvφ) ∂ (ρwφ) ∂  ∂ φ + + + = Γ  ∂t ∂x ∂y ∂z ∂x  ∂x  +

∂  ∂ φ ∂  ∂ φ Γ  + Γ  + S ∂y  ∂y  ∂z  ∂z 

(2.74)

Persamaan untuk turunan berkembang penuhnya adalah

aP φ P = aW φW +aE φ E +aS φ S +a N φ N +aB φ B +aT φT +aPo TPo +Su

Dimana

aP = aW +aE +aS +a N +aB +aT +a op +ΔF S P

dengan aPo =

dan

(2.75)

ρ op ΔV Δt

SΔV = Su + S p φ P

Dimana koeffisien-koeffisien untuk persamaan (2.75) dapat dilihat pada tabel 2.6dan 2.7 untuk kasus 1D,2D. dan 3D

48


Tabel 2.6 nilai-nilai koeffisien untuk persamaan (2.75) Aliran 1D aW

aE

Max Fw , Dw +

Aliran 2D Fw ,0 2

Max  Fe  − Fe ,  De − 2  

aS

aN

aB

-

-

-

Max Fw , Dw +

Aliran 3D Fw ,0 2

Max   ,0  

 Fe  − Fe ,  De − 2  

Max Fw , Dw +

Fw ,0 2

Max   ,0  

 Fe  − Fe ,  De − 2  

  ,0  

Max

Max

  Fs   ,0  Fs ,  D s + 2    

  Fs   ,0  Fs ,  D s + 2    

Max

Max

 Fn    ,0 − Fn ,  Dn + 2    

 Fn    ,0 − Fn ,  Dn + 2    

-

Max   Fb  Fb ,  Db + 2  

  ,0  

aT

-

-

 Ft    ,0 − Ft ,  Dt + 2    

ΔF

Fe – Fw

F-Fw+Fn-Fs

Fe-Fw+Fn-Fs+Ft-Fb

Untuk nilai-nilai koeffisien tersebut nilai F dan D dihitung berdasarkan formula berikut :

49


Tabel 2.7 nilai F dan D untuk tabel 2.10 Face

w

e

s

n

b

t

F

( ρu) w Aw

( ρu) e Ae

( ρu) s As

( ρu) n An

( ρu) b Ab

( ρu) t At

D

Γw Aw δx WP

Γe Ae δx PE

Γs As δy SP

Γn An δy PN

Γb Ab δz BP

Γt At δz PT

2.15

Kelebihan dan Kekurangan Relaksasi Dalam solusi pendekatan dari persamaan-persamaan aljabar atau skema

yang berulang untuk menganalisa masalah non linear, sering diinginkan untuk mempercepat atau memperlambat perubahan yang terjadi, dari suatu iterasi ke iterasi berikutnya. Proses ini disebut kelebihan atau kekurangan relaksasi (overrelaxation atau underrelaxation), bergantung pada variabel berubah untuk dipercepat atau diperlambat. Overrelaxation sering digunakan yang berhubungan dengan metode Gauss-Seidel yang hasil skemanya dikenal sebagai Successive Over-Relaxation (SOR). Underrelaxation merupakan hal yang berguna untuk masalah nonlinear. Hal ini sering digunakan untuk mencegah terjadinya perbedaan (divergen) dalam solusi pendekatan. Ada banyak cara memperkenalkan overrelaxation atau underrelaxation. Untuk masalah konduksi dengan persamaan diskretisasi umum dengan metode titik ke titik memiliki bentuk :

a P TP = ∑ a nbTnb + b

(2.76)

Dimana nilai Tnb memiliki nilai yang bersebelahan dengan titik yang sekarang Tp yang nilainya sudah diketahui dari iterasi sebelumnya. Dalam tiap kasus nilai Tnb merupakan nilai terakhir dari temperatur yang berada disebelah titik sekarang (Tp). Persamaan (2.76) juga dapat ditulis

TP =

∑a

T +b

nb nb

(2.77)

aP

50


Jika kita masukkan nilai Tp* yang nilainya diperoleh dari nilai Tp dari iterasi sebelumnya pada sisi kanan persamaan, diperoleh  ∑ a nbTnb + b  − TP *  TP = TP * +    aP  

(2.78)

Bentuk ini dapat dimodifikasi dengan menambahkan sebauh faktor relaksasi α, sehingga diperoleh :  ∑ a nbTnb + b  − TP *  TP = TP * +α    aP  

(2.79a)

Atau

aP

α

TP = ∑ a nbTnb + b + (1 − α )

aP

α

TP *

(2.79b)

Harus dicatat bahawa, ketika iterasi konvergen, nilai Tp menjadi sama dengan Tp*. Persamaan (2.79) menunjukkan hasil nilai T digunakan pada persamaan (2.76). nilai relaksasi factor antara 0 dan 1 yang berdampak pada underrelaxation sehingga nilai Tp semakin dekat dengan Tp*. Untuk nilai α yang kecil . perubahan Tp menjadi lebih lambat. Ketika α lebih besar dari 1 maka akan menghasilkan overrelaxation. Tidak ada peraturan khusus dalam memilih nilai α yang terbaik. Nilai yang paling optimum bergantung pada beberapa faktor, seperti kondisi masalah, banyaknya titik grid, jarak antar grid,dsb. Nilai α yang sesuai dapat ditentukan berdasarkan pengalaman dan ekplorasi perhitungan dari masalah yang diberikan.

2.16 Perhitungan Metode Elemen Hingga Pada kenyatannya sangat banyak permasalahan di bidang keteknikan yang tidak dapat diselesaikan dengan persamaan eksak. Untuk memperoleh solusinya dapat digunakan turunan-turunan dari beberapa persamaan yang kompleks ditambah lagi dengan kondisi awal benda dan kondisi batas yang digunakan. Untuk mengatasinya maka dilakukanlah pendekatan numerik. Metode elemen hingga lebih menggunakan persamaan yang berkaitan (integral formulations)

51


dibandingkan dengan menggunakan turunan dari suatu persamaan untuk membuat suatu sistem persamaan aljabar. Ada

beberapa

tahap

dalam

penyelesaian

elemen

hingga

yaitu,

preprocessing ; merupakan tahap pembuatan domain dari elemen hingga (node dan elemen) serta menentukan kondisi batas, beban, dan kondisi awal seperti yang dapat dilihat pada Gambar 2.22. Solution; merupakan tahap penyelesaian persamaan linear dan nonlinear untuk memperoleh hasil pada node seperti nilai perpindahan

atau

temperatur

pada

premasalahan

perpindahan

panas.

Postprocessing; merupakan tahapan untuk memperoleh informasi dari tahap solution, dalam tahap ini kita mendapatkan nilai-nilai tegangan, regangan, heat flux, dsb.

Elemen 1

1

Elemen 2

2 Elemen 1

Elemen 3

3 Elemen 2

Elemen n

n+1

n Elemen 3

Elemen n

Gambar 2.22 Menentukan node dan elemen termasuk dalam tahapan preprocessing[7]

Untuk menetukan regangan dan tegangan akibat beban termal dengan menggunakan metode elemen hingga kita, dapat dijabarkan dengan menggunakan persamaan dasar regangan dan tegangan termal pada persamaan (2.8) dan (2.9), maka dengan menggunakan persamaan tegangan rata-rata dan mensubtitusikannya dengan persamaan (2.8) dan (2.9) diperoleh

F A F = α∆TEA δEA F= = kδ L

σ=

(2.80)

52


Untuk benda dengan luas penampang berbeda pada arah-y. Dilakukan pendekatan dengan menggunakan prinsip pegas, disebabkan persamaan (2.80) mirip dengan persamaan pegas linear F=k.x f = k eq (Ti +1 − Ti ) = k eq

Aavg E l

(u i +1 − u i ) =

( Ai +1 + Ai ) E (u i +1 − u i ) 2l

( A + Ai ) E = i +1 2l

(2.81)

Ai dan Ai+1, adalah luas penampang permukaan pada node 1 dan i+1, dan l merupakan panjang elemen. Dengan memodelkan elemen, maka kita dapat membuat diagram benda bebas pada node seperti yang terlihat pada Gambar 2.23

R1

Node 1

K1(u2-u1) K1(u2-u1) Node 2 K2(u3-u2)-T K2(u3-u2)-T

Node 3 K3(u4-u3) K3(un-u3) Node n Kn(un+1-un) Kn(un+1-un) Node n+1 Rn+1

Gambar 2.23 Diagram benda bebas pada node[7]

Berdasarkan rumus kesetimbangan dimana total gaya yang bekerja pada tiap node adalah nol. Pernyataan tersebut dapat menghasilkan persamaan berikut Node 1: R1 − k1 (T2 − T1 ) = 0 Node 2: k1 (T2 − T1 ) − k 2 (T3 − T2 ) 2 = 0 Node 3: k 2 (T3 − T2 ) − k 3 (Tn − T3 ) − T = 0

53


Node n: k 3 (Tn − T3 ) − k n (Tn +1 − Tn ) = 0 Node n+1 : k n (Tn +1 − Tn ) − R n +1 = 0

(2.82)

Nilai beban eksternal T pada beberapa kasus juga tidak selamanya terdapat pada node-node tertentu bergantung pada kondisi permasalahan, mis: tidak ada temperatur

di

titik

3,

maka

persamaan

menjadi

Node

3:

k 2 (T3 − T2 ) − k 3 (Tn − T3 ) = 0 Dengan mengatur ulang persamaan (2.82) dengan memisahkan antara nilai k, T dan gaya luar P , maka diperoleh :

k 1u1 − k 1u1

− k 1u 2 k 1u 2

k 2u2 k 2u2

− k 2u3 k 2u3

− k 3u n k 3u n

k 3u 3 − k 3u3

k 3u n − k 3u n

=

R1

= =

0 T

(2.83)

− k n u n +1 = 0 k n u n +1 = Rn +1

Jika persamaan kesetimbangan yang diberikan oleh persamaan (2.83) dalam bentuk matriks, menjadi :  k1 − k  1  0   0  0

− k1 k1 + k 2 − k2 0 0

0 − k2 k2 + k3 − k3 0

0 0 − k3 k3 + kn − kn

  u1  − R1   u   0     2    u3  =  T   − k n  u n   0      k n  u n +1   Rn +1  0 0 0

(2.84)

Jika dipisahkan antara matriks pembebanan dengan matriks perlawanan (resistance), makabentuk matriks menjadi − R1   k1  0  − k   1   0 = 0  0   0     Rn +1   0

{R}

=

− k1 k1 + k 2 − k2 0 0

0 − k2 k2 + k3 − k3 0

0 0 − k3 k3 + kn − kn

[k]

  u1   0   u  0    2      u 3  − T   − k n  u n  0      k n  u n +1   0  0 0 0

(2.85)

{u} – {T}

54


Jika dilihat pada Gambar 2.22 , maka diketahui benda diikat ada pada node 1 dan Pn+1 maka, maka baris pertama dan kelima pada persamaan (2.85) haruslah u1=0 dan un+1=0. Maka dengan menggunakan kondisi batas tersebut didapatkan bentuk persamaan matriks  1 − k  1  0   0  0

0 k1 + k 2 − k2

0 − k2 k2 + k3

0 0

− k3 0

  u1   0   u  0    2      u 3  = T   − k n  u n  0      1  u n +1   0 

0 0

0 0 0

− k3 k3 + kn 0

(2.86)

Dikarenakan dalam suatu elemen terdapat dua node, dan setiap nodenya berhubungan dengan perubahan temperatur, maka perlu membuat dua persamaan untuk setiap elemennya seperti yang terlihat pada Gambar 2.24

fi=keq(ui -ui+1)

fi=keq(ui+1-ui)

Node i

Node i atau

Node i+1

Node i+1

fi+1=keq(ui+1-ui)

fi+1=keq(ui+1-ui)

Gambar 2.24 Gaya-gaya pada sebuah elemen[7]

Dengan prinsip kesetimbangan dimana fi + fi+1 haruslah nol. Maka dapat dituliskan gaya-gaya yang bekerja pada node I dan i+1 adalah f i = k eq (u i − u i +1 )

(2.87)

f i +1 = k eq (u i +1 − u i )

Persamaan (2.87) dapat dinyatakan dalam bentuk matriks  f i   k eq  =  f i +1  − k eq

− k eq   u i    k eq  u i +1 

(2.88)

55


Dengan menggunakan persamaan (2.88) maka untuk tiap elemen dan kemudian disatukan kedalam bentuk formasi matriks global. Matriks untuk elemen (1) adalah :

[K ](1) = 

k1 − k1

− k1  k1 

dan bentuknya dalam matriks global adalah

[K ](1)

 k1 − k  1 = 0   0  0

− k1 k1

0 0 0  u1 0 0 0 u 2 0 0 0 u 3  0 0 0 u n 0 0 0 u n +1

0 0 0

Perubahan temperatur ditunjukkan pada matriks global untuk membantu dalam melihat pengaruh node terhadap elemen disebelahnya. Sama dengan elemen (1) maka untuk elemen selanjutnya

[K ]( 2) = 

k2 − k 2

− k2  k 2 

Posisinya dalam matriks global

[K ]( 2)

0 0 0 k 2  = 0 − k 2  0 0 0 0

[K ](3) = 

k3 − k 3

0 − k2 k2

0 0 0 0 0

0 0

0  u1 0 u 2 0 u 3  0 u n 0 u n +1

− k3  k 3 

Dan posisinya dalam bentuk matriks global

[K ](3)

0 0  = 0  0 0

0 0 0 0 0 k3 0 − k3 0 0

0 0 − k3 k3 0

0  u1 0 u 2 0 u 3  0 u n 0 u n +1

Dan

56


[K ]( 4) = 

k4 − k 4

− k4  k 4 

Bentuknya dalam matriks global

[K ]( 4)

0 0  = 0  0 0

0 0 0 0

0 0

0 0 0 0

0

k4 0 0 − k4

0  u1 0  u 2 0  u3  − k4  un k 4  u n +1

Bentuk akhir dalam matriks global lokal jika digabungkan, atau ditambahkan, digabungkan untuk setiap elemennya dalam bentuk matriks global

[K ](G ) = [K ](1) + [K ]( 2) + [K ](3) + [K ]( 4) [K ]

(G )

 k1 − k  1 = 0   0  0

− k1 k1 + k 2 − k2 0 0

0

0

− k2 k2 + k3 − k3 0

0 − k3 k3 + k4 − k4

0  0  0   − k4  k 4 

(2.89)

Dapat dilihat bahwa persamaan matriks global yang diperoleh menggunakan deskripsi tiap elemen pada persamaan (2.89) identik dengan persamaan global yang ditentukan dengan diagram benda bebas pada persamaan (2.85). Jika didapat beban di node 3 , maka bentuk matriks globalnya menjadi

[K ]

(G )

 1 − k  1 = 0   0  0

0 k1 + k 2 − k2 0 0

0 − k2 k2 + k3 − k3

0 0 − k3 k3 + k4

0

0

0   u1   0  0   u1   0      0   u1  = T   − k 4  u n  0      1  u n +1   o 

(2.90)

Sehingga dapat disimpulkan, formula metode elemen hingga selalu mengacu kepada bentuk persamaan umum : [matriks global K]{matriks perpindahan} = {matriks beban}

57


BAB II TINJAUAN PUSTAKA

Recommend Documents