nonisothermtr.m

Matlab在化学工程中的应用

function NonIsothermTR
% 模拟计算非等温固定床管式反应器的轴向温度分布和转化率分布
% 在一列管反应器中进行邻二甲苯(A)氧化制邻苯二酸酐(B)
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/05/27 $

clear all
clc
global  Ct Ac rhoB Cp H1 H3 U dt TJ Cp H1 H3 U TJ rhog us E1 E2 E3 R yO2
L = 1;                  % 反应管长, m
G = 4684;               % 表观质量流速, kg/m2 hr
rhoB = 1300;            % 催化剂堆积密度, kg/m3
Mm = 29.48;      	    % 气体的平均分子量, kg/kmol
yA0 = 0.00924;  	    % 入口处邻二甲苯的摩尔分率
yO2 = 0.208;    		% 氧的摩尔分率(恒为常数)
Cp = 1.047;      	    % 比热, kJ/kmol K
U = 345.686; 		    % 传热系数, kJ/m2 hr K
P = 101.325;            % 假设总压为定值 = 1 atm = 101.325 kJ
dt = 0.0254;    	    % 管径, m  
TJ = 580;      		    % 冷却温度, K
T0 = 700;     	        % 物料进口温度(初始温度), K 
H1 = -1.285e+6;   	    % 反应A→B的反应热, kJ/kmol
H3 = -4.564e+6;         % 反应A→C的反应热, kJ/kmol

% 活化能, kJ/kmol
E1 = 1.1304e5;           
E2 = 1.315e5;		    
E3 = 1.197e5;		   

R = 8.314;  		    % 理想气体常数, kJ/kmol K

Ac = pi*(dt/2)^2;       % 反应管的横截面积, m2
Ft = G*Ac/Mm;		    % 总摩尔流率, moles/hr

us = 3600;              % 线速度，m/hr  
rhog = G/us;
Ct  = Ft/(Ac*us);
FA0 = yA0*Ft;           % A的进料摩尔流率, kmol/hr
CA0 = FA0/(Ac*us);
CB0 = 0;                % FB0 = 0
CC0 = 0;                % FC0 = 0
[z, y] = ode45(@Equations, [0 L], [CA0 CB0 CC0 T0])
CA = y(:, 1);
CB = y(:, 2);
CC = y(:, 3);
xA = (CA0-CA)./CA0;		            % A的转化率
xB = CB(2:end)./(CA0-CA(2:end));	% 生成的B/反应的A 
xB = [0; xB]
xC = CC(2:end)./(CA0-CA(2:end));	% 生成的C/反应的A
xC = [0; xC]

% 图形输出
plot(z, y(:, 4))        % 温度分布
xlabel('z')
ylabel('T (K)')
figure
plot(z, xA, 'r-')       % 转化率分布
xlabel('z')
ylabel('x_A')
figure
plot(z, CA, 'r-', z, CB, 'k--', z, CC, 'b:')    % 浓度分布
xlabel('z')
ylabel('C_A, C_B, C_C')
legend('C_A', 'C_B', 'C_C')

% ------------------------------------------------------------------
function dydz = Equations(z, y)                         % 模型方程组
global  yO2 Ct Ac rhoB Cp H1 H3 U dt TJ Cp H1 H3 U TJ rhog us
CA = y(1);
CB = y(2);
CC = y(3);
T = y(4);

% 摩尔分率
yA = CA/Ct;
yB = CB/Ct;
yC = CC/Ct;

% 反应速度
[rA, rB, rC, k1, k2, k3] = Rates(yA, yB, yC, T);

% 物料平衡
dCAdz = rhoB*rA/us;
dCBdz = rhoB*rB/us;
dCCdz = rhoB*rC/us;

% 热量衡算
dTdz = ( rhoB*(-H1*k1 -H3*k3)*yA*yO2-4*U*(T-TJ)/dt )/(us*rhog*Cp);
dydz = [dCAdz; dCBdz; dCCdz; dTdz];

% ------------------------------------------------------------------
function [rA, rB, rC, k1, k2, k3] = Rates(yA, yB, yC, T)    % 反应动力学
global E1 E2 E3 R yO2

% 速度常数, kmol/kg catalyst hr
k1 = exp(-E1/(R*T) + 19.837);
k2 = exp(-E2/(R*T) + 20.86);
k3 = exp(-E3/(R*T) + 18.97);

% 反应速度, kmol/kg catalyst hr
rA = -(k1+k3)*yA*yO2;               % A的总反应速度
rB = k1*yA*yO2 - k2*yB*yO2;         % B的净生成速率
rC = k2*yB*yO2 + k3*yA*yO2;         % C的总生成速率