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