#include "math.h"
#include "iostream"
#include "fstream"
#include<iostream>
using namespace std;
double F(double row,double R,double P,double T,double A0,double B0,double C0,
double D0,double E0,double a,double b,double c,double d,double alfa,
double gama)
{
double aaa=row*R*T+(B0*R*T-A0-C0/T/T+D0/pow(T,3.0)-E0/pow(T,4.0)*row*row
+(b*R*T-a-d/T)*row*row*row+alfa*(a+d/T)*pow(row,6.0)+c*pow(row,3.0)/T/T*exp(-gama*row*row)
-P);
return aaa;
}
double Density(double R,double P,double T,double A0,double B0,double C0,
double D0,double E0,double a,double b,double c,double d,double alfa,
double gama)
{
double ro1,ro2,ro0;
ro1=0.0,ro2=P/R/T;
while(fabs(ro1-ro2)>pow(10.0,-4.0))
{
ro0=ro2;
ro2=(ro1*F(ro2,R,P,T,A0,B0,C0,D0,E0,a,b,c,d,alfa,gama)
-ro2*F(ro1,R,P,T,A0,B0,C0,D0,E0,a,b,c,d,alfa,gama))
/(F(ro2,R,P,T,A0,B0,C0,D0,E0,a,b,c,d,alfa,gama)
-F(ro1,R,P,T,A0,B0,C0,D0,E0,a,b,c,d,alfa,gama));
ro1=ro0;
}
return ro2;
}
int main()
{
//通用常数
double A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11;
//临界温度,临界压力,临界密度,偏心因子,分子量
double Tc[9],Pc[9],roc[9],w[9];
//单组分
double Aa0[9],Bb0[9],Cc0[9],Dd0[9],Ee0[9],a0[9],b0[9],c0[9],d0[9],alfa0[9],gama0[9];
//混合物
double A0,B0,C0,D0,E0,a,b,c,d,alfa,gama;
//二元交互作用系数
double K[9][9];
//气体的摩尔分数
double y[9];
//PVT 方程中的其它变量
double P,T,R,Z;
int i,j,N;//N,组分数
double density;//密度
double Mw;//混合气体分子量
/*///////////////////////// 1. 输入数据 //////////////////////*/
N=9;
R=8.3143;//kJ/(kmol.K)
//通用常数默认
A1=0.443690; A2=1.284380; A3=0.356306; A4=0.544979; A5=0.528629;
A6=0.484011; A7=0.0705233; A8=0.504087; A9=0.0307452; A10=0.0732828;
A11=0.006450;
B1=0.115449; B2=-0.920731; B3=1.70871; B4=-0.270896; B5=0.349621;
B6=0.754130; B7=-0.044448; B8=1.32245; B9=0.179433; B10=0.463492;
B11=-0.022143;
//临界温度,临界压力,临界密度,偏心因子,可输入
Tc[1]=190.69; Tc[2]=305.38; Tc[3]=369.89; Tc[4]=425.18; Tc[5]=408.13;
Tc[6]=469.49; Tc[7]=460.37; Tc[8]=126.15; Tc[9]=304.09;
Pc[1]=4.604; Pc[2]=4.880; Pc[3]=4.250; Pc[4]=3.797; Pc[5] =3.648;
Pc[6]=3.369; Pc[7]=3.374; Pc[8]=3.394; Pc[9]=7.376;
roc[1]=10.050; roc[2]=6.7566; roc[3]=4.9994; roc[4]=3.9213; roc[5]=3.8012;
roc[6]=3.2149; roc[7]=3.2469; roc[8]=11.099; roc[9]=10.638;
w[1]=0.013; w[2]=0.1018; w[3]=0.157; w[4]=0.197; w[5]=0.183;
w[6]=0.252; w[7]=0.226; w[8]=0.035; w[9]=0.210;
//二元交互作用系数
K[1][1]=0.;
K[2][1]=0.010; K[2][2]=0.0;
K[3][1]=0.023; K[3][2]=0.0031; K[3][3]=0.0;
K[4][1]=0.031; K[4][2]=0.0045; K[4][3]=0.0035; K[4][4]=0.0;
K[5][1]=0.0275; K[5][2]=0.004; K[5][3]=0.003; K[5][4]=0.0; K[5][5]=0.0;
K[6][1]=0.041; K[6][2]=0.006; K[6][3]=0.0045; K[6][4]=0.001; K[6][5]=0.001; K[6][6]=0.0;
K[7][1]=0.036; K[7][2]=0.005; K[7][3]=0.004; K[7][4]=0.008; K[7][5]=0.008; K[7][6]=0.0; K[7][7]=0.0;
K[8][1]=0.025; K[8][2]=0.007; K[8][3]=0.10; K[8][4]=0.12; K[8][5]=0.11; K[8][6]=0.148; K[8][7]=0.134; K[8][8]=0.0;
K[9][1]=0.05; K[9][2]=0.048; K[9][3]=0.045; K[9][4]=0.05; K[9][5]=0.05; K[9][6]=0.05; K[9][7]=0.05; K[9][8]=0.0; K[9][9]=0.0;
for(i=1;i<=N;i++)
{
for(j=i+1;j<=N;j++)
K[i][j]=K[j][i];
}
//各组分分子量
static double u[9]={16.042,30.068,44.094,58.120,58.120,72.146,72.146,28.016,44.010};
//气体的摩尔分数,可输入
y[1]=0.866; y[2]=0.042; y[3]=0.035; y[4]=0.019;y[5]=0.007; y[6]=0.005; y[7]=0.006; y[8]=0.011;
y[9]=0.009;
//输入工况
// cout<<"请输入要计算的工况(压力P--MPa,温度T--K):"<<endl;
P=4.0;T=293;
/*///////////////////// 2. 单组分系数 ///////////////////////////////////////*/
//平均分子量
Mw=0;//混合气体分子量
for(i=1;i<=N;i++)
Mw+=y[i]*u[i];
for(i=1;i<=N;i++)
{
Bb0[i]=(A1+B1*w[i])/roc[i];
Aa0[i]=(A2+B2*w[i])/roc[i]*R*Tc[i];
Cc0[i]=(A3+B3*w[i])/roc[i]*R*pow(Tc[i],3.0);
gama0[i]=(A4+B4*w[i])/roc[i]/roc[i];
b0[i]=(A5+B5*w[i])/roc[i]/roc[i];
a0[i]=(A6+B6*w[i])/roc[i]/roc[i]*R*Tc[i];
alfa0[i]=(A7+B7*w[i])/pow(roc[i],3.0);
c0[i]=(A8+B8*w[i])/roc[i]/roc[i]*R*pow(Tc[i],3.0);
Dd0[i]=(A9+B9*w[i])/roc[i]*R*pow(Tc[i],4.0);
d0[i]=(A10+B10*w[i])/roc[i]/roc[i]*R*Tc[i]*Tc[i];
Ee0[i]=(A11+B11*w[i]*exp(-3.8*w[i]))*R*pow(Tc[i],5.0)/roc[i];
}
/*///////////////////// 3. 混合气体系数 ///////////////////////////////////////*/
A0=B0=C0=D0=E0=a=b=c=d=alfa=gama=0;
for(i=1;i<=N;i++)
{
B0+=y[i]*Bb0[i];
b+=y[i]*pow(b0[i],1.0/3.0);
a+=y[i]*pow(a0[i],1.0/3.0);
c+=y[i]*pow(c0[i],1.0/3.0);
d+=y[i]*pow(d0[i],1.0/3.0);
alfa+=y[i]*pow(alfa0[i],1.0/3.0);
gama+=y[i]*pow(gama0[i],0.5);
for(j=1;j<=N;j++)
{
A0+=y[i]*y[j]*sqrt(Aa0[i]*Aa0[j])*(1.0-K[i][j]);
C0+=y[i]*y[j]*sqrt(Cc0[i]*Cc0[j])*pow(1.0-K[i][j],3.0);
D0+=y[i]*y[j]*sqrt(Dd0[i]*Dd0[j])*pow(1.0-K[i][j],4.0);
E0+=y[i]*y[j]*sqrt(Ee0[i]*Ee0[j])*pow(1.0-K[i][j],5.0);
}
}
a=pow(a,3.0);
b=pow(b,3.0);
c=pow(c,3.0);
d=pow(d,3.0);
alfa=pow(alfa,3.0);
gama=pow(gama,2.0);
/*///////////////////// 5. 求流体密度和压缩系数Z ///////////////////////////////////////*/
density=Density(R,P,T,A0,B0,C0,D0,E0,a,b,c,d,alfa,gama);
Z=P/(density*R*T);
/*///////////////////// 6.输出计算结果 /////////////////////////////////////*/
cout<<"\n++++++++++++++++++++++++++计算结果++++++++++++++++++++++++"<<endl;
cout<<" 密度p="<<density<<" kmol/m^3,\t\t"<<"压缩系数Z="<<Z<<""<<endl<<endl;
}
- 1
- 2
前往页