连续搅拌反应釜(CSTR)代码和实验仿真2019-06-20

如需转载，请标记出处。

在之前的一篇文章连续搅拌反应釜(CSTR)2019-06-11中我介绍了CSTR的模型机理和数学微分方程，以及参数设置可参考的文献。接下来将给出模型求解的代码，给出输入输出数据。

1. 数学微分方程代码：

function [CA,T]=CSTR_DIS(CA_INI,T_INI,DeltaT,qc)

% Parameters:
% Process Flowrate (m^3/sec)
q = 100;
% Feed Concentration (mol/m^3)
CA0 = 1;
% Feed Temperature (K)
T0 = 350;
% Inlet Coolant Temperature (K)
Tc0 = 350;
% Volume of CSTR (m^3)
V = 100;
% Heat Transfer Term (cal/(min K))
hA = 7e5;
% Reaction Rate Constant (min^{-1})
k0 = 7.2e10;
% E - Activation energy in the Arrhenius Equation (J/mol)
% R - Universal Gas Constant = 8.31451 J/mol-K
EoverR = 1e4;
% Heat of Reaction (cal/mol)
DeltaH = 2e5;
% Liquid Density (g/L)
rou = 1e3;
rouc = 1e3;
% Specific Heats (cal/(g K))
Cp = 1;
Cpc = 1;
% INITIALIZATION OF THE PARAMETERS
TimeLen=length(qc);
CA=zeros(TimeLen+1,1);
T=zeros(TimeLen+1,1);
CA(1)=CA_INI;
T(1)=T_INI;

for i=1:TimeLen
    CA(i+1)=CA(i)+DeltaT*(q/V*(CA0-CA(i))-k0*exp(-EoverR/T(i))*CA(i));
    T(i+1)=T(i)+DeltaT*(q/V*(T0-T(i))+DeltaH/(rou*Cp)*k0*exp(-EoverR/T(i))*CA(i)+rouc*Cpc*qc(i)/V/rou/Cp*(1-exp(-hA/(rou*Cpc*qc(i))))*(Tc0-T(i)));
end
CA=CA(2:TimeLen+1);
T=T(2:TimeLen+1);
end

2. 通过给定初始的值CA,T，以及输入qc,带入函数即可求解得到数据。

clc;
clear;
N = 1000;% 样本量
qc1 = unifrnd(98,108,[N,1]); % 输入qc,均匀分布
qc2 = linspace(98,108,N);% 输入qc
Initi_CA1=0.0795; % mol/L C_A
Initi_T1=443.4566; %437.44; % K T
DeltaT=0.1;
[CA1,T1]=CSTR_DIS(Initi_CA1,Initi_T1,DeltaT,qc1);
[CA2,T2]=CSTR_DIS(Initi_CA1,Initi_T1,DeltaT,qc2);
subplot(4,1,1);
plot(CA1,'b-.');
xlabel('采样时刻');
ylabel('CA1');
subplot(4,1,2);
plot(CA2,'g-.');
xlabel('采样时刻');
ylabel('CA2');
subplot(4,1,3);
scatter(qc1,CA1,'m.');
xlabel('qc1');
ylabel('CA1');
subplot(4,1,4);
scatter(qc2,CA2,'c.');
xlabel('qc2');
ylabel('CA2');

结果如图所示： CSTR_该_大白.jpg