From 5f723ce5855e05294d670dfcf2d0484da58484d3 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Mon, 15 Sep 2014 14:19:22 +0200 Subject: [PATCH] Added refinement of the phase envelope, PQ and TQ work for saturation See https://github.com/CoolProp/CoolProp/issues/133 Closes https://github.com/CoolProp/CoolProp/issues/143 Signed-off-by: Ian Bell --- src/Backends/Helmholtz/FlashRoutines.cpp | 15 +++- .../Helmholtz/PhaseEnvelopeRoutines.cpp | 68 +++++++++++++++++-- .../Helmholtz/PhaseEnvelopeRoutines.h | 12 ++-- src/Backends/Helmholtz/VLERoutines.h | 14 ++-- 4 files changed, 89 insertions(+), 20 deletions(-) diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index 7f52c331..d549142c 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -427,6 +427,10 @@ void FlashRoutines::PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS, parame std::size_t &imax = solutions[0]; + // Shift the solution if needed to ensure that imax+2 and imax-1 are both in range + if (imax+2 >= env.T.size()){ imax--; } + else if (imax-1 < 0){ imax++; } + SaturationSolvers::newton_raphson_saturation NR; SaturationSolvers::newton_raphson_saturation_options IO; @@ -435,13 +439,16 @@ void FlashRoutines::PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS, parame IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::P_IMPOSED; // p -> rhomolar_vap IO.rhomolar_vap = CubicInterp(env.p, env.rhomolar_vap, imax-1, imax, imax+1, imax+2, static_cast(IO.p)); + IO.T = CubicInterp(env.rhomolar_vap, env.T, imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); } else if (other == iT){ IO.T = HEOS._T; IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::T_IMPOSED; - // p -> rhomolar_vap + // T -> rhomolar_vap IO.rhomolar_vap = CubicInterp(env.T, env.rhomolar_vap, imax-1, imax, imax+1, imax+2, static_cast(IO.T)); + IO.p = CubicInterp(env.rhomolar_vap, env.p, imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); } + IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); if (quality == SATURATED_VAPOR){ IO.bubble_point = false; @@ -458,9 +465,13 @@ void FlashRoutines::PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS, parame IO.bubble_point = true; IO.x = HEOS.get_mole_fractions(); // Because Q = 0 IO.y.resize(IO.x.size()); + // Phases are inverted, so "liquid" is actually the lighter phase + std::swap(IO.rhomolar_liq, IO.rhomolar_vap); for (std::size_t i = 0; i < IO.y.size()-1; ++i) // First N-1 elements { - IO.y[i] = CubicInterp(env.rhomolar_vap, env.x[i], imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); + // Phases are inverted, so liquid mole fraction (x) of phase envelope is actually the vapor phase mole fraction + // Use the liquid density as well + IO.y[i] = CubicInterp(env.rhomolar_vap, env.x[i], imax-1, imax, imax+1, imax+2, IO.rhomolar_liq); } IO.y[IO.y.size()-1] = 1 - std::accumulate(IO.y.begin(), IO.y.end()-1, 0.0); NR.call(HEOS, IO.x, IO.y, IO); diff --git a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp index fda697d6..f9a82963 100644 --- a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp +++ b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp @@ -145,8 +145,8 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) } } catch(std::exception &e){ - std::cout << e.what() << std::endl; - std::cout << IO.T << " " << IO.p << std::endl; + //std::cout << e.what() << std::endl; + //std::cout << IO.T << " " << IO.p << std::endl; // Try again, but with a smaller step IO.rhomolar_vap /= factor; factor = 1 + (factor-1)/2; @@ -199,14 +199,73 @@ 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]){ env.built = true; std::cout << format("envelope built.\n"); return; } + if (iter > 4 && IO.p < env.p[0]){ + env.built = true; + std::cout << format("envelope built.\n"); + + // Now we refine the phase envelope to add some points in places that are still pretty rough + // Can be re-enabled by just uncommenting the following line + refine(HEOS); + + return; + } // Reset the failure counter failure_count = 0; } - } +void PhaseEnvelopeRoutines::refine(HelmholtzEOSMixtureBackend &HEOS) +{ + PhaseEnvelopeData &env = HEOS.PhaseEnvelope; + SaturationSolvers::newton_raphson_saturation NR; + SaturationSolvers::newton_raphson_saturation_options IO; + IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::RHOV_IMPOSED; + IO.bubble_point = false; + IO.y = HEOS.get_mole_fractions(); + + std::size_t i = 0; + do{ + // Don't do anything if change in density is small enough + if (std::abs(env.rhomolar_vap[i]/env.rhomolar_vap[i+1]-1) < 0.2){ i++; continue; } + + double rhomolar_vap_start = env.rhomolar_vap[i], + rhomolar_vap_end = env.rhomolar_vap[i+1]; + + // Ok, now we are going to do some more refining in this step + for (double rhomolar_vap = rhomolar_vap_start*1.1; rhomolar_vap < rhomolar_vap_end; rhomolar_vap *= 1.1) + { + IO.rhomolar_vap = rhomolar_vap; + IO.x.resize(IO.y.size()); + if (i < env.T.size()-3){ + IO.T = CubicInterp(env.rhomolar_vap, env.T, i, i+1, i+2, i+3, IO.rhomolar_vap); + IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, i, i+1, i+2, i+3, IO.rhomolar_vap); + for (std::size_t j = 0; j < IO.x.size()-1; ++j){ // First N-1 elements + IO.x[j] = CubicInterp(env.rhomolar_vap, env.x[j], i, i+1, i+2, i+3, IO.rhomolar_vap); + } + } + else{ + IO.T = CubicInterp(env.rhomolar_vap, env.T, i, i-1, i-2, i-3, IO.rhomolar_vap); + IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, i, i-1, i-2, i-3, IO.rhomolar_vap); + for (std::size_t j = 0; j < IO.x.size()-1; ++j){ // First N-1 elements + IO.x[j] = CubicInterp(env.rhomolar_vap, env.x[j], i, i-1, i-2, i-3, IO.rhomolar_vap); + } + } + IO.x[IO.x.size()-1] = 1 - std::accumulate(IO.x.begin(), IO.x.end()-1, 0.0); + try{ + NR.call(HEOS, IO.y, IO.x, IO); + env.insert_variables(IO.T, IO.p, IO.rhomolar_liq, IO.rhomolar_vap, IO.hmolar_liq, + IO.hmolar_vap, IO.smolar_liq, IO.smolar_vap, IO.x, IO.y, i+1); + std::cout << "dv " << IO.rhomolar_vap << " dl " << IO.rhomolar_liq << " T " << IO.T << " p " << IO.p << " hl " << IO.hmolar_liq << " hv " << IO.hmolar_vap << " sl " << IO.smolar_liq << " sv " << IO.smolar_vap << " x " << vec_to_string(IO.x, "%0.10Lg") << " Ns " << IO.Nsteps << std::endl; + } + catch(std::exception &e){ + continue; + } + i++; + } + } + while (i < env.T.size()-1); +} void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS) { enum maxima_points {PMAX_SAT = 0, TMAX_SAT = 1}; @@ -354,7 +413,6 @@ std::vector > PhaseEnvelopeRoutines::find_in } if (matched){ - std::cout << i << " " << i+1 << std::endl; intersections.push_back(std::pair(i, i+1)); } } diff --git a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h index f1345cb2..b987b88d 100644 --- a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h +++ b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h @@ -10,18 +10,18 @@ class PhaseEnvelopeRoutines{ */ static void build(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 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 diff --git a/src/Backends/Helmholtz/VLERoutines.h b/src/Backends/Helmholtz/VLERoutines.h index d8d4d04b..0117cbce 100644 --- a/src/Backends/Helmholtz/VLERoutines.h +++ b/src/Backends/Helmholtz/VLERoutines.h @@ -347,7 +347,7 @@ namespace SaturationSolvers rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE; }; - /** Call the Newton-Raphson VLE Solver + /** \brief Call the Newton-Raphson VLE Solver * * This solver must be passed reasonable guess values for the mole fractions, * densities, etc. You may want to take a few steps of successive substitution @@ -360,15 +360,15 @@ namespace SaturationSolvers */ void call(HelmholtzEOSMixtureBackend &HEOS, const std::vector &z, std::vector &z_incipient, newton_raphson_saturation_options &IO); - /*! Build the arrays for the Newton-Raphson solve - - This method builds the Jacobian matrix, the sensitivity matrix, etc. + /** \brief Build the arrays for the Newton-Raphson solve * - */ + * This method builds the Jacobian matrix, the sensitivity matrix, etc. + * + */ void build_arrays(); - /** Check the derivatives in the Jacobian using numerical derivatives. - */ + /** \brief Check the derivatives in the Jacobian using numerical derivatives. + */ void check_Jacobian(); }; };