8 #if defined(ENABLE_CATCH)
9 # include <catch2/catch_all.hpp>
22 if (!twophase && HEOS.
_T > closest_state.
T) {
70 HEOS.
_Q = (o.
z[0] - o.
x[0]) / (o.
y[0] - o.
x[0]);
97 bool saturation_called =
false;
105 throw ValueError(
"twophase not implemented yet");
162 double omega, R, kappa, a, b, A, B, C, Tc, pc, V = 1 / rhomolar;
168 kappa = 0.37464 + 1.54226 * omega - 0.26992 * omega * omega;
169 a = 0.457235 * R * R * Tc * Tc / pc;
170 b = 0.077796 * R * Tc / pc;
171 double den = V * V + 2 * b * V - b * b;
174 A = R * Tc / (V - b) - a * kappa * kappa / (den);
175 B = +2 * a * kappa * (1 + kappa) / (den);
176 C = -a * (1 + 2 * kappa + kappa * kappa) / (den)-p;
180 double sqrt_Tr1 = (-B + sqrt(B * B - 4 * A * C)) / (2 * A);
182 return sqrt_Tr1 * sqrt_Tr1 * Tc;
193 bool saturation_called =
false;
200 if (saturation_called) {
217 Halley(resid, T0, 1e-10, 100);
260 double Tmin = HEOS.
Tmin() + 0.1;
263 const double eps = 1e-12;
291 if (std::abs(HEOS.
Q() - 1) > 1e-10) {
292 throw ValueError(
format(
"non-unity quality not currently allowed for HQ_flash"));
313 }
else if (std::abs(HEOS.
Q()) < 1e-10) {
324 }
else if (std::abs(HEOS.
Q() - 1) < 1e-10) {
336 throw ValueError(
format(
"non-zero or 1 quality not currently allowed for QS_flash"));
352 Tmin_sat = std::max(Tmin_satL, Tmin_satV) - 1e-13;
358 if ((
get_config_bool(CRITICAL_WITHIN_1UK) && std::abs(T - Tmax_sat) < 1e-6) || std::abs(T - Tmax_sat) < 1e-12) {
362 HEOS.
_p = 0.5 * HEOS.
SatV->p() + 0.5 * HEOS.
SatL->p();
364 throw ValueError(
format(
"Temperature to QT_flash [%0.8Lg K] must be in range [%0.8Lg K, %0.8Lg K]", T, Tmin_sat - 0.1, Tmax_sat));
366 double rhoL = _HUGE, rhoV = _HUGE;
371 HEOS.
_p = 0.5 * HEOS.
SatV->p() + 0.5 * HEOS.
SatL->p();
373 }
else if (!(HEOS.
components[0].EOS().pseudo_pure)) {
380 HEOS.
_p = 0.5 * HEOS.
SatV->p() + 0.5 * HEOS.
SatL->p();
384 CoolPropDbl rhoLanc = _HUGE, rhoVanc = _HUGE, rhoLsat = _HUGE, rhoVsat = _HUGE;
387 rhoLanc = HEOS.
components[0].ancillaries.rhoL.evaluate(HEOS.
_T);
388 HEOS.
SatL->update_TP_guessrho(HEOS.
_T, HEOS.
_p, rhoLanc);
392 rhoVanc = HEOS.
components[0].ancillaries.rhoV.evaluate(HEOS.
_T);
393 HEOS.
SatV->update_TP_guessrho(HEOS.
_T, HEOS.
_p, rhoVanc);
396 throw CoolProp::ValueError(
format(
"For pseudo-pure fluid, quality must be equal to 0 or 1. Two-phase quality is not defined"));
447 NR.
call(HEOS, IO.
x, IO.
y, IO);
450 NR.
call(HEOS, IO.
y, IO.
x, IO);
464 void get_Henrys_coeffs_FP(
const std::string& CAS,
double& A,
double& B,
double& C,
double& Tmin,
double& Tmax) {
466 if (CAS ==
"7440-59-7")
473 }
else if (CAS ==
"7440-01-9")
480 }
else if (CAS ==
"7440-37-1")
487 }
else if (CAS ==
"7439-90-9")
494 }
else if (CAS ==
"7440-63-3")
501 }
else if (CAS ==
"1333-74-0")
508 }
else if (CAS ==
"7727-37-9")
515 }
else if (CAS ==
"7782-44-7")
522 }
else if (CAS ==
"630-08-0")
529 }
else if (CAS ==
"124-38-9")
536 }
else if (CAS ==
"7783-06-4")
543 }
else if (CAS ==
"74-82-8")
550 }
else if (CAS ==
"74-84-0")
557 }
else if (CAS ==
"2551-62-4")
565 throw ValueError(
"Bad component in Henry's law constants");
579 HEOS.
SatL->update_TP_guessrho(HEOS.
_TLanc, HEOS.
_p, rhoL);
580 HEOS.
SatV->update_TP_guessrho(HEOS.
_TVanc, HEOS.
_p, rhoV);
584 HEOS.
_p = HEOS.
_Q * HEOS.
SatV->p() + (1 - HEOS.
_Q) * HEOS.
SatL->p();
585 HEOS.
_T = HEOS.
_Q * HEOS.
SatV->T() + (1 - HEOS.
_Q) * HEOS.
SatL->T();
594 pmin_sat = std::max(pmin_satL, pmin_satV);
609 throw ValueError(
format(
"Pressure to PQ_flash [%6g Pa] must be in range [%8Lg Pa, %8Lg Pa]", HEOS.
_p, pmin_sat, pmax_sat));
623 double increment = 0.4;
626 for (
double omega = 1.0; omega > 0; omega -= increment) {
628 options.
omega = omega;
636 if (omega < 1.1 * increment) {
649 HEOS.
_p = HEOS.
_Q * HEOS.
SatV->p() + (1 - HEOS.
_Q) * HEOS.
SatL->p();
669 std::vector<CoolPropDbl> K = HEOS.
K;
671 if (
get_config_bool(HENRYS_LAW_TO_GENERATE_VLE_GUESSES) && std::abs(HEOS.
_Q - 1) < 1e-10) {
672 const std::vector<CoolPropFluid>& components = HEOS.
get_components();
673 std::size_t iWater = 0;
674 double p1star =
PropsSI(
"P",
"T", Tguess,
"Q", 1,
"Water");
676 std::vector<CoolPropDbl> x(y.size());
677 for (std::size_t i = 0; i < components.size(); ++i) {
678 if (components[i].CAS ==
"7732-18-5") {
682 double A, B, C, Tmin, Tmax;
684 double T_R = Tguess / 647.096, tau = 1 - T_R;
685 double k_H = p1star * exp(A / T_R + B * pow(tau, 0.355) / T_R + C * pow(T_R, -0.41) * exp(tau));
686 x[i] = y[i] * HEOS.
_p / k_H;
693 for (std::size_t i = 0; i < y.size(); ++i) {
699 K[iWater] = y[iWater] / x[iWater];
724 NR.
call(HEOS, IO.
x, IO.
y, IO);
727 NR.
call(HEOS, IO.
y, IO.
x, IO);
744 IO.
x = std::vector<CoolPropDbl>(guess.
x.begin(), guess.
x.end());
745 IO.
y = std::vector<CoolPropDbl>(guess.
y.begin(), guess.
y.end());
751 if (std::abs(HEOS.
Q()) < 1e-10) {
753 NR.
call(HEOS, IO.
x, IO.
y, IO);
754 }
else if (std::abs(HEOS.
Q() - 1) < 1e-10) {
756 NR.
call(HEOS, IO.
y, IO.
x, IO);
771 IO.
x = std::vector<CoolPropDbl>(guess.
x.begin(), guess.
x.end());
772 IO.
y = std::vector<CoolPropDbl>(guess.
y.begin(), guess.
y.end());
783 if (std::abs(HEOS.
Q()) < 1e-10) {
785 NR.
call(HEOS, IO.
x, IO.
y, IO);
786 }
else if (std::abs(HEOS.
Q() - 1) < 1e-10) {
788 NR.
call(HEOS, IO.
y, IO.
x, IO);
827 std::vector<std::pair<std::size_t, std::size_t>> intersections =
838 quality_options quality;
840 quality = SATURATED_LIQUID;
842 quality = SATURATED_VAPOR;
843 }
else if (HEOS.
_Q > 0 && HEOS.
_Q < 1) {
846 throw ValueError(
"Quality is not within 0 and 1");
849 if (quality == SATURATED_LIQUID || quality == SATURATED_VAPOR) {
854 std::vector<std::size_t> solutions;
855 for (std::vector<std::pair<std::size_t, std::size_t>>::const_iterator it = intersections.begin(); it != intersections.end(); ++it) {
857 solutions.push_back(it->first);
861 if (solutions.size() == 1) {
863 std::size_t& imax = solutions[0];
866 if (imax + 2 >= env.T.size()) {
868 }
else if (imax == 0) {
885 }
else if (other ==
iT) {
896 if (quality == SATURATED_VAPOR) {
899 IO.
x.resize(IO.
y.size());
900 for (std::size_t i = 0; i < IO.
x.size() - 1; ++i)
904 IO.
x[IO.
x.size() - 1] = 1 - std::accumulate(IO.
x.begin(), IO.
x.end() - 1, 0.0);
905 NR.
call(HEOS, IO.
y, IO.
x, IO);
909 IO.
y.resize(IO.
x.size());
912 for (std::size_t i = 0; i < IO.
y.size() - 1; ++i)
918 IO.
y[IO.
y.size() - 1] = 1 - std::accumulate(IO.
y.begin(), IO.
y.end() - 1, 0.0);
919 NR.
call(HEOS, IO.
x, IO.
y, IO);
921 }
else if (solutions.size() == 0) {
922 throw ValueError(
"No solution was found in PQ_flash");
924 throw ValueError(
"More than 1 solution was found in PQ_flash");
932 std::vector<std::size_t> liquid_solutions, vapor_solutions;
933 for (std::vector<std::pair<std::size_t, std::size_t>>::const_iterator it = intersections.begin(); it != intersections.end(); ++it) {
934 if (std::abs(env.Q[it->first] - 0) < 10 *
DBL_EPSILON && std::abs(env.Q[it->second] - 0) < 10 *
DBL_EPSILON) {
935 liquid_solutions.push_back(it->first);
937 if (std::abs(env.Q[it->first] - 1) < 10 *
DBL_EPSILON && std::abs(env.Q[it->second] - 1) < 10 *
DBL_EPSILON) {
938 vapor_solutions.push_back(it->first);
942 if (liquid_solutions.size() != 1 || vapor_solutions.size() != 1) {
943 throw ValueError(
format(
"Number liquid solutions [%d] or vapor solutions [%d] != 1", liquid_solutions.size(), vapor_solutions.size()));
945 std::size_t iliq = liquid_solutions[0], ivap = vapor_solutions[0];
951 CoolPropDbl 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,
961 rhomolar_vap_sat_vap =
CubicInterp(env.p, env.rhomolar_vap, ivap - 1, ivap, ivap + 1, ivap + 2,
static_cast<CoolPropDbl>(IO.
p));
962 T_sat_vap =
CubicInterp(env.rhomolar_vap, env.T, ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
963 rhomolar_liq_sat_vap =
CubicInterp(env.rhomolar_vap, env.rhomolar_liq, ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
966 rhomolar_liq_sat_liq =
CubicInterp(env.p, env.rhomolar_vap, iliq - 1, iliq, iliq + 1, iliq + 2,
static_cast<CoolPropDbl>(IO.
p));
967 T_sat_liq =
CubicInterp(env.rhomolar_vap, env.T, iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
968 rhomolar_vap_sat_liq =
CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
969 }
else if (other ==
iT) {
976 rhomolar_vap_sat_vap =
CubicInterp(env.T, env.rhomolar_vap, ivap - 1, ivap, ivap + 1, ivap + 2,
static_cast<CoolPropDbl>(IO.
T));
977 p_sat_vap =
CubicInterp(env.rhomolar_vap, env.p, ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
978 rhomolar_liq_sat_vap =
CubicInterp(env.rhomolar_vap, env.rhomolar_liq, ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
981 rhomolar_liq_sat_liq =
CubicInterp(env.T, env.rhomolar_vap, iliq - 1, iliq, iliq + 1, iliq + 2,
static_cast<CoolPropDbl>(IO.
T));
982 p_sat_liq =
CubicInterp(env.rhomolar_vap, env.p, iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
983 rhomolar_vap_sat_liq =
CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
991 IO.
T = IO.
beta * T_sat_vap + (1 - IO.
beta) * T_sat_liq;
992 IO.
p = IO.
beta * p_sat_vap + (1 - IO.
beta) * p_sat_liq;
995 IO.
x.resize(IO.
z.size());
996 IO.
y.resize(IO.
z.size());
998 for (std::size_t i = 0; i < IO.
x.size() - 1; ++i)
1000 CoolPropDbl x_sat_vap =
CubicInterp(env.rhomolar_vap, env.x[i], ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
1001 CoolPropDbl y_sat_vap =
CubicInterp(env.rhomolar_vap, env.y[i], ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
1003 CoolPropDbl x_sat_liq =
CubicInterp(env.rhomolar_vap, env.y[i], iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
1004 CoolPropDbl y_sat_liq =
CubicInterp(env.rhomolar_vap, env.x[i], iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
1006 IO.
x[i] = IO.
beta * x_sat_vap + (1 - IO.
beta) * x_sat_liq;
1007 IO.
y[i] = IO.
beta * y_sat_vap + (1 - IO.
beta) * y_sat_liq;
1009 IO.
x[IO.
x.size() - 1] = 1 - std::accumulate(IO.
x.begin(), IO.
x.end() - 1, 0.0);
1010 IO.
y[IO.
y.size() - 1] = 1 - std::accumulate(IO.
y.begin(), IO.
y.end() - 1, 0.0);
1025 : HEOS(HEOS), rhomolar_spec(rhomolar_spec), other(other), value(value) {
1028 double call(
double T) {
1032 Qd = (1 / rhomolar_spec - 1 / SatL.
rhomolar()) / (1 / SatV.rhomolar() - 1 / SatL.
rhomolar());
1038 } resid(HEOS, rhomolar_spec, other, value);
1046 Tmin_sat = std::max(Tmin_satL, Tmin_satV) - 1e-13;
1064 : HEOS(HEOS), rhomolar(rhomolar), value(value), other(other), Tmin(Tmin), Tmax(Tmax) {
1068 double call(
double T) {
1073 return (eos - value) / value;
1079 double deriv(
double T) {
1085 double second_deriv(
double T) {
1091 bool input_not_in_range(
double T) {
1092 return (T < Tmin || T > Tmax);
1099 shared_ptr<HelmholtzEOSMixtureBackend> Sat;
1133 y_solid = (yV - yL) / (1 / rhoVtriple - 1 / rhoLtriple) * (1 / HEOS.
_rhomolar - 1 / rhoLtriple) + yL;
1135 if (value < y_solid) {
1144 for (
int i_try = 0; i_try < 7; i_try++)
1166 optionsD.
omega /= 2;
1168 if (i_try >= 6){
throw;}
1173 if (value > Sat->keyed_output(other)) {
1174 solver_resid resid(&HEOS, HEOS.
_rhomolar, value, other, Sat->keyed_output(
iT), HEOS.
Tmax() * 1.5);
1176 HEOS.
_T =
Halley(resid, 0.5 * (Sat->keyed_output(
iT) + HEOS.
Tmax() * 1.5), 1e-10, 100);
1202 HEOS.
_Q = (1 / HEOS.
_rhomolar - 1 / HEOS.
SatL->rhomolar()) / (1 / HEOS.
SatV->rhomolar() - 1 / HEOS.
SatL->rhomolar());
1203 HEOS.
_T = HEOS.
SatL->T();
1239 solver_resid resid(&HEOS, HEOS.
_rhomolar, value, other, TVtriple, HEOS.
Tmax() * 1.5);
1280 solver_resid resid(&HEOS, HEOS.
_rhomolar, value, other, TLtriple, HEOS.
Tmax() * 1.5);
1302 double A[2][2], B[2][2];
1318 throw ValueError(
"other is invalid in HSU_P_flash_singlephase_Newton");
1323 bool failed =
false;
1324 CoolPropDbl omega = 1.0, f2, df2_dtau, df2_ddelta;
1326 while (worst_error > 1e-6 && failed ==
false) {
1342 CoolPropDbl f1 = delta / tau * (1 + delta * dar_ddelta) - p / (rhoc * R * Tc);
1343 CoolPropDbl df1_dtau = (1 + delta * dar_ddelta) * (-delta / tau / tau) + delta / tau * (delta * d2ar_ddelta_dtau);
1344 CoolPropDbl df1_ddelta = (1.0 / tau) * (1 + 2.0 * delta * dar_ddelta + delta * delta * d2ar_ddelta2);
1347 f2 = (1 + delta * dar_ddelta) + tau * (da0_dtau + dar_dtau) - tau * y / (R * Tc);
1348 df2_dtau = delta * d2ar_ddelta_dtau + da0_dtau + dar_dtau + tau * (d2a0_dtau2 + d2ar_dtau2) - y / (R * Tc);
1349 df2_ddelta = (dar_ddelta + delta * d2ar_ddelta2) + tau * (d2a0_ddelta_dtau + d2ar_ddelta_dtau);
1353 f2 = tau * (da0_dtau + dar_dtau) - ar - a0 - y / R;
1354 df2_dtau = tau * (d2a0_dtau2 + d2ar_dtau2) + (da0_dtau + dar_dtau) - dar_dtau - da0_dtau;
1355 df2_ddelta = tau * (d2a0_ddelta_dtau + d2ar_ddelta_dtau) - dar_ddelta - da0_ddelta;
1359 throw ValueError(
"other variable in HSU_P_flash_singlephase_Newton is invalid");
1364 A[0][1] = df1_ddelta;
1366 A[1][1] = df2_ddelta;
1371 tau -= omega * (B[0][0] * f1 + B[0][1] * f2);
1372 delta -= omega * (B[1][0] * f1 + B[1][1] * f2);
1374 if (std::abs(f1) > std::abs(f2))
1375 worst_error = std::abs(f1);
1377 worst_error = std::abs(f2);
1380 throw SolutionError(
format(
"Invalid values for inputs p=%g y=%g for fluid %s in HSU_P_flash_singlephase", p, y, _HEOS.
name().c_str()));
1385 throw SolutionError(
format(
"HSU_P_flash_singlephase did not converge with inputs p=%g h=%g for fluid %s", p, y, _HEOS.
name().c_str()));
1394 throw ValueError(
"value for p in HSU_P_flash_singlephase_Brent is invalid");
1397 throw ValueError(
"value for other in HSU_P_flash_singlephase_Brent is invalid");
1406 CoolPropDbl eos0, eos1, rhomolar, rhomolar0, rhomolar1;
1432 double call(
double T) {
1434 if (iter < 2 || std::abs(rhomolar1 / rhomolar0 - 1) > 0.05) {
1454 rhomolar0 = rhomolar;
1455 }
else if (iter == 1) {
1457 rhomolar1 = rhomolar;
1461 rhomolar0 = rhomolar1;
1462 rhomolar1 = rhomolar;
1468 double deriv(
double T) {
1471 double second_deriv(
double T) {
1474 bool input_not_in_range(
double x) {
1475 return (x < Tmin || x > Tmax);
1478 solver_resid resid(&HEOS, HEOS.
_p, value, other, Tmin, Tmax);
1482 Halley(resid, Tmin, 1e-12, 100);
1484 throw ValueError(
"Halley's method was unable to find a solution in HSU_P_flash_singlephase_Brent");
1504 if (eos1 > eos0 && value > eos1) {
1506 format(
"HSU_P_flash_singlephase_Brent could not find a solution because %s [%Lg %s] is above the maximum value of %0.12Lg %s",
1507 name.c_str(), value, units.c_str(), eos1, units.c_str()));
1509 if (eos1 > eos0 && value < eos0) {
1511 format(
"HSU_P_flash_singlephase_Brent could not find a solution because %s [%Lg %s] is below the minimum value of %0.12Lg %s",
1512 name.c_str(), value, units.c_str(), eos0, units.c_str()));
1521 bool saturation_called =
false;
1549 Tmax = 1.5 * HEOS.
Tmax();
1553 if (saturation_called) {
1554 Tmin = HEOS.
SatV->T();
1562 if (saturation_called) {
1563 Tmax = HEOS.
SatL->T();
1572 Tmin = HEOS.
Tmin() - 1e-3;
1579 Tmax = 1.5 * HEOS.
Tmax();
1584 Tmin = HEOS.
Tmin() - 1e-3;
1594 }
catch (std::exception& e) {
1595 throw ValueError(
format(
"unable to solve 1phase PY flash with Tmin=%Lg, Tmax=%Lg due to error: %s", Tmin, Tmax, e.what()));
1605 std::size_t iclosest;
1614 throw ValueError(
"two-phase solution for Y");
1618 throw ValueError(
"phase envelope must be built to carry out HSU_P_flash for mixture");
1632 : HEOS(HEOS), T(T), value(value), other(other) {}
1633 double call(
double rhomolar) {
1638 double deriv(
double rhomolar) {
1641 double second_deriv(
double rhomolar) {
1645 solver_resid resid(&HEOS, T, value, other);
1677 Brent(resid, rhoc, rhomin, LDBL_EPSILON, 1e-9, 100);
1678 }
else if (y < yc) {
1696 if (step_count > 30) {
1697 throw ValueError(
format(
"Even by increasing rhoc, not able to bound input; input %Lg is not in range %Lg,%Lg", y, yc, ymin));
1701 Brent(resid, rhomin, rhoc, LDBL_EPSILON, 1e-9, 100);
1703 throw ValueError(
format(
"input %Lg is not in range %Lg,%Lg,%Lg", y, yc, ymin));
1741 CoolPropDbl rhomolar_guess = (rhomelt - rhoL) / (ymelt - yL) * (y - yL) + rhoL;
1744 Halley(resid, rhomolar_guess, 1e-8, 100);
1746 Secant(resid, rhomolar_guess, 0.0001 * rhomolar_guess, 1e-12, 100);
1755 Halley(resid, 0.5 * (rhomin + rhoV), 1e-8, 100);
1758 Brent(resid, rhomin, rhoV, LDBL_EPSILON, 1e-12, 100);
1764 throw ValueError(
format(
"phase to solver_for_rho_given_T_oneof_HSU is invalid"));
1803 Q = (1 / HEOS.
rhomolar() - 1 / HEOS1.
SatL->rhomolar()) / (1 / HEOS1.
SatV->rhomolar() - 1 / HEOS1.
SatL->rhomolar());
1806 Q = (HEOS.
smolar() - HEOS1.
SatL->smolar()) / (HEOS1.
SatV->smolar() - HEOS1.
SatL->smolar());
1809 Q = (HEOS.
hmolar() - HEOS1.
SatL->hmolar()) / (HEOS1.
SatV->hmolar() - HEOS1.
SatL->hmolar());
1812 Q = (HEOS.
umolar() - HEOS1.
SatL->umolar()) / (HEOS1.
SatV->umolar() - HEOS1.
SatL->umolar());
1824 HEOS.
_p = HEOS.
_Q * HEOS1.
SatV->p() + (1 - HEOS.
_Q) * HEOS1.
SatL->p();
1830 throw ValueError(
format(
"Temperature specified is not the imposed phase region."));
1895 : HEOS(HEOS), hmolar(hmolar_spec), smolar(smolar_spec), Qs(_HUGE){};
1896 double call(
double T) {
1900 Qs = (smolar - SatL.
smolar()) / (SatV.smolar() - SatL.
smolar());
1906 } resid(HEOS, hmolar_spec, smolar_spec);
1914 Tmin_sat = std::max(Tmin_satL, Tmin_satV) - 1e-13;
1923 double resid = 9e30, resid_old = 9e30;
1929 r(0) = HEOS.
hmolar() - hmolar_spec;
1930 r(1) = HEOS.
smolar() - smolar_spec;
1936 Eigen::Vector2d v = J.colPivHouseholderQr().solve(-r);
1937 bool good_solution =
false;
1938 double tau0 = HEOS.
tau(), delta0 = HEOS.
delta();
1941 for (
double frac = 1.0; frac > 0.001; frac /= 2) {
1944 double tau_new = tau0 + options.
omega * frac * v(0);
1945 double delta_new = delta0 + options.
omega * frac * v(1);
1946 double T_new = reducing.
T / tau_new;
1947 double rhomolar_new = delta_new * reducing.
rhomolar;
1951 if (resid > resid_old) {
1952 throw ValueError(
format(
"residual not decreasing; frac: %g, resid: %g, resid_old: %g", frac, resid, resid_old));
1954 good_solution =
true;
1961 if (!good_solution) {
1966 throw ValueError(
format(
"HS_flash_singlephase took too many iterations; residual is %g; prior was %g", resid, resid_old));
1968 }
while (std::abs(resid) > 1e-9);
1972 double logp = ((double)rand() / (double)RAND_MAX) * (log(HEOS.
pmax()) - log(HEOS.
p_triple())) + log(HEOS.
p_triple());
1973 T = ((double)rand() / (double)RAND_MAX) * (HEOS.
Tmax() - HEOS.
Ttriple()) + HEOS.
Ttriple();
1986 : HEOS(HEOS), hmolar(hmolar_spec), smolar(smolar_spec){};
1987 double call(
double T) {
1989 double r = HEOS.
hmolar() - hmolar;
1992 } resid(HEOS, hmolar, smolar);
1995 bool good_Tmin =
false;
2000 rmin = resid.call(Tmin);
2005 if (Tmin > HEOS.
Tmax()) {
2008 }
while (!good_Tmin);
2011 bool good_Tmax =
false;
2012 double Tmax = HEOS.
Tmax() * 1.01;
2016 rmax = resid.call(Tmax);
2024 }
while (!good_Tmax);
2025 if (rmin * rmax > 0 && std::abs(rmax) < std::abs(rmin)) {
2031 #if defined(ENABLE_CATCH)
2033 TEST_CASE(
"PD with T very large should yield error",
"[PDflash]") {
2035 double Tc = HEOS->T_critical();
2037 CHECK_THROWS(HEOS->update(
DmassP_INPUTS, 2, 5 * HEOS->p()));
2040 TEST_CASE(
"Stability testing",
"[stability]") {
2041 shared_ptr<HelmholtzEOSMixtureBackend> HEOS(
new HelmholtzEOSMixtureBackend(
strsplit(
"n-Propane&n-Butane&n-Pentane&n-Hexane",
'&')));
2042 std::vector<double> z(4);
2047 HEOS->set_mole_fractions(z);
2050 double TL = HEOS->T(), rhoL = HEOS->rhomolar();
2053 double TV = HEOS->T(), rhoV = HEOS->rhomolar();
2055 SECTION(
"Liquid (feed is stable)") {
2056 StabilityRoutines::StabilityEvaluationClass stability_tester(*HEOS);
2057 for (
double T = TL - 1; T >= 100; T -= 1) {
2058 stability_tester.set_TP(T, 101325);
2060 CHECK_NOTHROW(stability_tester.is_stable());
2063 SECTION(
"Vapor (feed is stable)") {
2064 StabilityRoutines::StabilityEvaluationClass stability_tester(*HEOS);
2065 for (
double T = TV + 1; T <= 500; T += 1) {
2066 stability_tester.set_TP(T, 101325);
2068 CHECK_NOTHROW(stability_tester.is_stable());
2071 SECTION(
"Two-phase (feed is unstable)") {
2072 StabilityRoutines::StabilityEvaluationClass stability_tester(*HEOS);
2073 stability_tester.set_TP((TV + TL) / 2.0, 101325);
2074 CHECK(stability_tester.is_stable() ==
false);
2078 TEST_CASE(
"Test critical points for methane + H2S",
"[critical_points]") {
2079 shared_ptr<HelmholtzEOSMixtureBackend> HEOS(
new HelmholtzEOSMixtureBackend(
strsplit(
"Methane&H2S",
'&')));
2081 double zz[] = {0.998, 0.97, 0.9475, 0.94, 0.93, 0.86, 0.85, 0.84, 0.75, 0.53, 0.52, 0.51, 0.49, 0.36, 0.24, 0.229, 0.09};
2082 int Npts[] = {2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 2, 2, 2, 1, 1, 1, 1};
2083 int imax =
sizeof(zz) /
sizeof(
double);
2085 for (
int i = 0; i < imax; ++i) {
2087 std::vector<double> z(2);
2090 HEOS->set_mole_fractions(z);
2092 std::vector<CriticalState> pts = HEOS->all_critical_points();
2093 CHECK(pts.size() == Npts[i]);
2097 TEST_CASE(
"Test critical points for nitrogen + ethane with HEOS",
"[critical_points]") {
2098 shared_ptr<HelmholtzEOSMixtureBackend> HEOS(
new HelmholtzEOSMixtureBackend(
strsplit(
"Nitrogen&Ethane",
'&')));
2099 std::vector<double> zz =
linspace(0.001, 0.999, 21);
2100 for (
int i = 0; i < static_cast<std::size_t>(zz.size()); ++i) {
2102 std::vector<double> z(2);
2105 HEOS->set_mole_fractions(z);
2107 std::vector<CriticalState> pts;
2108 CHECK_NOTHROW(pts = HEOS->all_critical_points());
2112 TEST_CASE(
"Test critical points for nitrogen + ethane with PR",
"[critical_points]") {
2113 shared_ptr<PengRobinsonBackend> HEOS(
new PengRobinsonBackend(
strsplit(
"Nitrogen&Ethane",
'&')));
2114 HEOS->set_binary_interaction_double(0, 1,
"kij", 0.0407);
2115 std::vector<double> zz =
linspace(0.001, 0.999, 21);
2116 for (
int i = 0; i < static_cast<std::size_t>(zz.size()); ++i) {
2118 std::vector<double> z(2);
2121 HEOS->set_mole_fractions(z);
2123 std::vector<CriticalState> pts;
2124 CHECK_NOTHROW(pts = HEOS->all_critical_points());