34 if (generate_SatL_and_SatV) {
35 bool SatLSatV =
false;
47 if (components.empty()) {
51 for (std::size_t i = 0; i < N; ++i) {
52 const std::string& alpha_type = components[i].alpha_type;
53 if (alpha_type !=
"default") {
54 const std::vector<double>& c = components[i].alpha_coeffs;
55 shared_ptr<AbstractCubicAlphaFunction> acaf;
56 if (alpha_type ==
"Twu") {
57 acaf.reset(
new TwuAlphaFunction(get_cubic()->a0_ii(i), c[0], c[1], c[2], get_cubic()->get_Tr() / get_cubic()->get_Tc()[i]));
58 }
else if (alpha_type ==
"MathiasCopeman" || alpha_type ==
"Mathias-Copeman") {
62 throw ValueError(
"alpha function is not understood");
64 cubic->set_alpha_function(i, acaf);
72 if (components.size() == 0) {
79 for (std::size_t i = 0; i < N; ++i) {
83 _components.push_back(fld);
89 if (!ParamName.compare(
"name")) {
91 }
else if (!ParamName.compare(
"aliases")) {
93 }
else if (!ParamName.compare(
"CAS") || !ParamName.compare(
"CAS_number")) {
95 }
else if (!ParamName.compare(
"formula")) {
97 }
else if (!ParamName.compare(
"ASHRAE34")) {
99 }
else if (!ParamName.compare(
"REFPROPName") || !ParamName.compare(
"REFPROP_name") || !ParamName.compare(
"REFPROPname")) {
101 }
else if (ParamName.find(
"BibTeX") == 0)
104 }
else if (ParamName.find(
"pure") == 0) {
110 }
else if (ParamName ==
"INCHI" || ParamName ==
"InChI" || ParamName ==
"INCHI_STRING") {
112 }
else if (ParamName ==
"INCHI_Key" || ParamName ==
"InChIKey" || ParamName ==
"INCHIKEY") {
114 }
else if (ParamName ==
"2DPNG_URL") {
116 }
else if (ParamName ==
"SMILES" || ParamName ==
"smiles") {
118 }
else if (ParamName ==
"CHEMSPIDER_ID") {
119 throw NotImplementedError(
"Parameter \"CHEMSPIDER_ID\" not available for cubic backends.");
120 }
else if (ParamName ==
"JSON") {
123 throw ValueError(
format(
"fluid parameter [%s] is invalid", ParamName.c_str()));
128 std::vector<std::string> out;
129 for (std::size_t i = 0; i < components.size(); ++i) {
130 out.push_back(components[i].name);
140 const std::vector<double> Tc = cubic->get_Tc(), pc = cubic->get_pc();
141 for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
142 T_r += mole_fractions[i] * Tc[i];
144 double v_c_Lmol = 2.14107171795 * (Tc[i] / pc[i] * 1000) + 0.00773144012514;
145 v_r += mole_fractions[i] * v_c_Lmol / 1000.0;
147 rhomolar_r = 1 / v_r;
155 double Tr_GERGlike, rhor_GERGlike;
156 get_linear_reducing_parameters(rhor_GERGlike, Tr_GERGlike);
157 R_delta *= rhor_GERGlike / rhomolar_reducing() * 5;
158 R_tau *= T_reducing() / Tr_GERGlike * 5;
164 double b = get_cubic()->bm_term(mole_fractions);
165 double v = 1 / (delta * rhomolar_reducing());
166 bool covolume_check = v < 1.1 * b;
168 return covolume_check;
187 double Tr_GERGlike, rhor_GERGlike;
188 get_linear_reducing_parameters(rhor_GERGlike, Tr_GERGlike);
189 delta0 *= rhor_GERGlike / rhomolar_reducing();
190 tau0 *= T_reducing() / Tr_GERGlike;
194 double tau = cubic->
get_Tr() / T;
195 double delta = rhomolar / cubic->
get_rhor();
196 return _rhomolar * gas_constant() * _T * (1 + delta * cubic->
alphar(tau, delta, this->get_mole_fractions_doubleref(), 0, 1));
201 _p = calc_pressure_nocache(_T, _rhomolar);
203 _phase = imposed_phase_index;
211 const std::vector<CoolPropDbl>& mole_fractions,
const CoolPropDbl& tau,
213 bool cache_values =
true;
214 HelmholtzDerivatives derivs = residual_helmholtz->all(*
this, mole_fractions, tau, delta, cache_values);
219 return derivs.alphar;
221 return derivs.dalphar_ddelta;
223 return derivs.d2alphar_ddelta2;
225 return derivs.d3alphar_ddelta3;
227 return derivs.d4alphar_ddelta4;
236 return derivs.dalphar_dtau;
238 return derivs.d2alphar_ddelta_dtau;
240 return derivs.d3alphar_ddelta2_dtau;
242 return derivs.d4alphar_ddelta3_dtau;
251 return derivs.d2alphar_dtau2;
253 return derivs.d3alphar_ddelta_dtau2;
255 return derivs.d4alphar_ddelta2_dtau2;
263 return derivs.d3alphar_dtau3;
265 return derivs.d4alphar_ddelta_dtau3;
273 return derivs.d4alphar_dtau4;
285 std::cout <<
format(
"%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
290 CoolPropDbl ld_value1 = value1, ld_value2 = value2;
291 pre_update(input_pair, ld_value1, ld_value2);
295 switch (input_pair) {
299 _rhomolar = solver_rho_Tp(value2 , value1 );
304 saturation(input_pair);
309 saturation(input_pair);
340 double am = cubic->
am_term(cubic->
get_Tr() / T, mole_fractions, 0);
341 double bm = cubic->
bm_term(mole_fractions);
351 double crho1 = R * T - p * (d1 + d2 + d3);
352 double crho2 = R * T * (d2 + d3) - p * (d1 * (d2 + d3) + d2 * d3) - am;
353 double crho3 = R * T * d2 * d3 - p * d1 * d2 * d3 - d1 * am;
356 solve_cubic(crho3, crho2, crho1, crho0, Nsolns, rho0, rho1, rho2);
357 sort3(rho0, rho1, rho2);
371 : ACB(ACB), inputs(inputs), imposed_variable(imposed_variable){};
375 double rho0 = -1, rho1 = -1, rho2 = -1, T, p;
379 p = imposed_variable;
382 T = imposed_variable;
395 double tau = T_r / T;
397 deltaL = rho2 / rho_r;
398 deltaV = rho0 / rho_r;
400 double DELTAgibbs = log(deltaV) - log(deltaL);
414 double tau = cubic->
get_Tr() / _T;
415 std::vector<double> x(1, 1);
416 double a = cubic->
am_term(tau, x, 0);
422 double crho4 = -
powInt(Delta_1 * Delta_2, 2) * R * _T *
powInt(b, 4) + a *
powInt(b, 3) * (Delta_1 + Delta_2);
424 -2 * ((Delta_1 * Delta_1 * Delta_2 + Delta_1 * Delta_2 * Delta_2) * R * _T *
powInt(b, 3) + a *
powInt(b, 2) * (Delta_1 + Delta_2 - 1));
425 double crho2 = ((-Delta_1 * Delta_1 - Delta_2 * Delta_2 - 4 * Delta_1 * Delta_2) * R * _T *
powInt(b, 2) + a * b * (Delta_1 + Delta_2 - 4));
426 double crho1 = -2 * (Delta_1 + Delta_2) * R * _T * b + 2 * a;
427 double crho0 = -R * _T;
428 double rho0, rho1, rho2, rho3;
430 solve_quartic(crho4, crho3, crho2, crho1, crho0, Nsoln, rho0, rho1, rho2, rho3);
431 std::vector<double> roots;
432 if (rho0 > 0 && 1 / rho0 > b) {
433 roots.push_back(rho0);
435 if (rho1 > 0 && 1 / rho1 > b) {
436 roots.push_back(rho1);
438 if (rho2 > 0 && 1 / rho2 > b) {
439 roots.push_back(rho2);
441 if (rho3 > 0 && 1 / rho3 > b) {
442 roots.push_back(rho3);
450 double rhoL = -1, rhoV = -1;
452 if (is_pure_or_pseudopure) {
454 double theta = -log10(_p / pc) * (1 / 0.7 - 1) / (acentric + 1);
455 double Ts_est = Tc / (theta + 1);
457 static std::string errstr;
469 if (is_pure_or_pseudopure) {
471 static std::string errstr;
474 std::vector<double> roots;
477 double neg_log10_pr = (acentric + 1) / (1 / 0.7 - 1) * (Tc / _T - 1);
478 double ps_est = pc * pow(10.0, -neg_log10_pr);
481 if (roots.size() == 2) {
482 double p0 = calc_pressure_nocache(_T, roots[0]);
483 double p1 = calc_pressure_nocache(_T, roots[1]);
488 if (p0 > 0 && p1 < pc) {
507 _rhomolar = 1 / (_Q / rhoV + (1 - _Q) / rhoL);
511 _rhomolar = solver_rho_Tp(T, p, 40000);
512 return static_cast<double>(_rhomolar);
516 double rho0 = 0, rho1 = 0, rho2 = 0, rho = -1;
517 rho_Tp_cubic(T, p, Nsoln, rho0, rho1, rho2);
520 }
else if (Nsoln == 3) {
526 }
else if (rho1 > 0) {
528 }
else if (rho2 > 0) {
536 throw ValueError(
"Specified phase is invalid");
539 if (p < p_critical()) {
540 add_transient_pure_state();
541 transient_pure_state->set_mole_fractions(this->mole_fractions);
542 transient_pure_state->update(
PQ_INPUTS, p, 0);
543 if (T > transient_pure_state->T()) {
544 double rhoV = transient_pure_state->saturated_vapor_keyed_output(
iDmolar);
546 if (rho0 > 0 && rho0 < rhoV) {
548 }
else if (rho1 > 0 && rho1 < rhoV) {
558 throw ValueError(
"Cubic has three roots, but phase not imposed and guess density not provided");
562 throw ValueError(
"Obtained neither 1 nor three roots");
564 if (is_pure_or_pseudopure) {
566 this->recalculate_singlephase_phase();
570 _phase = imposed_phase_index;
581 for (
unsigned int i = 0; i < N; ++i) {
582 summer += mole_fractions[i] * components[i].molemass;
588 const double value) {
590 if (i < 0 || i >= N) {
591 if (j < 0 || j >= N) {
592 throw ValueError(
format(
"Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N - 1));
594 throw ValueError(
format(
"Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
596 }
else if (j < 0 || j >= N) {
597 throw ValueError(
format(
"Index j [%d] is out of bounds. Must be between 0 and %d.", j, N - 1));
599 if (parameter ==
"kij" || parameter ==
"k_ij") {
600 get_cubic()->set_kij(i, j, value);
602 throw ValueError(
format(
"I don't know what to do with parameter [%s]", parameter.c_str()));
604 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
605 (*it)->set_binary_interaction_double(i, j, parameter, value);
610 if (i < 0 || i >= N) {
611 if (j < 0 || j >= N) {
612 throw ValueError(
format(
"Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N - 1));
614 throw ValueError(
format(
"Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
616 }
else if (j < 0 || j >= N) {
617 throw ValueError(
format(
"Index j [%d] is out of bounds. Must be between 0 and %d.", j, N - 1));
619 if (parameter ==
"kij" || parameter ==
"k_ij") {
620 return get_cubic()->get_kij(i, j);
622 throw ValueError(
format(
"I don't know what to do with parameter [%s]", parameter.c_str()));
627 get_cubic()->set_all_alpha_functions(donor->
get_cubic()->get_all_alpha_functions());
628 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
635 get_cubic()->set_kmat(donor->
get_cubic()->get_kmat());
636 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
643 this->copy_k(&donor);
646 this->set_alpha_from_components();
647 this->set_alpha0_from_components();
648 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
659 if (i < 0 || i >= N) {
660 throw ValueError(
format(
"Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
662 if (parameter ==
"MC" || parameter ==
"mc" || parameter ==
"Mathias-Copeman") {
663 get_cubic()->set_C_MC(i, c1, c2, c3);
664 }
else if (parameter ==
"TWU" || parameter ==
"Twu" || parameter ==
"twu") {
665 get_cubic()->set_C_Twu(i, c1, c2, c3);
667 throw ValueError(
format(
"I don't know what to do with parameter [%s]", parameter.c_str()));
669 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
677 if (i < 0 || i >= N) {
678 throw ValueError(
format(
"Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
681 if (parameter ==
"c" || parameter ==
"cm" || parameter ==
"c_m") {
682 get_cubic()->set_cm(value);
683 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
687 }
else if (parameter ==
"Q" || parameter ==
"Qk" || parameter ==
"Q_k") {
688 get_cubic()->set_Q_k(i, value);
689 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
694 throw ValueError(
format(
"I don't know what to do with parameter [%s]", parameter.c_str()));
699 if (i < 0 || i >= N) {
700 throw ValueError(
format(
"Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
703 if (parameter ==
"c" || parameter ==
"cm" || parameter ==
"c_m") {
704 return get_cubic()->get_cm();
705 }
else if (parameter ==
"Q" || parameter ==
"Qk" || parameter ==
"Q_k") {
706 return get_cubic()->get_Q_k(i);
708 throw ValueError(
format(
"I don't know what to do with parameter [%s]", parameter.c_str()));