Why my graph is not the same as paper!!

2 views (last 30 days)
clear all
rlength =10; %meter
dx = 0.0001;
Numdata = rlength/dx;
%constant value
Cpf=0.707; %kcal/kg.K
Cpg=0.719; %kcal/kg.K
f=1.0;
H=-26000; %kcal/kg.mol
R=1.987; %kcal/kg.mol.K
S1=10; %meter
S2=0.78; %square meter
U=500; %kcal/h.m^2.K
W=26400; %kg/h
N0=701.2; %kg.mol/m^2.h
%initial Temperature and mass flow rate of nitrogen
Tf(1)=650;
Tg(1)=650;
N2(1)=700;
%Find k values
k1=1.78954e4*exp(-20800/(R*Tg));
k2=2.5714e16*exp(-47400/(R*Tg));
%find partial pressure
pN2=286*N2/(2.598*N0 + 2*N2);
pH2=3*pN2;
pNH3=286*(2.23*N0-2*N2)/(2.598*N0 + 2*N2);
for n=1:1:Numdata
%Process
Tf(n+1)=Tf(n)+dx*(-((U*S1)/(W*Cpf))*(Tg(n)-Tf(n)));
Tg(n+1)=Tg(n)+dx*(-((U*S1)/(W*Cpg))*(Tg(n)-Tf(n))+((-H*S2)/(W*Cpg))*(f*(k1*pN2*(pH2^1.5)/pNH3 - k2*pNH3/(pH2^1.5))));
N2(n+1)=N2(n)+dx*(-f*(k1*pN2*(pH2^1.5)/pNH3 - k2*pNH3/(pH2^1.5)));
end;
for n = 1:1:Numdata+1
%convert length
xi(n) = (n-1)*dx;
end
%show figure
figure
plot(xi,Tf,'-r'),...
ylabel('Temperature(K)'),title('Feed gas'),grid
xlabel('Reactor Length (m)')
axis([0 10 100 900])
figure
plot(xi,Tg,'-b'),...
ylabel('Temperature(K)'),title('Product gas'),grid
xlabel('Reactor Length(m)')
axis([0 10 100 900])
figure
plot(xi,N2,'-g'),...
ylabel('Temperature(K)'),title('flow rate of nitrogen'),grid
xlabel('Reactor Length(m)')
axis([0 10 100 900])
from paper
  3 Comments
Torsten
Torsten on 12 Jan 2022
k values and partial pressures have to be recomputed during the course of integration because Tf, Tg and N2 change.

Sign in to comment.

Accepted Answer

Torsten
Torsten on 12 Jan 2022
clear all
rlength =10; %meter
dx = 0.0001;
Numdata = rlength/dx;
%constant value
Cpf=0.707; %kcal/kg.K
Cpg=0.719; %kcal/kg.K
f=1.0;
H=-26000; %kcal/kg.mol
R=1.987; %kcal/kg.mol.K
S1=10; %meter
S2=0.78; %square meter
U=500; %kcal/h.m^2.K
W=26400; %kg/h
N0=701.2; %kg.mol/m^2.h
%initial Temperature and mass flow rate of nitrogen
Tf(1)=650;
Tg(1)=650;
N2(1)=700;
for n=1:1:Numdata
%Find k values
k1=1.78954e4*exp(-20800/(R*Tg(n)));
k2=2.5714e16*exp(-47400/(R*Tg(n)));
%find partial pressure
pN2=286*N2(n)/(2.598*N0 + 2*N2(n));
pH2=3*pN2;
pNH3=286*(2.23*N0-2*N2(n))/(2.598*N0 + 2*N2(n));
%Process
Tf(n+1)=Tf(n)+dx*(-((U*S1)/(W*Cpf))*(Tg(n)-Tf(n)));
Tg(n+1)=Tg(n)+dx*(-((U*S1)/(W*Cpg))*(Tg(n)-Tf(n))+((-H*S2)/(W*Cpg))*(f*(k1*pN2*(pH2^1.5)/pNH3 - k2*pNH3/(pH2^1.5))));
N2(n+1)=N2(n)+dx*(-f*(k1*pN2*(pH2^1.5)/pNH3 - k2*pNH3/(pH2^1.5)));
end;
for n = 1:1:Numdata+1
%convert length
xi(n) = (n-1)*dx;
end
%show figure
figure
plot(xi,Tf,'-r'),...
ylabel('Temperature(K)'),title('Feed gas'),grid
xlabel('Reactor Length (m)')
axis([0 10 100 900])
figure
plot(xi,Tg,'-b'),...
ylabel('Temperature(K)'),title('Product gas'),grid
xlabel('Reactor Length(m)')
axis([0 10 100 900])
figure
plot(xi,N2,'-g'),...
ylabel('Temperature(K)'),title('flow rate of nitrogen'),grid
xlabel('Reactor Length(m)')
axis([0 10 100 900])

More Answers (0)

Categories

Find more on 2-D and 3-D Plots in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!