#include"math.h"
#include"stdio.h"
double b[9][4]; /*常数*/
double R; /*气体常数*/
double tc; /*临界温度*/
double pc; /*临界压力*/
double denc; /*临界密度*/
double Ri,Pa; /*饱和蒸汽压状态方程的准则数*/
double a1,a2; /*饱和液体密度方程的常数*/
double d[5]; /*定容比热容方程的常数*/
double hc,sc; /*根据基准态的焓值和熵值确定的常数*/
/*理想气体状态下的比内能u0*/
double u0(double tr)
{
double max,max1,max2;
int i;
max1=0.0;
for(i=0;i<=4;i++)
{
max2=d[i]*pow(tr,(i+1))/(i+1);
max1=max1+max2;
}
max=tc*max1;
return(max);
}
/*理想气体状态下的定容比熵s0*/
double s0(double tr)
{
double max,max1,max2;
int i;
max1=0.0;
for(i=1;i<=4;i++)
{
max2=d[i]*pow(tr,i)/i;
max1=max1+max2;
}
max=d[0]*log(tr)+max1;
return(max);
}
/*中间变量l*/
double l(int q,double tr,int n,double den)
{
int i,j;
double l1,l2,r;
l1=0.0;
l2=0.0;
for(i=1;i<=n;i++)
{
for(j=1;j<=3;j++)
{
{
if(q==1)
r=1;
else if(q==2)
r=double(1/i);
else
r=double(j/i);
}
l1=l1+r*b[i][j]*pow(den,i)*pow(tr,i);
}
l2=l2+l1;
}
return(l2);
}
/*气相区的状态方程*/
double z(double den,double tr,int n)
{
int i,j;
double z1,z2,z3;
z2=0.0;
for(i=1;i<=n;i++)
{
z1=0.0;
for(j=0;j<=3;j++)
{
z1=z1+b[i][j]/pow(tr,j);
}
z2=z2+z1*pow(den,i);
}
z3=1.0+z2;
return(z3);
}
/*饱和蒸汽压方程的str*/
double str(double tr)
{
double max;
max=(tr-1.0)*(0.2*pow((tr+1.0),2.0)+0.5);
return(max);
}
double bqp,bqden,bqh,bqs;
void main()
{
double t,T,tr,bq1,bqden1,z0;
int n=8;
printf("输入温度t:\n");
scanf("%lf",&t);
b[1][0]=4.53465,b[1][1]=-14.6817,b[1][2]=13.5886,b[1][3]=-9.91264;
b[2][0]=-12.0457647,b[2][1]=74.4227631,b[2][2]=-95.9087499,b[2][3]=50.8152618;
b[3][0]=118.985614,b[3][1]=-629.351137,b[3][2]=557.606974,b[3][3]=70.9903655;
b[4][0]=-109.451026,b[4][1]=2106.02214,b[4][2]=-1785.21688,b[4][3]=0;
b[5][0]=-1095.97522,b[5][1]=-2799.47726,b[5][2]=2461.60212,b[5][3]=0;
b[6][0]=3436.96088,b[6][1]=1525.29115,b[6][2]=-943.536139,b[6][3]=0;
b[7][0]=-4151.8895,b[7][1]=-839.346672,b[7][2]=0,b[7][3]=0;
b[8][0]=2350.2697,b[8][1]=0,b[8][2]=0,b[8][3]=0;
R=0.48816;
tc=405.55;
pc=113.97;
denc=0.2350;
Ri=7.0284,Pa=-0.3958;
a1=1.6389,a2=0.3859;
d[0]=1.7262,d[1]=-1.4470,d[2]=2.5387,d[3]=-1.2409,d[4]=0.21857;
hc=1144.764,sc=7.1147;
T=273.15+t;
tr=T/tc;
printf("%lf\n",tr);
bq1=4*(tr-1.0)/tr+str(tr)-5.3*log(tr);
printf("%lf\n",bq1);
bqp=pc*exp(Ri*log(tr)+(Ri-4.0+Pa)*bq1);
printf("饱和蒸汽压:%lf\n",bqp);
/*利用迭代法求出密度*/
bqden1=bqp/(10.0*R*T);
do
{
bqden=bqden1;
printf("密度:%lf\n",bqden);
z0=z(bqden,tr,n);
printf("z:%lf\n",z0);
bqden1=bqp/(10*R*T*z0);
}while(fabs(bqden1-bqden)>1e-6);
/*求出焓值h和熵值s*/
bqh=R*tc*tr*(1+l(1,tr,n,bqden)+l(3,tr,n,bqden))+u0(tr)+hc;
bqs=R*(l(3,tr,n,bqden)-l(2,tr,n,bqden)-log(bqden))+s0(tr)+sc;
}
