CoolProp  6.6.0
An open-source fluid property and humid air property database
FlashRoutines.cpp
Go to the documentation of this file.
1 #include "VLERoutines.h"
2 #include "FlashRoutines.h"
4 #include "HelmholtzEOSBackend.h"
6 #include "Configuration.h"
7 
8 #if defined(ENABLE_CATCH)
9 # include <catch2/catch_all.hpp>
11 #endif
12 
13 namespace CoolProp {
14 
16  if (HEOS.PhaseEnvelope.built) {
17  // Use the phase envelope if already constructed to determine phase boundary
18  // Determine whether you are inside (two-phase) or outside (single-phase)
19  SimpleState closest_state;
20  std::size_t i;
21  bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS.PhaseEnvelope, iP, HEOS._p, iT, HEOS._T, i, closest_state);
22  if (!twophase && HEOS._T > closest_state.T) {
23  // Gas solution - bounded between phase envelope temperature and very high temperature
24  //
25  // Start with a guess value from SRK
26  CoolPropDbl rhomolar_guess = HEOS.solver_rho_Tp_SRK(HEOS._T, HEOS._p, iphase_gas);
27 
28  solver_TP_resid resid(HEOS, HEOS._T, HEOS._p);
29  std::string errstr;
31  try {
32  // Try using Newton's method
33  CoolPropDbl rhomolar = Newton(resid, rhomolar_guess, 1e-10, 100);
34  // Make sure the solution is within the bounds
35  if (!is_in_closed_range(static_cast<CoolPropDbl>(closest_state.rhomolar), static_cast<CoolPropDbl>(0.0), rhomolar)) {
36  throw ValueError("out of range");
37  }
38  HEOS.update_DmolarT_direct(rhomolar, HEOS._T);
39  } catch (...) {
40  // If that fails, try a bounded solver
41  CoolPropDbl rhomolar = Brent(resid, closest_state.rhomolar, 1e-10, DBL_EPSILON, 1e-10, 100);
42  // Make sure the solution is within the bounds
43  if (!is_in_closed_range(static_cast<CoolPropDbl>(closest_state.rhomolar), static_cast<CoolPropDbl>(0.0), rhomolar)) {
44  throw ValueError("out of range");
45  }
46  }
47  HEOS.unspecify_phase();
48  HEOS._Q = -1;
49  } else {
50  // Liquid solution
51  throw ValueError();
52  }
53  } else {
55  // Blind flash call
56  // Following the strategy of Gernert, 2014
57  StabilityRoutines::StabilityEvaluationClass stability_tester(HEOS);
58  if (!stability_tester.is_stable()) {
59  // There is a phase split and liquid and vapor phases are formed
61  stability_tester.get_liq(o.x, o.rhomolar_liq);
62  stability_tester.get_vap(o.y, o.rhomolar_vap);
63  o.z = HEOS.get_mole_fractions();
64  o.T = HEOS.T();
65  o.p = HEOS.p();
66  o.omega = 1.0;
68  solver.solve();
69  HEOS._phase = iphase_twophase;
70  HEOS._Q = (o.z[0] - o.x[0]) / (o.y[0] - o.x[0]); // All vapor qualities are the same (these are the residuals in the solver)
71  HEOS._rhomolar = 1 / (HEOS._Q / HEOS.SatV->rhomolar() + (1 - HEOS._Q) / HEOS.SatL->rhomolar());
72  } else {
73  // It's single-phase
74  double rho = HEOS.solver_rho_Tp_global(HEOS.T(), HEOS.p(), 20000);
75  HEOS.update_DmolarT_direct(rho, HEOS.T());
76  HEOS._Q = -1;
77  HEOS._phase = iphase_liquid;
78  }
79  } else {
80  // It's single-phase, and phase is imposed
81  double rho = HEOS.solver_rho_Tp(HEOS.T(), HEOS.p());
82  HEOS.update_DmolarT_direct(rho, HEOS.T());
83  HEOS._Q = -1;
84  HEOS._phase = HEOS.imposed_phase_index;
85  }
86  }
87 }
89  if (HEOS.is_pure_or_pseudopure) {
90  if (HEOS.imposed_phase_index == iphase_not_imposed) // If no phase index is imposed (see set_components function)
91  {
92  // At very low temperature (near the triple point temp), the isotherms are VERY steep
93  // Thus it can be very difficult to determine state based on ps = f(T)
94  // So in this case, we do a phase determination based on p, generally it will be useful enough
95  if (HEOS._T < 0.9 * HEOS.Ttriple() + 0.1 * HEOS.calc_Tmax_sat()) {
96  // Find the phase, while updating all internal variables possible using the pressure
97  bool saturation_called = false;
98  HEOS.p_phase_determination_pure_or_pseudopure(iT, HEOS._T, saturation_called);
99  } else {
100  // Find the phase, while updating all internal variables possible using the temperature
102  }
103  // Check if twophase solution
104  if (!HEOS.isHomogeneousPhase()) {
105  throw ValueError("twophase not implemented yet");
106  }
107  } else {
108  // Phase is imposed. Update _phase in case it was reset elsewhere by another call
109  HEOS._phase = HEOS.imposed_phase_index;
110  }
111  // Find density
112  HEOS._rhomolar = HEOS.solver_rho_Tp(HEOS._T, HEOS._p);
113  HEOS._Q = -1;
114  } else {
115  PT_flash_mixtures(HEOS);
116  }
117 }
118 
119 // Define the residual to be driven to zero
121 {
122  public:
126  double call(double T) {
128  CoolPropDbl peos = HEOS->p();
129  CoolPropDbl r = (peos - p) / p;
130  return r;
131  };
132  double deriv(double T) {
133  // dp/dT|rho / pspecified
134  return HEOS->first_partial_deriv(iP, iT, iDmolar) / p;
135  };
136  double second_deriv(double T) {
137  // d2p/dT2|rho / pspecified
138  return HEOS->second_partial_deriv(iP, iT, iDmolar, iT, iDmolar) / p;
139  };
140 };
141 
142 /***
143 \f[
144 \begin{array}{l}
145 p = \frac{{RT}}{{v - b}} - \frac{{a\alpha }}{{v\left( {v + b} \right)}}\\
146 \alpha = \left( {1 + \kappa \left( {1 - \sqrt {{T_r}} } \right)} \right)\left( {1 + \kappa \left( {1 - \sqrt {{T_r}} } \right)} \right) = 1 + 2\kappa \left( {1 - \sqrt {{T_r}} } \right) + {\kappa ^2}{\left( {1 - \sqrt {{T_r}} } \right)^2}\\
147 \alpha = 1 + 2\kappa \left( {1 - \sqrt {{T_r}} } \right) + {\kappa ^2}{\left( {1 - \sqrt {{T_r}} } \right)^2}\\
148 \alpha = 1 + 2\kappa - 2\kappa \sqrt {{T_r}} + {\kappa ^2}\left[ {1 - 2\sqrt {{T_r}} + {T_r}} \right]\\
149 T = {T_r}{T_c}\\
150 p = \frac{{R{T_r}{T_c}}}{{v - b}} - \frac{{a\left( {1 + 2\kappa - 2\kappa \sqrt {{T_r}} + {\kappa ^2}\left[ {1 - 2\sqrt {{T_r}} + {T_r}} \right]} \right)}}{{v\left( {v + b} \right)}}\\
151 \\
152 {\rm{Factor in terms of }}\sqrt {{T_r}} \\
153 \\
154 p = \frac{{R{T_r}{T_c}}}{{v - b}} - \frac{{a\left( {1 + 2\kappa + {\kappa ^2} - 2\kappa \sqrt {{T_r}} + {\kappa ^2}\left[ { - 2\sqrt {{T_r}} + {T_r}} \right]} \right)}}{{v\left( {v + b} \right)}}\\
155 p = \frac{{R{T_r}{T_c}}}{{v - b}} - \frac{{a\left( {1 + 2\kappa + {\kappa ^2} - 2\kappa (1 + \kappa )\sqrt {{T_r}} + {\kappa ^2}{T_r}} \right)}}{{v\left( {v + b} \right)}}\\
156 p = \frac{{R{T_r}{T_c}}}{{v - b}} - \frac{{a\left( {1 + 2\kappa + {\kappa ^2}} \right)}}{{v\left( {v + b} \right)}} + \frac{{2a\kappa (1 + \kappa )}}{{v\left( {v + b} \right)}}\sqrt {{T_r}} - \frac{{a{\kappa ^2}}}{{v\left( {v + b} \right)}}{T_r}\\
157 0 = \left[ {\frac{{R{T_c}}}{{v - b}} - \frac{{a{\kappa ^2}}}{{v\left( {v + b} \right)}}} \right]{T_r} + \frac{{2a\kappa (1 + \kappa )}}{{v\left( {v + b} \right)}}\sqrt {{T_r}} - \frac{{a\left( {1 + 2\kappa + {\kappa ^2}} \right)}}{{v\left( {v + b} \right)}} - p
158 \end{array}
159 \f]
160  */
161 double FlashRoutines::T_DP_PengRobinson(HelmholtzEOSMixtureBackend& HEOS, double rhomolar, double p) {
162  double omega, R, kappa, a, b, A, B, C, Tc, pc, V = 1 / rhomolar;
163  omega = HEOS.acentric_factor();
164  Tc = HEOS.T_critical();
165  pc = HEOS.p_critical();
166  R = HEOS.gas_constant();
167 
168  kappa = 0.37464 + 1.54226 * omega - 0.26992 * omega * omega;
169  a = 0.457235 * R * R * Tc * Tc / pc;
170  b = 0.077796 * R * Tc / pc;
171  double den = V * V + 2 * b * V - b * b;
172 
173  // A sqrt(Tr)^2 + B sqrt(Tr) + C = 0
174  A = R * Tc / (V - b) - a * kappa * kappa / (den);
175  B = +2 * a * kappa * (1 + kappa) / (den);
176  C = -a * (1 + 2 * kappa + kappa * kappa) / (den)-p;
177 
178  //D = B*B-4*A*C;
179 
180  double sqrt_Tr1 = (-B + sqrt(B * B - 4 * A * C)) / (2 * A);
181  //double sqrt_Tr2 = (-B-sqrt(B*B-4*A*C))/(2*A);
182  return sqrt_Tr1 * sqrt_Tr1 * Tc;
183 };
184 
186  // Comment out the check for an imposed phase. There's no code to handle if it is!
187  // Solver below and flash calculations (if two phase) have to be called anyway.
188  //
189  // if (HEOS.imposed_phase_index == iphase_not_imposed) // If no phase index is imposed (see set_components function)
190  // {
191  if (HEOS.is_pure_or_pseudopure) {
192  // Find the phase, while updating all internal variables possible using the pressure
193  bool saturation_called = false;
194  HEOS.p_phase_determination_pure_or_pseudopure(iDmolar, HEOS._rhomolar, saturation_called);
195 
196  if (HEOS.isHomogeneousPhase()) {
197  CoolPropDbl T0;
198  if (HEOS._phase == iphase_liquid) {
199  // If it is a liquid, start off at the ancillary value
200  if (saturation_called) {
201  T0 = HEOS.SatL->T();
202  } else {
203  T0 = HEOS._TLanc.pt();
204  }
205  } else if (HEOS._phase == iphase_supercritical_liquid) {
206  // If it is a supercritical
207  T0 = 1.1 * HEOS.T_critical();
208  } else if (HEOS._phase == iphase_gas || HEOS._phase == iphase_supercritical_gas || HEOS._phase == iphase_supercritical) {
209  // First, get a guess for density from Peng-Robinson
210  T0 = T_DP_PengRobinson(HEOS, HEOS.rhomolar(), HEOS.p());
211  } else {
212  throw ValueError("I should never get here");
213  }
214  // Then, do the solver using the full EOS
215  solver_DP_resid resid(&HEOS, HEOS.rhomolar(), HEOS.p());
216  std::string errstr;
217  Halley(resid, T0, 1e-10, 100);
218  HEOS._Q = -1;
219  // Update the state for conditions where the state was guessed
221  } else {
222  // Nothing to do here; phase determination has handled this already
223  }
224  } else {
225  throw NotImplementedError("DP_flash not ready for mixtures");
226  }
227  // }
228  // TO DO: Put the imposed phase check back in
229  // and provide the else code here if it is imposed.
230 }
231 
233 {
234  public:
238  double call(double T) {
239  HEOS.update(QT_INPUTS, 0, T); // Doesn't matter whether liquid or vapor, we are just doing a full VLE call for given T
243  return (1 / rhomolar - 1 / rhoL) / (1 / rhoV - 1 / rhoL) - Q_target;
244  }
245  double deriv(double T) {
246  return _HUGE;
247  }
248  double second_deriv(double T) {
249  return _HUGE;
250  }
251 };
252 
255  options.use_logdelta = false;
257  if (HEOS.is_pure_or_pseudopure) {
258  // Bump the temperatures to hopefully yield more reliable results
259  double Tmax = HEOS.T_critical() - 0.1;
260  double Tmin = HEOS.Tmin() + 0.1;
261  double rhomolar = HEOS._rhomolar;
262  double Q = HEOS._Q;
263  const double eps = 1e-12; // small tolerance to allow for slop in iterative calculations
264  if (rhomolar >= (HEOS.rhomolar_critical() + eps) && Q > (0 + eps)){
265  throw CoolProp::OutOfRangeError(format("DQ inputs are not defined for density (%g) above critical density (%g) and Q>0", rhomolar, HEOS.rhomolar_critical()).c_str());
266  }
267  DQ_flash_residual resid(HEOS, rhomolar, Q);
268  Brent(resid, Tmin, Tmax, DBL_EPSILON, 1e-10, 100);
269  HEOS._p = HEOS.SatV->p();
270  HEOS._T = HEOS.SatV->T();
271  HEOS._rhomolar = rhomolar;
272  HEOS._Q = Q;
273  HEOS._phase = iphase_twophase;
274  } else {
275  throw NotImplementedError("DQ_flash not ready for mixtures");
276  }
277 }
280  options.use_logdelta = false;
282  if (Tguess < 0) {
283  options.use_guesses = true;
284  options.T = Tguess;
285  CoolProp::SaturationAncillaryFunction& rhoL = HEOS.get_components()[0].ancillaries.rhoL;
286  CoolProp::SaturationAncillaryFunction& rhoV = HEOS.get_components()[0].ancillaries.rhoV;
287  options.rhoL = rhoL.evaluate(Tguess);
288  options.rhoV = rhoV.evaluate(Tguess);
289  }
290  if (HEOS.is_pure_or_pseudopure) {
291  if (std::abs(HEOS.Q() - 1) > 1e-10) {
292  throw ValueError(format("non-unity quality not currently allowed for HQ_flash"));
293  }
294  // Do a saturation call for given h for vapor, first with ancillaries, then with full saturation call
296  SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS.hmolar(), options);
297  HEOS._p = HEOS.SatV->p();
298  HEOS._T = HEOS.SatV->T();
299  HEOS._rhomolar = HEOS.SatV->rhomolar();
300  HEOS._phase = iphase_twophase;
301  } else {
302  throw NotImplementedError("HQ_flash not ready for mixtures");
303  }
304 }
306  if (HEOS.is_pure_or_pseudopure) {
307 
308  if (std::abs(HEOS.smolar() - HEOS.get_state("reducing").smolar) < 0.001) {
309  HEOS._p = HEOS.p_critical();
310  HEOS._T = HEOS.T_critical();
311  HEOS._rhomolar = HEOS.rhomolar_critical();
313  } else if (std::abs(HEOS.Q()) < 1e-10) {
314  // Do a saturation call for given s for liquid, first with ancillaries, then with full saturation call
317  options.use_logdelta = false;
319  SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS.smolar(), options);
320  HEOS._p = HEOS.SatL->p();
321  HEOS._T = HEOS.SatL->T();
322  HEOS._rhomolar = HEOS.SatL->rhomolar();
323  HEOS._phase = iphase_twophase;
324  } else if (std::abs(HEOS.Q() - 1) < 1e-10) {
325  // Do a saturation call for given s for vapor, first with ancillaries, then with full saturation call
328  options.use_logdelta = false;
330  SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS.smolar(), options);
331  HEOS._p = HEOS.SatV->p();
332  HEOS._T = HEOS.SatV->T();
333  HEOS._rhomolar = HEOS.SatV->rhomolar();
334  HEOS._phase = iphase_twophase;
335  } else {
336  throw ValueError(format("non-zero or 1 quality not currently allowed for QS_flash"));
337  }
338  } else {
339  throw NotImplementedError("QS_flash not ready for mixtures");
340  }
341 }
343  CoolPropDbl T = HEOS._T;
344  if (HEOS.is_pure_or_pseudopure) {
345  // The maximum possible saturation temperature
346  // Critical point for pure fluids, slightly different for pseudo-pure, very different for mixtures
347  CoolPropDbl Tmax_sat = HEOS.calc_Tmax_sat() + 1e-13;
348 
349  // Check what the minimum limits for the equation of state are
350  CoolPropDbl Tmin_satL, Tmin_satV, Tmin_sat;
351  HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
352  Tmin_sat = std::max(Tmin_satL, Tmin_satV) - 1e-13;
353 
354  // Get a reference to keep the code a bit cleaner
355  const CriticalRegionSplines& splines = HEOS.components[0].EOS().critical_region_splines;
356 
357  // If exactly(ish) at the critical temperature, liquid and vapor have the critial density
358  if ((get_config_bool(CRITICAL_WITHIN_1UK) && std::abs(T - Tmax_sat) < 1e-6) || std::abs(T - Tmax_sat) < 1e-12) {
359  HEOS.SatL->update(DmolarT_INPUTS, HEOS.rhomolar_critical(), HEOS._T);
360  HEOS.SatV->update(DmolarT_INPUTS, HEOS.rhomolar_critical(), HEOS._T);
361  HEOS._rhomolar = HEOS.rhomolar_critical();
362  HEOS._p = 0.5 * HEOS.SatV->p() + 0.5 * HEOS.SatL->p();
363  } else if (!is_in_closed_range(Tmin_sat - 0.1, Tmax_sat, T) && (CoolProp::get_config_bool(DONT_CHECK_PROPERTY_LIMITS) == false)) {
364  throw ValueError(format("Temperature to QT_flash [%0.8Lg K] must be in range [%0.8Lg K, %0.8Lg K]", T, Tmin_sat - 0.1, Tmax_sat));
365  } else if (get_config_bool(CRITICAL_SPLINES_ENABLED) && splines.enabled && HEOS._T > splines.T_min) {
366  double rhoL = _HUGE, rhoV = _HUGE;
367  // Use critical region spline if it has it and temperature is in its range
368  splines.get_densities(T, splines.rhomolar_min, HEOS.rhomolar_critical(), splines.rhomolar_max, rhoL, rhoV);
369  HEOS.SatL->update(DmolarT_INPUTS, rhoL, HEOS._T);
370  HEOS.SatV->update(DmolarT_INPUTS, rhoV, HEOS._T);
371  HEOS._p = 0.5 * HEOS.SatV->p() + 0.5 * HEOS.SatL->p();
372  HEOS._rhomolar = 1 / (HEOS._Q / HEOS.SatV->rhomolar() + (1 - HEOS._Q) / HEOS.SatL->rhomolar());
373  } else if (!(HEOS.components[0].EOS().pseudo_pure)) {
374  // Set some input options
376 
377  // Actually call the solver
379 
380  HEOS._p = 0.5 * HEOS.SatV->p() + 0.5 * HEOS.SatL->p();
381  HEOS._rhomolar = 1 / (HEOS._Q / HEOS.SatV->rhomolar() + (1 - HEOS._Q) / HEOS.SatL->rhomolar());
382  } else {
383  // Pseudo-pure fluid
384  CoolPropDbl rhoLanc = _HUGE, rhoVanc = _HUGE, rhoLsat = _HUGE, rhoVsat = _HUGE;
385  if (std::abs(HEOS._Q) < DBL_EPSILON) {
386  HEOS._p = HEOS.components[0].ancillaries.pL.evaluate(HEOS._T); // These ancillaries are used explicitly
387  rhoLanc = HEOS.components[0].ancillaries.rhoL.evaluate(HEOS._T);
388  HEOS.SatL->update_TP_guessrho(HEOS._T, HEOS._p, rhoLanc);
389  HEOS._rhomolar = HEOS.SatL->rhomolar();
390  } else if (std::abs(HEOS._Q - 1) < DBL_EPSILON) {
391  HEOS._p = HEOS.components[0].ancillaries.pV.evaluate(HEOS._T); // These ancillaries are used explicitly
392  rhoVanc = HEOS.components[0].ancillaries.rhoV.evaluate(HEOS._T);
393  HEOS.SatV->update_TP_guessrho(HEOS._T, HEOS._p, rhoVanc);
394  HEOS._rhomolar = HEOS.SatV->rhomolar();
395  } else {
396  throw CoolProp::ValueError(format("For pseudo-pure fluid, quality must be equal to 0 or 1. Two-phase quality is not defined"));
397  }
398 
399  try {
400  } catch (...) {
401  // Near the critical point, the behavior is not very nice, so we will just use the ancillary
402  rhoLsat = rhoLanc;
403  rhoVsat = rhoVanc;
404  }
405  }
406  // Load the outputs
407  HEOS._phase = iphase_twophase;
408  } else {
409  if (HEOS.PhaseEnvelope.built) {
410  PT_Q_flash_mixtures(HEOS, iT, HEOS._T);
411  } else {
412  // Set some input options
415  options.Nstep_max = 20;
416 
417  // Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted
419 
420  // Use Wilson iteration to obtain updated guess for pressure
422 
423  // Actually call the successive substitution solver
424  SaturationSolvers::successive_substitution(HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options);
425 
426  // -----
427  // Newton-Raphson
428  // -----
429 
432 
433  IO.bubble_point = (HEOS._Q < 0.5);
434 
435  IO.x = options.x;
436  IO.y = options.y;
437  IO.rhomolar_liq = options.rhomolar_liq;
438  IO.rhomolar_vap = options.rhomolar_vap;
439  IO.T = options.T;
440  IO.p = options.p;
441  IO.Nstep_max = 30;
442 
444 
445  if (IO.bubble_point) {
446  // Compositions are z, z_incipient
447  NR.call(HEOS, IO.x, IO.y, IO);
448  } else {
449  // Compositions are z, z_incipient
450  NR.call(HEOS, IO.y, IO.x, IO);
451  }
452 
453  HEOS._p = IO.p;
454  HEOS._rhomolar = 1 / (HEOS._Q / IO.rhomolar_vap + (1 - HEOS._Q) / IO.rhomolar_liq);
455  }
456  // Load the outputs
457  HEOS._phase = iphase_twophase;
458  HEOS._p = HEOS.SatV->p();
459  HEOS._rhomolar = 1 / (HEOS._Q / HEOS.SatV->rhomolar() + (1 - HEOS._Q) / HEOS.SatL->rhomolar());
460  HEOS._T = HEOS.SatL->T();
461  }
462 }
463 
464 void get_Henrys_coeffs_FP(const std::string& CAS, double& A, double& B, double& C, double& Tmin, double& Tmax) {
465  // Coeffs from Fernandez-Prini JPCRD 2003 DOI: 10.1063/1.1564818
466  if (CAS == "7440-59-7") //Helium
467  {
468  A = -3.52839;
469  B = 7.12983;
470  C = 4.47770;
471  Tmin = 273.21;
472  Tmax = 553.18;
473  } else if (CAS == "7440-01-9") // Ne
474  {
475  A = -3.18301;
476  B = 5.31448;
477  C = 5.43774;
478  Tmin = 273.20;
479  Tmax = 543.36;
480  } else if (CAS == "7440-37-1") // Ar
481  {
482  A = -8.40954;
483  B = 4.29587;
484  C = 10.52779;
485  Tmin = 273.19;
486  Tmax = 568.36;
487  } else if (CAS == "7439-90-9") // Kr
488  {
489  A = -8.97358;
490  B = 3.61508;
491  C = 11.29963;
492  Tmin = 273.19;
493  Tmax = 525.56;
494  } else if (CAS == "7440-63-3") // Xe
495  {
496  A = -14.21635;
497  B = 4.00041;
498  C = 15.60999;
499  Tmin = 273.22;
500  Tmax = 574.85;
501  } else if (CAS == "1333-74-0") // H2
502  {
503  A = -4.73284;
504  B = 6.08954;
505  C = 6.06066;
506  Tmin = 273.15;
507  Tmax = 636.09;
508  } else if (CAS == "7727-37-9") // N2
509  {
510  A = -9.67578;
511  B = 4.72162;
512  C = 11.70585;
513  Tmin = 278.12;
514  Tmax = 636.46;
515  } else if (CAS == "7782-44-7") // O2
516  {
517  A = -9.44833;
518  B = 4.43822;
519  C = 11.42005;
520  Tmin = 274.15;
521  Tmax = 616.52;
522  } else if (CAS == "630-08-0") // CO
523  {
524  A = -10.52862;
525  B = 5.13259;
526  C = 12.01421;
527  Tmin = 278.15;
528  Tmax = 588.67;
529  } else if (CAS == "124-38-9") // CO2
530  {
531  A = -8.55445;
532  B = 4.01195;
533  C = 9.52345;
534  Tmin = 274.19;
535  Tmax = 642.66;
536  } else if (CAS == "7783-06-4") // H2S
537  {
538  A = -4.51499;
539  B = 5.23538;
540  C = 4.42126;
541  Tmin = 273.15;
542  Tmax = 533.09;
543  } else if (CAS == "74-82-8") // CH4
544  {
545  A = -10.44708;
546  B = 4.66491;
547  C = 12.12986;
548  Tmin = 275.46;
549  Tmax = 633.11;
550  } else if (CAS == "74-84-0") // C2H6
551  {
552  A = -19.67563;
553  B = 4.51222;
554  C = 20.62567;
555  Tmin = 275.44;
556  Tmax = 473.46;
557  } else if (CAS == "2551-62-4") // SF6
558  {
559  A = -16.56118;
560  B = 2.15289;
561  C = 20.35440;
562  Tmin = 283.14;
563  Tmax = 505.55;
564  } else {
565  throw ValueError("Bad component in Henry's law constants");
566  }
567 }
569  if (HEOS.is_pure_or_pseudopure) {
570  if (HEOS.components[0].EOS().pseudo_pure) {
571  // It is a pseudo-pure mixture
572 
573  HEOS._TLanc = HEOS.components[0].ancillaries.pL.invert(HEOS._p);
574  HEOS._TVanc = HEOS.components[0].ancillaries.pV.invert(HEOS._p);
575  // Get guesses for the ancillaries for density
576  CoolPropDbl rhoL = HEOS.components[0].ancillaries.rhoL.evaluate(HEOS._TLanc);
577  CoolPropDbl rhoV = HEOS.components[0].ancillaries.rhoV.evaluate(HEOS._TVanc);
578  // Solve for the density
579  HEOS.SatL->update_TP_guessrho(HEOS._TLanc, HEOS._p, rhoL);
580  HEOS.SatV->update_TP_guessrho(HEOS._TVanc, HEOS._p, rhoV);
581 
582  // Load the outputs
583  HEOS._phase = iphase_twophase;
584  HEOS._p = HEOS._Q * HEOS.SatV->p() + (1 - HEOS._Q) * HEOS.SatL->p();
585  HEOS._T = HEOS._Q * HEOS.SatV->T() + (1 - HEOS._Q) * HEOS.SatL->T();
586  HEOS._rhomolar = 1 / (HEOS._Q / HEOS.SatV->rhomolar() + (1 - HEOS._Q) / HEOS.SatL->rhomolar());
587  } else {
588  // Critical point for pure fluids, slightly different for pseudo-pure, very different for mixtures
589  CoolPropDbl pmax_sat = HEOS.calc_pmax_sat();
590 
591  // Check what the minimum limits for the equation of state are
592  CoolPropDbl pmin_satL, pmin_satV, pmin_sat;
593  HEOS.calc_pmin_sat(pmin_satL, pmin_satV);
594  pmin_sat = std::max(pmin_satL, pmin_satV);
595 
596  // Check for being AT the critical point
597  if (is_in_closed_range(pmax_sat * (1 - 1e-10), pmax_sat * (1 + 1e-10), static_cast<CoolPropDbl>(HEOS._p))) {
598  // Load the outputs
600  HEOS._p = HEOS.p_critical();
601  HEOS._rhomolar = HEOS.rhomolar_critical();
602  HEOS._T = HEOS.T_critical();
603  return;
604  }
605 
606  // Check limits
607  if (CoolProp::get_config_bool(DONT_CHECK_PROPERTY_LIMITS) == false) {
608  if (!is_in_closed_range(pmin_sat * 0.999999, pmax_sat * 1.000001, static_cast<CoolPropDbl>(HEOS._p))) {
609  throw ValueError(format("Pressure to PQ_flash [%6g Pa] must be in range [%8Lg Pa, %8Lg Pa]", HEOS._p, pmin_sat, pmax_sat));
610  }
611  }
612  // ------------------
613  // It is a pure fluid
614  // ------------------
615 
616  // Set some input options
618  // Specified variable is pressure
620  // Use logarithm of delta as independent variables
621  options.use_logdelta = false;
622 
623  double increment = 0.4;
624 
625  try {
626  for (double omega = 1.0; omega > 0; omega -= increment) {
627  try {
628  options.omega = omega;
629 
630  // Actually call the solver
631  SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS._p, options);
632 
633  // If you get here, there was no error, all is well
634  break;
635  } catch (...) {
636  if (omega < 1.1 * increment) {
637  throw;
638  }
639  // else we are going to try again with a smaller omega
640  }
641  }
642  } catch (...) {
643  // We may need to polish the solution at low pressure
644  SaturationSolvers::saturation_P_pure_1D_T(HEOS, HEOS._p, options);
645  }
646 
647  // Load the outputs
648  HEOS._phase = iphase_twophase;
649  HEOS._p = HEOS._Q * HEOS.SatV->p() + (1 - HEOS._Q) * HEOS.SatL->p();
650  HEOS._rhomolar = 1 / (HEOS._Q / HEOS.SatV->rhomolar() + (1 - HEOS._Q) / HEOS.SatL->rhomolar());
651  HEOS._T = HEOS.SatL->T();
652  }
653  } else {
654  if (HEOS.PhaseEnvelope.built) {
655  PT_Q_flash_mixtures(HEOS, iP, HEOS._p);
656  } else {
657 
658  // Set some input options
661  io.Nstep_max = 10;
662 
663  // Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted
665 
666  // Use Wilson iteration to obtain updated guess for temperature
668 
669  std::vector<CoolPropDbl> K = HEOS.K;
670 
671  if (get_config_bool(HENRYS_LAW_TO_GENERATE_VLE_GUESSES) && std::abs(HEOS._Q - 1) < 1e-10) {
672  const std::vector<CoolPropFluid>& components = HEOS.get_components();
673  std::size_t iWater = 0;
674  double p1star = PropsSI("P", "T", Tguess, "Q", 1, "Water");
675  const std::vector<CoolPropDbl> y = HEOS.mole_fractions;
676  std::vector<CoolPropDbl> x(y.size());
677  for (std::size_t i = 0; i < components.size(); ++i) {
678  if (components[i].CAS == "7732-18-5") {
679  iWater = i;
680  continue;
681  } else {
682  double A, B, C, Tmin, Tmax;
683  get_Henrys_coeffs_FP(components[i].CAS, A, B, C, Tmin, Tmax);
684  double T_R = Tguess / 647.096, tau = 1 - T_R;
685  double k_H = p1star * exp(A / T_R + B * pow(tau, 0.355) / T_R + C * pow(T_R, -0.41) * exp(tau));
686  x[i] = y[i] * HEOS._p / k_H;
687  //
688  K[i] = y[i] / x[i];
689  }
690  }
691  // Update water K factor
692  double summer = 0;
693  for (std::size_t i = 0; i < y.size(); ++i) {
694  if (i != iWater) {
695  summer += x[i];
696  }
697  }
698  x[iWater] = summer;
699  K[iWater] = y[iWater] / x[iWater];
700  }
701 
702  // Actually call the successive substitution solver
703  SaturationSolvers::successive_substitution(HEOS, HEOS._Q, Tguess, HEOS._p, HEOS.mole_fractions, K, io);
704 
705  // -----
706  // Newton-Raphson
707  // -----
708 
711 
712  IO.bubble_point = (HEOS._Q < 0.5);
713  IO.x = io.x;
714  IO.y = io.y;
715  IO.rhomolar_liq = io.rhomolar_liq;
716  IO.rhomolar_vap = io.rhomolar_vap;
717  IO.T = io.T;
718  IO.p = io.p;
719  IO.Nstep_max = 30;
721 
722  if (IO.bubble_point) {
723  // Compositions are z, z_incipient
724  NR.call(HEOS, IO.x, IO.y, IO);
725  } else {
726  // Compositions are z, z_incipient
727  NR.call(HEOS, IO.y, IO.x, IO);
728  }
729  }
730 
731  // Load the outputs
732  HEOS._phase = iphase_twophase;
733  HEOS._p = HEOS.SatV->p();
734  HEOS._rhomolar = 1 / (HEOS._Q / HEOS.SatV->rhomolar() + (1 - HEOS._Q) / HEOS.SatL->rhomolar());
735  HEOS._T = HEOS.SatL->T();
736  }
737 }
738 
742  IO.rhomolar_liq = guess.rhomolar_liq;
743  IO.rhomolar_vap = guess.rhomolar_vap;
744  IO.x = std::vector<CoolPropDbl>(guess.x.begin(), guess.x.end());
745  IO.y = std::vector<CoolPropDbl>(guess.y.begin(), guess.y.end());
746  IO.T = guess.T;
747  IO.p = HEOS._p;
748  IO.bubble_point = false;
750 
751  if (std::abs(HEOS.Q()) < 1e-10) {
752  IO.bubble_point = true;
753  NR.call(HEOS, IO.x, IO.y, IO);
754  } else if (std::abs(HEOS.Q() - 1) < 1e-10) {
755  IO.bubble_point = false;
756  NR.call(HEOS, IO.y, IO.x, IO);
757  } else {
758  throw ValueError(format("Quality must be 0 or 1"));
759  }
760 
761  // Load the other outputs
762  HEOS._phase = iphase_twophase;
763  HEOS._rhomolar = 1 / (HEOS._Q / IO.rhomolar_vap + (1 - HEOS._Q) / IO.rhomolar_liq);
764  HEOS._T = IO.T;
765 }
769  IO.rhomolar_liq = guess.rhomolar_liq;
770  IO.rhomolar_vap = guess.rhomolar_vap;
771  IO.x = std::vector<CoolPropDbl>(guess.x.begin(), guess.x.end());
772  IO.y = std::vector<CoolPropDbl>(guess.y.begin(), guess.y.end());
773  IO.T = HEOS._T;
774  IO.p = guess.p;
775  IO.bubble_point = false;
777 
778  if (get_debug_level() > 9) {
779  std::cout << format(" QT w/ guess p %g T %g dl %g dv %g x %s y %s\n", IO.p, IO.T, IO.rhomolar_liq, IO.rhomolar_vap,
780  vec_to_string(IO.x, "%g").c_str(), vec_to_string(IO.y, "%g").c_str());
781  }
782 
783  if (std::abs(HEOS.Q()) < 1e-10) {
784  IO.bubble_point = true;
785  NR.call(HEOS, IO.x, IO.y, IO);
786  } else if (std::abs(HEOS.Q() - 1) < 1e-10) {
787  IO.bubble_point = false;
788  NR.call(HEOS, IO.y, IO.x, IO);
789  } else {
790  throw ValueError(format("Quality must be 0 or 1"));
791  }
792 
793  // Load the other outputs
794  HEOS._p = IO.p;
795  HEOS._phase = iphase_twophase;
796  HEOS._rhomolar = 1 / (HEOS._Q / IO.rhomolar_vap + (1 - HEOS._Q) / IO.rhomolar_liq);
797 }
798 
800  HEOS.solver_rho_Tp(HEOS.T(), HEOS.p(), guess.rhomolar);
801  // Load the other outputs
802  HEOS._phase = iphase_gas; // Guessed for mixtures
803  if (HEOS.is_pure_or_pseudopure) {
804  if (HEOS._p > HEOS.p_critical()) {
805  if (HEOS._T > HEOS.T_critical()) {
807  } else {
809  }
810  } else {
811  if (HEOS._T > HEOS.T_critical()) {
813  } else if (HEOS._rhomolar > HEOS.rhomolar_critical()) {
814  HEOS._phase = iphase_liquid;
815  } else {
816  HEOS._phase = iphase_gas;
817  }
818  }
819  }
820 
821  HEOS._Q = -1;
822 }
823 
825 
826  // Find the intersections in the phase envelope
827  std::vector<std::pair<std::size_t, std::size_t>> intersections =
829 
830  PhaseEnvelopeData& env = HEOS.PhaseEnvelope;
831 
832  enum quality_options
833  {
834  SATURATED_LIQUID,
835  SATURATED_VAPOR,
836  TWO_PHASE
837  };
838  quality_options quality;
839  if (std::abs(HEOS._Q) < 100 * DBL_EPSILON) {
840  quality = SATURATED_LIQUID;
841  } else if (std::abs(HEOS._Q - 1) < 100 * DBL_EPSILON) {
842  quality = SATURATED_VAPOR;
843  } else if (HEOS._Q > 0 && HEOS._Q < 1) {
844  quality = TWO_PHASE;
845  } else {
846  throw ValueError("Quality is not within 0 and 1");
847  }
848 
849  if (quality == SATURATED_LIQUID || quality == SATURATED_VAPOR) {
850  // *********************************************************
851  // Bubble- or dew-point calculation
852  // *********************************************************
853  // Find the correct solution
854  std::vector<std::size_t> solutions;
855  for (std::vector<std::pair<std::size_t, std::size_t>>::const_iterator it = intersections.begin(); it != intersections.end(); ++it) {
856  if (std::abs(env.Q[it->first] - HEOS._Q) < 10 * DBL_EPSILON && std::abs(env.Q[it->second] - HEOS._Q) < 10 * DBL_EPSILON) {
857  solutions.push_back(it->first);
858  }
859  }
860 
861  if (solutions.size() == 1) {
862 
863  std::size_t& imax = solutions[0];
864 
865  // Shift the solution if needed to ensure that imax+2 and imax-1 are both in range
866  if (imax + 2 >= env.T.size()) {
867  imax--;
868  } else if (imax == 0) {
869  imax++;
870  }
871  // Here imax+2 or imax-1 is still possibly out of range:
872  // 1. If imax initially is 1, and env.T.size() <= 3, then imax will become 0.
873  // 2. If imax initially is 0, and env.T.size() <= 2, then imax will become MAX_UINT.
874  // 3. If imax+2 initially is more than env.T.size(), then single decrement will not bring it to range
875 
878 
879  if (other == iP) {
880  IO.p = HEOS._p;
882  // p -> rhomolar_vap
883  IO.rhomolar_vap = CubicInterp(env.p, env.rhomolar_vap, imax - 1, imax, imax + 1, imax + 2, static_cast<CoolPropDbl>(IO.p));
884  IO.T = CubicInterp(env.rhomolar_vap, env.T, imax - 1, imax, imax + 1, imax + 2, IO.rhomolar_vap);
885  } else if (other == iT) {
886  IO.T = HEOS._T;
888  // T -> rhomolar_vap
889  IO.rhomolar_vap = CubicInterp(env.T, env.rhomolar_vap, imax - 1, imax, imax + 1, imax + 2, static_cast<CoolPropDbl>(IO.T));
890  IO.p = CubicInterp(env.rhomolar_vap, env.p, imax - 1, imax, imax + 1, imax + 2, IO.rhomolar_vap);
891  } else {
892  throw ValueError();
893  }
894  IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, imax - 1, imax, imax + 1, imax + 2, IO.rhomolar_vap);
895 
896  if (quality == SATURATED_VAPOR) {
897  IO.bubble_point = false;
898  IO.y = HEOS.get_mole_fractions(); // Because Q = 1
899  IO.x.resize(IO.y.size());
900  for (std::size_t i = 0; i < IO.x.size() - 1; ++i) // First N-1 elements
901  {
902  IO.x[i] = CubicInterp(env.rhomolar_vap, env.x[i], imax - 1, imax, imax + 1, imax + 2, IO.rhomolar_vap);
903  }
904  IO.x[IO.x.size() - 1] = 1 - std::accumulate(IO.x.begin(), IO.x.end() - 1, 0.0);
905  NR.call(HEOS, IO.y, IO.x, IO);
906  } else {
907  IO.bubble_point = true;
908  IO.x = HEOS.get_mole_fractions(); // Because Q = 0
909  IO.y.resize(IO.x.size());
910  // Phases are inverted, so "liquid" is actually the lighter phase
911  std::swap(IO.rhomolar_liq, IO.rhomolar_vap);
912  for (std::size_t i = 0; i < IO.y.size() - 1; ++i) // First N-1 elements
913  {
914  // Phases are inverted, so liquid mole fraction (x) of phase envelope is actually the vapor phase mole fraction
915  // Use the liquid density as well
916  IO.y[i] = CubicInterp(env.rhomolar_vap, env.x[i], imax - 1, imax, imax + 1, imax + 2, IO.rhomolar_liq);
917  }
918  IO.y[IO.y.size() - 1] = 1 - std::accumulate(IO.y.begin(), IO.y.end() - 1, 0.0);
919  NR.call(HEOS, IO.x, IO.y, IO);
920  }
921  } else if (solutions.size() == 0) {
922  throw ValueError("No solution was found in PQ_flash");
923  } else {
924  throw ValueError("More than 1 solution was found in PQ_flash");
925  }
926  } else {
927  // *********************************************************
928  // Two-phase calculation for given vapor quality
929  // *********************************************************
930 
931  // Find the correct solution
932  std::vector<std::size_t> liquid_solutions, vapor_solutions;
933  for (std::vector<std::pair<std::size_t, std::size_t>>::const_iterator it = intersections.begin(); it != intersections.end(); ++it) {
934  if (std::abs(env.Q[it->first] - 0) < 10 * DBL_EPSILON && std::abs(env.Q[it->second] - 0) < 10 * DBL_EPSILON) {
935  liquid_solutions.push_back(it->first);
936  }
937  if (std::abs(env.Q[it->first] - 1) < 10 * DBL_EPSILON && std::abs(env.Q[it->second] - 1) < 10 * DBL_EPSILON) {
938  vapor_solutions.push_back(it->first);
939  }
940  }
941 
942  if (liquid_solutions.size() != 1 || vapor_solutions.size() != 1) {
943  throw ValueError(format("Number liquid solutions [%d] or vapor solutions [%d] != 1", liquid_solutions.size(), vapor_solutions.size()));
944  }
945  std::size_t iliq = liquid_solutions[0], ivap = vapor_solutions[0];
946 
949  IO.beta = HEOS._Q;
950 
951  CoolPropDbl rhomolar_vap_sat_vap, T_sat_vap, rhomolar_liq_sat_vap, rhomolar_liq_sat_liq, T_sat_liq, rhomolar_vap_sat_liq, p_sat_liq,
952  p_sat_vap;
953 
954  if (other == iP) {
955  IO.p = HEOS._p;
956  p_sat_liq = IO.p;
957  p_sat_vap = IO.p;
959 
960  // Calculate the interpolated values for beta = 0 and beta = 1
961  rhomolar_vap_sat_vap = CubicInterp(env.p, env.rhomolar_vap, ivap - 1, ivap, ivap + 1, ivap + 2, static_cast<CoolPropDbl>(IO.p));
962  T_sat_vap = CubicInterp(env.rhomolar_vap, env.T, ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
963  rhomolar_liq_sat_vap = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
964 
965  // Phase inversion for liquid solution (liquid is vapor and vice versa)
966  rhomolar_liq_sat_liq = CubicInterp(env.p, env.rhomolar_vap, iliq - 1, iliq, iliq + 1, iliq + 2, static_cast<CoolPropDbl>(IO.p));
967  T_sat_liq = CubicInterp(env.rhomolar_vap, env.T, iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
968  rhomolar_vap_sat_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
969  } else if (other == iT) {
970  IO.T = HEOS._T;
971  T_sat_liq = IO.T;
972  T_sat_vap = IO.T;
974 
975  // Calculate the interpolated values for beta = 0 and beta = 1
976  rhomolar_vap_sat_vap = CubicInterp(env.T, env.rhomolar_vap, ivap - 1, ivap, ivap + 1, ivap + 2, static_cast<CoolPropDbl>(IO.T));
977  p_sat_vap = CubicInterp(env.rhomolar_vap, env.p, ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
978  rhomolar_liq_sat_vap = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
979 
980  // Phase inversion for liquid solution (liquid is vapor and vice versa)
981  rhomolar_liq_sat_liq = CubicInterp(env.T, env.rhomolar_vap, iliq - 1, iliq, iliq + 1, iliq + 2, static_cast<CoolPropDbl>(IO.T));
982  p_sat_liq = CubicInterp(env.rhomolar_vap, env.p, iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
983  rhomolar_vap_sat_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
984  } else {
985  throw ValueError();
986  }
987 
988  // Weight the guesses by the vapor mole fraction
989  IO.rhomolar_vap = IO.beta * rhomolar_vap_sat_vap + (1 - IO.beta) * rhomolar_vap_sat_liq;
990  IO.rhomolar_liq = IO.beta * rhomolar_liq_sat_vap + (1 - IO.beta) * rhomolar_liq_sat_liq;
991  IO.T = IO.beta * T_sat_vap + (1 - IO.beta) * T_sat_liq;
992  IO.p = IO.beta * p_sat_vap + (1 - IO.beta) * p_sat_liq;
993 
994  IO.z = HEOS.get_mole_fractions();
995  IO.x.resize(IO.z.size());
996  IO.y.resize(IO.z.size());
997 
998  for (std::size_t i = 0; i < IO.x.size() - 1; ++i) // First N-1 elements
999  {
1000  CoolPropDbl x_sat_vap = CubicInterp(env.rhomolar_vap, env.x[i], ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
1001  CoolPropDbl y_sat_vap = CubicInterp(env.rhomolar_vap, env.y[i], ivap - 1, ivap, ivap + 1, ivap + 2, rhomolar_vap_sat_vap);
1002 
1003  CoolPropDbl x_sat_liq = CubicInterp(env.rhomolar_vap, env.y[i], iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
1004  CoolPropDbl y_sat_liq = CubicInterp(env.rhomolar_vap, env.x[i], iliq - 1, iliq, iliq + 1, iliq + 2, rhomolar_liq_sat_liq);
1005 
1006  IO.x[i] = IO.beta * x_sat_vap + (1 - IO.beta) * x_sat_liq;
1007  IO.y[i] = IO.beta * y_sat_vap + (1 - IO.beta) * y_sat_liq;
1008  }
1009  IO.x[IO.x.size() - 1] = 1 - std::accumulate(IO.x.begin(), IO.x.end() - 1, 0.0);
1010  IO.y[IO.y.size() - 1] = 1 - std::accumulate(IO.y.begin(), IO.y.end() - 1, 0.0);
1011  NR.call(HEOS, IO);
1012  }
1013 }
1015  class Residual : public FuncWrapper1D
1016  {
1017 
1018  public:
1020  CoolPropDbl rhomolar_spec; // Specified value for density
1021  parameters other; // Key for other value
1022  CoolPropDbl value; // value for S,H,U
1023  CoolPropDbl Qd; // Quality from density
1024  Residual(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl rhomolar_spec, parameters other, CoolPropDbl value)
1025  : HEOS(HEOS), rhomolar_spec(rhomolar_spec), other(other), value(value) {
1026  Qd = _HUGE;
1027  };
1028  double call(double T) {
1029  HEOS.update(QT_INPUTS, 0, T);
1030  HelmholtzEOSMixtureBackend &SatL = HEOS.get_SatL(), &SatV = HEOS.get_SatV();
1031  // Quality from density
1032  Qd = (1 / rhomolar_spec - 1 / SatL.rhomolar()) / (1 / SatV.rhomolar() - 1 / SatL.rhomolar());
1033  // Quality from other parameter (H,S,U)
1034  CoolPropDbl Qo = (value - SatL.keyed_output(other)) / (SatV.keyed_output(other) - SatL.keyed_output(other));
1035  // Residual is the difference between the two
1036  return Qo - Qd;
1037  }
1038  } resid(HEOS, rhomolar_spec, other, value);
1039 
1040  // Critical point for pure fluids, slightly different for pseudo-pure, very different for mixtures
1041  CoolPropDbl Tmax_sat = HEOS.calc_Tmax_sat() - 1e-13;
1042 
1043  // Check what the minimum limits for the equation of state are
1044  CoolPropDbl Tmin_satL, Tmin_satV, Tmin_sat;
1045  HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
1046  Tmin_sat = std::max(Tmin_satL, Tmin_satV) - 1e-13;
1047 
1048  Brent(resid, Tmin_sat, Tmax_sat - 0.01, DBL_EPSILON, 1e-12, 20);
1049  // Solve once more with the final vapor quality
1050  HEOS.update(QT_INPUTS, resid.Qd, HEOS.T());
1051 }
1052 // D given and one of P,H,S,U
1054  // Define the residual to be driven to zero
1055  class solver_resid : public FuncWrapper1DWithTwoDerivs
1056  {
1057  public:
1059  CoolPropDbl rhomolar, value;
1060  parameters other;
1061  CoolPropDbl Tmin, Tmax;
1062 
1063  solver_resid(HelmholtzEOSMixtureBackend* HEOS, CoolPropDbl rhomolar, CoolPropDbl value, parameters other, CoolPropDbl Tmin, CoolPropDbl Tmax)
1064  : HEOS(HEOS), rhomolar(rhomolar), value(value), other(other), Tmin(Tmin), Tmax(Tmax) {
1066  HEOS->specify_phase(iphase_gas);
1067  };
1068  double call(double T) {
1069  HEOS->update_DmolarT_direct(rhomolar, T);
1070  double eos = HEOS->keyed_output(other);
1071  if (other == iP) {
1072  // For p, should use fractional error
1073  return (eos - value) / value;
1074  } else {
1075  // For everything else, use absolute error
1076  return eos - value;
1077  }
1078  };
1079  double deriv(double T) {
1080  if (other == iP) {
1081  return HEOS->first_partial_deriv(other, iT, iDmolar) / value;
1082  }
1083  return HEOS->first_partial_deriv(other, iT, iDmolar);
1084  };
1085  double second_deriv(double T) {
1086  if (other == iP) {
1087  return HEOS->second_partial_deriv(other, iT, iDmolar, iT, iDmolar) / value;
1088  }
1089  return HEOS->second_partial_deriv(other, iT, iDmolar, iT, iDmolar);
1090  };
1091  bool input_not_in_range(double T) {
1092  return (T < Tmin || T > Tmax);
1093  }
1094  };
1095 
1096  if (HEOS.is_pure_or_pseudopure) {
1097  CoolPropFluid& component = HEOS.components[0];
1098 
1099  shared_ptr<HelmholtzEOSMixtureBackend> Sat;
1100  CoolPropDbl rhoLtriple = component.triple_liquid.rhomolar;
1101  CoolPropDbl rhoVtriple = component.triple_vapor.rhomolar;
1102  // Check if in the "normal" region
1103  if (HEOS._rhomolar >= rhoVtriple && HEOS._rhomolar <= rhoLtriple) {
1104  CoolPropDbl yL, yV, value, y_solid;
1105  CoolPropDbl TLtriple = component.triple_liquid.T;
1106  CoolPropDbl TVtriple = component.triple_vapor.T;
1107 
1108  // First check if solid (below the line connecting the triple point values) - this is an error for now
1109  switch (other) {
1110  case iSmolar:
1111  yL = HEOS.calc_smolar_nocache(TLtriple, rhoLtriple);
1112  yV = HEOS.calc_smolar_nocache(TVtriple, rhoVtriple);
1113  value = HEOS._smolar;
1114  break;
1115  case iHmolar:
1116  yL = HEOS.calc_hmolar_nocache(TLtriple, rhoLtriple);
1117  yV = HEOS.calc_hmolar_nocache(TVtriple, rhoVtriple);
1118  value = HEOS._hmolar;
1119  break;
1120  case iUmolar:
1121  yL = HEOS.calc_umolar_nocache(TLtriple, rhoLtriple);
1122  yV = HEOS.calc_umolar_nocache(TVtriple, rhoVtriple);
1123  value = HEOS._umolar;
1124  break;
1125  case iP:
1126  yL = HEOS.calc_pressure_nocache(TLtriple, rhoLtriple);
1127  yV = HEOS.calc_pressure_nocache(TVtriple, rhoVtriple);
1128  value = HEOS._p;
1129  break;
1130  default:
1131  throw ValueError(format("Input is invalid"));
1132  }
1133  y_solid = (yV - yL) / (1 / rhoVtriple - 1 / rhoLtriple) * (1 / HEOS._rhomolar - 1 / rhoLtriple) + yL;
1134 
1135  if (value < y_solid) {
1136  throw ValueError(format("Other input [%d:%g] is solid", other, value));
1137  }
1138 
1139  // Check if other is above the saturation value.
1141  optionsD.omega = 1;
1142  optionsD.use_logdelta = false;
1143  optionsD.max_iterations = 200;
1144  for (int i_try = 0; i_try < 7; i_try++)
1145  {
1146  try
1147  {
1148  if (HEOS._rhomolar > HEOS._crit.rhomolar)
1149  {
1151  SaturationSolvers::saturation_D_pure(HEOS, HEOS._rhomolar, optionsD);
1152  // SatL and SatV have the saturation values
1153  Sat = HEOS.SatL;
1154  }
1155  else
1156  {
1158  SaturationSolvers::saturation_D_pure(HEOS, HEOS._rhomolar, optionsD);
1159  // SatL and SatV have the saturation values
1160  Sat = HEOS.SatV;
1161  }
1162  break; // good solve
1163  }
1164  catch(CoolPropBaseError)
1165  {
1166  optionsD.omega /= 2;
1167  optionsD.max_iterations *= 2;
1168  if (i_try >= 6){throw;}
1169  }
1170  }
1171 
1172  // If it is above, it is not two-phase and either liquid, vapor or supercritical
1173  if (value > Sat->keyed_output(other)) {
1174  solver_resid resid(&HEOS, HEOS._rhomolar, value, other, Sat->keyed_output(iT), HEOS.Tmax() * 1.5);
1175  try {
1176  HEOS._T = Halley(resid, 0.5 * (Sat->keyed_output(iT) + HEOS.Tmax() * 1.5), 1e-10, 100);
1177  } catch (...) {
1178  HEOS._T = Brent(resid, Sat->keyed_output(iT), HEOS.Tmax() * 1.5, DBL_EPSILON, 1e-12, 100);
1179  }
1180  HEOS._Q = 10000;
1181  HEOS._p = HEOS.calc_pressure_nocache(HEOS.T(), HEOS.rhomolar());
1182  HEOS.unspecify_phase();
1183  // Update the phase flag
1185  } else {
1186  // Now we know that temperature is between Tsat(D) +- tolerance and the minimum temperature for the fluid
1187  if (other == iP) {
1188  // Iterate to find T(p), its just a saturation call
1189 
1190  // Set some input options
1192  // Specified variable is pressure
1194  // Use logarithm of delta as independent variables
1195  optionsPHSU.use_logdelta = false;
1196 
1197  // Actually call the solver
1198  SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS._p, optionsPHSU);
1199 
1200  // Load the outputs
1201  HEOS._phase = iphase_twophase;
1202  HEOS._Q = (1 / HEOS._rhomolar - 1 / HEOS.SatL->rhomolar()) / (1 / HEOS.SatV->rhomolar() - 1 / HEOS.SatL->rhomolar());
1203  HEOS._T = HEOS.SatL->T();
1204  } else {
1205  // Residual is difference in quality calculated from density and quality calculated from the other parameter
1206  // Iterate to find T
1207  HSU_D_flash_twophase(HEOS, HEOS._rhomolar, other, value);
1208  HEOS._phase = iphase_twophase;
1209  }
1210  }
1211  }
1212  // Check if vapor/solid region below triple point vapor density
1213  else if (HEOS._rhomolar < component.triple_vapor.rhomolar) {
1214  CoolPropDbl y, value;
1215  CoolPropDbl TVtriple = component.triple_vapor.T; //TODO: separate TL and TV for ppure
1216 
1217  // If value is above the value calculated from X(Ttriple, _rhomolar), it is vapor
1218  switch (other) {
1219  case iSmolar:
1220  y = HEOS.calc_smolar_nocache(TVtriple, HEOS._rhomolar);
1221  value = HEOS._smolar;
1222  break;
1223  case iHmolar:
1224  y = HEOS.calc_hmolar_nocache(TVtriple, HEOS._rhomolar);
1225  value = HEOS._hmolar;
1226  break;
1227  case iUmolar:
1228  y = HEOS.calc_umolar_nocache(TVtriple, HEOS._rhomolar);
1229  value = HEOS._umolar;
1230  break;
1231  case iP:
1232  y = HEOS.calc_pressure_nocache(TVtriple, HEOS._rhomolar);
1233  value = HEOS._p;
1234  break;
1235  default:
1236  throw ValueError(format("Input is invalid"));
1237  }
1238  if (value > y) {
1239  solver_resid resid(&HEOS, HEOS._rhomolar, value, other, TVtriple, HEOS.Tmax() * 1.5);
1240  HEOS._phase = iphase_gas;
1241  try {
1242  HEOS._T = Halley(resid, 0.5 * (TVtriple + HEOS.Tmax() * 1.5), DBL_EPSILON, 100);
1243  } catch (...) {
1244  HEOS._T = Brent(resid, TVtriple, HEOS.Tmax() * 1.5, DBL_EPSILON, 1e-12, 100);
1245  }
1246  HEOS._Q = 10000;
1247  HEOS.calc_pressure();
1248  } else {
1249  throw ValueError(format("D < DLtriple %g %g", value, y));
1250  }
1251 
1252  }
1253  // Check in the liquid/solid region above the triple point density
1254  else {
1255  CoolPropDbl y, value;
1256  CoolPropDbl TLtriple = component.EOS().Ttriple;
1257 
1258  // If value is above the value calculated from X(Ttriple, _rhomolar), it is vapor
1259  switch (other) {
1260  case iSmolar:
1261  y = HEOS.calc_smolar_nocache(TLtriple, HEOS._rhomolar);
1262  value = HEOS._smolar;
1263  break;
1264  case iHmolar:
1265  y = HEOS.calc_hmolar_nocache(TLtriple, HEOS._rhomolar);
1266  value = HEOS._hmolar;
1267  break;
1268  case iUmolar:
1269  y = HEOS.calc_umolar_nocache(TLtriple, HEOS._rhomolar);
1270  value = HEOS._umolar;
1271  break;
1272  case iP:
1273  y = HEOS.calc_pressure_nocache(TLtriple, HEOS._rhomolar);
1274  value = HEOS._p;
1275  break;
1276  default:
1277  throw ValueError(format("Input is invalid"));
1278  }
1279  if (value > y) {
1280  solver_resid resid(&HEOS, HEOS._rhomolar, value, other, TLtriple, HEOS.Tmax() * 1.5);
1281  HEOS._phase = iphase_liquid;
1282  try {
1283  HEOS._T = Halley(resid, 0.5 * (TLtriple + HEOS.Tmax() * 1.5), DBL_EPSILON, 100);
1284  } catch (...) {
1285  HEOS._T = Brent(resid, TLtriple, HEOS.Tmax() * 1.5, DBL_EPSILON, 1e-12, 100);
1286  }
1287  HEOS._Q = 10000;
1288  HEOS.calc_pressure();
1289  } else {
1290  throw ValueError(format("D < DLtriple %g %g", value, y));
1291  }
1292  }
1293  // Update the state for conditions where the state was guessed
1294  if (HEOS.phase() != iphase_twophase) {
1296  }
1297  } else
1298  throw NotImplementedError("PHSU_D_flash not ready for mixtures");
1299 }
1300 
1302  double A[2][2], B[2][2];
1303  CoolPropDbl y = _HUGE;
1305  _HEOS.update(DmolarT_INPUTS, rhomolar0, T0);
1306  CoolPropDbl Tc = HEOS.calc_T_critical();
1307  CoolPropDbl rhoc = HEOS.calc_rhomolar_critical();
1308  CoolPropDbl R = HEOS.gas_constant();
1309  CoolPropDbl p = HEOS.p();
1310  switch (other) {
1311  case iHmolar:
1312  y = HEOS.hmolar();
1313  break;
1314  case iSmolar:
1315  y = HEOS.smolar();
1316  break;
1317  default:
1318  throw ValueError("other is invalid in HSU_P_flash_singlephase_Newton");
1319  }
1320 
1321  CoolPropDbl worst_error = 999;
1322  int iter = 0;
1323  bool failed = false;
1324  CoolPropDbl omega = 1.0, f2, df2_dtau, df2_ddelta;
1325  CoolPropDbl tau = _HEOS.tau(), delta = _HEOS.delta();
1326  while (worst_error > 1e-6 && failed == false) {
1327 
1328  // All the required partial derivatives
1329  CoolPropDbl a0 = _HEOS.calc_alpha0_deriv_nocache(0, 0, HEOS.mole_fractions, tau, delta, Tc, rhoc);
1330  CoolPropDbl da0_ddelta = _HEOS.calc_alpha0_deriv_nocache(0, 1, HEOS.mole_fractions, tau, delta, Tc, rhoc);
1331  CoolPropDbl da0_dtau = _HEOS.calc_alpha0_deriv_nocache(1, 0, HEOS.mole_fractions, tau, delta, Tc, rhoc);
1332  CoolPropDbl d2a0_dtau2 = _HEOS.calc_alpha0_deriv_nocache(2, 0, HEOS.mole_fractions, tau, delta, Tc, rhoc);
1333  CoolPropDbl d2a0_ddelta_dtau = 0.0;
1334 
1335  CoolPropDbl ar = _HEOS.calc_alphar_deriv_nocache(0, 0, HEOS.mole_fractions, tau, delta);
1336  CoolPropDbl dar_dtau = _HEOS.calc_alphar_deriv_nocache(1, 0, HEOS.mole_fractions, tau, delta);
1337  CoolPropDbl dar_ddelta = _HEOS.calc_alphar_deriv_nocache(0, 1, HEOS.mole_fractions, tau, delta);
1338  CoolPropDbl d2ar_ddelta_dtau = _HEOS.calc_alphar_deriv_nocache(1, 1, HEOS.mole_fractions, tau, delta);
1339  CoolPropDbl d2ar_ddelta2 = _HEOS.calc_alphar_deriv_nocache(0, 2, HEOS.mole_fractions, tau, delta);
1340  CoolPropDbl d2ar_dtau2 = _HEOS.calc_alphar_deriv_nocache(2, 0, HEOS.mole_fractions, tau, delta);
1341 
1342  CoolPropDbl f1 = delta / tau * (1 + delta * dar_ddelta) - p / (rhoc * R * Tc);
1343  CoolPropDbl df1_dtau = (1 + delta * dar_ddelta) * (-delta / tau / tau) + delta / tau * (delta * d2ar_ddelta_dtau);
1344  CoolPropDbl df1_ddelta = (1.0 / tau) * (1 + 2.0 * delta * dar_ddelta + delta * delta * d2ar_ddelta2);
1345  switch (other) {
1346  case iHmolar: {
1347  f2 = (1 + delta * dar_ddelta) + tau * (da0_dtau + dar_dtau) - tau * y / (R * Tc);
1348  df2_dtau = delta * d2ar_ddelta_dtau + da0_dtau + dar_dtau + tau * (d2a0_dtau2 + d2ar_dtau2) - y / (R * Tc);
1349  df2_ddelta = (dar_ddelta + delta * d2ar_ddelta2) + tau * (d2a0_ddelta_dtau + d2ar_ddelta_dtau);
1350  break;
1351  }
1352  case iSmolar: {
1353  f2 = tau * (da0_dtau + dar_dtau) - ar - a0 - y / R;
1354  df2_dtau = tau * (d2a0_dtau2 + d2ar_dtau2) + (da0_dtau + dar_dtau) - dar_dtau - da0_dtau;
1355  df2_ddelta = tau * (d2a0_ddelta_dtau + d2ar_ddelta_dtau) - dar_ddelta - da0_ddelta;
1356  break;
1357  }
1358  default:
1359  throw ValueError("other variable in HSU_P_flash_singlephase_Newton is invalid");
1360  }
1361 
1362  //First index is the row, second index is the column
1363  A[0][0] = df1_dtau;
1364  A[0][1] = df1_ddelta;
1365  A[1][0] = df2_dtau;
1366  A[1][1] = df2_ddelta;
1367 
1368  //double det = A[0][0]*A[1][1]-A[1][0]*A[0][1];
1369 
1370  MatInv_2(A, B);
1371  tau -= omega * (B[0][0] * f1 + B[0][1] * f2);
1372  delta -= omega * (B[1][0] * f1 + B[1][1] * f2);
1373 
1374  if (std::abs(f1) > std::abs(f2))
1375  worst_error = std::abs(f1);
1376  else
1377  worst_error = std::abs(f2);
1378 
1379  if (!ValidNumber(f1) || !ValidNumber(f2)) {
1380  throw SolutionError(format("Invalid values for inputs p=%g y=%g for fluid %s in HSU_P_flash_singlephase", p, y, _HEOS.name().c_str()));
1381  }
1382 
1383  iter += 1;
1384  if (iter > 100) {
1385  throw SolutionError(format("HSU_P_flash_singlephase did not converge with inputs p=%g h=%g for fluid %s", p, y, _HEOS.name().c_str()));
1386  }
1387  }
1388 
1389  HEOS.update(DmolarT_INPUTS, rhoc * delta, Tc / tau);
1390 }
1392  CoolPropDbl Tmax, phases phase) {
1393  if (!ValidNumber(HEOS._p)) {
1394  throw ValueError("value for p in HSU_P_flash_singlephase_Brent is invalid");
1395  };
1396  if (!ValidNumber(value)) {
1397  throw ValueError("value for other in HSU_P_flash_singlephase_Brent is invalid");
1398  };
1399  class solver_resid : public FuncWrapper1DWithTwoDerivs
1400  {
1401  public:
1403  CoolPropDbl p, value;
1404  parameters other;
1405  int iter;
1406  CoolPropDbl eos0, eos1, rhomolar, rhomolar0, rhomolar1;
1407  CoolPropDbl Tmin, Tmax;
1408 
1409  solver_resid(HelmholtzEOSMixtureBackend* HEOS, CoolPropDbl p, CoolPropDbl value, parameters other, double Tmin, double Tmax)
1410  : HEOS(HEOS),
1411  p(p),
1412  value(value),
1413  other(other),
1414  iter(0),
1415  eos0(-_HUGE),
1416  eos1(-_HUGE),
1417  rhomolar(_HUGE),
1418  rhomolar0(_HUGE),
1419  rhomolar1(_HUGE),
1420  Tmin(Tmin),
1421  Tmax(Tmax) {
1422  // Specify the state to avoid saturation calls, but only if phase is subcritical
1423  switch (CoolProp::phases phase = HEOS->phase()) {
1424  case iphase_liquid:
1425  case iphase_gas:
1426  HEOS->specify_phase(phase);
1427  default:
1428  // Otherwise don't do anything (this is to make compiler happy)
1429  {}
1430  }
1431  }
1432  double call(double T) {
1433 
1434  if (iter < 2 || std::abs(rhomolar1 / rhomolar0 - 1) > 0.05) {
1435  // Run the solver with T,P as inputs; but only if the last change in density was greater than a few percent
1436  HEOS->update(PT_INPUTS, p, T);
1437  } else {
1438  // Run the solver with T,P as inputs; but use the guess value for density from before
1439  HEOS->update_TP_guessrho(T, p, rhomolar);
1440  }
1441 
1442  // Get the value of the desired variable
1443  CoolPropDbl eos = HEOS->keyed_output(other);
1444 
1445  // Store the value of density
1446  rhomolar = HEOS->rhomolar();
1447 
1448  // Difference between the two is to be driven to zero
1449  CoolPropDbl r = eos - value;
1450 
1451  // Store values for later use if there are errors
1452  if (iter == 0) {
1453  eos0 = eos;
1454  rhomolar0 = rhomolar;
1455  } else if (iter == 1) {
1456  eos1 = eos;
1457  rhomolar1 = rhomolar;
1458  } else {
1459  eos0 = eos1;
1460  eos1 = eos;
1461  rhomolar0 = rhomolar1;
1462  rhomolar1 = rhomolar;
1463  }
1464 
1465  iter++;
1466  return r;
1467  };
1468  double deriv(double T) {
1469  return HEOS->first_partial_deriv(other, iT, iP);
1470  }
1471  double second_deriv(double T) {
1472  return HEOS->second_partial_deriv(other, iT, iP, iT, iP);
1473  }
1474  bool input_not_in_range(double x) {
1475  return (x < Tmin || x > Tmax);
1476  }
1477  };
1478  solver_resid resid(&HEOS, HEOS._p, value, other, Tmin, Tmax);
1479 
1480  try {
1481  // First try to use Halley's method (including two derivatives)
1482  Halley(resid, Tmin, 1e-12, 100);
1483  if (!is_in_closed_range(Tmin, Tmax, static_cast<CoolPropDbl>(resid.HEOS->T())) || resid.HEOS->phase() != phase) {
1484  throw ValueError("Halley's method was unable to find a solution in HSU_P_flash_singlephase_Brent");
1485  }
1486  // Un-specify the phase of the fluid
1487  HEOS.unspecify_phase();
1488  } catch (...) {
1489  try {
1490  resid.iter = 0;
1491  // Halley's method failed, so now we try Brent's method
1492  Brent(resid, Tmin, Tmax, DBL_EPSILON, 1e-12, 100);
1493  // Un-specify the phase of the fluid
1494  HEOS.unspecify_phase();
1495  } catch (...) {
1496  // Un-specify the phase of the fluid
1497  HEOS.unspecify_phase();
1498 
1499  // Determine why you were out of range if you can
1500  //
1501  CoolPropDbl eos0 = resid.eos0, eos1 = resid.eos1;
1502  std::string name = get_parameter_information(other, "short");
1503  std::string units = get_parameter_information(other, "units");
1504  if (eos1 > eos0 && value > eos1) {
1505  throw ValueError(
1506  format("HSU_P_flash_singlephase_Brent could not find a solution because %s [%Lg %s] is above the maximum value of %0.12Lg %s",
1507  name.c_str(), value, units.c_str(), eos1, units.c_str()));
1508  }
1509  if (eos1 > eos0 && value < eos0) {
1510  throw ValueError(
1511  format("HSU_P_flash_singlephase_Brent could not find a solution because %s [%Lg %s] is below the minimum value of %0.12Lg %s",
1512  name.c_str(), value, units.c_str(), eos0, units.c_str()));
1513  }
1514  throw;
1515  }
1516  }
1517 }
1518 
1519 // P given and one of H, S, or U
1521  bool saturation_called = false;
1522  CoolPropDbl value;
1523 
1524  // Find the phase, while updating all internal variables possible
1525  switch (other) {
1526  case iSmolar:
1527  value = HEOS.smolar();
1528  break;
1529  case iHmolar:
1530  value = HEOS.hmolar();
1531  break;
1532  case iUmolar:
1533  value = HEOS.umolar();
1534  break;
1535  default:
1536  throw ValueError(format("Input for other [%s] is invalid", get_parameter_information(other, "long").c_str()));
1537  }
1538  if (HEOS.is_pure_or_pseudopure) {
1539 
1540  // Find the phase, while updating all internal variables possible
1541  HEOS.p_phase_determination_pure_or_pseudopure(other, value, saturation_called);
1542 
1543  if (HEOS.isHomogeneousPhase()) {
1544  // Now we use the single-phase solver to find T,rho given P,Y using a
1545  // bounded 1D solver by adjusting T and using given value of p
1546  CoolPropDbl Tmin, Tmax;
1547  switch (HEOS._phase) {
1548  case iphase_gas: {
1549  Tmax = 1.5 * HEOS.Tmax();
1550  if (HEOS._p < HEOS.p_triple()) {
1551  Tmin = std::max(HEOS.Tmin(), HEOS.Ttriple());
1552  } else {
1553  if (saturation_called) {
1554  Tmin = HEOS.SatV->T();
1555  } else {
1556  Tmin = HEOS._TVanc.pt() - 0.01;
1557  }
1558  }
1559  break;
1560  }
1561  case iphase_liquid: {
1562  if (saturation_called) {
1563  Tmax = HEOS.SatL->T();
1564  } else {
1565  Tmax = HEOS._TLanc.pt() + 0.01;
1566  }
1567 
1568  // Sometimes the minimum pressure for the melting line is a bit above the triple point pressure
1569  if (HEOS.has_melting_line() && HEOS._p > HEOS.calc_melting_line(iP_min, -1, -1)) {
1570  Tmin = HEOS.calc_melting_line(iT, iP, HEOS._p) - 1e-3;
1571  } else {
1572  Tmin = HEOS.Tmin() - 1e-3;
1573  }
1574  break;
1575  }
1578  case iphase_supercritical: {
1579  Tmax = 1.5 * HEOS.Tmax();
1580  // Sometimes the minimum pressure for the melting line is a bit above the triple point pressure
1581  if (HEOS.has_melting_line() && HEOS._p > HEOS.calc_melting_line(iP_min, -1, -1)) {
1582  Tmin = HEOS.calc_melting_line(iT, iP, HEOS._p) - 1e-3;
1583  } else {
1584  Tmin = HEOS.Tmin() - 1e-3;
1585  }
1586  break;
1587  }
1588  default: {
1589  throw ValueError(format("Not a valid homogeneous state"));
1590  }
1591  }
1592  try {
1593  HSU_P_flash_singlephase_Brent(HEOS, other, value, Tmin, Tmax, HEOS._phase);
1594  } catch (std::exception& e) {
1595  throw ValueError(format("unable to solve 1phase PY flash with Tmin=%Lg, Tmax=%Lg due to error: %s", Tmin, Tmax, e.what()));
1596  }
1597  HEOS._Q = -1;
1598  // Update the state for conditions where the state was guessed
1600  }
1601  } else {
1602  if (HEOS.PhaseEnvelope.built) {
1603  // Determine whether you are inside or outside
1604  SimpleState closest_state;
1605  std::size_t iclosest;
1606  bool twophase = PhaseEnvelopeRoutines::is_inside(HEOS.PhaseEnvelope, iP, HEOS._p, other, value, iclosest, closest_state);
1607 
1608  if (!twophase) {
1609  PY_singlephase_flash_resid resid(HEOS, HEOS._p, other, value);
1610  // If that fails, try a bounded solver
1611  Brent(resid, closest_state.T + 10, 1000, DBL_EPSILON, 1e-10, 100);
1612  HEOS.unspecify_phase();
1613  } else {
1614  throw ValueError("two-phase solution for Y");
1615  }
1616 
1617  } else {
1618  throw ValueError("phase envelope must be built to carry out HSU_P_flash for mixture");
1619  }
1620  }
1621 }
1623  // Define the residual to be driven to zero
1624  class solver_resid : public FuncWrapper1DWithTwoDerivs
1625  {
1626  public:
1628  CoolPropDbl T, value;
1629  parameters other;
1630 
1631  solver_resid(HelmholtzEOSMixtureBackend* HEOS, CoolPropDbl T, CoolPropDbl value, parameters other)
1632  : HEOS(HEOS), T(T), value(value), other(other) {}
1633  double call(double rhomolar) {
1634  HEOS->update_DmolarT_direct(rhomolar, T);
1635  double eos = HEOS->keyed_output(other);
1636  return eos - value;
1637  };
1638  double deriv(double rhomolar) {
1639  return HEOS->first_partial_deriv(other, iDmolar, iT);
1640  }
1641  double second_deriv(double rhomolar) {
1642  return HEOS->second_partial_deriv(other, iDmolar, iT, iDmolar, iT);
1643  }
1644  };
1645  solver_resid resid(&HEOS, T, value, other);
1646 
1647  // Supercritical temperature
1648  if (HEOS._T > HEOS._crit.T) {
1649  CoolPropDbl yc, ymin, y;
1650  CoolPropDbl rhoc = HEOS.components[0].crit.rhomolar;
1651  CoolPropDbl rhomin = 1e-10;
1652 
1653  // Determine limits for the other variable
1654  switch (other) {
1655  case iSmolar: {
1656  yc = HEOS.calc_smolar_nocache(HEOS._T, rhoc);
1657  ymin = HEOS.calc_smolar_nocache(HEOS._T, rhomin);
1658  y = HEOS._smolar;
1659  break;
1660  }
1661  case iHmolar: {
1662  yc = HEOS.calc_hmolar_nocache(HEOS._T, rhoc);
1663  ymin = HEOS.calc_hmolar_nocache(HEOS._T, rhomin);
1664  y = HEOS._hmolar;
1665  break;
1666  }
1667  case iUmolar: {
1668  yc = HEOS.calc_umolar_nocache(HEOS._T, rhoc);
1669  ymin = HEOS.calc_umolar_nocache(HEOS._T, rhomin);
1670  y = HEOS._umolar;
1671  break;
1672  }
1673  default:
1674  throw ValueError();
1675  }
1676  if (is_in_closed_range(yc, ymin, y)) {
1677  Brent(resid, rhoc, rhomin, LDBL_EPSILON, 1e-9, 100);
1678  } else if (y < yc) {
1679  // Increase rhomelt until it bounds the solution
1680  int step_count = 0;
1681  while (!is_in_closed_range(ymin, yc, y)) {
1682  rhoc *= 1.1; // Increase density by a few percent
1683  switch (other) {
1684  case iSmolar:
1685  yc = HEOS.calc_smolar_nocache(HEOS._T, rhoc);
1686  break;
1687  case iHmolar:
1688  yc = HEOS.calc_hmolar_nocache(HEOS._T, rhoc);
1689  break;
1690  case iUmolar:
1691  yc = HEOS.calc_umolar_nocache(HEOS._T, rhoc);
1692  break;
1693  default:
1694  throw ValueError(format("Input is invalid"));
1695  }
1696  if (step_count > 30) {
1697  throw ValueError(format("Even by increasing rhoc, not able to bound input; input %Lg is not in range %Lg,%Lg", y, yc, ymin));
1698  }
1699  step_count++;
1700  }
1701  Brent(resid, rhomin, rhoc, LDBL_EPSILON, 1e-9, 100);
1702  } else {
1703  throw ValueError(format("input %Lg is not in range %Lg,%Lg,%Lg", y, yc, ymin));
1704  }
1705  // Update the state (T > Tc)
1706  if (HEOS._p < HEOS.p_critical()) {
1708  } else {
1710  }
1711  }
1712  // Subcritical temperature liquid
1713  else if ((HEOS._phase == iphase_liquid) || (HEOS._phase == iphase_supercritical_liquid)) {
1714  CoolPropDbl ymelt, yL, y;
1715  CoolPropDbl rhomelt = HEOS.components[0].triple_liquid.rhomolar;
1716  CoolPropDbl rhoL = static_cast<double>(HEOS._rhoLanc);
1717 
1718  switch (other) {
1719  case iSmolar: {
1720  ymelt = HEOS.calc_smolar_nocache(HEOS._T, rhomelt);
1721  yL = HEOS.calc_smolar_nocache(HEOS._T, rhoL);
1722  y = HEOS._smolar;
1723  break;
1724  }
1725  case iHmolar: {
1726  ymelt = HEOS.calc_hmolar_nocache(HEOS._T, rhomelt);
1727  yL = HEOS.calc_hmolar_nocache(HEOS._T, rhoL);
1728  y = HEOS._hmolar;
1729  break;
1730  }
1731  case iUmolar: {
1732  ymelt = HEOS.calc_umolar_nocache(HEOS._T, rhomelt);
1733  yL = HEOS.calc_umolar_nocache(HEOS._T, rhoL);
1734  y = HEOS._umolar;
1735  break;
1736  }
1737  default:
1738  throw ValueError();
1739  }
1740 
1741  CoolPropDbl rhomolar_guess = (rhomelt - rhoL) / (ymelt - yL) * (y - yL) + rhoL;
1742 
1743  try {
1744  Halley(resid, rhomolar_guess, 1e-8, 100);
1745  } catch (...) {
1746  Secant(resid, rhomolar_guess, 0.0001 * rhomolar_guess, 1e-12, 100);
1747  }
1748  }
1749  // Subcritical temperature gas
1750  else if (HEOS._phase == iphase_gas) {
1751  CoolPropDbl rhomin = 1e-14;
1752  CoolPropDbl rhoV = static_cast<double>(HEOS._rhoVanc);
1753 
1754  try {
1755  Halley(resid, 0.5 * (rhomin + rhoV), 1e-8, 100);
1756  } catch (...) {
1757  try {
1758  Brent(resid, rhomin, rhoV, LDBL_EPSILON, 1e-12, 100);
1759  } catch (...) {
1760  throw ValueError();
1761  }
1762  }
1763  } else {
1764  throw ValueError(format("phase to solver_for_rho_given_T_oneof_HSU is invalid"));
1765  }
1766 };
1767 
1770  // Use the phase defined by the imposed phase
1771  HEOS._phase = HEOS.imposed_phase_index;
1772  // The remaining code in this branch was added to set some needed parameters if phase is imposed,
1773  // since HEOS.T_phase_determination_pure_or_pseudopure() is not being called.
1774  if (HEOS._T < HEOS._crit.T) //
1775  {
1776  HEOS._rhoVanc = HEOS.components[0].ancillaries.rhoV.evaluate(HEOS._T);
1777  HEOS._rhoLanc = HEOS.components[0].ancillaries.rhoL.evaluate(HEOS._T);
1778  if (HEOS._phase == iphase_liquid) {
1779  HEOS._Q = -1000;
1780  } else if (HEOS._phase == iphase_gas) {
1781  HEOS._Q = 1000;
1782  } else if (HEOS._phase == iphase_twophase) {
1783  // Actually have to use saturation information sadly
1784  // For the given temperature, find the saturation state
1785  // Run the saturation routines to determine the saturation densities and pressures
1788  SaturationSolvers::saturation_T_pure(HEOS1, HEOS._T, options);
1789 
1790  if (other != iDmolar) {
1791  // Update the states
1792  if (HEOS.SatL) HEOS.SatL->update(DmolarT_INPUTS, HEOS._rhoLanc, HEOS._T);
1793  if (HEOS.SatV) HEOS.SatV->update(DmolarT_INPUTS, HEOS._rhoVanc, HEOS._T);
1794  // Update the two-Phase variables
1795  HEOS._rhoLmolar = HEOS.SatL->rhomolar();
1796  HEOS._rhoVmolar = HEOS.SatV->rhomolar();
1797  }
1798 
1799  CoolPropDbl Q;
1800 
1801  switch (other) {
1802  case iDmolar:
1803  Q = (1 / HEOS.rhomolar() - 1 / HEOS1.SatL->rhomolar()) / (1 / HEOS1.SatV->rhomolar() - 1 / HEOS1.SatL->rhomolar());
1804  break;
1805  case iSmolar:
1806  Q = (HEOS.smolar() - HEOS1.SatL->smolar()) / (HEOS1.SatV->smolar() - HEOS1.SatL->smolar());
1807  break;
1808  case iHmolar:
1809  Q = (HEOS.hmolar() - HEOS1.SatL->hmolar()) / (HEOS1.SatV->hmolar() - HEOS1.SatL->hmolar());
1810  break;
1811  case iUmolar:
1812  Q = (HEOS.umolar() - HEOS1.SatL->umolar()) / (HEOS1.SatV->umolar() - HEOS1.SatL->umolar());
1813  break;
1814  default:
1815  throw ValueError(format("bad input for other"));
1816  }
1817  if (Q < 0) {
1818  HEOS._Q = -1;
1819  } else if (Q > 1) {
1820  HEOS._Q = 1;
1821  } else {
1822  HEOS._Q = Q;
1823  // Load the outputs
1824  HEOS._p = HEOS._Q * HEOS1.SatV->p() + (1 - HEOS._Q) * HEOS1.SatL->p();
1825  HEOS._rhomolar = 1 / (HEOS._Q / HEOS.SatV->rhomolar() + (1 - HEOS._Q) / HEOS.SatL->rhomolar());
1826  }
1827  } else if (HEOS._phase == iphase_supercritical_liquid) {
1828  HEOS._Q = -1000;
1829  } else
1830  throw ValueError(format("Temperature specified is not the imposed phase region."));
1831  } else if (HEOS._T > HEOS._crit.T && HEOS._T > HEOS.components[0].EOS().Ttriple) {
1832  HEOS._Q = 1e9;
1833  }
1834  } else {
1835  if (HEOS.is_pure_or_pseudopure) {
1836  // Find the phase, while updating all internal variables possible
1837  switch (other) {
1838  case iDmolar:
1840  break;
1841  case iSmolar:
1843  break;
1844  case iHmolar:
1846  break;
1847  case iUmolar:
1849  break;
1850  default:
1851  throw ValueError(format("Input is invalid"));
1852  }
1853  } else {
1854  HEOS._phase = iphase_gas;
1855  throw NotImplementedError("DHSU_T_flash does not support mixtures (yet)");
1856  }
1857  }
1858 
1859  //if (HEOS.isHomogeneousPhase() && !ValidNumber(HEOS._p)) // original, pre 1352
1860  // only the solver requires single phase
1861  if (((other == iDmolar) || HEOS.isHomogeneousPhase()) && !ValidNumber(HEOS._p)) // post 1352
1862  {
1863  switch (other) {
1864  case iDmolar:
1865  break;
1866  case iHmolar:
1868  break;
1869  case iSmolar:
1871  break;
1872  case iUmolar:
1874  break;
1875  default:
1876  break;
1877  }
1878  HEOS.calc_pressure();
1879  HEOS._Q = -1;
1880  }
1881  if (HEOS.is_pure_or_pseudopure && HEOS._phase != iphase_twophase) {
1882  // Update the state for conditions where the state was guessed
1884  }
1885 }
1887  HS_flash_twophaseOptions& options) {
1888  class Residual : public FuncWrapper1D
1889  {
1890 
1891  public:
1893  CoolPropDbl hmolar, smolar, Qs;
1894  Residual(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl hmolar_spec, CoolPropDbl smolar_spec)
1895  : HEOS(HEOS), hmolar(hmolar_spec), smolar(smolar_spec), Qs(_HUGE){};
1896  double call(double T) {
1897  HEOS.update(QT_INPUTS, 0, T);
1898  HelmholtzEOSMixtureBackend &SatL = HEOS.get_SatL(), &SatV = HEOS.get_SatV();
1899  // Quality from entropy
1900  Qs = (smolar - SatL.smolar()) / (SatV.smolar() - SatL.smolar());
1901  // Quality from enthalpy
1902  CoolPropDbl Qh = (hmolar - SatL.hmolar()) / (SatV.hmolar() - SatL.hmolar());
1903  // Residual is the difference between the two
1904  return Qh - Qs;
1905  }
1906  } resid(HEOS, hmolar_spec, smolar_spec);
1907 
1908  // Critical point for pure fluids, slightly different for pseudo-pure, very different for mixtures
1909  CoolPropDbl Tmax_sat = HEOS.calc_Tmax_sat() - 1e-13;
1910 
1911  // Check what the minimum limits for the equation of state are
1912  CoolPropDbl Tmin_satL, Tmin_satV, Tmin_sat;
1913  HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
1914  Tmin_sat = std::max(Tmin_satL, Tmin_satV) - 1e-13;
1915 
1916  Brent(resid, Tmin_sat, Tmax_sat - 0.01, DBL_EPSILON, 1e-12, 20);
1917  // Run once more with the final vapor quality
1918  HEOS.update(QT_INPUTS, resid.Qs, HEOS.T());
1919 }
1921  HS_flash_singlephaseOptions& options) {
1922  int iter = 0;
1923  double resid = 9e30, resid_old = 9e30;
1924  CoolProp::SimpleState reducing = HEOS.get_state("reducing");
1925  do {
1926  // Independent variables are T0 and rhomolar0, residuals are matching h and s
1927  Eigen::Vector2d r;
1928  Eigen::Matrix2d J;
1929  r(0) = HEOS.hmolar() - hmolar_spec;
1930  r(1) = HEOS.smolar() - smolar_spec;
1931  J(0, 0) = HEOS.first_partial_deriv(iHmolar, iTau, iDelta);
1932  J(0, 1) = HEOS.first_partial_deriv(iHmolar, iDelta, iTau);
1933  J(1, 0) = HEOS.first_partial_deriv(iSmolar, iTau, iDelta);
1934  J(1, 1) = HEOS.first_partial_deriv(iSmolar, iDelta, iTau);
1935  // Step in v obtained from Jv = -r
1936  Eigen::Vector2d v = J.colPivHouseholderQr().solve(-r);
1937  bool good_solution = false;
1938  double tau0 = HEOS.tau(), delta0 = HEOS.delta();
1939  // Calculate the old residual after the last step
1940  resid_old = sqrt(POW2(HEOS.hmolar() - hmolar_spec) + POW2(HEOS.smolar() - smolar_spec));
1941  for (double frac = 1.0; frac > 0.001; frac /= 2) {
1942  try {
1943  // Calculate new values
1944  double tau_new = tau0 + options.omega * frac * v(0);
1945  double delta_new = delta0 + options.omega * frac * v(1);
1946  double T_new = reducing.T / tau_new;
1947  double rhomolar_new = delta_new * reducing.rhomolar;
1948  // Update state with step
1949  HEOS.update(DmolarT_INPUTS, rhomolar_new, T_new);
1950  resid = sqrt(POW2(HEOS.hmolar() - hmolar_spec) + POW2(HEOS.smolar() - smolar_spec));
1951  if (resid > resid_old) {
1952  throw ValueError(format("residual not decreasing; frac: %g, resid: %g, resid_old: %g", frac, resid, resid_old));
1953  }
1954  good_solution = true;
1955  break;
1956  } catch (...) {
1957  HEOS.clear();
1958  continue;
1959  }
1960  }
1961  if (!good_solution) {
1962  throw ValueError(format("Not able to get a solution"));
1963  }
1964  iter++;
1965  if (iter > 50) {
1966  throw ValueError(format("HS_flash_singlephase took too many iterations; residual is %g; prior was %g", resid, resid_old));
1967  }
1968  } while (std::abs(resid) > 1e-9);
1969 }
1971  // Randomly obtain a starting value that is single-phase
1972  double logp = ((double)rand() / (double)RAND_MAX) * (log(HEOS.pmax()) - log(HEOS.p_triple())) + log(HEOS.p_triple());
1973  T = ((double)rand() / (double)RAND_MAX) * (HEOS.Tmax() - HEOS.Ttriple()) + HEOS.Ttriple();
1974  p = exp(logp);
1975 }
1977  // Use TS flash and iterate on T (known to be between Tmin and Tmax)
1978  // in order to find H
1979  double hmolar = HEOS.hmolar(), smolar = HEOS.smolar();
1980  class Residual : public FuncWrapper1D
1981  {
1982  public:
1984  CoolPropDbl hmolar, smolar;
1985  Residual(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl hmolar_spec, CoolPropDbl smolar_spec)
1986  : HEOS(HEOS), hmolar(hmolar_spec), smolar(smolar_spec){};
1987  double call(double T) {
1988  HEOS.update(SmolarT_INPUTS, smolar, T);
1989  double r = HEOS.hmolar() - hmolar;
1990  return r;
1991  }
1992  } resid(HEOS, hmolar, smolar);
1993  std::string errstr;
1994  // Find minimum temperature
1995  bool good_Tmin = false;
1996  double Tmin = HEOS.Ttriple();
1997  double rmin;
1998  do {
1999  try {
2000  rmin = resid.call(Tmin);
2001  good_Tmin = true;
2002  } catch (...) {
2003  Tmin += 0.5;
2004  }
2005  if (Tmin > HEOS.Tmax()) {
2006  throw ValueError("Cannot find good Tmin");
2007  }
2008  } while (!good_Tmin);
2009 
2010  // Find maximum temperature
2011  bool good_Tmax = false;
2012  double Tmax = HEOS.Tmax() * 1.01; // Just a little above, so if we use Tmax as input, it should still work
2013  double rmax;
2014  do {
2015  try {
2016  rmax = resid.call(Tmax);
2017  good_Tmax = true;
2018  } catch (...) {
2019  Tmax -= 0.1;
2020  }
2021  if (Tmax < Tmin) {
2022  throw ValueError("Cannot find good Tmax");
2023  }
2024  } while (!good_Tmax);
2025  if (rmin * rmax > 0 && std::abs(rmax) < std::abs(rmin)) {
2026  throw CoolProp::ValueError(format("HS inputs correspond to temperature above maximum temperature of EOS [%g K]", HEOS.Tmax()));
2027  }
2028  Brent(resid, Tmin, Tmax, DBL_EPSILON, 1e-10, 100);
2029 }
2030 
2031 #if defined(ENABLE_CATCH)
2032 
2033 TEST_CASE("PD with T very large should yield error", "[PDflash]") {
2034  shared_ptr<HelmholtzEOSBackend> HEOS(new HelmholtzEOSBackend("R134a"));
2035  double Tc = HEOS->T_critical();
2036  HEOS->update(DmassT_INPUTS, 1.1, 1.5 * Tc);
2037  CHECK_THROWS(HEOS->update(DmassP_INPUTS, 2, 5 * HEOS->p()));
2038 }
2039 
2040 TEST_CASE("Stability testing", "[stability]") {
2041  shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend(strsplit("n-Propane&n-Butane&n-Pentane&n-Hexane", '&')));
2042  std::vector<double> z(4);
2043  z[0] = 0.1;
2044  z[1] = 0.2;
2045  z[2] = 0.3;
2046  z[3] = 0.4;
2047  HEOS->set_mole_fractions(z);
2048 
2049  HEOS->update(PQ_INPUTS, 101325, 0);
2050  double TL = HEOS->T(), rhoL = HEOS->rhomolar();
2051 
2052  HEOS->update(PQ_INPUTS, 101325, 1);
2053  double TV = HEOS->T(), rhoV = HEOS->rhomolar();
2054 
2055  SECTION("Liquid (feed is stable)") {
2056  StabilityRoutines::StabilityEvaluationClass stability_tester(*HEOS);
2057  for (double T = TL - 1; T >= 100; T -= 1) {
2058  stability_tester.set_TP(T, 101325);
2059  CAPTURE(T);
2060  CHECK_NOTHROW(stability_tester.is_stable());
2061  }
2062  }
2063  SECTION("Vapor (feed is stable)") {
2064  StabilityRoutines::StabilityEvaluationClass stability_tester(*HEOS);
2065  for (double T = TV + 1; T <= 500; T += 1) {
2066  stability_tester.set_TP(T, 101325);
2067  CAPTURE(T);
2068  CHECK_NOTHROW(stability_tester.is_stable());
2069  }
2070  }
2071  SECTION("Two-phase (feed is unstable)") {
2072  StabilityRoutines::StabilityEvaluationClass stability_tester(*HEOS);
2073  stability_tester.set_TP((TV + TL) / 2.0, 101325);
2074  CHECK(stability_tester.is_stable() == false);
2075  }
2076 }
2077 
2078 TEST_CASE("Test critical points for methane + H2S", "[critical_points]") {
2079  shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend(strsplit("Methane&H2S", '&')));
2080 
2081  double zz[] = {0.998, 0.97, 0.9475, 0.94, 0.93, 0.86, 0.85, 0.84, 0.75, 0.53, 0.52, 0.51, 0.49, 0.36, 0.24, 0.229, 0.09};
2082  int Npts[] = {2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 2, 2, 2, 1, 1, 1, 1}; // Number of critical points that should be found
2083  int imax = sizeof(zz) / sizeof(double);
2084 
2085  for (int i = 0; i < imax; ++i) {
2086  double z0 = zz[i];
2087  std::vector<double> z(2);
2088  z[0] = z0;
2089  z[1] = 1 - z0;
2090  HEOS->set_mole_fractions(z);
2091  CAPTURE(z0);
2092  std::vector<CriticalState> pts = HEOS->all_critical_points();
2093  CHECK(pts.size() == Npts[i]);
2094  }
2095 }
2096 
2097 TEST_CASE("Test critical points for nitrogen + ethane with HEOS", "[critical_points]") {
2098  shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend(strsplit("Nitrogen&Ethane", '&')));
2099  std::vector<double> zz = linspace(0.001, 0.999, 21);
2100  for (int i = 0; i < static_cast<std::size_t>(zz.size()); ++i) {
2101  double z0 = zz[i];
2102  std::vector<double> z(2);
2103  z[0] = z0;
2104  z[1] = 1 - z0;
2105  HEOS->set_mole_fractions(z);
2106  CAPTURE(z0);
2107  std::vector<CriticalState> pts;
2108  CHECK_NOTHROW(pts = HEOS->all_critical_points());
2109  }
2110 }
2111 
2112 TEST_CASE("Test critical points for nitrogen + ethane with PR", "[critical_points]") {
2113  shared_ptr<PengRobinsonBackend> HEOS(new PengRobinsonBackend(strsplit("Nitrogen&Ethane", '&')));
2114  HEOS->set_binary_interaction_double(0, 1, "kij", 0.0407); // Ramırez-Jimenez et al.
2115  std::vector<double> zz = linspace(0.001, 0.999, 21);
2116  for (int i = 0; i < static_cast<std::size_t>(zz.size()); ++i) {
2117  double z0 = zz[i];
2118  std::vector<double> z(2);
2119  z[0] = z0;
2120  z[1] = 1 - z0;
2121  HEOS->set_mole_fractions(z);
2122  CAPTURE(z0);
2123  std::vector<CriticalState> pts;
2124  CHECK_NOTHROW(pts = HEOS->all_critical_points());
2125  }
2126 }
2127 
2128 #endif
2129 
2130 } /* namespace CoolProp */