CoolProp 6.8.1dev
An open-source fluid property and humid air property database
CubicBackend.cpp
Go to the documentation of this file.
1#include "CubicBackend.h"
2#include "Solvers.h"
3#include "Configuration.h"
4#include "Backends/Helmholtz/VLERoutines.h"
5
6void CoolProp::AbstractCubicBackend::setup(bool generate_SatL_and_SatV) {
7 N = cubic->get_Tc().size();
8
9 // Set the pure fluid flag
10 is_pure_or_pseudopure = (N == 1);
11
12 // Resize the vectors
13 resize(N);
14
15 // Reset the residual Helmholtz energy class
16 residual_helmholtz.reset(new CubicResidualHelmholtz(this));
17 // If pure, set the mole fractions to be unity
18 if (is_pure_or_pseudopure) {
19 mole_fractions = std::vector<CoolPropDbl>(1, 1.0);
20 } else {
21 mole_fractions.clear();
22 }
23 // Now set the reducing function for the mixture
24 Reducing.reset(new ConstantReducingFunction(cubic->get_Tr(), cubic->get_rhor()));
25
26 // Set the alpha function based on the components in use
27 set_alpha_from_components();
28
29 // Set the ideal-gas helmholtz energy based on the components in use;
30 set_alpha0_from_components();
31
32 // Top-level class can hold copies of the base saturation classes,
33 // saturation classes cannot hold copies of the saturation classes
34 if (generate_SatL_and_SatV) {
35 bool SatLSatV = false;
36 SatL.reset(this->get_copy(SatLSatV));
37 SatL->specify_phase(iphase_liquid);
38 linked_states.push_back(SatL);
39 SatV.reset(this->get_copy(SatLSatV));
40 SatV->specify_phase(iphase_gas);
41 linked_states.push_back(SatV);
42 }
43}
44
45void CoolProp::AbstractCubicBackend::set_alpha_from_components() {
47 if (components.empty()) {
48 return;
49 }
50
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") {
59 acaf.reset(
60 new MathiasCopemanAlphaFunction(get_cubic()->a0_ii(i), c[0], c[1], c[2], get_cubic()->get_Tr() / get_cubic()->get_Tc()[i]));
61 } else {
62 throw ValueError("alpha function is not understood");
63 }
64 cubic->set_alpha_function(i, acaf);
65 }
66 }
67}
68
69void CoolProp::AbstractCubicBackend::set_alpha0_from_components() {
70
71 // If empty so far (e.g., for SatL and SatV)
72 if (components.size() == 0) {
73 return;
74 }
75
76 // Get the vector of CoolProp fluids from the base class
77 std::vector<CoolPropFluid>& _components = HelmholtzEOSMixtureBackend::get_components();
78
79 for (std::size_t i = 0; i < N; ++i) {
80 CoolPropFluid fld;
81 fld.EOSVector.push_back(EquationOfState());
82 fld.EOS().alpha0 = components[i].alpha0;
83 _components.push_back(fld);
84 }
85}
86
87std::string CoolProp::AbstractCubicBackend::fluid_param_string(const std::string& ParamName) {
88 CoolProp::CubicLibrary::CubicsValues cpfluid = components[0];
89 if (!ParamName.compare("name")) {
90 return cpfluid.name;
91 } else if (!ParamName.compare("aliases")) {
92 return strjoin(cpfluid.aliases, get_config_string(LIST_STRING_DELIMITER));
93 } else if (!ParamName.compare("CAS") || !ParamName.compare("CAS_number")) {
94 return cpfluid.CAS;
95 } else if (!ParamName.compare("formula")) {
96 throw NotImplementedError("Parameter \"formula\" not available for cubic backends.");
97 } else if (!ParamName.compare("ASHRAE34")) {
98 throw NotImplementedError("Parameter \"ASHRAE34\" not available for cubic backends.");
99 } else if (!ParamName.compare("REFPROPName") || !ParamName.compare("REFPROP_name") || !ParamName.compare("REFPROPname")) {
100 throw NotImplementedError("Parameter \"REFPROPName\" not available for cubic backends.");
101 } else if (ParamName.find("BibTeX") == 0) // Starts with "BibTeX"
102 {
103 throw NotImplementedError("BibTeX parameters not available for cubic backends.");
104 } else if (ParamName.find("pure") == 0) {
105 if (is_pure()) {
106 return "true";
107 } else {
108 return "false";
109 }
110 } else if (ParamName == "INCHI" || ParamName == "InChI" || ParamName == "INCHI_STRING") {
111 throw NotImplementedError("Parameter \"INCHI\" not available for cubic backends.");
112 } else if (ParamName == "INCHI_Key" || ParamName == "InChIKey" || ParamName == "INCHIKEY") {
113 throw NotImplementedError("Parameter \"INCHI_Key\" not available for cubic backends.");
114 } else if (ParamName == "2DPNG_URL") {
115 throw NotImplementedError("Parameter \"2DPNG_URL\" not available for cubic backends.");
116 } else if (ParamName == "SMILES" || ParamName == "smiles") {
117 throw NotImplementedError("Parameter \"SMILES\" not available for cubic backends.");
118 } else if (ParamName == "CHEMSPIDER_ID") {
119 throw NotImplementedError("Parameter \"CHEMSPIDER_ID\" not available for cubic backends.");
120 } else if (ParamName == "JSON") {
121 return CubicLibrary::get_fluid_as_JSONstring(cpfluid.name);
122 } else {
123 throw ValueError(format("fluid parameter [%s] is invalid", ParamName.c_str()));
124 }
125}
126
127std::vector<std::string> CoolProp::AbstractCubicBackend::calc_fluid_names(void) {
128 std::vector<std::string> out;
129 for (std::size_t i = 0; i < components.size(); ++i) {
130 out.push_back(components[i].name);
131 }
132 return out;
133}
134
135void CoolProp::AbstractCubicBackend::get_linear_reducing_parameters(double& rhomolar_r, double& T_r) {
136 // In the case of models where the reducing temperature is not a function of composition (SRK, PR, etc.),
137 // we need to use an appropriate value for T_r and v_r, here we use a linear weighting
138 T_r = 0;
139 double v_r = 0;
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];
143 // Curve fit from all the pure fluids in CoolProp (thanks to recommendation of A. Kazakov)
144 double v_c_Lmol = 2.14107171795 * (Tc[i] / pc[i] * 1000) + 0.00773144012514; // [L/mol]
145 v_r += mole_fractions[i] * v_c_Lmol / 1000.0;
146 }
147 rhomolar_r = 1 / v_r;
148}
149
150void CoolProp::AbstractCubicBackend::get_critical_point_search_radii(double& R_delta, double& R_tau) {
151 // Get the starting values from the base class
152 CoolProp::HelmholtzEOSMixtureBackend::get_critical_point_search_radii(R_delta, R_tau);
153
154 // Now we scale them to get the appropriate search radii
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;
159}
160
161bool CoolProp::AbstractCubicBackend::get_critical_is_terminated(double& delta, double& tau) {
162 // If the volume is less than the mixture covolume, stop. The mixture covolume is the
163 // smallest volume that is physically allowed for a cubic EOS
164 double b = get_cubic()->bm_term(mole_fractions); // [m^3/mol]
165 double v = 1 / (delta * rhomolar_reducing()); //[m^3/mol]
166 bool covolume_check = v < 1.1 * b;
167
168 return covolume_check;
169}
170
171void CoolProp::AbstractCubicBackend::get_critical_point_starting_values(double& delta0, double& tau0) {
172
173 // Get the starting values from the base class
174 CoolProp::HelmholtzEOSMixtureBackend::get_critical_point_starting_values(delta0, tau0);
175
176 // The starting tau and delta values can be thought of as being given with the
177 // GERG reducing values, or tau0 = Tr_GERG/T = 0.xxx and delta0 = rho/rhor_GERG = 0.xxx
178 // Then we need to multiply by
179 //
180 // Tr_cubic/Tr_GERGlike & rhor_GERGlike/rhor_cubic
181 //
182 // to get shifted values:
183 //
184 // delta0 = rho/rhor_GERG*(rhor_GERGlike/rhor_cubic)
185 // tau0 = Tr_GERG/T*(Tr_cubic/Tr_GERGlike)
186 //
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;
191}
192CoolPropDbl CoolProp::AbstractCubicBackend::calc_pressure_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
193 AbstractCubic* cubic = get_cubic().get();
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));
197}
198void CoolProp::AbstractCubicBackend::update_DmolarT() {
199 // Only works for now when phase is specified
200 if (this->imposed_phase_index != iphase_not_imposed) {
201 _p = calc_pressure_nocache(_T, _rhomolar);
202 _Q = -1;
203 _phase = imposed_phase_index;
204 } else {
205 // Pass call to parent class
206 HelmholtzEOSMixtureBackend::update(DmolarT_INPUTS, this->_rhomolar, this->_T);
207 }
208};
209
210CoolPropDbl CoolProp::AbstractCubicBackend::calc_alphar_deriv_nocache(const int nTau, const int nDelta,
211 const std::vector<CoolPropDbl>& mole_fractions, const CoolPropDbl& tau,
212 const CoolPropDbl& delta) {
213 bool cache_values = true;
214 HelmholtzDerivatives derivs = residual_helmholtz->all(*this, mole_fractions, tau, delta, cache_values);
215 switch (nTau) {
216 case 0: {
217 switch (nDelta) {
218 case 0:
219 return derivs.alphar;
220 case 1:
221 return derivs.dalphar_ddelta;
222 case 2:
223 return derivs.d2alphar_ddelta2;
224 case 3:
225 return derivs.d3alphar_ddelta3;
226 case 4:
227 return derivs.d4alphar_ddelta4;
228 default:
229 throw ValueError(format("nDelta (%d) is invalid", nDelta));
230 }
231 break;
232 }
233 case 1: {
234 switch (nDelta) {
235 case 0:
236 return derivs.dalphar_dtau;
237 case 1:
238 return derivs.d2alphar_ddelta_dtau;
239 case 2:
240 return derivs.d3alphar_ddelta2_dtau;
241 case 3:
242 return derivs.d4alphar_ddelta3_dtau;
243 default:
244 throw ValueError(format("nDelta (%d) is invalid", nDelta));
245 }
246 break;
247 }
248 case 2: {
249 switch (nDelta) {
250 case 0:
251 return derivs.d2alphar_dtau2;
252 case 1:
253 return derivs.d3alphar_ddelta_dtau2;
254 case 2:
255 return derivs.d4alphar_ddelta2_dtau2;
256 default:
257 throw ValueError(format("nDelta (%d) is invalid", nDelta));
258 }
259 }
260 case 3: {
261 switch (nDelta) {
262 case 0:
263 return derivs.d3alphar_dtau3;
264 case 1:
265 return derivs.d4alphar_ddelta_dtau3;
266 default:
267 throw ValueError(format("nDelta (%d) is invalid", nDelta));
268 }
269 }
270 case 4: {
271 switch (nDelta) {
272 case 0:
273 return derivs.d4alphar_dtau4;
274 default:
275 throw ValueError(format("nDelta (%d) is invalid", nDelta));
276 }
277 }
278 default:
279 throw ValueError(format("nTau (%d) is invalid", nTau));
280 }
281}
282
283void CoolProp::AbstractCubicBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) {
284 if (get_debug_level() > 10) {
285 std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
286 get_input_pair_short_desc(input_pair).c_str(), value1, value2)
287 << std::endl;
288 }
289
290 CoolPropDbl ld_value1 = value1, ld_value2 = value2;
291 pre_update(input_pair, ld_value1, ld_value2);
292 value1 = ld_value1;
293 value2 = ld_value2;
294
295 switch (input_pair) {
296 case PT_INPUTS:
297 _p = value1;
298 _T = value2;
299 _rhomolar = solver_rho_Tp(value2 /*T*/, value1 /*p*/);
300 break;
301 case QT_INPUTS:
302 _Q = value1;
303 _T = value2;
304 saturation(input_pair);
305 break;
306 case PQ_INPUTS:
307 _p = value1;
308 _Q = value2;
309 saturation(input_pair);
310 break;
311 case DmolarT_INPUTS:
312 _rhomolar = value1;
313 _T = value2;
314 update_DmolarT();
315 break;
316 case SmolarT_INPUTS:
317 case DmolarP_INPUTS:
318 case DmolarHmolar_INPUTS:
319 case DmolarSmolar_INPUTS:
320 case DmolarUmolar_INPUTS:
321 case HmolarP_INPUTS:
322 case PSmolar_INPUTS:
323 case PUmolar_INPUTS:
324 case HmolarSmolar_INPUTS:
325 case QSmolar_INPUTS:
326 case HmolarQ_INPUTS:
327 case DmolarQ_INPUTS:
328 HelmholtzEOSMixtureBackend::update(input_pair, value1, value2);
329 break;
330 default:
331 throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
332 }
333
334 post_update();
335}
336
337void CoolProp::AbstractCubicBackend::rho_Tp_cubic(CoolPropDbl T, CoolPropDbl p, int& Nsolns, double& rho0, double& rho1, double& rho2) {
338 AbstractCubic* cubic = get_cubic().get();
339 double R = cubic->get_R_u();
340 double am = cubic->am_term(cubic->get_Tr() / T, mole_fractions, 0);
341 double bm = cubic->bm_term(mole_fractions);
342 double cm = cubic->cm_term();
343
344 // Introducing new variables to simplify the equation:
345 double d1 = cm - bm;
346 double d2 = cm + cubic->get_Delta_1() * bm;
347 double d3 = cm + cubic->get_Delta_2() * bm;
348
349 // Cubic coefficients:
350 double crho0 = -p;
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;
354
355 // Solving the cubic:
356 solve_cubic(crho3, crho2, crho1, crho0, Nsolns, rho0, rho1, rho2);
357 sort3(rho0, rho1, rho2);
358 return;
359}
360
361class SaturationResidual : public CoolProp::FuncWrapper1D
362{
363 public:
364 CoolProp::AbstractCubicBackend* ACB;
365 CoolProp::input_pairs inputs;
366 double imposed_variable;
367 double deltaL, deltaV;
368
369 SaturationResidual(){};
370 SaturationResidual(CoolProp::AbstractCubicBackend* ACB, CoolProp::input_pairs inputs, double imposed_variable)
371 : ACB(ACB), inputs(inputs), imposed_variable(imposed_variable){};
372
373 double call(double value) {
374 int Nsolns = 0;
375 double rho0 = -1, rho1 = -1, rho2 = -1, T, p;
376
377 if (inputs == CoolProp::PQ_INPUTS) {
378 T = value;
379 p = imposed_variable;
380 } else if (inputs == CoolProp::QT_INPUTS) {
381 p = value;
382 T = imposed_variable;
383 } else {
384 throw CoolProp::ValueError("Cannot have something other than PQ_INPUTS or QT_INPUTS here");
385 }
386
387 // Calculate the liquid and vapor densities
388 ACB->rho_Tp_cubic(T, p, Nsolns, rho0, rho1, rho2);
389
390 // -----------------------------------------------------
391 // Calculate the difference in Gibbs between the phases
392 // -----------------------------------------------------
393 AbstractCubic* cubic = ACB->get_cubic().get();
394 double rho_r = cubic->get_rhor(), T_r = cubic->get_Tr();
395 double tau = T_r / T;
396 // There are three density solutions, we know the highest is the liquid, the lowest is the vapor
397 deltaL = rho2 / rho_r;
398 deltaV = rho0 / rho_r;
399 // From alpha0; all terms that are only a function of temperature drop out since TL=TV
400 double DELTAgibbs = log(deltaV) - log(deltaL);
401 // From alphar;
402 DELTAgibbs += (cubic->alphar(tau, deltaV, ACB->get_mole_fractions_doubleref(), 0, 0)
403 - cubic->alphar(tau, deltaL, ACB->get_mole_fractions_doubleref(), 0, 0));
404 // From delta*dalphar_dDelta
405 DELTAgibbs += (deltaV * cubic->alphar(tau, deltaV, ACB->get_mole_fractions_doubleref(), 0, 1)
406 - deltaL * cubic->alphar(tau, deltaL, ACB->get_mole_fractions_doubleref(), 0, 1));
407 return DELTAgibbs;
408 };
409};
410
411std::vector<double> CoolProp::AbstractCubicBackend::spinodal_densities() {
413 AbstractCubic* cubic = get_cubic().get();
414 double tau = cubic->get_Tr() / _T;
415 std::vector<double> x(1, 1);
416 double a = cubic->am_term(tau, x, 0);
417 double b = cubic->bm_term(x);
418 double R = cubic->get_R_u();
419 double Delta_1 = cubic->get_Delta_1();
420 double Delta_2 = cubic->get_Delta_2();
421
422 double crho4 = -powInt(Delta_1 * Delta_2, 2) * R * _T * powInt(b, 4) + a * powInt(b, 3) * (Delta_1 + Delta_2);
423 double crho3 =
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;
429 int Nsoln;
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);
434 }
435 if (rho1 > 0 && 1 / rho1 > b) {
436 roots.push_back(rho1);
437 }
438 if (rho2 > 0 && 1 / rho2 > b) {
439 roots.push_back(rho2);
440 }
441 if (rho3 > 0 && 1 / rho3 > b) {
442 roots.push_back(rho3);
443 }
444 return roots;
445}
446
447void CoolProp::AbstractCubicBackend::saturation(CoolProp::input_pairs inputs) {
448 AbstractCubic* cubic = get_cubic().get();
449 double Tc = cubic->get_Tc()[0], pc = cubic->get_pc()[0], acentric = cubic->get_acentric()[0];
450 double rhoL = -1, rhoV = -1;
451 if (inputs == PQ_INPUTS) {
452 if (is_pure_or_pseudopure) {
453 // Estimate temperature from the acentric factor relationship
454 double theta = -log10(_p / pc) * (1 / 0.7 - 1) / (acentric + 1);
455 double Ts_est = Tc / (theta + 1);
456 SaturationResidual resid(this, inputs, _p);
457 static std::string errstr;
458 double Ts = CoolProp::Secant(resid, Ts_est, -0.1, 1e-10, 100);
459 _T = Ts;
460 rhoL = resid.deltaL * cubic->get_Tr();
461 rhoV = resid.deltaV * cubic->get_Tr();
462 this->SatL->update(DmolarT_INPUTS, rhoL, _T);
463 this->SatV->update(DmolarT_INPUTS, rhoV, _T);
464 } else {
465 HelmholtzEOSMixtureBackend::update(PQ_INPUTS, _p, _Q);
466 return;
467 }
468 } else if (inputs == QT_INPUTS) {
469 if (is_pure_or_pseudopure) {
470 SaturationResidual resid(this, inputs, _T);
471 static std::string errstr;
472 // ** Spinodal densities is disabled for now because it is VERY slow :(
473 // std::vector<double> roots = spinodal_densities();
474 std::vector<double> roots;
475
476 // Estimate pressure from the acentric factor relationship
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);
479
480 double ps;
481 if (roots.size() == 2) {
482 double p0 = calc_pressure_nocache(_T, roots[0]);
483 double p1 = calc_pressure_nocache(_T, roots[1]);
484 if (p1 < p0) {
485 std::swap(p0, p1);
486 }
487 //ps = CoolProp::BoundedSecant(resid, p0, p1, pc, -0.01*ps_est, 1e-5, 100);
488 if (p0 > 0 && p1 < pc) {
489 ps = CoolProp::Brent(resid, p0 * 1.0001, p1 * 0.9999, DBL_EPSILON, 1e-10, 100);
490 } else {
491 ps = CoolProp::BoundedSecant(resid, ps_est, 1e-10, pc, -0.01 * ps_est, 1e-5, 100);
492 }
493 } else {
494 ps = CoolProp::BoundedSecant(resid, ps_est, 1e-10, pc, -0.01 * ps_est, 1e-5, 100);
495 }
496
497 _p = ps;
498 rhoL = resid.deltaL * cubic->get_Tr();
499 rhoV = resid.deltaV * cubic->get_Tr();
500 this->SatL->update(DmolarT_INPUTS, rhoL, _T);
501 this->SatV->update(DmolarT_INPUTS, rhoV, _T);
502 } else {
503 HelmholtzEOSMixtureBackend::update(QT_INPUTS, _Q, _T);
504 return;
505 }
506 }
507 _rhomolar = 1 / (_Q / rhoV + (1 - _Q) / rhoL);
508 _phase = iphase_twophase;
509}
510CoolPropDbl CoolProp::AbstractCubicBackend::solver_rho_Tp_global(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomolar_max) {
511 _rhomolar = solver_rho_Tp(T, p, 40000);
512 return static_cast<double>(_rhomolar);
513}
514CoolPropDbl CoolProp::AbstractCubicBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rho_guess) {
515 int Nsoln = 0;
516 double rho0 = 0, rho1 = 0, rho2 = 0, rho = -1;
517 rho_Tp_cubic(T, p, Nsoln, rho0, rho1, rho2); // Densities are sorted in increasing order
518 if (Nsoln == 1) {
519 rho = rho0;
520 } else if (Nsoln == 3) {
521 if (imposed_phase_index != iphase_not_imposed) {
522 // Use imposed phase to select root
523 if (imposed_phase_index == iphase_gas || imposed_phase_index == iphase_supercritical_gas) {
524 if (rho0 > 0) {
525 rho = rho0;
526 } else if (rho1 > 0) {
527 rho = rho1;
528 } else if (rho2 > 0) {
529 rho = rho2;
530 } else {
531 throw CoolProp::ValueError(format("Unable to find gaseous density for T: %g K, p: %g Pa", T, p));
532 }
533 } else if (imposed_phase_index == iphase_liquid || imposed_phase_index == iphase_supercritical_liquid) {
534 rho = rho2;
535 } else {
536 throw ValueError("Specified phase is invalid");
537 }
538 } else {
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);
545 // Gas
546 if (rho0 > 0 && rho0 < rhoV) {
547 rho = rho0;
548 } else if (rho1 > 0 && rho1 < rhoV) {
549 rho = rho1;
550 } else {
551 throw CoolProp::ValueError(format("Unable to find gaseous density for T: %g K, p: %g Pa", T, p));
552 }
553 } else {
554 // Liquid
555 rho = rho2;
556 }
557 } else {
558 throw ValueError("Cubic has three roots, but phase not imposed and guess density not provided");
559 }
560 }
561 } else {
562 throw ValueError("Obtained neither 1 nor three roots");
563 }
564 if (is_pure_or_pseudopure) {
565 // Set some variables at the end
566 this->recalculate_singlephase_phase();
567 } else {
568 if (imposed_phase_index != iphase_not_imposed) {
569 // Use the imposed phase index
570 _phase = imposed_phase_index;
571 } else {
572 _phase = iphase_gas; // TODO: fix this
573 }
574 }
575 _Q = -1;
576 return rho;
577}
578
579CoolPropDbl CoolProp::AbstractCubicBackend::calc_molar_mass(void) {
580 double summer = 0;
581 for (unsigned int i = 0; i < N; ++i) {
582 summer += mole_fractions[i] * components[i].molemass;
583 }
584 return summer;
585}
586
587void CoolProp::AbstractCubicBackend::set_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter,
588 const double value) {
589 // bound-check indices
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));
593 } else {
594 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
595 }
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));
598 }
599 if (parameter == "kij" || parameter == "k_ij") {
600 get_cubic()->set_kij(i, j, value);
601 } else {
602 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
603 }
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);
606 }
607};
608double CoolProp::AbstractCubicBackend::get_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter) {
609 // bound-check indices
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));
613 } else {
614 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N - 1));
615 }
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));
618 }
619 if (parameter == "kij" || parameter == "k_ij") {
620 return get_cubic()->get_kij(i, j);
621 } else {
622 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
623 }
624};
625
626void CoolProp::AbstractCubicBackend::copy_all_alpha_functions(AbstractCubicBackend* donor) {
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) {
629 AbstractCubicBackend* ACB = static_cast<AbstractCubicBackend*>(it->get());
630 ACB->copy_all_alpha_functions(this);
631 }
632}
633
634void CoolProp::AbstractCubicBackend::copy_k(AbstractCubicBackend* donor) {
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) {
637 AbstractCubicBackend* ACB = static_cast<AbstractCubicBackend*>(it->get());
638 ACB->copy_k(this);
639 }
640}
641
642void CoolProp::AbstractCubicBackend::copy_internals(AbstractCubicBackend& donor) {
643 this->copy_k(&donor);
644
645 this->components = donor.components;
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) {
649 AbstractCubicBackend* ACB = static_cast<AbstractCubicBackend*>(it->get());
650 ACB->components = donor.components;
651 ACB->set_alpha_from_components();
652 ACB->set_alpha0_from_components();
653 }
654}
655
656void CoolProp::AbstractCubicBackend::set_cubic_alpha_C(const size_t i, const std::string& parameter, const double c1, const double c2,
657 const double c3) {
658 // bound-check indices
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));
661 }
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);
666 } else {
667 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
668 }
669 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
670 AbstractCubicBackend* ACB = static_cast<AbstractCubicBackend*>(it->get());
671 ACB->set_cubic_alpha_C(i, parameter, c1, c2, c3);
672 }
673}
674
675void CoolProp::AbstractCubicBackend::set_fluid_parameter_double(const size_t i, const std::string& parameter, const double value) {
676 // bound-check indices
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));
679 }
680 // Set the volume translation parrameter, currently applied to the whole fluid, not to components.
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) {
684 AbstractCubicBackend* ACB = static_cast<AbstractCubicBackend*>(it->get());
685 ACB->set_fluid_parameter_double(i, parameter, value);
686 }
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) {
690 AbstractCubicBackend* ACB = static_cast<AbstractCubicBackend*>(it->get());
691 ACB->set_fluid_parameter_double(i, parameter, value);
692 }
693 } else {
694 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
695 }
696}
697double CoolProp::AbstractCubicBackend::get_fluid_parameter_double(const size_t i, const std::string& parameter) {
698 // bound-check indices
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));
701 }
702 // Get the volume translation parrameter, currently applied to the whole fluid, not to components.
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);
707 } else {
708 throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
709 }
710}
Generated on Thu Jun 5 2025 04:10:24 for CoolProp by doxygen 1.9.4