bbgasmodel.cpp

气体热力性质计算程序

#include "..\\include\\BBGasModel.h"

namespace XZGas
{
    BBGasModel::BBGasModel(const GasData &gas):
    GasModel(gas)
    {
		A0=gas.BBEquData.A0;
		B0=gas.BBEquData.B0;
		a=gas.BBEquData.a;
		b=gas.BBEquData.b;
		c=gas.BBEquData.c;
    }
    
    BBGasModel::~BBGasModel()
    {
    }

	inline double BBGasModel::dpdvT(double TT,double vv)
	{
		double dpdvt = 0.0,TT2=TT*TT;
		dpdvt-=4*b*B0*Rg*c/TT2;
		dpdvt/=vv;		
		dpdvt-=3*(A0*a-B0*Rg*c/TT2-B0*b*Rg*TT);
		dpdvt/=vv;
		dpdvt-=2*((B0*TT-c/TT2)*Rg-A0);
		dpdvt/=vv;
        dpdvt-=Rg*TT;
		dpdvt/=vv;
		dpdvt/=vv;
		
		return dpdvt;
	}

	inline double BBGasModel::dpdTv(double TT,double vv)
	{
		double dpdtv=0.0,T3=TT*TT*TT;
        dpdtv-=2*b*B0*Rg*c/T3;
		dpdtv/=vv;
		dpdtv+=B0*Rg*(2*c/T3-b);
		dpdtv/=vv;
		dpdtv+=Rg*(B0+2*c/T3);
		dpdtv/=vv;
		dpdtv+=Rg;
		dpdtv/=vv;
		
		return dpdtv;
	}

    inline double BBGasModel::Cp(double TT,double pp)
    {
		double cp=0.0;
		double vv=v(TT,pp);
		double cv=Cv(TT,vv),dpdtv=dpdTv(TT,vv);
		//-dvdTp, N means negtive
		double NdvdTp=dpdtv/dpdvT(TT,vv);

		return cv-TT*dpdtv*NdvdTp;
    }
    
    inline double BBGasModel::Cv(double TT,double vv)
    {
		//Cv=Cv(T)+T*IdpdTv2*dv
		double cv1=_Cp(TT)-Rg;
	
		//the second part of Cv which is close ralated to the change of volume
		double cv2=0.0,v0=v(T0,p0);
		double v_=1/vv,v0_=1/v0;
		cv2+=v_-v0_;
		v_/=vv;
		v0_/=v0;
		cv2+=((v_-v0_)*B0/2);
		v_/=vv;
		v0_/=v0;
		cv2-=(b*B0/3*(v_-v0_));
		cv2*=6*Rg*c/TT/TT/TT;

		return cv1+cv2;
    }	

	inline double BBGasModel::Gama(double TT,double pp)
	{
		double vv=v(TT,pp);
		double cv=Cv(TT,vv),dpdtv=dpdTv(TT,vv);
		
		//-dvdTp, N means negtive
		double NdvdTp=dpdtv/dpdvT(TT,vv);
		double cp=cv-TT*dpdtv*NdvdTp;
		return cp/cv;
	}

	inline double BBGasModel::h(double TT,double pp)
	{
		//h=u+pv
		double vv=v(TT,pp),v0=v(T0,p0);
		double uu=u(TT,vv);
			
		return uu+pp*vv-p0*v0+h0;
	}

	inline double BBGasModel::u(double TT,double vv)
	{
		double v0=v(T0,p0);
		//the first part of integration above temperature: ICvdT
		double u1=_ICp(TT)-_ICp(T0)-Rg*(TT-T0);

		//the second part of calculus above special volume: I(T*dpdTv-p)dv
		
		//the expression of Rg*C/T^2 which is offen used in the calculation below
		double RcT2=Rg*c/TT/TT; 
		//the gadget used to reduce the times of multiplication
		double v_=1/vv,v0_=1/v0;
		double u2=0.0;
		u2-=(3.0*RcT2+A0)*(v_-v0_);

		v_/=vv;
		v0_/=v0;
		u2-=(1.5*B0*RcT2+a*A0/2)*(v_-v0_);

		v_/=vv;
		v0_/=v0;
		u2+=B0*b*RcT2*(v_-v0_);

		return u1+u2+u0;
	}

	inline double BBGasModel::s(double TT,double pp)
	{
		double vv=v(TT,pp);

		//the first part of entropy ICv/TdT
		double s1=_ICpT(TT)-_ICpT(T0)-Rg*log(TT/T0);
		
		//the seconde part of entropy, which is related to one variable v
		double s2=0.0;
		double v0=v(T0,p0);
        double T3=TT*TT*TT;
		s2+=2*c/T3*(vv-v0);
		s2+=(2*c*B0/T3+1)*log(vv/v0);
		s2-=(1-2*c*b/T3)*B0*(1/vv-1/v0);
		s2+=B0*b/2*(1/vv/vv-1/v0/v0);
		s2*=Rg;

		return s1+s2+s0;
	}

	inline double BBGasModel::p(double TT,double vv)
	{
		double _p=0.0;
		
		double T3=TT*TT*TT,tmp=Rg*TT;
		
		double BB=B0-A0/tmp-c/T3;
		double CC=A0*a/tmp-B0*b-c*B0/T3;
		double DD=B0*b*c/T3;

		_p=DD/vv;
		_p+=CC;
		_p/=vv;
		_p+=BB;
		_p/=vv;
		_p+=1;
		_p*=(tmp/vv);

		return _p;
	}

	inline double BBGasModel::v(double TT,double pp)
	{
		//the special volume is worked out by iteration, which converges very fast
		//the initial value of special volume is assigned to 1.0 m^3/kg
        if(TT<=0||pp<=0)
            return 0.0;
		//the context below is the Newtow Iteration on the density of gas:

        double _den1=2.0,_den2=2.0;	//the density of the gas
		double d;
        double T3=TT*TT*TT,tmp=Rg*TT;
		double BB=B0-A0/tmp-c/T3;
		double CC=A0*a/tmp-B0*b-c*B0/T3;
		double DD=B0*b*c/T3;
		double df=0.0,f=0.0;
        do
        {
            df=4*DD*_den2;
			df+=3*CC;
			df*=_den2;
			df+=2*BB;
			df*=_den2;
			df+=1;
			df*=tmp;

			f=DD*_den2;
			f+=CC;
			f*=_den2;
			f+=BB;
			f*=_den2*_den2;
			f+=_den2;
			f*=tmp;
			f-=pp;

			_den1=_den2-f/df;
			d=_den1-_den2;
			_den2=_den1;
        }while(d>1e-10||d<-1e-10);

        return 1/_den1;
	}

	inline double BBGasModel::T(double pp,double vv)
    {
		//the initial value of temperature for the iteration below
		//is assigned to 298.15K (25 centigrade)
    
		double _T1=298.15,_T2=298.15;
        double d;
		unsigned int n=0;
		//because the value of vv and pp will be constant in the process 
		//of iteration, the expression below will not change. So the final
		//value will be constant. The variable exprC1 can be used to reduce
		//the computation. And so is exprC2.
				
		//problem provoking!!!!!!!!!!!!!!!!!!!!!!! this iteration may deverge
		//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		double exprC1=(pp*vv*vv+A0*(1-a/vv))/Rg/(vv+B0*(1-b/vv));
		double exprC2=c/vv;
		do
        {
			_T1=exprC1+exprC2/_T2/_T2;
			_T1=(_T1+_T2)/2;
            d=_T1-_T2;
			_T2=_T1;
			n++;
			if(n>500)
				break;
        }while(d<-1e-10||d>1e-10);
        return _T1;
    }

	inline double BBGasModel::Alphav(double pp,double vv)
	{
		double TT=T(pp,vv);
		double dvdTp=dpdTv(TT,vv)/dpdvT(TT,vv);

		return dvdTp/vv;
	}

	inline double BBGasModel::KT(double TT,double vv)
	{
		return -1/vv/dpdvT(TT,vv);
	}
	inline double BBGasModel::Beta(double pp,double vv)
	{
		double TT=T(pp,vv);
		return dpdTv(TT,vv)/pp;
	}

	inline void BBGasModel::ChangeGas(const GasData &newgas)
	{
		//the same kind of gas
		if(strcmp(GasName,newgas.Name) == 0)
			return ;

		GasModel::ChangeGas(newgas);

		A0=newgas.BBEquData.A0;
		B0=newgas.BBEquData.B0;
		a=newgas.BBEquData.a;
		b=newgas.BBEquData.b;
		c=newgas.BBEquData.c;
	}
}