jeweled-望洋兴叹什么意思
高等工程热力学作业
姓名：
班级：
学号：
XX
XXXX
XXXXXXX
1
第一章
1．用PR方程计算制冷剂R32，R125，和混合制冷剂R410a（R32R125：5050
Wt%）的
pvT性质。
程序说明：进入程序后选择所要计算的制冷剂，输入p，T后可得其
比体积（两相区时分别
输出气液相比体积）
源程序：
#include
#include
#define R 8.31451
double
Newton(double A,double B,double x)
{
double x0;
double f,df;
do
{
x0=x;
f=x*x*x-(1-B)*x*x+(A-3*B*B-2*B)*x-(A*B-B*B-B*B*B);
df=3*x*x-2*(1-B)*x+(A-3*B*B-2*B);
x=x-fdf;
}while(fabs(x-x0)>1e-6);
return
x;
}
void R32(double T,double
p,double *a,double *b,double *M)
{
double
Tc,pc,w,k,a1,Tr;
*M=52.024e-3;
Tc=351.255;
pc=5780000;
w=0.277;
k=0.37464+1.54226*w-0.26992*w*w;
Tr=TTc;
a1=pow(1+k*(1-pow(Tr,0.5)),2);
*a=0.45727*a1*R*R*Tc*Tcpc;
*b=0.07780*R*Tcpc;
}
void R125(double T,double p,double
*a,double *b,double *M)
{
double
Tc,pc,w,k,a1,Tr;
2
*M=120.03e-3;
Tc=339.45;
pc=3630600;
w=0.299;
k=0.37464+1.54226*w-0.26992*w*w;
Tr=TTc;
a1=pow(1+k*(1-pow(Tr,0.5)),2);
*a=0.45727*a1*R*R*Tc*Tcpc;
*b=0.07780*R*Tcpc;
}
void R410a(double T,double p,double
*a,double *b,double *M)
{
double
a1,a2,b1,b2,x1,x2,k12,M1,M2;
k12=0.01;
R32(T,p,&a1,&b1,&M1);
R125(T,p,&a2,&b2,&M2);
x1=1(1+M1M2);
x2=1(1+M2M1);
*a=x1*x1*a1+x2*x2*a2+2*x1*x2*(1-k12)*sqrt(a1*a2);
*b=x1*b1+x2*b2;
*M=x1*M1+x2*M2;
}
void main()
{
double M,T,a,b,p,A,B;
int i;
N1:
cout<<
cin>>i;
if(i!=1&&i!=2&&i!=3)
{
cout<<
goto
N1;
}
cout<<
cin>>T;
cout<<
cin>>p;
p=p*1e6;
if(i==1)
{
R32(T,p,&a,&b,&M);
}
else
if(i==2)
3
{
R125(T,p,&a,&b,&M);
}
else if(i==3)
{
R410a(T,p,&a,&b,&M);
}
A=a*p(R*R*T*T);
B=b*p(R*T);
double
z1=Newton(A,B,1000);
double
z2=Newton(A,B,0.001);
if(fabs(z1-z2)<1e-4)
{
double v1=z1*R*TpM;
cout<<单位比体积为：
}
else
{
double
v1=z1*R*TpM;
double v2=z2*R*TpM;
cout<<气体单位比体积为：
cout<<液体单位比体积为：
}
}
（一）
please enter 1(R32),2(R125)or3(R410a)
1
please enter T(K)
300
please
enter p(Mpa)
1.7499
气体单位比体积为：0.0214228m^3kg
液体单位比体积为：0.00122247m^3kg
Press any key to
continue
（二）
please enter
1(R32),2(R125)or3(R410a)
2
please enter
T(K)
300
please enter p(Mpa)
1.7499
气体单位比体积为：0.00762949m^3kg
液体单位比体积为：0.000859649m^3kg
4
Press any key to continue
（三）
please enter 1(R32),2(R125)or3(R410a)
3
please enter T(K)
300
please enter
p(Mpa)
1.7499
气体单位比体积为：0.0147329m^3kg
液体单位比体积为：0.00105639m^3kg
Press any key to
continue
第二章
1. 利用热力学普遍关系式推导：
??
?s
??
?s
??
?s
??
?v
?
?
?T
?
?
?
??
?
??
???
p
?
?T
?
v
?
?v
?
T
?
?T
?
p
证明：
?c
p
?c
v
?R
g
?
R
g
c
p
R
g
c
v
R
g
R
g
p
?
v
?
v
?
p
?
v
?
p
由理想气体状态方程得:
?
?
?v
?
R
g
?
?p
?
R
g
?
?T
?
?
?
，
?
p
p
?
?T
?
?
?
v
v
??
?v
?
代入可得：
?
?T
?
?
?<
br>c
p
p
v
?
?
?
?p
?
?
?T
?
?
?
c
v
v
p
?
?
?
?p
??
?v
?
?
?T
?
?
?
?
v
?
?T
?
?
p
根据热力状态基本表达式得：
c
p
v
?
1
?
，
c
v
??
1
?
?T
?
p??T
?
?
?p
?
?
??
s
?
?v
?
s
?
?
?v
?
?
?p
?<
br>代入得：
?
?T
?
?
p
??
?
??
?
?T
?
v
?
?p
??
?v
?<
br>?T
??
?T
?
?
?
?
?T
??
?
?
?
?
v
?
?
?p
?<
br>?
?
s
?
?v
?
?T
?
p
?
s
利用麦克斯韦关系式:
?
?
?p
???
?s
?
?
??
?
?T
?
?
?
?v
?
?
T
?
?v
?
v
?
?
s
?
?
??
?
?
?T
?
?
T
?
?p
?
v
5
得：
?1
?
?
?1
?
?p
??
?v
?
?T
?
?
?p
?
?
?
? p
?
?
?
?
?
?T
??
?v
?< br>?
?
?
?T
?
?
?
?
?
s
?
?s
?
?
T
?
?v
?
?
V
?
?T
?
p
?
?
?
s
?s
?
?
?
T
带入倒数关系
?
?
?v
?
?
?1
?
，
?< br>?p
?
?
?1
?
?T
?
?
p
?
?p
?
?
?
?T
?
?
v
?< br>?
?s
?
?
?
?v
?
?
T
?
?s
?
T
?
?
?s
??
?s
?
?T
?
?
?
?
?
?
?
?
?
?p
?
?
?
?
?v
?
p?
?T
?
v
?
?T
?
?
?T
?
v
??
p
由麦克斯韦关系：
?
?
?s
??
?
?
?v
?
??
?
p
?
?
T
?
?T
?
v
得
?
?
?s?
?
?T
?
?
?
?
?
?s
?
?
?
?
?
?s
?
?
?
?
?
?v
?
?
p
?
?T
?
v
??v
?
T
?
?T
?
p
2. 利用热力学普遍关系式推导第三dh和ds方程：
dh?
?
?
v?c
?
?T
?
?
?
?T
?
?
?p
?
?
?
dp?c
?
p
?
p
??
dv
v
?
?
?v
?
p
ds?
cv
?
?T
?
c
p
?
?T
?
T
?
?
?p
?
?
dp?
??
dv
v
T
?
?v
?
p
解：若状态方程以p，v为独立变 量
(1)比焓的变化为:
dh?
?
?
?h< br>?
?
?p
?
?
dp?
?
?
?h?
v
?
?v
?
?
dv
p
式中:
?
?h
?
?
?
?h
?< br>??
?
?p
?
?
??
?T
?
??p
?
?
?
?
?T
?
v
?
?
?
T
?
?
?s
?
?
?v
???
?
?
?T
?
?
?c
v
??
?v
v
?
?
?p
?
?
?
?T
?< br>v
?
?T
?
v
?
?
?p
?
v
?
?p
?
v
?
?T
?
?
v
?
?
?
?h
?
?
?h
?
?
?T
?
?
p
?
?T
?
?v
?
?
??c
?
?
p
?
p
?
?
?p
?
?
?v
?
p
?
?T
?
?
p
代回得：
6
dh?
?
?
v
?c
?
?T
?
?
p
?
?
?p
?<
br>?
?
?
dp?c
?
?T
?
p
?v
?
?
?v
?
?
dv
p
(2)比熵的变化
ds?
?
?
?s
?
?
dp?<
br>?
?
?s
?
?
?
p
?
v
?
?v
?
?
dv
p
?
2.2
?
?
?s
?
?
??s
?
?
?
?T
?
?
v
c
v
?
?
?
?p
?
?
?
T
?
v
?
?
?p
?
?
T
?
?
?p?
?
v
?
?T
?
?
v
?
?s
?
?
?
?s
?
?
?
?p
?
?
?
?
?T
?
?
v
c
v
?
?T
?
v
?
?
?p
?
?T
?
?
?p
?
?
v
?
?T
?
?
v
代回得：
ds?
c
v
?
?
T
?
c
p
?
?T
?
T
?
?
?p
?
?
dp?
T
?
?
?v
?
?
dv
v
p
3.
推导PR方程的导出热力性质余函数
h
r
、
s
r
。
解：PR方程：
p?
RTa
v?b
?
v
?
v?b
?
?b
?
v?b
?
a
v
?
RT
?
v
r
?
?
?
?
?p?
v
?
?
dv?RTln
v
?
?
?
v
?
RTRTa
?
v
?
?
?<
br>v?b
?
v
?
v
?
v?b
?
?b<
br>?
v?b
?
?
dv?RTln
?
v
?
?RT
?
v
v
?
ln
?
v?b
?
?
?
?
?RT
?
lnv
?
v
a
?
?
?
?
v
?
v?b
?
?b<
br>?
v?b
?
dv?RTln
v
v
?
?RTl
n
v?b
v
?
?
v
a
?
v
?v?b
?
?b
?
v?b
?
dv?RTln
v<
br>v
?
上式中：
?
v
a
?
v
?
v?b
?
?b
?
v?b
?
dv?a
?
v
1
?
?
22b
?
1
?
v?0.414b
?
1
?
v?2.414b
?
?
dv
7
?
a
?
v
1
22b
?
?
?
?
v?0.414b
dv
?
?
v
1
?
?
v?0.414b
dv
?<
br>?
?
av?0.414b
22b
ln
v?2.414b
a
v?bav?0.41
r
?RTln
v?
22b
ln
b4
v?2.41b4
?RTl
v
v
n
?
所以：
s
?<
br>?
?
v?b
?
v?
r
?s?s??
?T?
a?a
?
v
??Rln
v
?
22
l
n
0.41b4
v?2.41b4
?Rl
v
b
v
n
?
其中：
?
?
?a
?T
。
h
?r
?h?h?a
r
??T
r
?s
?
p
?
?v
?
v
?
T
?<
br>?a
22b
ln
v?0.414b
?
v
?
a
v
v?2.414b
?RT
?
?
v?b
?1
??
?
v
?
v?b
?
?b
?
v?b?
?
T
?
?a
22bln
v?0.414b
v?2.414b
?RT
bav
v?b<
br>?
v
?
v?b
?
?b
?
v?b
?<
br>
4．用PR方程计算工质R32，R125，和混合工R32R125的导出热力性
质焓和熵。
（一）计算R32，R125的焓熵值
源程序
#include
#include
#define R 8.31
double get_a(double
w,double T,double Tc,double pc)
{
double
k=0.37464+1.54226*w-0.26992*w*w;
double
ar=(1+k*(1-sqrt(TTc)))*(1+k*(1-sqrt(TTc)));
double a=0.45724*ar*R*R*Tc*Tcpc;
return a;
}
double get_b(double Tc,double pc)
{
double b=0.0778*R*Tcpc;
return b;
}
8
double
Newton(double A,double B,double x)
{
double x0;
double f,df;
do
{
x0=x;
f=x*x*x-(1-B)*x*x+(A-3*B*B-2*B)*x-(A*B-B*B-B*B*B);
df=3*x*x-2*(1-B)*x+(A-3*B*B-2*B);
x=x-fdf;
}while(fabs(x-x0)>1e-6);
return
x;
}
double get_ar(double T,double
v,double vv,double a,double b)
{
double a
r=R*T*log((v-b)v)-a*log((v-0.414*b)(v+2.414*b))(2*
1.414*b)+R*T*log(vvv);
return ar;
}
double get_sr(double T,double v,double
vv,double a,double b,double bb)
{
double
sr=-1*R*log((v-b)v)+bb*log((v-0.414*b)(v+2.414*b))
(2*1.414*b)-R*log(vvv);
return sr;
}
void get_hr(double Tc,double pc,double
w,double T,double p,double *hr,double *sr,double
zz)
{
double a=get_a(w,T,Tc,pc);
double b=get_b(Tc,pc);
double bb=(get_a(w,T+0
.25,Tc,pc)-get_a(w,T-0.25,Tc,pc))0.5;
double
A=a*p(R*R*T*T);
double B=b*p(R*T);
double z=Newton(A,B,zz);
double v=z*R*Tp;
double vv=R*Tp;
double
ar=get_ar(T,v,vv,a,b);
*sr=get_sr(T,v,vv,a,b,bb);
*hr=ar+T*(*sr)+R*T*(1-z);
}
void main()
{
int i;
double M;
double
w;
9
double h0=200.0;
double s0=1.0;
double T0;
double p0;
double c0,c1,c2,c3;
double T,p,Tc,pc;
double hr0,sr0,hrv,hrl,srv,srl;
double pM[2]={52.024,120.03};
double
pTc[2]={351.255,339.45};
double
ppc[2]={5780000,3630600};
double
pT0[2]={273.15,273.15};
double
pp0[2]={813100,617500};
double
pw[2]={0.277,0.299};
double
pc0[2]={4.424901,2.838680};
double
pc1[2]={-2.661170,11.581633};
double
pc2[2]={5.580232,-1.704482};
double
pc3[2]={-1.680558,-0.266732};
N1:
cout<<
cin>>i;
if(i==1||i==2)
{
M=pM[i-1];
Tc=pTc[i-1];
pc=ppc[i-1];
T0=pT0[i-1];
p0=pp0[i-1];
w=pw[i-1];
c0=pc0[i-1];
c1=pc1[i-1];
c2=pc2[i-1];
c3=pc3[i-1];
}
else
{
cout<<
goto N1;
}
cout<<
cin>>T;
cout<<
cin>>p;
10
p=p*1e6;
get_hr(Tc,pc,w,T0,p0,&hr0,&sr0,0.001);
get_hr(Tc,pc,w,T,p,&hrv,&srv,1.1);
get_hr(Tc,pc,w,T,p,&hrl,&srl,0.001);
if(fabs(hrv-hrl)<1e-4)
{
double
h=h0*M+hr0-hrv+R*(c0*(T-T0)+c12Tc*(pow(T,2)
-pow(T0,2))+c23pow(Tc,2)*(pow(T,3)-pow(T0,3))+c34p
ow(Tc,3)
*(pow(T,4)-pow(T0,4)));
double s=s0*M+sr0-srv-R*log(pp0)+R*(c0*log(TT0)+c1
*(T-T0)Tc
+c22pow(Tc,2)*(pow(T,2)-pow(T0,2)
)+c33pow(Tc,3)*(pow(T,3)-pow(T0,3)));
h=hM;
s=sM;
cout<<
cout<<
}
else
{
double
h=h0*M+hr0-hrv+R*(c0*(T-T0)+c12Tc*(pow(T,2)
-pow(T0,2))+c23pow(Tc,2)*(pow(T,3)-pow(T0,3))+c34p
ow(Tc,3)
*(pow(T,4)-pow(T0,4)));
double s=s0*M+sr0-srv-R*log(pp0)+R*(c0*log(TT0)+c1
*(T-T0)Tc
+c22pow(Tc,2)*(pow(T,2)-pow(T0,2)
)+c33pow(Tc,3)*(pow(T,3)-pow(T0,3)));
h=hM;
s=sM;
cout<<气相
cout<<
cout<<
h=h0*M+hr0-hrl+R*(c0*(T-T0)+c12Tc*(pow(T,2)
-pow(T0,2))+c23pow(Tc,2)*(pow(T,3)-pow(T0,3))+c34p
ow(Tc,3)
*(pow(T,4)-pow(T0,4)));
s=s0
*M+sr0-srl-R*log(pp0)+R*(c0*log(TT0)+c1*(T-T0)Tc
+c22pow(Tc,2)*(pow(T,2)-pow(T0,2))+c33pow(T
c,3)*(pow(T,3)-pow(T0,3)));
h=hM;
s=sM;
cout<<液相
cout<<
cout<<
}
}
please enter 1(R32),2(R125)
1
please enter T(K)
11
300
please enter p(Mpa)
1.7749
气相
h=533.901kJkg
s=2.11853kJ(kg*K)
液相
h=254.453kJkg
s=1.18608kJ(kg*K)
Press
any key to continue
（二）计算混合工质的焓熵值
源程序
#include
#include
#define R 8.31451
#define k12 0.01
double get_a(double w1,double w2,double
T,double Tc1,double pc1,double Tc2,
double
pc2,double x1,double x2)
{
double
k=0.37464+1.54226*w1-0.26992*w1*w1;
double
ar=(1+k*(1-sqrt(TTc1)))*(1+k*(1-sqrt(TTc1)));
double a1=0.45724*ar*R*R*Tc1*Tc1pc1;
k=0.37464+1.54226*w2-0.26992*w2*w2;
ar=(1+k*(1-sqrt(TTc2)))*(1+k*(1-sqrt(TTc2)));
double a2=0.45724*ar*R*R*Tc2*Tc2pc2;
double
a=2*x1*x2*(1-k12)*sqrt(a1*a2)+x1*x1*a1+x2*x2*a2;
return a;
}
double get_b(double
Tc1,double pc1,double Tc2,double pc2,double
x1,double x2)
{
double
b1=0.0778*R*Tc1pc1;
double
b2=0.0778*R*Tc2pc2;
double b=x1*b1+x2*b2;
return b;
}
double get_bb(double
w1,double w2,double T,double Tc1,double pc1,double
Tc2,
double pc2,double x1,double x2)
{
double
a1=get_a(w1,w2,T+0.25,Tc1,pc1,Tc2,pc2,x1,x2);
double
a2=get_a(w1,w2,T-0.25,Tc1,pc1,Tc2,pc2,x1,x2);
12
double bb=(a1-a2)0.5;
return bb;
}
double Newton(double
A,double B,double x)
{
double x0;
double f,df;
do
{
x0=x;
f=x*x*x-(1-B)*x*x+(A-3*B*B-2*B)*x-(A*B-B*B-B*B*B);
df=3*x*x-2*(1-B)*x+(A-3*B*B-2*B);
x=x-fdf;
}while(fabs(x-x0)>1e-6);
return
x;
}
double get_ar(double T,double
v,double vv,double a,double b)
{
double a
r=R*T*log((v-b)v)-a*log((v-0.414*b)(v+2.414*b))(2*
1.414*b)+R*T*log(vvv);
return ar;
}
double get_sr(double T,double v,double
vv,double a,double b,double bb)
{
double
sr=-1*R*log((v-b)v)+bb*log((v-0.414*b)(v+2.414*b))
(2*1.414*b)-R*log(vvv);
return sr;
}
void main()
{
double M1=52.024;
double Tc1=351.255;
double pc1=5780000;
double Ts1=273.15;
double ps1=813100;
double w1=0.277;
double x1;
double
c01=4.424901;
double c11=-2.661170;
double c21=5.580232;
double c31=-1.680558;
double M2=120.03;
13
double Tc2=339.45;
double
pc2=3630600;
double Ts2=273.15;
double
ps2=671500;
double w2=0.299;
double x2;
double c02=2.838680;
double c12=11.581633;
double c22=-1.704482;
double
c32=-0.266732;
double T0=273.15;
double p0=799100;
double
hr0,sr0,hrv,srv,hrl,srl,h,s;
double p,T;
x1=(M2)(M1+M2);
x2=1-x1;
double
M=M1*x1+M2*x2;
double
a=get_a(w1,w2,T0,Tc1,pc1,Tc2,pc2,x1,x2);
double b=get_b(Tc1,pc1,Tc2,pc2,x1,x2);
double
bb=get_bb(w1,w2,T0,Tc1,pc1,Tc2,pc2,x1,x2);
double A=a*p0(R*R*T0*T0);
double
B=b*p0(R*T0);
double z=Newton(A,B,0.001);
double v=z*R*T0p0;
double vv=R*T0p0;
double ar=get_ar(T0,v,vv,a,b);
sr0=get_sr(T0,v,vv,a,b,bb);
hr0=ar+T0*sr0+R*T0*(1-z);
cout<<
cin>>T;
cout<<
cin>>p;
p=p*1e6;
a=get_a(w1,w2,T,Tc1,pc1,Tc2,pc2,x1,x2);
bb=get_bb(w1,w2,T,Tc1,pc1,Tc2,pc2,x1,x2);
A=a*p(R*R*T*T);
B=b*p(R*T);
z=Newton(A,B,0.001);
14
v=z*R*Tp;
vv=R*Tp;
ar=get_ar(T,v,vv,a,b);
srv=get_sr(T,v,vv,a,b,bb);
hrv=ar+T*srv+R*T*(1-z);
z=Newton(A,B,1.1);
v=z*R*Tp;
vv=R*Tp;
ar=get_ar(T,v,vv,a,b);
srl=get_sr(T,v,vv,a,b,bb);
hrl=ar+T*srl+R*T*(1-z);
double
dh1=R*(c01*(T-T0)+c112Tc1*(pow(T,2)-pow(T0,2))+
c213pow(Tc1,2)*(pow(T,3)-pow(T0,3))+c314pow(T
c1,3)*(pow(T,4)-pow(T0,4)));
double
dh2=R*(c02*(T-T0)+c122Tc2*(pow(T,2)-pow(T0,2))+
c223pow(Tc2,2)*(pow(T,3)-pow(T0,3))+c324pow(T
c2,3)*(pow(T,4)-pow(T0,4)));
double
dh=dh1*x1+dh2*x2;
double
ds1=-R*log(pp0)+R*(c01*log(TT0)+c11*(T-T0)Tc1+
c212pow(Tc1,2)*(pow(T,2)-pow(T0,2))+c313pow(T
c1,3)*(pow(T,3)-pow(T0,3)));
double
ds2=-R*log(pp0)+R*(c02*log(TT0)+c12*(T-T0)Tc2+
c222pow(Tc2,2)*(pow(T,2)-pow(T0,2))+c323pow(T
c2,3)*(pow(T,3)-pow(T0,3)));
double
ds=ds1*x1+ds2*x2;
if(fabs(hrv-hrl)<1e-4)
{
h=200+(hr0-hrv+dh)M;
s=1.0+(sr0-srv+ds)M;
cout<<
cout<<
}
else
{
h=200+(hr0-hrv+dh)M;
s=1.0+(sr0-srv+ds)M;
cout<<气相
cout<<
cout<<
h=200+(hr0-hrl+dh)M;
s=1.0+(sr0-srl+ds)M;
cout<<液相
cout<<
cout<<
}
15
}
please
enter T(K)
300
please enter p(Mpa)
1.7406
气相
h=246.288kJkg
s=1.1582kJ(kg*K)
液相
h=433.398kJkg
s=1.78313kJ(kg*K)
Press any key to
continue
第三章
1.
溶液的化学势、逸度与逸度系数定义是什么？物理意义是什么？
(1)化学势的定义是：溶液系统在某
些特定的条件下，作为该系统的特征函数的广延量参数
对该组成摩尔数的变化率。i组成化学势的定义式
为：
?
i
?
?
?
?U
?
?<
br>?n
?
i
?
S,V,n
j
?
j?i
?
物理意义为：化学势差是传递质量的驱动力，是由于i组成摩尔数改变而引起的相应特
征广延性质改变的势位。
(2)溶液逸度及逸度系数：
对于变成分溶液系统，i组成的逸度及逸度系数由下式定义：
?
?
d
?
i
?dG
i
?R
m
Tdlnf
i
T<
br>?
?
?
?c(T)
?
?G
i
?R
m
Tlnf
iii
?
T
?
?
?
lim(f
i
?
p?0
xp
)?1
i
?
??
f
i
?
?
?
i
?
x
i
p
?
?
??
??
上面两式中：
f
?
i
— 表示溶液中i组成的逸度；
?
—
表示溶液中i组成的逸度系数。
?
i
?
是溶液中组成的热力性质，并且是强
度性质，逸度和压力一样，表征物理意义为：
f
i
了物质的逃逸势，它表示组成的假想
分压力。
?
中
xp
为
i
组成的分压力，随着实际气体接近
理想气体，在数值
f
?
上接近于
p
，所
?
i
i
i
16
以
?
?
i
是度量气体非理想性的标尺之一。
?
?
i
是无量纲参数。
3．结合专业选用适当的方程并推导其i组成逸度的数学表达式。
解：选用PR方程
由：
dA??SdT?pdV?
?
r
?
i
dn
i
i?1
得：
?
?
?A
?
?<
br>?V
?
?
??p
T,n
A
?
T,
V,n
?
?A
*
?
T,V
*
,n
?
??
?
V
V
*
pdV
p
*
?
0
而
V
*
??
，对上式右边加减
?
V
nR
m
T
V
*
V
dV
，得
A
?T,V,n
?
?A
*
?
T,V
*
,n
?
?
?
V
?
nRT
?
V
nR
m<
br>T
V
*
?
m
?
V
?p
?
?
dV?
?
V
*
V
dV
当
T,V,n
j
为常数时，上式对
n
j
求导，则
?
?
?A
?
V
?
RT
?
?
?n
?
?
?
?
?A
*
?
?
?<
br>?
m
?
?
?
?p
?
?
?
d
V?R
V
m
Tln
i
?
T,V,n
?
?n
i
?
?
V
*
T,V
*
,n
j?
?
V
?
?n
i
?
T,V,n
j?
?
V
*
又因为
?
?
?A
?
?
?n
?
?
?
?
?G
?
?
?n
?
?G
f
i
?
T,V,n
i
?
T
,p,n
j
?
?A
*
?
*
?
?
?n
?
?
?
?
?G
?
?
?G*
f
i
?
T,V
*
,n
j
?
?n
i
?
T,p
*
,n
j
得
*
?
G?
?
V
R
m
T
?
?p
??
i
?G
i
V
*
?
?
?
V<
br>
?
V
?
?
?
?n
?
?dV?R
m
Tln
*
i
?
T,V,n
j
?
V
左侧可以表示成：
??
G
i
?G
*
i
?R
m
Tln
f
i
?
?R
m
Tln
f
i
f
*
xp
*
i
i
RTln
V
?
Z
n
R
m
T
nR
m
T
?
p
*
m
V
*
?Rm
Tln
?
?
p
?
p
*
?
?
?R
m
TlnZ?R
m
Tln
p
代回：
17
R
m
Tln<
br>?
?
?
又
p?
?
?
V
?
?
?p
?
R
m
T
?
?
?
?
dV?R
m
TlnZ
?
?
V
??
??n
i
?
T,V,n
j
??
RTa
?
v?bv
?
v?b
?
?b
?
v?b
?<
br>rrrr
V
其中
v?,n?
?
n
i
,a?<
br>??
x
i
x
j
(1?k
ij
)a
i
a
j
,b?
?
x
i
b
i
n
i?1i?1j?1i?1
?
?p
?
?p?v?x
i<
br>?p?a?x
i
?p?b?x
i
???
则：
?
?
?n?v?x?n?a?x?n?b?x?n
iiiiii
?
i
?
T,V,n
j
带回上式进行积分得：
?
2
?
x
i
a
ij
?
b
i
b
i
?
A
?
j
Z?2.414B
ln
?
i
?<
br>?
Z?1
?
?ln
?
Z?B
?
???ln<
br>
??
bab
?
Z?0.414B
22B
?
??
?
式中
A?
pV
m
ap
bp
B?Z?
22
R
m
TR
m
T
R
m
T
第四章
1．用Peng-Robinson方程计算纯质R32的p-
T图和溶液R32R125分别在p=1atm和p=12atm
下的T-x图。
（一）R32的p-T图
源程序：
#include
#include
double Newton(double
A,double B,double x)
{
double x0;
double f,df;
do
{
x0=x;
f=x*x*x-(1-B)*x*x+(A-3*B*B-2*B)*x-(A*B-B*B-B*B*B);
df=3*x*x-2*(1-B)*x+(A-3*B*B-2*B);
x=x-fdf;
}while(fabs(x-x0)>1e-6);
return
x;
}
void main()
{
18
double Tc=351.255;
double
pc=5780000;
double w=0.277;
double
R=8.31451;
double
T,p=1e6,a,b,A,B,faiv,fail,z;
cout<<
p(Mpa)
for(int i=0;i<=10;i++)
{
T=270;
do
{
double
k=0.37464+1.54226*w-0.26992*w*w;
double
ar=(1+k*(1-sqrt(TTc)))*(1+k*(1-sqrt(TTc)));
a=0.45724*ar*R*R*Tc*Tcpc;
b=0.0778*R*Tcpc;
A=a*p(R*R*T*T);
B=b*p(R*T);
z=Newton(A,B,1.1);
faiv=(z-1)-log(z-
B)-A(2*1.414*B)*log((z+2.414*B)(z-0.414*B));
faiv=exp(faiv);
z=Newton(A,B,0.001);
fail=(z-1)-log(z-B)-A(
2*1.414*B)*log((z+2.414*B)(z-0.414*B));
fail=exp(fail);
T=T+0.001;
}while(fabs(faiv-fail)>1e-4);
cout<
}
}
T(K) p(Mpa)
279.619 1
282.759 1.1
285.691 1.2
288.443 1.3
291.041 1.4
293.503 1.5
295.844 1.6
298.078 1.7
300.215 1.8
302.265 1.9
304.236 2
Press any key
to continue
19
使用matlab拟和图像为
（二）溶液R32R125分别在p=1atm和p=12atm下的T-x图
源程序：
#include
#include
#include
#define R 8.31451
double Newton(double A,double B,double x)
{
double x0;
double f,df;
do
{
x0=x;
f=x*x*x-(1-B)*x*x+(A-3*B*B-2*B)*x-(A*B-B*B-B*B*B);
df=3*x*x-2*(1-B)*x+(A-3*B*B-2*B);
x=x-fdf;
}while(fabs(x-x0)>1e-6);
return
x;
}
void main( )
20
{
int k=1;
double ysum;
double M1=52.024;
double Tc1=351.255;
double pc1=5780000;
double w1=0.277;
double M2=120.03;
double Tc2=339.45;
double pc2=3630600;
double w2=0.299;
double k12=0.01;
double p;
double P[2]={1.013e5,12*1.013e5};
double
TT[2]={200,260};
for(int j=0;j<=1;j++)
{
p=P[j];
cout<<
cout<< y(R32)
T
for(int i=0;i<=100;i++)
{
double x1=(double)i100,x2=1-x1;
double
y1=0.5,y2=0.5;
double T=TT[j];
N1:
double k1=0.37464+1.54226*w1-0.26992*w1*w1;
double Tr1=TTc1;
double e1=(1.0+k1*(1-sqrt(
Tr1)))*(1.0+k1*(1.0-sqrt(Tr1)));
double
a1=0.45724*e1*R*R*Tc1*Tc1pc1;
double
b1=0.0778*R*Tc1pc1;
double
k2=0.37464+1.54226*w2-0.26992*w2*w2;
double
Tr2=TTc2;
double e2=(1.0+k2*(1-sqrt(Tr2)))*
(1.0+k2*(1.0-sqrt(Tr2)));
double
a2=0.45724*e2*R*R*Tc2*Tc2pc2;
double
b2=0.0778*R*Tc2pc2;
double
a12=sqrt(a1*a2)*(1.0-k12);
double
ax=2.0*x1*x2*a12+x1*x1*a1+x2*x2*a2;
double
bx=b1*x1+b2*x2;
21
double A=(ax*p)(R*R*T*T);
double
B=(bx*p)(R*T);
double
z=Newton(A,B,0.001);
double
QL
1=(b1bx)*(z-1.0)-log(z-B)-(A(2.0*sqrt(2.0)*B))*((2
.0ax)*(x1*a1+x2*a12)-b1bx)*log(
(z+2.414*B)(z-0
.414*B));
QL1=exp(QL1);
double
QL2=(b2bx)*(z-1.0)-log(z-B)-(A(2.0*sqrt(2.0)*B))*(
(2.0ax)*(x1*a12+x2*a2)-b2bx)*log(
(z+2.414*B)(z
-0.414*B));
QL2=exp(QL2);
N2:
double ay=2.0*y1*y2*a12+y1*y1*a1+y2*y2*a2;
double by=b1*y1+b2*y2;
A=(ay*p)(R*R*T*T);
B=(by*p)(R*T);
z=Newton(A,B,1.1);
double
QV1=(b1by)*(z-1.0)-log(z-B)-(A(2.0*sqrt(2.0)*
B))*((2.0ay)*(x1*a1+x2*a12)-b1by)*log(
(z+2.414
*B)(z-0.414*B));
QV1=exp(QV1);
double
QV2=(b2by)*(z-1.0)-log(z-B)-(A(2.0*sqr
t(2.0)*B))*((2.0ay)*(x1*a12+x2*a2)-b2by)*log(
(
z+2.414*B)(z-0.414*B));
QV2=exp(QV2);
double K1=QL1QV1;
double
K2=QL2QV2;
if(k==1)
{
double y1=K1*x1(K1*x1+K2*x2);
double
y2=K2*x2(K1*x1+K2*x2);
ysum=K1*x1+K2*x2;
k++;
goto N2;
}
if(fabs(ysum-K1*x1-K2*x2)>1e-4)
22
{
ysum=K1*x1+K2*x2;
y1=K1*x1ysum;
y2=K2*x2ysum;
goto
N2;
}
else
if(fabs((K1*x1+K2*x2)-1.0)>1e-4)
{
T=T+0.001;
goto N1;
}
printf( %3.4f %3.4fn
}
cout<<
}
p=1atm
x(R32) y(R32) T
0.0000 0.0000 224.7630
0.0100 0.0155
224.6400
0.0200 0.0308 224.5200
0.0300
0.0458 224.4020
（x每隔0.01计算一次，数据较多，此处省略）
0.9500 0.9347 221.3590
0.9600 0.9469
221.4210
0.9700 0.9595 221.4880
0.9800
0.9725 221.5590
0.9900 0.9860 221.6360
1.0000 1.0000 221.7180
*****************
*************************************************
p=12atm
x(R32) y(R32) T
0.0000
0.0000 293.3800
0.0100 0.0132 293.1980
0.0200 0.0262 293.0190
0.0300 0.0392
292.8420
0.0400 0.0520 292.6680
0.0500
0.0647 292.4960
0.0600 0.0773 292.3270
（x每隔0.01计算一次，数据较多，此处省略）
0.9600 0.9577
285.9870
0.9700 0.9680 286.0130
0.9800
0.9785 286.0410
23
0.9900
0.9892 286.0730
1.000 1.0000 286.1080
Press any key to continue
使用matlab拟和图像为
24
25