diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index f2d4ef2c..2826cbc9 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -251,22 +251,70 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS) } else { - // Set some input options - SaturationSolvers::mixture_VLE_IO options; - options.sstype = SaturationSolvers::imposed_T; - options.Nstep_max = 20; + 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; + + // 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){ + if (std::abs(HEOS._Q-1) > 1e-6){ throw NotImplementedError("Q = 1 only for now"); } + std::size_t &imax = solutions[0]; + SaturationSolvers::newton_raphson_saturation NR; + SaturationSolvers::newton_raphson_saturation_options IO; + IO.bubble_point = false; // because Q = 1 + 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.y = HEOS.get_mole_fractions(); // because Q = 1 + IO.x = IO.y; // Just to give it good size + 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); + 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(); + } + else{ + // Set some input options + SaturationSolvers::mixture_VLE_IO options; + options.sstype = SaturationSolvers::imposed_T; + options.Nstep_max = 20; - // Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted - long double pguess = SaturationSolvers::saturation_preconditioner(HEOS, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions); + // Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted + long double pguess = SaturationSolvers::saturation_preconditioner(HEOS, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions); - // Use Wilson iteration to obtain updated guess for pressure - pguess = SaturationSolvers::saturation_Wilson(HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess); + // Use Wilson iteration to obtain updated guess for pressure + pguess = SaturationSolvers::saturation_Wilson(HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess); - // Actually call the successive substitution solver - SaturationSolvers::successive_substitution(HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options); + // Actually call the successive substitution solver + SaturationSolvers::successive_substitution(HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options); - HEOS._p = options.p; - HEOS._rhomolar = 1/(HEOS._Q/options.rhomolar_vap+(1-HEOS._Q)/options.rhomolar_liq); + HEOS._p = options.p; + HEOS._rhomolar = 1/(HEOS._Q/options.rhomolar_vap+(1-HEOS._Q)/options.rhomolar_liq); + } } } void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS) diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index 66d8902b..a511cd9a 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -805,7 +805,7 @@ void SaturationSolvers::newton_raphson_saturation::resize(unsigned int N) } err_rel.resize(N+1); } - else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){ + else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED){ r.resize(N); negative_r.resize(N); J.resize(N, std::vector(N, 0)); @@ -959,7 +959,12 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke } x[N-1] = 1 - std::accumulate(x.begin(), x.end()-1, 0.0); } - T += v[N-1]; err_rel[N-1] = v[N-1]/T; + if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){ + T += v[N-1]; err_rel[N-1] = v[N-1]/T; + } + else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED){ + p += v[N-1]; err_rel[N-1] = v[N-1]/p; + } if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){ rhomolar_liq += v[N]; err_rel[N] = v[N]/rhomolar_liq; } @@ -1008,9 +1013,9 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays() rSatL.update(DmolarT_INPUTS, rhomolar_liq, T); rSatV.update(DmolarT_INPUTS, rhomolar_vap, T); } - else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){ + else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED){ rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar(); - rSatV.update_TP_guessrho(T, p, rhomolar_vap); + rSatV.update_TP_guessrho(T, p, rhomolar_vap); rhomolar_vap = rSatV.rhomolar(); } else{ throw ValueError("imposed variable not set for NR VLE"); @@ -1084,6 +1089,30 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays() 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); } } + else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED){ + // Independent variables are N-1 mole fractions of incipient phase, p and rho' + + // 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_p_xi(rSatL, i, j, xN_flag); + } + } + J[i][N-1] = 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(); } diff --git a/src/Backends/Helmholtz/VLERoutines.h b/src/Backends/Helmholtz/VLERoutines.h index 1235b8ad..4a23b5a6 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 {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, RHOV_IMPOSED}; + enum imposed_variable_options {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, RHOV_IMPOSED, T_IMPOSED}; int Nstep_max; bool bubble_point; std::size_t Nsteps;