diff --git a/include/DataStructures.h b/include/DataStructures.h index c4077e7a..7d8c3564 100644 --- a/include/DataStructures.h +++ b/include/DataStructures.h @@ -22,6 +22,15 @@ struct SimpleState bool is_valid(){return ValidNumber(rhomolar) && ValidNumber(T) && ValidNumber(hmolar) && ValidNumber(p);} }; +/// A modified class for the state point at the maximum saturation entropy on the vapor curve +struct SsatSimpleState : public SimpleState +{ + enum SsatSimpleStateEnum {SSAT_MAX_NOT_SET=0, SSAT_MAX_DOESNT_EXIST, SSAT_MAX_DOES_EXIST}; + SsatSimpleStateEnum exists; + SsatSimpleState(){ SimpleState(); } +}; + + /// -------------------------------------------------- /// Define some constants that will be used throughout /// -------------------------------------------------- @@ -166,8 +175,10 @@ enum fluid_types{FLUID_TYPE_PURE, FLUID_TYPE_PSEUDOPURE, FLUID_TYPE_REFPROP, FLU enum input_pairs{ QT_INPUTS, ///< Molar quality, Temperature in K PQ_INPUTS, ///< Pressure in Pa, Molar quality - QS_INPUTS, ///< Molar quality, Entropy in J/mol/K - HQ_INPUTS, ///< Enthalpy in J/mol, Molar quality + QSmolar_INPUTS, ///< Molar quality, Entropy in J/mol/K + QSmass_INPUTS, ///< Molar quality, Entropy in J/kg/K + HmolarQ_INPUTS, ///< Enthalpy in J/mol, Molar quality + HmassQ_INPUTS, ///< Enthalpy in J/kg, Molar quality PT_INPUTS, ///< Pressure in Pa, Temperature in K diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index dd322736..f23e4cad 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -136,15 +136,23 @@ void FlashRoutines::PT_flash(HelmholtzEOSMixtureBackend &HEOS) HEOS._rhomolar = HEOS.solver_rho_Tp(HEOS._T, HEOS._p); HEOS._Q = -1; } -void FlashRoutines::HQ_flash(HelmholtzEOSMixtureBackend &HEOS) +void FlashRoutines::HQ_flash(HelmholtzEOSMixtureBackend &HEOS, long double Tguess) { + SaturationSolvers::saturation_PHSU_pure_options options; + options.use_logdelta = false; + HEOS.specify_phase(iphase_twophase); + if (Tguess < 0){ + options.use_guesses = true; + options.T = Tguess; + CoolProp::SaturationAncillaryFunction &rhoL = HEOS.get_components()[0]->ancillaries.rhoL; + CoolProp::SaturationAncillaryFunction &rhoV = HEOS.get_components()[0]->ancillaries.rhoV; + options.rhoL = rhoL.evaluate(Tguess); + options.rhoV = rhoV.evaluate(Tguess); + } if (HEOS.is_pure_or_pseudopure){ if (std::abs(HEOS.Q()-1) > 1e-10){throw ValueError(format("non-unity quality not currently allowed for HQ_flash"));} // Do a saturation call for given h for vapor, first with ancillaries, then with full saturation call - SaturationSolvers::saturation_PHSU_pure_options options; options.specified_variable = SaturationSolvers::saturation_PHSU_pure_options::IMPOSED_HV; - options.use_logdelta = false; - HEOS.specify_phase(iphase_twophase); SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS.hmolar(), options); HEOS._p = HEOS.SatV->p(); HEOS._T = HEOS.SatV->T(); @@ -159,7 +167,14 @@ void FlashRoutines::QS_flash(HelmholtzEOSMixtureBackend &HEOS) { if (HEOS.is_pure_or_pseudopure){ - if (std::abs(HEOS.Q()) < 1e-10){ + if (std::abs(HEOS.smolar() - HEOS.get_state("reducing").smolar) < 0.001) + { + HEOS._p = HEOS.p_critical(); + HEOS._T = HEOS.T_critical(); + HEOS._rhomolar = HEOS.rhomolar_critical(); + HEOS._phase = iphase_critical_point; + } + else if (std::abs(HEOS.Q()) < 1e-10){ // Do a saturation call for given s for liquid, first with ancillaries, then with full saturation call SaturationSolvers::saturation_PHSU_pure_options options; options.specified_variable = SaturationSolvers::saturation_PHSU_pure_options::IMPOSED_SL; @@ -1141,6 +1156,8 @@ void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, parameters ot HEOS.calc_pressure(); HEOS._Q = -1; } + // Update the state for conditions where the state was guessed + } void FlashRoutines::HS_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, long double hmolar_spec, long double smolar_spec, HS_flash_singlephaseOptions &options) @@ -1194,7 +1211,7 @@ void FlashRoutines::HS_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, long throw ValueError(format("HS_flash_singlephase took too many iterations; residual is %g; prior was %g", resid, resid_old)); } } - while(std::abs(resid) > 1e-10); + while(std::abs(resid) > 1e-9); } void FlashRoutines::HS_flash_generate_TP_singlephase_guess(HelmholtzEOSMixtureBackend &HEOS, double &T, double &p) { @@ -1219,47 +1236,71 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS) // Find maxima states if needed // Cache the maximum enthalpy saturation state; - //HEOS.calc_hsat_max(); + HEOS.calc_hsat_max(); // For weird fluids like the siloxanes, there can also be a maximum // entropy along the vapor saturation line. Try to find it if it has one - // HEOS.calc_ssat_max(); + HEOS.calc_ssat_max(); CoolProp::SimpleState crit = HEOS.get_state("reducing"); CoolProp::SimpleState &tripleL = HEOS.components[0]->triple_liquid, &tripleV = HEOS.components[0]->triple_vapor; - // Enthalpy at solid line + + double first_maxima_in_saturation_entropy; + if (HEOS.ssat_max.exists == SsatSimpleState::SSAT_MAX_DOES_EXIST){ + first_maxima_in_saturation_entropy = HEOS.ssat_max.smolar; + } + else{ + first_maxima_in_saturation_entropy = tripleV.smolar; + } + + double h1 = HEOS.hmolar(), s1 = HEOS.smolar(); + // Enthalpy at solid line for given entropy 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)); + if (HEOS.hmolar() < hsolid-0.1){ // -0.1 is for a buffer + throw ValueError(format("Enthalpy [%g J/mol] is below solid enthalpy [%g J/mol] for entropy [%g J/mol/K]", HEOS.hmolar(), hsolid-0.1, HEOS.smolar())); } - /*// 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)); + /* Now check up to the first maxima in saturated vapor entropy. + * For nicely behaved fluids, this means all the way up to the triple point vapor + * For not-nicely behaved fluids with a local maxima on the saturated vapor entropy curve, + * the entropy at the local maxima in entropy + * + * In this section, there can only be one solution for the saturation temperature, which is why the solver + * is divided up in this way. + */ + else if (HEOS.smolar() < first_maxima_in_saturation_entropy){ + double Q; + if (HEOS.smolar() < crit.smolar){ + Q = 0; // Liquid part + } + else{ + Q = 1; // Vapor part } - }*/ - // Part C - 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){ - - // Update the temporary instance with saturated liquid entropy - HEOS_copy->update(QS_INPUTS, 0, HEOS.smolar()); + // Update the temporary instance with saturated entropy + HEOS_copy->update(QSmolar_INPUTS, Q, HEOS.smolar()); + // Check if above the saturation enthalpy for given entropy + // If it is, the inputs are definitely single-phase. We are done here + double h1 = HEOS.hmolar(), h2 = HEOS_copy->hmolar(); if (HEOS.hmolar() > HEOS_copy->hmolar()){ solution = single_phase_solution; } else{ - // C2: It is below hsatL(ssatL) + // C2: It is below hsat(ssat) // 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 - // Update the temporary instance with saturated vapor enthalpy - HEOS_copy->update(HQ_INPUTS, HEOS.hmolar(), 1); - if (HEOS.smolar() > HEOS_copy->smolar()) - { + /* If enthalpy is between enthalpy at maxima in h and triple + * point enthalpy, search in enthalpy to find possible solutions + * that yield the correct enthalpy + * + * There should only be one solution since we have already bypassed the local maxima in enthalpy + */ + + HEOS_copy->update_HmolarQ_with_guessT(HEOS.hmolar(), 1, HEOS.hsat_max.T); + + if (HEOS.smolar() > HEOS_copy->smolar()){ solution = single_phase_solution; } else{ @@ -1268,22 +1309,6 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS) } } } - // Part C - if tripleV.s > s > sc - else if (HEOS.smolar() > crit.smolar && HEOS.smolar() < tripleV.smolar){ - // Do a saturation_s call for given s on the vapor line to determine whether two-phase or vapor - // Update the temporary instance with saturated vapor entropy - HEOS_copy->update(QS_INPUTS, 1, HEOS.smolar()); - - double h1 = HEOS.hmolar(), h2 = HEOS_copy->hmolar(); - if (HEOS.hmolar() > HEOS_copy->hmolar()){ - // D2b: It is above hsatV(ssatV) --> gas - solution = single_phase_solution; - } - else{ - // C2a: It is below ssatV(hsatV) --> two-phase - solution = two_phase_solution; - } - } // Part D - Check higher limit else if (HEOS.smolar() > tripleV.smolar){ solution = single_phase_solution; @@ -1311,12 +1336,22 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS) break; } catch(std::exception &e){ - HEOS_copy->update(DmolarT_INPUTS, HEOS.rhomolar_critical()*1.3, HEOS.Tmax()); - // Do the flash calcs starting from the guess value - HS_flash_singlephase(*HEOS_copy, HEOS.hmolar(), HEOS.smolar(), options); - // Copy the results - HEOS.update(DmolarT_INPUTS, HEOS_copy->rhomolar(), HEOS_copy->T()); - break; + try{ + // Trying again with another guessed value + HEOS_copy->update(DmolarT_INPUTS, HEOS.rhomolar_critical()*1.3, HEOS.Tmax()); + HS_flash_singlephase(*HEOS_copy, HEOS.hmolar(), HEOS.smolar(), options); + // Copy the results + HEOS.update(DmolarT_INPUTS, HEOS_copy->rhomolar(), HEOS_copy->T()); + break; + } + catch (std::exception &e){ + // Trying again with another guessed value + HEOS_copy->update(DmolarT_INPUTS, HEOS.rhomolar_critical(), 0.5*HEOS.Tmax() + 0.5*HEOS.T_critical()); + HS_flash_singlephase(*HEOS_copy, HEOS.hmolar(), HEOS.smolar(), options); + // Copy the results + HEOS.update(DmolarT_INPUTS, HEOS_copy->rhomolar(), HEOS_copy->T()); + break; + } } } case two_phase_solution: diff --git a/src/Backends/Helmholtz/FlashRoutines.h b/src/Backends/Helmholtz/FlashRoutines.h index 920c136b..2b53eb7c 100644 --- a/src/Backends/Helmholtz/FlashRoutines.h +++ b/src/Backends/Helmholtz/FlashRoutines.h @@ -42,7 +42,8 @@ public: /// Flash for given molar enthalpy and (molar) quality /// @param HEOS The HelmholtzEOSMixtureBackend to be used - static void HQ_flash(HelmholtzEOSMixtureBackend &HEOS); + /// @param Tguess (optional) The guess temperature in K to start from, ignored if < 0 + static void HQ_flash(HelmholtzEOSMixtureBackend &HEOS, long double Tguess = -1); /// Flash for mixture given temperature or pressure and (molar) quality /// @param HEOS The HelmholtzEOSMixtureBackend to be used diff --git a/src/Backends/Helmholtz/Fluids/Ancillaries.cpp b/src/Backends/Helmholtz/Fluids/Ancillaries.cpp index 5b78b48a..ca0d5053 100644 --- a/src/Backends/Helmholtz/Fluids/Ancillaries.cpp +++ b/src/Backends/Helmholtz/Fluids/Ancillaries.cpp @@ -114,7 +114,7 @@ double SaturationAncillaryFunction::invert(double value, double min_bound, doubl try{ // Safe to expand the domain a little bit to lower temperature, absolutely cannot exceed Tmax // because then you get (negative number)^(double) which is undefined. - return Brent(resid,min_bound,max_bound,DBL_EPSILON,1e-12,100,errstring); + return Brent(resid,min_bound,max_bound,DBL_EPSILON,1e-10,100,errstring); } catch(std::exception &e){ return Secant(resid,max_bound, -0.01, 1e-12, 100, errstring); diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 3255c6b2..65d4fa77 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -658,6 +658,20 @@ void HelmholtzEOSMixtureBackend::update_DmolarT_direct(long double rhomolar, lon post_update(); } +void HelmholtzEOSMixtureBackend::update_HmolarQ_with_guessT(long double hmolar, long double Q, long double Tguess) +{ + CoolProp::input_pairs pair = CoolProp::HmolarQ_INPUTS; + // Set up the state + pre_update(pair, hmolar, Q); + + _hmolar = hmolar; + _Q = Q; + FlashRoutines::HQ_flash(*this, Tguess); + + // Cleanup + post_update(); +} + void HelmholtzEOSMixtureBackend::update_TP_guessrho(long double T, long double p, long double rhomolar_guess) { CoolProp::input_pairs pair = PT_INPUTS; @@ -784,9 +798,9 @@ void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double _Q = value1; _T = value2; FlashRoutines::QT_flash(*this); break; case PQ_INPUTS: _p = value1; _Q = value2; FlashRoutines::PQ_flash(*this); break; - case QS_INPUTS: + case QSmolar_INPUTS: _Q = value1; _smolar = value2; FlashRoutines::QS_flash(*this); break; - case HQ_INPUTS: + case HmolarQ_INPUTS: _hmolar = value1; _Q = value2; FlashRoutines::HQ_flash(*this); break; default: throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str())); @@ -1112,16 +1126,27 @@ void HelmholtzEOSMixtureBackend::calc_ssat_max(void) return HEOS->SatV->first_partial_deriv(iSmolar,iT,iP)+HEOS->SatV->first_partial_deriv(iSmolar,iP,iT)/dTdp_along_sat; } }; - if (!ssat_max.is_valid()) + if (!ssat_max.is_valid() && ssat_max.exists != SsatSimpleState::SSAT_MAX_DOESNT_EXIST) { - Residual resid(*this); - std::string errstr; - Secant(resid, T_critical() - 1, 0.1, 1e-8, 30, errstr); - ssat_max.T = resid.HEOS->T(); - ssat_max.p = resid.HEOS->p(); - ssat_max.rhomolar = resid.HEOS->rhomolar(); - ssat_max.hmolar = resid.HEOS->hmolar(); - ssat_max.smolar = resid.HEOS->smolar(); + shared_ptr HEOS_copy(new CoolProp::HelmholtzEOSMixtureBackend(get_components())); + Residual resid(*HEOS_copy); + CoolProp::SimpleState &tripleV = HEOS_copy->get_components()[0]->triple_vapor; + double v1 = resid.call(hsat_max.T); + double v2 = resid.call(tripleV.T); + // If there is a sign change, there is a maxima, otherwise there is no local maxima/minima + if (v1*v2 < 0){ + std::string errstr; + Brent(resid, hsat_max.T, tripleV.T, DBL_EPSILON, 1e-8, 30, errstr); + ssat_max.T = resid.HEOS->T(); + ssat_max.p = resid.HEOS->p(); + ssat_max.rhomolar = resid.HEOS->rhomolar(); + ssat_max.hmolar = resid.HEOS->hmolar(); + ssat_max.smolar = resid.HEOS->smolar(); + ssat_max.exists = SsatSimpleState::SSAT_MAX_DOES_EXIST; + } + else{ + ssat_max.exists = SsatSimpleState::SSAT_MAX_DOESNT_EXIST; + } } } void HelmholtzEOSMixtureBackend::calc_hsat_max(void) @@ -1141,9 +1166,10 @@ void HelmholtzEOSMixtureBackend::calc_hsat_max(void) }; if (!hsat_max.is_valid()) { - Residualhmax residhmax(*this); + shared_ptr HEOS_copy(new CoolProp::HelmholtzEOSMixtureBackend(get_components())); + Residualhmax residhmax(*HEOS_copy); std::string errstrhmax; - Secant(residhmax, T_critical() - 1, 0.1, 1e-8, 30, errstrhmax); + Brent(residhmax, T_critical() - 0.1, HEOS_copy->Ttriple() + 1, DBL_EPSILON, 1e-8, 30, errstrhmax); hsat_max.T = residhmax.HEOS->T(); hsat_max.p = residhmax.HEOS->p(); hsat_max.rhomolar = residhmax.HEOS->rhomolar(); diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index bbd3d58f..99ab02e2 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -41,7 +41,8 @@ public: shared_ptr Reducing; ExcessTerm Excess; PhaseEnvelopeData PhaseEnvelope; - SimpleState ssat_max, hsat_max; + SimpleState hsat_max; + SsatSimpleState ssat_max; friend class FlashRoutines; // Allows the static methods in the FlashRoutines class to have access to all the protected members and methods of this class friend class TransportRoutines; // Allows the static methods in the TransportRoutines class to have access to all the protected members and methods of this class @@ -77,10 +78,21 @@ public: void resize(unsigned int N); shared_ptr SatL, SatV; ///< + /** \brief The standard update function + * @param input_pair The pair of inputs that will be provided + * @param value1 The first input value + * @param value2 The second input value + */ void update(CoolProp::input_pairs input_pair, double value1, double value2); + /** \brief Update with TP and a guess for rho + * @param T Temperature in K + * @param p Pressure in Pa + * @param rho_guess Density in mol/m^3 guessed + */ void update_TP_guessrho(long double T, long double p, long double rho_guess); void update_DmolarT_direct(long double rhomolar, long double T); + void update_HmolarQ_with_guessT(long double hmolar, long double Q, long double Tguess); /** \brief Set the components of the mixture * diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index 3f3f6426..961b54ab 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -178,6 +178,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l HEOS.calc_reducing_state(); const SimpleState & reduce = HEOS.get_reducing_state(); + CoolProp::SimpleState crit = HEOS.get_state("reducing"); shared_ptr SatL = HEOS.SatL, SatV = HEOS.SatV; @@ -235,27 +236,53 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l double Tmax = HEOS.calc_Tmax_sat(); try{ T = Brent(resid, Tmin-3, Tmax + 1, DBL_EPSILON, 1e-10, 50, errstr); } catch(std::exception &e){ - throw ValueError(format("Unable to invert ancillary equation for hV: %s",e.what())); + shared_ptr HEOS_copy(new HelmholtzEOSMixtureBackend(HEOS.get_components())); + HEOS_copy->update(QT_INPUTS, 1, Tmin); double hTmin = HEOS_copy->hmolar(); + HEOS_copy->update(QT_INPUTS, 1, Tmax); double hTmax = HEOS_copy->hmolar(); + T = (Tmax-Tmin)/(hTmax-hTmin)*(HEOS.hmolar()-hTmin) + Tmin; } } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SL) { - CoolProp::SaturationAncillaryFunction &anc = HEOS.get_components()[0]->ancillaries.sL; + CoolPropFluid &component = *(HEOS.get_components()[0]); + CoolProp::SaturationAncillaryFunction &anc = component.ancillaries.sL; CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor"); - // Ancillary is deltas = s - hs_anchor.s - try{ - T = anc.invert(specified_value - hs_anchor.smolar); + // If near the critical point, use a near critical guess value for T + if (std::abs(HEOS.smolar() - crit.smolar) < std::abs(component.ancillaries.sL.get_max_abs_error())) + { + T = std::max(0.99*crit.T, crit.T-0.1); + } + else{ + long double Tmin, Tmax, Tmin_satV; + HEOS.calc_Tmin_sat(Tmin, Tmin_satV); + Tmax = HEOS.calc_Tmax_sat(); + // Ancillary is deltas = s - hs_anchor.s + // First try a conventional call + try{ + T = anc.invert(specified_value - hs_anchor.smolar, Tmin, Tmax); + } + catch(std::exception &e){ + try{ + T = anc.invert(specified_value - hs_anchor.smolar, Tmin - 3, Tmax + 3); + } + catch(std::exception &e){ + double vmin = anc.evaluate(Tmin); + double vmax = anc.evaluate(Tmax); + if (std::abs(specified_value - hs_anchor.smolar) < std::abs(vmax)){ + T = Tmax - 0.1; + } + else{ + throw ValueError(format("Unable to invert ancillary equation for sL for value %Lg with Tminval %g and Tmaxval %g ", specified_value - hs_anchor.smolar, vmin, vmax)); + } + } + } } - catch(std::exception &e){ - double vmin = anc.evaluate(anc.get_Tmin()); - double vmax = anc.evaluate(anc.get_Tmax()); - throw ValueError(format("Unable to invert ancillary equation for sL for value %Lg with Tminval %g and Tmaxval %g ", specified_value - hs_anchor.smolar, vmin, vmax)); - } } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SV) { CoolPropFluid &component = *(HEOS.get_components()[0]); CoolProp::SimpleState crit = HEOS.get_state("reducing"); + CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor"); class Residual : public FuncWrapper1D { public: @@ -267,26 +294,33 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l } double call(double T){ long double s_liq = component->ancillaries.sL.evaluate(T) + component->pEOS->hs_anchor.smolar; - return s_liq + component->ancillaries.sLV.evaluate(T) - s; + long double resid = s_liq + component->ancillaries.sLV.evaluate(T) - s; + + return resid; }; }; Residual resid(component, HEOS.smolar()); - // If near the critical point, use a near critical guess value for T - if (std::abs(HEOS.smolar() - crit.smolar) < std::abs(component.ancillaries.sL.get_max_abs_error() + component.ancillaries.sLV.get_max_abs_error())) - { - T = std::max(0.99*crit.T, crit.T-0.1); + // Ancillary is deltas = s - hs_anchor.s + std::string errstr; + long double Tmin_satL, Tmin_satV; + HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV); + double Tmin = Tmin_satL; + double Tmax = HEOS.calc_Tmax_sat(); + try{ + T = Brent(resid, Tmin-3, Tmax, DBL_EPSILON, 1e-10, 50, errstr); } - else{ - // Ancillary is deltas = s - hs_anchor.s - std::string errstr; - long double Tmin_satL, Tmin_satV; - HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV); - double Tmin = Tmin_satL; - double Tmax = HEOS.calc_Tmax_sat(); - try{ T = Brent(resid, Tmin-3, Tmax+1, DBL_EPSILON, 1e-10, 50, errstr); } - catch(std::exception &e){ - throw ValueError(format("Unable to invert ancillary equation for sV: %s",e.what())); + catch(std::exception &e){ + long double vmax = resid.call(Tmax); + // If near the critical point, use a near critical guess value for T + if (std::abs(specified_value - hs_anchor.smolar) < std::abs(vmax)){ + T = std::max(0.99*crit.T, crit.T-0.1); + } + else{ + shared_ptr HEOS_copy(new HelmholtzEOSMixtureBackend(HEOS.get_components())); + HEOS_copy->update(QT_INPUTS, 1, Tmin); double sTmin = HEOS_copy->smolar(); + HEOS_copy->update(QT_INPUTS, 1, Tmax); double sTmax = HEOS_copy->smolar(); + T = (Tmax-Tmin)/(sTmax-sTmin)*(HEOS.smolar()-sTmin) + Tmin; } } } @@ -294,7 +328,9 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l { throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid",options.specified_variable)); } - if (T > HEOS.T_critical()-1){ T -= 1; } + // If T from the ancillaries is above the critical temp, this will cause failure + // in ancillaries for rhoV and rhoL, decrease if needed + T = std::min(T, static_cast(HEOS.T_critical()-0.1)); // Evaluate densities from the ancillary equations rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T); @@ -306,7 +342,10 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l // be way off, and often times negative SatL->update(DmolarT_INPUTS, rhoL, T); SatV->update(DmolarT_INPUTS, rhoV, T); - rhoL += -(SatL->p()-SatV->p())/SatL->first_partial_deriv(iP, iDmolar, iT); + double rhoL_updated = rhoL -(SatL->p()-SatV->p())/SatL->first_partial_deriv(iP, iDmolar, iT); + + // Accept the update if the liquid density is greater than the vapor density + if (rhoL_updated > rhoV){ rhoL = rhoL_updated; } // Update the state again with the better guess for the liquid density SatL->update(DmolarT_INPUTS, rhoL, T); @@ -480,16 +519,25 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l } v = linsolve(J, negativer); - - tau += options.omega*v[0]; - - if (options.use_logdelta){ - deltaL = exp(log(deltaL)+options.omega*v[1]); - deltaV = exp(log(deltaV)+options.omega*v[2]); - } - else{ - deltaL += options.omega*v[1]; - deltaV += options.omega*v[2]; + + // Conditions for an acceptable step are: + // a) tau > 1 + // b) rhoL > rhoV or deltaL > deltaV + double tau0 = tau, deltaL0 = deltaL, deltaV0 = deltaV; + for (double omega_local = 1.0; omega_local > 0.1; omega_local /= 1.1) + { + tau = tau0 + omega_local*options.omega*v[0]; + if (options.use_logdelta){ + deltaL = exp(log(deltaL0)+omega_local*options.omega*v[1]); + deltaV = exp(log(deltaV0)+omega_local*options.omega*v[2]); + } + else{ + deltaL = deltaL0 + omega_local*options.omega*v[1]; + deltaV = deltaV0 + omega_local*options.omega*v[2]; + } + if (tau > 1 && deltaL > deltaV){ + break; + } } rhoL = deltaL*reduce.rhomolar; @@ -808,13 +856,21 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HE stepL = options.omega/DELTA*( (KV-KL)*dJV-(JV-JL)*dKV); stepV = options.omega/DELTA*( (KV-KL)*dJL-(JV-JL)*dKL); - if (deltaL+stepL > 1 && deltaV+stepV < 1 && deltaV+stepV > 0){ - deltaL += stepL; deltaV += stepV; - rhoL = deltaL*reduce.rhomolar; rhoV = deltaV*reduce.rhomolar; - } - else{ - throw ValueError(format("rhosatPure_Akasaka densities are crossed")); + long double deltaL0 = deltaL, deltaV0 = deltaV; + // Conditions for an acceptable step are: + // a) rhoL > rhoV or deltaL > deltaV + for (double omega_local = 1.0; omega_local > 0.1; omega_local /= 1.1) + { + deltaL = deltaL0 + omega_local*stepL; + deltaV = deltaV0 + omega_local*stepV; + + if (deltaL > 1 && deltaV < 1 && deltaV > 0){ + break; + } } + + rhoL = deltaL*reduce.rhomolar; + rhoV = deltaV*reduce.rhomolar; iter++; if (iter > 100){ throw SolutionError(format("Akasaka solver did not converge after 100 iterations")); diff --git a/src/DataStructures.cpp b/src/DataStructures.cpp index 8acf0026..5ac0944d 100644 --- a/src/DataStructures.cpp +++ b/src/DataStructures.cpp @@ -372,8 +372,10 @@ public: input_pair_info input_pair_list[] = { input_pair_info(QT_INPUTS,"QT_INPUTS","Molar quality, Temperature in K"), - input_pair_info(QS_INPUTS,"QS_INPUTS","Molar quality, Entropy in J/mol/K"), - input_pair_info(HQ_INPUTS,"HQ_INPUTS","Enthalpy in J/mol, Molar quality"), + input_pair_info(QSmolar_INPUTS,"QS_INPUTS","Molar quality, Entropy in J/mol/K"), + input_pair_info(QSmass_INPUTS,"QS_INPUTS","Molar quality, Entropy in J/kg/K"), + input_pair_info(HmolarQ_INPUTS,"HQ_INPUTS","Enthalpy in J/mol, Molar quality"), + input_pair_info(HmassQ_INPUTS,"HQ_INPUTS","Enthalpy in J/kg, Molar quality"), input_pair_info(PQ_INPUTS,"PQ_INPUTS","Pressure in Pa, Molar quality"), input_pair_info(PT_INPUTS, "PT_INPUTS","Pressure in Pa, Temperature in K"), diff --git a/src/Tests/CoolProp-Tests.cpp b/src/Tests/CoolProp-Tests.cpp index b25e1614..908af647 100644 --- a/src/Tests/CoolProp-Tests.cpp +++ b/src/Tests/CoolProp-Tests.cpp @@ -1264,6 +1264,33 @@ TEST_CASE("Triple point checks", "[triple_point]") } } +TEST_CASE("Test that saturation solvers solve all the way to T = Tc", "[sat_T_to_Tc]") +{ + std::vector fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),','); + for (std::size_t i = 0; i < fluids.size(); ++i) + { + double Tc = Props1SI(fluids[i], "Tcrit"); + std::ostringstream ss1; + ss1 << "Check sat_T at Tc for " << fluids[i]; + SECTION(ss1.str(),"") + { + double pc = PropsSI("P","T",Tc,"Q",0,fluids[i]); + CAPTURE(pc); + CHECK(ValidNumber(pc)); + } + for (double j = 0.1; j > 1e-10; j /= 10) + { + std::ostringstream ss2; + ss2 << "Check sat_T for " << fluids[i] << " for Tc - T = " << j << " K"; + SECTION(ss2.str(),"") + { + CHECK(ValidNumber(PropsSI("D","T",Tc-j,"Q",0,fluids[i]))); + } + } + + } +} + TEST_CASE("Test that reference states are correct", "[reference_states]") { std::vector fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),',');