From 8f9c6ca4383c35bd9748edbc35fb6b1ee6ec85bc Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 12 Aug 2014 21:17:34 +0200 Subject: [PATCH] Added hs_anchor test (passing) and function get_state to allow to get access to states stored in the backends Signed-off-by: Ian Bell --- dev/codelite/coolprop.project | 1 - include/AbstractState.h | 23 ++++++++- src/AbstractState.cpp | 13 +++-- src/Backends/Helmholtz/Fluids/Ancillaries.cpp | 26 ++++++++++ src/Backends/Helmholtz/Fluids/FluidLibrary.h | 1 + .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 51 +++++++++++++++---- .../Helmholtz/HelmholtzEOSMixtureBackend.h | 3 ++ src/Backends/Helmholtz/TransportRoutines.cpp | 6 +-- src/Backends/Helmholtz/VLERoutines.cpp | 8 +-- src/CoolProp.cpp | 6 +-- 10 files changed, 111 insertions(+), 27 deletions(-) diff --git a/dev/codelite/coolprop.project b/dev/codelite/coolprop.project index 76968496..75e0c55a 100644 --- a/dev/codelite/coolprop.project +++ b/dev/codelite/coolprop.project @@ -121,7 +121,6 @@ - diff --git a/include/AbstractState.h b/include/AbstractState.h index c0e91e25..aa0327d9 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -204,6 +204,8 @@ protected: /// Using this backend, get the triple point temperature in K virtual long double calc_Ttriple(void){throw NotImplementedError("calc_Ttriple is not implemented for this backend");}; + /// Using this backend, get the triple point pressure in Pa + virtual long double calc_p_triple(void){throw NotImplementedError("calc_p_triple is not implemented for this backend");}; /// Using this backend, get the critical point temperature in K virtual long double calc_T_critical(void){throw NotImplementedError("calc_T_critical is not implemented for this backend");}; @@ -229,6 +231,10 @@ protected: virtual long double calc_melting_line(int param, int given, long double value){throw NotImplementedError("This backend does not implement calc_melting_line function");}; virtual int calc_phase(void){throw NotImplementedError("This backend does not implement calc_phase function");}; + + /// Using this backend, calculate a phase given by the state string + /// @param state A string that describes the state desired, one of "hs_anchor", "critical"/"crit", "reducing" + virtual const CoolProp::SimpleState calc_state(const std::string &state){throw NotImplementedError("calc_state is not implemented for this backend");}; public: @@ -278,7 +284,8 @@ public: void set_mass_fractions(const std::vector &mass_fractions){set_mass_fractions(std::vector(mass_fractions.begin(), mass_fractions.end()));}; void set_volu_fractions(const std::vector &volu_fractions){set_volu_fractions(std::vector(volu_fractions.begin(), volu_fractions.end()));}; - const CoolProp::SimpleState & get_reducing(){return _reducing;}; + const CoolProp::SimpleState & get_reducing_state(){return _reducing;}; + const CoolProp::SimpleState & get_state(const std::string &state){return calc_state(state);}; double keyed_output(int key); double trivial_keyed_output(int key); @@ -311,6 +318,9 @@ public: For mixtures, it is the exact critical point molar density calculated by the methods of Michelsen( \todo fill in reference) */ double rhomolar_critical(void); + + /// Return the triple point pressure + double p_triple(void); std::string name(){return calc_name();}; @@ -376,7 +386,18 @@ public: /// Return true if the fluid has a melting line - default is false, but can be re-implemented by derived class virtual bool has_melting_line(void){return false;}; + /// Return a value from the melting line + /// @param param The key for the parameter to be returned + /// @param given The key for the parameter that is given + /// @param value The value for the parameter that is given double melting_line(int param, int given, double value); + + /// Return the value from a saturation ancillary curve (if the backend implements it) + /// @param param The key for the parameter to be returned + /// @param Q The quality for the parameter that is given (0 = saturated liquid, 1 = saturated vapor) + /// @param given The key for the parameter that is given + /// @param value The value for the parameter that is given + double saturation_ancillary(parameters param, int Q, int given, double value); // ---------------------------------------- // Transport properties diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index aa71e645..45fb11fd 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -168,9 +168,9 @@ double AbstractState::trivial_keyed_output(int key) case iP_max: return pmax(); case iT_reducing: - return get_reducing().T; + return get_reducing_state().T; case irhomolar_reducing: - return get_reducing().rhomolar; + return get_reducing_state().rhomolar; case iP_critical: return this->p_critical(); case iT_critical: @@ -179,6 +179,8 @@ double AbstractState::trivial_keyed_output(int key) return this->rhomolar_critical(); case irhomass_critical: return this->rhomolar_critical()*molar_mass(); + case iP_triple: + return this->p_triple(); default: throw ValueError(format("This input [%d: \"%s\"] is not valid for trivial_keyed_output",key,get_parameter_information(key,"short").c_str())); } @@ -226,9 +228,9 @@ double AbstractState::keyed_output(int key) case imolar_mass: return molar_mass(); case iT_reducing: - return get_reducing().T; + return get_reducing_state().T; case irhomolar_reducing: - return get_reducing().rhomolar; + return get_reducing_state().rhomolar; case ispeed_sound: return speed_sound(); case ialpha0: @@ -282,6 +284,9 @@ double AbstractState::T_critical(void){ double AbstractState::p_critical(void){ return calc_p_critical(); } +double AbstractState::p_triple(void){ + return calc_p_triple(); +} double AbstractState::rhomolar_critical(void){ return calc_rhomolar_critical(); } diff --git a/src/Backends/Helmholtz/Fluids/Ancillaries.cpp b/src/Backends/Helmholtz/Fluids/Ancillaries.cpp index 95fc0bb5..69153bd7 100644 --- a/src/Backends/Helmholtz/Fluids/Ancillaries.cpp +++ b/src/Backends/Helmholtz/Fluids/Ancillaries.cpp @@ -377,4 +377,30 @@ TEST_CASE("Tests for values from melting lines", "[melting]") } } } + +TEST_CASE("Test that hs_anchor enthalpy/entropy agrees with EOS", "[ancillaries]") +{ + std::vector fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),','); + for (std::size_t i = 0; i < fluids.size(); ++i) + { + shared_ptr AS(CoolProp::AbstractState::factory("HEOS",fluids[i])); + + CoolProp::SimpleState hs_anchor = AS->get_state("hs_anchor"); + + // See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU + std::ostringstream ss1; + ss1 << "Check hs_anchor for " << fluids[i]; + SECTION(ss1.str(),"") + { + std::string note = "The enthalpy and entropy are hardcoded in the fluid JSON files. They MUST agree with the values calculated by the EOS"; + AS->update(CoolProp::DmolarT_INPUTS, hs_anchor.rhomolar, hs_anchor.T); + CAPTURE(hs_anchor.hmolar); + CAPTURE(hs_anchor.smolar); + double EOS_hmolar = AS->hmolar(); + double EOS_smolar = AS->smolar(); + CHECK( std::abs(EOS_hmolar - hs_anchor.hmolar) < 1e-3); + CHECK( std::abs(EOS_smolar - hs_anchor.smolar) < 1e-3); + } + } +} #endif \ No newline at end of file diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.h b/src/Backends/Helmholtz/Fluids/FluidLibrary.h index 05e672ba..378d4e64 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h +++ b/src/Backends/Helmholtz/Fluids/FluidLibrary.h @@ -320,6 +320,7 @@ protected: /// \todo: define limits of EOS better EOS.limits.Tmin = cpjson::get_double(satminL_state, "T"); + EOS.ptriple = cpjson::get_double(satminL_state, "p"); EOS.Ttriple = EOS.limits.Tmin; // BibTex keys diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index c7eba8e4..94978749 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -131,6 +131,27 @@ void HelmholtzEOSMixtureBackend::update_states(void) // Clear again just to be sure clear(); } +const CoolProp::SimpleState HelmholtzEOSMixtureBackend::calc_state(const std::string &state) +{ + if (is_pure_or_pseudopure) + { + if (!state.compare("hs_anchor")){ + return components[0]->pEOS->hs_anchor; + } + else if (!state.compare("max_sat_T")){ + return components[0]->pEOS->max_sat_T; + } + else if (!state.compare("max_sat_p")){ + return components[0]->pEOS->max_sat_p; + } + else{ + throw ValueError(format("This state [%s] is invalid to calc_state",state.c_str())); + } + } + else{ + throw ValueError(format("calc_state not supported for mixtures")); + } +}; long double HelmholtzEOSMixtureBackend::calc_gas_constant(void) { double summer = 0; @@ -376,6 +397,14 @@ long double HelmholtzEOSMixtureBackend::calc_Ttriple(void) } return summer; } +long double HelmholtzEOSMixtureBackend::calc_p_triple(void) +{ + double summer = 0; + for (unsigned int i = 0; i < components.size(); ++i){ + summer += mole_fractions[i]*components[i]->pEOS->ptriple; + } + return summer; +} std::string HelmholtzEOSMixtureBackend::calc_name(void) { if (components.size() != 1){ @@ -625,21 +654,21 @@ void HelmholtzEOSMixtureBackend::update(long input_pair, double value1, double v long double HelmholtzEOSMixtureBackend::calc_Bvirial() { - return 1/get_reducing().rhomolar*calc_alphar_deriv_nocache(0,1,mole_fractions,_tau,1e-12); + return 1/get_reducing_state().rhomolar*calc_alphar_deriv_nocache(0,1,mole_fractions,_tau,1e-12); } long double HelmholtzEOSMixtureBackend::calc_dBvirial_dT() { - long double dtau_dT =-get_reducing().T/pow(_T,2); - return 1/get_reducing().rhomolar*calc_alphar_deriv_nocache(1,1,mole_fractions,_tau,1e-12)*dtau_dT; + long double dtau_dT =-get_reducing_state().T/pow(_T,2); + return 1/get_reducing_state().rhomolar*calc_alphar_deriv_nocache(1,1,mole_fractions,_tau,1e-12)*dtau_dT; } long double HelmholtzEOSMixtureBackend::calc_Cvirial() { - return 1/pow(get_reducing().rhomolar,2)*calc_alphar_deriv_nocache(0,2,mole_fractions,_tau,1e-12); + return 1/pow(get_reducing_state().rhomolar,2)*calc_alphar_deriv_nocache(0,2,mole_fractions,_tau,1e-12); } long double HelmholtzEOSMixtureBackend::calc_dCvirial_dT() { - long double dtau_dT =-get_reducing().T/pow(_T,2); - return 1/pow(get_reducing().rhomolar,2)*calc_alphar_deriv_nocache(1,2,mole_fractions,_tau,1e-12)*dtau_dT; + long double dtau_dT =-get_reducing_state().T/pow(_T,2); + return 1/pow(get_reducing_state().rhomolar,2)*calc_alphar_deriv_nocache(1,2,mole_fractions,_tau,1e-12)*dtau_dT; } void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int other, long double value, bool &saturation_called) { @@ -738,7 +767,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot case iHmolar: { // Ancillaries are h-h_anchor, so add back h_anchor - long double h_liq = component.ancillaries.hL.evaluate(_TLanc) + component.EOSVector[0].hs_anchor.hmolar; + long double h_liq = component.ancillaries.hL.evaluate(_TLanc) + component.pEOS->hs_anchor.hmolar; long double h_liq_error_band = component.ancillaries.hL.get_max_abs_error(); long double h_vap = h_liq + component.ancillaries.hLV.evaluate(_TLanc); long double h_vap_error_band = h_liq_error_band + component.ancillaries.hLV.get_max_abs_error(); @@ -1211,8 +1240,8 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot void get_dtau_ddelta(HelmholtzEOSMixtureBackend *HEOS, long double T, long double rho, int index, long double &dtau, long double &ddelta) { - long double rhor = HEOS->get_reducing().rhomolar, - Tr = HEOS->get_reducing().T, + long double rhor = HEOS->get_reducing_state().rhomolar, + Tr = HEOS->get_reducing_state().T, dT_dtau = -pow(T, 2)/Tr, R = HEOS->gas_constant(), delta = rho/rhor, @@ -1441,8 +1470,8 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double HelmholtzEOSMixtureBackend *HEOS; solver_TP_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double p){ - this->HEOS = HEOS; this->T = T; this->p = p; this->rhor = HEOS->get_reducing().rhomolar; - this->tau = HEOS->get_reducing().T/T; this->R_u = HEOS->gas_constant(); + this->HEOS = HEOS; this->T = T; this->p = p; this->rhor = HEOS->get_reducing_state().rhomolar; + this->tau = HEOS->get_reducing_state().T/T; this->R_u = HEOS->gas_constant(); }; double call(double rhomolar){ this->rhomolar = rhomolar; diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index ab91f7e3..2f3a9538 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -48,6 +48,8 @@ public: bool has_melting_line(){ return is_pure_or_pseudopure && components[0]->ancillaries.melting_line.enabled();}; long double calc_melting_line(int param, int given, long double value); int calc_phase(void){return _phase;}; + + const CoolProp::SimpleState calc_state(const std::string &state); const std::vector &get_components(){return components;}; std::vector &get_K(){return K;}; @@ -141,6 +143,7 @@ public: long double calc_Tmax(void); long double calc_pmax(void); long double calc_Ttriple(void); + long double calc_p_triple(void); long double calc_pmax_sat(void); long double calc_Tmax_sat(void); void calc_Tmin_sat(long double &Tmin_satL, long double &Tmin_satV); diff --git a/src/Backends/Helmholtz/TransportRoutines.cpp b/src/Backends/Helmholtz/TransportRoutines.cpp index 66e29e1b..f7c09e95 100644 --- a/src/Backends/Helmholtz/TransportRoutines.cpp +++ b/src/Backends/Helmholtz/TransportRoutines.cpp @@ -527,9 +527,9 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers( GAMMA = data.GAMMA, zeta0 = data.zeta0, qD = data.qD, - Tc = HEOS.get_reducing().T, // [K] - rhoc = HEOS.get_reducing().rhomolar, // [mol/m^3] - Pcrit = HEOS.get_reducing().p, // [Pa] + Tc = HEOS.get_reducing_state().T, // [K] + rhoc = HEOS.get_reducing_state().rhomolar, // [mol/m^3] + Pcrit = HEOS.get_reducing_state().p, // [Pa] Tref, // [K] cp,cv,delta,num,zeta,mu,pi=M_PI,tau, OMEGA_tilde,OMEGA_tilde0; diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index a5ca1719..def74baf 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -20,7 +20,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l std::vector > J(3, std::vector(3,_HUGE)); HEOS->calc_reducing_state(); - const SimpleState & reduce = HEOS->get_reducing(); + const SimpleState & reduce = HEOS->get_reducing_state(); shared_ptr SatL = HEOS->SatL, SatV = HEOS->SatV; @@ -223,7 +223,7 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long std::vector > J(2, std::vector(2,_HUGE)); HEOS->calc_reducing_state(); - const SimpleState & reduce = HEOS->get_reducing(); + const SimpleState & reduce = HEOS->get_reducing_state(); shared_ptr SatL = HEOS->SatL, SatV = HEOS->SatV; @@ -392,7 +392,7 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE */ HEOS->calc_reducing_state(); - const SimpleState & reduce = HEOS->get_reducing(); + const SimpleState & reduce = HEOS->get_reducing_state(); long double R_u = HEOS->calc_gas_constant(); shared_ptr SatL = HEOS->SatL, SatV = HEOS->SatV; @@ -414,7 +414,7 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE // Use the density ancillary function as the starting point for the solver // If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature - if (T > 0.99*HEOS->get_reducing().T){ + if (T > 0.99*HEOS->get_reducing_state().T){ rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T-0.1); rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T-0.1); } diff --git a/src/CoolProp.cpp b/src/CoolProp.cpp index 7ebe8700..dd957412 100644 --- a/src/CoolProp.cpp +++ b/src/CoolProp.cpp @@ -729,7 +729,7 @@ void set_reference_stateS(std::string Ref, std::string reference_state) double deltah = HEOS->hmass() - 200000; // offset from 200000 J/kg enthalpy double deltas = HEOS->smass() - 1000; // offset from 1000 J/kg/K entropy double delta_a1 = deltas/(8.314472/HEOS->molar_mass()); - double delta_a2 = -deltah/(8.314472/HEOS->molar_mass()*HEOS->get_reducing().T); + double delta_a2 = -deltah/(8.314472/HEOS->molar_mass()*HEOS->get_reducing_state().T); HEOS->get_components()[0]->pEOS->alpha0.EnthalpyEntropyOffset.set(delta_a1, delta_a2, "IIR"); HEOS->update_states(); } @@ -741,7 +741,7 @@ void set_reference_stateS(std::string Ref, std::string reference_state) double deltah = HEOS->hmass() - 0; // offset from 0 J/kg enthalpy double deltas = HEOS->smass() - 0; // offset from 0 J/kg/K entropy double delta_a1 = deltas/(8.314472/HEOS->molar_mass()); - double delta_a2 = -deltah/(8.314472/HEOS->molar_mass()*HEOS->get_reducing().T); + double delta_a2 = -deltah/(8.314472/HEOS->molar_mass()*HEOS->get_reducing_state().T); HEOS->get_components()[0]->pEOS->alpha0.EnthalpyEntropyOffset.set(delta_a1, delta_a2, "ASHRAE"); HEOS->update_states(); } @@ -753,7 +753,7 @@ void set_reference_stateS(std::string Ref, std::string reference_state) double deltah = HEOS->hmass() - 0; // offset from 0 kJ/kg enthalpy double deltas = HEOS->smass() - 0; // offset from 0 kJ/kg/K entropy double delta_a1 = deltas/(8.314472/HEOS->molar_mass()); - double delta_a2 = -deltah/(8.314472/HEOS->molar_mass()*HEOS->get_reducing().T); + double delta_a2 = -deltah/(8.314472/HEOS->molar_mass()*HEOS->get_reducing_state().T); HEOS->get_components()[0]->pEOS->alpha0.EnthalpyEntropyOffset.set(delta_a1, delta_a2, "NBP"); HEOS->update_states(); }