matlab-物料衡算參考word

1、工程應(yīng)用數(shù)學(xué)作業(yè)5材研1401 642080502011 黃惠洪題目：1.用不同溫度下固體鉛的定壓熱容數(shù)據(jù)，求298K時(shí)的絕對(duì)熵。2.氣液平衡計(jì)算。第一題：實(shí)驗(yàn)數(shù)據(jù)的數(shù)值積分問(wèn)題用不同溫度下固體鉛的定壓熱容數(shù)據(jù)，求298K時(shí)的絕對(duì)熵。數(shù)據(jù)見(jiàn)下表：由熱力學(xué)方法可知，絕對(duì)熵與定壓熱容有如下關(guān)系：本題可用多種方法來(lái)解決，如Simpson數(shù)值積分法，或采用的形式，或先用數(shù)據(jù)點(diǎn)擬合方程得回歸參數(shù)，然后對(duì)回歸方程進(jìn)行積分（回歸方程可用，也可在此基礎(chǔ)上探索更好的回歸方程形式）。1方法與思路采用變型公式：求解給定區(qū)間數(shù)值積分問(wèn)題。2. 過(guò)程與結(jié)果1）積分區(qū)間的確定：根據(jù)網(wǎng)上教程，matlab對(duì)于(-, +)的

2、定義是x=-inf: inf，但是正常情況下是無(wú)需這樣定義的，因系統(tǒng)默認(rèn)變量的取值是從負(fù)無(wú)窮到正無(wú)窮的。2 ) 編寫程序Entropy_Int1.m并保存：這里將最小二乘B樣條擬合法和三次樣條插值函數(shù)計(jì)算作比較，命令行程序如下：function Entropy_Intclear all, clc% 讀入數(shù)據(jù)T = 5. 10. 15. 20. 25. 30. 50. 70. 100. 150. 200. 250. 298.;Cp = 0.305 2.8 7.0 10.8 14.1 16.5 21.4 23.3 24.5 25.4 25.8 26.2 26.5;lnT=log(T);% 用最小二

3、乘B樣條擬合法和三次樣條插值函數(shù)計(jì)算lnTi = linspace(lnT(1),lnT(end),299);knots = 3; K = 3; % 三次B樣條sp = spap2(knots,K,lnT,Cp);cs = csapi(lnT,Cp); % 生成三次樣條插值函數(shù)cs % 繪制濃度擬合曲線Ti = linspace(lnT(1),lnT(end),299);推薦精選Cp_B = fnval(sp,Ti);Cp_c = fnval(cs,Ti);plot(lnT,Cp,ro,Ti,Cp_B,b-,Ti,Cp_c,Bl-.)xlabel(LnT)ylabel(Cp)legend(實(shí)驗(yàn)

4、值,B樣條擬合,立方樣條擬合)% 進(jìn)行數(shù)值積分pp1= fnint(sp);s1= fnval(pp1,lnT)pp2= fnint(cs);s2= fnval(pp2,lnT)end命令窗口輸入Entropy_Int，結(jié)果如下。s1 = Columns 1 through 11 0 1.1069 3.3382 5.9047 8.5613 11.2145 20.9711 28.7027 37.3706 47.4362 54.7432 Columns 12 through 13 60.5284 65.1646s2 = Columns 1 through 11 0 0.6604 2.5965 5.

5、1400 7.9193 10.7146 20.5046 28.0510 36.5924 46.7244 54.0893 Columns 12 through 13 59.8897 64.51873.結(jié)果分析1）查閱CRC化學(xué)和物理手冊(cè)，可得到固體鉛298K時(shí)的絕對(duì)熵為64.8推薦精選，與計(jì)算結(jié)果之間的差別很小。s1=(s1-s)/s=0.5%，s2=(s2-s)/s=0.4%。2）對(duì)比B樣條擬合和立方樣條擬合，可以發(fā)現(xiàn)，立方樣條擬合和實(shí)驗(yàn)值的匹配度更高。最后擬合出的結(jié)果也是立方樣條更接近絕對(duì)熵的理論值。但是B樣條擬合可以適當(dāng)?shù)脑蠼Y(jié)點(diǎn)數(shù)knots來(lái)提高其擬合程度。3）其中，T=0K時(shí)的實(shí)驗(yàn)數(shù)據(jù)

6、由于等于-，如果采用-inf帶入求解會(huì)出現(xiàn)如下警告：Warning: All data points with NaN or Inf as their site will be ignored. 說(shuō)明-將被系統(tǒng)忽略。因此直接將此點(diǎn)舍去來(lái)讓程序可以正常運(yùn)行。便可計(jì)算出結(jié)果。第二題：氣液平衡計(jì)算準(zhǔn)確的汽液平衡數(shù)據(jù)可以用來(lái)計(jì)算液相活度系數(shù)和Gibbs過(guò)量自由能。對(duì)下表的雙組分汽液平衡數(shù)據(jù)，計(jì)算液相活度系數(shù)，并對(duì)結(jié)果進(jìn)行熱力學(xué)的一致性檢驗(yàn)，假定氣相是理想的，即 ，其中 為純組分的飽和蒸汽壓。數(shù)據(jù)來(lái)源：S. Weissman and S. E. Wood, Vapor-Liquid Equilibriu

7、m of Benzene-2,2,4-trimethylpentane Mixtures, J. Chem. Phys., vol.32, 1960, p.1153表51汽液平衡數(shù)據(jù)(液、汽相摩爾分率為組分1：苯)實(shí)驗(yàn)點(diǎn)計(jì)算值液相摩爾分率汽相摩爾分率平衡總壓P(mmHg)-0.00.0178.081.4411.0-0.3650.08190.1869202.741.4151.008-0.3386-0.41180.21920.4065236.861.3431.011-0.2840-0.48370.35840.5509266.041.2501.046-0.1787-0.54210.38310.574

8、8270.731.2421.048-0.1670-0.56950.52560.6786293.361.1581.116-0.0367-0.65180.84780.8741324.661.0231.5080.38760.57150.98720.9863327.391.0001.9680.67670.69471.01.0327.051.02.0260.706 對(duì)低壓汽液平衡，假定氣相為理想是合理的，此時(shí)有：由上表數(shù)據(jù)和公式，結(jié)合由數(shù)據(jù)手冊(cè)上得到的飽和蒸汽壓數(shù)據(jù)或計(jì)算數(shù)據(jù)，可得到數(shù)據(jù)點(diǎn)對(duì)應(yīng)的活度系數(shù)，見(jiàn)上表。對(duì)汽液平衡數(shù)據(jù)進(jìn)行熱力學(xué)一致性檢驗(yàn)要用Gibbs-Duhem方程。對(duì)雙組分情況，Gibbs-

9、Duhem方程的形式是：推薦精選如采用面積檢驗(yàn)法，有：積分的區(qū)間是0,1，因此還缺少兩端點(diǎn)，即無(wú)限稀釋的活度系數(shù)，需采用估計(jì)計(jì)算。無(wú)限稀釋的活度系數(shù)的定義為：和 在分離過(guò)程設(shè)計(jì)和汽液平衡活度系數(shù)模型的參數(shù)估計(jì)中都需要無(wú)限稀釋的活度系數(shù)，但往往很難用實(shí)驗(yàn)測(cè)定，因此它的計(jì)算是很有用的。這里介紹一個(gè)估計(jì)計(jì)算無(wú)限稀釋活度系數(shù)的方法。首先計(jì)算和-，然后將它們和一起對(duì)作圖，可得下圖。分別將-和曲線外推到，然后將它們和，其中用接近零的液相摩爾分率和活度系數(shù)計(jì)算，這樣計(jì)算得到的兩端點(diǎn)的也列于上表，并可作出整個(gè)區(qū)間0,1的待積分圖形。1方法與思路1）首先計(jì)算、-和，對(duì)、-進(jìn)行多項(xiàng)式擬合，得到端點(diǎn)值，從而確定邊值

10、條件。推薦精選2）對(duì)所有數(shù)據(jù)點(diǎn)進(jìn)行三次樣條插值，對(duì)汽液平衡數(shù)據(jù)進(jìn)行熱力學(xué)一致性檢驗(yàn)。2. 過(guò)程與結(jié)果編寫程序VLE.m并保存：function VELclear all,clc%確定邊界條件%導(dǎo)入數(shù)據(jù)，不含端點(diǎn)x1=0.0819,0.2192,0.3584,0.3831,0.5256,0.8478,0.9872;y1=0.1869,0.4065,0.5509,0.5748,0.6786,0.8741,0.9863;p=202.74,236.86,266.04,270.73,293.36,324.66,327.39;p10=327.05;p20=178.08;r1=(p.*y1)./(p10.*

11、x1);r2=(p.*(1-y1)./(p20.*(1-x1);t1=log(r2)./(x1).2);t2=-log(r1)./(1-x1).2);t3=log(r2./r1);x15=0.0819,0.2192,0.3584,0.3831,0.5256;T15=t2(1) t2(2) t2(3) t2(4) t2(5);x67=0.8478,0.9872;T67=t1(6) t1(7);%多項(xiàng)式擬合P1=polyfit(x15,T15,1);P2=polyfit(x67,T67,1);%端值計(jì)算t1end=polyval(P2,1-1e-100);t2int=polyval(P1,1e-1

12、00);%立方樣條差值擬合t1int=log(1)./(1e-100).2);t1=t1 t1end;t1=t1int t1;t2end=-log(1)./(1e-100).2);t2=t2int t2;t2=t2 t2end;t3=t2int t3 t1end;X1=1e-100 x1 1-1e-100;cs=csapi(X1,t3); X2=1e-100 x1;X3=x1 1-1e-100;X1all=linspace(X1(1),X1(end),100);cs_X1=fnval(cs,X1all);Y1=1e-100 x15;Z1=polyval(P1,Y1);Y2=x67 1-1e-1

13、00;推薦精選Z2=polyval(P2,Y2);%端點(diǎn)活度系數(shù)Vint=exp(-polyval(P1,0)Vend=exp(polyval(P2,1)t3int=log(1/Vint)t3end=log(Vend)%繪制氣液平衡曲線與多項(xiàng)式擬合曲線plot(X1all,cs_X1,b-,X1,t3,b*,Y1,Z1,r-o,Y2,Z2,g-o)xlabel(x1)ylabel(%)legend(ln(2/1)立方樣條插值擬合,ln(2/1)實(shí)驗(yàn)數(shù)據(jù)點(diǎn),ln2/x12,-ln1/x22)%進(jìn)行數(shù)值積分pp=fnint(cs);s=fnval(pp,X1)命令窗口輸入VEL，結(jié)果如下：t1 =

14、 0 1.2283 0.2280 0.3479 0.3185 0.3975 0.5716 0.6945 0.7058t2 = -0.3646 -0.4115 -0.4838 -0.5428 -0.5695 -0.6522 -1.0022 -0.7750 0t3 = -0.3646 -0.3386 -0.2840 -0.1787 -0.1700 -0.0370 0.3876 0.6767 0.7058Vint = 1.4400Vend = 2.0255t3int = -0.3646t3end = 0.7058s = 0 -0.0287 -0.0721 -0.1039 -0.1082 -0.124

15、2 -0.0733 0.0002 0.0090推薦精選3.結(jié)果分析1) 將對(duì)應(yīng)計(jì)算結(jié)果填入下表表51汽液平衡數(shù)據(jù)(液、汽相摩爾分率為組分1：苯)實(shí)驗(yàn)點(diǎn)計(jì)算值液相摩爾分率汽相摩爾分率平衡總壓P(mmHg)-0.00.0178.081.4411.0-0.3646-0.36460.08190.1869202.741.4151.008-0.3386-0.41151.22830.21920.4065236.861.3431.011-0.2840-0.48380.22800.35840.5509266.041.2501.046-0.1787-0.54280.34790.38310.5748270.731.2421.048-0.1700-0.56950.31850.52560.6786293.361.1581.116-0.0370-