Case 1
Case 1
Case 1
Disusun Oleh :
Himmah Sekar Eka Ayu Gustiana
(16/405832/PTK/11287)
Program Pascasarjana
Jurusan Teknik Kimia
Fakultas Teknik
Universitas Gadjah Mada
Yogyakarta
CASE 1
Unsteady state temperature distribution in a cylindrical rod in the air. A cylinder rod has a length of L and a
diameter D. One of its end is attaching to a wall having temperature of Ts while air temperature around is
constant at Ta. Diameter of the rod is much smaller than the length so that radial temperature gradient is
negligible coefficient of heat transfer between the rod surface to the air is h. The rod is made of material having
heat capacity of Cp, density and heat conductivity k. Please develope mathematical equations to determine
temperature distribution in the rod as a function of position time (unsteady state).
Method of line :
2 1 .
2
. . ( ) =
x=0 x=L
2 +1 2 + 1
=
2 2
+1 2 + 1 1
= [ . . ( )] ;2
. 2
At i=1 : 1 =
At i=N+1
. . 2 |= = . . 2 (|= )
4 4
| = (|= )
=
3+1 4 + 1
= (+1 )
2
2
4 1 +
+1 =
2
3+
for i=2:N
T0(i)=Ta;
end
time_span=[0 time];
[time,T] = ode15s(@ode_case1,time_span,T0,[],Ts,dx,h,k,Ta,rho,Cp,N,d);
%membuat grafik
for i=1:N+1
x(i)=dx*(i-1)
end
surf(time,x,T')
xlabel('waktu,detik')
ylabel('posisi,cm')
zlabel('suhu,C')
title('Grafik antara waktu,posisi dan suhu')
function ode=ode_case1(time,T,Ts,dx,h,k,Ta,rho,Cp,N,d)
%nilai T1
T(1)=Ts;
%nilai T(N+1)
T(N+1)=(4*T(N)-T(N-1)+(2*dx*h/k)*Ta)/(3+(2*dx*h/k));
for i=2:N
x(i)=dx*(i-1);
alpa=k/rho/Cp;
d2Tdx2(i)=(T(i+1)-2*T(i)+T(i-1))/(dx^2);
ode(i)=alpa*(d2Tdx2(i)-(h/k/d*(T(i)-Ta)));
end
ode=ode';
Run Program :
Grafik antara waktu,posisi dan suhu
400
300
200
suhu,C
100
10
20
30 15
waktu,detik 5 10
0
posisi,cm
CASE 2
Liquid containing component A with a concentration of CAin is flowing through a cylindrical pipe having a radius of
B and a length of L. The flow is in laminar regime which velocity following
2
= 2 (1 ( )
From a certain position in a length direction, the pipe material consists of dissolvable component B. In the liquid
body, the component A react with B to form a product following first order reaction which responds to A and B.
The liquid near containing B pipe surface is always saturated with B. Please determine steady state concentration
of A and B in the following liquid.
Method of Line:
2 1
+ =0
2
2 1
+ =0
2
2 1
= [ + ]
2
2 1
= [ 2 + ]
r=0 r=R
+1 1
=
2
2 +1 22 + 1
=
2 2
For i=1
=0
=0
3 +42 3 +42
1 = 3
1 = 3
For i=N+1
(1) +4
=0 (+1) = 3
(+1) =
for i=2:N
Ca(:,i)=C(:,i);
Cb(:,i)=C(:,N+i);
end
Ca(:,1)=(-Ca(:,3)+4*Ca(:,2))/3;
Cb(:,1)=(-Cb(:,3)+4*Cb(:,2))/3;
Ca(:,N+1)=(-Ca(:,N-1)+4*Ca(:,N))/3;
Cb(:,N+1)=Cbs;
for i=1:N+1
r(i)=(i-1)*dr;
end
function ode=ode_case2(z,C,N,Cbs,dr,Vavg,R,Da,k,Db)
for i=2:N
Ca(i)=C(i);
Cb(i)=C(N+i);
end
%nilai Ca1 dan Cb1
Ca(1)=(-Ca(3)+4*Ca(2))/3;
Cb(1)=(-Cb(3)+4*Cb(2))/3;
for i=2:N
r(i)=dr*(i-1);
Vz(i)=2*Vavg*(1-((r(i)/R)^2));
dCadr(i)=(Ca(i+1)-Ca(i-1))/2/dr;
d2Cadr2(i)=(Ca(i+1)-2*Ca(i)+Ca(i-1))/(dr^2);
odeA(i)=Da/Vz(i)*(d2Cadr2(i)+dCadr(i)/r(i)-k*Ca(i)*Cb(i)/Da);
dCbdr(i)=(Cb(i+1)-Cb(i-1))/2/dr;
d2Cbdr2(i)=(Cb(i+1)-2*Cb(i)+Cb(i-1))/(dr^2);
odeB(i)=Db/Vz(i)*(d2Cbdr2(i)+dCbdr(i)/r(i)-k*Ca(i)*Cb(i)/Db);
end
ode=[odeA odeB]';
Run Program :
Ca profile
0.12
0.1
3
0.08
Ca,mole/cm
0.06
0.04
0.02
2
1.5 100
80
1
60
0.5 40
20
0 0
r position,cm
z position,cm
Cb profile
0.15
0.1
3
Cb,mole/cm
0.05
-0.05
2
1.5 100
80
1
60
0.5 40
20
0 0
r position,cm
z position,cm
CASE 3
A solution containing A is flowing down on a vertical solid surface due to gravitational effect that form a film with
a thickness of . The surface contains component B that is soluble in the solution. In the solution A reacts with B,
A+Bproduct with a rate of reaction (amount of A is assumed to be unlimited): -rA=k.CB (mole/time/volume)
The solution near the surface is always in saturation of B (CBS). Please determine the concentration distribution of
B in the solution at unsteady state condition. The flow pattern of the solution in y direction is following laminar
flow (x is distance from the solid surface)
2
= ( )
2
Method of Line :
2
= [ 2 ]
In x direction x=0 x=
(,) (,)
2 (,) (+1,) 2(,) + (1,)
=
2 2
(,) (,+1) (,1)
=
2
(,) (+1,) 2(,) + (1,) (,+1) (,1)
= [ ( ) ]
2 2 (,)
2 ; 2
Nx=20; %increment of x
Ny=20; %increment of y
dx=delta/Nx;
dy=L/Ny;
% initital condition
Cb0=0;
%karena punya (Nx-1)*(Ny-1) ODES jadi
for j=2:Ny
for i=2:Nx
Cbcal0(i,j)=Cb0;
end
end
% penyelesaian ODE
tspan=[0 100];
[t,Cbcalc]=ode15s(@ode_case3,tspan,Cb0calc,[],Ny,Nx,Cbs,Cbin,dx,rho,g,delta,miu,dy,Db
,k);
size(t)
size(Cbcalc)
ij=1;
for j=1:Ny+1
y(ij)=(j-1)*dy;
ij=ij+1;
end
ij=1;
for i=1:Nx+1
x(ij)=(j-1)*dy;
ij=ij+1;
end
size(x)
size(y)
Cbhit(1,2:Nx+1)=Cbin;
Cbhit(2:Ny+1,1)=Cbs;
ij=1;
%for k=1:length(t)
k=2;
for i=2:Ny
for j=2:Nx
Cbhit(i,j)=Cbcalc(k,ij)
ij=ij+1;
end
end
size(Cbhit)
%Definisi nilai IC untuk plot grafik
surf(x,y,Cbhit)
xlabel('posisi x')
ylabel('posisi y')
zlabel('Cb')
function ode=ode_case3(t,Cbcalc,Ny,Nx,Cbs,Cbin,dx,rho,g,delta,miu,dy,Db,k)
ij=1;
%convert 1 coloumn dimension into 2 coloumn dimension
for j=2:Ny
for i=2:Nx
Cb(i,j)=Cbcalc(ij);
ij=ij+1;
end
end
for i=2:Nx
Cb(i,1)=Cbin; %at j=1
Cb(i,Ny+1)=(4*Cb(i,Ny)-Cb(i,Ny-1))/3; %at j=Ny+1
end
ij=1;
for j=2:Ny
for i=2:Nx
x=(i-1)*dx;
vy=(rho*g*delta/miu)*(x-x^2/(2*delta));
ddCbdx=(Cb(i+1,j)-2*Cb(i,j)+Cb(i-1,j))/dx^2;
dCbdy=(Cb(i,j+1)-Cb(i,j-1))/(2*dy);
ode(ij)=Db*ddCbdx-vy*dCbdy-k*Cb(i,j);
ij=ij+1;
end
end
ode=ode';
Run Program :
0.1
0.08
0.06
Cb
0.04
0.02
0
60
40
20 51
50.5
50
49.5
0
posisi y 49
posisi x
CASE 4
Heat Conduction and Mass Diffusion in Slab Catalyst
Method of Line :
2 1
( . )
= ( 2 ) ; =
2
= ( 2 + ( ) ) ; =
.
x=0 x=L
At point i
(+1) 2() + (1)
= (( ) )
2
(+1) 2() + (1)
= (( 2 ) + ( ) ) ;2
At point i=1
= 0 ()
= 0 ()
42 3
1 =
3
42 3
1 =
3
At point i=N+1
CN+1=CAS
TN+1=Tf
clear all
clc
global De delta eps rho Cp deltaHr k A E R T0 Tf CAs N dx
dx = delta/N;
% Initial Condition
for i=2:N
CA0(i)=0;
T0(i)=303;
end
D0=[CA0 T0]
tspan=[0 100]
for i=2:N
CA(:,i)=D(:,i);
T(:,i)=D(:,N+i);
end
CA(:,1)=(-CA(:,3)+4*CA(:,2))/3;
CA(:,N+1)=CAs;
T(:,1)=(-T(:,3)+4*T(:,2))/3;
T(:,N+1)=Tf;
for i=1:N+1
x(i)=(i-1)*dx
end
figure (1)
surf (t,x,CA')
title ('CA Profile')
xlabel ('time,minute')
ylabel ('thickness,cm')
zlabel ('concentration of A,mole/cm^3')
figure (2)
surf (t,x,T')
title ('T Profile')
xlabel ('time,minute')
ylabel ('thickness,cm')
zlabel ('temperature,K')
function ode=ode_func4(t,D)
global De delta eps rho Cp deltaHr k A E R T0 Tf CAs N dx
for i=2:N
CA(i)=D(i);
T(i)=D(N+i);
end
CA(1)=(-CA(3)+4*CA(2))/3;
CA(N+1)=CAs;
T(1)=(-T(3)+4*T(2))/3;
T(N+1)=Tf;
for i=2:N
kr(i) = A*exp(-E/R/T(i));
x(i)=(i-1)*dx;
ddCAdx2(i)=(CA(i+1)-2*CA(i)+CA(i-1))/dx^2;
odeCA1(i)=ddCAdx2(i)-kr(i)/De*CA(i);
odeCA(i)=De/eps*odeCA1(i);
ddTdx2(i)=(T(i+1)-2*T(i)+T(i-1))/dx^2;
odeT1(i)=ddTdx2(i)-kr(i)/k*(-deltaHr)*CA(i);
odeT(i)=k/(rho*Cp)*odeT1(i);
end
ode=[odeCA odeT]';
Run Program :
CA Profile
0.6
0.5
3
0.4
concentration of A,mole/cm
0.3
0.2
0.1
0
100
-0.1
50
0.5
0.4
0.3
0.2
0.1 0
0
time,minute
thickness,cm
T Profile
1600
1400
1200
temperature,K
1000
800
600
0.8
400
0.6
200 0.4
0
20 0.2
40
60
80 0
100 thickness,cm
time,minute