Why my graph is not the same as paper!!

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])

[Nattanit Benjangvisanu]

it's from my code.Not same as paper. Help me please! I think because N2 equation and I try to do everything but they are not the same.

[Torsten]

k values and partial pressures have to be recomputed during the course of integration because Tf, Tg and N2 change.

[Nattanit Benjangvisanu]

Sorry,I'm beginner. How to write?

Answer 1

[Torsten]

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])

[Nattanit Benjangvisanu]

Thank you so much!! You're amazing.