版块导航: 正在加载中...

登录注册

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

24小时热门版块排行榜

返回列表

apolloking

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1083.1
散金: 35
红花: 2
帖子: 156
在线: 60.5小时
虫号: 75702
注册: 2005-06-20
专业: 能源化工

[交流] 【求助】使用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 (幼儿园)
金币: 3024.6
红花: 3
帖子: 797
在线: 22.5小时
虫号: 522854
注册: 2008-03-12
性别: MM
专业: 化工系统工程

私下解答。

回复此楼

2楼2009-07-19 14:48:56

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

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

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师

应助: 44 (小学生)
金币: 17724.5
红花: 17
帖子: 2413
在线: 546.7小时
虫号: 496517
注册: 2008-01-19
专业: 化学反应工程

★
小木虫(金币+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

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1083.1
散金: 35
红花: 2
帖子: 156
在线: 60.5小时
虫号: 75702
注册: 2005-06-20
专业: 能源化工

改了这个也不对，运行还是有错误

赞一下

回复此楼

5楼2009-07-20 19:13:06

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

hitzhang

木虫 (正式写手)

仿真EPI: 1
应助: 0 (幼儿园)
贵宾: 2.15
金币: 1376.7
散金: 969
红花: 8
帖子: 862
在线: 226.4小时
虫号: 390575
注册: 2007-06-02
性别: GG
专业: 无机非金属类电介质与电解

★
小木虫(金币+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

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1083.1
散金: 35
红花: 2
帖子: 156
在线: 60.5小时
虫号: 75702
注册: 2005-06-20
专业: 能源化工

实话说这个程序是我从另外一个程序改过来的，我也不是很清楚这个 tspan具体的含义

赞一下

回复此楼

7楼2009-07-20 22:05:23

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

apolloking

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1083.1
散金: 35
红花: 2
帖子: 156
在线: 60.5小时
虫号: 75702
注册: 2005-06-20
专业: 能源化工

snipher950 网友说先固定一个温度，然后输入空速和转化率，但是之后我就不知道怎么办了。我看相关的文献好像都是一次性就做出来一个k0和E的

赞一下

回复此楼

8楼2009-07-20 22:07:53

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

hitzhang

木虫 (正式写手)

仿真EPI: 1
应助: 0 (幼儿园)
贵宾: 2.15
金币: 1376.7
散金: 969
红花: 8
帖子: 862
在线: 226.4小时
虫号: 390575
注册: 2007-06-02
性别: GG
专业: 无机非金属类电介质与电解

是这样，能把你程序里的函数调用关系以及各个函数的作用描述一下吗？

赞一下

回复此楼

9楼2009-07-20 22:21:52

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

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];

%以上是我要输入的参数，一共是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 的主题更新

返回列表

24小时热门版块排行榜

[交流] 【求助】使用Matlab预估动力学方程问题已有4人参与

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

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助: