| 5 | 1/1 | 返回列表 |
| 查看: 3206 | 回復(fù): 13 | ||||
| 當(dāng)前只顯示滿足指定條件的回帖,點(diǎn)擊這里查看本話題的所有回帖 | ||||
apolloking金蟲 (小有名氣)
|
[交流]
【求助】使用Matlab預(yù)估動(dòng)力學(xué)方程問(wèn)題 已有4人參與
|
|||
|
大家好,問(wèn)題如下: 使用固定床積分反應(yīng)器來(lái)測(cè)量甲烷催化燃燒的本征動(dòng)力學(xué)實(shí)驗(yàn)數(shù)據(jù)。 整個(gè)計(jì)算過(guò)程見(jiàn)圖片附件 本征動(dòng)力學(xué)方程為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 文獻(xiàn)使用了Matlab的半隱式龍格庫(kù)塔法和最小二乘法去約束優(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為流速 又因?yàn)镃a=C0*(1-xa) 所以上式轉(zhuǎn)化為 dxa/d(w/Fa0)=-k0*exp(-E/(R*T))*C0(1-xa) 因?yàn)閷?shí)驗(yàn)所用的甲烷體積濃度為1%,上式可轉(zhuǎn)化為 1%*22.4L/mol dxa/d(w/Fa0)=-k0*exp(-E/(R*T))*(1-xa)/2.24 (我的理解是w/Fa0是空速的倒數(shù),但文獻(xiàn)中給的空速值單位是ml/(g*h),我不知道需不需要在進(jìn)行單位轉(zhuǎn)化,還是直接取倒數(shù)就行了) 文獻(xiàn)中說(shuō),“采用Matlab程序?qū)Ρ菊鲃?dòng)力學(xué)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行搜索,選優(yōu)”求的k0和E 按照這個(gè)思路,是不是針對(duì)每一個(gè)溫度,認(rèn)為溫度值已知,空速和轉(zhuǎn)化率為變量求的一組k0和E。然后將所有溫度得到的k0和E值取平均值? 我自己寫了一個(gè)550度的計(jì)算程序,但是得不到正確的結(jié)果,希望大家能幫忙看一下。 function jifen %積分法動(dòng)力學(xué)參數(shù)估算 close all, clear, clc, %實(shí)驗(yàn)數(shù)據(jù) t=[1/5000 1/10000 1/15000 1/20000 1/25000]; ca=[0.989578 0.972126 0.956323 0.897261 0.870735]; %最小二乘法估計(jì)動(dòng)力學(xué)參數(shù) x0=[0 1]; x=lsqnonlin(@myfun,x0,[],[],optimset('Display','iter'),ca,t); k0=x(1); E=x(2); %給出估算結(jié)果 fprintf('動(dòng)力學(xué)參數(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í)驗(yàn)數(shù)據(jù)','模擬曲線') %實(shí)驗(yàn)導(dǎo)數(shù)與計(jì)算導(dǎo)數(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 ] |
動(dòng)力學(xué)擬合 |
鐵桿木蟲 (著名寫手)
方丈大師
金蟲 (小有名氣)
|
function jifen %積分法動(dòng)力學(xué)參數(shù)估算 close all, clear, clc, %實(shí)驗(yàn)數(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]; %最小二乘法估計(jì)動(dòng)力學(xué)參數(shù) x0=[3E8 20000]; x=lsqnonlin(@myfun,x0,[],[],optimset('Display','iter'),t,xa,w); k0=x(1); e=x(2); %給出估算結(jié)果 fprintf('動(dòng)力學(xué)參數(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í)驗(yàn)數(shù)據(jù)','模擬曲線') %實(shí)驗(yàn)導(dǎo)數(shù)與計(jì)算導(dǎo)數(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; 錯(cuò)誤提示 ??? 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 ] |
金蟲 (小有名氣)
木蟲 (正式寫手)
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 學(xué)碩274求調(diào)劑 +10 | Li李魚 2026-03-26 | 10/500 |
|
|---|---|---|---|---|
|
[考研] 生物學(xué)學(xué)碩,一志愿湖南大學(xué),初試成績(jī)338 +7 | YYYYYNNNNN 2026-03-26 | 9/450 |
|
|
[考研] 一志愿鄭州大學(xué),080500學(xué)碩,總分317分求調(diào)劑 +10 | 舉個(gè)栗子oi 2026-03-24 | 11/550 |
|
|
[考研] 085601材料工程找調(diào)劑 +17 | oatmealR 2026-03-29 | 18/900 |
|
|
[考研] 材料學(xué)碩333求調(diào)劑 +14 | 北道巷 2026-03-24 | 14/700 |
|
|
[考研] 317分 一志愿南理工材料工程 本科湖工大 求調(diào)劑 +12 | 芋泥小鈴鐺 2026-03-28 | 12/600 |
|
|
[考研] 310求調(diào)劑 +10 | 爭(zhēng)取九點(diǎn)睡 2026-03-30 | 10/500 |
|
|
[考研] 343求調(diào)劑 +6 | 愛(ài)羈絆 2026-03-29 | 6/300 |
|
|
[考研] 調(diào)劑考研 +3 | 王杰一 2026-03-29 | 3/150 |
|
|
[考研] 0703化學(xué)調(diào)劑,求導(dǎo)師收 +9 | 天天好運(yùn)來(lái)上岸?/a> 2026-03-24 | 10/500 |
|
|
[考研] 316求調(diào)劑 +7 | 江辭666 2026-03-26 | 7/350 |
|
|
[考研] 343求調(diào)劑 +5 | 愛(ài)羈絆 2026-03-28 | 5/250 |
|
|
[考研] 339求調(diào)劑,想調(diào)回江蘇 +6 | 烤麥芽 2026-03-27 | 8/400 |
|
|
[材料工程] 一志愿C9材料與化工專業(yè)總分300求調(diào)劑 +8 | 曼111 2026-03-24 | 9/450 |
|
|
[考研] 295求調(diào)劑 +5 | 1428151015 2026-03-27 | 6/300 |
|
|
[考研] 285求調(diào)劑 +4 | AZMK 2026-03-27 | 7/350 |
|
|
[考研] 機(jī)械學(xué)碩310分,數(shù)一英一,一志愿211本科雙非找調(diào)劑信息 +3 | @357 2026-03-25 | 3/150 |
|
|
[考研] 一志愿 南京郵電大學(xué) 288分 材料考研 求調(diào)劑 +3 | jl0720 2026-03-26 | 3/150 |
|
|
[考研] 機(jī)械學(xué)碩總分317求調(diào)劑!。。 +4 | Acaciad 2026-03-25 | 4/200 |
|
|
[考研] 網(wǎng)絡(luò)空間安全0839招調(diào)劑 +4 | w320357296 2026-03-25 | 6/300 |
|