bwrs.cpp

bwrs方程计算天然气水合物的物性参数值

#include "math.h"
#include "iostream.h"
#include "fstream.h"

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

void main()
{
	ofstream fout("bwrsresult.txt");
	//通用常数
	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],u[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;	ro[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=0;T=273;

/*/////////////////////   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;
	fout.close();
}