% ammonia PBR simulation
% basic data
global Rgas Rgas2;
Rgas = 8.314; % kJ/kgmol/K
Rgas2 = 0.082057; % atm*m3/kgmol/K
% molecular weights
MWN2 = 28.0134; % kg/kgmol
MWH2 = 2.016;
MWNH3 = 17.02;
MWAr = 39.948;
MWCH4 = 16.043;
% heat capacity coefficients from fit of Hysys properties at 150 atm
CpCoef = [ 40.442   -35.279      46.918     -19.41      0 ;
           28.75      1.86        0           0         0  ;
           1088.5    -5691.4     11787      -10900     3801.6 ;
           36.819     -53.01      63.973     -27.022    0 ;
           24.709     50.297       0           0        0 ];
% reaction kinetics from Hysys case
% forward reaction
k0f = 3.6e7; % kgmol/m3/h/atm^2
Ef = 9.1e4; % kJ/kgmol
% reverse reaction
k0r = 4.68e13; % kgmol/m3/h/atm
Er = 1.41e5; % kJ/kgmol
% reactor parameters
Dr = 3; % diameter, m
Lr = 1; % length, m
Ar = pi*Dr^2/4; % x-sectional area, m2
Vr = Ar*Lr; % volume, m3
% catalyst particles and packing
Dp = 1e-3; % particle diameter, m
eps = 0.4; % void fraction
% gas kinematic viscosity estimate from Hysys case
nu = 5.075e-7; % m2/s
% feed conditions from Hysys case
FeedN2 = 12348; % kgmol/h
FeedH2 = 37044;
FeedNH3 = 0;
FeedAr = 12391;
FeedCH4 = 5652;
FeedP = 150; % atm
FeedT = 270; % degC
% mass flux
M0 = FeedN2*MWN2 + FeedH2*MWH2 + FeedNH3*MWNH3 + FeedAr*MWAr + FeedCH4*MWCH4;
G0 = M0/Ar/3600;  % kg/s/m2
% initial conditions
T0 = FeedT;
H0 = (FeedN2*HtCap(1,T0,CpCoef)+FeedH2*HtCap(2,T0,CpCoef) ...
    +FeedNH3*HtCap(3,T0,CpCoef)+FeedAr*HtCap(4,T0,CpCoef) ...
    +FeedCH4*HtCap(5,T0,CpCoef))*T0;
y0 = [ FeedN2 ; FeedH2 ; FeedNH3 ; FeedP ; FeedT ; H0 ];
% mass matrix
M = [ 1 0 0 0 0 0 ;
      0 0 0 0 0 0 ;
      0 0 0 0 0 0 ;
      0 0 0 1 0 0 ;
      0 0 0 0 0 0 ;
      0 0 0 0 0 1 ] ;
% solution span
vspan = linspace(0,Vr,200);
% anonymous function
PFRanon = @(v,y) PBRsysFullModel(v,y,FeedN2,FeedH2,FeedAr,FeedCH4 ... 
    ,k0f,Ef,k0r,Er,eps,nu,G0,Dp,Ar,MWN2,MWH2,MWNH3,MWAr,MWCH4,CpCoef);
options = odeset('Mass',M);
[ v , ysoln ] = ode15s(PFRanon,vspan,y0,options);
% unpack solution
MflowN2 = ysoln(:,1);
MflowH2 = ysoln(:,2);
MflowNH3 = ysoln(:,3);
P = ysoln(:,4);
T = ysoln(:,5);
H = ysoln(:,6);
Conv = (FeedN2 - MflowN2)/FeedN2 *100; % Conversion, N2 basis, in %
% create plots
figure(1) % Temperature Profile
plot(v,T,'k-')
grid
xlabel('Reactor Volume - m3')
ylabel('Temperature - degC')
title('Temperature Profile')
figure(2) % Pressure Profile
plot(v,P,'m-')
grid
xlabel('Reactor Volume - m3')
ylabel('Presssure - atm')
title('Pressure Profile')
figure(3) % Molar Flows Profile
plot(v,MflowN2,'b-',v,MflowH2,'m-',v,MflowNH3,'g-')
grid
xlabel('Reactor Volume - m3')
ylabel('Molar Flow Rate - kgmol/h')
title('Molar Flow Profile')
legend('N2','H2','NH3','Location','northeastoutside')
figure(4)  % Conversion Profile
plot(v,Conv,'k-')
grid
xlabel('Reactor Volume - m3')
ylabel('Conversion - %')
title('Conversion Profile')
format long g
disp([ v(end) ysoln(end,:) ])
disp(Conv(end))
SaveVars = [v MflowN2 MflowH2 MflowNH3 P T H Conv];
save AmmoniaFullModel.txt SaveVars -ascii ;



