diff --git a/include/IncompressibleFluid.h b/include/IncompressibleFluid.h index 61d29e05..5e9f94c9 100644 --- a/include/IncompressibleFluid.h +++ b/include/IncompressibleFluid.h @@ -60,7 +60,6 @@ protected: composition_types xid; double TminPsat; - double Tref, pref, xref, href, sref; double xbase, Tbase; /// These are the objects that hold the coefficients @@ -144,11 +143,6 @@ public: double getxmin() const {return xmin;} composition_types getxid() const {return xid;} double getTminPsat() const {return TminPsat;} - double getTref() const {return Tref;} - double getpref() const {return pref;} - double getxref() const {return xref;} - double gethref() const {return href;} - double getsref() {return sref;} double getTbase() const {return Tbase;} double getxbase() const {return xbase;} @@ -161,10 +155,6 @@ public: void setxmin(double xmin) {this->xmin = xmin;} void setxid(composition_types xid) {this->xid = xid;} void setTminPsat(double TminPsat) {this->TminPsat = TminPsat;} - //void setTref(double Tref) {this->Tref = Tref;} - //void setpref(double pref) {this->pref = pref;} - //void setxref(double xref) {this->xref = xref;} - void set_reference_state(double T0, double p0, double x0, double h0, double s0); void setTbase(double Tbase) {this->Tbase = Tbase;} void setxbase(double xbase) {this->xbase = xbase;} @@ -203,14 +193,14 @@ public: double rho (double T, double p, double x); /// Heat capacities as a function of temperature, pressure and composition. double c (double T, double p, double x); - double cp (double T, double p, double x){return c(T,p,x);}; - double cv (double T, double p, double x){return c(T,p,x);}; + double cp (double T, double p, double x){throw ValueError(format("%s (%d): Please use the c-function instead.",__FILE__,__LINE__));} + double cv (double T, double p, double x){throw ValueError(format("%s (%d): Please use the c-function instead.",__FILE__,__LINE__));} /// Entropy as a function of temperature, pressure and composition. - double s (double T, double p, double x); + double s (double T, double p, double x){throw ValueError(format("%s (%d): The internal calculations have changed, use the backend to calculate entropy from the partial derivatives.",__FILE__,__LINE__));} /// Internal energy as a function of temperature, pressure and composition. - double u (double T, double p, double x); + double u (double T, double p, double x){throw ValueError(format("%s (%d): The internal calculations have changed, use the backend to calculate internal energy from enthalpy.",__FILE__,__LINE__));} /// Enthalpy as a function of temperature, pressure and composition. - double h (double T, double p, double x); + double h (double T, double p, double x){throw ValueError(format("%s (%d): The internal calculations have changed, use the backend to calculate enthalpy from the partial derivatives.",__FILE__,__LINE__));} /// Viscosity as a function of temperature, pressure and composition. double visc(double T, double p, double x); /// Thermal conductivity as a function of temperature, pressure and composition. @@ -220,29 +210,27 @@ public: /// Freezing temperature as a function of pressure and composition. double Tfreeze( double p, double x); - /* Below are direct calculations of the derivatives. Nothing - * special is going on, we simply use the polynomial class to - * derive the different functions with respect to temperature. - */ - /// Partial derivative of density with respect to temperature at constant pressure and composition - double drhodTatPx (double T, double p, double x); - /// Partial derivative of entropy with respect to temperature at constant pressure and composition - double dsdTatPx (double T, double p, double x){return c(T,p,x)/T;}; - /// Partial derivative of internal energy with respect to temperature at constant pressure and composition - double dudTatPx (double T, double p, double x){return c(T,p,x);}; - /// Partial derivative of enthalpy with respect to temperature at constant pressure and composition - double dhdTatPx (double T, double p, double x){return c(T,p,x);}; - - - /* Other useful derivatives - */ - /// Partial derivative of enthalpy with respect to pressure at constant temperature and composition - // \f[ \left( \frac{\partial h}{\partial p} \right)_T = v - T \left( \frac{\partial v}{\partial T} \right)_p = \rho^{-1} \left( 1 + T \rho^{-1} \left( \frac{\partial \rho}{\partial T} \right)_p \right) \f] - double dhdpatTx (double T, double p, double x); - /// Partial derivative of entropy with respect to pressure at constant temperature and composition - // \f[ \left( \frac{\partial s}{\partial p} \right)_T = - \left( \frac{\partial v}{\partial T} \right)_p = \rho^{-2} \left( \frac{\partial \rho}{\partial T} \right)_p \right) \f] - double dsdpatTx (double T, double p, double x); + * special is going on, we simply use the polynomial class to + * derive the different functions with respect to temperature. + */ + /// Partial derivative of density + // with respect to temperature at constant pressure and composition + double drhodTatPx(double T, double p, double x); + ///// Partial derivative of entropy + //// with respect to temperature at constant pressure and composition + //double dsdTatPx (double T, double p, double x){return c(T,p,x)/T;}; + ///// Partial derivative of enthalpy + //// with respect to temperature at constant pressure and composition + //double dhdTatPx (double T, double p, double x){return c(T,p,x);}; + /// Partial derivative of entropy + // with respect to temperature at constant pressure and composition + // integrated in temperature + double dsdTatPxdT(double T, double p, double x); + /// Partial derivative of enthalpy + // with respect to temperature at constant pressure and composition + // integrated in temperature + double dhdTatPxdT(double T, double p, double x); /// Mass fraction conversion function diff --git a/src/Backends/Incompressible/IncompressibleBackend.cpp b/src/Backends/Incompressible/IncompressibleBackend.cpp index 0531eaeb..7742741a 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.cpp +++ b/src/Backends/Incompressible/IncompressibleBackend.cpp @@ -20,27 +20,29 @@ namespace CoolProp { IncompressibleBackend::IncompressibleBackend() { - //this->_fractions_id = ifrac_undefined; + throw NotImplementedError("Empty constructor is not implemented for incompressible fluids"); + this->fluid = NULL; } IncompressibleBackend::IncompressibleBackend(IncompressibleFluid* fluid) { - //this->_fractions_id = fluid->getxid(); this->fluid = fluid; + this->set_reference_state(); if (this->fluid->is_pure()){ - this->set_fractions(std::vector(1,1.0)); + this->set_fractions(std::vector(1,0.0)); } } IncompressibleBackend::IncompressibleBackend(const std::string &fluid_name) { this->fluid = &get_incompressible_fluid(fluid_name); - //this->_fractions_id = this->fluid->getxid(); + this->set_reference_state(); if (this->fluid->is_pure()){ - this->set_fractions(std::vector(1,1.0)); + this->set_fractions(std::vector(1,0.0)); } } IncompressibleBackend::IncompressibleBackend(const std::vector &component_names) { throw NotImplementedError("Mixture-style constructor is not implemented yet for incompressible fluids"); + this->fluid = NULL; } void IncompressibleBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) { @@ -125,6 +127,32 @@ void IncompressibleBackend::update(CoolProp::input_pairs input_pair, double valu fluid->checkTPX(_T,_p,_fractions[0]); } +/// Clear all the cached values +bool IncompressibleBackend::clear() { + AbstractState::clear(); // Call the base class + /// Reference values, no need to calculate them each time + this->_hmass_ref.clear(); + this->_smass_ref.clear(); + /// Additional cached elements used for the partial derivatives + this->_cmass.clear(); + this->_hmass.clear(); + this->_rhomass.clear(); + this->_smass.clear(); + this->_umass.clear(); + // Done + return true; +} + +/// Update the reference values and clear the state +void IncompressibleBackend::set_reference_state(double T0, double p0, double x0, double h0, double s0){ + this->clear(); + this->_T_ref = T0; + this->_p_ref = p0; + this->_x_ref = x0; + this->_h_ref = h0; + this->_s_ref = s0; +} + /// Set the fractions /** @param fractions The vector of fractions of the components converted to the correct input @@ -136,7 +164,7 @@ void IncompressibleBackend::set_fractions(const std::vector &fracti ( this->_fractions[0]!=fractions[0] ) ) { // Change it! if (get_debug_level()>=20) std::cout << format("Incompressible backend: Updating the fractions triggered a change in reference state %s -> %s",vec_to_string(this->_fractions).c_str(),vec_to_string(fractions).c_str()) << std::endl; this->_fractions = fractions; - fluid->set_reference_state(fluid->getTref(), fluid->getpref(), this->_fractions[0], fluid->gethref(), fluid->getsref()); + set_reference_state(T_ref(), p_ref(), this->_fractions[0], h_ref(), s_ref()); } } @@ -208,6 +236,74 @@ void IncompressibleBackend::check_status() { throw NotImplementedError("Cannot check status for incompressible fluid"); } +/** We have to override some of the functions from the AbstractState. + * The incompressibles are only mass-based and do not support conversion + * from molar to specific quantities. + * We also have a few new chaced variables that we need. + */ +/// Return the mass density in kg/m^3 +double IncompressibleBackend::rhomass(void){ + if (!_rhomass) _rhomass = calc_rhomass(); + return _rhomass; +} +/// Return the mass enthalpy in J/kg +double IncompressibleBackend::hmass(void){ + if (!_hmass) _hmass = calc_hmass(); + return _hmass; +} +/// Return the molar entropy in J/mol/K +double IncompressibleBackend::smass(void){ + if (!_smass) _smass = calc_smass(); + return _smass; +} +/// Return the molar internal energy in J/mol +double IncompressibleBackend::umass(void){ + if (!_umass) _umass = calc_umass(); + return _umass; +} +/// Return the mass constant pressure specific heat in J/kg/K +double IncompressibleBackend::cmass(void){ + if (!_cmass) _cmass = calc_cmass(); + return _cmass; +} + +/// Return the temperature in K +double IncompressibleBackend::T_ref(void){ + if (!_T_ref) throw ValueError("Reference temperature is not set"); + return _T_ref; +} +/// Return the pressure in Pa +double IncompressibleBackend::p_ref(void){ + if (!_p_ref) throw ValueError("Reference pressure is not set"); + return _p_ref; +} +/// Return the composition +double IncompressibleBackend::x_ref(void){ + if (!_x_ref) throw ValueError("Reference composition is not set"); + return _x_ref; +} +/// Return the mass enthalpy in J/kg +double IncompressibleBackend::h_ref(void){ + if (!_h_ref) throw ValueError("Reference enthalpy is not set"); + return _h_ref; +} +/// Return the molar entropy in J/mol/K +double IncompressibleBackend::s_ref(void){ + if (!_s_ref) throw ValueError("Reference entropy is not set"); + return _s_ref; +} + +/// Return the mass enthalpy in J/kg +double IncompressibleBackend::hmass_ref(void){ + if (!_hmass_ref) _hmass_ref = raw_calc_hmass(T_ref(),p_ref(),x_ref()); + return _hmass_ref; +} +/// Return the molar entropy in J/mol/K +double IncompressibleBackend::smass_ref(void){ + if (!_smass_ref) _smass_ref = raw_calc_smass(T_ref(),p_ref(),x_ref()); + return _smass_ref; +} + /// Calculate T given pressure and density /** @param rhomass The mass density in kg/m^3 @@ -228,26 +324,24 @@ long double IncompressibleBackend::HmassP_flash(long double hmass, long double p class HmassP_residual : public FuncWrapper1D { protected: double p,x,h_in; - IncompressibleFluid* fluid; + IncompressibleBackend* backend; protected: HmassP_residual(){}; public: - HmassP_residual(IncompressibleFluid* fluid, const double &p, const double &x, const double &h_in){ - this->p = p; - this->x = x; + HmassP_residual(IncompressibleBackend* backend, const double &p, const double &x, const double &h_in){ + this->p=p; + this->x=x; this->h_in = h_in; - this->fluid = fluid; + this->backend=backend; } virtual ~HmassP_residual(){}; double call(double target){ - return fluid->h(target,p,x) - h_in; //fluid.u(target,p,x)+ p / fluid.rho(target,p,x) - h_in; + return backend->raw_calc_hmass(target,p,x) - h_in; //fluid.u(target,p,x)+ p / fluid.rho(target,p,x) - h_in; } //double deriv(double target); }; - //double T_tmp = this->PUmass_flash(p, hmass); // guess value from u=h - - HmassP_residual res = HmassP_residual(fluid, p, _fractions[0], hmass); + HmassP_residual res = HmassP_residual(this, p, _fractions[0], hmass-h_ref()+hmass_ref()); std::string errstring; double macheps = DBL_EPSILON; @@ -268,23 +362,23 @@ long double IncompressibleBackend::PSmass_flash(long double p, long double smass class PSmass_residual : public FuncWrapper1D { protected: double p,x,s_in; - IncompressibleFluid* fluid; + IncompressibleBackend* backend; protected: PSmass_residual(){}; public: - PSmass_residual(IncompressibleFluid* fluid, const double &p, const double &x, const double &s_in){ + PSmass_residual(IncompressibleBackend* backend, const double &p, const double &x, const double &s_in){ this->p = p; this->x = x; this->s_in = s_in; - this->fluid = fluid; + this->backend = backend; } virtual ~PSmass_residual(){}; double call(double target){ - return fluid->s(target,p,x) - s_in; + return backend->raw_calc_smass(target,p,x) - s_in; } }; - PSmass_residual res = PSmass_residual(fluid, p, _fractions[0], smass); + PSmass_residual res = PSmass_residual(this, p, _fractions[0], smass-s_ref()+smass_ref()); std::string errstring; double macheps = DBL_EPSILON; @@ -295,42 +389,84 @@ long double IncompressibleBackend::PSmass_flash(long double p, long double smass return result; } -/// Calculate T given pressure and internal energy -/** -@param umass The mass internal energy in J/kg -@param p The pressure in Pa -@returns T The temperature in K -*/ -long double IncompressibleBackend::PUmass_flash(long double p, long double umass){ +///// Calculate T given pressure and internal energy +///** +//@param umass The mass internal energy in J/kg +//@param p The pressure in Pa +//@returns T The temperature in K +//*/ +//long double IncompressibleBackend::PUmass_flash(long double p, long double umass){ +// +// class PUmass_residual : public FuncWrapper1D { +// protected: +// double p,x,u_in; +// IncompressibleBackend* fluid; +// protected: +// PUmass_residual(){}; +// public: +// PUmass_residual(IncompressibleBackend* fluid, const double &p, const double &x, const double &u_in){ +// this->p = p; +// this->x = x; +// this->u_in = u_in; +// this->fluid = fluid; +// } +// virtual ~PUmass_residual(){}; +// double call(double target){ +// return fluid->u(target,p,x) - u_in; +// } +// }; +// +// PUmass_residual res = PUmass_residual(*this, p, _fractions[0], umass); +// +// std::string errstring; +// double macheps = DBL_EPSILON; +// double tol = DBL_EPSILON*1e3; +// int maxiter = 10; +// double result = Brent(&res, fluid->getTmin(), fluid->getTmax(), macheps, tol, maxiter, errstring); +// //if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; +// return result; +//} - class PUmass_residual : public FuncWrapper1D { - protected: - double p,x,u_in; - IncompressibleFluid* fluid; - protected: - PUmass_residual(){}; - public: - PUmass_residual(IncompressibleFluid* fluid, const double &p, const double &x, const double &u_in){ - this->p = p; - this->x = x; - this->u_in = u_in; - this->fluid = fluid; - } - virtual ~PUmass_residual(){}; - double call(double target){ - return fluid->u(target,p,x) - u_in; - } - }; +/// We start with the functions that do not need a reference state +long double IncompressibleBackend::calc_melting_line(int param, int given, long double value){ + if (param == iT && given == iP){ + return fluid->Tfreeze(value, _fractions[0]); + } else { + throw ValueError("For incompressibles, the only valid inputs to calc_melting_line are T(p)"); + } +}; +long double IncompressibleBackend::calc_umass(void){ + return hmass()-_p/rhomass(); +}; - PUmass_residual res = PUmass_residual(fluid, p, _fractions[0], umass); +/// ... and continue with the ones that depend on reference conditions. +long double IncompressibleBackend::calc_hmass(void){ + return h_ref() + raw_calc_hmass(_T,_p,_fractions[0]) - hmass_ref(); +}; +long double IncompressibleBackend::calc_smass(void){ + return s_ref() + raw_calc_smass(_T,_p,_fractions[0]) - smass_ref(); +}; - std::string errstring; - double macheps = DBL_EPSILON; - double tol = DBL_EPSILON*1e3; - int maxiter = 10; - double result = Brent(&res, fluid->getTmin(), fluid->getTmax(), macheps, tol, maxiter, errstring); - //if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; - return result; + +/// Functions that can be used with the solver, they miss the reference values! +long double IncompressibleBackend::raw_calc_hmass(double T, double p, double x){ + return calc_dhdTatPxdT(T,p,x) + p * calc_dhdpatTx(T,fluid->rho(T, p, x),calc_drhodTatPx(T,p,x)); +}; +long double IncompressibleBackend::raw_calc_smass(double T, double p, double x){ + return calc_dsdTatPxdT(T,p,x) + p * calc_dsdpatTx( fluid->rho(T, p, x),calc_drhodTatPx(T,p,x)); +}; + +/* Other useful derivatives + */ +/// Partial derivative of entropy with respect to pressure at constant temperature and composition +// \f[ \left( \frac{\partial s}{\partial p} \right)_T = - \left( \frac{\partial v}{\partial T} \right)_p = \rho^{-2} \left( \frac{\partial \rho}{\partial T} \right)_p \right) \f] +double IncompressibleBackend::calc_dsdpatTx (double rho, double drhodTatPx){ + return 1/rho/rho * drhodTatPx; +} +/// Partial derivative of enthalpy with respect to pressure at constant temperature and composition +// \f[ \left( \frac{\partial h}{\partial p} \right)_T = v - T \left( \frac{\partial v}{\partial T} \right)_p = \rho^{-1} \left( 1 + T \rho^{-1} \left( \frac{\partial \rho}{\partial T} \right)_p \right) \f] +double IncompressibleBackend::calc_dhdpatTx (double T, double rho, double drhodTatPx){ + return 1/rho * ( 1 + T/rho * drhodTatPx); } @@ -395,47 +531,35 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CHECK( check_abs(T,res,acc) ); } - // Compare h - val = fluid.h(T, p, x); - res = backend.HmassP_flash(val, p); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(val); - CAPTURE(res); - CHECK( check_abs(T,res,acc) ); - } - - // Compare s - val = fluid.s(T, p, x); - res = backend.PSmass_flash(p, val); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(val); - CAPTURE(res); - CHECK( check_abs(T,res,acc) ); - } - - // Compare u - val = fluid.u(T, p, x); - res = backend.PUmass_flash(p, val); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(val); - CAPTURE(res); - CHECK( check_abs(T,res,acc) ); - } - // call the update function to set internal variables, // concentration has been set before. CHECK_THROWS( backend.update( CoolProp::DmassT_INPUTS, val, T ) ); // First with wrong parameters CHECK_NOTHROW( backend.update( CoolProp::PT_INPUTS, p, T ) ); // ... and then with the correct ones. + // Compare enthalpy flash + val = backend.hmass(); + res = backend.HmassP_flash(val, p); + { + CAPTURE(T); + CAPTURE(p); + CAPTURE(x); + CAPTURE(val); + CAPTURE(res); + CHECK( check_abs(T,res,acc) ); + } + + // Compare entropy flash + val = backend.smass(); + res = backend.PSmass_flash(p, val); + { + CAPTURE(T); + CAPTURE(p); + CAPTURE(x); + CAPTURE(val); + CAPTURE(res); + CHECK( check_abs(T,res,acc) ); + } + /// Get the viscosity [Pa-s] val = fluid.visc(T, p, x); res = backend.calc_viscosity(); @@ -461,7 +585,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi } val = fluid.rho(T, p, x); - res = backend.calc_rhomass(); + res = backend.rhomass(); { CAPTURE(T); CAPTURE(p); @@ -471,42 +595,31 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CHECK( check_abs(val,res,acc) ); } - val = fluid.h(T, p, x); - res = backend.calc_hmass(); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(val); - CAPTURE(res); - CHECK( check_abs(val,res,acc) ); - } - - val = fluid.s(T, p, x); - res = backend.calc_smass(); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(val); - CAPTURE(res); - CHECK( check_abs(val,res,acc) ); - } +// val = fluid.h(T, p, x); +// res = backend.hmass(); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( check_abs(val,res,acc) ); +// } +// +// val = fluid.s(T, p, x); +// res = backend.smass(); +// { +// CAPTURE(T); +// CAPTURE(p); +// CAPTURE(x); +// CAPTURE(val); +// CAPTURE(res); +// CHECK( check_abs(val,res,acc) ); +// } // Make sure the result does not change -> reference state... - val = backend.calc_smass(); + val = backend.smass(); backend.update( CoolProp::PT_INPUTS, p, T ); - res = backend.calc_smass(); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(val); - CAPTURE(res); - CHECK( check_abs(val,res,acc) ); - } - - val = fluid.u(T, p, x); - res = backend.calc_umass(); + res = backend.smass(); { CAPTURE(T); CAPTURE(p); @@ -517,7 +630,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi } val = fluid.c(T, p, x); - res = backend.calc_cpmass(); + res = backend.cpmass(); { CAPTURE(T); CAPTURE(p); @@ -527,7 +640,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CHECK( check_abs(val,res,acc) ); } - res = backend.calc_cvmass(); + res = backend.cvmass(); { CAPTURE(T); CAPTURE(p); @@ -537,10 +650,9 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CHECK( check_abs(val,res,acc) ); } - // Compare Tfreeze val = fluid.Tfreeze(p, x);//-20.02+273.15;// 253.1293105454671; - res = -20.02+273.15; + res = backend.calc_T_freeze();// -20.02+273.15; { CAPTURE(T); CAPTURE(p); diff --git a/src/Backends/Incompressible/IncompressibleBackend.h b/src/Backends/Incompressible/IncompressibleBackend.h index e8db03d5..e65e61bb 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.h +++ b/src/Backends/Incompressible/IncompressibleBackend.h @@ -12,11 +12,21 @@ namespace CoolProp { class IncompressibleBackend : public AbstractState { + protected: - //int Ncomp; - //static bool _REFPROP_supported; - //int _fractions_id; - std::vector _fractions; + + + /// Bulk values, state variables + //double _T, _p; // From AbstractState + std::vector _fractions; + + /// Reference values, no need to calculate them each time + CachedElement _T_ref,_p_ref,_x_ref,_h_ref,_s_ref; + CachedElement _hmass_ref, _smass_ref; + + /// Additional cached elements used for the partial derivatives + CachedElement _cmass, _hmass, _rhomass, _smass, _umass; + IncompressibleFluid *fluid; /// Set the fractions @@ -55,6 +65,12 @@ public: */ void update(CoolProp::input_pairs input_pair, double value1, double value2); + /// Clear all the cached values + bool clear(); + + /// Update the reference values and clear the state + void set_reference_state(double T0=20+273.15, double p0=101325, double x0=0.0, double h0=0.0, double s0=0.0); + /// Set the mole fractions /** @param mole_fractions The vector of mole fractions of the components @@ -77,6 +93,42 @@ public: /// Check if the mole fractions have been set, etc. void check_status(); + /** We have to override some of the functions from the AbstractState. + * The incompressibles are only mass-based and do not support conversion + * from molar to specific quantities. + * We also have a few new chaced variables that we need. + */ + /// Return the mass density in kg/m^3 + double rhomass(void); + /// Return the mass enthalpy in J/kg + double hmass(void); + /// Return the molar entropy in J/mol/K + double smass(void); + /// Return the molar internal energy in J/mol + double umass(void); + /// Return the mass constant pressure specific heat in J/kg/K + double cmass(void); + + /// Return the temperature in K + double T_ref(void); + /// Return the pressure in Pa + double p_ref(void); + /// Return the composition + double x_ref(void); + /// Return the mass enthalpy in J/kg + double h_ref(void); + /// Return the molar entropy in J/mol/K + double s_ref(void); + + /// Return the mass enthalpy in J/kg + double hmass_ref(void); + /// Return the molar entropy in J/mol/K + double smass_ref(void); + + + /** These functions should be protected, but that requires new tests. + * I'll leave that as a TODO item for now. + */ /// Calculate T given pressure and density /** @param rhomass The mass density in kg/m^3 @@ -99,42 +151,72 @@ public: */ long double PSmass_flash(long double p, long double smass); - /// Calculate T given pressure and internal energy - /** - @param umass The mass internal energy in J/kg - @param p The pressure in Pa - @returns T The temperature in K - */ - long double PUmass_flash(long double p, long double umass); - - /// Get the viscosity [Pa-s] - long double calc_viscosity(void){return fluid->visc(_T, _p, _fractions[0]);}; - /// Get the thermal conductivity [W/m/K] (based on the temperature and pressure in the state class) - long double calc_conductivity(void){return fluid->cond(_T, _p, _fractions[0]);}; +// /// Calculate T given pressure and internal energy +// /** +// @param umass The mass internal energy in J/kg +// @param p The pressure in Pa +// @returns T The temperature in K +// */ +// long double PUmass_flash(long double p, long double umass); + /// We start with the functions that do not need a reference state long double calc_rhomass(void){return fluid->rho(_T, _p, _fractions[0]);}; - long double calc_hmass(void){return fluid->h(_T, _p, _fractions[0]);}; - long double calc_smass(void){return fluid->s(_T, _p, _fractions[0]);}; - long double calc_umass(void){return fluid->u(_T, _p, _fractions[0]);}; - long double calc_cpmass(void){return fluid->cp(_T, _p, _fractions[0]);}; - long double calc_cvmass(void){return fluid->cv(_T, _p, _fractions[0]);}; + long double calc_cmass(void){return fluid->c(_T, _p, _fractions[0]);}; + long double calc_cpmass(void){return cmass();}; + long double calc_cvmass(void){return cmass();}; + long double calc_viscosity(void){return fluid->visc(_T, _p, _fractions[0]);}; + long double calc_conductivity(void){return fluid->cond(_T, _p, _fractions[0]);}; + long double calc_T_freeze(void){return fluid->Tfreeze(_p, _fractions[0]);}; + long double calc_melting_line(int param, int given, long double value); + long double calc_umass(void); + + /// ... and continue with the ones that depend on reference conditions. + long double calc_hmass(void); + long double calc_smass(void); + +public: + /// Functions that can be used with the solver, they miss the reference values! + long double raw_calc_hmass(double T, double p, double x); + long double raw_calc_smass(double T, double p, double x); + + +protected: + /* Below are direct calculations of the derivatives. Nothing + * special is going on, we simply use the polynomial class to + * derive the different functions with respect to temperature. + */ + /// Partial derivative of density with respect to temperature at constant pressure and composition + double calc_drhodTatPx(double T, double p, double x){return fluid->drhodTatPx(T,p,x);}; + /// Partial derivative of entropy with respect to temperature at constant pressure and composition + double calc_dsdTatPx (double T, double p, double x){return fluid->c(T,p,x)/T;}; + /// Partial derivative of enthalpy with respect to temperature at constant pressure and composition + double calc_dhdTatPx (double T, double p, double x){return fluid->c(T,p,x);}; + /// Partial derivative of entropy + // with respect to temperature at constant pressure and composition + // integrated in temperature + double calc_dsdTatPxdT(double T, double p, double x){return fluid->dsdTatPxdT(T,p,x);}; + /// Partial derivative of enthalpy + // with respect to temperature at constant pressure and composition + // integrated in temperature + double calc_dhdTatPxdT(double T, double p, double x){return fluid->dhdTatPxdT(T,p,x);}; + + + /* Other useful derivatives + */ + /// Partial derivative of entropy with respect to pressure at constant temperature and composition + // \f[ \left( \frac{\partial s}{\partial p} \right)_T = - \left( \frac{\partial v}{\partial T} \right)_p = \rho^{-2} \left( \frac{\partial \rho}{\partial T} \right)_p \right) \f] + double calc_dsdpatTx (double rho, double drhodTatPx); + /// Partial derivative of enthalpy with respect to pressure at constant temperature and composition + // \f[ \left( \frac{\partial h}{\partial p} \right)_T = v - T \left( \frac{\partial v}{\partial T} \right)_p = \rho^{-1} \left( 1 + T \rho^{-1} \left( \frac{\partial \rho}{\partial T} \right)_p \right) \f] + double calc_dhdpatTx (double T, double rho, double drhodTatPx); + + +public: + /// Constants from the fluid object long double calc_Tmax(void){return fluid->getTmax();}; long double calc_Tmin(void){return fluid->getTmin();}; - long double calc_fraction_min(void){return fluid->getxmin();}; long double calc_fraction_max(void){return fluid->getxmax();}; - long double calc_T_freeze(void){ - return fluid->Tfreeze(_p, _fractions[0]);}; - - long double calc_melting_line(int param, int given, long double value){ - if (param == iT && given == iP){ - return fluid->Tfreeze(value, _fractions[0]); - } - else{ - throw ValueError("For incompressibles, the only valid inputs to calc_melting_line are T(p)"); - } - }; - std::string calc_name(void){return fluid->getDescription();}; }; diff --git a/src/Backends/Incompressible/IncompressibleFluid.cpp b/src/Backends/Incompressible/IncompressibleFluid.cpp index 9d5b5b8d..2672b90f 100644 --- a/src/Backends/Incompressible/IncompressibleFluid.cpp +++ b/src/Backends/Incompressible/IncompressibleFluid.cpp @@ -18,14 +18,6 @@ and transport properties. */ //IncompressibleFluid::IncompressibleFluid(); -void IncompressibleFluid::set_reference_state(double T0, double p0, double x0, double h0, double s0){ - this->Tref = T0; - this->pref = p0; - this->xref = x0; - this->href = h0; - this->sref = s0; -} - void IncompressibleFluid::validate(){ return; // TODO: Implement validation function @@ -122,52 +114,6 @@ double IncompressibleFluid::c (double T, double p, double x){ return _HUGE; } -/// Entropy as a function of temperature, pressure and composition. -double IncompressibleFluid::s (double T, double p, double x){ - double dsdTatp = _HUGE; - double dsdpatT = _HUGE; - switch (specific_heat.type) { - case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: - //TODO: Optimise this, add a cached value or direct calculation of definite integral - dsdTatp = poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase) - poly.integral(specific_heat.coeffs, Tref, xref, 0, -1, 0, Tbase, xbase); - dsdpatT = this->rho(T,p,x); - dsdpatT = this->drhodTatPx(T,p,x)/dsdpatT/dsdpatT * (p-pref); - return sref + dsdTatp + dsdpatT; - break; - case IncompressibleData::INCOMPRESSIBLE_NOT_SET: - throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); - break; - default: - throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for entropy.",__FILE__,__LINE__,specific_heat.type)); - break; - } - return _HUGE; -} - -/// Internal energy as a function of temperature, pressure and composition. -double IncompressibleFluid::u (double T, double p, double x){return u_h(T,p,x);}; - -/// Enthalpy as a function of temperature, pressure and composition. -double IncompressibleFluid::h (double T, double p, double x){ - double dhdTatP = _HUGE; - double dhdpatT = _HUGE; - switch (specific_heat.type) { - case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: - //TODO: Optimise this, add a cached value or direct calculation of definite integral - dhdTatP = poly.integral(specific_heat.coeffs, T, x, 0, 0, 0, Tbase, xbase) - poly.integral(specific_heat.coeffs, Tref, xref, 0, 0, 0, Tbase, xbase); - dhdpatT = this->dhdpatTx(T, p, x) * (p-pref); - return href + dhdTatP + dhdpatT; - break; - case IncompressibleData::INCOMPRESSIBLE_NOT_SET: - throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); - break; - default: - throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for internal energy.",__FILE__,__LINE__,specific_heat.type)); - break; - } - return _HUGE; - } - /// Viscosity as a function of temperature, pressure and composition. double IncompressibleFluid::visc(double T, double p, double x){ switch (viscosity.type) { @@ -298,21 +244,42 @@ double IncompressibleFluid::drhodTatPx (double T, double p, double x){ } return _HUGE; } +/// Partial derivative of entropy +// with respect to temperature at constant pressure and composition +// integrated in temperature +double IncompressibleFluid::dsdTatPxdT(double T, double p, double x){ + switch (specific_heat.type) { + case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: + return poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase); + break; + case IncompressibleData::INCOMPRESSIBLE_NOT_SET: + throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); + break; + default: + throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for entropy.",__FILE__,__LINE__,specific_heat.type)); + break; + } + return _HUGE; +} +/// Partial derivative of enthalpy +// with respect to temperature at constant pressure and composition +// integrated in temperature +double IncompressibleFluid::dhdTatPxdT(double T, double p, double x){ + switch (specific_heat.type) { + case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: + return poly.integral(specific_heat.coeffs, T, x, 0, 0, 0, Tbase, xbase); + break; + case IncompressibleData::INCOMPRESSIBLE_NOT_SET: + throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); + break; + default: + throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for entropy.",__FILE__,__LINE__,specific_heat.type)); + break; + } + return _HUGE; +} + -/* Other useful derivatives - */ -/// Partial derivative of enthalpy with respect to temperature at constant pressure and composition -double IncompressibleFluid::dhdpatTx (double T, double p, double x){ - double rho = 0; - rho = this->rho(T,p,x); - return 1/rho * ( 1 + T/rho * this->drhodTatPx(T,p,x)); -} -/// Partial derivative of entropy with respect to temperature at constant pressure and composition -double IncompressibleFluid::dsdpatTx (double T, double p, double x){ - double rho = 0; - rho = this->rho(T,p,x); - return 1/rho/rho * this->drhodTatPx(T,p,x); -} /// Mass fraction conversion function @@ -647,66 +614,16 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi //XLT.setVolToMass(parse_coefficients(fluid_json, "volume2mass", false)); //XLT.setMassToMole(parse_coefficients(fluid_json, "mass2mole", false)); - //XLT.set_reference_state(25+273.15, 1.01325e5, 0.0, 0.0, 0.0); - double Tref = 25+273.15; - double pref = 2e5; - double xref = 0.0; - double href = 127.0; - double sref = 23.0; - XLT.set_reference_state(Tref, pref, xref, href, sref); - /// A function to check coefficients and equation types. //XLT.validate(); double acc = 0.0001; double val = 0; double res = 0; - // Compare reference state - { - res = XLT.h(Tref,pref,xref); - CAPTURE(Tref); - CAPTURE(pref); - CAPTURE(xref); - CAPTURE(href); - CAPTURE(res); - CHECK( check_abs(href,res,acc) ); - } - { - res = XLT.s(Tref,pref,xref); - CAPTURE(Tref); - CAPTURE(pref); - CAPTURE(xref); - CAPTURE(sref); - CAPTURE(res); - CHECK( check_abs(sref,res,acc) ); - double res2 = XLT.s(Tref,pref,xref); - CAPTURE(Tref); - CAPTURE(pref); - CAPTURE(xref); - CAPTURE(res); - CAPTURE(res2); - CHECK( check_abs(res2,res,acc) ); - res = XLT.s(Tref,pref,xref); - CAPTURE(Tref); - CAPTURE(pref); - CAPTURE(xref); - CAPTURE(res); - CAPTURE(res2); - CHECK( check_abs(res2,res,acc) ); - } - - - - Tref = 25+273.15; - pref = 0.0; - xref = 0.0; - href = 0.0; - sref = 0.0; - XLT.set_reference_state(Tref, pref, xref, href, sref); // Prepare the results and compare them to the calculated values double T = 273.15+50; double p = 10e5; - double x = xref; + double x = 0.0; // Compare density val = 824.4615702148608; @@ -728,62 +645,10 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CHECK( check_abs(val,res,acc) ); } - // Compare s - val = 144.08; - res = XLT.s(T,p,x); - { - CAPTURE(T); - CAPTURE(val); - CAPTURE(res); - CHECK( check_abs(val,res,acc) ); - } - - val = 0.0; - res = XLT.s(Tref,pref,xref); - { - CAPTURE(T); - CAPTURE(val); - CAPTURE(res); - CHECK( val==res ); - } - - // Compare u - val = 44724.1; - res = XLT.u(T,p,x); - { - CAPTURE(T); - CAPTURE(val); - CAPTURE(res); - CHECK( check_abs(val,res,acc) ); - } - - val = href - pref/XLT.rho(Tref,pref,xref); - res = XLT.u(Tref,pref,xref); - { - CAPTURE(T); - CAPTURE(val); - CAPTURE(res); - CHECK( val==res ); - } - - // Compare h - val = 45937; - res = XLT.h(T,p,x); - { - CAPTURE(T); - CAPTURE(val); - CAPTURE(res); - CHECK( check_abs(val,res,acc) ); - } - - val = 0.0; - res = XLT.h(Tref,pref,xref); - { - CAPTURE(T); - CAPTURE(val); - CAPTURE(res); - CHECK( val==res ); - } + // Check property functions + CHECK_THROWS(XLT.s(T,p,x)); + CHECK_THROWS(XLT.h(T,p,x)); + CHECK_THROWS(XLT.u(T,p,x)); // Compare v val = 0.0008931435169681835; @@ -811,17 +676,6 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CoolProp::IncompressibleFluid CH3OH = CoolPropTesting::incompressibleFluidObject(); - //XLT.set_reference_state(25+273.15, 1.01325e5, 0.0, 0.0, 0.0); - double Tref = 25+273.15; - double pref = 0.0; - double xref = 0.25; - double href = 0.0; - double sref = 0.0; - CH3OH.set_reference_state(Tref, pref, xref, href, sref); - - /// A function to check coefficients and equation types. - //CH3OH.validate(); - // Prepare the results and compare them to the calculated values double acc = 0.0001; double T = 273.15+10; @@ -854,94 +708,10 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CHECK( check_abs(expected,actual,acc) ); } - // Compare s - expected = -207.027; - actual = CH3OH.s(T,p,x); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(expected); - CAPTURE(actual); - CHECK( check_abs(expected,actual,acc) ); - } - expected = CH3OH.s(T,p,x); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(expected); - CAPTURE(actual); - CHECK( check_abs(expected,actual,acc) ); - } - - expected = 0.0; - actual = CH3OH.s(Tref,pref,xref); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(expected); - CAPTURE(actual); - CHECK( expected==actual ); - } - expected = CH3OH.s(Tref,pref,xref); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(expected); - CAPTURE(actual); - CHECK( expected==actual ); - } - - // Compare u - expected = -60157.1; - actual = CH3OH.u(T,p,x); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(expected); - CAPTURE(actual); - CHECK( check_abs(expected,actual,acc) ); - } - - expected = href - pref/CH3OH.rho(Tref,pref,xref); - actual = CH3OH.u(Tref,pref,xref); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(expected); - CAPTURE(actual); - CHECK( expected==actual ); - } - - // Compare h - expected = -59119; - actual = CH3OH.h(T,p,x); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(expected); - CAPTURE(actual); - CHECK( check_abs(expected,actual,acc) ); - } - - expected = 0.0; - actual = CH3OH.h(Tref,pref,xref); - { - CAPTURE(T); - CAPTURE(p); - CAPTURE(x); - CAPTURE(expected); - CAPTURE(actual); - std::string errmsg = CoolProp::get_global_param_string("errstring"); - CAPTURE(errmsg); - CHECK( expected==actual ); - } + // Check property functions + CHECK_THROWS(CH3OH.s(T,p,x)); + CHECK_THROWS(CH3OH.h(T,p,x)); + CHECK_THROWS(CH3OH.u(T,p,x)); // Compare v expected = 0.0023970245009602097; diff --git a/src/Backends/Incompressible/IncompressibleLibrary.cpp b/src/Backends/Incompressible/IncompressibleLibrary.cpp index 8fc1a6a0..9cd90cc2 100644 --- a/src/Backends/Incompressible/IncompressibleLibrary.cpp +++ b/src/Backends/Incompressible/IncompressibleLibrary.cpp @@ -7,326 +7,321 @@ namespace CoolProp{ -/// Class to access Lithium-Bromide solutions -/** Employs some basic wrapper-like functionality - * to bridge the gap between the solution functions - * used in the paper by Pátek and Klomfar: - * http://dx.doi.org/10.1016/j.ijrefrig.2005.10.007 - * - * We owe gratitude to the authors for providing - * both access to the paper as well as the equations - * in the form of C source code. */ - -double const LiBrSolution::M_H2O = 0.018015268; /* kg/mol, molar mass of H2O */ -double const LiBrSolution::M_LiBr = 0.08685; /* kg/mol, molar mass of LiBr */ -double const LiBrSolution::T0 = 221; /* K, constant */ - -/* Critical point of H2O */ -double const LiBrSolution::Tc_H2O = 647.096; /* K, temperature */ -double const LiBrSolution::pc_H2O = 22.064; /* MPa, pressure */ -double const LiBrSolution::rhoc_H2O = 17873; /* mol/m^3 (322 kg/m^3), molar density */ -double const LiBrSolution::hc_H2O = 37548.5; /* J/mol, molar enthalpy */ -double const LiBrSolution::sc_H2O = 79.3933; /* J/(mol.K) molar entropy*/ - -/*Triple point of H2O */ -double const LiBrSolution::Tt_H2O = 273.16; /* K, temperature */ -double const LiBrSolution::cpt_H2O = 76.0226; /* J/(mol.K), molar isobaric heat capacity*/ - -double LiBrSolution::ps_mix(double T, double x) -/* Equation (1) */ -{ - static double m[8] = { 3.0, 4.0, 4.0, 8.0, 1.0, 1.0, 4.0, 6.0 }; - static double n[8] = { 0.0, 5.0, 6.0, 3.0, 0.0, 2.0, 6.0, 0.0 }; - static double t[8] = { 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 }; - static double a[8] = { -2.41303e2, 1.91750e7, -1.75521e8, 3.25430e7, - 3.92571e2, -2.12626e3, 1.85127e8, 1.91216e3 }; - double tau, suma = 0.0; - int i; - - tau = T / Tc_H2O; - for (i = 0; i <= 7; i++) - suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]); - return (ps_H2O(T - suma)); - -} /* end function ps_mix */ - -double LiBrSolution::rho_mix(double T, double x) -/* Equation (2) */ -{ - static double m[2] = { 1.0, 1.0 }; - static double n[2] = { 0.0, 6.0 }; - static double a[2] = { 1.746, 4.709 }; - - double tau, suma = 0.0; - int i; - - tau = T / Tc_H2O; - for (i = 0; i <= 1; i++) - suma += a[i] * pow(x, m[i]) * pow(tau, n[i]); - - return ((1.0 - x) * rho_H2O(T) + rhoc_H2O * suma); - -} /* end function rho_mix */ - -double LiBrSolution::cp_mix(double T, double x) -/* Equation (3) */ -{ - static double m[8] = { 2.0, 3.0, 3.0, 3.0, 3.0, 2.0, 1.0, 1.0 }; - static double n[8] = { 0.0, 0.0, 1.0, 2.0, 3.0, 0.0, 3.0, 2.0 }; - static double t[8] = { 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 3.0, 4.0 }; - static double a[8] = { -1.42094e1, 4.04943e1, 1.11135e2, 2.29980e2, - 1.34526e3, -1.41010e-2, 1.24977e-2, -6.83209e-4 }; - - double tau, suma = 0.0; - int i; - - tau = Tc_H2O / (T - T0); - for (i = 0; i <= 7; i++) - suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]); - - return ((1.0 - x) * cp_H2O(T) + cpt_H2O * suma); - -} /* end function cp_mix */ - -double LiBrSolution::h_mix(double T, double x) -/* Equation (4) */ -{ - static double m[30] = { 1.0, 1.0, 2.0, 3.0, 6.0, 1.0, 3.0, 5.0, 4.0, - 5.0, 5.0, 6.0, 6.0, 1.0, 2.0, 2.0, 2.0, 5.0, 6.0, 7.0, 1.0, 1.0, - 2.0, 2.0, 2.0, 3.0, 1.0, 1.0, 1.0, 1.0 }; - - static double n[30] = { 0.0, 1.0, 6.0, 6.0, 2.0, 0.0, 0.0, 4.0, 0.0, - 4.0, 5.0, 5.0, 6.0, 0.0, 3.0, 5.0, 7.0, 0.0, 3.0, 1.0, 0.0, 4.0, - 2.0, 6.0, 7.0, 0.0, 0.0, 1.0, 2.0, 3.0 }; - - static double t[30] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, - 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0, - 4.0, 4.0, 4.0, 4.0, 5.0, 5.0, 5.0, 5.0 }; - - static double a[30] = { 2.27431, -7.99511, 3.85239e2, -1.63940e4, - -4.22562e2, 1.13314e-1, -8.33474, -1.73833e4, 6.49763, - 3.24552e3, -1.34643e4, 3.99322e4, -2.58877e5, -1.93046e-3, - 2.80616, -4.04479e1, 1.45342e2, -2.74873, -4.49743e2, - -1.21794e1, -5.83739e-3, 2.33910e-1, 3.41888e-1, 8.85259, - -1.78731e1, 7.35179e-2, -1.79430e-4, 1.84261e-3, -6.24282e-3, - 6.84765e-3 }; - - double tau, suma = 0.0; - int i; - - tau = Tc_H2O / (T - T0); - for (i = 0; i <= 29; i++) - suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]); - - return ((1.0 - x) * h_H2O(T) + hc_H2O * suma); - -} /* end function h_mix */ - -double LiBrSolution::s_mix(double T, double x) -/* Equation (5) */ -{ - static double m[29] = { 1.0, 1.0, 2.0, 3.0, 6.0, 1.0, 3.0, 5.0, 1.0, - 2.0, 2.0, 4.0, 5.0, 5.0, 6.0, 6.0, 1.0, 3.0, 5.0, 7.0, 1.0, 1.0, - 1.0, 2.0, 3.0, 1.0, 1.0, 1.0, 1.0 }; - - static double n[29] = { 0.0, 1.0, 6.0, 6.0, 2.0, 0.0, 0.0, 4.0, 0.0, - 0.0, 4.0, 0.0, 4.0, 5.0, 2.0, 5.0, 0.0, 4.0, 0.0, 1.0, 0.0, 2.0, - 4.0, 7.0, 1.0, 0.0, 1.0, 2.0, 3.0 }; - - static double t[29] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, - 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0, - 4.0, 4.0, 4.0, 5.0, 5.0, 5.0, 5.0 }; - - static double a[29] = { 1.53091, -4.52564, 6.98302e+2, -2.16664e+4, - -1.47533e+3, 8.47012e-2, -6.59523, -2.95331e+4, 9.56314e-3, - -1.88679e-1, 9.31752, 5.78104, 1.38931e+4, -1.71762e+4, - 4.15108e+2, -5.55647e+4, -4.23409e-3, 3.05242e+1, -1.67620, - 1.48283e+1, 3.03055e-3, -4.01810e-2, 1.49252e-1, 2.59240, - -1.77421e-1, -6.99650e-5, 6.05007e-4, -1.65228e-3, 1.22966e-3 }; - - double tau, suma = 0.0; - int i; - - tau = Tc_H2O / (T - T0); - for (i = 0; i <= 28; i++) - suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]); - - return ((1.0 - x) * s_H2O(T) + sc_H2O * suma); - -} /* end function s_mix */ - -double LiBrSolution::ps_H2O(double T) -/* Equation (28) */ -{ - static double a[7] = { 0.0, -7.85951783, 1.84408259, -11.7866497, - 22.6807411, -15.9618719, 1.80122502 }; - - double tau, ps; - - tau = 1 - T / Tc_H2O; - - ps = pc_H2O - * exp( - Tc_H2O / T - * (a[1] * tau + a[2] * pow(tau, 1.5) - + a[3] * pow(tau, 3.0) - + a[4] * pow(tau, 3.5) - + a[5] * pow(tau, 4.0) - + a[6] * pow(tau, 7.5))); - - return (ps * 1.0e6); - -} /* end function ps_H2O */ - -double LiBrSolution::rho_H2O(double T) -/* Equation (29) */ -{ - static double b[7] = { 0.0, 1.99274064, 1.09965342, -0.510839303, - -1.75493479, -45.5170352, -6.7469445e5 }; - double theta, rho; - - theta = 1.0 - T / Tc_H2O; - - rho = rhoc_H2O - * (1.0 + b[1] * pow(theta, 1.0 / 3.0) - + b[2] * pow(theta, 2.0 / 3.0) - + b[3] * pow(theta, 5.0 / 3.0) - + b[4] * pow(theta, 16.0 / 3.0) - + b[5] * pow(theta, 43.0 / 3.0) - + b[6] * pow(theta, 110.0 / 3.0)); - - return (rho); -} /* end function rho_H2O */ - -double LiBrSolution::cp_H2O(double T) -/* Equation (30) */ -{ - static double a[5] = - { 1.38801, -2.95318, 3.18721, -0.645473, 9.18946e5 }; - static double b[5] = { 0.0, 2.0, 3.0, 6.0, 34.0 }; - static double c[5] = { 0.0, 2.0, 3.0, 5.0, 0.0 }; - - double suma = 0.0; - int i; - - for (i = 0; i <= 4; i++) - suma += a[i] * exp(b[i] * log(1.0 - T / Tc_H2O)) - * exp(c[i] * log(T / Tt_H2O)); - - return (cpt_H2O * suma); - -} /* end function cp_H2O */ - -double LiBrSolution::h_H2O(double T) -/* Equation (31) */ -{ - static double a[4] = { -4.37196e-1, 3.03440e-1, -1.29582, -1.76410e-1 }; - static double alpha[4] = { 1.0 / 3.0, 2.0 / 3.0, 5.0 / 6.0, 21.0 / 6.0 }; - - double suma = 0.0; - int i; - - for (i = 0; i <= 3; i++) - suma += a[i] * exp(alpha[i] * log(1.0 - T / Tc_H2O)); - - return (hc_H2O * (1.0 + suma)); - -} /* end function h_H2O */ - -double LiBrSolution::s_H2O(double T) -/* Equation (32) */ -{ - static double a[4] = { -3.34112e-1, -8.47987e-1, -9.11980e-1, -1.64046 }; - static double alpha[4] = { 1.0 / 3.0, 3.0 / 3.0, 8.0 / 3.0, 24.0 / 3.0 }; - - double suma = 0.0; - int i; - - for (i = 0; i <= 3; i++) - suma += a[i] * exp(alpha[i] * log(1.0 - T / Tc_H2O)); - - return (sc_H2O * (1.0 + suma)); - -} /* end function s_H2O */ - - -/** Finished with the code from the paper. Now we need to - * convert the molar values to mass-based units. */ -double LiBrSolution::massToMole(double w) -/* Equation (7) */ -{ - return (w/M_LiBr)/(w/M_LiBr+(1.-w)/M_H2O); - //return (w*M_LiBr)/(w*M_LiBr+(1.-w)*M_H2O); -} - -double LiBrSolution::molarToSpecific(double w, double value) -/* Equation (7,8) */ -{ - double x = massToMole(w); - //return w/(x*M_LiBr) * value; - return 1. / ( x*M_LiBr + (1.-x)*M_H2O ) * value; -} - -bool const LiBrSolution::debug = false; - - - -LiBrSolution::LiBrSolution():IncompressibleFluid(){ - name = std::string("LiBr"); - description = std::string("Lithium-Bromide solution from Patek2006"); - reference = std::string("Patek2006"); - - Tmin = 273.00; - Tmax = 500.00; - TminPsat = Tmin; - - xmin = 0.0; - xmax = 1.0; - - xref = 0.0; - Tref = 0.0; - pref = 0.0; - href = 0.0; - sref = 0.0; - xbase = 0.0; - Tbase = 0.0; - -}; - -double LiBrSolution::rho(double T, double p, double x){ - checkTPX(T, p, x); - return 1./molarToSpecific(x, 1./rho_mix(T,massToMole(x))); -} -double LiBrSolution::c(double T, double p, double x){ - checkTPX(T, p, x); - return molarToSpecific(x, cp_mix(T,massToMole(x))); -} -//double h(double T, double p, double x){ -// return h_u(T,p,x); +///// Class to access Lithium-Bromide solutions +///** Employs some basic wrapper-like functionality +// * to bridge the gap between the solution functions +// * used in the paper by Pátek and Klomfar: +// * http://dx.doi.org/10.1016/j.ijrefrig.2005.10.007 +// * +// * We owe gratitude to the authors for providing +// * both access to the paper as well as the equations +// * in the form of C source code. */ +// +//double const LiBrSolution::M_H2O = 0.018015268; /* kg/mol, molar mass of H2O */ +//double const LiBrSolution::M_LiBr = 0.08685; /* kg/mol, molar mass of LiBr */ +//double const LiBrSolution::T0 = 221; /* K, constant */ +// +///* Critical point of H2O */ +//double const LiBrSolution::Tc_H2O = 647.096; /* K, temperature */ +//double const LiBrSolution::pc_H2O = 22.064; /* MPa, pressure */ +//double const LiBrSolution::rhoc_H2O = 17873; /* mol/m^3 (322 kg/m^3), molar density */ +//double const LiBrSolution::hc_H2O = 37548.5; /* J/mol, molar enthalpy */ +//double const LiBrSolution::sc_H2O = 79.3933; /* J/(mol.K) molar entropy*/ +// +///*Triple point of H2O */ +//double const LiBrSolution::Tt_H2O = 273.16; /* K, temperature */ +//double const LiBrSolution::cpt_H2O = 76.0226; /* J/(mol.K), molar isobaric heat capacity*/ +// +//double LiBrSolution::ps_mix(double T, double x) +///* Equation (1) */ +//{ +// static double m[8] = { 3.0, 4.0, 4.0, 8.0, 1.0, 1.0, 4.0, 6.0 }; +// static double n[8] = { 0.0, 5.0, 6.0, 3.0, 0.0, 2.0, 6.0, 0.0 }; +// static double t[8] = { 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 }; +// static double a[8] = { -2.41303e2, 1.91750e7, -1.75521e8, 3.25430e7, +// 3.92571e2, -2.12626e3, 1.85127e8, 1.91216e3 }; +// double tau, suma = 0.0; +// int i; +// +// tau = T / Tc_H2O; +// for (i = 0; i <= 7; i++) +// suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]); +// return (ps_H2O(T - suma)); +// +//} /* end function ps_mix */ +// +//double LiBrSolution::rho_mix(double T, double x) +///* Equation (2) */ +//{ +// static double m[2] = { 1.0, 1.0 }; +// static double n[2] = { 0.0, 6.0 }; +// static double a[2] = { 1.746, 4.709 }; +// +// double tau, suma = 0.0; +// int i; +// +// tau = T / Tc_H2O; +// for (i = 0; i <= 1; i++) +// suma += a[i] * pow(x, m[i]) * pow(tau, n[i]); +// +// return ((1.0 - x) * rho_H2O(T) + rhoc_H2O * suma); +// +//} /* end function rho_mix */ +// +//double LiBrSolution::cp_mix(double T, double x) +///* Equation (3) */ +//{ +// static double m[8] = { 2.0, 3.0, 3.0, 3.0, 3.0, 2.0, 1.0, 1.0 }; +// static double n[8] = { 0.0, 0.0, 1.0, 2.0, 3.0, 0.0, 3.0, 2.0 }; +// static double t[8] = { 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 3.0, 4.0 }; +// static double a[8] = { -1.42094e1, 4.04943e1, 1.11135e2, 2.29980e2, +// 1.34526e3, -1.41010e-2, 1.24977e-2, -6.83209e-4 }; +// +// double tau, suma = 0.0; +// int i; +// +// tau = Tc_H2O / (T - T0); +// for (i = 0; i <= 7; i++) +// suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]); +// +// return ((1.0 - x) * cp_H2O(T) + cpt_H2O * suma); +// +//} /* end function cp_mix */ +// +//double LiBrSolution::h_mix(double T, double x) +///* Equation (4) */ +//{ +// static double m[30] = { 1.0, 1.0, 2.0, 3.0, 6.0, 1.0, 3.0, 5.0, 4.0, +// 5.0, 5.0, 6.0, 6.0, 1.0, 2.0, 2.0, 2.0, 5.0, 6.0, 7.0, 1.0, 1.0, +// 2.0, 2.0, 2.0, 3.0, 1.0, 1.0, 1.0, 1.0 }; +// +// static double n[30] = { 0.0, 1.0, 6.0, 6.0, 2.0, 0.0, 0.0, 4.0, 0.0, +// 4.0, 5.0, 5.0, 6.0, 0.0, 3.0, 5.0, 7.0, 0.0, 3.0, 1.0, 0.0, 4.0, +// 2.0, 6.0, 7.0, 0.0, 0.0, 1.0, 2.0, 3.0 }; +// +// static double t[30] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, +// 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0, +// 4.0, 4.0, 4.0, 4.0, 5.0, 5.0, 5.0, 5.0 }; +// +// static double a[30] = { 2.27431, -7.99511, 3.85239e2, -1.63940e4, +// -4.22562e2, 1.13314e-1, -8.33474, -1.73833e4, 6.49763, +// 3.24552e3, -1.34643e4, 3.99322e4, -2.58877e5, -1.93046e-3, +// 2.80616, -4.04479e1, 1.45342e2, -2.74873, -4.49743e2, +// -1.21794e1, -5.83739e-3, 2.33910e-1, 3.41888e-1, 8.85259, +// -1.78731e1, 7.35179e-2, -1.79430e-4, 1.84261e-3, -6.24282e-3, +// 6.84765e-3 }; +// +// double tau, suma = 0.0; +// int i; +// +// tau = Tc_H2O / (T - T0); +// for (i = 0; i <= 29; i++) +// suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]); +// +// return ((1.0 - x) * h_H2O(T) + hc_H2O * suma); +// +//} /* end function h_mix */ +// +//double LiBrSolution::s_mix(double T, double x) +///* Equation (5) */ +//{ +// static double m[29] = { 1.0, 1.0, 2.0, 3.0, 6.0, 1.0, 3.0, 5.0, 1.0, +// 2.0, 2.0, 4.0, 5.0, 5.0, 6.0, 6.0, 1.0, 3.0, 5.0, 7.0, 1.0, 1.0, +// 1.0, 2.0, 3.0, 1.0, 1.0, 1.0, 1.0 }; +// +// static double n[29] = { 0.0, 1.0, 6.0, 6.0, 2.0, 0.0, 0.0, 4.0, 0.0, +// 0.0, 4.0, 0.0, 4.0, 5.0, 2.0, 5.0, 0.0, 4.0, 0.0, 1.0, 0.0, 2.0, +// 4.0, 7.0, 1.0, 0.0, 1.0, 2.0, 3.0 }; +// +// static double t[29] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, +// 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0, +// 4.0, 4.0, 4.0, 5.0, 5.0, 5.0, 5.0 }; +// +// static double a[29] = { 1.53091, -4.52564, 6.98302e+2, -2.16664e+4, +// -1.47533e+3, 8.47012e-2, -6.59523, -2.95331e+4, 9.56314e-3, +// -1.88679e-1, 9.31752, 5.78104, 1.38931e+4, -1.71762e+4, +// 4.15108e+2, -5.55647e+4, -4.23409e-3, 3.05242e+1, -1.67620, +// 1.48283e+1, 3.03055e-3, -4.01810e-2, 1.49252e-1, 2.59240, +// -1.77421e-1, -6.99650e-5, 6.05007e-4, -1.65228e-3, 1.22966e-3 }; +// +// double tau, suma = 0.0; +// int i; +// +// tau = Tc_H2O / (T - T0); +// for (i = 0; i <= 28; i++) +// suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]); +// +// return ((1.0 - x) * s_H2O(T) + sc_H2O * suma); +// +//} /* end function s_mix */ +// +//double LiBrSolution::ps_H2O(double T) +///* Equation (28) */ +//{ +// static double a[7] = { 0.0, -7.85951783, 1.84408259, -11.7866497, +// 22.6807411, -15.9618719, 1.80122502 }; +// +// double tau, ps; +// +// tau = 1 - T / Tc_H2O; +// +// ps = pc_H2O +// * exp( +// Tc_H2O / T +// * (a[1] * tau + a[2] * pow(tau, 1.5) +// + a[3] * pow(tau, 3.0) +// + a[4] * pow(tau, 3.5) +// + a[5] * pow(tau, 4.0) +// + a[6] * pow(tau, 7.5))); +// +// return (ps * 1.0e6); +// +//} /* end function ps_H2O */ +// +//double LiBrSolution::rho_H2O(double T) +///* Equation (29) */ +//{ +// static double b[7] = { 0.0, 1.99274064, 1.09965342, -0.510839303, +// -1.75493479, -45.5170352, -6.7469445e5 }; +// double theta, rho; +// +// theta = 1.0 - T / Tc_H2O; +// +// rho = rhoc_H2O +// * (1.0 + b[1] * pow(theta, 1.0 / 3.0) +// + b[2] * pow(theta, 2.0 / 3.0) +// + b[3] * pow(theta, 5.0 / 3.0) +// + b[4] * pow(theta, 16.0 / 3.0) +// + b[5] * pow(theta, 43.0 / 3.0) +// + b[6] * pow(theta, 110.0 / 3.0)); +// +// return (rho); +//} /* end function rho_H2O */ +// +//double LiBrSolution::cp_H2O(double T) +///* Equation (30) */ +//{ +// static double a[5] = +// { 1.38801, -2.95318, 3.18721, -0.645473, 9.18946e5 }; +// static double b[5] = { 0.0, 2.0, 3.0, 6.0, 34.0 }; +// static double c[5] = { 0.0, 2.0, 3.0, 5.0, 0.0 }; +// +// double suma = 0.0; +// int i; +// +// for (i = 0; i <= 4; i++) +// suma += a[i] * exp(b[i] * log(1.0 - T / Tc_H2O)) +// * exp(c[i] * log(T / Tt_H2O)); +// +// return (cpt_H2O * suma); +// +//} /* end function cp_H2O */ +// +//double LiBrSolution::h_H2O(double T) +///* Equation (31) */ +//{ +// static double a[4] = { -4.37196e-1, 3.03440e-1, -1.29582, -1.76410e-1 }; +// static double alpha[4] = { 1.0 / 3.0, 2.0 / 3.0, 5.0 / 6.0, 21.0 / 6.0 }; +// +// double suma = 0.0; +// int i; +// +// for (i = 0; i <= 3; i++) +// suma += a[i] * exp(alpha[i] * log(1.0 - T / Tc_H2O)); +// +// return (hc_H2O * (1.0 + suma)); +// +//} /* end function h_H2O */ +// +//double LiBrSolution::s_H2O(double T) +///* Equation (32) */ +//{ +// static double a[4] = { -3.34112e-1, -8.47987e-1, -9.11980e-1, -1.64046 }; +// static double alpha[4] = { 1.0 / 3.0, 3.0 / 3.0, 8.0 / 3.0, 24.0 / 3.0 }; +// +// double suma = 0.0; +// int i; +// +// for (i = 0; i <= 3; i++) +// suma += a[i] * exp(alpha[i] * log(1.0 - T / Tc_H2O)); +// +// return (sc_H2O * (1.0 + suma)); +// +//} /* end function s_H2O */ +// +// +///** Finished with the code from the paper. Now we need to +// * convert the molar values to mass-based units. */ +//double LiBrSolution::massToMole(double w) +///* Equation (7) */ +//{ +// return (w/M_LiBr)/(w/M_LiBr+(1.-w)/M_H2O); +// //return (w*M_LiBr)/(w*M_LiBr+(1.-w)*M_H2O); +//} +// +//double LiBrSolution::molarToSpecific(double w, double value) +///* Equation (7,8) */ +//{ +// double x = massToMole(w); +// //return w/(x*M_LiBr) * value; +// return 1. / ( x*M_LiBr + (1.-x)*M_H2O ) * value; +//} +// +//bool const LiBrSolution::debug = false; +// +// +// +//LiBrSolution::LiBrSolution():IncompressibleFluid(){ +// name = std::string("LiBr"); +// description = std::string("Lithium-Bromide solution from Patek2006"); +// reference = std::string("Patek2006"); +// +// Tmin = 273.00; +// Tmax = 500.00; +// TminPsat = Tmin; +// +// xmin = 0.0; +// xmax = 1.0; +// +// xbase = 0.0; +// Tbase = 0.0; +// +//}; +// +//double LiBrSolution::rho(double T, double p, double x){ +// checkTPX(T, p, x); +// return 1./molarToSpecific(x, 1./rho_mix(T,massToMole(x))); +//} +//double LiBrSolution::c(double T, double p, double x){ +// checkTPX(T, p, x); +// return molarToSpecific(x, cp_mix(T,massToMole(x))); +//} +////double h(double T, double p, double x){ +//// return h_u(T,p,x); +////} +//double LiBrSolution::s(double T, double p, double x){ +// checkTPX(T, p, x); +// return molarToSpecific(x, s_mix(T,massToMole(x))); +//} +//double LiBrSolution::visc(double T, double p, double x){ +// throw ValueError("Viscosity is not defined for LiBr-solutions."); +//} +//double LiBrSolution::cond(double T, double p, double x){ +// throw ValueError("Thermal conductivity is not defined for LiBr-solutions."); +//} +//double LiBrSolution::u(double T, double p, double x){ +// checkTPX(T, p, x); +// return molarToSpecific(x, h_mix(T,massToMole(x))); +//} +//double LiBrSolution::psat(double T, double x){ +// //checkT(T,p,x); +// if (debug) throw ValueError(format("Your concentration is %f in kg/kg and %f in mol/mol.",x,massToMole(x))); +// return ps_mix(T,massToMole(x)); +//}; +//double LiBrSolution::Tfreeze(double p, double x){ +// if (debug) throw ValueError(format("No freezing point data available for Lithium-Bromide: p=%f, x=%f",p,x)); +// return Tmin; //} -double LiBrSolution::s(double T, double p, double x){ - checkTPX(T, p, x); - return molarToSpecific(x, s_mix(T,massToMole(x))); -} -double LiBrSolution::visc(double T, double p, double x){ - throw ValueError("Viscosity is not defined for LiBr-solutions."); -} -double LiBrSolution::cond(double T, double p, double x){ - throw ValueError("Thermal conductivity is not defined for LiBr-solutions."); -} -double LiBrSolution::u(double T, double p, double x){ - checkTPX(T, p, x); - return molarToSpecific(x, h_mix(T,massToMole(x))); -} -double LiBrSolution::psat(double T, double x){ - //checkT(T,p,x); - if (debug) throw ValueError(format("Your concentration is %f in kg/kg and %f in mol/mol.",x,massToMole(x))); - return ps_mix(T,massToMole(x)); -}; -double LiBrSolution::Tfreeze(double p, double x){ - if (debug) throw ValueError(format("No freezing point data available for Lithium-Bromide: p=%f, x=%f",p,x)); - return Tmin; -} /// Default constructor @@ -482,14 +477,14 @@ void JSONIncompressibleLibrary::add_one(rapidjson::Value &fluid_json) { fluid.setVolume2input(parse_coefficients(fluid_json, "volume2input", false)); fluid.setMole2input(parse_coefficients(fluid_json, "mole2input", false)); - if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading reference state for %s ",fluid.getName().c_str()) << std::endl; - fluid.set_reference_state( - parse_value(fluid_json, "Tref", false, 20+273.15) , - parse_value(fluid_json, "pref", false, 1.01325e5) , - parse_value(fluid_json, "xref", false, 0.0) , - parse_value(fluid_json, "href", false, 0.0) , - parse_value(fluid_json, "sref", false, 0.0) - ); + //if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading reference state for %s ",fluid.getName().c_str()) << std::endl; + //fluid.set_reference_state( + // parse_value(fluid_json, "Tref", false, 20+273.15) , + // parse_value(fluid_json, "pref", false, 1.01325e5) , + // parse_value(fluid_json, "xref", false, 0.0) , + // parse_value(fluid_json, "href", false, 0.0) , + // parse_value(fluid_json, "sref", false, 0.0) + // ); /// A function to check coefficients and equation types. fluid.validate(); diff --git a/src/Backends/Incompressible/IncompressibleLibrary.h b/src/Backends/Incompressible/IncompressibleLibrary.h index 89cfcfb6..205de8eb 100644 --- a/src/Backends/Incompressible/IncompressibleLibrary.h +++ b/src/Backends/Incompressible/IncompressibleLibrary.h @@ -17,116 +17,116 @@ namespace CoolProp{ extern int get_debug_level(); -/// Class to access Lithium-Bromide solutions -/** Employs some basic wrapper-like functionality - * to bridge the gap between the solution functions - * used in the paper by Pátek and Klomfar: - * http://dx.doi.org/10.1016/j.ijrefrig.2005.10.007 - * - * We owe gratitude to the authors for providing - * both access to the paper as well as the equations - * in the form of C source code. */ -class LiBrSolution : public IncompressibleFluid{ - -protected: - static double const M_H2O; /* kg/mol, molar mass of H2O */ - static double const M_LiBr; /* kg/mol, molar mass of LiBr */ - static double const T0; /* K, constant */ - - /* Critical point of H2O */ - static double const Tc_H2O; /* K, temperature */ - static double const pc_H2O; /* MPa, pressure */ - static double const rhoc_H2O; /* mol/m^3 (322 kg/m^3), molar density */ - static double const hc_H2O; /* J/mol, molar enthalpy */ - static double const sc_H2O; /* J/(mol.K) molar entropy*/ - - /*Triple point of H2O */ - static double const Tt_H2O; /* K, temperature */ - static double const cpt_H2O; /* J/(mol.K), molar isobaric heat capacity*/ - - double ps_mix(double T, double x); - double rho_mix(double T, double x); - double cp_mix(double T, double x); - double h_mix(double T, double x); - double s_mix(double T, double x); - double ps_H2O(double T); - double rho_H2O(double T); - double cp_H2O(double T); - double h_H2O(double T); - double s_H2O(double T); - - /** Finished with the code from the paper. Now we need to - * convert the molar values to mass-based units. */ - double massToMole(double w); - double molarToSpecific(double w, double value); - - static const bool debug; - -public: - - LiBrSolution(); - - double rho(double T, double p, double x); - double c(double T, double p, double x); - //double h(double T_K, double p, double x); - double s(double T, double p, double x); - double visc(double T, double p, double x); - double cond(double T, double p, double x); - double u(double T, double p, double x); - double psat(double T, double x); - double Tfreeze(double p, double x); - - /* Some functions can be inverted directly, those are listed - * here. It is also possible to solve for other quantities, but - * that involves some more sophisticated processing and is not - * done here, but in the backend, T(h,p) for example. - */ - /// Temperature as a function of density, pressure and composition. - double T_rho (double Dmass, double p, double x){throw NotImplementedError(format("%s (%d): T from density is not implemented for LiBr.",__FILE__,__LINE__));} - /// Temperature as a function of heat capacities as a function of temperature, pressure and composition. - double T_c (double Cmass, double p, double x){throw NotImplementedError(format("%s (%d): T from heat capacity is not implemented for LiBr.",__FILE__,__LINE__));} - /// Temperature as a function of entropy as a function of temperature, pressure and composition. - double T_s (double Smass, double p, double x){throw NotImplementedError(format("%s (%d): T from entropy is not implemented for LiBr.",__FILE__,__LINE__));} - /// Temperature as a function of internal energy as a function of temperature, pressure and composition. - double T_u (double Umass, double p, double x){throw NotImplementedError(format("%s (%d): T from internal energy is not implemented for LiBr.",__FILE__,__LINE__));} - /// Temperature as a function of enthalpy, pressure and composition. - //double T_h (double Hmass, double p, double x){throw NotImplementedError(format("%s (%d): T from enthalpy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));} - /// Viscosity as a function of temperature, pressure and composition. - double T_visc(double visc, double p, double x){throw NotImplementedError(format("%s (%d): T from viscosity is not implemented.",__FILE__,__LINE__));} - /// Thermal conductivity as a function of temperature, pressure and composition. - double T_cond(double cond, double p, double x){throw NotImplementedError(format("%s (%d): T from conductivity is not implemented.",__FILE__,__LINE__));} - /// Saturation pressure as a function of temperature and composition. - double T_psat(double psat, double x){throw NotImplementedError(format("%s (%d): T from psat is not implemented.",__FILE__,__LINE__));} - /// Composition as a function of freezing temperature and pressure. - double x_Tfreeze( double Tfreeze, double p){throw NotImplementedError(format("%s (%d): x from T_freeze is not implemented.",__FILE__,__LINE__));} - - - /// Overwrite some standard functions that cannot be used with LiBr - void setName(std::string name){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setDescription(std::string description){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setReference(std::string reference){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setTmax(double Tmax){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setTmin(double Tmin){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setxmax(double xmax){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setxmin(double xmin){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setTminPsat(double TminPsat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - - void setTbase(double Tbase){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setxbase(double xbase){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - - void setDensity(IncompressibleData density){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setSpecificHeat(IncompressibleData specific_heat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setViscosity(IncompressibleData viscosity){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setConductivity(IncompressibleData conductivity){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setPsat(IncompressibleData p_sat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setTfreeze(IncompressibleData T_freeze){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setVolToMass(IncompressibleData volToMass){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - void setMassToMole(IncompressibleData massToMole){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} - - bool is_pure() {return false;}; - -}; - +///// Class to access Lithium-Bromide solutions +///** Employs some basic wrapper-like functionality +// * to bridge the gap between the solution functions +// * used in the paper by Pátek and Klomfar: +// * http://dx.doi.org/10.1016/j.ijrefrig.2005.10.007 +// * +// * We owe gratitude to the authors for providing +// * both access to the paper as well as the equations +// * in the form of C source code. */ +//class LiBrSolution : public IncompressibleFluid{ +// +//protected: +// static double const M_H2O; /* kg/mol, molar mass of H2O */ +// static double const M_LiBr; /* kg/mol, molar mass of LiBr */ +// static double const T0; /* K, constant */ +// +// /* Critical point of H2O */ +// static double const Tc_H2O; /* K, temperature */ +// static double const pc_H2O; /* MPa, pressure */ +// static double const rhoc_H2O; /* mol/m^3 (322 kg/m^3), molar density */ +// static double const hc_H2O; /* J/mol, molar enthalpy */ +// static double const sc_H2O; /* J/(mol.K) molar entropy*/ +// +// /*Triple point of H2O */ +// static double const Tt_H2O; /* K, temperature */ +// static double const cpt_H2O; /* J/(mol.K), molar isobaric heat capacity*/ +// +// double ps_mix(double T, double x); +// double rho_mix(double T, double x); +// double cp_mix(double T, double x); +// double h_mix(double T, double x); +// double s_mix(double T, double x); +// double ps_H2O(double T); +// double rho_H2O(double T); +// double cp_H2O(double T); +// double h_H2O(double T); +// double s_H2O(double T); +// +// /** Finished with the code from the paper. Now we need to +// * convert the molar values to mass-based units. */ +// double massToMole(double w); +// double molarToSpecific(double w, double value); +// +// static const bool debug; +// +//public: +// +// LiBrSolution(); +// +// double rho(double T, double p, double x); +// double c(double T, double p, double x); +// //double h(double T_K, double p, double x); +// double s(double T, double p, double x); +// double visc(double T, double p, double x); +// double cond(double T, double p, double x); +// double u(double T, double p, double x); +// double psat(double T, double x); +// double Tfreeze(double p, double x); +// +// /* Some functions can be inverted directly, those are listed +// * here. It is also possible to solve for other quantities, but +// * that involves some more sophisticated processing and is not +// * done here, but in the backend, T(h,p) for example. +// */ +// /// Temperature as a function of density, pressure and composition. +// double T_rho (double Dmass, double p, double x){throw NotImplementedError(format("%s (%d): T from density is not implemented for LiBr.",__FILE__,__LINE__));} +// /// Temperature as a function of heat capacities as a function of temperature, pressure and composition. +// double T_c (double Cmass, double p, double x){throw NotImplementedError(format("%s (%d): T from heat capacity is not implemented for LiBr.",__FILE__,__LINE__));} +// /// Temperature as a function of entropy as a function of temperature, pressure and composition. +// double T_s (double Smass, double p, double x){throw NotImplementedError(format("%s (%d): T from entropy is not implemented for LiBr.",__FILE__,__LINE__));} +// /// Temperature as a function of internal energy as a function of temperature, pressure and composition. +// double T_u (double Umass, double p, double x){throw NotImplementedError(format("%s (%d): T from internal energy is not implemented for LiBr.",__FILE__,__LINE__));} +// /// Temperature as a function of enthalpy, pressure and composition. +// //double T_h (double Hmass, double p, double x){throw NotImplementedError(format("%s (%d): T from enthalpy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));} +// /// Viscosity as a function of temperature, pressure and composition. +// double T_visc(double visc, double p, double x){throw NotImplementedError(format("%s (%d): T from viscosity is not implemented.",__FILE__,__LINE__));} +// /// Thermal conductivity as a function of temperature, pressure and composition. +// double T_cond(double cond, double p, double x){throw NotImplementedError(format("%s (%d): T from conductivity is not implemented.",__FILE__,__LINE__));} +// /// Saturation pressure as a function of temperature and composition. +// double T_psat(double psat, double x){throw NotImplementedError(format("%s (%d): T from psat is not implemented.",__FILE__,__LINE__));} +// /// Composition as a function of freezing temperature and pressure. +// double x_Tfreeze( double Tfreeze, double p){throw NotImplementedError(format("%s (%d): x from T_freeze is not implemented.",__FILE__,__LINE__));} +// +// +// /// Overwrite some standard functions that cannot be used with LiBr +// void setName(std::string name){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setDescription(std::string description){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setReference(std::string reference){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setTmax(double Tmax){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setTmin(double Tmin){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setxmax(double xmax){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setxmin(double xmin){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setTminPsat(double TminPsat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// +// void setTbase(double Tbase){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setxbase(double xbase){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// +// void setDensity(IncompressibleData density){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setSpecificHeat(IncompressibleData specific_heat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setViscosity(IncompressibleData viscosity){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setConductivity(IncompressibleData conductivity){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setPsat(IncompressibleData p_sat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setTfreeze(IncompressibleData T_freeze){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setVolToMass(IncompressibleData volToMass){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// void setMassToMole(IncompressibleData massToMole){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));} +// +// bool is_pure() {return false;}; +// +//}; +// @@ -170,13 +170,13 @@ public: void add_obj(IncompressibleFluid fluid_obj); /** \brief Get an IncompressibleFluid instance stored in this library - * + * * @param name Name of the fluid */ IncompressibleFluid& get(std::string name); /** \brief Get a CoolPropFluid instance stored in this library - * + * * @param key The index of the fluid in the map */ IncompressibleFluid& get(std::size_t key); diff --git a/src/Tests/TestObjects.cpp b/src/Tests/TestObjects.cpp index 9a875e73..cca144cc 100644 --- a/src/Tests/TestObjects.cpp +++ b/src/Tests/TestObjects.cpp @@ -198,14 +198,6 @@ CoolProp::IncompressibleFluid CoolPropTesting::incompressibleFluidObject(){ //CH3OH.setVolToMass(volume2mass); //CH3OH.setMassToMole(mass2mole); - //XLT.set_reference_state(25+273.15, 1.01325e5, 0.0, 0.0, 0.0); - double Tref = 20+273.15; - double pref = 101325; - double xref = 0.0; - double href = 0.0; - double sref = 0.0; - CH3OH.set_reference_state(Tref, pref, xref, href, sref); - /// A function to check coefficients and equation types. CH3OH.validate();