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