CoolProp  6.6.1dev
An open-source fluid property and humid air property database
VLERoutines.h
Go to the documentation of this file.
1 
2 #ifndef VLEROUTINES_H
3 #define VLEROUTINES_H
4 
6 
7 namespace CoolProp {
8 
9 namespace SaturationSolvers {
11 {
12  bool use_guesses;
15  : use_guesses(use_guesses), omega(_HUGE), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE) {}
16 };
18 {
20  saturation_T_pure_options() : omega(_HUGE), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE), p(_HUGE), T(_HUGE) {}
21 };
22 
24 {
26  {
29  };
30  bool use_guesses,
35  saturation_D_pure_options() : use_guesses(false), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE), imposed_rho(0), max_iterations(200) {
36  use_logdelta = true;
37  omega = 1.0;
38  } // Defaults
39 };
40 
42 {
44  imposed_p
45 };
47 {
49  int Nstep_max;
51  std::vector<CoolPropDbl> x, y, K;
52 };
53 
63 static CoolPropDbl Wilson_lnK_factor(const HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, CoolPropDbl p, std::size_t i) {
64  const double pci = HEOS.get_fluid_constant(i, iP_critical);
65  const double Tci = HEOS.get_fluid_constant(i, iT_critical);
66  const double omegai = HEOS.get_fluid_constant(i, iacentric_factor);
67  return log(pci / p) + 5.373 * (1 + omegai) * (1 - Tci / T);
68 };
69 
70 void saturation_D_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl rhomolar, saturation_D_pure_options& options);
71 void saturation_T_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options);
72 void saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options);
73 void saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options);
74 
76 {
78  {
88  };
89  bool use_guesses,
93  saturation_PHSU_pure_options() : use_logdelta(true), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE), T(_HUGE), p(_HUGE) {
95  use_guesses = true;
96  omega = 1.0;
97  }
98 };
99 
100 #if !defined(NO_FMTLIB) && FMT_VERSION >= 90000
102  return fmt::underlying(option);
103 }
104 #endif
105 
106 void saturation_PHSU_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl specified_value, saturation_PHSU_pure_options& options);
107 
108 /* \brief This is a backup saturation_p solver for the case where the Newton solver cannot approach closely enough the solution
109  *
110  * This is especially a problem at low pressures where catastrophic truncation error occurs, especially in the saturated vapor side
111  *
112  * @param HEOS The Helmholtz EOS backend instance to be used
113  * @param p Imposed pressure in kPa
114  * @param options Options to be passed to the function (at least T, rhoL and rhoV must be provided)
115  */
116 void saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl p, saturation_PHSU_pure_options& options);
117 
118 /* \brief This is a backup saturation_T solver for the case where the Newton solver cannot approach closely enough the solution
119  *
120  * This is especially a problem at low pressures where catastrophic truncation error occurs, especially in the saturated vapor side
121  *
122  * @param HEOS The Helmholtz EOS backend instance to be used
123  * @param T Imposed temperature in K
124  * @param options Options to be passed to the function (at least p, rhoL and rhoV must be provided)
125  */
126 void saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options);
127 
128 /* \brief A robust but slow solver in the very-near-critical region
129  *
130  * This solver operates in the following fashion:
131  * 1. Using a bounded interval for rho'':[rhoc, rhoc-??], guess a value for rho''
132  * 2. For guessed value of rho'' and given value of T, calculate p
133  * 3. Using a Brent solver on the other co-existing phase (rho'), calculate the (bounded) value of rho' that yields the same pressure
134  * 4. Use another outer Brent solver on rho'' to enforce the same Gibbs function between liquid and vapor
135  * 5. Fin.
136  *
137  * @param HEOS The Helmholtz EOS backend instance to be used
138  * @param ykey The CoolProp::parameters key to be imposed - one of iT or iP
139  * @param y The value for the imposed variable
140  */
141 void saturation_critical(HelmholtzEOSMixtureBackend& HEOS, CoolProp::parameters ykey, CoolPropDbl y);
142 
143 void successive_substitution(HelmholtzEOSMixtureBackend& HEOS, const CoolPropDbl beta, CoolPropDbl T, CoolPropDbl p,
144  const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& K, mixture_VLE_IO& options);
152 void x_and_y_from_K(CoolPropDbl beta, const std::vector<CoolPropDbl>& K, const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& x,
153  std::vector<CoolPropDbl>& y);
154 
161 {
162  public:
164  double T, p, beta;
165  const std::vector<CoolPropDbl>& z;
166  std::vector<CoolPropDbl>& K;
168 
170  const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& K)
172  T(imposed_value),
173  p(imposed_value),
174  beta(beta),
175  z(z),
176  K(K),
177  HEOS(HEOS) {} // if input_type == imposed_T -> use T, else use p; init both
178  double call(double input_value) {
179  double summer = 0;
180  if (input_type == imposed_T) {
181  p = input_value; // Iterate on pressure
182  } else {
183  T = input_value; // Iterate on temperature, pressure imposed
184  }
185  for (unsigned int i = 0; i < z.size(); i++) {
186  K[i] = exp(Wilson_lnK_factor(HEOS, T, p, i));
187  summer += z[i] * (K[i] - 1) / (1 - beta + beta * K[i]);
188  }
189  return summer;
190  }
191 };
192 inline double saturation_preconditioner(HelmholtzEOSMixtureBackend& HEOS, double input_value, sstype_enum input_type,
193  const std::vector<CoolPropDbl>& z) {
194  double ptriple = 0, pcrit = 0, Ttriple = 0, Tcrit = 0;
195 
196  if (HEOS.get_components().empty()) {
197  return -1;
198  }
199 
200  for (unsigned int i = 0; i < z.size(); i++) {
201  Tcrit += HEOS.get_fluid_constant(i, iT_critical) * z[i];
202  pcrit += HEOS.get_fluid_constant(i, iP_critical) * z[i];
203  Ttriple += HEOS.get_fluid_constant(i, iT_triple) * z[i];
204  ptriple += HEOS.get_fluid_constant(i, iP_triple) * z[i];
205  }
206  // Return an invalid number if either triple point temperature or pressure are not available
207  if (!ValidNumber(Ttriple) || !ValidNumber(ptriple)) {
208  return _HUGE;
209  }
210 
211  if (input_type == imposed_T) {
212  return exp(log(pcrit / ptriple) / (Tcrit - Ttriple) * (input_value - Ttriple) + log(ptriple));
213  } else if (input_type == imposed_p) {
214  return 1 / (1 / Tcrit - (1 / Ttriple - 1 / Tcrit) / log(pcrit / ptriple) * log(input_value / pcrit));
215  } else {
216  throw ValueError();
217  }
218 }
246 inline double saturation_Wilson(HelmholtzEOSMixtureBackend& HEOS, double beta, double input_value, sstype_enum input_type,
247  const std::vector<CoolPropDbl>& z, double guess) {
248  double out = 0;
249  std::string errstr;
250 
251  // If T is input and beta = 0 or beta = 1, explicit solution for p is possible
252  if (input_type == imposed_T && (std::abs(beta) < 1e-12 || std::abs(beta - 1) < 1e-12)) {
253  const std::vector<double> z = HEOS.get_mole_fractions_ref();
254  bool beta0 = std::abs(beta) < 1e-12; // True is beta is approx. zero
255  for (int i = 0; i < static_cast<int>(z.size()); ++i) {
256  double pci = HEOS.get_fluid_constant(i, iP_critical);
257  double Tci = HEOS.get_fluid_constant(i, iT_critical);
258  double omegai = HEOS.get_fluid_constant(i, iacentric_factor);
259  if (beta0) {
260  out += z[i] * pci * exp(5.373 * (1 + omegai) * (1 - Tci / input_value));
261  } else {
262  out += z[i] / (pci * exp(5.373 * (1 + omegai) * (1 - Tci / input_value)));
263  }
264  }
265  if (!beta0) { // beta = 1
266  out = 1 / out; // summation is for 1/p, take reciprocal to get p
267  }
268  std::vector<CoolPropDbl>& K = HEOS.get_K();
269  for (int i = 0; i < static_cast<int>(z.size()); ++i) {
270  double pci = HEOS.get_fluid_constant(i, iP_critical);
271  double Tci = HEOS.get_fluid_constant(i, iT_critical);
272  double omegai = HEOS.get_fluid_constant(i, iacentric_factor);
273  K[i] = pci / out * exp(5.373 * (1 + omegai) * (1 - Tci / input_value));
274  }
275  } else {
276  // Find first guess for output variable using Wilson K-factors
277  WilsonK_resid Resid(HEOS, beta, input_value, input_type, z, HEOS.get_K());
278  if (guess < 0 || !ValidNumber(guess))
279  out = Brent(Resid, 50, 10000, 1e-10, 1e-10, 100);
280  else
281  out = Secant(Resid, guess, 0.001, 1e-10, 100);
282  if (!ValidNumber(out)) {
283  throw ValueError("saturation_p_Wilson failed to get good output value");
284  }
285  }
286  return out;
287 }
289 {
291 };
292 
294 {
296  {
299  T_IMPOSED
300  };
302  std::size_t Nsteps;
305  std::vector<CoolPropDbl> x, y, z;
307  : Nstep_max(30),
308  Nsteps(0),
309  beta(-1),
310  omega(1),
311  rhomolar_liq(_HUGE),
312  rhomolar_vap(_HUGE),
313  pL(_HUGE),
314  pV(_HUGE),
315  p(_HUGE),
316  T(_HUGE),
317  hmolar_liq(_HUGE),
318  hmolar_vap(_HUGE),
319  smolar_liq(_HUGE),
320  smolar_vap(_HUGE),
321  imposed_variable(NO_VARIABLE_IMPOSED) {} // Defaults
322 };
323 
355 {
356  public:
360  std::size_t N;
361  bool logging;
362  int Nsteps;
363  Eigen::MatrixXd J;
364  Eigen::Vector2d r, err_rel;
365  std::vector<CoolPropDbl> K, x, y, z;
366  std::vector<SuccessiveSubstitutionStep> step_logger;
367 
369  : HEOS(NULL),
371  error_rms(_HUGE),
372  rhomolar_liq(_HUGE),
373  rhomolar_vap(_HUGE),
374  T(_HUGE),
375  p(_HUGE),
376  min_rel_change(_HUGE),
377  beta(_HUGE),
378  N(0),
379  logging(false),
380  Nsteps(0){};
381 
382  void resize(unsigned int N);
383 
384  // Reset the state of all the internal variables
385  void pre_call() {
386  K.clear();
387  x.clear();
388  y.clear();
389  step_logger.clear();
390  error_rms = 1e99;
391  Nsteps = 0;
392  rhomolar_liq = _HUGE;
393  rhomolar_vap = _HUGE;
394  T = _HUGE;
395  p = _HUGE;
396  };
397 
408 
409  /* \brief Build the arrays for the Newton-Raphson solve
410  *
411  */
412  void build_arrays();
413 };
414 
416 {
418  {
422  T_IMPOSED
423  };
426  std::size_t Nsteps;
429  std::vector<CoolPropDbl> x, y;
431  : bubble_point(false),
432  omega(_HUGE),
433  rhomolar_liq(_HUGE),
434  rhomolar_vap(_HUGE),
435  pL(_HUGE),
436  pV(_HUGE),
437  p(_HUGE),
438  T(_HUGE),
439  hmolar_liq(_HUGE),
440  hmolar_vap(_HUGE),
441  smolar_liq(_HUGE),
442  smolar_vap(_HUGE),
444  Nstep_max = 30;
445  Nsteps = 0;
446  } // Defaults
447 };
448 
477 {
478  public:
481  std::size_t N;
482  bool logging;
484  int Nsteps;
485  Eigen::MatrixXd J;
488  std::vector<CoolPropDbl> K, x, y;
489  Eigen::VectorXd r, err_rel;
490  std::vector<SuccessiveSubstitutionStep> step_logger;
491 
493 
494  void resize(std::size_t N);
495 
496  // Reset the state of all the internal variables
497  void pre_call() {
498  step_logger.clear();
499  error_rms = 1e99;
500  Nsteps = 0;
501  rhomolar_liq = _HUGE;
502  rhomolar_vap = _HUGE;
503  T = _HUGE;
504  p = _HUGE;
505  };
506 
518  void call(HelmholtzEOSMixtureBackend& HEOS, const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& z_incipient,
520 
526  void build_arrays();
527 
530  void check_Jacobian();
531 };
532 
534 {
536  std::size_t Nsteps;
538  std::vector<CoolPropDbl> x,
539  y,
540  z;
541  PTflash_twophase_options() : omega(_HUGE), rhomolar_liq(_HUGE), rhomolar_vap(_HUGE), pL(_HUGE), pV(_HUGE), p(_HUGE), T(_HUGE) {
542  Nstep_max = 30;
543  Nsteps = 0; // Defaults
544  }
545 };
546 
548 {
549  public:
550  double error_rms;
551  bool logging;
552  int Nsteps;
553  Eigen::MatrixXd J;
554  Eigen::VectorXd r, err_rel;
557  std::vector<SuccessiveSubstitutionStep> step_logger;
558 
560 
568  void solve();
569 
575  void build_arrays();
576 };
577 }; // namespace SaturationSolvers
578 
579 namespace StabilityRoutines {
580 
585 {
586  protected:
588  std::vector<double> lnK, K, K0, x, y, xL, xH;
589  const std::vector<double>& z;
591  double m_T,
592  m_p;
593  private:
594  bool _stable;
595  bool debug;
596 
597  public:
599  : HEOS(HEOS),
600  z(HEOS.get_mole_fractions_doubleref()),
601  rhomolar_liq(-1),
602  rhomolar_vap(-1),
603  beta(-1),
604  tpd_liq(10000),
605  tpd_vap(100000),
606  DELTAG_nRT(10000),
607  m_T(-1),
608  m_p(-1),
609  _stable(false),
610  debug(false){};
613  void set_TP(double T, double p) {
614  m_T = T;
615  m_p = p;
616  };
619  void rho_TP_w_guesses();
622  void rho_TP_global();
625  void rho_TP_SRK_translated();
626 
629  void trial_compositions();
632  void successive_substitution(int num_steps);
637  void check_stability();
640  bool is_stable() {
643  check_stability();
644  return _stable;
645  }
647  void get_liq(std::vector<double>& x, double& rhomolar) {
648  x = this->x;
649  rhomolar = rhomolar_liq;
650  }
652  void get_vap(std::vector<double>& y, double& rhomolar) {
653  y = this->y;
654  rhomolar = rhomolar_vap;
655  }
656 };
657 
658 }; /* namespace StabilityRoutines*/
659 
660 }; /* namespace CoolProp*/
661 
662 #endif