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