isothermtr.m

《实用化工计算机模拟：MATLAB在化学工程中的应用 》这本书光盘里的程序~

function IsothermTR
% 等温管式反应器（苯脱氢反应）
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center,
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/07/21 $

clear all
clc
global keq1 keq2 k1 k2 P

T = 573.9;            % 反应温度, K
keq1 = 0.312;       % 平衡常数
keq2 = 0.480;       % 平衡常数
P = 1;
FA0 = 128.2;        % 纯苯的进料流量, lbmol/h
CA0 = 1/22.4;       % 标准状态下的进料浓度，mol/L
CA0 = CA0/453.6/(1e-3 * 3.2808^3);     % lbmol/ft^3  （1 mol = 1/453.6 
lbmol, 1 m = 3.2808 ft）
v0 = FA0/CA0;               % 标准状态下的进料体积流量，ft^3/h

% 初始转化率
x0 = [0 0];

% 反应速度常数, 1/atm2 hr
k1 = 1.496e7 *exp(-8.44e3/T);
k2 = 8.67e6*exp(-8.44e3/T);

%  In the following, t represents (V / FA0), its unit: ft^3.hr/lbmol
t = [0  0.005  0.010  0.020  0.040  0.060  0.080  0.100  0.120  0.140  0.180 
  0.220  0.260  0.300  0.400];
[t, x] = ode45(@Euqations, t, x0);          % 解ODE方程组
tau = CA0*t;                        % 空时，hr
SV = 1./tau;                         % 空速，1/hr
xBt = x(:,1) + x(:,2);                % xBt: 
苯转化为二苯基和三苯基的总转化率
t42p = spline(xBt, t, 0.42);   % 总转化率为42%时的(VR / FA0)
VR = t42p * FA0;                 % 总转化率为42%时所需的反应器体积, ft^3
PB = 1 - x(:,1) - x(:,2);          % 苯的分压, atm
PD = x(:,1)/2 - x(:,2);            % 二苯基的分压, atm
PT = x(:,2);                          % 三苯基的分压, atm
tmax = t(find(PD==max(PD)));             % 
使二苯基的分压(浓度)最大所对应的(V / FA0)
tau_max = tau(find(PD==max(PD)));   % 使二苯基的分压(浓度)最大所对应的空时
SV_max = SV(find(PD==max(PD)));    % 使二苯基的分压(浓度)最大所对应的空速

% 输出结果
disp('  Results:')
fprintf('\n    所需的反应器体积为：%.1f %s\n\n', VR, 'ft^3 
(苯的总转化率为42%时)')
fprintf('    在V/FA0为%.3f %s\n', tmax, 
'ft^3·hr/lbmol下，可使二苯基的浓度达最大值')
fprintf('    即在空速SV为%4.1f %s\n\n', SV_max, 
'1/hr下，可使二苯基的浓度达最大值')
disp('       t         x1       x2       xBt        PB        PD        PT')
disp([t  x(:,1:2) xBt PB PD PT])
plot(t, xBt)
xlabel('V/F_A_0,  ft^3·hr/lbmol')
ylabel('x_B_t')
figure
plot(t, PB, 'b--', t, PD, 'm-.', t, PT, 'r-')
xlabel('V/F_A_0,   ft^3·hr/lbmol')
ylabel('P,  atm')
legend('P_B', 'P_D', 'P_T')


% ------------------------------------------------------------------
function dxdt = Euqations(t, x)        % here t = V / FA0
global keq1 keq2 k1 k2 P
% 分压, atm
PB = P*(1 - x(1) - x(2));
PD = P*(x(1)/2 - x(2));
PH = P*(x(1)/2 + x(2));
PT = P*x(2);

% 反应速度, lbmol/ft^3.hr
r1 = k1*(PB*PB - PD*PH/keq1);
r2 = k2*(PB*PD - PH*PT/keq2);

% 物料平衡方程
dx1dt = r1;
dx2dt = r2;
dxdt = [dx1dt; dx2dt];