Pemodelan Matematika dan Metode Numerik

Bab 3 Pemodelan Matematika dan Metode Numerik

3.1 Model Keadaan Tunak Model keadaan tunak hanya tergantung pada jarak saja. Oleh karena itu, distribusi temperatur gas sepanjang pipa sebagai fungsi dari jarak dapat dihitung dengan menggunakan hukum kekekalan energi yang menyatakan laju perubahan energi suatu sistem sama dengan jumlah panas dikurangi jumlah kerja pada sistem tersebut, dengan mengasumsikan proses perpindahan panas yang terjadi bersifat tunak dan tidak ada kerja yang dilakukan oleh sistem, sehingga diperoleh persamaan sebagai berikut : q=

dE . dt

(3.1)

Perhitungan laju perubahan energi dalam kasus ini dijelaskan dengan menggunakan konsep fluks. Akan diperhatikan proses aliran energi dalam suatu segmen x sampai dengan x + ∆x pada kontrol volum (Gambar 3.1), untuk mendapatkan laju 29

BAB 3. PEMODELAN MATEMATIKA DAN METODE NUMERIK

30

perubahan energi dengan konsep fluks. Misalkan energi yang masuk melalui titik x adalah f (x) sedangkan energi yang masuk melalui titik x + ∆x adalah f (x + ∆x). Jika f bernilai positif, maka energi mengalir masuk ke dalam segmen melalui sebelah kiri titik ujung x, sedangkan penulisan tanda minus untuk f (x + ∆x) dibutuhkan karena f (x + ∆x) > 0 menunjukkan energi mengalir ke sebelah kanan titik ujung x + ∆x. Sehingga laju perubahan energi adalah : f (x) − f (x + ∆x) .

(3.2)

Berdasarkan aproksimasi Taylor, Persamaan (3.2) merupakan turunan pertama f (x), yaitu : −

∂f ∆x . ∂x

(3.3)

Dengan mengabaikan energi potensial dan kinetik, dan mengasumsikan tidak ada efek nuklir, listrik, dan magnetik, maka dengan menguraikan energi yang terjadi pada sistem yaitu hanya energi panas (MCv T ) dan energi yang menyebabkan kehilangan tekanan (M ρp ) , dengan M sebagai mass flow, yaitu aliran massa yang melewati luas penampang A, M = ρvA, maka f adalah : p M(Cv T + ) . ρ

Gambar 3.1: Segmen x sampai dengan x + ∆x pada Kontrol Volum.

(3.4)


31

Dengan mensubstitusikan f yang berbentuk Persamaan (3.4) ke dalam Persamaan (3.3) dan menotasikan V sebagai volume sehingga V = 1/ρ, maka laju perubahan energi menjadi, −M

d(Cv T + pV) dx . dx

(3.5)

Dengan mengasumsikan gas yang dialirkan bersifat ideal, sehingga pV = RT dan C p = Cv + R, maka Persamaan (3.5) menjadi, − MC p dT .

(3.6)

Dengan mensubstitusikan Persamaan (3.6) ke dalam Persamaan (3.1), maka diperoleh, q = −MC p dT .

(3.7)

Mengenai panas (q) yang terjadi di sistem, dengan sistem pada kasus ini adalah sebuah pipa (Gambar 3.2), akan dibahas sebagai berikut. Terjadi aliran gas sepanjang pipa dari ujung pipa di kiri x sampai dengan ujung pipa di kanan x+∆x, dengan temperatur gas (T ) yang lebih besar dibandingkan dengan temperatur lingkungan (T amb ). Hal ini yang mengakibatkan terjadinya perpindahan panas sepanjang dx dari gas ke lingkungan sekitarnya secara konveksi, sehingga temperatur gas terus berkurang sampai mendekati temperatur sekitarnya. Konveksi adalah perpindahan panas yang terjadi antara permukaan dan media bergerak (fluida) yang mempunyai temperatur berbeda melalui proses difusi ataupun dengan cara mengalirnya fluida tersebut dari molekul dengan temperatur yang lebih tinggi ke molekul dengan temperatur yang lebih rendah. Dengan demikian, panas yang terjadi akibat proses konveksi menurut Newton’s law of cooling, dapat dituliskan ke dalam bentuk persamaan seperti, q = kL (T − T amb )dx ,

(3.8)


32

dengan kL adalah koefisien perpindahan panas secara konveksi yang bergantung pada kondisi batas yang dipengaruhi oleh geometri permukaannya dan gerakan fluida.

Gambar 3.2: Perpindahan Panas secara Konveksi di dalam Pipa. Dengan mensubstitusikan Persamaan (3.8) ke dalam Persamaan (3.7), maka diperoleh persamaan, MC p dT = −kL (T − T amb )dx .

(3.9)

Dengan membentuk Persamaan (3.9) seperti berikut, dT −kL = dx . (T − T amb ) MC p

akan diperoleh persamaan untuk menghitung distribusi temperatur gas sepanjang pipa sebagai fungsi dari jarak, yaitu dengan cara mengintegralkan kedua suku per-


33

samaan tersebut. Langkah - langkahnya adalah sebagai berikut : T (x) Z

T (0)

dT −kL = (T − T amb ) MC p

Zx dx . 0

Hasil pengintegralannya adalah sebagai berikut, ln(T − T amb )|TT (x) (0) =

−kL x x| . MC p 0

Dengan mensubstitusi batas integralnya akan menghasilkan persamaan, ! (T (x) − T amb ) −kL ln = x. (T (0) − T amb ) MC p

Dengan memberikan eksponen di kedua ruas persamaan di atas, maka akan diperoleh persamaan untuk menghitung distribusi temperatur gas sepanjang pipa sebagai fungsi dari jarak, yaitu: T (x) = T amb + (T (0) − T amb ) exp−αx ,

(3.10)

dengan α = kL /MC p Sedangkan untuk mendapatkan persamaan yang dapat menghitung distribusi tekanan sepanjang pipa sebagai fungsi dari jarak, dapat diperoleh dengan cara berikut : • Dari persamaan untuk mendapatkan mass flow, bisa diperoleh v = M/ρA, persamaan ini akan di substitusikan ke dalam persamaan penurunan tekanan aki-

BAB 3. PEMODELAN MATEMATIKA DAN METODE NUMERIK bat gaya gesek, Persamaan (2.37) ( ddxp =

−2 fg0 ρv2 D ),

34

sehingga akan diperoleh,

0 2 d p −2 fg M = . dx ρA2 D

• Persamaan mencari rapat massa, Persamaan (2.15) (ρg =

(3.11)

pMg ZRT ),

akan disubs-

titusikan ke dalam Persamaan (3.11), sehingga diperoleh, 0 2 d p −2 fg M ZRT = . dx pMg A2 D

• Dengan menggunakan turunan parsial,

d p2 dx

(3.12)

yang dapat dituliskan menjadi

bentuk 2p ddxp , sehingga dengan mensubstitusi Persamaan (3.12) ke dalam bentuk turunan parsial tersebut, akan diperoleh, 0 2 d p2 −4 fg M ZRT = . dx Mg A2 D

(3.13)

• Dengan mensubstitusi persamaan untuk menghitung distribusi temperatur gas sepanjang pipa sebagai fungsi dari jarak, Persamaan(3.10) ke dalam Persamaan (3.13), maka akan diperoleh, d p2 =

−4 fg0 M 2 ZR(T amb dx + (T (0) − T amb ) exp−αx dx) Mg A2 D

.

(3.14)

• Dengan mengintegralkan Persamaan (3.14), suku kiri terhadap p2 dengan batas p dari p(0) sampai dengan p(x) dan suku kiri terhadap x dengan batas x dari 0 sampai dengan x, maka akan diperoleh, p(x) p2 | p(0)

" i ! # −ZR 4 f 0 T (0) − T amb h −αx x x = 2 ) (exp ) M 2 . T amb x|0 − ( 0 α A Mg D

(3.15)


35

• Dengan memasukkan nilai batas integralnya, maka akan diperoleh persamaan akhir untuk menghitung distribusi tekanan yang hanya bergantung pada jarak saja, yaitu :

q p(x) =

p(0)2 − K M 2 ,

(3.16)

dengan " !# ZR 4 f 0 T (0) − T amb −αx K= 2 T amb x + ( )(1 − exp ) . α A Mg D

3.2

Model Keadaan Transien Berbeda dengan keadaan tunak, keadaan transien merupakan keadaan yang

bergantung pada jarak dan waktu. Pada Persamaan energi (2.48), akan dijabarkan panas per unit massa per unit satuan luas (qρA) yang terjadi pada sistem dengan kondisi transien. Dengan mengasumsikan perpindahan panas yang terjadi hanya proses konduksi yang melewati dinding pipa dan konveksi yang terjadi antara partikel fluida di dalam pipa, yang dapat mengakibatkan perpindahan panas dari gas ke lingkungan sekitar (Gambar 3.3). Dengan menggunakan konsep fluks pada proses konduksi dan asumsi seperti yang telah disebutkan di atas, maka didapatkan persamaan, qρAdx = qkonduksi | x − qkonduksi | x+dx − qkonveksi .

(3.17)

Konduksi adalah perpindahan panas melalui media diam yang diakibatkan oleh aktivitas partikel dan energi yang berpindah dari partikel dengan temperatur yang lebih tinggi ke partikel dengan temperatur yang lebih rendah. Dengan demikian, panas yang terjadi akibat proses konduksi menurut Fourier’s law of cooling, dapat


36

Gambar 3.3: Perpindahan Panas Konveksi dan Konduksi di dalam Pipa. dituliskan ke dalam bentuk persamaan, q = −λ

∂T , ∂x

(3.18)

dengan λ adalah konduktivitas bahan yang dilalui panas. Lalu, Persamaan (3.18) dan (3.8) akan disubstitusikan ke dalam Persaman (3.17), sehingga menjadi, ! ∂T ∂ ∂T ∂T + λ + (λ ) − kL (T − T amb )dx , qρAdx = −λ ∂x ∂x ∂x ∂x

dengan

∂ ∂T ∂x (λ ∂x )

sebagai perubahan panas akibat konduksi sepanjang dx. Dengan

menyederhanakan persamaan tersebut, akan diperoleh persamaan akhir yaitu : qρAdx =

∂ ∂T (λ ) − kL (T − T amb )dx . ∂x ∂x

(3.19)

Dengan mengkombinasikan antara Persamaan (2.48) dan (3.19), akan diperoleh


37

persamaan, ! ∂ ∂ ρvAp ∂ ∂T (ρvACv T dx) + ( dx ) − λA dx ∂x ρ ∂x ∂x |∂x {z } | {z } | {z } 1

2

(3.20)

3

∂ + kL (T − T amb ) dx + (ρACv T dx) = 0 . | {z } |∂t {z } 4

5

Akan diintegralkan tiap suku Persamaan (3.20) terhadap x, dengan x dari x = 0 sampai x = L. Lalu dengan memasukkan data yang dibutuhkan pada hasil pengintegralan, akan dilakukan analisis dimensi untuk mendapatkan model yang lebih sederhana. Data masukan yang dibutuhkan adalah,

Besaran γg P0 T0 TL T amb R Cv Cp kL L D λ Q

Keterangan Specific Gravity gas Tekanan di inlet Temperatur di inlet Temperatur di outlet Temperatur lingkungan Konstanta gas universal Specific Heat Specific Heat Koef. perpindahan panas konveksi Panjang pipa Diameter pipa Koef. kekasaran pipa Konduktivitas bahan Laju alir

Nilai 0, 6538 1146, 17 306, 48 285, 7 284, 7 518, 8 1, 759x103 2, 278x103 25 369000 0, 67945 0, 00001968 3, 4x10−2 8508791, 67

Satuan − psia 0K 0K 0K J/kg 0 K J/kg0 K J/kg0 K W/m0 K m m − W/m0 K m3 /h

Tabel 3.1: Data Masukan.

Dari data masukan di atas, dapat dicari rapat massa, faktor deviasi, kecepatan suara, dan faktor gesekan dengan korelasi pada bab 2, yaitu :


Besaran ρ Z c fg

Keterangan Rapat massa Faktor deviasi Kecepatan suara Faktor gesekan

Nilai 69, 51853 0, 8445 337, 188 0, 008

38

Satuan kg/m3 − m/s −

Tabel 3.2: Hasil Perhitungan ρ, Z, c dan fg .

Proses pengintegralan Persamaan (3.20) dan proses memberikan data yang ada di Tabel 3.1 dan 3.2 pada hasil pengintegralan adalah sebagai berikut : 1. ZL 0

∂ (ρvACv T ) dx ≈ ρQCv (T | x=L − T | x=0 ) ∂x 8508791, 67 × 1, 759x103 × −20, 78 3600 9 ≈ −6x10 . ≈ 69, 5 ×

2. ZL 0

! ∂ ρvAp dx ≈ ρQZR(T | x=L − T | x=0 ) ∂x ρ 8508791, 67 × 0, 8445 × 518, 8 × −20, 78 3600 ≈ −1, 5x109 .

≈ 69, 5 ×


39

3. ZL 0

! ∂ ∂T ∆T | x=L − ∆T | x=0 −λA dx ≈ −λA( ) ∂x ∂x ∆x π ≈ −3, 4x10−2 × 0, 679452 × (0, 56 − 23, 93) 4 −5 ≈ 4, 36x10 ,

Dengan ∆T | x=0 dan ∆T | x=L diperoleh dari selisih nilai temperatur pada keadaan tunak di titik x = 0 dan x = L. 4. ZL kL (T − T amb ) dx ≈ kL L (T − T amb ) 0

≈ 25 × 369000 × (306, 48 − 284, 7) ≈ 2x108 .

5. ZL 0

∂ ∆T (ρACv T ) dx ≈ ρACv L ∂t ∆t π 5 ≈ 69, 5 × 0, 679452 × 1759 × 369000 × 4 3600 7 ≈ 2, 2x10 ,

Dengan memisalkan penambahan temperatur yang terjadi adalah 5 0 K dalam waktu 1 jam. Dengan melihat hasil pengintegralan di atas, akan dilakukan analisis dimensi, yaitu dengan mengabaikan suku yang nilainya sangat kecil dibandingkan


40

dengan nilai yang lainnya, dalam arti nilai tersebut sangat kecil pengaruhnya. Dengan demikian, suku ke empat akan diabaikan, karena nilainya terlalu kecil dibandingkan dengan nilai suku lainnya, sehingga akan didapatkan penyederhanaan dari Persamaan energi (3.20) menjadi, " !# ∂ ∂ p ρAdx (Cv T ) + ρvAdx Cv T + = −kL (T − T amb ) dx . ∂t ∂x ρ

(3.21)

Dengan mensubstitusikan p = c2 ρ dan m(x, t) = ρ(x, t)v(x, t) ke dalam Persamaan (3.21), maka model aliran transien dengan pipa horizontal dapat direpresentasikan oleh persamaaan berikut :    ∂ρ ∂(m)   =0,  ∂t + ∂x   2     ∂ mρ +c2 ρ − fg m|m| ∂m  + = 2Dρ ,   ∂t ∂x     h i  ∂    ∂t∂ ρACv T + ∂x mA Cv T + c2 = −kL (T − T amb ) .

(3.22)

3.3 Metode Numerik Pada subbab ini akan dibahas skema numerik yang akan digunakan pada Persamaan (3.22) untuk mengetahui distribusi aliran yang bersifat transien sepanjang pipa. Sebelumnya akan dibahas terlebih dahulu mengenai analisis dimensi, syarat awal dan syarat batas.

3.3.1 Analisis Dimensi Analisis dimensi yang dilakukan disini adalah mengubah besaran menjadi besaran tidak berdimensi dengan tujuan menyederhanakan model yang akan diselesaikan secara numerik. Akan dilakukan analisis dimensi pada Persamaan (3.22),


41

dengan memilih beberapa besaran sebagai acuan. Untuk besaran panjang, dipilih panjang pipa (L) sebagai acuan, lalu untuk rapat massa dipilih rapat massa di inlet (ρ0 ) sebagai acuan, sedangkan untuk temperatur dipilih temperatur lingkungan (T amb ) sebagai acuan dan untuk kecepatan dipilih kecepatan suara (c) sebagai acuan. Selain itu, ada besaran yang dibuat tak berdimensi terhadap besaran acuan yang telah ditentukan di atas, seperti fluks massa, besaran ini akan dibuat tak berdimensi terhadap fluks massa di inlet (m0 ) dengan m0 = ρ0 c. Selain itu, besaran waktu, besaran ini akan dibuat tak berdimensi terhadap t0 dengan t0 = L/c. Apabila ruas kanan t0 = L/c dikalikan dengan

ρ0 ρ0 ,

maka t0 =

Lρ0 m0 .

Secara umum, analisis dimensi

dapat diringkas sebagai berikut : e x=

x ρ T m t e= ,e ρ = , Te = ,m ,e t= . L ρ0 T amb m0 t0

Dengan mensubstitusi besaran yang telah dibuat tak berdimensi pada Persamaan (3.22) yang pertama, diperoleh ∂e ρ ∂e m + =0. x ∂e t ∂e

(3.23)

Sedangkan dengan mensubstitusikan besaran yang telah dibuat tak berdimensi untuk Persamaan (3.22) yang kedua, diperoleh 2 e m ρ e |e −L fg m m| ∂e m ∂ eρ +e + = . ∂e x 2De ρ ∂e t

(3.24)

Dan yang terakhir, akan disubstitusikan besaran yang telah dibuat tak berdimensi untuk Persamaan (3.22) yang ketiga, dengan sebelumnya membuat persamaan tersebut menjadi lebih sederhana, yaitu dengan membagi persamaan tersebut dengan


42

ACv , sehingga menjadi, " !# ∂ ∂ c2 −kL (T − T amb ) . ρT + m T+ = ∂t ∂x Cv ACv

Akan dibuat pemisalan, yaitu λ1 =

c2 Cv

dan λ2 =

kL ACv .

Dengan mensubstitusikan λ1 ,

λ2 , dan besaran yang telah dibuat tak berdimensi, maka akan diperoleh λ1 e ∂e m T + e ∂e ρT λ2 t0 e T amb + =− T −1 . ∂e x ρ0 ∂e t

(3.25)

Persamaan (3.23), (3.24) dan (3.25) memuat semua variabel yang sudah tidak berdimensi lagi. Ketiga persamaan tersebut yang akan digunakan dalam skema numerik. Namun, untuk kemudahan notasi, tandae. akan dihilangkan, sehingga penulisannya menjadi,

   ∂ρ ∂m     ∂t + ∂x = 0 ,    2    ∂m ∂ mρ +ρ −L fg m|m|  + ∂x = 2Dρ ,  ∂t     λ   ∂m T + T 1   amb  + = − λρ20t0 (T − 1) .  ∂ρT ∂t ∂x

(3.26)

3.3.2 Syarat Awal Dalam melakukan proses numerik, dibutuhkan syarat awal. Pada kasus ini, proses aliran bersifat tunak pada kondisi awalnya. Oleh karena itu, keadaan tunak digunakan sebagai syarat awal. Proses aliran bersifat tunak dalam arti sifat fluidanya tidak mengalami perubahan terhadap waktu. Apabila direpresentasikan ke dalam bentuk persamaan, maka menjadi ∂ρ ∂m ∂(ρT ) = 0, = 0, =0. ∂t ∂t ∂t

(3.27)


43

Akan disubstitusikan Persamaan (3.27) ke dalam Persamaan (3.26), yaitu menjadi    ∂m   =0,   ∂x  2      ∂ mρ +ρ −L fg m|m|  , = 2Dρ  ∂x     λ   ∂m T + T 1   amb  = − λρ20t0 (T − 1) .  ∂x

Dari Persamaan (3.28) yang pertama,

∂m ∂x

(3.28)

= 0 memberi arti bahwa fluks

massa bernilai konstan sepanjang pipa, karena yang diketahui adalah nilai fluks massa di inlet yaitu m0 , maka untuk keadaan tunak nilai fluks massa sepanjang pipa konstan sebesar nilai fluks massa di inlet yaitu m0 . Dengan hubungan antara fluks massa dan laju alir, akan diperoleh distribusi laju alir sepanjang pipa untuk keadaan tunak. Distribusi laju alir untuk keadaan tunak bersifat konstan sepanjang pipa, sama seperti nilai fluks massa untuk keadaan tunak.

Sedangkan untuk Persamaan (3.28) kedua,

∂

m2 ρ +ρ

∂x

=

−L fg m|m| 2Dρ

dengan meng-

ganti m dengan m0 diperoleh 1 2 −L fg m0 2 x − m0 ln ρ + ρ = +C . 2 2D 2

(3.29)

Dengan mensubstitusi ρ0 = 1, maka akan diperoleh C = 1/2, sehingga apabila disubstitusikan ke dalam Persamaan(3.29) diperoleh persamaan akhir, yaitu : f (ρ) =

2D ln ρ D 2 − ρ −1 − x = 0 . L fg L f g m0 2

(3.30)

Masalah di atas sama dengan mencari akar fungsi terhadap ρ. Dalam tugas akhir ini, digunakan metode Newton Raphson untuk mencari akar fungsi terhadap ρ, dengan x yang merupakan panjang pipa akan dibagi menjadi beberapa segmen, misal-


44

kan J segmen, dengan tiap segmen mempunyai panjang ∆x. Dengan proses tersebut dan hubungan antara rapat massa dan tekanan yang direpresentasikan melalui persamaan keadaan, maka akan diperoleh distribusi tekanan sepanjang pipa untuk keadaan tunak (Gambar 3.4). Distribusi tekanan untuk keadaan tunak bersifat menurun sepanjang pipa.

Gambar 3.4: Tekanan pada Keadaan Tunak.

Dan terakhir untuk persamaan (3.28) ketiga,

λ ∂m T + T 1 amb

∂x

= − λρ20t0 (T − 1) de-

ngan mengganti m dengan m0 , diperoleh m0 ln(T − 1) = −

λ2 t0 x+C . ρ0

(3.31)

Dengan mensubstitusikan temperatur di inlet yang telah dibuat tak berdimensi yaitu 0 0 T = TTamb , maka diperoleh C = m0 ln TTamb − 1 , sehingga apabila disubstitusikan ke dalam Persamaan (3.31) dan membentuk kedalam fungsi temperatur T terhadap x,


45

diperoleh persamaan akhir, yaitu : λ t

− ρ 2m0 x

T (x) = 1 + exp

0 0

! T0 −1 . T amb

(3.32)

Dengan Persamaan (3.32) akan diperoleh distribusi temperatur sepanjang pipa pada keadaan tunak (Gambar 3.5). Distribusi temperatur bersifat menurun menuju temperatur lingkungan, setelah mencapai temperatur lingkungan, nilai temperatur tidak dapat turun lagi.

Gambar 3.5: Temperatur pada Keadaan Tunak.

3.3.3 Syarat Batas Pada dasarnya syarat batas diperoleh dari masalah di lapangan. Dalam tugas akhir ini, syarat batas yang diketahui adalah nilai laju alir di inlet dan di outlet (Gambar 3.6) dan (Gambar 3.7) , yang akan dikonversikan ke dalam fluks massa.


46

Gambar 3.6: Laju Alir di Inlet Waktu Simulasi 7 jam. Data diberikan pada Tabel 3.3. Selain itu, diketahui juga syarat batas untuk temperatur di inlet yaitu bernilai konstan. Waktu 0 − 1 jam 1 − 2 jam 2 − 3 jam 3 − 4 jam 4 − 5 jam 5 − 6 jam 6 − 7 jam

Laju Alir di Inlet 204.211 204.780 Turun secara linear 0 0 Naik secara linear 285.620

Laju ALir di Outlet 191.42 175.192 166.53 176.556 163.224 166.564 160.842

Satuan MMS CF/D MMS CF/D MMS CF/D MMS CF/D MMS CF/D MMS CF/D MMS CF/D

Tabel 3.3: Syarat Batas.

Untuk nilai yang tidak diketahui pada batasnya, dalam kasus ini adalah rapat massa (ρ) dan (ρT ) di outlet, dapat diperoleh dengan cara diskritisasi Persamaan (3.26) untuk persamaan yang pertama dan ketiga. Sebelumnya L yang merupakan panjang pipa akan dibagi menjadi beberapa segmen, misalkan J segmen, dengan


47

tiap segmen mempunyai panjang ∆x dan t yang merupakan waktu proses terjadinya transien akan dibagi menjadi beberapa segmen, misalkan N segmen, dengan tiap segmen mempunyai panjang ∆t. Maka notasi ρn0 dan ρnJ yang berturut - turut adalah rapat massa gas di inlet dan outlet pada waktu ke-n. Proses diskritisasi persamaan (3.26) untuk persamaan yang pertama dan ketiga, yaitu :

1. Diskritisasi untuk rapat massa (ρ) di inlet : n ρn+1 0 = ρ0 +

∆t n m0 − mn1 . ∆x

(3.33)

2. Diskritisasi untuk rapat massa (ρ) di outlet : n ρn+1 J = ρJ +

∆t n m J−1 − mnJ . ∆x

(3.34)

3. Diskritisasi untuk (ρT ) di outlet : ρn+1 T Jn+1 = ρnJ T Jn − J

n 1 mnJ T Jn + Tλamb − mnJ−1 T J−1 + − λρ20t0 T Jn − 1 ∆t .

∆t ∆x

h

λ1 T amb

i

(3.35)

3.3.4 Skema Lax-Wendroff Pada sub bab ini, akan dijelaskan mengenai skema numerik yang digunakan dalam penyelesaian. Sebelumnya perhatikan Persamaan (3.26). Persamaan tersebut dapat ditulis ke dalam bentuk vektor seperti, ~ ∂F( ~ U) ~ ∂U ~ , + = ~r(U) ∂t ∂x

(3.36)


48

Gambar 3.7: Laju Alir di Outlet Waktu Simulasi 7 jam. dengan     ρ   m       2 m ~ =  m  , F( ~ U) ~ =  U    ρ + ρ    1 m T + Tλamb ρT

        ~  , ~r(U) =     

0 −L fg m|m| 2ρD −λ2 t0 ρ0 (T −

1)

     .  

Skema numerik yang akan digunakan merupakan skema Lax-Wendroff dua langkah, dengan stencil (Gambar 3.8). Adapun skema Lax-Wendroff dua langkah yaitu : U~n+11 = j+ 2

∆t 1 ~n F~ U~nj+1 − F~ U~nj +~r(U~nj )∆t , U j+1 + U~nj − 2 2∆x

U~n+1 = U~nj − j

∆t ~ ~n+1 ~ ~n+1 F U 1 − F U 1 +~r(U~nj )∆t , j+ 2 j− 2 ∆x

(3.37)

(3.38)


49

dengan j = 1, 2, ..., J − 1, dan   ρn j   n n ~ U j =  m j   n n ρj T j

       ~n ~n   , F j (U j ) =     

mnj 2 ρnj

mnj + ρnj

mnj T nj +

λ1 T amb

        ~ =   , ~r(U)     

  0   −L fg mnj mnj  n  2ρ j D  −λ2 t0  n − 1 .  T ρ0 j

Gambar 3.8: Stencil Skema Lax-Wendroff Dua Langkah. Notasi U nj menyatakan rapat massa (ρ), fluks massa (m) dan (ρT ) di segmen ke- j pada step waktu ke-n.

Pemodelan Matematika dan Metode Numerik

Recommend Documents