CoolPropFluid.h

Go to the documentation of this file.

/*
 * 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 <variant>
#include <optional>
#include <cassert>
#include <iterator>
#include "Eigen/Core"
#include "PolyMath.h"
#include "Ancillaries.h"
#include "superancillary/superancillary.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 = CONDUCTIVITY_DILUTE_NOT_SET;
    ConductivityDiluteRatioPolynomialsData ratio_polynomials;
    ConductivityDiluteEta0AndPolyData eta0_and_poly;
};
 
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 = CONDUCTIVITY_RESIDUAL_NOT_SET;
    ConductivityResidualPolynomialData polynomials;
    ConductivityResidualPolynomialAndExponentialData polynomial_and_exponential;
};
 
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 = CONDUCTIVITY_CRITICAL_NOT_SET;
    ConductivityCriticalSimplifiedOlchowySengersData Olchowy_Sengers;
};
 
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 = VISCOSITY_DILUTE_NOT_SET;
    ViscosityDiluteGasCollisionIntegralData collision_integral;  
    ViscosityDiluteCollisionIntegralPowersOfTstarData
      collision_integral_powers_of_Tstar;       
    ViscosityDiluteGasPowersOfT powers_of_T;    
    ViscosityDiluteGasPowersOfTr powers_of_Tr;  
};
 
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 = VISCOSITY_INITIAL_DENSITY_NOT_SET;
    ViscosityRainWaterFriendData rainwater_friend;   
    ViscosityInitialDensityEmpiricalData empirical;  
};
 
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 = VISCOSITY_HIGHER_ORDER_NOT_SET;
    ViscosityModifiedBatschinskiHildebrandData
      modified_Batschinski_Hildebrand;            
    ViscosityFrictionTheoryData friction_theory;  
};
 
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:
    using SuperAncillary_t = superancillary::SuperAncillary<std::vector<double>>;
    
private:
    std::string superancillaries_str;
    std::optional<SuperAncillary_t> superancillaries = std::nullopt; 
    
   public:
    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;  
    
    std::optional<SuperAncillary_t> & get_superanc_optional(){
        
        if (!superancillaries){
            if (!superancillaries_str.empty()){
                auto start = std::chrono::high_resolution_clock::now(); // Start time
                // Now do the parsing pass and replace with the actual superancillary
                superancillaries.emplace(SuperAncillary_t(superancillaries_str));
                auto end = std::chrono::high_resolution_clock::now();  // End time
                auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
                //std::cout << "Execution time: " << duration.count() << " microseconds for " << BibTeX_EOS << std::endl;
            }
        }
        return superancillaries;
    }
    
    void set_superancillaries_str(const std::string &s){
        superancillaries_str = s;
        // Do the construction greedily by default, but allow it to be lazy if you want
#if !defined(LAZY_LOAD_SUPERANCILLARIES)
        get_superanc_optional();
#endif
    }
    
    
 
    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_ */