西安交通大学工程期末编程大作业完整版.docx

西安交通大学工程期末编程大作业完整版.docx

《西安交通大学工程期末编程大作业完整版.docx》由会员分享，可在线阅读，更多相关《西安交通大学工程期末编程大作业完整版.docx（33页珍藏版）》请在冰豆网上搜索。

西安交通大学工程期末编程大作业完整版

高等工程热力学作业

姓名：

XX

班级：

XXXX

学号：

XXXXXXX

第一章

1．用PR方程计算制冷剂R32，R125，和混合制冷剂R410a（R32/R125：

50/50Wt%）的pvT性质。

程序说明：

进入程序后选择所要计算的制冷剂，输入p，T后可得其比体积（两相区时分别输出气液相比体积）

源程序：

#include"iostream.h"

#include"math.h"

#defineR8.31451

doubleNewton（doubleA,doubleB,doublex）

{

doublex0;

doublef,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-f/df;

}while（fabs（x-x0）>1e-6）;

returnx;

}

voidR32（doubleT,doublep,double*a,double*b,double*M）

{

doubleTc,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=T/Tc;

a1=pow（1+k*（1-pow（Tr,0.5））,2）;

*a=0.45727*a1*R*R*Tc*Tc/pc;

*b=0.07780*R*Tc/pc;

}

voidR125（doubleT,doublep,double*a,double*b,double*M）

{

doubleTc,pc,w,k,a1,Tr;

*M=120.03e-3;

Tc=339.45;

pc=3630600;

w=0.299;

k=0.37464+1.54226*w-0.26992*w*w;

Tr=T/Tc;

a1=pow（1+k*（1-pow（Tr,0.5））,2）;

*a=0.45727*a1*R*R*Tc*Tc/pc;

*b=0.07780*R*Tc/pc;

}

voidR410a（doubleT,doublep,double*a,double*b,double*M）

{

doublea1,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+M1/M2）;

x2=1/（1+M2/M1）;

*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;

}

voidmain（）

{

doubleM,T,a,b,p,A,B;

inti;

N1:

cout<<"pleaseenter1（R32）,2（R125）or3（R410a）"<

cin>>i;

if（i!

=1&&i!

=2&&i!

=3）

{

cout<<"Thenumberiswrong"<

gotoN1;

}

cout<<"pleaseenterT（K）"<

cin>>T;

cout<<"pleaseenterp（Mpa）"<

cin>>p;

p=p*1e6;

if（i==1）

{

R32（T,p,&a,&b,&M）;

}

elseif（i==2）

{

R125（T,p,&a,&b,&M）;

}

elseif（i==3）

{

R410a（T,p,&a,&b,&M）;

}

A=a*p/（R*R*T*T）;

B=b*p/（R*T）;

doublez1=Newton（A,B,1000）;

doublez2=Newton（A,B,0.001）;

if（fabs（z1-z2）<1e-4）

{

doublev1=z1*R*T/p/M;

cout<<"单位比体积为：

"<

}

else

{

doublev1=z1*R*T/p/M;

doublev2=z2*R*T/p/M;

cout<<"气体单位比体积为：

"<

cout<<"液体单位比体积为：

"<

}

}

（一）

pleaseenter1（R32）,2（R125）or3（R410a）

1

pleaseenterT（K）

300

pleaseenterp（Mpa）

1.7499

气体单位比体积为：

0.0214228m^3/kg

液体单位比体积为：

0.00122247m^3/kg

Pressanykeytocontinue

（二）

pleaseenter1（R32）,2（R125）or3（R410a）

2

pleaseenterT（K）

300

pleaseenterp（Mpa）

1.7499

气体单位比体积为：

0.00762949m^3/kg

液体单位比体积为：

0.000859649m^3/kg

Pressanykeytocontinue

（三）

pleaseenter1（R32）,2（R125）or3（R410a）

3

pleaseenterT（K）

300

pleaseenterp（Mpa）

1.7499

气体单位比体积为：

0.0147329m^3/kg

液体单位比体积为：

0.00105639m^3/kg

Pressanykeytocontinue

第二章

1.利用热力学普遍关系式推导：

证明：

由理想气体状态方程得:

，

代入可得：

根据热力状态基本表达式得：

，

代入得：

利用麦克斯韦关系式:

得：

带入倒数关系

，

由麦克斯韦关系：

得

2.利用热力学普遍关系式推导第三dh和ds方程：

解：

若状态方程以p，v为独立变量

（1）比焓的变化为:

式中:

代回得：

（2）比熵的变化

代回得：

3.推导PR方程的导出热力性质余函数

、

。

解：

PR方程：

上式中：

所以：

其中：

。

4．用PR方程计算工质R32，R125，和混合工R32/R125的导出热力性质焓和熵。

（一）计算R32，R125的焓熵值

源程序

#include

#include

#defineR8.31

doubleget_a（doublew,doubleT,doubleTc,doublepc）

{

doublek=0.37464+1.54226*w-0.26992*w*w;

doublear=（1+k*（1-sqrt（T/Tc）））*（1+k*（1-sqrt（T/Tc）））;

doublea=0.45724*ar*R*R*Tc*Tc/pc;

returna;

}

doubleget_b（doubleTc,doublepc）

{

doubleb=0.0778*R*Tc/pc;

returnb;

}

doubleNewton（doubleA,doubleB,doublex）

{

doublex0;

doublef,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-f/df;

}while（fabs（x-x0）>1e-6）;

returnx;

}

doubleget_ar（doubleT,doublev,doublevv,doublea,doubleb）

{

doublear=R*T*log（（v-b）/v）-a*log（（v-0.414*b）/（v+2.414*b））/（2*1.414*b）+R*T*log（v/vv）;

returnar;

}

doubleget_sr（doubleT,doublev,doublevv,doublea,doubleb,doublebb）

{

doublesr=-1*R*log（（v-b）/v）+bb*log（（v-0.414*b）/（v+2.414*b））/（2*1.414*b）-R*log（v/vv）;

returnsr;

}

voidget_hr（doubleTc,doublepc,doublew,doubleT,doublep,double*hr,double*sr,doublezz）

{

doublea=get_a（w,T,Tc,pc）;

doubleb=get_b（Tc,pc）;

doublebb=（get_a（w,T+0.25,Tc,pc）-get_a（w,T-0.25,Tc,pc））/0.5;

doubleA=a*p/（R*R*T*T）;

doubleB=b*p/（R*T）;

doublez=Newton（A,B,zz）;

doublev=z*R*T/p;

doublevv=R*T/p;

doublear=get_ar（T,v,vv,a,b）;

*sr=get_sr（T,v,vv,a,b,bb）;

*hr=ar+T*（*sr）+R*T*（1-z）;

}

voidmain（）

{

inti;

doubleM;

doublew;

doubleh0=200.0;

doubles0=1.0;

doubleT0;

doublep0;

doublec0,c1,c2,c3;

doubleT,p,Tc,pc;

doublehr0,sr0,hrv,hrl,srv,srl;

doublepM[2]={52.024,120.03};

doublepTc[2]={351.255,339.45};

doubleppc[2]={5780000,3630600};

doublepT0[2]={273.15,273.15};

doublepp0[2]={813100,617500};

doublepw[2]={0.277,0.299};

doublepc0[2]={4.424901,2.838680};

doublepc1[2]={-2.661170,11.581633};

doublepc2[2]={5.580232,-1.704482};

doublepc3[2]={-1.680558,-0.266732};

N1:

cout<<"pleaseenter1（R32）,2（R125）"<

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<<"Thenumberiswrong"<

gotoN1;

}

cout<<"pleaseenterT（K）"<

cin>>T;

cout<<"pleaseenterp（Mpa）"<

cin>>p;

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）

{

doubleh=h0*M+hr0-hrv+R*（c0*（T-T0）+c1/2/Tc*（pow（T,2）

-pow（T0,2））+c2/3/pow（Tc,2）*（pow（T,3）-pow（T0,3））+c3/4/pow（Tc,3）

*（pow（T,4）-pow（T0,4）））;

doubles=s0*M+sr0-srv-R*log（p/p0）+R*（c0*log（T/T0）+c1*（T-T0）/Tc

+c2/2/pow（Tc,2）*（pow（T,2）-pow（T0,2））+c3/3/pow（Tc,3）*（pow（T,3）-pow（T0,3）））;

h=h/M;

s=s/M;

cout<<"h="<

cout<<"s="<

}

else

{

doubleh=h0*M+hr0-hrv+R*（c0*（T-T0）+c1/2/Tc*（pow（T,2）

-pow（T0,2））+c2/3/pow（Tc,2）*（pow（T,3）-pow（T0,3））+c3/4/pow（Tc,3）

*（pow（T,4）-pow（T0,4）））;

doubles=s0*M+sr0-srv-R*log（p/p0）+R*（c0*log（T/T0）+c1*（T-T0）/Tc

+c2/2/pow（Tc,2）*（pow（T,2）-pow（T0,2））+c3/3/pow（Tc,3）*（pow（T,3）-pow（T0,3）））;

h=h/M;

s=s/M;

cout<<"气相"<

cout<<"h="<

cout<<"s="<

h=h0*M+hr0-hrl+R*（c0*（T-T0）+c1/2/Tc*（pow（T,2）

-pow（T0,2））+c2/3/pow（Tc,2）*（pow（T,3）-pow（T0,3））+c3/4/pow（Tc,3）

*（pow（T,4）-pow（T0,4）））;

s=s0*M+sr0-srl-R*log（p/p0）+R*（c0*log（T/T0）+c1*（T-T0）/Tc

+c2/2/pow（Tc,2）*（pow（T,2）-pow（T0,2））+c3/3/pow（Tc,3）*（pow（T,3）-pow（T0,3）））;

h=h/M;

s=s/M;

cout<<"液相"<

cout<<"h="<

cout<<"s="<

}

}

pleaseenter1（R32）,2（R125）

1

pleaseenterT（K）

300

pleaseenterp（Mpa）

1.7749

气相

h=533.901kJ/kg

s=2.11853kJ/（kg*K）

液相

h=254.453kJ/kg

s=1.18608kJ/（kg*K）

Pressanykeytocontinue

（二）计算混合工质的焓熵值

源程序

#include

#include

#defineR8.31451

#definek120.01

doubleget_a（doublew1,doublew2,doubleT,doubleTc1,doublepc1,doubleTc2,

doublepc2,doublex1,doublex2）

{

doublek=0.37464+1.54226*w1-0.26992*w1*w1;

doublear=（1+k*（1-sqrt（T/Tc1）））*（1+k*（1-sqrt（T/Tc1）））;

doublea1=0.45724*ar*R*R*Tc1*Tc1/pc1;

k=0.37464+1.54226*w2-0.26992*w2*w2;

ar=（1+k*（1-sqrt（T/Tc2）））*（1+k*（1-sqrt（T/Tc2）））;

doublea2=0.45724*ar*R*R*Tc2*Tc2/pc2;

doublea=2*x1*x2*（1-k12）*sqrt（a1*a2）+x1*x1*a1+x2*x2*a2;

returna;

}

doubleget_b（doubleTc1,doublepc1,doubleTc2,doublepc2,doublex1,doublex2）

{

doubleb1=0.0778*R*Tc1/pc1;

doubleb2=0.0778*R*Tc2/pc2;

doubleb=x1*b1+x2*b2;

returnb;

}

doubleget_bb（doublew1,doublew2,doubleT,doubleTc1,doublepc1,doubleTc2,

doublepc2,doublex1,doublex2）

{

doublea1=get_a（w1,w2,T+0.25,Tc1,pc1,Tc2,pc2,x1,x2）;

doublea2=get_a（w1,w2,T-0.25,Tc1,pc1,Tc2,pc2,x1,x2）;

doublebb=（a1-a2）/0.5;

returnbb;

}

doubleNewton（doubleA,doubleB,doublex）

{

doublex0;

doublef,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-f/df;

}while（fabs（x-x0）>1e-6）;

returnx;

}

doubleget_ar（doubleT,doublev,doublevv,doublea,doubleb）

{

doublear=R*T*log（（v-b）/v）-a*log（（v-0.414*b）/（v+2.414*b））/（2*1.414*b）+R*T*log（v/vv）;

returnar;

}

doubleget_sr（doubleT,doublev,doublevv,doublea,doubleb,doublebb）

{

doublesr=-1*R*log（（v-b）/v）+bb*log（（v-0.414*b）/（v+2.414*b））/（2*1.414*b）-R*log（v/vv）;

returnsr;

}

voidmain（）

{

doubleM1=52.024;

doubleTc1=351.255;

doublepc1=5780000;

doubleTs1=273.15;

doubleps1=813100;

doublew1=0.277;

doublex1;

doublec01=4.424901;

doublec11=-2.661170;

doublec21=5.580232;

doublec31=-1.680558;

doubleM2=120.03;

doubleTc2=339.45;

doublepc2=3630600;

doubleTs2=273.15;

doubleps2=671500;

doublew2=0.299;

doublex2;

doublec02=2.838680;

doublec12=11.581633;

doublec22=-1.704482;

doublec32=-0.266732;

doubleT0=273.15;

doublep0=799100;

doublehr0,sr0,hrv,srv,hrl,srl,h,s;

doublep,T;

x1=（M2）/（M1+M2）;

x2=1-x1;

doubleM=M1*x1+M2*x2;

doublea=get_a（w1,w2,T0,Tc1,pc1,Tc2,pc2,x1,x2）;

doubleb=get_b（Tc1,pc1,Tc2,pc2,x1,x2）;

doublebb=get_bb（w1,w2,T0,Tc1,pc1,Tc2,pc2,x1,x2）;

doubleA=a*p0/（R*R*T0*T0）;

doubleB=b*p0/（R*T0）;

doublez=Newton（A,B,0.001）;

doublev=z*R*T0/p0;

doublevv=R*T0/p0;

doublear=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<<"pleaseenterT（K）"<

cin>>T;

cout<<"pleaseenterp（Mpa）"<

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）;

v=z*R*T/p;

vv=R*T/p;

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*T/p;

vv=R*T/p;

ar=get_ar（T,v,vv,a,b）;

srl=get_sr（T,v,vv,a,b,