高等工程热力学作业（编程）第三章实际气体状态方程第四章实际气体导出热力学性质与过程题目：一、用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:1000f=Z^3-(1-B)*Z^2+Z*(A-2*B-3*B^2)-(A*B-B^2-B^3);Z=Z-f/(3*Z^2-2*(1-B)*Z+(A-2*B-3*B^2));if(abs(f)<err)breakendend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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*R^2*N1(2)^2/N1(3)/10^6;aa1=0.45724*R^2*N1(2)^2/N1(3)/10^6*2*sqrt(alpha1)*(-k1/(2*sqrt(N1(2)*T))); b1=0.07780*R*N1(2)/N1(3)/10^6;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*R^2*N2(2)^2/N2(3)/10^6;aa2=0.45724*R^2*N2(2)^2/N2(3)/10^6*2*sqrt(alpha2)*(-k2/(2*sqrt(N2(2)*T))); b2=0.07780*R*N2(2)/N2(3)/10^6;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*10^6/(R^2*T^2);B(i)=b(i)*p*10^6/(R*T);Z(i)=newton(A(i),B(i),1);vv(i)=R*T*Z(i)/p/10^6;digits(5);v(i)=vpa(vv(i),5);i=i+1;enda=[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/10^6));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/10^6));hr(i)=ar(i)+T*sr(i)+R*T*(1-Z(i));endend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%4求绝对焓熵（以0℃饱和液体为标准)/(1satl ℃,0K kg kJ s ⋅=）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/(m^3/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*w1^2;alpha1=(1+k1*(1-Tr1^0.5))^2;a1=0.45724*alpha1*(R^2)*(Tc1^2)/Pc1;b1=0.07780*R*Tc1/Pc1;A1=a1*P1/((R^2)*(T1^2));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.89for n=1:Nphi1v=phi(T1,p1,1.1);phi1L=phi(T1,p1,0.001);if abs(phi1v-phi1L)<=errbreakelsep1=p1+dp;endendif n==N+1fprintf('error!')break;elseplot(T1,p1/10^6,'r-');hold on;endendgrid;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*R^2*N1(2)^2/N1(3)/10^6;b1=0.07780*R*N1(2)/N1(3)/10^6;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*R^2*N2(2)^2/N2(3)/10^6;b2=0.07780*R*N2(2)/N2(3)/10^6;a=x1*x1*a1+2*x1*x2*(1-0.01)*sqrt(a1*a2)+x2*x2*a2; b=x1*b1+x2*b2;A=a*p*10^6/(R^2*T^2);B=b*p*10^6/(R*T);Z=newton(A,B,Z);if type==1bi=b1;sai=2*(x1*a1+x2*0.99*sqrt(a1*a2));else if type==2bi=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:tT=220;while 1n=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==1y1(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);endwhile 1if abs((sumy-sumy1)/sumy1)<1e-4breakendsumy1=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);endif abs(sumy-1)<=1e-4q(i)=T;breakendT=T+0.01;endR(:,i)=[x1(i),y1(i),T];ends=0;for i=1:tif R(3,i)<265s=s+1;L(:,s)=R(:,i);endendL(:,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:tT=290;while 1n=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==1y1(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);endwhile 1if abs((sumy-sumy1)/sumy1)<1e-4breakendsumy1=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);endif abs(sumy-1)<=1e-4q(i)=T;breakendT=T+0.01;endR(:,i)=[x1(i),y1(i),T];endplot(R(1,:),R(3,:),'r'); hold on;plot(R(2,:),R(3,:));legend('泡点线'，’露点线’);xlabel('R290摩尔分数');ylabel('混合工质温度/K');title('p=10atm下，R290/R600a的T-x图');grid on %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%运行结果：。