cantera/ext/tpx/Sub.cpp

/*
 * The Substance class
 * D. Goodwin, Caltech Nov. 1996
 */
#include "Sub.h"
#include <math.h>
#include <fstream>
#include <stdio.h>

using std::string;

namespace tpx {

    static string fp2str(double x, string fmt="%g") {
        char buf[30];
        sprintf(buf, fmt.c_str(), x);
        return string(buf);
    }

/*
    static string int2str(int n, string fmt="%d") {
        char buf[30];
        sprintf(buf, fmt.c_str(), n);
        return string(buf);
    }
*/

    string TPX_Error::ErrorMessage = "";
    string TPX_Error::ErrorProcedure = "";

    string errorMsg(int flag) {
        switch(flag) {
        case NoConverge: return "no convergence";
        case GenError:   return "general error";
        case InvalidInput:   return "invalid input";
        case TempError:   return "temperature error";
        case PresError:   return "pressure error";
        default: return "(unknown error)";
        }
    }

    //-------------- Public Member Functions --------------

    /// Pressure [Pa]. If two phases are present, return the
    /// saturation pressure; otherwise return the pressure
    /// computed directly from the underlying eos.
    double Substance::P() {
	double ppp = (TwoPhase() ? Ps() : Pp());
	return (Err ? Undef : ppp);
    }

    const double DeltaT = 0.000001;

    /// The derivative of the saturation pressure 
    /// with respect to temperature.
    double Substance::dPsdT() {
	double ps1, tsave, dpdt;
	tsave = T;
	ps1 = Ps();
	set_T(T + DeltaT);
	dpdt = (Ps() - ps1)/DeltaT;
	set_T(tsave);
	return (Err ? Undef : dpdt);
    }

    /// true if a liquid/vapor mixture, false otherwise
    int Substance::TwoPhase() {
	if (T >= Tcrit()) return 0;
	update_sat();
	return ((Rho < Rhf) && (Rho > Rhv) ? 1 : 0);
}

    /// Vapor fraction. 
    /// If T >= Tcrit, 0 is returned for v < Vcrit, and 1 is
    /// returned if v > Vcrit.
    double Substance::x() {
	double xx, vv, vl;
	if (T >= Tcrit())
            return (1.0/Rho < Vcrit() ? 0.0 : 1.0);
	else {
            update_sat();
            if (Rho <= Rhv)
                xx = 1.0;
            else if (Rho >= Rhf)
                xx = 0.0;
            else {
                vv = 1.0/Rhv;
                vl = 1.0/Rhf;
                xx = (1.0/Rho - vl)/(vv - vl);
            }
            return (Err ? Undef : xx);
	}
    }
    
    /// Saturation temperature at pressure p.
    double Substance::Tsat(double p) {
	double Tsave, p_here, dp, dt, dpdt, dta,
            dtm, tsat;
	if (Err || (p <= 0.0) || (p > Pcrit())) {
            throw TPX_Error("Substance::Tsat","illegal pressure value");
            set_Err(PresError);
            return Undef;
	}
	int LoopCount = 0;
	double tol = 1.e-6*p;
	Tsave = T;
	if (T < Tmin()) T = 0.5*(Tcrit() - Tmin());
	if (T >= Tcrit()) T = 0.5*(Tcrit() - Tmin());
        do {
            if (Err) break;
            if (T > Tcrit()) T = Tcrit() - 0.001;
            if (T < Tmin()) T = Tmin() + 0.001;
            p_here = Ps();
            dpdt = dPsdT();
            dp = p - p_here;
            dt = dp/dpdt;
            dta = fabs(dt);
            dtm = 0.1*T;
            if (dta > dtm) dt = dt*dtm/dta;
            T = T + dt;
            LoopCount++;
            if (LoopCount > 100) {
                T = Tsave;
                set_Err(NoConverge);
                return Undef;
            }
	} while (fabs(dp) > tol);
	tsat = T;
	T = Tsave;
	return (Err ? Undef : tsat);
    }  
    
    
    // absolute tolerances
    
    double TolAbsH = 0.0001;                       // J/kg
    double TolAbsU = 0.0001;
    double TolAbsS = 1.e-7;
    double TolAbsP = 0.000;                      // Pa
    double TolAbsV = 1.e-8;
    double TolAbsT = 1.e-3;
    double TolRel = 3.e-8;
    
    void Substance::Set(int XY, double x0, double y0) {
	double temp;
	
	clear_Err();  // clear error flag

	/* if inverted (PT) switch order and change sign of XY (TP = -PT) */
	if (XY < 0) {
            double tmp = x0;
            x0 = y0;
            y0 = tmp;
            XY *= -1;
	}
	
	switch(XY) {
	case TV:
            set_T(x0); set_v(y0); break;
            
	case HP:
            if (Lever(Pgiven, y0, x0, EvalH)) return;
            set_xy(EvalH, EvalP, x0, y0, TolAbsH, TolAbsP, TolRel, TolRel);
            break;
            
	case SP:
            if (Lever(Pgiven, y0, x0, EvalS)) return;
            set_xy(EvalS, EvalP, x0, y0, TolAbsS, TolAbsP, TolRel, TolRel);
            break;
            
	case PV:
            if (Lever(Pgiven, x0, y0, EvalV)) return;
            set_xy(EvalP, EvalV, x0, y0, TolAbsP, TolAbsV, TolRel, TolRel);
            break;
            
	case TP:
            if (x0 < Tcrit()) {
                set_T(x0);
                if (y0 < Ps()) {
                    Set(TX, x0, Vapor);
                }
                else {
                    Set(TX, x0, Liquid);
                }
            }
            else
                set_T(x0);
            set_xy(EvalT, EvalP, x0, y0, TolAbsT, TolAbsP, TolRel, TolRel);
            break;    
            
	case UV:
            set_xy(EvalU, EvalV, x0, y0, TolAbsU, TolAbsV, TolRel, TolRel);
            break;
            
	case ST:
            if (Lever(Tgiven, y0, x0, EvalS)) return;
            set_xy(EvalS, EvalT, x0, y0, TolAbsS, TolAbsT, TolRel, TolRel);
            break;
            
	case SV:
            set_xy(EvalS, EvalV, x0, y0, TolAbsS, TolAbsV, TolRel, TolRel);
            break;
            
	case UP:
            if (Lever(Pgiven, y0, x0, EvalU)) return;
            set_xy(EvalU, EvalP, x0, y0, TolAbsU, TolAbsP, TolRel, TolRel);
            break;
            
	case VH:
            set_xy(EvalV, EvalH, x0, y0, TolAbsV, TolAbsH, TolRel, TolRel);
            break;
            
	case TH:
            set_xy(EvalT, EvalH, x0, y0, TolAbsT, TolAbsH, TolRel, TolRel);
            break;
            
	case SH:
            set_xy(EvalS, EvalH, x0, y0, TolAbsS, TolAbsH, TolRel, TolRel);
            break;
            
	case PX:  
            temp = Tsat(x0);
            if ((y0 >= 0.0) && (y0 <= 1.0) && (temp < Tcrit())) {
                set_T(temp);
                update_sat();
                Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
            }
            else 
                set_Err(InvalidInput);
            break;
            
	case TX:  
            if ((y0 >= 0.0) && (y0 <= 1.0) && (x0 < Tcrit())) {
                set_T(x0);
                update_sat();
                Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
            }
            else {
                set_Err(InvalidInput);
            }
            break;
            
	default:
            set_Err(InvalidInput);
	}
	if (Err) {
            T = Undef; 
            Rho = Undef;
            Tslast = Undef; 
            Rhf = Undef; 
            Rhv = Undef; 
	}
    }

    //------------------ Protected and Private Functions -------------------

    void Substance::set_Rho(double r0) {
	if (r0 > 0.0) {
            Rho = r0;
	}
	else {
            set_Err(InvalidInput);
        }
    }
    
    void Substance::set_T(double t0) {
	if ((t0 >= Tmin()) && (t0 <= Tmax())) {
            T = t0;
	}
	else {
            throw TPX_Error("Substance::set_T",
                "illegal temperature value "+fp2str(t0));
            set_Err(TempError);
	}
    }
    
    void Substance::set_v(double v0) {
	if (v0 > 0) { 
            Rho = 1.0/v0;
	}
	else {
            throw TPX_Error("Substance::set_v",
                "negative specific volume: "+fp2str(v0));
            set_Err(InvalidInput);
        }
    }
    
    double Substance::Ps() {
	if (T < Tmin() || T > Tcrit()) {
            throw TPX_Error("Substance::Ps",
                "illegal temperature value "+fp2str(T));
            set_Err(TempError);
            return Undef;
	}
	update_sat();
        return Pst;
    }

    // update saturated liquid and vapor densities and saturation pressure
    
    void Substance::Set_meta(double phase, double pp) {
	if (phase == Liquid) 
            Rho = ldens();                       // trial value = liquid dens
	else
            Rho = pp*MolWt()/(8314.0*T);         // trial value = ideal gas
	set_TPp(T, pp);	
    }

    void Substance::update_sat() {
	if ((T != Tslast) && (T < Tcrit())) {
            double Rho_save = Rho;
            double gf, gv, dg, dp, dlp, psold;
            
            double pp = Psat();
            double lps = log(pp);
            // trial value = Psat from correlation
            int i;
            
            for (i = 0; i<20; i++) {
                Rho = ldens();                // trial value = liquid density
                set_TPp(T,pp);
                Rhf = Rho;                    // sat liquid density

                gf = hp() - T*sp();
                Rho = pp*MolWt()/(8314.0*T);  // trial value = ideal gas
                set_TPp(T,pp);

                Rhv = Rho;                    // sat vapor density
                gv = hp() - T*sp();
                dg = gv - gf;
                
                if (fabs(dg) < 0.001) break;
                dp = dg/(1.0/Rhv - 1.0/Rhf);
                psold = pp;
                
                if (fabs(dp) > pp) {
                    dlp = dg/(pp*(1.0/Rhv - 1.0/Rhf));
                    lps -= dlp;
                    pp = exp(lps);
                }
                else {
                    pp -= dp;
                    lps = log(pp); // added 10/5/04
                }
                if (pp > Pcrit()) {
                    pp = psold + 0.5*(Pcrit() - psold);
                    lps = log(pp); // added 10/5/04
                }
                else if (pp < 0.0) {
                    pp = psold/2.0;
                    lps = log(pp); // added 10/5/04
                }
            }	
            if (Rhf <= Rhv) {
                throw TPX_Error("Substance::update_sat",
                    "wrong root found for sat. liquid or vapor at P = "+fp2str(pp));
            }

            if (i >= 20) {
                Pst = Undef;
                Rhv = Undef;
                Rhf = Undef;
                Tslast = Undef;
                throw TPX_Error("substance::update_sat","no convergence");
                set_Err(NoConverge);
            }
            else {
                Pst = pp;
                Tslast = T;
            }
            Rho = Rho_save;
	}
    }

    double Substance::vprop(int ijob) {
	switch (ijob) {
	case EvalH: return hp();
	case EvalS: return sp();
	case EvalU: return up();
	case EvalV: return vp();
	case EvalP: return Pp();
	default: return Undef;
	}
    }

    int Substance::Lever(int itp, double sat, double val, int ifunc) {
        /*
         * uses lever rule to set state in the dome. Returns 1 if in dome,
         * 0 if not, in which case state not set.
         */
	double Valf, Valg, Tsave, Rhosave, xx, vv, psat;
	Tsave = T;
	Rhosave = Rho;
	if (itp == Tgiven) {
            if (sat >= Tcrit()) return 0;
            set_T(sat);
            psat = Ps();
	}
	else if (itp == Pgiven) {
            if (sat >= Pcrit()) return 0;
            psat = sat;
            T = Tsat(psat);
	    if (T == Undef) {
	      Err = 0;
	      T = Tsave;
	      Rho = Rhosave;
	      return 0;
	    }
	}
	else {
            throw TPX_Error("Substance::Lever","general error");
            set_Err(GenError);
            return GenError;
	}
	Set(TX, T, Vapor);
	Valg = vprop(ifunc);
	Set(TX, T, Liquid);
	Valf = vprop(ifunc);
	if (Err)
            return Err;
	else if ((val >= Valf) && (val <= Valg)) {
            xx = (val - Valf)/(Valg - Valf);
            vv = (1.0 - xx)/Rhf + xx/Rhv;
            set_v(vv);
            return 1;
	}
	else {
            T = Tsave;
            Rho = Rhosave;
            return 0;
	}
    }	


    void Substance::set_xy(int ifx, int ify, double X, double Y, 
        double atx, double aty, 
        double rtx, double rty) {
        
	double v_here, t_here, dv, dt, dxdt, dydt, dxdv, dydv,
            det, x_here, y_here, dvm, dtm, dva, dta;
	double Xa, Ya, err_x, err_y;
	double dvs1 = 2.0*Vcrit();
	double dvs2 = 0.7*Vcrit();
	int LoopCount = 0;
	
	double v_save = 1.0/Rho;
	double t_save = T;
	
	if (Err) return;
	
	if ((T == Undef) && (Rho == Undef)) {  // new object, try to pick
            Set(TV,Tcrit()*1.1,Vcrit()*1.1);   // "reasonable" starting point
            t_here = T;
            v_here = 1.0/Rho;
	}
	else if (Rho == Undef) {  // new object, try to pick
            Set(TV,T,Vcrit()*1.1);   // "reasonable" starting point
            t_here = T;
            v_here = 1.0/Rho;
	}
	else {
            v_here = v_save;
            t_here = t_save;
	}
        
	Xa = fabs(X);
	Ya = fabs(Y);

      
        // loop
	do {
            if (Err) break;
            x_here = prop(ifx);
            y_here = prop(ify);
            err_x = fabs(X - x_here);
            err_y = fabs(Y - y_here);

            if ((err_x < atx + rtx*Xa) && (err_y < aty + rty*Ya)) break;
            
            /* perturb t */
            dt = 0.001*t_here;
            if (t_here + dt > Tmax())
                dt *= -1.0;
            
            /* perturb v */
            dv = 0.001*v_here;
            if (v_here <= Vcrit()) dv *= -1.0;
            
            /* derivatives with respect to T */
            Set(TV, t_here + dt, v_here);
            dxdt = (prop(ifx) - x_here)/dt;
            dydt = (prop(ify) - y_here)/dt;
            
            /* derivatives with respect to v */
            Set(TV, t_here, v_here + dv);
            dxdv = (prop(ifx) - x_here)/dv;
            dydv = (prop(ify) - y_here)/dv;
            
            det = dxdt*dydv - dydt*dxdv;
            dt = ((X - x_here)*dydv - (Y - y_here)*dxdv)/det;
            dv = ((Y - y_here)*dxdt - (X - x_here)*dydt)/det;
            
            dvm = 0.2*v_here;
            if (v_here < dvs1) dvm *= 0.5;
            if (v_here < dvs2) dvm *= 0.5;
            dtm = 0.1*t_here;
            dva = fabs(dv);
            dta = fabs(dt);
            if (dva > dvm) dv *= dvm/dva;
            if (dta > dtm) dt *= dtm/dta;
            v_here += dv;
            t_here += dt;
            if (t_here >= Tmax())
                t_here = Tmax() - 0.001;
            else if (t_here <= Tmin())
                t_here = Tmin() + 0.001;
            if (v_here <= 0.0)
                v_here = 0.0001;
            Set(TV, t_here, v_here);
            LoopCount++;
            if (LoopCount > 200) {
                throw TPX_Error("Substance::set_xy","no convergence");
                set_Err(NoConverge); 
                break;
            }
	} while(1);
    }


    double Substance::prop(int ijob) {
	double xx, pp, lp, vp, Rho_save;
	if (ijob == EvalP) return P();
	if (ijob == EvalT) return T;
	xx = x();
	if ((xx > 0.0) && (xx < 1.0)) {
            Rho_save = Rho;
            Rho = Rhv;
            vp = vprop(ijob);
            Rho = Rhf;
            lp = vprop(ijob);
            pp = (1.0 - xx)*lp + xx*vp;
            Rho = Rho_save;
            return pp;
	}
	else
            return vprop(ijob);
    }	
    
    int kbr;
    const double ErrP = 1.e-7;
    const double Big = 1.e30;
    
    double Vmin, Vmax, Pmin, Pmax, dvs1, dvs2, dpdv, dvbf, dv, dva, dvm, 
        dt, v_here, P_here, vt;
    
    void Substance::BracketSlope(double Pressure) {
	if (kbr == 0) {
            dv = (v_here < Vcrit() ? -0.05*v_here : 0.2*v_here);
            if (Vmin > 0.0)
                dv = 0.2*v_here;
            if (Vmax < Big)
                dv = -0.05*v_here;
	}
	else {
            double dpdv = (Pmax - Pmin)/(Vmax - Vmin);
            v_here = Vmax;
            P_here = Pmax;
            dv = dvbf*(Pressure - P_here)/dpdv;
            dvbf = 0.5*dvbf;
	}
    }
    
    void Substance::set_TPp(double Temp, double Pressure) {
	kbr = 0;
	dvbf = 1.0;
	Vmin = 0.0;
	Vmax = Big;
	Pmin = Big;
	Pmax = 0.0;
	dvs1 = 2.0*Vcrit();
	dvs2 = 0.7*Vcrit();
	int LoopCount = 0;
	
	double v_save = 1.0/Rho;
	set_T(Temp);
	v_here = vp();
	
	// loop
	
	while(P_here = Pp(), fabs(Pressure - P_here) >= ErrP*Pressure) {
            if (P_here < 0.0) {
                BracketSlope(Pressure);
            }
            else {
                dv = 0.001*v_here;
                if (v_here <= Vcrit()) {
                    dv *= -1.0;
                }
                Set(TV, Temp, v_here+dv);
                dpdv = (Pp() - P_here)/dv;
                if (dpdv > 0.0) {
                    BracketSlope(Pressure);
                }
                else {
                    if ((P_here > Pressure) && (v_here > Vmin))
                        Vmin = v_here;
                    else if ((P_here < Pressure) && (v_here < Vmax))
                        Vmax = v_here;
                    if (v_here == Vmin) Pmin = P_here;
                    if (v_here == Vmax) Pmax = P_here;
                    if (Vmin >= Vmax) {
                        throw TPX_Error("Substance::set_TPp","Vmin >= Vmax");
                        set_Err(GenError); 
                    }
                    else if ((Vmin > 0.0) && (Vmax < Big)) kbr = 1;
                    dvbf = 1.0;
                    if (dpdv == 0.0) {
                        dvbf = 0.5;
                        BracketSlope(Pressure);
                    }
                    else {
                        dv = (Pressure - P_here)/dpdv;
                    }
                }
            }
            dvm = 0.2*v_here;
            if (v_here < dvs1) dvm *= 0.5;
            if (v_here < dvs2) dvm *= 0.5;
            if (kbr != 0) {
                vt = v_here + dv;
                if ((vt < Vmin) || (vt > Vmax))
                    dv = Vmin + (Pressure - Pmin)*(Vmax - Vmin)/(Pmax - Pmin) - v_here;
            }
            dva = fabs(dv);
            if (dva > dvm) dv *= dvm/dva;
            v_here += dv;
            if (dv == 0.0) {
                throw TPX_Error("Substance::set_TPp","dv = 0 and no convergence");
                set_Err(NoConverge); 
                return;
            }
            Set(TV, Temp, v_here);
            LoopCount++;
            if (LoopCount > 100) {       
                Set(TV, Temp, v_save);
                throw TPX_Error("Substance::set_TPp",string("no convergence for ")
                    +"P* = "+fp2str(Pressure/Pcrit())+". V* = "
                    +fp2str(v_save/Vcrit()));
                set_Err(NoConverge);
                return;
            }
	}
	Set(TV, Temp,v_here);
    }
}