From 20cae192b08b74edcd44ecb09439431e628b4183 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Thu, 11 Sep 2014 02:15:19 +0200 Subject: [PATCH] Added necessary code to get sat_p working for dewpoint calculations, phase envelopes build for R32/R125 - was getting bad density roots Signed-off-by: Ian Bell --- include/PhaseEnvelope.h | 16 +- .../Helmholtz/HelmholtzEOSMixtureBackend.h | 3 + src/Backends/Helmholtz/MixtureDerivatives.cpp | 65 +++++++- src/Backends/Helmholtz/MixtureDerivatives.h | 4 + .../Helmholtz/PhaseEnvelopeRoutines.cpp | 44 ++++-- src/Backends/Helmholtz/VLERoutines.cpp | 146 ++++++++++++------ src/Backends/Helmholtz/VLERoutines.h | 9 +- 7 files changed, 222 insertions(+), 65 deletions(-) diff --git a/include/PhaseEnvelope.h b/include/PhaseEnvelope.h index e7f5bb68..b1103c01 100644 --- a/include/PhaseEnvelope.h +++ b/include/PhaseEnvelope.h @@ -18,7 +18,7 @@ public: icrit; ///< The index of the point corresponding to the critical point std::vector< std::vector > K, lnK, x, y; - std::vector T, p, lnT, lnp, rhomolar_liq, rhomolar_vap, lnrhomolar_liq, lnrhomolar_vap, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap; + std::vector T, p, lnT, lnp, rhomolar_liq, rhomolar_vap, lnrhomolar_liq, lnrhomolar_vap, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap, Q; PhaseEnvelopeData(){ built = false; TypeI = false; }; @@ -32,7 +32,7 @@ public: void clear(){ T.clear(); p.clear(); lnT.clear(); lnp.clear(); rhomolar_liq.clear(); rhomolar_vap.clear(); lnrhomolar_liq.clear(); lnrhomolar_vap.clear(); hmolar_liq.clear(); hmolar_vap.clear(); smolar_liq.clear(); smolar_vap.clear(); - K.clear(); lnK.clear(); x.clear(); y.clear(); + K.clear(); lnK.clear(); x.clear(); y.clear(); Q.clear(); } void insert_variables(const long double T, const long double p, @@ -67,6 +67,12 @@ public: this->x[j].insert(this->x[j].begin() + i, x[j]); this->y[j].insert(this->y[j].begin() + i, y[j]); } + if (rhomolar_liq > rhomolar_vap){ + this->Q.insert(this->Q.begin(), 1); + } + else{ + this->Q.insert(this->Q.begin(), 0); + } }; void store_variables(const long double T, const long double p, @@ -100,6 +106,12 @@ public: this->x[i].push_back(x[i]); this->y[i].push_back(y[i]); } + if (rhomolar_liq > rhomolar_vap){ + this->Q.push_back(1); + } + else{ + this->Q.push_back(0); + } }; }; diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index 7b13c73f..9df97cb0 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -64,6 +64,9 @@ public: std::vector &get_K(){return K;}; std::vector &get_lnK(){return lnK;}; + std::vector calc_mole_fractions_liquid(void){return SatL->get_mole_fractions();}; + std::vector calc_mole_fractions_vapor(void){return SatV->get_mole_fractions();}; + const CoolProp::PhaseEnvelopeData &calc_phase_envelope_data(){return PhaseEnvelope;}; void resize(unsigned int N); diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index ab999dad..1311ae34 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -96,6 +96,14 @@ long double MixtureDerivatives::d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, } } +long double MixtureDerivatives::dln_fugacity_i_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag) +{ + return dln_fugacity_coefficient_dT__constp_n(HEOS, i, xN_flag); +} +long double MixtureDerivatives::dln_fugacity_i_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag) +{ + return dln_fugacity_coefficient_dp__constT_n(HEOS, i, xN_flag) + 1/HEOS.p(); +} long double MixtureDerivatives::fugacity_i(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){ return HEOS.mole_fractions[i]*HEOS.rhomolar()*HEOS.gas_constant()*HEOS.T()*exp( dnalphar_dni__constT_V_nj(HEOS, i, xN_flag)); } @@ -148,6 +156,20 @@ long double MixtureDerivatives::dln_fugacity_i_ddelta__consttau_x(HelmholtzEOSMi { return 1 + HEOS.delta()*HEOS.dalphar_dDelta() + HEOS.delta()*d_ndalphardni_dDelta(HEOS, i, xN_flag); } +long double MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) +{ + // This is a term to which some more might be added depending on i and j + long double val = dln_fugacity_coefficient_dxj__constT_p_xi(HEOS, i, j, xN_flag); + const std::vector &x = HEOS.get_const_mole_fractions(); + std::size_t N = x.size(); + if (i == N-1){ + val += -1/x[N-1]; + } + else if (i == j){ + val += 1/x[j]; + } + return val; +} long double MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) { long double rhor = HEOS.Reducing->rhormolar(HEOS.get_const_mole_fractions()); @@ -363,6 +385,8 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss00 << "Mixture with " << Ncomp << " components"; SECTION(ss00.str(),"") { + x_N_dependency_flag xN_flag = XN_INDEPENDENT; + std::vector names; std::vector z; shared_ptr HEOS, HEOS_plusT_constrho, HEOS_plusrho_constT, HEOS_minusT_constrho, HEOS_minusrho_constT; @@ -414,7 +438,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") rHEOS.update(DmolarT_INPUTS, rho1, T1); - x_N_dependency_flag xN_flag = XN_DEPENDENT; // These ones only require the i index for (std::size_t i = 0; i< z.size();++i) { @@ -423,7 +446,9 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") HEOS_pluszi->specify_phase(iphase_gas); HelmholtzEOSMixtureBackend &rHEOS_pluszi = *(HEOS_pluszi.get()); std::vector zp = z; /// Copy base composition - zp[i] += dz; zp[z.size()-1] -= dz; + zp[i] += dz; + if (xN_flag = XN_DEPENDENT) + zp[z.size()-1] -= dz; rHEOS_pluszi.set_mole_fractions(zp); rHEOS_pluszi.update(DmolarT_INPUTS, rho1, T1); @@ -431,7 +456,9 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") HEOS_minuszi->specify_phase(iphase_gas); HelmholtzEOSMixtureBackend &rHEOS_minuszi = *(HEOS_minuszi.get()); std::vector zm = z; /// Copy base composition - zm[i] -= dz; zm[z.size()-1] += dz; + zm[i] -= dz; + if (xN_flag = XN_DEPENDENT) + zm[z.size()-1] += dz; rHEOS_minuszi.set_mole_fractions(zm); rHEOS_minuszi.update(DmolarT_INPUTS, rho1, T1); @@ -448,6 +475,19 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-8); } + std::ostringstream ss0a; + ss0a << "dln_fugacity_i_dp__constT_n, i=" << i; + SECTION(ss0a.str(),"") + { + double analytic = MixtureDerivatives::dln_fugacity_i_dp__constT_n(rHEOS, i, xN_flag); + double v1 = log(MixtureDerivatives::fugacity_i(rHEOS_plusrho_constT, i, xN_flag)), p1 = rHEOS_plusrho_constT.p(); + double v2 = log(MixtureDerivatives::fugacity_i(rHEOS_minusrho_constT, i, xN_flag)), p2 = rHEOS_minusrho_constT.p(); + double numeric = (v1 - v2)/(p1 - p2); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } std::ostringstream ss0b; ss0b << "dln_fugacity_i_drho__constT_n, i=" << i; SECTION(ss0b.str(), "") @@ -491,6 +531,25 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double err = std::abs((numeric-analytic)/analytic); CHECK(err < 1e-6); } + std::ostringstream ss1a; + ss1a << "dln_fugacity_i_dT__constp_n, i=" << i; + SECTION(ss1a.str(), "") + { + double T1 = 300, dT = 1e-3; + rHEOS.specify_phase(iphase_gas); + + rHEOS.update(PT_INPUTS, 101325, T1); + double analytic = MixtureDerivatives::dln_fugacity_i_dT__constp_n(rHEOS, i, xN_flag); + + rHEOS.update(PT_INPUTS, 101325, T1 + dT); + double v1 = log(MixtureDerivatives::fugacity_i(rHEOS, i, xN_flag)); + rHEOS.update(PT_INPUTS, 101325, T1 - dT); + double v2 = log(MixtureDerivatives::fugacity_i(rHEOS, i, xN_flag)); + + double numeric = (v1 - v2)/(2*dT); + double err = std::abs((numeric-analytic)/analytic); + CHECK(err < 1e-6); + } std::ostringstream ss3; ss3 << "d_ndalphardni_dDelta, i=" << i; diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index 521357c7..94588343 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.h +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -124,9 +124,13 @@ class MixtureDerivatives{ */ static long double dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag); + static long double dln_fugacity_i_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag); + static long double dln_fugacity_i_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag); + static long double dln_fugacity_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); static long double dln_fugacity_i_dtau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag); static long double dln_fugacity_i_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag); static long double dln_fugacity_dxj__constT_rho_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + /** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions * * The derivative term diff --git a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp index 3cfada43..4ecb3e36 100644 --- a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp +++ b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp @@ -15,7 +15,7 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) // Set some imput options SaturationSolvers::mixture_VLE_IO io; io.sstype = SaturationSolvers::imposed_p; - io.Nstep_max = 2; + io.Nstep_max = 20; // Set the pressure to a low pressure HEOS._p = 100000; @@ -44,6 +44,22 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) IO.p = io.p; IO.Nstep_max = 30; + /* + IO.p = 1e5; + IO.rhomolar_liq = 17257.17130; + IO.rhomolar_vap = 56.80022884; + IO.T = 219.5200523; + IO.x[0] = 0.6689704673; + IO.x[1] = 0.3310295327; + */ + IO.rhomolar_liq *= 1.2; + IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::P_IMPOSED; + + NR.call(HEOS, IO.y, IO.x, IO); + + // Switch to density imposed + IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::RHOV_IMPOSED; + bool dont_extrapolate = false; PhaseEnvelopeData &env = HEOS.PhaseEnvelope; @@ -63,7 +79,6 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) } if (iter - iter0 > 0){ IO.rhomolar_vap *= factor;} - long double x = IO.rhomolar_vap; if (dont_extrapolate) { // Reset the step to a reasonably small size @@ -71,29 +86,29 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) } else if (iter - iter0 == 2) { - IO.T = LinearInterp(env.rhomolar_vap, env.T, iter-2, iter-1, x); - IO.rhomolar_liq = LinearInterp(env.rhomolar_vap, env.rhomolar_liq, iter-2, iter-1, x); + IO.T = LinearInterp(env.rhomolar_vap, env.T, iter-2, iter-1, IO.rhomolar_vap); + IO.rhomolar_liq = LinearInterp(env.rhomolar_vap, env.rhomolar_liq, iter-2, iter-1, IO.rhomolar_vap); for (std::size_t i = 0; i < IO.x.size()-1; ++i) // First N-1 elements { - IO.x[i] = LinearInterp(env.rhomolar_vap, env.x[i], iter-2, iter-1, x); + IO.x[i] = LinearInterp(env.rhomolar_vap, env.x[i], iter-2, iter-1, IO.rhomolar_vap); } } else if (iter - iter0 == 3) { - IO.T = QuadInterp(env.rhomolar_vap, env.T, iter-3, iter-2, iter-1, x); - IO.rhomolar_liq = QuadInterp(env.rhomolar_vap, env.rhomolar_liq, iter-3, iter-2, iter-1, x); + IO.T = QuadInterp(env.rhomolar_vap, env.T, iter-3, iter-2, iter-1, IO.rhomolar_vap); + IO.rhomolar_liq = QuadInterp(env.rhomolar_vap, env.rhomolar_liq, iter-3, iter-2, iter-1, IO.rhomolar_vap); for (std::size_t i = 0; i < IO.x.size()-1; ++i) // First N-1 elements { - IO.x[i] = QuadInterp(env.rhomolar_vap, env.x[i], iter-3, iter-2, iter-1, x); + IO.x[i] = QuadInterp(env.rhomolar_vap, env.x[i], iter-3, iter-2, iter-1, IO.rhomolar_vap); } } else if (iter - iter0 > 3) { - IO.T = CubicInterp(env.rhomolar_vap, env.T, iter-4, iter-3, iter-2, iter-1, x); - IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iter-4, iter-3, iter-2, iter-1, x); + IO.T = CubicInterp(env.rhomolar_vap, env.T, iter-4, iter-3, iter-2, iter-1, IO.rhomolar_vap); + IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iter-4, iter-3, iter-2, iter-1, IO.rhomolar_vap); // Check if there is a large deviation from linear interpolation - this suggests a step size that is so large that a minima or maxima of the interpolation function is crossed - long double T_linear = LinearInterp(env.rhomolar_vap, env.T, iter-2, iter-1, x); + long double T_linear = LinearInterp(env.rhomolar_vap, env.T, iter-2, iter-1, IO.rhomolar_vap); if (std::abs((T_linear-IO.T)/IO.T) > 0.1){ // Try again, but with a smaller step IO.rhomolar_vap /= factor; @@ -103,7 +118,7 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) } for (std::size_t i = 0; i < IO.x.size()-1; ++i) // First N-1 elements { - IO.x[i] = CubicInterp(env.rhomolar_vap, env.x[i], iter-4, iter-3, iter-2, iter-1, x); + IO.x[i] = CubicInterp(env.rhomolar_vap, env.x[i], iter-4, iter-3, iter-2, iter-1, IO.rhomolar_vap); if (IO.x[i] < 0 || IO.x[i] > 1){ // Try again, but with a smaller step IO.rhomolar_vap /= factor; @@ -123,13 +138,18 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) // Dewpoint calculation, liquid (x) is incipient phase try{ NR.call(HEOS, IO.y, IO.x, IO); + if (!ValidNumber(IO.rhomolar_liq) || !ValidNumber(IO.p) || !ValidNumber(IO.T)){ + throw ValueError("Invalid number"); + } } catch(std::exception &e){ 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; failure_count++; + continue; } diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index 142ae08f..75291cec 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -722,7 +722,7 @@ void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend &HEOS long double rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(T, p, iphase_liquid); // [mol/m^3] long double rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(T, p, iphase_gas); // [mol/m^3] - HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq*1.3); + HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq*1.5); HEOS.SatV->update_TP_guessrho(T, p, rhomolar_vap); do @@ -796,9 +796,24 @@ void SaturationSolvers::newton_raphson_saturation::resize(unsigned int N) x.resize(N); y.resize(N); - r.resize(N+1); - negative_r.resize(N+1); - J.resize(N+1, std::vector(N+1, 0)); + if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){ + r.resize(N+1); + negative_r.resize(N+1); + J.resize(N+1); + for (std::size_t i = 0; i < N+1; ++i){ + J[i].resize(N+1); + } + err_rel.resize(N+1); + } + else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){ + r.resize(N); + negative_r.resize(N); + J.resize(N, std::vector(N, 0)); + err_rel.resize(N); + } + else{ + throw ValueError(); + } } void SaturationSolvers::newton_raphson_saturation::check_Jacobian() { @@ -893,16 +908,18 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke int iter = 0; if (get_debug_level() > 9){std::cout << " NRsat::call: p" << IO.p << " T" << IO.T << " dl" << IO.rhomolar_liq << " dv" << IO.rhomolar_vap << std::endl;} - + // Reset all the variables and resize pre_call(); - resize(z.size()); this->bubble_point = IO.bubble_point; rhomolar_liq = IO.rhomolar_liq; rhomolar_vap = IO.rhomolar_vap; T = IO.T; p = IO.p; + imposed_variable = IO.imposed_variable; + + resize(z.size()); if (bubble_point){ // Bubblepoint, vapor (y) is the incipient phase @@ -928,9 +945,6 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke // Solve for the step; v is the step with the contents // [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)] std::vector v = linsolve(J, negative_r); - - max_rel_change = max_abs_value(v); - min_abs_change = min_abs_value(v); if (bubble_point){ for (unsigned int i = 0; i < N-1; ++i){ @@ -940,21 +954,26 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke } else{ for (unsigned int i = 0; i < N-1; ++i){ + err_rel[i] = v[i]/y[i]; x[i] += v[i]; } x[N-1] = 1 - std::accumulate(x.begin(), x.end()-1, 0.0); } - T += v[N-1]; - rhomolar_liq += v[N]; - //std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " " << vec_to_string(v, "%0.10Lg") << " " << vec_to_string(x, "%0.10Lg") << std::endl; - + T += v[N-1]; err_rel[N-1] = v[N-1]/T; + if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){ + rhomolar_liq += v[N]; err_rel[N] = v[N]/rhomolar_liq; + } + //std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " v " << vec_to_string(v, "%0.10Lg") << " x " << vec_to_string(x, "%0.10Lg") << " r " << vec_to_string(r, "%0.10Lg") << std::endl; + + min_rel_change = min_abs_value(err_rel); iter++; if (iter == IO.Nstep_max){ throw ValueError(format("newton_raphson_saturation::call reached max number of iterations [%d]",IO.Nstep_max)); } } - while(this->error_rms > 1e-7 && max_rel_change > 1000*LDBL_EPSILON && min_abs_change > 100*DBL_EPSILON && iter < IO.Nstep_max); + while(this->error_rms > 1e-9 && min_rel_change > 10*DBL_EPSILON && iter < IO.Nstep_max); + IO.Nsteps = iter; IO.p = p; IO.x = x; // Mole fractions in liquid @@ -983,54 +1002,93 @@ 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); } - rSatL.update(DmolarT_INPUTS, rhomolar_liq, T); - rSatV.update(DmolarT_INPUTS, rhomolar_vap, T); - - //std::cout << format("\t\t%Lg %Lg %Lg %s\n", T, rhomolar_liq, rhomolar_vap, vec_to_string(x,"%0.12Lg").c_str()); + if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){ + rSatL.update(DmolarT_INPUTS, rhomolar_liq, T); + rSatV.update(DmolarT_INPUTS, rhomolar_vap, T); + } + else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){ + rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar(); + rSatV.update_TP_guessrho(T, p, rhomolar_vap); + } + else{ + throw ValueError(); + } // For diagnostic purposes calculate the pressures (no derivatives are evaluated) long double p_liq = rSatL.p(); long double p_vap = rSatV.p(); p = 0.5*(p_liq + p_vap); - + // Step 2: // ------- // Build the residual vector and the Jacobian matrix 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) - { - // Equate the liquid and vapor fugacities - long double ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag)); - long double ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag)); - r[i] = ln_f_liq - ln_f_vap; + + 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) + { + // Equate the liquid and vapor fugacities + long double ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag)); + long double ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag)); + r[i] = ln_f_liq - ln_f_vap; + + if (bubble_point){ + for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2 + J[i][j] = -MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatV, i, j, xN_flag); + } + } + else{ + for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2 + J[i][j] = MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatL, i, j, xN_flag); + } + } + J[i][N-1] = MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatV, i, xN_flag); + J[i][N] = MixtureDerivatives::dln_fugacity_i_drho__constT_n(rSatL, i, xN_flag); + } + // --------------------------------------------------------------- + // Derivatives of pL(T,rho',x)-p(T,rho'',y) with respect to inputs + // --------------------------------------------------------------- + r[N] = p_liq - p_vap; + for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2 + J[N][j] = MixtureDerivatives::dpdxj__constT_V_xi(rSatL, j, xN_flag); // p'' not a function of x0 + } + // Fixed composition derivatives + J[N][N-1] = rSatL.first_partial_deriv(iP, iT, iDmolar)-rSatV.first_partial_deriv(iP, iT, iDmolar); + J[N][N] = rSatL.first_partial_deriv(iP, iDmolar, iT); + } + else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){ + // Independent variables are N-1 mole fractions of incipient phase, T and rho' - if (bubble_point){ - for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2 - J[i][j] = -MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatV, i, j, xN_flag); + // For the residuals F_i (equality of fugacities) + for (std::size_t i = 0; i < N; ++i) + { + // Equate the liquid and vapor fugacities + long double ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag)); + long double ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag)); + r[i] = ln_f_liq - ln_f_vap; + + if (bubble_point){ + for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2 + J[i][j] = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag); + } } - } - else{ - for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2 - J[i][j] = MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatL, i, j, xN_flag); + else{ + for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2 + J[i][j] = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag); + } } + J[i][N-1] = MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag); } - J[i][N-1] = MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatV, i, xN_flag); - J[i][N] = MixtureDerivatives::dln_fugacity_i_drho__constT_n(rSatL, i, xN_flag); } - // --------------------------------------------------------------- - // Derivatives of pL(T,rho',x)-p(T,rho'',y) with respect to inputs - // --------------------------------------------------------------- - r[N] = p_liq - p_vap; - for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2 - J[N][j] = MixtureDerivatives::dpdxj__constT_V_xi(rSatL, j, xN_flag); // p'' not a function of x0 + else{ + throw ValueError(); } - // Fixed composition derivatives - J[N][N-1] = rSatL.first_partial_deriv(iP, iT, iDmolar)-rSatV.first_partial_deriv(iP, iT, iDmolar); - J[N][N] = rSatL.first_partial_deriv(iP, iDmolar, iT); // Flip all the signs of the entries in the residual vector since we are solving Jv = -r, not Jv=r // Also calculate the rms error of the residual vector at this step diff --git a/src/Backends/Helmholtz/VLERoutines.h b/src/Backends/Helmholtz/VLERoutines.h index 10868d32..1235b8ad 100644 --- a/src/Backends/Helmholtz/VLERoutines.h +++ b/src/Backends/Helmholtz/VLERoutines.h @@ -201,7 +201,7 @@ namespace SaturationSolvers }; struct newton_raphson_saturation_options{ - enum imposed_variable_options {IMPOSED_P, IMPOSED_T}; + enum imposed_variable_options {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, RHOV_IMPOSED}; int Nstep_max; bool bubble_point; std::size_t Nsteps; @@ -219,8 +219,9 @@ namespace SaturationSolvers */ class newton_raphson_saturation { - public: - long double error_rms, rhomolar_liq, rhomolar_vap, T, p, max_rel_change, min_abs_change; + public: + newton_raphson_saturation_options::imposed_variable_options imposed_variable; + long double error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change; unsigned int N; bool logging; bool bubble_point; @@ -228,7 +229,7 @@ namespace SaturationSolvers long double dTsat_dPsat, dPsat_dTsat; STLMatrix J; HelmholtzEOSMixtureBackend *HEOS; - std::vector K, x, y, phi_ij_liq, phi_ij_vap, dlnphi_drho_liq, dlnphi_drho_vap, r, negative_r, dXdS, neg_dFdS; + std::vector K, x, y, phi_ij_liq, phi_ij_vap, dlnphi_drho_liq, dlnphi_drho_vap, r, negative_r, dXdS, neg_dFdS, err_rel; std::vector step_logger; newton_raphson_saturation(){};