| 查看: 831 | 回复: 2 | ||
happytan新虫 (小有名气)
|
[求助]
请那位帮我看看这个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]; |
回复此楼» 猜你喜欢
欢迎发来filecode的Mz6后的代码验证其规律
已经有123人回复
-
今天基金会出结果吗？20260819
已经有13人回复
-
时间戳变了，能看出什么问题？
已经有16人回复
-
filecode，4个jtjc了
已经有13人回复
-
为什么余额宝的年化利率越来越低？主要原因有哪些？
已经有5人回复
-
今天放榜没戏了吧
已经有5人回复
-
今天放榜吗？
已经有15人回复
-
filecode=后面第一个是大写字母
已经有12人回复
-
2027广东省杰青
已经有10人回复
-
朋友圈看到的
已经有8人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 帮忙看看这个平面图用matlab怎么画
已经有19人回复
- 师哥师姐帮我看看我这个邀请函符合要求吗~~先谢谢了
已经有12人回复
- 请大家帮我看看这个核磁谱图吧，实在是没有办法了
已经有30人回复
- 【求助】请教一个matlab程序，谢谢帮忙！！
已经有15人回复
- 【请教】大家帮我看看这种薄膜样品的截面SEM照片是怎么获得的，谢谢了
已经有18人回复
- 【求助】急急！！！哪位懂xps，帮我看一下我的数据吧。谢谢
已经有5人回复
- 【求助】大家帮我看看这个MATLAB程序？
已经有6人回复
2楼2011-09-28 16:59:12
信彼南山
木虫 (著名写手)
- 应助: 33 (小学生)
- 金币: 4142.9
- 散金: 1221
- 红花: 16
- 帖子: 1178
- 在线: 233.5小时
- 虫号: 1133529
- 注册: 2010-10-27
- 专业: 导航、制导与传感技术
3楼2011-09-28 19:44:47
