Buku Pemrograman Komputer
Buku Pemrograman Komputer
Buku Pemrograman Komputer
Y.C. Danarto
Bregas Siswahjono T.S.
Paryanto
LEMBAR PENGESAHAN
1. Judul buku ajar
2. Tim Penyusun
3. Waktu penulisan
: 4 ( empat ) bulan.
4. Biaya
Mengetahui/Menyetujui,
Pembantu Rektor I UNS
Adrian Nur
Y.C. Danarto
Bregas Siswahjono T.S.
Paryanto
Adrian Nur
Y.C.. Danarto
Bregas Siswahjono T.S.
Paryanto
Edisi Pertama
Edisi pertama, Cetakan pertama, November 2005
KATA PENGANTAR
Tim Penulis
iv
DAFTAR ISI
Kata Pengantar
Daftar Isi
I.
II.
III.
IV.
iii
v
1
2
9
16
23
2.1.
2.2.
2.3.
23
31
36
43
3.1.
3.2.
3.3.
43
48
52
INTEGRAL NUMERIS
75
4.1.
4.2.
4.3.
76
80
85
Aturan Trapesium
Aturan Simpson
Kuadratur Gauss
vi
V.
PENCOCOKAN KURVA
5.1.
5.2.
VI.
VII.
100
101
112
119
121
122
127
131
OPTIMASI
133
6.1.
6.2.
6.3.
133
136
146
Pemograman Linier
Metode Golden Section
Optimasi Multivariabel
169
7.1.
170
170
173
174
176
206
207
212
227
8.1.
228
228
236
237
237
247
7.2.
VIII.
99
8.2.
Daftar Pustaka
251
PENGANTAR
PEMROGRAMAN MATLAB
2
telah memiliki pengalaman dalam pemrograman bahasa lain seperti bahasa C,
PASCAL, atau FORTRAN.
MATLAB merupakan perangkat lunak (software) yang dikembangkan
oleh Mathworks, Inc. (lihat http://www.mathworks.com) dan merupakan
perangkat lunak yang paling efisien untuk perhitungan numerik berbasis matriks.
Dengan demikian jika di dalam perhitungan kita dapat memformulasikan
masalah ke dalam format matriks, maka MATLAB merupakan pilihan terbaik
untuk penyelesaian numeriknya.
Dalam penggunaannya MATLAB hanya bersaing dengan MATCAD
atau MAPLE. Namun demikian sifatnya yang fleksibel menjadikannya lebih
banyak digunakan oleh para praktisi maupun akademisi dari pada yang lain.
MATLAB merupakan bahasa pemrograman tingkat tinggi berbasis
matriks yang sering digunakan untuk teknik komputasi numerik dalam
menyelesaikan masalah-masalah yang melibatkan operasi matematika elemen,
matrik, optimasi, aproksimasi, dan lain-lain. Secara umum MATLAB banyak
digunakan pada bidang-bidang :
3
MATLAB Command Window/Editor merupakan jendela yang dibuka
pertama kali setiap kali MATLAB dijalankan.
Pada jendela ini dapat dilakukan akses-akses ke perintah-perintah
MATLAB dengan mengetikkan ekspresi MATLAB, seperti mengakses help
window, helpdesk, demo, dan lain-lainnya.
Jika variabel-variabel dan perintah-perintah pada layar command
window yang sedang aktif akan disimpan, maka dapat dilakukan dengan
menggunakan perintah diary. Sebagai contoh, jika ingin menyimpan output
m =
1
4
2
5
3
6
di direktori c:\coba dengan nama file data1.txt, maka dapat dilakukan dengan
mengetik :
m=[1 2 3;4 5 6]
m =
1
2
3
4
5
6
diary 'c:\coba\data1.txt'
Apapun yang tertulis di command window setelah perintah ini akan disimpan
sampai kita menutup file ini dengan mengetikkan :
diary off
pwd
4
mkdir
what
: digunakan untuk melihat nama m-file (yang akan dibahas kemudian) dalam direktori aktif.
who
whos
delete
clear
clc
doc
demo
>> edit
5
Coba anda ketikkan script berikut ini pada command window dan
perhatikan hasilnya tergambar pada figure window.
x=0:2:360;
y=sin(x*pi/180);
plot(x,y)
dan kernudian menekan enter. Di layar akan muncul informasi dalam bentuk teks
pada layar MATLAB yaitu:
SIN
Sine.
SIN(X) is the sine of the elements of X.
Overloaded methods
help sym/sin.m
Bentuk help yang lain adalah Help Window yaitu dengan mengetikkan
helpwin
Informasi selengkapnya dapat diperoleh dengan mengklik dua kali salah satu
baris yang ada di MATLAB Help Window atau dapat juga dengan mengetikkan
informasi yang diinginkan pada kotak di sebelah kiri atas MATLAB Help
Window.
Help yang lain lagi adalah dalam bentuk Help Desk.
Anda bisa
dan anda akan dapatkan berbagai kemudahan dalam memperoleh informasi yang
diinginkan.
OPERATOR MATEMATIKA
MATLAB dapat melakukan operasi-operasi aritmatika dasar berikut:
No.
OPERASI
Pemangkatan
Perkalian; Pembagian
Penambahan; Pengurangan
SIMBOL
^
* ; / atau \
+ ; -
7
help elfun
menyimpan variabel x dan y dalam format biner di file data dalam format ASCII.
Untuk membuka data gunakan perintah load.
Contoh:
load data.mat
8
sernua ekspresi relasi dan logika adalah satu untuk benar dan nol untuk salah
dengan tipe array logika, yaitu hasilnya mernuat bilangan 1 dan 0 yang tidak saja
dapat digunakan untuk statemen matematika akan tetapi dapat juga digunakan
untuk pengalamatan.
Operator relasi MATLAB menunjukkan perbandingan:
OPERATOR
RELASI
<
DISKRIPSI
Kurang dari
>
Lebih dari
<=
>=
Sama dengan
~=
DISKRIPSI
&
AND
OR
NOT
9
1.
http://www.mathworks.com/
3.
http://dir.yahoo.com/science/mathematics/software/matlab/
http://www.cse.uiuc.edu/cse301/matlab.html
String
String dalam MATLAB adalah tipe data yang terdiri dari huruf-huruf
dan/atau nilai-nilai ASCII yang ditampilkan representasinya. String adalah teks
yang diawali dan diakhiri dengan tanda (......).
Contoh:
g='selamat datang'
g =
selamat datang
size(g)
ans =
1
14
whos
Name
Size
Bytes Class
ans
1x2
16 double array
g
1x14
28 char array
Grand total is 16 elements using 44 bytes
10
Setiap karakter dalam suatu string adalah satu elernen dalam array, dengan setiap
elemennya berukuran 2 byte.
Untuk melihat representasi ASCII karakter string dapat dilakukan
dengan
operasi
aritmetik
terhadap
string
atau
mengkonversikannya
109
97
116
32
100
97
116
109
97
116
32
100
97
116
Karena string merupakan array numerik dengan atribut khusus, string dapat
dimanipulasi menggunakan sernua metode manipulasi array yang tersedia, dalam
MATLAB.
t=g(1:5)
t =
selam
11
Fungsi-fungsi string yang lain adalah:
Fungsi input
nama=input('Nama :','s')
Nama : jojon
nama =
jojon
Fungsi disp
Fungsi
disp
memungkinkan
untuk
menampilkan
string
tanpa
atau
disp('jojon')
jojon
num2stro
52
54
Jelas bahwa h+2 bukan 236, karena h mewakili nilai ASCII dari string 2, 3 dan 4
yaitu masing-masing 50, 51 dan 52
Skalar
Skalar adalah sebuah bilangan tunggal yang secara matriks berukuran
1x1. Operasi skalar merupakan dasar matematika. Kalkulator biasa melakukan
perhitungan skalar. Tentu saja MATLAB bisa melakukan dengan lebih baik.
12
Namun jika dalam sekali waktu kita ingin melakukan operasi yang sama pada
beberapa bilangan, perulangan operasi skalar akan menghabiskan waktu dan
tidak praktis. Untuk itu MATLAB menyediakan operasi pada tipe data array
yang akan dibahas berikutnya.
Array
Untuk menghitung nilai, y=sin(x); 0 x . Tidak mungkin
menghitung semua nilai pada interval tersebut. Kita bisa mengevaluasi nilainya
pada setiap jarak 0,1. Untuk mudahnya kita membuat tabel:
x
0.1
0..2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0.31
0.59
0.81
0.95
1.0
0.95
0.81
0.59
0.31
atau
x=[0:0.1:1]*pi
1.8850
13
atau
x=linspace(0,pi,11)
Perkalian skalar
Penambahan array
a+b=[a1+b1 a2+b2
Perkalian array
Pemangkatan array
...
an+bn]
Matriks
Matriks merupakan bentuk utama MATLAB. Seperti haInya array,
matriks juga didefenisikan elemen demi elemen. Untuk memahami operasi
matriks dasar, perhatikan contoh-contoh berikut ini.
a=[1 2 3] %penentuan matriks a
a =
1
2
3
b=[1 2 1; 3 2 1; 2 1 3] %penentuan matriks b
b =
1
2
1
3
2
1
2
1
3
14
Tanda titik koma memisahkan antar baris matriks, sedangkan antar kolom cukup
dengan spasi.
c=a*b %perkalian a dan b
c =
13
9
12
Pada perkalian matriks di atas yang perlu diingat adalah jumlah kolom a harus
sama dengan jumlah baris b. Jika tidak demikian maka MATLAB akan
memberikan pesan kesalahan:
d=b*a
??? Error using ==> *
Inner matrix dimensions must agree.
Jika b bukan matriks bujur sangkar, maka akan ada pesan kesalahan
det(a)
??? Error using ==> det
Matrix must be square.
inv(b) %invers b, b harus matriks bujursangkar
ans =
-0.5000
0.5000
0
0.7000
-0.1000
-0.2000
0.1000
-0.3000
0.4000
f=a' %transpose matriks a
f =
1
2
3
15
3
Elemen pertama pada variabel ans adalah jumlah baris dan kedua adalah jumlah
kolom.
Terdapat berbagai fungsi matriks yang tersedia dalam MATLAB yang
bisa anda pelajari dengan mengetikkan:
help matfun
Tentu saja ini bisa dipahami jika anda sudah memahami matematika matriks.
Selain operasi matriks dasar MATLAB juga menyediakan operasi
antar elemen yang sangat berguna. Dot atau titik (.) adalah tanda yang
digunakan unutk operasi-operasi tersebut. Sebagai contoh operasi matriks X dan
matriks Y meliputi X.^Y; X.*Y dan X./Y. Pada operasi ini tiap elemen
dioperasikan dengan elemen yang bersesuaian posisinya dalam matriks. Dengan
demikian ukuran matriks X dan Y harus sama.
Bentuk ini tidak perlu digunakan dalam penjumlahan atau pengurangan
karena sudah merupakan operasi antar elemen.
X=[1 2;3 4]
X =
1
2
3
4
Y=[2 3;1 4]
Y =
2
3
1
4
X.*Y
% perkalian antar elemen
ans =
2
6
3
16
X*Y
% perkalian standar matriks
ans =
4
11
10
25
Z=[2 3]
Z =
2
3
X.*Z
??? Error using ==> .*
Matrix dimensions must agree.
1.3. M-FILE
16
Untuk menyelesaikan masalah yang cukup sederhana, mengetikkan
beberapa perintah langsung di command window memang cukup cepat dan
efektif. Tetapi jika jumlah perintahnya sangat banyak atau jika anda ingin
mengubah beberapa variabel dan mengulang kembali perhitungannya, maka
mengetikkan perintah langsung di command window akan sangat menjemukan.
Untuk itu, MATLAB menyediakan suatu struktur untuk membuat fungsi anda
sendiri atau suatu teknik pemrograman dalam bentuk M-File. M-File berasal dari
aturan bahwa nama file ini harus dikhiri dengan ekstensi .m; misalnya
gambar1.m
M-File juga memudahkan pemakai untuk menyimpan dan memanggil
kembali perhitungan yang sering digunakan. Perintah diary yang biasanya
digunakan untuk menyimpan file yang berisi barisan perintah dan keluaran pada
command window hanya sekedar merekam tetapi tidak bisa langsung digunakan
dalam perhitungan jika dipanggil kembali.
Membentuk M-File
Untuk membuat M-File, klik File pada menu pull-down di MATLAB
command window, lalu pilih New dan klik M-File. Di layar akan muncul
MATLAB Editor/Debugger. Selanjutnya di layar ini dapat ditulis argumenargumen yang diinginkan, dapat diedit penulisannya dan sebagainya. Setelah
selesai melakukan pengetikan, klik File pada layar MATLAB Editor/Debugger
dan pilih Save As... . Beri nama yang Anda inginkan untuk file tersebut,
misalnya contoh11.m kemudian klik Save. Pastikan file yang disimpan dalam
direktori yang mudah untuk dipanggil.
Berikut ini contoh sederhana script file.
Contoh 1.1
% Script file contoh11.
% Kurva fungsi y=(sin x)/x
x = pi/100:pi/100:10*pi;
y = sin(x)./x;
plot(x,y)
grid
17
Untuk mengesekusi script file di atas maka pada Command Window ketikkan
nama file yang telah disimpan sebelumnya misal
contoh11.
Dengan syarat
0.8
0.6
0.4
0.2
-0.2
-0.4
0
10
15
20
25
30
35
Narna fungsi dan narna file harus identik. Contoh: fhitung disimpan dalam
file yang bernama fhitung.m
2.
18
melibatkan pemanggilan ke fungsi M-File yang lain, fungsi M-File yang
dipanggil itu juga akan dikompilasi ke dalam memori.
3.
Setiap fungsi memiliki ruang kerja sendiri yang berbeda dengan ruang kerja
MATLAB. Satu-satunya hubungan antara ruang kerja MATLAB dengan
variabel-variabel dalam fungsi adalah variabel-variabel input dan output
fungsi. Jika suatu fungsi mengubah nilai dalam bentuk suatu variabel input,
perubahan itu hanya tampak dalam fungsi dan tidak mempengaruhi ruang
kerja MATLAB.
5.
Jika suatu fungsi dipanggil, jumlah argument input dan output yang
digunakan hanya terdapat dalam fungsi tersebut.
6.
Fungsi dapat berbagi variabel dengan fungsi lain, ruang ke.ja MATLAB,
dan
pemanggilan
rekursi
untuk
dirinya
sendiri
jika
variabelnya
8.
19
berfungsi sebaliknya. Perintah input memungkinkan untuk meminta input dari
pemakai saat M-File dijalankan.
Secara umum fungsi M-File didefinisikan menggunakan command
function. Sintaks standar untuk command function adalah:
function[output1,output2,...]=NamaFunction(inputl,input2,..)
Jika kita menginginkan keterangan tentang M-file tersebut, kita ketikkan pada
Command Window:
help gerak_parabola
Menjalankan M-File
Untuk memanggil atau mengeksekusi M File, yang perlu dilakukan
adalah dengan memindahkan path search dari MATLAB compiler. Pada
dasarnya, proses eksekusi compiler MATLAB adalah dengan mencan command
atau definisi operator yang ada dan mengeksekusi definisi script atau operator
pertama yang ditulis dan ditemui direktori MATLAB (di direkton bin atau
toolbox).
20
Set Path dapat dilakukan dengan cara membuka menu pull down di
MATLAB command editor, pilih File kemudian pilih Set Path. Dari menu
tersebut arahkan current directory ke direktori tempat di mana disimpan script
yang ingin dieksekusi dengan cara mengetikkan atau browsing directory ke
tempat penyimpanan script yang akan dieksekusi.
Ada beberapa cara menjalankan fungsi M-File. Misalkan fungsi
gerak_parabola (contoh 12) dapat dijalankan dengan cara, antara lain:
[jarak,tinggi]=gerak_parabola(pi/4,40,9.8,1)
x =
28.2843
y =
23.3843
Bukan
y=V0*sin(sudut)*t-g/2*t^2;
21
Contoh 13
g=inline('6-x^2')
g =
Inline function:
g(x) = 6-x^2
g(2)
ans =
2
Contoh 14
radius=inline('sqrt(x^2+y^2)','x','y')
radius =
Inline function:
radius(x,y) = sqrt(x^2+y^2)
radius(2,3)
ans =
3.6056
PERSAMAAN ALJABAR
LINIER
24
a 11
a 21
a
31
a 11
x3 =
x2 =
x1 =
a 12
a 22
a 32
a 13 | b1
a '23 | b '2
''
a 33
| b 3''
a 12
a '22
b 3''
''
a 33
(b
'
2
a 13 | b 1
a 23 | b 2
a 33 | b 3
,
a '23 x 3
a '22
(b
),
a 12 x 2 a 13 x 3
.
a 11
Persamaan dasar
a11x1 + a12x2 + a13x3 = b1
.... (2.1)
.... (2.2)
.... (2.3)
1
1
1
a12x2 +
a13x3 =
b1
a 11
a 11
a 11
.... (2.4)
a
a
a
a 22 21 a 12 x 2 + a 23 21 a 13 x 3 = b 2 21 b1
a 11
a 11
a 11
a
a
a
a 32 31 a 12 x 2 + a 33 31 a 13 x 3 = b 3 31 b1
a 11
a 11
a 11
.... (2.5)
25
'
'
a 32
x 2 + a 33
x 3 = b 3'
.... (2.6)
x2 +
a '22
a '23 x 3 =
1
a '22
b '2
.... (2.7)
[Pers.(2.6) - a 32 x Pers.(2.7)]
'
'
'
a '22
a '22
''
a 33
x 3 = b 3''
.... (2.8)
b 3''
''
a 33
(b
'
2
a '23 x 3
a '22
(b
a 12 x 2 a 13 x 3
a 11
26
%
Eliminasi Gauss
for k=1:n
for i=k+1:n
m(i,k)=A(i,k)/A(k,k);
for j=k:n
A(i,j)=A(i,j)-m(i,k)*A(k,j);
end
B(i)=B(i)-m(i,k)*B(k);
end
end
%
Substitusi Balik
x(n)=B(n)/A(n,n);
for i=n-1:-1:1
sum=B(i);
for k=i+1:n
sum=sum-A(i,k)*x(k);
end
x(i)=sum/A(i,i);
end
%S
%P
% abu
0,2
0,05
1,0
0,06
0,5
0,03
0,7
0,03
x2 +
x3 +
x4
0,2x1 +
x2 +
0,5x3 +
0,7x4
0,61
0,05x1 +
0,06x2 +
0,03x3 +
0,03x4
0,043
2x1 +
3x2 +
x3 +
x4
1,8
x1 +
x2 +
x3 +
x4
2x1 +
10x2 +
5x3 +
7x4
6,1
5x1 +
6x2 +
3x3 +
3x4
4,3
Proses penskalaan
27
2x1 +
3x2 +
x3 +
x4
1,8
Campuran Batubara
Hasil Penyusunan Persamaan Aljabar Linier
x1 +
x2 + x3 + x4 = 1
2x1 + 10x2 + 5x3 + 7x4 = 6.1
5x1 + 6x2 + 3x3 + 3x4 = 4.3
2x1 + 3x2 + x3 + x4 = 1.8
Penyelesaian dilakukan dengan Eliminasi Gauss
-------------------------------------------------------Surakarta, Oktober 2005
Jurusan Teknik Kimia, Fak. Teknik
Universitas Sebelas Maret
********************************************************
clear all
clc
%
Penyusunan Matriks A
A = [1 1 1 1
2 10 5 7
5 6 3 3
2 3 1 1];
%
Penyusunan Matriks B
B = [1 6.1 4.3 1.8]';
n = 4;
Orde matriks
x = elimgauss(A,B,n);
disp(' ')
disp('Fraksi masing-masing type batubara')
disp('----------------------------------')
x1 = x(1)
x2 = x(2)
x3 = x(3)
x4 = x(4)
Keluaran program
Fraksi masing-masing type batubara
---------------------------------x1 =
0.2000
x2 =
0.3000
x3 =
0.4000
x4 =
0.1000
Sehingga batubara
type 1 = 20 %
type 2 = 30 %
type 3 = 40 %
type 4 = 10 %
28
(2.2.1)
yi = 1
(2.2.2)
yi PT =
xi Pio
H1 x1 + H2 x2 + H3 x3 = H
(2.2.3)
(2.2.4)
Penyelesaian :
Persamaan (2.2.1)
x 1 + x2 + x3 = 1
Substitusi persamaan (2.2.3) ke persamaan (2.2.2)
PO
P1O
PO
x1 + 2 x2 + 3 x3 = 1
PT
PT
PT
(2.2.5)
Persamaan (2.2.4)
20 x1 + 30 x2 + 55 x3 = 30
Persamaan aljabar simultan
(2.2.6)
29
x 1 + x2 + x3 = 1
(2.2.1)
30 x1 + 20 x2 + 15 x1 = 24
(2.2.5)
20 x1 + 30 x2 + 55 x3 = 30
(2.2.6)
clear all
clc
%
Data-data
P = [5/2 5/3 5/4];
H = [20 30 55];
Pt = 2;
Ht = 30;
%
%
%
%
atm
kcal
atm
kcal
%
Penyusunan Matriks A
A = [1 1 1
P(1)/Pt P(2)/Pt P(3)/Pt
H(1) H(2) H(3)];
%
Penyusunan Matriks B
B = [1 1 Ht]';
n = 3;
Orde matriks
x = elimgauss(A,B,n);
disp(' ')
disp('Fraksi Mol masing-masing komponen')
disp('---------------------------------')
x1 = x(1)
x2 = x(2)
x3 = x(3)
30
.
a n1
a 12
a 22
a 13
a 23
...
...
a 32
a 33
...
a n2
a n3
...
a 1n
a 2 n
a 3n
a nn
x 1 b1
x b
2 2
x 3 = b3
. .
x n b n
. (2.9)
x = B
Suatu matriks jika dikalikan dengan matriks identitas maka akan dihasilkan
matriks itu sendiri. Matriks identitas adalah matriks bujursangkar yang terdiri
dari bilangan 0 pada setiap elemennya kecuali pada diagonal matriks yang
mempunyai bilangan 1.
Contoh matriks identitas 3 x 3
1 0 0
0 1 0 = I
0 0 1
Jika matriks identitas tersebut dikalikan dengan matriks lain maka akan
dihasilkan matriks itu sendiri. Contohnya jika matriks tersebut dikalikan dengan
B yaitu suatu matriks yang berukuran 1 x 3.
1 0 0 b 1
b1
0 1 0 b = b
2
2
0 0 1 b 3
b 3
= B
31
A x
= B
-1
= A-1 B
A A x
I x = A-1 B
= A-1 B
(2.10)
a
Inversi matriks A ukuran 2 x 2 11
a 21
A-1 =
a 11a 22
a 12
dapat ditentukan dengan cara
a 22
a 22
1
a 12 a 21 a 21
a 12
a 11
inv.
Contoh 2.3. Konstanta Kecepatan Reaksi Phase Gas Nitrit Oksida dengan
Oksigen
Data konstanta kecepatan reaksi phase gas nitrit oksida dengan oksigen
T, K
300
7,1
413
4,0
564
2,8
k = A T exp (-E/T)
Tentukan A, m, dan E
Penyelesaian
k = A Tm exp (-E/T)
ln k = ln A + m ln T E/T
Dari data-data di atas
22,68336 = ln A + 5,70378 m 0,00333 E
22,10956 = ln A + 6,02345 m 0,00242 E
21,75289 = ln A + 6,33505 m 0,00177 E
untuk
ln A
= x1
= x2
32
E
= x3
clear all
clc
%
Data-data
T = [ 300 413 564];
k = [ 7.1 4.0 2.8]*10^9;
% K
% cc/(mol.det)
%
Penyusunan Matriks AA
AA = [ 1 log(T(1)) 1/T(1)
1 log(T(2)) 1/T(2)
1 log(T(3)) 1/T(3)];
%
Penyusunan Matriks BB
BB = log(k);
CC = inv(AA)*BB';
disp(' ')
disp('Parameter konstanta kecepatan reaksi')
disp('------------------------------------')
A = exp(CC(1))
m = CC(2)
E = -CC(3)
33
A =
1.3509e+007
m =
0.6064
E =
-841.6950
= K 1
kf a dp
De
K2
D
e
K3
Sh = K1 (Re)K2 (Sc)K3
Tentukan K1, K2, dan K3 dengan data-data berikut:
Sh
Re
Sc
(Sherwood)
(Reynold)
(Schmidt)
43,7
10800
0,6
21,5
5290
0,6
24,2
3120
1,8
Penyelesaian
Sh = K1 (Re)K2 (Sc)K3
log (Sh) = log (K1) + K2 log (Re) + K3 log (Sc)
Dari tabel-tabel di atas, diperoleh
1,64048 = x1 + 4,03342 x2 0,22185 x3
1,33244 = x1 + 3,72346 x2 0,22185 x3
1,38382 = x1 + 3,49415 x2 + 0,25527 x3
disusun dalam bentuk matriks
34
1 4,03342 0,22185 x 1 1,64048
1 3,72346 0,22185 x = 1,33244
2
1 3,49415 0,25527 x 3 1,38382
Program Penyelesaian
%
%
%
%
%
%
%
%
%
%
clear all
clc
%
Data-data
Sh = [43.7 21.5 24.2];
Re = [10800 5290 3120];
Sc = [0.6 0.6 1.8];
%
Penyusunan Matriks A
A = [1 log10(Re(1)) log10(Sc(1))
1 log10(Re(2)) log10(Sc(2))
1 log10(Re(3)) log10(Sc(3))];
%
Penyusunan Matriks B
B = [log10(Sh(1))
log10(Sh(2))
log10(Sh(3))]
k = A\B;
disp(' ')
disp('Konstanta Parameter Transfer Massa ')
disp('-----------------------------------')
k1 = 10^k(1)
k2 = k(2)
k3 = k(3)
Keluaran Program
Konstanta Parameter Transfer Massa
----------------------------------k1 =
0.0058
k2 =
0.9938
k3 =
0.5853
35
(2.11)
(2.12)
Contoh 2.5. Neraca Massa Steady State pada Kolom Distilasi Bertingkat
Xylena, styrena, toluena, dan benzena dipisahkan dengan kolom distilasi seperti
ditunjukkan pada gambar 2.1.
Tentukan kecepatan alir D1, D2, B1, B2, D, dan B.
36
D1
7 % Xylena
4 % Styrena
54 % Toluena
35 % Benzena
D
#2
B1
F = 70 mol/menit
18 % Xylena
24 % Styrena
42 % Toluena
16 % Benzena
#1
D2
15 % Xylena
25 % Styrena
40 % Toluena
20 % Benzena
15 % Xylena
10 % Styrena
54 % Toluena
21 % Benzena
B
#3
B2
24 % Xylena
65 % Styrena
10 % Toluena
1 % Benzena
clear all
clc
37
%
Penyusunan Matriks A
A = [0.07 0.18 0.15 0.24
0.04 0.24 0.10 0.65
0.54 0.42 0.54 0.1
0.35 0.16 0.21 0.01];
%
Penyusunan Matriks B
B = [0.15*70 0.25*70 0.40*70 0.20*70]';
disp(' ')
disp('Penyelesaian untuk D1 B1 D2 B2 adalah : ')
%
Dekomposisi Matriks A
[L,U] = lu(A);
%
Membentuk Matriks B*
Bstar = inv(L)*B;
%
Menentukan Matriks x
X = inv(U)*Bstar
disp('Umpan kolom 2')
D1 = X(1);
B1 = X(2);
D = D1 + B1
X_Dx = (A(1,1)*D1+A(1,2)*B1)/D
X_Ds = (A(2,1)*D1+A(2,2)*B1)/D
X_Dt = (A(3,1)*D1+A(3,2)*B1)/D
X_Db = (A(4,1)*D1+A(4,2)*B1)/D
disp('Umpan kolom 3')
D2 = X(3);
B2 = X(4);
B = D2 + B2
X_Dx = (A(1,3)*D2+A(1,4)*B2)/B
X_Ds = (A(2,3)*D2+A(2,4)*B2)/B
X_Dt = (A(3,3)*D2+A(3,4)*B2)/B
X_Db = (A(4,3)*D2+A(4,4)*B2)/B
Keluaran program :
Penyelesaian untuk D1 B1 D2 B2 adalah :
X =
26.2500
17.5000
8.7500
17.5000
Umpan kolom 2
D =
43.7500
X_Dx =
0.1140
X_Ds =
0.1200
X_Dt =
0.4920
X_Db =
0.2740
Umpan kolom 3
38
B =
26.2500
X_Dx =
0.2100
X_Ds =
0.4667
X_Dt =
0.2467
X_Db =
0.0767
V (L/mol)
9/2
16/3
Program Matlab
%
%
%
%
%
%
%
%
%
%
Persamaan Virial
PV = a + bP + cP^2
Untuk mengevaluasi konstanta-konstanta a, b, dan c,
telah dilakukan percobaan dengan data-data di bawah ini
--------------------------------------------------------------Surakarta, Oktober 2005
Jurusan Teknik Kimia, Fak. Teknik
Universitas Sebelas Maret
***************************************************************
clear all
clc
%
Data-data
P = [1 2 3];
V = [4 9/2 16/3];
%
Matriks A
A = [1 P(1) (P(1))^2
1 P(2)
(P(2))^2
1 P(3) (P(3))^2];
%
Matriks B
B = [P(1)*V(1)
% atm
% L/mol
39
P(2)*V(2)
P(3)*V(3)];
C = A\B;
disp(' ')
disp('Konstanta persamaan virial adalah : ')
disp('------------------------------------')
a = C(1)
b = C(2)
c = C(3)
Keluaran Program
Konstanta persamaan virial adalah :
-----------------------------------a =
1
b =
2
c =
1
tabel.
Tentukan persamaan kecepatan reaksi tersebut !
Ca
Cb
0,7
0,2
0,4567
0,9
0,3
1,0707
1,6
0,2
2,5917
40
%
Data-data
Ca = [0.7 0.9 1.6];
Cb = [0.2 0.3 0.2];
ra = [0.4567 1.0707 2.5917];
% Matriks A
A = [1 log10(Ca(1)) log10(Cb(1))
1 log10(Ca(2)) log10(Cb(2))
1 log10(Ca(3)) log10(Cb(3))];
% Matriks B
B = [log10(ra(1))
log10(ra(2))
log10(ra(3))];
konst = A\B;
disp(' ')
disp('Konstanta adalah : ')
disp('-------------------')
k = 10^konst(1)
a = konst(2)
b = konst(3)
Keluaran Program
%
Persamaan Virial
%
ra=k*(Ca)^a*(Cb)^b
%
log(ra) = log(k) + a*log(Ca) + b*log(Cb)
%
Untuk mengevaluasi konstanta-konstanta a, b, dan k,
%
telah dilakukan percobaan dengan data-data di bawah ini
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
echo off
Konstanta adalah :
------------------k =
3.4990
a =
2.1000
b =
0.7998
PERSAMAAN ALJABAR
NON LINIER
44
dicari nilai f(x) dari kedua titik dan dari titik tengah interval. Dengan
menganalisa nilai f(x) ketiga titik maka kita dapat membuang setengah interval
dalam setiap langkahnya.
Langkah penyelesaian :
1.
2.
3.
Interval yang benar akan menghasilkan f(xA) dan f(xB) pada daerah yang
berbeda (sebelah kanan dan kiri akar persamaan yang dicari)
Akar persamaan
xA
xB
5.
xA + xB
2
Hitung f(xM)
a.
. (3.1)
45
xA
=
xM
xA
xB
baru
x M = xA
Baru
xA
xB
Uji apakah f(xm) < Toleransi. Jika tidak, ulangi langkah (4) untuk
menentukan xm yang baru.
2p
(1 x ) (2 + x )
46
Program matlab untuk metode bagidua :
%
Persamaan aljabar non linier
%
%
Penyelesaian dilakukan dengan Metode Bagidua (bisection)
%
%
-------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
**************************************************************
clc
clear all
disp 'Masukkan Perkiraan Batas Bawah dan Batas Atas '
xa=input('batas bawah = ');
xb=input('batas atas = ');
fa = contoh31(xa);
fb = contoh31(xb);
x=[xa:(xb-xa)/50:xb];
plot(x,contoh31(x));
xlabel('x')
ylabel('f(x)')
grid on
while fa*fb > 0
disp 'Tebakan salah, ganti batas atas dan/atau batas bawah'
disp 'Perhatikan grafik!! Perkirakan nilai x untuk f(x)
mendekati nol'
xa=input('batas bawah = ');
xb=input('batas atas = ');
fa=contoh31(xa);
fb=contoh31(xb);
end
tol=0.000001;
xm=(xa+xb)/2;
fm=contoh31(xm);
while abs(xa-xb)>tol
xm=(xa+xb)/2;
fm=contoh31(xm);
if fa*fm<0
xb=xm;
fb=contoh31(xb);
else
xa=xm;
fa=contoh31(xa);
end
end
t =['f(x)=0 untuk x = 'num2str(xm)];
disp (t)
47
p = 0.2;
kp = 0.4568;
f_x = x./(1-x).*sqrt(2*p./(2+x))-kp;
%atm
Dalam memberikan masukan untuk batas bawah dan atas, secara teori nilai
konversi berkisar antara 0 0,99. Untuk x = 1, nilai kp akan tak terdefinisi.
Keluaran program :
Masukkan Perkiraan Batas Bawah dan Batas Atas
batas bawah = 0
batas atas = 0.99
f(x)=0 untuk x = 0.53487
1 y + y 2 y3
1 y3
dengan y = b/4v, untuk b adalah koreksi van der Waals dan v adalah volum
molar. Jika z = 0,892 berapakah y ?
Penyelesaian dilakukan dengan program bisection yaitu dengan mengganti
fungsi disosiasi H2O dengan fungsi untuk persoalan faktor kompresibilitas gas
ideal.
function f_y = contoh32(y)
%
Faktor Kompresibilitas Gas Ideal
%
1 - y + y^2 - y^3
%
z = ------------------%
1 - y^3
%
z = 0,896
%
%
%
z = 0.892;
f_y = (1-y+y.^2-y.^3)-z.*(1-y.^3);
Secara umum ada 3 nilai yang memenuhi persamaan di atas karena f(y)
merupakan persamaan orde 3. Tetapi nilai untuk y yang tepat adalah antara 0 1.
Masukkan Perkiraan Batas Bawah dan Batas Atas
48
batas bawah = 0
batas atas = 1
f(x)=0 untuk x = 0.1229
f(xi)
Akar persamaan
Garis
singgung
f(xi+1)
xi+1
xi
xi+1 = xi - f ( x i )
f ' (x i )
Langkah penyelesaian
1.
2.
3.
4.
.(3.2)
49
| f(xi) | < toleransi,
Jika ya, maka selesai.
Jika tidak, xi baru = xi+1, ulangi langkah (2)
Adakalanya f(x) sulit dicari dengan analitis. Untuk itu dapat dilakukan
pendekatan secara numeris
f(x)= f ( x 0 + ) f ( x 0 )
. (3.3)
B
untuk pO dalam mmHg
T+C
dan T dalam OC
Benzena
Toluena
6,89745
6,95334
1206,35
1343,94
220,237
219,377
50
%
%
%
%
%
%
%
%
%
clc
clear all
xnew=input(' Nilai untuk tebakan awal =
xold=0;
tol=0.0001;
eps=0.0001;
');
Proram penyelesaian
function fT=contoh33(T)
%
Temperatur Dew Point untuk Campuran Benzena dan Toluena
%
Komposisi uap y1 = 0,77 dan y2 = 0,33
%
Tekanan uap murni
%
o
(
B
)
%
Pi = 10^( A + ------- )
%
(
(T + C) )
%
Hukum Roult-Dalton
%
yi * P
%
xi = -------%
o
%
Pi
%
x1 + x2 = 1
%
%
%
P=760;
% Konversi dari atm ke mmHg
y1=0.77;
y2=0.23;
p1o=10^(6.89745-1206.35/(T+220.237));
p2o=10^(6.95334-1343.94/(T+219.377));
x1=y1*P/p1o;
x2=y2*P/p2o;
fT=x1+x2-1;
Keluaran program
Nilai untuk tebakan awal =
x-old
f(x-old)
85
51
85
0.17029
89.639
0.014053
90.0946
0.0001176
90.0985
8.4065e-009
akar persamaan, x = 90.0985
f(x) = 8.4065e-009
Sehingga T = 90,0985 OC
Contoh 3.4. Hubungan Faktor Friksi suatu Pelarut dengan Bilangan
Reynolds
Hubungan faktor friksi untuk aliran suatu pelarut dengan bilangan Reynolds (Re)
secara empiris adalah
1
f
5,6
1
= ln Re f + 14
k
k
Re = 3750;
k = 0.28;
f_f = ((1/k)*log(Re*sqrt(f))+(14-5.6/k))*sqrt(f)-1;
Keluaran program
Nilai untuk tebakan awal = 0.01
x-old
f(x-old)
0.01
0.51676
0.0044847
-0.080152
0.005105
-0.0020744
0.0051219
-1.4377e-006
akar persamaan, x = 0.0051219 f(x) = -1.4377e-006
52
sehingga f = 0,0051219
fzero.
dengan :
adalah nama suatu fungsi yang berisi persamaan yang
akan dicari nilai nol-nya.
xo adalah nilai perkiraan awal
tol adalah ketepatan penyelesaian jika trace mempunyai nilai lebih
dari 1.
trace adalah jumlah iterasi yang dilakukan.
Dua parameter terdepan harus diisi, sehingga alternatif untuk menggunakan
nama_fungsi
= fzero(nama_fungsi,xo)
RT
a
V b V (V + b ) + b( V b)
%
%
%
%
%
%
%
%
Persamaan Peng-Robinson
Penyelesaian dilakukan dengan fungsi fzero
--------------------------------------------------------------Surakarta, Oktober 2005
Jurusan Teknik Kimia, Fak. Teknik
Universitas Sebelas Maret
53
%
***************************************************************
% Data-data
clear all
clc
global R T P
R=8.3137;
T=340;
P=10^4;
%
%
%
m^3.kPa/(kgmol.K)
K
kPa
R.T
P = -----V - b
Parameter utk CO2
a
------------------V(V + b) + b(V - b)
%
Nama File : F35.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global R T P
%Parameter Peng-Robinson untuk CO2
a=364.61;
b=0.02664;
fV=R*T/(V-b)-a/(V*(V+b)+b*(V-b))-P;
Keluaran program
V =
1.6792e-001
54
T
q=
Cp
dt
To
q = -2616;
To = 550;
% BTU/lbmol
% oF
%
Integral secara analitis
fT=7.053*(T-To)+1.2242*10^-3/2*(T^2-To^2)...
-2.6124*10^-7/3*(T^3-To^3)-q;
Contoh-contoh lain
Contoh 3.7. Debit Aliran dengan Pompa pada Pipa
Suatu cairan akan dialirkan dari tangki 1 ke tangki 2 melalui pipa berdiameter D,
dengan bantuan pompa. Panjang ekuivalen pipa, Le. Dari persamaan Bernoulli
antara titik 1 dan titik 2 diperoleh persamaan berikut :
z 2 z1 +
f.Le.v 2
Hm = 0
2.g.D
f=
0,0596
Re 0,215
Re =
.v.D
dengan
55
Karakteristik pompa sentrifugal yang dapat dipakai berupa hubungan antara head
pompa ( Hm, cm) dengan debit ( Q, cm3/dtk) dapat didekati dengan persamaan :
Hm = 3718,5 2,34967.Q + 7,8474.10-4.Q2 9,5812.10-8.Q3
Debit aliran dihitung dengan persamaan :
Q=
.D 2 .v
%
%
%
%
%
%
%
Densitas, g/cm^3
Viskositas, g/cm.dtk
grafitasi, cm/dtk^2
tinggi titik 1, cm
tinggi titik 2, cm
diameter pipa, cm
panjang ekivalen, cm
% Perhitungan
vhit=fzero(@F37,200)
Qhit=pi/4*D^2*vhit
Program terkait
function fv = F37(v)
%
Langkah perhitungan
%
tebak v ---> hitung bil. Reynold ---> hitung faktor friksi
%
---> hitung debit ---> hitung head pompa --->
%
hitung pers. Bernoulli
%
f Le v^2
%
z2 - z1 + -------- - Hm = 0
%
2 g D
%
%
%
56
Q=pi/4*D^2*v;
% Perhitungan Debit
Hm=3718.5-2.3496*Q+7.8474e-4*Q^2-9.5812e-8*Q^3;
% Perhitungan head pompa
fv=z2-z1+f*Le*v^2/(2*g*D)-Hm; % Persamaan Bernoulli
Keluaran program
vhit =
227.6735
Qhit =
2.8610e+003
T = TF -
0,25.(1 - z).C pA
0,25.z. H FR
+ (0,75 0,5.z).C pB + 0,25.C pC
dengan :
B
K
ln
A
dengan K =
0,25.(1 z)
.P
(1 0,5.z)
(0,75 0,5.z)
PB =
.P
(1 0,5.z)
0,25.z
.P
PC =
(1 0,5.z)
PA =
PC
PA .PB
57
Data-data yang diketahui :
Tref = 298 K; A = 8.10
-6
HRo
global dHr Tf
% Data-data
dHr=-20000;
% Panas reaksi, kal/gmol
Tf=400; P=20;
% Temperatur umpan, K
Tref=298;
% Temperatur referensi, K
% Parameter temperatur kesetimbangan
A=8e-6;
% atm
B=4500;
% K
% Kapasitas panas
Cpa=7;
% kal/mol.K
Cpb=8;
% kal/mol.K
Cpc=12;
% kal/mol.K
% Perhitungan konversi kesetimbangan
konversi=fzero(@F38,0.6)
Program terkait
function fz = F38(z)
%
Langkah perhitungan
%
tebak konversi kesetimbangan z ---> hitung tekanan parsial
%
komponen a, b, dan c ---> hitung panas reaksi
%
---> hitung temperatur reaktor ---> hitung temperatur
%
kesetimbangan ---> bandingkan T reaktor dan T ktmbgn
%
T kesetimbangan - T reaktor = 0
%
%
%
%
58
fz=Tstb-T;
Keluaran program
konversi =
0.4143
3764
PAo = exp14,95
4497
PBo = exp16,07
4934
PCo = exp16,27
global ya yb yc Pt
% Data-data tekanan parsial
ya=0.4; yb=0.3; yc=0.3;
% Tekanan total, cmHg
Pt=76;
% Perhitungan temperatur
T=fzero(@F39,300)
Program terkait
59
function fT=F39(T)
%
Tekanan uap murni
%
o
(
B )
%
Pi = 10^( A + --- )
%
(
T )
%
Hukum Roult-Dalton
%
yi * P
%
xi = -------%
o
%
Pi
%
x1 + x2 + x3 = 1
%
%
%
global ya yb yc Pt
Pao=exp(14.95-3764/T);
Pbo=exp(16.07-4497/T);
Pco=exp(16.27-4934/T);
fT=Pt*(ya/Pao+yb/Pbo+yc/Pco)-1;
Keluaran program
T =
390.2252
4g p D p
v t =
3C D
(3.10.1)
CD =
24
Re
24
1 + 0,14 Re 0,7
Re
untuk
0,1 Re 1000
60
CD = 0,44
untuk
4
1000< Re 35000
untuk
35000 Re
Tentukan terminal velocity untuk partikel batubara dengan P = 1800 kg/m3 dan
Dp = 0,208.10-3 m yang jatuh dalam air pada T = 298,15 K dengan = 994,6
kg/m3 dan = 8,931.10-4 kg/mdetik.
Penyelesaian
Tebak vt hitung Re Tentukan CD Cek vt dengan persamaan (3.10.1 )
Program Matlab
function fv_t = contoh310(v_t)
%
Terminal Velocity untuk Partikel Jatuh dalam Fluida
%
Tebak vt --> hitung Re --> Tentukan CD --> Cek vt
%
/ 4.g.(rhop - rho).Dp
%
vt = / ------------------%
\/
3.CD.rho
%
%
%
%
Data-data
rho_p = 1800;
D_p = 0.208*10^-3;
T = 298.15;
rho = 994.6;
mu = 8.931*10^-4;
g = 9.80665;
%kg/m^3
%m
%K
%kg/m^3
%kg/m/s
%m/s^2
%
Menghitung bilangan Reynold
Re = D_p*v_t*rho/mu;
%
Menentukan koef. drag (CD)
if Re < 0.1
C_D = 24/Re;
elseif Re < 1000
C_D = 24*(1+0.14*Re^0.7)/Re;
elseif Re < 350000;
C_D = 0.44;
else
C_D = 0.19-80000/Re;
end
%
Membandingkan vt awal dan vt perhitungan
fv_t = v_t^2*(3*C_D*rho)-4*g*(rho_p-rho)*D_p;
61
1.5782e-002
Hukum gas ideal dapat menunjukkan hubungan tekanan, volum, dan temperatur
(PVT) hanya pada tekanan gas rendah (mendekati tekanan atmosfer). Untuk
tekanan tinggi persamaan yang harus digunakan lebih kompleks daripada hukum
gas ideal tersebut. Perhitungan volum molar dan faktor kompressibilitas
menggunakan persamaan yang lebih kompleks yang memerlukan penyelesaian
secara numeris. Salah satu persamaan yang dapat digunakan adalah persamaan
van der Waals
a
P + 2 (V b ) = RT
V
dengan
a=
27 R 2TC 2
64 PC
b=
RTC
8PC
dengan
P
= tekanan (atm)
= temperatur (K)
TC
PC
P
PC
PV
RT
Tentukan volum molar dan faktor kompressibilitas untuk gas ammonia pada
tekanan P = 56 atm dan temperatur T = 450 K menggunakan persamaan van
der Waals !
b.
62
c.
Penyelesaian
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
Volum Molar dan Faktor Kompressibilitas utk Pers van der Waals
Varibel berubah
P
Pr = ----Pr
Dari data-data yang ada
Hitung P
Hitung V dengan terlebih dahulu menebak V,
lalu hitung dengan pers. van der Waals
Hitung Z
P.V
Z = -----R.T
--------------------------------------------------------------Surakarta, Oktober 2005
Jurusan Teknik Kimia, Fak. Teknik
Universitas Sebelas Maret
***************************************************************
clear all
format short e
global P a b R T
%
Data-data
Pc = 111.3;
%
Tc = 405.5;
%
R = 0.08206;
%
T = 450;
%
P = 56;
%
Pr = [P/Pc 1 2 4 10 20];
a = 27/64*R^2*Tc^2/Pc;
b = R*Tc/(8*Pc);
atm
K
atm.L/gmol.K
K
atm
%
Perhitungan
for j = 1:6
P = Pc*Pr(j);
Vtebak = R*T/P;
V = fzero('waals',Vtebak);
z = P*V/(R*T);
hasil(j,1)=Pr(j);
hasil(j,2)=V;
hasil(j,3)=P*V/(R*T);
end
%
Menampilkan hasil perhitungan
disp('
Pr
Volum molar
faktor Z')
disp(hasil)
plot(hasil(:,1),hasil(:,3))
title('Faktor Kompresibilitas vs Tekanan Tereduksi')
xlabel('Tekanan Tereduksi')
ylabel('Faktor Kompresibilitas')
= RT
63
%
Nama File : waals.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global P a b R T
% Modifikasi persamaan van der Waals
f_V = P*V^3-P*b*V^2-R*T*V^2+a*V-a*b;
Faktor kompresibilitas
Hasil perhitungan
Volum molar
5.7489e-001
2.3351e-001
7.7268e-002
6.0654e-002
5.0875e-002
4.6175e-002
faktor Z
8.7183e-001
7.0381e-001
4.6578e-001
7.3126e-001
1.5334e+000
2.7835e+000
64
Neraca energi pada batang logam
Eacc = Eg Eout
Ecc = kecepatan energi tersimpan pada batang logam
Eg = kecepatan pembangkitan energi dalam batang karena panas hambatan
Eout = kecepatan aliran energi total keluar dari batang karena konveksi dan radiasi
T lingkungan
Panas yang dilepaskan
I, kuat arus
D, diameter
L, panjang pipa
dT
dt
Eg = I2RLL
Dengan RL adalah hambatan per satuan panjang dan satuan persamaan di atas
adalah (amp)2ohm.
Eout = hA(T - T) + A(T4 T4sur)
Dengan A adalah luas permukaan perpindahan panas
T adalah temperatur fluida
Tsur adalah temperatur lingkungan karena perpindahan panas radiasi
= emisivitas material
= konstanta Stefan boltzmann = 5,67.108 W/m2-K4
Substitusi ke persamaan
cV
dT 2
= I RLL - hA(T - T) + A(T4 T4sur)
dt
65
dengan
A = DL
V = D2L/4
Kondisi awal T(0) = T0
Jika arus mula-mula adalah nol, kemudian temperatur awal T0 maka didapatkan
persamaan
hA(T - T) + A(T04 T4sur) = 0
jika T dan Tsur sama, maka T0 = T = Tsur.
D = 1 mm
= 0,8
RL = 0,4 ohm/m
h = 100 W/m2K
T = Tsur = 25 oC
%
%
%
%
%
%
%
clear all
global D emiss R_L sig h Tinf Tsur Isq
D = 0.001;
% Diameter kawat (m)
emiss = 0.8;
% Emisivitas
R_L = 0.4;
% Risistansi per satuan panjang (ohm/m)
sig = 5.67e-8;
% Konstanta Stefan-Bolzmann (W/m^2K^4)
h = 100;
% Koefisien transfer panas (W/m^2K)
Tinf = 25+273.15;
% Temperatur fluida (K)
Tsur = Tinf;
% Temperatur lingkungan untuk radiasi (K)
NI = 21;
I = linspace(0,10,NI);
Tss = zeros(size(I));
Isq = I(1)*I(1);
Tg = Tinf;
Tss(1) = fzero('neracaenergi',Tg);
%
Perhitungan Temperatur
for n=2:NI
Isq = I(n)*I(n);
Tss(n) = fzero('neracaenergi',Tss(n-1));
end
Tssa = Tinf +I.*I*R_L/(h*pi*D);
%
Plot hasil perhitungan
plot(I,Tss-273.15,'k',I,Tssa-273.15,'k-.','LineWidth',2), grid on
rr = axis; rr(3) = 0; axis(rr);
title('Temperatur Konduktor Steady State vs Arus')
66
xlabel('Arus (amps)'), ylabel('Temperatur Steady State (C)')
legend('T(ss) dengan radiasi','T(ss) tanpa radiasi')
%
Perhitungan neraca energi
A = pi*D;
gen = I.*I.*R_L;
conv = h*A*(Tss-Tinf);
rad = emiss*sig*A*(Tss.^4-Tsur^4);
bal = (gen-conv-rad);
%
Menampilkan hasil
fprintf(1,'\n\n')
fprintf(1,'
HASIL PERHITUNGAN \n')
fprintf(1,' Parameter perancangan \n')
fprintf(1,'
Diameter kawat (m)
= %8.3f \n',D)
fprintf(1,'
Emisivitas
= %8.3f
\n',emiss)
fprintf(1,'
Resistansi (ohm/m)
= %8.3f
\n',R_L)
fprintf(1,'
Konstanta Stefan-Bolzmann (W/m^2K^4)= %9.2e
\n',sig)
fprintf(1,'
Koefisien Transfer Panas (W/m^2K)
= %8.3f \n',h)
fprintf(1,'
Temperatur Fluida (C)
= %8f
\n',Tinf-273.15)
fprintf(1,'
Temperatur Lingkungan (C)
= %8f
\n',Tsur-273.15)
fprintf(1,'\n')
fprintf(1,' Komponen Neraca Energi vs Arus \n')
fprintf(1,'
I
Energi yg Muncul
Konveksi
Radiasi Neraca
Temperatur \n')
fprintf(1,' (A)
(W/m)
(W/m)
(W/m)
(W/m)
C
\n')
for n = 1:NI
fprintf(1,' %4.1f
%8.2f
%8.2f
%8.2f %8.4f %8.2f
\n',...
I(n),gen(n),conv(n),rad(n),bal(n),Tss(n)-273)
end
67
Keluaran program
HASIL PERHITUNGAN
Parameter perancangan
Diameter kawat (m)
=
0.001
Emisivitas
=
0.800
Resistansi (ohm/m)
=
0.400
Konstanta Stefan-Bolzmann (W/m^2K^4)= 5.67e-008
Koefisien Transfer Panas (W/m^2K)
= 100.000
Temperatur Fluida (C)
= 25.000000
Temperatur Lingkungan (C)
= 25.000000
Komponen Neraca Energi vs Arus
I
Energi yg Muncul
Konveksi
(A)
(W/m)
(W/m)
0.0
0.00
0.00
0.5
0.10
0.10
1.0
0.40
0.38
1.5
0.90
0.86
2.0
1.60
1.52
2.5
2.50
2.38
3.0
3.60
3.43
3.5
4.90
4.66
4.0
6.40
6.08
4.5
8.10
7.68
5.0
10.00
9.47
5.5
12.10
11.44
6.0
14.40
13.59
6.5
16.90
15.92
7.0
19.60
18.42
7.5
22.50
21.09
8.0
25.60
23.93
8.5
28.90
26.93
9.0
32.40
30.09
9.5
36.10
33.41
10.0
40.00
36.88
Radiasi
(W/m)
0.00
0.00
0.02
0.04
0.08
0.12
0.17
0.24
0.32
0.42
0.53
0.66
0.81
0.98
1.18
1.41
1.67
1.97
2.31
2.69
3.12
Neraca Temperatur
(W/m)
C
0.0000
25.15
0.0000
25.45
0.0000
26.36
-0.0000
27.88
0.0000
30.00
0.0000
32.73
0.0000
36.06
0.0000
39.98
-0.0000
44.50
-0.0000
49.60
-0.0000
55.30
-0.0000
61.57
-0.0000
68.41
0.0000
75.81
0.0000
83.77
0.0000
92.28
0.0000
101.32
-0.0000
110.88
-0.0000
120.94
-0.0000
131.50
0.0000
142.53
68
69
Neraca energi pada sistem pada keadaan steady :
p1
v2
p
v2
+ z1 + 1 + h A h R h L = 2 + z1 + 2
2g
2g
LV 2
D 2g
Persamaan menjadi
V2
LV 2
= h+Lf
2g
D 2g
V2 =
2g(h + L )
L
1 + f
D
0,25
/ D 5,74
+ 0,9
log10
3,7 Re
Pipa 1
Pipa 1 D3 = 0,1342 ft
Pemograman Matlab
70
%
%
%
%
%
%
%
clear all
global DD h LL g visc rough
%
Data-data, D dan L bervariasi
D = [0.0518 0.0874 0.1342];
%
h = 5;
%
g = 32.2;
%
visc = 1.21e-5;
%
rough = 1.5e-4;
%
NL = 50;
L = linspace(1/12,20,NL)';
ND = length(D);
V = zeros(NL,ND);
%
Perhitungan
for n=1:ND
DD = D(n);
Vo = sqrt(2*g*h);
LL = L(1);
V(1,n) = fzero('friksi_pipa',Vo);
for i=2:NL
LL = L(i);
V(i,n) = fzero('friksi_pipa',V(i-1,n));
end
end
%
Menampilkan hasil perhitungan
plot(L,V(:,1),'k -',L,V(:,2),'k+',L,V(:,3),'k.','lineWidth',2)
title('Friksi Pipa : Velocity vs Panjang Pipa')
xlabel('Panjang pipa (ft)')
ylabel('Velocity Aliran Keluar (ft/detik)')
legend('pipa Sch 40 1/2"','pipa Sch 40 1"','pipa Sch 40 1 1/2"')
Keluaran program
71
INTEGRAL NUMERIS
. (4.1)
I = f ( x ) dx
a
x
a
76
x
a
Dx
f (x) dx
x
(f (a ) + f ( x 1 ) ) + x (f ( x 1 ) + f ( x 2 ) ) + ...... + x (f ( x n 1 ) + f (b) )
2
2
2
x
[f (a ) + 2f (x 1 ) + 2f ( x 2 ) + ...... + 2f (x n 1 ) + f (b)]
2
(b a )
[f (a ) + 2f (x 1 ) + 2f (x 2 ) + ...... + 2f (x n 1 ) + f (b)]
2n
.(4.2)
Dalam Matlab, metode perhitungan integrasi dengan aturan trapesium dapat
menggunakan fungsi trapz.
77
Penggunaan fungsi trapz adalah
Z = trapz(x,y)
D = F = 10 lbmol/jam
78
Neraca massa benzena
F xF D yD = M
dt =
dx w
dt
M
dx w
F x F D yD
x w
1 ( 1) x w
adalah
xw
t=
M
dx w
F
x
D yD
F
xwo
79
% Nilai t pada berbagai xw
yd = alpha.*xw./(1+(alpha-1).*xw);
ti=-M./(F.*xf-D.*yd);
% Jumlah t pada seluruh bagian
t=trapz(xw,ti)
10
11
17
24
32
41
48
51
v(10) =
10
a(t) dt + v(0) =
a(t) dt
Program Matlab
%
%
%
%
%
%
%
%
%
%
%
Kecepatan accelerometer
Hubungan waktu (t) dengan percepatan a (m^2/detik)
ditunjukkan pada tabel
Kecepatan secara numeris
v = int ( a dt)
--------------------------------------------------------------Surakarta, Oktober 2005
Jurusan Teknik Kimia, Fak. Teknik
Universitas Sebelas Maret
***************************************************************
% Data-data
t=[0:10];
a=[0 2 4 7 11 17 24 32 41 48 51];
% detik
% m^2/detik
v(1)=0;
% keadaan awal
80
for k=[2:11]
v(k)=trapz(t(1:k),a(1:k));
end
% Cetak Hasil
disp(' ')
disp('
t(detik) v(m/detik)')
disp([t' v' ])
Keluaran Program
t(detik)
0
1.0000
2.0000
3.0000
4.0000
5.0000
6.0000
7.0000
8.0000
9.0000
10.0000
v(m/detik)
0
1.0000
4.0000
9.5000
18.5000
32.5000
53.0000
81.0000
117.5000
162.0000
211.5000
f (x) dx = 3 (f 0 + 4f1 + f 2 )
x0
.(4.3.)
81
f (x)
f(x0)
f(x1)
f(x2)
h
x
x0
x1
x2
f (x) dx
h
(f 0 + 4f1 + f 2 ) + h (f 2 + 4f3 + f 4 ) + ...... + h (f 2n 2 + 4fx 2n 1 + f 2n )
3
3
3
h
[f 0 + 4(f1 + f 3 + ...... + f 2n 1 ) + 2(f 2 + f 4 + ...... + f 2n 2 ) + f 2n )]
3
.(4.4)
Untuk menggambarkan penggunaan aturan Simpson, akan ditampilkan dua
alternatif yaitu fungsi simp1 dan simp2.
simp1
82
%
---------------------------------------------------------------
if (m/2)~= floor(m/2)
disp('m harus genap '); break
end
h = (b-a)/m;
x = [a:h:b];
y = feval(func,x);
v = 2*ones(m+1,1);
v2 = 2*ones(m/2,1);
v(2:2:m) = v(2:2:m) + v2;
v(1) = 1;
v(m+1) = 1;
q = y*v;
q = q*h/3;
Alternatif 2
function q = simp1(func,a,b,m)
%
Integrasi Numeris
%
%
Penyelesaian dilakukan dengan Aturan Simpson
%
Menggunakan loop
%
q = simp2(func,a,b,m)
%
integrasi fungsi func dr a sampai b dg m pembagian
%
%
Nama File : simp2.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------if (m/2)~= floor(m/2)
disp('m harus genap '); break
end
h = (b-a)/m;
s = 0;
y1 = feval(func,a);
for j = 2:2:m
x = a+(j-1)*h;
ym = feval(func,x);
x = a+j*h;
yh = feval(func,x);
s = s+y1+4*ym+yh;
y1=yh;
end
q = s*h/3;
y = x 7 dx
0
83
function fv = f43(x)
% Fungsi pangkat sederhana
%
f(x) = x^7
%
%
Nama File : F43.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------fv = x.^7;
n = 2; i = 1;
t = clock;
disp(' n
nilai integral')
while n<51200
simpval=simp1('f43',0,1,n);
fprintf('%3.0f%14.9f\n',n,simpval);
n=2*n; i=i+1;
end
fprintf('\nwaktu = %4.2f detik ',etime(clock,t));
Keluaran program
n
2
4
8
16
32
64
128
256
512
1024
2048
4096
8192
16384
32768
nilai integral
0.171875000
0.129150391
0.125278473
0.125017703
0.125001111
0.125000070
0.125000004
0.125000000
0.125000000
0.125000000
0.125000000
0.125000000
0.125000000
0.125000000
0.125000000
84
%
***************************************************************
n = 2; i = 1;
t = clock;
disp(' n
nilai integral')
while n<51200
simpval=simp2('F43',0,1,n);
fprintf('%3.0f%14.9f\n',n,simpval);
n=2*n; i=i+1;
end
fprintf('\nwaktu = %4.2f detik ',etime(clock,t));
Keluaran program
n
2
4
8
16
32
64
128
256
512
1024
2048
4096
8192
16384
32768
nilai integral
0.171875000
0.129150391
0.125278473
0.125017703
0.125001111
0.125000070
0.125000004
0.125000000
0.125000000
0.125000000
0.125000000
0.125000000
0.125000000
0.125000000
0.125000000
85
Misalkan kendala titik-titik basis yang tetap ini diperbaiki dengan
menentukan dua titik pada tertentu kurva. Dengan menempatkan titik-titik ini
dengan bijaksana, dapat dibuat suatu garis lurus yang mengimbangi kesalahan
positif dan negatif, sehingga perkiraan integral dapat diperbaiki (Gambar 4.5.b).
Kuadratur Gauss adalah suatu teknik untuk melaksanakan strategi ini.
Integral dinyatakan dengan
1
i =1
I = f ( x ) dx = c i f ( x i )
.(4.5)
untuk n = 2, kita harus menentukan 4 parameter yaitu c1, c2, x1, dan x2.
1
I = f ( x ) dx = c1f(x1) + c2f(x2)
1
Aturan integrasi akan tepat untuk fungsi polinomial 1, x, x2, dan x3.
f(x) = 1
memberikan
=2
= c1 + c2
x dx
=0
= c1x1 + c2x2
= 2/3
= c1x12 + c2x22
=0
= c1x13 + c2x23
1 dx
f(x) = x
memberikan
f(x) = x2
memberikan
dx
f(x) = x3
memberikan
dx
86
f (x)
x
(a)
f (x)
x
(b)
sehingga
x1 = -
1
3
x2 =
1
3
1
1
1
I = f ( x ) dx = f + f
3
3
Prosedur yang sama dapat dilakukan dengan penggunaan titik yang lebih banyak.
87
Fungsi Matlab
quad
dan
quad8
Faktor Bobot
Argumen Fungsi
c1 =
1,000000000 x1 =
c2 =
1,000000000 x2 =
-0,577350269
0,577350269
c1 =
0,555555556 x1 =
-0,774596669
c2 =
0,888888889 x2 =
c3 =
0,555555556 x3 =
0,774596669
c1 =
0,347854845 x1 =
-0,861136312
c2 =
0,652145155 x2 =
-0,339981044
c3 =
0,652145155 x3 =
0,339981044
c4 =
0,347854845 x4 =
0,861136312
c1 =
0,236926885 x1 =
-0,906179846
c2 =
0,478628670 x2 =
-0,538469310
c3 =
0,568888889 x3 =
c4 =
0,478628670 x4 =
0,538469310
c5 =
0,236926885 x5 =
0,906179846
c1 =
0,171324492 x1 =
-0,932469514
c2 =
0,360761573 x2 =
-0,661209386
c3 =
0,467913935 x3 =
-0,238619186
c4 =
0,467913935 x4 =
0,238619186
c5 =
0,360761573 x5 =
0,661209386
c6 =
0,171324492 x6 =
0,932469514
88
%
% Variabel lain
%
p
... exponent/pangkat
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
global p
disp(' ')
disp('Perhitungan integral dari a sampai b (x^p dx)')
disp(' ')
% Input parameter
p=input('Masukan nilai pangkat p: ');
a=input('Masukan batas integrasi bawah a : ');
b=input('Masukan batas integrasi atas b : ');
% integral
s = quad8('F44', a, b);
% cetak hasil
disp(['Hasil perhitungan secara numeris : ', num2str(s)])
disp(['Hasil
perhitungan
secara
analitis
:
',
num2str(b^(p+1)/(p+1)-a^(p+1)/(p+1))])
Keluaran program
Perhitungan integral dari a sampai b (x^p dx)
Masukan nilai pangkat p: 2.5
Masukan batas integrasi bawah a :
Masukan batas integrasi atas b :
Hasil perhitungan secara numeris
Hasil perhitungan secara analitis
0
1
:
:
0.28571
0.28571
Contoh-contoh lain
Contoh 4.5. Kecepatan Reaksi pada Katalis Porous
Suatu proses dengan katalis porous mempunyai kecepatan reaksi
89
dC
= 0,7C 2
dt
1,0357 + 0,3173
1 + 0,4172
= 12 C
% mol/gr katalis
% mol/gr katalis
C = linspace(Cn,Co,101) ;
% Perhitungan integral
phi = 12*sqrt(C);
eta = (1.0357+0.3173*phi)./(1+0.4172*phi);
x=1./(0.7.*eta.*C.^2);
% karena kondisi batas dibalik,
% tanda negatif hilang
t = trapz(C,x);
% Cetak Hasil
disp(['Waktu yang diperlukan menurunkan konsentrasi dari Co = ',
num2str(Co), ' mol/gr katalis'])
disp(['
menjadi Cn = ', num2str(Cn), ' mol/gr katalis'])
disp(['
adalah ', num2str(t), ' detik'])
Keluaran program
Waktu yang diperlukan menurunkan konsentrasi dari Co = 2 mol/gr
katalis
menjadi Cn = 1 mol/gr katalis
adalah 0.89188 detik
90
T dalam oF dan Cp dalam Btu/lbmol oF. Jika panas yang dilepaskan
untuk menurunkan temperatur campuran gas panas tersebut dari 550 oF adalah
2616 Btu/lbmol gas sampai temperatur berapakah campuran gas tersebut dapat
didinginkan.
T
q=
Cp
dt
To
q = -2616;
To = 550;
% BTU/lbmol
% oF
%
Integral secara numeris aturan trapesium
T=linspace(Tn,To,1000);
cp=-(7.053+1.2242.*10^-3.*T-2.6124.*10^-7.*T.^2);
qtebak=trapz(T,cp);
fq=qtebak-q;
91
Fv, Cao, To
A ---> hasil
xout
Fv
C Ao
xout
xin
1
k (1 x ) 2
dx
k = A.exp E
RT
92
%
---------------------------------------------------------------
% Data-data
Fv=200;
Cao=5;
rho=1.1;
Cp=0.8;
A=3.12*10^8;
E=18600;
Hr=-15;
R=1.987;
Vol=8000;
xin=0;
xout=0.8;
%
%
%
%
%
%
%
%
%
%
%
% Perhitungan Volume
x=linspace(xin,xout,1000);
T=To-Cao*Hr/rho/Cp*x
;
k=A*exp(-E/R./T);
eq=1./k./(1-x).^2;
V=Fv/Cao*trapz(x,eq);
% Dibandingkan dg Volume sesungguhnya
FV=Vol-V;
x S = exp 8,8053
T
t=
kX
Ro
(x
0
dr
S x)
93
3
m = 4 r Nb / 3
T = To +
(mo - m)
x
(W + m) Cp
%
%
%
%
%
%
%
%
jari-jari awal, cm
densitas, g/mL
jumlah padatan
Temperatur awal, K
panas pelarutan, kal/g
gram solvent
kapasitas panas solvent
konstanta pelarutan, g/(cm2.dtk)
% Perhitungan
waktu=quad(@F48,0,Ro)
Program terkait
function area = F48(r)
%
Integral
%
Hitung berat mula-mula ---> hitung berat pada r
%
---> hitung T ---> hitung x ---> hitung t
%
%
rho
[Ro
dr
%
t = ----- int [ -------%
kx
[0 (xs - x)
%
%
%
Keluaran program
94
waktu =
180.0717
3333
x S = exp 8,8053
T
t=
kX
Ro
(x
0
dr
S x)
T = To +
(mo - m)
(W + m) Cp
Jika diketahui W = 100.000 g, kx = 0,01 g/(cm2.dtk) dan tdata = 200 dtk; tentukan
suhu masuk To. Suhu masuk To dapat diketahui dengan persamaan berikut,
t tdata = 0
Program Penyelesaian
%
%
%
%
%
%
%
95
tdat=200;
% waktu pelarutan
To=fzero(@F49,400)
Program terkait
function fTo = F49(To)
%
tebak To ---> Hitung berat mula-mula
%
---> hitung berat pada r ---> hitung T
%
---> hitung x ---> hitung t
%
rho
[Ro
dr
%
t = ----- int [ -------%
kx
[0 (xs - x)
%
---> bandingkan dg t data
%
%
%
Keluaran program
To =
347.2318
Contoh 4.10.
Suatu cairan akan dialirkan dari tangki dengan grafitasi melalui pipa berdiameter
D. Panjang ekuivalen pipa, Le. Tinggi permukaan cairan pada tangki mula-mula
za.
Waktu pengosongan tangki dinyatakan dengan persamaan berikut :
2 z
D 2f 1
t = t2 . .dz
D z 2i v
Persamaan Bernoulli:
f.Le.v 2 v 2
=0
2.g.D 2g
f=
0,0596
Re 0,215
dengan
96
Re =
.v.D
.D 2 .v
Program terkait
function fz=F410(z)
%
%
Nama File : F410.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global rho D Dt miu Le g
vo=0.4*sqrt(2*g*z);
v=fsolve(@F410F,vo,[],z);
fz=1./v;
Program lain
function fv=F410F(v,z)
%
%
Nama File : F410F.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global rho D Dt miu Le g
Re=rho*v*D/miu;
97
f=0.0596./Re.^0.215;
fv=z-f*Le.*v.^2/(2*g*D)-v.^2/(2*g);
Keluaran program
waktu =
1.9688e+004
PENCOCOKAN KURVA
Datadata yang dimiliki sering berupa harga-harga yang diskrit dari suatu
garis yang kontinyu. Sering kita membutuhkan data-data yang terletak di antara
harga-harga yang diskrit tersebut. Bagian ini akan mempelajari bagaimana cara
kita memperkirakan data-data untuk harga-harga yang terletak di antara hargaharga yang diskrit.
Ada dua pendekatan yang umum dilakukan, yaitu :
1.
2.
Interpolasi
Pendekatan interpolasi dilakukan jika data-data yang ada telah diketahui
dengan teliti dan valid ( misal data dari buku-buku referensi/handbook ).
Persamaan yang dibuat harus melalui titik-titik yang telah diketahui
tersebut.
100
(5.1)
n =1
i =1
Ei =
(yi ao a1xi)
(5.2)
i =1
Tetapi kriteria ini kurang baik. Gambar 5.1 yang melukiskan pencocokan
kurva terhadap dua titik. Sembarang garis lurus yang melalui titik tengah
dari garis penghubung itu akan meminimalkan harga pers (3.2) menjadi
nol, karena kesalahan akan saling meniadakan.
2.
101
n
Ei =
i =1
|yi ao a1xi|
(5.3)
i =1
o
o
3.
SSE =
Ei2 = (y i a o a1x i )
n
n =1
i=1
(5.4)
(5.5)
)2
SSE = E i2 = (y i a 0 a1x i )2
n
i =1
n =1
(5.6)
(5.7)
102
Minimasi SSE, jika SSE = 0 dan SSE = 0.
a 0
a 1
SSE = -2
a 0
( yi ao a1xi ) = 0
yi - ao - a1xi = 0
yi - n ao - a1 xi = 0
n ao + a1 xi = yi
SSE = -2
a 1
(5.8)
[( yi ao a1xi )xi] = 0
xiyi - xiao - a1xi2 = 0
xiyi - ao xi - a1 xi2 = 0
ao xi + a1 xi2 = xiyi
(5.9)
(5.10)
a0 = y a 1 x
(5.11)
n x i2 ( x i ) 2
KUANTIFIKASI KESALAHAN
Kuantifikasi kesalahan dalam regresi kuadrat terkecil dapat dinyatakan dengan
persamaan
R2 =
St Sr
St
...(5.12)
dan
St =
(yi - y )2
...(5.13)
Sr =
(yi a0 a1xi)2
...(5.14)
103
LINIERISASI HUBUNGAN YANG TIDAK LINIER
Hubungan variabel bebas dan variabel tak bebas tidak selalu merupakan
hubungan yang linier, tetapi dapat juga hubungan yang tidak linier namun dapat
dibawa ke bentuk hubungan yang linier. Proses ini disebut linierisasi.
Contoh linierisasi hubungan yang tidak linier
1.
y=axb
ln (y) = ln (a x b)
ln (y) = ln (a) + ln (xb)
ln (y) = ln (a) + b ln (x)
Hubungan variabel bebas dan variabel tak bebas menjadi linier dengan
ln (y) = y
ln (x) = x
ln (a) = a0
b
= a1
sehingga, y = a1x + a0
2.
y = a ebx
ln (y) = ln (a ebx)
ln (y) = ln (a) + b x
Hubungan linier, dengan
ln (y) = y
x
=x
ln (a) = a0
b
= a1
sehingga, y = a1x + a0
3.
y=a+
b
x
1
=x
x
a = a0
104
b = a1
sehingga, y = a1x + a0
4.
y=
x
a + bx
1 a + bx
=
y
x
1 a
=
+b
y x
Hubungan linier, dengan
1
=y
y
1
=x
x
b = a0
a = a1
sehingga, y = a1x + a0
5.
y=
ax
b+x
1 b+x
=
y
ax
1
a
1
=
+
y bx a
Hubungan linier, dengan
1
=y
y
1
=x
x
1
= a0
a
105
a
= a1
b
sehingga, y = a1x + a0
Contoh 5.1. Persamaan BET untuk proses adsorpsi
Dalam proses adsorpsi, hubungan logam yang dijerap, q(mg/g adsorben), dengan
konsentrasi kesetimbangan C (mg/L) dapat dihubungkan dengan persamaan BET
q=
q maks d C
(CS C )1 + (d 1) C
C
S
0,8
4,88
8,04
16,1
19,25
17,23
32,09
46,75
62,90
64,73
b.
Penyelesaian
q=
q maks d C
(C S C )1 + (d 1) C
CS
1
=
q
C S
(C S C )1 + (d 1) C
q maks dC
C
1 + (d 1)
CS
1
=
q(C S C )
q maks dC
C
1
(d 1)C
=
+
q(C S C ) q maks d q maks dC S
y = a0 + a1x
106
dengan
C
q (C S C )
=C
a0
1
q maks d
a1
(d 1)
q maks d C S
Program Matlab
% Persamaan BET untuk proses adsorpsi
%
qmaks d C
%
q = ------------------------------%
(Cs - C)[ 1 + (d - 1)*C/Cs ]
%
Linierisasi
%
C
1
(d-1)C
%
--------- = ------- + ---------%
q(Cs - C)
qmaks d
qmaks d Cs
%
C
%
y = --------x = C
%
q(Cs - C)
%
1
(d - 1)
%
a0 = ------a1 = ---------%
qmaks d
qmaks d Cs
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
% Data - data
C = [0.8 4.88 8.04 16.1 19.25];
q = [17.23 32.09 46.75 62.90 64.73];
adsorben
% Konsentrasi, mg/L
% Logam terjerap, mg/g
n=5;
Cs=49.46;
% jumlah data
% konstanta BET, mg/L
% Linierisasi
y = C./(Cs-C)./q;
x = C;
% Perhitungan regresi linier
xkuad=0;
xy=0;
for i=1:n
xkuad=xkuad+(x(i))^2;
xy=xy+x(i)*y(i);
end
% Konstanta regresi linier
a1=(n*xy-sum(x)*sum(y))/(n*xkuad-(sum(x))^2);
a0=mean(y)-a1*mean(x);
% Parameter BET
d=a1/a0*Cs+1
qmaks=1/a0/d
% Plot hasil
q=a0+a1*x;
plot(x,y,'k*',x,q,'k-','lineWidth', 2)
107
title('Perbandingan Data Percobaan dan Perhitungan')
Keluaran program
d =
32.0593
qmaks =
42.9100
polyfit
polyfit
polyfit
dipilih n = 1, maka akan dihasilkan persamaan garis lurus yaitu regresi linier.
Program Matlab dengan penggunaan fungsi polyfit untuk contoh yang sama
% Persamaan BET untuk proses adsorpsi
%
qmaks d C
%
q = ------------------------------%
(Cs - C)[ 1 + (d - 1)*C/Cs ]
%
Linierisasi
%
C
1
(d-1)C
%
--------- = ------- + ---------%
q(Cs - C)
qmaks d
qmaks d Cs
%
C
%
y = --------x = C
108
%
%
%
%
%
%
%
%
%
%
q(Cs - C)
1
(d - 1)
a0 = ------a1 = ---------qmaks d
qmaks d Cs
Penyelesaian dengan fungsi polyfit
--------------------------------------------------------------Surakarta, Oktober 2005
Jurusan Teknik Kimia, Fak. Teknik
Universitas Sebelas Maret
***************************************************************
% Data - data
C = [0.8 4.88 8.04 16.1 19.25];
q = [17.23 32.09 46.75 62.90 64.73];
adsorben
% Konsentrasi, mg/L
% Logam terjerap, mg/g
Cs=49.46;
% Linierisasi
y = C./(Cs-C)./q;
x = C;
% fungsi polyfit
p=polyfit(x,y,1);
% Parameter BET
d=p(1)*Cs/p(2)+1
qmaks=1/p(2)/d
Keluaran program
d =
32.0593
qmaks =
42.9100
-26,7
-4,4
6,4
18,4
31,8
40,3
p (torr)
10
20
40
60
p = exp a +
T
+
273
,
2
Penyelesaiannya
b
p = exp a +
T + 273,2
ln(p) = a +
b
T + 273,2
109
y=a+bx
untuk y = ln(p) dan x =
1
T + 273,2
% dalam C
% dalam Torr
% Linierisasi
y = log(p);
x = 1./(T+273.15);
kons = polyfit(x,y,1);
a = kons(2)
b = kons(1)
Keluaran program
a =
19.1905
b =
-4.7260e+003
Sehingga
4726,0
= exp19,1905
T
+ 273,2
110
15
12
= konstanta
Dari persamaan Torricelli terlihat hubungan laju alir dan volume cairan dalam
tangki adalah pangkat . Sehingga jika diplotkan log10(f) dan log10(V) akan
dihasilkan persamaan garis lurus.
111
Program Matlab
%
%
%
%
%
%
%
%
%
%
%
Resistansi Hidrolik
persamaan Torricelli adalah
f = r*V^(1/2)
f = laju alir cairan melalui katup (L/detik)
V = volum cairan di dalam tangki (L)
r = konstanta
--------------------------------------------------------------Surakarta, Oktober 2005
Jurusan Teknik Kimia, Fak. Teknik
Universitas Sebelas Maret
***************************************************************
% Data - data
V = [6 9 12 15];
waktu = [9 8 7 6 ];
laju_alir = 1./waktu;
%
p
m
r
% volume, L
% waktu pengisian, detik
% laju alir, 1/detik
Penggunaan polyfit
= polyfit(log10(V),log10(laju_alir),1);
= p(1);
= 10^p(2)
Keluaran program
waktu
Laju alir
112
....(5.15)
sehingga
SSE = (y data y pers )2
n
....(5.16)
n =1
i =1
n =1
( yi ao a1xi - a2 x i - .... - am x im ) = 0
....(5.17)
113
SSE = -2
a 1
SSE = -2
a 2
x i ( yi ao a1xi - a2 x i - .... - am x im ) = 0
.
.
.
.
SSE
= -2
a m
xi ( yi ao a1xi - a2 x i - .... - am x im ) = 0
a0 xi + a1 xi + a3 xi + .... + am xim + 2 = xi yi
2
.
.
.
.
2m = m y
a0 xim + a1 xim + 1 + a2 xim + 2 + .... + am xi
xi i
xi
x i2
x 3i
xi
x2
i
....
m
x i
x im +1
x i2
x 3i
x i4
x im + 2
....
....
x im a 0 yi
x im +1 a1 = x2i yi
xim + 2 a 2 x i yi
x i2 m
....
a
m
....
xmy
i i
114
Fungsi
polyfit
polyfit
koefisien-koefisien polinomial.
Selain fungsi
polyfit,
polyval
untuk
Tekanan
Temperatur
Tekanan
(OC)
(mmHg)
(OC)
(mmHg)
-36,7
15,4
60
-19,6
26,1
100
-11,5
10
42,2
200
-2,6
20
60,6
400
7,6
40
80,1
760
B
T + 273,15
dengan P adalah tekanan uap dalam mmHg dan T adalah temperatur dalam OC.
Parameter A dan B ditentukan dengan regresi.
Program pertama untuk penyelesaian persamaan empiris polinomial
% Persamaan Fit untuk Data Tekanan Uap Benzena
%
Hubungan tekanan sebagai fungsi temperatur
%
sebagai persamaan empiris polinomial sederhana
%
P = a0 + a1*T + a2*T^2 + a3*T^3 + ... + an*T^n
%
Penyelesaian dengan menggunakan fungsi polyfit
%
Persamaan polinomial tersebut kemudian dievaluasi
115
%
%
%
%
%
%
% Data - data
% tekanan, mmHg
vp = [ 1 5 10 20 40 60 100 200 400 760];
% temperatur, oC
T =[-36.7 -19.6 -11.5 -2.6 7.6 15.4 26.1 42.2 60.6 80.1];
% Penggunaan fungsi polyfit
m = 3;
%orde polinomial
p = polyfit(T,vp,m);
% Evaluasi persamaan polinomial
z = polyval(p,T);
% nilai vp dengan persamaan polinomial
% Cetak hasil
plot(T,z,'k-',T,vp,'ko','linewidth',2)
xlabel('T(C)')
ylabel('vp(mmHg)')
vp(mmHg)
Keluaran program
116
% Persamaan Fit untuk Data Tekanan Uap Benzena
%
Hubungan tekanan sebagai fungsi temperatur
%
sebagai persamaan empiris polinomial sederhana
%
P = a0 + a1*T + a2*T^2 + a3*T^3 + ... + an*T^n
%
Penyelesaian dengan menggunakan fungsi polyfit
%
Persamaan polinomial tersebut kemudian dievaluasi
%
dengan fungsi polyval
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
% Data - data
% tekanan, mmHg
vp = [ 1 5 10 20 40 60 100 200 400 760];
% temperatur, oC
T =[-36.7 -19.6 -11.5 -2.6 7.6 15.4 26.1 42.2 60.6 80.1];
% Linierisasi
y = log10(vp);
x = 1./(T+273.15);
% Penggunaan fungsi polyfit
p = polyfit(x,y,1);
% Evaluasi persamaan
z = 10.^(p(2)+(p(1)./(T+273.15)));
% Cetak hasil
plot(T,z,'k-',T,vp,'ko','linewidth',2)
xlabel('T(C)')
ylabel('vp(mmHg)')
Keluaran program
vp(mmHg)
117
Dari kedua grafik terlihat persamaan polinomial memberikan hasil yang lebih
baik.
0,0
0,1
0,2
0,3
0,4
-ra
0,0053
0,0052
0,0050
0,0045
0,0040
0,5
0,6
0,7
0,8
0,85
-ra
0,0033
0,0025
0,0018
0,00125
0,001
Program Matlab
%Regresi Polinomial
118
%
%
%
%
%
%
%
%
%
%
%
% Data-data
X=[0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.85];
ra=[0.0053 0.0052 0.005 0.0045 0.004 0.0033 ...
0.0025 0.0018 0.00125 0.001];
% Fungsi polyfit orde 3
p=polyfit(X,ra,3)
% Evaluasi persamaan polinomial
f=polyval(p,X);
% Plot data percobaan dan perhitungan
plot(X,ra,'ko',X,f,'k-','linewidth',2)
xlabel('x'); ylabel('y');
axis([0 0.85 0 0.006]);
title('Regresi Polinomial');
Keluaran program
119
....(5.18)
sehingga
SSE = (y data y pers )2
n
....(5.19)
n =1
SSE =
2
2
E i = (yi a 0 a1x1,i a 2 x 2,i .... a k x k ,i )
i =1
n =1
SSE
SSE
SSE
SSE
= 0,
= 0,
= 0, ....,
= 0.
a 0
a 1
a 2
a k
a0n + a1x1,i + a2x2,i + .... + akxk,i
= yi
x1,i
x 2, i
2
x1,i x1,i x 2,i
x 2,i x1,i x 22,i
x1,i
x
2 ,i
....
x k ,i
x k ,i x1,i
....
x k ,i x 2,i ....
x k,i a 0 yi
x1,i x k a 1 x1,i yi
x 2,i x k a 2 = x 2,i yi
x 2k ,i
....
a k
....
x y
m,i i
....(5.20)
120
Penyelesaian dapat dilakukan dengan metode-metode yang telah dijelaskan pada
Bab II.
unsur 1 (%)
unsur 2 (%)
x1
x2
7,1
19,2
31
45
11
Keluaran program
121
a =
0.8000
10.2429
1.2143
max_persen_error =
3.2193
5.2. INTERPOLASI
Interpolasi berguna untuk memperkirakan harga-harga antara dari datadata yang telah diketahui dengan tepat.
Interpolasi yang paling populer adalah polinomial interpolasi diferensiasi
terbagi Newton. Sebelum pembahasan interpolasi polinomial secara umum, akan
dibahas interpolasi linier (versi pertama) dan interpolasi kuadratik (versi kedua).
f(x1)
f(x)
f(x0)
x0
x1
122
f (x) f (x 0 ) = f (x1 ) f (x 0 )
x1 x 0
x x0
....(5.21)
INTERPOLASI KUADRATIK
Interpolasi kuadratik adalah interpolasi polinomial versi 2, yang
membutuhkan 3 titik yang diketahui.
Bentuk umum
f(x) = b0 + b1(x-x0) + b2(x-x0)(x-x1)
....(5.22)
dengan
b0 = f(x0)
b1 = f ( x 1 ) f ( x 0 )
( x1 x 0 )
f ( x 2 ) f ( x1 ) f ( x1 ) f ( x 0 )
( x 2 x1 )
( x1 x 0 )
b2 =
(x 2 x 0 )
INTERPOLASI POLINOMIAL
Bentuk umum interpolasi polinomial Newton
f(x) = b0 + b1(x-x0) + b2(x-x0)(x-x1) + .... + bn(x-x0)(x-x1)(x-x2)....(x-xn) ....(5.23)
dengan
b0 = f(x0)
b1 = f[x1,x0]
b2 = f[x2,x1,x0]
.
.
.
bn = f[xn,xn-1, ....,x1,x0]
Differensiasi terbagi pertama
123
f[xi,xj] = f ( x i ) f ( x j )
(x i x j )
.
.
.
Differensiasi terbagi ke-n
f[xn,xn-1,....,x1,x0] = f [ x n , x n 1 ,..., x1 ] f [ x n 1 , x n 2 ,..., x 0 ]
xn x0
sehingga
f(x) = f(x0)+(x-x0)f[x1,x0]+(x-x0)(x-x1)f[x2,x1,x0]+....+
(x-x0)(x-x1)(x-x2)....(x-xn)f[xn,xn-1,....,x1,x0]
....(5.24)
xi
f(xi)
x0
f(xi)
x1
f(x1)
II
III
f[x1,x0]
f[x2,x1,x0]
f[x2,x1]
2
3
x2
f(x2)
x3
f(x3)
f[x3,x2]
f[x3,x2,x1,x0]
f[x3,x2,x1]
124
(%)
(g/cm3)
0,997
10
1,011
20
1,024
30
1,035
40
1,045
50
1,053
60
1,06
70
1,064
80
1,065
90
1,061
100
1,044
Penyelesaian
% Densitas Asam Asetat pada 25 oC
%
Interpolasi Polinomial Newton
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
format long
% Data-data
x=[60 70 80 90 100];
rho = [1.06 1.064 1.065 1.061 1.044];
% konsentrasi, %
% Densitas, g/cm^3
125
disp(t)
disp '----------------------------'
fa=1;
p=0;
for j=1:n
p=p+frho(1,j)*fa;
t=['Nilai densitas dengan Polinomial Newton orde ' num2str(j-1)
' = ' num2str(p)];
disp(t)
fa=fa*(xi-x(j));
end
disp ' '
Keluaran program
rho =
1.0600
0.0004
-0.0000
-0.0000
-0.0000
1.0640
0.0001
-0.0000
-0.0000
1.0650
-0.0004
-0.0001
1.0610
-0.0017
1.0440
1.06
1.062
1.0624
1.0623
1.0625
Dalam matlab dapat digunakan fungsi interp1 untuk interpolasi satu dimensi.
YI = interp1 (X,Y,XI,metode) adalah interpolasi untuk menentukan nilai
YI jika X = X1, dengan menggunakan data-data hubungan Y terhadap X.
Metode yang dapat digunakan antara lain :
nearest
- interpolasi dengan tetangga terdekat
linier
- interpolasi linier
spline
- interpolasi spline kubik
cubic
- interpolasi kubik
126
k
[BTU/jam.ft.OF]
T (OF)
0,0057
32
0,0074
115
0,0099
212
0,0147
363
Program penyelesaian
% Konduktivitas termal aseton
%
Hubungan ln(k) dan ln(T) mendekati garis lurus
%
Tentukan k pada T = 300 oF
%
Diselesaikan dengan fungsi interp1
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
% Data-data
k = [0.0057 0.0074 0.0099 0.0147];
T = [32 115 212 363]+460;
% BTU/jam.ft.oF
% temperatur, R
Keluaran program
k_pd_300 =
0.0126
127
Polinomial orde ke-n telah dapat digunakan untuk interpolasi di antara n
+ 1 titik data. Misalnya untuk delapan titik dapat diturunkan suatu polinomial
orde ke tujuh yang sempurna. Kurva ini akan melalui semua titik-titik data yang
ada. Tetapi ada kasus di mana fungsi-fungsi ini dapat membawa hasil yang
keliru. Suatu pendekatan alternatif adalah menerapkan polinomial yang lebih
rendah terhadap sub kumpulan titik data. Polinomial penyambungan demikian
disebut fungsi spline
Kurva orde tiga fungsi spline dilaksanakan untuk menyambung setiap pasang
titik-titik data yang dinamakan
spline kubik
SPLINE LINIER
Sambungan paling sederhana adalah sebuah garis lurus. Spline orde
pertama untuk sekumpulan susunan titik data dapat didefinisikan sebagai
sekumpulan fungsi linier yang menghubungkan titik-titik :
f(x) = f(x0) + m0(x x0)
x0 # x # x1
x1 # x # x2
xn-1 # x # xn
f ( x i +1 ) f ( x i )
x i +1 x i
128
SPLINE KUADRATIK
a1x2 + b1x + c
interval 1
x0
i=0
a2x2 + b2x + c
interval 2
a3x2 + b3x + c
interval 3
x1
x2
i=1
i=2
Gambar 5.10. Konsep Dasar Spline Kuadratik
x3
i=3
B2
+ B3T
T
129
Dengan dalam poise, dan T dalam K. Bi adalah parameter dengan ketentuan
sebagai berikut:
Fraksi berat
B1
B2.10-3
B3.102
-19,52
3,913
2,112
0,1
-22,14
4,475
2,470
0,2
-25,16
5,157
2,859
0,3
-28,38
5,908
3,255
0,4
-31,52
6,678
3,634
0,5
-34,51
7,417
MDEA
3,972
O
Penyelesaian
% Viskositas larutan MDEA
%
B2
%
ln miu = B1 + ----- + B3*T
%
T
%
dengan miu = viskositas, poise
%
T
= temperatur, K
%
Bi = parameter viskositas sebagai fungsi konsentrasi, x
%
(data dalam bentuk tabel)
%
viskositas MDE pada konsentrasi x = 0,43 dan T = 45 oC
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
% Data parameter viskositas sebagai fungsi x
x = [0 0.1 0.2 0.3 0.4 0.5];
B1 = [-19.52 -22.14 -25.16 -28.38 -31.52 -34.51];
B2 = [3.913 4.475 5.157 5.908 6.678 7.417]*10^3;
B3 = [2.112 2.470 2.859 3.255 3.634 3.972]*10^-2;
% Data x dan T yang ingin dicari viskositasnya
x_i = 0.43;
% Konsentrasi
T_i = 45+273.15;
% Temperatur, K
% Konstanta B1, B2, dan B3 pada 0.43 ditentukan terlebih dahulu
% Interpolasi nilai B pada konsentrasi 0.43
B1_i = interp1(x,B1,x_i,'spline');
B2_i = interp1(x,B2,x_i,'spline');
B3_i = interp1(x,B3,x_i,'spline');
% Viskositas dihitung
viscA = exp(B1_i + B2_i/T_i + B3_i*T_i);
% Viscositas pada 45 C ditentukan terlebih dahulu
% Viskositas pada 45 C dihitung
visc = exp(B1 + B2./T_i + B3.*T_i)
% Interpolasi viskositas pada konsentrasi 0.43
130
viscB = interp1(x,visc,x_i,'spline');
% Cetak Hasil
disp(' ')
disp(['Viskositas MDEA pada x = ',num2str(x_i), ' dan T =
',num2str(T_i)])
disp(['jika parameter viskositas pada x = ',num2str(x_i), '
dihitung dulu'])
viskositas = viscA
disp(['jika viskositas pada T = ',num2str(T_i),' dihitung dulu'])
viskositas = viscB
Keluaran program
Viskositas MDEA pada x = 0.43 dan T = 318.15
jika parameter viskositas pada x = 0.43 dihitung dulu
viskositas =
3.2302
jika viskositas pada T = 318.15 dihitung dulu
viskositas =
3.2093
ZI,
dari fungsi Z pada titik XI dan YI dari matriks X dan Y. Metode yang ditampilkan
pada interpolasi 1 dimensi dapat juga digunakan pada interpolasi 2 dimensi.
50
55
1258,2
1282,2
90
1328,7
1378,1
p (psia)
700
-
131
95
1328,4
1377,8
Program penyelesaian
% Beda entalpi uap yang didinginkan
%
uap mula-mula adalah 640 oF dan 92 psia
%
menjadi 480 oF dan 52 psia
%
Digunakan interpolasi 2 dimensi
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
% Kondisi awal
To = [600 700];
po = [90 95];
Ho = [1328.7 1378.1
1328.4 1377.8];
H1 = interp2(To,po,Ho,640,92)
% Kondisi akhir
Tn = [400 500];
pn = [50 55];
Hn = [1258.7 1282.6
1258.2 1282.2];
H2 = interp2(Tn,pn,Hn,480,52)
% Perhitungan beda entalpi
del_H = H2 - H1
Keluaran program
H1 =
1.3483e+003
H2 =
1.2777e+003
del_H =
-70.6880
% Temperatur, oF
% Tekanan, psia
% Entalphi, BTU/lb
% Interpolasi 2 dimensi
% Temperatur, oF
% Tekanan, psia
% Entalphi, BTU/lb
% Interpolasi 2 dimensi
OPTIMASI
Optimasi adalah suatu proses untuk mencari kondisi optimum dalam arti
yang paling menguntungkan. Optimasi bisa berupa proses mencari nilai
maksimum (maksimasi) atau proses mencari nilai minimum (minimasi).
Variabel yang dimaksimumkan atau diminimimkan disebut objective
function, sedangkan variabel yang dicari nilainya sehingga objective function
menjadi maksimum atau minimum disebut design variabel.
Proses optimasi secara analitis dilakukan dengan cara melihat turunan
(derivatif) objective function terhadap design variabel. Dalam prakteknya, proses
optimasi sangat sulit diselesaikan dengan cara analitis (kalkulus). Jika demikian
maka dapat dilakukan proses optimasi dengan cara direct search, yaitu
melakukan perhitungan harga objective function untuk berbagai nilai design
variabel.
134
dengan x adalah vektor kolom n komponen yang ingin ditentukan nilainya.
Konstanta yang diberikan sistem disediakan oleh vektor kolom b dengan m
komponen, matriks A ukuran m x n, dan vektor kolom c dengan n komponen.
Seluruh
persamaan
dan
fungsi
berada
dalam
bentuk
linier.
Tujuan
keseluruhannya adalah minimasi suatu fungsi linier cTx yang disebut sebagai
objective function. Penyelesaian secara detail tidak dijelaskan disini. Tetapi
program berikut dapat digunakan untuk menyelesaikan persoalan pemograman
linier.
function [xsol,basic]=proglin(A,b,c,tol)
%
Program ini bertujuan mencari nilai minimum
%
untuk pemograman linier c'x
%
subyek Ax = b
%
%
A adalah matriks koefisien untuk batasan-batasan
%
b adalah vektor kolom sisi kanan
%
c adalah vektor baris koefisien biaya
%
xsol adalah vektor hasil penyelesaian
%
basic adalah urutan variabel dasar (basic)
%
%
%
135
end
if basiscount==m
for k=basic
B=[B A(:,k)]; cb=[cb c(k)];
end
primalsol=b'/B'; xsol=primalsol;
break
end
iter=iter+1;
end
objective=c*x
Minimasi - keuntungan = - 3A 4B
Subyek 4A + 2B # 80
2A + 5B # 120
A, B 0
136
b=[80 120]';
% Perhitungan
[xsol,ind]=proglin(a,b,c,1e-10);
% Cetak hasil
i=1;fprintf('\Penyelesaian adalah : ');
for j=ind
fprintf('\nx(%1.0f)=%8.4f\',j,xsol(i));
i=i+1;
end
fprintf('\nNilai variabel lain adalah nol\n')
Keluaran program
objective =
-110.0000
x(1)= 10.0000
x(2)= 20.0000
Nilai variabel lain adalah nol
XA baru = XP lama
XB tetap
XP baru = XQ lama
XQ dicari
XA tetap
137
XB baru = XQ lama
XP dicari
XQ baru = XP lama
YP > YQ maka
XA baru = XP lama
XB tetap
XP baru = XQ lama
XQ dicari
if p(1)<p(2)
a=p(1); b=p(2);
else
a=p(2); b=p(1);
end
g=(-1+sqrt(5))/2;
r=b-a;
iter=0;
while r>tol
x=[a+(1-g)*r a+g*r];
y=feval(func,x);
if y(1)<y(2)
b=x(2);
else
a=x(1);
end
r=b-a;
iter=iter+1;
end
iter
f=feval(func,a);
138
l = 0,618
l
l
XA
XP
XQ
XB
XA
XP
XQ
XB
XA
XP XQ
XB
F xo
(4W + 0,79F)
a.
Buktikan bahwa: x
b.
c.
139
xo = 0,002 lb as. benzoat/ lb air
F = 10.000 lb air/jam
W = lb benzena / jam
Phase benzena
y lb as. benzoat/lb
benzena
Phase air
x lb as.benzoat / (lb air + lb benzena larut)
Gambar 6.2. Proses Ekstraksi Asam Benzoat dengan Benzena
Penyelesaian
a. Fraksi massa Asam Benzoat pada fase air
Neraca massa komponen Asam Benzoat
input output = akumulasi
As. Benzoat umpan As. Benzoat fase benzena As. Benzoat fase air = 0
As. Benzoat umpan (jumlah benzena x fraksi As. Benzoat pada fase
benzena) (jumlah benzena dan air x fraksi As. Benzoat pada fase
air) = 0
F x0 (W 0,07 F) y (F + 0,07 F) x = 0
F x0 (W 0,07 F) 4 x (1,07 F) x = 0
F x0 4 W x 0,79 F x = 0
x=
Fxo
(4W + 0,79F)
b. Keuntungan (P)
Keuntungan ( P ) = harga jual As. Benzoat fase benzena harga beli benzena
= As. Benzoat fase benzena x $ 0,4 Benzena x $0,001
140
= (x0 F x F) $ 0,4 (W) $0,001
= 0,4F(xo x) 0,001W
c. Jumlah benzena (W) optimal untuk mendapatkan keuntungan (P) maksimal
P
= 0,4F(xo x) 0,001W
= 0,4F(xo
Fxo
) 0,001W
(4W + 0,79F)
= 0,4.10000 (0,002 -
=8-
80000
(4W + 7900)
10000.0,002
(4W + 0,79.10000)
) 0,001.W
- 0,001W
function fw=contoh62(w)
% Ekstraksi asam benzoat dari limbah industri
%
Program bertujuan menentukan jumlah benzena yang ditambahkan
%
untuk pengambilan asam benzoat dari limbah industri
%
w adalah benzena yang ditambahkan
%
%
%
% Data-data
xo=0.002;
F=10000;
141
Matlab menyediakan fungsi untuk menentukan nilai minimasi satu variabel yaitu
fungsi fminbnd. Cara memanggil fungsi ini adalah
x = fminbnd(fun,x1,x2)
dari fungsi
fun
pada interval
dimulai pada
x0
x
x
xA
, dengan CA=CAO(1 - xA).
kC AO (1 x A )
a.
b.
CB
0,1.x A (1 x A )
=
ts
x A + 2,25(1 x A )
c.
d.
Penyelesaian
a.
dC A
= k.C 2A
dt
Ca
dC A
2
Cao C A
= k dt
0
1 - 1 = kt
C A C Ao
142
1 = kt
1
C Ao (1 x A ) C Ao
xA
= kt
C Ao (1 x A )
tr =
b.
xA
kC AO (1 x A )
CB
=
ts
C Ao x A
1
kC Ao
xA
1 xA
+ tp
= C Ao x A kC Ao (1 x A )
x A + tp kC Ao (1 x A )
=
c.
0,1.x A (1 x A )
x A + 2,25(1 x A )
Program matlab
% gmol/L
% L/(gmol.menit)
% menit
% Interval konversi
xo=0.01;
% konversi minimal yang mungkin terjadi
xn=0.99;
% konversi maksimal yang mungkin terjadi
% Perhitungan konversi optimal
x_opt=fminbnd('F63',xo,xn)
Fungsi matlab
function Cb_ts_min=F63(x)
%
Program ini berisi fungsi untuk menghitung
143
%
%
%
%
global k Cao tp
tr = x/(k*Cao*(1-x));
ts = tp + tr;
Cb_ts = Cao*x/ts;
Cb_ts_min = 1/Cb_ts;
%
%
%
%
Waktu reaksi
Waktu siklus
Maksimum B terbentuk per waktu siklus
Program minimasi
Keluaran program
x_opt =
0.6000
tr =
xA
kC AO (1 x A )
= 30 menit
Cara klasik untuk mendapatkan suatu nilai yang optimal adalah dengan mencoba
semua titik yang mungkin. Cara ini lebih lama tetapi cukup menarik untuk
optimasi proses yang sederhana.
144
% Produksi optimal
%
Perhitungan jumlah unit barang yang diproduksi setiap tahun
%
Untuk mendapatkan keuntungan maksimal
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
n = 5000;
%Jumlah simulasi random
level = [25:50];
%vektor tingkat produksi awal
cost = 30000+2000*level;
%biaya
for k = 1:26
cum_profit = 0;
for m = 1:n
demand = floor(rand*(50-25)+26);
if demand >= level(k)
income = 4000*level(k);
else
income = 4000*demand+1000*(level(k)-demand);
end
profit = income-cost(k);
cum_profit = cum_profit+profit;
end
expected_profit = cum_profit/n;
p(k,1) = level(k);
p(k,2) = expected_profit;
end
plot(p(:,1),p(:,2),'o',p(:,1),p(:,2),'-')
xlabel('Jumlah unit'),ylabel('Profit ($)')
yaitu antara 25 50. loop m menunjukkan jumlah simulasi secara acak. Variabel
demand menunjukkan jumlah permintaan untuk produksi. Jika permintaan lebih
besar daripada jumlah produksi, maka seluruh barang akan dijual dan incomenya
adalah $ 4.000 kali jumlah barang. Jika income lebih sedikit daripada produksi
barang, maka income adalah $ 4.000 kali jumlah barang terjual ditambah $ 1.000
kali jumlah barang yang tidak terjual. Pada setiap permintaan yang acak, maka
profit adalah selisih pdapatn dan biaya. profit_kum adalah profit kumulatif
Profit ($)
145
optimal pada fungsi fun. fun menerima input x, dan keluar adalah F yang
dievaluasi pada x. x0dapat berupa scalar, vektor, atau matriks.
Contoh 5.4 tentang hubungan tekanan uap benzena dan temperatur, dapat
pula digunakan persamaan Antoine
log(P) = A
B
T+C
146
Program
% Persamaan Fit untuk Data Tekanan Uap Benzena
%
Regresi non linier
%
B
%
log (P) = A - ------%
T + C
%
Penyelesaian dengan fungsi fminsearch
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
global vp T
% Data-data
vp = [ 1 5 10 20 40 60 100 200 400 760];
T = [-36.7 -19.6 -11.5 -2.6 7.6 15.4 26.1 42.2 60.6 80.1];
%Tebakan untuk nilai parameter
Ko = [10 2000 273] ;
% Perhitungan optimasi
K = fminsearch('F65',Ko)
% Evaluasi parameter
z = 10.^(K(1)-(K(2)./(T+K(3))));
% Menampilkan hasil
plot(T,z,'k-',T,vp,'ko','linewidth',2)
set(gca,'fontsize',14)
xlabel('T(C)')
ylabel('vp(mmHg)')
Fungsi terkait
function ff=F65(K)
%
Program ini bertujuan mencari nilai minimum
%
program non linier
%
Persamaan Antoine
%
K(1) = A
%
K(2) = B
%
K(3) = C
%
%
Nama File : F65.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global vp T
% Perhitungan minimasi SSE
f = log10(vp)-K(1)+K(2)./(T+K(3));
ff=sum(f.*f);
Keluaran program
K =
5.7673
677.0928
Artinya A = 5,7673
153.8852
147
B = 677,0928
vp(mmHg)
C = 153,8852
PE
PEA
PH
(mol/kgkat.detik)
(atm)
(atm)
(atm)
1,04
3,13
5,21
3,82
4,19
2,391
0,5
3,867
0,5
0,5
148
2,199
0,5
0,75
0,5
kPE PH
1 + K A PEA + K E PE
Penyelesaian
% Hidrogenasi Etilen menjadi Etana
%
Reaksi hidrogenasi etilen menjadi etana
%
H2 + C2H4 --> C2H6
%
Regresi Non Linear
%
kPePh
%
ra = -----------------%
1 + KaPea + KePe
%
Ingin menghitung nilai k, Ka, dan Ke
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
clear all
clc
global ra Pe Ph Pea n
% Data-data Percobaan
ra=[1.04 3.13 5.21 3.82 4.19 2.391
%
%
Pe=[1 1 1 3 5 0.5 0.5 0.5 0.5]; %
Pea=[1 1 1 1 1 1 0.5 3 5];
%
Ph=[1 3 5 3 3 3 5 3 1];
%
% Jumlah data
n=9;
% Tebakan awal k, Ka, Ke masing-masing 1
xo=[1 1 1];
%Minimasi di sub program F66
disp 'Harga k, Ka, dan Ke hasil minimasi'
fminsearch('F66',xo)
149
%
%
%
%
%
%
%
% Data-data
global ra Pe Ph Pea n
% Minimasi SSE
f=0;
for i=1:n
f=f+(ra(i)-(x(1)*Pe(i)*Ph(i))/(1+x(3)*Pea(i)+x(2)*Pe(i)))^2;
end
Hasil program
>> Harga k, Ka, dan Ke hasil minimasi
ans =
3.3479
2.2111
0.0429
3,3479PE PH
1 + 2,2111PEA + 0,0429PE
Gambar 6.5 menunjukkan cross section dari kanal saluran air. Analisa awal
menunjukkan bahwa luas penampang saluran harus 100 ft3 untuk dapat
menampung kecepatan aliran air. Untuk meminimasi biaya digunakan kanal
lurus dengan cara meminimasi panjang perimeter kanal. Tentukan d, b, dan
untuk meminimasi panjang kanal ini.
d
b
150
L = b + 2d
sin
Luas penampang saluran
2
= 100
A = db + d
tan
b = 100 d
d
tan
100
d
2d
+
d
tan sin
Program Matlab
% Optimasi Kanal Saluran Air
%
menentukan d dan sudut tetha
%
100
d
2d
%
L = ---- - ---------- + ---------%
d
tan(tetha)
sin(tetha)
%
Ingin menghitung nilai k, Ka, dan Ke
%
---------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
*****************************************************************
%
clear all
clc
% Tebakan awal
d_tebakan = 20;
teta_tebakan = 1;
% ft
% Perhitungan optimasi
x =fminsearch('F67',[d_tebakan teta_tebakan]);
% Hasil perhitungan
d = x(1)
teta = x(2)*180/pi
% konversi ke derajat
Program perhitungan L
function L = F67(x)
%
Program ini berisi fungsi untuk menghitung L
%
sebagai fungsi d dan tetha
%
%
Nama File : F67.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------d = x(1);
teta = x(2);
% Persamaan yang diminimasi
L = 100./d-d./tan(teta)+2*d./sin(teta);
151
Keluaran program
>>
d =
7.5983
teta =
60.0004
clear all
global xDataz yDataz
% Header Program
disp('Demonstrasi regresi nonlinear')
disp(' y = alpha*exp(beta*x)')
% Data-data
xDataz = [ 0 2 3 5 7 8 9 ]';
yDataz = [ 10 8 5 3 2 1 1 ]';
% Tebakan
guess = [15 -0.5];
% Perhitungan
param = fminsearch('F68', guess);
% Cetak hasil perhitungan
alpha = param(1);
beta = param(2);
disp('Koefisien adalah: ')
disp([ 'alpha = ', num2str(alpha) ])
disp([ 'beta = ', num2str(beta) ])
% Evaluasi hasil
yfit = alpha.*exp(beta.*xDataz);
% Plot dan membandingkan hasil
plot(xDataz, yDataz, 'o', xDataz, yfit, '-', 'linewidth',2)
set(gca,'fontsize',16)
152
global xDataz yDataz
% Mengubah variabel
alpha = param(1);
beta = param(2);
% Perhitungan SSE
error = yDataz - alpha.*exp(beta.*xDataz);
sse = sum(error.^2);
Keluaran Program
y = k.x2
153
S1
S3
S2
1
x0
x1
F
x2
F
x3=xR
S1
S2
S3
y1
y2
y3
x1 =
x2 =
F + F 2 + 4F.S1 .k.x o
2S1 .k
F + F 2 + 4F.S 2 .k.x 1
xR = x3 =
2S 2 .k
F + F 2 + 4F.S3 .k.x 2
2S3 .k
S3 = S S1 S2
Data-data yang diketahui adalah F = 100 kg A/mnt; S = 180 kg B/mnt; xo = 0,2
kg C/ kg A; dan k = 20
Program Penyelesaian
% Ekstraktor Tiga Tingkat
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
global Stot xo F S3 k
% Data-data
Stot=180;
xo=0.2;
F=100;
k=20;
%
%
%
%
154
% Hasil perhitungan
S1=S(1)
S2=S(2)
S3=Stot-S1-S2
xR_min=xR
Program terkait
function xR=F69(S)
%
Program ini berisi fungsi untuk menghitung XR
%
fungsi konstanta-konstanta yang ingin dicari
%
S(1) & S(2) ---> S(3) akan terhitung
%
%
Nama File : F69.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global Stot xo F S3 k
S3=Stot-S(1)-S(2);
x1= (-F+sqrt(F^2+4*F*S(1)*k*xo))/(2*S(1)*k);
x2= (-F+sqrt(F^2+4*F*S(2)*k*x1))/(2*S(2)*k);
xR= (-F+sqrt(F^2+4*F*S3*k*x2))/(2*S3*k);
Keluaran program
S1 =
47.3733
S2 =
59.8036
S3 =
72.8231
xR_min =
0.0382
Pada suhu 350 K, campuran biner bahan A dan B, fasa uapnya mengikuti
hukum-hukum gas ideal sedangkan fasa cairnya non-ideal dengan koefisien
aktivitas mengikuti persamaan :
A = exp ( .xB2 )
dan
B = exp ( .xA2 )
Pada 350 K, tekanan uap murni A dan B masing-masing 80 cmHg dan 60 cmHg.
Data tekanan total kesetimbangan pada 350 K pd berbagai fraksi mol A adalah :
xA
Pt
(cmHg)
0,1
0,2
0,3
0,4
0,6
0,7
0,8
0,9
65,8
70,3
73,9
76,5
79,9
80,9
81,3
81
155
Berdasarkan data-data tersebut, perkirakan nilai dengan membandingkan Pt
data dengan Pt hasil perhitungan sampai diperoleh SSE minimum.
Program terkait
function fP=F610(B)
%
Program ini berisi fungsi untuk menghitung SSE
%
%
SSE = sigma [ P data - P hit ]^2
%
%
Nama File : F610.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global xa Ptdat Pao Pbo
Pt=xa*Pao.*exp(B*(1-xa).^2)+(1-xa)*Pbo.*exp(B*xa.^2);
fP=sum((Ptdat-Pt).^2);
Keluaran program
beta =
0.4584
Contoh 6.11.
156
Reaksi fasa gas bolak-balik eksotermis SO2 + O2 SO3 dijalankan dalam
reaktor tabung adiabatik dengan bantuan katalisator padat berbentuk butir-butir.
Aliran gas dianggap plug-flow dan tekanan system dianggap tetap Pt = 1 atm.
Proses ditunjukkan pada gambar 6.8.
Kecepatan reaksi mengikuti persamaan :
PSO3
k 1 .PO 2 .PSO 2 .1
0
,
5
P
.
P
.
K
gmol SO 2
SO 2
O2
P
=
r
2
jam . kg katalisator 22,414. 1 + K 2 .PSO 2 + K 3 .PSO 3
dengan :
5473
k 1 = exp12,16
8619
K 2 = exp 9,953 +
52956
K 3 = exp 71,745
11300
K P = exp 10,68 +
PSO 2 =
(1 x)
.Pt
(12,82 0,5.x)
PO 2 =
(1 x)
.Pt
(12,82 0,5.x)
PSO3 =
x
.Pt
(12,82 0,5.x)
xF
Tin
T = TF + 241,7.(x xF)
xin
Tin
xout
Tout
W kg katalis
F gmol N2 / j
SO2 = 7,8%
O2 = 10,8%
N2 = 81,4%
R gmol N2 / j
157
out
R 1
W = 137,3.1 + . . dx
F x in r
R
x F + .x out
R
x in =
R
1+
F
Diketahui data-data sebagai berikut, xF = 0; xout = 0,64; TF = 660 K
Program penyelesaian
% Reaktor Adiabatis untuk Reaksi Fase Gas Bolak-balik Eksotermis
%
%
SO2 + O2 ---> SO3
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
global xf xout Tf
% Data-data
xf=0;
% Konversi umpan
xout=0.64; % Konversi keluar reaktor
Tf=660;
% Temperatur umpan, K
% Perhitungan
[RperF,Wmin]=fminsearch(@F611,4)
Program terkait
158
function W=F611(RperF)
%
Program ini berisi fungsi untuk menghitung berat katalis
%
%
Nama File : F611.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global xf xout Tf
% Perhitungan konversi masuk reaktor
xin=(xf+RperF*xout)/(1+RperF);
% Distribusi nilai konversi dari konversi masuk sampai keluar
x=linspace(xin,xout,100);
% Perhitungan temperatur reaktor
T=Tf+241.7*(x-xf);
% Perhitungan tekanan parsial
Pso2=(1-x)./(12.82-0.5*x);
Po2=(1.3846-0.5*x)./(12.82-0.5*x);
Pso3=x./(12.82-0.5*x);
% Perhitungan kecepatan reaksi
k1=exp(12.16-5473./T);
K2=exp(-9.953+8619./T);
K3=exp(-71.745+52596./T);
Kp=exp(-10.68+11300./T);
Ra=k1.*Po2.*Pso2.*(1-Pso3./(Pso2.*sqrt(Po2).*Kp))./...
(22.414*(1+K2.*Pso2+K3.*Pso3).^2);
Ra1=1./Ra;
% Perhitungan berat katalis
area=trapz(x,Ra1);
W=137.3*(1+RperF)*area;
Keluaran program
RperF =
2.2554
Wmin =
4.4758e+003
Sehingga pada R/F = 2,25 akan diperoleh berat katalis minimum yaitu 4475,8 kg.
Reaktor batch beroperasi secara adiabatis untuk reaksi fasa cair order 2 : A B.
Perubahan entalpi reaksi HR, volume reaktor VR, dan kapasitas panas larutan Cp
dianggap tetap. Waktu bongkar pasang tp. Umpan masuk pada suhu TF dan
konsentrasi A CAo. Reaksi dihentikan pada konversi xR. Konversi A mula-mula
xRo. Ingin dicari konversi A yang memberikan kecepatan produksi B tiap waktu
maksimum.
Persamaan-persamaan yang dipakai adalah :
159
tR =
1
C Ao
x out
x in
1
dx A
k.(1 x A ) 2
E
k = A. exp
RT
T = TF B = VR
C Ao H R
xA
.Cp
C Ao .x R
(t R + t p )
%
%
%
%
%
%
%
%
%
%
Temperatur umpan, K
Konsentrasi umpan, gmol/L
Beda entalphi reaksi, kkal/gmol
Densitas larutan, kg/L
Kapasitas panas larutan, kcal/kg.K
Konstanta Arrhenius, L/gmol.menit
Konstanta Arrhenius, kkal/gmol
Konstanta gas ideal, kkal/gmol.K
Konversi mula-mula
Volume reaktor, L
% Perhitungan konversi
konversi=fminsearch(@F612,0.5)
Program terkait
function fxr=soal_15b(xr)
%
Program ini berisi fungsi untuk menghitung konversi
%
%
Nama File : F612.m
%
Surakarta, Oktober 2005
%
---------------------------------------------------------------
160
global Tf Cao dHr rho Cp A Ea R xro Vr tp
tr=quad(@F612F,xro,xr);
fxr=(tr+tp)/(Vr*Cao*xr);
Program lain
function fx=F612F(x)
%
Program ini berisi fungsi untuk menghitung integral
%
waktu reaksi
%
Nama File : F612F.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global Tf Cao dHr rho Cp A Ea R xro Vr tp
T=Tf-Cao*dHr*x/(rho*Cp);
k=A*exp(-Ea./(R*T));
fx=1./(Cao*k.*(1-x).^2);
Keluaran program
konversi =
0.6267
t=
xN
dx
k[(1 x) K.x ]
xo
E
k = A. exp
RT
K = . exp
T
T = To +
.C Ao
(x - x o )
.Cp
Data-data yang diketahui adalah : A = 25; E = 10.000; R = 1,987; = 1750; = 5.000; = 1100; Cp = 1,2; CAo = 1; xo = 0; xN = 0,3.
161
Program penyelesaian
% Reaktor Batch Adiabatis Reaksi Fase Cair Bolak-balik
%
A <---> B
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
global alpa rho Cp Cao xo xn beta A Ea R lamda
% Data-data
alpa=1750;
rho=1100;
Cp=1.2;
Cao=1;
xo=0;
xn=0.3;
beta=-5e3;
A=25;
Ea=1e4;
R=1.987;
lamda=-8e4;
% Perhitungan temperatur awal optimum
[To,waktu]=fminsearch(@F613,700)
Program terkait
function ft=F613(To)
%
Program ini berisi fungsi untuk menghitung waktu reaksi
%
%
Nama File : F613.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global alpa rho Cp Cao xo xn beta A Ea R lamda
% Distribusi konversi di reaktor
x=linspace(xo,xn,100);
% Perhitungan temperatur reaktor
T=To+lamda*Cao*(x-xo)/rho/Cp;
% Perhitungan waktu reaksi
k=A*exp(-Ea./(R*T));
K=alpa*exp(beta./T) ;
y=1./(k.*((1-x)-K.*x));
ft=trapz(x,y);
Keluaran program
To =
736.3329
waktu =
25.6426
162
Reaksi endotermis order 1 fasa gas : A B + C, dijalankan dengan katalisator
padat dalam 2 reaktor adiabatis yang disusun seri seperti terlihat pada gambar
6.9. Suhu umpan masuk ke reaktor-reaktor dijaga Tin dengan preheater dan interheater. Umpan reaktor 1 berupa A murni sebanyak FAo. Tekanan reaktor-reaktor
dianggap tetap P. Aliran gas dianggap plug flow. Carilah konversi antara yang
memberikan berat katalisator W minimum.
Berat katalisator dinyatakan dengan persamaan,
x out
Wi = FAo .
x in
dx
, ( i =1 , 2 )
rA
rA =
k.P 1 x
R .T 1 + x
T = TR
x0, Tin
x1, T1
x1, Tin
x2, T2
x0, T0
W2
W1
HE
HE
Program penyelesaian
% Dua Reaktor Adiabatis Tersusun Seri
%
Reaksi endotermis order 1 fasa gas
%
A ---> B + C
%
dengan katalisator padat dalam 2 reaktor adiabatis
%
yang disusun seri
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
global xin xout Tr Tin Cpa Cpb Cpc dHr R P Fao
% Data-data
xin=0;
xout=0.6;
163
Tr=298;
Tin=1200;
Cpa=4;
Cpb=5.5;
Cpc=7;
dHr=2000;
P=2;
R=0.082;
Fao=10;
%
%
%
%
%
%
%
%
%
Temperatur referensi, K
Temperatur masuk reaktor, K
Kapasitas panas A, kal/gmol.K
Kapasitas panas B, kal/gmol.K
Kapasitas panas C, kal/gmol.K
Beda entalphi reaktsi, kal/gmol
Tekanan reaktor, atm
Konstanta gas ideal, L.atm/gmol.K
Kecepatan umpan masuk,gmol/detik
Program terkait
function fw=F614(x1)
%
Program ini berisi fungsi untuk menghitung berat katalis
%
%
Nama File : F614.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global xin xout Tr Tin Cpa Cpb Cpc dHr R P Fao
% Distribusi nilai x dari xin sampai x1
x=linspace(xin,x1,40);
% Perhitungan temperatur keluar reaktor 1
T=Tr-((Tr-Tin)*((1-xin)*Cpa+xin*(Cpb+Cpc))+(x-xin)*dHr)./...
((1-x)*Cpa+x*(Cpb+Cpc));
% Perhitungan berat katalis reaktor 1
k=2.5e16*exp(-33000./T);
ya=Fao./(k*P.*(1-x)./(R*T.*(1+x)));
wa=trapz(x,ya);
% Distribusi nilai x dari x1 sampai xout
x=linspace(x1,xout,40);
% Perhitungan temperatur keluar reaktor 2
T=Tr-((Tr-Tin)*((1-x1)*Cpa+x1*(Cpb+Cpc))+(x-x1)*dHr)./...
((1-x)*Cpa+x*(Cpb+Cpc));
% Perhitungan berat katalis reaktor 2
k=2.5e16*exp(-33000./T);
yb=Fao./(k*P.*(1-x)./(R*T.*(1+x)));
wb=trapz(x,yb);
% Berat katalis total
fw=wa+wb;
Keluaran program
x1 =
0.2506
W =
560.4560
Sehingga konversi dari reaktor 1 yang optimal adalah 0,2506 dengan berat
katalis total minimal 560,456 kg.
164
Konversi A
(menit)
Suhu
(K)
340
30
0,11
344
60
0,28
350
90
0,53
359
120
0,78
369
150
0,90
373
t=
1
C Ao
dx
k.(1 x)
E
k = A. exp
RT
T = To + CAo..x
Tentukan nilai dengan cara regresi linear kemudian tentukan nilai A dan B
dengan mencari SSE minimum, SSE dinyatakan dengan persamaan berikut,
Penyelesaian
% Reaktor Batch untuk Reaksi Gas Cair
%
Reaksi fasa cair eksotermis
%
A ---> B
%
reaktor tangki berpengaduk yang bekerja
%
secara batch dan adiabatis
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
% Tebakan konstanta
165
Ao=[2e18 1.8e4];
% Penentuan konstanta
[A]=fminsearch(@F615,Ao);
% Menampilkan hasil perhitungan
Ar=A(1)
B=A(2)
Program terkait
function ft=F615(A)
%
Program ini berisi fungsi untuk menghitung SEE
%
%
Nama File : F615.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------% Data-data
Cao=1.6; To=340;
tdat=[0 30 60 90 120 150];
Tdt=[340 344 350 359 369 373];
xdt=[0 0.11 0.28 0.53 0.78 0.9];
% Penentuan Beta
p=polyfit(xdt,Tdt,1);
beta=p(1)/Cao;
% Perhitungan waktu
for i=1:6
x=linspace(0,xdt(i),100);
T=To+beta*Cao*x;
k=A(1)*exp(-A(2)./T);
y=1./(Cao*k.*(1-x).^2);
wkt(i)=trapz(x,y);
end
% Perhitungan SEE
ft=sum((tdat-wkt).^2);
Keluaran program
Ar =
9.9683e+017
B =
1.6204e+004
PERSAMAAN
DIFFERENSIAL ORDINER
170
maka persoalan tersebut disebut masalah harga batas (boundary value problem).
Pada bagian awal bab ini akan dibahas masalah harga awal terlebih dahulu yang
akan dilanjutkan penyelesaian untuk masalah harga batas pada bagian akhir bab.
Ada dua kategori utama metode numerik dalam penyelesaian persamaan
differensial ordiner, yaitu metode satu langkah dan metode multi langkah. Untuk
metode satu langkah akan dibahas metode Euler dan metode Runge-Kutta.
Sedangkan metode multi langkah akan dibahas metode predictor-korektor.
slope
xi
xi+1
h
(7.1)
171
yn+1 = yn + f(tn,yn) h
(7.2)
Contoh 7.1.
Selesaikan dengan Metode Euler
dy/dt = y 20
Program Matlab
% Demonstrasi penggunaan metode euler
%
Program ini bertujuan menyelesaikan
%
dy
%
---- = y - 20
y(0) = 100
%
dt
%
Tentukan y sebagai fungsi t
%
Penyelesaian dengan Metode Euler
%
Dicoba untuk berbagai besaran step yang berbeda
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
[t y] = euler('F71',[0 5],100,0.5); % step = 0.5
plot(t,y,'ko','LineWidth',2);
axis([0 5 0 12000])
pause(3)
hold on;
[t y] = euler('F71',[0 5],100,0.2); % step = 0.2
plot(t,y,'kx','LineWidth',2);
pause(3)
[t y] = euler('F71',[0 5],100,0.1); % step = 0.1
plot(t,y,'k+','LineWidth',2);
pause(3)
xlabel('time','FontSize',14);ylabel('y','FontSize',14);
plot(t,80*exp(t)+20,'k-')
hold off;
172
Fungsi terkait
function v=F71(t,y)
%
Program ini berisi fungsi untuk menghitung v
%
dy
%
v = ---- = y - 20
%
dt
%
Nama File : F71.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------v= y - 20;
Keluaran program
(7.3)
f(tn, yn, h) disebut fungsi inkremen yang dapat diinterpretasikan sebagai sebuah
slope rata-rata sepanjang interval. Fungsi inkremen dapat ditulis sebagai
173
f(tn, yn, h) = a1k1 + a2k2 + + ankn
Ada beberapa jenis metode Runge-Kutta yang dapat digunakan tergantung dari
pengembangan jumlah suku-suku yang berbeda pada fungsi inkremen. Metode
Runge-Kutta yang paling popular adalah Metode Runge-Kutta klasikal yang
menggunakan orde 4 ( n = 4 ).
yn+1 = yn + (k1 + 2k2 + 2k3 + k4)/6
dengan
(7.4)
k1 = h(tn, yn)
k2 = h(tn + h/2, yn + k1/2)
k3 = h(tn + h/2, yn + k2/2)
k4 = h(tn + h, yn + k3)
174
Metode Euler dan metode Runge-Kutta menggunakan harga variabel tak
bebas yn pada sebuah titik berikutnya tn+1. Pendekatan alternatif lain adalah
dengan mendasarkan pengertian bahwa sekali perhitungan dimulai informasi
yang berharga dari titik-titik terdahulu berada dalam pengaturan. Pendekatan ini
disebut metode multi langkah. Metode prediktor-korektor adalah metode multi
langkah yang menggunakan beberapa persamaan untuk memperkirakan harga
variabel tak bebas ke-n+1. Hasilnya kemudian koreksi dengan persamaan
korektor. Sebagai contoh dapat digunakan metode Adam-Bashforth-Moulton.
yn+1 = yn + h(55yn 59yn-1 + 37yn-2 9yn-3)/24
(7.5)
y n+1 = f(tn+1,yn+1)
yn+1 = yn + h(9yn+1 19yn 5yn-1 + yn-2)/24
(7.6)
y n+1 = f(tn+1,yn+1)
tn+1 = tn + h
Terlihat dibutuhkan tiga titik awal untuk dapat menggunakan metode AdamBashforth-Moulton. Tiga titik awal ini dapat diperoleh dengan menggunakan
metode Runge-Kutta.
Program Matlab untuk metode Adam-Bashforth-Moulton
function [tvals,yvals]=abm(func,tspan,startval,step)
%
Program ini bertujuan untuk penyelesaian
%
persamaan differensial ordiner
%
Metode Prediktor-Korektor (Adams, Bashforth, Moulton)
%
Cara menggunakan fungsi ini
%
tvals adalah variabel bebas (t)
%
yvals adalah variabel tak bebas hasil perhitungan (y)
%
func adalah nama fungsi yang ingin diselesaikan
%
tspan adalah variabel bebas
%
starval adalah nilai awal variabel bebas (yo)
%
step adalah besaran langkah untuk interval tspan
%
%
Nama File : abm.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------b=[ ];c=[ ];d=[ ];order = 4;
b=[1/6 1/3 1/3 1/6]; d=[0 0.5 0.5 1];
c=[0 0 0 0; 0.5 0 0 0; 0 0.5 0 0; 0 0 1 0];
steps=(tspan(2)-tspan(1))/step+1;
y=startval;t=tspan(1);fval(1)=feval(func,t,y);
ys(1)=startval;yvals=startval;tvals=tspan(1);
% Tiga langkah awal dengan Metode Runge-Kutta
for j=2:4
k(1)=step*feval(func,t,y);
for I=2:order
k(i)=step*feval(func,t+step*d(i),y+c(I,1:I-1)*k(1:I-1));
end;
y1=y+b*k;ys(j)=y1;t1=t+step;
175
fval(j)=feval(func,t1,y1);
tvals=[tvals,t1]; yvals=[yvals,y1];
t=t1;y=y1;
end;
% Langkah selanjutnya dg ABM orde 4
for I=5:steps
y1=ys(4)+step*(55*fval(4)-59*fval(3)+37*fval(2)-9*fval(1))/24;
t1=t+step;fval(5)=feval(func,t1,y1);
yc=ys(4)+step*(9*fval(5)+19*fval(4)-5*fval(3)+fval(2))/24;
fval(5)=feval(func,t1,yc);
fval(1:4)=fval(2:5);
ys(4)=yc;
tvals=[tvals,t1];yvals=[yvals,yc];
t=t1; y=y1;
end
Contoh 7.2.
Selesaikan
dy/dt = 2yt
dengan y = 2 untuk t = 0
Program terkait
function v=F72(t,y)
%
Program ini berisi fungsi untuk menghitung v
%
dy
%
v = ---- = 2*y*t
%
dt
%
Nama File : F72.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------v=2*t*y;
Keluaran program
176
t
0.00
0.25
0.50
0.75
1.00
1.25
1.50
1.75
2.00
Runge-Kutta
2.0000000
2.1289876
2.5680329
3.5099767
5.4357436
9.5369365
18.9519740
42.6424234
108.5814979
Prediktor-Korektor
Eksak
2.0000000
2.0000000
2.1289876
2.1289889
2.5680329
2.5680508
3.5099767
3.5101093
5.4340314
5.4365637
9.5206761
9.5414664
18.8575896
18.9754717
42.1631012
42.7618855
106.2068597
109.1963001
ode
mendekati penyelesaiannya. Fungsi ode45 adalah fungsi ode pertama yang dapat
dicoba pada masalah yang baru. Untuk melakukan hal tersebut perlu ditulis
sebuah fungsi M-File yang menghasilkan derivatif untuk nilai variabel tak bebas.
177
% Plot hasil
plot(t, y, 'k-','linewidth',2)
Program perhitungan
function dydt = F73(time, y)
%
Program ini berisi fungsi untuk menghitung v
%
dc/dt = 1 - c - xk*c
... neraca massa
%
dT/dt = 1 - T - alpha*xk*c
... neraca entalpi
%
Nama File : F73.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------global Alpha Beta
% Variabel sesungguhnya
c = y(1);
T = y(2);
% Konstanta kecepatan reaksi
xk = Beta*T*T;
% Persamaan Differensial Ordiner
dydt(1) = 1. - c - xk*c
;
dydt(2) = 1. - T - Alpha*xk*c ;
dydt = dydt';
Keluaran program
178
Recycle
Tangki I
Tangki II
Produk
179
Co = [C1o C2o];
% Batasan waktu
to = 0;
% waktu awal, menit
tn = 30;
% waktu akhir, menit
n = 31;
% jumlah data
datat = linspace(to,tn,n);
% Perhitungan
[t,C]= ode45('F74',datat,Co);
% Plot hasil
plot(t,C(:,1),'k-',t,C(:,2),'k+','LineWidth',2);
legend('C1','C2')
% Menampilkan hasil
t = t(n)
C1 = C(n,1)
C2 = C(n,2)
Program terkait
function dCdt = F74(t,C)
%
Program ini berisi fungsi untuk menghitung
%
konsentrasi garam tangki I dan II
%
sebagai fungsi waktu
%
Neraca massa tangki I
%
dC1/dt = (q0*C0 + qR*C2 - q1*C1)/V1
%
Neraca massa tangki II
%
dC2/dt = (q1*C1 - q2*C2)/V2
%
Nama File : F74.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------% Data-data
V1 = 100;
V2 = 100;
q0 = 5;
q1 = 8;
q2 = 8;
qR = 3;
qproduk = 5;
C0 = 0;
%
%
%
%
%
%
%
%
Volume tangki I, L
Volume tangki II, L
Laju alir umpan, L/menit
Laju alir keluar tangki I, L/menit
Laju alir keluar tangki II, L/menit
Laju alir recycle, L/menit
Laju alir produk, L/menit
konsentrasi umpan, L/menit
% Persamaan differensial
dCdt1 = (q0*C0 + qR*C(2) - q1*C(1))/V1;
dCdt2 = (q1*C(1) - q2*C(2))/V2;
dCdt = [dCdt1
dCdt2];
Keluaran program
180
181
O
dengan UA = 10 kJ/menit C yaitu koefisien transfer panas dan luas oil setiap
tangki. T adalah temperatur minyak dalam tangki (OC) dan Q adalah kecepatan
panas ditransfer dalam kJ/menit.
To = 20 C
W = 100 kg/menit
T1
T1
steam
T2
T2
T3
steam
T3
steam
dT1
= WCpTo + UA(Tsteam T1) WCpT1
dt
dT1
= [WCp(To T1) + UA(Tsteam T1)]/ (MCp)
dt
analog untuk tangki 2 dan 3
dT2
= [WCp(T1 T2) + UA(Tsteam T2)]/ (MCp)
dt
dT3
= [WCp(T2 T3) + UA(Tsteam T3)]/ (MCp)
dt
Program Matlab
% Penukar Panas dalam Tangki Seri
%
Larutan minyak multikomponen dipanaskan dalam
%
tiga tangki tersusun seri.
%
Sebagai pemanas digunakan steam jenuh yang
%
terkondensasi dalam koil
%
Tentukan suhu di dalam tangki sebagai fungsi waktu
%
Penyelesaian dengan fungsi ode45
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
182
clear all
clc
global W UA M Cp Tsteam Jml_Tangki To
% Data-data
Jml_Tangki = 3;
W = 100;
UA = 10;
M = 1000;
Cp = 2.0;
Tsteam = 250;
To = 20;
%
%
%
%
%
%
%
Jumlah tangki
Laju alir minyak, kg/menit
Koef. transfer panas & luas koil, kJ/menit.C
massa minyak dalam tangki, kg
kapasitas panas minyak, kj/kg
temperatur steam jenuh, C
temperatur minyak masuk, C
Program terkait
function dTdt = F75(t,T)
%
Program ini berisi fungsi untuk menghitung
%
temperatur di dalam tangki sebagai fungsi waktu
%
Neraca panas di setiap tangki
%
dTi
W Cp (To - Ti) + UA (Tsteam - Ti)
%
----- = ----------------------------------%
dt
M Cp
%
Nama File : F75.m
%
Surakarta, Oktober 2005
%
--------------------------------------------------------------global W UA M Cp Tsteam Jml_Tangki To
for j =1:Jml_Tangki
if j==1
dTdt(j)=(W*Cp*(To-T(j))+UA*(Tsteam-T(j)))/(M*Cp); % Tangki 1
else
dTdt(j)=(W*Cp*(T(j-1)-T(j))+UA*(Tsteam-T(j)))/(M*Cp);
% Tangki 2 & 3
end
end
dTdt=dTdt';
Keluaran program
183
1
B
A
2
A + C
B + C
3
C + B
2B
dC A
= -k1CA + k2CBCC
dt
dC B
= k1CA k2CBCC k3CB2
dt
dC C
= k3CB2
dt
184
CA(0) =1, CB(0) = CC(0) = 0.
Tentukan CA(10), CB(10), dan CC(10), jika k1 = 0,08, k2 =2 x 104, k3 = 6.107 !
Program Matlab
% Reaksi-reaksi Simultan
%
k1
k2
k3
%
A ----> B
B + C ----> A + C
2B ----> C + B
%
ra = -k1*Ca
%
rb = k1*Ca - k2*Cb*Cc - k3*Cb^2
%
rc = k3*Cb^2
%
Tentukan Ca, Cb, dan Cc sebagai fungsi waktu
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
clear all
clc
global k1 k2 k3
% Data konstanta kecepatan reaksi
k1 = 8e-2;
k2 = 2e-4;
k3 = 6e-2;
% Data kondisi awal
Cao = 1.0;
Cbo = 0;
Cco = 0;
Co = [Cao Cbo Cco];
% Batasan waktu
to = 0;
tn = 10;
tspan = [to tn];
% Perhitungan
[t C] = ode45('F76',tspan,Co);
% Plot Hasil
plot(t,C(:,1),'k-',t,C(:,2),'k+',t,C(:,3),'ko','LineWidth',2)
title('Tiga Reaksi Simultan','FontSize',14)
xlabel('waktu','FontSize',14)
ylabel('konsentrasi','FontSize',14)
legend('Ca','Cb','Cc',14)
Program terkait
function dC_dt=F76(t,C)
%
Program ini berisi fungsi untuk menghitung Ca, Cb, Cc
%
ra = -k1*Ca
%
rb = k1*Ca - k2*Cb*Cc - k3*Cb^2
%
rc = k3*Cb^2
%
Nama File : F76.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------global k1 k2 k3
% Variabel awal
Ca = C(1);
Cb = C(2);
Cc = C(3);
185
dC_dt(1) = -k1*Ca + k2*Cb*Cc;
dC_dt(2) = k1*Ca - k2*Cb*Cc - k3*Cb^2;
dC_dt(3) = k3*Cb^2;
dC_dt=dC_dt';
konsentrasi
Keluaran program
O
1/2 O 2
H2C
CH2
Etilen oksida (A) dan O2 (B) masuk reaktor pada temperatur 260 OC. Etilen
masuk dengan kecepatan 0,3 lbmol/detik dan tekanan 10 atm. Reaksi dijalankan
pada sepuluh unit pipa yang masing-masing unit terdiri 100 pipa diameter 1 in
schedule 40.
186
Akibatnya laju alir 3.10-4 lbmol/detik. Sifat fisis fluida bereaksi dipertimbangkan
sebagai udara pada tekanan dan temperatur yang sama. Densitas katalis yang
berdiameter in 120 lb/ft3 dan bed void fraction 0,45.
Kecepatan reaksi rA = kPA1/3PB2/3 dengan k = 0,0141 lbmol/(atm.lbkatalis.jam)
pada 260 OC. Tentukan konversi dan penurunan tekanan sebagai fungsi berat
katalis ! Diketahui o = 0,0775 atm/ft.
Penyelesaian
FAo = 3.10-4 lbmol/detik = 1,08 lbmol/jam
FBo = 1,5.10-4 lbmol/detik = 0,54 lbmol/jam
FI = 1,5.10-4 lbmol/detik x 0,79/0,21 = 5,64.10-4 lbmol/detik = 2,03 lbmol/jam
FTo = 3,65 lbmol/jam
PAo = yAoPo = (FAo / FTo)Po
PAo = (1,08/2,03).10 atm = 3 atm
-rA = kPA1/3PB2/3
-rA = k(CART)1/3(CBRT)2/3 = kRT CA1/3CB2/3
v = vo (1+X)(Po/P)
CA = FA = C Ao (1 x ) P
v
1 + x Po
CB = FB = C Ao (1 / 2 x / 2) P
v
1 + x
Po
-rA = kC Ao RTo P (1- x)1/3(1/2 - x)2/3
1 + x Po
-rA = k (1 x ) y
(1 + x )
k = 0,0141lbmol/(atm.lbkat.jam).3atm.0,63
k = 0,0266 lbmol/(lbkat.jam)
= yAo = (FAo / FTo) = (1,08/2,03)(1 1 ) = -0,15
FAo
dx = -rA
dW
187
dx
dW
= -rA / FAo
Po
dP
=
(1 + x )
dW
2 P / Po
dy
= (1 + x )
dW
2y
dengan =
2 o
Ac(1 ) kat Po
= 0,0166 / lbkat
o = G (1 )
c g c D P 3
150(1 )
+ 1,75G
Dp
188
figure(1)
plot(W,x,'k+',W,f,'k-','linewidth',2)
title('Reaksi Pembentukan Etilen Oksida','fontsize',14)
xlabel('W (lb)','fontsize',14)
legend('X','y','f','fontsize',14)
rate=kprime*((1-x(:,1))/(1+eps*x(:,1)))*x(:,2);
figure(2)
plot(W,rate,'k-','linewidth',2)
title('Kecepatan Reaksi Pembentukan Etilen Oksida','fontsize',14)
xlabel('W (lb)','fontsize',14)
Program lain
function dx_dW=F77(W,x)
%
Nama File : F77.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------global eps kprime
X = x(1);
y = x(2);
kprime=0.0266;
eps=-0.15;
alpha=0.0166;
rate=kprime*((1-X)/(1+eps*X))*y;
fa0=1;
dx_dW(1,:)=rate/fa0;
dx_dW(2,:)=-alpha*(1+eps*X)/(2*y);
Keluaran program
189
Ta
Fao
To
Ta
Berbagai nilai parameter untuk reaktor ini ditunjukkan pada tabel di bawah ini :
190
CPa = 40,0 J/gmolK
R = 8,314 j/gmolK
Hr = -40.000 J/gmol
Ua = 0,8 J/kgmenitK
Ea = 41.800 J/gmolK
Ta = 500 K
= 0,015 kg-1
Po = 10 atm
yao = 1,0
To = 450 K
Plot konversi (X), penurunan tekanan (y) dan temperatur (T.10-3) sepanjang
reaktor untuk W = 0 sampai W = 20 kg. Plot juga profil konsentrasi untuk
reaktan A dan produk C untuk harga W yang sama !
Penyelesaian
Neraca massa untuk reaktor katalitik
Fao
dX = -r
a
dW
-ra = k Ca 2 Cc
KC
191
d(P / Po) (1 + X) Po T
=
dW
2
P To
dy (1 + X) T
=
dW
2y
To
Secara umum neraca energi dapat ditulis sebagai berikut
dT Ua (Ta T) + r ' a Hr
=
dW
FaoCpa
% Reaksi Phase Gas Reversibel dan Eksotermis dlm Reaktor Katalitik
%
Reaksi elementari phase gas 2A --> C
%
dijalankan pada reaktor packed bed
%
Reaktor dilengkapi dengan penukar panas
%
terjadi penurunan tekanan sepanjang reaktor
%
Tentukan konveri, penurunan tekanan, dan temperatur
%
sebagai fungsi berat katalis
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
global Ta delH CPA FAO
% Data-data
Ta = 500;
delH = -40000;
CPA = 40;
FAO = 5;
%
%
%
%
% Kondisi awal
var0(1) = 0;
var0(2) = 450;
var0(3) = 1.0;
Wspan = [0 20];
% konversi
% temperatur
% penurunan tekanan
temperatur, K
panas reaksi, J/gmol
kapasitas panas A, J/(gmol.K)
laju alir A, gmol/menit
% Perhitungan
[W var]=ode45('F78',Wspan,var0);
% Plot hasil
tt = var(:,2)/1000;
x=var(:,1);
T=var(:,2);
y=var(:,3);
figure(1)
plot(W,x,'k-',W,y,'k+',W,tt,'ko','LineWidth',2)
title('Model Reaktor','FontSize',14)
xlabel('W','FontSize',14)
legend('konversi A', 'Tekanan ternormalisasi', 'T/1000')
temp=0.271*(450./T).*y./(1-0.5*x);
CA=temp.*(1-x);
CC=temp.*0.5.*x;
figure(2)
plot(W,CA,'k-',W,CC,'k+','LineWidth',2)
title('Konsentrasi','FontSize',14)
xlabel('W','FontSize',14)
legend('Konsentrasi A','Konsentrasi C')
192
Program terkait
function der=F78(W,var)
%
Program ini berisi fungsi untuk menghitung x, T, y
%
dx/dW = -ra'/Fao
%
dT
Ua*(Ta-T) + ra'*delH
%
---- = ---------------------%
dW
Fao*Cpa
%
dy
- alpha*(1+eps*x)
T
%
---- = ------------------ * ---%
dW
2*y
To
%
Nama File : F78.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------global Ta delH CPA FAO
% variabel awal
x = var(1);
T = var(2);
y = var(3);
% perhitungan
k = 0.5*exp(5032*(1/450-1/T));
temp = 0.271*(450/T)*y/(1-0.5*x);
CA = temp*(1-x);
CC = temp*0.5*x;
Kc = 25000*exp(delH/8.314*(1/450-1/T));
rA = -k*(CA*CA-CC/Kc);
% persamaan differensial
der(1) = -rA/FAO;
der(2) = (0.8*(Ta-T)+rA*delH)/(CPA*FAO);
der(3) = -0.015*(1-0.5*x)*(T/450)/(2*y);
der=der';
Keluaran program
193
194
Gambar 7.11.b. Keluaran Program Contoh 7.8
Tangki berisi air berbentuk bola diisi dari suatu lubang pada bagian puncak dan
dikeluarkan melalui lubang yang lain pada bagian dasar kolom. Jika jari-jari
tangki adalah r, maka volume air dalam tangki dapat dinyatakan sebagai fungsi h
(tinggi cairan) sebagai berikut
V(h) = rh2 h3
Laju alir volum melalui lubang (q) merupakan fungsi tinggi cairan (h) dan
dinyatakan dengan :
q = CdA 2gh
dengan
= luas lubang
Cd
= gravitasi
Tentukan waktu tangki menjadi kosong jika tinggi cairan mula-mula 9 ft. Tangki
mempunyai jari-jari r = 5 ft dan lubang mempunyai diameter 1 in di dasar tangki.
Gunakan g = 32,2 ft/second2.
dV
= -q
dt
195
2rh
dh
dh
- h2
= - CdA
dt
dt
2gh
dh = C d A 2gh
dt
( 2rh h 2 )
% Tinggi Cairan dalam Tangki Sperikal
%
Menentukan tinggi cairan dalam tangki sebagai fungsi waktu
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
clear all
% kondisi awal
ho = 9;
% ft
% Batasan waktu
to = 0;
tn = 2475;
tspan = [to tn];
% menit
% menit
% Perhitungan
[t,h]=ode45('F79',tspan,ho);
% Plot hasil
plot(t,h,'k-','LineWidth',2)
xlabel('Waktu (detik)','FontSize',14)
ylabel('Tinggi (ft)','FontSize',14)
Program terkait
function hdot = F79(t,h)
%
dh
Cd*A*(2*g*h)^0.5
%
---- = - ------------------%
dt
phi*(2*r*h - h^2)
%
Nama File : F63.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------%Data-data
Cd = 0.6;
r = 5;
dhole = 1/12;
A = pi*(dhole)^2;
g = 32.2;
%
%
%
%
jari-jari tangki, ft
diameter lubang, konversi ke ft
luas lubang pengeluaran
ft/det^2
hdot = -Cd*A*(2*g*h)^0.5/(pi*(2*r*h-h^2));
Tinggi(ft)
Keluaran program
196
Sebanyak 20.000 barrel minyak diproses per hari untuk meningkatkan angka
oktannya dalam reaktor reforming. Laju alir umpan adalah 44 kg per detik atau
440 mol per detik.
Reaksi dehidrogenasi
Parafin olefin + H2
Reaksi orde satu terhadap parafin. Asumsi parafin masih murni pada tekanan
2000 kPa dan konsentrasi 0,32 mol/dm3. Tentukan konsentrasi dan penurunan
tekanan ketika reaksi ini dijalankan dalam reaktor bola berdiameter 27 dm !
-rA
= kCA
-rA
= B (-rA) = C (1 - ) kCA
Berat katalis
= 173,870 kg
= 0,032 kg/dm3
Dp
= 0,02 dm
= 0,4
197
3
= 0,02 dm /(kgkatdetik)
= 1,5.10-6 kg/dmdetik
L = L
= 27 dm
= 2,6 kg/dm3
dFA
= rA Ac
dz
dx rA Ac
=
dz
FAo
rA = - kCA = - k CA C (1 - )
CA = CAo 1 x y
1 + x
= yAo = 1 x (1+1-1) = 1
y = P/Po
Persamaan untuk penurunan tekanan
dy
= o (1 + x )
dz
Po y
o = G (1 ) 150(1 ) + 1,75G
c g c D P 3
Dp
G = m / Ac
Reaktor sperikal
Ac = [R2 (z L)2]
o =
0,02
0,032.0,02.0,4 3
o = [98,87G + 25630G2]0,01
w = C (1 - ) [R2z 1/3(z L)3 - 1/3 L3]
% Reaksi Dehidrogenasi dalam Reaktor Sperikal
%
Reaksi dehidrogenasi
%
Parafin --> olefin + H2
%
Reaksi orde satu terhadap parafin.
%
Tentukan konsentrasi dan penurunan tekanan
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
global eps kprime
198
% Data awal
Xo = 0;
yo = 1;
xo = [Xo yo];
% batasan jarak
zspan = [0 54];
%Perhitungan
[z,x]=ode45('F710',zspan,xo);
% Plot hasil
plot(z,x,'k-','linewidth',2)
title('Reaksi Dehidrogenasi pada Reaktor Bola')
xlabel('z (dm)','FontSize',12)
legend('X','y',12)
axis([0 54 0 1.2])
Program terkait
function dx_dz=F710(z,x)
%
Nama File : F710.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------% Variabel awal
X = x(1);
y = x(2);
% Data-data
Fa0 = 440;
P0 = 2000;
Ca0 = 0.32;
R = 30;
phi = 0.4;
kprime = 0.02;
L = 27;
rhocat = 2.6;
m = 44;
% Perhitungan
Ca = Ca0*(1-X)*y/(1+X);
Ac = pi*(R^2-(z-L)^2);
V=pi*(z*R^2-1/3*(z-L)^3-1/3*L^3);
G=m/Ac;
ra=-kprime*Ca*rhocat*(1-phi);
beta=(98.87*G+25630*G^2)*0.01;
W=rhocat*(1-phi)*V;
% Differnsial
dx_dz(1,:)=-ra*Ac/Fa0;
dx_dz(2,:)=-beta/P0/y*(1+X);
Keluaran program
199
15200
dx1
T p2 pDpH
= 14,96.10 6 e
B
d(V F)
K1
r2 =
15200
dx 2
T p p pT pH
= 8,67.10 6 e
B D
d(V F)
K2
T = temperatur, 1033 K
K1 = 0,312
200
x2 = konversi pada reaksi 2
K2 = 0,480
dx1
= 6,089 (1 x1 x 2 )2
d (V F)
(12 x1 x 2 )(12 x1 + x 2 )
x 2 12 x1 + x 2
dx 2
= 3,529 (1 x1 x 2 ) 12 x1 x 2
0,480
d(V F)
0,312
% Temperatur, K
% konstanta kesetimbangan reaksi 1
% konstanta kesetimbangan reaksi 2
% Kondisi awal
x1o = 0;
x2o = 0;
xo = [x1o x2o];
% Batasan V/F
V_Fo = 0;
% V/F awal, ft^3.jam/lbmol
V_Ff = 0.5;
% V/F awal, ft^3.jam/lbmol
V_Fspan = [V_Fo V_Ff];
% Perhitungan
[V_F,x]=ode45('F711',V_Fspan,xo);
% Plot hasil
plot(V_F,x(:,1),'k-',V_F,x(:,2),'k+','Linewidth',2)
title('Dehidrogenasi Benzena','FontSize',12)
xlabel('V/F dlm ft^3.jam/lbmol','FontSize',12)
ylabel('konversi Benzena','FontSize',12)
legend('reaksi 1', 'reaksi 2')
disp(['Pada V/F ',num2str(V_F(end))])
konversi = x(end,:)
Program terkait
function dx_dVF = F711(V_F,x)
%
Program ini berisi fungsi untuk menghitung
201
%
%
%
%
%
%
%
%
%
%
global T K1 K2
dx_dVF(1) = 14.96*10^6*exp(-15200/T)*((1-x(1)-x(2))^2-...
(0.5*x(1)-x(2))*(0.5*x(1)+x(2))/K1);
dx_dVF(2) = 8.67*10^6*exp(-15200/T)*((1-x(1)-x(2))*...
(0.5*x(1)-x(2))-x(2)*(0.5*x(1)+x(2))/K2);
dx_dVF=dx_dVF';
konversi Benzena
Keluaran program
Proses distilasi secara biner melibatkan komponen 1 dan 2. Mol liquid tersisa (L)
dinyatakan sebagai fungsi fraksi mol komponen ke-2, x2 sebagai berikut :
dL
L
=
dx 2 x 2 (k 2 1)
202
dengan k2 adalah rasio kesetimbangan uap cair komponen ke-2. Jika sistem
dianggap berada pada keadaan ideal, rasio kesetimbangan uap cair dapat dihitung
sebagi ki = Pi /P dengan Pi adalah tekanan uap komponen i dan P adalah tekanan
total.
Secara umum model tekanan uap menggunakan persamaan Antoine dengan T
adalah temperatur (OC).
Pi = 10^ A B
T +C
k1x1 + k2x2 = 1
Untuk sistem biner dengan komponen benzena (komponen ke-1) dan toluena
(komponen ke-2) diasumsikan berada pada kesetimbangan. Konstanta Antoine
untuk benzena A1 = 6,90565, B1 = 1211,033, dan C1 = 220,79. Sedang untuk
toluena A2 = 6,95464, B2 = 1344,8, dan C2 = 219,482. P adalah tekanan (mmHg)
dan T adalah temperatur (OC).
Hitunglah jumlah liquid tersisa
mencapai 80 %, jika diketahui 100 mol liquid umpan terdiri 60% benzena dan
40% toluena (fraksi mol) pada tekanan 1,2 atm.
Penyelesaian
% Distilasi Secera Batch Komponen Biner
%
Proses distilasi secara biner melibatkan komponen 1 dan 2
%
Hitunglah jumlah liquid tersisa dalam distilasi
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
clear all
clc
global A B C P
%
A
B
C
% Data-data proses
P = 1.2*760;
Lo = 100;
x2_awal = 0.4;
x2_akhir = 0.8;
%
%
%
%
mmHg
mol
mol toluena awal
mol toluena akhir
203
% Penyelesaian PD Ordiner
x2span = [x2_awal x2_akhir];
[x2 L] = ode45('F712',x2span,Lo);
% Plot hasil
plot(x2,L,'k-','linewidth',2)
title('Distillasi Batch','fontsize',14)
xlabel('Fraksi mol toluena','fontsize',14)
ylabel('Mol cairan','fontsize',14)
% Save hasil
output = [x2 L];
save batch.dat output -ascii
Program terkait 1
function dLdx2 = F712(x2,L)
%
Program ini berisi fungsi untuk menghitung PD ordiner
%
dL
L
%
----- = --------------%
dx2
x2 ( k2 - 1 )
%
Nama File : F712.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------global A B C P x2
T_tebak = (80.1+110.6)/2;
T = fzero('F712F', T_tebak);
%
%
% T
%
P_i = 10.^(A-B./(T+C));
k = P_i/P;
dLdx2 = L/x2/(k(2)-1);
Program terkait 2
function f = F712F(T)
%
Program ini berisi fungsi untuk menentukan T
%
(
B
)
%
Pi = 10^(A - -------)
%
(
T + C )
%
Pi
%
ki = ---%
P
%
k1*x1 + k2*x2 = 1
%
Nama File : F712F.m
%
Surakarta, Oktober 2005
%
-------------------------------------------------------------global A B C P x2
x1 = 1-x2;
P_i = 10.^(A-B./(T+C));
k = P_i/P;
f = 1-k(1)*x1-k(2)*x2;
Keluaran program
Mol cairan
204
205
end
w1 = zbc(1,3);
206
a2 = zbc(2,1);
%
%
%
%
%
%
%
%
%
%
%
%
b2 = zbc(2,2);
&
w2 = zbc(2,3);
icnt = 1;
icntmax = 25;
Contoh 7.13
+ 3xy + 7 y = cos(2 x )
207
Program matlab
% Demonstrasi penyelesaian Persamaan Differensial Ordiner
%
Masalah Nilai Batas
%
Program ini bertujuan menyelesaikan
%
y'' + 3xy' + 7y = cos(2x)
%
y(0) = 1 and y(pi) = 0
%
Penyelesaian dengan Metode Shooting
%
File prepared by J. R. White, UMass-Lowell (Aug. 2003)
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
clear all, close all, nfig = 0;
%
Nilai ALF tebakan adalah 1, -1, dan -5.4713
gs = [1 -1 -5.4715];
%
%
set kondisi batas dan penyelesaian IVP
xo = 0;
xf = pi;
tol = 1.0e-6;
yxo = 1.0;
options = odeset('RelTol',tol);
[x1,z1] = ode23('F713',[xo xf],[yxo gs(1)],options);
[x2,z2] = ode23('F713',[xo xf],[yxo gs(2)],options);
[x3,z3] = ode23('F713',[xo xf],[yxo gs(3)],options);
%
%
plot hasil
nfig = nfig+1; figure(nfig)
plot(x1,z1(:,1),'k-',x2,z2(:,1),'k-.',x3,z3(:,1),'k-','LineWidth',2)
title('Metode Shooting (Iterasi Manual)')
xlabel('x'),ylabel('y'),grid
legend('alf = 1','alf = -1','alf = -5.4715')
Program terkait
function zp = odefile(x,z)
%
Persamaan Differensial Ordiner
%
y'' + 3xy' + 7y = cos(2x)
%
Nama File : F713.m
%
Surakarta, Oktober 2005
%
File prepared by J. R. White, UMass-Lowell (March 2003)
%
--------------------------------------------------------------zp = zeros(length(z),1);
zp(1) = z(2);
zp(2) = -3*x*z(2) - 7*z(1) + cos(2*x);
Keluaran program
208
209
% Plot hasil
plot(xs,zs(:,1),'k','LineWidth',2),grid
title('Metode Shooting untuk PDO orde 2')
Keluaran program
210
dy yi +1 yi
dx
x
(forward)
dy yi yi 1
dx
x
(backward)
dy yi +1 yi 1
dx
2x
(central)
d 2 y d dy
=
dx 2 dx dx
y i+1 y i y i y i1
x = y i+1 2 y i + y i1
x
x
(x )2
Contoh 7.14
+ 3xy + 7 y = cos(2 x )
xf xo
N +1
y i 1 2 y i + y i +1
x 2
yi = yi +1 yi 1
2x
+ 3xy + 7 y = cos(2 x )
y i 1 2 y i + y i +1
y y i 1 + 7y = cos(2x )
+3xi i +1
i
i
2
2x
x
211
3x i x
2
2
3x x
1
yi-1 + (-2 + 7x )yi + 1 + i yi+1 = x cos(2xi)
2
2
selanjutnya
A(i,i-1) = 1
3x i x
2
A(i,i) = -2 + 7x2
A(i,i+1) = 1 + 3x i x
2
b(i) = x cos(2xi)
untuk i = 1, yi-1 = yo = yxo
3x x
(-2 + 7x2)y1 + 1 + i y2 = x2cos(2x1)
2
3x 1 x yxo
1
Kasus B
Pendekatan y =
d2y
dt 2
y i +1 y i
x
+ 3xy + 7 y = cos(2 x )
y i 1 2 y i + y i +1
yi +1 yi + 7y = cos(2x )
i
i
+3xi x
2
x
212
%
%
%
%
%
%
%
%
%
% Kondisi batas
xo = 0;
xf = pi;
yxo = 1;
yxf = 0;
% jumlah titik yang berbeda
NN = [20 40 60 160];
for n=1:4
N = NN(n);
% substitusi finite difference
dx = (xf-xo)/(N+1);
dx2 = dx*dx;
x = (xo+dx):dx:(xf-dx);
% Metode 1 : central
a = zeros(N,N); b = zeros(N,1);
% titik awal
a(1,1) = -2+7*dx2;
a(1,2) = 1+3*x(1)*dx/2;
b(1) = dx2*cos(2*x(1))-(1-3*x(1)*dx/2)*yxo;
% titik ke 2 sampai ke N-1
for i = 2:N-1
a(i,i-1) = 1-3*x(i)*dx/2;
a(i,i) = -2+7*dx2;
a(i,i+1) = 1+3*x(i)*dx/2;
b(i) = dx2*cos(2*x(i));
end
% titik ke N
a(N,N-1) = 1-3*x(N)*dx/2;
a(N,N) = -2+7*dx2;
b(N) = dx2*cos(2*x(N))-(1-3*x(N)*dx/2)*yxf;
% penyelesaian persamaan aljabar simultan
y = a\b;
za = [yxo y' yxf];
xa = [xo x xf];
% Metode 2 : forward
a = zeros(N,N); b = zeros(N,1);
% titik awal
a(1,1) = -(2+3*x(i)*dx-7*dx2);
a(1,2) = 1+3*x(i)*dx/2;
b(1) = dx2*cos(2*x(1))-yxo;
% titik ke 2 sampai ke N-1
for i = 2:N-1
a(i,i-1) = 1;
a(i,i) = -(2+3*x(i)*dx-7*dx2);
a(i,i+1) = 1+3*x(i)*dx;
b(i) = dx2*cos(2*x(i));
end
% titik ke N
a(N,N-1) = 1;
a(N,N) = -(2+3*x(N)*dx+7*dx2);
b(N) = dx2*cos(2*x(N))-(1-3*x(N)*dx)*yxf;
% penyelesaian persamaan aljabar simultan
y = a\b;
213
zb = [yxo y' yxf];
xb = [xo x xf];
% plot hasil
t = '220+n';
subplot(eval(t)),plot(xa,za,xb,zb,'LineWidth',2)
axis([0 3.2 -2 1]);
title(['Metode Finite Difference (',num2str(N),' titik)'])
legend('Metode Central','Metode Forward')
end
Keluaran hasil
qr
qr+ r
qc = 0
214
Dalam keadaan steady
qr qr+r qc = 0
- kAr
k 2r
dT
dT
+ kAr
dr r
dr
dT
dr
k 2r
r + r
hAc (T - T) = 0
r + r
dT
dr
2[2r+r]h(T - T) = 0
limit r 0
dT
d
k 2r 4rh(T T) = 0
dr
dr
2
d T
dT
2h r2(T T) = 0
r2 dr 2 + r
dr
k
T(rw) = Tw
dan
dT
dr
=0
rs
T T
Tw T
dan
x=
r
rs
rs
du du dr
dT
=
=
dx dr dx Tw T dr
rs2
d 2 u d du dr
d 2T
=
=
dx 2 dr dx dx Tw T dr 2
Tw T
(rsx)
2
rs
2
d2u
Tw T du
2h
2
du
d2u
+x
2x2u = 0
x
2
dx
dx
2
x2 u + x u 2x2u = 0
Kondisi batas
2hrs2
dengan =
k
2
215
x=
rw
=a
rs
u(a) = 1
x=
rs
=b=1
rs
u(b) = 0
rw = 1 in
rs = 1,5 in
O
= 0,0625 in
Tw = 200 F
T = 70 F
h = 20 BTU/jam.ft2.oF
h = 75 BTU/jam.ft.oF
Tentukan distribusi panas, total hilang panas fin, dan hitung effisiensi, , yaitu
=
Penyelesaian
x2 u + x u 2x2u = 0
Kondisi batas
x=
rw
=a
rs
u(a) = 1
x=
rs
=b=1
rs
u(b) = 0
Shooting Method
z1 = z2
z2 =
1
z2 + 2z1
x
1
z(a) = u (a ) =
u ' (a )
216
Program penyelesaian
% Transfer Panas dalam Circular Fin
%
Program ini bertujuan menyelesaikan
%
x2 u" + x u' - a2x2u = 0
%
z1' = z2
%
z2' = - (1/x) z2 + a2z1
%
u(a) = 1
u'(b) = 0
%
Penyelesaian dengan Metode Shooting
%
File prepared by J. R. White, UMass-Lowell (Aug. 2003)
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
clear all
clc
global alpha
rw = 1/12;
rs = 1.5/12;
thk = 0.0625/12;
Tw = 200;
Tinf = 70;
h = 20;
k = 75;
%
%
%
%
%
%
%
a = rw/rs;
b = rs/rs;
alf2 = (2*h*rs*rs)/(k*thk);
alpha = sqrt(alf2);
Qideal = 2*pi*h*(rs*rs-rw*rw)*(Tw-Tinf);
fid = fopen('cylfinsh.out','w');
fprintf(fid,'\n Data dan Hasil Perhitungan \n');
fprintf(fid,'\n \n Data \n');
fprintf(fid,'Jari-jari inside, rw
= %6.2f ft \n',rw);
fprintf(fid,'Jari-jari outside, rs
= %6.2f ft \n',rs);
fprintf(fid,'Tebal fin, thk
= %6.2f ft \n',thk);
fprintf(fid,'Temperatur inside, Tw
= %6.2f F \n',Tw);
fprintf(fid,'Temperatur ambient, Tinf
= %6.2f F \n',Tinf);
fprintf(fid,'Koef. transfer panas, h
= %6.2f BTU/jamft^2F \n',h);
fprintf(fid,'Konduktivitas termal, k
= %6.2f BTU/jamftF \n',k);
zbc = [1 0 1.0
0 1 0.0];
[x,z] = bvp2sh('F715',[a b],zbc);
[nr,nc] = size(z);
Qactual = -k*(2*pi*a*thk)*(Tw-Tinf)*z(1,2);
Ttip = Tinf+(Tw-Tinf)*z(nr,1);
fineff = Qactual/Qideal;
fprintf(fid,'\n \n-------------------------------\n');
fprintf(fid, 'Temperatur dinding (F)
= %8.3f \n',Tw);
fprintf(fid, 'Temperatur ambient (F)
= %8.3f \n',Tinf);
fprintf(fid, 'Temperatur tip (F)
= %8.3f \n',Ttip);
fprintf(fid, 'Q ideal (BTU/jam)
= %8.3f \n',Qideal);
fprintf(fid, 'Q actual (BTU/jam)
= %8.3f \n',Qactual);
fprintf(fid, 'Eff fin
= %8.3f \n',fineff);
subplot(2,1,1)
plot(x,z(:,1),'k','LineWidth',2)
title('Profil Temperatur pd Fin Lingkaran dg Metode Shooting')
grid, ylabel('Temperatur')
217
subplot(2,1,2)
plot(x,z(:,2),'k','LineWidth',2)
title('Gradien Temperatur pd Fin Lingkaran dg Metode Shooting')
grid, xlabel('jarak'), ylabel('grandien temperatur')
fclose(fid)
Program terkait
function zp = odefile(x,z)
%
Persamaan Differensial Ordiner
%
x2 u" + x u' - a2x2u = 0
%
Nama File : F715.m
%
Surakarta, Oktober 2005
%
File prepared by J. R. White, UMass-Lowell (March 2003)
%
--------------------------------------------------------------global alpha
zp = zeros(length(z),1);
zp(1) = z(2);
zp(2) = -z(2)/x+alpha*alpha*z(1);
Gradien temperatur
Temperatur
Keluaran program
218
xi+1 = xi + x
x xi
u(x) u(xi) = ui
u(x) u(xi) = ui
x =
ba
N
ui =
u i 1 2u i + u i +1
x 2
ui =
u i +1 u i 1
2x
u i 1 2u i + u i +1 + x u i +1 u i 1 2x 2u = 0
i i
i
x 2
2 x
xi2
dengan 2 =
x
1
2x i
2hrs2
k
ui-1 (2 + 2x2)ui +
A(i,i-1) = 1
x u = 0
i+1
1 +
2x i
x
2x i
A(i,i) = (2 + 2x2)
A(i,i+1) = 1 +
x
2x i
b(i) = 0
untuk i = 1, u0 = u(a) = 1
(2 + 2x2)u1 +
x
1 +
u2 =
2x 1
x
1
u0
2x 1
dT
dr
hAc(TN - T) = 0
r = rs r / 2
qxN = 0
219
Ar = 2(rs -
r
)
2
Ac =2 rs rs
2
dT
dr
=
r = rs r / 2
r 2
=
2
r
r
hr
f N (TN T )
k
dengan fN =
du
dx
du
dx
=
xN x / 2
rs
dT
Tw T dr
(rs r / 4)
(rs r / 2)
2
= 2hrs 1 r f (TN T )
N
k 2 rs (Tw T )
rs r / 2
xN x / 2
= 2 2 fNuN
atau
du
dx
=
xN x / 2
x
2
u N u N 1
x
x 2
1 + 2
f N uN + uN-1 = 0
2
Program Matlab
% Transfer Panas dalam Circular Fin
%
Program ini bertujuan menyelesaikan
%
x2 u" + x u' - a2x2u = 0
%
z1' = z2
%
z2' = - (1/x) z2 + a2z1
%
u(a) = 1
u'(b) = 0
%
Penyelesaian dengan Metode Finite Difference
%
File prepared by J. R. White, UMass-Lowell (Aug. 2003)
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
clear all
clc
rw = 1/12;
rs = 1.5/12;
220
thk = 0.0625/12;
Tw = 200;
Tinf = 70;
h = 20;
k = 75;
%
%
%
%
%
a = rw/rs;
b = rs/rs;
alf2 = (2*h*rs*rs)/(k*thk);
alpha = sqrt(alf2);
Qideal = 2*pi*h*(rs*rs-rw*rw)*(Tw-Tinf);
fid = fopen('cylfinfd.out','w');
fprintf(fid,'\n Data dan Hasil Perhitungan \n');
fprintf(fid,'\n \n Data \n');
fprintf(fid,'Jari-jari inside, rw
= %6.2f ft \n',rw);
fprintf(fid,'Jari-jari outside, rs
= %6.2f ft \n',rs);
fprintf(fid,'Tebal fin, thk
= %6.2f ft \n',thk);
fprintf(fid,'Temperatur inside, Tw
= %6.2f F \n',Tw);
fprintf(fid,'Temperatur ambient, Tinf
= %6.2f F \n',Tinf);
fprintf(fid,'Koef. transfer panas, h
= %6.2f BTU/jamft^2F \n',h);
fprintf(fid,'Konduktivitas termal, k
= %6.2f BTU/jamftF \n',k);
u0 = 1;
N = input('Enter nilai N yg diinginkan (tekan 0 untuk keluar) = ');
while N ~= 0
clear A B u x xx zz zzp
A = zeros(N,N);
B = zeros(N,1);
dx = (b-a)/N;
dx2 = dx*dx;
x = linspace(a+dx,b,N);
fn = (b-dx/4)/(b-dx/2);
for n = 2:N-1
A(n,n-1) = 1-dx/(2*x(n));
A(n,n) = -(2+alf2*dx2);
A(n,n+1) = 1+dx/(2*x(n));
B(n) = 0;
end
A(1,1) = -(2+alf2*dx2);
A(1,2) = 1+dx/(2*x(1));
B(1) = -(1-dx/(2*x(1)))*u0;
A(N,N-1) = 1;
A(N,N) = -(1+alf2*dx2*fn/2);
B(N) = 0;
u = A\B;
xx = [a x];
zz = [u0 u'];
zzp(1) = (zz(2)-zz(1))/dx;
zzp(N+1) =0;
for n=2:N
zzp(n) = (zz(n+1)-zz(n-1))/(2*dx);
end
Qactual1 = -k*(2*pi*a*thk)*(Tw-Tinf)*zzp(1);
dr=dx*rs;
Qactual2 = h*2*pi*((rw+dr/2)^2-rw^2)*(Tw-Tinf)*zz(1) + ...
h*2*pi*(rs^2-(rs-dr/2)^2)*(Tw-Tinf)*zz(N+1);
for n = 2:N
r2 = (xx(n)+dx/2)*rs;
r1 = (xx(n)-dx/2)*rs;
area = 2*pi*(r2^2-r1^2);
Qactual2=Qactual2+h*area*(Tw-Tinf)*zz(n);
221
end
Ttip = Tinf+(Tw-Tinf)*zz(N+1);
fineff = Qactual2/Qideal;
fprintf(fid,'\n \n-------------------------------\n');
fprintf(fid, 'Temperatur dinding (F)
= %8.3f \n',Tw);
fprintf(fid, 'Temperatur ambient (F)
= %8.3f \n',Tinf);
fprintf(fid, 'Temperatur tip (F)
= %8.3f \n',Ttip);
fprintf(fid, 'Q ideal (BTU/jam)
= %8.3f \n',Qideal);
fprintf(fid, 'Q actual (BTU/jam)
= %8.3f (konduksi)\n',
Qactual1);
fprintf(fid, 'Q actual (BTU/jam)
= %8.3f (konveksi)\n',
Qactual2);
fprintf(fid, 'Eff fin
= %8.3f \n',fineff);
subplot(2,1,1)
plot(xx,zz,'LineWidth',2)
title(['Profil Temperatur pd Fin Lingkaran dg Metode FD . . .
untuk N = ', num2str(N)])
grid, ylabel('Temperatur')
subplot(2,1,2)
plot(xx,zzp,'LineWidth',2)
title(['Gradien Temperatur pd Fin Lingkaran dg Metode FD . . .
untuk N = ', num2str(N)])
grid, xlabel('jarak'), ylabel('grandien temperatur')
N = input('Enter nilai N lain yang diinginkan (tekan 0 untuk
keluar) = ');
end
fclose(fid)
Gradien temperatur
Temperatur
222
PERSAMAAN
DIFFERENSIAL PARSIAL
2u
x 2
+B
u u
2u
2u
+ C 2 + D x , y, u , , = 0
x y
xy
y
(8.1)
228
Persamaan Laplace
A
2u
x 2
+C
2u
=0
y 2
(8.2)
Persamaan Poisson
A
2u
x
+C
u u
+ D x, y, u, , = 0
x y
2u
y 2
(8.3)
2u
y 2
=0
(8.4)
2u
y 2
u i +1, j 2u i, j + i 1, j
(x )
u i, j+1 2u i, j + i, j1
(y )2
(8.5)
Jika diambil x = y = h
2u
x
2u
y
1
h2
(8.6)
229
Program Matlab
% Profil Temperatur pada Plat Rektanguler
%
Persamaan Laplace
%
d^2 U
d^2 U
%
----- + ----- = 0
%
d x^2
d y^2
%
Temperatur salah satu sisi = 100*sin(Pi*y)
%
sedang ketiga sisi yang lain = 0
%
Penyelesaian dengan Metode Finite Difference
%
%
Author's Data: Housam BINOUS
%
Department of Chemical Engineering
%
National Institute of Applied Sciences and Technology
%
Tunis, TUNISIA
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
% Jumlah titik adalah 29*29
% Penyelesaian sistem AU = X dengan U adalah temperatur
% yang tidak diketahui pada titik interior.
clear;
N=30;
for j=0:N-2
for i=1:N-2
X(i+j*(N-1))=0;
end
end
for i=1:N-1
X(i*(N-1))=-100*sin(i*pi/N);
end
% Penyusunan matriks A
for i=1:(N-1)*(N-1)
A(i,i)=-4;
end
for i=1:N-2
for k=0:N-2
A(i+k*(N-1),i+1+k*(N-1))=1;
end
end
for k=0:N-3
for i=1:N-1
A(i+k*(N-1),i+(k+1)*(N-1))=1;
end
end
for i=1:(N-1)*(N-1)
for j=1:i
A(i,j)=A(j,i);
end
end
% Inversi Matriks dan Perhitungan temperatur
M=inv(A);
U=M*X';
% Plot hasil pada bentuk contour
for i=1:N-1
for j=1:N-1
x(i,j)=U(j+(i-1)*(N-1));
end
230
end
[i,j]=meshgrid(1:1:N-1,1:1:N-1);
[c,h]=contourf(i,j,x);
Penyelesaian
2u
x
2u
y 2
=0
dengan u(x,0) = 0,
u(x,10) = 0,
u(0,y) = 0,
231
u(20,y) = 100.
Program Matlab
% Profil Temperatur pada Plat Rektanguler
%
Persamaan Laplace
%
d^2 U
d^2 U
%
----- + ----- = 0
%
d x^2
d y^2
%
Plat ukuran 10cm x 20cm
%
Temperatur salah satu sisi = 100
%
sedang ketiga sisi yang lain = 0
%
Penyelesaian dengan Metode Finite Difference
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
% Jumlah titik adalah 9*19
% Penyelesaian sistem AU = X dengan U adalah temperatur
% yang tidak diketahui pada titik interior.
clear;
N=20;
M=10;
for j=0:N-2
for i=1:M-2
X(i+j*(M-1))=0;
end
end
for i=1:M-1
X(i*(N-1))=-100;
end
% Penyusunan matriks A
for i=1:(N-1)*(M-1)
A(i,i)=-4;
end
for i=1:N-2
for k=0:M-2
A(i+k*(N-1),i+1+k*(N-1))=1;
end
end
for k=0:M-3
for i=1:N-1
A(i+k*(N-1),i+(k+1)*(N-1))=1;
end
end
for i=1:(N-1)*(M-1)
for j=1:i
A(i,j)=A(j,i);
end
end
% Inversi Matriks dan Perhitungan Temperatur
G=inv(A);
U=G*X';
% Plot hasil bentuk contour
for i=1:M-1
for j=1:N-1
x(i,j)=U(j+(i-1)*(N-1));
end
232
end
T = x
[i,j]=meshgrid(1:1:N-1,1:1:M-1);
[c,h]=contourf(i,j,x);
Keluaran program
233
test=0;
if F==zeros(ny+1,nx+1), test=1; end
if bx0==zeros(1,ny+1), test=test+1; end
if bxn==zeros(1,ny+1), test=test+1; end
if by0==zeros(1,nx+1), test=test+1; end
if byn==zeros(1,nx+1), test=test+1; end
if test==5
disp(' WARNING ')
disp(' ')
break
end
bx0=bx0(1,ny+1:-1:1); bxn=bxn(1,ny+1:-1:1);
a(1,:)=byn; a(ny+1,:)=by0;
a(:,1)=bx0'; a(:,nx+1)=bxn';ncase=1;
end
for i=2:ny
for j=2:nx
nn=(i-2)*(nx-1)+(j-1);
q(nn,1)=i; q(nn,2)=j; p(i,j)=nn;
end
end
C=zeros(nmax,nmax); e=zeros(nmax,1); om=zeros(nmax,1);
if ncase==1, g=zeros(nmax,1); end
for i=2:ny
for j=2:nx
nn=p(i,j); C(nn,nn)=-(2+2*r^2); e(nn)=hy^2*G(i,j);
if ncase==1, g(nn)=g(nn)+hy^2*F(i,j); end
if p(i+1,j)~=0
np=p(i+1,j); C(nn,np)=1;
else
if ncase==1, g(nn)=g(nn)-by0(j); end
end
if p(i-1,j)~=0
np=p(i-1,j); C(nn,np)=1;
else
if ncase==1, g(nn)=g(nn)-byn(j); end
end
if p(i,j+1)~=0
np=p(i,j+1); C(nn,np)=r^2;
else
if ncase==1, g(nn)=g(nn)-r^2*bxn(i); end
end
if p(i, j-1)~=0
np=p(i,j-1); C(nn,np)=r^2;
else
if ncase==1, g(nn)=g(nn)-r^2*bx0(i); end
end
end
end
if ncase==1
C=C+diag(e); z=C\g;
for nn=1:nmax
i=q(nn,1); j=q(nn,2); a(i,j)=z(nn);
end
else
[u,lam]=eig(C,-diag(e));
[om,k]=sort(diag(lam)); u=u(:,k);
for nn=1:nmax
i=q(nn,1); j=q(nn,2);
a(i,j)=u(nn,mode);
end
end
234
Tentukan distribusi temperatur dalam suatu plat rektanguler, dengan kondisi
batas sebagai berikut
x = 0, T = 100y
x = 3, T = 250y
y = 0, T = 0
y = 2, T = 200 + (100/3)x2
Penyelesaian untuk ukuran 6 x 6.
Program Matlab
% Profil Temperatur pada Plat Rektanguler
%
Persamaan Laplace
%
Plat ukuran 2 x 3
%
pada x = 0, T = 100y
%
pada x = 3, T = 250y
%
pada y = 0, T = 0
%
pada y = 2, T = 200 + (100/3)x^2
%
Penyelesaian dengan Metode Finite Difference
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
clear all
% Jumlah titik arah x, y
nx=6; ny=6;
% Ukuran titik arah x, y
hx=0.5; hy=0.3333;
% Input data pada kondisi batas
by0=[0 0 0 0 0 0 0];
byn=[200 208.33 233.33 275 333.33 408.33 500];
bx0=[0 33.33 66.67 100 133.33 166.67 200];
bxn=[0 83.33 166.67 250 333.33 416.67 500];
% Penyelesaian dg fungsi ellipgen
F=zeros(ny+1,nx+1); G=F;
% PD Laplace
a=ellipgen(nx, hx, ny, hy, G, F, bx0, bxn, by0, byn);
% Plot hasil
contourf(a)
xlabel('Titik-titik dalam arah x');
ylabel('Titik-titik dalam arah y');
Keluaran program
235
2u
x
+C
u u
= F x, y, u , ,
x y
y
2u
2
Permasalahan
ini
mengikuti
persamaan
Poisson
dengan
F(x,y)
236
% Ukuran titik arah x, y
hx=1/6; hy=1/6;
% Input data
by0=[0 0 0 0
byn=[0 0 0 0
bx0=[0 0 0 0
bxn=[0 0 0 0
Keluaran program
237
2u
x
c u
k t
(8.7)
u ij+1 + 2u ij + u ij1
(8.8)
(x )2
j+1
j
u u i u i
=
t
t
(8.9)
Substitusi persamaan (8.8) dan (8.9) pada persamaan (8.7) dan penyelesaian
untuk uij+1
u ij+1 =
kt
c(x )
Penyederhanaan jika
u ij+1 =
(u
j
i +1
kt
c(x )2
1 j
u i +1 + u ij1
2
2kt
+ u ij1 + 1
c(x )2
j
u
i
(8.10)
(8.11)
Contoh 8.5
238
uin = 0.5;
uA =1;
uB = 0.2;
Keluaran program
>> kondisi awal t = 0
x =
0
0.0400
0.0800
0.1200
0.1600
0.2000
0.2400
0.2800
0.3200
0.3600
0.4000
0.4400
0.4800
0.5200
0.5600
0.6000
239
0.6400
0.6800
0.7200
0.7600
0.8000
0.8400
0.8800
0.9200
0.9600
1.0000
u0 =
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
0.5000
kondisi t >= 0
x =
0
0.0400
0.0800
0.1200
0.1600
0.2000
0.2400
0.2800
0.3200
0.3600
0.4000
0.4400
0.4800
0.5200
0.5600
0.6000
0.6400
0.6800
0.7200
0.7600
0.8000
0.8400
0.8800
0.9200
0.9600
1.0000
240
t =
0.0080
0.0160
0.0240
0.0320
0.0400
1.0000
0.8714
0.7557
0.6620
0.5939
0.5494
0.5235
0.5100
0.5038
0.5013
0.5004
0.5001
0.5000
0.5000
0.4999
0.4998
0.4992
0.4977
0.4940
0.4859
0.4703
0.4437
0.4028
0.3466
0.2772
0.2000
1.0000
0.9099
0.8244
0.7472
0.6810
0.6272
0.5855
0.5550
0.5337
0.5197
0.5109
0.5055
0.5021
0.4996
0.4971
0.4937
0.4882
0.4798
0.4670
0.4487
0.4237
0.3914
0.3517
0.3054
0.2540
0.2000
1.0000
0.9267
0.8558
0.7896
0.7298
0.6776
0.6335
0.5974
0.5689
0.5470
0.5304
0.5179
0.5083
0.5001
0.4922
0.4835
0.4728
0.4592
0.4418
0.4200
0.3935
0.3621
0.3262
0.2865
0.2440
0.2000
1.0000
0.9366
0.8748
0.8160
0.7614
0.7121
0.6684
0.6307
0.5988
0.5722
0.5502
0.5319
0.5163
0.5023
0.4889
0.4751
0.4600
0.4428
0.4228
0.3997
0.3732
0.3434
0.3105
0.2752
0.2381
0.2000
1.0000
0.9432
0.8876
0.8341
0.7836
0.7369
0.6943
0.6562
0.6226
0.5931
0.5673
0.5447
0.5244
0.5057
0.4879
0.4700
0.4513
0.4313
0.4095
0.3855
0.3592
0.3307
0.3001
0.2677
0.2342
0.2000
un =
241
Contoh 8.6. Distribusi Temperatur sebagai Fungsi Waktu pada Plat Tipis
Plat besi yang sangat luas mempunyai tebal 2 cm. Temperatur mula-mula dalam
plat merupakan fungsi jarak dari salah satu sisinya sebagai berikut :
u = 100x
untuk 0 # x # 1,
u = 100( 2 x)
untuk 0 # x # 1.
Tentukan temperatur tebal plat sebagai fungsi x dan t, jika kedua permukaan
tetap dijaga 0 OC. Untuk besi k = 0,13 kal/detik.cm. OC, c = 0,11 cal/g. OC, =
7,8 g/cm3.
Penyelesaian
kt
c(x )
1 1 ,
=
2 M
Program Matlab
% Persamaan Differensial Parsial
%
Persamaan Parabolik
%
d^2 u
c*rho
d u
%
------- = ----- ----%
d x^2
k
d t
%
Plat tebal 2 cm
%
Kondisi awal
%
u(x,0) = 100x
utk 0 <= x <= 1
%
u(x,0) = 100(2-x)
utk 1 <= x <= 2
%
Kondisi batas
%
u(0,t) = 0
%
u(2,t) = 0
%
Penyelesaian dengan Finite Difference
%
Metode Eksplisit
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
format short
clc
clear all
% Data-data
L=2;
k=0.13;
c=0.11;
rho=7.8;
% Interval
N=8;
M=0.5;
delx=L/N;
%
%
%
%
tebal plat, cm
konduktivitas panas, kal/det.cm.oC
kapasitas panas, kal/g.oC
densitas, g/cm^3
% modulus
242
delt=M*c*rho*delx^2/k;
xo=0;
Jend=10;
for i=1:N+1
x(i)=xo+delx*(i-1);
end
% Kondisi awal
for i=1:ceil((N+1)/2)
u(1,i)=100*x(i);
end
for i=N+1:-1:ceil((N+1)/2)
u(1,i)=100*(2-x(i));
end
for i=1:Jend
t(i)=delt*i;
end
for j=1:Jend
for i=2:ceil((N+1)/2)+1
u(j+1,i)=(k*delt/(c*rho*delx^2))*(u(j,i-1)+u(j,i+1))+...
(1-2*k*delt/(c*rho*delx^2))*u(j,i);
end
for i=N+1:-1:ceil((N+1)/2)+1
u(j+1,i)=u(j+1,ceil((N+1)/2)*2-i);
end
end
x
t=t'
u
Keluaran program
x =
0
0.2500
0.5000
0.7500
1.0000
1.2500
1.5000
1.7500
2.0000
25.0000
25.0000
25.0000
25.0000
21.8750
21.8750
18.7500
18.7500
16.0156
16.0156
13.6719
50.0000
50.0000
50.0000
43.7500
43.7500
37.5000
37.5000
32.0313
32.0313
27.3438
27.3438
75.0000
75.0000
62.5000
62.5000
53.1250
53.1250
45.3125
45.3125
38.6719
38.6719
33.0078
100.000
75.0000
75.0000
62.5000
62.5000
53.1250
53.1250
45.3125
45.3125
38.6719
38.6719
75.0000
75.0000
62.5000
62.5000
53.1250
53.1250
45.3125
45.3125
38.6719
38.6719
33.0078
50.0000
50.0000
50.0000
43.7500
43.7500
37.5000
37.5000
32.0313
32.0313
27.3438
27.3438
25.0000
25.0000
25.0000
25.0000
21.8750
21.8750
18.7500
18.7500
16.0156
16.0156
13.6719
0
0
0
0
0
0
0
0
0
0
0
t =
0.2062
0.4125
0.6187
0.8250
1.0313
1.2375
1.4437
1.6500
1.8563
2.0625
u =
0
0
0
0
0
0
0
0
0
0
0
243
Suatu tabung panjang 20 cm mula-mula berisi udara dengan 2 % uap alkohol.
Pada bagian bawah tabung berhubungan dengan bejana berisi alkohol sehingga
alkohol tersebut menguap melalui tabung yang mula-mula berisi udara diam
tersebut. Pada bagian ini konsentrasi alkohol dijaga tetap 10 %. Pada bagian atas
(puncak) tabung uap alkohol di permukaan atas tabung dapat dianggap selalu
nol.
20 cm
Tentukan distribusi konsentrasi alkohol pada tabung sampai minimal 100 detik.
Diketahui = 0,119 cm2 / detik. Ambil r = t/(x)2 =1/2 dan x = 4 cm.
Penyelesaian
t = 0.5(x)2/ = 67,2 detik
Program Matlab
% Persamaan Differensial Parsial
%
Persamaan Parabolik
%
d^2 c
d c
%
D ------- = ----%
d x^2
d t
%
Kondisi awal
%
c(0,t) = 2
%
Kondisi batas
%
c(0,t) = 0
c(20,t) = 10
%
Penyelesaian dengan Finite Difference
%
Metode Eksplisit
%
%
--------------------------------------------------------------%
Surakarta, Oktober 2005
%
Jurusan Teknik Kimia, Fak. Teknik
%
Universitas Sebelas Maret
%
***************************************************************
%
244
format short
clc
clear all
% Data-data
L=20;
D=0.119;
% panjang tube, cm
% koefisien difusivitas, cm^2/detik
N=5;
M=0.5;
delx=L/N;
delt=M/D*delx^2;
xo=0;
Jend=16;
for i=1:N+1
x(i)=xo+delx*(i-1);
end
x
%Kondisi awal
u(1,1)=0.0;
for i=2:N
u(1,i)=2;
end
u(1,N+1)=10.0;
% interval waktu
for i=1:Jend+1
t(i) = delt*i-delt;
end
t=t'
for j=1:Jend
u(j+1,1)=0.0;
for i=2:N
u(j+1,i)=(delt*D/delx^2)*(u(j,i-1)+u(j,i+1))+...
(1-2*delt*D/delx^2)*u(j,i);
end
u(j+1,N+1)=10.0;
end
u
Keluaran program
x =
0
t =
1.0e+003 *
0
0.0672
0.1345
0.2017
0.2689
0.3361
0.4034
0.4706
0.5378
0.6050
0.6723
0.7395
0.8067
0.8739
0.9412
1.0084
1.0756
12
16
20
245
u =
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
2.0000
1.0000
1.0000
0.7500
1.2500
1.1250
1.5000
1.4219
1.6719
1.6211
1.7852
1.7520
1.8594
1.8376
1.9080
1.8937
1.9398
2.0000
2.0000
1.5000
2.5000
2.2500
3.0000
2.8438
3.3438
3.2422
3.5703
3.5039
3.7188
3.6753
3.8159
3.7875
3.8795
3.8609
2.0000
2.0000
4.0000
3.7500
4.7500
4.5625
5.1875
5.0625
5.4688
5.3867
5.6523
5.5986
5.7725
5.7373
5.8511
5.8281
5.9025
2.0000
6.0000
6.0000
7.0000
6.8750
7.3750
7.2813
7.5938
7.5313
7.7344
7.6934
7.8262
7.7993
7.8862
7.8687
7.9255
7.9140
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
10.0000
heat
246
y=L\b;x=U\y;
u(j,2:nx)=x';b=x;
b(1,1)=b(1,1)+lowb*alpha;
b(nx-1,1)=b(nx-1,1)+hib*alpha;
end
Contoh 8.8. Distribusi Temperatur pada Dinding Batu Tahan Api sebagai
Fungsi Waktu
Dinding batu tahan api setebal 0,3 m mula-mula mempunyai temperatur seragam
100 OC. Temperatur kedua permukaan diturunkan secara tiba-tiba menjadi 20 OC
dan tetap dijaga pada temperatur ini. Tentukan distribusi temperatur pada
dinding pada interval 440 detik sampai 22000 detik. Untuk batu tahan api K =
5.10-7 m2/detik.
heat,
%
%
%
%
%
m^2/detik
ukuran interval x
jumlah interval x
detik, ukuran interval t
jumlah interval t
% Kondisi awal
init = 100*ones(1,nx+1);
% Kondisi batas
lowb=20; hib=20;
% Perhitungan
u=heat(nx,hx,nt,ht,init,lowb,hib,K);
% Plot Hasil
surfl(u)
axis([0 16 0 50 0 120])
view([-217 30])
xlabel('x'); ylabel('waktu');
zlabel('temperatur')
247
Keluaran Program
DAFTAR PUSTAKA
---, 2005, MATLAB Central > File Exchange > Chemistry and Physics,
www.mathworks.com/matlabcentral/fileexchange/load Category.do
Chapra, S.C. dan Canale, R.P., 1991, Metode Numerik untuk Teknik, Penerbit
UI Jakarta
Cutlip, M.B. dan Schacham, M., 2002, Application of Mathematical Software
Packages in Chemical Engineering Education, Workshop Presenters,
ASSE Chem. Eng. Division Summer School, University of Colorado,
Boulder
Davis, M.E., 1965, Numerical Analysis for Chemical Engineering, John Wiley
and Sons., New York
Fogler, H.S., 1999, Elements of Chemical Reaction Engineering, 3 ed.,
Prentice Hall International Ed. Englewood Cliffs, New Jersey.
Gerald, C.F. dan Wheatley, P.O., 1985, Applied Numerical Analysis, AdisonWesley Publising Company, Inc., USA
Hanna, O.T. dan Sandall, O.C., 1995, Computational Methods in Chemical
Engineering, Prentice Hall PTR, New Jersey
Hanselman, D. dan Littlefield, B., 2000, Matlab Bahasa Komputasi Teknik,
Penerbit Andi, Yoyakarta
Jenson, V.G. dan Jeffreys, G.V., 1977, Mathematical Methods in Chemical
Engineering, Academic Press, Inc., London
Levenspiel, O., 1999, Chemical Reaction Engineering, 3 ed., John Wiley and
Sons, New York
Mickley, H.S., Sherwood, T.S., dan Reed, C.E., 1975, Applied Mathematical in
Chemical Engineering, Tata Mc. Graw Hill Pub. Ca., New York
Palm, W.J., 1999, Matlab for Engineering Applications, Mc.Graw-Hill, New
York
252
Penny, J. dan Lindfield, G., 2000, Numerical Methods using Matlab, Prentice
Hall, Inc., New Jersey
Pike, R.W., 1986, Optimization for Engineering System, Van Nostrand
Reinhold Company Inc., New York
Rudd, D.F. dan Watson, C.C., 1968, Strategy of Process Engineering, John
Wiley & Sons, Inc., New York
Sediawan, W.B. dan Prasetya, A., 1977, Pemodelan Matematis dan
Penyelesaian Numeris dalam Teknik Kimia, Penerbit Andi Yogyakarta
Smith, J.M., 1981, Chemical Engineering Kinetics, 3rded., Mc. Graw Hill Inc.,
Tokyo
Wang, N. S., 2002, Computer Methods in Chemical Engineering, Chemical
Engineering Dept., University of Maryland, www.glue.umd.edu/
White, J.R., 2005, Mathematical Methods for Engineering, Chemical &
Nuclear Engineering Dept., University of Massachussetts Lowell,
www.profjrwhite.com/math_methods.htm.
White, J.R., 2005, Problem Solving with Matlab, Chemical & Nuclear
Engineering Dept., University of Massachussetts Lowell,
www.profjrwhite.com/matlab_course.htm.