diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index aef35924..2ddbb779 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -1086,4 +1086,76 @@ void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, parameters ot } } +void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS) +{ + if (HEOS.imposed_phase_index != iphase_not_imposed) + { + // Use the phase defined by the imposed phase + HEOS._phase = HEOS.imposed_phase_index; + } + else + { + CoolProp::SimpleState &crit = HEOS.components[0]->pEOS->reduce; + CoolProp::SimpleState &tripleL = HEOS.components[0]->triple_liquid, + &tripleV = HEOS.components[0]->triple_vapor; + // Enthalpy at solid line + double hsolid = (tripleV.hmolar-tripleL.hmolar)/(tripleV.smolar-tripleL.smolar)*(HEOS.smolar()-tripleL.smolar) + tripleL.hmolar; + // Part A - first check if HS is below triple line formed by connecting the triple point states + // If so, it is solid, and not supported + if (HEOS.hmolar() < hsolid){ + throw ValueError(format("Enthalpy [%g J/mol] is below solid enthalpy [%g J/mol]", HEOS.hmolar(), hsolid)); + } + // Part B - Check lower limit + else if (HEOS.smolar() < tripleL.smolar){ + // If fluid is other than water (which can have solutions below tripleL), cannot have any solution, fail + if (upper(HEOS.name()) != "Water"){ + throw ValueError(format("Entropy [%g J/mol/K] is below triple point liquid entropy [%g J/mol/K]", HEOS.smolar(), tripleL.smolar)); + } + } + // Part C - Check higher limit + else if (HEOS.smolar() > tripleV.smolar){ + throw ValueError(format("Entropy [%g J/mol/K] is above triple point vapor entropy [%g J/mol/K]", HEOS.smolar(), tripleV.smolar)); + } + // Part D - if s < sc, a few options are possible. It could be two-phase, or liquid (more likely), or gas (less likely) + else if (HEOS.smolar() < crit.smolar){ + // Do a saturation call for given s for liquid, first with ancillaries, then with full saturation call + // Check if above the saturation enthalpy each time + // D1: It is above hsat + // It is liquid, supercritical liquid, or out of range + // Come up with some guess values, maybe start at the saturation curve with under-relaxation? + // D2: It is below hsat + // Either two-phase, or vapor (for funky saturation curves like the siloxanes) + // Do a saturation_h call for given h on the vapor line to determine whether two-phase or vapor + // D2a: It is below ssat + // two-phase + // D2b: It is above ssat + // vapor + throw ValueError("partD"); + } + // Part E - HEOS.smolar() > crit.hmolar > tripleV.smolar + else{ + throw ValueError("partE"); + } + } + + /*if (HEOS.isHomogeneousPhase() && !ValidNumber(HEOS._p)) + { + switch (other) + { + case iDmolar: + break; + case iHmolar: + HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._hmolar, iHmolar); break; + case iSmolar: + HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._smolar, iSmolar); break; + case iUmolar: + HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._umolar, iUmolar); break; + default: + break; + } + HEOS.calc_pressure(); + HEOS._Q = -1; + }*/ +} + } /* namespace CoolProp */ diff --git a/src/Backends/Helmholtz/FlashRoutines.h b/src/Backends/Helmholtz/FlashRoutines.h index 2fd16fdd..03c17dec 100644 --- a/src/Backends/Helmholtz/FlashRoutines.h +++ b/src/Backends/Helmholtz/FlashRoutines.h @@ -79,6 +79,10 @@ public: /// @param HEOS The HelmholtzEOSMixtureBackend to be used /// @param other The index for the other input from CoolProp::parameters; allowed values are iP, iHmolar, iSmolar, iUmolar static void PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, parameters other); + + /// A flash routine for (H,S) + /// @param HEOS The HelmholtzEOSMixtureBackend to be used + static void HS_flash(HelmholtzEOSMixtureBackend &HEOS); }; diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index cb11a6cb..1c5f68b9 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -778,6 +778,8 @@ void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double _p = value1; _smolar = value2; FlashRoutines::HSU_P_flash(*this, iSmolar); break; case PUmolar_INPUTS: _p = value1; _umolar = value2; FlashRoutines::HSU_P_flash(*this, iUmolar); break; + case HmolarSmolar_INPUTS: + _hmolar = value1; _smolar = value2; FlashRoutines::HS_flash(*this); break; case QT_INPUTS: _Q = value1; _T = value2; FlashRoutines::QT_flash(*this); break; case PQ_INPUTS: diff --git a/src/Tests/CoolProp-Tests.cpp b/src/Tests/CoolProp-Tests.cpp index f4aac67c..20605038 100644 --- a/src/Tests/CoolProp-Tests.cpp +++ b/src/Tests/CoolProp-Tests.cpp @@ -1264,7 +1264,7 @@ TEST_CASE("Triple point checks", "[triple_point]") } } -TEST_CASE("Test that states agree with CoolProp", "[reference_states]") +TEST_CASE("Test that reference states are correct", "[reference_states]") { std::vector fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),','); for (std::size_t i = 0; i < fluids.size(); ++i) @@ -1338,6 +1338,40 @@ TEST_CASE("Test that states agree with CoolProp", "[reference_states]") } } } - } + +TEST_CASE("Test that HS solver works for a few fluids", "[HS_solver]") +{ + std::vector fluids; fluids.push_back("Propane"); fluids.push_back("D4"); fluids.push_back("Water"); + for (std::size_t i = 0; i < fluids.size(); ++i) + { + std::vector fl(1,fluids[i]); + shared_ptr HEOS(new CoolProp::HelmholtzEOSMixtureBackend(fl)); + for (double p = HEOS->p_triple()*10; p < HEOS->pmax(); p *= 10) + { + double Tmin = HEOS->Ttriple(); + double Tmax = HEOS->Tmax(); + for (double T = Tmin; T < Tmax; T += 100) + { + std::ostringstream ss; + ss << "Check HS for " << fluids[i] << " for T=" << T << ", p=" << p; + SECTION(ss.str(),"") + { + CHECK_NOTHROW(HEOS->update(PT_INPUTS, p, T)); + std::ostringstream ss1; + ss1 << "h=" << HEOS->hmolar() << ", s=" << HEOS->smolar(); + SECTION(ss1.str(),"") + { + CAPTURE(T); + CAPTURE(p); + CAPTURE(HEOS->hmolar()); + CAPTURE(HEOS->smolar()); + CHECK_NOTHROW(HEOS->update(HmolarSmolar_INPUTS, HEOS->hmolar(), HEOS->smolar())); + } + } + } + } + } +} + #endif