Figure 4.20 (page 152):

Benzene conversion versus reactor volume.

Code for Figure 4.20

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.

global k1 K1 k2 K2 R T P NBf

% Example in material balance chapter.
% Data are: 
%           NBf = 60e3 mol/hr
%           R  = 0.08205 lit atm/mol K
%           T  = 1033 K
%           P  = 1 atm
%           k1 = 7e5 lit/mol hr
%           k2 = 4e5 lit/mol hr
%           K1 = 0.31
%           K2 = 0.48
%           Reactor volume: 0 to 1600 lit

NBf = 60e3;
R  = 0.08205;
T  = 1033;
P  = 1;
k1 = 7e5;
k2 = 4e5;
K1 = 0.31;
K2 = 0.48;

x0=[0;0];

vout = linspace(0,1600,50)';


opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, x] = ode15s(@benzene_pyrol_rhs, vout, x0, opts);

conv = ( 2*x(:,1) + x(:,2) ) /NBf;
yB   = ( NBf - 2*x(:,1) - x(:,2) ) /NBf;
yD   = ( x(:,1) - x(:,2) ) /NBf;
yT   = ( x(:,2) ) /NBf;
yH   = ( x(:,1) + x(:,2) ) /NBf;

data = [vout conv yB yD yT yH];

subplot(2,1,1);
plot (data(:,1),data(:,2));
title ('Figure 4.20')

subplot(2,1,2);
plot (data(:,1),data(:,3:6));
title ('Figure 4.20')

benzene_pyrol_rhs.m

function xdot = benzene_pyrol_rhs(volume,x)
     global k1 K1 k2 K2 R T P NBf

     % Example in material balance chapter.
     % Data are: 
     %           NBf = 60e3 mol/hr
     %           R  = 0.08205 lit atm/mol K
     %           T  = 1033 K
     %           P  = 1 atm
     %           k1 = 7e5 lit/mol hr
     %           k2 = 4e5 lit/mol hr
     %           K1 = 0.31
     %           K2 = 0.48
     %

     ext1 = x(1);
     ext2 = x(2);

     NB = NBf - 2*ext1 - ext2;
     ND = ext1 - ext2;
     NH = ext1 + ext2;
     NT = ext2;

     Q = NBf*R*T/P;

     cB = NB/Q;
     cD = ND/Q;
     cT = NT/Q;
     cH = NH/Q;

     xdot = [k1*(cB*cB - cD*cH/K1); k2*(cB*cD - cT*cH/K2)];