D-A-K方法求天然气偏差系数 以及等温压缩系数、天然气体积系数、天然气黏度的计算。牛顿迭代法。

D-A-K方法求天然气偏差系数 以及等温压缩系数、天然气体积系数、天然气黏度。牛顿迭代法。

这是第一次发文章，由于对该方法原理了解不是很透彻，没有考虑很多条件，仅给大家提供一些思路。

#include<stdio.h>
#include<math.h>
double Fcount(double a,double b,double pc,double tc,double pp ,double tt,double pg,double mg,double G,double B[])
	{
	    double P=G;
		double S,T,U,M,Z,Cp,Cg,Bg,Ug,K,X,Y;
		S=-0.27*a/b+P+(B[0]+B[1]/b+B[2]/pow(b,3)+B[3]/pow(b,4)+B[4]*pow(b,5))*pow(P,2)+(B[5]+B[6]/b+B[7]/pow(b,2))*pow(P,3)-
			B[8]*(B[6]/b+B[7]/pow(b,2))*pow(P,6)+B[9]*(1+B[10]*pow(P,2))*(pow(P,3)/pow(b,3))*exp(-B[10]*pow(P,2));
		T=1+2*(B[0]+B[1]/b+B[2]/pow(b,3)+B[3]/pow(b,4)+B[4]/pow(b,5))*P+3*(B[5]+B[6]/b+B[7]/pow(b,2))*pow(P,2)-6*B[8]*(B[6]/b+
			B[7]/pow(b,2))*pow(P,5)+(B[9]/pow(b,3))*(3*pow(P,2)+B[10]*(3*pow(P,4)-2*B[10]*pow(P,6)))*exp(-B[10]*pow(P,2));
		U=P-S/T;
		if(fabs(U-P)<0.0001)
		{
            Z=(0.27*a)/(P*b);
			M=(B[0]+B[1]/b+B[2]/pow(b,3)+B[3]/pow(b,4)+B[4]/pow(b,5))+2*(B[5]+B[6]/b+B[7]/pow(b,2))*P-5*B[8]*(B[6]/b+
				B[7]/pow(b,2))*pow(P,4)+2*B[9]*(P+B[10]*pow(P,3)-B[10]*B[10]*pow(P,5))/pow(b,3);
			Cp=1/a-0.27*(M/(1+b*M/Z))/(Z*Z*b);
			Cg=Cp/pc;
			Bg=3.458*0.0001*Z*tt/pp;
			K=(22.650+0.0388*mg)*pow(tt,1.5)/(209.2+19.26*mg+1.8*tt);
            X=3.448+548/tt+0.01*mg;
			Y=2.447-0.224*X;
            Ug=0.0001*K*exp(X*pow(pg,Y));
            printf("偏差系数Z的值=%lf\n天然气黏度Ug的值=%lf\n天然气气体系数Bg的值=%lf\n等温压缩系数Cg的值=%lf\n",Z,Ug,Bg,Cg);
		}
		else
			Fcount(a,b,pc,tc,pp,tt,pg,mg,U,B);
		return 0;
}
		
void main()
{

	double tr,pr,pc,tc,pp,tt,pg,mg;   
	double A[11]={0.3265,-1.01,-0.5339,0.01569,-0.05165,0.5475,-0.7361,0.1844,0.1056,0.6134,0.7210};
	printf("请输入拟对比压力pr:");
		scanf("%lf",&pr);
	printf("请输入拟对比温度tr:");
	    scanf("%lf",&tr);
	printf("请输入拟临界温度tc:");
	    scanf("%lf",&tc);
	printf("请输入拟临界压力pc:");
	    scanf("%lf",&pc);
	printf("请输入温度tt:");
	    scanf("%lf",&tt);
	printf("请输入压力pp:");
	    scanf("%lf",&pp);
    printf("请输入混合气体密度pg:");
	    scanf("%lf",&pg);
    printf("请输入混合气体视相对分子质量mg:");
	    scanf("%lf",&mg);
    double R=0.27*pr/tr;
    Fcount(pr,tr,pc,tc,pp,tt,pg,mg,R,A);


	
}

下面是简单的参数验证：

大家可以根据具体的条件继续修改。

有错误也请大家指正。