Matlab code

File: solve_proj.m

% filename solve_proj.m

%------------------------------------------------------------------

% ********************************************

% * Oxidization of Ammonia to Nitric Oxide *

% ********************************************

% f(1:6)= Flow rate of: 1-NH3 2-O2 3-NO 4-H2O 5-N2 6-H2

% f(7)=pressure f(8)=temperature

%------------------------------------------------------------------

clear all

global H cpT P

format long

%

% Reactor Parameters, in "cm"

Length=1.;

Height=4.e-3;

Width=3.e-2;

P=2*(Height+Width); % perimeter

A=Height*Width; % area

%

% Inlet Information

Q=100.; % Volume flow rate, "L/hour"

T0=673; % temperature, "K"

p0=[0.046 0.069 0 0 0 0]; % pressure, "torr"

p0T=sum(p0);

f0T=p0T*101325/760 * Q/(1000*3600) / (8.3145*T0); % mole flow rate, "mol/sec"

%

% Heat Capacity & Reaction Heat

cp(1:6)=[8.38 7.016 7.133 8.025 6.961 6.889]; % "cal/(K.mol)"

dHf(1:6)=[-11.02 0 21.57 -57.796 0 0]; % "kcal/mol"

dHf=dHf*1000;

cpT=0;

for i=1:6

cpT=cpT+p0(i)*cp(i)/p0T;

end

H(1)=dHf(3)+1.5*dHf(4)-dHf(1)-1.25*dHf(2);

H(2)=1.25*dHf(5)+1.5*dHf(4)-dHf(1)-1.5*dHf(3);

H(3)=0.5*dHf(5)+1.5*dHf(6)-dHf(1);

%

% Calculating

for i=1:6

f0(i)=f0T*p0(i)/p0T;

end

f0(7)=p0T;

f0(8)=T0;

options=odeset('AbsTol',[1e-15 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-10]);

[z,f]=ode45('rhs_proj',[0 Length], f0, options);

%

% Plotting

% If index=1; plot flow rate (all components)~z

% `` 2; `` flow rate (NH3 & NO only )~z

% `` 3; `` T ~ z

% `` else; `` p ~ z

index=2;

if index==1

plot(z,f(:,1),'k-',z,f(:,2),'m:',z,f(:,3),'c-.',z,f(:,4),'r--',z,f(:,5),'y-',z,f(:,6),'b-.','LineWidth',1)

title('Flow rate vs. z')

xlabel('z (cm)')

ylabel('flow rate (mol/sec)')

legend('NH3','O2','NO','H2O','N2','H2')

else if index==2

plot(z,f(:,1),'m-',z,f(:,3),'b:','LineWidth',1)

title('Flow rate vs. z')

xlabel('z (cm)')

ylabel('flow rate (mol/sec)')

legend('NH3','NO')

else if index==3

plot(z,f(:,8),'LineWidth',2)

xlabel('z (cm)')

ylabel('T (K)')

title('T ~ z')

else

plot(z,f(:,7),'LineWidth',2)

xlabel('z (cm)')

ylabel('total p (torr)')

title('Total p ~ z')

end

end

end

%

% Conversion, NO yield & Selectivity

m=size(z);

mm=m(1);

format compact

conversion=1.-f(mm,1)/f0(1)

selectivity=4/5*(f0(2)-f(mm,2))/(f0(1)-f(mm,1))

yield_NO=f(mm,3)/f0(1) % in "mol/(mol NH3 in)"

Texit=f(mm,8)

File: rhs_proj.m

% filename rhs_proj

function fdot=rhs_proj(z,f)

% 1-NH3 2-O2 3-NO 4-H2O 5-N2 6-H2 7-pressure 8-temperature

global H cpT P

R=1.987;

pT=f(7);

T=f(8);

fT=0;

for i=1:6

fT=fT+f(i);

end

for i=1:6

y(i)=f(i)/fT;

p(i)=pT*y(i);

end

RT=R*T;

%

% Reaction rate

r11=3.4e-8*exp(21700/RT)*p(1)*p(2)^0.5;

r12=1+0.08*exp(4400/RT)*p(2)^0.5;

r13=1+1.6e-3*exp(25500/RT)*p(1);

r(1)=r11/(r12*r13);

r21=3.36e-9*exp(26000/RT)*p(1)*p(3);

r22=1+5.4e-2*exp(10900/RT)*p(2)+4.5e-6*exp(20900/RT)*p(1);

r(2)=r21/r22^2;

r31=4.4e-6*exp(-570/RT)*p(1);

r(3)=r31/r22;

%r(2)=0; % use this if neglecting 2nd reaction

%r(3)=0; % use this if neglecting 3rd reaction

HT=0;

for i=1:3

HT=HT-H(i)*r(i);

end

%

% Rate of change for f

fdot(1)=P*(-r(1)-r(2)-r(3));

fdot(2)=P*(-1.25*r(1));

fdot(3)=P*(r(1)-1.5*r(2));

fdot(4)=P*(1.5*r(1)+1.5*r(2));

fdot(5)=P*(1.25*r(2)+0.5*r(3));

fdot(6)=P*(1.5*r(3));

fdot(7)=-0.115/2;

fdot(8)=P*HT/cpT/fT;

%fdot(7)=0; % use this if neglecting pressure change

%fdot(8)=0; % use this if neglecting temperature change

fdot=fdot';