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的總生成速率