Tugas Pemodelan Matematis Kelas B
Disusun Oleh :
1.
Atholloh Adkha Falakh Karim
15/378993/TK/42935
2.
Arief Fathur Rahman
15/378991/TK/42933
DEPARTEMEN TEKNIK KIMIA FAKULTAS TEKNIK UNIVERSITAS GADJAH MADA YOGYAKARTA 2017
1. Reaksi Heterogen A ----> Produk terjadi dalam katalisator padat bentuk bola berjari-jari R, dengan laju reaksi: (-rA)=k Ca
.
Luas Total bahan aktif dalam katalis sebesar a (luas/volum). Difusivitas efektif A dalam katalis De. Bila Konsentrasi A fasa gas dipermukaan luar katalis Caf. Laju Transfer massa Gas0Permukaan sbb: .
NA=kc (Caf-Cas)
Tentukan Distribusi Konsentrasi A (Ca) didalam katalis pada keadaan Steady state. (Kasus 9A: Difusi dan Reaksi pada Katalisator Padat bentuk bola)
ANSWER: First step : Mencari Persamaan Pemodelan Matematis
Neraca Massa Katalis di elemen volume:
NA . A| r - NA . A| r +Δr - k Ca=0
-De.4r2
|r +
De.4(r+ )2
(r+ )2 − r2 lim →0
2 (r )
r2
² + ²
2r
-
ka
-
ka
Ca r 2 =0
k a
Ca r 2 =0
|r +Δr -k a Ca 4 ̅ 2 Δr =0
Ca ̅
2
=0
² 2 k a + ²
Ca =0
(1)
dengan Boundary Condition
r=0
=0
r=R
kc a = (Ca-Caf)
Second step: Merubah persamaan ke FDA
(+1 )−() =
()−(−1) =
² ²
=
(forward)
(2)
(backward)
(3)
(+1)−2()+(−1) ²
(4)
sehingga, dengan mensubstitusi 2 dan 3 persamaan diatas dengan persamaan 1, didapat persamaan berikut:
(+1)−2()+(−1) 2 (+1)−() k a + ²
Ca(i) =0
Third Step : Diselesaikan dengan MATLAB
Hasil dan Pembahasan
Reaksi Heterogen A ----> Produk terjadi dalam katalisator padat bentuk bola berjari-jari R, dengan laju reaksi: (-rA)=k Ca
.
Dapat dihitung menggunakan fsolve atau bvp4c, jika menggunakan fsolve harus diubah bentuk menjadi persamaan aljabar menggunakan FDA, sedangkan jika menggunakan bvp4c persaman langsung dapat diselesaikan. (kedua script terlampir). Hasil kedua perhitungan, baik menggunakan fsolve maupun bvp4c menunjukan bahwa semakin besar r (jarak dari pusat katalis) maka Ca (konsentrasi komponen A dalam katalis semakin besar, dapat dilihat di grafik berikut.
Bvp4c merupakan toolbox yang tersedia dalam matlab yang digunakan untuk menyelesaikan persamaan diferensial dengan Boundary value problem. dalam bvp4c, program dibagi menjadi 3 bagian, yaitu file berisi persamaan diferensial, BVP, dan Main program. Bvp4c ini tidak dapat menyelesaikan persamaan diferensial orde 2 secara langsung oleh karena itu, persamaan diferensial orde dua dimisalkan agar bentuknya menjadi orde satu. Sehingga kita memperoleh 1 variabel independent (r) dan 2 variabel dependent ( Persamaan Diferensial orde dua dan Ca). Boundary value problem pada bvp4c ada 2 yakni initial condition dan boundary condition. main program pada bvp4c berisi persamaan differensial, BVP,dan tebakan awal, dapat dituliskan sol=(@fun1,@fun2,solinit), solinit=(linspace(r=0,r=R), [tebakan_awal tebakan awal]). Untuk menekstrak hasil digunakan yint=deval(sol,xinit), dengan xinit=linspace(r=0, r=R). Semisal jika ingin mengekstrak Ca maka Ca=yint(1,:), jika ingin mengekstrak Ca persamaan orde 2, Ca=yint(2,:).
Fungsi fsolve merupakan fungsi yang digunakan untuk menyelesaikan persamaan aljabar, sehingga persamaan diferensial orde dua yang didapat harus diubah dahulu menjadi persamaan aljabar menggunakan metode finite difference approximation (FDA). Metode ini terdiri dari 2 bagian, pertama berisi main program dan yang kedua berisi persamaan aljabar, persamaan aljabar ini terdiri dari BC dan PD.
Kesimpulan
1.Kenaikan nilai jari jari berbanding lurus dengan kenaikan nilai konsentrasi 2. Persamaan diferensial orde dua dapat diselesaikan dengan fungsi fsolve atau toolbox bvp4c
Lampiran
clc clear global Cas De a Rout k alfa Caf kc Cas=.5; Rout=.1; De=.1; a=200; k=.1; alfa=k*a/De; Caf=0.2; kc=.01; solinit=bvpinit(linspace(1e-4,Rout),[Cas 0]); sol=bvp4c(@f12,@f13,solinit);
xint=linspace(1e-4, Rout); yint=deval(sol,xint) %extract solution
Caplot=yint(1,:); xint' Caplot' plot(xint,Caplot); ylabel('Ca'); xlabel('r'); legend('zzz');
function dCadr=f12(r,var) global De a kc Ca=var(1); u=var(2); dCadr(1)=u; dCadr(2)=kc*a*Ca/De-2*u/r; dCadr=dCadr'; end
Error using f12 (line 3) Not enough input arguments. Published with MATLAB® R2012b
function BC=f13(a,b) global kc De Caf Ca1=a(1); u1=a(2); Ca2=b(1); u2=b(2); BC=zeros(2,1); BC(1)=u1; BC(2)=u2+kc/De*(Ca2-Caf); end
Error using f13 (line 3) Not enough input arguments. Published with MATLAB® R2012b
yint =
Columns 1 through 7
0.0255 0.0255 0.0255 0.0255 0.0255 0.0255 0.0255 0 0.0002 0.0004 0.0005 0.0007 0.0009 0.0010
Columns 8 through 14
0.0255
0.0255
0.0255
0.0255
0.0255
0.0255
0.0255
0.0012 0.0014 0.0016 0.0017 0.0019 0.0021 0.0023
Columns 15 through 21
0.0255 0.0255 0.0256 0.0256 0.0256 0.0256 0.0256 0.0024 0.0026 0.0028 0.0029 0.0031 0.0033 0.0035
Columns 22 through 28
0.0256 0.0256 0.0256 0.0256 0.0256 0.0256 0.0256 0.0036 0.0038 0.0040 0.0041 0.0043 0.0045 0.0047
Columns 29 through 35
0.0256 0.0256 0.0256 0.0256 0.0256 0.0256 0.0256 0.0048 0.0050 0.0052 0.0054 0.0055 0.0057 0.0059
Columns 36 through 42
0.0256 0.0256 0.0256 0.0257 0.0257 0.0257 0.0257 0.0060 0.0062 0.0064 0.0066 0.0067 0.0069 0.0071
Columns 43 through 49
0.0257 0.0257 0.0257 0.0257 0.0257 0.0257 0.0257 0.0073 0.0074 0.0076 0.0078 0.0080 0.0081 0.0083
Columns 50 through 56
0.0257 0.0257 0.0258 0.0258 0.0258 0.0258 0.0258 0.0085 0.0086 0.0088 0.0090 0.0092 0.0093 0.0095
Columns 57 through 63
0.0258 0.0258 0.0258 0.0258 0.0258 0.0259 0.0259 0.0097 0.0099 0.0100 0.0102 0.0104 0.0106 0.0107
Columns 64 through 70
0.0259 0.0259 0.0259 0.0259 0.0259 0.0259 0.0259 0.0109 0.0111 0.0113 0.0115 0.0116 0.0118 0.0120
Columns 71 through 77
0.0260 0.0260 0.0260 0.0260 0.0260 0.0260 0.0260 0.0122 0.0123 0.0125 0.0127 0.0129 0.0130 0.0132
Columns 78 through 84
0.0260 0.0261 0.0261 0.0261 0.0261 0.0261 0.0261 0.0134 0.0136 0.0138 0.0139 0.0141 0.0143 0.0145
Columns 85 through 91
0.0261 0.0262 0.0262 0.0262 0.0262 0.0262 0.0262 0.0147 0.0148 0.0150 0.0152 0.0154 0.0155 0.0157
Columns 92 through 98
0.0263 0.0263 0.0263 0.0263 0.0263 0.0263 0.0264 0.0159 0.0161 0.0163 0.0165 0.0166 0.0168 0.0170
Columns 99 through 100
0.0264 0.0264 0.0172 0.0174
ans =
0.0162
0.0344
0.0526
0.0707
0.0173
0.0354
0.0536
0.0717
0.0001
0.0183
0.0364
0.0546
0.0728
0.0011
0.0193
0.0374
0.0556
0.0738
0.0021
0.0203
0.0384
0.0566
0.0748
0.0031
0.0213
0.0395
0.0576
0.0758
0.0041
0.0223
0.0405
0.0586
0.0768
0.0051
0.0233
0.0415
0.0596
0.0778
0.0062
0.0243
0.0425
0.0606
0.0788
0.0072
0.0253
0.0435
0.0617
0.0798
0.0082
0.0263
0.0445
0.0627
0.0808
0.0092
0.0273
0.0455
0.0637
0.0818
0.0102
0.0284
0.0465
0.0647
0.0828
0.0112
0.0294
0.0475
0.0657
0.0839
0.0122
0.0304
0.0485
0.0667
0.0849
0.0132
0.0314
0.0495
0.0677
0.0859
0.0142
0.0324
0.0506
0.0687
0.0869
0.0152
0.0334
0.0516
0.0697
0.0879
0.0889
0.0255
0.0257
0.0260
0.0899
0.0255
0.0257
0.0260
0.0909
0.0255
0.0257
0.0260
0.0919
0.0256
0.0257
0.0260
0.0929
0.0256
0.0257
0.0260
0.0939
0.0256
0.0257
0.0260
0.0950
0.0256
0.0257
0.0260
0.0960
0.0256
0.0257
0.0261
0.0970
0.0256
0.0257
0.0261
0.0980
0.0256
0.0258
0.0261
0.0990
0.0256
0.0258
0.0261
0.1000
0.0256
0.0258
0.0261
0.0256
0.0258
0.0261
0.0256
0.0258
0.0261
0.0256
0.0258
0.0262
0.0256
0.0258
0.0262
0.0255
0.0256
0.0258
0.0262
0.0255
0.0256
0.0258
0.0262
0.0255
0.0256
0.0258
0.0262
0.0255
0.0256
0.0259
0.0262
0.0255
0.0256
0.0259
0.0263
0.0255
0.0256
0.0259
0.0263
0.0255
0.0256
0.0259
0.0263
0.0255
0.0256
0.0259
0.0263
0.0255
0.0256
0.0259
0.0263
0.0255
0.0257
0.0259
0.0263
0.0255
0.0257
0.0259
0.0264
0.0255
0.0257
0.0259
0.0264
0.0255
0.0257
0.0260
0.0264
ans =
Published with MATLAB® R2012b
clc clear all global k kc dp De a CAs n dr Caf %Data yang digunakan untuk perhitungan k kc
=0.1; %Ketetapan reaksi(cm/s) =0.01; %Diameter katalis(cm)
dp De
=0.2; % =1e-2;
a =200; CAs =0.5; Caf =0.2; n =21; dr =dp/(n-1); %Solver CA0(1,1:n)=0.3; CA = fsolve(@odefun9a,CA0,[]); %Plotting
r = linspace(0,dp/2,n); plot(r,CA) title('Profil Konsentrasi A terhadap jari-jari katalisator') xlabel('Jari-Jari Katalisator, cm') ylabel('Konsentrasi A, mol/cm3') legend('Konsentrasi A') %Mencetak hasil perhitungan fprintf('%5.4f \t\t %5.4f \n',[r;CA]) function fy=odefun9a(CA) global k kc dp De a CAs n dr Caf fy=zeros(n,1); r = linspace(0,dp/2,n); for i=1 fy(i)= (CA(i+1)-CA(i))/dr; end for i=2:n-1 fy(i)= (CA(i+1)-2*CA(i)+CA(i-1))/dr^2 +... 2/r(i)*(CA(i+1)-CA(i-1))/(2*dr) - ... k*a/De*CA(i); end for i=n fy(i)=De*(CA(i)-CA(i-1))/dr + kc*(Caf-CAs); end end
Error using odefun9a (line 7) Not enough input arguments. Published with MATLAB® R2012b
Equation solved. fsolve completed because the vector of function values is near zero as measured by the default value of the function tolerance, and the problem appears regular as measured by the gradient.
0.0000 0.0050 0.0100
0.0001 0.0001 0.0001
0.0150 0.0200
0.0001 0.0001
0.0250 0.0300
0.0001 0.0002
0.0350 0.0400 0.0450
0.0002 0.0003 0.0003
0.0500 0.0550
0.0004 0.0006
0.0600 0.0650 0.0700
0.0007 0.0010 0.0014
0.0750 0.0800
0.0019 0.0026
0.0850 0.0900
0.0036 0.0051
0.0950 0.1000
0.0072 0.0102
Published with MATLAB® R2012b
