应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

Znn3bq.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】使用Matlab预估动力学方程问题
北京石油化工学院2026年研究生招生接收调剂公告
141/2返回列表12下一页
查看: 3222  |  回复: 13
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

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 ]
回复此楼

» 收录本帖的淘帖专辑推荐

动力学拟合

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2009-07-18 17:24:51
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snipher950

木虫 (正式写手)
私下解答。
回复此楼
2楼2009-07-19 14:48:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

apolloking

金虫 (小有名气)
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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
r=-k0*exp(-e/8.314/(t+273.15))*(1-xa)/2.24;
赞一下(1人) 回复此楼
4楼2009-07-20 13:01:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

apolloking

金虫 (小有名气)
改了这个也不对，运行还是有错误
赞一下 回复此楼
5楼2009-07-20 19:13:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hitzhang

木虫 (正式写手)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
这句x=lsqnonlin(@myfun,x0,[],[],optimset('Display','iter'),t,xa,w);，t, xa, w应该是myfun的额外输入参数，而你的myfun函数有x, xa, tspan，不懂作者是什么意思
赞一下(1人) 回复此楼
6楼2009-07-20 22:03:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

apolloking

金虫 (小有名气)
实话说这个程序是我从另外一个程序改过来的，我也不是很清楚这个 tspan具体的含义
赞一下 回复此楼
7楼2009-07-20 22:05:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

apolloking

金虫 (小有名气)
snipher950  网友说先固定一个温度，然后输入空速和转化率，但是之后我就不知道怎么办了。我看相关的文献好像都是一次性就做出来一个k0和E的
赞一下 回复此楼
8楼2009-07-20 22:07:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hitzhang

木虫 (正式写手)
是这样，能把你程序里的函数调用关系以及各个函数的作用描述一下吗？
赞一下 回复此楼
9楼2009-07-20 22:21:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

apolloking

金虫 (小有名气)
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];

%以上是我要输入的参数，一共是45组，每组参数分别为温度t 空速的倒数w 转化率xa
%最小二乘法估计动力学参数
x0=[3E8 20000];
x=lsqnonlin(@myfun,x0,[],[],optimset('Display','iter'),t,xa,w);

%这个是非线性最小二乘法来调用myfun这个函数来判断我计算出来的转化率和实际测得的转化率之差
k0=x(1);
e=x(2);
%给出估算结果
fprintf('动力学参数为：\n'),
fprintf('k=%f\ta=%f\n',k0,e),
%实验导数与计算导数之差
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;


%这个是我的动力学方程
赞一下 回复此楼
10楼2009-07-20 22:30:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 apolloking 的主题更新
141/2返回列表12下一页
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 290求调剂085701 +13 1314捧花 2026-04-02 13/650 2026-04-07 13:14 by guanxin1001
[考研] 085404 293求调剂 +8 勇远库爱314 2026-04-06 9/450 2026-04-07 13:05 by flydream1314
[考研] 284求调剂 +11 梵@@ 2026-04-06 11/550 2026-04-07 09:33 by lovelychen
[考研] 327求调剂 +3 Xxjc1107. 2026-04-06 3/150 2026-04-06 22:34 by qlm5820
[考研] 一志愿C9的化学工程（085602） 340分，感觉校内调剂无望，求调剂 +12 万事宜臻 2026-04-04 12/600 2026-04-06 07:46 by 无际的草原
[考研] 化学调剂 +17 艾志恒 2026-04-03 18/900 2026-04-06 07:10 by jj987
[考研] 求调剂 +10 Hll胡 2026-04-04 10/500 2026-04-05 20:09 by nepu_uu
[考研] 生物学308分求调剂（一志愿华东师大） +8 相信必会光芒万� 2026-04-05 10/500 2026-04-05 12:19 by Hdyxbekcb
[考研] 男生，一志愿沪9生物学071000，初试308求调剂 +3 刘墨墨 2026-04-04 3/150 2026-04-05 08:26 by barlinike
[考研] 22408，264求调剂 +3 ywh729 2026-04-03 4/200 2026-04-04 11:04 by ywh729
[考研] 26调剂 086003 +6 失活的细胞 2026-04-04 6/300 2026-04-04 09:50 by zhangdingwa
[考研] 297求调剂 +11 ljy20040718！ 2026-04-03 13/650 2026-04-04 09:23 by 来看流星雨10
[考研] 265求调剂 +17 林深温澜 2026-04-01 20/1000 2026-04-04 01:09 by userper
[考研] 336求调剂 +8 kiyy 2026-04-01 8/400 2026-04-03 19:41 by lijunpoly
[考研] 一志愿华中农业071010，总分320求调剂 +7 困困困困坤坤 2026-04-02 7/350 2026-04-03 17:26 by Yuena_Wang
[考研] 机械专硕297 +3 Afksy 2026-04-03 3/150 2026-04-03 14:24 by 1753564080
[考研] 266求调剂 +3 08电气工程 2026-04-03 3/150 2026-04-03 14:05 by 1753564080
[考研] 372分材料与化工（085600）一志愿湖南大学求调剂 +5 蓝笺片 2026-04-02 6/300 2026-04-02 21:37 by dongzh2009
[考研] 一志愿北京科技大学材料学硕328分求调剂 +6 1段时间 2026-03-31 7/350 2026-04-02 13:57 by 3041
[考研] 314求调剂 +11 1xiaojun23 2026-03-31 12/600 2026-04-02 12:31 by 1xiaojun23
信息提示
关闭
请填处理意见
关闭