|
本帖最后由 wangchuanlei 于 2012-4-27 11:37 编辑
#define N 14
#define O 0.45742
#define R 8.31451
#include<stdio.h>
#include<math.h>
void main()
{
double Tc[N]={190.56,305.32,369.83,425.12,407.8,469.7,460.4,433.8,507.6,497.7,504.6,489.0,500.0,540.2};
double Pc[N]={4599e3,4872e3,4248e3,3796e3,3640e3,3370e3,3380e3,3196e3,3025e3,3040e3,3120e3,3100e3,3150e3,2740e3};
double w[N]={0.011,0.099,0.152,0.199,0.177,0.249,0.228,0.196,0.305,0.278,0.274,0.234,0.248,0.351};
int i,j;
float p,T;
double aij,A,B,b=0,a=0,b1,c1,d1,p1,q1,r,y0,y1,y2,z0,z1,z2,x1,m,n;
/* printf("按顺序输入天然气摩尔组成:\n");
for(i=0;i<N;i++)
scanf("%f",&x);*/
double x[N]={1.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00};
printf("\n绝对压力/Mpa:");
scanf("%f",&p);
printf("\n绝对温度/K:");
scanf("%f",&T);
printf("\n");
/* K1=1-1.41263+0.02461*T-1.40764e-4*T*T+2.55299e-7*T*T*T;
printf("结果:K1=%.4f\n",K1);*/
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
aij=pow(pow(O,2)*pow(R*Tc,2)*pow(R*Tc[j],2)*pow(1+(0.37464+1.54226*w-0.26992*w*w)*(1-pow(T/Tc,0.5)),2)*pow(1+(0.37464+1.54226*w[j]-0.26992*w[j]*w[j])*(1-pow(T/Tc[j],0.5)),2)/(Pc*Pc[j]),0.5);
a+=x*x[j]*aij;
printf("结果:aij=%.4f\n",aij);
}
}
printf("结果:a=%.4f\n",a);
for(i=0;i<N;i++)
{
b+=x*0.0778*R*Tc/Pc;
}
printf("结果:b=%.4f\n",b);
B=b*p*1e6/(R*T);
A=a*p*1e6/pow(R*T,2);
printf("结果:B=%.4f\n",B);
printf("结果:A=%.4f\n",A);
b1=B-1;
c1=A-2*B-3*B*B;
d1=B*B*B+B*B-A*B;
p1=c1-b1*b1/3;
q1=d1-b1*c1/3+2*b1*b1*b1/27;
r=q1*q1/4+p1*p1*p1/27;
printf("结果:q1=%.4f\n",q1);
printf("结果:r=%.6f\n",r);
printf("结果:d1=%.4f\n",d1);
printf("结果:p1=%.4f\n",p1);
if(r>0)
{
m=-q1/2+pow(r,0.5);
n=-q1/2-pow(r,0.5);
printf("结果:m=%.4f\n",m);
printf("结果:n=%.4f\n",n);
if(n>=0)
{
y0=pow(m,1.0/3)+pow(n,1.0/3);
}
else
{
y0=pow(m,1.0/3)-pow(-n,1.0/3);
}
printf("结果:y0=%.6f\n",y0);
z0=y0-b1/3;
printf("结果:y0=%.4f\n",y0);
printf("结果:Z0=%.4f\n",z0);
}
else if(r=0)
{
y0=2*pow(-q1/2,1.0/3.0);
y1=pow(q1/2,1.0/3.0);
z0=y0-b1/3;
z1=y1-b1/3;
printf("结果:Z0=%.4f\n",z0);
printf("结果:Z1=%.4f\n",z1);
}
else
{
m=-q1/(2*pow(-p1*p1*p1/27,0.5));
x1=acos(m)/3;
y0=2*pow(-p1/3,0.5)*cos(x1);
y1=2*pow(-p1/3,0.5)*cos(x1+2*3.14159/3);
y2=2*pow(-p1/3,0.5)*cos(x1+4*3.14159/3);
z0=y0-b1/3;
z1=y1-b1/3;
z2=y2-b1/3;
printf("结果:m=%.6f\n",m);
printf("结果:x1=%.6f\n",x1);
printf("结果:Z0=%.4f\n",z0);
printf("结果:Z1=%.4f\n",z1);
printf("结果:Z2=%.4f\n",z2);
}
}
|
|