From 2376baffddfc5fbb843fd4e5a342827ac8ac7700 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Thu, 7 Apr 2016 21:06:47 -0600 Subject: [PATCH] Finish stability calculations for critical point calculations! --- .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 3 +-- src/Backends/Helmholtz/VLERoutines.cpp | 23 ++++++++++++------- src/Backends/Helmholtz/VLERoutines.h | 5 ++-- 3 files changed, 19 insertions(+), 12 deletions(-) diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 5f246b47..9a8db743 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -3374,8 +3374,7 @@ double HelmholtzEOSMixtureBackend::calc_tangent_plane_distance(const double T, c throw ValueError(format("Trial composition vector size [%d] is not the same as bulk composition [%d]", w.size(), z.size())); } update_TPD_state(); - TPD_state->set_mole_fractions(std::vector(w.begin(), w.end())); - TPD_state->specify_phase(iphase_liquid); // Something homogeneous + TPD_state->set_mole_fractions(w); if (rhomolar_guess < 0){ TPD_state->update(PT_INPUTS, p, T); } diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index 224b4150..d45a2d64 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -1722,6 +1722,8 @@ void StabilityRoutines::StabilityEvaluationClass::trial_compositions(){ g0 += z[i]*(K[i]-1); // The summation for beta = 0 g1 += z[i]*(1-1/K[i]); // The summation for beta = 1 } + // Copy K-factors for later use + K0 = K; // Now see what to do about the g(0) and g(1) values // ----- // @@ -1740,16 +1742,15 @@ void StabilityRoutines::StabilityEvaluationClass::trial_compositions(){ SaturationSolvers::x_and_y_from_K(beta, K, z, x, y); normalize_vector(x); normalize_vector(y); + if (debug){ fmt::printf("1) T: %g p: %g beta: %g\n", HEOS.T(), HEOS.p(), beta); } } void StabilityRoutines::StabilityEvaluationClass::successive_substitution(int num_steps){ // ---- // Do a few steps of successive substitution // ---- - HEOS.SatL->set_mole_fractions(x); - HEOS.SatV->set_mole_fractions(y); - HEOS.SatL->calc_reducing_state(); - HEOS.SatV->calc_reducing_state(); + HEOS.SatL->set_mole_fractions(x); HEOS.SatL->calc_reducing_state(); + HEOS.SatV->set_mole_fractions(y); HEOS.SatV->calc_reducing_state(); // First use cubic as a guess for the density of liquid and vapor phases rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(HEOS.T(), HEOS.p(), iphase_liquid); // [mol/m^3] @@ -1770,6 +1771,7 @@ void StabilityRoutines::StabilityEvaluationClass::successive_substitution(int nu rhomolar_liq = 1/(v_SRK - summer_c); } + if (debug){ fmt::printf("2) SS1: i beta K rho' rho''\n"); } for (int step_count = 0; step_count < num_steps; ++step_count){ // Set the composition HEOS.SatL->set_mole_fractions(x); HEOS.SatV->set_mole_fractions(y); @@ -1805,6 +1807,7 @@ void StabilityRoutines::StabilityEvaluationClass::successive_substitution(int nu SaturationSolvers::x_and_y_from_K(beta, K, z, x, y); normalize_vector(x); normalize_vector(y); + if (debug){ fmt::printf("2) %d %g %s %g %g\n", step_count, beta, vec_to_string(K, "%0.6f"), rhomolar_liq, rhomolar_vap); } } } void StabilityRoutines::StabilityEvaluationClass::check_stability(){ @@ -1815,6 +1818,7 @@ void StabilityRoutines::StabilityEvaluationClass::check_stability(){ this->tpd_liq = HEOS.tangent_plane_distance(HEOS.T(), HEOS.p(), x, rhomolar_liq); this->tpd_vap = HEOS.tangent_plane_distance(HEOS.T(), HEOS.p(), y, rhomolar_vap); this->DELTAG_nRT = (1-beta)*tpd_liq + beta*(tpd_vap); + if (debug){ fmt::printf("3) tpd': %g tpd'': %g DELTAG/nRT: %g\n", tpd_liq, tpd_vap, DELTAG_nRT); } // If any of these cases are met, conclusively unstable, stop! if (this->tpd_liq < 0 || this->tpd_vap < 0 || this->DELTAG_nRT < 0){ @@ -1827,13 +1831,14 @@ void StabilityRoutines::StabilityEvaluationClass::check_stability(){ // Generate light and heavy test compositions (Gernert, 2014, Eq. 23) std::vector xL(z.size()), xH(z.size()); for (std::size_t i = 0; i < z.size(); ++i){ - xL[i] = z[i]*K[i]; - xH[i] = z[i]/K[i]; + xL[i] = z[i]*K0[i]; + xH[i] = z[i]/K0[i]; } normalize_vector(xL); normalize_vector(xH); // For each composition, use successive substitution to try to evaluate stability + if (debug){ fmt::printf("3) SS2: i x' x'' rho' rho'' tpd' tpd''\n"); } for (int step_count = 0; step_count < 100; ++step_count){ // Set the composition HEOS.SatL->set_mole_fractions(xH); HEOS.SatV->set_mole_fractions(xL); @@ -1849,16 +1854,18 @@ void StabilityRoutines::StabilityEvaluationClass::check_stability(){ double tpd_L = HEOS.tangent_plane_distance(HEOS.T(), HEOS.p(), xL, rhomolar_vap); double tpd_H = HEOS.tangent_plane_distance(HEOS.T(), HEOS.p(), xH, rhomolar_liq); tpdL.push_back(tpd_L); tpdH.push_back(tpd_H); - if (tpd_L < 0 || tpd_H < 0){_stable = false; return;} // Calculate the new composition from the fugacity coefficients for (std::size_t i = 0; i < z.size(); ++i){ - // TODO: this seems weird; what fugacity coefficient is meant to be used in numerator of Gernert, 2014, Eq. 24 ? xL[i] = z[i]*HEOS.fugacity_coefficient(i)/HEOS.SatV->fugacity_coefficient(i); xH[i] = z[i]*HEOS.fugacity_coefficient(i)/HEOS.SatL->fugacity_coefficient(i); } normalize_vector(xL); normalize_vector(xH); + if (debug){ fmt::printf("2) %d %s %s %g %g %g %g\n", vec_to_string(xL, "%0.6f"), vec_to_string(xH, "%0.6f"), rhomolar_liq, rhomolar_vap, tpd_L, tpd_H); } + + if (tpd_L < 0 || tpd_H < 0){_stable = false; return;} + } } diff --git a/src/Backends/Helmholtz/VLERoutines.h b/src/Backends/Helmholtz/VLERoutines.h index a803b5ee..4e877c8d 100644 --- a/src/Backends/Helmholtz/VLERoutines.h +++ b/src/Backends/Helmholtz/VLERoutines.h @@ -441,14 +441,15 @@ namespace StabilityRoutines{ class StabilityEvaluationClass{ protected: HelmholtzEOSMixtureBackend &HEOS; - std::vector lnK, K, x, y; + std::vector lnK, K, K0, x, y; const std::vector &z; double rhomolar_liq, rhomolar_vap, beta, tpd_liq, tpd_vap, DELTAG_nRT; private: bool _stable; + bool debug; public: StabilityEvaluationClass(HelmholtzEOSMixtureBackend &HEOS) - : HEOS(HEOS), z(HEOS.get_mole_fractions_doubleref()), rhomolar_liq(-1), rhomolar_vap(-1), beta(-1), tpd_liq(10000), tpd_vap(100000), DELTAG_nRT(10000), _stable(false) {}; + : HEOS(HEOS), z(HEOS.get_mole_fractions_doubleref()), rhomolar_liq(-1), rhomolar_vap(-1), beta(-1), tpd_liq(10000), tpd_vap(100000), DELTAG_nRT(10000), _stable(false), debug(false) {}; /** \brief Calculate trial compositions */ void trial_compositions();