高等工程热力学作业-.doc

高等工程热力学作业（编程）第三章 实际气体状态方程第四章 实际气体导出热力学性质与过程题目：一、用PR方程计算制冷剂R290、R600a和混合制冷剂R290/R600a：50/50wt%的PVT性质。二、用PR方程计算制冷剂R290、R600a和混合制冷剂R290/R600a的导出热力学性质焓和熵。源程序：1、牛顿迭代法求Zfunction Z=newton(A,B,Z)err=1e-6;for n=0:1000 f=Z3-(1-B)*Z2+Z*(A-2*B-3*B2)-(A*B-B2-B3); Z=Z-f/(3*Z2-2*(1-B)*Z+(A-2*B-3*B2); if(abs(f)err) break endend%2、求a、b、Z、v等参数函数function v,Z,a,b,beta=vv(p,T)R=8.31451;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.26992*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); 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 a=a1 a2 a3;b=b1 b2 b3;beta=aa1 aa2 aa3;end%3、余函数法求ar、sr、hrfunction ar,sr,hr=as(p,T)v,Z,a,b,beta=vv(p,T);R=8.31451;for i=1:3;sr(i)=-R*log(v(i)-b(i)/v(i)+beta(i)/(2*sqrt(2)*b(i)*log(v(i)-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);endend%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); Mm=M1*x1+M2*x2; M=M1 M2 Mm;ps=0.015696 0.32979 0.47446; T0=273.15;R=8.31451;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.981*10(-7)*t.2);cps3=inline(-64.79./t+6.798-3.415*10(-3)*t+6.294*10(-7)*t.2);cph1=inline(-95.80+6.945*t-3.597*10(-3)*t.2+7.290*10(-7)*t.3);cph2=inline(-23.91+6.605*t-3.176*10(-3)*t.2+4.981*10(-7)*t.3);cph3=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.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);for i=1: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); h(i)=200*M(i)+hr0(i)+Ih(i)*M(i)-hr(i);endend%5、主程序求v、h、sclearP=input(输入R600a工质压力:P/MPa:n);T=input(输入R600a工质温度:T/K:n);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分别在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=newton(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 break else p1=p1+dp; end end if n=N+1 fprintf(error!) break; else plot(T1,p1/106,r-); hold on; endend grid;title(R290工质p-T相图);xlabel(T/K);ylabel(p/MPa);运行结果：溶液R290/600a分别在P=1atm和10atm下的T-X相图源程序：1、求R290/R600a逸度系数函数function phimix=phimix(type,x1,T,p,Z)R=8.3145; x2=1-x1; 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;b1=0.07780*R*N1(2)/N1(3)/106; k2=0.37464+1.54226*N2(4)-0.26992*N2(4)2;alpha2=(1+k2*(1-(T/N2(2)0.5)2;a2=0.45724*alpha2*R2*N2(2)2/N2(3)/106;b2=0.07780*R*N2(2)/N2(3)/106; 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); else if type=2 bi=b2; sai=2*(x2*a2+x1*0.99*sqrt(a1*a2); endendphimix=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);end%2、1atm下R290/R600a的T-x图clearx1=0:0.01:1;x2=1-x1;t=length(x1); y1=x1;y2=1-y1;P=0.1; n=0;for i=1: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); fug2_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; fug1_v=phimix(1,y1(i),T,P,1.1); 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); end while 1 if abs(sumy-sumy1)/sumy1)1e-4 break end sumy1=sumy; y1(i)=k1*x1(i)/sumy; y2(i)=k2*x2(i)/sumy; fug1_v=phimix(1,y1(i),T,P,1.1); 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); end if abs(sumy-1)=1e-4 q(i)=T; break end T=T+0.01; end R(:,i)=x1(i),y1(i),T;end s=0;for i=1:t if R(3,i)265 s=s+1; L(:,s)=R(:,i); end end L(:,1)=0,0,261;L(:,s+1)=1,1,230.61;plot(L(1,:),L(3,:),r); hold on;plot(L(2,:),L(3,:); legend(泡点线，露点线);xlabel(R290摩尔分数);ylabel(混合工质温度/K);title(p=1atm下，R290/R600a的T-x图);grid on %3、10atm下R290/R600a的T-x图clearx1=0:0.01:1;x2=1-x1;t=length(x1);y1=x1;y2=1-y1;P=1;n=0;for i=1:t T=290; 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); fug2_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)