CoolProp  6.6.0
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 }
558 void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl rhomolar, saturation_D_pure_options& options)
559 {
560  /*
561  This function is inspired by the method of Akasaka:
562 
563  R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
564  Helmholtz Energy Equations of State",
565  Journal of Thermal Science and Technology v3 n3,2008
566 
567  Ancillary equations are used to get a sensible starting point
568  */
569  std::vector<CoolPropDbl> r(2, _HUGE), v;
570  std::vector<std::vector<CoolPropDbl>> J(2, std::vector<CoolPropDbl>(2, _HUGE));
571 
572  HEOS.calc_reducing_state();
573  const SimpleState& reduce = HEOS.get_reducing_state();
574  shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
575 
576  CoolPropDbl T, rhoL, rhoV;
577  CoolPropDbl deltaL = 0, deltaV = 0, tau = 0, error, p_error;
578  int iter = 0;
579 
580  // Use the density ancillary function as the starting point for the solver
581  try {
582  if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL) {
583  // Invert liquid density ancillary to get temperature
584  // TODO: fit inverse ancillaries too
585  T = HEOS.get_components()[0].ancillaries.rhoL.invert(rhomolar);
586  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
587  rhoL = rhomolar;
588  } else if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOV) {
589  // Invert vapor density ancillary to get temperature
590  // TODO: fit inverse ancillaries too
591  T = HEOS.get_components()[0].ancillaries.rhoV.invert(rhomolar);
592  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
593  rhoV = rhomolar;
594  } else {
595  throw ValueError(format("imposed rho to saturation_D_pure [%d%] is invalid", options.imposed_rho));
596  }
597 
598  deltaL = rhoL / reduce.rhomolar;
599  deltaV = rhoV / reduce.rhomolar;
600  tau = reduce.T / T;
601  } catch (NotImplementedError&) {
602  throw; // ??? What is this try...catch for?
603  }
604 
605  do {
606  /*if (get_debug_level()>8){
607  std::cout << format("%s:%d: right before the derivs with deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
608  }*/
609 
610  // Calculate once to save on calls to EOS
611  SatL->update(DmolarT_INPUTS, rhoL, T);
612  SatV->update(DmolarT_INPUTS, rhoV, T);
613 
614  CoolPropDbl pL = SatL->p();
615  CoolPropDbl pV = SatV->p();
616 
617  // These derivatives are needed for both cases
618  CoolPropDbl dalphar_dtauL = SatL->dalphar_dTau();
619  CoolPropDbl dalphar_dtauV = SatV->dalphar_dTau();
620  CoolPropDbl d2alphar_ddelta_dtauL = SatL->d2alphar_dDelta_dTau();
621  CoolPropDbl d2alphar_ddelta_dtauV = SatV->d2alphar_dDelta_dTau();
622  CoolPropDbl alpharL = SatL->alphar();
623  CoolPropDbl alpharV = SatV->alphar();
624  CoolPropDbl dalphar_ddeltaL = SatL->dalphar_dDelta();
625  CoolPropDbl dalphar_ddeltaV = SatV->dalphar_dDelta();
626 
627  // -r_1
628  r[0] = -(deltaV * (1 + deltaV * dalphar_ddeltaV) - deltaL * (1 + deltaL * dalphar_ddeltaL));
629  // -r_2
630  r[1] = -(deltaV * dalphar_ddeltaV + alpharV + log(deltaV) - deltaL * dalphar_ddeltaL - alpharL - log(deltaL));
631 
632  // dr1_dtau
633  J[0][0] = pow(deltaV, 2) * d2alphar_ddelta_dtauV - pow(deltaL, 2) * d2alphar_ddelta_dtauL;
634  // dr2_dtau
635  J[1][0] = deltaV * d2alphar_ddelta_dtauV + dalphar_dtauV - deltaL * d2alphar_ddelta_dtauL - dalphar_dtauL;
636 
637  if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL) {
638  CoolPropDbl d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
639  if (options.use_logdelta) {
640  J[0][1] = deltaV + 2 * pow(deltaV, 2) * dalphar_ddeltaV + pow(deltaV, 3) * d2alphar_ddelta2V;
641  J[1][1] = pow(deltaV, 2) * d2alphar_ddelta2V + 2 * deltaV * dalphar_ddeltaV + 1;
642  } else {
643  J[0][1] = 1 + 2 * deltaV * dalphar_ddeltaV + pow(deltaV, 2) * d2alphar_ddelta2V;
644  J[1][1] = deltaV * d2alphar_ddelta2V + 2 * dalphar_ddeltaV + 1 / deltaV;
645  }
646  } else if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOV) {
647  CoolPropDbl d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
648  if (options.use_logdelta) {
649  J[0][1] = -deltaL - 2 * pow(deltaL, 2) * dalphar_ddeltaL - pow(deltaL, 3) * d2alphar_ddelta2L;
650  J[1][1] = -pow(deltaL, 2) * d2alphar_ddelta2L - 2 * deltaL * dalphar_ddeltaL - 1;
651  } else {
652  J[0][1] = -1 - 2 * deltaL * dalphar_ddeltaL - pow(deltaL, 2) * d2alphar_ddelta2L;
653  J[1][1] = -deltaL * d2alphar_ddelta2L - 2 * dalphar_ddeltaL - 1 / deltaL;
654  }
655  }
656 
657  //double DET = J[0][0]*J[1][1]-J[0][1]*J[1][0];
658 
659  v = linsolve(J, r);
660 
661  double omega = options.omega;
662 
663  if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL) {
664  if (options.use_logdelta)
665  deltaV = exp(log(deltaV)+omega*v[1]);
666  else
667  {
668  if (deltaV + omega*v[1] <= 0) {omega = 0.5*deltaV/v[1];} // gone off track, take a smaller step
669  if (tau + omega*v[0] <= 0) {omega = 0.5*tau/v[0];}
670  deltaV += omega*v[1];
671  }
672  }
673  else
674  {
675  if (options.use_logdelta)
676  deltaL = exp(log(deltaL)+omega*v[1]);
677  else
678  {
679  if (deltaL + omega*v[1] <= 0) {omega = 0.5*deltaL/v[1];} // gone off track, take a smaller step
680  if (tau + omega*v[0] <= 0) {omega = 0.5*tau/v[0];}
681  deltaL += omega*v[1];
682  }
683  }
684 
685  tau += omega*v[0];
686 
687  rhoL = deltaL*reduce.rhomolar;
688  rhoV = deltaV*reduce.rhomolar;
689  T = reduce.T/tau;
690 
691  p_error = (pL-pV)/pL;
692 
693  error = sqrt(pow(r[0], 2) + pow(r[1], 2));
694  iter++;
695  if (T < 0) {
696  throw SolutionError(format("saturation_D_pure solver T < 0"));
697  }
698  if (iter > options.max_iterations){
699  throw SolutionError(format("saturation_D_pure solver did not converge after %d iterations with rho: %g mol/m^3",options.max_iterations,rhomolar));
700  }
701  } while (error > 1e-9);
702  CoolPropDbl p_error_limit = 1e-3;
703  if (std::abs(p_error) > p_error_limit) {
704  throw SolutionError(format("saturation_D_pure solver abs error on p [%Lg] > limit [%Lg]", p_error, p_error_limit));
705  }
706 }
707 void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options) {
708  // Set some input options
709  SaturationSolvers::saturation_T_pure_Akasaka_options _options(false);
710  _options.omega = 1.0;
711  try {
712  // Actually call the solver
713  SaturationSolvers::saturation_T_pure_Maxwell(HEOS, T, _options);
714  } catch (...) {
715  try {
716  // Actually call the solver
717  SaturationSolvers::saturation_T_pure_Akasaka(HEOS, T, _options);
718  } catch (...) {
719  // If there was an error, store values for use in later solvers
720  options.pL = _options.pL;
721  options.pV = _options.pV;
722  options.rhoL = _options.rhoL;
723  options.rhoV = _options.rhoV;
724  options.p = _options.pL;
725  SaturationSolvers::saturation_T_pure_1D_P(HEOS, T, options);
726  }
727  }
728 }
729 void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options) {
730  // Start with the method of Akasaka
731 
732  /*
733  This function implements the method of Akasaka
734 
735  R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
736  Helmholtz Energy Equations of State",
737  Journal of Thermal Science and Technology v3 n3,2008
738 
739  Ancillary equations are used to get a sensible starting point
740  */
741 
742  HEOS.calc_reducing_state();
743  const SimpleState& reduce = HEOS.get_reducing_state();
744  CoolPropDbl R_u = HEOS.gas_constant();
745  shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
746 
747  CoolPropDbl rhoL = _HUGE, rhoV = _HUGE, JL, JV, KL, KV, dJL, dJV, dKL, dKV;
748  CoolPropDbl DELTA, deltaL = 0, deltaV = 0, error, PL, PV, stepL, stepV;
749  int iter = 0;
750 
751  try {
752  if (options.use_guesses) {
753  // Use the guesses provided in the options structure
754  rhoL = options.rhoL;
755  rhoV = options.rhoV;
756  } else {
757  // Use the density ancillary function as the starting point for the solver
758 
759  // If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
760  if (T > 0.99 * HEOS.get_reducing_state().T) {
761  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T - 0.1);
762  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T - 0.1);
763  } else {
764  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
765  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
766 
767  // Apply a single step of Newton's method to improve guess value for liquid
768  // based on the error between the gas pressure (which is usually very close already)
769  // and the liquid pressure, which can sometimes (especially at low pressure),
770  // be way off, and often times negative
771  SatL->update(DmolarT_INPUTS, rhoL, T);
772  SatV->update(DmolarT_INPUTS, rhoV, T);
773 
774  // Update the guess for liquid density using density solver with vapor pressure
775  // and liquid density guess from ancillaries
776  HEOS.specify_phase(iphase_liquid);
777  rhoL = HEOS.solver_rho_Tp(T, SatV->p(), rhoL);
778  HEOS.unspecify_phase();
779  }
780  }
781 
782  deltaL = rhoL / reduce.rhomolar;
783  deltaV = rhoV / reduce.rhomolar;
784  } catch (NotImplementedError&) {
785  /*double Tc = crit.T;
786  double pc = crit.p.Pa;
787  double w = 6.67228479e-09*Tc*Tc*Tc-7.20464352e-06*Tc*Tc+3.16947758e-03*Tc-2.88760012e-01;
788  double q = -6.08930221451*w -5.42477887222;
789  double pt = exp(q*(Tc/T-1))*pc;*/
790 
791  //double rhoL = density_Tp_Soave(T, pt, 0), rhoV = density_Tp_Soave(T, pt, 1);
792 
793  //deltaL = rhoL/reduce.rhomolar;
794  //deltaV = rhoV/reduce.rhomolar;
795  //tau = reduce.T/T;
796  }
797  //if (get_debug_level()>5){
798  // std::cout << format("%s:%d: Akasaka guess values deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
799  // }
800 
801  do {
802  /*if (get_debug_level()>8){
803  std::cout << format("%s:%d: right before the derivs with deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
804  }*/
805 
806  // Calculate once to save on calls to EOS
807  SatL->update(DmolarT_INPUTS, rhoL, T);
808  SatV->update(DmolarT_INPUTS, rhoV, T);
809  CoolPropDbl alpharL = SatL->alphar();
810  CoolPropDbl alpharV = SatV->alphar();
811  CoolPropDbl dalphar_ddeltaL = SatL->dalphar_dDelta();
812  CoolPropDbl dalphar_ddeltaV = SatV->dalphar_dDelta();
813  CoolPropDbl d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
814  CoolPropDbl d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
815 
816  JL = deltaL * (1 + deltaL * dalphar_ddeltaL);
817  JV = deltaV * (1 + deltaV * dalphar_ddeltaV);
818  KL = deltaL * dalphar_ddeltaL + alpharL + log(deltaL);
819  KV = deltaV * dalphar_ddeltaV + alpharV + log(deltaV);
820 
821  PL = R_u * reduce.rhomolar * T * JL;
822  PV = R_u * reduce.rhomolar * T * JV;
823 
824  // At low pressure, the magnitude of d2alphar_ddelta2L and d2alphar_ddelta2V are enormous, truncation problems arise for all the partials
825  dJL = 1 + 2 * deltaL * dalphar_ddeltaL + deltaL * deltaL * d2alphar_ddelta2L;
826  dJV = 1 + 2 * deltaV * dalphar_ddeltaV + deltaV * deltaV * d2alphar_ddelta2V;
827  dKL = 2 * dalphar_ddeltaL + deltaL * d2alphar_ddelta2L + 1 / deltaL;
828  dKV = 2 * dalphar_ddeltaV + deltaV * d2alphar_ddelta2V + 1 / deltaV;
829 
830  DELTA = dJV * dKL - dJL * dKV;
831 
832  error = sqrt((KL - KV) * (KL - KV) + (JL - JV) * (JL - JV));
833 
834  // Get the predicted step
835  stepL = options.omega / DELTA * ((KV - KL) * dJV - (JV - JL) * dKV);
836  stepV = options.omega / DELTA * ((KV - KL) * dJL - (JV - JL) * dKL);
837 
838  CoolPropDbl deltaL0 = deltaL, deltaV0 = deltaV;
839  // Conditions for an acceptable step are:
840  // a) rhoL > rhoV or deltaL > deltaV
841  for (double omega_local = 1.0; omega_local > 0.1; omega_local /= 1.1) {
842  deltaL = deltaL0 + omega_local * stepL;
843  deltaV = deltaV0 + omega_local * stepV;
844 
845  if (deltaL > 1 && deltaV < 1 && deltaV > 0) {
846  break;
847  }
848  }
849 
850  rhoL = deltaL * reduce.rhomolar;
851  rhoV = deltaV * reduce.rhomolar;
852  iter++;
853  if (iter > 100) {
854  throw SolutionError(format("Akasaka solver did not converge after 100 iterations"));
855  }
856  } while (error > 1e-10 && std::abs(stepL) > 10 * DBL_EPSILON * std::abs(stepL) && std::abs(stepV) > 10 * DBL_EPSILON * std::abs(stepV));
857 
858  CoolPropDbl p_error_limit = 1e-3;
859  CoolPropDbl p_error = (PL - PV) / PL;
860  if (std::abs(p_error) > p_error_limit) {
861  options.pL = PL;
862  options.pV = PV;
863  options.rhoL = rhoL;
864  options.rhoV = rhoV;
865  throw SolutionError(format("saturation_T_pure_Akasaka solver abs error on p [%g] > limit [%g]", std::abs(p_error), p_error_limit));
866  }
867 }
868 
869 CoolPropDbl sign(CoolPropDbl x) {
870  if (x > 0) {
871  return 1;
872  } else {
873  return -1;
874  }
875 }
876 
877 void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options) {
878 
879  /*
880  This function implements the method of
881 
882  Ancillary equations are used to get a sensible starting point
883  */
884 
885  HEOS.calc_reducing_state();
886  shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
887  CoolProp::SimpleState& crit = HEOS.get_components()[0].crit;
888  CoolPropDbl rhoL = _HUGE, rhoV = _HUGE, error = 999, DeltavL, DeltavV, pL, pV, p, last_error;
889  int iter = 0, small_step_count = 0,
890  backwards_step_count = 0; // Counter for the number of times you have taken a step that increases error
891 
892  try {
893  if (options.use_guesses) {
894  // Use the guesses provided in the options structure
895  rhoL = options.rhoL;
896  rhoV = options.rhoV;
897  } else {
898  // Use the density ancillary function as the starting point for the solver
899 
900  // If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
901  if (T > 0.9999 * HEOS.get_reducing_state().T) {
902  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T - 0.1);
903  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T - 0.1);
904  } else {
905  rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
906  rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
907  p = HEOS.get_components()[0].ancillaries.pV.evaluate(T);
908 
909  CoolProp::SimpleState& tripleL = HEOS.get_components()[0].triple_liquid;
910  CoolProp::SimpleState& tripleV = HEOS.get_components()[0].triple_vapor;
911 
912  // If the guesses are terrible, apply a simple correction
913  // but only if the limits are being checked
914  if ((rhoL < crit.rhomolar * 0.8 || rhoL > tripleL.rhomolar * 1.2 || rhoV > crit.rhomolar * 1.2 || rhoV < tripleV.rhomolar * 0.8)
915  && !get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
916  // Lets assume that liquid density is more or less linear with T
917  rhoL = (crit.rhomolar - tripleL.rhomolar) / (crit.T - tripleL.T) * (T - tripleL.T) + tripleL.rhomolar;
918  // Then we calculate pressure from this density
919  SatL->update_DmolarT_direct(rhoL, T);
920  // Then we assume vapor to be ideal gas
921  if (SatL->p() > 0) {
922  rhoV = SatL->p() / (SatL->gas_constant() * T);
923  } else {
924  rhoV = p / (SatL->gas_constant() * T);
925  }
926  // Update the vapor state
927  SatV->update_DmolarT_direct(rhoV, T);
928  } else {
929  SatL->update_DmolarT_direct(rhoL, T);
930  SatV->update_DmolarT_direct(rhoV, T);
931  }
932  if (get_debug_level() > 5) {
933  std::cout << format("[Maxwell] ancillaries T: %0.16Lg rhoL: %0.16Lg rhoV: %0.16Lg pL: %g pV: %g\n", T, rhoL, rhoV, SatL->p(),
934  SatV->p());
935  }
936 
937  // Update the guess for liquid density using density solver with vapor pressure
938  // and liquid density guess from ancillaries, but only if the pressures are not
939  // close to each other
940  if (std::abs((SatL->p() - p) / p) > 0.1) {
941  for (int iii = 0; iii < 6; ++iii) {
942  // Use Halley's method to update the liquid density (http://en.wikipedia.org/wiki/Halley%27s_method)
943  CoolPropDbl f = SatL->p() - SatV->p();
944  CoolPropDbl dpdrho = SatL->first_partial_deriv(iP, iDmolar, iT);
945  CoolPropDbl d2pdrho2 = SatL->second_partial_deriv(iP, iDmolar, iT, iDmolar, iT);
946  CoolPropDbl deltarhoLHalley = -(2 * f * dpdrho) / (2 * POW2(dpdrho) - f * d2pdrho2);
947  rhoL += deltarhoLHalley;
948  if (std::abs(deltarhoLHalley / rhoL) < DBL_EPSILON) {
949  break;
950  }
951  SatL->update_DmolarT_direct(rhoL, T);
952  }
953  }
954 
955  SatL->update_DmolarT_direct(rhoL, T);
956  SatV->update_DmolarT_direct(rhoV, T);
957 
958  // Update the guess for vapor density using density solver with vapor pressure
959  // and density guess from ancillaries, but only if the pressures are not
960  // close to each other
961  if (std::abs((SatV->p() - p) / p) > 0.1) {
962  HEOS.specify_phase(iphase_gas);
963  rhoV = SatV->solver_rho_Tp(T, p, rhoV);
964  HEOS.unspecify_phase();
965  }
966  }
967  }
968  } catch (NotImplementedError&) {
969  }
970 
971  if (rhoL < crit.rhomolar) {
972  rhoL = 1.01 * crit.rhomolar;
973  }
974  if (rhoV > crit.rhomolar) {
975  rhoV = 0.99 * crit.rhomolar;
976  }
977  last_error = _HUGE;
978  SatL->update_DmolarT_direct(rhoL, T);
979  SatV->update_DmolarT_direct(rhoV, T);
980  if (get_debug_level() > 5) {
981  std::cout << format("[Maxwell] starting T: %0.16Lg rhoL: %Lg rhoV: %Lg pL: %Lg pV: %g\n", T, rhoL, rhoV, SatL->p(), SatV->p());
982  }
983  do {
984  pL = SatL->p();
985  pV = SatV->p();
986  CoolPropDbl vL = 1 / SatL->rhomolar(), vV = 1 / SatV->rhomolar();
987  // Get alpha, the pressure derivative with volume at constant T
988  // Given by (dp/drho|T)*drhodv
989  CoolPropDbl alphaL = SatL->first_partial_deriv(iP, iDmolar, iT) * (-POW2(SatL->rhomolar()));
990  CoolPropDbl alphaV = SatV->first_partial_deriv(iP, iDmolar, iT) * (-POW2(SatV->rhomolar()));
991 
992  // Total helmholtz energy for liquid and vapor
993  CoolPropDbl RT = SatL->gas_constant() * T;
994  CoolPropDbl helmholtzL = (SatL->calc_alpha0() + SatL->calc_alphar()) * RT;
995  CoolPropDbl helmholtzV = (SatV->calc_alpha0() + SatV->calc_alphar()) * RT;
996 
997  // Calculate the mean pressure
998  CoolPropDbl pM = (helmholtzL - helmholtzV) / (vV - vL);
999 
1000  // Coefficients for the quadratic in the step
1001  CoolPropDbl A = 0.5 * alphaL * (alphaL - alphaV);
1002  CoolPropDbl B = alphaL * (pL - pV - alphaV * (vL - vV));
1003  CoolPropDbl C = alphaV * (vL - vV) * (pM - pL) + 0.5 * POW2(pL - pV);
1004 
1005  // Argument to square root
1006  CoolPropDbl sqrt_arg = std::abs(B * B / (4 * A * A) - C / A);
1007 
1008  // If the argument to sqrt is very small, we multiply it by a large factor to make it
1009  // larger, and then also divide the sqrt by the sqrt of the factor
1010  if (std::abs(sqrt_arg) > 1e-10) {
1011  DeltavL = -0.5 * B / A + sign((alphaL - alphaV) / alphaV) * sqrt(sqrt_arg);
1012  } else {
1013  // Scale the argument to sqrt() function to make it about 1.0, and divide by sqrt(factor) to yield a factor of 1
1014  CoolPropDbl powerfactor = -log10(sqrt_arg);
1015  DeltavL = -0.5 * B / A + sign((alphaL - alphaV) / alphaV) * sqrt(sqrt_arg * powerfactor) / sqrt(powerfactor);
1016  }
1017  DeltavV = (pL - pV + alphaL * DeltavL) / alphaV;
1018 
1019  // Update the densities of liquid and vapor
1020  rhoL = 1 / (vL + DeltavL);
1021  rhoV = 1 / (vV + DeltavV);
1022  if (B * B / (4 * A * A) - C / A < -10 * DBL_EPSILON) {
1023  rhoL *= 1.01;
1024  rhoV /= 1.01;
1025  }
1026 
1027  // Update the states again
1028  SatL->update_DmolarT_direct(rhoL, T);
1029  SatV->update_DmolarT_direct(rhoV, T);
1030 
1031  // Calculate the error (here the relative error in pressure)
1032  error = std::abs((SatL->p() - SatV->p()) / SatL->p());
1033 
1034  if (get_debug_level() > 5) {
1035  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,
1036  DeltavL / vL, DeltavV / vV, pL, pV);
1037  }
1038 
1039  // If the step size is small, start a counter to allow the other density
1040  // to be corrected a few times
1041  if (std::abs(DeltavL * rhoL) < 1e-13 || std::abs(DeltavV * rhoV) < 1e-13) {
1042  small_step_count++;
1043  }
1044  // If you are not continuing to march towards the solution, after a couple of times, stop
1045  // This is especially a problem for water
1046  if (std::abs(error) > std::abs(last_error)) {
1047  backwards_step_count++;
1048  }
1049 
1050  iter++;
1051  last_error = error;
1052  if (iter > 30) {
1053  throw SolutionError(format("Maxwell solver did not converge after 30 iterations; rhoL: %0.16Lg rhoV: %0.16Lg error: %Lg dvL/vL: %Lg "
1054  "dvV/vV: %Lg pL: %Lg pV: %Lg\n",
1055  rhoL, rhoV, error, DeltavL / vL, DeltavV / vV, pL, pV));
1056  }
1057  } while ((SatL->p() < 0) || (error > 1e-10 && small_step_count < 4 && backwards_step_count < 6));
1058  if (get_debug_level() > 5) {
1059  std::cout << format("[Maxwell] pL: %g pV: %g\n", SatL->p(), SatV->p());
1060  }
1061 }
1062 
1063 void SaturationSolvers::x_and_y_from_K(CoolPropDbl beta, const std::vector<CoolPropDbl>& K, const std::vector<CoolPropDbl>& z,
1064  std::vector<CoolPropDbl>& x, std::vector<CoolPropDbl>& y) {
1065  for (unsigned int i = 0; i < K.size(); i++) {
1066  double denominator = (1 - beta + beta * K[i]); // Common denominator
1067  x[i] = z[i] / denominator;
1068  y[i] = K[i] * z[i] / denominator;
1069  }
1070 }
1071 
1072 void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend& HEOS, const CoolPropDbl beta, CoolPropDbl T, CoolPropDbl p,
1073  const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& K, mixture_VLE_IO& options) {
1074  int iter = 1;
1075  CoolPropDbl change, f, df, deriv_liq, deriv_vap;
1076  std::size_t N = z.size();
1077  std::vector<CoolPropDbl> ln_phi_liq, ln_phi_vap;
1078  ln_phi_liq.resize(N);
1079  ln_phi_vap.resize(N);
1080 
1081  std::vector<CoolPropDbl>&x = HEOS.SatL->get_mole_fractions_ref(), &y = HEOS.SatV->get_mole_fractions_ref();
1082  x_and_y_from_K(beta, K, z, x, y);
1083 
1084  HEOS.SatL->specify_phase(iphase_liquid);
1085  HEOS.SatV->specify_phase(iphase_gas);
1086 
1087  normalize_vector(x);
1088  normalize_vector(y);
1089 
1090  HEOS.SatL->set_mole_fractions(x);
1091  HEOS.SatV->set_mole_fractions(y);
1092  HEOS.SatL->calc_reducing_state();
1093  HEOS.SatV->calc_reducing_state();
1094  CoolPropDbl rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(T, p, iphase_liquid); // [mol/m^3]
1095  CoolPropDbl rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(T, p, iphase_gas); // [mol/m^3]
1096 
1097  // Use Peneloux volume translation to shift liquid volume
1098  // As in Horstmann :: doi:10.1016/j.fluid.2004.11.002
1099  double summer_c = 0, v_SRK = 1 / rhomolar_liq;
1100  const std::vector<CoolPropFluid>& components = HEOS.get_components();
1101  for (std::size_t i = 0; i < components.size(); ++i) {
1102  // Get the parameters for the cubic EOS
1103  CoolPropDbl Tc = HEOS.get_fluid_constant(i, iT_critical);
1104  CoolPropDbl pc = HEOS.get_fluid_constant(i, iP_critical);
1105  CoolPropDbl rhomolarc = HEOS.get_fluid_constant(i, irhomolar_critical);
1106  CoolPropDbl R = 8.3144598;
1107 
1108  summer_c += z[i] * (0.40768 * R * Tc / pc * (0.29441 - pc / (rhomolarc * R * Tc)));
1109  }
1110  rhomolar_liq = 1 / (v_SRK - summer_c);
1111  HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq);
1112  HEOS.SatV->update_TP_guessrho(T, p, rhomolar_vap);
1113 
1114  do {
1115  HEOS.SatL->update_TP_guessrho(T, p, HEOS.SatL->rhomolar());
1116  HEOS.SatV->update_TP_guessrho(T, p, HEOS.SatV->rhomolar());
1117 
1118  f = 0;
1119  df = 0;
1120 
1121  x_N_dependency_flag xN_flag = XN_INDEPENDENT;
1122  for (std::size_t i = 0; i < N; ++i) {
1123  ln_phi_liq[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatL.get()), i, xN_flag);
1124  ln_phi_vap[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatV.get()), i, xN_flag);
1125 
1126  if (options.sstype == imposed_p) {
1127  deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatL.get()), i, xN_flag);
1128  deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatV.get()), i, xN_flag);
1129  } else if (options.sstype == imposed_T) {
1130  deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatL.get()), i, xN_flag);
1131  deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatV.get()), i, xN_flag);
1132  } else {
1133  throw ValueError();
1134  }
1135 
1136  K[i] = exp(ln_phi_liq[i] - ln_phi_vap[i]);
1137 
1138  f += z[i] * (K[i] - 1) / (1 - beta + beta * K[i]);
1139 
1140  double dfdK = K[i] * z[i] / pow(1 - beta + beta * K[i], (int)2);
1141  df += dfdK * (deriv_liq - deriv_vap);
1142  }
1143 
1144  if (std::abs(df) <= 1e-14) { // To avoid dividing by 0
1145  if (std::abs(f) <= 1e-12) // 1e-12 is the loop convergence criterion
1146  {
1147  change = -f; // Should be converged. f <= e-12, so change will have nearly no impact.
1148  } else {
1149  throw ValueError(format("df very small (df = %g) in successive_substitution but f is not converged (f = %g > 1e-12).", df, f));
1150  }
1151  } else {
1152  change = -f / df;
1153  }
1154 
1155  double omega = 1.0;
1156  if (options.sstype == imposed_p) {
1157  T += change;
1158  } else if (options.sstype == imposed_T) {
1159  if (std::abs(change) > 0.05 * p) {
1160  omega = 0.1;
1161  }
1162  p += omega * change;
1163  }
1164 
1165  x_and_y_from_K(beta, K, z, x, y);
1166  normalize_vector(x);
1167  normalize_vector(y);
1168  HEOS.SatL->set_mole_fractions(x);
1169  HEOS.SatV->set_mole_fractions(y);
1170 
1171  iter += 1;
1172  if (iter > 50) {
1173  throw ValueError(format("saturation_p was unable to reach a solution within 50 iterations"));
1174  }
1175  } while (std::abs(f) > 1e-12 && iter < options.Nstep_max);
1176 
1177  HEOS.SatL->update_TP_guessrho(T, p, HEOS.SatL->rhomolar());
1178  HEOS.SatV->update_TP_guessrho(T, p, HEOS.SatV->rhomolar());
1179 
1180  options.p = HEOS.SatL->p();
1181  options.T = HEOS.SatL->T();
1182  options.rhomolar_liq = HEOS.SatL->rhomolar();
1183  options.rhomolar_vap = HEOS.SatV->rhomolar();
1184  options.x = x;
1185  options.y = y;
1186 }
1187 void SaturationSolvers::newton_raphson_saturation::resize(std::size_t N) {
1188  this->N = N;
1189  x.resize(N);
1190  y.resize(N);
1191 
1192  if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1193  r.resize(N + 1);
1194  err_rel.resize(N + 1);
1195  J.resize(N + 1, N + 1);
1196  } else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1197  r.resize(N);
1198  err_rel.resize(N);
1199  J.resize(N, N);
1200  } else {
1201  throw ValueError();
1202  }
1203 }
1204 void SaturationSolvers::newton_raphson_saturation::check_Jacobian() {
1205  // References to the classes for concision
1206  HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
1207 
1208  // Build the Jacobian and residual vectors
1209  build_arrays();
1210 
1211  // Make copies of the base
1212  CoolPropDbl T0 = T;
1213  std::vector<CoolPropDbl> x0 = x;
1214  Eigen::VectorXd r0 = r;
1215  Eigen::MatrixXd J0 = J;
1216  CoolPropDbl rhomolar_liq0 = rSatL.rhomolar();
1217  CoolPropDbl rhomolar_vap0 = rSatV.rhomolar();
1218 
1219  {
1220  // Derivatives with respect to T
1221  double dT = 1e-3, T1 = T + dT, T2 = T - dT;
1222  this->T = T1;
1223  this->rhomolar_liq = rhomolar_liq0;
1224  this->rhomolar_vap = rhomolar_vap0;
1225  build_arrays();
1226  Eigen::VectorXd r1 = r;
1227  this->T = T2;
1228  this->rhomolar_liq = rhomolar_liq0;
1229  this->rhomolar_vap = rhomolar_vap0;
1230  build_arrays();
1231  Eigen::VectorXd r2 = r;
1232 
1233  Eigen::VectorXd diffn = (r1 - r2) / (2 * dT);
1234  std::cout << format("For T\n");
1235  //std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
1236  //std::cout << "analytic: " << vec_to_string(J0.col(N-1), "%0.11Lg") << std::endl;
1237  }
1238  {
1239  // Derivatives with respect to rho'
1240  double drho = 1;
1241  this->T = T0;
1242  this->rhomolar_liq = rhomolar_liq0 + drho;
1243  this->rhomolar_vap = rhomolar_vap0;
1244  build_arrays();
1245  Eigen::VectorXd rr1 = r;
1246  this->T = T0;
1247  this->rhomolar_liq = rhomolar_liq0 - drho;
1248  this->rhomolar_vap = rhomolar_vap0;
1249  build_arrays();
1250  Eigen::VectorXd rr2 = r;
1251 
1252  Eigen::VectorXd diffn = (rr1 - rr2) / (2 * drho);
1253  std::cout << format("For rho\n");
1254  //std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
1255  //std::cout << "analytic: " << vec_to_string(J0.col(N-1), "%0.11Lg") << std::endl;
1256  }
1257  for (std::size_t i = 0; i < x.size() - 1; ++i) {
1258  // Derivatives with respect to x[i]
1259  double dx = 1e-5;
1260  this->x = x0;
1261  this->x[i] += dx;
1262  this->x[x.size() - 1] -= dx;
1263  this->T = T0;
1264  this->rhomolar_liq = rhomolar_liq0;
1265  this->rhomolar_vap = rhomolar_vap0;
1266  build_arrays();
1267  Eigen::VectorXd r1 = this->r;
1268 
1269  this->x = x0;
1270  this->x[i] -= dx;
1271  this->x[x.size() - 1] += dx;
1272  this->T = T0;
1273  this->rhomolar_liq = rhomolar_liq0;
1274  this->rhomolar_vap = rhomolar_vap0;
1275  build_arrays();
1276  Eigen::VectorXd r2 = this->r;
1277 
1278  Eigen::VectorXd diffn = (r1 - r2) / (2 * dx);
1279  std::cout << format("For x%d N %d\n", i, N);
1280  //std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
1281  //std::cout << "analytic: " << vec_to_string(J0.col(i), "%0.11Lg") << std::endl;
1282  }
1283 }
1284 void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBackend& HEOS, const std::vector<CoolPropDbl>& z,
1285  std::vector<CoolPropDbl>& z_incipient, newton_raphson_saturation_options& IO) {
1286  int iter = 0;
1287  bool debug = get_debug_level() > 9 || false;
1288 
1289  if (debug) {
1290  std::cout << " NRsat::call: p " << IO.p << " T " << IO.T << " dl " << IO.rhomolar_liq << " dv " << IO.rhomolar_vap << std::endl;
1291  }
1292 
1293  // Reset all the variables and resize
1294  pre_call();
1295 
1296  this->bubble_point = IO.bubble_point;
1297  rhomolar_liq = IO.rhomolar_liq;
1298  rhomolar_vap = IO.rhomolar_vap;
1299  T = IO.T;
1300  p = IO.p;
1301  imposed_variable = IO.imposed_variable;
1302 
1303  resize(z.size());
1304 
1305  if (bubble_point) {
1306  // Bubblepoint, vapor (y) is the incipient phase
1307  x = z;
1308  y = z_incipient;
1309  } else {
1310  // Dewpoint, liquid (x) is the incipient phase
1311  y = z;
1312  x = z_incipient;
1313  }
1314 
1315  // Hold a pointer to the backend
1316  this->HEOS = &HEOS;
1317 
1318  //check_Jacobian();
1319 
1320  do {
1321  // Build the Jacobian and residual vectors
1322  build_arrays();
1323 
1324  // Solve for the step; v is the step with the contents
1325  // [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]
1326  Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
1327 
1328  if (bubble_point) {
1329  for (unsigned int i = 0; i < N - 1; ++i) {
1330  err_rel[i] = v[i] / y[i];
1331  y[i] += v[i];
1332  }
1333  y[N - 1] = 1 - std::accumulate(y.begin(), y.end() - 1, 0.0);
1334  } else {
1335  for (unsigned int i = 0; i < N - 1; ++i) {
1336  err_rel[i] = v[i] / x[i];
1337  x[i] += v[i];
1338  }
1339  x[N - 1] = 1 - std::accumulate(x.begin(), x.end() - 1, 0.0);
1340  }
1341  if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED) {
1342  T += v[N - 1];
1343  err_rel[N - 1] = v[N - 1] / T;
1344  } else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1345  p += v[N - 1];
1346  err_rel[N - 1] = v[N - 1] / p;
1347  } else if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1348  T += v[N - 1];
1349  err_rel[N - 1] = v[N - 1] / T;
1350  rhomolar_liq += v[N];
1351  err_rel[N] = v[N] / rhomolar_liq;
1352  } else {
1353  throw ValueError("invalid imposed_variable");
1354  }
1355  if (debug) {
1356  //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;
1357  }
1358 
1359  min_rel_change = err_rel.cwiseAbs().minCoeff();
1360  iter++;
1361 
1362  if (iter == IO.Nstep_max) {
1363  throw ValueError(format("newton_raphson_saturation::call reached max number of iterations [%d]", IO.Nstep_max));
1364  }
1365  } while (this->error_rms > 1e-7 && min_rel_change > 1000 * DBL_EPSILON && iter < IO.Nstep_max);
1366 
1367  IO.Nsteps = iter;
1368  IO.p = p;
1369  IO.x = x; // Mole fractions in liquid
1370  IO.y = y; // Mole fractions in vapor
1371  IO.T = T;
1372  IO.rhomolar_liq = rhomolar_liq;
1373  IO.rhomolar_vap = rhomolar_vap;
1374  const std::vector<CoolPropFluid>& fluidsL = HEOS.SatL->get_components();
1375  const std::vector<CoolPropFluid>& fluidsV = HEOS.SatV->get_components();
1376  if (!fluidsL.empty() && !fluidsV.empty()) {
1377  IO.hmolar_liq = HEOS.SatL->hmolar();
1378  IO.hmolar_vap = HEOS.SatV->hmolar();
1379  IO.smolar_liq = HEOS.SatL->smolar();
1380  IO.smolar_vap = HEOS.SatV->smolar();
1381  }
1382 }
1383 
1384 void SaturationSolvers::newton_raphson_saturation::build_arrays() {
1385  // References to the classes for concision
1386  HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
1387 
1388  // Step 0:
1389  // -------
1390  // Set mole fractions for the incipient phase
1391  if (bubble_point) {
1392  // Vapor is incipient phase, set its composition
1393  rSatV.set_mole_fractions(y);
1394  rSatL.set_mole_fractions(x);
1395  } else {
1396  // Liquid is incipient phase, set its composition
1397  rSatL.set_mole_fractions(x);
1398  rSatV.set_mole_fractions(y);
1399  }
1400 
1401  if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1402  rSatL.update(DmolarT_INPUTS, rhomolar_liq, T);
1403  rSatV.update(DmolarT_INPUTS, rhomolar_vap, T);
1404  } else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1405  rSatL.update_TP_guessrho(T, p, rhomolar_liq);
1406  rhomolar_liq = rSatL.rhomolar();
1407  rSatV.update_TP_guessrho(T, p, rhomolar_vap);
1408  rhomolar_vap = rSatV.rhomolar();
1409  } else {
1410  throw ValueError("imposed variable not set for NR VLE");
1411  }
1412 
1413  // For diagnostic purposes calculate the pressures (no derivatives are evaluated)
1414  CoolPropDbl p_liq = rSatL.p();
1415  CoolPropDbl p_vap = rSatV.p();
1416  p = 0.5 * (p_liq + p_vap);
1417 
1418  // Step 2:
1419  // -------
1420  // Build the residual vector and the Jacobian matrix
1421 
1422  x_N_dependency_flag xN_flag = XN_DEPENDENT;
1423 
1424  if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1425  // For the residuals F_i (equality of fugacities)
1426  for (std::size_t i = 0; i < N; ++i) {
1427  // Equate the liquid and vapor fugacities
1428  CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1429  CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1430  r(i) = ln_f_liq - ln_f_vap;
1431 
1432  for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1433  if (bubble_point) {
1434  J(i, j) = -MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatV, i, j, xN_flag);
1435  } else {
1436  J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatL, i, j, xN_flag);
1437  }
1438  }
1439  J(i, N - 1) = MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatL, i, xN_flag)
1440  - MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatV, i, xN_flag);
1441  J(i, N) = MixtureDerivatives::dln_fugacity_i_drho__constT_n(rSatL, i, xN_flag);
1442  }
1443  // ---------------------------------------------------------------
1444  // Derivatives of pL(T,rho',x)-p(T,rho'',y) with respect to inputs
1445  // ---------------------------------------------------------------
1446  r(N) = p_liq - p_vap;
1447  for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1448  J(N, j) = MixtureDerivatives::dpdxj__constT_V_xi(rSatL, j, xN_flag); // p'' not a function of x0
1449  }
1450  // Fixed composition derivatives
1451  J(N, N - 1) = rSatL.first_partial_deriv(iP, iT, iDmolar) - rSatV.first_partial_deriv(iP, iT, iDmolar);
1452  J(N, N) = rSatL.first_partial_deriv(iP, iDmolar, iT);
1453  } else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED) {
1454  // Independent variables are N-1 mole fractions of incipient phase and T
1455 
1456  // For the residuals F_i (equality of fugacities)
1457  for (std::size_t i = 0; i < N; ++i) {
1458  // Equate the liquid and vapor fugacities
1459  CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1460  CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1461  r(i) = ln_f_liq - ln_f_vap;
1462 
1463  for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1464  if (bubble_point) {
1465  J(i, j) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
1466  } else {
1467  J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
1468  }
1469  }
1470  J(i, N - 1) =
1471  MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag);
1472  }
1473  } else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1474  // Independent variables are N-1 mole fractions of incipient phase and p
1475 
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_p_xi(rSatV, i, j, xN_flag);
1486  } else {
1487  J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
1488  }
1489  }
1490  J(i, N - 1) =
1491  MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatV, i, xN_flag);
1492  }
1493  } else {
1494  throw ValueError();
1495  }
1496 
1497  error_rms = r.norm();
1498 
1499  // Calculate derivatives along phase boundary;
1500  // Gernert thesis 3.96 and 3.97
1501  CoolPropDbl dQ_dPsat = 0, dQ_dTsat = 0;
1502  for (std::size_t i = 0; i < N; ++i) {
1503  dQ_dPsat += x[i]
1504  * (MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rSatL, i, xN_flag)
1505  - MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rSatV, i, xN_flag));
1506  dQ_dTsat += x[i]
1507  * (MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatL, i, xN_flag)
1508  - MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatV, i, xN_flag));
1509  }
1510  dTsat_dPsat = -dQ_dPsat / dQ_dTsat;
1511  dPsat_dTsat = -dQ_dTsat / dQ_dPsat;
1512 }
1513 
1514 void SaturationSolvers::newton_raphson_twophase::call(HelmholtzEOSMixtureBackend& HEOS, newton_raphson_twophase_options& IO) {
1515  int iter = 0;
1516 
1517  if (get_debug_level() > 9) {
1518  std::cout << " NRsat::call: p" << IO.p << " T" << IO.T << " dl" << IO.rhomolar_liq << " dv" << IO.rhomolar_vap << std::endl;
1519  }
1520 
1521  // Reset all the variables and resize
1522  pre_call();
1523 
1524  rhomolar_liq = IO.rhomolar_liq;
1525  rhomolar_vap = IO.rhomolar_vap;
1526  T = IO.T;
1527  p = IO.p;
1528  imposed_variable = IO.imposed_variable;
1529  x = IO.x;
1530  y = IO.y;
1531  z = IO.z;
1532  beta = IO.beta;
1533 
1534  this->N = z.size();
1535  x.resize(N);
1536  y.resize(N);
1537  r.resize(2 * N - 1);
1538  J.resize(2 * N - 1, 2 * N - 1);
1539  err_rel.resize(2 * N - 1);
1540 
1541  // Hold a pointer to the backend
1542  this->HEOS = &HEOS;
1543 
1544  do {
1545  // Build the Jacobian and residual vectors
1546  build_arrays();
1547 
1548  // Solve for the step; v is the step with the contents
1549  // [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]
1550 
1551  // Uncomment to see Jacobian and residual at every step
1552  // std::cout << vec_to_string(J, "%0.12Lg") << std::endl;
1553  // std::cout << vec_to_string(negative_r, "%0.12Lg") << std::endl;
1554 
1555  Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
1556 
1557  for (unsigned int i = 0; i < N - 1; ++i) {
1558  err_rel[i] = v[i] / x[i];
1559  x[i] += v[i];
1560  err_rel[i + (N - 1)] = v[i + (N - 1)] / y[i];
1561  y[i] += v[i + (N - 1)];
1562  }
1563  x[N - 1] = 1 - std::accumulate(x.begin(), x.end() - 1, 0.0);
1564  y[N - 1] = 1 - std::accumulate(y.begin(), y.end() - 1, 0.0);
1565 
1566  if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED) {
1567  T += v[2 * N - 2];
1568  err_rel[2 * N - 2] = v[2 * N - 2] / T;
1569  } else if (imposed_variable == newton_raphson_twophase_options::T_IMPOSED) {
1570  p += v[2 * N - 2];
1571  err_rel[2 * N - 2] = v[2 * N - 2] / p;
1572  } else {
1573  throw ValueError("invalid imposed_variable");
1574  }
1575  //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;
1576 
1577  min_rel_change = err_rel.cwiseAbs().minCoeff();
1578  iter++;
1579 
1580  if (iter == IO.Nstep_max) {
1581  throw ValueError(format("newton_raphson_saturation::call reached max number of iterations [%d]", IO.Nstep_max));
1582  }
1583  } while (this->error_rms > 1e-9 && min_rel_change > 1000 * DBL_EPSILON && iter < IO.Nstep_max);
1584 
1585  IO.Nsteps = iter;
1586  IO.p = p;
1587  IO.x = x; // Mole fractions in liquid
1588  IO.y = y; // Mole fractions in vapor
1589  IO.T = T;
1590  IO.rhomolar_liq = rhomolar_liq;
1591  IO.rhomolar_vap = rhomolar_vap;
1592  IO.hmolar_liq = HEOS.SatL.get()->hmolar();
1593  IO.hmolar_vap = HEOS.SatV.get()->hmolar();
1594  IO.smolar_liq = HEOS.SatL.get()->smolar();
1595  IO.smolar_vap = HEOS.SatV.get()->smolar();
1596 }
1597 
1598 void SaturationSolvers::newton_raphson_twophase::build_arrays() {
1599  // References to the classes for concision
1600  HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
1601 
1602  // Step 0:
1603  // -------
1604  // Set mole fractions
1605  rSatL.set_mole_fractions(x);
1606  rSatV.set_mole_fractions(y);
1607 
1608  //std::vector<CoolPropDbl> &x = rSatL.get_mole_fractions();
1609  //std::vector<CoolPropDbl> &y = rSatV.get_mole_fractions();
1610 
1611  rSatL.update_TP_guessrho(T, p, rhomolar_liq);
1612  rhomolar_liq = rSatL.rhomolar();
1613  rSatV.update_TP_guessrho(T, p, rhomolar_vap);
1614  rhomolar_vap = rSatV.rhomolar();
1615 
1616  // For diagnostic purposes calculate the pressures (no derivatives are evaluated)
1617  CoolPropDbl p_liq = rSatL.p();
1618  CoolPropDbl p_vap = rSatV.p();
1619  p = 0.5 * (p_liq + p_vap);
1620 
1621  // Step 2:
1622  // -------
1623  // Build the residual vector and the Jacobian matrix
1624 
1625  x_N_dependency_flag xN_flag = XN_DEPENDENT;
1626 
1627  // Form of residuals do not depend on which variable is imposed
1628  for (std::size_t i = 0; i < N; ++i) {
1629  // Equate the liquid and vapor fugacities
1630  CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1631  CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1632  r[i] = ln_f_liq - ln_f_vap; // N of these
1633 
1634  if (i != N - 1) {
1635  // Equate the specified vapor mole fraction and that given defined by the ith component
1636  r[i + N] = (z[i] - x[i]) / (y[i] - x[i]) - beta; // N-1 of these
1637  }
1638  }
1639 
1640  // First part of derivatives with respect to ln f_i
1641  for (std::size_t i = 0; i < N; ++i) {
1642  for (std::size_t j = 0; j < N - 1; ++j) {
1643  J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
1644  J(i, j + N - 1) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
1645  }
1646 
1647  // Last derivative with respect to either T or p depending on what is imposed
1648  if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED) {
1649  J(i, 2 * N - 2) =
1650  MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag);
1651  } else if (imposed_variable == newton_raphson_twophase_options::T_IMPOSED) {
1652  J(i, 2 * N - 2) =
1653  MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatV, i, xN_flag);
1654  } else {
1655  throw ValueError();
1656  }
1657  }
1658  // Derivatives with respect to the vapor mole fractions residual
1659  for (std::size_t i = 0; i < N - 1; ++i) {
1660  std::size_t k = i + N; // N ln f_i residuals
1661  J(k, i) = (z[i] - y[i]) / pow(y[i] - x[i], 2);
1662  J(k, i + (N - 1)) = -(z[i] - x[i]) / pow(y[i] - x[i], 2);
1663  }
1664 
1665  error_rms = r.norm(); // Square-root (The R in RMS)
1666 }
1667 
1668 class RachfordRiceResidual : public FuncWrapper1DWithDeriv
1669 {
1670  private:
1671  const std::vector<double>&z, &lnK;
1672 
1673  public:
1674  RachfordRiceResidual(const std::vector<double>& z, const std::vector<double>& lnK) : z(z), lnK(lnK){};
1675  double call(double beta) {
1676  return FlashRoutines::g_RachfordRice(z, lnK, beta);
1677  }
1678  double deriv(double beta) {
1679  return FlashRoutines::dgdbeta_RachfordRice(z, lnK, beta);
1680  }
1681 };
1682 
1683 void StabilityRoutines::StabilityEvaluationClass::trial_compositions() {
1684 
1685  x.resize(z.size());
1686  y.resize(z.size());
1687  lnK.resize(z.size());
1688  K.resize(z.size());
1689  double g0 = 0, g1 = 0, beta = -1;
1690 
1691  for (int i = 0; i < static_cast<int>(z.size()); ++i) {
1692  // Calculate the K-factor
1693  if (m_T < 0 && m_p < 0) {
1694  // Using T&P from the class
1695  lnK[i] = SaturationSolvers::Wilson_lnK_factor(HEOS, HEOS.T(), HEOS.p(), i);
1696  } else {
1697  // Using specified T&P
1698  lnK[i] = SaturationSolvers::Wilson_lnK_factor(HEOS, m_T, m_p, i);
1699  }
1700  K[i] = exp(lnK[i]);
1701  g0 += z[i] * (K[i] - 1); // The summation for beta = 0
1702  g1 += z[i] * (1 - 1 / K[i]); // The summation for beta = 1
1703  }
1704  // Copy K-factors for later use
1705  K0 = K;
1706  // Now see what to do about the g(0) and g(1) values
1707  // -----
1708  //
1709  if (g0 < 0) {
1710  beta = 0; // Assumed to be at bubble-point temperature
1711  } else if (g1 > 0) {
1712  beta = 1; // Assumed to be at the dew-point temperature
1713  } else {
1714  // Need to iterate to find beta that makes g of Rachford-Rice zero
1715  RachfordRiceResidual resid(z, lnK);
1716  beta = Brent(resid, 0, 1, DBL_EPSILON, 1e-10, 100);
1717  }
1718  // Get the compositions from given value for beta, K, z
1719  SaturationSolvers::x_and_y_from_K(beta, K, z, x, y);
1720  normalize_vector(x);
1721  normalize_vector(y);
1722  if (debug) {
1723  std::cout << format("1) T: %g p: %g beta: %g\n", HEOS.T(), HEOS.p(), beta);
1724  }
1725 }
1726 void StabilityRoutines::StabilityEvaluationClass::successive_substitution(int num_steps) {
1727  // ----
1728  // Do a few steps of successive substitution
1729  // ----
1730 
1731  HEOS.SatL->set_mole_fractions(x);
1732  HEOS.SatL->calc_reducing_state();
1733  HEOS.SatV->set_mole_fractions(y);
1734  HEOS.SatV->calc_reducing_state();
1735 
1736  if (debug) {
1737  std::cout << format("2) SS1: i beta K x y rho' rho''\n");
1738  }
1739  for (int step_count = 0; step_count < num_steps; ++step_count) {
1740  // Set the composition
1741  HEOS.SatL->set_mole_fractions(x);
1742  HEOS.SatV->set_mole_fractions(y);
1743  HEOS.SatL->calc_reducing_state();
1744  HEOS.SatV->calc_reducing_state();
1745 
1746  this->rho_TP_global();
1747 
1748  // Calculate the new K-factors from the fugacity coefficients
1749  double g0 = 0, g1 = 0;
1750  for (std::size_t i = 0; i < z.size(); ++i) {
1751  lnK[i] = log(HEOS.SatL->fugacity_coefficient(i) / HEOS.SatV->fugacity_coefficient(i));
1752  K[i] = exp(lnK[i]);
1753  g0 += z[i] * (K[i] - 1); // The summation for beta = 0
1754  g1 += z[i] * (1 - 1 / K[i]); // The summation for beta = 1
1755  }
1756  RachfordRiceResidual resid(z, lnK);
1757  if (g0 < 0) {
1758  beta = 0;
1759  } else if (g1 > 0) {
1760  beta = 1;
1761  } else {
1762  // Need to iterate to find beta that makes g of Rachford-Rice zero
1763  beta = Brent(resid, 0, 1, DBL_EPSILON, 1e-10, 100);
1764  }
1765 
1766  // Get the compositions from given values for beta, K, z
1767  SaturationSolvers::x_and_y_from_K(beta, K, z, x, y);
1768  normalize_vector(x);
1769  normalize_vector(y);
1770  if (debug) {
1771  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(),
1772  vec_to_string(y, "%0.6f").c_str(), rhomolar_liq, rhomolar_vap);
1773  }
1774  }
1775 }
1776 void StabilityRoutines::StabilityEvaluationClass::check_stability() {
1777  std::vector<double> tpdL, tpdH;
1778 
1779  // Calculate the temperature and pressure to be used
1780  double the_T = (m_T > 0 && m_p > 0) ? m_T : HEOS.T();
1781  double the_p = (m_T > 0 && m_p > 0) ? m_p : HEOS.p();
1782 
1783  // If beta value is between epsilon and 1-epsilon, check the TPD
1784  if (beta > DBL_EPSILON && beta < 1 - DBL_EPSILON) {
1785 
1786  // Set the composition back to the bulk composition for both liquid and vapor phases
1787  HEOS.SatL->set_mole_fractions(z);
1788  HEOS.SatV->set_mole_fractions(z);
1789  HEOS.SatL->calc_reducing_state();
1790  HEOS.SatV->calc_reducing_state();
1791 
1792  // Update the densities in each class
1793  double rhoL = HEOS.SatL->solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SatL->SRK_covolume());
1794  double rhoV = HEOS.SatV->solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SatV->SRK_covolume());
1795  HEOS.SatL->update_DmolarT_direct(rhoL, the_T);
1796  HEOS.SatV->update_DmolarT_direct(rhoV, the_T);
1797 
1798  // Calculate the tpd and the Gibbs energy difference (Gernert, 2014, Eqs. 20-22)
1799  // The trial compositions are the phase compositions from before
1800  this->tpd_liq = HEOS.SatL->tangent_plane_distance(the_T, the_p, x, rhomolar_liq);
1801  this->tpd_vap = HEOS.SatV->tangent_plane_distance(the_T, the_p, y, rhomolar_vap);
1802 
1803  this->DELTAG_nRT = (1 - beta) * tpd_liq + beta * (tpd_vap);
1804  if (debug) {
1805  std::cout << format("3) tpd': %g tpd'': %g DELTAG/nRT: %g\n", tpd_liq, tpd_vap, DELTAG_nRT);
1806  }
1807 
1808  // If any of these cases are met, feed is conclusively unstable, stop!
1809  if (this->tpd_liq < -DBL_EPSILON || this->tpd_vap < -DBL_EPSILON || this->DELTAG_nRT < -DBL_EPSILON) {
1810  if (debug) {
1811  std::cout << format("3) PHASE SPLIT beta in (eps,1-eps) \n");
1812  }
1813  _stable = false;
1814  return;
1815  }
1816  }
1817 
1818  // Ok, we aren't sure about stability, need to keep going with the full tpd analysis
1819 
1820  // Use the global density solver to obtain the density root (or the lowest Gibbs energy root if more than one)
1821  CoolPropDbl rho_bulk = HEOS.solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SRK_covolume());
1822  HEOS.update_DmolarT_direct(rho_bulk, the_T);
1823 
1824  // Calculate the fugacity coefficient at initial composition of the bulk phase
1825  std::vector<double> fugacity_coefficient0(z.size()), fugacity0(z.size());
1826  for (std::size_t i = 0; i < z.size(); ++i) {
1827  fugacity_coefficient0[i] = HEOS.fugacity_coefficient(i);
1828  fugacity0[i] = HEOS.fugacity(i);
1829  }
1830 
1831  // Generate light and heavy test compositions (Gernert, 2014, Eq. 23)
1832  xL.resize(z.size());
1833  xH.resize(z.size());
1834  for (std::size_t i = 0; i < z.size(); ++i) {
1835  xL[i] = z[i] * K0[i]; // Light-phase composition
1836  xH[i] = z[i] / K0[i]; // Heavy-phase composition
1837  }
1838  normalize_vector(xL);
1839  normalize_vector(xH);
1840 
1841  // For each composition, use successive substitution to try to evaluate stability
1842  if (debug) {
1843  std::cout << format("3) SS2: i x' x'' rho' rho'' tpd' tpd''\n");
1844  }
1845 
1846  // We got this far, we assume stable phases
1847  _stable = true;
1848 
1849  double diffbulkL = 0, diffbulkH = 0;
1850  for (int step_count = 0; step_count < 100; ++step_count) {
1851 
1852  // Set the composition
1853  HEOS.SatL->set_mole_fractions(xH);
1854  HEOS.SatV->set_mole_fractions(xL);
1855  HEOS.SatL->calc_reducing_state();
1856  HEOS.SatV->calc_reducing_state();
1857 
1858  // Do the global density solver for both phases
1859  rho_TP_global();
1860 
1861  double tpd_L = 0, tpd_H = 0;
1862  for (std::size_t i = 0; i < xL.size(); ++i) {
1863  tpd_L += xL[i] * (log(MixtureDerivatives::fugacity_i(*HEOS.SatV, i, XN_DEPENDENT)) - log(fugacity0[i]));
1864  tpd_H += xH[i] * (log(MixtureDerivatives::fugacity_i(*HEOS.SatL, i, XN_DEPENDENT)) - log(fugacity0[i]));
1865  }
1866  tpdL.push_back(tpd_L);
1867  tpdH.push_back(tpd_H);
1868 
1869  // Calculate the new composition from the fugacity coefficients
1870  diffbulkL = 0, diffbulkH = 0;
1871  for (std::size_t i = 0; i < z.size(); ++i) {
1872  xL[i] = z[i] * fugacity_coefficient0[i] / HEOS.SatV->fugacity_coefficient(i);
1873  diffbulkL += std::abs(xL[i] - z[i]);
1874  xH[i] = z[i] * fugacity_coefficient0[i] / HEOS.SatL->fugacity_coefficient(i);
1875  diffbulkH += std::abs(xH[i] - z[i]);
1876  }
1877  normalize_vector(xL);
1878  normalize_vector(xH);
1879  if (debug) {
1880  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(),
1881  rhomolar_liq, rhomolar_vap, tpd_L, tpd_H);
1882  }
1883 
1884  // Check if either of the phases have the bulk composition. If so, no phase split
1885  if (diffbulkL < 1e-2 || diffbulkH < 1e-2) {
1886  _stable = true;
1887  return;
1888  }
1889 
1890  // Check if either tpd is negative, if so, phases definitively split, quit
1891  if (tpd_L < -1e-12 || tpd_H < -1e-12) {
1892  _stable = false;
1893  return;
1894  }
1895  }
1896  if (diffbulkH > 0.25 || diffbulkL > 0.25) {
1897  // At least one test phase is definitely not the bulk composition, so phase split predicted
1898  _stable = false;
1899  }
1900 }
1901 
1902 void StabilityRoutines::StabilityEvaluationClass::rho_TP_global() {
1903 
1904  // Calculate the temperature and pressure to be used
1905  double the_T = (m_T > 0 && m_p > 0) ? m_T : HEOS.T();
1906  double the_p = (m_T > 0 && m_p > 0) ? m_p : HEOS.p();
1907 
1908  // Calculate covolume of SRK, use it as the maximum density
1909  double rhoL = HEOS.SatL->solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SatL->SRK_covolume());
1910  double rhoV = HEOS.SatV->solver_rho_Tp_global(the_T, the_p, 0.9 / HEOS.SatV->SRK_covolume());
1911  HEOS.SatL->update_DmolarT_direct(rhoL, the_T);
1912  HEOS.SatV->update_DmolarT_direct(rhoV, the_T);
1913 
1914  rhomolar_liq = HEOS.SatL->rhomolar();
1915  rhomolar_vap = HEOS.SatV->rhomolar();
1916 }
1917 
1918 void StabilityRoutines::StabilityEvaluationClass::rho_TP_w_guesses() {
1919 
1920  // Re-calculate the density
1921  if (m_T > 0 && m_p > 0) {
1922  HEOS.SatL->update_TP_guessrho(m_T, m_p, rhomolar_liq);
1923  HEOS.SatV->update_TP_guessrho(m_T, m_p, rhomolar_vap);
1924  } else {
1925  HEOS.SatL->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_liq);
1926  HEOS.SatV->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_vap);
1927  }
1928  rhomolar_liq = HEOS.SatL->rhomolar();
1929  rhomolar_vap = HEOS.SatV->rhomolar();
1930 }
1931 
1932 void StabilityRoutines::StabilityEvaluationClass::rho_TP_SRK_translated() {
1933 
1934  // First use cubic as a guess for the density of liquid and vapor phases
1935  if (m_T > 0 && m_p > 0) {
1936  rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(m_T, m_p, iphase_liquid); // [mol/m^3]
1937  rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(m_T, m_p, iphase_gas); // [mol/m^3]
1938  } else {
1939  rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(HEOS.T(), HEOS.p(), iphase_liquid); // [mol/m^3]
1940  rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(HEOS.T(), HEOS.p(), iphase_gas); // [mol/m^3]
1941  }
1942 
1943  // Apply volume translation to liquid density only
1944  if (HEOS.backend_name().find("Helmholtz") == 0) {
1945  // Use Peneloux volume translation to shift liquid volume
1946  // As in Horstmann :: doi:10.1016/j.fluid.2004.11.002
1947  double summer_c = 0, v_SRK = 1 / rhomolar_liq;
1948  for (std::size_t i = 0; i < z.size(); ++i) {
1949  // Get the parameters for the cubic EOS
1950  CoolPropDbl Tc = HEOS.get_fluid_constant(i, iT_critical), pc = HEOS.get_fluid_constant(i, iP_critical),
1951  rhomolarc = HEOS.get_fluid_constant(i, irhomolar_critical);
1952  CoolPropDbl R = 8.3144598;
1953  summer_c += z[i] * (0.40768 * R * Tc / pc * (0.29441 - pc / (rhomolarc * R * Tc)));
1954  }
1955  rhomolar_liq = 1 / (v_SRK - summer_c);
1956  }
1957 }
1958 
1959 void SaturationSolvers::PTflash_twophase::solve() {
1960  const std::size_t N = IO.x.size();
1961  int iter = 0;
1962  double min_rel_change;
1963  do {
1964  // Build the Jacobian and residual vectors
1965  build_arrays();
1966 
1967  // Solve for the step; v is the step with the contents
1968  // [delta(x'_0), delta(x'_1), ..., delta(x'_{N-1}), delta(x''_0), delta(x''_1), ..., delta(x''_{N-1})]
1969 
1970  // Uncomment to see Jacobian and residual at every step
1971  //std::cout << vec_to_string(J, "%0.12Lg") << std::endl;
1972  //std::cout << vec_to_string(-r, "%0.12Lg") << std::endl;
1973 
1974  Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
1975 
1976  for (unsigned int i = 0; i < N - 1; ++i) {
1977  err_rel[i] = v[i] / IO.x[i];
1978  IO.x[i] += v[i];
1979  err_rel[i + (N - 1)] = v[i + (N - 1)] / IO.y[i];
1980  IO.y[i] += v[i + (N - 1)];
1981  }
1982  IO.x[N - 1] = 1 - std::accumulate(IO.x.begin(), IO.x.end() - 1, 0.0);
1983  IO.y[N - 1] = 1 - std::accumulate(IO.y.begin(), IO.y.end() - 1, 0.0);
1984 
1985  //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;
1986 
1987  min_rel_change = err_rel.cwiseAbs().minCoeff();
1988  iter++;
1989 
1990  if (iter == IO.Nstep_max) {
1991  throw ValueError(format("PTflash_twophase::call reached max number of iterations [%d]", IO.Nstep_max));
1992  }
1993  } while (this->error_rms > 1e-9 && min_rel_change > 1000 * DBL_EPSILON && iter < IO.Nstep_max);
1994 }
1995 void SaturationSolvers::PTflash_twophase::build_arrays() {
1996  const std::size_t N = IO.x.size();
1997 
1998  r.resize(2 * N - 2);
1999  J.resize(2 * N - 2, 2 * N - 2);
2000  err_rel.resize(2 * N - 2);
2001 
2002  HEOS.SatL->set_mole_fractions(IO.x);
2003  HEOS.SatL->update_TP_guessrho(IO.T, IO.p, IO.rhomolar_liq);
2004 
2005  HEOS.SatV->set_mole_fractions(IO.y);
2006  HEOS.SatV->update_TP_guessrho(IO.T, IO.p, IO.rhomolar_vap);
2007 
2008  // Independent variables are
2009  // [delta(x'_0), delta(x'_1), ..., delta(x'_{N-1}), delta(x''_0), delta(x''_1), ..., delta(x''_{N-1})]
2010 
2011  // First N residuals are the iso-fugacity condition
2012  for (std::size_t k = 0; k < N; ++k) {
2013  r(k) = log(HEOS.SatL->fugacity(k) / HEOS.SatV->fugacity(k));
2014  for (std::size_t j = 0; j < N - 1; ++j) {
2015  if (k == N - 1) {
2016  J(k, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatL.get()), k, j, XN_DEPENDENT);
2017  J(k, j + N - 1) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatV.get()), k, j, XN_DEPENDENT);
2018  } else {
2019  J(k, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatL.get()), k, j, XN_DEPENDENT);
2020  J(k, j + N - 1) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatV.get()), k, j, XN_DEPENDENT);
2021  }
2022  }
2023  }
2024  // Next N-2 residuals are amount of substance balances
2025  for (std::size_t i = 0; i < N - 2; ++i) {
2026  std::size_t k = i + N;
2027  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]);
2028  for (std::size_t j = 0; j < N - 2; ++j) {
2029  J(k, j) = (IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
2030  J(k, j + N - 1) = -(IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
2031  }
2032  std::size_t j = N - 2;
2033  J(k, j) = -(IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
2034  J(k, j + N - 1) = +(IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
2035  }
2036  this->error_rms = r.norm();
2037 }
2038 } /* namespace CoolProp*/
2039 
2040 #if defined(ENABLE_CATCH)
2041 # include <catch2/catch_all.hpp>
2042 
2043 TEST_CASE("Check the PT flash calculation for two-phase inputs", "[PTflash_twophase]") {
2044  shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", "Propane&Ethane"));
2045  AS->set_mole_fractions(std::vector<double>(2, 0.5));
2046  // Dewpoint calculation
2047  AS->update(CoolProp::PQ_INPUTS, 101325, 1);
2048 
2049  // Do a dummy calculation at the dewpoint state - make sure all the residuals are zero as they should be
2050  CoolProp::SaturationSolvers::PTflash_twophase_options o;
2051  o.x = AS->mole_fractions_liquid();
2052  o.y = AS->mole_fractions_vapor();
2053  o.z = AS->get_mole_fractions();
2054  o.rhomolar_liq = AS->saturated_liquid_keyed_output(CoolProp::iDmolar);
2055  o.rhomolar_vap = AS->saturated_vapor_keyed_output(CoolProp::iDmolar);
2056  o.T = AS->T();
2057  o.p = AS->p();
2058  o.omega = 1.0;
2059  CoolProp::SaturationSolvers::PTflash_twophase solver(*static_cast<CoolProp::HelmholtzEOSMixtureBackend*>(AS.get()), o);
2060  solver.build_arrays();
2061  double err = solver.r.norm();
2062  REQUIRE(std::abs(err) < 1e-10);
2063 
2064  // Now, perturb the solution a little bit and actually do the solve
2065  std::vector<double> x0 = o.x;
2066  o.x[0] *= 1.1;
2067  o.x[1] = 1 - o.x[0];
2068  solver.solve();
2069  // Make sure we end up with the same liquid composition
2070  double diffx0 = o.x[0] - x0[0];
2071  REQUIRE(std::abs(diffx0) < 1e-10);
2072 
2073  // Now do the blind flash call with PT as inputs
2074  AS->update(CoolProp::PT_INPUTS, AS->p(), AS->T() - 2);
2075  REQUIRE(AS->phase() == CoolProp::iphase_twophase);
2076 }
2077 
2078 #endif
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1