高等工程热力学作业Word文档下载推荐.docx

1、N1=44.096 369.89 4.2512 0.1521;N2=58.122 407.81 3.6290 0.1840;k1=0.37464+1.54226*N1（4）-0.26992*N1（4）2;alpha1=（1+k1*（1-（T/N1（2）0.5）2;a1=0.45724*alpha1*R2*N1（2）2/N1（3）/106;aa1=0.45724*R2*N1（2）2/N1（3）/106*2*sqrt（alpha1）*（-k1/（2*sqrt（N1（2）*T）;b1=0.07780*R*N1（2）/N1（3）/106; k2=0.37464+1.54226*N2（4）-0.2699

2、2*N2（4）2;alpha2=（1+k2*（1-（T/N2（2）0.5）2;a2=0.45724*alpha2*R2*N2（2）2/N2（3）/106;aa2=0.45724*R2*N2（2）2/N2（3）/106*2*sqrt（alpha2）*（-k2/（2*sqrt（N2（2）*T）;b2=0.07780*R*N2（2）/N2（3）/106;a3=0.25*a1+0.5*（1-0.01）*sqrt（a1*a2）+0.25*a2;aa3=0.25*aa1+0.5*（1-0.01）*1/2/sqrt（a1*a2）*（a1*aa2+a2*aa1）+0.25*aa2;b3=0.5*（b1+b2）

3、;a=a1 a2 a3;b=b1 b2 b3;beta=aa1 aa2 aa3;for i=1:3;A（i）=a（i）*p*106/（R2*T2）;B（i）=b（i）*p*106/（R*T）;Z（i）=newton（A（i）,B（i）,1）;vv（i）=R*T*Z（i）/p/106;digits（5）;v（i）=vpa（vv（i）,5）;i=i+1;end 3、余函数法求ar、sr、hrfunction ar,sr,hr=as（p,T）v,Z,a,b,beta=vv（p,T）;sr（i）=-R*log（v（i）-b（i）/v（i）+beta（i）/（2*sqrt（2）*b（i）*log（v（i

4、）-0.414*b（i）/（v（i）+2.414*b（i）-R*log（v（i）/（R*T/p/106）;ar（i）=R*T*log（v（i）-b（i）/v（i）-a（i）/（2*sqrt（2）*b（i）*log（v（i）-0.414*b（i）/（v（i）+2.414*b（i）+R*T*log（v（i）/（R*T/p/106）;hr（i）=ar（i）+T*sr（i）+R*T*（1-Z（i）;4求绝对焓熵（以0饱和液体为标准）function h,s=hs（p,T） M1=44.096;M2=58.122;x1=（1/M1）/（1/M1+1/M2）;x2=（1/M2）/（1/M1+1/M2）;M

5、m=M1*x1+M2*x2;M=M1 M2 Mm;ps=0.015696 0.32979 0.47446;T0=273.15;c1=-95.80 6.945 -3.597*10（-3） 7.290*10（-7）;c2=-23.91 6.605 -3.176*10（-3） 4.981*10（-7）;c3=-64.79 6.798 -3.415*10（-3） 6.294*10（-7）;cps1=inline（-95.80./t+6.945-3.597*10（-3）*t+7.290*10（-7）*t.2）;cps2=inline（-23.91./t+6.605-3.176*10（-3）*t+4.98

6、1*10（-7）*t.2cps3=inline（-64.79./t+6.798-3.415*10（-3）*t+6.294*10（-7）*t.2cph1=inline（-95.80+6.945*t-3.597*10（-3）*t.2+7.290*10（-7）*t.3cph2=inline（-23.91+6.605*t-3.176*10（-3）*t.2+4.981*10（-7）*t.3cph3=inline（-64.79+6.798*t-3.415*10（-3）*t.2+6.294*10（-7）*t.3 Is1=quad（cps1,273.15,T）/1000; Is2=quad（cps2,273.

7、15,T）/1000; Is3=quad（cps3,273.15,T）/1000;Ih1=quad（cph1,273.15,T）/1000; Ih2=quad（cph2,273.15,T）/1000; Ih3=quad（cph3,273.15,T）/1000;Is=Is1 Is2 Is3; Ih=Ih1 Ih2 Ih3;ar,sr,hr=as（p,T）;3 ar1,sr1,hr1=as（ps（i）,T0）; ar0（i）=ar1（i）; sr0（i）=sr1（i）; hr0（i）=hr1（i）; s（i）=1*M（i）+sr0（i）+Is（i）*M（i）-R*log（p/ps（i）-sr（i）

8、; h（i）=200*M（i）+hr0（i）+Ih（i）*M（i）-hr（i）;5、主程序求v、h、sclearP=input（输入R600a工质压力:P/MPa:nT=input（输入R600a工质温度:T/K:v=vv（p,T） h,s=hs（p,T）R290、R600a、R290/R600a的比体积v/（m3/mol）;R290、R600a、R290/R600a的焓h/（J/mol）;R290、R600a、R290/R600a的熵s/（J/（mol.K）;运行结果：单位制：SI第六章 气液相平衡试用Peng-Robinson方程计算纯质R290 P-T相图和溶液R290/R600a分别在

9、p=1atm和p=10atm下的T-X相图。纯质R290 P-T相图：1、求纯质R290逸度系数函数function phi1=phi（T1,P1,Z）;R=8.3145;M1=44.096e-3;Tc1=369.89;Pc1=4.2512e6;w1=0.1512;Tr1=T1/Tc1;k1=0.37464+1.54226*w1-0.26992*w12;alpha1=（1+k1*（1-Tr10.5）2;a1=0.45724*alpha1*（R2）*（Tc12）/Pc1;b1=0.07780*R*Tc1/Pc1;A1=a1*P1/（R2）*（T12）;B1=b1*P1/（R*T1）;Z=newt

10、on（A1,B1,Z）;phi1=exp（Z-1-log（Z-B1）-A1*log（Z+2.414*B1）/（Z-0.414*B1）/（2*sqrt（2）*B1）;2、主程序p1=3e5;dp=100;N=20000;err=1e-8;for T1=200:0.1:369.89 for n=1:N phi1v=phi（T1,p1,1.1）; phi1L=phi（T1,p1,0.001）; if abs（phi1v-phi1L）=err else p1=p1+dp; if n=N+1 fprintf（error! break; plot（T1,p1/106,r- hold on;grid;tit

11、le（R290工质p-T相图xlabel（T/Kylabel（p/MPa溶液R290/600a分别在P=1atm和10atm下的T-X相图1、求R290/R600a逸度系数函数function phimix=phimix（type,x1,T,p,Z）x2=1-x1;a=x1*x1*a1+2*x1*x2*（1-0.01）*sqrt（a1*a2）+x2*x2*a2;b=x1*b1+x2*b2;A=a*p*106/（R2*T2）;B=b*p*106/（R*T）;Z=newton（A,B,Z）;if type=1 bi=b1; sai=2*（x1*a1+x2*0.99*sqrt（a1*a2）; els

12、e if type=2 bi=b2; sai=2*（x2*a2+x1*0.99*sqrt（a1*a2）;phimix=exp（bi/b*（Z-1）-log（Z-B）-A*（sai/a-bi/b）*log（Z+2.414*B）/（Z-0.414*B）/（2*sqrt（2）*B）;2、1atm下R290/R600a的T-x图x1=0:0.01:1;t=length（x1）;y1=x1;y2=1-y1;P=0.1;n=0;t T=220; while 1 n=n+1; fug1_l=phimix（1,x1（i）,T,P,0.01）; fug1_v=phimix（1,y1（i）,T,P,1.1）; f

13、ug2_l=phimix（2,x1（i）,T,P,0.01）; fug2_v=phimix（2,y1（i）,T,P,1.1）; k1=fug1_l/fug1_v; k2=fug2_l/fug2_v; y1（i）=k1*x1（i）; y2（i）=k2*x2（i）; sumy=y1（i）+y2（i）; sumy1=sumy; if n=1 y1（i）=k1*x1（i）/sumy; y2（i）=k2*x2（i）/sumy; if abs（sumy-sumy1）/sumy1）1e-4 if abs（sumy-1）=1e-4 q（i）=T; T=T+0.01; R（:,i）=x1（i）,y1（i）,T;s=0; if R（3,i）265 s=s+1; L（:,s）=R（:,i）;L（:,1）=0,0,261;,s+1）=1,1,230.61;plot（L（1,:）,L（3,:）,rplot（L（2,:）;legend（泡点线，露点线）;R290摩尔分数混合工质温度/Kp=1atm下，R290/R600a的T-x图grid on 3、10atm下R290/R600a的T-x图P=1; T=290;plot（R（1,:）,R（3,:plot（R（2,:p=10atm下，R290/R600a的T-x图grid on