cstr.m

本代码为黄华江编著《实用化工计算机模拟—MATLAB在化学工程中的应用》的配套车程序

function CSTR
% CSTR反应器的热稳定性分析（Thermal stability analysis of a CSTR）--相平面图
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/08/03 $

global E R k0
global F CA0 V T0 UA TJ rho Cp Hr HG HL tau
STOPTIME = 5000;
F = 1.0e-8;         % 液体体积流量, m3/s 
CA0 = 5.0;          % 进料浓度, kmol/m3
T0 = 300;           % 进料温度, K   
TJ = 305;           % 夹套温度, K
V = 2.0e-6;         % 反应器体积,m3 
rho = 1000.0;       % 液体密度, kg/m3
Cp = 4.187;         % 液体热容, kJ/kg K 
UA = 5.68e-6;       % 传热系数与传热面积的乘积, kJ/K s
Hr = -4.19e4;       % 反应热, kJ/kmol  
k0 = 8.03E+12;      % 指前因子, 1/s	
E = 9.42e+4;        % 活化能, kJ/kmol 
R = 8.317;          % 气体常数, kJ/kmol K 
tau = V/F;

% 绘制Qg (Qr)~T的关系图以便观察稳态点
% -----------------------------------
T = (290:380);
k = RateConstant(T)
CA = CA0./(1+tau*k);                % STEADY STATE MASS BALANCE
[Qg,Qr] = Heat(k,CA,T);
figure
plot(T,Qg,'b-',T,Qr,'k--');
xlabel('T (K)');
ylabel('Q_g, Q_r (kJ/s)');
legend('Q_g','Q_r')

% 动态计算
% --------
tspan = [0  STOPTIME];
CAI = linspace(0,5,2)       % Initial temp, K	
TI = linspace(300,340,10)   % Initial conc., kmol/m3
% 绘制CSTR反应器的相平面图
figure
hold on
for i=1:length(CAI)
    for j=1:length(TI)
        [t,y] = ode45(@DynamicModel,tspan,[CAI(i),TI(j)]);
        CA = y(:,1);
        T = y(:,2);
        xA = 1 - CA/CA0;
        plot(T,xA);
    end
end
xlabel('T (K)');
ylabel('x_A');
hold off


% ------------------------------------------------------------------
function dydt = DynamicModel(t,y)           % Dynamic Model
global F CA0 V T0 UA TJ rho Cp Hr tau
CA = y(1);
T = y(2);
k = RateConstant(T);
% Dynamic mass balance, first order kinetics
dCAdt = (CA0-CA)/tau - k*CA;         

% Dynamic energy balance 
[Qg,Qr] = Heat(k,CA,T);
dTdt = (Qg - Qr)/(V*rho*Cp);                % 热量平衡, K/s
dydt = [dCAdt; dTdt];


% ------------------------------------------------------------------
function k = RateConstant(T)
global k0 E R
% 以下两式由式k=k0exp(-E/R/T)变换而成，通过此变换以避免溢出。
temp = k0*exp(-E/(R*273));                  % Rate constant at 273 K, 1/s
k = temp*exp(-(E/R)*(1./T-1/273));          % This form avoids math overflow

% ------------------------------------------------------------------
function [Qg,Qr] = Heat(k,CA,T)
global F T0 V UA TJ rho Cp Hr
Qr = F*rho*Cp*(T-T0) + UA*(T-TJ);           % 移走的热量, kJ/s
Qg = k.*CA*V*(-Hr);                         % 产生的热量, kJ/s