24小时热门版块排行榜

返回列表

happytan

新虫 (小有名气)

应助: 0 (幼儿园)
金币: 76.8
散金: 50
帖子: 214
在线: 33.2小时
虫号: 71431
注册: 2005-05-31
性别: GG
专业: 有机合成

[求助] 请那位帮我看看这个MATLAB程序那里有问题！先谢谢

请那位帮我看看这个MATLAB程序那里有问题！先谢谢

function RH_Cat_Kinetics %用数值积分法进行反应速率分析得到速率常数，并将结果与实验值对比分析
clear all;clc
global C kk0

t=[0:300:2400];    % t:时间，s
C=[0.00 0.11 0.19 0.25 0.31 0.44 0.54 0.62 0.65;
0.00 0.10 0.18 0.24 0.28 0.39 0.48 0.52 0.53;
0.10 3.20 1.70 1.30 2.50 3.90 5.20 8.20 9.7];
C=C';%温度150℃的动力学实验数据C(1) &C(2): kmol/m^3, C(3):%
%非线性拟合
kk0=[5.1095e-004; 2.8939;9.5453e+003];
                  %利用suresh模型算出同样温度条件下的速率常数作为其初始值
tspan=[0:300:2400]; %时间阶梯
C0=[0;0;0.1];    % 0.1:氧气在反应体系气相中的初始分率，考虑了环己烷的饱和蒸
                  %汽压，在实验条件下为50%左右，除去该蒸汽压以及惰性气体所
                  %占的比例后得出的上述数值
kk=lsqnonlin(@myfunc,kk0,0,inf,[ ],tspan,C0,C); %非线性最小二乘法拟合。调用函数myfunc
ci=nlparci(beta,resid,jacobian);    %回归系数，残差，雅阁比矩阵
%拟合效果图(实验值与拟合值的比较)
M=[ 1 0 0 0; %解微分代数方程组时的系数矩阵
0 1 0 0;
0 0 0 0;
0 0 0 0];
options=odeset;  %产生/改变参数结构
options.Mass=M;
options.TolRol=le-7;%定义精度
tspan=[5:40];
C0=[0 0 0.1 ];%初值C
P=1.2;%MPa

G=1.5;%L/min
kla=0.14;%1/s;
apxl=input('Pls.input gas holdup you got according reaction conditions using hydro_mass_ c:','s');
[t_sim,C_sim]=ode45(@lj_ssm,tspan,CO,options,kk,P,G,kla,apxl);
plot(t,C_exp(:,1),'-ks',t_sim,C_sim(:,1),':ko',t,C_exp(:,2),'-kd',t_sim,C_sim(:,2),':k+');
legend('Exp.CP','Model.CP','Exp.CI','Mode1.CI'),xlabel('时间t/s'),ylabel('浓度C/(kmol/m^3)');

%定义非线性优化的目标函数
function f=myfunc(kk,tspan,C0,C)
P=12;
G=1.5;
kla=0.14;
apxl=0.2;
[t,CC]=ode45(@lj_ssm,tspan,C0,[],kk,P,G,kla,apxl);
f=CC-C;

%定义待求的动力学方程
function dCdt=lj_ssm(t,C,kk,P,G,kla,apxl)
G=G/(60*22.4);%改变气体流量的单位，便于后面的等式成立，,mol/s
yin=0.21;
GI=G*(1-yin);%按照进气中氧气含量21%计算惰性气体的体积流量
Vm=22.4;%气体的摩尔体积，,L
Vl=0.3;          % VL:液体体积，按300mL计算
VG1=Vl/(1-apxl); %反应体系的总体积，L
P=P*10;          % P原来的单位是MPa，为了能够利用后面的公式，将单位转变为bar
H=1.09e-2;       % H:氧气的Henry系数
kk=kk0.*exp(-Ea/(R*T)); % T:K k0:指前因子Ea:活化能
dCdt=zeros(1,4);%先预分配空间给浓度C，这样可以加快调用速度
dCPdt=kk(1)*kk(3)*C(1)*kk(4)/(kk(1)+kk(2)*C(1)+kk(3)*kk(4));
dC1dt=kk(3)*C(1)*kk(4)*(kk(1)-kk(2)*C(3))/(kk(1)+kk(2)*C(1)+kk(3)*kk(4));
gas_mass=Vm*(1-apxl)/(Vl*apxl)*(GI*(yin/(1-yin)-C(4)/(1-C(4)))-kla*(H*P*C(4)-kk(4))*VG1);
liquid_mass=kla*(H*P*C(4)-C(2))-(kk(3)*C(1)*C(2)*(kk(1)+kk(2)*C(1))/(kk(1)+kk(2)*C(1)+kk(3)*C(2))); % gas mass & liquidse mass分别达到稳态时氧气在气液相中的物料平衡方程左端
dCdt=[dCPdt;dC1dt;gas_mass;liquid_mass];

回复此楼