Chemical Reaction Engineering Assignment

Matlab function for Runge Kutta method, Chemical Reaction Engineering

In this sample matlab assignment the subject matter expert is demonstrating to develop a Matlab function to solve problems using 4th order Runge Kutta method. The expert has also tested the implementation of this function in two problems. Both the problems are from chemical reaction engineering domain that involves finding the steady state concentration and use of mass balance principles, etc

SOLUTION : –

clear all

clc

y0=0;

V=10;

F=200;

k=0.1;

F_xy=@(t,c) 200/10-1/10-0.1*c^2;

h=0.1;  %1/10 minute

ntmax=5/h;  %5 minutes

y(1)=y0;

t=0:h:5;

for i=1:ntmax

k1 = h*F_xy(t(i),y(i));

k2 =  h*F_xy(t(i)+h/2,y(i)+k1/2);

k3 = h*F_xy(t(i)+h/2,y(i)+k2/2);

k4 = h*F_xy(t(i)+h,y(i)+k3) ;

y(i+1) = y(i) + (1/6)*(k1+2*k2+2*k3+k4);  % main equation

end

[tt,yy]=ode45(F_xy,t,y0)

plot(t,y,’o’,’Linewidth’,3)

grid on

hold on

plot(t,y,’k’,’Linewidth’,3)

xlabel(‘time (minutes)’)

ylabel(‘concentration g/m^3’)

legend(‘4th order Runge-Kutta’,’ODE45′)

print(‘File1R_1′,’-dpng’)

clear all

clc

c1(1)=0;

c2(1)=0;

c3(1)=0;

c4(1)=0;

cd1(1)=0;

cd2(1)=0;

cd3(1)=0;

cd4(1)=0;

ca0=10;

tau=5;

k=0.1;

h=0.01;

ntmax=20/h;  %

t=0:h:20;

f1=@(t,c1,c2,c3,c4) 1/tau*(ca0-c1)-k*c1;

f2=@(t,c1,c2,c3,c4)  -1/tau*(c2)+k*c1;

f3=@(t,c1,c2,c3,c4)  1/tau*(c1-c3)-k*c3;

f4=@(t,c1,c2,c3,c4) 1/tau*(c2-c4)+k*c3;

for i=1:ntmax

k1=h*f1(t(i),c1(i),c2(i),c3(i),c4(i));

l1=h*f2(t(i),c1(i),c2(i),c3(i),c4(i));

m1=h*f3(t(i),c1(i),c2(i),c3(i),c4(i));

n1=h*f4(t(i),c1(i),c2(i),c3(i),c4(i));

k2=h*f1(t(i)+h/2,c1(i)+k1/2,c2(i)+l1/2,c3(i)+m1/2,c4(i)+n1/2);

l2=h*f2(t(i)+h/2,c1(i)+k1/2,c2(i)+l1/2,c3(i)+m1/2,c4(i)+n1/2);

m2=h*f3(t(i)+h/2,c1(i)+k1/2,c2(i)+l1/2,c3(i)+m1/2,c4(i)+n1/2);

n2=h*f4(t(i)+h/2,c1(i)+k1/2,c2(i)+l1/2,c3(i)+m1/2,c4(i)+n1/2);

k3=h*f1(t(i)+h/2,c1(i)+k2/2,c2(i)+l2/2,c3(i)+m2/2,c4(i)+n2/2);

l3=h*f2(t(i)+h/2,c1(i)+k2/2,c2(i)+l2/2,c3(i)+m2/2,c4(i)+n2/2);

m3=h*f3(t(i)+h/2,c1(i)+k2/2,c2(i)+l2/2,c3(i)+m2/2,c4(i)+n2/2);

n3=h*f4(t(i)+h/2,c1(i)+k2/2,c2(i)+l2/2,c3(i)+m2/2,c4(i)+n2/2);

k4=h*f1(t(i)+h,c1(i)+k3,c2(i)+l3,c3(i)+m3,c4(i)+n3);

l4=h*f2(t(i)+h,c1(i)+k3,c2(i)+l3,c3(i)+m3,c4(i)+n3);

m4=h*f3(t(i)+h,c1(i)+k3,c2(i)+l3,c3(i)+m3,c4(i)+n3);

n4=h*f4(t(i)+h,c1(i)+k3,c2(i)+l3,c3(i)+m3,c4(i)+n3);

c1(i+1)=c1(i)+(1/6)*(k1+2*k2+2*k3+k4);

c2(i+1)=c2(i)+(1/6)*(l1+2*l2+2*l3+l4);

c3(i+1)=c3(i)+(1/6)*(m1+2*m2+2*m3+m4);

c4(i+1)=c4(i)+(1/6)*(n1+2*n2+2*n3+n4);

cd1(i+1)=cd1(i)+h*(1/tau*(ca0-cd1(i))-k*cd1(i));

cd2(i+1)=cd2(i)+h*(-1/tau*(cd2(i))+k*cd1(i));

cd3(i+1)=cd3(i)+h*(1/tau*(cd1(i)-cd3(i))-k*cd3(i));

cd4(i+1)=cd4(i)+h*(1/tau*(cd2(i)-cd4(i))+k*cd3(i));

end

%[tt,cc]=ode45(@rhsFile1R2,[0 20], [0 0 0 0] );

subplot(2,2,1)

plot(t,cd1,’Linewidth’,3)

grid on

xlabel(‘time (minutes)’)

ylabel(‘concentration mg/L’)

%hold on

%plot(tt,cc(:,1),’or’)

%legend(‘R-K 4th order’,’ODE45′)

title(‘Conc. at A outlet’)

subplot(2,2,2)

plot(t,cd2,’b’,’Linewidth’,3)

grid on

%hold on

%plot(tt,cc(:,2),’or’)

%legend(‘R-K 4th order’,’ODE45′)

title(‘Conc. at B outlet’)

xlabel(‘time (minutes)’)

ylabel(‘concentration mg/L’)

subplot(2,2,3)

plot(t,cd3,’b’,’Linewidth’,3)

grid on

%hold on

%plot(tt,cc(:,3),’or’)

%legend(‘R-K 4th order’,’ODE45′)

title(‘Conc. in A ‘)

xlabel(‘time (minutes)’)

ylabel(‘concentration mg/L’)

subplot(2,2,4)

plot(t,cd4,’b’,’Linewidth’,3)

grid on

%hold on

%plot(tt,cc(:,4),’or’)

%legend(‘R-K 4th order’,’ODE45′)

title(‘Conc. in B ‘)

xlabel(‘time (minutes)’)

ylabel(‘concentration mg/L’)

print(‘File1R_2′,’-dpng’)