Nonisothermal Plug Flow Reactor
کد:
% rate1.m
% This function gives the right-hand side for a simple
% reactor problem which is isothermal.
function ydot=rate1(VR,y)
% y(1) is CA, y(2) is CB, y(3) is CC
% k=0.3 and u=0.5
CA=y(1);
rate=0.3*CA*CA;
ydot(1)=-2.*rate/0.5;
ydot(2) =+rate/0.5;
ydot(3)=0.;
ydot=ydot’;
plot(z,y(:,1),‘*-’,z,y(:,2),‘+-’,z,y(:,3),‘x-’)
% rateSO2.m for SO2 reaction
function ydot=rateSO2(z,y)
% X is concentration SO2 divided by the inlet concentration
X=y(1);
% T is temperature in degrees K
T=y(2);
k1=exp(-14.96+11070/T);
k2=exp(-1.331+2331/T);
Keq=exp(-11.02+11570/T);
term1=X*sqrt(1-0.167*(1-X));
term2=2.2*(1-X)/Keq;
denom=(k1+k2*(1-X))^2;
rate=(term1-term2)/denom;
ydot(1)=-50*rate;
ydot(2)=-4.1*(T-673.2)+1.02e4*rate;
ydot=ydot’;
y(1)=0.2; y(2)=573.2;
rateSO2(0.1,y)
ans=-0.001148 410.235
% run_SO2.m
% set the dimensionless initial conditions
y0=[1 673.2]
% set the integration range
zspan=[0 1]
% call the solver
[z y]=ode45(@rateSO2,zspan,y0)
% plot the result
plot(z,y(:,1),‘r-’)
xlabel(‘dimensionless axial position’)
legend(‘fraction converted’)
pause
plot(z,y(:,2),‘g-’)
xlabel(‘dimensionless axial position’)
legend(‘temperature (K)’)