| 查看: 540 | 回复: 0 | ||
happytan新虫 (小有名气)
|
[求助]
请那位帮我看看这个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]; |
回复此楼» 猜你喜欢
能否退出参与的面上项目解除限项
已经有21人回复
-
没有任何消息-是不是就凉了
已经有3人回复
-
建议基金发布提前给出明确的时间点
已经有16人回复
-
如果此刻你正在为国基感到焦虑，不妨来听听这首《基金之外》
已经有4人回复
-
今日不放榜？网传国自然预计 8 月 27 日可查结果
已经有21人回复
-
估计是周四
已经有4人回复
-
我面上完蛋了
已经有9人回复
-
2026国自然函评费到账
已经有18人回复
-
范进中举一文的中心思想
已经有7人回复
-
朋友圈看到的
已经有9人回复
找到一些相关的精华帖子，希望有用哦~
帮忙看看这个平面图用matlab怎么画
已经有19人回复
- 师哥师姐帮我看看我这个邀请函符合要求吗~~先谢谢了
已经有12人回复
- 请大家帮我看看这个核磁谱图吧，实在是没有办法了
已经有30人回复
- 【求助】请教一个matlab程序，谢谢帮忙！！
已经有15人回复
- 【请教】大家帮我看看这种薄膜样品的截面SEM照片是怎么获得的，谢谢了
已经有18人回复
- 【求助】急急！！！哪位懂xps，帮我看一下我的数据吧。谢谢
已经有5人回复
- 【求助】大家帮我看看这个MATLAB程序？
已经有6人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
