| 查看: 3000 | 回复: 13 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
apolloking金虫 (小有名气)
|
[交流]
【求助】使用Matlab预估动力学方程问题已有4人参与
|
||
|
大家好，问题如下： 使用固定床积分反应器来测量甲烷催化燃烧的本征动力学实验数据。 整个计算过程见图片附件 本征动力学方程为ra=k0*exp（-E/(R*T)）*Ca ra:反应速率 mol/(kgcat*s) k0:指前因子 m3/(kgcat*s) E:活化能 kJ/mol Ca：甲烷摩尔浓度 mol/m3 待定参数为k0和E 文献使用了Matlab的半隐式龙格库塔法和最小二乘法去约束优化 最后得出的k0为3.77*10E8 E为25kJ/mol 在积分反应器上有; -ra=dxa/d(w/Fa0)=-k0*exp（-E/(R*T)）*Ca xa为转化率 w为催化剂的质量 Fa0为流速 又因为Ca=C0*(1-xa) 所以上式转化为 dxa/d(w/Fa0)=-k0*exp（-E/(R*T)）*C0（1-xa） 因为实验所用的甲烷体积浓度为1%，上式可转化为 1%*22.4L/mol dxa/d(w/Fa0)=-k0*exp（-E/(R*T)）*（1-xa）/2.24 （我的理解是w/Fa0是空速的倒数,但文献中给的空速值单位是ml/(g*h)，我不知道需不需要在进行单位转化，还是直接取倒数就行了） 文献中说，“采用Matlab程序对本征动力学实验数据进行搜索，选优”求的k0和E 按照这个思路，是不是针对每一个温度，认为温度值已知，空速和转化率为变量求的一组k0和E。然后将所有温度得到的k0和E值取平均值？ 我自己写了一个550度的计算程序，但是得不到正确的结果，希望大家能帮忙看一下。 function jifen %积分法动力学参数估算 close all, clear, clc, %实验数据 t=[1/5000 1/10000 1/15000 1/20000 1/25000]; ca=[0.989578 0.972126 0.956323 0.897261 0.870735]; %最小二乘法估计动力学参数 x0=[0 1]; x=lsqnonlin(@myfun,x0,[],[],optimset('Display','iter'),ca,t); k0=x(1); E=x(2); %给出估算结果 fprintf('动力学参数为：\n'), fprintf('k=%f\ta=%f\n',k0,E), %绘制反应速率拟合曲线 [T,CA]=ode45(@reaction_rate,[0,1/5000],ca(1),[],k0,a); plot(t,ca,'*k'), hold on, plot(T,CA,'-k'), xlabel('t/min') ylabel('C_A/mol/L') legend('实验数据','模拟曲线') %实验导数与计算导数之差 function F=myfun(x,ca,tspan) k0=x(1); E=x(2); [t,cac]=ode45(@reaction_rate,tspan,ca(1),[],k0,E); F=cac-ca'; %反应速率方程 function r=reaction_rate(t,ca,k0,E) R=8.314; T=823.15; r=k0*exp(-E/R/T)*(1-ca)/2.24; [ Last edited by apolloking on 2009-7-18 at 18:11 ] |
回复此楼» 猜你喜欢
博士读完未来一定会好吗
已经有23人回复
-
导师想让我从独立一作变成了共一第一
已经有7人回复
-
到新单位后，换了新的研究方向，没有团队，持续积累2区以上论文，能申请到面上吗
已经有11人回复
-
读博
已经有4人回复
-
JMPT 期刊投稿流程
已经有4人回复
-
心脉受损
已经有5人回复
-
Springer期刊投稿求助
已经有4人回复
-
小论文投稿
已经有3人回复
-
Bioresource Technology期刊，第一次返修的时候被退回好几次了
已经有9人回复
-
申请2026年博士
已经有6人回复
feng198800
金虫 (正式写手)
- 应助: 15 (小学生)
- 金币: 860.2
- 散金: 234
- 红花: 1
- 帖子: 415
- 在线: 242.1小时
- 虫号: 700353
- 注册: 2009-02-12
- 性别: GG
- 专业: 化学反应工程
13楼2012-03-20 19:53:11
apolloking
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1083.1
- 散金: 35
- 红花: 2
- 帖子: 156
- 在线: 60.5小时
- 虫号: 75702
- 注册: 2005-06-20
- 专业: 能源化工
|
function jifen %积分法动力学参数估算 close all, clear, clc, %实验数据 t=[350 375 400 425 450 475 500 525 550 350 375 400 425 450 475 500 525 550 350 375 400 425 450 475 500 525 550 350 375 400 425 450 475 500 525 550 350 375 400 425 450 475 500 525 550]; w=[22.4*3600/5000 22.4*3600/5000 22.4*3600/5000 22.4*3600/5000 22.4*3600/5000 22.4*3600/5000 22.4*3600/5000 22.4*3600/5000 22.4*3600/5000 22.4*3600/10000 22.4*3600/10000 22.4*3600/10000 22.4*3600/10000 22.4*3600/10000 22.4*3600/10000 22.4*3600/10000 22.4*3600/10000 22.4*3600/10000 22.4*3600/15000 22.4*3600/15000 22.4*3600/15000 22.4*3600/15000 22.4*3600/15000 22.4*3600/15000 22.4*3600/15000 22.4*3600/15000 22.4*3600/15000 22.4*3600/20000 22.4*3600/20000 22.4*3600/20000 22.4*3600/20000 22.4*3600/20000 22.4*3600/20000 22.4*3600/20000 22.4*3600/20000 22.4*3600/20000 22.4*3600/25000 22.4*3600/25000 22.4*3600/25000 22.4*3600/25000 22.4*3600/25000 22.4*3600/25000 22.4*3600/25000 22.4*3600/25000 22.4*3600/25000]; xa=[0.224067 0.308978 0.418433 0.583276 0.712133 0.834063 0.927528 0.971279 0.989578 0.167333 0.26823 0.359565 0.492762 0.605102 0.706325 0.870419 0.936043 0.972126 0.161258 0.229717 0.305041 0.445914 0.562024 0.702027 0.847299 0.901723 0.956323 0.703751 0.413094 0.227758 0.297666 0.456384 0.592814 0.743624 0.825191 0.897261 0.089297 0.149173 0.213542 0.304452 0.433972 0.563372 0.703383 0.801953 0.870735]; %最小二乘法估计动力学参数 x0=[3E8 20000]; x=lsqnonlin(@myfun,x0,[],[],optimset('Display','iter'),t,xa,w); k0=x(1); e=x(2); %给出估算结果 fprintf('动力学参数为：\n'), fprintf('k=%f\ta=%f\n',k0,e), %绘制反应速率拟合曲线 [W,XA]=ode45(@reaction_rate,[0,22.4*3600/5000],xa(1),[],k0,e); plot(w,xa,'*r'), hold on, plot(W,XA,'-r'), xlabel('t/min') ylabel('C_A/mol/L') legend('实验数据','模拟曲线') %实验导数与计算导数之差 function F=myfun(x,xa,tspan) k0=x(1); e=x(2); [w,xac]=ode45(@reaction_rate,tspan,xa(1),[],k0,e); F=xac-xa'; %反应速率方程 function r=reaction_rate(t,w,xa,k0,e) r=-k0*exp(-e/8.314/(t+273.15))(1-xa)/2.24; 错误提示 ??? Error: File: D:\My Documents\Untitled3.m Line: 34 Column: 7 ()-indexing must appear last in an index expression. 第34行是 r=-k0*exp(-e/8.314/(t+273.15))(1-xa)/2.24; [ Last edited by apolloking on 2009-7-20 at 11:53 ] |
3楼2009-07-20 11:51:15
change0618
铁杆木虫 (著名写手)
方丈大师
- 应助: 44 (小学生)
- 金币: 17724.5
- 红花: 17
- 帖子: 2413
- 在线: 546.7小时
- 虫号: 496517
- 注册: 2008-01-19
- 专业: 化学反应工程
4楼2009-07-20 13:01:01
apolloking
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1083.1
- 散金: 35
- 红花: 2
- 帖子: 156
- 在线: 60.5小时
- 虫号: 75702
- 注册: 2005-06-20
- 专业: 能源化工
5楼2009-07-20 19:13:06