Figure 6.40 (page 329):

Ammonia mole fraction versus reactor length.

Code for Figure 6.40

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 K0 Tstd delH0 P R Rg Nn0 Nh0 Nnh0 A k0 E delH0 K0 Tstd deltaH ...
    UA Ta0 len  xnh0 xh0 xn0 Q0 tau rxrate
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,250)'; % 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.
%%
% = 2.4e12/2*(3/4)^(1.5); % reverse rate constant
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
Nnh0 = xnh0*Nt0;
Nh0 = xh0*Nt0;
Nn0 = xn0*Nt0;



%
%% Solve for the multiple steady states using either T0guess=360 or 600
%% or 800 K
%

T0guess=[300;600;900];
% T0guess=[600];

table1 = [len];
for i=1:length(T0guess)
  opts = optimset ('MaxFunEvals', 2000, ...
		   'MaxIter', 500);
  [T0,fval,info] = fsolve('feed_condition',T0guess(i),opts);
%  T0, info

%% Solve for the exact profiles for no. of moles of Ammonia, x(:,1), and
%% temperatures x(:,2) and x(:,3)

  x0=[Nnh0;T0;T0];
  opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
  [tsolver, x] = ode15s(@ivpjbr,len,x0,opts);
  NT = Nn0 + Nh0 + 2*Nnh0 - x(:,1);
  xnh = x(:,1) ./ NT;
  table1 = [table1 x xnh];
end




subplot(2,1,1);
plot (table1(:,1),table1(:,[3,4,7,8,11,12]));
title ('Figure 6.40')

subplot(2,1,2);
plot (table1(:,1),table1(:,[5,9,13]));
title ('Figure 6.40')

ivpjbr.m

function xdot=ivpjbr(t,x)

  global P R Rg Nn0 Nh0 Nnh0 A k0 E delH0 Tstd K0 deltaH UA Q0 tau rxrate
  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;
  xdot(2) = qadd - rxrate*deltaH;
  xdot(3) = qadd;

feed_condition.m

function error = feed_condition (T0)
  global Ta0 len Nnh0 Nn0 Nh0 P R xnh0 xh0 xn0 Q0
  x0 = [Nnh0;T0;T0];
  opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
  [tsolver, x]  = ode15s(@ivpjbr,len,x0,opts);
  Taend = x((length(x)),3);
  error = Ta0 - Taend;