CoolPropFluid.h

/*
 * CoolPropFluid.h
 *
 *  Created on: 20 Dec 2013
 *      Author: jowr
 */

#ifndef COOLPROPFLUID_H_
#define COOLPROPFLUID_H_

#include "DataStructures.h"
#include "Helmholtz.h"
#include "Solvers.h"

#include <numeric>
#include <string>
#include <vector>
#include <map>
#include <cassert>
#include <iterator>
#include "Eigen/Core"
#include "PolyMath.h"
#include "Ancillaries.h"

namespace CoolProp {

struct BibTeXKeysStruct
{
    std::string EOS, CP0, VISCOSITY, CONDUCTIVITY, ECS_LENNARD_JONES, ECS_FITS, SURFACE_TENSION;
};

struct EnvironmentalFactorsStruct
{
    double GWP20, GWP100, GWP500, ODP, HH, PH, FH;
    std::string ASHRAE34;
};
struct CriticalRegionSplines
{
    double T_min, T_max, rhomolar_min, rhomolar_max;
    std::vector<double> cL, cV;
    bool enabled;
    CriticalRegionSplines() : T_min(_HUGE), T_max(_HUGE), rhomolar_min(_HUGE), rhomolar_max(_HUGE), enabled(false){};

    const void get_densities(double T, double rho_min, double rho_crit, double rho_max, double& rhoL, double& rhoV) const {
        int Nsoln = -1, Ngood = 0;
        double rho1 = 0, rho2 = 0, rho3 = 0;

        // -----------
        // Liquid part
        // -----------
        Ngood = 0;
        solve_cubic(cL[0], cL[1], cL[2], cL[3] - T, Nsoln, rho1, rho2, rho3);
        if (Nsoln == 1 && rho1 < rho_max && rho1 > rho_crit) {
            rhoL = rho1;
        } else {
            if (rho1 < rho_max && rho1 > rho_crit) {
                Ngood++;
                rhoL = rho1;
            }
            if (rho2 < rho_max && rho2 > rho_crit) {
                Ngood++;
                rhoL = rho2;
            }
            if (Nsoln > 2 && rho3 < rho_max && rho3 > rho_crit) {
                Ngood++;
                rhoL = rho3;
            }
            if (Ngood > 1) {
                throw ValueError(format("More than one liquid solution found for critical spline for T=%0.12g", T));
            };
            if (Ngood < 1) {
                throw ValueError(format("No liquid solution found for critical spline for T=%0.12g", T));
            };
        }

        // ----------
        // Vapor part
        // ----------
        Ngood = 0;
        Nsoln = 0;
        solve_cubic(cV[0], cV[1], cV[2], cV[3] - T, Nsoln, rho1, rho2, rho3);
        if (Nsoln == 1 && rho1 > rho_min && rho1 < rho_crit) {
            rhoV = rho1;
        } else {
            if (rho1 > rho_min && rho1 < rho_crit) {
                Ngood++;
                rhoV = rho1;
            }
            if (rho2 > rho_min && rho2 < rho_crit) {
                Ngood++;
                rhoV = rho2;
            }
            if (Nsoln > 2 && rho3 > rho_min && rho3 < rho_crit) {
                Ngood++;
                rhoV = rho3;
            }
            if (Ngood > 1) {
                throw ValueError(format("More than one vapor solution found for critical spline for T=%0.12g", T));
            };
            if (Ngood < 1) {
                throw ValueError(format("No vapor solution found for critical spline for T=%0.12g", T));
            };
        }
    };
};

struct EOSLimits
{
    double Tmin, Tmax, rhomax, pmax;
};

struct ConductivityECSVariables
{
    std::string reference_fluid;
    CoolPropDbl psi_rhomolar_reducing, f_int_T_reducing;
    std::vector<CoolPropDbl> psi_a, psi_t, f_int_a, f_int_t;
};

struct ConductivityDiluteEta0AndPolyData
{
    std::vector<CoolPropDbl> A, t;
};

struct ConductivityDiluteRatioPolynomialsData
{
    CoolPropDbl T_reducing, p_reducing;
    std::vector<CoolPropDbl> A, B, n, m;
};
struct ConductivityDiluteVariables
{
    enum ConductivityDiluteEnum
    {
        CONDUCTIVITY_DILUTE_RATIO_POLYNOMIALS,
        CONDUCTIVITY_DILUTE_ETA0_AND_POLY,
        CONDUCTIVITY_DILUTE_CO2,
        CONDUCTIVITY_DILUTE_CO2_HUBER_JPCRD_2016,
        CONDUCTIVITY_DILUTE_ETHANE,
        CONDUCTIVITY_DILUTE_NONE,
        CONDUCTIVITY_DILUTE_NOT_SET
    };
    int type;
    ConductivityDiluteRatioPolynomialsData ratio_polynomials;
    ConductivityDiluteEta0AndPolyData eta0_and_poly;

    ConductivityDiluteVariables() {
        type = CONDUCTIVITY_DILUTE_NOT_SET;
    }
};

struct ConductivityResidualPolynomialAndExponentialData
{
    CoolPropDbl T_reducing, rhomass_reducing;
    std::vector<CoolPropDbl> A, t, d, gamma, l;
};

struct ConductivityResidualPolynomialData
{
    CoolPropDbl T_reducing, rhomass_reducing;
    std::vector<CoolPropDbl> B, t, d;
};
struct ConductivityResidualVariables
{
    enum ConductivityResidualEnum
    {
        CONDUCTIVITY_RESIDUAL_POLYNOMIAL,
        CONDUCTIVITY_RESIDUAL_POLYNOMIAL_AND_EXPONENTIAL,
        CONDUCTIVITY_RESIDUAL_R123,
        CONDUCTIVITY_RESIDUAL_CO2,
        CONDUCTIVITY_RESIDUAL_NOT_SET
    };
    int type;
    ConductivityResidualPolynomialData polynomials;
    ConductivityResidualPolynomialAndExponentialData polynomial_and_exponential;

    ConductivityResidualVariables() {
        type = CONDUCTIVITY_RESIDUAL_NOT_SET;
    }
};

struct ConductivityCriticalSimplifiedOlchowySengersData
{
    CoolPropDbl k, R0, gamma, nu, GAMMA, zeta0, qD, T_reducing, p_reducing, T_ref;
    ConductivityCriticalSimplifiedOlchowySengersData()
      :                    // Universal constants - can still be adjusted if need be
        k(1.3806488e-23),  //[J/K]
        R0(1.03),          //[-]
        gamma(1.239),      //[-]
        nu(0.63),          //[-]
        // Suggested default values - can be over-written
        GAMMA(0.0496),    //[-]
        zeta0(1.94e-10),  //[m]
        qD(2e9),          //[m]
        // Set to invalid number, can be provided in the JSON file
        // T_ref default is 1.5*Tc
        T_reducing(_HUGE),
        p_reducing(_HUGE),
        T_ref(_HUGE) {}
};
struct ConductivityCriticalVariables
{
    enum ConductivityResidualEnum
    {
        CONDUCTIVITY_CRITICAL_SIMPLIFIED_OLCHOWY_SENGERS,
        CONDUCTIVITY_CRITICAL_R123,
        CONDUCTIVITY_CRITICAL_AMMONIA,
        CONDUCTIVITY_CRITICAL_NONE,
        CONDUCTIVITY_CRITICAL_CARBONDIOXIDE_SCALABRIN_JPCRD_2006,
        CONDUCTIVITY_CRITICAL_NOT_SET
    };
    int type;
    ConductivityCriticalSimplifiedOlchowySengersData Olchowy_Sengers;

    ConductivityCriticalVariables() {
        type = CONDUCTIVITY_CRITICAL_NOT_SET;
    }
};

struct ViscosityDiluteGasCollisionIntegralData
{
    CoolPropDbl molar_mass, C;
    std::vector<CoolPropDbl> a, t;
};
struct ViscosityDiluteCollisionIntegralPowersOfTstarData
{
    CoolPropDbl T_reducing,  
      C;                     
    std::vector<CoolPropDbl> a, t;
};
struct ViscosityDiluteGasPowersOfT
{
    std::vector<CoolPropDbl> a, t;
};
struct ViscosityDiluteGasPowersOfTr
{
    std::vector<CoolPropDbl> a, t;
    CoolPropDbl T_reducing;
};
struct ViscosityDiluteVariables
{
    enum ViscosityDiluteType
    {
        VISCOSITY_DILUTE_COLLISION_INTEGRAL,                  
        VISCOSITY_DILUTE_COLLISION_INTEGRAL_POWERS_OF_TSTAR,  
        VISCOSITY_DILUTE_KINETIC_THEORY,                      
        VISCOSITY_DILUTE_ETHANE,                              
        VISCOSITY_DILUTE_CYCLOHEXANE,                         
        VISCOSITY_DILUTE_CO2_LAESECKE_JPCRD_2017,             
        VISCOSITY_DILUTE_POWERS_OF_T,                         
        VISCOSITY_DILUTE_POWERS_OF_TR,                        
        VISCOSITY_DILUTE_NOT_SET
    };
    ViscosityDiluteType type;
    ViscosityDiluteGasCollisionIntegralData collision_integral;  
    ViscosityDiluteCollisionIntegralPowersOfTstarData
      collision_integral_powers_of_Tstar;       
    ViscosityDiluteGasPowersOfT powers_of_T;    
    ViscosityDiluteGasPowersOfTr powers_of_Tr;  
    ViscosityDiluteVariables() {
        type = VISCOSITY_DILUTE_NOT_SET;
    }
};

struct ViscosityRainWaterFriendData
{
    std::vector<CoolPropDbl> b, t;
};
struct ViscosityInitialDensityEmpiricalData
{
    std::vector<CoolPropDbl> n, d, t;
    CoolPropDbl T_reducing, rhomolar_reducing;
};

struct ViscosityInitialDensityVariables
{
    enum ViscosityInitialDensityEnum
    {
        VISCOSITY_INITIAL_DENSITY_RAINWATER_FRIEND,  
        VISCOSITY_INITIAL_DENSITY_EMPIRICAL,         
        VISCOSITY_INITIAL_DENSITY_NOT_SET
    };
    ViscosityInitialDensityEnum type;
    ViscosityRainWaterFriendData rainwater_friend;   
    ViscosityInitialDensityEmpiricalData empirical;  
    ViscosityInitialDensityVariables() {
        type = VISCOSITY_INITIAL_DENSITY_NOT_SET;
    }
};

struct ViscosityModifiedBatschinskiHildebrandData
{
    std::vector<CoolPropDbl> a, d1, d2, t1, t2, f, g, h, p, q, gamma, l;
    CoolPropDbl T_reduce, rhomolar_reduce;
};
struct ViscosityFrictionTheoryData
{
    std::vector<CoolPropDbl> Aa, Aaa, Aaaa, Ar, Arr, Adrdr, Arrr, Ai, Aii, AdrAdr;
    int Na, Naa, Naaa, Nr, Nrr, Nrrr, Nii;
    CoolPropDbl c1, c2, T_reduce, rhomolar_reduce;
};
struct ViscosityHigherOrderVariables
{
    enum ViscosityHigherOrderEnum
    {
        VISCOSITY_HIGHER_ORDER_BATSCHINKI_HILDEBRAND,  
        VISCOSITY_HIGHER_ORDER_HYDROGEN,               
        VISCOSITY_HIGHER_ORDER_HEXANE,                 
        VISCOSITY_HIGHER_ORDER_HEPTANE,                
        VISCOSITY_HIGHER_ORDER_ETHANE,                 
        VISCOSITY_HIGHER_ORDER_BENZENE,                
        VISCOSITY_HIGHER_ORDER_TOLUENE,                
        VISCOSITY_HIGHER_ORDER_CO2_LAESECKE_JPCRD_2017,
        VISCOSITY_HIGHER_ORDER_FRICTION_THEORY,        
        VISCOSITY_HIGHER_ORDER_NOT_SET
    };
    ViscosityHigherOrderEnum type;
    ViscosityModifiedBatschinskiHildebrandData
      modified_Batschinski_Hildebrand;            
    ViscosityFrictionTheoryData friction_theory;  
    ViscosityHigherOrderVariables() {
        type = VISCOSITY_HIGHER_ORDER_NOT_SET;
    };
};

struct ViscosityRhoSrVariables
{
    std::vector<double> c_liq, c_vap;
    double C, x_crossover, rhosr_critical;
};
struct ViscosityECSVariables
{
    std::string reference_fluid;
    CoolPropDbl psi_rhomolar_reducing;
    std::vector<CoolPropDbl> psi_a, psi_t;
};
struct ViscosityChungData
{
    CoolPropDbl rhomolar_critical, acentric, T_critical, molar_mass, dipole_moment_D;
};

class TransportPropertyData
{
   public:
    enum ViscosityHardcodedEnum
    {
        VISCOSITY_HARDCODED_WATER,       
        VISCOSITY_HARDCODED_HEAVYWATER,  
        VISCOSITY_HARDCODED_HELIUM,      
        VISCOSITY_HARDCODED_R23,         
        VISCOSITY_HARDCODED_METHANOL,    
        VISCOSITY_HARDCODED_M_XYLENE,    
        VISCOSITY_HARDCODED_O_XYLENE,    
        VISCOSITY_HARDCODED_P_XYLENE,    
        VISCOSITY_NOT_HARDCODED
    };
    enum ConductivityHardcodedEnum
    {
        CONDUCTIVITY_HARDCODED_WATER,       
        CONDUCTIVITY_HARDCODED_HEAVYWATER,  
        CONDUCTIVITY_HARDCODED_R23,         
        CONDUCTIVITY_HARDCODED_HELIUM,      
        CONDUCTIVITY_HARDCODED_METHANE,     
        CONDUCTIVITY_NOT_HARDCODED
    };
    ViscosityDiluteVariables viscosity_dilute;
    ViscosityInitialDensityVariables viscosity_initial;
    ViscosityHigherOrderVariables viscosity_higher_order;
    ViscosityECSVariables viscosity_ecs;
    ViscosityRhoSrVariables viscosity_rhosr;
    ViscosityChungData viscosity_Chung;

    ConductivityDiluteVariables conductivity_dilute;
    ConductivityResidualVariables conductivity_residual;
    ConductivityCriticalVariables conductivity_critical;
    ConductivityECSVariables conductivity_ecs;

    std::string BibTeX_viscosity,                      
      BibTeX_conductivity;                             
    bool viscosity_using_ECS;                          
    bool conductivity_using_ECS;                       
    bool viscosity_using_Chung;                        
    bool viscosity_using_rhosr;                        
    bool viscosity_model_provided;                     
    bool conductivity_model_provided;                  
    CoolPropDbl sigma_eta,                             
      epsilon_over_k;                                  
    ViscosityHardcodedEnum hardcoded_viscosity;        
    ConductivityHardcodedEnum hardcoded_conductivity;  
    TransportPropertyData()
      : viscosity_using_ECS(false),
        conductivity_using_ECS(false),
        viscosity_using_Chung(false),
        viscosity_using_rhosr(false),
        viscosity_model_provided(false),
        conductivity_model_provided(false),
        sigma_eta(_HUGE),
        epsilon_over_k(_HUGE),
        hardcoded_viscosity(VISCOSITY_NOT_HARDCODED),
        hardcoded_conductivity(CONDUCTIVITY_NOT_HARDCODED) {}
};

struct Ancillaries
{
    SaturationAncillaryFunction pL, pV, rhoL, rhoV, hL, hLV, sL, sLV;
    MeltingLineVariables melting_line;
    SurfaceTensionCorrelation surface_tension;
};


class EquationOfState
{
   public:
    EquationOfState(){};
    ~EquationOfState(){};
    SimpleState reduce,  
      sat_min_liquid,    
      sat_min_vapor,     
      hs_anchor,         
      max_sat_T,         
      max_sat_p;         
    EOSLimits limits;    
    double R_u,          
      molar_mass,        
      acentric,          
      Ttriple,           
      ptriple;           
    bool pseudo_pure;    
    ResidualHelmholtzContainer alphar;  
    IdealHelmholtzContainer alpha0;     
    std::string BibTeX_EOS,             
      BibTeX_CP0;                       
    CriticalRegionSplines
      critical_region_splines;  

    void validate() {
        assert(R_u < 9 && R_u > 8);
        assert(molar_mass > 0.001 && molar_mass < 1);
    };
    CoolPropDbl baser(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.base(tau, delta);
    };
    // First partials
    CoolPropDbl dalphar_dDelta(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.dDelta(tau, delta);
    };
    CoolPropDbl dalphar_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.dTau(tau, delta);
    };
    // Second partials
    CoolPropDbl d2alphar_dDelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.dDelta2(tau, delta);
    };
    CoolPropDbl d2alphar_dDelta_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.dDelta_dTau(tau, delta);
    };
    CoolPropDbl d2alphar_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.dTau2(tau, delta);
    };
    // Third partials
    CoolPropDbl d3alphar_dDelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.dDelta3(tau, delta);
    };
    CoolPropDbl d3alphar_dDelta2_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.dDelta2_dTau(tau, delta);
    };
    CoolPropDbl d3alphar_dDelta_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.dDelta_dTau2(tau, delta);
    };
    CoolPropDbl d3alphar_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alphar.dTau3(tau, delta);
    };

    CoolPropDbl base0(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.base(tau, delta);
    };
    // First partials
    CoolPropDbl dalpha0_dDelta(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.dDelta(tau, delta);
    };
    CoolPropDbl dalpha0_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.dTau(tau, delta);
    };
    // Second partials
    CoolPropDbl d2alpha0_dDelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.dDelta2(tau, delta);
    };
    CoolPropDbl d2alpha0_dDelta_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.dDelta_dTau(tau, delta);
    };
    CoolPropDbl d2alpha0_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.dTau2(tau, delta);
    };
    // Third partials
    CoolPropDbl d3alpha0_dDelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.dDelta3(tau, delta);
    };
    CoolPropDbl d3alpha0_dDelta2_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.dDelta2_dTau(tau, delta);
    };
    CoolPropDbl d3alpha0_dDelta_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.dDelta_dTau2(tau, delta);
    };
    CoolPropDbl d3alpha0_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) {
        return alpha0.dTau3(tau, delta);
    };
};


class CoolPropFluid
{
   protected:
    // Transport property data
    std::string ECSReferenceFluid;  
    double ECS_qd;                  
   public:
    CoolPropFluid() : ECS_qd(-_HUGE) {
        this->ChemSpider_id = -1;
    };
    ~CoolPropFluid(){};
    const EquationOfState& EOS() const {
        return EOSVector[0];
    }  
    EquationOfState& EOS() {
        return EOSVector[0];
    }                                        
    std::vector<EquationOfState> EOSVector;  

    std::string name;  
    std::string
      REFPROPname;  
    std::string CAS;                   
    std::string formula;               
    std::vector<std::string> aliases;  
    std::string InChI;                 
    std::string InChIKey;              
    std::string smiles;                
    int ChemSpider_id;                 
    std::string TwoDPNG_URL;           

    BibTeXKeysStruct BibTeXKeys;             
    EnvironmentalFactorsStruct environment;  
    Ancillaries ancillaries;                 
    TransportPropertyData transport;
    SimpleState crit,  
      triple_liquid,   
      triple_vapor;    

    double gas_constant() {
        return EOS().R_u;
    };
    double molar_mass() {
        return EOS().molar_mass;
    };
};

#if !defined(NO_FMTLIB) && FMT_VERSION >= 90000
inline int format_as(ViscosityDiluteVariables::ViscosityDiluteType type) {
    return fmt::underlying(type);
}

inline int format_as(TransportPropertyData::ViscosityHardcodedEnum viscosity) {
    return fmt::underlying(viscosity);
}

inline int format_as(TransportPropertyData::ConductivityHardcodedEnum conductivity) {
    return fmt::underlying(conductivity);
}
#endif

} /* namespace CoolProp */
#endif /* COOLPROPFLUID_H_ */