From 20e55df16e83045ed17a9f71aa840693d8882ace Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sun, 14 Sep 2014 22:11:59 +0200 Subject: [PATCH] Added solver for mixtures with vapor quality set and one of P or T. Q no longer needs to be 0 or 1 Signed-off-by: Ian Bell --- src/Backends/Helmholtz/FlashRoutines.cpp | 324 +++++++++++++---------- src/Backends/Helmholtz/FlashRoutines.h | 6 + src/Backends/Helmholtz/VLERoutines.cpp | 165 ++++++++++++ src/Backends/Helmholtz/VLERoutines.h | 69 ++++- 4 files changed, 424 insertions(+), 140 deletions(-) diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index bf24cfbf..7f52c331 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -252,76 +252,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS) else { if(HEOS.PhaseEnvelope.built){ - // Find the intersections in the phase envelope - std::vector< std::pair > intersections = PhaseEnvelopeRoutines::find_intersections(HEOS, iT, HEOS._T); - - PhaseEnvelopeData &env = HEOS.PhaseEnvelope; - - std::size_t quality; - if (std::abs(HEOS._Q) < 100*DBL_EPSILON){ - quality = 0; - } - else if (std::abs(HEOS._Q - 1) < 100*DBL_EPSILON){ - quality = 1; - } - else{ - throw ValueError("Quality is neither 0 or 1"); - } - - // Find the correct solution - std::vector solutions; - for (std::vector< std::pair >::iterator it = intersections.begin(); it != intersections.end(); ++it){ - if (std::abs(env.Q[it->first] - HEOS._Q) < 10*DBL_EPSILON && std::abs(env.Q[it->second] - HEOS._Q) < 10*DBL_EPSILON ){ - solutions.push_back(it->first); - } - } - if (solutions.size() == 1){ - std::size_t &imax = solutions[0]; - SaturationSolvers::newton_raphson_saturation NR; - SaturationSolvers::newton_raphson_saturation_options IO; - IO.T = HEOS._T; - IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::T_IMPOSED; - // p -> 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 == 1){ - IO.bubble_point = false; - IO.y = HEOS.get_mole_fractions(); - IO.x.resize(IO.y.size()); - 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], imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); - } - IO.x[IO.x.size()-1] = 1 - std::accumulate(IO.x.begin(), IO.x.end()-1, 0.0); - NR.call(HEOS, IO.y, IO.x, IO); - } - else{ - IO.bubble_point = true; - IO.x = HEOS.get_mole_fractions(); - IO.y.resize(IO.x.size()); - 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); - } - 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); - } - - - } - else if (solutions.size() == 0){ - throw ValueError("No solution was found in PQ_flash"); - } - else{ - throw ValueError("More than 1 solution was "); - } - // Load the outputs - HEOS._phase = iphase_twophase; - HEOS._p = HEOS.SatV->p(); - HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar()); - HEOS._T = HEOS.SatL->T(); + PT_Q_flash_mixtures(HEOS, iT, HEOS._T); } else{ // Set some input options @@ -341,6 +272,11 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS) HEOS._p = options.p; HEOS._rhomolar = 1/(HEOS._Q/options.rhomolar_vap+(1-HEOS._Q)/options.rhomolar_liq); } + // Load the outputs + HEOS._phase = iphase_twophase; + HEOS._p = HEOS.SatV->p(); + HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar()); + HEOS._T = HEOS.SatL->T(); } } void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS) @@ -425,65 +361,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS) else { if (HEOS.PhaseEnvelope.built){ - // Find the intersections in the phase envelope - std::vector< std::pair > intersections = PhaseEnvelopeRoutines::find_intersections(HEOS, iP, HEOS._p); - - PhaseEnvelopeData &env = HEOS.PhaseEnvelope; - - std::size_t quality; - if (std::abs(HEOS._Q) < 100*DBL_EPSILON){ - quality = 0; - } - else if (std::abs(HEOS._Q - 1) < 100*DBL_EPSILON){ - quality = 1; - } - else{ - throw ValueError("Quality is neither 0 or 1"); - } - - // Find the correct solution - std::vector solutions; - for (std::vector< std::pair >::iterator it = intersections.begin(); it != intersections.end(); ++it){ - if (std::abs(env.Q[it->first] - HEOS._Q) < 10*DBL_EPSILON && std::abs(env.Q[it->second] - HEOS._Q) < 10*DBL_EPSILON ){ - solutions.push_back(it->first); - } - } - if (solutions.size() == 1){ - std::size_t &imax = solutions[0]; - SaturationSolvers::newton_raphson_saturation NR; - SaturationSolvers::newton_raphson_saturation_options IO; - if (quality == 1){ - IO.bubble_point = false; - } - else{ - IO.bubble_point = true; - } - IO.p = HEOS._p; - 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.y = HEOS.get_mole_fractions(); // because Q = 1 - IO.x.resize(IO.y.size()); - IO.T = CubicInterp(env.rhomolar_vap, env.T, 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); - 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], imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); - } - IO.x[IO.x.size()-1] = 1 - std::accumulate(IO.x.begin(), IO.x.end()-1, 0.0); - NR.call(HEOS, IO.y, IO.x, IO); - } - else if (solutions.size() == 0){ - throw ValueError("No solution was found in PQ_flash"); - } - else{ - throw ValueError("More than 1 solution was "); - } - // Load the outputs - HEOS._phase = iphase_twophase; - HEOS._p = HEOS.SatV->p(); - HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar()); - HEOS._T = HEOS.SatL->T(); + PT_Q_flash_mixtures(HEOS, iP, HEOS._p); } else{ // Set some imput options @@ -499,13 +377,189 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS) // Actually call the successive substitution solver SaturationSolvers::successive_substitution(HEOS, HEOS._Q, Tguess, HEOS._p, HEOS.mole_fractions, HEOS.K, io); - - // Load the outputs - HEOS._phase = iphase_twophase; - HEOS._p = HEOS.SatV->p(); - HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar()); - HEOS._T = HEOS.SatL->T(); } + + // Load the outputs + HEOS._phase = iphase_twophase; + HEOS._p = HEOS.SatV->p(); + HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar()); + HEOS._T = HEOS.SatL->T(); + } +} + +void FlashRoutines::PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS, parameters other, long double value) +{ + + // Find the intersections in the phase envelope + std::vector< std::pair > intersections = PhaseEnvelopeRoutines::find_intersections(HEOS, other, value); + + PhaseEnvelopeData &env = HEOS.PhaseEnvelope; + + enum quality_options{SATURATED_LIQUID, SATURATED_VAPOR, TWO_PHASE}; + quality_options quality; + if (std::abs(HEOS._Q) < 100*DBL_EPSILON){ + quality = SATURATED_LIQUID; + } + else if (std::abs(HEOS._Q - 1) < 100*DBL_EPSILON){ + quality = SATURATED_VAPOR; + } + else if (HEOS._Q > 0 && HEOS._Q < 1){ + quality = TWO_PHASE; + } + else{ + throw ValueError("Quality is not within 0 and 1"); + } + + if (quality == SATURATED_LIQUID || quality == SATURATED_VAPOR) + { + // ********************************************************* + // Bubble- or dew-point calculation + // ********************************************************* + // Find the correct solution + std::vector solutions; + for (std::vector< std::pair >::iterator it = intersections.begin(); it != intersections.end(); ++it){ + if (std::abs(env.Q[it->first] - HEOS._Q) < 10*DBL_EPSILON && std::abs(env.Q[it->second] - HEOS._Q) < 10*DBL_EPSILON ){ + solutions.push_back(it->first); + } + } + + if (solutions.size() == 1){ + + std::size_t &imax = solutions[0]; + + SaturationSolvers::newton_raphson_saturation NR; + SaturationSolvers::newton_raphson_saturation_options IO; + + if (other == iP){ + IO.p = HEOS._p; + 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)); + } + else if (other == iT){ + IO.T = HEOS._T; + IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::T_IMPOSED; + // p -> rhomolar_vap + IO.rhomolar_vap = CubicInterp(env.T, env.rhomolar_vap, imax-1, imax, imax+1, imax+2, static_cast(IO.T)); + } + + if (quality == SATURATED_VAPOR){ + IO.bubble_point = false; + IO.y = HEOS.get_mole_fractions(); // Because Q = 1 + IO.x.resize(IO.y.size()); + 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], imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); + } + IO.x[IO.x.size()-1] = 1 - std::accumulate(IO.x.begin(), IO.x.end()-1, 0.0); + NR.call(HEOS, IO.y, IO.x, IO); + } + else{ + IO.bubble_point = true; + IO.x = HEOS.get_mole_fractions(); // Because Q = 0 + IO.y.resize(IO.x.size()); + 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); + } + 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); + } + } + else if (solutions.size() == 0){ + throw ValueError("No solution was found in PQ_flash"); + } + else{ + throw ValueError("More than 1 solution was found in PQ_flash"); + } + } + else{ + // ********************************************************* + // Two-phase calculation for given vapor quality + // ********************************************************* + + // Find the correct solution + std::vector liquid_solutions, vapor_solutions; + for (std::vector< std::pair >::iterator it = intersections.begin(); it != intersections.end(); ++it){ + if (std::abs(env.Q[it->first] - 0) < 10*DBL_EPSILON && std::abs(env.Q[it->second] - 0) < 10*DBL_EPSILON ){ + liquid_solutions.push_back(it->first); + } + if (std::abs(env.Q[it->first] - 1) < 10*DBL_EPSILON && std::abs(env.Q[it->second] - 1) < 10*DBL_EPSILON ){ + vapor_solutions.push_back(it->first); + } + } + + if (liquid_solutions.size() != 1 || vapor_solutions.size() != 1){ + throw ValueError(format("Number liquid solutions [%d] or vapor solutions [%d] != 1", liquid_solutions.size(), vapor_solutions.size() )); + } + std::size_t iliq = liquid_solutions[0], ivap = vapor_solutions[0]; + + SaturationSolvers::newton_raphson_twophase NR; + SaturationSolvers::newton_raphson_twophase_options IO; + IO.beta = HEOS._Q; + + long double rhomolar_vap_sat_vap, T_sat_vap, rhomolar_liq_sat_vap, rhomolar_liq_sat_liq, T_sat_liq, rhomolar_vap_sat_liq, p_sat_liq, p_sat_vap; + + if (other == iP){ + IO.p = HEOS._p; + p_sat_liq = IO.p; p_sat_vap = IO.p; + IO.imposed_variable = SaturationSolvers::newton_raphson_twophase_options::P_IMPOSED; + + // Calculate the interpolated values for beta = 0 and beta = 1 + rhomolar_vap_sat_vap = CubicInterp(env.p, env.rhomolar_vap, ivap-1, ivap, ivap+1, ivap+2, static_cast(IO.p)); + T_sat_vap = CubicInterp(env.rhomolar_vap, env.T, ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap); + rhomolar_liq_sat_vap = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap); + + // Phase inversion for liquid solution (liquid is vapor and vice versa) + rhomolar_liq_sat_liq = CubicInterp(env.p, env.rhomolar_vap, iliq-1, iliq, iliq+1, iliq+2, static_cast(IO.p)); + T_sat_liq = CubicInterp(env.rhomolar_vap, env.T, iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq); + rhomolar_vap_sat_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq); + } + else if (other == iT){ + IO.T = HEOS._T; + T_sat_liq = IO.T; T_sat_vap = IO.T; + IO.imposed_variable = SaturationSolvers::newton_raphson_twophase_options::T_IMPOSED; + + // Calculate the interpolated values for beta = 0 and beta = 1 + rhomolar_vap_sat_vap = CubicInterp(env.T, env.rhomolar_vap, ivap-1, ivap, ivap+1, ivap+2, static_cast(IO.T)); + p_sat_vap = CubicInterp(env.rhomolar_vap, env.p, ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap); + rhomolar_liq_sat_vap = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap); + + // Phase inversion for liquid solution (liquid is vapor and vice versa) + rhomolar_liq_sat_liq = CubicInterp(env.T, env.rhomolar_vap, iliq-1, iliq, iliq+1, iliq+2, static_cast(IO.T)); + p_sat_liq = CubicInterp(env.rhomolar_vap, env.p, iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq); + rhomolar_vap_sat_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq); + } + else{ + throw ValueError(); + } + + // Weight the guesses by the vapor mole fraction + IO.rhomolar_vap = IO.beta*rhomolar_vap_sat_vap + (1-IO.beta)*rhomolar_vap_sat_liq; + IO.rhomolar_liq = IO.beta*rhomolar_liq_sat_vap + (1-IO.beta)*rhomolar_liq_sat_liq; + IO.T = IO.beta*T_sat_vap + (1-IO.beta)*T_sat_liq; + IO.p = IO.beta*p_sat_vap + (1-IO.beta)*p_sat_liq; + + IO.z = HEOS.get_mole_fractions(); + IO.x.resize(IO.z.size()); + IO.y.resize(IO.z.size()); + + for (std::size_t i = 0; i < IO.x.size()-1; ++i) // First N-1 elements + { + long double x_sat_vap = CubicInterp(env.rhomolar_vap, env.x[i], ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap); + long double y_sat_vap = CubicInterp(env.rhomolar_vap, env.y[i], ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap); + + long double x_sat_liq = CubicInterp(env.rhomolar_vap, env.y[i], iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq); + long double y_sat_liq = CubicInterp(env.rhomolar_vap, env.x[i], iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq); + + IO.x[i] = IO.beta*x_sat_vap + (1-IO.beta)*x_sat_liq; + IO.y[i] = IO.beta*y_sat_vap + (1-IO.beta)*y_sat_liq; + } + IO.x[IO.x.size()-1] = 1 - std::accumulate(IO.x.begin(), IO.x.end()-1, 0.0); + IO.y[IO.y.size()-1] = 1 - std::accumulate(IO.y.begin(), IO.y.end()-1, 0.0); + std::vector &XX = IO.x; + std::vector &YY = IO.y; + NR.call(HEOS, IO); } } // D given and one of P,H,S,U diff --git a/src/Backends/Helmholtz/FlashRoutines.h b/src/Backends/Helmholtz/FlashRoutines.h index 86e20342..2fd16fdd 100644 --- a/src/Backends/Helmholtz/FlashRoutines.h +++ b/src/Backends/Helmholtz/FlashRoutines.h @@ -36,6 +36,12 @@ public: /// @param HEOS The HelmholtzEOSMixtureBackend to be used static void QT_flash(HelmholtzEOSMixtureBackend &HEOS); + /// Flash for mixture given temperature or pressure and (molar) quality + /// @param HEOS The HelmholtzEOSMixtureBackend to be used + /// @param other The parameter that is imposed, either iT or iP + /// @param value The value for the imposed parameter + static void PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS, parameters other, long double value); + /// Flash for given pressure and temperature /// @param HEOS The HelmholtzEOSMixtureBackend to be used static void PT_flash(HelmholtzEOSMixtureBackend &HEOS); diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index ca1bcfc0..05bd722b 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -1141,4 +1141,169 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays() dPsat_dTsat = -dQ_dTsat/dQ_dPsat; } +void SaturationSolvers::newton_raphson_twophase::call(HelmholtzEOSMixtureBackend &HEOS, newton_raphson_twophase_options &IO) +{ + 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(); + + rhomolar_liq = IO.rhomolar_liq; + rhomolar_vap = IO.rhomolar_vap; + T = IO.T; + p = IO.p; + imposed_variable = IO.imposed_variable; + x = IO.x; + y = IO.y; + z = IO.z; + beta = IO.beta; + + this->N = z.size(); + x.resize(N); + y.resize(N); + r.resize(2*N-1); + negative_r.resize(2*N-1); + J.resize(2*N-1, std::vector(2*N-1, 0)); + err_rel.resize(2*N-1); + + // Hold a pointer to the backend + this->HEOS = &HEOS; + + do + { + // Build the Jacobian and residual vectors + build_arrays(); + + // Solve for the step; v is the step with the contents + // [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)] + + // Uncomment to see Jacobian and residual at every step + // std::cout << vec_to_string(J, "%0.12Lg") << std::endl; + // std::cout << vec_to_string(negative_r, "%0.12Lg") << std::endl; + + std::vector v = linsolve(J, negative_r); + for (unsigned int i = 0; i < N-1; ++i){ + err_rel[i] = v[i]/x[i]; + x[i] += v[i]; + err_rel[i+(N-1)] = v[i+(N-1)]/y[i]; + y[i] += v[i+(N-1)]; + } + x[N-1] = 1 - std::accumulate(x.begin(), x.end()-1, 0.0); + y[N-1] = 1 - std::accumulate(y.begin(), y.end()-1, 0.0); + + if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED){ + T += v[2*N-2]; err_rel[2*N-2] = v[2*N-2]/T; + } + else if (imposed_variable == newton_raphson_twophase_options::T_IMPOSED){ + p += v[2*N-2]; err_rel[2*N-2] = v[2*N-2]/p; + } + else{ + throw ValueError("invalid imposed_variable"); + } + //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-9 && min_rel_change > 1000*DBL_EPSILON && iter < IO.Nstep_max); + + IO.Nsteps = iter; + IO.p = p; + IO.x = x; // Mole fractions in liquid + IO.y = y; // Mole fractions in vapor + IO.T = T; + IO.rhomolar_liq = rhomolar_liq; + IO.rhomolar_vap = rhomolar_vap; + IO.hmolar_liq = HEOS.SatL.get()->hmolar(); + IO.hmolar_vap = HEOS.SatV.get()->hmolar(); + IO.smolar_liq = HEOS.SatL.get()->smolar(); + IO.smolar_vap = HEOS.SatV.get()->smolar(); +} + +void SaturationSolvers::newton_raphson_twophase::build_arrays() +{ + // References to the classes for concision + HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get()); + + // Step 0: + // ------- + // Set mole fractions + rSatL.set_mole_fractions(x); + rSatV.set_mole_fractions(y); + + //std::vector &x = rSatL.get_mole_fractions(); + //std::vector &y = rSatV.get_mole_fractions(); + + rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar(); + rSatV.update_TP_guessrho(T, p, rhomolar_vap); rhomolar_vap = rSatV.rhomolar(); + + // 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; + + // Form of residuals do not depend on which variable is imposed + 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; // N of these + + if (i != N-1){ + // Equate the specified vapor mole fraction and that given defined by the ith component + r[i+N] = (z[i]-x[i])/(y[i]-x[i]) - beta; // N-1 of these + } + } + + // First part of derivatives with respect to ln f_i + for (std::size_t i = 0; i < N; ++i) + { + for (std::size_t j = 0; j < N-1; ++j){ + J[i][j] = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag); + J[i][j+N-1] = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag); + } + + // Last derivative with respect to either T or p depending on what is imposed + if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED){ + J[i][2*N-2] = MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag); + } + else if (imposed_variable == newton_raphson_twophase_options::T_IMPOSED){ + J[i][2*N-2] = MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatV, i, xN_flag); + } + else{ + throw ValueError(); + } + } + // Derivatives with respect to the vapor mole fractions residual + for (std::size_t i = 0; i < N-1; ++i) + { + std::size_t k = i + N; // N ln f_i residuals + J[k][i] = (z[i]-y[i])/pow(y[i]-x[i], 2); + J[k][i+(N-1)] = -(z[i]-x[i])/pow(y[i]-x[i], 2); + } + + // 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 + error_rms = 0; + for (unsigned int i = 0; i < J.size(); ++i) + { + negative_r[i] = -r[i]; + error_rms += r[i]*r[i]; // Sum the squares + } + error_rms = sqrt(error_rms); // Square-root (The R in RMS) +} + } /* namespace CoolProp*/ diff --git a/src/Backends/Helmholtz/VLERoutines.h b/src/Backends/Helmholtz/VLERoutines.h index 4a23b5a6..c2343e95 100644 --- a/src/Backends/Helmholtz/VLERoutines.h +++ b/src/Backends/Helmholtz/VLERoutines.h @@ -199,6 +199,66 @@ namespace SaturationSolvers { long double T,p; }; + + struct newton_raphson_twophase_options{ + enum imposed_variable_options {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, T_IMPOSED}; + int Nstep_max; + std::size_t Nsteps; + long double beta, omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap; + imposed_variable_options imposed_variable; + std::vector x, y, z; + newton_raphson_twophase_options(){ Nstep_max = 30; Nsteps = 0; beta = -1; omega =1;} // Defaults + }; + + /** \brief A class to do newton raphson solver for VLE for p,q or T,q + * + * A class is used rather than a function so that it is easier to store iteration histories, additional output values, etc. + * + * This class only handles bubble and dew lines since the independent variables are N-1 of the mole fractions in the incipient phase along with one of T, p, or rho + */ + class newton_raphson_twophase + { + public: + newton_raphson_twophase_options::imposed_variable_options imposed_variable; + long double error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change, beta; + unsigned int N; + bool logging; + int Nsteps; + STLMatrix J; + HelmholtzEOSMixtureBackend *HEOS; + std::vector K, x, y, z, r, negative_r, err_rel; + std::vector step_logger; + + newton_raphson_twophase(){}; + + void resize(unsigned int N); + + // Reset the state of all the internal variables + void pre_call() + { + K.clear(); x.clear(); y.clear(); step_logger.clear(); error_rms = 1e99; Nsteps = 0; + rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE; + }; + + /** 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 + * before you start with Newton Raphson. + * + * @param HEOS HelmholtzEOSMixtureBackend instance + * @param IO The input/output data structure + */ + void call(HelmholtzEOSMixtureBackend &HEOS, newton_raphson_twophase_options &IO); + + /*! Build the arrays for the Newton-Raphson solve + + This method builds the Jacobian matrix, the sensitivity matrix, etc. + * + */ + void build_arrays(); + }; + struct newton_raphson_saturation_options{ enum imposed_variable_options {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, RHOV_IMPOSED, T_IMPOSED}; @@ -208,7 +268,7 @@ namespace SaturationSolvers long double omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap; imposed_variable_options imposed_variable; std::vector x, y; - newton_raphson_saturation_options(){ Nsteps = 0;} // Defaults + newton_raphson_saturation_options(){ Nstep_max = 30; Nsteps = 0;} // Defaults }; /** \brief A class to do newton raphson solver for VLE given guess values for vapor-liquid equilibria. This class will then be included in the Mixture class @@ -226,10 +286,10 @@ namespace SaturationSolvers bool logging; bool bubble_point; int Nsteps; - 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, err_rel; + long double dTsat_dPsat, dPsat_dTsat; + std::vector K, x, y, r, negative_r, err_rel; std::vector step_logger; newton_raphson_saturation(){}; @@ -239,8 +299,7 @@ namespace SaturationSolvers // Reset the state of all the internal variables void pre_call() { - K.clear(); x.clear(); y.clear(); phi_ij_liq.clear(); - phi_ij_vap.clear(); dlnphi_drho_liq.clear(), dlnphi_drho_vap.clear(), + K.clear(); x.clear(); y.clear(); step_logger.clear(); error_rms = 1e99; Nsteps = 0; rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE; };