rksgasmodel.cpp

气体热力性质计算程序

#include "..\\includes\\RKSGasModel.h"

namespace XZGas
{    
	inline double RKSGasModel::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 RKSGasModel::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);
	}

    RKSGasModel::RKSGasModel(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;
    }
    
    RKSGasModel::~RKSGasModel()
    {
    }

    inline double RKSGasModel::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 RKSGasModel::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 RKSGasModel::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 RKSGasModel::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;
	}

	inline double RKSGasModel::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;
	}

	inline double RKSGasModel::s(double x,double y,Pid id)
	{		
		//the formula used here is ds=Cv/T*dT+dpdTv*dv
		double vv,TT;
        switch(id)
        {   
			//transform parameters according to the ideal gas equation -- pv=RgT
            case(T_P):
                TT=x;
                vv=v(TT,y);
            break;

            case(T_V):
                TT=x;
                vv=y;  
            break;
            
			case(P_V):
                TT=T(x,y);
                vv=y;
            break;
        }
		//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;
	}
		
    inline double RKSGasModel::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 RKSGasModel::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 RKSGasModel::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 w=pp*1e-10,d;
        do
        {
            _t1=a/vv/(vv+b)/(Rg*_t1/(vv-b)-pp);
            _t1*=_t1;
            d=_t1-_t2;
			_t2=_t1;
        }while(d<-1e-10||d>1e-10);
        return _t1;
    }

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

		return dvdTp/vv;
	}

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

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