CoolProp 8.0.0
An open-source fluid property and humid air property database
CubicBackend.cpp
Go to the documentation of this file.
1#include "CubicBackend.h"
2
3#include <cmath>
4#include "CoolProp/numerics/Solvers.h"
5#include "CoolProp/Configuration.h"
6#include "Backends/Helmholtz/VLERoutines.h"
7#include "Backends/Helmholtz/MixtureDerivatives.h"
8#include <Eigen/Dense>
9
10void CoolProp::AbstractCubicBackend::setup(bool generate_SatL_and_SatV) {
11 N = cubic->get_Tc().size();
12
13 // Set the pure fluid flag
14 is_pure_or_pseudopure = (N == 1);
15
16 // Resize the vectors
17 resize(N);
18
19 // Reset the residual Helmholtz energy class
20 residual_helmholtz = std::make_shared<CubicResidualHelmholtz>(this);
21 // If pure, set the mole fractions to be unity
22 if (is_pure_or_pseudopure) {
23 mole_fractions = std::vector<CoolPropDbl>(1, 1.0);
24 } else {
25 mole_fractions.clear();
26 }
27 // Now set the reducing function for the mixture
28 Reducing = std::make_shared<ConstantReducingFunction>(cubic->get_Tr(), cubic->get_rhor());
29
30 // Set the alpha function based on the components in use
31 set_alpha_from_components();
32
33 // Set the ideal-gas helmholtz energy based on the components in use;
34 set_alpha0_from_components();
35
36 // Top-level class can hold copies of the base saturation classes,
37 // saturation classes cannot hold copies of the saturation classes
38 if (generate_SatL_and_SatV) {
39 bool SatLSatV = false;
40 SatL.reset(this->get_copy(SatLSatV));
41 SatL->specify_phase(iphase_liquid);
42 linked_states.push_back(SatL);
43 SatV.reset(this->get_copy(SatLSatV));
44 SatV->specify_phase(iphase_gas);
45 linked_states.push_back(SatV);
46 }
47}
48
49void CoolProp::AbstractCubicBackend::set_alpha_from_components() {
51 if (components.empty()) {
52 return;
53 }
54
55 for (std::size_t i = 0; i < N; ++i) {
56 const std::string& alpha_type = components[i].alpha_type;
57 if (alpha_type != "default") {
58 const std::vector<double>& c = components[i].alpha_coeffs;
59 shared_ptr<AbstractCubicAlphaFunction> acaf;
60 if (alpha_type == "Twu") {
61 acaf = std::make_shared<TwuAlphaFunction>(get_cubic()->a0_ii(i), c[0], c[1], c[2], get_cubic()->get_Tr() / get_cubic()->get_Tc()[i]);
62 } else if (alpha_type == "MathiasCopeman" || alpha_type == "Mathias-Copeman") {
63 acaf.reset(
64 new MathiasCopemanAlphaFunction(get_cubic()->a0_ii(i), c[0], c[1], c[2], get_cubic()->get_Tr() / get_cubic()->get_Tc()[i]));
65 } else {
66 throw ValueError("alpha function is not understood");
67 }
68 cubic->set_alpha_function(i, acaf);
69 }
70 }
71}
72
73void CoolProp::AbstractCubicBackend::set_alpha0_from_components() {
74
75 // If empty so far (e.g., for SatL and SatV)
76 if (components.size() == 0) {
77 return;
78 }
79
80 // Get the vector of CoolProp fluids from the base class
81 std::vector<CoolPropFluid>& _components = HelmholtzEOSMixtureBackend::get_components();
82
83 for (std::size_t i = 0; i < N; ++i) {
84 CoolPropFluid fld;
85 fld.EOSVector.emplace_back();
86 fld.EOS().alpha0 = components[i].alpha0;
87 _components.push_back(fld);
88 }
89}
90
91std::string CoolProp::AbstractCubicBackend::fluid_param_string(const std::string& ParamName) {
92 CoolProp::CubicLibrary::CubicsValues cpfluid = components[0];
93 if (!ParamName.compare("name")) {
94 return cpfluid.name;
95 } else if (!ParamName.compare("aliases")) {
96 return strjoin(cpfluid.aliases, get_config_string(LIST_STRING_DELIMITER));
97 } else if (!ParamName.compare("CAS") || !ParamName.compare("CAS_number")) {
98 return cpfluid.CAS;
99 } else if (!ParamName.compare("formula")) {
100 throw NotImplementedError("Parameter \"formula\" not available for cubic backends.");
101 } else if (!ParamName.compare("ASHRAE34")) {
102 throw NotImplementedError("Parameter \"ASHRAE34\" not available for cubic backends.");
103 } else if (!ParamName.compare("REFPROPName") || !ParamName.compare("REFPROP_name") || !ParamName.compare("REFPROPname")) {
104 throw NotImplementedError("Parameter \"REFPROPName\" not available for cubic backends.");
105 } else if (ParamName.find("BibTeX") == 0) // Starts with "BibTeX"
106 {
107 throw NotImplementedError("BibTeX parameters not available for cubic backends.");
108 } else if (ParamName.find("pure") == 0) {
109 if (is_pure()) {
110 return "true";
111 } else {
112 return "false";
113 }
114 } else if (ParamName == "INCHI" || ParamName == "InChI" || ParamName == "INCHI_STRING") {
115 throw NotImplementedError("Parameter \"INCHI\" not available for cubic backends.");
116 } else if (ParamName == "INCHI_Key" || ParamName == "InChIKey" || ParamName == "INCHIKEY") {
117 throw NotImplementedError("Parameter \"INCHI_Key\" not available for cubic backends.");
118 } else if (ParamName == "2DPNG_URL") {
119 throw NotImplementedError("Parameter \"2DPNG_URL\" not available for cubic backends.");
120 } else if (ParamName == "SMILES" || ParamName == "smiles") {
121 throw NotImplementedError("Parameter \"SMILES\" not available for cubic backends.");
122 } else if (ParamName == "CHEMSPIDER_ID") {
123 throw NotImplementedError("Parameter \"CHEMSPIDER_ID\" not available for cubic backends.");
124 } else if (ParamName == "JSON") {
125 return CubicLibrary::get_fluid_as_JSONstring(cpfluid.name);
126 } else {
127 throw ValueError(format("fluid parameter [%s] is invalid", ParamName.c_str()));
128 }
129}
130
131std::vector<std::string> CoolProp::AbstractCubicBackend::calc_fluid_names() {
132 std::vector<std::string> out;
133 out.reserve(components.size());
134 for (auto& component : components) {
135 out.push_back(component.name);
136 }
137 return out;
138}
139
140void CoolProp::AbstractCubicBackend::get_linear_reducing_parameters(double& rhomolar_r, double& T_r) {
141 // In the case of models where the reducing temperature is not a function of composition (SRK, PR, etc.),
142 // we need to use an appropriate value for T_r and v_r, here we use a linear weighting
143 T_r = 0;
144 double v_r = 0;
145 const std::vector<double> Tc = cubic->get_Tc(), pc = cubic->get_pc();
146 for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
147 T_r += mole_fractions[i] * Tc[i];
148 // Curve fit from all the pure fluids in CoolProp (thanks to recommendation of A. Kazakov)
149 double v_c_Lmol = 2.14107171795 * (Tc[i] / pc[i] * 1000) + 0.00773144012514; // [L/mol]
150 v_r += mole_fractions[i] * v_c_Lmol / 1000.0;
151 }
152 rhomolar_r = 1 / v_r;
153}
154
155void CoolProp::AbstractCubicBackend::get_critical_point_search_radii(double& R_delta, double& R_tau) {
156 // Get the starting values from the base class
157 CoolProp::HelmholtzEOSMixtureBackend::get_critical_point_search_radii(R_delta, R_tau);
158
159 // Now we scale them to get the appropriate search radii
160 double Tr_GERGlike = NAN, rhor_GERGlike = NAN;
161 get_linear_reducing_parameters(rhor_GERGlike, Tr_GERGlike);
162 R_delta *= rhor_GERGlike / rhomolar_reducing() * 5;
163 R_tau *= T_reducing() / Tr_GERGlike * 5;
164}
165
166bool CoolProp::AbstractCubicBackend::get_critical_is_terminated(double& delta, double& tau) {
167 // If the volume is less than the mixture covolume, stop. The mixture covolume is the
168 // smallest volume that is physically allowed for a cubic EOS
169 double b = get_cubic()->bm_term(mole_fractions); // [m^3/mol]
170 double v = 1 / (delta * rhomolar_reducing()); //[m^3/mol]
171 bool covolume_check = v < 1.1 * b;
172
173 return covolume_check;
174}
175
176void CoolProp::AbstractCubicBackend::get_critical_point_starting_values(double& delta0, double& tau0) {
177
178 // Get the starting values from the base class
179 CoolProp::HelmholtzEOSMixtureBackend::get_critical_point_starting_values(delta0, tau0);
180
181 // The starting tau and delta values can be thought of as being given with the
182 // GERG reducing values, or tau0 = Tr_GERG/T = 0.xxx and delta0 = rho/rhor_GERG = 0.xxx
183 // Then we need to multiply by
184 //
185 // Tr_cubic/Tr_GERGlike & rhor_GERGlike/rhor_cubic
186 //
187 // to get shifted values:
188 //
189 // delta0 = rho/rhor_GERG*(rhor_GERGlike/rhor_cubic)
190 // tau0 = Tr_GERG/T*(Tr_cubic/Tr_GERGlike)
191 //
192 double Tr_GERGlike = NAN, rhor_GERGlike = NAN;
193 get_linear_reducing_parameters(rhor_GERGlike, Tr_GERGlike);
194 delta0 *= rhor_GERGlike / rhomolar_reducing();
195 tau0 *= T_reducing() / Tr_GERGlike;
196}
197CoolPropDbl CoolProp::AbstractCubicBackend::calc_pressure_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
198 AbstractCubic* cubic = get_cubic().get();
199 double tau = cubic->get_Tr() / T;
200 double delta = rhomolar / cubic->get_rhor();
201 return _rhomolar * gas_constant() * _T * (1 + delta * cubic->alphar(tau, delta, this->get_mole_fractions_doubleref(), 0, 1));
202}
203void CoolProp::AbstractCubicBackend::update_DmolarT() {
204 if (this->imposed_phase_index != iphase_not_imposed) {
205 _p = calc_pressure_nocache(_T, _rhomolar);
206 _Q = -1;
207 _phase = imposed_phase_index;
208 return;
209 }
210 // Mixtures don't have a cubic superancillary; fall through to the HEOS path so the
211 // user sees the same diagnostic as before for unsupported mixture DmolarT input.
212 if (!is_pure_or_pseudopure || get_superanc_eos_code() == CubicSuperAncillary::UNKNOWN_CODE) {
213 HelmholtzEOSMixtureBackend::update(DmolarT_INPUTS, this->_rhomolar, this->_T);
214 return;
215 }
216 // Pure SRK/PR: use the Chebyshev superancillary to bracket the dome at T, then pick
217 // the matching EOS root. Above the superancillary Tmax we are unambiguously single-
218 // phase and just evaluate the EOS directly.
219 const double Tmax = calc_superanc_Tmax();
220 if (_T >= Tmax) {
221 _p = calc_pressure_nocache(_T, _rhomolar);
222 _Q = -1;
223 recalculate_singlephase_phase();
224 return;
225 }
226 const double rhoL_sat = calc_saturation_ancillary(iDmolar, 0, iT, _T);
227 const double rhoV_sat = calc_saturation_ancillary(iDmolar, 1, iT, _T);
228 if (_rhomolar >= rhoL_sat) {
229 _p = calc_pressure_nocache(_T, _rhomolar);
230 _Q = -1;
231 _phase = iphase_liquid;
232 } else if (_rhomolar <= rhoV_sat) {
233 _p = calc_pressure_nocache(_T, _rhomolar);
234 _Q = -1;
235 _phase = iphase_gas;
236 } else {
237 // Inside the dome: pressure is the saturation pressure and quality follows from
238 // the lever rule, 1/rho = (1-Q)/rhoL + Q/rhoV.
239 const double rho = _rhomolar;
240 _p = calc_saturation_ancillary(iP, 0, iT, _T);
241 _Q = (rhoV_sat * (rhoL_sat - rho)) / (rho * (rhoL_sat - rhoV_sat));
242 _phase = iphase_twophase;
243 SatL->update_TDmolarP_unchecked(_T, rhoL_sat, _p);
244 SatV->update_TDmolarP_unchecked(_T, rhoV_sat, _p);
245 }
246};
247
248CoolPropDbl CoolProp::AbstractCubicBackend::calc_alphar_deriv_nocache(const int nTau, const int nDelta,
249 const std::vector<CoolPropDbl>& mole_fractions, const CoolPropDbl& tau,
250 const CoolPropDbl& delta) {
251 bool cache_values = true;
252 HelmholtzDerivatives derivs = residual_helmholtz->all(*this, mole_fractions, tau, delta, cache_values);
253 switch (nTau) {
254 case 0: {
255 switch (nDelta) {
256 case 0:
257 return derivs.alphar;
258 case 1:
259 return derivs.dalphar_ddelta;
260 case 2:
261 return derivs.d2alphar_ddelta2;
262 case 3:
263 return derivs.d3alphar_ddelta3;
264 case 4:
265 return derivs.d4alphar_ddelta4;
266 default:
267 throw ValueError(format("nDelta (%d) is invalid", nDelta));
268 }
269 break;
270 }
271 case 1: {
272 switch (nDelta) {
273 case 0:
274 return derivs.dalphar_dtau;
275 case 1:
276 return derivs.d2alphar_ddelta_dtau;
277 case 2:
278 return derivs.d3alphar_ddelta2_dtau;
279 case 3:
280 return derivs.d4alphar_ddelta3_dtau;
281 default:
282 throw ValueError(format("nDelta (%d) is invalid", nDelta));
283 }
284 break;
285 }
286 case 2: {
287 switch (nDelta) {
288 case 0:
289 return derivs.d2alphar_dtau2;
290 case 1:
291 return derivs.d3alphar_ddelta_dtau2;
292 case 2:
293 return derivs.d4alphar_ddelta2_dtau2;
294 default:
295 throw ValueError(format("nDelta (%d) is invalid", nDelta));
296 }
297 }
298 case 3: {
299 switch (nDelta) {
300 case 0:
301 return derivs.d3alphar_dtau3;
302 case 1:
303 return derivs.d4alphar_ddelta_dtau3;
304 default:
305 throw ValueError(format("nDelta (%d) is invalid", nDelta));
306 }
307 }
308 case 4: {
309 switch (nDelta) {
310 case 0:
311 return derivs.d4alphar_dtau4;
312 default:
313 throw ValueError(format("nDelta (%d) is invalid", nDelta));
314 }
315 }
316 default:
317 throw ValueError(format("nTau (%d) is invalid", nTau));
318 }
319}
320
321void CoolProp::AbstractCubicBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) {
322 if (get_debug_level() > 10) {
323 std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
324 get_input_pair_short_desc(input_pair).c_str(), value1, value2)
325 << '\n';
326 }
327
328 // Mass-quality input pair on a true mixture: solve iteratively for Qmolar
329 // before delegating to the molar-pair flash. The inherited HEOS override of
330 // calc_phase_molar_masses (using SatL/SatV->mole_fractions and components[i].molar_mass())
331 // works for cubic backends since they share the same SatL/SatV machinery.
332 if (CoolProp::is_Qmass_pair(input_pair) && mole_fractions.size() > 1) {
333 update_Qmass_pair(input_pair, value1, value2);
334 return;
335 }
336
337 CoolPropDbl ld_value1 = value1, ld_value2 = value2;
338 pre_update(input_pair, ld_value1, ld_value2);
339 value1 = ld_value1;
340 value2 = ld_value2;
341
342 switch (input_pair) {
343 case PT_INPUTS:
344 _p = value1;
345 _T = value2;
346 if (is_pure_or_pseudopure || imposed_phase_index != iphase_not_imposed) {
347 _rhomolar = solver_rho_Tp(value2 /*T*/, value1 /*p*/);
348 } else {
349 HelmholtzEOSMixtureBackend::update(PT_INPUTS, value1, value2);
350 }
351 break;
352 case QT_INPUTS:
353 _Q = value1;
354 _T = value2;
355 saturation(input_pair);
356 break;
357 case PQ_INPUTS:
358 _p = value1;
359 _Q = value2;
360 saturation(input_pair);
361 break;
362 case DmolarT_INPUTS:
363 _rhomolar = value1;
364 _T = value2;
365 update_DmolarT();
366 break;
367 case HmolarT_INPUTS:
368 case SmolarT_INPUTS:
369 case TUmolar_INPUTS:
370 case DmolarP_INPUTS:
371 case DmolarHmolar_INPUTS:
372 case DmolarSmolar_INPUTS:
373 case DmolarUmolar_INPUTS:
374 case HmolarP_INPUTS:
375 case PSmolar_INPUTS:
376 case PUmolar_INPUTS:
377 case HmolarSmolar_INPUTS:
378 case QSmolar_INPUTS:
379 case HmolarQ_INPUTS:
380 case DmolarQ_INPUTS:
381 HelmholtzEOSMixtureBackend::update(input_pair, value1, value2);
382 break;
383 default:
384 throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
385 }
386
387 post_update();
388}
389
390void CoolProp::AbstractCubicBackend::rho_Tp_cubic(CoolPropDbl T, CoolPropDbl p, int& Nsolns, double& rho0, double& rho1, double& rho2) {
391 AbstractCubic* cubic = get_cubic().get();
392 double R = cubic->get_R_u();
393 double am = cubic->am_term(cubic->get_Tr() / T, mole_fractions, 0);
394 double bm = cubic->bm_term(mole_fractions);
395 double cm = cubic->cm_term();
396
397 // Introducing new variables to simplify the equation:
398 double d1 = cm - bm;
399 double d2 = cm + cubic->get_Delta_1() * bm;
400 double d3 = cm + cubic->get_Delta_2() * bm;
401
402 // Cubic coefficients:
403 double crho0 = -p;
404 double crho1 = R * T - p * (d1 + d2 + d3);
405 double crho2 = R * T * (d2 + d3) - p * (d1 * (d2 + d3) + d2 * d3) - am;
406 double crho3 = R * T * d2 * d3 - p * d1 * d2 * d3 - d1 * am;
407
408 // Solving the cubic:
409 solve_cubic(crho3, crho2, crho1, crho0, Nsolns, rho0, rho1, rho2);
410 sort3(rho0, rho1, rho2);
411 return;
412}
413
414class SaturationResidual : public CoolProp::FuncWrapper1D
415{
416 public:
417 CoolProp::AbstractCubicBackend* ACB = nullptr;
418 CoolProp::input_pairs inputs = CoolProp::INPUT_PAIR_INVALID;
419 double imposed_variable = _HUGE;
420 double deltaL = _HUGE, deltaV = _HUGE;
421
422 SaturationResidual() = default;
423 SaturationResidual(CoolProp::AbstractCubicBackend* ACB, CoolProp::input_pairs inputs, double imposed_variable)
424 : ACB(ACB), inputs(inputs), imposed_variable(imposed_variable) {};
425
426 double call(double value) override {
427 int Nsolns = 0;
428 double rho0 = -1, rho1 = -1, rho2 = -1, T = NAN, p = NAN;
429
430 if (inputs == CoolProp::PQ_INPUTS) {
431 T = value;
432 p = imposed_variable;
433 } else if (inputs == CoolProp::QT_INPUTS) {
434 p = value;
435 T = imposed_variable;
436 } else {
437 throw CoolProp::ValueError("Cannot have something other than PQ_INPUTS or QT_INPUTS here");
438 }
439
440 // Calculate the liquid and vapor densities
441 ACB->rho_Tp_cubic(T, p, Nsolns, rho0, rho1, rho2);
442
443 // -----------------------------------------------------
444 // Calculate the difference in Gibbs between the phases
445 // -----------------------------------------------------
446 AbstractCubic* cubic = ACB->get_cubic().get();
447 double rho_r = cubic->get_rhor(), T_r = cubic->get_Tr();
448 double tau = T_r / T;
449 // There are three density solutions, we know the highest is the liquid, the lowest is the vapor
450 deltaL = rho2 / rho_r;
451 deltaV = rho0 / rho_r;
452 // From alpha0; all terms that are only a function of temperature drop out since TL=TV
453 double DELTAgibbs = log(deltaV) - log(deltaL);
454 // From alphar;
455 DELTAgibbs += (cubic->alphar(tau, deltaV, ACB->get_mole_fractions_doubleref(), 0, 0)
456 - cubic->alphar(tau, deltaL, ACB->get_mole_fractions_doubleref(), 0, 0));
457 // From delta*dalphar_dDelta
458 DELTAgibbs += (deltaV * cubic->alphar(tau, deltaV, ACB->get_mole_fractions_doubleref(), 0, 1)
459 - deltaL * cubic->alphar(tau, deltaL, ACB->get_mole_fractions_doubleref(), 0, 1));
460 return DELTAgibbs;
461 };
462};
463
464std::vector<double> CoolProp::AbstractCubicBackend::spinodal_densities() {
466 AbstractCubic* cubic = get_cubic().get();
467 double tau = cubic->get_Tr() / _T;
468 std::vector<double> x(1, 1);
469 double a = cubic->am_term(tau, x, 0);
470 double b = cubic->bm_term(x);
471 double R = cubic->get_R_u();
472 double Delta_1 = cubic->get_Delta_1();
473 double Delta_2 = cubic->get_Delta_2();
474
475 double crho4 = -powInt(Delta_1 * Delta_2, 2) * R * _T * powInt(b, 4) + a * powInt(b, 3) * (Delta_1 + Delta_2);
476 double crho3 =
477 -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));
478 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));
479 double crho1 = -2 * (Delta_1 + Delta_2) * R * _T * b + 2 * a;
480 double crho0 = -R * _T;
481 double rho0 = NAN, rho1 = NAN, rho2 = NAN, rho3 = NAN;
482 int Nsoln = 0;
483 solve_quartic(crho4, crho3, crho2, crho1, crho0, Nsoln, rho0, rho1, rho2, rho3);
484 std::vector<double> roots;
485 if (rho0 > 0 && 1 / rho0 > b) {
486 roots.push_back(rho0);
487 }
488 if (rho1 > 0 && 1 / rho1 > b) {
489 roots.push_back(rho1);
490 }
491 if (rho2 > 0 && 1 / rho2 > b) {
492 roots.push_back(rho2);
493 }
494 if (rho3 > 0 && 1 / rho3 > b) {
495 roots.push_back(rho3);
496 }
497 return roots;
498}
499
500void CoolProp::AbstractCubicBackend::saturation(CoolProp::input_pairs inputs) {
501 AbstractCubic* cubic = get_cubic().get();
502 double Tc = cubic->get_Tc()[0], pc = cubic->get_pc()[0], acentric = cubic->get_acentric()[0];
503 double rhoL = -1, rhoV = -1;
504 if (inputs == PQ_INPUTS) {
505 if (is_pure_or_pseudopure) {
506 // Estimate temperature from the acentric factor relationship
507 double theta = -log10(_p / pc) * (1 / 0.7 - 1) / (acentric + 1);
508 double Ts_est = Tc / (theta + 1);
509 SaturationResidual resid(this, inputs, _p);
510 static std::string errstr;
511 double Ts = CoolProp::Secant(resid, Ts_est, -0.1, 1e-10, 100);
512 _T = Ts;
513 rhoL = resid.deltaL * cubic->get_Tr();
514 rhoV = resid.deltaV * cubic->get_Tr();
515 this->SatL->update(DmolarT_INPUTS, rhoL, _T);
516 this->SatV->update(DmolarT_INPUTS, rhoV, _T);
517 } else {
518 HelmholtzEOSMixtureBackend::update(PQ_INPUTS, _p, _Q);
519 return;
520 }
521 } else if (inputs == QT_INPUTS) {
522 if (is_pure_or_pseudopure) {
523 SaturationResidual resid(this, inputs, _T);
524 static std::string errstr;
525 // ** Spinodal densities is disabled for now because it is VERY slow :(
526 // std::vector<double> roots = spinodal_densities();
527 std::vector<double> roots;
528
529 // Estimate pressure from the acentric factor relationship
530 double neg_log10_pr = (acentric + 1) / (1 / 0.7 - 1) * (Tc / _T - 1);
531 double ps_est = pc * pow(10.0, -neg_log10_pr);
532
533 double ps = NAN;
534 if (roots.size() == 2) {
535 double p0 = calc_pressure_nocache(_T, roots[0]);
536 double p1 = calc_pressure_nocache(_T, roots[1]);
537 if (p1 < p0) {
538 std::swap(p0, p1);
539 }
540 //ps = CoolProp::BoundedSecant(resid, p0, p1, pc, -0.01*ps_est, 1e-5, 100);
541 if (p0 > 0 && p1 < pc) {
542 ps = CoolProp::Brent(resid, p0 * 1.0001, p1 * 0.9999, DBL_EPSILON, 1e-10, 100);
543 } else {
544 ps = CoolProp::BoundedSecant(resid, ps_est, 1e-10, pc, -0.01 * ps_est, 1e-5, 100);
545 }
546 } else {
547 ps = CoolProp::BoundedSecant(resid, ps_est, 1e-10, pc, -0.01 * ps_est, 1e-5, 100);
548 }
549
550 _p = ps;
551 rhoL = resid.deltaL * cubic->get_Tr();
552 rhoV = resid.deltaV * cubic->get_Tr();
553 this->SatL->update(DmolarT_INPUTS, rhoL, _T);
554 this->SatV->update(DmolarT_INPUTS, rhoV, _T);
555 } else {
556 HelmholtzEOSMixtureBackend::update(QT_INPUTS, _Q, _T);
557 return;
558 }
559 }
560 _rhomolar = 1 / (_Q / rhoV + (1 - _Q) / rhoL);
561 _phase = iphase_twophase;
562}
563CoolPropDbl CoolProp::AbstractCubicBackend::solver_rho_Tp_global(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomolar_max) {
564 int Nsoln = 0;
565 double rho0 = 0, rho1 = 0, rho2 = 0;
566 rho_Tp_cubic(T, p, Nsoln, rho0, rho1, rho2);
567
568 double rho;
569 if (Nsoln == 1) {
570 rho = rho0;
571 } else if (Nsoln == 3) {
572 if (imposed_phase_index != iphase_not_imposed) {
573 // Delegate to solver_rho_Tp which handles imposed phase
574 return solver_rho_Tp(T, p);
575 }
576 // Three roots, no imposed phase: pick the root with the lowest
577 // Gibbs energy. This avoids the expensive p_critical() call that
578 // triggers all_critical_points() for many-component mixtures.
579 double rho_vap = (rho0 > 0) ? rho0 : rho1;
580 double rho_liq = rho2;
581 // Guard against a non-physical (<= 0) candidate root: its Gibbs energy is
582 // garbage/NaN, and `g_liq < g_vap` is false for a NaN g_vap, which would
583 // silently select the non-physical root (CoolProp-1tbe.8 finding 3A).
584 bool vap_ok = rho_vap > 0;
585 bool liq_ok = rho_liq > 0;
586 if (vap_ok && liq_ok) {
587 double g_vap = calc_gibbsmolar_nocache(T, rho_vap);
588 double g_liq = calc_gibbsmolar_nocache(T, rho_liq);
589 // Prefer the lower-Gibbs root; if one Gibbs is non-finite, keep the finite one.
590 rho = (ValidNumber(g_liq) && (!ValidNumber(g_vap) || g_liq < g_vap)) ? rho_liq : rho_vap;
591 } else if (vap_ok || liq_ok) {
592 rho = vap_ok ? rho_vap : rho_liq;
593 } else {
594 throw ValueError("solver_rho_Tp_global: no positive density root among the three cubic roots");
595 }
596 } else {
597 throw ValueError("Obtained neither 1 nor three roots");
598 }
599
600 _rhomolar = rho;
601 update_DmolarT_direct(rho, T);
602 _Q = -1;
603 if (imposed_phase_index != iphase_not_imposed) {
604 _phase = imposed_phase_index;
605 } else {
606 _phase = (rho < rhomolar_reducing()) ? iphase_gas : iphase_liquid;
607 }
608 return rho;
609}
610CoolPropDbl CoolProp::AbstractCubicBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, phases phase) {
611 int Nsoln = 0;
612 double rho0 = 0, rho1 = 0, rho2 = 0, rho = -1;
613 rho_Tp_cubic(T, p, Nsoln, rho0, rho1, rho2); // Densities are sorted in increasing order
614 if (Nsoln == 1) {
615 rho = rho0;
616 } else if (Nsoln == 3) {
617 if (phase == iphase_gas || phase == iphase_supercritical_gas) {
618 if (rho0 > 0) {
619 rho = rho0;
620 } else if (rho1 > 0) {
621 rho = rho1;
622 } else if (rho2 > 0) {
623 rho = rho2;
624 } else {
625 throw CoolProp::ValueError(format("Unable to find gaseous density for T: %g K, p: %g Pa", T, p));
626 }
627 } else if (phase == iphase_liquid || phase == iphase_supercritical_liquid) {
628 rho = rho2;
629 } else {
630 throw ValueError("Specified phase is invalid");
631 }
632 } else {
633 throw ValueError("Obtained neither 1 nor three roots");
634 }
635 _phase = phase;
636 _Q = -1;
637 return rho;
638}
639
640CoolPropDbl CoolProp::AbstractCubicBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rho_guess) {
641 if (imposed_phase_index != iphase_not_imposed) {
642 return solver_rho_Tp(T, p, imposed_phase_index);
643 }
644 int Nsoln = 0;
645 double rho0 = 0, rho1 = 0, rho2 = 0, rho = -1;
646 rho_Tp_cubic(T, p, Nsoln, rho0, rho1, rho2); // Densities are sorted in increasing order
647 if (Nsoln == 1) {
648 rho = rho0;
649 } else if (Nsoln == 3) {
650 if (p < p_critical()) {
651 add_transient_pure_state();
652 transient_pure_state->set_mole_fractions(this->mole_fractions);
653 transient_pure_state->update(PQ_INPUTS, p, 0);
654 if (T > transient_pure_state->T()) {
655 double rhoV = transient_pure_state->saturated_vapor_keyed_output(iDmolar);
656 // Gas
657 if (rho0 > 0 && rho0 < rhoV) {
658 rho = rho0;
659 } else if (rho1 > 0 && rho1 < rhoV) {
660 rho = rho1;
661 } else {
662 throw CoolProp::ValueError(format("Unable to find gaseous density for T: %g K, p: %g Pa", T, p));
663 }
664 } else {
665 // Liquid
666 rho = rho2;
667 }
668 } else {
669 throw ValueError("Cubic has three roots, but phase not imposed and guess density not provided");
670 }
671 } else {
672 throw ValueError("Obtained neither 1 nor three roots");
673 }
674 if (is_pure_or_pseudopure) {
675 this->recalculate_singlephase_phase();
676 } else {
677 if (rho < rhomolar_reducing()) {
678 _phase = iphase_gas;
679 } else {
680 _phase = iphase_liquid;
681 }
682 }
683 _Q = -1;
684 return rho;
685}
686
687CoolPropDbl CoolProp::AbstractCubicBackend::calc_molar_mass() {
688 double summer = 0;
689 for (unsigned int i = 0; i < N; ++i) {
690 summer += mole_fractions[i] * components[i].molemass;
691 }
692 return summer;
693}
694
695void CoolProp::AbstractCubicBackend::set_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter,
696 const double value) {
697 // bound-check indices
698 if (i >= N) {
699 if (j >= N) {
700 throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N - 1));
701 } else {
702 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
703 }
704 } else if (j >= N) {
705 throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N - 1));
706 }
707 if (parameter == "kij" || parameter == "k_ij") {
708 get_cubic()->set_kij(i, j, value);
709 } else {
710 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
711 }
712 for (auto& state : linked_states) {
713 state->set_binary_interaction_double(i, j, parameter, value);
714 }
715};
716double CoolProp::AbstractCubicBackend::get_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter) {
717 // bound-check indices
718 if (i >= N) {
719 if (j >= N) {
720 throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N - 1));
721 } else {
722 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
723 }
724 } else if (j >= N) {
725 throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N - 1));
726 }
727 if (parameter == "kij" || parameter == "k_ij") {
728 return get_cubic()->get_kij(i, j);
729 } else {
730 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
731 }
732};
733
734void CoolProp::AbstractCubicBackend::copy_all_alpha_functions(AbstractCubicBackend* donor) {
735 get_cubic()->set_all_alpha_functions(donor->get_cubic()->get_all_alpha_functions());
736 for (auto& state : linked_states) {
737 auto* ACB = static_cast<AbstractCubicBackend*>(state.get());
738 ACB->copy_all_alpha_functions(this);
739 }
740}
741
742void CoolProp::AbstractCubicBackend::copy_k(AbstractCubicBackend* donor) {
743 get_cubic()->set_kmat(donor->get_cubic()->get_kmat());
744 for (auto& state : linked_states) {
745 auto* ACB = static_cast<AbstractCubicBackend*>(state.get());
746 ACB->copy_k(this);
747 }
748}
749
750void CoolProp::AbstractCubicBackend::copy_internals(AbstractCubicBackend& donor) {
751 this->copy_k(&donor);
752
753 this->components = donor.components;
754 this->set_alpha_from_components();
755 this->set_alpha0_from_components();
756 for (auto& state : linked_states) {
757 auto* ACB = static_cast<AbstractCubicBackend*>(state.get());
758 ACB->components = donor.components;
759 ACB->set_alpha_from_components();
760 ACB->set_alpha0_from_components();
761 }
762}
763
764void CoolProp::AbstractCubicBackend::set_cubic_alpha_C(const size_t i, const std::string& parameter, const double c1, const double c2,
765 const double c3) {
766 // bound-check indices
767 if (i >= N) {
768 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
769 }
770 if (parameter == "MC" || parameter == "mc" || parameter == "Mathias-Copeman") {
771 get_cubic()->set_C_MC(i, c1, c2, c3);
772 if (!components.empty() && i < components.size()) {
773 components[i].alpha_type = "MathiasCopeman";
774 components[i].alpha_coeffs = {c1, c2, c3};
775 }
776 } else if (parameter == "TWU" || parameter == "Twu" || parameter == "twu") {
777 get_cubic()->set_C_Twu(i, c1, c2, c3);
778 if (!components.empty() && i < components.size()) {
779 components[i].alpha_type = "Twu";
780 components[i].alpha_coeffs = {c1, c2, c3};
781 }
782 } else {
783 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
784 }
785 for (auto& state : linked_states) {
786 auto* ACB = static_cast<AbstractCubicBackend*>(state.get());
787 ACB->set_cubic_alpha_C(i, parameter, c1, c2, c3);
788 }
789}
790
791void CoolProp::AbstractCubicBackend::set_fluid_parameter_double(const size_t i, const std::string& parameter, const double value) {
792 // bound-check indices
793 if (i >= N) {
794 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
795 }
796 // Set the volume translation parrameter, currently applied to the whole fluid, not to components.
797 if (parameter == "c" || parameter == "cm" || parameter == "c_m") {
798 get_cubic()->set_cm(value);
799 for (auto& state : linked_states) {
800 auto* ACB = static_cast<AbstractCubicBackend*>(state.get());
801 ACB->set_fluid_parameter_double(i, parameter, value);
802 }
803 } else if (parameter == "Q" || parameter == "Qk" || parameter == "Q_k") {
804 get_cubic()->set_Q_k(i, value);
805 for (auto& state : linked_states) {
806 auto* ACB = static_cast<AbstractCubicBackend*>(state.get());
807 ACB->set_fluid_parameter_double(i, parameter, value);
808 }
809 } else {
810 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
811 }
812}
813double CoolProp::AbstractCubicBackend::get_fluid_parameter_double(const size_t i, const std::string& parameter) {
814 // bound-check indices
815 if (i >= N) {
816 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
817 }
818 // Get the volume translation parrameter, currently applied to the whole fluid, not to components.
819 if (parameter == "c" || parameter == "cm" || parameter == "c_m") {
820 return get_cubic()->get_cm();
821 } else if (parameter == "Q" || parameter == "Qk" || parameter == "Q_k") {
822 return get_cubic()->get_Q_k(i);
823 } else {
824 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
825 }
826}
827
828double CoolProp::AbstractCubicBackend::calc_superanc_Tmax() {
829 if (!is_pure_or_pseudopure) {
830 throw ValueError("calc_superanc_Tmax: fluid must be pure");
831 }
832 int eos_code = get_superanc_eos_code();
833 if (eos_code == CubicSuperAncillary::UNKNOWN_CODE) {
834 throw NotImplementedError("calc_superanc_Tmax: no superancillary available for this cubic EOS");
835 }
836 double Ttilde_max = CubicSuperAncillary::get_Ttilde_max(eos_code);
837 std::vector<double> x = {1.0};
838 // tau = T_r/Tc gives alpha=1 regardless of whether T_r is 1.0 or Tc
839 double tau_at_Tc = cubic->get_Tr() / cubic->get_Tc()[0];
840 double a0 = cubic->am_term(tau_at_Tc, x, 0);
841 double bm = cubic->bm_term(x);
842 return Ttilde_max * a0 / (cubic->get_R_u() * bm);
843}
844
845CoolPropDbl CoolProp::AbstractCubicBackend::calc_saturation_ancillary(parameters param, int Q, parameters given, double value) {
846 if (!is_pure_or_pseudopure) {
847 throw NotImplementedError("calc_saturation_ancillary is not implemented for mixtures in the cubic backend");
848 }
849 int eos_code = get_superanc_eos_code();
850 if (eos_code == CubicSuperAncillary::UNKNOWN_CODE) {
851 throw NotImplementedError("calc_saturation_ancillary: no superancillary available for this cubic EOS");
852 }
853 if (given != iT) {
854 throw NotImplementedError(format("calc_saturation_ancillary: only T-given is supported for cubic EOS; got given=%s",
855 get_parameter_information(given, "short").c_str()));
856 }
857 double T = value;
858 double Tmax = calc_superanc_Tmax();
859 if (T > Tmax) {
860 throw ValueError(format("calc_saturation_ancillary: T (%g K) exceeds superancillary Tmax (%g K)", T, Tmax));
861 }
862 std::vector<double> x = {1.0};
863 double tau = cubic->get_Tr() / T;
864 double am = cubic->am_term(tau, x, 0);
865 double bm = cubic->bm_term(x);
866 double Ttilde = cubic->get_R_u() * T * bm / am;
867
868 using namespace CubicSuperAncillary;
869 if (param == iP) {
870 double p_tilde = supercubic(eos_code, P_CODE, Ttilde);
871 return p_tilde * am / (bm * bm);
872 } else if (param == iDmolar) {
873 if (Q == 0) {
874 return supercubic(eos_code, RHOL_CODE, Ttilde) / bm;
875 } else {
876 return supercubic(eos_code, RHOV_CODE, Ttilde) / bm;
877 }
878 } else {
879 throw NotImplementedError(
880 format("calc_saturation_ancillary: unsupported param=%s for cubic EOS", get_parameter_information(param, "short").c_str()));
881 }
882}
883
884void CoolProp::AbstractCubicBackend::update_QT_pure_superanc(CoolPropDbl Q, CoolPropDbl T) {
885 if (!is_pure_or_pseudopure) {
886 throw ValueError("update_QT_pure_superanc: fluid must be pure");
887 }
888 int eos_code = get_superanc_eos_code();
889 if (eos_code == CubicSuperAncillary::UNKNOWN_CODE) {
890 throw NotImplementedError("update_QT_pure_superanc: no superancillary available for this cubic EOS");
891 }
892 double Tmax = calc_superanc_Tmax();
893 if (T > Tmax) {
894 throw ValueError(format("update_QT_pure_superanc: T (%g K) exceeds superancillary Tmax (%g K)", T, Tmax));
895 }
896
897 std::vector<double> x = {1.0};
898 double tau = cubic->get_Tr() / T;
899 double am = cubic->am_term(tau, x, 0);
900 double bm = cubic->bm_term(x);
901 double Ttilde = cubic->get_R_u() * T * bm / am;
902
903 using namespace CubicSuperAncillary;
904 CoolPropDbl rhoL = supercubic(eos_code, RHOL_CODE, Ttilde) / bm;
905 CoolPropDbl rhoV = supercubic(eos_code, RHOV_CODE, Ttilde) / bm;
906 CoolPropDbl p = supercubic(eos_code, P_CODE, Ttilde) * am / (bm * bm);
907
908 clear();
909 _Q = Q;
910 _T = T;
911 _p = p;
912 _rhomolar = 1.0 / (Q / rhoV + (1.0 - Q) / rhoL);
913 _phase = iphase_twophase;
914
915 SatL->update_TDmolarP_unchecked(T, rhoL, p);
916 SatV->update_TDmolarP_unchecked(T, rhoV, p);
917
918 post_update(false);
919}
920
Generated by doxygen 1.9.4