matlab工程热力学大程序,西安交通大学工程期末编程大作业(完整版

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

高等工程热力学作业姓名：XX班级：XXXX学号：XXXXXXX第一章1用PR方程计算制冷剂R32，R125，和混合制冷剂R410a(R32/R125：50/50 Wt%)的pvT性质。程序说明：进入程序后选择所要计算的制冷剂，输入p，T后可得其比体积(两相区时分别输出气液相比体积)源程序：#include iostream.h#include math.h#define R 8.31451double Newton(double A,double B,double x)double x0;double f,df;dox0=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);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=;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;void R125(double T,double p,double *a,double *b,double *M)double Tc,pc,w,k,a1,Tr;*M=120.03e-3;Tc=339.45;pc=;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;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+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;void main()double M,T,a,b,p,A,B;int i;N1:coutplease enter 1(R32),2(R125)or3(R410a)endl;cini;if(i!=1&i!=2&i!=3)coutThe number is wrongendl;goto N1;coutplease enter T(K)endl;cinT;coutplease enter p(Mpa)endl;cinp;p=p*1e6;if(i=1)R32(T,p,&a,&b,&M);else if(i=2)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*T/p/M;cout单位比体积为：v1m3/kgendl;elsedouble v1=z1*R*T/p/M;double v2=z2*R*T/p/M;cout气体单位比体积为：v1m3/kgendl;cout液体单位比体积为：v2m3/kgendl;(一)please enter 1(R32),2(R125)or3(R410a)1please enter T(K)300please enter p(Mpa)1.7499气体单位比体积为：0.m3/kg液体单位比体积为：0.m3/kgPress any key to continue(二)please enter 1(R32),2(R125)or3(R410a)2please enter T(K)300please enter p(Mpa)1.7499气体单位比体积为：0.m3/kg液体单位比体积为：0.m3/kgPress any key to continue(三)please enter 1(R32),2(R125)or3(R410a)3please enter T(K)300please enter p(Mpa)1.7499气体单位比体积为：0.m3/kg液体单位比体积为：0.m3/kgPress any key to continue第二章1. 利用热力学普遍关系式推导：证明： 由理想气体状态方程得:， 代入可得： 根据热力状态基本表达式得： ， 代入得：利用麦克斯韦关系式:得：带入倒数关系， 由麦克斯韦关系：得2. 利用热力学普遍关系式推导第三dh和ds方程：解：若状态方程以p，v为独立变量 (1)比焓的变化为:式中:代回得：(2)比熵的变化代回得：3. 推导PR方程的导出热力性质余函数、。解：PR方程：上式中：所以：其中：。4用PR方程计算工质R32，R125，和混合工R32/R125的导出热力性质焓和熵。(一)计算R32，R125的焓熵值源程序#includeiostream.h#includemath.h#define R 8.31double 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(T/Tc)*(1+k*(1-sqrt(T/Tc);double a=0.45724*ar*R*R*Tc*Tc/pc;return a;double get_b(double Tc,double pc)double b=0.0778*R*Tc/pc;return b;double Newton(double A,double B,double x)double x0;double f,df;dox0=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);return x;double get_ar(double T,double v,double vv,double a,double b)double ar=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);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(v/vv);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*T/p;double vv=R*T/p;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;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;doublepM2=52.024,120.03;doublepTc2=351.255,339.45;doubleppc2=,;doublepT02=273.15,273.15;doublepp02=,;doublepw2=0.277,0.299;doublepc02=4.,2.;doublepc12=-2.,11.;doublepc22=5.,-1.;doublepc32=-1.,-0.;N1:coutplease enter 1(R32),2(R125)endl;cini;if(i=1|i=2)M=pMi-1;Tc=pTci-1;pc=ppci-1;T0=pT0i-1;p0=pp0i-1;w=pwi-1;c0=pc0i-1;c1=pc1i-1;c2=pc2i-1;c3=pc3i-1;elsecoutThe number is wrongendl;goto N1;coutplease enter T(K)endl;cinT;coutplease enter p(Mpa)endl;cinp;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)+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);double s=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;couth=hkJ/kgendl;couts=skJ/(kg*K)endl;elsedouble h=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);double s=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气相endl;couth=hkJ/kgendl;couts=skJ/(kg*K)endl;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液相endl;couth=hkJ/kgendl;couts=skJ/(kg*K)endl;please enter 1(R32),2(R125)1please enter T(K)300please enter p(Mpa)1.7749气相h=533.901kJ/kgs=2.11853kJ/(kg*K)液相h=254.453kJ/kgs=1.18608kJ/(kg*K)Press any key to continue(二)计算混合工质的焓熵值源程序#includeiostream.h#includemath.h#define R 8.31451#define k12 0.01double 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(T/Tc1)*(1+k*(1-sqrt(T/Tc1);double a1=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);double a2=0.45724*ar*R*R*Tc2*Tc2/pc2;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*Tc1/pc1;double b2=0.0778*R*Tc2/pc2;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);double bb=(a1-a2)/0.5;return bb;double Newton(double A,double B,double x)double x0;double f,df;dox0=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);return x;double get_ar(double T,double v,double vv,double a,double b)double ar=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);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(v/vv);return sr;void main()double M1=52.024;double Tc1=351.255;double pc1=;double Ts1=273.15;double ps1=;double w1=0.277;double x1;double c01=4.;double c11=-2.;double c21=5.;double c31=-1.;double M2=120.03;double Tc2=339.45;double pc2=;double Ts2=273.15;double ps2=;double w2=0.299;double x2;doublec02=2.;doublec12=11.;doublec22=-1.;doublec32=-0.;double T0=273.15;double p0=;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*T0/p0;double vv=R*T0/p0;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);coutplease enter T(K)endl;cinT;coutplease enter p(Mpa)endl;cinp;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,bb);hrl=ar+T*srl+R*T*(1-z);double dh1=R*(c01*(T-T0)+c11/2/Tc1*(pow(T,2)-pow(T0,2)+c21/3/pow(Tc1,2)*(pow(T,3)-pow(T0,3)+c31/4/pow(Tc1,3)*(pow(T,4)-pow(T0,4);double dh2=R*(c02*(T-T0)+c12/2/Tc2*(pow(T,2)-pow(T0,2)+c22/3/pow(Tc2,2)*(pow(T,3)-pow(T0,3)+c32/4/pow(Tc2,3)*(pow(T,4)-pow(T0,4);double dh=dh1*x1+dh2*x2;double ds1=-R*log(p/p0)+R*(c01*log(T/T0)+c11*(T-T0)/Tc1+c21/2/pow(Tc1,2)*(pow(T,2)-pow(T0,2)+c31/3/pow(Tc1,3)*(pow(T,3)-pow(T0,3);double ds2=-R*log(p/p0)+R*(c02*log(T/T0)+c12*(T-T0)/Tc2+c22/2/pow(Tc2,2)*(pow(T,2)-pow(T0,2)+c32/3/pow(Tc2,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;couth=hkJ/kgendl;couts=skJ/(kg*K)endl;elseh=200+(hr0-hrv+dh)/M;s=1.0+(sr0-srv+ds)/M;cout气相endl;couth=hkJ/kgendl;couts=skJ/(kg*K)endl;h=200+(hr0-hrl+dh)/M;s=1.0+(sr0-srl+ds)/M;cout液相endl;couth=hkJ/kgendl;couts=skJ/(kg*K)endl;please enter T(K)300please enter p(Mpa)1.7406气相h=246.288kJ/kgs=1.1582kJ/(kg*K)液相h=433.398kJ/kgs=1.78313kJ/(kg*K)Press any key to continue第三章1. 溶液的化学势、逸度与逸度系数定义是什么？物理意义是什么？(1)化学势的定义是：溶液系统在某些特定的条件下，作为该系统的特征函数的广延量参数对该组成摩尔数的变化率。i组成化学势的定义式为：物理意义为：化学势差是传递质量的驱动力，是由于i组成摩尔数改变而引起的相应特征广延性质改变的势位。(2)溶液逸度及逸度系数：对于变成分溶液系统，i组成的逸度及逸度系数由下式定义：上面两式中： 表示溶液中i组成的逸度； 表示溶液中i组成的逸度系数。物理意义为：是溶液中组成的热力性质，并且是强度性质，逸度和压力一样，表征了物质的逃逸势，它表示组成的假想分压力。中为组成的分压力，随着实际气体接近理想气体，在数值上接近于，所以是度量气体非理想性的标尺之一。是无量纲参数。3结合专业选用适当的方程并推导其i组成逸度的数学表达式。解：选用PR方程由： 得：而，对上式右边加减，得当为常数时，上式对求导，则又因为得左侧可以表示成：代回：又其中则：带回上式进行积分得：式中第四章1用Peng-Robinson方程计算纯质R32的p-T图和溶液R32/R125分别在p=1atm和p=12atm下的T-x图。(一)R32的p-T图源程序：#includeiostream.h#includemath.hdouble Newton(double A,double B,double x)double x0;double f,df;dox0=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);return x;void main()double Tc=351.255;double pc=;double w=0.277;double R=8.31451;double T,p=1e6,a,b,A,B,faiv,fail,z;coutT(K) p(Mpa)endl;for(int i=0;i=10;i+)T=270;dodouble k=0.37464+1.54226*w-0.26992*w*w;double ar=(1+k*(1-sqrt(T/Tc)*(1+k*(1-sqrt(T/Tc);a=0.45724*ar*R*R*Tc*Tc/pc;b=0.0778*R*Tc/pc;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);coutT p/1e6endl;p=p+0.1e6;T(K) p(Mpa)279.619 1282.759 1.1285.691 1.2288.443 1.3291.041 1.4293.503 1.5295.844 1.6298.078 1.7300.215 1.8302.265 1.9304.236 2Press any key to continue使用matlab拟和图像为(二)溶液R32/R125分别在p=1atm和p=12atm下的T-x图源程序：#includeiostream.h#includemath.h#includestdio.h#define R 8.31451double Newton(double A,double B,double x)double x0;double f,df;dox0=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);return x;void main( )int k=1;double ysum;double M1=52.024;double Tc1=351.255;double pc1=;double w1=0.277;double M2=120.03;double Tc2=339.45;double pc2=;double w2=0.299;double k12=0.01;double p;double P2=1.013e5,12*1.013e5;double TT2=200,260;for(int j=0;j=1;j+)p=Pj;coutp=p/1.013e5atmendl;coutx(R32) y(R32) Tendl;for(int i=0;i=100;i+)double x1=(double)i/100,x2=1-x1;double y1=0.5,y2=0.5;double T=TTj;N1:double k1=0.37464+1.54226*w1-0.26992*w1*w1;double Tr1=T/Tc1;double e1=(1.0+k1*(1-sqrt(Tr1)*(1.0+k1*(1.0-sqrt(Tr1);double a1=0.45724*e1*R*R*Tc1*Tc1/pc1;double b1=0.0778*R*Tc1/pc1;double k2=0.37464+1.54226*w2-0.26992*w2*w2;double Tr2=T/Tc2;double e2=(1.0+k2*(1-sqrt(Tr2)*(1.0+k2*(1.0-sqrt(Tr2);double a2=0.45724*e2*R*R*Tc2*Tc2/pc2;double b2=0.0778*R*Tc2/pc2;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;double A=(ax*p)/(R*R*T*T);double B=(bx*p)/(R*T);double z=Newton(A,B,0.001);double QL1=(b1/bx)*(z-1.0)-log(z-B)-(A/(2.0*sqrt(2.0)*B)*(2.0/ax)*(x1*a1+x2*a12)-b1/bx)*log(z+2.414*B)/(z-0.414*B);QL1=exp(QL1);double QL2=(b2/bx)*(z-1.0)-log(z-B)-(A/(2.0*sqrt(2.0)*B)*(2.0/ax)*(x1*a12+x2*a2)-b2/bx)*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=(b1/by)*(z-1.0)-log(z-B)-(A/(2.0*sqrt(2.0)*B)*(2.0/ay)*(x1*a1+x2*a12)-b1/by)*log(z+2.414*B)/(z-0.414*B);QV1=exp(QV1);double QV2=(b2/by)*(z-1.0)-log(z-B)-(A/(2.0*sqrt(2.0)*B)*(2.0/ay)*(x1*a12+x2*a2)-b2/by)*log(z+2.414*B)/(z-0.414*B);QV2=exp(QV2);double K1=QL1/QV1;double K2=QL2/QV2;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)ysum=K1*x1+K2*x2;y1=K1*x1/ysum;y2=K2*x2/ysum;goto N2;else if(fabs(K1*x1+K2*x2)-1.0)1e-4)T=T+0.001;goto N1;printf(%3.4f %3.4f %3.4fn,x1,y1,T);cout*endl;p=1atmx(R32) y(R32) T0.0000 0.0000 224.76300.0100 0.0155 224.64000.0200 0.0308 224.52000.0300 0.0458 224.4020(x每隔0.01计算一次，数据较多，此处省略)0.9500 0.9347 221.35900.9600 0.9469 221.42100.9700 0.9595 221.48800.9800 0.9725 221.55900.9900 0.9860 221.63601.0000 1.0000 221.7180*p=12atmx(R32) y(R32) T0.0000 0.0000 293.38000.0100 0.0132 293.19800.0200 0.0262 293.01900.0300 0.0392 292.84200.0400 0.0520 292.66800.0500 0.0647 292.49600.0600 0.0773 292.3270(x每隔0.01计算一次，数据较多，此处省略)0.9600 0.9577 285.98700.9700 0.9680 286.01300.9800 0.9785 286.04100.9900 0.9892 286.07301.000 1.0000 286.1080 Press any key to continue使用matlab拟和图像为