diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index f4d3844e..6ef32277 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -1,6 +1,7 @@ #include "VLERoutines.h" #include "FlashRoutines.h" #include "HelmholtzEOSMixtureBackend.h" +#include "PhaseEnvelopeRoutines.h" namespace CoolProp{ @@ -29,8 +30,45 @@ template T dgdbeta_RachfordRice(const std::vector &z, const std::vec void FlashRoutines::PT_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS) { - if (false){//HEOS.PhaseEnvelope.built){ + if (HEOS.PhaseEnvelope.built){ // Use the phase envelope if already constructed to determine phase boundary + // Determine whether you are inside (two-phase) or outside (single-phase) + SimpleState closest_state; + std::size_t i; + bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS, iP, HEOS._p, iT, HEOS._T, i, closest_state); + if (!twophase && HEOS._T > closest_state.T){ + // Gas solution - bounded between phase envelope temperature and very high temperature + // + // Start with a guess value from SRK + long double rhomolar_guess = HEOS.solver_rho_Tp_SRK(HEOS._T, HEOS._p, iphase_gas); + + solver_TP_resid resid(HEOS, HEOS._T, HEOS._p); + std::string errstr; + HEOS.specify_phase(iphase_gas); + try{ + // Try using Newton's method + long double rhomolar = Newton(resid, rhomolar_guess, 1e-10, 100, errstr); + // Make sure the solution is within the bounds + if (!is_in_closed_range(static_cast(closest_state.rhomolar), 0.0L, rhomolar)){ + throw ValueError("out of range"); + } + HEOS.update_DmolarT_direct(rhomolar, HEOS._T); + } + catch(std::exception &e){ + // If that fails, try a bounded solver + long double rhomolar = Brent(resid, closest_state.rhomolar, 1e-10, DBL_EPSILON, 1e-10, 100, errstr); + // Make sure the solution is within the bounds + if (!is_in_closed_range(static_cast(closest_state.rhomolar), 0.0L, rhomolar)){ + throw ValueError("out of range"); + } + } + HEOS.unspecify_phase(); + + } + else{ + // Liquid solution + throw ValueError(); + } } else{ // Following the strategy of Gernert, 2014 @@ -643,7 +681,6 @@ void FlashRoutines::HSU_P_flash_singlephase_Brent(HelmholtzEOSMixtureBackend &HE // Get the value of the desired variable eos = HEOS->keyed_output(other); pp = HEOS->p(); - // Difference between the two is to be driven to zero r = eos - value; @@ -694,6 +731,7 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth { bool saturation_called = false; long double value; + std::cout << "PHSU" << std::endl; if (HEOS.imposed_phase_index != iphase_not_imposed) { // Use the phase defined by the imposed phase @@ -701,21 +739,25 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth } else { + std::cout << "PHSU not imposed" << std::endl; + // Find the phase, while updating all internal variables possible + switch (other) + { + case iSmolar: + value = HEOS.smolar(); break; + case iHmolar: + value = HEOS.hmolar(); break; + case iUmolar: + value = HEOS.umolar(); break; + default: + throw ValueError(format("Input for other [%s] is invalid", get_parameter_information(other, "long").c_str())); + } if (HEOS.is_pure_or_pseudopure) { // Find the phase, while updating all internal variables possible - switch (other) - { - case iSmolar: - value = HEOS.smolar(); HEOS.p_phase_determination_pure_or_pseudopure(iSmolar, value, saturation_called); break; - case iHmolar: - value = HEOS.hmolar(); HEOS.p_phase_determination_pure_or_pseudopure(iHmolar, value, saturation_called); break; - case iUmolar: - value = HEOS.umolar(); HEOS.p_phase_determination_pure_or_pseudopure(iUmolar, value, saturation_called); break; - default: - throw ValueError(format("Input for other [%s] is invalid", get_parameter_information(other, "long").c_str())); - } + HEOS.p_phase_determination_pure_or_pseudopure(other, value, saturation_called); + if (HEOS.isHomogeneousPhase()) { // Now we use the single-phase solver to find T,rho given P,Y using a @@ -765,7 +807,30 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth } else { - throw NotImplementedError("HSU_P_flash does not support mixtures (yet)"); + std::cout << format("PHSU flash for mixture\n"); + if (HEOS.PhaseEnvelope.built){ + // Determine whether you are inside or outside + SimpleState closest_state; + std::size_t iclosest; + std::cout << format("pre is inside\n"); + bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS, iP, HEOS._p, other, value, iclosest, closest_state); + std::cout << format("post is inside\n"); + + std::string errstr; + if (!twophase){ + PY_singlephase_flash_resid resid(HEOS, HEOS._p, other, value); + // If that fails, try a bounded solver + long double rhomolar = Brent(resid, closest_state.T+10, 1000, DBL_EPSILON, 1e-10, 100, errstr); + HEOS.unspecify_phase(); + } + else{ + throw ValueError("two-phase solution for Y"); + } + + } + else{ + throw ValueError("phase envelope must be built"); + } } } } diff --git a/src/Backends/Helmholtz/FlashRoutines.h b/src/Backends/Helmholtz/FlashRoutines.h index fb95b67f..c8d6ddc7 100644 --- a/src/Backends/Helmholtz/FlashRoutines.h +++ b/src/Backends/Helmholtz/FlashRoutines.h @@ -13,6 +13,7 @@ state is based. #define FLASHROUTINES_H #include "HelmholtzEOSMixtureBackend.h" +#include "Solvers.h" namespace CoolProp{ @@ -74,5 +75,87 @@ public: static void PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, parameters other); }; + +/** A residual function for the rho(T,P) solver + */ +class solver_TP_resid : public FuncWrapper1D +{ +public: + long double T, p, r, peos, rhomolar, rhor, tau, R_u, delta, dalphar_dDelta; + 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_state().rhomolar; + this->tau = HEOS.get_reducing_state().T/T; this->R_u = HEOS.gas_constant(); + }; + double call(double rhomolar){ + this->rhomolar = rhomolar; + delta = rhomolar/rhor; // needed for derivative + HEOS->update_DmolarT_direct(rhomolar, T); + peos = HEOS->p(); + r = (peos-p)/p; + return r; + }; + double deriv(double rhomolar){ + // dp/drho|T / pspecified + return R_u*T*(1+2*delta*HEOS->dalphar_dDelta()+pow(delta, 2)*HEOS->d2alphar_dDelta2())/p; + }; +}; + +/** A residual function for the f(P, Y) solver + */ +class PY_singlephase_flash_resid : public FuncWrapper1D +{ +public: + + HelmholtzEOSMixtureBackend *HEOS; + long double p; + parameters other; + long double r, eos, value, T, rhomolar; + + int iter; + long double r0, r1, T1, T0, eos0, eos1, pp; + PY_singlephase_flash_resid(HelmholtzEOSMixtureBackend &HEOS, long double p, parameters other, long double value) : + HEOS(&HEOS), p(p), other(other), value(value) + { + iter = 0; + // Specify the state to avoid saturation calls, but only if phase is subcritical + if (HEOS.phase() == iphase_liquid || HEOS.phase() == iphase_gas ){ + HEOS.specify_phase(HEOS.phase()); + } + }; + double call(double T){ + + this->T = T; + + // Run the solver with T,P as inputs; + HEOS->update(PT_INPUTS, p, T); + + rhomolar = HEOS->rhomolar(); + HEOS->update(DmolarT_INPUTS, rhomolar, T); + // Get the value of the desired variable + eos = HEOS->keyed_output(other); + pp = HEOS->p(); + + // Difference between the two is to be driven to zero + r = eos - value; + + // Store values for later use if there are errors + if (iter == 0){ + r0 = r; T0 = T; eos0 = eos; + } + else if (iter == 1){ + r1 = r; T1 = T; eos1 = eos; + } + else{ + r0 = r1; T0 = T1; eos0 = eos1; + r1 = r; T1 = T; eos1 = eos; + } + + iter++; + return r; + }; +}; + } /* namespace CoolProp */ #endif /* FLASHROUTINES_H */ diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index fa195ac2..1460683b 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -1696,11 +1696,19 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double } try{ + // First we try with Newton's method with analytic derivative double rhomolar = Newton(resid, rhomolar_guess, 1e-8, 100, errstring); if (!ValidNumber(rhomolar)){ throw ValueError(); } +// double rr = this->rhomolar(); +// for (double rho = rr*0.1; rho < 1.1*rr; rho += 100){ +// specify_phase(iphase_gas); +// this->update(DmolarT_INPUTS, rho, T); +// unspecify_phase(); +// std::cout << format("%g %g\n", rho, this->p()); +// } return rhomolar; } catch(std::exception &) @@ -1710,6 +1718,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double double rhomolar = Secant(resid, rhomolar_guess, 1.1*rhomolar_guess, 1e-8, 100, errstring); if (!ValidNumber(rhomolar)){throw ValueError();} return rhomolar; + } catch(std::exception &) { diff --git a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp index 4ecb3e36..fda697d6 100644 --- a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp +++ b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp @@ -52,7 +52,9 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) IO.x[0] = 0.6689704673; IO.x[1] = 0.3310295327; */ - IO.rhomolar_liq *= 1.2; + + //IO.rhomolar_liq *= 1.2; + IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::P_IMPOSED; NR.call(HEOS, IO.y, IO.x, IO); @@ -197,11 +199,12 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) factor = std::max(factor, static_cast(1.01)); // Stop if the pressure is below the starting pressure - if (iter > 4 && IO.p < env.p[0]){ return; } + if (iter > 4 && IO.p < env.p[0]){ env.built = true; std::cout << format("envelope built.\n"); return; } // Reset the failure counter failure_count = 0; } + } void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS) @@ -282,6 +285,7 @@ void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS) }; double call(double rhomolar_vap){ PhaseEnvelopeData &env = HEOS->PhaseEnvelope; + IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::RHOV_IMPOSED; IO.bubble_point = false; IO.rhomolar_vap = rhomolar_vap; IO.y = HEOS->get_mole_fractions(); @@ -329,6 +333,109 @@ void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS) } } +std::vector > PhaseEnvelopeRoutines::find_intersections(HelmholtzEOSMixtureBackend &HEOS, parameters iInput, long double value) +{ + std::vector > intersections; + + PhaseEnvelopeData &env = HEOS.PhaseEnvelope; + for (std::size_t i = 0; i < env.p.size()-1; ++i){ + bool matched = false; + switch(iInput){ + case iP: + if (is_in_closed_range(env.p[i], env.p[i+1], value)){ matched = true; } break; + case iT: + if (is_in_closed_range(env.T[i], env.T[i+1], value)){ matched = true; } break; + case iHmolar: + if (is_in_closed_range(env.hmolar_vap[i], env.hmolar_vap[i+1], value)){ matched = true; } break; + case iSmolar: + if (is_in_closed_range(env.smolar_vap[i], env.smolar_vap[i+1], value)){ matched = true; } break; + default: + throw ValueError(format("bad index to find_intersections")); + } + + if (matched){ + std::cout << i << " " << i+1 << std::endl; + intersections.push_back(std::pair(i, i+1)); + } + } + return intersections; +} +bool PhaseEnvelopeRoutines::is_inside(HelmholtzEOSMixtureBackend &HEOS, parameters iInput1, long double value1, parameters iInput2, long double value2, std::size_t &iclosest, SimpleState &closest_state) +{ + PhaseEnvelopeData &env = HEOS.PhaseEnvelope; + std::vector > intersections = find_intersections(HEOS, iInput1, value1); + + // For now, first input must be p + if (iInput1 != iP){throw ValueError("For now, first input must be p in is_inside");} + + // If number of intersections is 0, input is out of range, quit + if (intersections.size() == 0){ throw ValueError("Input is out of range; no intersections found"); } + + // If number of intersections is 1, input will be determined based on the single intersection + // Need to know if values increase or decrease to the right of the intersection point + if (intersections.size()%2 != 0){ throw ValueError("Input is weird; odd number of intersections found"); } + + // If number of intersections is even, might be a bound + if (intersections.size()%2 == 0){ + if (intersections.size() != 2){throw ValueError("for now only even value accepted is 2"); } + std::vector other_indices(4, 0); + std::vector *y; + std::vector other_values(4, 0); + other_indices[0] = intersections[0].first; other_indices[1] = intersections[0].second; + other_indices[2] = intersections[1].first; other_indices[3] = intersections[1].second; + + switch(iInput2){ + case iT: y = &(env.T); break; + case iP: y = &(env.p); break; + case iHmolar: y = &(env.hmolar_vap); break; + case iSmolar: y = &(env.smolar_vap); break; + default: break; + } + + other_values[0] = (*y)[other_indices[0]]; other_values[1] = (*y)[other_indices[1]]; + other_values[2] = (*y)[other_indices[2]]; other_values[3] = (*y)[other_indices[3]]; + + long double min_other = *(std::min_element(other_values.begin(), other_values.end())); + long double max_other = *(std::max_element(other_values.begin(), other_values.end())); + + std::cout << format("min: %Lg max: %Lg val: %Lg\n", min_other, max_other, value2); + + // If by using the outer bounds of the second variable, we are outside the range, + // then the value is definitely not inside the phase envelope and we don't need to + // do any more analysis. + if (!is_in_closed_range(min_other, max_other, value2)){ + std::vector d(4, 0); + d[0] = std::abs(other_values[0]-value2); d[1] = std::abs(other_values[1]-value2); + d[2] = std::abs(other_values[2]-value2); d[3] = std::abs(other_values[3]-value2); + + // Index of minimum distance in the other_values vector + std::size_t idist = std::distance(d.begin(), std::min_element(d.begin(), d.end())); + // Index of closest point in the phase envelope + std::size_t iclosest = other_indices[idist]; + + // Get the state for the point which is closest to the desired value - this + // can be used as a bounding value in the outer single-phase flash routine + // since you know (100%) that it is a good bound + closest_state.T = env.T[iclosest]; + closest_state.p = env.p[iclosest]; + closest_state.rhomolar = env.rhomolar_vap[iclosest]; + closest_state.hmolar = env.hmolar_vap[iclosest]; + closest_state.smolar = env.smolar_vap[iclosest]; + closest_state.Q = env.Q[iclosest]; + + std::cout << format("it is not inside") << std::endl; + return false; + } + else{ + // Now we have to do a saturation flash call in order to determine whether or not we are inside the phase envelope or not + + // First we can interpolate using the phase envelope to get good guesses for the necessary values + + throw ValueError("For now can't be inside"); + } + } +} + } /* namespace CoolProp */ #endif \ No newline at end of file diff --git a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h index 8a1cc4b0..f1345cb2 100644 --- a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h +++ b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h @@ -3,12 +3,49 @@ namespace CoolProp{ class PhaseEnvelopeRoutines{ -public: + public: + /** \brief Build the phase envelope + * + * @param HEOS The HelmholtzEOSMixtureBackend instance to be used + */ static void build(HelmholtzEOSMixtureBackend &HEOS); /** \brief Finalize the phase envelope and calculate maxima values, critical point, etc. + * + * @param HEOS The HelmholtzEOSMixtureBackend instance to be used */ static void finalize(HelmholtzEOSMixtureBackend &HEOS); + + /** \brief Refine the phase envelope, adding points in places that are sparse + * + * @param HEOS The HelmholtzEOSMixtureBackend instance to be used + */ + //static void refine(HelmholtzEOSMixtureBackend &HEOS); + + /** \brief Determine which indices bound a given value + * + * If you provide pressure for instance, it will return each of the indices + * that bound crossings in the pressure versus rhov curve. Thus this information + * can be used to determine whether another input is "inside" or "outside" the phase + * boundary. + * + * @param HEOS The HelmholtzEOSMixtureBackend instance to be used + * @param iInput The key for the variable type that is to be checked + * @param value The value associated with iInput + */ + static std::vector > find_intersections(HelmholtzEOSMixtureBackend &HEOS, parameters iInput, long double value); + + /** \brief Determine whether a pair of inputs is inside or outside the phase envelope + * + * @param HEOS The HelmholtzEOSMixtureBackend instance to be used + * @param iInput1 The key for the first input + * @param value1 The value of the first input + * @param iInput2 The key for the second input + * @param value2 The value of the second input + * @param iclosest The index of the phase envelope for the closest point + * @param closest_state A SimpleState corresponding to the closest point found on the phase envelope + */ + static bool is_inside(HelmholtzEOSMixtureBackend &HEOS, parameters iInput1, long double value1, parameters iInput2, long double value2, std::size_t &iclosest, SimpleState &closest_state); }; } /* namespace CoolProp */ \ No newline at end of file diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index 75291cec..66d8902b 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -972,7 +972,7 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke throw ValueError(format("newton_raphson_saturation::call reached max number of iterations [%d]",IO.Nstep_max)); } } - while(this->error_rms > 1e-9 && min_rel_change > 10*DBL_EPSILON && iter < IO.Nstep_max); + while(this->error_rms > 1e-9 && min_rel_change > 1000*DBL_EPSILON && iter < IO.Nstep_max); IO.Nsteps = iter; IO.p = p; @@ -1002,7 +1002,6 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays() else{ // Liquid is incipient phase, set its composition rSatL.set_mole_fractions(x); - rSatV.set_mole_fractions(y); } if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){ @@ -1014,9 +1013,9 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays() rSatV.update_TP_guessrho(T, p, rhomolar_vap); } else{ - throw ValueError(); + throw ValueError("imposed variable not set for NR VLE"); } - + // For diagnostic purposes calculate the pressures (no derivatives are evaluated) long double p_liq = rSatL.p(); long double p_vap = rSatV.p(); @@ -1029,7 +1028,6 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays() x_N_dependency_flag xN_flag = XN_DEPENDENT; if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){ - x_N_dependency_flag xN_flag = XN_DEPENDENT; // For the residuals F_i (equality of fugacities) for (std::size_t i = 0; i < N; ++i) { @@ -1101,6 +1099,7 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays() error_rms = sqrt(error_rms); // Square-root (The R in RMS) // Calculate derivatives along phase boundary; + // Gernert thesis 3.96 and 3.97 long double dQ_dPsat = 0, dQ_dTsat = 0; for (std::size_t i = 0; i < N; ++i) { diff --git a/src/main.cxx b/src/main.cxx index d487c37e..7ebeebcb 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -423,9 +423,9 @@ int main() #endif #if 0 { - ::set_debug_level(11); + ::set_debug_level(0); std::vector tags; - tags.push_back("[REFPROP_backwards_compatibility]"); + tags.push_back("[mixture_derivs]"); run_user_defined_tests(tags); char c; std::cin >> c; @@ -440,8 +440,8 @@ int main() #endif #if 0 { - double TTT0 = PropsSI("T","P",1e6,"Q",1,"REFPROP::Ethane[0.5]&Propane[0.5]"); - double TTT1 = PropsSI("T","P",1e6,"Q",1,"HEOS::Ethane[0.5]&Propane[0.5]"); + double TTT0 = PropsSI("T","Q",1,"P",3e6,"REFPROP::R32[0.3]&R125[0.7]"); + double TTT1 = PropsSI("T","Q",1,"P",3e6,"HEOS::R125[0.7]&R32[0.3]"); int rr =0; } #endif @@ -450,8 +450,8 @@ int main() ::set_debug_level(0); - shared_ptr HEOS(AbstractState::factory("HEOS","Methane&Propane")); - std::vector z(2, 0.3); z[1] = 0.7; + shared_ptr HEOS(AbstractState::factory("HEOS","Methane&Ethane")); + std::vector z(2, 0.8); z[1] = 1-z[0]; //shared_ptr HEOS(AbstractState::factory("HEOS","Methane&Propane&Ethane&n-Butane")); //std::vector z(4, 0.1); z[1] = 0.35; z[2] = 0.35, z[3] = 0.2; HEOS->set_mole_fractions(z); @@ -465,6 +465,16 @@ int main() std::cout << get_global_param_string("errstring") << std::endl; } t2 = clock(); + HEOS->update(PSmolar_INPUTS, 4e6, 79.1048486373); + long double TT = HEOS->T(); + HEOS->update(PQ_INPUTS, 1.3e5, 1); + double ssat = HEOS->smolar(); + double hsat = HEOS->hmolar(); + double dsat = HEOS->rhomolar(); + for (long double s = ssat + 60; s > ssat; s -= 5){ + HEOS->update(PSmolar_INPUTS, 1.3e5, s); + std::cout << s << " " << HEOS->rhomolar() << " " << dsat << std::endl; + } std::cout << format("time: %g s/call\n", ((double)(t2-t1))/CLOCKS_PER_SEC); exit(EXIT_SUCCESS);