Figure 6.38 (page 328):

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

Code for Figure 6.38

Text of the GNU GPL.

main.m

%% 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 
%

global P R Rg Nn0 Nh0 Nnh0 A k0 E delH0 Tstd K0 deltaH UA tau

P  = 300;      % atm - pressure of the reactor
R  = 82.05e-6; %(m^3 atm)/(mol K) - Gas constant
% v0 = 0.16;      % m/sec
%% Repairing bug.
%% Change velocity to get profile C in Figures 6.40, 6.41
%% jbr, 3/5/07
v0 = 0.05713;   % m/sec
A  = 1;        % m^2 - Cross sectional area available for reaction
%UA = 1/3;        % 1/m 
UA = 1/2;        % 1/m 
delG0 = -4250;   % cal/mol ammonia
delH0 = -12000; % cal/mol ammonia
Tstd  = 298;   % K
Rg = 1.987;     % cal/mol K
K0 = exp(-delG0/(Rg*Tstd)); % equilibrium constant at Tstd=298
%deltaH = 495;  % K/mol
deltaH = - 2.342 ;  % s K/mol
tau    = 1500; % K
Ta0 = 323;  % K  - entering coolant (feed) temperature of the reacting gas
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.
%%
k0 = 7.794e11;            % reverse rate constant
E = 20000;   % K  reverse rate activation energy
%
%% Initial flows
%
Q0   = v0*A; % m^3/s
Nt0  = Q0*P/(R*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);
  Nnh0 = xnh0*Nt0;
  Nh0 = xh0*Nt0;
  Nn0 = xn0*Nt0;
  x0 = [Nnh0;T;T];
  opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
  [tsolver, x]  = ode15s(@ivpjbr,len,x0,opts);
  NT = Nn0 + Nh0 + 2*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 Ta0*ones(nTs,1) Taend Trend xnhvec];

plot (table(:,1),table(:,2:3));
title ('Figure 6.38')

ivpjbr.m

function xdot=ivpjbr(t,x)

  global P R Rg Nn0 Nh0 Nnh0 A k0 E delH0 Tstd K0 deltaH UA tau
  Nnh = x(1);
  T   = x(2);
  Ta  = x(3);
  RT  = R*T;
  Nn=Nn0 - 1/2*(Nnh - Nnh0);
  Nh=Nh0 - 3/2*(Nnh - Nnh0);
  Nt = Nn + Nh + Nnh;
  pn = (Nn/Nt)*P;
  ph = (Nh/Nt)*P;
  pnh = (Nnh/Nt)*P;
  qadd = (Ta - T)*UA;
  k    = k0*exp(-E/T);
  K = K0*exp(-delH0/Rg*(1/T - 1/Tstd));
  rxrate  = k/RT* ( K^2 * pn*ph^(1.5)/pnh - pnh/ph^(1.5) );
  xdot = zeros (3, 1);
  xdot(1) = A * 2* rxrate;
%%% Van Heerden fudge factor; for comparison purposes only
%    xnh = Nnh/Nt;
%    taufix = - 2*tau*A*(1+xnh)/Nt;
%    xdot(2) = qadd - rxrate*taufix;
  xdot(2) = qadd - rxrate*deltaH;
  xdot(3) = qadd;