rkgasmodel.cpp

气体热力性质计算程序

#include "..\\include\\RKGasModel.h"

namespace XZGas
{    
	inline double RKGasModel::dpdTv(double pp,double vT,double TT)
	{
		//the simplest form of dpdTv is 3*Rg/2/(v-b)-p/2/TT
		if(TT<0)	//the default value of TT is -1
			if(vT<20)	//if vT is too small, it carries the special volume instead of temperature
                TT=T(pp,vT);
			else	//in this range, vT cannot be the special volume of gas
			{
				TT=vT;
				vT=v(TT,pp);
			}
		
		return Rg*3/2/(vT-b)-pp/2/TT;
	}

	inline double RKGasModel::dpdvT(double TT,double vv)
	{
		double vb=vv+b,vvb=vv*vb;
		
		return (vb+b)*a/sqrt(TT)/vvb/vvb-Rg*TT/(vv-b)/(vv-b);
	}

    RKGasModel::RKGasModel(const GasData &gas):
    GasModel(gas)
    {
		a=0.42748*Rg*Rg/gas.pc*pow(gas.Tc,2.5);
        b=0.08664*Rg*gas.Tc/gas.pc;
    }
    
    RKGasModel::~RKGasModel()
    {
    }

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

		GasModel::ChangeGas(newgas);

		a=0.42748*Rg*Rg/newgas.pc*pow(newgas.Tc,2.5);
        b=0.08664*Rg*newgas.Tc/newgas.pc;
	}

    inline double RKGasModel::Cp(double TT,double pp)
    {
		//from the relationship between Cp and Cv: Cp-Cv=T*dvdTp*dpdTv
		double vv=v(TT,pp);
		double cv=Cv(TT,vv),dpdtv=dpdTv(pp,vv,TT);
		//-dvdTp, N means negtive
		double NdvdTp=dpdtv/dpdvT(TT,vv);

		return cv-TT*dpdtv*NdvdTp;
    }

	inline double RKGasModel::Cv(double TT,double vv)
	{
		//Cv=Cv(T)+T*IdpdTv2*dv
		double cv1=_Cp(TT)-Rg;
		
		double v0=v(T0,p0);
		double vvb=vv/(vv+b);
		vvb/=v0/(v0+b);
		double cv2=a*0.75/b*pow(TT,-1.5)*log(vvb);
		
		return cv1+cv2;
	}

	inline double RKGasModel::Gama(double TT,double pp)
	{
		double vv=v(TT,pp);
		double cv=Cv(TT,vv),dpdtv=dpdTv(pp,vv,TT);
		
		//-dvdTp, N means negtive
		double NdvdTp=dpdtv/dpdvT(TT,vv);
		double cp=cv-TT*dpdtv*NdvdTp;
		return cp/cv;
	}
	
	inline double RKGasModel::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 RKGasModel::u(double TT,double vv)
	{
		//from the formula: du=CvdT+(T*dpdTv-p)dv

		//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
		double v0=v(T0,p0);
		double vvb=vv/(vv+b);
		vvb/=v0/(v0+b);
		double u2=a*1.5/b/sqrt(TT)*log(vvb);
		
		return u1+u2+u0;
	}

	inline double RKGasModel::s(double TT,double pp)
	{		
		//the formula used here is ds=Cv/T*dT+dpdTv*dv
		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 has only one variable v
		double v0=v(T0,p0);
		double sv=a/2/pow(TT,1.5)/b*log(vv/(vv+b))+Rg*log(vv-b);
	
		double sv0=a/2/pow(TT,1.5)/b*log(v0/(v0+b))+Rg*log(v0-b);
		
		return s1+sv-sv0+s0;
	}
		
    inline double RKGasModel::p(double TT,double vv)
    {
        if(TT<=0||vv<=0)
            return 0.0;
        return (Rg*TT/(vv-b)-a/(sqrt(TT)*vv*(vv+b)));
    }

	inline double RKGasModel::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;
        double _v1=1.0,_v2=1.0,d;
        //double 1e-10;
		short n=0;
        do
        {
            _v1=pp+a/sqrt(TT)/_v1/(_v1+b);
            _v1=b+Rg*TT/_v1;
            d=_v1-_v2;
			++n;
			if(n>7)
				break;
			_v2=_v1;
        }while(d<-1e-10||d>1e10);

        return _v1;
    }
    
    inline double RKGasModel::T(double pp,double vv)
    {
		//the initial value of temperature for the iteration below
		//is assigned to 273.15K (0 celcius centigrade)
        double _T1=273.15,_T2=273.15;
        double d;

		//the calculation below is the Newton iteration on T:

		//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.
		
		double exprC1=a/vv/(vv+b);
		double exprC2=Rg/(vv-b);
		double f,df;
        do
        {
            f=pp-exprC2*_T2+exprC1/sqrt(_T2);
			df=-exprC2-exprC1/2/pow(_T2,1.5);
			_T1=_T2-f/df;
            d=_T1-_T2;
			_T2=_T1;
        }while(d<-1e-10||d>1e-10);
        return _T1;
    }

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

		return dvdTp/vv;
	}

	inline double RKGasModel::KT(double TT,double vv)
	{
		return -1/vv/dpdvT(TT,vv);
	}

	inline double RKGasModel::Beta(double pp,double vv)
	{
		return dpdTv(pp,vv)/pp;
	}
}