apolloking

金蟲 (小有名氣)

應(yīng)助: 0 (幼兒園)
金幣: 1083.1
散金: 35
紅花: 2
帖子: 156
在線: 60.5小時
蟲號: 75702
注冊: 2005-06-20
專業(yè): 能源化工

[交流] 【求助】使用Matlab預估動力學方程問題已有4人參與

大家好，問題如下：
使用固定床積分反應(yīng)器來測量甲烷催化燃燒的本征動力學實驗數(shù)據(jù)。
整個計算過程見圖片附件

本征動力學方程為ra=k0*exp（-E/(R*T)）*Ca
ra:反應(yīng)速率 mol/(kgcat*s)
k0:指前因子 m3/(kgcat*s)
E:活化能 kJ/mol
Ca：甲烷摩爾濃度 mol/m3

待定參數(shù)為k0和E

文獻使用了Matlab的半隱式龍格庫塔法和最小二乘法去約束優(yōu)化
最后得出的k0為3.77*10E8
E為25kJ/mol

在積分反應(yīng)器上有;
-ra=dxa/d(w/Fa0)=-k0*exp（-E/(R*T)）*Ca
xa為轉(zhuǎn)化率 w為催化劑的質(zhì)量 Fa0為流速
又因為Ca=C0*(1-xa)
所以上式轉(zhuǎn)化為
dxa/d(w/Fa0)=-k0*exp（-E/(R*T)）*C0（1-xa）
因為實驗所用的甲烷體積濃度為1%，上式可轉(zhuǎn)化為
1%*22.4L/mol

dxa/d(w/Fa0)=-k0*exp（-E/(R*T)）*（1-xa）/2.24
（我的理解是w/Fa0是空速的倒數(shù),但文獻中給的空速值單位是ml/(g*h)，我不知道需不需要在進行單位轉(zhuǎn)化，還是直接取倒數(shù)就行了）

文獻中說，“采用Matlab程序?qū)Ρ菊鲃恿W實驗數(shù)據(jù)進行搜索，選優(yōu)”求的k0和E
按照這個思路，是不是針對每一個溫度，認為溫度值已知，空速和轉(zhuǎn)化率為變量求的一組k0和E。然后將所有溫度得到的k0和E值取平均值？
我自己寫了一個550度的計算程序，但是得不到正確的結(jié)果，希望大家能幫忙看一下。

function jifen
%積分法動力學參數(shù)估算
close all,
clear,
clc,
%實驗數(shù)據(jù)
t=[1/5000 1/10000 1/15000 1/20000 1/25000];
ca=[0.989578 0.972126 0.956323 0.897261 0.870735];
%最小二乘法估計動力學參數(shù)
x0=[0 1];
x=lsqnonlin(@myfun,x0,[],[],optimset('Display','iter'),ca,t);
k0=x(1);
E=x(2);
%給出估算結(jié)果
fprintf('動力學參數(shù)為：\n'),
fprintf('k=%f\ta=%f\n',k0,E),
%繪制反應(yīng)速率擬合曲線
[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('實驗數(shù)據(jù)','模擬曲線')
%實驗導數(shù)與計算導數(shù)之差
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';
%反應(yīng)速率方程
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

已閱回復此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

apolloking

金蟲 (小有名氣)

應(yīng)助: 0 (幼兒園)
金幣: 1083.1
散金: 35
紅花: 2
帖子: 156
在線: 60.5小時
蟲號: 75702
注冊: 2005-06-20
專業(yè): 能源化工

function jifen
%積分法動力學參數(shù)估算
close all,
clear,
clc,
%實驗數(shù)據(jù)
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];
%最小二乘法估計動力學參數(shù)
x0=[3E8 20000];
x=lsqnonlin(@myfun,x0,[],[],optimset('Display','iter'),t,xa,w);
k0=x(1);
e=x(2);
%給出估算結(jié)果
fprintf('動力學參數(shù)為：\n'),
fprintf('k=%f\ta=%f\n',k0,e),
%繪制反應(yīng)速率擬合曲線
[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('實驗數(shù)據(jù)','模擬曲線')
%實驗導數(shù)與計算導數(shù)之差
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';
%反應(yīng)速率方程
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

已閱回復此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

查看全部 14 個回答

change0618

鐵桿木蟲 (著名寫手)

方丈大師

應(yīng)助: 44 (小學生)
金幣: 17724.5
紅花: 17
帖子: 2413
在線: 546.7小時
蟲號: 496517
注冊: 2008-01-19
專業(yè): 化學反應(yīng)工程

★
小木蟲(金幣+0.5):給個紅包，謝謝回帖交流

r=-k0*exp(-e/8.314/(t+273.15))*(1-xa)/2.24;

贊一下(1人)

回復此樓

4樓2009-07-20 13:01:01

已閱回復此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

apolloking

金蟲 (小有名氣)

應(yīng)助: 0 (幼兒園)
金幣: 1083.1
散金: 35
紅花: 2
帖子: 156
在線: 60.5小時
蟲號: 75702
注冊: 2005-06-20
專業(yè): 能源化工

改了這個也不對，運行還是有錯誤

贊一下

回復此樓

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

已閱回復此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

hitzhang

木蟲 (正式寫手)

仿真EPI: 1
應(yīng)助: 0 (幼兒園)
貴賓: 2.15
金幣: 1376.7
散金: 969
紅花: 8
帖子: 862
在線: 226.4小時
蟲號: 390575
注冊: 2007-06-02
性別: GG
專業(yè): 無機非金屬類電介質(zhì)與電解

★
小木蟲(金幣+0.5):給個紅包，謝謝回帖交流

這句x=lsqnonlin(@myfun,x0,[],[],optimset('Display','iter'),t,xa,w);，t, xa, w應(yīng)該是myfun的額外輸入?yún)?shù)，而你的myfun函數(shù)有x, xa, tspan，不懂作者是什么意思

贊一下(1人)

回復此樓

6樓2009-07-20 22:03:01

已閱回復此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

最具人氣熱帖推薦 [查看全部]	作者	回/看	最后發(fā)表

[考研] 一志愿大連理工大學材料求調(diào)劑 +5	Gymno 2026-03-30	5/250	2026-03-31 04:18 by fmesaito
[考研] 一志愿南京航空航天大學，080500材料科學與工程學碩 +4	@taotao 2026-03-30	4/200	2026-03-31 04:04 by zhyzzh
[考研] 277跪求調(diào)劑 +8	1915668 2026-03-27	12/600	2026-03-30 21:01 by dophin1985
[考研] 085600 286分材料求調(diào)劑 +11	麻辣魷魚 2026-03-27	12/600	2026-03-30 19:33 by Wang200018
[考研] 303求調(diào)劑 +7	DLkz1314. 2026-03-30	7/350	2026-03-30 16:05 by shuang5186
[考研] 324求調(diào)劑 +9	hanamiko 2026-03-26	11/550	2026-03-30 14:27 by JourneyLucky
[考研] 本科211生物醫(yī)學工程085409求調(diào)劑339分 +3	里子木yy 2026-03-29	3/150	2026-03-30 13:29 by gyzj2026
[考研] 材料與化工272求調(diào)劑 +21	阿斯蒂芬2004 2026-03-28	21/1050	2026-03-30 10:52 by 晴空210210
[考研] 311求調(diào)劑 +10	lin0039 2026-03-26	10/500	2026-03-30 10:26 by herarysara
[考研] 332求92調(diào)劑 +8	蕉蕉123 2026-03-28	8/400	2026-03-29 10:46 by 周梓丹
[考研] 321求調(diào)劑 +7	璞玉~~ 2026-03-25	8/400	2026-03-29 06:41 by 544594351
[考研] 275求調(diào)劑 +10	jjjjjjjjjjl 2026-03-27	10/500	2026-03-27 23:47 by barnett0632
[考研] 266分求材料化工冶金礦業(yè)等專業(yè)的調(diào)劑 +4	哇呼哼呼哼 2026-03-26	4/200	2026-03-27 17:02 by zhyzzh
[考研] 安徽大學專碩生物與醫(yī)藥專業(yè)(086000)324分，英語已過四六級，六級521，求調(diào)劑 +4	美味可樂雞翅 2026-03-26	4/200	2026-03-27 15:27 by 星空星月
[考研] 274求調(diào)劑 +17	顧九笙要謙虛 2026-03-24	23/1150	2026-03-27 15:16 by caszguilin
[考研] 化學調(diào)劑一志愿上海交通大學336分-本科上海211 +4	小魚愛有機 2026-03-25	4/200	2026-03-26 10:19 by aa331100
[考研] 考研一志愿蘇州大學初始315（英一）求調(diào)劑 +3	sbdksD 2026-03-24	4/200	2026-03-25 18:16 by xcjcqu
[考研] 【2026考研調(diào)劑】制藥工程 284分求相關(guān)專業(yè)調(diào)劑名額 +4	袁奐奐 2026-03-25	8/400	2026-03-25 14:32 by lbsjt
[考研] 一志愿武理085500機械專業(yè)總分300求調(diào)劑 +3	an10101 2026-03-24	7/350	2026-03-25 00:00 by 山鬼0-
[考研] 材料專碩331求調(diào)劑 +4	鮮當牛 2026-03-24	4/200	2026-03-24 15:58 by JourneyLucky