Figure 6.38

% 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
Figure 6.38:

Coolant temperature at reactor outlet versus temperature at reactor inlet; three steady-state solutions.

Code for Figure 6.38

main.m

ivpjbr.m
