CoolProp  6.6.0
An open-source fluid property and humid air property database
HelmholtzEOSMixtureBackend.cpp
Go to the documentation of this file.
1 /*
2  * AbstractBackend.cpp
3  *
4  * Created on: 20 Dec 2013
5  * Author: jowr
6  */
7 
8 #include <memory>
9 
10 #if defined(_MSC_VER)
11 # define _CRTDBG_MAP_ALLOC
12 # ifndef _CRT_SECURE_NO_WARNINGS
13 # define _CRT_SECURE_NO_WARNINGS
14 # endif
15 # include <crtdbg.h>
16 # include <sys/stat.h>
17 #else
18 # include <sys/stat.h>
19 #endif
20 
21 #include <string>
22 //#include "CoolProp.h"
23 
24 #include "HelmholtzEOSMixtureBackend.h"
25 #include "HelmholtzEOSBackend.h"
26 #include "Fluids/FluidLibrary.h"
27 #include "Solvers.h"
28 #include "MatrixMath.h"
29 #include "VLERoutines.h"
30 #include "FlashRoutines.h"
31 #include "TransportRoutines.h"
32 #include "MixtureDerivatives.h"
33 #include "PhaseEnvelopeRoutines.h"
34 #include "ReducingFunctions.h"
35 #include "MixtureParameters.h"
36 #include "IdealCurves.h"
37 #include "MixtureParameters.h"
38 #include <stdlib.h>
39 
40 static int deriv_counter = 0;
41 
42 namespace CoolProp {
43 
44 class HEOSGenerator : public AbstractStateGenerator
45 {
46  public:
47  AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
48  if (fluid_names.size() == 1) {
49  return new HelmholtzEOSBackend(fluid_names[0]);
50  } else {
51  return new HelmholtzEOSMixtureBackend(fluid_names);
52  }
53  };
54 };
55 // This static initialization will cause the generator to register
56 static CoolProp::GeneratorInitializer<HEOSGenerator> heos_gen(CoolProp::HEOS_BACKEND_FAMILY);
57 
58 HelmholtzEOSMixtureBackend::HelmholtzEOSMixtureBackend() {
59  imposed_phase_index = iphase_not_imposed;
60  is_pure_or_pseudopure = false;
61  N = 0;
62  _phase = iphase_unknown;
63  // Reset the residual Helmholtz energy class
64  residual_helmholtz.reset(new ResidualHelmholtz());
65 }
66 HelmholtzEOSMixtureBackend::HelmholtzEOSMixtureBackend(const std::vector<std::string>& component_names, bool generate_SatL_and_SatV) {
67  std::vector<CoolPropFluid> components(component_names.size());
68  for (unsigned int i = 0; i < components.size(); ++i) {
69  components[i] = get_library().get(component_names[i]);
70  }
71 
72  // Reset the residual Helmholtz energy class
73  residual_helmholtz.reset(new ResidualHelmholtz());
74 
75  // Set the components and associated flags
76  set_components(components, generate_SatL_and_SatV);
77 
78  // Set the phase to default unknown value
79  _phase = iphase_unknown;
80 }
81 HelmholtzEOSMixtureBackend::HelmholtzEOSMixtureBackend(const std::vector<CoolPropFluid>& components, bool generate_SatL_and_SatV) {
82 
83  // Reset the residual Helmholtz energy class
84  residual_helmholtz.reset(new ResidualHelmholtz());
85 
86  // Set the components and associated flags
87  set_components(components, generate_SatL_and_SatV);
88 
89  // Set the phase to default unknown value
90  _phase = iphase_unknown;
91 }
92 void HelmholtzEOSMixtureBackend::set_components(const std::vector<CoolPropFluid>& components, bool generate_SatL_and_SatV) {
93 
94  // Copy the components
95  this->components = components;
96  this->N = components.size();
97 
98  is_pure_or_pseudopure = (components.size() == 1);
99  if (is_pure_or_pseudopure) {
100  mole_fractions = std::vector<CoolPropDbl>(1, 1);
101  std::vector<std::vector<double>> ones(1, std::vector<double>(1, 1));
102  Reducing = shared_ptr<ReducingFunction>(new GERG2008ReducingFunction(components, ones, ones, ones, ones));
103  } else {
104  // Set the mixture parameters - binary pair reducing functions, departure functions, F_ij, etc.
105  set_mixture_parameters();
106  }
107 
108  imposed_phase_index = iphase_not_imposed;
109 
110  // Top-level class can hold copies of the base saturation classes,
111  // saturation classes cannot hold copies of the saturation classes
112  if (generate_SatL_and_SatV) {
113  SatL.reset(get_copy(false));
114  SatL->specify_phase(iphase_liquid);
115  linked_states.push_back(SatL);
116  SatV.reset(get_copy(false));
117  SatV->specify_phase(iphase_gas);
118  linked_states.push_back(SatV);
119  }
120 }
121 void HelmholtzEOSMixtureBackend::set_mole_fractions(const std::vector<CoolPropDbl>& mf) {
122  if (mf.size() != N) {
123  throw ValueError(format("size of mole fraction vector [%d] does not equal that of component vector [%d]", mf.size(), N));
124  }
125  // Copy values without reallocating memory
126  this->mole_fractions = mf; // Most effective copy
127  this->resize(N); // No reallocation of this->mole_fractions happens
128  clear_comp_change();
129 };
130 void HelmholtzEOSMixtureBackend::sync_linked_states(const HelmholtzEOSMixtureBackend* const source) {
131  residual_helmholtz.reset(source->residual_helmholtz->copy_ptr());
132  if (source->Reducing) {
133  Reducing.reset(source->Reducing->copy());
134  }
135  // Recurse into linked states of the class
136  for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
137  it->get()->sync_linked_states(source);
138  }
139 }
140 HelmholtzEOSMixtureBackend* HelmholtzEOSMixtureBackend::get_copy(bool generate_SatL_and_SatV) {
141  // Set up the class with these components
142  HelmholtzEOSMixtureBackend* ptr = new HelmholtzEOSMixtureBackend(components, generate_SatL_and_SatV);
143  // Recursively walk into linked states, setting the departure and reducing terms
144  // to be equal to the parent (this instance)
145  ptr->sync_linked_states(this);
146  return ptr;
147 };
148 void HelmholtzEOSMixtureBackend::set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions) {
149  if (mass_fractions.size() != N) {
150  throw ValueError(format("size of mass fraction vector [%d] does not equal that of component vector [%d]", mass_fractions.size(), N));
151  }
152  std::vector<CoolPropDbl> moles;
153  CoolPropDbl sum_moles = 0.0;
154  CoolPropDbl tmp = 0.0;
155  for (unsigned int i = 0; i < components.size(); ++i) {
156  tmp = mass_fractions[i] / components[i].molar_mass();
157  moles.push_back(tmp);
158  sum_moles += tmp;
159  }
160  std::vector<CoolPropDbl> mole_fractions;
161  for (std::vector<CoolPropDbl>::iterator it = moles.begin(); it != moles.end(); ++it) {
162  mole_fractions.push_back(*it / sum_moles);
163  }
164  this->set_mole_fractions(mole_fractions);
165 };
166 void HelmholtzEOSMixtureBackend::resize(std::size_t N) {
167  this->mole_fractions.resize(N);
168  this->K.resize(N);
169  this->lnK.resize(N);
170  for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
171  it->get()->N = N;
172  it->get()->resize(N);
173  }
174 }
175 void HelmholtzEOSMixtureBackend::recalculate_singlephase_phase() {
176  if (p() > p_critical()) {
177  if (T() > T_critical()) {
178  _phase = iphase_supercritical;
179  } else {
180  _phase = iphase_supercritical_liquid;
181  }
182  } else {
183  if (T() > T_critical()) {
184  _phase = iphase_supercritical_gas;
185  } else {
186  // Liquid or vapor
187  if (rhomolar() > rhomolar_critical()) {
188  _phase = iphase_liquid;
189  } else {
190  _phase = iphase_gas;
191  }
192  }
193  }
194 }
195 std::string HelmholtzEOSMixtureBackend::fluid_param_string(const std::string& ParamName) {
196  CoolProp::CoolPropFluid cpfluid = get_components()[0];
197  if (!ParamName.compare("name")) {
198  return cpfluid.name;
199  } else if (!ParamName.compare("aliases")) {
200  return strjoin(cpfluid.aliases, get_config_string(LIST_STRING_DELIMITER));
201  } else if (!ParamName.compare("CAS") || !ParamName.compare("CAS_number")) {
202  return cpfluid.CAS;
203  } else if (!ParamName.compare("formula")) {
204  return cpfluid.formula;
205  } else if (!ParamName.compare("ASHRAE34")) {
206  return cpfluid.environment.ASHRAE34;
207  } else if (!ParamName.compare("REFPROPName") || !ParamName.compare("REFPROP_name") || !ParamName.compare("REFPROPname")) {
208  return cpfluid.REFPROPname;
209  } else if (ParamName.find("BibTeX") == 0) // Starts with "BibTeX"
210  {
211  std::vector<std::string> parts = strsplit(ParamName, '-');
212  if (parts.size() != 2) {
213  throw ValueError(format("Unable to parse BibTeX string %s", ParamName.c_str()));
214  }
215  std::string key = parts[1];
216  if (!key.compare("EOS")) {
217  return cpfluid.EOS().BibTeX_EOS;
218  } else if (!key.compare("CP0")) {
219  return cpfluid.EOS().BibTeX_CP0;
220  } else if (!key.compare("VISCOSITY")) {
221  return cpfluid.transport.BibTeX_viscosity;
222  } else if (!key.compare("CONDUCTIVITY")) {
223  return cpfluid.transport.BibTeX_conductivity;
224  } else if (!key.compare("ECS_LENNARD_JONES")) {
225  throw NotImplementedError();
226  } else if (!key.compare("ECS_VISCOSITY_FITS")) {
227  throw NotImplementedError();
228  } else if (!key.compare("ECS_CONDUCTIVITY_FITS")) {
229  throw NotImplementedError();
230  } else if (!key.compare("SURFACE_TENSION")) {
231  return cpfluid.ancillaries.surface_tension.BibTeX;
232  } else if (!key.compare("MELTING_LINE")) {
233  return cpfluid.ancillaries.melting_line.BibTeX;
234  } else {
235  throw CoolProp::KeyError(format("Bad key to get_BibTeXKey [%s]", key.c_str()));
236  }
237  } else if (ParamName.find("pure") == 0) {
238  if (is_pure()) {
239  return "true";
240  } else {
241  return "false";
242  }
243  } else if (ParamName == "INCHI" || ParamName == "InChI" || ParamName == "INCHI_STRING") {
244  return cpfluid.InChI;
245  } else if (ParamName == "INCHI_Key" || ParamName == "InChIKey" || ParamName == "INCHIKEY") {
246  return cpfluid.InChIKey;
247  } else if (ParamName == "2DPNG_URL") {
248  return cpfluid.TwoDPNG_URL;
249  } else if (ParamName == "SMILES" || ParamName == "smiles") {
250  return cpfluid.smiles;
251  } else if (ParamName == "CHEMSPIDER_ID") {
252  return format("%d", cpfluid.ChemSpider_id);
253  } else if (ParamName == "JSON") {
254  return get_fluid_as_JSONstring(cpfluid.CAS);
255  } else {
256  throw ValueError(format("fluid parameter [%s] is invalid", ParamName.c_str()));
257  }
258 }
259 
260 void HelmholtzEOSMixtureBackend::apply_simple_mixing_rule(std::size_t i, std::size_t j, const std::string& model) {
261  // bound-check indices
262  if (i < 0 || i >= N) {
263  if (j < 0 || j >= N) {
264  throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N-1));
265  } else {
266  throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N-1));
267  }
268  } else if (j < 0 || j >= N) {
269  throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N-1));
270  }
271  if (model == "linear") {
272  double Tc1 = get_fluid_constant(i, iT_critical), Tc2 = get_fluid_constant(j, iT_critical);
273  double gammaT = 0.5 * (Tc1 + Tc2) / sqrt(Tc1 * Tc2);
274  double rhoc1 = get_fluid_constant(i, irhomolar_critical), rhoc2 = get_fluid_constant(j, irhomolar_critical);
275  double gammaV = 4.0 * (1 / rhoc1 + 1 / rhoc2) / pow(pow(rhoc1, -1.0 / 3.0) + pow(rhoc2, -1.0 / 3.0), 3);
276  set_binary_interaction_double(i, j, "betaT", 1.0);
277  set_binary_interaction_double(i, j, "gammaT", gammaT);
278  set_binary_interaction_double(i, j, "betaV", 1.0);
279  set_binary_interaction_double(i, j, "gammaV", gammaV);
280  } else if (model == "Lorentz-Berthelot") {
281  set_binary_interaction_double(i, j, "betaT", 1.0);
282  set_binary_interaction_double(i, j, "gammaT", 1.0);
283  set_binary_interaction_double(i, j, "betaV", 1.0);
284  set_binary_interaction_double(i, j, "gammaV", 1.0);
285  } else {
286  throw ValueError(format("mixing rule [%s] is not understood", model.c_str()));
287  }
288 }
290 void HelmholtzEOSMixtureBackend::set_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter,
291  const double value) {
292  // bound-check indices
293  if (i < 0 || i >= N) {
294  if (j < 0 || j >= N) {
295  throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N-1));
296  } else {
297  throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N-1));
298  }
299  } else if (j < 0 || j >= N) {
300  throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N-1));
301  }
302  if (parameter == "Fij") {
303  residual_helmholtz->Excess.F[i][j] = value;
304  residual_helmholtz->Excess.F[j][i] = value;
305  } else {
306  Reducing->set_binary_interaction_double(i, j, parameter, value);
307  }
309  for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
310  it->get()->set_binary_interaction_double(i, j, parameter, value);
311  }
312 };
314 double HelmholtzEOSMixtureBackend::get_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter) {
315  // bound-check indices
316  if (i < 0 || i >= N) {
317  if (j < 0 || j >= N) {
318  throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N-1));
319  } else {
320  throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N-1));
321  }
322  } else if (j < 0 || j >= N) {
323  throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N-1));
324  }
325  if (parameter == "Fij") {
326  return residual_helmholtz->Excess.F[i][j];
327  } else {
328  return Reducing->get_binary_interaction_double(i, j, parameter);
329  }
330 };
332 //std::string HelmholtzEOSMixtureBackend::get_binary_interaction_string(const std::string &CAS1, const std::string &CAS2, const std::string &parameter){
333 // return CoolProp::get_mixture_binary_pair_data(CAS1, CAS2, parameter);
334 //}
336 void HelmholtzEOSMixtureBackend::set_binary_interaction_string(const std::size_t i, const std::size_t j, const std::string& parameter,
337  const std::string& value) {
338  // bound-check indices
339  if (i < 0 || i >= N) {
340  if (j < 0 || j >= N) {
341  throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N-1));
342  } else {
343  throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N-1));
344  }
345  } else if (j < 0 || j >= N) {
346  throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N-1));
347  }
348  if (parameter == "function") {
349  residual_helmholtz->Excess.DepartureFunctionMatrix[i][j].reset(get_departure_function(value));
350  residual_helmholtz->Excess.DepartureFunctionMatrix[j][i].reset(get_departure_function(value));
351  } else {
352  throw ValueError(format("Cannot process this string parameter [%s] in set_binary_interaction_string", parameter.c_str()));
353  }
355  for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
356  it->get()->set_binary_interaction_string(i, j, parameter, value);
357  }
358 };
359 
360 void HelmholtzEOSMixtureBackend::calc_change_EOS(const std::size_t i, const std::string& EOS_name) {
361 
362  if (i < components.size()) {
363  CoolPropFluid& fluid = components[i];
364  EquationOfState& EOS = fluid.EOSVector[0];
365 
366  if (EOS_name == "SRK" || EOS_name == "Peng-Robinson") {
367 
368  // Get the parameters for the cubic EOS
369  CoolPropDbl Tc = EOS.reduce.T;
370  CoolPropDbl pc = EOS.reduce.p;
371  CoolPropDbl rhomolarc = EOS.reduce.rhomolar;
372  CoolPropDbl acentric = EOS.acentric;
373  CoolPropDbl R = 8.3144598;
374 
375  // Remove the residual part
376  EOS.alphar.empty_the_EOS();
377  // Set the contribution
378  shared_ptr<AbstractCubic> ac;
379  if (EOS_name == "SRK") {
380  ac.reset(new SRK(Tc, pc, acentric, R));
381  } else {
382  ac.reset(new PengRobinson(Tc, pc, acentric, R));
383  }
384  ac->set_Tr(Tc);
385  ac->set_rhor(rhomolarc);
386  EOS.alphar.cubic = ResidualHelmholtzGeneralizedCubic(ac);
387  } else if (EOS_name == "XiangDeiters") {
388 
389  // Get the parameters for the EOS
390  CoolPropDbl Tc = EOS.reduce.T;
391  CoolPropDbl pc = EOS.reduce.p;
392  CoolPropDbl rhomolarc = EOS.reduce.rhomolar;
393  CoolPropDbl acentric = EOS.acentric;
394  CoolPropDbl R = 8.3144598;
395 
396  // Remove the residual part
397  EOS.alphar.empty_the_EOS();
398  // Set the Xiang & Deiters contribution
399  EOS.alphar.XiangDeiters = ResidualHelmholtzXiangDeiters(Tc, pc, rhomolarc, acentric, R);
400  }
401  } else {
402  throw ValueError(format("Index [%d] is invalid", i));
403  }
404  // Now do the same thing to the saturated liquid and vapor instances if possible
405  if (this->SatL) SatL->change_EOS(i, EOS_name);
406  if (this->SatV) SatV->change_EOS(i, EOS_name);
407 }
408 void HelmholtzEOSMixtureBackend::calc_phase_envelope(const std::string& type) {
409  // Clear the phase envelope data
410  PhaseEnvelope = PhaseEnvelopeData();
411  // Build the phase envelope
412  PhaseEnvelopeRoutines::build(*this, type);
413  // Finalize the phase envelope
414  PhaseEnvelopeRoutines::finalize(*this);
415 };
416 void HelmholtzEOSMixtureBackend::set_mixture_parameters() {
417  // Build the matrix of binary-pair reducing functions
418  MixtureParameters::set_mixture_parameters(*this);
419 }
420 void HelmholtzEOSMixtureBackend::update_states(void) {
421  CoolPropFluid& component = components[0];
422  EquationOfState& EOS = component.EOSVector[0];
423 
424  // Clear the state class
425  clear();
426 
427  // Calculate the new enthalpy and entropy values
428  update(DmolarT_INPUTS, EOS.hs_anchor.rhomolar, EOS.hs_anchor.T);
429  EOS.hs_anchor.hmolar = hmolar();
430  EOS.hs_anchor.smolar = smolar();
431 
432  // Calculate the new enthalpy and entropy values at the reducing state
433  update(DmolarT_INPUTS, EOS.reduce.rhomolar, EOS.reduce.T);
434  EOS.reduce.hmolar = hmolar();
435  EOS.reduce.smolar = smolar();
436 
437  // Clear again just to be sure
438  clear();
439 }
440 const CoolProp::SimpleState& HelmholtzEOSMixtureBackend::calc_state(const std::string& state) {
441  if (is_pure_or_pseudopure) {
442  if (!state.compare("hs_anchor")) {
443  return components[0].EOS().hs_anchor;
444  } else if (!state.compare("max_sat_T")) {
445  return components[0].EOS().max_sat_T;
446  } else if (!state.compare("max_sat_p")) {
447  return components[0].EOS().max_sat_p;
448  } else if (!state.compare("reducing")) {
449  return components[0].EOS().reduce;
450  } else if (!state.compare("critical")) {
451  return components[0].crit;
452  } else if (!state.compare("triple_liquid")) {
453  return components[0].triple_liquid;
454  } else if (!state.compare("triple_vapor")) {
455  return components[0].triple_vapor;
456  } else {
457  throw ValueError(format("This state [%s] is invalid to calc_state", state.c_str()));
458  }
459  } else {
460  if (!state.compare("critical")) {
461  return _critical;
462  } else {
463  throw ValueError(format("calc_state not supported for mixtures"));
464  }
465  }
466 };
467 CoolPropDbl HelmholtzEOSMixtureBackend::calc_acentric_factor(void) {
468  if (is_pure_or_pseudopure) {
469  return components[0].EOS().acentric;
470  } else {
471  throw ValueError("acentric factor cannot be calculated for mixtures");
472  }
473 }
474 CoolPropDbl HelmholtzEOSMixtureBackend::calc_gas_constant(void) {
475  if (is_pure_or_pseudopure) {
476  return components[0].gas_constant();
477  } else {
478  if (get_config_bool(NORMALIZE_GAS_CONSTANTS)) {
479  return get_config_double(R_U_CODATA);
480  } else {
481  // mass fraction weighted average of the components
482  double summer = 0;
483  for (unsigned int i = 0; i < components.size(); ++i) {
484  summer += mole_fractions[i] * components[i].gas_constant();
485  }
486  return summer;
487  }
488  }
489 }
490 CoolPropDbl HelmholtzEOSMixtureBackend::calc_molar_mass(void) {
491  double summer = 0;
492  for (unsigned int i = 0; i < components.size(); ++i) {
493  summer += mole_fractions[i] * components[i].molar_mass();
494  }
495  return summer;
496 }
497 CoolPropDbl HelmholtzEOSMixtureBackend::calc_saturation_ancillary(parameters param, int Q, parameters given, double value) {
498  if (is_pure_or_pseudopure) {
499  if (param == iP && given == iT) {
500  // p = f(T), direct evaluation
501  switch (Q) {
502  case 0:
503  return components[0].ancillaries.pL.evaluate(value);
504  case 1:
505  return components[0].ancillaries.pV.evaluate(value);
506  }
507  } else if (param == iT && given == iP) {
508  // T = f(p), inverse evaluation
509  switch (Q) {
510  case 0:
511  return components[0].ancillaries.pL.invert(value);
512  case 1:
513  return components[0].ancillaries.pV.invert(value);
514  }
515  } else if (param == iDmolar && given == iT) {
516  // rho = f(T), inverse evaluation
517  switch (Q) {
518  case 0:
519  return components[0].ancillaries.rhoL.evaluate(value);
520  case 1:
521  return components[0].ancillaries.rhoV.evaluate(value);
522  }
523  } else if (param == iT && given == iDmolar) {
524  // T = f(rho), inverse evaluation
525  switch (Q) {
526  case 0:
527  return components[0].ancillaries.rhoL.invert(value);
528  case 1:
529  return components[0].ancillaries.rhoV.invert(value);
530  }
531  } else if (param == isurface_tension && given == iT) {
532  return components[0].ancillaries.surface_tension.evaluate(value);
533  } else {
534  throw ValueError(format("calc of %s given %s is invalid in calc_saturation_ancillary", get_parameter_information(param, "short").c_str(),
535  get_parameter_information(given, "short").c_str()));
536  }
537 
538  throw ValueError(format("Q [%d] is invalid in calc_saturation_ancillary", Q));
539  } else {
540  throw NotImplementedError(format("calc_saturation_ancillary not implemented for mixtures"));
541  }
542 }
543 
544 CoolPropDbl HelmholtzEOSMixtureBackend::calc_melting_line(int param, int given, CoolPropDbl value) {
545  if (is_pure_or_pseudopure) {
546  return components[0].ancillaries.melting_line.evaluate(param, given, value);
547  } else {
548  throw NotImplementedError(format("calc_melting_line not implemented for mixtures"));
549  }
550 }
551 CoolPropDbl HelmholtzEOSMixtureBackend::calc_surface_tension(void) {
552  if (is_pure_or_pseudopure) {
553  if ((_phase == iphase_twophase) || (_phase == iphase_critical_point)) { // if within the two phase region or at critical point
554  return components[0].ancillaries.surface_tension.evaluate(T()); // calculate surface tension and return
555  } else { // else state point not in the two phase region
556  throw ValueError(format("surface tension is only defined within the two-phase region; Try PQ or QT inputs")); // throw error
557  }
558  } else {
559  throw NotImplementedError(format("surface tension not implemented for mixtures"));
560  }
561 }
562 CoolPropDbl HelmholtzEOSMixtureBackend::calc_viscosity_dilute(void) {
563  if (is_pure_or_pseudopure) {
564  CoolPropDbl eta_dilute;
565  switch (components[0].transport.viscosity_dilute.type) {
566  case ViscosityDiluteVariables::VISCOSITY_DILUTE_KINETIC_THEORY:
567  eta_dilute = TransportRoutines::viscosity_dilute_kinetic_theory(*this);
568  break;
569  case ViscosityDiluteVariables::VISCOSITY_DILUTE_COLLISION_INTEGRAL:
570  eta_dilute = TransportRoutines::viscosity_dilute_collision_integral(*this);
571  break;
572  case ViscosityDiluteVariables::VISCOSITY_DILUTE_POWERS_OF_T:
573  eta_dilute = TransportRoutines::viscosity_dilute_powers_of_T(*this);
574  break;
575  case ViscosityDiluteVariables::VISCOSITY_DILUTE_POWERS_OF_TR:
576  eta_dilute = TransportRoutines::viscosity_dilute_powers_of_Tr(*this);
577  break;
578  case ViscosityDiluteVariables::VISCOSITY_DILUTE_COLLISION_INTEGRAL_POWERS_OF_TSTAR:
579  eta_dilute = TransportRoutines::viscosity_dilute_collision_integral_powers_of_T(*this);
580  break;
581  case ViscosityDiluteVariables::VISCOSITY_DILUTE_ETHANE:
582  eta_dilute = TransportRoutines::viscosity_dilute_ethane(*this);
583  break;
584  case ViscosityDiluteVariables::VISCOSITY_DILUTE_CYCLOHEXANE:
585  eta_dilute = TransportRoutines::viscosity_dilute_cyclohexane(*this);
586  break;
587  case ViscosityDiluteVariables::VISCOSITY_DILUTE_CO2_LAESECKE_JPCRD_2017:
588  eta_dilute = TransportRoutines::viscosity_dilute_CO2_LaeseckeJPCRD2017(*this);
589  break;
590  default:
591  throw ValueError(
592  format("dilute viscosity type [%d] is invalid for fluid %s", components[0].transport.viscosity_dilute.type, name().c_str()));
593  }
594  return eta_dilute;
595  } else {
596  throw NotImplementedError(format("dilute viscosity not implemented for mixtures"));
597  }
598 }
599 CoolPropDbl HelmholtzEOSMixtureBackend::calc_viscosity_background() {
600  CoolPropDbl eta_dilute = calc_viscosity_dilute(), initial_density = 0, residual = 0;
601  return calc_viscosity_background(eta_dilute, initial_density, residual);
602 }
603 CoolPropDbl HelmholtzEOSMixtureBackend::calc_viscosity_background(CoolPropDbl eta_dilute, CoolPropDbl& initial_density, CoolPropDbl& residual) {
604 
605  switch (components[0].transport.viscosity_initial.type) {
606  case ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_RAINWATER_FRIEND: {
607  CoolPropDbl B_eta_initial = TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(*this);
608  CoolPropDbl rho = rhomolar();
609  initial_density = eta_dilute * B_eta_initial * rho; //TODO: Check units once AMTG
610  break;
611  }
612  case ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_EMPIRICAL: {
613  initial_density = TransportRoutines::viscosity_initial_density_dependence_empirical(*this);
614  break;
615  }
616  case ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_NOT_SET: {
617  break;
618  }
619  }
620 
621  // Higher order terms
622  switch (components[0].transport.viscosity_higher_order.type) {
623  case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_BATSCHINKI_HILDEBRAND:
624  residual = TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand(*this);
625  break;
626  case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_FRICTION_THEORY:
627  residual = TransportRoutines::viscosity_higher_order_friction_theory(*this);
628  break;
629  case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HYDROGEN:
630  residual = TransportRoutines::viscosity_hydrogen_higher_order_hardcoded(*this);
631  break;
632  case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_TOLUENE:
633  residual = TransportRoutines::viscosity_toluene_higher_order_hardcoded(*this);
634  break;
635  case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HEXANE:
636  residual = TransportRoutines::viscosity_hexane_higher_order_hardcoded(*this);
637  break;
638  case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HEPTANE:
639  residual = TransportRoutines::viscosity_heptane_higher_order_hardcoded(*this);
640  break;
641  case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_ETHANE:
642  residual = TransportRoutines::viscosity_ethane_higher_order_hardcoded(*this);
643  break;
644  case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_BENZENE:
645  residual = TransportRoutines::viscosity_benzene_higher_order_hardcoded(*this);
646  break;
647  case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_CO2_LAESECKE_JPCRD_2017:
648  residual = TransportRoutines::viscosity_CO2_higher_order_hardcoded_LaeseckeJPCRD2017(*this);
649  break;
650  default:
651  throw ValueError(
652  format("higher order viscosity type [%d] is invalid for fluid %s", components[0].transport.viscosity_dilute.type, name().c_str()));
653  }
654 
655  return initial_density + residual;
656 }
657 
658 CoolPropDbl HelmholtzEOSMixtureBackend::calc_viscosity(void) {
659  if (is_pure_or_pseudopure) {
660  CoolPropDbl dilute = 0, initial_density = 0, residual = 0, critical = 0;
661  calc_viscosity_contributions(dilute, initial_density, residual, critical);
662  return dilute + initial_density + residual + critical;
663  } else {
664  set_warning_string("Mixture model for viscosity is highly approximate");
665  CoolPropDbl summer = 0;
666  for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
667  shared_ptr<HelmholtzEOSBackend> HEOS(new HelmholtzEOSBackend(components[i]));
668  HEOS->update(DmolarT_INPUTS, _rhomolar, _T);
669  summer += mole_fractions[i] * log(HEOS->viscosity());
670  }
671  return exp(summer);
672  }
673 }
674 void HelmholtzEOSMixtureBackend::calc_viscosity_contributions(CoolPropDbl& dilute, CoolPropDbl& initial_density, CoolPropDbl& residual,
675  CoolPropDbl& critical) {
676  if (is_pure_or_pseudopure) {
677  // Reset the variables
678  dilute = 0;
679  initial_density = 0;
680  residual = 0;
681  critical = 0;
682 
683  // Get a reference for code cleanness
684  CoolPropFluid& component = components[0];
685 
686  if (!component.transport.viscosity_model_provided) {
687  throw ValueError(format("Viscosity model is not available for this fluid"));
688  }
689 
690  // Check if using ECS
691  if (component.transport.viscosity_using_ECS) {
692  // Get reference fluid name
693  std::string fluid_name = component.transport.viscosity_ecs.reference_fluid;
694  std::vector<std::string> names(1, fluid_name);
695  // Get a managed pointer to the reference fluid for ECS
696  shared_ptr<HelmholtzEOSMixtureBackend> ref_fluid(new HelmholtzEOSMixtureBackend(names));
697  // Get the viscosity using ECS and stick in the critical value
698  critical = TransportRoutines::viscosity_ECS(*this, *ref_fluid);
699  return;
700  }
701 
702  // Check if using Chung model
703  if (component.transport.viscosity_using_Chung) {
704  // Get the viscosity using ECS and stick in the critical value
705  critical = TransportRoutines::viscosity_Chung(*this);
706  return;
707  }
708 
709  // Check if using rho*sr model
710  if (component.transport.viscosity_using_rhosr) {
711  // Get the viscosity using rho*sr model and stick in the critical value
712  critical = TransportRoutines::viscosity_rhosr(*this);
713  return;
714  }
715 
716  if (component.transport.hardcoded_viscosity != CoolProp::TransportPropertyData::VISCOSITY_NOT_HARDCODED) {
717  // Evaluate hardcoded model and stick in the critical value
718  switch (component.transport.hardcoded_viscosity) {
719  case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_WATER:
720  critical = TransportRoutines::viscosity_water_hardcoded(*this);
721  break;
722  case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_HEAVYWATER:
723  critical = TransportRoutines::viscosity_heavywater_hardcoded(*this);
724  break;
725  case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_HELIUM:
726  critical = TransportRoutines::viscosity_helium_hardcoded(*this);
727  break;
728  case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_R23:
729  critical = TransportRoutines::viscosity_R23_hardcoded(*this);
730  break;
731  case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_METHANOL:
732  critical = TransportRoutines::viscosity_methanol_hardcoded(*this);
733  break;
734  case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_M_XYLENE:
735  critical = TransportRoutines::viscosity_m_xylene_hardcoded(*this);
736  break;
737  case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_O_XYLENE:
738  critical = TransportRoutines::viscosity_o_xylene_hardcoded(*this);
739  break;
740  case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_P_XYLENE:
741  critical = TransportRoutines::viscosity_p_xylene_hardcoded(*this);
742  break;
743  default:
744  throw ValueError(
745  format("hardcoded viscosity type [%d] is invalid for fluid %s", component.transport.hardcoded_viscosity, name().c_str()));
746  }
747  return;
748  }
749 
750  // -------------------------
751  // Normal evaluation
752  // -------------------------
753 
754  // Dilute part
755  dilute = calc_viscosity_dilute();
756 
757  // Background viscosity given by the sum of the initial density dependence and higher order terms
758  calc_viscosity_background(dilute, initial_density, residual);
759 
760  // Critical part (no fluids have critical enhancement for viscosity currently)
761  critical = 0;
762  } else {
763  throw ValueError("calc_viscosity_contributions invalid for mixtures");
764  }
765 }
766 void HelmholtzEOSMixtureBackend::calc_conductivity_contributions(CoolPropDbl& dilute, CoolPropDbl& initial_density, CoolPropDbl& residual,
767  CoolPropDbl& critical) {
768  if (is_pure_or_pseudopure) {
769  // Reset the variables
770  dilute = 0;
771  initial_density = 0;
772  residual = 0;
773  critical = 0;
774 
775  // Get a reference for code cleanness
776  CoolPropFluid& component = components[0];
777 
778  if (!component.transport.conductivity_model_provided) {
779  throw ValueError(format("Thermal conductivity model is not available for this fluid"));
780  }
781 
782  // Check if using ECS
783  if (component.transport.conductivity_using_ECS) {
784  // Get reference fluid name
785  std::string fluid_name = component.transport.conductivity_ecs.reference_fluid;
786  std::vector<std::string> name(1, fluid_name);
787  // Get a managed pointer to the reference fluid for ECS
788  shared_ptr<HelmholtzEOSMixtureBackend> ref_fluid(new HelmholtzEOSMixtureBackend(name));
789  // Get the viscosity using ECS and store in initial_density (not normally used);
790  initial_density = TransportRoutines::conductivity_ECS(*this, *ref_fluid); // Warning: not actually initial_density
791  return;
792  }
793 
794  if (component.transport.hardcoded_conductivity != CoolProp::TransportPropertyData::CONDUCTIVITY_NOT_HARDCODED) {
795  // Evaluate hardcoded model and deposit in initial_density variable
796  // Warning: not actually initial_density
797  switch (component.transport.hardcoded_conductivity) {
798  case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_WATER:
799  initial_density = TransportRoutines::conductivity_hardcoded_water(*this);
800  break;
801  case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_HEAVYWATER:
802  initial_density = TransportRoutines::conductivity_hardcoded_heavywater(*this);
803  break;
804  case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_R23:
805  initial_density = TransportRoutines::conductivity_hardcoded_R23(*this);
806  break;
807  case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_HELIUM:
808  initial_density = TransportRoutines::conductivity_hardcoded_helium(*this);
809  break;
810  case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_METHANE:
811  initial_density = TransportRoutines::conductivity_hardcoded_methane(*this);
812  break;
813  default:
814  throw ValueError(format("hardcoded conductivity type [%d] is invalid for fluid %s",
815  components[0].transport.hardcoded_conductivity, name().c_str()));
816  }
817  return;
818  }
819 
820  // -------------------------
821  // Normal evaluation
822  // -------------------------
823 
824  // Dilute part
825  switch (component.transport.conductivity_dilute.type) {
826  case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_RATIO_POLYNOMIALS:
827  dilute = TransportRoutines::conductivity_dilute_ratio_polynomials(*this);
828  break;
829  case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_ETA0_AND_POLY:
830  dilute = TransportRoutines::conductivity_dilute_eta0_and_poly(*this);
831  break;
832  case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_CO2:
833  dilute = TransportRoutines::conductivity_dilute_hardcoded_CO2(*this);
834  break;
835  case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_CO2_HUBER_JPCRD_2016:
836  dilute = TransportRoutines::conductivity_dilute_hardcoded_CO2_HuberJPCRD2016(*this);
837  break;
838  case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_ETHANE:
839  dilute = TransportRoutines::conductivity_dilute_hardcoded_ethane(*this);
840  break;
841  case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_NONE:
842  dilute = 0.0;
843  break;
844  default:
845  throw ValueError(
846  format("dilute conductivity type [%d] is invalid for fluid %s", components[0].transport.conductivity_dilute.type, name().c_str()));
847  }
848 
849  // Residual part
850  residual = calc_conductivity_background();
851 
852  // Critical part
853  switch (component.transport.conductivity_critical.type) {
854  case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_SIMPLIFIED_OLCHOWY_SENGERS:
855  critical = TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(*this);
856  break;
857  case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_R123:
858  critical = TransportRoutines::conductivity_critical_hardcoded_R123(*this);
859  break;
860  case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_AMMONIA:
861  critical = TransportRoutines::conductivity_critical_hardcoded_ammonia(*this);
862  break;
863  case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_NONE:
864  critical = 0.0;
865  break;
866  case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_CARBONDIOXIDE_SCALABRIN_JPCRD_2006:
867  critical = TransportRoutines::conductivity_critical_hardcoded_CO2_ScalabrinJPCRD2006(*this);
868  break;
869  default:
870  throw ValueError(
871  format("critical conductivity type [%d] is invalid for fluid %s", components[0].transport.viscosity_dilute.type, name().c_str()));
872  }
873  } else {
874  throw ValueError("calc_conductivity_contributions invalid for mixtures");
875  }
876 };
877 
878 CoolPropDbl HelmholtzEOSMixtureBackend::calc_conductivity_background(void) {
879  // Residual part
880  CoolPropDbl lambda_residual = _HUGE;
881  switch (components[0].transport.conductivity_residual.type) {
882  case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL:
883  lambda_residual = TransportRoutines::conductivity_residual_polynomial(*this);
884  break;
885  case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL_AND_EXPONENTIAL:
886  lambda_residual = TransportRoutines::conductivity_residual_polynomial_and_exponential(*this);
887  break;
888  default:
889  throw ValueError(
890  format("residual conductivity type [%d] is invalid for fluid %s", components[0].transport.conductivity_residual.type, name().c_str()));
891  }
892  return lambda_residual;
893 }
894 CoolPropDbl HelmholtzEOSMixtureBackend::calc_conductivity(void) {
895  if (is_pure_or_pseudopure) {
896  CoolPropDbl dilute = 0, initial_density = 0, residual = 0, critical = 0;
897  calc_conductivity_contributions(dilute, initial_density, residual, critical);
898  return dilute + initial_density + residual + critical;
899  } else {
900  set_warning_string("Mixture model for conductivity is highly approximate");
901  CoolPropDbl summer = 0;
902  for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
903  shared_ptr<HelmholtzEOSBackend> HEOS(new HelmholtzEOSBackend(components[i]));
904  HEOS->update(DmolarT_INPUTS, _rhomolar, _T);
905  summer += mole_fractions[i] * HEOS->conductivity();
906  }
907  return summer;
908  }
909 }
910 void HelmholtzEOSMixtureBackend::calc_conformal_state(const std::string& reference_fluid, CoolPropDbl& T, CoolPropDbl& rhomolar) {
911  shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> REF(new CoolProp::HelmholtzEOSBackend(reference_fluid));
912 
913  if (T < 0 && rhomolar < 0) {
914  // Collect some parameters
915  CoolPropDbl Tc = T_critical(), Tc0 = REF->T_critical(), rhocmolar = rhomolar_critical(), rhocmolar0 = REF->rhomolar_critical();
916 
917  // Starting guess values for shape factors
918  CoolPropDbl theta = 1;
919  CoolPropDbl phi = 1;
920 
921  // The equivalent substance reducing ratios
922  CoolPropDbl f = Tc / Tc0 * theta;
923  CoolPropDbl h = rhocmolar0 / rhocmolar * phi; // Must be the ratio of MOLAR densities!!
924 
925  // Starting guesses for conformal state
926  T = this->T() / f;
927  rhomolar = this->rhomolar() * h;
928  }
929 
930  TransportRoutines::conformal_state_solver(*this, *REF, T, rhomolar);
931 }
932 CoolPropDbl HelmholtzEOSMixtureBackend::calc_Ttriple(void) {
933  double summer = 0;
934  for (unsigned int i = 0; i < components.size(); ++i) {
935  summer += mole_fractions[i] * components[i].EOS().Ttriple;
936  }
937  return summer;
938 }
939 CoolPropDbl HelmholtzEOSMixtureBackend::calc_p_triple(void) {
940  double summer = 0;
941  for (unsigned int i = 0; i < components.size(); ++i) {
942  summer += mole_fractions[i] * components[i].EOS().ptriple;
943  }
944  return summer;
945 }
946 std::string HelmholtzEOSMixtureBackend::calc_name(void) {
947  if (components.size() != 1) {
948  throw ValueError(format("calc_name is only valid for pure and pseudo-pure fluids, %d components", components.size()));
949  } else {
950  return components[0].name;
951  }
952 }
953 
954 CoolPropDbl HelmholtzEOSMixtureBackend::calc_saturated_liquid_keyed_output(parameters key) {
955  if ((key == iDmolar) && _rhoLmolar) return _rhoLmolar;
956  if (!SatL) throw ValueError("The saturated liquid state has not been set.");
957  return SatL->keyed_output(key);
958 }
959 CoolPropDbl HelmholtzEOSMixtureBackend::calc_saturated_vapor_keyed_output(parameters key) {
960  if ((key == iDmolar) && _rhoVmolar) return _rhoVmolar;
961  if (!SatV) throw ValueError("The saturated vapor state has not been set.");
962  return SatV->keyed_output(key);
963 }
964 
965 void HelmholtzEOSMixtureBackend::calc_ideal_curve(const std::string& type, std::vector<double>& T, std::vector<double>& p) {
966  if (type == "Joule-Thomson") {
967  JouleThomsonCurveTracer JTCT(this, 1e5, 800);
968  JTCT.trace(T, p);
969  } else if (type == "Joule-Inversion") {
970  JouleInversionCurveTracer JICT(this, 1e5, 800);
971  JICT.trace(T, p);
972  } else if (type == "Ideal") {
973  IdealCurveTracer ICT(this, 1e5, 800);
974  ICT.trace(T, p);
975  } else if (type == "Boyle") {
976  BoyleCurveTracer BCT(this, 1e5, 800);
977  BCT.trace(T, p);
978  } else {
979  throw ValueError(format("Invalid ideal curve type: %s", type.c_str()));
980  }
981 };
982 std::vector<std::string> HelmholtzEOSMixtureBackend::calc_fluid_names(void) {
983  std::vector<std::string> out;
984  for (std::size_t i = 0; i < components.size(); ++i) {
985  out.push_back(components[i].name);
986  }
987  return out;
988 }
989 CoolPropDbl HelmholtzEOSMixtureBackend::calc_ODP(void) {
990  if (components.size() != 1) {
991  throw ValueError(format("For now, calc_ODP is only valid for pure and pseudo-pure fluids, %d components", components.size()));
992  } else {
993  CoolPropDbl v = components[0].environment.ODP;
994  if (!ValidNumber(v) || v < 0) {
995  throw ValueError(format("ODP value is not specified or invalid"));
996  }
997  return v;
998  }
999 }
1000 CoolPropDbl HelmholtzEOSMixtureBackend::calc_GWP20(void) {
1001  if (components.size() != 1) {
1002  throw ValueError(format("For now, calc_GWP20 is only valid for pure and pseudo-pure fluids, %d components", components.size()));
1003  } else {
1004  CoolPropDbl v = components[0].environment.GWP20;
1005  if (!ValidNumber(v) || v < 0) {
1006  throw ValueError(format("GWP20 value is not specified or invalid"));
1007  }
1008  return v;
1009  }
1010 }
1011 CoolPropDbl HelmholtzEOSMixtureBackend::calc_GWP100(void) {
1012  if (components.size() != 1) {
1013  throw ValueError(format("For now, calc_GWP100 is only valid for pure and pseudo-pure fluids, %d components", components.size()));
1014  } else {
1015  CoolPropDbl v = components[0].environment.GWP100;
1016  if (!ValidNumber(v) || v < 0) {
1017  throw ValueError(format("GWP100 value is not specified or invalid"));
1018  }
1019  return v;
1020  }
1021 }
1022 CoolPropDbl HelmholtzEOSMixtureBackend::calc_GWP500(void) {
1023  if (components.size() != 1) {
1024  throw ValueError(format("For now, calc_GWP500 is only valid for pure and pseudo-pure fluids, %d components", components.size()));
1025  } else {
1026  CoolPropDbl v = components[0].environment.GWP500;
1027  if (!ValidNumber(v) || v < 0) {
1028  throw ValueError(format("GWP500 value is not specified or invalid"));
1029  }
1030  return v;
1031  }
1032 }
1033 CoolPropDbl HelmholtzEOSMixtureBackend::calc_T_critical(void) {
1034  if (components.size() != 1) {
1035  std::vector<CriticalState> critpts = calc_all_critical_points();
1036  if (critpts.size() == 1) {
1037  //if (!critpts[0].stable){ throw ValueError(format("found one critical point but critical point is not stable")); }
1038  return critpts[0].T;
1039  } else {
1040  throw ValueError(format("critical point finding routine found %d critical points", critpts.size()));
1041  }
1042  } else {
1043  return components[0].crit.T;
1044  }
1045 }
1046 CoolPropDbl HelmholtzEOSMixtureBackend::calc_p_critical(void) {
1047  if (components.size() != 1) {
1048  std::vector<CriticalState> critpts = calc_all_critical_points();
1049  if (critpts.size() == 1) {
1050  //if (!critpts[0].stable){ throw ValueError(format("found one critical point but critical point is not stable")); }
1051  return critpts[0].p;
1052  } else {
1053  throw ValueError(format("critical point finding routine found %d critical points", critpts.size()));
1054  }
1055  } else {
1056  return components[0].crit.p;
1057  }
1058 }
1059 CoolPropDbl HelmholtzEOSMixtureBackend::calc_rhomolar_critical(void) {
1060  if (components.size() != 1) {
1061  std::vector<CriticalState> critpts = calc_all_critical_points();
1062  if (critpts.size() == 1) {
1063  //if (!critpts[0].stable){ throw ValueError(format("found one critical point but critical point is not stable")); }
1064  return critpts[0].rhomolar;
1065  } else {
1066  throw ValueError(format("critical point finding routine found %d critical points", critpts.size()));
1067  }
1068  } else {
1069  return components[0].crit.rhomolar;
1070  }
1071 }
1072 CoolPropDbl HelmholtzEOSMixtureBackend::calc_pmax_sat(void) {
1073  if (is_pure_or_pseudopure) {
1074  if (components[0].EOS().pseudo_pure) {
1075  return components[0].EOS().max_sat_p.p;
1076  } else {
1077  return p_critical();
1078  }
1079  } else {
1080  throw ValueError("calc_pmax_sat not yet defined for mixtures");
1081  }
1082 }
1083 CoolPropDbl HelmholtzEOSMixtureBackend::calc_Tmax_sat(void) {
1084  if (is_pure_or_pseudopure) {
1085  if (components[0].EOS().pseudo_pure) {
1086  double Tmax_sat = components[0].EOS().max_sat_T.T;
1087  if (!ValidNumber(Tmax_sat)) {
1088  return T_critical();
1089  } else {
1090  return Tmax_sat;
1091  }
1092  } else {
1093  return T_critical();
1094  }
1095  } else {
1096  throw ValueError("calc_Tmax_sat not yet defined for mixtures");
1097  }
1098 }
1099 
1100 void HelmholtzEOSMixtureBackend::calc_Tmin_sat(CoolPropDbl& Tmin_satL, CoolPropDbl& Tmin_satV) {
1101  if (is_pure_or_pseudopure) {
1102  Tmin_satL = components[0].EOS().sat_min_liquid.T;
1103  Tmin_satV = components[0].EOS().sat_min_vapor.T;
1104  return;
1105  } else {
1106  throw ValueError("calc_Tmin_sat not yet defined for mixtures");
1107  }
1108 }
1109 
1110 void HelmholtzEOSMixtureBackend::calc_pmin_sat(CoolPropDbl& pmin_satL, CoolPropDbl& pmin_satV) {
1111  if (is_pure_or_pseudopure) {
1112  pmin_satL = components[0].EOS().sat_min_liquid.p;
1113  pmin_satV = components[0].EOS().sat_min_vapor.p;
1114  return;
1115  } else {
1116  throw ValueError("calc_pmin_sat not yet defined for mixtures");
1117  }
1118 }
1119 
1120 // Minimum allowed saturation temperature the maximum of the saturation temperatures of liquid and vapor
1121 // For pure fluids, both values are the same, for pseudo-pure they are probably the same, for mixtures they are definitely not the same
1122 
1123 CoolPropDbl HelmholtzEOSMixtureBackend::calc_Tmax(void) {
1124  double summer = 0;
1125  for (unsigned int i = 0; i < components.size(); ++i) {
1126  summer += mole_fractions[i] * components[i].EOS().limits.Tmax;
1127  }
1128  return summer;
1129 }
1130 CoolPropDbl HelmholtzEOSMixtureBackend::calc_Tmin(void) {
1131  double summer = 0;
1132  for (unsigned int i = 0; i < components.size(); ++i) {
1133  summer += mole_fractions[i] * components[i].EOS().limits.Tmin;
1134  }
1135  return summer;
1136 }
1137 CoolPropDbl HelmholtzEOSMixtureBackend::calc_pmax(void) {
1138  double summer = 0;
1139  for (unsigned int i = 0; i < components.size(); ++i) {
1140  summer += mole_fractions[i] * components[i].EOS().limits.pmax;
1141  }
1142  return summer;
1143 }
1144 
1145 void HelmholtzEOSMixtureBackend::update_DmolarT_direct(CoolPropDbl rhomolar, CoolPropDbl T) {
1146  // TODO: This is just a quick fix for #878 - should be done more systematically
1147  const CoolPropDbl rhomolar_min = 0;
1148  const CoolPropDbl T_min = 0;
1149 
1150  if (rhomolar < rhomolar_min) {
1151  throw ValueError(format("The molar density of %f mol/m3 is below the minimum of %f mol/m3", rhomolar, rhomolar_min));
1152  }
1153 
1154  if (T < T_min) {
1155  throw ValueError(format("The temperature of %f K is below the minimum of %f K", T, T_min));
1156  }
1157 
1158  CoolProp::input_pairs pair = DmolarT_INPUTS;
1159  // Set up the state
1160  pre_update(pair, rhomolar, T);
1161 
1162  _rhomolar = rhomolar;
1163  _T = T;
1164  _p = calc_pressure();
1165 
1166  // Cleanup
1167  bool optional_checks = false;
1168  post_update(optional_checks);
1169 }
1170 
1171 void HelmholtzEOSMixtureBackend::update_HmolarQ_with_guessT(CoolPropDbl hmolar, CoolPropDbl Q, CoolPropDbl Tguess) {
1172  CoolProp::input_pairs pair = CoolProp::HmolarQ_INPUTS;
1173  // Set up the state
1174  pre_update(pair, hmolar, Q);
1175 
1176  _hmolar = hmolar;
1177  _Q = Q;
1178  FlashRoutines::HQ_flash(*this, Tguess);
1179 
1180  // Cleanup
1181  post_update();
1182 }
1183 void HelmholtzEOSMixtureBackend::update_internal(HelmholtzEOSMixtureBackend& HEOS) {
1184  this->_hmolar = HEOS.hmolar();
1185  this->_smolar = HEOS.smolar();
1186  this->_T = HEOS.T();
1187  this->_umolar = HEOS.umolar();
1188  this->_p = HEOS.p();
1189  this->_rhomolar = HEOS.rhomolar();
1190  this->_Q = HEOS.Q();
1191  this->_phase = HEOS.phase();
1192 
1193  // Copy the derivatives as well
1194 }
1195 void HelmholtzEOSMixtureBackend::update_TP_guessrho(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomolar_guess) {
1196  CoolProp::input_pairs pair = PT_INPUTS;
1197  // Set up the state
1198  pre_update(pair, p, T);
1199 
1200  // Do the flash call to find rho = f(T,p)
1201  CoolPropDbl rhomolar = solver_rho_Tp(T, p, rhomolar_guess);
1202 
1203  // Update the class with the new calculated density
1204  update_DmolarT_direct(rhomolar, T);
1205 
1206  // Skip the cleanup, already done in update_DmolarT_direct
1207 }
1208 
1209 void HelmholtzEOSMixtureBackend::pre_update(CoolProp::input_pairs& input_pair, CoolPropDbl& value1, CoolPropDbl& value2) {
1210  // Clear the state
1211  clear();
1212 
1213  if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) {
1214  throw ValueError("Mole fractions must be set");
1215  }
1216 
1217  // If the inputs are in mass units, convert them to molar units
1218  mass_to_molar_inputs(input_pair, value1, value2);
1219 
1220  // Set the mole-fraction weighted gas constant for the mixture
1221  // (or the pure/pseudo-pure fluid) if it hasn't been set yet
1222  gas_constant();
1223 
1224  // Calculate and cache the reducing state
1225  calc_reducing_state();
1226 }
1227 
1228 void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) {
1229  if (get_debug_level() > 10) {
1230  std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
1231  get_input_pair_short_desc(input_pair).c_str(), value1, value2)
1232  << std::endl;
1233  }
1234 
1235  CoolPropDbl ld_value1 = value1, ld_value2 = value2;
1236  pre_update(input_pair, ld_value1, ld_value2);
1237  value1 = ld_value1;
1238  value2 = ld_value2;
1239 
1240  switch (input_pair) {
1241  case PT_INPUTS:
1242  _p = value1;
1243  _T = value2;
1244  FlashRoutines::PT_flash(*this);
1245  break;
1246  case DmolarT_INPUTS:
1247  _rhomolar = value1;
1248  _T = value2;
1249  FlashRoutines::DHSU_T_flash(*this, iDmolar);
1250  break;
1251  case SmolarT_INPUTS:
1252  _smolar = value1;
1253  _T = value2;
1254  FlashRoutines::DHSU_T_flash(*this, iSmolar);
1255  break;
1256  //case HmolarT_INPUTS:
1257  // _hmolar = value1; _T = value2; FlashRoutines::DHSU_T_flash(*this, iHmolar); break;
1258  //case TUmolar_INPUTS:
1259  // _T = value1; _umolar = value2; FlashRoutines::DHSU_T_flash(*this, iUmolar); break;
1260  case DmolarP_INPUTS:
1261  _rhomolar = value1;
1262  _p = value2;
1263  FlashRoutines::DP_flash(*this);
1264  break;
1265  case DmolarHmolar_INPUTS:
1266  _rhomolar = value1;
1267  _hmolar = value2;
1268  FlashRoutines::HSU_D_flash(*this, iHmolar);
1269  break;
1270  case DmolarSmolar_INPUTS:
1271  _rhomolar = value1;
1272  _smolar = value2;
1273  FlashRoutines::HSU_D_flash(*this, iSmolar);
1274  break;
1275  case DmolarUmolar_INPUTS:
1276  _rhomolar = value1;
1277  _umolar = value2;
1278  FlashRoutines::HSU_D_flash(*this, iUmolar);
1279  break;
1280  case HmolarP_INPUTS:
1281  _hmolar = value1;
1282  _p = value2;
1283  FlashRoutines::HSU_P_flash(*this, iHmolar);
1284  break;
1285  case PSmolar_INPUTS:
1286  _p = value1;
1287  _smolar = value2;
1288  FlashRoutines::HSU_P_flash(*this, iSmolar);
1289  break;
1290  case PUmolar_INPUTS:
1291  _p = value1;
1292  _umolar = value2;
1293  FlashRoutines::HSU_P_flash(*this, iUmolar);
1294  break;
1295  case HmolarSmolar_INPUTS:
1296  _hmolar = value1;
1297  _smolar = value2;
1298  FlashRoutines::HS_flash(*this);
1299  break;
1300  case QT_INPUTS:
1301  _Q = value1;
1302  _T = value2;
1303  if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1304  FlashRoutines::QT_flash(*this);
1305  break;
1306  case PQ_INPUTS:
1307  _p = value1;
1308  _Q = value2;
1309  if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1310  FlashRoutines::PQ_flash(*this);
1311  break;
1312  case QSmolar_INPUTS:
1313  _Q = value1;
1314  _smolar = value2;
1315  if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1316  FlashRoutines::QS_flash(*this);
1317  break;
1318  case HmolarQ_INPUTS:
1319  _hmolar = value1;
1320  _Q = value2;
1321  if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1322  FlashRoutines::HQ_flash(*this);
1323  break;
1324  case DmolarQ_INPUTS:
1325  _rhomolar = value1;
1326  _Q = value2;
1327  if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1328  FlashRoutines::DQ_flash(*this);
1329  break;
1330  default:
1331  throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
1332  }
1333 
1334  post_update();
1335 }
1336 const std::vector<CoolPropDbl> HelmholtzEOSMixtureBackend::calc_mass_fractions() {
1337  // mass fraction is mass_i/total_mass;
1338  CoolPropDbl mm = molar_mass();
1339  std::vector<CoolPropDbl>& mole_fractions = get_mole_fractions_ref();
1340  std::vector<CoolPropDbl> mass_fractions(mole_fractions.size());
1341  for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
1342  double mmi = get_fluid_constant(i, imolar_mass);
1343  mass_fractions[i] = mmi * (mole_fractions[i]) / mm;
1344  }
1345  return mass_fractions;
1346 }
1347 
1348 void HelmholtzEOSMixtureBackend::update_with_guesses(CoolProp::input_pairs input_pair, double value1, double value2,
1349  const GuessesStructure& guesses) {
1350  if (get_debug_level() > 10) {
1351  std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
1352  get_input_pair_short_desc(input_pair).c_str(), value1, value2)
1353  << std::endl;
1354  }
1355 
1356  CoolPropDbl ld_value1 = value1, ld_value2 = value2;
1357  pre_update(input_pair, ld_value1, ld_value2);
1358  value1 = ld_value1;
1359  value2 = ld_value2;
1360 
1361  switch (input_pair) {
1362  case PQ_INPUTS:
1363  _p = value1;
1364  _Q = value2;
1365  FlashRoutines::PQ_flash_with_guesses(*this, guesses);
1366  break;
1367  case QT_INPUTS:
1368  _Q = value1;
1369  _T = value2;
1370  FlashRoutines::QT_flash_with_guesses(*this, guesses);
1371  break;
1372  case PT_INPUTS:
1373  _p = value1;
1374  _T = value2;
1375  FlashRoutines::PT_flash_with_guesses(*this, guesses);
1376  break;
1377  default:
1378  throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
1379  }
1380  post_update();
1381 }
1382 
1383 void HelmholtzEOSMixtureBackend::post_update(bool optional_checks) {
1384  // Check the values that must always be set
1385  //if (_p < 0){ throw ValueError("p is less than zero");}
1386  if (!ValidNumber(_p)) {
1387  throw ValueError("p is not a valid number");
1388  }
1389  //if (_T < 0){ throw ValueError("T is less than zero");}
1390  if (!ValidNumber(_T)) {
1391  throw ValueError("T is not a valid number");
1392  }
1393  if (_rhomolar < 0) {
1394  throw ValueError("rhomolar is less than zero");
1395  }
1396  if (!ValidNumber(_rhomolar)) {
1397  throw ValueError("rhomolar is not a valid number");
1398  }
1399 
1400  if (optional_checks) {
1401  if (!ValidNumber(_Q)) {
1402  throw ValueError("Q is not a valid number");
1403  }
1404  if (_phase == iphase_unknown) {
1405  throw ValueError("_phase is unknown");
1406  }
1407  }
1408 
1409  // Set the reduced variables
1410  _tau = _reducing.T / _T;
1411  _delta = _rhomolar / _reducing.rhomolar;
1412 
1413  // Update the terms in the excess contribution
1414  residual_helmholtz->Excess.update(_tau, _delta);
1415 }
1416 
1417 CoolPropDbl HelmholtzEOSMixtureBackend::calc_Bvirial() {
1418  return 1 / rhomolar_reducing() * calc_alphar_deriv_nocache(0, 1, mole_fractions, _tau, 1e-12);
1419 }
1420 CoolPropDbl HelmholtzEOSMixtureBackend::calc_dBvirial_dT() {
1421  SimpleState red = get_reducing_state();
1422  CoolPropDbl dtau_dT = -red.T / pow(_T, 2);
1423  return 1 / red.rhomolar * calc_alphar_deriv_nocache(1, 1, mole_fractions, _tau, 1e-12) * dtau_dT;
1424 }
1425 CoolPropDbl HelmholtzEOSMixtureBackend::calc_Cvirial() {
1426  return 1 / pow(rhomolar_reducing(), 2) * calc_alphar_deriv_nocache(0, 2, mole_fractions, _tau, 1e-12);
1427 }
1428 CoolPropDbl HelmholtzEOSMixtureBackend::calc_dCvirial_dT() {
1429  SimpleState red = get_reducing_state();
1430  CoolPropDbl dtau_dT = -red.T / pow(_T, 2);
1431  return 1 / pow(red.rhomolar, 2) * calc_alphar_deriv_nocache(1, 2, mole_fractions, _tau, 1e-12) * dtau_dT;
1432 }
1433 void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int other, CoolPropDbl value, bool& saturation_called) {
1434  /*
1435  Determine the phase given p and one other state variable
1436  */
1437  saturation_called = false;
1438 
1439  // Reference declaration to save indexing
1440  CoolPropFluid& component = components[0];
1441 
1442  // Maximum saturation temperature - Equal to critical pressure for pure fluids
1443  CoolPropDbl psat_max = calc_pmax_sat();
1444 
1445  // Check supercritical pressure
1446  if (_p > psat_max) {
1447  _Q = 1e9;
1448  switch (other) {
1449  case iT: {
1450  if (_T > _crit.T) {
1451  this->_phase = iphase_supercritical;
1452  return;
1453  } else {
1454  this->_phase = iphase_supercritical_liquid;
1455  return;
1456  }
1457  }
1458  case iDmolar: {
1459  if (_rhomolar < _crit.rhomolar) {
1460  this->_phase = iphase_supercritical_gas;
1461  return;
1462  } else {
1463  this->_phase = iphase_supercritical_liquid;
1464  return;
1465  }
1466  }
1467  case iSmolar: {
1468  if (_smolar.pt() > _crit.smolar) {
1469  this->_phase = iphase_supercritical_gas;
1470  return;
1471  } else {
1472  this->_phase = iphase_supercritical_liquid;
1473  return;
1474  }
1475  }
1476  case iHmolar: {
1477  if (_hmolar.pt() > _crit.hmolar) {
1478  this->_phase = iphase_supercritical_gas;
1479  return;
1480  } else {
1481  this->_phase = iphase_supercritical_liquid;
1482  return;
1483  }
1484  }
1485  case iUmolar: {
1486  if (_umolar.pt() > _crit.umolar) {
1487  this->_phase = iphase_supercritical_gas;
1488  return;
1489  } else {
1490  this->_phase = iphase_supercritical_liquid;
1491  return;
1492  }
1493  }
1494  default: {
1495  throw ValueError("supercritical pressure but other invalid for now");
1496  }
1497  }
1498  }
1499  // Check between triple point pressure and psat_max
1500  else if (_p >= components[0].EOS().ptriple * 0.9999 && _p <= psat_max) {
1501  // First try the ancillaries, use them to determine the state if you can
1502 
1503  // Calculate dew and bubble temps from the ancillaries (everything needs them)
1504  _TLanc = components[0].ancillaries.pL.invert(_p);
1505  _TVanc = components[0].ancillaries.pV.invert(_p);
1506 
1507  bool definitely_two_phase = false;
1508 
1509  // Try using the ancillaries for P,H,S if they are there
1510  switch (other) {
1511  case iT: {
1512 
1513  if (has_melting_line()) {
1514  double Tm = melting_line(iT, iP, _p);
1515  if (get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
1516  _phase = iphase_liquid;
1517  } else {
1518  if (_T < Tm - 0.001) {
1519  throw ValueError(format("For now, we don't support T [%g K] below Tmelt(p) [%g K]", _T, Tm));
1520  }
1521  }
1522  } else {
1523  if (get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
1524  _phase = iphase_liquid;
1525  } else {
1526  if (_T < Tmin() - 0.001) {
1527  throw ValueError(format("For now, we don't support T [%g K] below Tmin(saturation) [%g K]", _T, Tmin()));
1528  }
1529  }
1530  }
1531 
1532  CoolPropDbl T_vap = 0.1 + static_cast<double>(_TVanc);
1533  CoolPropDbl T_liq = -0.1 + static_cast<double>(_TLanc);
1534 
1535  if (value > T_vap) {
1536  this->_phase = iphase_gas;
1537  _Q = -1000;
1538  return;
1539  } else if (value < T_liq) {
1540  this->_phase = iphase_liquid;
1541  _Q = 1000;
1542  return;
1543  }
1544  break;
1545  }
1546  case iHmolar: {
1547  if (!component.ancillaries.hL.enabled()) {
1548  break;
1549  }
1550  // Ancillaries are h-h_anchor, so add back h_anchor
1551  CoolPropDbl h_liq = component.ancillaries.hL.evaluate(_TLanc) + component.EOS().hs_anchor.hmolar;
1552  CoolPropDbl h_liq_error_band = component.ancillaries.hL.get_max_abs_error();
1553  CoolPropDbl h_vap = h_liq + component.ancillaries.hLV.evaluate(_TLanc);
1554  CoolPropDbl h_vap_error_band = h_liq_error_band + component.ancillaries.hLV.get_max_abs_error();
1555 
1556  // HelmholtzEOSMixtureBackend HEOS(components);
1557  // HEOS.update(QT_INPUTS, 1, _TLanc);
1558  // double h1 = HEOS.hmolar();
1559  // HEOS.update(QT_INPUTS, 0, _TLanc);
1560  // double h0 = HEOS.hmolar();
1561 
1562  // Check if in range given the accuracy of the fit
1563  if (value > h_vap + h_vap_error_band) {
1564  this->_phase = iphase_gas;
1565  _Q = -1000;
1566  return;
1567  } else if (value < h_liq - h_liq_error_band) {
1568  this->_phase = iphase_liquid;
1569  _Q = 1000;
1570  return;
1571  } else if (value > h_liq + h_liq_error_band && value < h_vap - h_vap_error_band) {
1572  definitely_two_phase = true;
1573  }
1574  break;
1575  }
1576  case iSmolar: {
1577  if (!component.ancillaries.sL.enabled()) {
1578  break;
1579  }
1580  // Ancillaries are s-s_anchor, so add back s_anchor
1581  CoolPropDbl s_anchor = component.EOS().hs_anchor.smolar;
1582  CoolPropDbl s_liq = component.ancillaries.sL.evaluate(_TLanc) + s_anchor;
1583  CoolPropDbl s_liq_error_band = component.ancillaries.sL.get_max_abs_error();
1584  CoolPropDbl s_vap = s_liq + component.ancillaries.sLV.evaluate(_TVanc);
1585  CoolPropDbl s_vap_error_band = s_liq_error_band + component.ancillaries.sLV.get_max_abs_error();
1586 
1587  // Check if in range given the accuracy of the fit
1588  if (value > s_vap + s_vap_error_band) {
1589  this->_phase = iphase_gas;
1590  _Q = -1000;
1591  return;
1592  } else if (value < s_liq - s_liq_error_band) {
1593  this->_phase = iphase_liquid;
1594  _Q = 1000;
1595  return;
1596  } else if (value > s_liq + s_liq_error_band && value < s_vap - s_vap_error_band) {
1597  definitely_two_phase = true;
1598  }
1599  break;
1600  }
1601  case iUmolar: {
1602  if (!component.ancillaries.hL.enabled()) {
1603  break;
1604  }
1605  // u = h-p/rho
1606 
1607  // Ancillaries are h-h_anchor, so add back h_anchor
1608  CoolPropDbl h_liq = component.ancillaries.hL.evaluate(_TLanc) + component.EOS().hs_anchor.hmolar;
1609  CoolPropDbl h_liq_error_band = component.ancillaries.hL.get_max_abs_error();
1610  CoolPropDbl h_vap = h_liq + component.ancillaries.hLV.evaluate(_TLanc);
1611  CoolPropDbl h_vap_error_band = h_liq_error_band + component.ancillaries.hLV.get_max_abs_error();
1612  CoolPropDbl rho_vap = component.ancillaries.rhoV.evaluate(_TVanc);
1613  CoolPropDbl rho_liq = component.ancillaries.rhoL.evaluate(_TLanc);
1614  CoolPropDbl u_liq = h_liq - _p / rho_liq;
1615  CoolPropDbl u_vap = h_vap - _p / rho_vap;
1616  CoolPropDbl u_liq_error_band = 1.5 * h_liq_error_band; // Most of error is in enthalpy
1617  CoolPropDbl u_vap_error_band = 1.5 * h_vap_error_band; // Most of error is in enthalpy
1618 
1619  // Check if in range given the accuracy of the fit
1620  if (value > u_vap + u_vap_error_band) {
1621  this->_phase = iphase_gas;
1622  _Q = -1000;
1623  return;
1624  } else if (value < u_liq - u_liq_error_band) {
1625  this->_phase = iphase_liquid;
1626  _Q = 1000;
1627  return;
1628  } else if (value > u_liq + u_liq_error_band && value < u_vap - u_vap_error_band) {
1629  definitely_two_phase = true;
1630  }
1631  break;
1632  }
1633  default: {
1634  }
1635  }
1636 
1637  // Now either density is an input, or an ancillary for h,s,u is missing
1638  // Always calculate the densities using the ancillaries
1639  if (!definitely_two_phase) {
1640  _rhoVanc = component.ancillaries.rhoV.evaluate(_TVanc);
1641  _rhoLanc = component.ancillaries.rhoL.evaluate(_TLanc);
1642  CoolPropDbl rho_vap = 0.95 * static_cast<double>(_rhoVanc);
1643  CoolPropDbl rho_liq = 1.05 * static_cast<double>(_rhoLanc);
1644  switch (other) {
1645  case iDmolar: {
1646  if (value < rho_vap) {
1647  this->_phase = iphase_gas;
1648  return;
1649  } else if (value > rho_liq) {
1650  this->_phase = iphase_liquid;
1651  return;
1652  }
1653  break;
1654  }
1655  }
1656  }
1657 
1658  if (!is_pure_or_pseudopure) {
1659  throw ValueError("possibly two-phase inputs not supported for mixtures for now");
1660  }
1661 
1662  // Actually have to use saturation information sadly
1663  // For the given pressure, find the saturation state
1664  // Run the saturation routines to determine the saturation densities and pressures
1665  HelmholtzEOSMixtureBackend HEOS(components);
1666  HEOS._p = this->_p;
1667  HEOS._Q = 0; // ?? What is the best to do here? Doesn't matter for our purposes since pure fluid
1668  FlashRoutines::PQ_flash(HEOS);
1669 
1670  // We called the saturation routines, so HEOS.SatL and HEOS.SatV are now updated
1671  // with the saturated liquid and vapor values, which can therefore be used in
1672  // the other solvers
1673  saturation_called = true;
1674 
1675  CoolPropDbl Q;
1676 
1677  if (other == iT) {
1678  if (value < HEOS.SatL->T() - 100 * DBL_EPSILON) {
1679  this->_phase = iphase_liquid;
1680  _Q = -1000;
1681  return;
1682  } else if (value > HEOS.SatV->T() + 100 * DBL_EPSILON) {
1683  this->_phase = iphase_gas;
1684  _Q = 1000;
1685  return;
1686  } else {
1687  this->_phase = iphase_twophase;
1688  }
1689  }
1690  switch (other) {
1691  case iDmolar:
1692  Q = (1 / value - 1 / HEOS.SatL->rhomolar()) / (1 / HEOS.SatV->rhomolar() - 1 / HEOS.SatL->rhomolar());
1693  break;
1694  case iSmolar:
1695  Q = (value - HEOS.SatL->smolar()) / (HEOS.SatV->smolar() - HEOS.SatL->smolar());
1696  break;
1697  case iHmolar:
1698  Q = (value - HEOS.SatL->hmolar()) / (HEOS.SatV->hmolar() - HEOS.SatL->hmolar());
1699  break;
1700  case iUmolar:
1701  Q = (value - HEOS.SatL->umolar()) / (HEOS.SatV->umolar() - HEOS.SatL->umolar());
1702  break;
1703  default:
1704  throw ValueError(format("bad input for other"));
1705  }
1706  // TODO: Check the speed penalty of these calls
1707  // Update the states
1708  if (this->SatL) this->SatL->update(DmolarT_INPUTS, HEOS.SatL->rhomolar(), HEOS.SatL->T());
1709  if (this->SatV) this->SatV->update(DmolarT_INPUTS, HEOS.SatV->rhomolar(), HEOS.SatV->T());
1710  // Update the two-Phase variables
1711  _rhoLmolar = HEOS.SatL->rhomolar();
1712  _rhoVmolar = HEOS.SatV->rhomolar();
1713 
1714  //
1715  if (Q < -1e-9) {
1716  this->_phase = iphase_liquid;
1717  _Q = -1000;
1718  return;
1719  } else if (Q > 1 + 1e-9) {
1720  this->_phase = iphase_gas;
1721  _Q = 1000;
1722  return;
1723  } else {
1724  this->_phase = iphase_twophase;
1725  }
1726 
1727  _Q = Q;
1728  // Load the outputs
1729  _T = _Q * HEOS.SatV->T() + (1 - _Q) * HEOS.SatL->T();
1730  _rhomolar = 1 / (_Q / HEOS.SatV->rhomolar() + (1 - _Q) / HEOS.SatL->rhomolar());
1731  return;
1732  } else if (_p < components[0].EOS().ptriple * 0.9999) {
1733  if (other == iT) {
1734  if (_T > std::max(Tmin(), Ttriple())) {
1735  _phase = iphase_gas;
1736  } else {
1737  if (get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
1738  _phase = iphase_gas;
1739  } else {
1740  throw NotImplementedError(format("For now, we don't support p [%g Pa] below ptriple [%g Pa] when T [%g] is less than Tmin [%g]",
1741  _p, components[0].EOS().ptriple, _T, std::max(Tmin(), Ttriple())));
1742  }
1743  }
1744  } else {
1745  _phase = iphase_gas;
1746  }
1747  } else {
1748  throw ValueError(format("The pressure [%g Pa] cannot be used in p_phase_determination", _p));
1749  }
1750 }
1751 void HelmholtzEOSMixtureBackend::calc_ssat_max(void) {
1752  class Residual : public FuncWrapper1D
1753  {
1754  public:
1755  HelmholtzEOSMixtureBackend* HEOS;
1756  Residual(HelmholtzEOSMixtureBackend& HEOS) : HEOS(&HEOS){};
1757  double call(double T) {
1758  HEOS->update(QT_INPUTS, 1, T);
1759  // dTdp_along_sat
1760  double dTdp_along_sat =
1761  HEOS->T() * (1 / HEOS->SatV->rhomolar() - 1 / HEOS->SatL->rhomolar()) / (HEOS->SatV->hmolar() - HEOS->SatL->hmolar());
1762  // dsdT_along_sat;
1763  return HEOS->SatV->first_partial_deriv(iSmolar, iT, iP) + HEOS->SatV->first_partial_deriv(iSmolar, iP, iT) / dTdp_along_sat;
1764  }
1765  };
1766  if (!ssat_max.is_valid() && ssat_max.exists != SsatSimpleState::SSAT_MAX_DOESNT_EXIST) {
1767  shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS_copy(new CoolProp::HelmholtzEOSMixtureBackend(get_components()));
1768  Residual resid(*HEOS_copy);
1769  const CoolProp::SimpleState& tripleV = HEOS_copy->get_components()[0].triple_vapor;
1770  double v1 = resid.call(hsat_max.T);
1771  double v2 = resid.call(tripleV.T);
1772  // If there is a sign change, there is a maxima, otherwise there is no local maxima/minima
1773  if (v1 * v2 < 0) {
1774  Brent(resid, hsat_max.T, tripleV.T, DBL_EPSILON, 1e-8, 30);
1775  ssat_max.T = resid.HEOS->T();
1776  ssat_max.p = resid.HEOS->p();
1777  ssat_max.rhomolar = resid.HEOS->rhomolar();
1778  ssat_max.hmolar = resid.HEOS->hmolar();
1779  ssat_max.smolar = resid.HEOS->smolar();
1780  ssat_max.exists = SsatSimpleState::SSAT_MAX_DOES_EXIST;
1781  } else {
1782  ssat_max.exists = SsatSimpleState::SSAT_MAX_DOESNT_EXIST;
1783  }
1784  }
1785 }
1786 void HelmholtzEOSMixtureBackend::calc_hsat_max(void) {
1787  class Residualhmax : public FuncWrapper1D
1788  {
1789  public:
1790  HelmholtzEOSMixtureBackend* HEOS;
1791  Residualhmax(HelmholtzEOSMixtureBackend& HEOS) : HEOS(&HEOS){};
1792  double call(double T) {
1793  HEOS->update(QT_INPUTS, 1, T);
1794  // dTdp_along_sat
1795  double dTdp_along_sat =
1796  HEOS->T() * (1 / HEOS->SatV->rhomolar() - 1 / HEOS->SatL->rhomolar()) / (HEOS->SatV->hmolar() - HEOS->SatL->hmolar());
1797  // dhdT_along_sat;
1798  return HEOS->SatV->first_partial_deriv(iHmolar, iT, iP) + HEOS->SatV->first_partial_deriv(iHmolar, iP, iT) / dTdp_along_sat;
1799  }
1800  };
1801  if (!hsat_max.is_valid()) {
1802  shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS_copy(new CoolProp::HelmholtzEOSMixtureBackend(get_components()));
1803  Residualhmax residhmax(*HEOS_copy);
1804  Brent(residhmax, T_critical() - 0.1, HEOS_copy->Ttriple() + 1, DBL_EPSILON, 1e-8, 30);
1805  hsat_max.T = residhmax.HEOS->T();
1806  hsat_max.p = residhmax.HEOS->p();
1807  hsat_max.rhomolar = residhmax.HEOS->rhomolar();
1808  hsat_max.hmolar = residhmax.HEOS->hmolar();
1809  hsat_max.smolar = residhmax.HEOS->smolar();
1810  }
1811 }
1812 void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int other, CoolPropDbl value) {
1813  if (!ValidNumber(value)) {
1814  throw ValueError(format("value to T_phase_determination_pure_or_pseudopure is invalid"));
1815  };
1816 
1817  // T is known, another input P, T, H, S, U is given (all molar)
1818  if (_T < _crit.T && _p > _crit.p) {
1819  // Only ever true if (other = iP); otherwise _p = -HUGE
1820  _phase = iphase_supercritical_liquid;
1821  } else if (std::abs(_T - _crit.T) < 10 * DBL_EPSILON) // Exactly at Tcrit
1822  {
1823  switch (other) {
1824  case iDmolar:
1825  if (std::abs(_rhomolar - _crit.rhomolar) < 10 * DBL_EPSILON) {
1826  _phase = iphase_critical_point;
1827  break;
1828  } else if (_rhomolar > _crit.rhomolar) {
1829  _phase = iphase_supercritical_liquid;
1830  break;
1831  } else {
1832  _phase = iphase_supercritical_gas;
1833  break;
1834  }
1835  case iP: {
1836  if (std::abs(_p - _crit.p) < 10 * DBL_EPSILON) {
1837  _phase = iphase_critical_point;
1838  break;
1839  } else if (_p > _crit.p) {
1840  _phase = iphase_supercritical_liquid;
1841  break;
1842  } else {
1843  _phase = iphase_supercritical_gas;
1844  break;
1845  }
1846  }
1847  default:
1848  throw ValueError(format("T=Tcrit; invalid input for other to T_phase_determination_pure_or_pseudopure"));
1849  }
1850  } else if (_T < _crit.T) // Gas, 2-Phase, Liquid, or Supercritical Liquid Region
1851  {
1852  // Start to think about the saturation stuff
1853  // First try to use the ancillary equations if you are far enough away
1854  // You know how accurate the ancillary equations are thanks to using CoolProp code to refit them
1855  switch (other) {
1856  case iP: {
1857  _pLanc = components[0].ancillaries.pL.evaluate(_T);
1858  _pVanc = components[0].ancillaries.pV.evaluate(_T);
1859  CoolPropDbl p_vap = 0.98 * static_cast<double>(_pVanc);
1860  CoolPropDbl p_liq = 1.02 * static_cast<double>(_pLanc);
1861 
1862  if (value < p_vap) {
1863  this->_phase = iphase_gas;
1864  _Q = -1000;
1865  return;
1866  } else if (value > p_liq) {
1867  this->_phase = iphase_liquid;
1868  _Q = 1000;
1869  return;
1870  } else if (!is_pure()) // pseudo-pure
1871  {
1872  // For pseudo-pure fluids, the ancillary pressure curves are the official
1873  // arbiter of the phase
1874  if (value > static_cast<CoolPropDbl>(_pLanc)) {
1875  this->_phase = iphase_liquid;
1876  _Q = 1000;
1877  return;
1878  } else if (value < static_cast<CoolPropDbl>(_pVanc)) {
1879  this->_phase = iphase_gas;
1880  _Q = -1000;
1881  return;
1882  } else {
1883  throw ValueError("Two-phase inputs not supported for pseudo-pure for now");
1884  }
1885  }
1886  break;
1887  }
1888  default: {
1889  // Always calculate the densities using the ancillaries
1890  _rhoVanc = components[0].ancillaries.rhoV.evaluate(_T);
1891  _rhoLanc = components[0].ancillaries.rhoL.evaluate(_T);
1892  CoolPropDbl rho_vap = 0.95 * static_cast<double>(_rhoVanc);
1893  CoolPropDbl rho_liq = 1.05 * static_cast<double>(_rhoLanc);
1894  switch (other) {
1895  case iDmolar: {
1896  if (value < rho_vap) {
1897  this->_phase = iphase_gas;
1898  return;
1899  } else if (value > rho_liq) {
1900  this->_phase = iphase_liquid;
1901  return;
1902  } else {
1903  // Next we check the vapor quality based on the ancillary values
1904  double Qanc = (1 / value - 1 / static_cast<double>(_rhoLanc))
1905  / (1 / static_cast<double>(_rhoVanc) - 1 / static_cast<double>(_rhoLanc));
1906  // If the vapor quality is significantly inside the two-phase zone, stop, we are definitely two-phase
1907  if (value > 0.95 * rho_liq || value < 1.05 * rho_vap) {
1908  // Definitely single-phase
1909  _phase = iphase_liquid; // Needed for direct update call
1910  _Q = -1000; // Needed for direct update call
1911  update_DmolarT_direct(value, _T);
1912  CoolPropDbl pL = components[0].ancillaries.pL.evaluate(_T);
1913  if (Qanc < 0.01 && _p > pL * 1.05 && first_partial_deriv(iP, iDmolar, iT) > 0
1914  && second_partial_deriv(iP, iDmolar, iT, iDmolar, iT) > 0) {
1915  _phase = iphase_liquid;
1916  _Q = -1000;
1917  return;
1918  } else if (Qanc > 1.01) {
1919  break;
1920  } else {
1921  _phase = iphase_unknown;
1922  _p = _HUGE;
1923  }
1924  }
1925  }
1926  break;
1927  }
1928  default: {
1929  if (!this->SatL || !this->SatV) {
1930  throw ValueError(format("The saturation properties are needed in T_phase_determination_pure_or_pseudopure"));
1931  }
1932  // If it is not density, update the states
1933  SatV->update(DmolarT_INPUTS, rho_vap, _T);
1934  SatL->update(DmolarT_INPUTS, rho_liq, _T);
1935 
1936  // First we check ancillaries
1937  switch (other) {
1938  case iSmolar: {
1939  if (value > SatV->calc_smolar()) {
1940  this->_phase = iphase_gas;
1941  return;
1942  }
1943  if (value < SatL->calc_smolar()) {
1944  this->_phase = iphase_liquid;
1945  return;
1946  }
1947  break;
1948  }
1949  case iHmolar: {
1950  if (value > SatV->calc_hmolar()) {
1951  this->_phase = iphase_gas;
1952  return;
1953  } else if (value < SatL->calc_hmolar()) {
1954  this->_phase = iphase_liquid;
1955  return;
1956  }
1957  break;
1958  }
1959  case iUmolar: {
1960  if (value > SatV->calc_umolar()) {
1961  this->_phase = iphase_gas;
1962  return;
1963  } else if (value < SatL->calc_umolar()) {
1964  this->_phase = iphase_liquid;
1965  return;
1966  }
1967  break;
1968  }
1969  default:
1970  throw ValueError(format("invalid input for other to T_phase_determination_pure_or_pseudopure"));
1971  }
1972  }
1973  }
1974  }
1975  }
1976 
1977  // Actually have to use saturation information sadly
1978  // For the given temperature, find the saturation state
1979  // Run the saturation routines to determine the saturation densities and pressures
1980  HelmholtzEOSMixtureBackend HEOS(components);
1981  SaturationSolvers::saturation_T_pure_options options;
1982  SaturationSolvers::saturation_T_pure(HEOS, _T, options);
1983 
1984  CoolPropDbl Q;
1985 
1986  if (other == iP) {
1987  if (value > HEOS.SatL->p() * (1e-6 + 1)) {
1988  this->_phase = iphase_liquid;
1989  _Q = -1000;
1990  return;
1991  } else if (value < HEOS.SatV->p() * (1 - 1e-6)) {
1992  this->_phase = iphase_gas;
1993  _Q = 1000;
1994  return;
1995  } else {
1996  throw ValueError(
1997  format("Saturation pressure [%g Pa] corresponding to T [%g K] is within 1e-4 %% of given p [%Lg Pa]", HEOS.SatL->p(), _T, value));
1998  }
1999  }
2000 
2001  switch (other) {
2002  case iDmolar:
2003  Q = (1 / value - 1 / HEOS.SatL->rhomolar()) / (1 / HEOS.SatV->rhomolar() - 1 / HEOS.SatL->rhomolar());
2004  break;
2005  case iSmolar:
2006  Q = (value - HEOS.SatL->smolar()) / (HEOS.SatV->smolar() - HEOS.SatL->smolar());
2007  break;
2008  case iHmolar:
2009  Q = (value - HEOS.SatL->hmolar()) / (HEOS.SatV->hmolar() - HEOS.SatL->hmolar());
2010  break;
2011  case iUmolar:
2012  Q = (value - HEOS.SatL->umolar()) / (HEOS.SatV->umolar() - HEOS.SatL->umolar());
2013  break;
2014  default:
2015  throw ValueError(format("bad input for other"));
2016  }
2017 
2018  // Update the states
2019  if (this->SatL) this->SatL->update(DmolarT_INPUTS, HEOS.SatL->rhomolar(), HEOS.SatL->T());
2020  if (this->SatV) this->SatV->update(DmolarT_INPUTS, HEOS.SatV->rhomolar(), HEOS.SatV->T());
2021  // Update the two-Phase variables
2022  _rhoLmolar = HEOS.SatL->rhomolar();
2023  _rhoVmolar = HEOS.SatV->rhomolar();
2024 
2025  if (Q < 0) {
2026  this->_phase = iphase_liquid;
2027  _Q = -1;
2028  return;
2029  } else if (Q > 1) {
2030  this->_phase = iphase_gas;
2031  _Q = 1;
2032  return;
2033  } else {
2034  this->_phase = iphase_twophase;
2035  }
2036  _Q = Q;
2037  // Load the outputs
2038  _p = _Q * HEOS.SatV->p() + (1 - _Q) * HEOS.SatL->p();
2039  _rhomolar = 1 / (_Q / HEOS.SatV->rhomolar() + (1 - _Q) / HEOS.SatL->rhomolar());
2040  return;
2041  } else if (_T > _crit.T && _T > components[0].EOS().Ttriple) // Supercritical or Supercritical Gas Region
2042  {
2043  _Q = 1e9;
2044  switch (other) {
2045  case iP: {
2046  if (_p > _crit.p) {
2047  this->_phase = iphase_supercritical;
2048  return;
2049  } else {
2050  this->_phase = iphase_supercritical_gas;
2051  return;
2052  }
2053  }
2054  case iDmolar: {
2055  if (_rhomolar > _crit.rhomolar) {
2056  this->_phase = iphase_supercritical_liquid;
2057  return;
2058  } else {
2059  this->_phase = iphase_supercritical_gas;
2060  return;
2061  }
2062  }
2063  case iSmolar: {
2064  if (_smolar.pt() > _crit.smolar) {
2065  this->_phase = iphase_supercritical_gas;
2066  return;
2067  } else {
2068  this->_phase = iphase_supercritical_liquid;
2069  return;
2070  }
2071  }
2072  case iHmolar: {
2073  if (_hmolar.pt() > _crit.hmolar) {
2074  this->_phase = iphase_supercritical_gas;
2075  return;
2076  } else {
2077  this->_phase = iphase_supercritical_liquid;
2078  return;
2079  }
2080  }
2081  case iUmolar: {
2082  if (_umolar.pt() > _crit.umolar) {
2083  this->_phase = iphase_supercritical_gas;
2084  return;
2085  } else {
2086  this->_phase = iphase_supercritical_liquid;
2087  return;
2088  }
2089  }
2090  default: {
2091  throw ValueError("supercritical temp but other invalid for now");
2092  }
2093  }
2094  } else {
2095  throw ValueError(format("For now, we don't support T [%g K] below Ttriple [%g K]", _T, components[0].EOS().Ttriple));
2096  }
2097 }
2098 void get_dT_drho(HelmholtzEOSMixtureBackend* HEOS, parameters index, CoolPropDbl& dT, CoolPropDbl& drho) {
2099  CoolPropDbl T = HEOS->T(), rho = HEOS->rhomolar(), rhor = HEOS->get_reducing_state().rhomolar, Tr = HEOS->get_reducing_state().T,
2100  dT_dtau = -pow(T, 2) / Tr, R = HEOS->gas_constant(), delta = rho / rhor, tau = Tr / T;
2101 
2102  switch (index) {
2103  case iT:
2104  dT = 1;
2105  drho = 0;
2106  break;
2107  case iDmolar:
2108  dT = 0;
2109  drho = 1;
2110  break;
2111  case iDmass:
2112  dT = 0;
2113  drho = HEOS->molar_mass();
2114  break;
2115  case iP: {
2116  // dp/drho|T
2117  drho = R * T * (1 + 2 * delta * HEOS->dalphar_dDelta() + pow(delta, 2) * HEOS->d2alphar_dDelta2());
2118  // dp/dT|rho
2119  dT = rho * R * (1 + delta * HEOS->dalphar_dDelta() - tau * delta * HEOS->d2alphar_dDelta_dTau());
2120  break;
2121  }
2122  case iHmass:
2123  case iHmolar: {
2124  // dh/dT|rho
2125  dT = R
2126  * (-pow(tau, 2) * (HEOS->d2alpha0_dTau2() + HEOS->d2alphar_dTau2())
2127  + (1 + delta * HEOS->dalphar_dDelta() - tau * delta * HEOS->d2alphar_dDelta_dTau()));
2128  // dh/drhomolar|T
2129  drho =
2130  T * R / rho * (tau * delta * HEOS->d2alphar_dDelta_dTau() + delta * HEOS->dalphar_dDelta() + pow(delta, 2) * HEOS->d2alphar_dDelta2());
2131  if (index == iHmass) {
2132  // dhmolar/drhomolar|T * dhmass/dhmolar where dhmass/dhmolar = 1/mole mass
2133  drho /= HEOS->molar_mass();
2134  dT /= HEOS->molar_mass();
2135  }
2136  break;
2137  }
2138  case iSmass:
2139  case iSmolar: {
2140  // ds/dT|rho
2141  dT = R / T * (-pow(tau, 2) * (HEOS->d2alpha0_dTau2() + HEOS->d2alphar_dTau2()));
2142  // ds/drho|T
2143  drho = R / rho * (-(1 + delta * HEOS->dalphar_dDelta() - tau * delta * HEOS->d2alphar_dDelta_dTau()));
2144  if (index == iSmass) {
2145  // ds/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
2146  drho /= HEOS->molar_mass();
2147  dT /= HEOS->molar_mass();
2148  }
2149  break;
2150  }
2151  case iUmass:
2152  case iUmolar: {
2153  // du/dT|rho
2154  dT = R * (-pow(tau, 2) * (HEOS->d2alpha0_dTau2() + HEOS->d2alphar_dTau2()));
2155  // du/drho|T
2156  drho = HEOS->T() * R / rho * (tau * delta * HEOS->d2alphar_dDelta_dTau());
2157  if (index == iUmass) {
2158  // du/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
2159  drho /= HEOS->molar_mass();
2160  dT /= HEOS->molar_mass();
2161  }
2162  break;
2163  }
2164  case iTau:
2165  dT = 1 / dT_dtau;
2166  drho = 0;
2167  break;
2168  case iDelta:
2169  dT = 0;
2170  drho = 1 / rhor;
2171  break;
2172  default:
2173  throw ValueError(format("input to get_dT_drho[%s] is invalid", get_parameter_information(index, "short").c_str()));
2174  }
2175 }
2176 
2177 CoolPropDbl HelmholtzEOSMixtureBackend::calc_pressure_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2178  SimpleState reducing = calc_reducing_state_nocache(mole_fractions);
2179  CoolPropDbl delta = rhomolar / reducing.rhomolar;
2180  CoolPropDbl tau = reducing.T / T;
2181 
2182  // Calculate derivative
2183  int nTau = 0, nDelta = 1;
2184  CoolPropDbl dalphar_dDelta = calc_alphar_deriv_nocache(nTau, nDelta, mole_fractions, tau, delta);
2185 
2186  // Get pressure
2187  return rhomolar * gas_constant() * T * (1 + delta * dalphar_dDelta);
2188 }
2189 HelmholtzEOSBackend::StationaryPointReturnFlag HelmholtzEOSMixtureBackend::solver_dpdrho0_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomax,
2190  CoolPropDbl& light, CoolPropDbl& heavy) {
2191 
2193  class dpdrho_resid : public FuncWrapper1DWithTwoDerivs
2194  {
2195  public:
2196  HelmholtzEOSMixtureBackend* HEOS;
2197  CoolPropDbl T, p, delta, rhor, tau, R_u;
2198 
2199  dpdrho_resid(HelmholtzEOSMixtureBackend* HEOS, CoolPropDbl T, CoolPropDbl p)
2200  : HEOS(HEOS),
2201  T(T),
2202  p(p),
2203  delta(_HUGE),
2204  rhor(HEOS->get_reducing_state().rhomolar),
2205  tau(HEOS->get_reducing_state().T / T),
2206  R_u(HEOS->gas_constant()) {}
2207  double call(double rhomolar) {
2208  delta = rhomolar / rhor; // needed for derivative
2209  HEOS->update_DmolarT_direct(rhomolar, T);
2210  // dp/drho|T
2211  return R_u * T * (1 + 2 * delta * HEOS->dalphar_dDelta() + POW2(delta) * HEOS->d2alphar_dDelta2());
2212  };
2213  double deriv(double rhomolar) {
2214  // d2p/drho2|T
2215  return R_u * T / rhor * (2 * HEOS->dalphar_dDelta() + 4 * delta * HEOS->d2alphar_dDelta2() + POW2(delta) * HEOS->calc_d3alphar_dDelta3());
2216  };
2217  double second_deriv(double rhomolar) {
2218  // d3p/drho3|T
2219  return R_u * T / POW2(rhor)
2220  * (6 * HEOS->d2alphar_dDelta2() + 6 * delta * HEOS->d3alphar_dDelta3() + POW2(delta) * HEOS->calc_d4alphar_dDelta4());
2221  };
2222  };
2223  dpdrho_resid resid(this, T, p);
2224  light = -1;
2225  heavy = -1;
2226  try {
2227  light = Halley(resid, 1e-6, 1e-8, 100);
2228  double d2pdrho2__constT = resid.deriv(light);
2229  if (d2pdrho2__constT > 0) {
2230  // Not possible since curvature should be negative
2231  throw CoolProp::ValueError("curvature cannot be positive");
2232  }
2233  } catch (std::exception& e) {
2234  if (get_debug_level() > 5) {
2235  std::cout << e.what() << std::endl;
2236  };
2237  light = -1;
2238  }
2239 
2240  if (light < 0) {
2241  try {
2242  // Now we are going to do something VERY slow - increase density until curvature is positive
2243  double rho = 1e-6;
2244  for (std::size_t counter = 0; counter <= 100; counter++) {
2245  resid.call(rho); // Updates the state
2246  double curvature = resid.deriv(rho);
2247  if (curvature > 0) {
2248  light = rho;
2249  break;
2250  }
2251  rho *= 2;
2252  }
2253  } catch (...) {
2254  }
2255  }
2256 
2257  // First try a "normal" calculation of the stationary point on the liquid side
2258  for (double omega = 0.7; omega > 0; omega -= 0.2) {
2259  try {
2260  resid.options.add_number("omega", omega);
2261  heavy = Halley(resid, rhomax, 1e-8, 100);
2262  double d2pdrho2__constT = resid.deriv(heavy);
2263  if (d2pdrho2__constT < 0) {
2264  // Not possible since curvature should be positive
2265  throw CoolProp::ValueError("curvature cannot be negative");
2266  }
2267  break; // Jump out, we got a good solution
2268  } catch (std::exception& e) {
2269  if (get_debug_level() > 5) {
2270  std::cout << e.what() << std::endl;
2271  };
2272  heavy = -1;
2273  }
2274  }
2275 
2276  if (heavy < 0) {
2277  try {
2278  // Now we are going to do something VERY slow - decrease density until curvature is negative or pressure is negative
2279  double rho = rhomax;
2280  for (std::size_t counter = 0; counter <= 100; counter++) {
2281  resid.call(rho); // Updates the state
2282  double curvature = resid.deriv(rho);
2283  if (curvature < 0 || this->p() < 0) {
2284  heavy = rho;
2285  break;
2286  }
2287  rho /= 1.1;
2288  }
2289  } catch (...) {
2290  }
2291  }
2292 
2293  if (light > 0 && heavy > 0) {
2294  // Found two stationary points, done!
2295  return TWO_STATIONARY_POINTS_FOUND;
2296  }
2297  // If no solution is found for dpdrho|T=0 starting at high and low densities,
2298  // then try to do a bounded solver to see if you can find any solutions. If you
2299  // can't, p = f(rho) is probably monotonic (supercritical?), and the bounds are
2300  else if (light < 0 && heavy < 0) {
2301  double dpdrho_min = resid.call(1e-10);
2302  double dpdrho_max = resid.call(rhomax);
2303  if (dpdrho_max * dpdrho_min > 0) {
2304  return ZERO_STATIONARY_POINTS;
2305  } else {
2306  throw CoolProp::ValueError("zero stationary points -- does this make sense?");
2307  }
2308  } else {
2309  return ONE_STATIONARY_POINT_FOUND;
2310  }
2311 }
2312 // Define the residual to be driven to zero
2313 class SolverTPResid : public FuncWrapper1DWithThreeDerivs
2314 {
2315  public:
2316  HelmholtzEOSMixtureBackend* HEOS;
2317  CoolPropDbl T, p, delta, rhor, tau, R_u;
2318 
2319  SolverTPResid(HelmholtzEOSMixtureBackend* HEOS, CoolPropDbl T, CoolPropDbl p)
2320  : HEOS(HEOS),
2321  T(T),
2322  p(p),
2323  delta(_HUGE),
2324  rhor(HEOS->get_reducing_state().rhomolar),
2325  tau(HEOS->get_reducing_state().T / T),
2326  R_u(HEOS->gas_constant()) {}
2327  double call(double rhomolar) {
2328  delta = rhomolar / rhor; // needed for derivative
2329  HEOS->update_DmolarT_direct(rhomolar, T);
2330  CoolPropDbl peos = HEOS->p();
2331  return (peos - p) / p;
2332  };
2333  double deriv(double rhomolar) {
2334  // dp/drho|T / pspecified
2335  return R_u * T * (1 + 2 * delta * HEOS->dalphar_dDelta() + POW2(delta) * HEOS->d2alphar_dDelta2()) / p;
2336  };
2337  double second_deriv(double rhomolar) {
2338  // d2p/drho2|T / pspecified
2339  return R_u * T / rhor * (2 * HEOS->dalphar_dDelta() + 4 * delta * HEOS->d2alphar_dDelta2() + POW2(delta) * HEOS->calc_d3alphar_dDelta3()) / p;
2340  };
2341  double third_deriv(double rhomolar) {
2342  // d3p/drho3|T / pspecified
2343  return R_u * T / POW2(rhor)
2344  * (6 * HEOS->d2alphar_dDelta2() + 6 * delta * HEOS->d3alphar_dDelta3() + POW2(delta) * HEOS->calc_d4alphar_dDelta4()) / p;
2345  };
2346 };
2347 CoolPropDbl HelmholtzEOSMixtureBackend::SRK_covolume() {
2348  double b = 0;
2349  for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
2350  // Get the parameters for the cubic EOS
2351  CoolPropDbl Tc = get_fluid_constant(i, iT_critical), pc = get_fluid_constant(i, iP_critical);
2352  CoolPropDbl R = 8.3144598;
2353  b += mole_fractions[i] * 0.08664 * R * Tc / pc;
2354  }
2355  return b;
2356 }
2357 CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp_global(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomolar_max) {
2358  // Find the densities along the isotherm where dpdrho|T = 0 (if you can)
2359  CoolPropDbl light = -1, heavy = -1;
2360  StationaryPointReturnFlag retval = solver_dpdrho0_Tp(T, p, rhomolar_max, light, heavy);
2361 
2362  // Define the solver class
2363  SolverTPResid resid(this, T, p);
2364 
2365  if (retval == ZERO_STATIONARY_POINTS) {
2366  // It's monotonic (no stationary points found), so do the full bounded solver
2367  // for the density
2368  double rho = Brent(resid, 1e-10, rhomolar_max, DBL_EPSILON, 1e-8, 100);
2369  return rho;
2370  } else if (retval == TWO_STATIONARY_POINTS_FOUND) {
2371 
2372  // Calculate the pressures at the min and max densities where dpdrho|T = 0
2373  double p_at_rhomin_stationary = calc_pressure_nocache(T, light);
2374  double p_at_rhomax_stationary = calc_pressure_nocache(T, heavy);
2375 
2376  double rho_liq = -1, rho_vap = -1;
2377  if (p > p_at_rhomax_stationary) {
2378  int counter = 0;
2379  for (/* init above, for debugging */; counter <= 10; counter++) {
2380  // Bump up rhomax if needed to bound the given pressure
2381  double p_at_rhomax = calc_pressure_nocache(T, rhomolar_max);
2382  if (p_at_rhomax < p) {
2383  rhomolar_max *= 1.05;
2384  } else {
2385  break;
2386  }
2387  }
2388  // Look for liquid root starting at stationary point density
2389  rho_liq = Brent(resid, heavy, rhomolar_max, DBL_EPSILON, 1e-8, 100);
2390  }
2391 
2392  if (p < p_at_rhomin_stationary) {
2393  // Look for vapor root starting at stationary point density
2394  rho_vap = Brent(resid, light, 1e-10, DBL_EPSILON, 1e-8, 100);
2395  }
2396 
2397  if (rho_vap > 0 && rho_liq > 0) {
2398  // Both densities are the same
2399  if (std::abs(rho_vap - rho_liq) < 1e-10) {
2400  // return one of them
2401  return rho_vap;
2402  } else {
2403  // Two solutions found, keep the one with lower Gibbs energy
2404  double gibbsmolar_vap = calc_gibbsmolar_nocache(T, rho_vap);
2405  double gibbsmolar_liq = calc_gibbsmolar_nocache(T, rho_liq);
2406  if (gibbsmolar_liq < gibbsmolar_vap) {
2407  return rho_liq;
2408  } else {
2409  return rho_vap;
2410  }
2411  }
2412  } else if (rho_vap < 0 && rho_liq > 0) {
2413  // Liquid root found, return it
2414  return rho_liq;
2415  } else if (rho_vap > 0 && rho_liq < 0) {
2416  // Vapor root found, return it
2417  return rho_vap;
2418  } else {
2419  throw CoolProp::ValueError(format("No density solutions for T=%g,p=%g,z=%s", T, p, vec_to_string(static_cast<std::vector<double>>(mole_fractions), "%0.12g").c_str()));
2420  }
2421  } else {
2422  throw CoolProp::ValueError(
2423  format("One stationary point (not good) for T=%g,p=%g,z=%s", T, p, vec_to_string(static_cast<std::vector<double>>(mole_fractions), "%0.12g").c_str()));
2424  }
2425 };
2426 
2427 CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomolar_guess) {
2428  phases phase;
2429 
2430  SolverTPResid resid(this, T, p);
2431 
2432  // Check if the phase is imposed
2433  if (imposed_phase_index != iphase_not_imposed)
2434  // Use the imposed phase index
2435  phase = imposed_phase_index;
2436  else
2437  // Use the phase index in the class
2438  phase = _phase;
2439 
2440  if (rhomolar_guess < 0) // Not provided
2441  {
2442  // Calculate a guess value using SRK equation of state
2443  rhomolar_guess = solver_rho_Tp_SRK(T, p, phase);
2444 
2445  // A gas-like phase, ideal gas might not be the perfect model, but probably good enough
2446  if (phase == iphase_gas || phase == iphase_supercritical_gas || phase == iphase_supercritical) {
2447  if (rhomolar_guess < 0 || !ValidNumber(rhomolar_guess)) // If the guess is bad, probably high temperature, use ideal gas
2448  {
2449  rhomolar_guess = p / (gas_constant() * T);
2450  }
2451  } else if (phase == iphase_liquid) {
2452  double rhomolar;
2453  if (is_pure_or_pseudopure) {
2454  // It's liquid at subcritical pressure, we can use ancillaries as guess value
2455  CoolPropDbl _rhoLancval = static_cast<CoolPropDbl>(components[0].ancillaries.rhoL.evaluate(T));
2456  try {
2457  // First we try with Halley's method starting at saturated liquid
2458  rhomolar = Halley(resid, _rhoLancval, 1e-8, 100);
2459  if (!ValidNumber(rhomolar) || first_partial_deriv(iP, iDmolar, iT) < 0
2460  || second_partial_deriv(iP, iDmolar, iT, iDmolar, iT) < 0) {
2461  throw ValueError("Liquid density is invalid");
2462  }
2463  } catch (std::exception&) {
2464  // Next we try with a Brent method bounded solver since the function is 1-1
2465  rhomolar = Brent(resid, _rhoLancval * 0.9, _rhoLancval * 1.3, DBL_EPSILON, 1e-8, 100);
2466  if (!ValidNumber(rhomolar)) {
2467  throw ValueError();
2468  }
2469  }
2470  } else {
2471  // Try with 4th order Householder method starting at a very high density
2472  rhomolar = Householder4(&resid, 3 * rhomolar_reducing(), 1e-8, 100);
2473  }
2474  return rhomolar;
2475  } else if (phase == iphase_supercritical_liquid) {
2476  CoolPropDbl rhoLancval = static_cast<CoolPropDbl>(components[0].ancillaries.rhoL.evaluate(T));
2477  // Next we try with a Brent method bounded solver since the function is 1-1
2478  double rhomolar = Brent(resid, rhoLancval * 0.99, rhomolar_critical() * 4, DBL_EPSILON, 1e-8, 100);
2479  if (!ValidNumber(rhomolar)) {
2480  throw ValueError();
2481  }
2482  return rhomolar;
2483  }
2484  }
2485 
2486  try {
2487  double rhomolar = Householder4(resid, rhomolar_guess, 1e-8, 20);
2488  if (!ValidNumber(rhomolar) || rhomolar < 0) {
2489  throw ValueError();
2490  }
2491  if (phase == iphase_liquid) {
2492  double dpdrho = first_partial_deriv(iP, iDmolar, iT);
2493  double d2pdrho2 = second_partial_deriv(iP, iDmolar, iT, iDmolar, iT);
2494  if (dpdrho < 0 || d2pdrho2 < 0) {
2495  // Try again with a larger density in order to end up at the right solution
2496  rhomolar = Householder4(resid, 3 * rhomolar_reducing(), 1e-8, 100);
2497  return rhomolar;
2498  }
2499  } else if (phase == iphase_gas) {
2500  double dpdrho = first_partial_deriv(iP, iDmolar, iT);
2501  double d2pdrho2 = second_partial_deriv(iP, iDmolar, iT, iDmolar, iT);
2502  if (dpdrho < 0 || d2pdrho2 > 0) {
2503  // Try again with a tiny density in order to end up at the right solution
2504  rhomolar = Householder4(resid, 1e-6, 1e-8, 100);
2505  return rhomolar;
2506  }
2507  }
2508  return rhomolar;
2509  } catch (std::exception& e) {
2510  if (phase == iphase_supercritical || phase == iphase_supercritical_gas) {
2511  double rhomolar = Brent(resid, 1e-10, 3 * rhomolar_reducing(), DBL_EPSILON, 1e-8, 100);
2512  return rhomolar;
2513  } else if (is_pure_or_pseudopure && T > T_critical()) {
2514  try {
2515  double rhomolar = Brent(resid, 1e-10, 5 * rhomolar_reducing(), DBL_EPSILON, 1e-8, 100);
2516  return rhomolar;
2517 
2518  } catch (...) {
2519  double rhomolar = Householder4(resid, 3 * rhomolar_reducing(), 1e-8, 100);
2520  return rhomolar;
2521  }
2522  }
2523  throw ValueError(format("solver_rho_Tp was unable to find a solution for T=%10Lg, p=%10Lg, with guess value %10Lg with error: %s", T, p,
2524  rhomolar_guess, e.what()));
2525  }
2526 }
2527 CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp_SRK(CoolPropDbl T, CoolPropDbl p, phases phase) {
2528  CoolPropDbl rhomolar, R_u = gas_constant(), a = 0, b = 0, k_ij = 0;
2529 
2530  for (std::size_t i = 0; i < components.size(); ++i) {
2531  CoolPropDbl Tci = components[i].EOS().reduce.T, pci = components[i].EOS().reduce.p, acentric_i = components[i].EOS().acentric;
2532  CoolPropDbl m_i = 0.480 + 1.574 * acentric_i - 0.176 * pow(acentric_i, 2);
2533  CoolPropDbl b_i = 0.08664 * R_u * Tci / pci;
2534  b += mole_fractions[i] * b_i;
2535 
2536  CoolPropDbl a_i = 0.42747 * pow(R_u * Tci, 2) / pci * pow(1 + m_i * (1 - sqrt(T / Tci)), 2);
2537 
2538  for (std::size_t j = 0; j < components.size(); ++j) {
2539  CoolPropDbl Tcj = components[j].EOS().reduce.T, pcj = components[j].EOS().reduce.p, acentric_j = components[j].EOS().acentric;
2540  CoolPropDbl m_j = 0.480 + 1.574 * acentric_j - 0.176 * pow(acentric_j, 2);
2541 
2542  CoolPropDbl a_j = 0.42747 * pow(R_u * Tcj, 2) / pcj * pow(1 + m_j * (1 - sqrt(T / Tcj)), 2);
2543 
2544  k_ij = 0;
2545  //if (i == j){
2546  // k_ij = 0;
2547  //}
2548  //else{
2549  // k_ij = 0;
2550  //}
2551 
2552  a += mole_fractions[i] * mole_fractions[j] * sqrt(a_i * a_j) * (1 - k_ij);
2553  }
2554  }
2555 
2556  CoolPropDbl A = a * p / pow(R_u * T, 2);
2557  CoolPropDbl B = b * p / (R_u * T);
2558 
2559  //Solve the cubic for solutions for Z = p/(rho*R*T)
2560  double Z0, Z1, Z2;
2561  int Nsolns;
2562  solve_cubic(1, -1, A - B - B * B, -A * B, Nsolns, Z0, Z1, Z2);
2563 
2564  // Determine the guess value
2565  if (Nsolns == 1) {
2566  rhomolar = p / (Z0 * R_u * T);
2567  } else {
2568  CoolPropDbl rhomolar0 = p / (Z0 * R_u * T);
2569  CoolPropDbl rhomolar1 = p / (Z1 * R_u * T);
2570  CoolPropDbl rhomolar2 = p / (Z2 * R_u * T);
2571 
2572  // Check if only one solution is positive, return the solution if that is the case
2573  if (rhomolar0 > 0 && rhomolar1 <= 0 && rhomolar2 <= 0) {
2574  return rhomolar0;
2575  }
2576  if (rhomolar0 <= 0 && rhomolar1 > 0 && rhomolar2 <= 0) {
2577  return rhomolar1;
2578  }
2579  if (rhomolar0 <= 0 && rhomolar1 <= 0 && rhomolar2 > 0) {
2580  return rhomolar2;
2581  }
2582 
2583  switch (phase) {
2584  case iphase_liquid:
2585  case iphase_supercritical_liquid:
2586  rhomolar = max3(rhomolar0, rhomolar1, rhomolar2);
2587  break;
2588  case iphase_gas:
2589  case iphase_supercritical_gas:
2590  case iphase_supercritical:
2591  rhomolar = min3(rhomolar0, rhomolar1, rhomolar2);
2592  break;
2593  default:
2594  throw ValueError("Bad phase to solver_rho_Tp_SRK");
2595  };
2596  }
2597  return rhomolar;
2598 }
2599 
2600 CoolPropDbl HelmholtzEOSMixtureBackend::calc_pressure(void) {
2601  // Calculate the reducing parameters
2602  _delta = _rhomolar / _reducing.rhomolar;
2603  _tau = _reducing.T / _T;
2604 
2605  // Calculate derivative if needed
2606  CoolPropDbl dar_dDelta = dalphar_dDelta();
2607  CoolPropDbl R_u = gas_constant();
2608 
2609  // Get pressure
2610  _p = _rhomolar * R_u * _T * (1 + _delta.pt() * dar_dDelta);
2611 
2612  //std::cout << format("p: %13.12f %13.12f %10.9f %10.9f %10.9f %10.9f %g\n",_T,_rhomolar,_tau,_delta,mole_fractions[0],dar_dDelta,_p);
2613  //if (_p < 0){
2614  // throw ValueError("Pressure is less than zero");
2615  //}
2616 
2617  return static_cast<CoolPropDbl>(_p);
2618 }
2619 CoolPropDbl HelmholtzEOSMixtureBackend::calc_hmolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2620  // Calculate the reducing parameters
2621  CoolPropDbl delta = rhomolar / _reducing.rhomolar;
2622  CoolPropDbl tau = _reducing.T / T;
2623 
2624  // Calculate derivatives if needed, or just use cached values
2625  // Calculate derivative if needed
2626  CoolPropDbl dar_dDelta = calc_alphar_deriv_nocache(0, 1, mole_fractions, tau, delta);
2627  CoolPropDbl dar_dTau = calc_alphar_deriv_nocache(1, 0, mole_fractions, tau, delta);
2628  CoolPropDbl da0_dTau = calc_alpha0_deriv_nocache(1, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2629  CoolPropDbl R_u = gas_constant();
2630 
2631  // Get molar enthalpy
2632  return R_u * T * (1 + tau * (da0_dTau + dar_dTau) + delta * dar_dDelta);
2633 }
2634 CoolPropDbl HelmholtzEOSMixtureBackend::calc_hmolar(void) {
2635  if (get_debug_level() >= 50)
2636  std::cout << format("HelmholtzEOSMixtureBackend::calc_hmolar: 2phase: %d T: %g rhomomolar: %g", isTwoPhase(), _T, _rhomolar) << std::endl;
2637  if (isTwoPhase()) {
2638  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
2639  if (std::abs(_Q) < DBL_EPSILON) {
2640  _hmolar = SatL->hmolar();
2641  } else if (std::abs(_Q - 1) < DBL_EPSILON) {
2642  _hmolar = SatV->hmolar();
2643  } else {
2644  _hmolar = _Q * SatV->hmolar() + (1 - _Q) * SatL->hmolar();
2645  }
2646  return static_cast<CoolPropDbl>(_hmolar);
2647  } else if (isHomogeneousPhase()) {
2648  // Calculate the reducing parameters
2649  _delta = _rhomolar / _reducing.rhomolar;
2650  _tau = _reducing.T / _T;
2651 
2652  // Calculate derivatives if needed, or just use cached values
2653  CoolPropDbl da0_dTau = dalpha0_dTau();
2654  CoolPropDbl dar_dTau = dalphar_dTau();
2655  CoolPropDbl dar_dDelta = dalphar_dDelta();
2656  CoolPropDbl R_u = gas_constant();
2657 
2658  // Get molar enthalpy
2659  _hmolar = R_u * _T * (1 + _tau.pt() * (da0_dTau + dar_dTau) + _delta.pt() * dar_dDelta);
2660 
2661  return static_cast<CoolPropDbl>(_hmolar);
2662  } else {
2663  throw ValueError(format("phase is invalid in calc_hmolar"));
2664  }
2665 }
2666 CoolPropDbl HelmholtzEOSMixtureBackend::calc_smolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2667  // Calculate the reducing parameters
2668  CoolPropDbl delta = rhomolar / _reducing.rhomolar;
2669  CoolPropDbl tau = _reducing.T / T;
2670 
2671  // Calculate derivatives if needed, or just use cached values
2672  // Calculate derivative if needed
2673  CoolPropDbl dar_dTau = calc_alphar_deriv_nocache(1, 0, mole_fractions, tau, delta);
2674  CoolPropDbl ar = calc_alphar_deriv_nocache(0, 0, mole_fractions, tau, delta);
2675  CoolPropDbl da0_dTau = calc_alpha0_deriv_nocache(1, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2676  CoolPropDbl a0 = calc_alpha0_deriv_nocache(0, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2677  CoolPropDbl R_u = gas_constant();
2678 
2679  // Get molar entropy
2680  return R_u * (tau * (da0_dTau + dar_dTau) - a0 - ar);
2681 }
2682 CoolPropDbl HelmholtzEOSMixtureBackend::calc_smolar(void) {
2683  if (isTwoPhase()) {
2684  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
2685  if (std::abs(_Q) < DBL_EPSILON) {
2686  _smolar = SatL->smolar();
2687  } else if (std::abs(_Q - 1) < DBL_EPSILON) {
2688  _smolar = SatV->smolar();
2689  } else {
2690  _smolar = _Q * SatV->smolar() + (1 - _Q) * SatL->smolar();
2691  }
2692  return static_cast<CoolPropDbl>(_smolar);
2693  } else if (isHomogeneousPhase()) {
2694  // Calculate the reducing parameters
2695  _delta = _rhomolar / _reducing.rhomolar;
2696  _tau = _reducing.T / _T;
2697 
2698  // Calculate derivatives if needed, or just use cached values
2699  CoolPropDbl da0_dTau = dalpha0_dTau();
2700  CoolPropDbl ar = alphar();
2701  CoolPropDbl a0 = alpha0();
2702  CoolPropDbl dar_dTau = dalphar_dTau();
2703  CoolPropDbl R_u = gas_constant();
2704 
2705  // Get molar entropy
2706  _smolar = R_u * (_tau.pt() * (da0_dTau + dar_dTau) - a0 - ar);
2707 
2708  return static_cast<CoolPropDbl>(_smolar);
2709  } else {
2710  throw ValueError(format("phase is invalid in calc_smolar"));
2711  }
2712 }
2713 CoolPropDbl HelmholtzEOSMixtureBackend::calc_umolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2714  // Calculate the reducing parameters
2715  CoolPropDbl delta = rhomolar / _reducing.rhomolar;
2716  CoolPropDbl tau = _reducing.T / T;
2717 
2718  // Calculate derivatives
2719  CoolPropDbl dar_dTau = calc_alphar_deriv_nocache(1, 0, mole_fractions, tau, delta);
2720  CoolPropDbl da0_dTau = calc_alpha0_deriv_nocache(1, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2721  CoolPropDbl R_u = gas_constant();
2722 
2723  // Get molar internal energy
2724  return R_u * T * tau * (da0_dTau + dar_dTau);
2725 }
2726 CoolPropDbl HelmholtzEOSMixtureBackend::calc_umolar(void) {
2727  if (isTwoPhase()) {
2728  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
2729  if (std::abs(_Q) < DBL_EPSILON) {
2730  _umolar = SatL->umolar();
2731  } else if (std::abs(_Q - 1) < DBL_EPSILON) {
2732  _umolar = SatV->umolar();
2733  } else {
2734  _umolar = _Q * SatV->umolar() + (1 - _Q) * SatL->umolar();
2735  }
2736  return static_cast<CoolPropDbl>(_umolar);
2737  } else if (isHomogeneousPhase()) {
2738  // Calculate the reducing parameters
2739  _delta = _rhomolar / _reducing.rhomolar;
2740  _tau = _reducing.T / _T;
2741 
2742  // Calculate derivatives if needed, or just use cached values
2743  CoolPropDbl da0_dTau = dalpha0_dTau();
2744  CoolPropDbl dar_dTau = dalphar_dTau();
2745  CoolPropDbl R_u = gas_constant();
2746 
2747  // Get molar internal energy
2748  _umolar = R_u * _T * _tau.pt() * (da0_dTau + dar_dTau);
2749 
2750  return static_cast<CoolPropDbl>(_umolar);
2751  } else {
2752  throw ValueError(format("phase is invalid in calc_umolar"));
2753  }
2754 }
2755 CoolPropDbl HelmholtzEOSMixtureBackend::calc_cvmolar(void) {
2756  // Calculate the reducing parameters
2757  _delta = _rhomolar / _reducing.rhomolar;
2758  _tau = _reducing.T / _T;
2759 
2760  // Calculate derivatives if needed, or just use cached values
2761  CoolPropDbl d2ar_dTau2 = d2alphar_dTau2();
2762  CoolPropDbl d2a0_dTau2 = d2alpha0_dTau2();
2763  CoolPropDbl R_u = gas_constant();
2764 
2765  // Get cv
2766  _cvmolar = -R_u * pow(_tau.pt(), 2) * (d2ar_dTau2 + d2a0_dTau2);
2767 
2768  return static_cast<double>(_cvmolar);
2769 }
2770 CoolPropDbl HelmholtzEOSMixtureBackend::calc_cpmolar(void) {
2771  // Calculate the reducing parameters
2772  _delta = _rhomolar / _reducing.rhomolar;
2773  _tau = _reducing.T / _T;
2774 
2775  // Calculate derivatives if needed, or just use cached values
2776  CoolPropDbl d2a0_dTau2 = d2alpha0_dTau2();
2777  CoolPropDbl dar_dDelta = dalphar_dDelta();
2778  CoolPropDbl d2ar_dDelta2 = d2alphar_dDelta2();
2779  CoolPropDbl d2ar_dDelta_dTau = d2alphar_dDelta_dTau();
2780  CoolPropDbl d2ar_dTau2 = d2alphar_dTau2();
2781  CoolPropDbl R_u = gas_constant();
2782 
2783  // Get cp
2784  _cpmolar = R_u
2785  * (-pow(_tau.pt(), 2) * (d2ar_dTau2 + d2a0_dTau2)
2786  + pow(1 + _delta.pt() * dar_dDelta - _delta.pt() * _tau.pt() * d2ar_dDelta_dTau, 2)
2787  / (1 + 2 * _delta.pt() * dar_dDelta + pow(_delta.pt(), 2) * d2ar_dDelta2));
2788 
2789  return static_cast<double>(_cpmolar);
2790 }
2791 CoolPropDbl HelmholtzEOSMixtureBackend::calc_cpmolar_idealgas(void) {
2792  // Calculate the reducing parameters
2793  _delta = _rhomolar / _reducing.rhomolar;
2794  _tau = _reducing.T / _T;
2795 
2796  // Calculate derivatives if needed, or just use cached values
2797  CoolPropDbl d2a0_dTau2 = d2alpha0_dTau2();
2798  CoolPropDbl R_u = gas_constant();
2799 
2800  // Get cp of the ideal gas
2801  return R_u * (1 + (-pow(_tau.pt(), 2)) * d2a0_dTau2);
2802 }
2803 CoolPropDbl HelmholtzEOSMixtureBackend::calc_speed_sound(void) {
2804  if (isTwoPhase()) {
2805  if (std::abs(_Q) < DBL_EPSILON) {
2806  return SatL->speed_sound();
2807  } else if (std::abs(_Q - 1) < DBL_EPSILON) {
2808  return SatV->speed_sound();
2809  } else {
2810  throw ValueError(format("Speed of sound is not defined for two-phase states because it depends on the distribution of phases."));
2811  }
2812  } else if (isHomogeneousPhase()) {
2813  // Calculate the reducing parameters
2814  _delta = _rhomolar / _reducing.rhomolar;
2815  _tau = _reducing.T / _T;
2816 
2817  // Calculate derivatives if needed, or just use cached values
2818  CoolPropDbl d2a0_dTau2 = d2alpha0_dTau2();
2819  CoolPropDbl dar_dDelta = dalphar_dDelta();
2820  CoolPropDbl d2ar_dDelta2 = d2alphar_dDelta2();
2821  CoolPropDbl d2ar_dDelta_dTau = d2alphar_dDelta_dTau();
2822  CoolPropDbl d2ar_dTau2 = d2alphar_dTau2();
2823  CoolPropDbl R_u = gas_constant();
2824  CoolPropDbl mm = molar_mass();
2825 
2826  // Get speed of sound
2827  _speed_sound = sqrt(
2828  R_u * _T / mm
2829  * (1 + 2 * _delta.pt() * dar_dDelta + pow(_delta.pt(), 2) * d2ar_dDelta2
2830  - pow(1 + _delta.pt() * dar_dDelta - _delta.pt() * _tau.pt() * d2ar_dDelta_dTau, 2) / (pow(_tau.pt(), 2) * (d2ar_dTau2 + d2a0_dTau2))));
2831 
2832  return static_cast<CoolPropDbl>(_speed_sound);
2833  } else {
2834  throw ValueError(format("phase is invalid in calc_speed_sound"));
2835  }
2836 }
2837 
2838 CoolPropDbl HelmholtzEOSMixtureBackend::calc_gibbsmolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2839  // Calculate the reducing parameters
2840  CoolPropDbl delta = rhomolar / _reducing.rhomolar;
2841  CoolPropDbl tau = _reducing.T / T;
2842 
2843  // Calculate derivatives
2844  CoolPropDbl dar_dDelta = calc_alphar_deriv_nocache(0, 1, mole_fractions, tau, delta);
2845  CoolPropDbl ar = calc_alphar_deriv_nocache(0, 0, mole_fractions, tau, delta);
2846  CoolPropDbl a0 = calc_alpha0_deriv_nocache(0, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2847 
2848  // Get molar gibbs function
2849  return gas_constant() * T * (1 + a0 + ar + delta * dar_dDelta);
2850 }
2851 CoolPropDbl HelmholtzEOSMixtureBackend::calc_gibbsmolar(void) {
2852  if (isTwoPhase()) {
2853  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
2854  _gibbsmolar = _Q * SatV->gibbsmolar() + (1 - _Q) * SatL->gibbsmolar();
2855  return static_cast<CoolPropDbl>(_gibbsmolar);
2856  } else if (isHomogeneousPhase()) {
2857  // Calculate the reducing parameters
2858  _delta = _rhomolar / _reducing.rhomolar;
2859  _tau = _reducing.T / _T;
2860 
2861  // Calculate derivatives if needed, or just use cached values
2862  CoolPropDbl ar = alphar();
2863  CoolPropDbl a0 = alpha0();
2864  CoolPropDbl dar_dDelta = dalphar_dDelta();
2865  CoolPropDbl R_u = gas_constant();
2866 
2867  // Get molar gibbs function
2868  _gibbsmolar = R_u * _T * (1 + a0 + ar + _delta.pt() * dar_dDelta);
2869 
2870  return static_cast<CoolPropDbl>(_gibbsmolar);
2871  } else {
2872  throw ValueError(format("phase is invalid in calc_gibbsmolar"));
2873  }
2874 }
2875 void HelmholtzEOSMixtureBackend::calc_excess_properties(void) {
2876  _gibbsmolar_excess = this->gibbsmolar(), _smolar_excess = this->smolar(), _hmolar_excess = this->hmolar();
2877  _umolar_excess = this->umolar();
2878  _volumemolar_excess = 1 / this->rhomolar();
2879  for (std::size_t i = 0; i < components.size(); ++i) {
2880  transient_pure_state.reset(new HelmholtzEOSBackend(components[i].name));
2881  transient_pure_state->update(PT_INPUTS, p(), T());
2882  double x_i = mole_fractions[i];
2883  double R = gas_constant();
2884  _gibbsmolar_excess = static_cast<double>(_gibbsmolar_excess) - x_i * (transient_pure_state->gibbsmolar() + R * T() * log(x_i));
2885  _hmolar_excess = static_cast<double>(_hmolar_excess) - x_i * transient_pure_state->hmolar();
2886  _umolar_excess = static_cast<double>(_umolar_excess) - x_i * transient_pure_state->umolar();
2887  _smolar_excess = static_cast<double>(_smolar_excess) - x_i * (transient_pure_state->smolar() - R * log(x_i));
2888  _volumemolar_excess = static_cast<double>(_volumemolar_excess) - x_i / transient_pure_state->rhomolar();
2889  }
2890  _helmholtzmolar_excess = static_cast<double>(_umolar_excess) - _T * static_cast<double>(_smolar_excess);
2891 }
2892 
2893 CoolPropDbl HelmholtzEOSMixtureBackend::calc_helmholtzmolar(void) {
2894  if (isTwoPhase()) {
2895  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
2896  _helmholtzmolar = _Q * SatV->helmholtzmolar() + (1 - _Q) * SatL->helmholtzmolar();
2897  return static_cast<CoolPropDbl>(_helmholtzmolar);
2898  } else if (isHomogeneousPhase()) {
2899  // Calculate the reducing parameters
2900  _delta = _rhomolar / _reducing.rhomolar;
2901  _tau = _reducing.T / _T;
2902 
2903  // Calculate derivatives if needed, or just use cached values
2904  CoolPropDbl ar = alphar();
2905  CoolPropDbl a0 = alpha0();
2906  CoolPropDbl R_u = gas_constant();
2907 
2908  // Get molar Helmholtz energy
2909  _helmholtzmolar = R_u * _T * (a0 + ar);
2910 
2911  return static_cast<CoolPropDbl>(_helmholtzmolar);
2912  } else {
2913  throw ValueError(format("phase is invalid in calc_helmholtzmolar"));
2914  }
2915 }
2916 CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(std::size_t i) {
2917  x_N_dependency_flag xN_flag = XN_DEPENDENT;
2918  return exp(MixtureDerivatives::ln_fugacity_coefficient(*this, i, xN_flag));
2919 }
2920 CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity(std::size_t i) {
2921  x_N_dependency_flag xN_flag = XN_DEPENDENT;
2922  return MixtureDerivatives::fugacity_i(*this, i, xN_flag);
2923 }
2924 CoolPropDbl HelmholtzEOSMixtureBackend::calc_chemical_potential(std::size_t i) {
2925  x_N_dependency_flag xN_flag = XN_DEPENDENT;
2926  double Tci = get_fluid_constant(i, iT_critical);
2927  double rhoci = get_fluid_constant(i, irhomolar_critical);
2928  double dnar_dni__const_T_V_nj = MixtureDerivatives::dnalphar_dni__constT_V_nj(*this, i, xN_flag);
2929  double dna0_dni__const_T_V_nj =
2930  components[i].EOS().alpha0.base(tau() * (Tci / T_reducing()), delta() / (rhoci / rhomolar_reducing())) + 1 + log(mole_fractions[i]);
2931  return gas_constant() * T() * (dna0_dni__const_T_V_nj + dnar_dni__const_T_V_nj);
2932 }
2933 CoolPropDbl HelmholtzEOSMixtureBackend::calc_phase_identification_parameter(void) {
2934  return 2
2935  - rhomolar()
2936  * (second_partial_deriv(iP, iDmolar, iT, iT, iDmolar) / first_partial_deriv(iP, iT, iDmolar)
2937  - second_partial_deriv(iP, iDmolar, iT, iDmolar, iT) / first_partial_deriv(iP, iDmolar, iT));
2938 }
2939 
2940 SimpleState HelmholtzEOSMixtureBackend::calc_reducing_state_nocache(const std::vector<CoolPropDbl>& mole_fractions) {
2941  SimpleState reducing;
2942  if (is_pure_or_pseudopure) {
2943  reducing = components[0].EOS().reduce;
2944  } else {
2945  reducing.T = Reducing->Tr(mole_fractions);
2946  reducing.rhomolar = Reducing->rhormolar(mole_fractions);
2947  }
2948  return reducing;
2949 }
2950 void HelmholtzEOSMixtureBackend::calc_reducing_state(void) {
2951  if (get_mole_fractions_ref().empty()) {
2952  throw ValueError("Mole fractions must be set before calling calc_reducing_state");
2953  }
2955  _reducing = calc_reducing_state_nocache(mole_fractions);
2956  _crit = _reducing;
2957 }
2958 void HelmholtzEOSMixtureBackend::calc_all_alphar_deriv_cache(const std::vector<CoolPropDbl>& mole_fractions, const CoolPropDbl& tau,
2959  const CoolPropDbl& delta) {
2960  deriv_counter++;
2961  bool cache_values = true;
2962  HelmholtzDerivatives derivs = residual_helmholtz->all(*this, get_mole_fractions_ref(), tau, delta, cache_values);
2963  _alphar = derivs.alphar;
2964  _dalphar_dDelta = derivs.dalphar_ddelta;
2965  _dalphar_dTau = derivs.dalphar_dtau;
2966  _d2alphar_dDelta2 = derivs.d2alphar_ddelta2;
2967  _d2alphar_dDelta_dTau = derivs.d2alphar_ddelta_dtau;
2968  _d2alphar_dTau2 = derivs.d2alphar_dtau2;
2969  _d3alphar_dDelta3 = derivs.d3alphar_ddelta3;
2970  _d3alphar_dDelta2_dTau = derivs.d3alphar_ddelta2_dtau;
2971  _d3alphar_dDelta_dTau2 = derivs.d3alphar_ddelta_dtau2;
2972  _d3alphar_dTau3 = derivs.d3alphar_dtau3;
2973  _d4alphar_dDelta4 = derivs.d4alphar_ddelta4;
2974  _d4alphar_dDelta3_dTau = derivs.d4alphar_ddelta3_dtau;
2975  _d4alphar_dDelta2_dTau2 = derivs.d4alphar_ddelta2_dtau2;
2976  _d4alphar_dDelta_dTau3 = derivs.d4alphar_ddelta_dtau3;
2977  _d4alphar_dTau4 = derivs.d4alphar_dtau4;
2978 }
2979 
2980 CoolPropDbl HelmholtzEOSMixtureBackend::calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl>& mole_fractions,
2981  const CoolPropDbl& tau, const CoolPropDbl& delta) {
2982  bool cache_values = false;
2983  HelmholtzDerivatives derivs = residual_helmholtz->all(*this, mole_fractions, tau, delta, cache_values);
2984  return derivs.get(nTau, nDelta);
2985 }
2986 CoolPropDbl HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl>& mole_fractions,
2987  const CoolPropDbl& tau, const CoolPropDbl& delta, const CoolPropDbl& Tr,
2988  const CoolPropDbl& rhor) {
2989  CoolPropDbl val;
2990  if (components.size() == 0) {
2991  throw ValueError("No alpha0 derivatives are available");
2992  }
2993  if (is_pure_or_pseudopure) {
2994  EquationOfState& E = components[0].EOS();
2995  // In the case of cubics, we need to use the shifted tau^*=Tc/T and delta^*=rho/rhoc
2996  // rather than tau=Tr/T and delta=rho/rhor
2997  // For multiparameter EOS, this changes nothing because Tc/Tr = 1 and rhoc/rhor = 1
2998  double Tc = get_fluid_constant(0, iT_reducing), rhomolarc = get_fluid_constant(0, irhomolar_reducing);
2999 
3000  // Cache the reducing temperature in some terms that need it (GERG-2004 models)
3001  E.alpha0.set_Tred(Tc);
3002  double taustar = Tc / Tr * tau, deltastar = rhor / rhomolarc * delta;
3003  if (nTau == 0 && nDelta == 0) {
3004  val = E.base0(taustar, deltastar);
3005  } else if (nTau == 0 && nDelta == 1) {
3006  val = E.dalpha0_dDelta(taustar, deltastar);
3007  } else if (nTau == 1 && nDelta == 0) {
3008  val = E.dalpha0_dTau(taustar, deltastar);
3009  } else if (nTau == 0 && nDelta == 2) {
3010  val = E.d2alpha0_dDelta2(taustar, deltastar);
3011  } else if (nTau == 1 && nDelta == 1) {
3012  val = E.d2alpha0_dDelta_dTau(taustar, deltastar);
3013  } else if (nTau == 2 && nDelta == 0) {
3014  val = E.d2alpha0_dTau2(taustar, deltastar);
3015  } else if (nTau == 0 && nDelta == 3) {
3016  val = E.d3alpha0_dDelta3(taustar, deltastar);
3017  } else if (nTau == 1 && nDelta == 2) {
3018  val = E.d3alpha0_dDelta2_dTau(taustar, deltastar);
3019  } else if (nTau == 2 && nDelta == 1) {
3020  val = E.d3alpha0_dDelta_dTau2(taustar, deltastar);
3021  } else if (nTau == 3 && nDelta == 0) {
3022  val = E.d3alpha0_dTau3(taustar, deltastar);
3023  } else {
3024  throw ValueError();
3025  }
3026  val *= pow(rhor / rhomolarc, nDelta);
3027  val /= pow(Tr / Tc, nTau);
3028  if (!ValidNumber(val)) {
3029  //calc_alpha0_deriv_nocache(nTau,nDelta,mole_fractions,tau,delta,Tr,rhor);
3030  throw ValueError(format("calc_alpha0_deriv_nocache returned invalid number with inputs nTau: %d, nDelta: %d, tau: %Lg, delta: %Lg", nTau,
3031  nDelta, tau, delta));
3032  } else {
3033  return val;
3034  }
3035  } else {
3036  // See Table B5, GERG 2008 from Kunz Wagner, JCED, 2012
3037  std::size_t N = mole_fractions.size();
3038  CoolPropDbl summer = 0;
3039  CoolPropDbl tau_i, delta_i, rho_ci, T_ci;
3040  CoolPropDbl Rmix = gas_constant();
3041  for (unsigned int i = 0; i < N; ++i) {
3042 
3043  rho_ci = get_fluid_constant(i, irhomolar_critical);
3044  T_ci = get_fluid_constant(i, iT_critical);
3045  CoolPropDbl Rcomponent = get_fluid_constant(i, igas_constant);
3046  tau_i = T_ci * tau / Tr;
3047  delta_i = delta * rhor / rho_ci;
3048  CoolPropDbl Rratio = Rcomponent / Rmix;
3049 
3050  // Cache the reducing temperature in some terms that need it (GERG-2004 models)
3051  components[i].EOS().alpha0.set_Tred(Tr);
3052 
3053  if (nTau == 0 && nDelta == 0) {
3054  double logxi = (std::abs(mole_fractions[i]) > DBL_EPSILON) ? log(mole_fractions[i]) : 0;
3055  summer += mole_fractions[i] * Rratio * (components[i].EOS().base0(tau_i, delta_i) + logxi);
3056  } else if (nTau == 0 && nDelta == 1) {
3057  summer += mole_fractions[i] * Rratio * rhor / rho_ci * components[i].EOS().dalpha0_dDelta(tau_i, delta_i);
3058  } else if (nTau == 1 && nDelta == 0) {
3059  summer += mole_fractions[i] * Rratio * T_ci / Tr * components[i].EOS().dalpha0_dTau(tau_i, delta_i);
3060  } else if (nTau == 0 && nDelta == 2) {
3061  summer += mole_fractions[i] * Rratio * pow(rhor / rho_ci, 2) * components[i].EOS().d2alpha0_dDelta2(tau_i, delta_i);
3062  } else if (nTau == 1 && nDelta == 1) {
3063  summer += mole_fractions[i] * Rratio * rhor / rho_ci * T_ci / Tr * components[i].EOS().d2alpha0_dDelta_dTau(tau_i, delta_i);
3064  } else if (nTau == 2 && nDelta == 0) {
3065  summer += mole_fractions[i] * Rratio * pow(T_ci / Tr, 2) * components[i].EOS().d2alpha0_dTau2(tau_i, delta_i);
3066  } else {
3067  throw ValueError();
3068  }
3069  }
3070  return summer;
3071  }
3072 }
3073 CoolPropDbl HelmholtzEOSMixtureBackend::calc_alphar(void) {
3074  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3075  return static_cast<CoolPropDbl>(_alphar);
3076 }
3077 CoolPropDbl HelmholtzEOSMixtureBackend::calc_dalphar_dDelta(void) {
3078  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3079  return static_cast<CoolPropDbl>(_dalphar_dDelta);
3080 }
3081 CoolPropDbl HelmholtzEOSMixtureBackend::calc_dalphar_dTau(void) {
3082  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3083  return static_cast<CoolPropDbl>(_dalphar_dTau);
3084 }
3085 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alphar_dTau2(void) {
3086  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3087  return static_cast<CoolPropDbl>(_d2alphar_dTau2);
3088 }
3089 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alphar_dDelta_dTau(void) {
3090  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3091  return static_cast<CoolPropDbl>(_d2alphar_dDelta_dTau);
3092 }
3093 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alphar_dDelta2(void) {
3094  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3095  return static_cast<CoolPropDbl>(_d2alphar_dDelta2);
3096 }
3097 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alphar_dDelta3(void) {
3098  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3099  return static_cast<CoolPropDbl>(_d3alphar_dDelta3);
3100 }
3101 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alphar_dDelta2_dTau(void) {
3102  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3103  return static_cast<CoolPropDbl>(_d3alphar_dDelta2_dTau);
3104 }
3105 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alphar_dDelta_dTau2(void) {
3106  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3107  return static_cast<CoolPropDbl>(_d3alphar_dDelta_dTau2);
3108 }
3109 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alphar_dTau3(void) {
3110  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3111  return static_cast<CoolPropDbl>(_d3alphar_dTau3);
3112 }
3113 
3114 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dDelta4(void) {
3115  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3116  return static_cast<CoolPropDbl>(_d4alphar_dDelta4);
3117 }
3118 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dDelta3_dTau(void) {
3119  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3120  return static_cast<CoolPropDbl>(_d4alphar_dDelta3_dTau);
3121 }
3122 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dDelta2_dTau2(void) {
3123  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3124  return static_cast<CoolPropDbl>(_d4alphar_dDelta2_dTau2);
3125 }
3126 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dDelta_dTau3(void) {
3127  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3128  return static_cast<CoolPropDbl>(_d4alphar_dDelta_dTau3);
3129 }
3130 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dTau4(void) {
3131  calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3132  return static_cast<CoolPropDbl>(_d4alphar_dTau4);
3133 }
3134 
3135 CoolPropDbl HelmholtzEOSMixtureBackend::calc_alpha0(void) {
3136  const int nTau = 0, nDelta = 0;
3137  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3138 }
3139 CoolPropDbl HelmholtzEOSMixtureBackend::calc_dalpha0_dDelta(void) {
3140  const int nTau = 0, nDelta = 1;
3141  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3142 }
3143 CoolPropDbl HelmholtzEOSMixtureBackend::calc_dalpha0_dTau(void) {
3144  const int nTau = 1, nDelta = 0;
3145  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3146 }
3147 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alpha0_dDelta2(void) {
3148  const int nTau = 0, nDelta = 2;
3149  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3150 }
3151 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alpha0_dDelta_dTau(void) {
3152  const int nTau = 1, nDelta = 1;
3153  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3154 }
3155 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alpha0_dTau2(void) {
3156  const int nTau = 2, nDelta = 0;
3157  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3158 }
3159 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alpha0_dDelta3(void) {
3160  const int nTau = 0, nDelta = 3;
3161  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3162 }
3163 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alpha0_dDelta2_dTau(void) {
3164  const int nTau = 1, nDelta = 2;
3165  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3166 }
3167 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alpha0_dDelta_dTau2(void) {
3168  const int nTau = 2, nDelta = 1;
3169  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3170 }
3171 CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alpha0_dTau3(void) {
3172  const int nTau = 3, nDelta = 0;
3173  return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3174 }
3175 CoolPropDbl HelmholtzEOSMixtureBackend::calc_first_saturation_deriv(parameters Of1, parameters Wrt1, HelmholtzEOSMixtureBackend& SatL,
3176  HelmholtzEOSMixtureBackend& SatV) {
3177  // Derivative of temperature w.r.t. pressure ALONG the saturation curve
3178  CoolPropDbl dTdP_sat = T() * (1 / SatV.rhomolar() - 1 / SatL.rhomolar()) / (SatV.hmolar() - SatL.hmolar());
3179 
3180  // "Trivial" inputs
3181  if (Of1 == iT && Wrt1 == iP) {
3182  return dTdP_sat;
3183  } else if (Of1 == iP && Wrt1 == iT) {
3184  return 1 / dTdP_sat;
3185  }
3186  // Derivative taken with respect to T
3187  else if (Wrt1 == iT) {
3188  return first_partial_deriv(Of1, iT, iP) + first_partial_deriv(Of1, iP, iT) / dTdP_sat;
3189  }
3190  // Derivative taken with respect to p
3191  else if (Wrt1 == iP) {
3192  return first_partial_deriv(Of1, iP, iT) + first_partial_deriv(Of1, iT, iP) * dTdP_sat;
3193  } else {
3194  throw ValueError(
3195  format("Not possible to take first saturation derivative with respect to %s", get_parameter_information(Wrt1, "short").c_str()));
3196  }
3197 }
3198 CoolPropDbl HelmholtzEOSMixtureBackend::calc_first_saturation_deriv(parameters Of1, parameters Wrt1) {
3199  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_first_saturation_deriv"));
3200 
3201  // Derivative of temperature w.r.t. pressure ALONG the saturation curve
3202  CoolPropDbl dTdP_sat = T() * (1 / SatV->rhomolar() - 1 / SatL->rhomolar()) / (SatV->hmolar() - SatL->hmolar());
3203 
3204  // "Trivial" inputs
3205  if (Of1 == iT && Wrt1 == iP) {
3206  return dTdP_sat;
3207  } else if (Of1 == iP && Wrt1 == iT) {
3208  return 1 / dTdP_sat;
3209  }
3210  // Derivative taken with respect to T
3211  else if (Wrt1 == iT) {
3212  return first_partial_deriv(Of1, iT, iP) + first_partial_deriv(Of1, iP, iT) / dTdP_sat;
3213  }
3214  // Derivative taken with respect to p
3215  else if (Wrt1 == iP) {
3216  return first_partial_deriv(Of1, iP, iT) + first_partial_deriv(Of1, iT, iP) * dTdP_sat;
3217  } else {
3218  throw ValueError(
3219  format("Not possible to take first saturation derivative with respect to %s", get_parameter_information(Wrt1, "short").c_str()));
3220  }
3221 }
3222 CoolPropDbl HelmholtzEOSMixtureBackend::calc_second_saturation_deriv(parameters Of1, parameters Wrt1, parameters Wrt2) {
3223  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_second_saturation_deriv"));
3224  if (Wrt1 == iP && Wrt2 == iP) {
3225  CoolPropDbl dydT_constp = this->first_partial_deriv(Of1, iT, iP);
3226  CoolPropDbl d2ydTdp = this->second_partial_deriv(Of1, iT, iP, iP, iT);
3227  CoolPropDbl d2ydp2_constT = this->second_partial_deriv(Of1, iP, iT, iP, iT);
3228  CoolPropDbl d2ydT2_constp = this->second_partial_deriv(Of1, iT, iP, iT, iP);
3229 
3230  CoolPropDbl dTdp_along_sat = calc_first_saturation_deriv(iT, iP);
3231  CoolPropDbl dvdrhoL = -1 / POW2(SatL->rhomolar());
3232  CoolPropDbl dvdrhoV = -1 / POW2(SatV->rhomolar());
3233  CoolPropDbl DELTAv = 1 / SatV->rhomolar() - 1 / SatL->rhomolar();
3234  CoolPropDbl dDELTAv_dT_constp = dvdrhoV * SatV->first_partial_deriv(iDmolar, iT, iP) - dvdrhoL * SatL->first_partial_deriv(iDmolar, iT, iP);
3235  CoolPropDbl dDELTAv_dp_constT = dvdrhoV * SatV->first_partial_deriv(iDmolar, iP, iT) - dvdrhoL * SatL->first_partial_deriv(iDmolar, iP, iT);
3236  CoolPropDbl DELTAh = SatV->hmolar() - SatL->hmolar();
3237  CoolPropDbl dDELTAh_dT_constp = SatV->first_partial_deriv(iHmolar, iT, iP) - SatL->first_partial_deriv(iHmolar, iT, iP);
3238  CoolPropDbl dDELTAh_dp_constT = SatV->first_partial_deriv(iHmolar, iP, iT) - SatL->first_partial_deriv(iHmolar, iP, iT);
3239  CoolPropDbl ddT_dTdp_along_sat_constp = (DELTAh * (_T * dDELTAv_dT_constp + DELTAv) - _T * DELTAv * dDELTAh_dT_constp) / POW2(DELTAh);
3240  CoolPropDbl ddp_dTdp_along_sat_constT = (DELTAh * (_T * dDELTAv_dp_constT) - _T * DELTAv * dDELTAh_dp_constT) / POW2(DELTAh);
3241 
3242  double ddp_dydpsigma = d2ydp2_constT + dydT_constp * ddp_dTdp_along_sat_constT + d2ydTdp * dTdp_along_sat;
3243  double ddT_dydpsigma = d2ydTdp + dydT_constp * ddT_dTdp_along_sat_constp + d2ydT2_constp * dTdp_along_sat;
3244  return ddp_dydpsigma + ddT_dydpsigma * dTdp_along_sat;
3245  } else {
3246  throw ValueError(format("Currently, only possible to take second saturation derivative w.r.t. P (both times)"));
3247  }
3248 }
3249 
3250 CoolPropDbl HelmholtzEOSMixtureBackend::calc_second_two_phase_deriv(parameters Of, parameters Wrt1, parameters Constant1, parameters Wrt2,
3251  parameters Constant2) {
3252  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_second_two_phase_deriv"));
3253 
3254  if (Of == iDmolar
3255  && ((Wrt1 == iHmolar && Constant1 == iP && Wrt2 == iP && Constant2 == iHmolar)
3256  || (Wrt2 == iHmolar && Constant2 == iP && Wrt1 == iP && Constant1 == iHmolar))) {
3257  parameters h_key = iHmolar, rho_key = iDmolar, p_key = iP;
3258  // taking the derivative of (drho/dv)*(dv/dh|p) with respect to p with h constant
3259  CoolPropDbl dv_dh_constp = calc_first_two_phase_deriv(rho_key, h_key, p_key) / (-POW2(rhomolar()));
3260  CoolPropDbl drhomolar_dp__consth = first_two_phase_deriv(rho_key, p_key, h_key);
3261 
3262  // Calculate the derivative of dvdh|p with respect to p at constant h
3263  CoolPropDbl dhL_dp_sat = SatL->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
3264  CoolPropDbl dhV_dp_sat = SatV->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
3265  CoolPropDbl drhoL_dp_sat = SatL->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
3266  CoolPropDbl drhoV_dp_sat = SatV->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
3267  CoolPropDbl numerator = 1 / SatV->keyed_output(rho_key) - 1 / SatL->keyed_output(rho_key);
3268  CoolPropDbl denominator = SatV->keyed_output(h_key) - SatL->keyed_output(h_key);
3269  CoolPropDbl dnumerator = -1 / POW2(SatV->keyed_output(rho_key)) * drhoV_dp_sat + 1 / POW2(SatL->keyed_output(rho_key)) * drhoL_dp_sat;
3270  CoolPropDbl ddenominator = dhV_dp_sat - dhL_dp_sat;
3271  CoolPropDbl d_dvdh_dp__consth = (denominator * dnumerator - numerator * ddenominator) / POW2(denominator);
3272  return -POW2(rhomolar()) * d_dvdh_dp__consth + dv_dh_constp * (-2 * rhomolar()) * drhomolar_dp__consth;
3273  } else if (Of == iDmass
3274  && ((Wrt1 == iHmass && Constant1 == iP && Wrt2 == iP && Constant2 == iHmass)
3275  || (Wrt2 == iHmass && Constant2 == iP && Wrt1 == iP && Constant1 == iHmass))) {
3276  parameters h_key = iHmass, rho_key = iDmass, p_key = iP;
3277  CoolPropDbl rho = keyed_output(rho_key);
3278  // taking the derivative of (drho/dv)*(dv/dh|p) with respect to p with h constant
3279  CoolPropDbl dv_dh_constp = calc_first_two_phase_deriv(rho_key, h_key, p_key) / (-POW2(rho));
3280  CoolPropDbl drho_dp__consth = first_two_phase_deriv(rho_key, p_key, h_key);
3281 
3282  // Calculate the derivative of dvdh|p with respect to p at constant h
3283  CoolPropDbl dhL_dp_sat = SatL->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
3284  CoolPropDbl dhV_dp_sat = SatV->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
3285  CoolPropDbl drhoL_dp_sat = SatL->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
3286  CoolPropDbl drhoV_dp_sat = SatV->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
3287  CoolPropDbl numerator = 1 / SatV->keyed_output(rho_key) - 1 / SatL->keyed_output(rho_key);
3288  CoolPropDbl denominator = SatV->keyed_output(h_key) - SatL->keyed_output(h_key);
3289  CoolPropDbl dnumerator = -1 / POW2(SatV->keyed_output(rho_key)) * drhoV_dp_sat + 1 / POW2(SatL->keyed_output(rho_key)) * drhoL_dp_sat;
3290  CoolPropDbl ddenominator = dhV_dp_sat - dhL_dp_sat;
3291  CoolPropDbl d_dvdh_dp__consth = (denominator * dnumerator - numerator * ddenominator) / POW2(denominator);
3292  return -POW2(rho) * d_dvdh_dp__consth + dv_dh_constp * (-2 * rho) * drho_dp__consth;
3293  } else {
3294  throw ValueError();
3295  }
3296 }
3297 CoolPropDbl HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant) {
3298  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_first_two_phase_deriv"));
3299  if (Of == iDmolar && Wrt == iHmolar && Constant == iP) {
3300  return -POW2(rhomolar()) * (1 / SatV->rhomolar() - 1 / SatL->rhomolar()) / (SatV->hmolar() - SatL->hmolar());
3301  } else if (Of == iDmass && Wrt == iHmass && Constant == iP) {
3302  return -POW2(rhomass()) * (1 / SatV->rhomass() - 1 / SatL->rhomass()) / (SatV->hmass() - SatL->hmass());
3303  } else if (Of == iDmolar && Wrt == iP && Constant == iHmolar) {
3304  // v = 1/rho; drhodv = -rho^2; dvdrho = -1/rho^2
3305  CoolPropDbl dvdrhoL = -1 / POW2(SatL->rhomolar());
3306  CoolPropDbl dvdrhoV = -1 / POW2(SatV->rhomolar());
3307  CoolPropDbl dvL_dp = dvdrhoL * SatL->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
3308  CoolPropDbl dvV_dp = dvdrhoV * SatV->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
3309  CoolPropDbl dhL_dp = SatL->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
3310  CoolPropDbl dhV_dp = SatV->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
3311  CoolPropDbl dxdp_h = (Q() * dhV_dp + (1 - Q()) * dhL_dp) / (SatL->hmolar() - SatV->hmolar());
3312  CoolPropDbl dvdp_h = dvL_dp + dxdp_h * (1 / SatV->rhomolar() - 1 / SatL->rhomolar()) + Q() * (dvV_dp - dvL_dp);
3313  return -POW2(rhomolar()) * dvdp_h;
3314  } else if (Of == iDmass && Wrt == iP && Constant == iHmass) {
3315  // v = 1/rho; drhodv = -rho^2; dvdrho = -1/rho^2
3316  CoolPropDbl dvdrhoL = -1 / POW2(SatL->rhomass());
3317  CoolPropDbl dvdrhoV = -1 / POW2(SatV->rhomass());
3318  CoolPropDbl dvL_dp = dvdrhoL * SatL->calc_first_saturation_deriv(iDmass, iP, *SatL, *SatV);
3319  CoolPropDbl dvV_dp = dvdrhoV * SatV->calc_first_saturation_deriv(iDmass, iP, *SatL, *SatV);
3320  CoolPropDbl dhL_dp = SatL->calc_first_saturation_deriv(iHmass, iP, *SatL, *SatV);
3321  CoolPropDbl dhV_dp = SatV->calc_first_saturation_deriv(iHmass, iP, *SatL, *SatV);
3322  CoolPropDbl dxdp_h = (Q() * dhV_dp + (1 - Q()) * dhL_dp) / (SatL->hmass() - SatV->hmass());
3323  CoolPropDbl dvdp_h = dvL_dp + dxdp_h * (1 / SatV->rhomass() - 1 / SatL->rhomass()) + Q() * (dvV_dp - dvL_dp);
3324  return -POW2(rhomass()) * dvdp_h;
3325  } else {
3326  throw ValueError("These inputs are not supported to calc_first_two_phase_deriv");
3327  }
3328 }
3329 CoolPropDbl HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, CoolPropDbl x_end) {
3330  // Note: If you need all three values (drho_dh__p, drho_dp__h and rho_spline),
3331  // you should calculate drho_dp__h first to avoid duplicate calculations.
3332 
3333  bool drho_dh__p = false;
3334  bool drho_dp__h = false;
3335  bool rho_spline = false;
3336 
3337  if (Of == iDmolar && Wrt == iHmolar && Constant == iP) {
3338  drho_dh__p = true;
3339  if (_drho_spline_dh__constp) return _drho_spline_dh__constp;
3340  } else if (Of == iDmass && Wrt == iHmass && Constant == iP) {
3341  return first_two_phase_deriv_splined(iDmolar, iHmolar, iP, x_end) * POW2(molar_mass());
3342  } else if (Of == iDmolar && Wrt == iP && Constant == iHmolar) {
3343  drho_dp__h = true;
3344  if (_drho_spline_dp__consth) return _drho_spline_dp__consth;
3345  } else if (Of == iDmass && Wrt == iP && Constant == iHmass) {
3346  return first_two_phase_deriv_splined(iDmolar, iP, iHmolar, x_end) * molar_mass();
3347  }
3348  // Add the special case for the splined density
3349  else if (Of == iDmolar && Wrt == iDmolar && Constant == iDmolar) {
3350  rho_spline = true;
3351  if (_rho_spline) return _rho_spline;
3352  } else if (Of == iDmass && Wrt == iDmass && Constant == iDmass) {
3353  return first_two_phase_deriv_splined(iDmolar, iDmolar, iDmolar, x_end) * molar_mass();
3354  } else {
3355  throw ValueError("These inputs are not supported to calc_first_two_phase_deriv");
3356  }
3357 
3358  if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_first_two_phase_deriv_splined"));
3359  if (_Q > x_end) {
3360  throw ValueError(format("Q [%g] is greater than x_end [%Lg]", _Q, x_end).c_str());
3361  }
3362  if (_phase != iphase_twophase) {
3363  throw ValueError(format("state is not two-phase"));
3364  }
3365 
3366  shared_ptr<HelmholtzEOSMixtureBackend> Liq(new HelmholtzEOSMixtureBackend(this->get_components())),
3367  End(new HelmholtzEOSMixtureBackend(this->get_components()));
3368 
3369  Liq->specify_phase(iphase_liquid);
3370  Liq->_Q = -1;
3371  Liq->update_DmolarT_direct(SatL->rhomolar(), SatL->T());
3372  End->update(QT_INPUTS, x_end, SatL->T());
3373 
3374  CoolPropDbl Delta = Q() * (SatV->keyed_output(iHmolar) - SatL->keyed_output(iHmolar));
3375  CoolPropDbl Delta_end = End->keyed_output(iHmolar) - SatL->keyed_output(iHmolar);
3376 
3377  // At the end of the zone to which spline is applied
3378  CoolPropDbl drho_dh_end = End->calc_first_two_phase_deriv(iDmolar, iHmolar, iP);
3379  CoolPropDbl rho_end = End->keyed_output(iDmolar);
3380 
3381  // Faking single-phase
3382  CoolPropDbl rho_liq = Liq->keyed_output(iDmolar);
3383  CoolPropDbl drho_dh_liq__constp = Liq->first_partial_deriv(iDmolar, iHmolar, iP);
3384 
3385  // Spline coordinates a, b, c, d
3386  CoolPropDbl Abracket = (2 * rho_liq - 2 * rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end));
3387  CoolPropDbl a = 1 / POW3(Delta_end) * Abracket;
3388  CoolPropDbl b = 3 / POW2(Delta_end) * (-rho_liq + rho_end) - 1 / Delta_end * (drho_dh_end + 2 * drho_dh_liq__constp);
3389  CoolPropDbl c = drho_dh_liq__constp;
3390  CoolPropDbl d = rho_liq;
3391 
3392  // Either the spline value or drho/dh|p can be directly evaluated now
3393  _rho_spline = a * POW3(Delta) + b * POW2(Delta) + c * Delta + d;
3394  _drho_spline_dh__constp = 3 * a * POW2(Delta) + 2 * b * Delta + c;
3395  if (rho_spline) return _rho_spline;
3396  if (drho_dh__p) return _drho_spline_dh__constp;
3397 
3398  // It's drho/dp|h
3399  // ... calculate some more things
3400 
3401  // Derivatives *along* the saturation curve using the special internal method
3402  CoolPropDbl dhL_dp_sat = SatL->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
3403  CoolPropDbl dhV_dp_sat = SatV->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
3404  CoolPropDbl drhoL_dp_sat = SatL->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
3405  CoolPropDbl drhoV_dp_sat = SatV->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
3406  CoolPropDbl rhoV = SatV->keyed_output(iDmolar);
3407  CoolPropDbl rhoL = SatL->keyed_output(iDmolar);
3408  CoolPropDbl drho_dp_end = POW2(End->keyed_output(iDmolar)) * (x_end / POW2(rhoV) * drhoV_dp_sat + (1 - x_end) / POW2(rhoL) * drhoL_dp_sat);
3409 
3410  // Faking single-phase
3411  //CoolPropDbl drho_dp__consth_liq = Liq->first_partial_deriv(iDmolar, iP, iHmolar);
3412  CoolPropDbl d2rhodhdp_liq = Liq->second_partial_deriv(iDmolar, iHmolar, iP, iP, iHmolar); // ?
3413 
3414  // Derivatives at the end point
3415  // CoolPropDbl drho_dp__consth_end = End->calc_first_two_phase_deriv(iDmolar, iP, iHmolar);
3416  CoolPropDbl d2rhodhdp_end = End->calc_second_two_phase_deriv(iDmolar, iHmolar, iP, iP, iHmolar);
3417 
3418  // Reminder:
3419  // Delta = Q()*(hV-hL) = h-hL
3420  // Delta_end = x_end*(hV-hL);
3421  CoolPropDbl d_Delta_dp__consth = -dhL_dp_sat;
3422  CoolPropDbl d_Delta_end_dp__consth = x_end * (dhV_dp_sat - dhL_dp_sat);
3423 
3424  // First pressure derivative at constant h of the coefficients a,b,c,d
3425  // CoolPropDbl Abracket = (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end));
3426  CoolPropDbl d_Abracket_dp_consth = (2 * drhoL_dp_sat - 2 * drho_dp_end + Delta_end * (d2rhodhdp_liq + d2rhodhdp_end)
3427  + d_Delta_end_dp__consth * (drho_dh_liq__constp + drho_dh_end));
3428  CoolPropDbl da_dp = 1 / POW3(Delta_end) * d_Abracket_dp_consth + Abracket * (-3 / POW4(Delta_end) * d_Delta_end_dp__consth);
3429  CoolPropDbl db_dp = -6 / POW3(Delta_end) * d_Delta_end_dp__consth * (rho_end - rho_liq) + 3 / POW2(Delta_end) * (drho_dp_end - drhoL_dp_sat)
3430  + (1 / POW2(Delta_end) * d_Delta_end_dp__consth) * (drho_dh_end + 2 * drho_dh_liq__constp)
3431  - (1 / Delta_end) * (d2rhodhdp_end + 2 * d2rhodhdp_liq);
3432  CoolPropDbl dc_dp = d2rhodhdp_liq;
3433  CoolPropDbl dd_dp = drhoL_dp_sat;
3434 
3435  _drho_spline_dp__consth =
3436  (3 * a * POW2(Delta) + 2 * b * Delta + c) * d_Delta_dp__consth + POW3(Delta) * da_dp + POW2(Delta) * db_dp + Delta * dc_dp + dd_dp;
3437  if (drho_dp__h) return _drho_spline_dp__consth;
3438 
3439  throw ValueError("Something went wrong in HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv_splined");
3440  return _HUGE;
3441 }
3442 
3443 CoolProp::CriticalState HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0) {
3444  class Resid : public FuncWrapperND
3445  {
3446  public:
3447  HelmholtzEOSMixtureBackend& HEOS;
3448  double L1, M1;
3449  Eigen::MatrixXd Lstar, Mstar;
3450  Resid(HelmholtzEOSMixtureBackend& HEOS) : HEOS(HEOS), L1(_HUGE), M1(_HUGE){};
3451  std::vector<double> call(const std::vector<double>& tau_delta) {
3452  double rhomolar = tau_delta[1] * HEOS.rhomolar_reducing();
3453  double T = HEOS.T_reducing() / tau_delta[0];
3454  HEOS.update(DmolarT_INPUTS, rhomolar, T);
3455  Lstar = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT);
3456  Mstar = MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT, Lstar);
3457  std::vector<double> o(2);
3458  o[0] = Lstar.determinant();
3459  o[1] = Mstar.determinant();
3460  return o;
3461  };
3462  std::vector<std::vector<double>> Jacobian(const std::vector<double>& x) {
3463  std::size_t N = x.size();
3464  std::vector<std::vector<double>> J(N, std::vector<double>(N, 0));
3465  Eigen::MatrixXd adjL = adjugate(Lstar), adjM = adjugate(Mstar), dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau),
3466  dLdDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta),
3467  dMdTau = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iTau, Lstar, dLdTau),
3468  dMdDelta = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iDelta, Lstar, dLdDelta);
3469 
3470  J[0][0] = (adjL * dLdTau).trace();
3471  J[0][1] = (adjL * dLdDelta).trace();
3472  J[1][0] = (adjM * dMdTau).trace();
3473  J[1][1] = (adjM * dMdDelta).trace();
3474  return J;
3475  }
3477  std::vector<std::vector<double>> numerical_Jacobian(const std::vector<double>& x) {
3478  std::size_t N = x.size();
3479  std::vector<double> rplus, rminus, xp;
3480  std::vector<std::vector<double>> J(N, std::vector<double>(N, 0));
3481 
3482  double epsilon = 0.0001;
3483 
3484  // Build the Jacobian by column
3485  for (std::size_t i = 0; i < N; ++i) {
3486  xp = x;
3487  xp[i] += epsilon;
3488  rplus = call(xp);
3489  xp[i] -= 2 * epsilon;
3490  rminus = call(xp);
3491 
3492  for (std::size_t j = 0; j < N; ++j) {
3493  J[j][i] = (rplus[j] - rminus[j]) / (2 * epsilon);
3494  }
3495  }
3496  std::cout << J[0][0] << " " << J[0][1] << std::endl;
3497  std::cout << J[1][0] << " " << J[1][1] << std::endl;
3498  return J;
3499  };
3500  };
3501  Resid resid(*this);
3502  std::vector<double> x, tau_delta(2);
3503  tau_delta[0] = T_reducing() / T0;
3504  tau_delta[1] = rho0 / rhomolar_reducing();
3505  x = NDNewtonRaphson_Jacobian(&resid, tau_delta, 1e-10, 100);
3506  _critical.T = T_reducing() / x[0];
3507  _critical.rhomolar = x[1] * rhomolar_reducing();
3508  _critical.p = calc_pressure_nocache(_critical.T, _critical.rhomolar);
3509 
3510  CriticalState critical;
3511  critical.T = _critical.T;
3512  critical.p = _critical.p;
3513  critical.rhomolar = _critical.rhomolar;
3514  if (_critical.p < 0) {
3515  critical.stable = false;
3516  } else {
3517  if (get_config_bool(ASSUME_CRITICAL_POINT_STABLE)) {
3518  critical.stable = true;
3519  } else {
3520  // Otherwise we try to check stability with TPD-based analysis
3521  StabilityRoutines::StabilityEvaluationClass stability_tester(*this);
3522  critical.stable = stability_tester.is_stable();
3523  }
3524  }
3525  return critical;
3526 }
3527 
3531 class OneDimObjective : public FuncWrapper1DWithTwoDerivs
3532 {
3533  public:
3534  CoolProp::HelmholtzEOSMixtureBackend& HEOS;
3535  const double delta;
3536  double _call, _deriv, _second_deriv;
3537  OneDimObjective(HelmholtzEOSMixtureBackend& HEOS, double delta0) : HEOS(HEOS), delta(delta0), _call(_HUGE), _deriv(_HUGE), _second_deriv(_HUGE){};
3538  double call(double tau) {
3539  double rhomolar = HEOS.rhomolar_reducing() * delta, T = HEOS.T_reducing() / tau;
3540  HEOS.update_DmolarT_direct(rhomolar, T);
3541  _call = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT).determinant();
3542  return _call;
3543  }
3544  double deriv(double tau) {
3545  Eigen::MatrixXd adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)),
3546  dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau);
3547  _deriv = (adjL * dLdTau).trace();
3548  return _deriv;
3549  };
3550  double second_deriv(double tau) {
3551  Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT),
3552  dLstardTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau),
3553  d2LstardTau2 = MixtureDerivatives::d2Lstar_dX2(HEOS, XN_INDEPENDENT, iTau, iTau), adjL = adjugate(Lstar),
3554  dadjLstardTau = adjugate_derivative(Lstar, dLstardTau);
3555  _second_deriv = (dLstardTau * dadjLstardTau + adjL * d2LstardTau2).trace();
3556  return _second_deriv;
3557  };
3558 };
3559 
3562 class L0CurveTracer : public FuncWrapper1DWithDeriv
3563 {
3564  public:
3565  CoolProp::HelmholtzEOSMixtureBackend& HEOS;
3566  double delta, tau,
3567  M1_last,
3568  theta_last,
3569  R_tau,
3570  R_delta,
3571  R_tau_tracer,
3572  R_delta_tracer;
3573  std::vector<CoolProp::CriticalState> critical_points;
3574  int N_critical_points;
3575  Eigen::MatrixXd Lstar, adjLstar, dLstardTau, d2LstardTau2, dLstardDelta;
3576  SpinodalData spinodal_values;
3577  bool
3578  find_critical_points;
3579  L0CurveTracer(HelmholtzEOSMixtureBackend& HEOS, double tau0, double delta0)
3580  : HEOS(HEOS), delta(delta0), tau(tau0), M1_last(_HUGE), N_critical_points(0), find_critical_points(true) {
3581  R_delta_tracer = 0.1;
3582  R_delta = R_delta_tracer;
3583  R_tau_tracer = 0.1;
3584  R_tau = R_tau_tracer;
3585  };
3586  /***
3587  \brief Update values for tau, delta
3588  @param theta The angle
3589  @param tau The old value of tau
3590  @param delta The old value of delta
3591  @param tau_new The new value of tau
3592  @param delta_new The new value of delta
3593  */
3594  void get_tau_delta(const double theta, const double tau, const double delta, double& tau_new, double& delta_new) {
3595  tau_new = tau + R_tau * cos(theta);
3596  delta_new = delta + R_delta * sin(theta);
3597  };
3598  /***
3599  \brief Calculate the value of L1
3600  @param theta The angle
3601  */
3602  double call(double theta) {
3603  double tau_new, delta_new;
3604  this->get_tau_delta(theta, tau, delta, tau_new, delta_new);
3605  double rhomolar = HEOS.rhomolar_reducing() * delta_new, T = HEOS.T_reducing() / tau_new;
3606  HEOS.update_DmolarT_direct(rhomolar, T);
3607  Lstar = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT);
3608  adjLstar = adjugate(Lstar);
3609  dLstardTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau);
3610  dLstardDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta);
3611  double L1 = Lstar.determinant();
3612  return L1;
3613  };
3614  /***
3615  \brief Calculate the first partial derivative of L1 with respect to the angle
3616  @param theta The angle
3617  */
3618  double deriv(double theta) {
3619  double dL1_dtau = (adjLstar * dLstardTau).trace(), dL1_ddelta = (adjLstar * dLstardDelta).trace();
3620  return -R_tau * sin(theta) * dL1_dtau + R_delta * cos(theta) * dL1_ddelta;
3621  };
3622 
3623  void trace() {
3624  bool debug = (get_debug_level() > 0) | false;
3625  double theta;
3626  for (int i = 0; i < 300; ++i) {
3627  if (i == 0) {
3628  // In the first iteration, search all angles in the positive delta direction using a
3629  // bounded solver with a very small radius in order to not hit other L1*=0 contours
3630  // that are in the vicinity
3631  R_tau = 0.001;
3632  R_delta = 0.001;
3633  theta = Brent(this, 0, M_PI, DBL_EPSILON, 1e-10, 100);
3634  } else {
3635  // In subsequent iterations, you already have an excellent guess for the direction to
3636  // be searching in, use Newton's method to refine the solution since we also
3637  // have an analytic solution for the derivative
3638  R_tau = R_tau_tracer;
3639  R_delta = R_delta_tracer;
3640  try {
3641  theta = Newton(this, theta_last, 1e-10, 100);
3642  } catch (...) {
3643  if (N_critical_points > 0 && delta > 1.5 * critical_points[0].rhomolar / HEOS.rhomolar_reducing()) {
3644  // Stopping the search - probably we have a kink in the L1*=0 contour
3645  // caused by poor low-temperature behavior
3646  continue;
3647  }
3648  }
3649 
3650  // If the solver takes a U-turn, going in the opposite direction of travel
3651  // this is not a good thing, and force a Brent's method solver call to make sure we keep
3652  // tracing in the same direction
3653  if (std::abs(angle_difference(theta, theta_last)) > M_PI / 2.0) {
3654  // We have found at least one critical point and we are now well above the density
3655  // associated with the first critical point
3656  if (N_critical_points > 0 && delta > 1.2 * critical_points[0].rhomolar / HEOS.rhomolar_reducing()) {
3657  // Stopping the search - probably we have a kink in the L1*=0 contour
3658  // caused by poor low-temperature behavior
3659  continue;
3660  } else {
3661  theta = Brent(this, theta_last - M_PI / 3.5, theta_last + M_PI / 3.5, DBL_EPSILON, 1e-10, 100);
3662  }
3663  }
3664  }
3665 
3666  // Calculate the second criticality condition
3667  double M1 = MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT, Lstar).determinant();
3668  double p_MPa = HEOS.p() / 1e6;
3669 
3670  // Calculate the new tau and delta at the new point
3671  double tau_new, delta_new;
3672  this->get_tau_delta(theta, tau, delta, tau_new, delta_new);
3673 
3674  // Stop if bounds are exceeded
3675  if (p_MPa > 500 || HEOS.get_critical_is_terminated(delta_new, tau_new)) {
3676  break;
3677  }
3678 
3679  // If the sign of M1 and the previous value of M1 have different signs, it means that
3680  // you have bracketed a critical point, run the full critical point solver to
3681  // find the critical point and store it
3682  // Only enabled if find_critical_points is true (the default)
3683  if (i > 0 && M1 * M1_last < 0 && find_critical_points) {
3684  double rhomolar = HEOS.rhomolar_reducing() * (delta + delta_new) / 2.0, T = HEOS.T_reducing() / ((tau + tau_new) / 2.0);
3685  CoolProp::CriticalState crit = HEOS.calc_critical_point(rhomolar, T);
3686  critical_points.push_back(crit);
3687  N_critical_points++;
3688  if (debug) {
3689  std::cout << HEOS.get_mole_fractions()[0] << " " << crit.rhomolar << " " << crit.T << " " << p_MPa << std::endl;
3690  }
3691  }
3692 
3693  // Update the storage values
3694  this->tau = tau_new;
3695  this->delta = delta_new;
3696  this->M1_last = M1;
3697  this->theta_last = theta;
3698 
3699  this->spinodal_values.tau.push_back(tau_new);
3700  this->spinodal_values.delta.push_back(delta_new);
3701  this->spinodal_values.M1.push_back(M1);
3702  }
3703  };
3704 };
3705 
3706 void HelmholtzEOSMixtureBackend::calc_criticality_contour_values(double& L1star, double& M1star) {
3707  Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(*this, XN_INDEPENDENT);
3708  Eigen::MatrixXd Mstar = MixtureDerivatives::Mstar(*this, XN_INDEPENDENT, Lstar);
3709  L1star = Lstar.determinant();
3710  M1star = Mstar.determinant();
3711 };
3712 
3713 void HelmholtzEOSMixtureBackend::get_critical_point_search_radii(double& R_delta, double& R_tau) {
3714  R_delta = 0.025;
3715  R_tau = 0.1;
3716 }
3717 std::vector<CoolProp::CriticalState> HelmholtzEOSMixtureBackend::_calc_all_critical_points(bool find_critical_points) {
3718  // Populate the temporary class used to calculate the critical point(s)
3719  add_critical_state();
3720  if (get_debug_level() > 10) {
3721  rapidjson::Document doc;
3722  doc.SetObject();
3723  rapidjson::Value& val = doc;
3724  std::vector<std::vector<DepartureFunctionPointer>>& mat = critical_state->residual_helmholtz->Excess.DepartureFunctionMatrix;
3725  if (mat.size() > 0) {
3726  mat[0][1]->phi.to_json(val, doc);
3727  std::cout << cpjson::to_string(doc);
3728  }
3729  }
3730  critical_state->set_mole_fractions(this->get_mole_fractions_ref());
3731 
3732  // Specify state to be something homogeneous to shortcut phase evaluation
3733  critical_state->specify_phase(iphase_gas);
3734 
3735  double delta0 = _HUGE, tau0 = _HUGE;
3736  critical_state->get_critical_point_starting_values(delta0, tau0);
3737 
3738  OneDimObjective resid_L0(*critical_state, delta0);
3739 
3740  // If the derivative of L1star with respect to tau is positive,
3741  // tau needs to be increased such that we sit on the other
3742  // side of the d(L1star)/dtau = 0 contour
3743  resid_L0.call(tau0);
3744  int bump_count = 0;
3745  while (resid_L0.deriv(tau0) > 0 && bump_count < 3) {
3746  tau0 *= 1.1;
3747  bump_count++;
3748  }
3749  double tau_L0 = Halley(resid_L0, tau0, 1e-10, 100);
3750  //double T0 = T_reducing()/tau_L0;
3751  //double rho0 = delta0*rhomolar_reducing();
3752 
3753  L0CurveTracer tracer(*critical_state, tau_L0, delta0);
3754  tracer.find_critical_points = find_critical_points;
3755 
3756  double R_delta = 0, R_tau = 0;
3757  critical_state->get_critical_point_search_radii(R_delta, R_tau);
3758  tracer.R_delta_tracer = R_delta;
3759  tracer.R_tau_tracer = R_tau;
3760  tracer.trace();
3761 
3762  this->spinodal_values = tracer.spinodal_values;
3763 
3764  return tracer.critical_points;
3765 }
3766 
3767 double HelmholtzEOSMixtureBackend::calc_tangent_plane_distance(const double T, const double p, const std::vector<double>& w,
3768  const double rhomolar_guess) {
3769 
3770  const std::vector<CoolPropDbl>& z = this->get_mole_fractions_ref();
3771  if (w.size() != z.size()) {
3772  throw ValueError(format("Trial composition vector size [%d] is not the same as bulk composition [%d]", w.size(), z.size()));
3773  }
3774  add_TPD_state();
3775  TPD_state->set_mole_fractions(w);
3776 
3777  CoolPropDbl rho = TPD_state->solver_rho_Tp_global(T, p, 0.9 / TPD_state->SRK_covolume());
3778  TPD_state->update_DmolarT_direct(rho, T);
3779 
3780  double summer = 0;
3781  for (std::size_t i = 0; i < w.size(); ++i) {
3782  summer +=
3783  w[i] * (log(MixtureDerivatives::fugacity_i(*TPD_state, i, XN_DEPENDENT)) - log(MixtureDerivatives::fugacity_i(*this, i, XN_DEPENDENT)));
3784  }
3785  return summer;
3786 }
3787 
3788 void HelmholtzEOSMixtureBackend::calc_build_spinodal() {
3789  // Ok, we are faking a little bit here, hijacking the code for critical points, but skipping evaluation of critical points
3790  bool find_critical_points = false;
3791  _calc_all_critical_points(find_critical_points);
3792 }
3793 
3794 void HelmholtzEOSMixtureBackend::set_reference_stateS(const std::string& reference_state) {
3795  for (std::size_t i = 0; i < components.size(); ++i) {
3796  CoolProp::HelmholtzEOSMixtureBackend HEOS(std::vector<CoolPropFluid>(1, components[i]));
3797  if (!reference_state.compare("IIR")) {
3798  if (HEOS.Ttriple() > 273.15) {
3799  throw ValueError(format("Cannot use IIR reference state; Ttriple [%Lg] is greater than 273.15 K", HEOS.Ttriple()));
3800  }
3801  HEOS.update(QT_INPUTS, 0, 273.15);
3802 
3803  // Get current values for the enthalpy and entropy
3804  double deltah = HEOS.hmass() - 200000; // offset from 200000 J/kg enthalpy
3805  double deltas = HEOS.smass() - 1000; // offset from 1000 J/kg/K entropy
3806  double delta_a1 = deltas / (HEOS.gas_constant() / HEOS.molar_mass());
3807  double delta_a2 = -deltah / (HEOS.gas_constant() / HEOS.molar_mass() * HEOS.get_reducing_state().T);
3808  // Change the value in the library for the given fluid
3809  set_fluid_enthalpy_entropy_offset(components[i], delta_a1, delta_a2, "IIR");
3810  if (get_debug_level() > 0) {
3811  std::cout << format("set offsets to %0.15g and %0.15g\n", delta_a1, delta_a2);
3812  }
3813  } else if (!reference_state.compare("ASHRAE")) {
3814  if (HEOS.Ttriple() > 233.15) {
3815  throw ValueError(format("Cannot use ASHRAE reference state; Ttriple [%Lg] is greater than than 233.15 K", HEOS.Ttriple()));
3816  }
3817  HEOS.update(QT_INPUTS, 0, 233.15);
3818 
3819  // Get current values for the enthalpy and entropy
3820  double deltah = HEOS.hmass() - 0; // offset from 0 J/kg enthalpy
3821  double deltas = HEOS.smass() - 0; // offset from 0 J/kg/K entropy
3822  double delta_a1 = deltas / (HEOS.gas_constant() / HEOS.molar_mass());
3823  double delta_a2 = -deltah / (HEOS.gas_constant() / HEOS.molar_mass() * HEOS.get_reducing_state().T);
3824  // Change the value in the library for the given fluid
3825  set_fluid_enthalpy_entropy_offset(components[i], delta_a1, delta_a2, "ASHRAE");
3826  if (get_debug_level() > 0) {
3827  std::cout << format("set offsets to %0.15g and %0.15g\n", delta_a1, delta_a2);
3828  }
3829  } else if (!reference_state.compare("NBP")) {
3830  if (HEOS.p_triple() > 101325) {
3831  throw ValueError(format("Cannot use NBP reference state; p_triple [%Lg Pa] is greater than than 101325 Pa", HEOS.p_triple()));
3832  }
3833  // Saturated liquid boiling point at 1 atmosphere
3834  HEOS.update(PQ_INPUTS, 101325, 0);
3835 
3836  double deltah = HEOS.hmass() - 0; // offset from 0 kJ/kg enthalpy
3837  double deltas = HEOS.smass() - 0; // offset from 0 kJ/kg/K entropy
3838  double delta_a1 = deltas / (HEOS.gas_constant() / HEOS.molar_mass());
3839  double delta_a2 = -deltah / (HEOS.gas_constant() / HEOS.molar_mass() * HEOS.get_reducing_state().T);
3840  // Change the value in the library for the given fluid
3841  set_fluid_enthalpy_entropy_offset(components[i], delta_a1, delta_a2, "NBP");
3842  if (get_debug_level() > 0) {
3843  std::cout << format("set offsets to %0.15g and %0.15g\n", delta_a1, delta_a2);
3844  }
3845  } else if (!reference_state.compare("DEF")) {
3846  set_fluid_enthalpy_entropy_offset(components[i], 0, 0, "DEF");
3847  } else if (!reference_state.compare("RESET")) {
3848  set_fluid_enthalpy_entropy_offset(components[i], 0, 0, "RESET");
3849  } else {
3850  throw ValueError(format("reference state string is invalid: [%s]", reference_state.c_str()));
3851  }
3852  }
3853 }
3854 
3860 void HelmholtzEOSMixtureBackend::set_reference_stateD(double T, double rhomolar, double hmolar0, double smolar0) {
3861  for (std::size_t i = 0; i < components.size(); ++i) {
3862  CoolProp::HelmholtzEOSMixtureBackend HEOS(std::vector<CoolPropFluid>(1, components[i]));
3863 
3864  HEOS.update(DmolarT_INPUTS, rhomolar, T);
3865 
3866  // Get current values for the enthalpy and entropy
3867  double deltah = HEOS.hmolar() - hmolar0; // offset from specified enthalpy in J/mol
3868  double deltas = HEOS.smolar() - smolar0; // offset from specified entropy in J/mol/K
3869  double delta_a1 = deltas / (HEOS.gas_constant());
3870  double delta_a2 = -deltah / (HEOS.gas_constant() * HEOS.get_reducing_state().T);
3871  set_fluid_enthalpy_entropy_offset(components[i], delta_a1, delta_a2, "custom");
3872  }
3873 }
3874 
3875 void HelmholtzEOSMixtureBackend::set_fluid_enthalpy_entropy_offset(CoolPropFluid& component, double delta_a1, double delta_a2,
3876  const std::string& ref) {
3877  component.EOS().alpha0.EnthalpyEntropyOffset.set(delta_a1, delta_a2, ref);
3878 
3879  shared_ptr<CoolProp::HelmholtzEOSBackend> HEOS(new CoolProp::HelmholtzEOSBackend(component));
3880  HEOS->specify_phase(iphase_gas); // Something homogeneous;
3881  // Calculate the new enthalpy and entropy values
3882  HEOS->update(DmolarT_INPUTS, component.EOS().hs_anchor.rhomolar, component.EOS().hs_anchor.T);
3883  component.EOS().hs_anchor.hmolar = HEOS->hmolar();
3884  component.EOS().hs_anchor.smolar = HEOS->smolar();
3885 
3886  double f = (HEOS->name() == "Water" || HEOS->name() == "CarbonDioxide") ? 1.00001 : 1.0;
3887 
3888  // Calculate the new enthalpy and entropy values at the reducing state
3889  HEOS->update(DmolarT_INPUTS, component.EOS().reduce.rhomolar * f, component.EOS().reduce.T * f);
3890  component.EOS().reduce.hmolar = HEOS->hmolar();
3891  component.EOS().reduce.smolar = HEOS->smolar();
3892 
3893  // Calculate the new enthalpy and entropy values at the critical state
3894  HEOS->update(DmolarT_INPUTS, component.crit.rhomolar * f, component.crit.T * f);
3895  component.crit.hmolar = HEOS->hmolar();
3896  component.crit.smolar = HEOS->smolar();
3897 
3898  // Calculate the new enthalpy and entropy values
3899  HEOS->update(DmolarT_INPUTS, component.triple_liquid.rhomolar, component.triple_liquid.T);
3900  component.triple_liquid.hmolar = HEOS->hmolar();
3901  component.triple_liquid.smolar = HEOS->smolar();
3902 
3903  // Calculate the new enthalpy and entropy values
3904  HEOS->update(DmolarT_INPUTS, component.triple_vapor.rhomolar, component.triple_vapor.T);
3905  component.triple_vapor.hmolar = HEOS->hmolar();
3906  component.triple_vapor.smolar = HEOS->smolar();
3907 
3908  if (!HEOS->is_pure()) {
3909  // Calculate the new enthalpy and entropy values
3910  HEOS->update(DmolarT_INPUTS, component.EOS().max_sat_T.rhomolar, component.EOS().max_sat_T.T);
3911  component.EOS().max_sat_T.hmolar = HEOS->hmolar();
3912  component.EOS().max_sat_T.smolar = HEOS->smolar();
3913  // Calculate the new enthalpy and entropy values
3914  HEOS->update(DmolarT_INPUTS, component.EOS().max_sat_p.rhomolar, component.EOS().max_sat_p.T);
3915  component.EOS().max_sat_p.hmolar = HEOS->hmolar();
3916  component.EOS().max_sat_p.smolar = HEOS->smolar();
3917  }
3918 }
3919 
3920 } /* namespace CoolProp */
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1