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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】使用Matlab预估动力学方程问题
141/2返回列表12下一页
查看: 3175  |  回复: 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下一页
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 085601材料工程专硕求调剂 +5 慕寒mio 2026-03-16 5/250 2026-03-17 21:31 by hmn_wj
[考研] 277调剂 +5 自由煎饼果子 2026-03-16 6/300 2026-03-17 19:26 by 李leezz
[考研] 梁成伟老师课题组欢迎你的加入 +8 一鸭鸭哟 2026-03-14 10/500 2026-03-17 15:07 by 一鸭鸭哟
[考研] 材料与化工专硕调剂 +5 heming3743 2026-03-16 5/250 2026-03-17 14:03 by 勇敢太监王公公
[考研] 285化工学硕求调剂（081700） +9 柴郡猫_ 2026-03-12 9/450 2026-03-17 10:18 by Sammy2
[考研] 药学383 求调剂 +3 药学chy 2026-03-15 4/200 2026-03-16 20:51 by 元子^0^
[考研] 326求调剂 +4 诺贝尔化学奖觊� 2026-03-15 7/350 2026-03-16 17:11 by 诺贝尔化学奖觊�
[考研] 085600调剂 +5 漾漾123sun 2026-03-12 6/300 2026-03-16 15:58 by 漾漾123sun
[考研] 277材料科学与工程080500求调剂 +3 自由煎饼果子 2026-03-16 3/150 2026-03-16 14:10 by 运气yunqi
[考研] 0856求调剂 +3 刘梦微 2026-03-15 3/150 2026-03-16 10:00 by houyaoxu
[考研] 22408总分284求调剂 +3 InAspic 2026-03-13 3/150 2026-03-15 11:10 by zhq0425
[考研] 294求调剂 +3 Zys010410@ 2026-03-13 4/200 2026-03-15 10:59 by zhq0425
[考研] 265求调剂 +4 威化饼07 2026-03-12 4/200 2026-03-14 17:23 by userper
[考研] 337一志愿华南理工0805材料求调剂 +7 mysdl 2026-03-11 9/450 2026-03-13 22:43 by JourneyLucky
[考研] 工科，求调剂 +3 我887 2026-03-11 3/150 2026-03-13 21:39 by JourneyLucky
[考研] 295求调剂 +3 小匕仔汁 2026-03-12 3/150 2026-03-13 15:17 by vgtyfty
[考研] 274求调剂 +3 S.H1 2026-03-12 3/150 2026-03-13 15:15 by JourneyLucky
[考研] 308求调剂 +3 是Lupa啊 2026-03-12 3/150 2026-03-13 14:30 by 求调剂zz
[考研] 321求调剂（食品/专硕） +3 xc321 2026-03-12 6/300 2026-03-13 08:45 by xc321
[考博] 26读博 +4 Rui135246 2026-03-12 10/500 2026-03-13 07:15 by gaobiao
信息提示
关闭
请填处理意见
关闭