高等工程热力学作业-

高等工程热力学作业（编程）

第三章实际气体状态方程

第四章实际气体导出热力学性质与过程

题目：

一、用PR方程计算制冷剂R290、R600a和混合制冷剂R290/R600a：50/50wt%的PVT性质。

二、用PR方程计算制冷剂R290、R600a和混合制冷剂R290/R600a的导出热力学性质焓和熵。

源程序：

1、牛顿迭代法求Z

function Z=newton(A,B,Z)

err=1e-6;

for n=0:1000

f=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)

break

end

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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;

end

a=[a1 a2 a3];

b=[b1 b2 b3];

beta=[aa1 aa2 aa3];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3、余函数法求ar 、sr 、hr

function [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));

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

4求绝对焓熵（以0℃饱和液体为标准)/(1sat

l ℃,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);

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

5、主程序求v、h、s

clear

P=input('输入R600a工质压力:P/MPa:\n');

T=input('输入R600a工质温度:T/K:\n');

[v]=vv(p,T)