nonisothermcstr.m

这是一些关于MATLAB的小程序

function NonIsothermCSTR
% 非等温CSTR反应器(在t=0 hr，进料浓度从0.5阶跃变为0.55的情况)
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center,
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/07/20 $

clear all
clc

global F0 F CA0 Kc T0 k0 dH rho rhoJ Cp VJ CJ U Ah FJ TJ0

% initial conditions
CA = 0.245;
T = 600;
TJ = 594.6;
V = 48;
VC = V*CA;
VT = V*T;

% parameter values
TJ0 = 530;
F0 = 40;
T0 = 530;
CA0 = 0.5;
Kc = 4;
delta = 0.01;

% disturbance
CA0 = 0.55;

VJ = 3.85;
E = 30000;
U = 150;
Cp = 0.75;
rho = 50;
k0 = 7.08e10;
R = 1.99;
Ah = 250;
dH = -30000;
CJ = 1.0;
rhoJ = 62.3;
Tset = 600;
[t,y] = ode23s(@Equations,[0:0.2:4],[V VC VT TJ]);  % solve ODEs
V = y(:,1);
CA = y(:,2)./V;
T = y(:,3)./V;
TJ = y(:,4);
FJ = 49.9-Kc*(600-T);
F = 40-10*(48-V);
Q = U*Ah*(T-TJ)

disp('  结果:')
disp('       t         CA      T          V         F        TJ        FJ')
disp( [t CA T V F TJ FJ] )
plot(t,CA)
xlabel('t，h')
ylabel('C_A，lbmol/ft^3')
axis([0 4 0.22 0.26])
figure
plot(t,T)
xlabel('t，h')
ylabel('T，^oR')
axis([0 4 600 610])
figure
plot(t,V)
xlabel('t，h')
ylabel('V，ft^3')
axis([0 4 40 50])
figure
plot(t,F)
xlabel(' t, h')
ylabel('F，ft^3/h')
axis([0 4 20 45])
figure
plot(t,TJ)
xlabel('t，h')
ylabel('T_J，^oR')
axis([0 4 590 600])
figure
plot(t,FJ)
xlabel('t，h')
ylabel('F_J，ft^3/h')
axis([0 4 20 80])


% ------------------------------------------------------------------
function dydt = Equations(t,y)
global F0 F CA0 Kc T0 k0 dH rho rhoJ Cp VJ CJ U Ah FJ TJ0

V  = y(1);
VC = y(2);
VT = y(3);
TJ = y(4);
CA = VC/V;
T = VT/V;

% feedback controllers
FJ = 49.9-Kc*(600-T);
F = 40-10*(48-V);

k = k0*exp(-30000/(1.99*T));   % reaction rate

Q = U*Ah*(T-TJ);

% simultaneous differential equations
dVdt = F0 - F;
dVCdt = F0*CA0 - F*CA - V*k*CA;
dVTdt = F0*T0-F*T-(dH*V*k*CA+Q)/(rho*Cp);   % dH:反应热
dTJdt = FJ*(TJ0-TJ)/VJ+Q/(rhoJ*VJ*CJ);
dydt = [dVdt;dVCdt;dVTdt;dTJdt];