CoolProp  6.6.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 #include "Solvers.h"
3 #include "Configuration.h"
4 #include "Backends/Helmholtz/VLERoutines.h"
5 
6 void 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 
45 void 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 
69 void 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 
87 std::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 
127 std::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 
135 void 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 
150 void 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 
161 bool 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 
171 void 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 }
192 CoolPropDbl 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 }
198 void 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 
210 CoolPropDbl 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 
283 void 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 
337 void 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 
361 class 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 
411 std::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 
447 void 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 }
510 CoolPropDbl 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 }
514 CoolPropDbl 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 
579 CoolPropDbl 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 
587 void 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 };
608 double 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 
626 void 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 
634 void 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 
642 void 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 
656 void 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 
675 void 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 }
697 double 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 Wed Nov 29 2023 14:26:23 for CoolProp by doxygen 1.9.1