% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING. If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.
% This program simulates the tubular reactor model for ammonia
% synthesis
%
% Revised 8/14/2018
p.P = 300; % atm - pressure of the reactor
p.R = 82.05e-6; %(m^3 atm)/(mol K) - Gas constant
% p.v0 = 0.16; % m/sec
% Repairing bug.
% Change velocity to get profile C in Figures 6.40, 6.41
% jbr, 3/5/07
p.v0 = 0.05713; % m/sec
p.A = 1; % m^2 - Cross sectional area available for reaction
%p.UA = 1/3; % 1/m
p.UA = 1/2; % 1/m
p.delG0 = -4250; % cal/mol ammonia
p.delH0 = -12000; % cal/mol ammonia
p.Tstd = 298; % K
p.Rg = 1.987; % cal/mol K
p.K0 = exp(-p.delG0/(p.Rg*p.Tstd)); % equilibrium constant at Tstd=298
%p.deltaH = 495; % K/mol
p.deltaH = - 2.342 ; % s K/mol
p.tau = 1500; % K
p.Ta0 = 323; % K - entering coolant (feed) temperature of the reacting gas
p.len = linspace(0,12,100)'; % m - length of reactor
%
% UA corresponds to Ua/C where U is the heat transfer coefficient, a is
% the total heat exchanging surface per unit length of the converter
% and C is the feed rate expressed as heat capacity per unit time.
%
% deltaH is the temperature rise of the reacting gas corresponding to
% the adiabatic formation of 1 mole of Ammonia.
%
% The following set of constants correspond to the rate constants(k1m
% and k2m) and activations energies/gas constants(E1m,E2m) of the
% forward and reverse reactions. Tm is the
% mean temperature which has been used to reparameterise.
%
p.k0 = 7.794e11; % reverse rate constant
p.E = 20000; % K reverse rate activation energy
%
% Initial flows
%
initial.Q0 = p.v0*p.A; % m^3/s
initial.Nt0 = initial.Q0*p.P/(p.R*p.Ta0);
xn0 = 0.985*0.25; % N_2 feed mole fraction
xh0 = 0.985*0.75; % H_2 feed mole fraction
xnh0 = 0.015; % Ammonia feed mole fraction
%
% Plot a graph showing multiple intersections of coolant feed
% temeperature Ta0 and solution to ODE on coolant side T(:,3)
%
nTs = 100;
Tvec= linspace(273,1000,nTs)';
Taend = zeros(nTs,1);
Trend = zeros(nTs,1);
xnhvec = zeros(nTs,1);
for i = 1:nTs
T = Tvec(i);
inital.Nnh0 = xnh0*initial.Nt0;
inital.Nh0 = xh0*initial.Nt0;
inital.Nn0 = xn0*initial.Nt0;
x0 = [inital.Nnh0; T; T];
opts = odeset('AbsTol', sqrt(eps), 'RelTol', sqrt(eps));
[tsolver, x] = ode15s(@(t, x) ivpjbr(t, x, p, inital), p.len, x0, opts);
NT = inital.Nn0 + inital.Nh0 + 2*inital.Nnh0 - x(:,1);
Taend(i) = x((length(x)),3);
Trend(i) = x((length(x)),2);
xnhvec(i) = x((length(x)),1)./ NT(length(x));
end
table = [Tvec p.Ta0*ones(nTs,1) Taend Trend xnhvec];
save -ascii ammonia_shooting.dat table;
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
plot (table(:,1),table(:,2:3));
% TITLE
end % PLOTTING
