CoolProp  6.6.1dev
An open-source fluid property and humid air property database
VLERoutines.cpp
Go to the documentation of this file.
1 
2 #include "HelmholtzEOSMixtureBackend.h"
3 #include "VLERoutines.h"
4 #include "MatrixMath.h"
5 #include "MixtureDerivatives.h"
6 #include "Configuration.h"
7 #include "FlashRoutines.h"
8 
9 namespace CoolProp {
10 
11 void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend& HEOS, parameters ykey, CoolPropDbl y) {
12 
13  class inner_resid : public FuncWrapper1D
14  {
15  public:
16  HelmholtzEOSMixtureBackend* HEOS;
17  CoolPropDbl T, desired_p;
18 
19  inner_resid(HelmholtzEOSMixtureBackend* HEOS, CoolPropDbl T, CoolPropDbl desired_p) : HEOS(HEOS), T(T), desired_p(desired_p){};
20  double call(double rhomolar_liq) {
21  HEOS->SatL->update(DmolarT_INPUTS, rhomolar_liq, T);
22  CoolPropDbl calc_p = HEOS->SatL->p();
23  std::cout << format("inner p: %0.16Lg; res: %0.16Lg", calc_p, calc_p - desired_p) << std::endl;
24  return calc_p - desired_p;
25  }
26  };
27 
28  // Define the outer residual to be driven to zero - this is the equality of
29  // Gibbs function for both co-existing phases
30  class outer_resid : public FuncWrapper1D
31  {
32  public:
33  HelmholtzEOSMixtureBackend* HEOS;
34  parameters ykey;
35  CoolPropDbl y;
36  CoolPropDbl rhomolar_crit;
37 
38  outer_resid(HelmholtzEOSMixtureBackend& HEOS, CoolProp::parameters ykey, CoolPropDbl y)
39  : HEOS(&HEOS), ykey(ykey), y(y), rhomolar_crit(HEOS.rhomolar_critical()){};
40  double call(double rhomolar_vap) {
41  // Calculate the other variable (T->p or p->T) for given vapor density
42  CoolPropDbl T, p, rhomolar_liq;
43  switch (ykey) {
44  case iT: {
45  T = y;
46  HEOS->SatV->update(DmolarT_INPUTS, rhomolar_vap, y);
47  p = HEOS->SatV->p();
48  std::cout << format("outer p: %0.16Lg", p) << std::endl;
49  inner_resid inner(HEOS, T, p);
50  rhomolar_liq = Brent(inner, rhomolar_crit * 1.5, rhomolar_crit * (1 + 1e-8), LDBL_EPSILON, 1e-10, 100);
51  break;
52  }
53  default:
54  throw ValueError("Wrong input for outer_resid");
55  }
56  HEOS->SatL->update(DmolarT_INPUTS, rhomolar_liq, T);
57  HEOS->SatV->update(DmolarT_INPUTS, rhomolar_vap, T);
58 
59  // Calculate the Gibbs functions for liquid and vapor
60  //CoolPropDbl gL = HEOS->SatL->gibbsmolar();
61  //CoolPropDbl gV = HEOS->SatV->gibbsmolar();
62 
63  // Residual is difference in Gibbs function
64  // r = gL - gV;
65 
66  return p;
67  };
68  };
69  outer_resid resid(HEOS, iT, y);
70 
71  double rhomolar_crit = HEOS.rhomolar_critical();
72 
73  Brent(&resid, rhomolar_crit * (1 - 1e-8), rhomolar_crit * 0.5, DBL_EPSILON, 1e-9, 20);
74 }
75 
76 void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options) {
77 
78  // Define the residual to be driven to zero
79  class solver_resid : public FuncWrapper1D
80  {
81  public:
82  HelmholtzEOSMixtureBackend* HEOS;
83  CoolPropDbl T, rhomolar_liq, rhomolar_vap;
84 
85  solver_resid(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, CoolPropDbl rhomolar_liq_guess, CoolPropDbl rhomolar_vap_guess)
86  : HEOS(&HEOS), T(T), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){};
87  double call(double p) {
88  // Recalculate the densities using the current guess values
89  HEOS->SatL->update_TP_guessrho(T, p, rhomolar_liq);
90  HEOS->SatV->update_TP_guessrho(T, p, rhomolar_vap);
91 
92  // Calculate the Gibbs functions for liquid and vapor
93  CoolPropDbl gL = HEOS->SatL->gibbsmolar();
94  CoolPropDbl gV = HEOS->SatV->gibbsmolar();
95 
96  // Residual is difference in Gibbs function
97  return gL - gV;
98  };
99  };
100  solver_resid resid(HEOS, T, options.rhoL, options.rhoV);
101 
102  if (!ValidNumber(options.p)) {
103  throw ValueError(format("options.p is not valid in saturation_T_pure_1D_P for T = %Lg", T));
104  };
105  if (!ValidNumber(options.rhoL)) {
106  throw ValueError(format("options.rhoL is not valid in saturation_T_pure_1D_P for T = %Lg", T));
107  };
108  if (!ValidNumber(options.rhoV)) {
109  throw ValueError(format("options.rhoV is not valid in saturation_T_pure_1D_P for T = %Lg", T));
110  };
111 
112  try {
113  Secant(resid, options.p, options.p * 1.1, 1e-10, 100);
114  } catch (...) {
115  CoolPropDbl pmax = std::min(options.p * 1.03, static_cast<CoolPropDbl>(HEOS.p_critical() + 1e-6));
116  CoolPropDbl pmin = std::max(options.p * 0.97, static_cast<CoolPropDbl>(HEOS.p_triple() - 1e-6));
117  Brent(resid, pmin, pmax, LDBL_EPSILON, 1e-8, 100);
118  }
119 }
120 
121 void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl p, saturation_PHSU_pure_options& options) {
122 
123  // Define the residual to be driven to zero
124  class solver_resid : public FuncWrapper1D
125  {
126  public:
127  HelmholtzEOSMixtureBackend* HEOS;
128  CoolPropDbl p, rhomolar_liq, rhomolar_vap;
129 
130  solver_resid(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl p, CoolPropDbl rhomolar_liq_guess, CoolPropDbl rhomolar_vap_guess)
131  : HEOS(&HEOS), p(p), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){};
132  double call(double T) {
133  // Recalculate the densities using the current guess values
134  HEOS->SatL->update_TP_guessrho(T, p, rhomolar_liq);
135  HEOS->SatV->update_TP_guessrho(T, p, rhomolar_vap);
136 
137  // Calculate the Gibbs functions for liquid and vapor
138  CoolPropDbl gL = HEOS->SatL->gibbsmolar();
139  CoolPropDbl gV = HEOS->SatV->gibbsmolar();
140 
141  // Residual is difference in Gibbs function
142  return gL - gV;
143  };
144  };
145  solver_resid resid(HEOS, p, options.rhoL, options.rhoV);
146 
147  if (!ValidNumber(options.T)) {
148  throw ValueError("options.T is not valid in saturation_P_pure_1D_T");
149  };
150  if (!ValidNumber(options.rhoL)) {
151  throw ValueError("options.rhoL is not valid in saturation_P_pure_1D_T");
152  };
153  if (!ValidNumber(options.rhoV)) {
154  throw ValueError("options.rhoV is not valid in saturation_P_pure_1D_T");
155  };
156 
157  CoolPropDbl Tmax = std::min(options.T + 2, static_cast<CoolPropDbl>(HEOS.T_critical() - 1e-6));
158  CoolPropDbl Tmin = std::max(options.T - 2, static_cast<CoolPropDbl>(HEOS.Ttriple() + 1e-6));
159  Brent(resid, Tmin, Tmax, LDBL_EPSILON, 1e-11, 100);
160 }
161 
162 void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl specified_value, saturation_PHSU_pure_options& options) {
163  /*
164  This function is inspired by the method of Akasaka:
165 
166  R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
167  Helmholtz Energy Equations of State",
168  Journal of Thermal Science and Technology v3 n3,2008
169 
170  Ancillary equations are used to get a sensible starting point
171  */
172  std::vector<CoolPropDbl> negativer(3, _HUGE), v;
173  std::vector<std::vector<CoolPropDbl>> J(3, std::vector<CoolPropDbl>(3, _HUGE));
174 
175  HEOS.calc_reducing_state();
176  const SimpleState& reduce = HEOS.get_reducing_state();
177  CoolProp::SimpleState crit = HEOS.get_state("reducing");
178  shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
179 
180  CoolPropDbl T, rhoL, rhoV, pL, pV, hL, sL, hV, sV;
181  CoolPropDbl deltaL = 0, deltaV = 0, tau = 0, error;
182  int iter = 0, specified_parameter;
183 
184  // Use the density ancillary function as the starting point for the solver
185  try {
186  if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_PL
187  || options.specified_variable == saturation_PHSU_pure_options::IMPOSED_PV) {
188  // Invert liquid density ancillary to get temperature
189  // TODO: fit inverse ancillaries too
190  try {
191  T = HEOS.get_components()[0].ancillaries.pL.invert(specified_value);
192  } catch (...) {
193  throw ValueError("Unable to invert ancillary equation");
194  }
195  } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_HL) {
196  CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
197  // Ancillary is deltah = h - hs_anchor.h
198  try {
199  T = HEOS.get_components()[0].ancillaries.hL.invert(specified_value - hs_anchor.hmolar);
200  } catch (...) {
201  throw ValueError("Unable to invert ancillary equation for hL");
202  }
203  } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_HV) {
204  class Residual : public FuncWrapper1D
205  {
206  public:
207  CoolPropFluid* component;
208  double h;
209  Residual(CoolPropFluid& component, double h) {
210  this->component = &component;
211  this->h = h;
212  }
213  double call(double T) {
214  CoolPropDbl h_liq = component->ancillaries.hL.evaluate(T) + component->EOS().hs_anchor.hmolar;
215  return h_liq + component->ancillaries.hLV.evaluate(T) - h;
216  };
217  };
218  Residual resid(HEOS.get_components()[0], HEOS.hmolar());
219 
220  // Ancillary is deltah = h - hs_anchor.h
221  CoolPropDbl Tmin_satL, Tmin_satV;
222  HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
223  double Tmin = Tmin_satL;
224  double Tmax = HEOS.calc_Tmax_sat();
225  try {
226  T = Brent(resid, Tmin - 3, Tmax + 1, DBL_EPSILON, 1e-10, 50);
227  } catch (...) {
228  shared_ptr<HelmholtzEOSMixtureBackend> HEOS_copy(new HelmholtzEOSMixtureBackend(HEOS.get_components()));
229  HEOS_copy->update(QT_INPUTS, 1, Tmin);
230  double hTmin = HEOS_copy->hmolar();
231  HEOS_copy->update(QT_INPUTS, 1, Tmax);
232  double hTmax = HEOS_copy->hmolar();
233  T = (Tmax - Tmin) / (hTmax - hTmin) * (HEOS.hmolar() - hTmin) + Tmin;
234  }
235  } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SL) {
236  CoolPropFluid& component = HEOS.get_components()[0];
237  CoolProp::SaturationAncillaryFunction& anc = component.ancillaries.sL;
238  CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
239  // If near the critical point, use a near critical guess value for T
240  if (std::abs(HEOS.smolar() - crit.smolar) < std::abs(component.ancillaries.sL.get_max_abs_error())) {
241  T = std::max(0.99 * crit.T, crit.T - 0.1);
242  } else {
243  CoolPropDbl Tmin, Tmax, Tmin_satV;
244  HEOS.calc_Tmin_sat(Tmin, Tmin_satV);
245  Tmax = HEOS.calc_Tmax_sat();
246  // Ancillary is deltas = s - hs_anchor.s
247  // First try a conventional call
248  try {
249  T = anc.invert(specified_value - hs_anchor.smolar, Tmin, Tmax);
250  } catch (...) {
251  try {
252  T = anc.invert(specified_value - hs_anchor.smolar, Tmin - 3, Tmax + 3);
253  } catch (...) {
254  double vmin = anc.evaluate(Tmin);
255  double vmax = anc.evaluate(Tmax);
256  if (std::abs(specified_value - hs_anchor.smolar) < std::abs(vmax)) {
257  T = Tmax - 0.1;
258  } else {
259  throw ValueError(format("Unable to invert ancillary equation for sL for value %Lg with Tminval %g and Tmaxval %g ",
260  specified_value - hs_anchor.smolar, vmin, vmax));
261  }
262  }
263  }
264  }
265  } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SV) {
266  CoolPropFluid& component = HEOS.get_components()[0];
267  CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
268  class Residual : public FuncWrapper1D
269  {
270  public:
271  CoolPropFluid* component;
272  double s;
273  Residual(CoolPropFluid& component, double s) {
274  this->component = &component;
275  this->s = s;
276  }
277  double call(double T) {
278  CoolPropDbl s_liq = component->ancillaries.sL.evaluate(T) + component->EOS().hs_anchor.smolar;
279  CoolPropDbl resid = s_liq + component->ancillaries.sLV.evaluate(T) - s;
280 
281  return resid;
282  };
283  };
284  Residual resid(component, HEOS.smolar());
285 
286  // Ancillary is deltas = s - hs_anchor.s
287  CoolPropDbl Tmin_satL, Tmin_satV;
288  HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
289  double Tmin = Tmin_satL;
290  double Tmax = HEOS.calc_Tmax_sat();
291  try {
292  T = Brent(resid, Tmin - 3, Tmax, DBL_EPSILON, 1e-10, 50);
293  } catch (...) {
294  CoolPropDbl vmax = resid.call(Tmax);
295  // If near the critical point, use a near critical guess value for T
296  if (std::abs(specified_value - hs_anchor.smolar) < std::abs(vmax)) {
297  T = std::max(0.99 * crit.T, crit.T - 0.1);
298  } else {
299  shared_ptr<HelmholtzEOSMixtureBackend> HEOS_copy(new HelmholtzEOSMixtureBackend(HEOS.get_components()));
300  HEOS_copy->update(QT_INPUTS, 1, Tmin);
301  double sTmin = HEOS_copy->smolar();
302  HEOS_copy->update(QT_INPUTS, 1, Tmax);
303  double sTmax = HEOS_copy->smolar();
304  T = (Tmax - Tmin) / (sTmax - sTmin) * (HEOS.smolar() - sTmin) + Tmin;
305  }
306  }
307  } else {
308  throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid", options.specified_variable));
309  }
310  // If T from the ancillaries is above the critical temp, this will cause failure
311  // in ancillaries for rhoV and rhoL, decrease if needed
312  T = std::min(T, static_cast<CoolPropDbl>(HEOS.T_critical() - 0.1));
313 
314  // Evaluate densities from the ancillary equations
315  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
316  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
317 
318  // Apply a single step of Newton's method to improve guess value for liquid
319  // based on the error between the gas pressure (which is usually very close already)
320  // and the liquid pressure, which can sometimes (especially at low pressure),
321  // be way off, and often times negative
322  SatL->update(DmolarT_INPUTS, rhoL, T);
323  SatV->update(DmolarT_INPUTS, rhoV, T);
324  double rhoL_updated = rhoL - (SatL->p() - SatV->p()) / SatL->first_partial_deriv(iP, iDmolar, iT);
325 
326  // Accept the update if the liquid density is greater than the vapor density
327  if (rhoL_updated > rhoV) {
328  rhoL = rhoL_updated;
329  }
330 
331  // Update the state again with the better guess for the liquid density
332  SatL->update(DmolarT_INPUTS, rhoL, T);
333  SatV->update(DmolarT_INPUTS, rhoV, T);
334 
335  deltaL = rhoL / reduce.rhomolar;
336  deltaV = rhoV / reduce.rhomolar;
337  tau = reduce.T / T;
338  } catch (NotImplementedError&) {
339  throw; // ??? What is this try...catch for?
340  }
341 
342  do {
343  /*if (get_debug_level()>8){
344  std::cout << format("%s:%d: right before the derivs with deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
345  }*/
346 
347  pL = SatL->p();
348  hL = SatL->hmolar();
349  sL = SatL->smolar();
350  pV = SatV->p();
351  hV = SatV->hmolar();
352  sV = SatV->smolar();
353 
354  // These derivatives are needed for both cases
355  CoolPropDbl alpharL = SatL->alphar();
356  CoolPropDbl alpharV = SatV->alphar();
357  CoolPropDbl dalphar_dtauL = SatL->dalphar_dTau();
358  CoolPropDbl dalphar_dtauV = SatV->dalphar_dTau();
359  CoolPropDbl d2alphar_ddelta_dtauL = SatL->d2alphar_dDelta_dTau();
360  CoolPropDbl d2alphar_ddelta_dtauV = SatV->d2alphar_dDelta_dTau();
361  CoolPropDbl dalphar_ddeltaL = SatL->dalphar_dDelta();
362  CoolPropDbl dalphar_ddeltaV = SatV->dalphar_dDelta();
363  CoolPropDbl d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
364  CoolPropDbl d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
365 
366  // -r_1 (equate the pressures)
367  negativer[0] = -(deltaV * (1 + deltaV * dalphar_ddeltaV) - deltaL * (1 + deltaL * dalphar_ddeltaL));
368  // -r_2 (equate the gibbs energy)
369  negativer[1] = -(deltaV * dalphar_ddeltaV + alpharV + log(deltaV) - deltaL * dalphar_ddeltaL - alpharL - log(deltaL));
370  switch (options.specified_variable) {
371  case saturation_PHSU_pure_options::IMPOSED_PL:
372  // -r_3 (equate calculated pressure and specified liquid pressure)
373  negativer[2] = -(pL / specified_value - 1);
374  break;
375  case saturation_PHSU_pure_options::IMPOSED_PV:
376  // -r_3 (equate calculated pressure and specified vapor pressure)
377  negativer[2] = -(pV / specified_value - 1);
378  break;
379  case saturation_PHSU_pure_options::IMPOSED_HL:
380  // -r_3 (equate calculated liquid enthalpy and specified liquid enthalpy)
381  negativer[2] = -(hL - specified_value);
382  break;
383  case saturation_PHSU_pure_options::IMPOSED_HV:
384  // -r_3 (equate calculated vapor enthalpy and specified vapor enthalpy)
385  negativer[2] = -(hV - specified_value);
386  break;
387  case saturation_PHSU_pure_options::IMPOSED_SL:
388  // -r_3 (equate calculated liquid entropy and specified liquid entropy)
389  negativer[2] = -(sL - specified_value);
390  break;
391  case saturation_PHSU_pure_options::IMPOSED_SV:
392  // -r_3 (equate calculated vapor entropy and specified vapor entropy)
393  negativer[2] = -(sV - specified_value);
394  break;
395  default:
396  throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid", options.specified_variable));
397  }
398 
399  // dr1_dtau
400  J[0][0] = pow(deltaV, 2) * d2alphar_ddelta_dtauV - pow(deltaL, 2) * d2alphar_ddelta_dtauL;
401  // dr2_dtau
402  J[1][0] = deltaV * d2alphar_ddelta_dtauV + dalphar_dtauV - deltaL * d2alphar_ddelta_dtauL - dalphar_dtauL;
403 
404  if (options.use_logdelta) {
405  // dr_1/d_log(delta'')
406  J[0][1] = -deltaL - 2 * pow(deltaL, 2) * dalphar_ddeltaL - pow(deltaL, 3) * d2alphar_ddelta2L;
407  // dr_2/d_log(delta'')
408  J[1][1] = -pow(deltaL, 2) * d2alphar_ddelta2L - 2 * deltaL * dalphar_ddeltaL - 1;
409  } else {
410  // dr_1/ddelta''
411  J[0][1] = -1 - 2 * deltaL * dalphar_ddeltaL - pow(deltaL, 2) * d2alphar_ddelta2L;
412  // dr_2/ddelta''
413  J[1][1] = -1 / deltaL - 2 * dalphar_ddeltaL - deltaL * d2alphar_ddelta2L;
414  }
415 
416  if (options.use_logdelta) {
417  // dr_1/d_log(delta'')
418  J[0][2] = deltaV + 2 * pow(deltaV, 2) * dalphar_ddeltaV + pow(deltaV, 3) * d2alphar_ddelta2V;
419  // dr_2/d_log(delta'')
420  J[1][2] = 1 + 2 * deltaV * dalphar_ddeltaV + 1 + pow(deltaV, 2) * d2alphar_ddelta2V;
421  } else {
422  // dr_1/ddelta''
423  J[0][2] = 1 + 2 * deltaV * dalphar_ddeltaV + pow(deltaV, 2) * d2alphar_ddelta2V;
424  // dr_2/ddelta''
425  J[1][2] = deltaV * d2alphar_ddelta2V + 2 * dalphar_ddeltaV + 1 / deltaV;
426  }
427 
428  // Derivatives of the specification equation
429  switch (options.specified_variable) {
430  case saturation_PHSU_pure_options::IMPOSED_PL:
431  // dr_3/dtau
432  J[2][0] = SatL->first_partial_deriv(iP, iTau, iDelta) / specified_value;
433  if (options.use_logdelta) {
434  // dr_3/d(log(delta'))
435  J[2][1] = deltaL * SatL->first_partial_deriv(iP, iDelta, iTau) / specified_value;
436  } else {
437  // dr_3/ddelta'
438  J[2][1] = SatL->first_partial_deriv(iP, iDelta, iTau) / specified_value;
439  }
440  // dr_3/ddelta'' (liquid pressure not a function of vapor density)
441  J[2][2] = 0;
442  specified_parameter = CoolProp::iP;
443  break;
444  case saturation_PHSU_pure_options::IMPOSED_PV:
445  // dr_3/dtau
446  J[2][0] = SatV->first_partial_deriv(iP, iTau, iDelta) / specified_value;
447  // dr_3/ddelta' (vapor pressure not a function of liquid density)
448  J[2][1] = 0;
449  if (options.use_logdelta) {
450  // dr_3/d(log(delta'')
451  J[2][2] = deltaV * SatV->first_partial_deriv(iP, iDelta, iTau) / specified_value;
452  } else {
453  // dr_3/ddelta''
454  J[2][2] = SatV->first_partial_deriv(iP, iDelta, iTau) / specified_value;
455  }
456  specified_parameter = CoolProp::iP;
457  break;
458  case saturation_PHSU_pure_options::IMPOSED_HL:
459  // dr_3/dtau
460  J[2][0] = SatL->first_partial_deriv(iHmolar, iTau, iDelta);
461  // dr_3/ddelta'
462  J[2][1] = SatL->first_partial_deriv(iHmolar, iDelta, iTau);
463  if (options.use_logdelta) {
464  J[2][1] *= deltaL;
465  }
466  // dr_3/ddelta''
467  J[2][2] = 0; //(liquid enthalpy not a function of vapor density)
468  specified_parameter = CoolProp::iHmolar;
469  break;
470  case saturation_PHSU_pure_options::IMPOSED_HV:
471  // dr_3/dtau
472  J[2][0] = SatV->first_partial_deriv(iHmolar, iTau, iDelta);
473  // dr_3/ddelta'
474  J[2][1] = 0; //(vapor enthalpy not a function of liquid density)
475  // dr_3/ddelta''
476  J[2][2] = SatV->first_partial_deriv(iHmolar, iDelta, iTau);
477  if (options.use_logdelta) {
478  J[2][2] *= deltaV;
479  }
480  specified_parameter = CoolProp::iHmolar;
481  break;
482  case saturation_PHSU_pure_options::IMPOSED_SL:
483  // dr_3/dtau
484  J[2][0] = SatL->first_partial_deriv(iSmolar, iTau, iDelta);
485  // dr_3/ddelta'
486  J[2][1] = SatL->first_partial_deriv(iSmolar, iDelta, iTau);
487  if (options.use_logdelta) {
488  J[2][1] *= deltaL;
489  }
490  // dr_3/ddelta''
491  J[2][2] = 0; //(liquid entropy not a function of vapor density)
492  specified_parameter = CoolProp::iSmolar;
493  break;
494  case saturation_PHSU_pure_options::IMPOSED_SV:
495  // dr_3/dtau
496  J[2][0] = SatV->first_partial_deriv(iSmolar, iTau, iDelta);
497  // dr_3/ddelta'
498  J[2][1] = 0; //(vapor enthalpy not a function of liquid density)
499  // dr_3/ddelta''
500  J[2][2] = SatV->first_partial_deriv(iSmolar, iDelta, iTau);
501  if (options.use_logdelta) {
502  J[2][2] *= deltaV;
503  }
504  specified_parameter = CoolProp::iSmolar;
505  break;
506  default:
507  throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid", options.specified_variable));
508  }
509 
510  v = linsolve(J, negativer);
511 
512  // Conditions for an acceptable step are:
513  // a) tau > 1
514  // b) rhoL > rhoV or deltaL > deltaV
515  double tau0 = tau, deltaL0 = deltaL, deltaV0 = deltaV;
516  for (double omega_local = 1.0; omega_local > 0.1; omega_local /= 1.1) {
517  tau = tau0 + omega_local * options.omega * v[0];
518  if (options.use_logdelta) {
519  deltaL = exp(log(deltaL0) + omega_local * options.omega * v[1]);
520  deltaV = exp(log(deltaV0) + omega_local * options.omega * v[2]);
521  } else {
522  deltaL = deltaL0 + omega_local * options.omega * v[1];
523  deltaV = deltaV0 + omega_local * options.omega * v[2];
524  }
525  if (tau > 1 && deltaL > deltaV) {
526  break;
527  }
528  }
529 
530  rhoL = deltaL * reduce.rhomolar;
531  rhoV = deltaV * reduce.rhomolar;
532  T = reduce.T / tau;
533 
534  SatL->update(DmolarT_INPUTS, rhoL, T);
535  SatV->update(DmolarT_INPUTS, rhoV, T);
536 
537  error = sqrt(pow(negativer[0], 2) + pow(negativer[1], 2) + pow(negativer[2], 2));
538  iter++;
539  if (T < 0) {
540  throw SolutionError(format("saturation_PHSU_pure solver T < 0"));
541  }
542  // If the change is very small, stop
543  if (max_abs_value(v) < 1e-10) {
544  break;
545  }
546  if (iter > 50) {
547  // Set values back into the options structure for use in next solver
548  options.rhoL = rhoL;
549  options.rhoV = rhoV;
550  options.T = T;
551  // Error out
552  std::string info = get_parameter_information(specified_parameter, "short");
553  throw SolutionError(format("saturation_PHSU_pure solver did not converge after 50 iterations for %s=%Lg current error is %Lg",
554  info.c_str(), specified_value, error));
555  }
556  } while (error > 1e-9);
557  // Recalculate error
558  // The result has changed since the last error calculation.
559  // In rare scenarios, the final step can become unstable due to solving a singular
560  // J matrix. This final error check verifies that the solution is still good.
561  // Furthermore, the forced phase of SatL and SatV may have caused errors. We will recalculate them without this assumption.
562  SatL->unspecify_phase();
563  SatV->unspecify_phase();
564  SatL->update(DmolarT_INPUTS, rhoL, T);
565  SatV->update(DmolarT_INPUTS, rhoV, T);
566  negativer[0] = -(deltaV * (1 + deltaV * SatV->dalphar_dDelta()) - deltaL * (1 + deltaL * SatL->dalphar_dDelta()));
567  negativer[1] = -(deltaV * SatV->dalphar_dDelta() + SatV->alphar() + log(deltaV) - deltaL * SatL->dalphar_dDelta() - SatL->alphar() - log(deltaL));
568  switch (options.specified_variable) {
569  case saturation_PHSU_pure_options::IMPOSED_PL:
570  // -r_3 (equate calculated pressure and specified liquid pressure)
571  negativer[2] = -(SatL->p() / specified_value - 1);
572  break;
573  case saturation_PHSU_pure_options::IMPOSED_PV:
574  // -r_3 (equate calculated pressure and specified vapor pressure)
575  negativer[2] = -(SatV->p() / specified_value - 1);
576  break;
577  case saturation_PHSU_pure_options::IMPOSED_HL:
578  // -r_3 (equate calculated liquid enthalpy and specified liquid enthalpy)
579  negativer[2] = -(SatL->hmolar() - specified_value);
580  break;
581  case saturation_PHSU_pure_options::IMPOSED_HV:
582  // -r_3 (equate calculated vapor enthalpy and specified vapor enthalpy)
583  negativer[2] = -(SatV->hmolar() - specified_value);
584  break;
585  case saturation_PHSU_pure_options::IMPOSED_SL:
586  // -r_3 (equate calculated liquid entropy and specified liquid entropy)
587  negativer[2] = -(SatL->smolar() - specified_value);
588  break;
589  case saturation_PHSU_pure_options::IMPOSED_SV:
590  // -r_3 (equate calculated vapor entropy and specified vapor entropy)
591  negativer[2] = -(SatV->smolar() - specified_value);
592  break;
593  default:
594  throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid", options.specified_variable));
595  }
596  error = sqrt(pow(negativer[0], 2) + pow(negativer[1], 2) + pow(negativer[2], 2));
597  // reset the phase for the next update.
598  SatL->specify_phase(iphase_liquid);
599  SatV->specify_phase(iphase_gas);
600  if (error > 1e-8 && max_abs_value(v) > 1e-9){
601  throw SolutionError(format("saturation_PHSU_pure solver was good, but went bad. Current error is %Lg", error));
602  }
603 }
604 void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl rhomolar, saturation_D_pure_options& options)
605 {
606  /*
607  This function is inspired by the method of Akasaka:
608 
609  R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
610  Helmholtz Energy Equations of State",
611  Journal of Thermal Science and Technology v3 n3,2008
612 
613  Ancillary equations are used to get a sensible starting point
614  */
615  std::vector<CoolPropDbl> r(2, _HUGE), v;
616  std::vector<std::vector<CoolPropDbl>> J(2, std::vector<CoolPropDbl>(2, _HUGE));
617 
618  HEOS.calc_reducing_state();
619  const SimpleState& reduce = HEOS.get_reducing_state();
620  shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
621 
622  CoolPropDbl T, rhoL, rhoV;
623  CoolPropDbl deltaL = 0, deltaV = 0, tau = 0, error, p_error;
624  int iter = 0;
625 
626  // Use the density ancillary function as the starting point for the solver
627  try {
628  if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL) {
629  // Invert liquid density ancillary to get temperature
630  // TODO: fit inverse ancillaries too
631  T = HEOS.get_components()[0].ancillaries.rhoL.invert(rhomolar);
632  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
633  rhoL = rhomolar;
634  } else if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOV) {
635  // Invert vapor density ancillary to get temperature
636  // TODO: fit inverse ancillaries too
637  T = HEOS.get_components()[0].ancillaries.rhoV.invert(rhomolar);
638  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
639  rhoV = rhomolar;
640  } else {
641  throw ValueError(format("imposed rho to saturation_D_pure [%d%] is invalid", options.imposed_rho));
642  }
643 
644  deltaL = rhoL / reduce.rhomolar;
645  deltaV = rhoV / reduce.rhomolar;
646  tau = reduce.T / T;
647  } catch (NotImplementedError&) {
648  throw; // ??? What is this try...catch for?
649  }
650 
651  do {
652  /*if (get_debug_level()>8){
653  std::cout << format("%s:%d: right before the derivs with deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
654  }*/
655 
656  // Calculate once to save on calls to EOS
657  SatL->update(DmolarT_INPUTS, rhoL, T);
658  SatV->update(DmolarT_INPUTS, rhoV, T);
659 
660  CoolPropDbl pL = SatL->p();
661  CoolPropDbl pV = SatV->p();
662 
663  // These derivatives are needed for both cases
664  CoolPropDbl dalphar_dtauL = SatL->dalphar_dTau();
665  CoolPropDbl dalphar_dtauV = SatV->dalphar_dTau();
666  CoolPropDbl d2alphar_ddelta_dtauL = SatL->d2alphar_dDelta_dTau();
667  CoolPropDbl d2alphar_ddelta_dtauV = SatV->d2alphar_dDelta_dTau();
668  CoolPropDbl alpharL = SatL->alphar();
669  CoolPropDbl alpharV = SatV->alphar();
670  CoolPropDbl dalphar_ddeltaL = SatL->dalphar_dDelta();
671  CoolPropDbl dalphar_ddeltaV = SatV->dalphar_dDelta();
672 
673  // -r_1
674  r[0] = -(deltaV * (1 + deltaV * dalphar_ddeltaV) - deltaL * (1 + deltaL * dalphar_ddeltaL));
675  // -r_2
676  r[1] = -(deltaV * dalphar_ddeltaV + alpharV + log(deltaV) - deltaL * dalphar_ddeltaL - alpharL - log(deltaL));
677 
678  // dr1_dtau
679  J[0][0] = pow(deltaV, 2) * d2alphar_ddelta_dtauV - pow(deltaL, 2) * d2alphar_ddelta_dtauL;
680  // dr2_dtau
681  J[1][0] = deltaV * d2alphar_ddelta_dtauV + dalphar_dtauV - deltaL * d2alphar_ddelta_dtauL - dalphar_dtauL;
682 
683  if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL) {
684  CoolPropDbl d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
685  if (options.use_logdelta) {
686  J[0][1] = deltaV + 2 * pow(deltaV, 2) * dalphar_ddeltaV + pow(deltaV, 3) * d2alphar_ddelta2V;
687  J[1][1] = pow(deltaV, 2) * d2alphar_ddelta2V + 2 * deltaV * dalphar_ddeltaV + 1;
688  } else {
689  J[0][1] = 1 + 2 * deltaV * dalphar_ddeltaV + pow(deltaV, 2) * d2alphar_ddelta2V;
690  J[1][1] = deltaV * d2alphar_ddelta2V + 2 * dalphar_ddeltaV + 1 / deltaV;
691  }
692  } else if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOV) {
693  CoolPropDbl d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
694  if (options.use_logdelta) {
695  J[0][1] = -deltaL - 2 * pow(deltaL, 2) * dalphar_ddeltaL - pow(deltaL, 3) * d2alphar_ddelta2L;
696  J[1][1] = -pow(deltaL, 2) * d2alphar_ddelta2L - 2 * deltaL * dalphar_ddeltaL - 1;
697  } else {
698  J[0][1] = -1 - 2 * deltaL * dalphar_ddeltaL - pow(deltaL, 2) * d2alphar_ddelta2L;
699  J[1][1] = -deltaL * d2alphar_ddelta2L - 2 * dalphar_ddeltaL - 1 / deltaL;
700  }
701  }
702 
703  //double DET = J[0][0]*J[1][1]-J[0][1]*J[1][0];
704 
705  v = linsolve(J, r);
706 
707  double omega = options.omega;
708 
709  if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL) {
710  if (options.use_logdelta)
711  deltaV = exp(log(deltaV)+omega*v[1]);
712  else
713  {
714  if (deltaV + omega*v[1] <= 0) {omega = 0.5*deltaV/v[1];} // gone off track, take a smaller step
715  if (tau + omega*v[0] <= 0) {omega = 0.5*tau/v[0];}
716  deltaV += omega*v[1];
717  }
718  }
719  else
720  {
721  if (options.use_logdelta)
722  deltaL = exp(log(deltaL)+omega*v[1]);
723  else
724  {
725  if (deltaL + omega*v[1] <= 0) {omega = 0.5*deltaL/v[1];} // gone off track, take a smaller step
726  if (tau + omega*v[0] <= 0) {omega = 0.5*tau/v[0];}
727  deltaL += omega*v[1];
728  }
729  }
730 
731  tau += omega*v[0];
732 
733  rhoL = deltaL*reduce.rhomolar;
734  rhoV = deltaV*reduce.rhomolar;
735  T = reduce.T/tau;
736 
737  p_error = (pL-pV)/pL;
738 
739  error = sqrt(pow(r[0], 2) + pow(r[1], 2));
740  iter++;
741  if (T < 0) {
742  throw SolutionError(format("saturation_D_pure solver T < 0"));
743  }
744  if (iter > options.max_iterations){
745  throw SolutionError(format("saturation_D_pure solver did not converge after %d iterations with rho: %g mol/m^3",options.max_iterations,rhomolar));
746  }
747  } while (error > 1e-9);
748  CoolPropDbl p_error_limit = 1e-3;
749  if (std::abs(p_error) > p_error_limit) {
750  throw SolutionError(format("saturation_D_pure solver abs error on p [%Lg] > limit [%Lg]", p_error, p_error_limit));
751  }
752 }
753 void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options) {
754  // Set some input options
755  SaturationSolvers::saturation_T_pure_Akasaka_options _options(false);
756  _options.omega = 1.0;
757  try {
758  // Actually call the solver
759  SaturationSolvers::saturation_T_pure_Maxwell(HEOS, T, _options);
760  } catch (...) {
761  try {
762  // Actually call the solver
763  SaturationSolvers::saturation_T_pure_Akasaka(HEOS, T, _options);
764  } catch (...) {
765  // If there was an error, store values for use in later solvers
766  options.pL = _options.pL;
767  options.pV = _options.pV;
768  options.rhoL = _options.rhoL;
769  options.rhoV = _options.rhoV;
770  options.p = _options.pL;
771  SaturationSolvers::saturation_T_pure_1D_P(HEOS, T, options);
772  }
773  }
774 }
775 void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options) {
776  // Start with the method of Akasaka
777 
778  /*
779  This function implements the method of Akasaka
780 
781  R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
782  Helmholtz Energy Equations of State",
783  Journal of Thermal Science and Technology v3 n3,2008
784 
785  Ancillary equations are used to get a sensible starting point
786  */
787 
788  HEOS.calc_reducing_state();
789  const SimpleState& reduce = HEOS.get_reducing_state();
790  CoolPropDbl R_u = HEOS.gas_constant();
791  shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
792 
793  CoolPropDbl rhoL = _HUGE, rhoV = _HUGE, JL, JV, KL, KV, dJL, dJV, dKL, dKV;
794  CoolPropDbl DELTA, deltaL = 0, deltaV = 0, error, PL, PV, stepL, stepV;
795  int iter = 0;
796 
797  try {
798  if (options.use_guesses) {
799  // Use the guesses provided in the options structure
800  rhoL = options.rhoL;
801  rhoV = options.rhoV;
802  } else {
803  // Use the density ancillary function as the starting point for the solver
804 
805  // If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
806  if (T > 0.99 * HEOS.get_reducing_state().T) {
807  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T - 0.1);
808  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T - 0.1);
809  } else {
810  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
811  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
812 
813  // Apply a single step of Newton's method to improve guess value for liquid
814  // based on the error between the gas pressure (which is usually very close already)
815  // and the liquid pressure, which can sometimes (especially at low pressure),
816  // be way off, and often times negative
817  SatL->update(DmolarT_INPUTS, rhoL, T);
818  SatV->update(DmolarT_INPUTS, rhoV, T);
819 
820  // Update the guess for liquid density using density solver with vapor pressure
821  // and liquid density guess from ancillaries
822  HEOS.specify_phase(iphase_liquid);
823  rhoL = HEOS.solver_rho_Tp(T, SatV->p(), rhoL);
824  HEOS.unspecify_phase();
825  }
826  }
827 
828  deltaL = rhoL / reduce.rhomolar;
829  deltaV = rhoV / reduce.rhomolar;
830  } catch (NotImplementedError&) {
831  /*double Tc = crit.T;
832  double pc = crit.p.Pa;
833  double w = 6.67228479e-09*Tc*Tc*Tc-7.20464352e-06*Tc*Tc+3.16947758e-03*Tc-2.88760012e-01;
834  double q = -6.08930221451*w -5.42477887222;
835  double pt = exp(q*(Tc/T-1))*pc;*/
836 
837  //double rhoL = density_Tp_Soave(T, pt, 0), rhoV = density_Tp_Soave(T, pt, 1);
838 
839  //deltaL = rhoL/reduce.rhomolar;
840  //deltaV = rhoV/reduce.rhomolar;
841  //tau = reduce.T/T;
842  }
843  //if (get_debug_level()>5){
844  // std::cout << format("%s:%d: Akasaka guess values deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
845  // }
846 
847  do {
848  /*if (get_debug_level()>8){
849  std::cout << format("%s:%d: right before the derivs with deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
850  }*/
851 
852  // Calculate once to save on calls to EOS
853  SatL->update(DmolarT_INPUTS, rhoL, T);
854  SatV->update(DmolarT_INPUTS, rhoV, T);
855  CoolPropDbl alpharL = SatL->alphar();
856  CoolPropDbl alpharV = SatV->alphar();
857  CoolPropDbl dalphar_ddeltaL = SatL->dalphar_dDelta();
858  CoolPropDbl dalphar_ddeltaV = SatV->dalphar_dDelta();
859  CoolPropDbl d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
860  CoolPropDbl d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
861 
862  JL = deltaL * (1 + deltaL * dalphar_ddeltaL);
863  JV = deltaV * (1 + deltaV * dalphar_ddeltaV);
864  KL = deltaL * dalphar_ddeltaL + alpharL + log(deltaL);
865  KV = deltaV * dalphar_ddeltaV + alpharV + log(deltaV);
866 
867  PL = R_u * reduce.rhomolar * T * JL;
868  PV = R_u * reduce.rhomolar * T * JV;
869 
870  // At low pressure, the magnitude of d2alphar_ddelta2L and d2alphar_ddelta2V are enormous, truncation problems arise for all the partials
871  dJL = 1 + 2 * deltaL * dalphar_ddeltaL + deltaL * deltaL * d2alphar_ddelta2L;
872  dJV = 1 + 2 * deltaV * dalphar_ddeltaV + deltaV * deltaV * d2alphar_ddelta2V;
873  dKL = 2 * dalphar_ddeltaL + deltaL * d2alphar_ddelta2L + 1 / deltaL;
874  dKV = 2 * dalphar_ddeltaV + deltaV * d2alphar_ddelta2V + 1 / deltaV;
875 
876  DELTA = dJV * dKL - dJL * dKV;
877 
878  error = sqrt((KL - KV) * (KL - KV) + (JL - JV) * (JL - JV));
879 
880  // Get the predicted step
881  stepL = options.omega / DELTA * ((KV - KL) * dJV - (JV - JL) * dKV);
882  stepV = options.omega / DELTA * ((KV - KL) * dJL - (JV - JL) * dKL);
883 
884  CoolPropDbl deltaL0 = deltaL, deltaV0 = deltaV;
885  // Conditions for an acceptable step are:
886  // a) rhoL > rhoV or deltaL > deltaV
887  for (double omega_local = 1.0; omega_local > 0.1; omega_local /= 1.1) {
888  deltaL = deltaL0 + omega_local * stepL;
889  deltaV = deltaV0 + omega_local * stepV;
890 
891  if (deltaL > 1 && deltaV < 1 && deltaV > 0) {
892  break;
893  }
894  }
895 
896  rhoL = deltaL * reduce.rhomolar;
897  rhoV = deltaV * reduce.rhomolar;
898  iter++;
899  if (iter > 100) {
900  throw SolutionError(format("Akasaka solver did not converge after 100 iterations"));
901  }
902  } while (error > 1e-10 && std::abs(stepL) > 10 * DBL_EPSILON * std::abs(stepL) && std::abs(stepV) > 10 * DBL_EPSILON * std::abs(stepV));
903 
904  CoolPropDbl p_error_limit = 1e-3;
905  CoolPropDbl p_error = (PL - PV) / PL;
906  if (std::abs(p_error) > p_error_limit) {
907  options.pL = PL;
908  options.pV = PV;
909  options.rhoL = rhoL;
910  options.rhoV = rhoV;
911  throw SolutionError(format("saturation_T_pure_Akasaka solver abs error on p [%g] > limit [%g]", std::abs(p_error), p_error_limit));
912  }
913 }
914 
915 CoolPropDbl sign(CoolPropDbl x) {
916  if (x > 0) {
917  return 1;
918  } else {
919  return -1;
920  }
921 }
922 
923 void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options) {
924 
925  /*
926  This function implements the method of
927 
928  Ancillary equations are used to get a sensible starting point
929  */
930 
931  HEOS.calc_reducing_state();
932  shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
933  CoolProp::SimpleState& crit = HEOS.get_components()[0].crit;
934  CoolPropDbl rhoL = _HUGE, rhoV = _HUGE, error = 999, DeltavL, DeltavV, pL, pV, p, last_error;
935  int iter = 0, small_step_count = 0,
936  backwards_step_count = 0; // Counter for the number of times you have taken a step that increases error
937 
938  try {
939  if (options.use_guesses) {
940  // Use the guesses provided in the options structure
941  rhoL = options.rhoL;
942  rhoV = options.rhoV;
943  } else {
944  // Use the density ancillary function as the starting point for the solver
945 
946  // If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
947  if (T > 0.9999 * HEOS.get_reducing_state().T) {
948  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T - 0.1);
949  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T - 0.1);
950  } else {
951  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
952  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
953  p = HEOS.get_components()[0].ancillaries.pV.evaluate(T);
954 
955  CoolProp::SimpleState& tripleL = HEOS.get_components()[0].triple_liquid;
956  CoolProp::SimpleState& tripleV = HEOS.get_components()[0].triple_vapor;
957 
958  // If the guesses are terrible, apply a simple correction
959  // but only if the limits are being checked
960  if ((rhoL < crit.rhomolar * 0.8 || rhoL > tripleL.rhomolar * 1.2 || rhoV > crit.rhomolar * 1.2 || rhoV < tripleV.rhomolar * 0.8)
961  && !get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
962 
963  if (get_debug_level() > 5) {
964  std::cout << format("[Maxwell] ancillaries correction T: %0.16Lg rhoL: %0.16Lg rhoV: %0.16Lg rhoc: %g rhoLtrip: %g rhoVtrip: %g\n", T, rhoL, rhoV, crit.rhomolar, tripleL.rhomolar, tripleV.rhomolar );
965  }
966 
967  // Lets assume that liquid density is more or less linear with T
968  rhoL = (crit.rhomolar - tripleL.rhomolar) / (crit.T - tripleL.T) * (T - tripleL.T) + tripleL.rhomolar;
969  // Then we calculate pressure from this density
970  SatL->update_DmolarT_direct(rhoL, T);
971  // Then we assume vapor to be ideal gas
972  if (SatL->p() > 0) {
973  rhoV = SatL->p() / (SatL->gas_constant() * T);
974  } else {
975  rhoV = p / (SatL->gas_constant() * T);
976  }
977  // Update the vapor state
978  SatV->update_DmolarT_direct(rhoV, T);
979  } else {
980  SatL->update_DmolarT_direct(rhoL, T);
981  SatV->update_DmolarT_direct(rhoV, T);
982  }
983  if (get_debug_level() > 5) {
984  std::cout << format("[Maxwell] ancillaries T: %0.16Lg rhoL: %0.16Lg rhoV: %0.16Lg pL: %g pV: %g\n", T, rhoL, rhoV, SatL->p(),
985  SatV->p());
986  }
987 
988  // Update the guess for liquid density using density solver with vapor pressure
989  // and liquid density guess from ancillaries, but only if the pressures are not
990  // close to each other
991  if (std::abs((SatL->p() - p) / p) > 0.1) {
992  for (int iii = 0; iii < 6; ++iii) {
993  // Use Halley's method to update the liquid density (http://en.wikipedia.org/wiki/Halley%27s_method)
994  CoolPropDbl f = SatL->p() - SatV->p();
995  CoolPropDbl dpdrho = SatL->first_partial_deriv(iP, iDmolar, iT);
996  CoolPropDbl d2pdrho2 = SatL->second_partial_deriv(iP, iDmolar, iT, iDmolar, iT);
997  CoolPropDbl deltarhoLHalley = -(2 * f * dpdrho) / (2 * POW2(dpdrho) - f * d2pdrho2);
998  rhoL += deltarhoLHalley;
999  if (std::abs(deltarhoLHalley / rhoL) < DBL_EPSILON) {
1000  break;
1001  }
1002  SatL->update_DmolarT_direct(rhoL, T);
1003  }
1004  }
1005 
1006  SatL->update_DmolarT_direct(rhoL, T);
1007  SatV->update_DmolarT_direct(rhoV, T);
1008 
1009  // Update the guess for vapor density using density solver with vapor pressure
1010  // and density guess from ancillaries, but only if the pressures are not
1011  // close to each other
1012  if (std::abs((SatV->p() - p) / p) > 0.1) {
1013  HEOS.specify_phase(iphase_gas);
1014  rhoV = SatV->solver_rho_Tp(T, p, rhoV);
1015  HEOS.unspecify_phase();
1016  }
1017  }
1018  }
1019  } catch (NotImplementedError&) {
1020  }
1021 
1022  if (rhoL < crit.rhomolar) {
1023  rhoL = 1.01 * crit.rhomolar;
1024  }
1025  if (rhoV > crit.rhomolar) {
1026  rhoV = 0.99 * crit.rhomolar;
1027  }
1028  last_error = _HUGE;
1029  SatL->update_DmolarT_direct(rhoL, T);
1030  SatV->update_DmolarT_direct(rhoV, T);
1031  if (get_debug_level() > 5) {
1032  std::cout << format("[Maxwell] starting T: %0.16Lg rhoL: %Lg rhoV: %Lg pL: %Lg pV: %g\n", T, rhoL, rhoV, SatL->p(), SatV->p());
1033  }
1034  do {
1035  pL = SatL->p();
1036  pV = SatV->p();
1037  CoolPropDbl vL = 1 / SatL->rhomolar(), vV = 1 / SatV->rhomolar();
1038  // Get alpha, the pressure derivative with volume at constant T
1039  // Given by (dp/drho|T)*drhodv
1040  CoolPropDbl alphaL = SatL->first_partial_deriv(iP, iDmolar, iT) * (-POW2(SatL->rhomolar()));
1041  CoolPropDbl alphaV = SatV->first_partial_deriv(iP, iDmolar, iT) * (-POW2(SatV->rhomolar()));
1042 
1043  // Total helmholtz energy for liquid and vapor
1044  CoolPropDbl RT = SatL->gas_constant() * T;
1045  CoolPropDbl helmholtzL = (SatL->calc_alpha0() + SatL->calc_alphar()) * RT;
1046  CoolPropDbl helmholtzV = (SatV->calc_alpha0() + SatV->calc_alphar()) * RT;
1047 
1048  // Calculate the mean pressure
1049  CoolPropDbl pM = (helmholtzL - helmholtzV) / (vV - vL);
1050 
1051  // Coefficients for the quadratic in the step
1052  CoolPropDbl A = 0.5 * alphaL * (alphaL - alphaV);
1053  CoolPropDbl B = alphaL * (pL - pV - alphaV * (vL - vV));
1054  CoolPropDbl C = alphaV * (vL - vV) * (pM - pL) + 0.5 * POW2(pL - pV);
1055 
1056  // Argument to square root
1057  CoolPropDbl sqrt_arg = std::abs(B * B / (4 * A * A) - C / A);
1058 
1059  // If the argument to sqrt is very small, we multiply it by a large factor to make it
1060  // larger, and then also divide the sqrt by the sqrt of the factor
1061  if (std::abs(sqrt_arg) > 1e-10) {
1062  DeltavL = -0.5 * B / A + sign((alphaL - alphaV) / alphaV) * sqrt(sqrt_arg);
1063  } else {
1064  // Scale the argument to sqrt() function to make it about 1.0, and divide by sqrt(factor) to yield a factor of 1
1065  CoolPropDbl powerfactor = -log10(sqrt_arg);
1066  DeltavL = -0.5 * B / A + sign((alphaL - alphaV) / alphaV) * sqrt(sqrt_arg * powerfactor) / sqrt(powerfactor);
1067  }
1068  DeltavV = (pL - pV + alphaL * DeltavL) / alphaV;
1069 
1070  // Update the densities of liquid and vapor
1071  rhoL = 1 / (vL + DeltavL);
1072  rhoV = 1 / (vV + DeltavV);
1073  if (B * B / (4 * A * A) - C / A < -10 * DBL_EPSILON) {
1074  rhoL *= 1.01;
1075  rhoV /= 1.01;
1076  }
1077 
1078  // Update the states again
1079  SatL->update_DmolarT_direct(rhoL, T);
1080  SatV->update_DmolarT_direct(rhoV, T);
1081 
1082  // Calculate the error (here the relative error in pressure)
1083  error = std::abs((SatL->p() - SatV->p()) / SatL->p());
1084 
1085  if (get_debug_level() > 5) {
1086  std::cout << format("[Maxwell] rhoL: %0.16Lg rhoV: %0.16Lg error: %Lg dvL/vL: %Lg dvV/vV: %Lg pL: %Lg pV: %Lg\n", rhoL, rhoV, error,
1087  DeltavL / vL, DeltavV / vV, pL, pV);
1088  }
1089 
1090  // If the step size is small, start a counter to allow the other density
1091  // to be corrected a few times
1092  if (std::abs(DeltavL * rhoL) < 1e-13 || std::abs(DeltavV * rhoV) < 1e-13) {
1093  small_step_count++;
1094  }
1095  // If you are not continuing to march towards the solution, after a couple of times, stop
1096  // This is especially a problem for water
1097  if (std::abs(error) > std::abs(last_error)) {
1098  backwards_step_count++;
1099  }
1100 
1101  iter++;
1102  last_error = error;
1103  if (iter > 30) {
1104  throw SolutionError(format("Maxwell solver did not converge after 30 iterations; rhoL: %0.16Lg rhoV: %0.16Lg error: %Lg dvL/vL: %Lg "
1105  "dvV/vV: %Lg pL: %Lg pV: %Lg\n",
1106  rhoL, rhoV, error, DeltavL / vL, DeltavV / vV, pL, pV));
1107  }
1108  } while ((SatL->p() < 0) || (error > 1e-10 && small_step_count < 4 && backwards_step_count < 6));
1109  if (get_debug_level() > 5) {
1110  std::cout << format("[Maxwell] pL: %g pV: %g\n", SatL->p(), SatV->p());
1111  }
1112 }
1113 
1114 void SaturationSolvers::x_and_y_from_K(CoolPropDbl beta, const std::vector<CoolPropDbl>& K, const std::vector<CoolPropDbl>& z,
1115  std::vector<CoolPropDbl>& x, std::vector<CoolPropDbl>& y) {
1116  for (unsigned int i = 0; i < K.size(); i++) {
1117  double denominator = (1 - beta + beta * K[i]); // Common denominator
1118  x[i] = z[i] / denominator;
1119  y[i] = K[i] * z[i] / denominator;
1120  }
1121 }
1122 
1123 void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend& HEOS, const CoolPropDbl beta, CoolPropDbl T, CoolPropDbl p,
1124  const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& K, mixture_VLE_IO& options) {
1125  int iter = 1;
1126  CoolPropDbl change, f, df, deriv_liq, deriv_vap;
1127  std::size_t N = z.size();
1128  std::vector<CoolPropDbl> ln_phi_liq, ln_phi_vap;
1129  ln_phi_liq.resize(N);
1130  ln_phi_vap.resize(N);
1131 
1132  std::vector<CoolPropDbl>&x = HEOS.SatL->get_mole_fractions_ref(), &y = HEOS.SatV->get_mole_fractions_ref();
1133  x_and_y_from_K(beta, K, z, x, y);
1134 
1135  HEOS.SatL->specify_phase(iphase_liquid);
1136  HEOS.SatV->specify_phase(iphase_gas);
1137 
1138  normalize_vector(x);
1139  normalize_vector(y);
1140 
1141  HEOS.SatL->set_mole_fractions(x);
1142  HEOS.SatV->set_mole_fractions(y);
1143  HEOS.SatL->calc_reducing_state();
1144  HEOS.SatV->calc_reducing_state();
1145  CoolPropDbl rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(T, p, iphase_liquid); // [mol/m^3]
1146  CoolPropDbl rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(T, p, iphase_gas); // [mol/m^3]
1147 
1148  // Use Peneloux volume translation to shift liquid volume
1149  // As in Horstmann :: doi:10.1016/j.fluid.2004.11.002
1150  double summer_c = 0, v_SRK = 1 / rhomolar_liq;
1151  const std::vector<CoolPropFluid>& components = HEOS.get_components();
1152  for (std::size_t i = 0; i < components.size(); ++i) {
1153  // Get the parameters for the cubic EOS
1154  CoolPropDbl Tc = HEOS.get_fluid_constant(i, iT_critical);
1155  CoolPropDbl pc = HEOS.get_fluid_constant(i, iP_critical);
1156  CoolPropDbl rhomolarc = HEOS.get_fluid_constant(i, irhomolar_critical);
1157  CoolPropDbl R = 8.3144598;
1158 
1159  summer_c += z[i] * (0.40768 * R * Tc / pc * (0.29441 - pc / (rhomolarc * R * Tc)));
1160  }
1161  rhomolar_liq = 1 / (v_SRK - summer_c);
1162  HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq);
1163  HEOS.SatV->update_TP_guessrho(T, p, rhomolar_vap);
1164 
1165  do {
1166  HEOS.SatL->update_TP_guessrho(T, p, HEOS.SatL->rhomolar());
1167  HEOS.SatV->update_TP_guessrho(T, p, HEOS.SatV->rhomolar());
1168 
1169  f = 0;
1170  df = 0;
1171 
1172  x_N_dependency_flag xN_flag = XN_INDEPENDENT;
1173  for (std::size_t i = 0; i < N; ++i) {
1174  ln_phi_liq[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatL.get()), i, xN_flag);
1175  ln_phi_vap[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatV.get()), i, xN_flag);
1176 
1177  if (options.sstype == imposed_p) {
1178  deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatL.get()), i, xN_flag);
1179  deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatV.get()), i, xN_flag);
1180  } else if (options.sstype == imposed_T) {
1181  deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatL.get()), i, xN_flag);
1182  deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatV.get()), i, xN_flag);
1183  } else {
1184  throw ValueError();
1185  }
1186 
1187  K[i] = exp(ln_phi_liq[i] - ln_phi_vap[i]);
1188 
1189  f += z[i] * (K[i] - 1) / (1 - beta + beta * K[i]);
1190 
1191  double dfdK = K[i] * z[i] / pow(1 - beta + beta * K[i], (int)2);
1192  df += dfdK * (deriv_liq - deriv_vap);
1193  }
1194 
1195  if (std::abs(df) <= 1e-14) { // To avoid dividing by 0
1196  if (std::abs(f) <= 1e-12) // 1e-12 is the loop convergence criterion
1197  {
1198  change = -f; // Should be converged. f <= e-12, so change will have nearly no impact.
1199  } else {
1200  throw ValueError(format("df very small (df = %g) in successive_substitution but f is not converged (f = %g > 1e-12).", df, f));
1201  }
1202  } else {
1203  change = -f / df;
1204  }
1205 
1206  double omega = 1.0;
1207  if (options.sstype == imposed_p) {
1208  T += change;
1209  } else if (options.sstype == imposed_T) {
1210  if (std::abs(change) > 0.05 * p) {
1211  omega = 0.1;
1212  }
1213  p += omega * change;
1214  }
1215 
1216  x_and_y_from_K(beta, K, z, x, y);
1217  normalize_vector(x);
1218  normalize_vector(y);
1219  HEOS.SatL->set_mole_fractions(x);
1220  HEOS.SatV->set_mole_fractions(y);
1221 
1222  iter += 1;
1223  if (iter > 50) {
1224  throw ValueError(format("saturation_p was unable to reach a solution within 50 iterations"));
1225  }
1226  } while (std::abs(f) > 1e-12 && iter < options.Nstep_max);
1227 
1228  HEOS.SatL->update_TP_guessrho(T, p, HEOS.SatL->rhomolar());
1229  HEOS.SatV->update_TP_guessrho(T, p, HEOS.SatV->rhomolar());
1230 
1231  options.p = HEOS.SatL->p();
1232  options.T = HEOS.SatL->T();
1233  options.rhomolar_liq = HEOS.SatL->rhomolar();
1234  options.rhomolar_vap = HEOS.SatV->rhomolar();
1235  options.x = x;
1236  options.y = y;
1237 }
1238 void SaturationSolvers::newton_raphson_saturation::resize(std::size_t N) {
1239  this->N = N;
1240  x.resize(N);
1241  y.resize(N);
1242 
1243  if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1244  r.resize(N + 1);
1245  err_rel.resize(N + 1);
1246  J.resize(N + 1, N + 1);
1247  } else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1248  r.resize(N);
1249  err_rel.resize(N);
1250  J.resize(N, N);
1251  } else {
1252  throw ValueError();
1253  }
1254 }
1255 void SaturationSolvers::newton_raphson_saturation::check_Jacobian() {
1256  // References to the classes for concision
1257  HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
1258 
1259  // Build the Jacobian and residual vectors
1260  build_arrays();
1261 
1262  // Make copies of the base
1263  CoolPropDbl T0 = T;
1264  std::vector<CoolPropDbl> x0 = x;
1265  Eigen::VectorXd r0 = r;
1266  Eigen::MatrixXd J0 = J;
1267  CoolPropDbl rhomolar_liq0 = rSatL.rhomolar();
1268  CoolPropDbl rhomolar_vap0 = rSatV.rhomolar();
1269 
1270  {
1271  // Derivatives with respect to T
1272  double dT = 1e-3, T1 = T + dT, T2 = T - dT;
1273  this->T = T1;
1274  this->rhomolar_liq = rhomolar_liq0;
1275  this->rhomolar_vap = rhomolar_vap0;
1276  build_arrays();
1277  Eigen::VectorXd r1 = r;
1278  this->T = T2;
1279  this->rhomolar_liq = rhomolar_liq0;
1280  this->rhomolar_vap = rhomolar_vap0;
1281  build_arrays();
1282  Eigen::VectorXd r2 = r;
1283 
1284  Eigen::VectorXd diffn = (r1 - r2) / (2 * dT);
1285  std::cout << format("For T\n");
1286  //std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
1287  //std::cout << "analytic: " << vec_to_string(J0.col(N-1), "%0.11Lg") << std::endl;
1288  }
1289  {
1290  // Derivatives with respect to rho'
1291  double drho = 1;
1292  this->T = T0;
1293  this->rhomolar_liq = rhomolar_liq0 + drho;
1294  this->rhomolar_vap = rhomolar_vap0;
1295  build_arrays();
1296  Eigen::VectorXd rr1 = r;
1297  this->T = T0;
1298  this->rhomolar_liq = rhomolar_liq0 - drho;
1299  this->rhomolar_vap = rhomolar_vap0;
1300  build_arrays();
1301  Eigen::VectorXd rr2 = r;
1302 
1303  Eigen::VectorXd diffn = (rr1 - rr2) / (2 * drho);
1304  std::cout << format("For rho\n");
1305  //std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
1306  //std::cout << "analytic: " << vec_to_string(J0.col(N-1), "%0.11Lg") << std::endl;
1307  }
1308  for (std::size_t i = 0; i < x.size() - 1; ++i) {
1309  // Derivatives with respect to x[i]
1310  double dx = 1e-5;
1311  this->x = x0;
1312  this->x[i] += dx;
1313  this->x[x.size() - 1] -= dx;
1314  this->T = T0;
1315  this->rhomolar_liq = rhomolar_liq0;
1316  this->rhomolar_vap = rhomolar_vap0;
1317  build_arrays();
1318  Eigen::VectorXd r1 = this->r;
1319 
1320  this->x = x0;
1321  this->x[i] -= dx;
1322  this->x[x.size() - 1] += dx;
1323  this->T = T0;
1324  this->rhomolar_liq = rhomolar_liq0;
1325  this->rhomolar_vap = rhomolar_vap0;
1326  build_arrays();
1327  Eigen::VectorXd r2 = this->r;
1328 
1329  Eigen::VectorXd diffn = (r1 - r2) / (2 * dx);
1330  std::cout << format("For x%d N %d\n", i, N);
1331  //std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
1332  //std::cout << "analytic: " << vec_to_string(J0.col(i), "%0.11Lg") << std::endl;
1333  }
1334 }
1335 void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBackend& HEOS, const std::vector<CoolPropDbl>& z,
1336  std::vector<CoolPropDbl>& z_incipient, newton_raphson_saturation_options& IO) {
1337  int iter = 0;
1338  bool debug = get_debug_level() > 9 || false;
1339 
1340  if (debug) {
1341  std::cout << " NRsat::call: p " << IO.p << " T " << IO.T << " dl " << IO.rhomolar_liq << " dv " << IO.rhomolar_vap << std::endl;
1342  }
1343 
1344  // Reset all the variables and resize
1345  pre_call();
1346 
1347  this->bubble_point = IO.bubble_point;
1348  rhomolar_liq = IO.rhomolar_liq;
1349  rhomolar_vap = IO.rhomolar_vap;
1350  T = IO.T;
1351  p = IO.p;
1352  imposed_variable = IO.imposed_variable;
1353 
1354  resize(z.size());
1355 
1356  if (bubble_point) {
1357  // Bubblepoint, vapor (y) is the incipient phase
1358  x = z;
1359  y = z_incipient;
1360  } else {
1361  // Dewpoint, liquid (x) is the incipient phase
1362  y = z;
1363  x = z_incipient;
1364  }
1365 
1366  // Hold a pointer to the backend
1367  this->HEOS = &HEOS;
1368 
1369  //check_Jacobian();
1370 
1371  do {
1372  // Build the Jacobian and residual vectors
1373  build_arrays();
1374 
1375  // Solve for the step; v is the step with the contents
1376  // [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]
1377  Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
1378 
1379  if (bubble_point) {
1380  for (unsigned int i = 0; i < N - 1; ++i) {
1381  err_rel[i] = v[i] / y[i];
1382  y[i] += v[i];
1383  }
1384  y[N - 1] = 1 - std::accumulate(y.begin(), y.end() - 1, 0.0);
1385  } else {
1386  for (unsigned int i = 0; i < N - 1; ++i) {
1387  err_rel[i] = v[i] / x[i];
1388  x[i] += v[i];
1389  }
1390  x[N - 1] = 1 - std::accumulate(x.begin(), x.end() - 1, 0.0);
1391  }
1392  if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED) {
1393  T += v[N - 1];
1394  err_rel[N - 1] = v[N - 1] / T;
1395  } else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1396  p += v[N - 1];
1397  err_rel[N - 1] = v[N - 1] / p;
1398  } else if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1399  T += v[N - 1];
1400  err_rel[N - 1] = v[N - 1] / T;
1401  rhomolar_liq += v[N];
1402  err_rel[N] = v[N] / rhomolar_liq;
1403  } else {
1404  throw ValueError("invalid imposed_variable");
1405  }
1406  if (debug) {
1407  //std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " v " << vec_to_string(v, "%0.10Lg") << " x " << vec_to_string(x, "%0.10Lg") << " r " << vec_to_string(r, "%0.10Lg") << std::endl;
1408  }
1409 
1410  min_rel_change = err_rel.cwiseAbs().minCoeff();
1411  iter++;
1412 
1413  if (iter == IO.Nstep_max) {
1414  throw ValueError(format("newton_raphson_saturation::call reached max number of iterations [%d]", IO.Nstep_max));
1415  }
1416  } while (this->error_rms > 1e-7 && min_rel_change > 1000 * DBL_EPSILON && iter < IO.Nstep_max);
1417 
1418  IO.Nsteps = iter;
1419  IO.p = p;
1420  IO.x = x; // Mole fractions in liquid
1421  IO.y = y; // Mole fractions in vapor
1422  IO.T = T;
1423  IO.rhomolar_liq = rhomolar_liq;
1424  IO.rhomolar_vap = rhomolar_vap;
1425  const std::vector<CoolPropFluid>& fluidsL = HEOS.SatL->get_components();
1426  const std::vector<CoolPropFluid>& fluidsV = HEOS.SatV->get_components();
1427  if (!fluidsL.empty() && !fluidsV.empty()) {
1428  IO.hmolar_liq = HEOS.SatL->hmolar();
1429  IO.hmolar_vap = HEOS.SatV->hmolar();
1430  IO.smolar_liq = HEOS.SatL->smolar();
1431  IO.smolar_vap = HEOS.SatV->smolar();
1432  }
1433 }
1434 
1435 void SaturationSolvers::newton_raphson_saturation::build_arrays() {
1436  // References to the classes for concision
1437  HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
1438 
1439  // Step 0:
1440  // -------
1441  // Set mole fractions for the incipient phase
1442  if (bubble_point) {
1443  // Vapor is incipient phase, set its composition
1444  rSatV.set_mole_fractions(y);
1445  rSatL.set_mole_fractions(x);
1446  } else {
1447  // Liquid is incipient phase, set its composition
1448  rSatL.set_mole_fractions(x);
1449  rSatV.set_mole_fractions(y);
1450  }
1451 
1452  if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1453  rSatL.update(DmolarT_INPUTS, rhomolar_liq, T);
1454  rSatV.update(DmolarT_INPUTS, rhomolar_vap, T);
1455  } else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1456  rSatL.update_TP_guessrho(T, p, rhomolar_liq);
1457  rhomolar_liq = rSatL.rhomolar();
1458  rSatV.update_TP_guessrho(T, p, rhomolar_vap);
1459  rhomolar_vap = rSatV.rhomolar();
1460  } else {
1461  throw ValueError("imposed variable not set for NR VLE");
1462  }
1463 
1464  // For diagnostic purposes calculate the pressures (no derivatives are evaluated)
1465  CoolPropDbl p_liq = rSatL.p();
1466  CoolPropDbl p_vap = rSatV.p();
1467  p = 0.5 * (p_liq + p_vap);
1468 
1469  // Step 2:
1470  // -------
1471  // Build the residual vector and the Jacobian matrix
1472 
1473  x_N_dependency_flag xN_flag = XN_DEPENDENT;
1474 
1475  if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1476  // For the residuals F_i (equality of fugacities)
1477  for (std::size_t i = 0; i < N; ++i) {
1478  // Equate the liquid and vapor fugacities
1479  CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1480  CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1481  r(i) = ln_f_liq - ln_f_vap;
1482 
1483  for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1484  if (bubble_point) {
1485  J(i, j) = -MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatV, i, j, xN_flag);
1486  } else {
1487  J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatL, i, j, xN_flag);
1488  }
1489  }
1490  J(i, N - 1) = MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatL, i, xN_flag)
1491  - MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatV, i, xN_flag);
1492  J(i, N) = MixtureDerivatives::dln_fugacity_i_drho__constT_n(rSatL, i, xN_flag);
1493  }
1494  // ---------------------------------------------------------------
1495  // Derivatives of pL(T,rho',x)-p(T,rho'',y) with respect to inputs
1496  // ---------------------------------------------------------------
1497  r(N) = p_liq - p_vap;
1498  for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1499  J(N, j) = MixtureDerivatives::dpdxj__constT_V_xi(rSatL, j, xN_flag); // p'' not a function of x0
1500  }
1501  // Fixed composition derivatives
1502  J(N, N - 1) = rSatL.first_partial_deriv(iP, iT, iDmolar) - rSatV.first_partial_deriv(iP, iT, iDmolar);
1503  J(N, N) = rSatL.first_partial_deriv(iP, iDmolar, iT);
1504  } else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED) {
1505  // Independent variables are N-1 mole fractions of incipient phase and T
1506 
1507  // For the residuals F_i (equality of fugacities)
1508  for (std::size_t i = 0; i < N; ++i) {
1509  // Equate the liquid and vapor fugacities
1510  CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1511  CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1512  r(i) = ln_f_liq - ln_f_vap;
1513 
1514  for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1515  if (bubble_point) {
1516  J(i, j) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
1517  } else {
1518  J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
1519  }
1520  }
1521  J(i, N - 1) =
1522  MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag);
1523  }
1524  } else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1525  // Independent variables are N-1 mole fractions of incipient phase and p
1526 
1527  // For the residuals F_i (equality of fugacities)
1528  for (std::size_t i = 0; i < N; ++i) {
1529  // Equate the liquid and vapor fugacities
1530  CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1531  CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1532  r(i) = ln_f_liq - ln_f_vap;
1533 
1534  for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1535  if (bubble_point) {
1536  J(i, j) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
1537  } else {
1538  J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
1539  }
1540  }
1541  J(i, N - 1) =
1542  MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatV, i, xN_flag);
1543  }
1544  } else {
1545  throw ValueError();
1546  }
1547 
1548  error_rms = r.norm();
1549 
1550  // Calculate derivatives along phase boundary;
1551  // Gernert thesis 3.96 and 3.97
1552  CoolPropDbl dQ_dPsat = 0, dQ_dTsat = 0;
1553  for (std::size_t i = 0; i < N; ++i) {
1554  dQ_dPsat += x[i]
1555  * (MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rSatL, i, xN_flag)
1556  - MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rSatV, i, xN_flag));
1557  dQ_dTsat += x[i]
1558  * (MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatL, i, xN_flag)
1559  - MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatV, i, xN_flag));
1560  }
1561  dTsat_dPsat = -dQ_dPsat / dQ_dTsat;
1562  dPsat_dTsat = -dQ_dTsat / dQ_dPsat;
1563 }
1564 
1565 void SaturationSolvers::newton_raphson_twophase::call(HelmholtzEOSMixtureBackend& HEOS, newton_raphson_twophase_options& IO) {
1566  int iter = 0;
1567 
1568  if (get_debug_level() > 9) {
1569  std::cout << " NRsat::call: p" << IO.p << " T" << IO.T << " dl" << IO.rhomolar_liq << " dv" << IO.rhomolar_vap << std::endl;
1570  }
1571 
1572  // Reset all the variables and resize
1573  pre_call();
1574 
1575  rhomolar_liq = IO.rhomolar_liq;
1576  rhomolar_vap = IO.rhomolar_vap;
1577  T = IO.T;
1578  p = IO.p;
1579  imposed_variable = IO.imposed_variable;
1580  x = IO.x;
1581  y = IO.y;
1582  z = IO.z;
1583  beta = IO.beta;
1584 
1585  this->N = z.size();
1586  x.resize(N);
1587  y.resize(N);
1588  r.resize(2 * N - 1);
1589  J.resize(2 * N - 1, 2 * N - 1);
1590  err_rel.resize(2 * N - 1);
1591 
1592  // Hold a pointer to the backend
1593  this->HEOS = &HEOS;
1594 
1595  do {
1596  // Build the Jacobian and residual vectors
1597  build_arrays();
1598 
1599  // Solve for the step; v is the step with the contents
1600  // [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]
1601 
1602  // Uncomment to see Jacobian and residual at every step
1603  // std::cout << vec_to_string(J, "%0.12Lg") << std::endl;
1604  // std::cout << vec_to_string(negative_r, "%0.12Lg") << std::endl;
1605 
1606  Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
1607 
1608  for (unsigned int i = 0; i < N - 1; ++i) {
1609  err_rel[i] = v[i] / x[i];
1610  x[i] += v[i];
1611  err_rel[i + (N - 1)] = v[i + (N - 1)] / y[i];
1612  y[i] += v[i + (N - 1)];
1613  }
1614  x[N - 1] = 1 - std::accumulate(x.begin(), x.end() - 1, 0.0);
1615  y[N - 1] = 1 - std::accumulate(y.begin(), y.end() - 1, 0.0);
1616 
1617  if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED) {
1618  T += v[2 * N - 2];
1619  err_rel[2 * N - 2] = v[2 * N - 2] / T;
1620  } else if (imposed_variable == newton_raphson_twophase_options::T_IMPOSED) {
1621  p += v[2 * N - 2];
1622  err_rel[2 * N - 2] = v[2 * N - 2] / p;
1623  } else {
1624  throw ValueError("invalid imposed_variable");
1625  }
1626  //std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " v " << vec_to_string(v, "%0.10Lg") << " x " << vec_to_string(x, "%0.10Lg") << " r " << vec_to_string(r, "%0.10Lg") << std::endl;
1627 
1628  min_rel_change = err_rel.cwiseAbs().minCoeff();
1629  iter++;
1630 
1631  if (iter == IO.Nstep_max) {
1632  throw ValueError(format("newton_raphson_saturation::call reached max number of iterations [%d]", IO.Nstep_max));
1633  }
1634  } while (this->error_rms > 1e-9 && min_rel_change > 1000 * DBL_EPSILON && iter < IO.Nstep_max);
1635 
1636  IO.Nsteps = iter;
1637  IO.p = p;
1638  IO.x = x; // Mole fractions in liquid
1639  IO.y = y; // Mole fractions in vapor
1640  IO.T = T;
1641  IO.rhomolar_liq = rhomolar_liq;
1642  IO.rhomolar_vap = rhomolar_vap;
1643  IO.hmolar_liq = HEOS.SatL.get()->hmolar();
1644  IO.hmolar_vap = HEOS.SatV.get()->hmolar();
1645  IO.smolar_liq = HEOS.SatL.get()->smolar();
1646  IO.smolar_vap = HEOS.SatV.get()->smolar();
1647 }
1648 
1649 void SaturationSolvers::newton_raphson_twophase::build_arrays() {
1650  // References to the classes for concision
1651  HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
1652 
1653  // Step 0:
1654  // -------
1655  // Set mole fractions
1656  rSatL.set_mole_fractions(x);
1657  rSatV.set_mole_fractions(y);
1658 
1659  //std::vector<CoolPropDbl> &x = rSatL.get_mole_fractions();
1660  //std::vector<CoolPropDbl> &y = rSatV.get_mole_fractions();
1661 
1662  rSatL.update_TP_guessrho(T, p, rhomolar_liq);
1663  rhomolar_liq = rSatL.rhomolar();
1664  rSatV.update_TP_guessrho(T, p, rhomolar_vap);
1665  rhomolar_vap = rSatV.rhomolar();
1666 
1667  // For diagnostic purposes calculate the pressures (no derivatives are evaluated)
1668  CoolPropDbl p_liq = rSatL.p();
1669  CoolPropDbl p_vap = rSatV.p();
1670  p = 0.5 * (p_liq + p_vap);
1671 
1672  // Step 2:
1673  // -------
1674  // Build the residual vector and the Jacobian matrix
1675 
1676  x_N_dependency_flag xN_flag = XN_DEPENDENT;
1677 
1678  // Form of residuals do not depend on which variable is imposed
1679  for (std::size_t i = 0; i < N; ++i) {
1680  // Equate the liquid and vapor fugacities
1681  CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1682  CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1683  r[i] = ln_f_liq - ln_f_vap; // N of these
1684 
1685  if (i != N - 1) {
1686  // Equate the specified vapor mole fraction and that given defined by the ith component
1687  r[i + N] = (z[i] - x[i]) / (y[i] - x[i]) - beta; // N-1 of these
1688  }
1689  }
1690 
1691  // First part of derivatives with respect to ln f_i
1692  for (std::size_t i = 0; i < N; ++i) {
1693  for (std::size_t j = 0; j < N - 1; ++j) {
1694  J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
1695  J(i, j + N - 1) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
1696  }
1697 
1698  // Last derivative with respect to either T or p depending on what is imposed
1699  if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED) {
1700  J(i, 2 * N - 2) =
1701  MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag);
1702  } else if (imposed_variable == newton_raphson_twophase_options::T_IMPOSED) {
1703  J(i, 2 * N - 2) =
1704  MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatV, i, xN_flag);
1705  } else {
1706  throw ValueError();
1707  }
1708  }
1709  // Derivatives with respect to the vapor mole fractions residual
1710  for (std::size_t i = 0; i < N - 1; ++i) {
1711  std::size_t k = i + N; // N ln f_i residuals
1712  J(k, i) = (z[i] - y[i]) / pow(y[i] - x[i], 2);
1713  J(k, i + (N - 1)) = -(z[i] - x[i]) / pow(y[i] - x[i], 2);
1714  }
1715 
1716  error_rms = r.norm(); // Square-root (The R in RMS)
1717 }
1718 
1719 class RachfordRiceResidual : public FuncWrapper1DWithDeriv
1720 {
1721  private:
1722  const std::vector<double>&z, &lnK;
1723 
1724  public:
1725  RachfordRiceResidual(const std::vector<double>& z, const std::vector<double>& lnK) : z(z), lnK(lnK){};
1726  double call(double beta) {
1727  return FlashRoutines::g_RachfordRice(z, lnK, beta);
1728  }
1729  double deriv(double beta) {
1730  return FlashRoutines::dgdbeta_RachfordRice(z, lnK, beta);
1731  }
1732 };
1733 
1734 void StabilityRoutines::StabilityEvaluationClass::trial_compositions() {
1735 
1736  x.resize(z.size());
1737  y.resize(z.size());
1738  lnK.resize(z.size());
1739  K.resize(z.size());
1740  double g0 = 0, g1 = 0, beta = -1;
1741 
1742  for (int i = 0; i < static_cast<int>(z.size()); ++i) {
1743  // Calculate the K-factor
1744  if (m_T < 0 && m_p < 0) {
1745  // Using T&P from the class
1746  lnK[i] = SaturationSolvers::Wilson_lnK_factor(HEOS, HEOS.T(), HEOS.p(), i);
1747  } else {
1748  // Using specified T&P
1749  lnK[i] = SaturationSolvers::Wilson_lnK_factor(HEOS, m_T, m_p, i);
1750  }
1751  K[i] = exp(lnK[i]);
1752  g0 += z[i] * (K[i] - 1); // The summation for beta = 0
1753  g1 += z[i] * (1 - 1 / K[i]); // The summation for beta = 1
1754  }
1755  // Copy K-factors for later use
1756  K0 = K;
1757  // Now see what to do about the g(0) and g(1) values
1758  // -----
1759  //
1760  if (g0 < 0) {
1761  beta = 0; // Assumed to be at bubble-point temperature
1762  } else if (g1 > 0) {
1763  beta = 1; // Assumed to be at the dew-point temperature
1764  } else {
1765  // Need to iterate to find beta that makes g of Rachford-Rice zero
1766  RachfordRiceResidual resid(z, lnK);
1767  beta = Brent(resid, 0, 1, DBL_EPSILON, 1e-10, 100);
1768  }
1769  // Get the compositions from given value for beta, K, z
1770  SaturationSolvers::x_and_y_from_K(beta, K, z, x, y);
1771  normalize_vector(x);
1772  normalize_vector(y);
1773  if (debug) {
1774  std::cout << format("1) T: %g p: %g beta: %g\n", HEOS.T(), HEOS.p(), beta);
1775  }
1776 }
1777 void StabilityRoutines::StabilityEvaluationClass::successive_substitution(int num_steps) {
1778  // ----
1779  // Do a few steps of successive substitution
1780  // ----
1781 
1782  HEOS.SatL->set_mole_fractions(x);
1783  HEOS.SatL->calc_reducing_state();
1784  HEOS.SatV->set_mole_fractions(y);
1785  HEOS.SatV->calc_reducing_state();
1786 
1787  if (debug) {
1788  std::cout << format("2) SS1: i beta K x y rho' rho''\n");
1789  }
1790  for (int step_count = 0; step_count < num_steps; ++step_count) {
1791  // Set the composition
1792  HEOS.SatL->set_mole_fractions(x);
1793  HEOS.SatV->set_mole_fractions(y);
1794  HEOS.SatL->calc_reducing_state();
1795  HEOS.SatV->calc_reducing_state();
1796 
1797  this->rho_TP_global();
1798 
1799  // Calculate the new K-factors from the fugacity coefficients
1800  double g0 = 0, g1 = 0;
1801  for (std::size_t i = 0; i < z.size(); ++i) {
1802  lnK[i] = log(HEOS.SatL->fugacity_coefficient(i) / HEOS.SatV->fugacity_coefficient(i));
1803  K[i] = exp(lnK[i]);
1804  g0 += z[i] * (K[i] - 1); // The summation for beta = 0
1805  g1 += z[i] * (1 - 1 / K[i]); // The summation for beta = 1
1806  }
1807  RachfordRiceResidual resid(z, lnK);
1808  if (g0 < 0) {
1809  beta = 0;
1810  } else if (g1 > 0) {
1811  beta = 1;
1812  } else {
1813  // Need to iterate to find beta that makes g of Rachford-Rice zero
1814  beta = Brent(resid, 0, 1, DBL_EPSILON, 1e-10, 100);
1815  }
1816 
1817  // Get the compositions from given values for beta, K, z
1818  SaturationSolvers::x_and_y_from_K(beta, K, z, x, y);
1819  normalize_vector(x);
1820  normalize_vector(y);
1821  if (debug) {
1822  std::cout << format("2) %d %g %s %s %s %g %g\n", step_count, beta, vec_to_string(K, "%0.6f").c_str(), vec_to_string(x, "%0.6f").c_str(),
1823  vec_to_string(y, "%0.6f").c_str(), rhomolar_liq, rhomolar_vap);
1824  }
1825  }
1826 }
1827 void StabilityRoutines::StabilityEvaluationClass::check_stability() {
1828  std::vector<double> tpdL, tpdH;
1829 
1830  // Calculate the temperature and pressure to be used
1831  double the_T = (m_T > 0 && m_p > 0) ? m_T : HEOS.T();
1832  double the_p = (m_T > 0 && m_p > 0) ? m_p : HEOS.p();
1833 
1834  // If beta value is between epsilon and 1-epsilon, check the TPD
1835  if (beta > DBL_EPSILON && beta < 1 - DBL_EPSILON) {
1836 
1837  // Set the composition back to the bulk composition for both liquid and vapor phases
1838  HEOS.SatL->set_mole_fractions(z);
1839  HEOS.SatV->set_mole_fractions(z);
1840  HEOS.SatL->calc_reducing_state();
1841  HEOS.SatV->calc_reducing_state();
1842 
1843  // Update the densities in each class
1844  double rhoL = HEOS.SatL->solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SatL->SRK_covolume());
1845  double rhoV = HEOS.SatV->solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SatV->SRK_covolume());
1846  HEOS.SatL->update_DmolarT_direct(rhoL, the_T);
1847  HEOS.SatV->update_DmolarT_direct(rhoV, the_T);
1848 
1849  // Calculate the tpd and the Gibbs energy difference (Gernert, 2014, Eqs. 20-22)
1850  // The trial compositions are the phase compositions from before
1851  this->tpd_liq = HEOS.SatL->tangent_plane_distance(the_T, the_p, x, rhomolar_liq);
1852  this->tpd_vap = HEOS.SatV->tangent_plane_distance(the_T, the_p, y, rhomolar_vap);
1853 
1854  this->DELTAG_nRT = (1 - beta) * tpd_liq + beta * (tpd_vap);
1855  if (debug) {
1856  std::cout << format("3) tpd': %g tpd'': %g DELTAG/nRT: %g\n", tpd_liq, tpd_vap, DELTAG_nRT);
1857  }
1858 
1859  // If any of these cases are met, feed is conclusively unstable, stop!
1860  if (this->tpd_liq < -DBL_EPSILON || this->tpd_vap < -DBL_EPSILON || this->DELTAG_nRT < -DBL_EPSILON) {
1861  if (debug) {
1862  std::cout << format("3) PHASE SPLIT beta in (eps,1-eps) \n");
1863  }
1864  _stable = false;
1865  return;
1866  }
1867  }
1868 
1869  // Ok, we aren't sure about stability, need to keep going with the full tpd analysis
1870 
1871  // Use the global density solver to obtain the density root (or the lowest Gibbs energy root if more than one)
1872  CoolPropDbl rho_bulk = HEOS.solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SRK_covolume());
1873  HEOS.update_DmolarT_direct(rho_bulk, the_T);
1874 
1875  // Calculate the fugacity coefficient at initial composition of the bulk phase
1876  std::vector<double> fugacity_coefficient0(z.size()), fugacity0(z.size());
1877  for (std::size_t i = 0; i < z.size(); ++i) {
1878  fugacity_coefficient0[i] = HEOS.fugacity_coefficient(i);
1879  fugacity0[i] = HEOS.fugacity(i);
1880  }
1881 
1882  // Generate light and heavy test compositions (Gernert, 2014, Eq. 23)
1883  xL.resize(z.size());
1884  xH.resize(z.size());
1885  for (std::size_t i = 0; i < z.size(); ++i) {
1886  xL[i] = z[i] * K0[i]; // Light-phase composition
1887  xH[i] = z[i] / K0[i]; // Heavy-phase composition
1888  }
1889  normalize_vector(xL);
1890  normalize_vector(xH);
1891 
1892  // For each composition, use successive substitution to try to evaluate stability
1893  if (debug) {
1894  std::cout << format("3) SS2: i x' x'' rho' rho'' tpd' tpd''\n");
1895  }
1896 
1897  // We got this far, we assume stable phases
1898  _stable = true;
1899 
1900  double diffbulkL = 0, diffbulkH = 0;
1901  for (int step_count = 0; step_count < 100; ++step_count) {
1902 
1903  // Set the composition
1904  HEOS.SatL->set_mole_fractions(xH);
1905  HEOS.SatV->set_mole_fractions(xL);
1906  HEOS.SatL->calc_reducing_state();
1907  HEOS.SatV->calc_reducing_state();
1908 
1909  // Do the global density solver for both phases
1910  rho_TP_global();
1911 
1912  double tpd_L = 0, tpd_H = 0;
1913  for (std::size_t i = 0; i < xL.size(); ++i) {
1914  tpd_L += xL[i] * (log(MixtureDerivatives::fugacity_i(*HEOS.SatV, i, XN_DEPENDENT)) - log(fugacity0[i]));
1915  tpd_H += xH[i] * (log(MixtureDerivatives::fugacity_i(*HEOS.SatL, i, XN_DEPENDENT)) - log(fugacity0[i]));
1916  }
1917  tpdL.push_back(tpd_L);
1918  tpdH.push_back(tpd_H);
1919 
1920  // Calculate the new composition from the fugacity coefficients
1921  diffbulkL = 0, diffbulkH = 0;
1922  for (std::size_t i = 0; i < z.size(); ++i) {
1923  xL[i] = z[i] * fugacity_coefficient0[i] / HEOS.SatV->fugacity_coefficient(i);
1924  diffbulkL += std::abs(xL[i] - z[i]);
1925  xH[i] = z[i] * fugacity_coefficient0[i] / HEOS.SatL->fugacity_coefficient(i);
1926  diffbulkH += std::abs(xH[i] - z[i]);
1927  }
1928  normalize_vector(xL);
1929  normalize_vector(xH);
1930  if (debug) {
1931  std::cout << format("3) %d %s %s %g %g %g %g\n", step_count, vec_to_string(xL, "%0.6f").c_str(), vec_to_string(xH, "%0.6f").c_str(),
1932  rhomolar_liq, rhomolar_vap, tpd_L, tpd_H);
1933  }
1934 
1935  // Check if either of the phases have the bulk composition. If so, no phase split
1936  if (diffbulkL < 1e-2 || diffbulkH < 1e-2) {
1937  _stable = true;
1938  return;
1939  }
1940 
1941  // Check if either tpd is negative, if so, phases definitively split, quit
1942  if (tpd_L < -1e-12 || tpd_H < -1e-12) {
1943  _stable = false;
1944  return;
1945  }
1946  }
1947  if (diffbulkH > 0.25 || diffbulkL > 0.25) {
1948  // At least one test phase is definitely not the bulk composition, so phase split predicted
1949  _stable = false;
1950  }
1951 }
1952 
1953 void StabilityRoutines::StabilityEvaluationClass::rho_TP_global() {
1954 
1955  // Calculate the temperature and pressure to be used
1956  double the_T = (m_T > 0 && m_p > 0) ? m_T : HEOS.T();
1957  double the_p = (m_T > 0 && m_p > 0) ? m_p : HEOS.p();
1958 
1959  // Calculate covolume of SRK, use it as the maximum density
1960  double rhoL = HEOS.SatL->solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SatL->SRK_covolume());
1961  double rhoV = HEOS.SatV->solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SatV->SRK_covolume());
1962  HEOS.SatL->update_DmolarT_direct(rhoL, the_T);
1963  HEOS.SatV->update_DmolarT_direct(rhoV, the_T);
1964 
1965  rhomolar_liq = HEOS.SatL->rhomolar();
1966  rhomolar_vap = HEOS.SatV->rhomolar();
1967 }
1968 
1969 void StabilityRoutines::StabilityEvaluationClass::rho_TP_w_guesses() {
1970 
1971  // Re-calculate the density
1972  if (m_T > 0 && m_p > 0) {
1973  HEOS.SatL->update_TP_guessrho(m_T, m_p, rhomolar_liq);
1974  HEOS.SatV->update_TP_guessrho(m_T, m_p, rhomolar_vap);
1975  } else {
1976  HEOS.SatL->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_liq);
1977  HEOS.SatV->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_vap);
1978  }
1979  rhomolar_liq = HEOS.SatL->rhomolar();
1980  rhomolar_vap = HEOS.SatV->rhomolar();
1981 }
1982 
1983 void StabilityRoutines::StabilityEvaluationClass::rho_TP_SRK_translated() {
1984 
1985  // First use cubic as a guess for the density of liquid and vapor phases
1986  if (m_T > 0 && m_p > 0) {
1987  rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(m_T, m_p, iphase_liquid); // [mol/m^3]
1988  rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(m_T, m_p, iphase_gas); // [mol/m^3]
1989  } else {
1990  rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(HEOS.T(), HEOS.p(), iphase_liquid); // [mol/m^3]
1991  rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(HEOS.T(), HEOS.p(), iphase_gas); // [mol/m^3]
1992  }
1993 
1994  // Apply volume translation to liquid density only
1995  if (HEOS.backend_name().find("Helmholtz") == 0) {
1996  // Use Peneloux volume translation to shift liquid volume
1997  // As in Horstmann :: doi:10.1016/j.fluid.2004.11.002
1998  double summer_c = 0, v_SRK = 1 / rhomolar_liq;
1999  for (std::size_t i = 0; i < z.size(); ++i) {
2000  // Get the parameters for the cubic EOS
2001  CoolPropDbl Tc = HEOS.get_fluid_constant(i, iT_critical), pc = HEOS.get_fluid_constant(i, iP_critical),
2002  rhomolarc = HEOS.get_fluid_constant(i, irhomolar_critical);
2003  CoolPropDbl R = 8.3144598;
2004  summer_c += z[i] * (0.40768 * R * Tc / pc * (0.29441 - pc / (rhomolarc * R * Tc)));
2005  }
2006  rhomolar_liq = 1 / (v_SRK - summer_c);
2007  }
2008 }
2009 
2010 void SaturationSolvers::PTflash_twophase::solve() {
2011  const std::size_t N = IO.x.size();
2012  int iter = 0;
2013  double min_rel_change;
2014  do {
2015  // Build the Jacobian and residual vectors
2016  build_arrays();
2017 
2018  // Solve for the step; v is the step with the contents
2019  // [delta(x'_0), delta(x'_1), ..., delta(x'_{N-1}), delta(x''_0), delta(x''_1), ..., delta(x''_{N-1})]
2020 
2021  // Uncomment to see Jacobian and residual at every step
2022  //std::cout << vec_to_string(J, "%0.12Lg") << std::endl;
2023  //std::cout << vec_to_string(-r, "%0.12Lg") << std::endl;
2024 
2025  Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
2026 
2027  for (unsigned int i = 0; i < N - 1; ++i) {
2028  err_rel[i] = v[i] / IO.x[i];
2029  IO.x[i] += v[i];
2030  err_rel[i + (N - 1)] = v[i + (N - 1)] / IO.y[i];
2031  IO.y[i] += v[i + (N - 1)];
2032  }
2033  IO.x[N - 1] = 1 - std::accumulate(IO.x.begin(), IO.x.end() - 1, 0.0);
2034  IO.y[N - 1] = 1 - std::accumulate(IO.y.begin(), IO.y.end() - 1, 0.0);
2035 
2036  //std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " v " << vec_to_string(v, "%0.10Lg") << " x " << vec_to_string(x, "%0.10Lg") << " r " << vec_to_string(r, "%0.10Lg") << std::endl;
2037 
2038  min_rel_change = err_rel.cwiseAbs().minCoeff();
2039  iter++;
2040 
2041  if (iter == IO.Nstep_max) {
2042  throw ValueError(format("PTflash_twophase::call reached max number of iterations [%d]", IO.Nstep_max));
2043  }
2044  } while (this->error_rms > 1e-9 && min_rel_change > 1000 * DBL_EPSILON && iter < IO.Nstep_max);
2045 }
2046 void SaturationSolvers::PTflash_twophase::build_arrays() {
2047  const std::size_t N = IO.x.size();
2048 
2049  r.resize(2 * N - 2);
2050  J.resize(2 * N - 2, 2 * N - 2);
2051  err_rel.resize(2 * N - 2);
2052 
2053  HEOS.SatL->set_mole_fractions(IO.x);
2054  HEOS.SatL->update_TP_guessrho(IO.T, IO.p, IO.rhomolar_liq);
2055 
2056  HEOS.SatV->set_mole_fractions(IO.y);
2057  HEOS.SatV->update_TP_guessrho(IO.T, IO.p, IO.rhomolar_vap);
2058 
2059  // Independent variables are
2060  // [delta(x'_0), delta(x'_1), ..., delta(x'_{N-1}), delta(x''_0), delta(x''_1), ..., delta(x''_{N-1})]
2061 
2062  // First N residuals are the iso-fugacity condition
2063  for (std::size_t k = 0; k < N; ++k) {
2064  r(k) = log(HEOS.SatL->fugacity(k) / HEOS.SatV->fugacity(k));
2065  for (std::size_t j = 0; j < N - 1; ++j) {
2066  if (k == N - 1) {
2067  J(k, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatL.get()), k, j, XN_DEPENDENT);
2068  J(k, j + N - 1) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatV.get()), k, j, XN_DEPENDENT);
2069  } else {
2070  J(k, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatL.get()), k, j, XN_DEPENDENT);
2071  J(k, j + N - 1) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatV.get()), k, j, XN_DEPENDENT);
2072  }
2073  }
2074  }
2075  // Next N-2 residuals are amount of substance balances
2076  for (std::size_t i = 0; i < N - 2; ++i) {
2077  std::size_t k = i + N;
2078  r(k) = (IO.z[i] - IO.x[i]) / (IO.y[i] - IO.x[i]) - (IO.z[N - 1] - IO.x[N - 1]) / (IO.y[N - 1] - IO.x[N - 1]);
2079  for (std::size_t j = 0; j < N - 2; ++j) {
2080  J(k, j) = (IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
2081  J(k, j + N - 1) = -(IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
2082  }
2083  std::size_t j = N - 2;
2084  J(k, j) = -(IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
2085  J(k, j + N - 1) = +(IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
2086  }
2087  this->error_rms = r.norm();
2088 }
2089 } /* namespace CoolProp*/
2090 
2091 #if defined(ENABLE_CATCH)
2092 # include <catch2/catch_all.hpp>
2093 
2094 TEST_CASE("Check the PT flash calculation for two-phase inputs", "[PTflash_twophase]") {
2095  shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", "Propane&Ethane"));
2096  AS->set_mole_fractions(std::vector<double>(2, 0.5));
2097  // Dewpoint calculation
2098  AS->update(CoolProp::PQ_INPUTS, 101325, 1);
2099 
2100  // Do a dummy calculation at the dewpoint state - make sure all the residuals are zero as they should be
2101  CoolProp::SaturationSolvers::PTflash_twophase_options o;
2102  o.x = AS->mole_fractions_liquid();
2103  o.y = AS->mole_fractions_vapor();
2104  o.z = AS->get_mole_fractions();
2105  o.rhomolar_liq = AS->saturated_liquid_keyed_output(CoolProp::iDmolar);
2106  o.rhomolar_vap = AS->saturated_vapor_keyed_output(CoolProp::iDmolar);
2107  o.T = AS->T();
2108  o.p = AS->p();
2109  o.omega = 1.0;
2110  CoolProp::SaturationSolvers::PTflash_twophase solver(*static_cast<CoolProp::HelmholtzEOSMixtureBackend*>(AS.get()), o);
2111  solver.build_arrays();
2112  double err = solver.r.norm();
2113  REQUIRE(std::abs(err) < 1e-10);
2114 
2115  // Now, perturb the solution a little bit and actually do the solve
2116  std::vector<double> x0 = o.x;
2117  o.x[0] *= 1.1;
2118  o.x[1] = 1 - o.x[0];
2119  solver.solve();
2120  // Make sure we end up with the same liquid composition
2121  double diffx0 = o.x[0] - x0[0];
2122  REQUIRE(std::abs(diffx0) < 1e-10);
2123 
2124  // Now do the blind flash call with PT as inputs
2125  AS->update(CoolProp::PT_INPUTS, AS->p(), AS->T() - 2);
2126  REQUIRE(AS->phase() == CoolProp::iphase_twophase);
2127 }
2128 
2129 #endif
Generated on Fri Oct 4 2024 09:06:33 for CoolProp by doxygen 1.9.1