diff --git a/include/AbstractState.h b/include/AbstractState.h index db7ac522..feb51bba 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -232,6 +232,12 @@ protected: virtual void update_states(void){throw NotImplementedError("This backend does not implement update_states function");}; virtual long double calc_melting_line(int param, int given, long double value){throw NotImplementedError("This backend does not implement calc_melting_line function");}; + + /// @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 + virtual long double calc_saturation_ancillary(parameters param, int Q, parameters given, double value){throw NotImplementedError("This backend does not implement calc_saturation_ancillary");}; virtual phases calc_phase(void){throw NotImplementedError("This backend does not implement calc_phase function");}; @@ -430,7 +436,7 @@ public: /// @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(int param, int Q, int given, double value); + double saturation_ancillary(parameters param, int Q, parameters given, double value); // ---------------------------------------- // Transport properties diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 241c8b37..144c2eae 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -330,6 +330,9 @@ double AbstractState::conductivity(void){ double AbstractState::melting_line(int param, int given, double value){ return calc_melting_line(param, given, value); } +double AbstractState::saturation_ancillary(parameters param, int Q, parameters given, double value){ + return calc_saturation_ancillary(param, Q, given, value); +} double AbstractState::surface_tension(void){ if (!_surface_tension) _surface_tension = calc_surface_tension(); return _surface_tension; diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 0f88086e..ed3bd384 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -174,6 +174,64 @@ long double HelmholtzEOSMixtureBackend::calc_molar_mass(void) } return summer; } +long double HelmholtzEOSMixtureBackend::calc_saturation_ancillary(parameters param, int Q, parameters given, double value) +{ + if (is_pure_or_pseudopure) + { + if (param == iP && given == iT){ + // p = f(T), direct evaluation + switch (Q) + { + case 0: + return components[0]->ancillaries.pL.evaluate(value); + case 1: + return components[0]->ancillaries.pV.evaluate(value); + } + } + else if (param == iT && given == iP){ + // T = f(p), inverse evaluation + switch (Q) + { + case 0: + return components[0]->ancillaries.pL.invert(value); + case 1: + return components[0]->ancillaries.pV.invert(value); + } + } + else if (param == iDmolar && given == iT){ + // rho = f(T), inverse evaluation + switch (Q) + { + case 0: + return components[0]->ancillaries.rhoL.evaluate(value); + case 1: + return components[0]->ancillaries.rhoV.evaluate(value); + } + } + else if (param == iT && given == iDmolar){ + // T = f(rho), inverse evaluation + switch (Q) + { + case 0: + return components[0]->ancillaries.rhoL.invert(value); + case 1: + return components[0]->ancillaries.rhoV.invert(value); + } + } + else{ + throw ValueError(format("calc of %s given %s is invalid in calc_saturation_ancillary", + get_parameter_information(param,"short").c_str(), + get_parameter_information(given,"short").c_str())); + } + + throw ValueError(format("Q [%d] is invalid in calc_saturation_ancillary", Q)); + } + else + { + throw NotImplementedError(format("calc_saturation_ancillary not implemented for mixtures")); + } +} + long double HelmholtzEOSMixtureBackend::calc_melting_line(int param, int given, long double value) { if (is_pure_or_pseudopure) @@ -1624,23 +1682,26 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double // Calculate a guess value using SRK equation of state rhomolar_guess = solver_rho_Tp_SRK(T, p, phase); - if (phase == iphase_gas || phase == iphase_supercritical_gas) + // A gas-like phase, ideal gas might not be the perfect model, but probably good enough + if (phase == iphase_gas || phase == iphase_supercritical_gas || iphase_supercritical) { if (rhomolar_guess < 0 || !ValidNumber(rhomolar_guess)) // If the guess is bad, probably high temperature, use ideal gas { rhomolar_guess = p/(gas_constant()*T); } } - else + // It's liquid at subcritical pressure, we can use ancillaries as a backup + else if (phase == iphase_liquid) { - if (phase == iphase_liquid) - { - long double _rhoLancval = static_cast(components[0]->ancillaries.rhoL.evaluate(T)); - if (!ValidNumber(rhomolar_guess) || rhomolar_guess < _rhoLancval){ - rhomolar_guess = _rhoLancval; - } + long double _rhoLancval = static_cast(components[0]->ancillaries.rhoL.evaluate(T)); + if (!ValidNumber(rhomolar_guess) || rhomolar_guess < _rhoLancval){ + rhomolar_guess = _rhoLancval; } } + else if (phase == iphase_supercritical_liquid){ + long double T = static_cast(components[0]->ancillaries.pL.invert(_p)); + rhomolar_guess = components[0]->ancillaries.rhoL.evaluate(T); + } } try{ diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index 21dfb6d4..15b07c18 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -53,6 +53,7 @@ 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); phases calc_phase(void){return _phase;}; + long double calc_saturation_ancillary(parameters param, int Q, parameters given, double value); const CoolProp::SimpleState &calc_state(const std::string &state); diff --git a/src/Tests/CoolProp-Tests.cpp b/src/Tests/CoolProp-Tests.cpp index 6eab8428..c415c2d3 100644 --- a/src/Tests/CoolProp-Tests.cpp +++ b/src/Tests/CoolProp-Tests.cpp @@ -916,6 +916,64 @@ TEST_CASE("Test partial derivatives using PropsSI", "[derivatives]") } } + +TEST_CASE("Ancillary functions", "[ancillary]") +{ + 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])); + + double Tc = AS->T_critical(); + double Tt = AS->Ttriple(); + + for (double f = 0.1; f < 1; f += 0.2) + { + double T = f*Tc + (1-f)*Tt; + + + std::ostringstream ss1; + ss1 << "Pressure error < 2% for fluid " << fluids[i] << " at " << T << " K"; + SECTION(ss1.str(), "") + { + AS->update(CoolProp::QT_INPUTS, 0, T); + double p_EOS = AS->p(); + double p_anc = AS->saturation_ancillary(CoolProp::iP, 0, CoolProp::iT, T); + double err = std::abs(p_EOS-p_anc)/p_anc; + CAPTURE(p_EOS); + CAPTURE(p_anc); + CAPTURE(T); + CHECK(err < 0.02); + } + std::ostringstream ss2; + ss2 << "Liquid density error < 3% for fluid " << fluids[i] << " at " << T << " K"; + SECTION(ss2.str(), "") + { + AS->update(CoolProp::QT_INPUTS, 0, T); + double rho_EOS = AS->rhomolar(); + double rho_anc = AS->saturation_ancillary(CoolProp::iDmolar, 0, CoolProp::iT, T); + double err = std::abs(rho_EOS-rho_anc)/rho_anc; + CAPTURE(rho_EOS); + CAPTURE(rho_anc); + CAPTURE(T); + CHECK(err < 0.03); + } + std::ostringstream ss3; + ss3 << "Vapor density error < 3% for fluid " << fluids[i] << " at " << T << " K"; + SECTION(ss3.str(), "") + { + AS->update(CoolProp::QT_INPUTS, 1, T); + double rho_EOS = AS->rhomolar(); + double rho_anc = AS->saturation_ancillary(CoolProp::iDmolar, 1, CoolProp::iT, T); + double err = std::abs(rho_EOS-rho_anc)/rho_anc; + CAPTURE(rho_EOS); + CAPTURE(rho_anc); + CAPTURE(T); + CHECK(err < 0.03); + } + } + } +} + //TEST_CASE("Test that states agree with CoolProp", "[states]") //{ // std::vector fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),',');