PR方程计算相平衡

#include 
#include 
#include
#include 

double secantmethod（doubledoubledoubledouble）;	///正割法求解
double pr（doubledoubledoubledouble）;			///关于压缩因子的PR方程
double newton（doubledoubledoubledouble）;		///牛顿法求解
double pr1（doubledoubledoubledouble）;	///求PR方程的导数值

void main（）
{
	double hydrocarbon[5][3]={{190.446.00.01155}{305.448.80.099}{369.842.50.153}{408.236.50.183}{425.238.00.199}};
									 ///甲烷methane、乙烷ethane、丙烷propane、异丁烷isobutane、正戊烷nbutane
									///临界温度K、临界压力bar、偏心因子omiga
	double component[5]={0.8420.1140.0280.0090.007};///天然气的组成

	double acom[5]={0};			///PR方程中纯组分的系数a
	double bcom[5]={0};			///PR方程中纯组分的系数b
	double omigaa=0.45724;		///计算PR方程中系数a和b的中间参数
	double omigab=0.07780;
	double r=8.314;
	double tr=0;
	double ac=0;
	double alf=0;
	double q=0;
	double a=0;
	double b=0;
	double k=0;
	double z=0;
	double aij[5][5]={0};		///混合规则中计算a
	double p=0;
	double t=0;
	double a1=0;
	double b1=0;
	double a2=0;
	double b2=0;
	double c2=0;
	double moldensity=0;
	double tc=0;
	int i=0;
	int j=0;
	double form=0;             ///迭代法求解Z时Z的初值

	printf（“Please input the pressure（bar） and temperature（℃）!\n“）;	///输入压力和温度
	scanf（“%lf“&p）;
	getchar（）;
	scanf（“%lf“&t）;
	t=t+273.15;
	printf（“The pressure:%.2fbar\n“p）;
	printf（“The temperature:%.2fK\n“t）;
	for（i=0;i<=4;i++）
		tc=tc+hydrocarbon[i][0]*component[i];

	for（i=0;i<=4;i++）	///计算每种纯组分的a和b
	{
		ac=omigaa*pow（r2）*pow（hydrocarbon[i][0]2）/hydrocarbon[i][1];
		tr=t/hydrocarbon[i][0];
		k=0.37464+1.54226*hydrocarbon[i][2]-0.26992*pow（hydrocarbon[i][2]2）;
		q=pow（tr0.5）;
		q=1-q;
		q=1+k*q;
		q=pow（q2）;
		acom[i]=ac*q;
		bcom[i]=omigab*r*hydrocarbon[i][0]/hydrocarbon[i][1];
	}

	for（i=0;i<=4;i++）	///利用混合规则求解混合物的a和b
		b=b+bcom[i]*component[i];
	for（i=0;i<=4;i++）
		for（j=i;j<=4;j++）
			aij[i][j]=pow（acom[i]*acom[j]0.5）;
	for（i=0;i<=4;i++）
		for（j=0;j			aij[i][j]=aij[j][i];
	for（i=0;i<=4;i++）
		for（j=0;j<=4;j++）
			a=a+component[i]*component[j]*aij[i][j];
	a1=a*p/pow（r*t2）;
	b1=b*p/（r*t）;
	a2=b1-1;
	b2=a1-2*b1-3*b1*b1;
	c2=b1*b1+b1*b1*b1-a1*b1;
										///判断所