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();
}