11 # define _CRTDBG_MAP_ALLOC
12 # ifndef _CRT_SECURE_NO_WARNINGS
13 # define _CRT_SECURE_NO_WARNINGS
16 # include <sys/stat.h>
18 # include <sys/stat.h>
40 static int deriv_counter = 0;
48 if (fluid_names.size() == 1) {
67 std::vector<CoolPropFluid>
components(component_names.size());
68 for (
unsigned int i = 0; i <
components.size(); ++i) {
96 this->
N = components.size();
101 std::vector<std::vector<double>> ones(1, std::vector<double>(1, 1));
112 if (generate_SatL_and_SatV) {
122 if (mf.size() !=
N) {
123 throw ValueError(
format(
"size of mole fraction vector [%d] does not equal that of component vector [%d]", mf.size(),
N));
136 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it =
linked_states.begin(); it !=
linked_states.end(); ++it) {
137 it->get()->sync_linked_states(source);
149 if (mass_fractions.size() !=
N) {
150 throw ValueError(
format(
"size of mass fraction vector [%d] does not equal that of component vector [%d]", mass_fractions.size(),
N));
152 std::vector<CoolPropDbl> moles;
155 for (
unsigned int i = 0; i <
components.size(); ++i) {
156 tmp = mass_fractions[i] /
components[i].molar_mass();
157 moles.push_back(tmp);
161 for (std::vector<CoolPropDbl>::iterator it = moles.begin(); it != moles.end(); ++it) {
170 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it =
linked_states.begin(); it !=
linked_states.end(); ++it) {
172 it->get()->resize(
N);
197 if (!ParamName.compare(
"name")) {
199 }
else if (!ParamName.compare(
"aliases")) {
201 }
else if (!ParamName.compare(
"CAS") || !ParamName.compare(
"CAS_number")) {
203 }
else if (!ParamName.compare(
"formula")) {
205 }
else if (!ParamName.compare(
"ASHRAE34")) {
207 }
else if (!ParamName.compare(
"REFPROPName") || !ParamName.compare(
"REFPROP_name") || !ParamName.compare(
"REFPROPname")) {
209 }
else if (ParamName.find(
"BibTeX") == 0)
211 std::vector<std::string> parts =
strsplit(ParamName,
'-');
212 if (parts.size() != 2) {
213 throw ValueError(
format(
"Unable to parse BibTeX string %s", ParamName.c_str()));
215 std::string key = parts[1];
216 if (!key.compare(
"EOS")) {
218 }
else if (!key.compare(
"CP0")) {
220 }
else if (!key.compare(
"VISCOSITY")) {
222 }
else if (!key.compare(
"CONDUCTIVITY")) {
224 }
else if (!key.compare(
"ECS_LENNARD_JONES")) {
226 }
else if (!key.compare(
"ECS_VISCOSITY_FITS")) {
228 }
else if (!key.compare(
"ECS_CONDUCTIVITY_FITS")) {
230 }
else if (!key.compare(
"SURFACE_TENSION")) {
232 }
else if (!key.compare(
"MELTING_LINE")) {
237 }
else if (ParamName.find(
"pure") == 0) {
243 }
else if (ParamName ==
"INCHI" || ParamName ==
"InChI" || ParamName ==
"INCHI_STRING") {
244 return cpfluid.
InChI;
245 }
else if (ParamName ==
"INCHI_Key" || ParamName ==
"InChIKey" || ParamName ==
"INCHIKEY") {
247 }
else if (ParamName ==
"2DPNG_URL") {
249 }
else if (ParamName ==
"SMILES" || ParamName ==
"smiles") {
251 }
else if (ParamName ==
"CHEMSPIDER_ID") {
253 }
else if (ParamName ==
"JSON") {
256 throw ValueError(
format(
"fluid parameter [%s] is invalid", ParamName.c_str()));
262 if (i < 0 || i >=
N) {
263 if (j < 0 || j >=
N) {
264 throw ValueError(
format(
"Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j,
N-1));
266 throw ValueError(
format(
"Index i [%d] is out of bounds. Must be between 0 and %d.", i,
N-1));
268 }
else if (j < 0 || j >=
N) {
269 throw ValueError(
format(
"Index j [%d] is out of bounds. Must be between 0 and %d.", j,
N-1));
271 if (model ==
"linear") {
273 double gammaT = 0.5 * (Tc1 + Tc2) / sqrt(Tc1 * Tc2);
275 double gammaV = 4.0 * (1 / rhoc1 + 1 / rhoc2) / pow(pow(rhoc1, -1.0 / 3.0) + pow(rhoc2, -1.0 / 3.0), 3);
280 }
else if (model ==
"Lorentz-Berthelot") {
286 throw ValueError(
format(
"mixing rule [%s] is not understood", model.c_str()));
291 const double value) {
293 if (i < 0 || i >=
N) {
294 if (j < 0 || j >=
N) {
295 throw ValueError(
format(
"Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j,
N-1));
297 throw ValueError(
format(
"Index i [%d] is out of bounds. Must be between 0 and %d.", i,
N-1));
299 }
else if (j < 0 || j >=
N) {
300 throw ValueError(
format(
"Index j [%d] is out of bounds. Must be between 0 and %d.", j,
N-1));
302 if (parameter ==
"Fij") {
306 Reducing->set_binary_interaction_double(i, j, parameter, value);
309 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it =
linked_states.begin(); it !=
linked_states.end(); ++it) {
310 it->get()->set_binary_interaction_double(i, j, parameter, value);
316 if (i < 0 || i >=
N) {
317 if (j < 0 || j >=
N) {
318 throw ValueError(
format(
"Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j,
N-1));
320 throw ValueError(
format(
"Index i [%d] is out of bounds. Must be between 0 and %d.", i,
N-1));
322 }
else if (j < 0 || j >=
N) {
323 throw ValueError(
format(
"Index j [%d] is out of bounds. Must be between 0 and %d.", j,
N-1));
325 if (parameter ==
"Fij") {
328 return Reducing->get_binary_interaction_double(i, j, parameter);
337 const std::string& value) {
339 if (i < 0 || i >=
N) {
340 if (j < 0 || j >=
N) {
341 throw ValueError(
format(
"Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j,
N-1));
343 throw ValueError(
format(
"Index i [%d] is out of bounds. Must be between 0 and %d.", i,
N-1));
345 }
else if (j < 0 || j >=
N) {
346 throw ValueError(
format(
"Index j [%d] is out of bounds. Must be between 0 and %d.", j,
N-1));
348 if (parameter ==
"function") {
352 throw ValueError(
format(
"Cannot process this string parameter [%s] in set_binary_interaction_string", parameter.c_str()));
355 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it =
linked_states.begin(); it !=
linked_states.end(); ++it) {
356 it->get()->set_binary_interaction_string(i, j, parameter, value);
366 if (EOS_name ==
"SRK" || EOS_name ==
"Peng-Robinson") {
378 shared_ptr<AbstractCubic> ac;
379 if (EOS_name ==
"SRK") {
380 ac.reset(
new SRK(Tc, pc, acentric, R));
385 ac->set_rhor(rhomolarc);
387 }
else if (EOS_name ==
"XiangDeiters") {
405 if (this->
SatL)
SatL->change_EOS(i, EOS_name);
406 if (this->
SatV)
SatV->change_EOS(i, EOS_name);
442 if (!state.compare(
"hs_anchor")) {
444 }
else if (!state.compare(
"max_sat_T")) {
446 }
else if (!state.compare(
"max_sat_p")) {
448 }
else if (!state.compare(
"reducing")) {
450 }
else if (!state.compare(
"critical")) {
452 }
else if (!state.compare(
"triple_liquid")) {
454 }
else if (!state.compare(
"triple_vapor")) {
457 throw ValueError(
format(
"This state [%s] is invalid to calc_state", state.c_str()));
460 if (!state.compare(
"critical")) {
471 throw ValueError(
"acentric factor cannot be calculated for mixtures");
483 for (
unsigned int i = 0; i <
components.size(); ++i) {
492 for (
unsigned int i = 0; i <
components.size(); ++i) {
499 if (param ==
iP && given ==
iT) {
503 return components[0].ancillaries.pL.evaluate(value);
505 return components[0].ancillaries.pV.evaluate(value);
507 }
else if (param ==
iT && given ==
iP) {
511 return components[0].ancillaries.pL.invert(value);
513 return components[0].ancillaries.pV.invert(value);
515 }
else if (param ==
iDmolar && given ==
iT) {
519 return components[0].ancillaries.rhoL.evaluate(value);
521 return components[0].ancillaries.rhoV.evaluate(value);
523 }
else if (param ==
iT && given ==
iDmolar) {
527 return components[0].ancillaries.rhoL.invert(value);
529 return components[0].ancillaries.rhoV.invert(value);
532 return components[0].ancillaries.surface_tension.evaluate(value);
546 return components[0].ancillaries.melting_line.evaluate(param, given, value);
554 return components[0].ancillaries.surface_tension.evaluate(
T());
556 throw ValueError(
format(
"surface tension is only defined within the two-phase region; Try PQ or QT inputs"));
565 switch (
components[0].transport.viscosity_dilute.type) {
592 format(
"dilute viscosity type [%d] is invalid for fluid %s",
components[0].transport.viscosity_dilute.type,
name().c_str()));
605 switch (
components[0].transport.viscosity_initial.type) {
609 initial_density = eta_dilute * B_eta_initial * rho;
622 switch (
components[0].transport.viscosity_higher_order.type) {
652 format(
"higher order viscosity type [%d] is invalid for fluid %s",
components[0].transport.viscosity_dilute.type,
name().c_str()));
655 return initial_density + residual;
660 CoolPropDbl dilute = 0, initial_density = 0, residual = 0, critical = 0;
662 return dilute + initial_density + residual + critical;
687 throw ValueError(
format(
"Viscosity model is not available for this fluid"));
694 std::vector<std::string> names(1, fluid_name);
763 throw ValueError(
"calc_viscosity_contributions invalid for mixtures");
779 throw ValueError(
format(
"Thermal conductivity model is not available for this fluid"));
786 std::vector<std::string>
name(1, fluid_name);
814 throw ValueError(
format(
"hardcoded conductivity type [%d] is invalid for fluid %s",
846 format(
"dilute conductivity type [%d] is invalid for fluid %s",
components[0].transport.conductivity_dilute.type,
name().c_str()));
871 format(
"critical conductivity type [%d] is invalid for fluid %s",
components[0].transport.viscosity_dilute.type,
name().c_str()));
874 throw ValueError(
"calc_conductivity_contributions invalid for mixtures");
881 switch (
components[0].transport.conductivity_residual.type) {
890 format(
"residual conductivity type [%d] is invalid for fluid %s",
components[0].transport.conductivity_residual.type,
name().c_str()));
892 return lambda_residual;
896 CoolPropDbl dilute = 0, initial_density = 0, residual = 0, critical = 0;
898 return dilute + initial_density + residual + critical;
934 for (
unsigned int i = 0; i <
components.size(); ++i) {
941 for (
unsigned int i = 0; i <
components.size(); ++i) {
956 if (!
SatL)
throw ValueError(
"The saturated liquid state has not been set.");
957 return SatL->keyed_output(key);
961 if (!
SatV)
throw ValueError(
"The saturated vapor state has not been set.");
962 return SatV->keyed_output(key);
966 if (type ==
"Joule-Thomson") {
969 }
else if (type ==
"Joule-Inversion") {
972 }
else if (type ==
"Ideal") {
975 }
else if (type ==
"Boyle") {
983 std::vector<std::string> out;
984 for (std::size_t i = 0; i <
components.size(); ++i) {
991 throw ValueError(
format(
"For now, calc_ODP is only valid for pure and pseudo-pure fluids, %d components",
components.size()));
1002 throw ValueError(
format(
"For now, calc_GWP20 is only valid for pure and pseudo-pure fluids, %d components",
components.size()));
1013 throw ValueError(
format(
"For now, calc_GWP100 is only valid for pure and pseudo-pure fluids, %d components",
components.size()));
1024 throw ValueError(
format(
"For now, calc_GWP500 is only valid for pure and pseudo-pure fluids, %d components",
components.size()));
1036 if (critpts.size() == 1) {
1038 return critpts[0].T;
1040 throw ValueError(
format(
"critical point finding routine found %d critical points", critpts.size()));
1049 if (critpts.size() == 1) {
1051 return critpts[0].p;
1053 throw ValueError(
format(
"critical point finding routine found %d critical points", critpts.size()));
1062 if (critpts.size() == 1) {
1064 return critpts[0].rhomolar;
1066 throw ValueError(
format(
"critical point finding routine found %d critical points", critpts.size()));
1080 throw ValueError(
"calc_pmax_sat not yet defined for mixtures");
1086 double Tmax_sat =
components[0].EOS().max_sat_T.T;
1096 throw ValueError(
"calc_Tmax_sat not yet defined for mixtures");
1102 Tmin_satL =
components[0].EOS().sat_min_liquid.T;
1103 Tmin_satV =
components[0].EOS().sat_min_vapor.T;
1106 throw ValueError(
"calc_Tmin_sat not yet defined for mixtures");
1112 pmin_satL =
components[0].EOS().sat_min_liquid.p;
1113 pmin_satV =
components[0].EOS().sat_min_vapor.p;
1116 throw ValueError(
"calc_pmin_sat not yet defined for mixtures");
1125 for (
unsigned int i = 0; i <
components.size(); ++i) {
1132 for (
unsigned int i = 0; i <
components.size(); ++i) {
1139 for (
unsigned int i = 0; i <
components.size(); ++i) {
1151 throw ValueError(
format(
"The molar density of %f mol/m3 is below the minimum of %f mol/m3",
rhomolar, rhomolar_min));
1155 throw ValueError(
format(
"The temperature of %f K is below the minimum of %f K",
T, T_min));
1167 bool optional_checks =
false;
1186 this->
_T = HEOS.
T();
1188 this->
_p = HEOS.
p();
1190 this->
_Q = HEOS.
Q();
1214 throw ValueError(
"Mole fractions must be set");
1230 std::cout <<
format(
"%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
1235 CoolPropDbl ld_value1 = value1, ld_value2 = value2;
1236 pre_update(input_pair, ld_value1, ld_value2);
1240 switch (input_pair) {
1345 return mass_fractions;
1351 std::cout <<
format(
"%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
1356 CoolPropDbl ld_value1 = value1, ld_value2 = value2;
1357 pre_update(input_pair, ld_value1, ld_value2);
1361 switch (input_pair) {
1394 throw ValueError(
"rhomolar is less than zero");
1397 throw ValueError(
"rhomolar is not a valid number");
1400 if (optional_checks) {
1437 saturation_called =
false;
1446 if (
_p > psat_max) {
1495 throw ValueError(
"supercritical pressure but other invalid for now");
1500 else if (
_p >=
components[0].EOS().ptriple * 0.9999 &&
_p <= psat_max) {
1507 bool definitely_two_phase =
false;
1518 if (
_T < Tm - 0.001) {
1519 throw ValueError(
format(
"For now, we don't support T [%g K] below Tmelt(p) [%g K]",
_T, Tm));
1526 if (
_T <
Tmin() - 0.001) {
1527 throw ValueError(
format(
"For now, we don't support T [%g K] below Tmin(saturation) [%g K]",
_T,
Tmin()));
1535 if (value > T_vap) {
1539 }
else if (value < T_liq) {
1563 if (value > h_vap + h_vap_error_band) {
1567 }
else if (value < h_liq - h_liq_error_band) {
1571 }
else if (value > h_liq + h_liq_error_band && value < h_vap - h_vap_error_band) {
1572 definitely_two_phase =
true;
1588 if (value > s_vap + s_vap_error_band) {
1592 }
else if (value < s_liq - s_liq_error_band) {
1596 }
else if (value > s_liq + s_liq_error_band && value < s_vap - s_vap_error_band) {
1597 definitely_two_phase =
true;
1616 CoolPropDbl u_liq_error_band = 1.5 * h_liq_error_band;
1617 CoolPropDbl u_vap_error_band = 1.5 * h_vap_error_band;
1620 if (value > u_vap + u_vap_error_band) {
1624 }
else if (value < u_liq - u_liq_error_band) {
1628 }
else if (value > u_liq + u_liq_error_band && value < u_vap - u_vap_error_band) {
1629 definitely_two_phase =
true;
1639 if (!definitely_two_phase) {
1646 if (value < rho_vap) {
1649 }
else if (value > rho_liq) {
1659 throw ValueError(
"possibly two-phase inputs not supported for mixtures for now");
1673 saturation_called =
true;
1692 Q = (1 / value - 1 / HEOS.
SatL->rhomolar()) / (1 / HEOS.
SatV->rhomolar() - 1 / HEOS.
SatL->rhomolar());
1695 Q = (value - HEOS.
SatL->smolar()) / (HEOS.
SatV->smolar() - HEOS.
SatL->smolar());
1698 Q = (value - HEOS.
SatL->hmolar()) / (HEOS.
SatV->hmolar() - HEOS.
SatL->hmolar());
1701 Q = (value - HEOS.
SatL->umolar()) / (HEOS.
SatV->umolar() - HEOS.
SatL->umolar());
1719 }
else if (Q > 1 + 1e-9) {
1732 }
else if (
_p <
components[0].EOS().ptriple * 0.9999) {
1740 throw NotImplementedError(
format(
"For now, we don't support p [%g Pa] below ptriple [%g Pa] when T [%g] is less than Tmin [%g]",
1748 throw ValueError(
format(
"The pressure [%g Pa] cannot be used in p_phase_determination",
_p));
1757 double call(
double T) {
1760 double dTdp_along_sat =
1761 HEOS->
T() * (1 / HEOS->
SatV->rhomolar() - 1 / HEOS->
SatL->rhomolar()) / (HEOS->
SatV->hmolar() - HEOS->
SatL->hmolar());
1768 Residual resid(*HEOS_copy);
1771 double v2 = resid.call(tripleV.
T);
1792 double call(
double T) {
1795 double dTdp_along_sat =
1796 HEOS->
T() * (1 / HEOS->
SatV->rhomolar() - 1 / HEOS->
SatL->rhomolar()) / (HEOS->
SatV->hmolar() - HEOS->
SatL->hmolar());
1803 Residualhmax residhmax(*HEOS_copy);
1814 throw ValueError(
format(
"value to T_phase_determination_pure_or_pseudopure is invalid"));
1818 if (_T < _crit.T && _p >
_crit.
p) {
1848 throw ValueError(
format(
"T=Tcrit; invalid input for other to T_phase_determination_pure_or_pseudopure"));
1862 if (value < p_vap) {
1866 }
else if (value > p_liq) {
1883 throw ValueError(
"Two-phase inputs not supported for pseudo-pure for now");
1896 if (value < rho_vap) {
1899 }
else if (value > rho_liq) {
1904 double Qanc = (1 / value - 1 /
static_cast<double>(
_rhoLanc))
1907 if (value > 0.95 * rho_liq || value < 1.05 * rho_vap) {
1918 }
else if (Qanc > 1.01) {
1930 throw ValueError(
format(
"The saturation properties are needed in T_phase_determination_pure_or_pseudopure"));
1939 if (value >
SatV->calc_smolar()) {
1950 if (value >
SatV->calc_hmolar()) {
1960 if (value >
SatV->calc_umolar()) {
1970 throw ValueError(
format(
"invalid input for other to T_phase_determination_pure_or_pseudopure"));
1987 if (value > HEOS.
SatL->p() * (1e-6 + 1)) {
1991 }
else if (value < HEOS.SatV->
p() * (1 - 1e-6)) {
1997 format(
"Saturation pressure [%g Pa] corresponding to T [%g K] is within 1e-4 %% of given p [%Lg Pa]", HEOS.
SatL->p(),
_T, value));
2003 Q = (1 / value - 1 / HEOS.
SatL->rhomolar()) / (1 / HEOS.
SatV->rhomolar() - 1 / HEOS.
SatL->rhomolar());
2006 Q = (value - HEOS.
SatL->smolar()) / (HEOS.
SatV->smolar() - HEOS.
SatL->smolar());
2009 Q = (value - HEOS.
SatL->hmolar()) / (HEOS.
SatV->hmolar() - HEOS.
SatL->hmolar());
2012 Q = (value - HEOS.
SatL->umolar()) / (HEOS.
SatV->umolar() - HEOS.
SatL->umolar());
2091 throw ValueError(
"supercritical temp but other invalid for now");
2100 dT_dtau = -pow(T, 2) / Tr, R = HEOS->
gas_constant(), delta = rho / rhor, tau = Tr / T;
2183 int nTau = 0, nDelta = 1;
2217 double second_deriv(
double rhomolar) {
2219 return R_u *
T /
POW2(rhor)
2223 dpdrho_resid resid(
this,
T,
p);
2227 light =
Halley(resid, 1e-6, 1e-8, 100);
2228 double d2pdrho2__constT = resid.deriv(light);
2229 if (d2pdrho2__constT > 0) {
2233 }
catch (std::exception& e) {
2235 std::cout << e.what() << std::endl;
2244 for (std::size_t counter = 0; counter <= 100; counter++) {
2246 double curvature = resid.deriv(rho);
2247 if (curvature > 0) {
2258 for (
double omega = 0.7; omega > 0; omega -= 0.2) {
2260 resid.options.add_number(
"omega", omega);
2261 heavy =
Halley(resid, rhomax, 1e-8, 100);
2262 double d2pdrho2__constT = resid.deriv(heavy);
2263 if (d2pdrho2__constT < 0) {
2268 }
catch (std::exception& e) {
2270 std::cout << e.what() << std::endl;
2279 double rho = rhomax;
2280 for (std::size_t counter = 0; counter <= 100; counter++) {
2282 double curvature = resid.deriv(rho);
2283 if (curvature < 0 || this->
p() < 0) {
2293 if (light > 0 && heavy > 0) {
2300 else if (light < 0 && heavy < 0) {
2301 double dpdrho_min = resid.call(1e-10);
2302 double dpdrho_max = resid.call(rhomax);
2303 if (dpdrho_max * dpdrho_min > 0) {
2324 rhor(
HEOS->get_reducing_state().rhomolar),
2331 return (peos -
p) /
p;
2376 double rho_liq = -1, rho_vap = -1;
2377 if (
p > p_at_rhomax_stationary) {
2379 for (; counter <= 10; counter++) {
2382 if (p_at_rhomax <
p) {
2383 rhomolar_max *= 1.05;
2392 if (
p < p_at_rhomin_stationary) {
2397 if (rho_vap > 0 && rho_liq > 0) {
2399 if (std::abs(rho_vap - rho_liq) < 1e-10) {
2406 if (gibbsmolar_liq < gibbsmolar_vap) {
2412 }
else if (rho_vap < 0 && rho_liq > 0) {
2415 }
else if (rho_vap > 0 && rho_liq < 0) {
2440 if (rhomolar_guess < 0)
2447 if (rhomolar_guess < 0 || !
ValidNumber(rhomolar_guess))
2461 throw ValueError(
"Liquid density is invalid");
2463 }
catch (std::exception&) {
2494 if (dpdrho < 0 || d2pdrho2 < 0) {
2502 if (dpdrho < 0 || d2pdrho2 > 0) {
2509 }
catch (std::exception& e) {
2523 throw ValueError(
format(
"solver_rho_Tp was unable to find a solution for T=%10Lg, p=%10Lg, with guess value %10Lg with error: %s",
T,
p,
2524 rhomolar_guess, e.what()));
2530 for (std::size_t i = 0; i <
components.size(); ++i) {
2532 CoolPropDbl m_i = 0.480 + 1.574 * acentric_i - 0.176 * pow(acentric_i, 2);
2536 CoolPropDbl a_i = 0.42747 * pow(R_u * Tci, 2) / pci * pow(1 + m_i * (1 - sqrt(
T / Tci)), 2);
2538 for (std::size_t j = 0; j <
components.size(); ++j) {
2540 CoolPropDbl m_j = 0.480 + 1.574 * acentric_j - 0.176 * pow(acentric_j, 2);
2542 CoolPropDbl a_j = 0.42747 * pow(R_u * Tcj, 2) / pcj * pow(1 + m_j * (1 - sqrt(
T / Tcj)), 2);
2562 solve_cubic(1, -1, A - B - B * B, -A * B, Nsolns, Z0, Z1, Z2);
2573 if (rhomolar0 > 0 && rhomolar1 <= 0 && rhomolar2 <= 0) {
2576 if (rhomolar0 <= 0 && rhomolar1 > 0 && rhomolar2 <= 0) {
2579 if (rhomolar0 <= 0 && rhomolar1 <= 0 && rhomolar2 > 0) {
2594 throw ValueError(
"Bad phase to solver_rho_Tp_SRK");
2632 return R_u *
T * (1 +
tau * (da0_dTau + dar_dTau) +
delta * dar_dDelta);
2636 std::cout <<
format(
"HelmholtzEOSMixtureBackend::calc_hmolar: 2phase: %d T: %g rhomomolar: %g",
isTwoPhase(),
_T,
_rhomolar) << std::endl;
2638 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for the two-phase properties"));
2680 return R_u * (
tau * (da0_dTau + dar_dTau) - a0 - ar);
2684 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for the two-phase properties"));
2706 _smolar = R_u * (
_tau.
pt() * (da0_dTau + dar_dTau) - a0 - ar);
2724 return R_u *
T *
tau * (da0_dTau + dar_dTau);
2728 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for the two-phase properties"));
2768 return static_cast<double>(
_cvmolar);
2785 * (-pow(
_tau.
pt(), 2) * (d2ar_dTau2 + d2a0_dTau2)
2787 / (1 + 2 *
_delta.
pt() * dar_dDelta + pow(
_delta.
pt(), 2) * d2ar_dDelta2));
2789 return static_cast<double>(
_cpmolar);
2801 return R_u * (1 + (-pow(
_tau.
pt(), 2)) * d2a0_dTau2);
2806 return SatL->speed_sound();
2808 return SatV->speed_sound();
2810 throw ValueError(
format(
"Speed of sound is not defined for two-phase states because it depends on the distribution of phases."));
2830 - pow(1 +
_delta.
pt() * dar_dDelta -
_delta.
pt() *
_tau.
pt() * d2ar_dDelta_dTau, 2) / (pow(
_tau.
pt(), 2) * (d2ar_dTau2 + d2a0_dTau2))));
2853 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for the two-phase properties"));
2879 for (std::size_t i = 0; i <
components.size(); ++i) {
2895 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for the two-phase properties"));
2929 double dna0_dni__const_T_V_nj =
2931 return gas_constant() *
T() * (dna0_dni__const_T_V_nj + dnar_dni__const_T_V_nj);
2952 throw ValueError(
"Mole fractions must be set before calling calc_reducing_state");
2961 bool cache_values =
true;
2982 bool cache_values =
false;
2984 return derivs.
get(nTau, nDelta);
2991 throw ValueError(
"No alpha0 derivatives are available");
3002 double taustar = Tc / Tr *
tau, deltastar = rhor / rhomolarc *
delta;
3003 if (nTau == 0 && nDelta == 0) {
3004 val = E.
base0(taustar, deltastar);
3005 }
else if (nTau == 0 && nDelta == 1) {
3007 }
else if (nTau == 1 && nDelta == 0) {
3009 }
else if (nTau == 0 && nDelta == 2) {
3011 }
else if (nTau == 1 && nDelta == 1) {
3013 }
else if (nTau == 2 && nDelta == 0) {
3015 }
else if (nTau == 0 && nDelta == 3) {
3017 }
else if (nTau == 1 && nDelta == 2) {
3019 }
else if (nTau == 2 && nDelta == 1) {
3021 }
else if (nTau == 3 && nDelta == 0) {
3026 val *= pow(rhor / rhomolarc, nDelta);
3027 val /= pow(Tr / Tc, nTau);
3030 throw ValueError(
format(
"calc_alpha0_deriv_nocache returned invalid number with inputs nTau: %d, nDelta: %d, tau: %Lg, delta: %Lg", nTau,
3041 for (
unsigned int i = 0; i <
N; ++i) {
3046 tau_i = T_ci *
tau / Tr;
3047 delta_i =
delta * rhor / rho_ci;
3053 if (nTau == 0 && nDelta == 0) {
3056 }
else if (nTau == 0 && nDelta == 1) {
3058 }
else if (nTau == 1 && nDelta == 0) {
3060 }
else if (nTau == 0 && nDelta == 2) {
3061 summer +=
mole_fractions[i] * Rratio * pow(rhor / rho_ci, 2) *
components[i].EOS().d2alpha0_dDelta2(tau_i, delta_i);
3062 }
else if (nTau == 1 && nDelta == 1) {
3063 summer +=
mole_fractions[i] * Rratio * rhor / rho_ci * T_ci / Tr *
components[i].EOS().d2alpha0_dDelta_dTau(tau_i, delta_i);
3064 }
else if (nTau == 2 && nDelta == 0) {
3136 const int nTau = 0, nDelta = 0;
3140 const int nTau = 0, nDelta = 1;
3144 const int nTau = 1, nDelta = 0;
3148 const int nTau = 0, nDelta = 2;
3152 const int nTau = 1, nDelta = 1;
3156 const int nTau = 2, nDelta = 0;
3160 const int nTau = 0, nDelta = 3;
3164 const int nTau = 1, nDelta = 2;
3168 const int nTau = 2, nDelta = 1;
3172 const int nTau = 3, nDelta = 0;
3181 if (Of1 ==
iT && Wrt1 ==
iP) {
3183 }
else if (Of1 ==
iP && Wrt1 ==
iT) {
3184 return 1 / dTdP_sat;
3187 else if (Wrt1 ==
iT) {
3191 else if (Wrt1 ==
iP) {
3199 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for calc_first_saturation_deriv"));
3205 if (Of1 ==
iT && Wrt1 ==
iP) {
3207 }
else if (Of1 ==
iP && Wrt1 ==
iT) {
3208 return 1 / dTdP_sat;
3211 else if (Wrt1 ==
iT) {
3215 else if (Wrt1 ==
iP) {
3223 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for calc_second_saturation_deriv"));
3224 if (Wrt1 ==
iP && Wrt2 ==
iP) {
3239 CoolPropDbl ddT_dTdp_along_sat_constp = (DELTAh * (
_T * dDELTAv_dT_constp + DELTAv) -
_T * DELTAv * dDELTAh_dT_constp) /
POW2(DELTAh);
3240 CoolPropDbl ddp_dTdp_along_sat_constT = (DELTAh * (
_T * dDELTAv_dp_constT) -
_T * DELTAv * dDELTAh_dp_constT) /
POW2(DELTAh);
3242 double ddp_dydpsigma = d2ydp2_constT + dydT_constp * ddp_dTdp_along_sat_constT + d2ydTdp * dTdp_along_sat;
3243 double ddT_dydpsigma = d2ydTdp + dydT_constp * ddT_dTdp_along_sat_constp + d2ydT2_constp * dTdp_along_sat;
3244 return ddp_dydpsigma + ddT_dydpsigma * dTdp_along_sat;
3246 throw ValueError(
format(
"Currently, only possible to take second saturation derivative w.r.t. P (both times)"));
3252 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for calc_second_two_phase_deriv"));
3267 CoolPropDbl numerator = 1 /
SatV->keyed_output(rho_key) - 1 /
SatL->keyed_output(rho_key);
3269 CoolPropDbl dnumerator = -1 /
POW2(
SatV->keyed_output(rho_key)) * drhoV_dp_sat + 1 /
POW2(
SatL->keyed_output(rho_key)) * drhoL_dp_sat;
3270 CoolPropDbl ddenominator = dhV_dp_sat - dhL_dp_sat;
3271 CoolPropDbl d_dvdh_dp__consth = (denominator * dnumerator - numerator * ddenominator) /
POW2(denominator);
3272 return -
POW2(
rhomolar()) * d_dvdh_dp__consth + dv_dh_constp * (-2 *
rhomolar()) * drhomolar_dp__consth;
3274 && ((Wrt1 ==
iHmass && Constant1 ==
iP && Wrt2 ==
iP && Constant2 ==
iHmass)
3275 || (Wrt2 ==
iHmass && Constant2 ==
iP && Wrt1 ==
iP && Constant1 ==
iHmass))) {
3287 CoolPropDbl numerator = 1 /
SatV->keyed_output(rho_key) - 1 /
SatL->keyed_output(rho_key);
3289 CoolPropDbl dnumerator = -1 /
POW2(
SatV->keyed_output(rho_key)) * drhoV_dp_sat + 1 /
POW2(
SatL->keyed_output(rho_key)) * drhoL_dp_sat;
3290 CoolPropDbl ddenominator = dhV_dp_sat - dhL_dp_sat;
3291 CoolPropDbl d_dvdh_dp__consth = (denominator * dnumerator - numerator * ddenominator) /
POW2(denominator);
3292 return -
POW2(rho) * d_dvdh_dp__consth + dv_dh_constp * (-2 * rho) * drho_dp__consth;
3298 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for calc_first_two_phase_deriv"));
3312 CoolPropDbl dvdp_h = dvL_dp + dxdp_h * (1 /
SatV->rhomolar() - 1 /
SatL->rhomolar()) +
Q() * (dvV_dp - dvL_dp);
3323 CoolPropDbl dvdp_h = dvL_dp + dxdp_h * (1 /
SatV->rhomass() - 1 /
SatL->rhomass()) +
Q() * (dvV_dp - dvL_dp);
3326 throw ValueError(
"These inputs are not supported to calc_first_two_phase_deriv");
3333 bool drho_dh__p =
false;
3334 bool drho_dp__h =
false;
3335 bool rho_spline =
false;
3355 throw ValueError(
"These inputs are not supported to calc_first_two_phase_deriv");
3358 if (!this->
SatL || !this->
SatV)
throw ValueError(
format(
"The saturation properties are needed for calc_first_two_phase_deriv_splined"));
3371 Liq->update_DmolarT_direct(
SatL->rhomolar(),
SatL->T());
3386 CoolPropDbl Abracket = (2 * rho_liq - 2 * rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end));
3388 CoolPropDbl b = 3 /
POW2(Delta_end) * (-rho_liq + rho_end) - 1 / Delta_end * (drho_dh_end + 2 * drho_dh_liq__constp);
3422 CoolPropDbl d_Delta_end_dp__consth = x_end * (dhV_dp_sat - dhL_dp_sat);
3426 CoolPropDbl d_Abracket_dp_consth = (2 * drhoL_dp_sat - 2 * drho_dp_end + Delta_end * (d2rhodhdp_liq + d2rhodhdp_end)
3427 + d_Delta_end_dp__consth * (drho_dh_liq__constp + drho_dh_end));
3428 CoolPropDbl da_dp = 1 /
POW3(Delta_end) * d_Abracket_dp_consth + Abracket * (-3 /
POW4(Delta_end) * d_Delta_end_dp__consth);
3429 CoolPropDbl db_dp = -6 /
POW3(Delta_end) * d_Delta_end_dp__consth * (rho_end - rho_liq) + 3 /
POW2(Delta_end) * (drho_dp_end - drhoL_dp_sat)
3430 + (1 /
POW2(Delta_end) * d_Delta_end_dp__consth) * (drho_dh_end + 2 * drho_dh_liq__constp)
3431 - (1 / Delta_end) * (d2rhodhdp_end + 2 * d2rhodhdp_liq);
3436 (3 * a *
POW2(Delta) + 2 * b * Delta + c) * d_Delta_dp__consth +
POW3(Delta) * da_dp +
POW2(Delta) * db_dp + Delta * dc_dp + dd_dp;
3439 throw ValueError(
"Something went wrong in HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv_splined");
3449 Eigen::MatrixXd Lstar, Mstar;
3451 std::vector<double> call(
const std::vector<double>& tau_delta) {
3457 std::vector<double> o(2);
3458 o[0] = Lstar.determinant();
3459 o[1] = Mstar.determinant();
3462 std::vector<std::vector<double>> Jacobian(
const std::vector<double>& x) {
3463 std::size_t
N = x.size();
3464 std::vector<std::vector<double>> J(
N, std::vector<double>(
N, 0));
3470 J[0][0] = (adjL * dLdTau).trace();
3471 J[0][1] = (adjL * dLdDelta).trace();
3472 J[1][0] = (adjM * dMdTau).trace();
3473 J[1][1] = (adjM * dMdDelta).trace();
3477 std::vector<std::vector<double>> numerical_Jacobian(
const std::vector<double>& x) {
3478 std::size_t
N = x.size();
3479 std::vector<double> rplus, rminus, xp;
3480 std::vector<std::vector<double>> J(
N, std::vector<double>(
N, 0));
3482 double epsilon = 0.0001;
3485 for (std::size_t i = 0; i <
N; ++i) {
3489 xp[i] -= 2 * epsilon;
3492 for (std::size_t j = 0; j <
N; ++j) {
3493 J[j][i] = (rplus[j] - rminus[j]) / (2 * epsilon);
3496 std::cout << J[0][0] <<
" " << J[0][1] << std::endl;
3497 std::cout << J[1][0] <<
" " << J[1][1] << std::endl;
3502 std::vector<double> x, tau_delta(2);
3547 _deriv = (adjL * dLdTau).trace();
3554 dadjLstardTau = adjugate_derivative(Lstar, dLstardTau);
3555 _second_deriv = (dLstardTau * dadjLstardTau + adjL * d2LstardTau2).trace();
3595 tau_new =
tau +
R_tau * cos(theta);
3603 double tau_new, delta_new;
3611 double L1 =
Lstar.determinant();
3620 return -
R_tau * sin(theta) * dL1_dtau +
R_delta * cos(theta) * dL1_ddelta;
3626 for (
int i = 0; i < 300; ++i) {
3668 double p_MPa =
HEOS.
p() / 1e6;
3671 double tau_new, delta_new;
3694 this->tau = tau_new;
3695 this->delta = delta_new;
3699 this->spinodal_values.
tau.push_back(tau_new);
3700 this->spinodal_values.
delta.push_back(delta_new);
3701 this->spinodal_values.
M1.push_back(M1);
3709 L1star = Lstar.determinant();
3710 M1star = Mstar.determinant();
3721 rapidjson::Document doc;
3723 rapidjson::Value& val = doc;
3724 std::vector<std::vector<DepartureFunctionPointer>>&
mat =
critical_state->residual_helmholtz->Excess.DepartureFunctionMatrix;
3725 if (
mat.size() > 0) {
3726 mat[0][1]->phi.to_json(val, doc);
3735 double delta0 = _HUGE, tau0 = _HUGE;
3736 critical_state->get_critical_point_starting_values(delta0, tau0);
3743 resid_L0.
call(tau0);
3745 while (resid_L0.
deriv(tau0) > 0 && bump_count < 3) {
3749 double tau_L0 =
Halley(resid_L0, tau0, 1e-10, 100);
3756 double R_delta = 0, R_tau = 0;
3768 const double rhomolar_guess) {
3771 if (w.size() != z.size()) {
3772 throw ValueError(
format(
"Trial composition vector size [%d] is not the same as bulk composition [%d]", w.size(), z.size()));
3781 for (std::size_t i = 0; i < w.size(); ++i) {
3790 bool find_critical_points =
false;
3795 for (std::size_t i = 0; i <
components.size(); ++i) {
3797 if (!reference_state.compare(
"IIR")) {
3798 if (HEOS.
Ttriple() > 273.15) {
3799 throw ValueError(
format(
"Cannot use IIR reference state; Ttriple [%Lg] is greater than 273.15 K", HEOS.
Ttriple()));
3804 double deltah = HEOS.
hmass() - 200000;
3805 double deltas = HEOS.
smass() - 1000;
3811 std::cout <<
format(
"set offsets to %0.15g and %0.15g\n", delta_a1, delta_a2);
3813 }
else if (!reference_state.compare(
"ASHRAE")) {
3814 if (HEOS.
Ttriple() > 233.15) {
3815 throw ValueError(
format(
"Cannot use ASHRAE reference state; Ttriple [%Lg] is greater than than 233.15 K", HEOS.
Ttriple()));
3820 double deltah = HEOS.
hmass() - 0;
3821 double deltas = HEOS.
smass() - 0;
3827 std::cout <<
format(
"set offsets to %0.15g and %0.15g\n", delta_a1, delta_a2);
3829 }
else if (!reference_state.compare(
"NBP")) {
3831 throw ValueError(
format(
"Cannot use NBP reference state; p_triple [%Lg Pa] is greater than than 101325 Pa", HEOS.
p_triple()));
3836 double deltah = HEOS.
hmass() - 0;
3837 double deltas = HEOS.
smass() - 0;
3843 std::cout <<
format(
"set offsets to %0.15g and %0.15g\n", delta_a1, delta_a2);
3845 }
else if (!reference_state.compare(
"DEF")) {
3847 }
else if (!reference_state.compare(
"RESET")) {
3850 throw ValueError(
format(
"reference state string is invalid: [%s]", reference_state.c_str()));
3861 for (std::size_t i = 0; i <
components.size(); ++i) {
3867 double deltah = HEOS.
hmolar() - hmolar0;
3868 double deltas = HEOS.
smolar() - smolar0;
3876 const std::string& ref) {
3886 double f = (HEOS->name() ==
"Water" || HEOS->name() ==
"CarbonDioxide") ? 1.00001 : 1.0;
3908 if (!HEOS->is_pure()) {