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
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1