CoolProp 8.0.0
An open-source fluid property and humid air property database
PCSAFTBackend.cpp
Go to the documentation of this file.
1#include <cmath>
2#include <utility>
3#include <vector>
4#include <string>
5#include <Eigen/Dense>
6#include <cstdlib>
7
8#include "PCSAFTBackend.h"
9#include "Backends/Helmholtz/VLERoutines.h"
10
11using std::vector;
12
13/*
14References
15----------
16* J. Gross and G. Sadowski, “Perturbed-Chain SAFT:  An Equation of State
17Based on a Perturbation Theory for Chain Molecules,” Ind. Eng. Chem.
18Res., vol. 40, no. 4, pp. 1244–1260, Feb. 2001.
19* M. Kleiner and G. Sadowski, “Modeling of Polar Systems Using PCP-SAFT: 
20An Approach to Account for Induced-Association Interactions,” J. Phys.
21Chem. C, vol. 111, no. 43, pp. 15544–15553, Nov. 2007.
22* Gross Joachim and Vrabec Jadran, “An equation‐of‐state contribution
23for polar components: Dipolar molecules,” AIChE J., vol. 52, no. 3,
24pp. 1194–1204, Feb. 2006.
25* A. J. de Villiers, C. E. Schwarz, and A. J. Burger, “Improving
26vapour–liquid-equilibria predictions for mixtures with non-associating polar
27components using sPC-SAFT extended with two dipolar terms,” Fluid Phase
28Equilibria, vol. 305, no. 2, pp. 174–184, Jun. 2011.
29* S. H. Huang and M. Radosz, “Equation of state for small, large,
30polydisperse, and associating molecules,” Ind. Eng. Chem. Res., vol. 29,
31no. 11, pp. 2284–2294, Nov. 1990.
32* S. H. Huang and M. Radosz, “Equation of state for small, large,
33polydisperse, and associating molecules: extension to fluid mixtures,”
34Ind. Eng. Chem. Res., vol. 30, no. 8, pp. 1994–2005, Aug. 1991.
35* S. H. Huang and M. Radosz, “Equation of state for small, large,
36polydisperse, and associating molecules: extension to fluid mixtures.
37[Erratum to document cited in CA115(8):79950j],” Ind. Eng. Chem. Res.,
38vol. 32, no. 4, pp. 762–762, Apr. 1993.
39* J. Gross and G. Sadowski, “Application of the Perturbed-Chain SAFT
40Equation of State to Associating Systems,” Ind. Eng. Chem. Res., vol.
4141, no. 22, pp. 5510–5515, Oct. 2002.
42* L. F. Cameretti, G. Sadowski, and J. M. Mollerup, “Modeling of Aqueous
43Electrolyte Solutions with Perturbed-Chain Statistical Associated Fluid
44Theory,” Ind. Eng. Chem. Res., vol. 44, no. 9, pp. 3355–3362, Apr. 2005.
45* L. F. Cameretti, G. Sadowski, and J. M. Mollerup, “Modeling of Aqueous
46Electrolyte Solutions with Perturbed-Chain Statistical Association Fluid
47Theory,” Ind. Eng. Chem. Res., vol. 44, no. 23, pp. 8944–8944, Nov. 2005.
48* C. Held, L. F. Cameretti, and G. Sadowski, “Modeling aqueous
49electrolyte solutions: Part 1. Fully dissociated electrolytes,” Fluid
50Phase Equilibria, vol. 270, no. 1, pp. 87–96, Aug. 2008.
51* C. Held, T. Reschke, S. Mohammad, A. Luza, and G. Sadowski, “ePC-SAFT
52revised,” Chem. Eng. Res. Des., vol. 92, no. 12, pp. 2884–2897, Dec. 2014.
53*/
54
55namespace CoolProp {
56
57PCSAFTBackend::PCSAFTBackend(const std::vector<std::string>& component_names, bool generate_SatL_and_SatV)
58 : N(component_names.size()), ion_term(false), polar_term(false), assoc_term(false), water_present(false), water_idx(0) {
59
60 components.resize(N);
61
62 for (unsigned int i = 0; i < N; ++i) {
63 components[i] = PCSAFTLibrary::get_library().get(component_names[i]);
64 // Determining which PC-SAFT terms should be used
65 if (components[i].getZ() != 0) {
66 ion_term = true;
67 }
68 if (components[i].getDipm() != 0) {
69 polar_term = true;
70 }
71 if (components[i].getVolA() != 0) {
72 assoc_term = true;
73 }
74 if (components[i].getCAS() == "7732-18-5") {
75 water_present = true;
76 water_idx = i;
77 }
78 }
79
80 // Set up association scheme
81 if (assoc_term) {
82 set_assoc_matrix();
83 }
84
85 // Set the components and associated flags
86 is_pure_or_pseudopure = (N == 1);
87
88 // loading interaction parameters
89 std::string kij_string;
90 std::string kijT_string;
91 if (is_pure_or_pseudopure) {
92 this->mole_fractions = std::vector<CoolPropDbl>(1, 1);
93 } else {
94 k_ij.resize(N * N, 0.0);
95 k_ijT.resize(N * N, 0.0);
96 for (unsigned int i = 0; i < N; ++i) {
97 for (unsigned int j = 0; j < N; ++j) {
98 if (i != j) {
99 kij_string = PCSAFTLibrary::get_library().get_binary_interaction_pcsaft(components[i].getCAS(), components[j].getCAS(), "kij");
100 kijT_string = PCSAFTLibrary::get_library().get_binary_interaction_pcsaft(components[i].getCAS(), components[j].getCAS(), "kijT");
101 k_ij[i * N + j] = atof(kij_string.c_str());
102 k_ijT[i * N + j] = atof(kijT_string.c_str());
103 }
104 }
105 }
106 }
107
108 if (generate_SatL_and_SatV) {
109 bool SatLSatV = false;
110 // Explicit scope: this runs inside the ctor; derived overrides
111 // of get_copy aren't active yet, so dispatch would resolve to
112 // PCSAFTBackend::get_copy regardless. Make that explicit and
113 // silence clang-analyzer-optin.cplusplus.VirtualCall.
114 SatL.reset(PCSAFTBackend::get_copy(SatLSatV));
115 SatL->specify_phase(iphase_liquid);
116 SatV.reset(PCSAFTBackend::get_copy(SatLSatV));
117 SatV->specify_phase(iphase_gas);
118 }
119
120 // Set the phase to default unknown value
121 imposed_phase_index = iphase_not_imposed;
122 _phase = iphase_unknown;
123}
124
125PCSAFTBackend::PCSAFTBackend(const std::vector<PCSAFTFluid>& components_in, bool generate_SatL_and_SatV)
126 : components(components_in), N(components.size()), ion_term(false), polar_term(false), assoc_term(false), water_present(false), water_idx(0) {
127
128 // Determining which PC-SAFT terms should be used
129
130 for (unsigned int i = 0; i < N; ++i) {
131 if (components[i].getZ() != 0) {
132 ion_term = true;
133 }
134 if (components[i].getDipm() != 0) {
135 polar_term = true;
136 }
137 if (components[i].getVolA() != 0) {
138 assoc_term = true;
139 }
140 if (components[i].getCAS() == "7732-18-5") {
141 water_present = true;
142 water_idx = i;
143 }
144 }
145
146 // Set up association scheme
147 if (assoc_term) {
148 set_assoc_matrix();
149 }
150
151 // Set the components and associated flags
152 is_pure_or_pseudopure = (N == 1);
153
154 // loading interaction parameters
155 std::string kij_string;
156 std::string kijT_string;
157 if (is_pure_or_pseudopure) {
158 this->mole_fractions = std::vector<CoolPropDbl>(1, 1);
159 } else {
160 k_ij.resize(N * N, 0.0);
161 k_ijT.resize(N * N, 0.0);
162 for (unsigned int i = 0; i < N; ++i) {
163 for (unsigned int j = 0; j < N; ++j) {
164 if (i != j) {
165 kij_string = PCSAFTLibrary::get_library().get_binary_interaction_pcsaft(components[i].getCAS(), components[j].getCAS(), "kij");
166 kijT_string = PCSAFTLibrary::get_library().get_binary_interaction_pcsaft(components[i].getCAS(), components[j].getCAS(), "kijT");
167 k_ij[i * N + j] = atof(kij_string.c_str());
168 k_ijT[i * N + j] = atof(kijT_string.c_str());
169 }
170 }
171 }
172 }
173
174 if (generate_SatL_and_SatV) {
175 bool SatLSatV = false;
176 // Explicit scope: this runs inside the ctor; derived overrides
177 // of get_copy aren't active yet, so dispatch would resolve to
178 // PCSAFTBackend::get_copy regardless. Make that explicit and
179 // silence clang-analyzer-optin.cplusplus.VirtualCall.
180 SatL.reset(PCSAFTBackend::get_copy(SatLSatV));
181 SatL->specify_phase(iphase_liquid);
182 SatV.reset(PCSAFTBackend::get_copy(SatLSatV));
183 SatV->specify_phase(iphase_gas);
184 }
185
186 // Set the phase to default unknown value
187 imposed_phase_index = iphase_not_imposed;
188 _phase = iphase_unknown;
189}
190
191void PCSAFTBackend::set_mole_fractions(const std::vector<CoolPropDbl>& mole_fractions) {
192 if (mole_fractions.size() != N) {
193 throw ValueError(format("size of mole fraction vector [%d] does not equal that of component vector [%d]", mole_fractions.size(), N));
194 }
195 // Copy values without reallocating memory
196 this->mole_fractions = mole_fractions; // Most effective copy
197 this->resize(N); // No reallocation of this->mole_fractions happens
198 // Also store the mole fractions as doubles
199 this->mole_fractions_double = std::vector<double>(mole_fractions.begin(), mole_fractions.end());
200};
201
202void PCSAFTBackend::set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions) {
203 if (mass_fractions.size() != N) {
204 throw ValueError(format("size of mass fraction vector [%d] does not equal that of component vector [%d]", mass_fractions.size(), N));
205 }
206 std::vector<CoolPropDbl> moles;
207 CoolPropDbl sum_moles = 0.0;
208 CoolPropDbl tmp = 0.0;
209 for (unsigned int i = 0; i < components.size(); ++i) {
210 tmp = mass_fractions[i] / components[i].molar_mass();
211 moles.push_back(tmp);
212 sum_moles += tmp;
213 }
214 std::vector<CoolPropDbl> mole_fractions;
215 for (const auto& mole : moles) {
216 mole_fractions.push_back(mole / sum_moles);
217 }
218 this->set_mole_fractions(mole_fractions);
219};
220
221PCSAFTBackend* PCSAFTBackend::get_copy(bool generate_SatL_and_SatV) {
222 // Set up the class with these components
223 auto* ptr = new PCSAFTBackend(components, generate_SatL_and_SatV);
224 return ptr;
225};
226
227void PCSAFTBackend::resize(std::size_t N) {
228 this->mole_fractions.resize(N);
229 this->mole_fractions_double.resize(N);
230 this->K.resize(N);
231 this->lnK.resize(N);
232}
233
234CoolPropDbl PCSAFTBackend::update_DmolarT(CoolPropDbl rho) {
235 _rhomolar = rho;
236 return this->calc_pressure();
237}
238
239CoolPropDbl PCSAFTBackend::calc_pressure() {
240 double den = _rhomolar * N_AV / 1.0e30;
241
242 CoolPropDbl Z = this->calc_compressibility_factor();
243 CoolPropDbl P = Z * kb * _T * den * 1.0e30; // Pa
244 return P;
245}
246
247CoolPropDbl PCSAFTBackend::calc_alphar() {
248 auto ncomp = N; // number of components
249 vector<double> d(ncomp);
250 for (auto i = 0U; i < ncomp; i++) {
251 d[i] = components[i].getSigma() * (1 - 0.12 * exp(-3 * components[i].getU() / _T));
252 }
253 if (ion_term) {
254 for (auto i = 0U; i < ncomp; i++) {
255 if (components[i].getZ() != 0) {
256 d[i] =
257 components[i].getSigma() * (1 - 0.12); // for ions the diameter is assumed to be temperature independent (see Held et al. 2014)
258 }
259 }
260 }
261
262 double den = _rhomolar * N_AV / 1.0e30;
263
264 vector<double> zeta(4, 0);
265 double summ = NAN;
266 for (int i = 0; i < 4; i++) {
267 summ = 0;
268 for (size_t j = 0; j < ncomp; j++) {
269 summ += mole_fractions[j] * components[j].getM() * pow(d[j], i);
270 }
271 zeta[i] = PI / 6 * den * summ;
272 }
273
274 double eta = zeta[3];
275 double m_avg = 0;
276 for (size_t i = 0; i < ncomp; i++) {
277 m_avg += mole_fractions[i] * components[i].getM();
278 }
279
280 vector<double> ghs(ncomp * ncomp, 0);
281 vector<double> e_ij(ncomp * ncomp, 0);
282 vector<double> s_ij(ncomp * ncomp, 0);
283 double m2es3 = 0.;
284 double m2e2s3 = 0.;
285 int idx = -1;
286 for (size_t i = 0; i < ncomp; i++) {
287 for (size_t j = 0; j < ncomp; j++) {
288 idx += 1;
289 s_ij[idx] = (components[i].getSigma() + components[j].getSigma()) / 2.;
290 if (ion_term) {
291 if (components[i].getZ() * components[j].getZ()
292 <= 0) { // for two cations or two anions e_ij is kept at zero to avoid dispersion between like ions (see Held et al. 2014)
293 if (k_ij.empty()) {
294 e_ij[idx] = sqrt(components[i].getU() * components[j].getU());
295 } else {
296 e_ij[idx] = sqrt(components[i].getU() * components[j].getU()) * (1 - (k_ij[idx] + k_ijT[idx] * _T));
297 }
298 }
299 } else {
300 if (k_ij.empty()) {
301 e_ij[idx] = sqrt(components[i].getU() * components[j].getU());
302 } else {
303 e_ij[idx] = sqrt(components[i].getU() * components[j].getU()) * (1 - (k_ij[idx] + k_ijT[idx] * _T));
304 }
305 }
306 m2es3 = m2es3 + mole_fractions[i] * mole_fractions[j] * components[i].getM() * components[j].getM() * e_ij[idx] / _T * pow(s_ij[idx], 3);
307 m2e2s3 =
308 m2e2s3
309 + mole_fractions[i] * mole_fractions[j] * components[i].getM() * components[j].getM() * pow(e_ij[idx] / _T, 2) * pow(s_ij[idx], 3);
310 ghs[idx] = 1 / (1 - zeta[3]) + (d[i] * d[j] / (d[i] + d[j])) * 3 * zeta[2] / (1 - zeta[3]) / (1 - zeta[3])
311 + pow(d[i] * d[j] / (d[i] + d[j]), 2) * 2 * zeta[2] * zeta[2] / pow(1 - zeta[3], 3);
312 }
313 }
314
315 double ares_hs = 1 / zeta[0]
316 * (3 * zeta[1] * zeta[2] / (1 - zeta[3]) + pow(zeta[2], 3.) / (zeta[3] * pow(1 - zeta[3], 2))
317 + (pow(zeta[2], 3.) / pow(zeta[3], 2.) - zeta[0]) * log(1 - zeta[3]));
318
319 static double a0[7] = {0.9105631445, 0.6361281449, 2.6861347891, -26.547362491, 97.759208784, -159.59154087, 91.297774084};
320 static double a1[7] = {-0.3084016918, 0.1860531159, -2.5030047259, 21.419793629, -65.255885330, 83.318680481, -33.746922930};
321 static double a2[7] = {-0.0906148351, 0.4527842806, 0.5962700728, -1.7241829131, -4.1302112531, 13.776631870, -8.6728470368};
322 static double b0[7] = {0.7240946941, 2.2382791861, -4.0025849485, -21.003576815, 26.855641363, 206.55133841, -355.60235612};
323 static double b1[7] = {-0.5755498075, 0.6995095521, 3.8925673390, -17.215471648, 192.67226447, -161.82646165, -165.20769346};
324 static double b2[7] = {0.0976883116, -0.2557574982, -9.1558561530, 20.642075974, -38.804430052, 93.626774077, -29.666905585};
325
326 vector<double> a(7, 0);
327 vector<double> b(7, 0);
328 for (int i = 0; i < 7; i++) {
329 a[i] = a0[i] + (m_avg - 1.) / m_avg * a1[i] + (m_avg - 1.) / m_avg * (m_avg - 2.) / m_avg * a2[i];
330 b[i] = b0[i] + (m_avg - 1.) / m_avg * b1[i] + (m_avg - 1.) / m_avg * (m_avg - 2.) / m_avg * b2[i];
331 }
332
333 double I1 = 0.0;
334 double I2 = 0.0;
335 for (int i = 0; i < 7; i++) {
336 I1 += a[i] * pow(eta, i);
337 I2 += b[i] * pow(eta, i);
338 }
339 double C1 = 1.
340 / (1. + m_avg * (8 * eta - 2 * eta * eta) / pow(1 - eta, 4)
341 + (1 - m_avg) * (20 * eta - 27 * eta * eta + 12 * pow(eta, 3) - 2 * pow(eta, 4)) / pow((1 - eta) * (2 - eta), 2.0));
342
343 summ = 0.0;
344 for (size_t i = 0; i < ncomp; i++) {
345 summ += mole_fractions[i] * (components[i].getM() - 1) * log(ghs[i * ncomp + i]);
346 }
347
348 double ares_hc = m_avg * ares_hs - summ;
349 double ares_disp = -2 * PI * den * I1 * m2es3 - PI * den * m_avg * C1 * I2 * m2e2s3;
350
351 // Dipole term (Gross and Vrabec term) --------------------------------------
352 double ares_polar = 0.;
353 if (polar_term) {
354 double A2 = 0.;
355 double A3 = 0.;
356 vector<double> dipmSQ(ncomp, 0);
357
358 static double a0dip[5] = {0.3043504, -0.1358588, 1.4493329, 0.3556977, -2.0653308};
359 static double a1dip[5] = {0.9534641, -1.8396383, 2.0131180, -7.3724958, 8.2374135};
360 static double a2dip[5] = {-1.1610080, 4.5258607, 0.9751222, -12.281038, 5.9397575};
361 static double b0dip[5] = {0.2187939, -1.1896431, 1.1626889, 0, 0};
362 static double b1dip[5] = {-0.5873164, 1.2489132, -0.5085280, 0, 0};
363 static double b2dip[5] = {3.4869576, -14.915974, 15.372022, 0, 0};
364 static double c0dip[5] = {-0.0646774, 0.1975882, -0.8087562, 0.6902849, 0};
365 static double c1dip[5] = {-0.9520876, 2.9924258, -2.3802636, -0.2701261, 0};
366 static double c2dip[5] = {-0.6260979, 1.2924686, 1.6542783, -3.4396744, 0};
367
368 const static double conv = 7242.702976750923; // conversion factor, see the note below Table 2 in Gross and Vrabec 2006
369
370 for (size_t i = 0; i < ncomp; i++) {
371 dipmSQ[i] = pow(components[i].getDipm(), 2.) / (components[i].getM() * components[i].getU() * pow(components[i].getSigma(), 3.)) * conv;
372 }
373
374 vector<double> adip(5, 0);
375 vector<double> bdip(5, 0);
376 vector<double> cdip(5, 0);
377 double J2 = NAN, J3 = NAN;
378 double m_ij = NAN;
379 double m_ijk = NAN;
380 for (size_t i = 0; i < ncomp; i++) {
381 for (size_t j = 0; j < ncomp; j++) {
382 m_ij = sqrt(components[i].getM() * components[j].getM());
383 if (m_ij > 2) {
384 m_ij = 2;
385 }
386 J2 = 0.;
387 for (int l = 0; l < 5; l++) {
388 adip[l] = a0dip[l] + (m_ij - 1) / m_ij * a1dip[l] + (m_ij - 1) / m_ij * (m_ij - 2) / m_ij * a2dip[l];
389 bdip[l] = b0dip[l] + (m_ij - 1) / m_ij * b1dip[l] + (m_ij - 1) / m_ij * (m_ij - 2) / m_ij * b2dip[l];
390 J2 += (adip[l] + bdip[l] * e_ij[i * ncomp + j] / _T)
391 * pow(eta, l); // i*ncomp+j needs to be used for e_ij because it is formatted as a 1D vector
392 }
393 A2 += mole_fractions[i] * mole_fractions[j] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * pow(s_ij[i * ncomp + i], 3)
394 * pow(s_ij[j * ncomp + j], 3) / pow(s_ij[i * ncomp + j], 3) * components[i].getDipnum() * components[j].getDipnum() * dipmSQ[i]
395 * dipmSQ[j] * J2;
396
397 for (size_t k = 0; k < ncomp; k++) {
398 m_ijk = pow((components[i].getM() * components[j].getM() * components[k].getM()), 1 / 3.);
399 if (m_ijk > 2) {
400 m_ijk = 2;
401 }
402 J3 = 0.;
403 for (int l = 0; l < 5; l++) {
404 cdip[l] = c0dip[l] + (m_ijk - 1) / m_ijk * c1dip[l] + (m_ijk - 1) / m_ijk * (m_ijk - 2) / m_ijk * c2dip[l];
405 J3 += cdip[l] * pow(eta, l);
406 }
407 A3 += mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T
408 * e_ij[k * ncomp + k] / _T * pow(s_ij[i * ncomp + i], 3) * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3)
409 / s_ij[i * ncomp + j] / s_ij[i * ncomp + k] / s_ij[j * ncomp + k] * components[i].getDipnum() * components[j].getDipnum()
410 * components[k].getDipnum() * dipmSQ[i] * dipmSQ[j] * dipmSQ[k] * J3;
411 }
412 }
413 }
414
415 A2 = -PI * den * A2;
416 A3 = -4 / 3. * PI * PI * den * den * A3;
417
418 if (A2 != 0) { // when the mole fraction of the polar compounds is 0 then A2 = 0 and division by 0 occurs
419 ares_polar = A2 / (1 - A3 / A2);
420 }
421 }
422
423 // Association term -------------------------------------------------------
424 double ares_assoc = 0.;
425 if (assoc_term) {
426 int num_sites = 0;
427 vector<int> iA; //indices of associating compounds
428 for (auto it = assoc_num.begin(); it != assoc_num.end(); ++it) {
429 num_sites += *it;
430 for (int i = 0; i < *it; i++) {
431 iA.push_back(static_cast<int>(it - assoc_num.begin()));
432 }
433 }
434
435 vector<double> x_assoc(num_sites); // mole fractions of only the associating compounds
436 for (int i = 0; i < num_sites; i++) {
437 x_assoc[i] = mole_fractions[iA[i]];
438 }
439
440 // these indices are necessary because we are only using 1D vectors
441 vector<double> XA(num_sites, 0);
442 vector<double> delta_ij(static_cast<std::size_t>(num_sites) * num_sites, 0);
443 auto idxa = 0ULL;
444 auto idxi = 0ULL; // index for the ii-th compound
445 auto idxj = 0ULL; // index for the jj-th compound
446 for (auto i = 0; i < num_sites; i++) {
447 idxi = iA[i] * ncomp + iA[i];
448 for (int j = 0; j < num_sites; j++) {
449 idxj = iA[j] * ncomp + iA[j];
450 if (assoc_matrix[idxa] != 0) {
451 double eABij = (components[iA[i]].getUAB() + components[iA[j]].getUAB()) / 2.;
452 double volABij = sqrt(components[iA[i]].getVolA() * components[iA[j]].getVolA())
453 * pow(sqrt(s_ij[idxi] * s_ij[idxj]) / (0.5 * (s_ij[idxi] + s_ij[idxj])), 3);
454 delta_ij[idxa] = ghs[iA[i] * ncomp + iA[j]] * (exp(eABij / _T) - 1) * pow(s_ij[iA[i] * ncomp + iA[j]], 3) * volABij;
455 }
456 idxa += 1;
457 }
458 XA[i] = (-1 + sqrt(1 + 8 * den * delta_ij[i * num_sites + i])) / (4 * den * delta_ij[i * num_sites + i]);
459 if (!std::isfinite(XA[i])) {
460 XA[i] = 0.02;
461 }
462 }
463
464 int ctr = 0;
465 double dif = 1000.;
466 vector<double> XA_old = XA;
467 while ((ctr < 100) && (dif > 1e-15)) {
468 ctr += 1;
469 XA = XA_find(XA_old, delta_ij, den, x_assoc);
470 dif = 0.;
471 for (int i = 0; i < num_sites; i++) {
472 dif += abs(XA[i] - XA_old[i]);
473 }
474 for (int i = 0; i < num_sites; i++) {
475 XA_old[i] = (XA[i] + XA_old[i]) / 2.0;
476 }
477 }
478
479 ares_assoc = 0.;
480 for (int i = 0; i < num_sites; i++) {
481 ares_assoc += mole_fractions[iA[i]] * (log(XA[i]) - 0.5 * XA[i] + 0.5);
482 }
483 }
484
485 // Ion term ---------------------------------------------------------------
486 double ares_ion = 0.;
487 if (ion_term) {
488 vector<double> q(ncomp);
489 for (size_t i = 0; i < ncomp; i++) {
490 q[i] = components[i].getZ() * E_CHRG;
491 }
492
493 summ = 0.;
494 for (size_t i = 0; i < ncomp; i++) {
495 summ += components[i].getZ() * components[i].getZ() * mole_fractions[i];
496 }
497 double kappa =
498 sqrt(den * E_CHRG * E_CHRG / kb / _T / (dielc * perm_vac) * summ); // the inverse Debye screening length. Equation 4 in Held et al. 2008.
499
500 if (kappa != 0) {
501 vector<double> chi(ncomp);
502 vector<double> sigma_k(ncomp);
503 summ = 0.;
504 for (size_t i = 0; i < ncomp; i++) {
505 chi[i] = 3 / pow(kappa * d[i], 3) * (1.5 + log(1 + kappa * d[i]) - 2 * (1 + kappa * d[i]) + 0.5 * pow(1 + kappa * d[i], 2));
506 summ += mole_fractions[i] * q[i] * q[i] * chi[i] * kappa;
507 }
508
509 ares_ion = -1 / 12. / PI / kb / _T / (dielc * perm_vac) * summ;
510 }
511 }
512
513 CoolPropDbl ares = ares_hc + ares_disp + ares_polar + ares_assoc + ares_ion;
514 return ares;
515}
516
517CoolPropDbl PCSAFTBackend::calc_dadt() {
518 auto ncomp = N; // number of components
519 vector<double> d(ncomp), dd_dt(ncomp);
520 for (auto i = 0U; i < ncomp; i++) {
521 d[i] = components[i].getSigma() * (1 - 0.12 * exp(-3 * components[i].getU() / _T));
522 dd_dt[i] = components[i].getSigma() * -3 * components[i].getU() / _T / _T * 0.12 * exp(-3 * components[i].getU() / _T);
523 }
524 if (ion_term) {
525 for (auto i = 0U; i < ncomp; i++) {
526 if (components[i].getZ() != 0) {
527 d[i] =
528 components[i].getSigma() * (1 - 0.12); // for ions the diameter is assumed to be temperature independent (see Held et al. 2014)
529 dd_dt[i] = 0.;
530 }
531 }
532 }
533
534 double den = _rhomolar * N_AV / 1.0e30;
535
536 vector<double> zeta(4, 0);
537 double summ = NAN;
538 for (int i = 0; i < 4; i++) {
539 summ = 0;
540 for (size_t j = 0; j < ncomp; j++) {
541 summ += mole_fractions[j] * components[j].getM() * pow(d[j], i);
542 }
543 zeta[i] = PI / 6 * den * summ;
544 }
545
546 vector<double> dzeta_dt(4, 0);
547 for (int i = 1; i < 4; i++) {
548 summ = 0;
549 for (size_t j = 0; j < ncomp; j++) {
550 summ += mole_fractions[j] * components[j].getM() * i * dd_dt[j] * pow(d[j], (i - 1));
551 }
552 dzeta_dt[i] = PI / 6 * den * summ;
553 }
554
555 double eta = zeta[3];
556 double m_avg = 0;
557 for (size_t i = 0; i < ncomp; i++) {
558 m_avg += mole_fractions[i] * components[i].getM();
559 }
560
561 vector<double> ghs(ncomp * ncomp, 0);
562 vector<double> dghs_dt(ncomp * ncomp, 0);
563 vector<double> e_ij(ncomp * ncomp, 0);
564 vector<double> s_ij(ncomp * ncomp, 0);
565 double m2es3 = 0.;
566 double m2e2s3 = 0.;
567 double ddij_dt = NAN;
568 int idx = -1;
569 for (size_t i = 0; i < ncomp; i++) {
570 for (size_t j = 0; j < ncomp; j++) {
571 idx += 1;
572 s_ij[idx] = (components[i].getSigma() + components[j].getSigma()) / 2.;
573 if (ion_term) {
574 if (components[i].getZ() * components[j].getZ()
575 <= 0) { // for two cations or two anions e_ij is kept at zero to avoid dispersion between like ions (see Held et al. 2014)
576 if (k_ij.empty()) {
577 e_ij[idx] = sqrt(components[i].getU() * components[j].getU());
578 } else {
579 e_ij[idx] = sqrt(components[i].getU() * components[j].getU()) * (1 - (k_ij[idx] + k_ijT[idx] * _T));
580 }
581 }
582 } else {
583 if (k_ij.empty()) {
584 e_ij[idx] = sqrt(components[i].getU() * components[j].getU());
585 } else {
586 e_ij[idx] = sqrt(components[i].getU() * components[j].getU()) * (1 - (k_ij[idx] + k_ijT[idx] * _T));
587 }
588 }
589 m2es3 = m2es3 + mole_fractions[i] * mole_fractions[j] * components[i].getM() * components[j].getM() * e_ij[idx] / _T * pow(s_ij[idx], 3);
590 m2e2s3 =
591 m2e2s3
592 + mole_fractions[i] * mole_fractions[j] * components[i].getM() * components[j].getM() * pow(e_ij[idx] / _T, 2) * pow(s_ij[idx], 3);
593 ghs[idx] = 1 / (1 - zeta[3]) + (d[i] * d[j] / (d[i] + d[j])) * 3 * zeta[2] / (1 - zeta[3]) / (1 - zeta[3])
594 + pow(d[i] * d[j] / (d[i] + d[j]), 2) * 2 * zeta[2] * zeta[2] / pow(1 - zeta[3], 3);
595 ddij_dt = (d[i] * d[j] / (d[i] + d[j])) * (dd_dt[i] / d[i] + dd_dt[j] / d[j] - (dd_dt[i] + dd_dt[j]) / (d[i] + d[j]));
596 dghs_dt[idx] = dzeta_dt[3] / pow(1 - zeta[3], 2.)
597 + 3 * (ddij_dt * zeta[2] + (d[i] * d[j] / (d[i] + d[j])) * dzeta_dt[2]) / pow(1 - zeta[3], 2.)
598 + 4 * (d[i] * d[j] / (d[i] + d[j])) * zeta[2]
599 * (1.5 * dzeta_dt[3] + ddij_dt * zeta[2] + (d[i] * d[j] / (d[i] + d[j])) * dzeta_dt[2]) / pow(1 - zeta[3], 3.)
600 + 6 * pow((d[i] * d[j] / (d[i] + d[j])) * zeta[2], 2.) * dzeta_dt[3] / pow(1 - zeta[3], 4.);
601 }
602 }
603
604 double dadt_hs =
605 1 / zeta[0]
606 * (3 * (dzeta_dt[1] * zeta[2] + zeta[1] * dzeta_dt[2]) / (1 - zeta[3]) + 3 * zeta[1] * zeta[2] * dzeta_dt[3] / pow(1 - zeta[3], 2.)
607 + 3 * pow(zeta[2], 2.) * dzeta_dt[2] / zeta[3] / pow(1 - zeta[3], 2.)
608 + pow(zeta[2], 3.) * dzeta_dt[3] * (3 * zeta[3] - 1) / pow(zeta[3], 2.) / pow(1 - zeta[3], 3.)
609 + (3 * pow(zeta[2], 2.) * dzeta_dt[2] * zeta[3] - 2 * pow(zeta[2], 3.) * dzeta_dt[3]) / pow(zeta[3], 3.) * log(1 - zeta[3])
610 + (zeta[0] - pow(zeta[2], 3) / pow(zeta[3], 2.)) * dzeta_dt[3] / (1 - zeta[3]));
611
612 static double a0[7] = {0.9105631445, 0.6361281449, 2.6861347891, -26.547362491, 97.759208784, -159.59154087, 91.297774084};
613 static double a1[7] = {-0.3084016918, 0.1860531159, -2.5030047259, 21.419793629, -65.255885330, 83.318680481, -33.746922930};
614 static double a2[7] = {-0.0906148351, 0.4527842806, 0.5962700728, -1.7241829131, -4.1302112531, 13.776631870, -8.6728470368};
615 static double b0[7] = {0.7240946941, 2.2382791861, -4.0025849485, -21.003576815, 26.855641363, 206.55133841, -355.60235612};
616 static double b1[7] = {-0.5755498075, 0.6995095521, 3.8925673390, -17.215471648, 192.67226447, -161.82646165, -165.20769346};
617 static double b2[7] = {0.0976883116, -0.2557574982, -9.1558561530, 20.642075974, -38.804430052, 93.626774077, -29.666905585};
618
619 vector<double> a(7, 0);
620 vector<double> b(7, 0);
621 for (int i = 0; i < 7; i++) {
622 a[i] = a0[i] + (m_avg - 1.) / m_avg * a1[i] + (m_avg - 1.) / m_avg * (m_avg - 2.) / m_avg * a2[i];
623 b[i] = b0[i] + (m_avg - 1.) / m_avg * b1[i] + (m_avg - 1.) / m_avg * (m_avg - 2.) / m_avg * b2[i];
624 }
625
626 double I1 = 0.0;
627 double I2 = 0.0;
628 double dI1_dt = 0.0, dI2_dt = 0.;
629 for (int i = 0; i < 7; i++) {
630 I1 += a[i] * pow(eta, i);
631 I2 += b[i] * pow(eta, i);
632 dI1_dt += a[i] * dzeta_dt[3] * i * pow(eta, i - 1);
633 dI2_dt += b[i] * dzeta_dt[3] * i * pow(eta, i - 1);
634 }
635 double C1 = 1.
636 / (1. + m_avg * (8 * eta - 2 * eta * eta) / pow(1 - eta, 4)
637 + (1 - m_avg) * (20 * eta - 27 * eta * eta + 12 * pow(eta, 3) - 2 * pow(eta, 4)) / pow((1 - eta) * (2 - eta), 2.0));
638 double C2 = -1 * C1 * C1
639 * (m_avg * (-4 * eta * eta + 20 * eta + 8) / pow(1 - eta, 5.)
640 + (1 - m_avg) * (2 * pow(eta, 3) + 12 * eta * eta - 48 * eta + 40) / pow((1 - eta) * (2 - eta), 3));
641 double dC1_dt = C2 * dzeta_dt[3];
642
643 summ = 0.;
644 for (size_t i = 0; i < ncomp; i++) {
645 summ += mole_fractions[i] * (components[i].getM() - 1) * dghs_dt[i * ncomp + i] / ghs[i * ncomp + i];
646 }
647
648 double dadt_hc = m_avg * dadt_hs - summ;
649 double dadt_disp = -2 * PI * den * (dI1_dt - I1 / _T) * m2es3 - PI * den * m_avg * (dC1_dt * I2 + C1 * dI2_dt - 2 * C1 * I2 / _T) * m2e2s3;
650
651 // Dipole term (Gross and Vrabec term) --------------------------------------
652 double dadt_polar = 0.;
653 if (polar_term) {
654 double A2 = 0.;
655 double A3 = 0.;
656 double dA2_dt = 0.;
657 double dA3_dt = 0.;
658 vector<double> dipmSQ(ncomp, 0);
659
660 static double a0dip[5] = {0.3043504, -0.1358588, 1.4493329, 0.3556977, -2.0653308};
661 static double a1dip[5] = {0.9534641, -1.8396383, 2.0131180, -7.3724958, 8.2374135};
662 static double a2dip[5] = {-1.1610080, 4.5258607, 0.9751222, -12.281038, 5.9397575};
663 static double b0dip[5] = {0.2187939, -1.1896431, 1.1626889, 0, 0};
664 static double b1dip[5] = {-0.5873164, 1.2489132, -0.5085280, 0, 0};
665 static double b2dip[5] = {3.4869576, -14.915974, 15.372022, 0, 0};
666 static double c0dip[5] = {-0.0646774, 0.1975882, -0.8087562, 0.6902849, 0};
667 static double c1dip[5] = {-0.9520876, 2.9924258, -2.3802636, -0.2701261, 0};
668 static double c2dip[5] = {-0.6260979, 1.2924686, 1.6542783, -3.4396744, 0};
669
670 const static double conv = 7242.702976750923; // conversion factor, see the note below Table 2 in Gross and Vrabec 2006
671
672 for (size_t i = 0; i < ncomp; i++) {
673 dipmSQ[i] = pow(components[i].getDipm(), 2.) / (components[i].getM() * components[i].getU() * pow(components[i].getSigma(), 3.)) * conv;
674 }
675
676 vector<double> adip(5, 0);
677 vector<double> bdip(5, 0);
678 vector<double> cdip(5, 0);
679 double J2 = NAN, J3 = NAN, dJ2_dt = NAN, dJ3_dt = NAN;
680 double m_ij = NAN;
681 double m_ijk = NAN;
682 for (size_t i = 0; i < ncomp; i++) {
683 for (size_t j = 0; j < ncomp; j++) {
684 m_ij = sqrt(components[i].getM() * components[j].getM());
685 if (m_ij > 2) {
686 m_ij = 2;
687 }
688 J2 = 0.;
689 dJ2_dt = 0.;
690 for (int l = 0; l < 5; l++) {
691 adip[l] = a0dip[l] + (m_ij - 1) / m_ij * a1dip[l] + (m_ij - 1) / m_ij * (m_ij - 2) / m_ij * a2dip[l];
692 bdip[l] = b0dip[l] + (m_ij - 1) / m_ij * b1dip[l] + (m_ij - 1) / m_ij * (m_ij - 2) / m_ij * b2dip[l];
693 J2 += (adip[l] + bdip[l] * e_ij[i * ncomp + j] / _T)
694 * pow(eta, l); // i*ncomp+j needs to be used for e_ij because it is formatted as a 1D vector
695 dJ2_dt += adip[l] * l * pow(eta, l - 1) * dzeta_dt[3]
696 + bdip[l] * e_ij[j * ncomp + j] * (1 / _T * l * pow(eta, l - 1) * dzeta_dt[3] - 1 / pow(_T, 2.) * pow(eta, l));
697 }
698 A2 += mole_fractions[i] * mole_fractions[j] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * pow(s_ij[i * ncomp + i], 3)
699 * pow(s_ij[j * ncomp + j], 3) / pow(s_ij[i * ncomp + j], 3) * components[i].getDipnum() * components[j].getDipnum() * dipmSQ[i]
700 * dipmSQ[j] * J2;
701 dA2_dt += mole_fractions[i] * mole_fractions[j] * e_ij[i * ncomp + i] * e_ij[j * ncomp + j] * pow(s_ij[i * ncomp + i], 3)
702 * pow(s_ij[j * ncomp + j], 3) / pow(s_ij[i * ncomp + j], 3) * components[i].getDipnum() * components[j].getDipnum()
703 * dipmSQ[i] * dipmSQ[j] * (dJ2_dt / pow(_T, 2) - 2 * J2 / pow(_T, 3));
704
705 for (size_t k = 0; k < ncomp; k++) {
706 m_ijk = pow((components[i].getM() * components[j].getM() * components[k].getM()), 1 / 3.);
707 if (m_ijk > 2) {
708 m_ijk = 2;
709 }
710 J3 = 0.;
711 dJ3_dt = 0.;
712 for (int l = 0; l < 5; l++) {
713 cdip[l] = c0dip[l] + (m_ijk - 1) / m_ijk * c1dip[l] + (m_ijk - 1) / m_ijk * (m_ijk - 2) / m_ijk * c2dip[l];
714 J3 += cdip[l] * pow(eta, l);
715 dJ3_dt += cdip[l] * l * pow(eta, l - 1) * dzeta_dt[3];
716 }
717 A3 += mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T
718 * e_ij[k * ncomp + k] / _T * pow(s_ij[i * ncomp + i], 3) * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3)
719 / s_ij[i * ncomp + j] / s_ij[i * ncomp + k] / s_ij[j * ncomp + k] * components[i].getDipnum() * components[j].getDipnum()
720 * components[k].getDipnum() * dipmSQ[i] * dipmSQ[j] * dipmSQ[k] * J3;
721 dA3_dt += mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * e_ij[i * ncomp + i] * e_ij[j * ncomp + j]
722 * e_ij[k * ncomp + k] * pow(s_ij[i * ncomp + i], 3) * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3)
723 / s_ij[i * ncomp + j] / s_ij[i * ncomp + k] / s_ij[j * ncomp + k] * components[i].getDipnum()
724 * components[j].getDipnum() * components[k].getDipnum() * dipmSQ[i] * dipmSQ[j] * dipmSQ[k]
725 * (-3 * J3 / pow(_T, 4) + dJ3_dt / pow(_T, 3));
726 }
727 }
728 }
729
730 A2 = -PI * den * A2;
731 A3 = -4 / 3. * PI * PI * den * den * A3;
732 dA2_dt = -PI * den * dA2_dt;
733 dA3_dt = -4 / 3. * PI * PI * den * den * dA3_dt;
734
735 if (A2 != 0) { // when the mole fraction of the polar compounds is 0 then A2 = 0 and division by 0 occurs
736 dadt_polar = (dA2_dt - 2 * A3 / A2 * dA2_dt + dA3_dt) / pow(1 - A3 / A2, 2.);
737 }
738 }
739
740 // Association term -------------------------------------------------------
741 double dadt_assoc = 0.;
742 if (assoc_term) {
743 int num_sites = 0;
744 vector<int> iA; //indices of associating compounds
745 for (auto it = assoc_num.begin(); it != assoc_num.end(); ++it) {
746 num_sites += *it;
747 for (int i = 0; i < *it; i++) {
748 iA.push_back(static_cast<int>(it - assoc_num.begin()));
749 }
750 }
751
752 vector<double> x_assoc(num_sites); // mole fractions of only the associating compounds
753 for (int i = 0; i < num_sites; i++) {
754 x_assoc[i] = mole_fractions[iA[i]];
755 }
756
757 // these indices are necessary because we are only using 1D vectors
758 vector<double> XA(num_sites, 0);
759 vector<double> delta_ij(static_cast<std::size_t>(num_sites) * num_sites, 0);
760 vector<double> ddelta_dt(static_cast<std::size_t>(num_sites) * num_sites, 0);
761 std::size_t idxa = 0UL;
762 std::size_t idxi = 0UL; // index for the ii-th compound
763 std::size_t idxj = 0UL; // index for the jj-th compound
764 for (auto i = 0; i < num_sites; i++) {
765 idxi = iA[i] * ncomp + iA[i];
766 for (auto j = 0; j < num_sites; j++) {
767 idxj = iA[j] * ncomp + iA[j];
768 if (assoc_matrix[idxa] != 0) {
769 double eABij = (components[iA[i]].getUAB() + components[iA[j]].getUAB()) / 2.;
770 double volABij = sqrt(components[iA[i]].getVolA() * components[iA[j]].getVolA())
771 * pow(sqrt(s_ij[idxi] * s_ij[idxj]) / (0.5 * (s_ij[idxi] + s_ij[idxj])), 3);
772 delta_ij[idxa] = ghs[iA[i] * ncomp + iA[j]] * (exp(eABij / _T) - 1) * pow(s_ij[iA[i] * ncomp + iA[j]], 3) * volABij;
773 ddelta_dt[idxa] =
774 pow(s_ij[idxj], 3) * volABij
775 * (-eABij / pow(_T, 2) * exp(eABij / _T) * ghs[iA[i] * ncomp + iA[j]] + dghs_dt[iA[i] * ncomp + iA[j]] * (exp(eABij / _T) - 1));
776 }
777 idxa += 1;
778 }
779 XA[i] = (-1 + sqrt(1 + 8 * den * delta_ij[i * num_sites + i])) / (4 * den * delta_ij[i * num_sites + i]);
780 if (!std::isfinite(XA[i])) {
781 XA[i] = 0.02;
782 }
783 }
784
785 int ctr = 0;
786 double dif = 1000.;
787 vector<double> XA_old = XA;
788 while ((ctr < 100) && (dif > 1e-15)) {
789 ctr += 1;
790 XA = XA_find(XA_old, delta_ij, den, x_assoc);
791 dif = 0.;
792 for (int i = 0; i < num_sites; i++) {
793 dif += abs(XA[i] - XA_old[i]);
794 }
795 for (int i = 0; i < num_sites; i++) {
796 XA_old[i] = (XA[i] + XA_old[i]) / 2.0;
797 }
798 }
799
800 vector<double> dXA_dt(num_sites, 0);
801 dXA_dt = dXAdt_find(delta_ij, den, XA, ddelta_dt, x_assoc);
802
803 for (int i = 0; i < num_sites; i++) {
804 dadt_assoc += mole_fractions[iA[i]] * (1 / XA[i] - 0.5) * dXA_dt[i];
805 }
806 }
807
808 // Ion term ---------------------------------------------------------------
809 double dadt_ion = 0.;
810 if (ion_term) {
811 vector<double> q(ncomp);
812 for (size_t i = 0; i < ncomp; i++) {
813 q[i] = components[i].getZ() * E_CHRG;
814 }
815
816 summ = 0.;
817 for (size_t i = 0; i < ncomp; i++) {
818 summ += components[i].getZ() * components[i].getZ() * mole_fractions[i];
819 }
820 double kappa =
821 sqrt(den * E_CHRG * E_CHRG / kb / _T / (dielc * perm_vac) * summ); // the inverse Debye screening length. Equation 4 in Held et al. 2008.
822
823 double dkappa_dt = NAN;
824 if (kappa != 0) {
825 vector<double> chi(ncomp);
826 vector<double> dchikap_dk(ncomp);
827 summ = 0.;
828 for (size_t i = 0; i < ncomp; i++) {
829 chi[i] = 3 / pow(kappa * d[i], 3) * (1.5 + log(1 + kappa * d[i]) - 2 * (1 + kappa * d[i]) + 0.5 * pow(1 + kappa * d[i], 2));
830 dchikap_dk[i] = -2 * chi[i] + 3 / (1 + kappa * d[i]);
831 summ += mole_fractions[i] * components[i].getZ() * components[i].getZ();
832 }
833 dkappa_dt = -0.5 * den * E_CHRG * E_CHRG / kb / _T / _T / (dielc * perm_vac) * summ / kappa;
834
835 summ = 0.;
836 for (size_t i = 0; i < ncomp; i++) {
837 summ += mole_fractions[i] * q[i] * q[i] * (dchikap_dk[i] * dkappa_dt / _T - kappa * chi[i] / _T / _T);
838 }
839 dadt_ion = -1 / 12. / PI / kb / (dielc * perm_vac) * summ;
840 }
841 }
842
843 double dadt = dadt_hc + dadt_disp + dadt_assoc + dadt_polar + dadt_ion;
844 return dadt;
845}
846
847CoolPropDbl PCSAFTBackend::calc_hmolar_residual() {
848 CoolPropDbl Z = calc_compressibility_factor();
849 CoolPropDbl dares_dt = calc_dadt();
850
851 CoolPropDbl hres = (-_T * dares_dt + (Z - 1)) * kb * N_AV * _T; // Equation A.46 from Gross and Sadowski 2001
852 return hres;
853}
854
855CoolPropDbl PCSAFTBackend::calc_smolar_residual() {
856 CoolPropDbl dares_dt = calc_dadt();
857 CoolPropDbl ares = calc_alphar();
858
859 CoolPropDbl sres = kb * N_AV * (-_T * dares_dt - ares);
860 return sres;
861}
862
863vector<CoolPropDbl> PCSAFTBackend::calc_fugacity_coefficients() {
864 auto ncomp = N; // number of components
865 vector<double> d(ncomp);
866 for (auto i = 0U; i < ncomp; i++) {
867 d[i] = components[i].getSigma() * (1 - 0.12 * exp(-3 * components[i].getU() / _T));
868 }
869 if (ion_term) {
870 for (auto i = 0U; i < ncomp; i++) {
871 if (components[i].getZ() != 0) {
872 d[i] =
873 components[i].getSigma() * (1 - 0.12); // for ions the diameter is assumed to be temperature independent (see Held et al. 2014)
874 }
875 }
876 }
877
878 double den = _rhomolar * N_AV / 1.0e30;
879
880 vector<double> zeta(4, 0);
881 double summ = NAN;
882 for (int i = 0; i < 4; i++) {
883 summ = 0;
884 for (size_t j = 0; j < ncomp; j++) {
885 summ += mole_fractions[j] * components[j].getM() * pow(d[j], i);
886 }
887 zeta[i] = PI / 6 * den * summ;
888 }
889
890 double eta = zeta[3];
891 double m_avg = 0;
892 for (size_t i = 0; i < ncomp; i++) {
893 m_avg += mole_fractions[i] * components[i].getM();
894 }
895
896 vector<double> ghs(ncomp * ncomp, 0);
897 vector<double> denghs(ncomp * ncomp, 0);
898 vector<double> e_ij(ncomp * ncomp, 0);
899 vector<double> s_ij(ncomp * ncomp, 0);
900 double m2es3 = 0.;
901 double m2e2s3 = 0.;
902 int idx = -1;
903 for (size_t i = 0; i < ncomp; i++) {
904 for (size_t j = 0; j < ncomp; j++) {
905 idx += 1;
906 s_ij[idx] = (components[i].getSigma() + components[j].getSigma()) / 2.;
907 if (ion_term) {
908 if (components[i].getZ() * components[j].getZ()
909 <= 0) { // for two cations or two anions e_ij is kept at zero to avoid dispersion between like ions (see Held et al. 2014)
910 if (k_ij.empty()) {
911 e_ij[idx] = sqrt(components[i].getU() * components[j].getU());
912 } else {
913 e_ij[idx] = sqrt(components[i].getU() * components[j].getU()) * (1 - (k_ij[idx] + k_ijT[idx] * _T));
914 }
915 }
916 } else {
917 if (k_ij.empty()) {
918 e_ij[idx] = sqrt(components[i].getU() * components[j].getU());
919 } else {
920 e_ij[idx] = sqrt(components[i].getU() * components[j].getU()) * (1 - (k_ij[idx] + k_ijT[idx] * _T));
921 }
922 }
923 m2es3 = m2es3 + mole_fractions[i] * mole_fractions[j] * components[i].getM() * components[j].getM() * e_ij[idx] / _T * pow(s_ij[idx], 3);
924 m2e2s3 =
925 m2e2s3
926 + mole_fractions[i] * mole_fractions[j] * components[i].getM() * components[j].getM() * pow(e_ij[idx] / _T, 2) * pow(s_ij[idx], 3);
927 ghs[idx] = 1 / (1 - zeta[3]) + (d[i] * d[j] / (d[i] + d[j])) * 3 * zeta[2] / (1 - zeta[3]) / (1 - zeta[3])
928 + pow(d[i] * d[j] / (d[i] + d[j]), 2) * 2 * zeta[2] * zeta[2] / pow(1 - zeta[3], 3);
929 denghs[idx] =
930 zeta[3] / (1 - zeta[3]) / (1 - zeta[3])
931 + (d[i] * d[j] / (d[i] + d[j])) * (3 * zeta[2] / (1 - zeta[3]) / (1 - zeta[3]) + 6 * zeta[2] * zeta[3] / pow(1 - zeta[3], 3))
932 + pow(d[i] * d[j] / (d[i] + d[j]), 2)
933 * (4 * zeta[2] * zeta[2] / pow(1 - zeta[3], 3) + 6 * zeta[2] * zeta[2] * zeta[3] / pow(1 - zeta[3], 4));
934 }
935 }
936
937 double ares_hs = 1 / zeta[0]
938 * (3 * zeta[1] * zeta[2] / (1 - zeta[3]) + pow(zeta[2], 3.) / (zeta[3] * pow(1 - zeta[3], 2))
939 + (pow(zeta[2], 3.) / pow(zeta[3], 2.) - zeta[0]) * log(1 - zeta[3]));
940 double Zhs = zeta[3] / (1 - zeta[3]) + 3. * zeta[1] * zeta[2] / zeta[0] / (1. - zeta[3]) / (1. - zeta[3])
941 + (3. * pow(zeta[2], 3.) - zeta[3] * pow(zeta[2], 3.)) / zeta[0] / pow(1. - zeta[3], 3.);
942
943 static double a0[7] = {0.9105631445, 0.6361281449, 2.6861347891, -26.547362491, 97.759208784, -159.59154087, 91.297774084};
944 static double a1[7] = {-0.3084016918, 0.1860531159, -2.5030047259, 21.419793629, -65.255885330, 83.318680481, -33.746922930};
945 static double a2[7] = {-0.0906148351, 0.4527842806, 0.5962700728, -1.7241829131, -4.1302112531, 13.776631870, -8.6728470368};
946 static double b0[7] = {0.7240946941, 2.2382791861, -4.0025849485, -21.003576815, 26.855641363, 206.55133841, -355.60235612};
947 static double b1[7] = {-0.5755498075, 0.6995095521, 3.8925673390, -17.215471648, 192.67226447, -161.82646165, -165.20769346};
948 static double b2[7] = {0.0976883116, -0.2557574982, -9.1558561530, 20.642075974, -38.804430052, 93.626774077, -29.666905585};
949
950 vector<double> a(7, 0);
951 vector<double> b(7, 0);
952 for (int i = 0; i < 7; i++) {
953 a[i] = a0[i] + (m_avg - 1.) / m_avg * a1[i] + (m_avg - 1.) / m_avg * (m_avg - 2.) / m_avg * a2[i];
954 b[i] = b0[i] + (m_avg - 1.) / m_avg * b1[i] + (m_avg - 1.) / m_avg * (m_avg - 2.) / m_avg * b2[i];
955 }
956
957 double detI1_det = 0.0;
958 double detI2_det = 0.0;
959 double I1 = 0.0;
960 double I2 = 0.0;
961 for (int i = 0; i < 7; i++) {
962 detI1_det += a[i] * (i + 1) * pow(eta, i);
963 detI2_det += b[i] * (i + 1) * pow(eta, i);
964 I2 += b[i] * pow(eta, i);
965 I1 += a[i] * pow(eta, i);
966 }
967 double C1 = 1.
968 / (1. + m_avg * (8 * eta - 2 * eta * eta) / pow(1 - eta, 4)
969 + (1 - m_avg) * (20 * eta - 27 * eta * eta + 12 * pow(eta, 3) - 2 * pow(eta, 4)) / pow((1 - eta) * (2 - eta), 2.0));
970 double C2 = -1. * C1 * C1
971 * (m_avg * (-4 * eta * eta + 20 * eta + 8) / pow(1 - eta, 5)
972 + (1 - m_avg) * (2 * pow(eta, 3) + 12 * eta * eta - 48 * eta + 40) / pow((1 - eta) * (2 - eta), 3.0));
973
974 summ = 0.0;
975 for (size_t i = 0; i < ncomp; i++) {
976 summ += mole_fractions[i] * (components[i].getM() - 1) * log(ghs[i * ncomp + i]);
977 }
978
979 double ares_hc = m_avg * ares_hs - summ;
980 double ares_disp = -2 * PI * den * I1 * m2es3 - PI * den * m_avg * C1 * I2 * m2e2s3;
981
982 summ = 0.0;
983 for (size_t i = 0; i < ncomp; i++) {
984 summ += mole_fractions[i] * (components[i].getM() - 1) / ghs[i * ncomp + i] * denghs[i * ncomp + i];
985 }
986
987 double Zhc = m_avg * Zhs - summ;
988 double Zdisp = -2 * PI * den * detI1_det * m2es3 - PI * den * m_avg * (C1 * detI2_det + C2 * eta * I2) * m2e2s3;
989
990 vector<double> dghsii_dx(ncomp * ncomp, 0);
991 vector<double> dahs_dx(ncomp, 0);
992 vector<double> dzeta_dx(4, 0);
993 idx = -1;
994 for (size_t i = 0; i < ncomp; i++) {
995 for (int l = 0; l < 4; l++) {
996 dzeta_dx[l] = PI / 6. * den * components[i].getM() * pow(d[i], l);
997 }
998 for (size_t j = 0; j < ncomp; j++) {
999 idx += 1;
1000 dghsii_dx[idx] =
1001 dzeta_dx[3] / (1 - zeta[3]) / (1 - zeta[3])
1002 + (d[j] * d[j] / (d[j] + d[j])) * (3 * dzeta_dx[2] / (1 - zeta[3]) / (1 - zeta[3]) + 6 * zeta[2] * dzeta_dx[3] / pow(1 - zeta[3], 3))
1003 + pow(d[j] * d[j] / (d[j] + d[j]), 2)
1004 * (4 * zeta[2] * dzeta_dx[2] / pow(1 - zeta[3], 3) + 6 * zeta[2] * zeta[2] * dzeta_dx[3] / pow(1 - zeta[3], 4));
1005 }
1006 dahs_dx[i] =
1007 -dzeta_dx[0] / zeta[0] * ares_hs
1008 + 1 / zeta[0]
1009 * (3 * (dzeta_dx[1] * zeta[2] + zeta[1] * dzeta_dx[2]) / (1 - zeta[3])
1010 + 3 * zeta[1] * zeta[2] * dzeta_dx[3] / (1 - zeta[3]) / (1 - zeta[3])
1011 + 3 * zeta[2] * zeta[2] * dzeta_dx[2] / zeta[3] / (1 - zeta[3]) / (1 - zeta[3])
1012 + pow(zeta[2], 3) * dzeta_dx[3] * (3 * zeta[3] - 1) / zeta[3] / zeta[3] / pow(1 - zeta[3], 3)
1013 + log(1 - zeta[3])
1014 * ((3 * zeta[2] * zeta[2] * dzeta_dx[2] * zeta[3] - 2 * pow(zeta[2], 3) * dzeta_dx[3]) / pow(zeta[3], 3) - dzeta_dx[0])
1015 + (zeta[0] - pow(zeta[2], 3) / zeta[3] / zeta[3]) * dzeta_dx[3] / (1 - zeta[3]));
1016 }
1017
1018 vector<double> dadisp_dx(ncomp, 0);
1019 vector<double> dahc_dx(ncomp, 0);
1020 double dzeta3_dx = NAN, daa_dx = NAN, db_dx = NAN, dI1_dx = NAN, dI2_dx = NAN, dm2es3_dx = NAN, dm2e2s3_dx = NAN, dC1_dx = NAN;
1021 for (size_t i = 0; i < ncomp; i++) {
1022 dzeta3_dx = PI / 6. * den * components[i].getM() * pow(d[i], 3);
1023 dI1_dx = 0.0;
1024 dI2_dx = 0.0;
1025 dm2es3_dx = 0.0;
1026 dm2e2s3_dx = 0.0;
1027 for (int l = 0; l < 7; l++) {
1028 daa_dx = components[i].getM() / m_avg / m_avg * a1[l] + components[i].getM() / m_avg / m_avg * (3 - 4 / m_avg) * a2[l];
1029 db_dx = components[i].getM() / m_avg / m_avg * b1[l] + components[i].getM() / m_avg / m_avg * (3 - 4 / m_avg) * b2[l];
1030 dI1_dx += a[l] * l * dzeta3_dx * pow(eta, l - 1) + daa_dx * pow(eta, l);
1031 dI2_dx += b[l] * l * dzeta3_dx * pow(eta, l - 1) + db_dx * pow(eta, l);
1032 }
1033 for (size_t j = 0; j < ncomp; j++) {
1034 dm2es3_dx += mole_fractions[j] * components[j].getM() * (e_ij[i * ncomp + j] / _T) * pow(s_ij[i * ncomp + j], 3);
1035 dm2e2s3_dx += mole_fractions[j] * components[j].getM() * pow(e_ij[i * ncomp + j] / _T, 2) * pow(s_ij[i * ncomp + j], 3);
1036 dahc_dx[i] += mole_fractions[j] * (components[j].getM() - 1) / ghs[j * ncomp + j] * dghsii_dx[i * ncomp + j];
1037 }
1038 dm2es3_dx = dm2es3_dx * 2 * components[i].getM();
1039 dm2e2s3_dx = dm2e2s3_dx * 2 * components[i].getM();
1040 dahc_dx[i] = components[i].getM() * ares_hs + m_avg * dahs_dx[i] - dahc_dx[i] - (components[i].getM() - 1) * log(ghs[i * ncomp + i]);
1041 dC1_dx = C2 * dzeta3_dx
1042 - C1 * C1
1043 * (components[i].getM() * (8 * eta - 2 * eta * eta) / pow(1 - eta, 4)
1044 - components[i].getM() * (20 * eta - 27 * eta * eta + 12 * pow(eta, 3) - 2 * pow(eta, 4)) / pow((1 - eta) * (2 - eta), 2));
1045
1046 dadisp_dx[i] =
1047 -2 * PI * den * (dI1_dx * m2es3 + I1 * dm2es3_dx)
1048 - PI * den * ((components[i].getM() * C1 * I2 + m_avg * dC1_dx * I2 + m_avg * C1 * dI2_dx) * m2e2s3 + m_avg * C1 * I2 * dm2e2s3_dx);
1049 }
1050
1051 vector<double> mu_hc(ncomp, 0);
1052 vector<double> mu_disp(ncomp, 0);
1053 for (size_t i = 0; i < ncomp; i++) {
1054 for (size_t j = 0; j < ncomp; j++) {
1055 mu_hc[i] += mole_fractions[j] * dahc_dx[j];
1056 mu_disp[i] += mole_fractions[j] * dadisp_dx[j];
1057 }
1058 mu_hc[i] = ares_hc + Zhc + dahc_dx[i] - mu_hc[i];
1059 mu_disp[i] = ares_disp + Zdisp + dadisp_dx[i] - mu_disp[i];
1060 }
1061
1062 // Dipole term (Gross and Vrabec term) --------------------------------------
1063 vector<double> mu_polar(ncomp, 0);
1064 if (polar_term) {
1065 double A2 = 0.;
1066 double A3 = 0.;
1067 double dA2_det = 0.;
1068 double dA3_det = 0.;
1069 vector<double> dA2_dx(ncomp, 0);
1070 vector<double> dA3_dx(ncomp, 0);
1071
1072 static double a0dip[5] = {0.3043504, -0.1358588, 1.4493329, 0.3556977, -2.0653308};
1073 static double a1dip[5] = {0.9534641, -1.8396383, 2.0131180, -7.3724958, 8.2374135};
1074 static double a2dip[5] = {-1.1610080, 4.5258607, 0.9751222, -12.281038, 5.9397575};
1075 static double b0dip[5] = {0.2187939, -1.1896431, 1.1626889, 0, 0};
1076 static double b1dip[5] = {-0.5873164, 1.2489132, -0.5085280, 0, 0};
1077 static double b2dip[5] = {3.4869576, -14.915974, 15.372022, 0, 0};
1078 static double c0dip[5] = {-0.0646774, 0.1975882, -0.8087562, 0.6902849, 0};
1079 static double c1dip[5] = {-0.9520876, 2.9924258, -2.3802636, -0.2701261, 0};
1080 static double c2dip[5] = {-0.6260979, 1.2924686, 1.6542783, -3.4396744, 0};
1081
1082 const static double conv = 7242.702976750923; // conversion factor, see the note below Table 2 in Gross and Vrabec 2006
1083
1084 vector<double> dipmSQ(ncomp, 0);
1085 for (size_t i = 0; i < ncomp; i++) {
1086 dipmSQ[i] = pow(components[i].getDipm(), 2.) / (components[i].getM() * components[i].getU() * pow(components[i].getSigma(), 3.)) * conv;
1087 }
1088
1089 vector<double> adip(5, 0);
1090 vector<double> bdip(5, 0);
1091 vector<double> cdip(5, 0);
1092 double J2 = NAN, dJ2_det = NAN, detJ2_det = NAN, J3 = NAN, dJ3_det = NAN, detJ3_det = NAN;
1093 double m_ij = NAN;
1094 double m_ijk = NAN;
1095 for (size_t i = 0; i < ncomp; i++) {
1096 for (size_t j = 0; j < ncomp; j++) {
1097 m_ij = sqrt(components[i].getM() * components[j].getM());
1098 if (m_ij > 2) {
1099 m_ij = 2;
1100 }
1101 J2 = 0.;
1102 dJ2_det = 0.;
1103 detJ2_det = 0.;
1104 for (int l = 0; l < 5; l++) {
1105 adip[l] = a0dip[l] + (m_ij - 1) / m_ij * a1dip[l] + (m_ij - 1) / m_ij * (m_ij - 2) / m_ij * a2dip[l];
1106 bdip[l] = b0dip[l] + (m_ij - 1) / m_ij * b1dip[l] + (m_ij - 1) / m_ij * (m_ij - 2) / m_ij * b2dip[l];
1107 J2 += (adip[l] + bdip[l] * e_ij[i * ncomp + j] / _T)
1108 * pow(eta, l); // i*ncomp+j needs to be used for e_ij because it is formatted as a 1D vector
1109 dJ2_det += (adip[l] + bdip[l] * e_ij[i * ncomp + j] / _T) * l * pow(eta, l - 1);
1110 detJ2_det += (adip[l] + bdip[l] * e_ij[i * ncomp + j] / _T) * (l + 1) * pow(eta, l);
1111 }
1112 A2 += mole_fractions[i] * mole_fractions[j] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * pow(s_ij[i * ncomp + i], 3)
1113 * pow(s_ij[j * ncomp + j], 3) / pow(s_ij[i * ncomp + j], 3) * components[i].getDipnum() * components[j].getDipnum() * dipmSQ[i]
1114 * dipmSQ[j] * J2;
1115 dA2_det += mole_fractions[i] * mole_fractions[j] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * pow(s_ij[i * ncomp + i], 3)
1116 * pow(s_ij[j * ncomp + j], 3) / pow(s_ij[i * ncomp + j], 3) * components[i].getDipnum() * components[j].getDipnum()
1117 * dipmSQ[i] * dipmSQ[j] * detJ2_det;
1118 if (i == j) {
1119 dA2_dx[i] += e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * pow(s_ij[i * ncomp + i], 3) * pow(s_ij[j * ncomp + j], 3)
1120 / pow(s_ij[i * ncomp + j], 3) * components[i].getDipnum() * components[j].getDipnum() * dipmSQ[i] * dipmSQ[j]
1121 * (mole_fractions[i] * mole_fractions[j] * dJ2_det * PI / 6. * den * components[i].getM() * pow(d[i], 3)
1122 + 2 * mole_fractions[j] * J2);
1123 } else {
1124 dA2_dx[i] += e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * pow(s_ij[i * ncomp + i], 3) * pow(s_ij[j * ncomp + j], 3)
1125 / pow(s_ij[i * ncomp + j], 3) * components[i].getDipnum() * components[j].getDipnum() * dipmSQ[i] * dipmSQ[j]
1126 * (mole_fractions[i] * mole_fractions[j] * dJ2_det * PI / 6. * den * components[i].getM() * pow(d[i], 3)
1127 + mole_fractions[j] * J2);
1128 }
1129
1130 for (size_t k = 0; k < ncomp; k++) {
1131 m_ijk = pow((components[i].getM() * components[j].getM() * components[k].getM()), 1 / 3.);
1132 if (m_ijk > 2) {
1133 m_ijk = 2;
1134 }
1135 J3 = 0.;
1136 dJ3_det = 0.;
1137 detJ3_det = 0.;
1138 for (int l = 0; l < 5; l++) {
1139 cdip[l] = c0dip[l] + (m_ijk - 1) / m_ijk * c1dip[l] + (m_ijk - 1) / m_ijk * (m_ijk - 2) / m_ijk * c2dip[l];
1140 J3 += cdip[l] * pow(eta, l);
1141 dJ3_det += cdip[l] * l * pow(eta, (l - 1));
1142 detJ3_det += cdip[l] * (l + 2) * pow(eta, (l + 1));
1143 }
1144 A3 += mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T
1145 * e_ij[k * ncomp + k] / _T * pow(s_ij[i * ncomp + i], 3) * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3)
1146 / s_ij[i * ncomp + j] / s_ij[i * ncomp + k] / s_ij[j * ncomp + k] * components[i].getDipnum() * components[j].getDipnum()
1147 * components[k].getDipnum() * dipmSQ[i] * dipmSQ[j] * dipmSQ[k] * J3;
1148 dA3_det += mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T
1149 * e_ij[k * ncomp + k] / _T * pow(s_ij[i * ncomp + i], 3) * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3)
1150 / s_ij[i * ncomp + j] / s_ij[i * ncomp + k] / s_ij[j * ncomp + k] * components[i].getDipnum()
1151 * components[j].getDipnum() * components[k].getDipnum() * dipmSQ[i] * dipmSQ[j] * dipmSQ[k] * detJ3_det;
1152 if ((i == j) && (i == k)) {
1153 dA3_dx[i] +=
1154 e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * e_ij[k * ncomp + k] / _T * pow(s_ij[i * ncomp + i], 3)
1155 * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3) / s_ij[i * ncomp + j] / s_ij[i * ncomp + k]
1156 / s_ij[j * ncomp + k] * components[i].getDipnum() * components[j].getDipnum() * components[k].getDipnum() * dipmSQ[i]
1157 * dipmSQ[j] * dipmSQ[k]
1158 * (mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * dJ3_det * PI / 6. * den * components[i].getM() * pow(d[i], 3)
1159 + 3 * mole_fractions[j] * mole_fractions[k] * J3);
1160 } else if ((i == j) || (i == k)) {
1161 dA3_dx[i] +=
1162 e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * e_ij[k * ncomp + k] / _T * pow(s_ij[i * ncomp + i], 3)
1163 * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3) / s_ij[i * ncomp + j] / s_ij[i * ncomp + k]
1164 / s_ij[j * ncomp + k] * components[i].getDipnum() * components[j].getDipnum() * components[k].getDipnum() * dipmSQ[i]
1165 * dipmSQ[j] * dipmSQ[k]
1166 * (mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * dJ3_det * PI / 6. * den * components[i].getM() * pow(d[i], 3)
1167 + 2 * mole_fractions[j] * mole_fractions[k] * J3);
1168 } else {
1169 dA3_dx[i] +=
1170 e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * e_ij[k * ncomp + k] / _T * pow(s_ij[i * ncomp + i], 3)
1171 * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3) / s_ij[i * ncomp + j] / s_ij[i * ncomp + k]
1172 / s_ij[j * ncomp + k] * components[i].getDipnum() * components[j].getDipnum() * components[k].getDipnum() * dipmSQ[i]
1173 * dipmSQ[j] * dipmSQ[k]
1174 * (mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * dJ3_det * PI / 6. * den * components[i].getM() * pow(d[i], 3)
1175 + mole_fractions[j] * mole_fractions[k] * J3);
1176 }
1177 }
1178 }
1179 }
1180
1181 A2 = -PI * den * A2;
1182 A3 = -4 / 3. * PI * PI * den * den * A3;
1183 dA2_det = -PI * den / eta * dA2_det;
1184 dA3_det = -4 / 3. * PI * PI * den / eta * den / eta * dA3_det;
1185 for (size_t i = 0; i < ncomp; i++) {
1186 dA2_dx[i] = -PI * den * dA2_dx[i];
1187 dA3_dx[i] = -4 / 3. * PI * PI * den * den * dA3_dx[i];
1188 }
1189
1190 vector<double> dapolar_dx(ncomp);
1191 for (size_t i = 0; i < ncomp; i++) {
1192 dapolar_dx[i] = (dA2_dx[i] * (1 - A3 / A2) + (dA3_dx[i] * A2 - A3 * dA2_dx[i]) / A2) / pow(1 - A3 / A2, 2);
1193 }
1194
1195 if (A2 != 0) { // when the mole fraction of the polar compounds is 0 then A2 = 0 and division by 0 occurs
1196 double ares_polar = A2 / (1 - A3 / A2);
1197 double Zpolar = eta * ((dA2_det * (1 - A3 / A2) + (dA3_det * A2 - A3 * dA2_det) / A2) / (1 - A3 / A2) / (1 - A3 / A2));
1198 for (size_t i = 0; i < ncomp; i++) {
1199 for (size_t j = 0; j < ncomp; j++) {
1200 mu_polar[i] += mole_fractions[j] * dapolar_dx[j];
1201 }
1202 mu_polar[i] = ares_polar + Zpolar + dapolar_dx[i] - mu_polar[i];
1203 }
1204 }
1205 }
1206
1207 // Association term -------------------------------------------------------
1208 vector<double> mu_assoc(ncomp, 0);
1209 if (assoc_term) {
1210 int num_sites = 0;
1211 vector<int> iA; //indices of associating compounds
1212 for (auto it = assoc_num.begin(); it != assoc_num.end(); ++it) {
1213 num_sites += *it;
1214 for (int i = 0; i < *it; i++) {
1215 iA.push_back(static_cast<int>(it - assoc_num.begin()));
1216 }
1217 }
1218
1219 vector<double> x_assoc(num_sites); // mole fractions of only the associating compounds
1220 for (int i = 0; i < num_sites; i++) {
1221 x_assoc[i] = mole_fractions[iA[i]];
1222 }
1223
1224 // these indices are necessary because we are only using 1D vectors
1225 vector<double> XA(num_sites, 0);
1226 vector<double> delta_ij(static_cast<std::size_t>(num_sites) * num_sites, 0);
1227 std::size_t idxa = 0UL;
1228 std::size_t idxi = 0UL; // index for the ii-th compound
1229 std::size_t idxj = 0UL; // index for the jj-th compound
1230 for (auto i = 0UL; i < static_cast<std::size_t>(num_sites); i++) {
1231 idxi = iA[i] * ncomp + iA[i];
1232 for (auto j = 0UL; j < static_cast<std::size_t>(num_sites); j++) {
1233 idxj = iA[j] * ncomp + iA[j];
1234 if (assoc_matrix[idxa] != 0) {
1235 double eABij = (components[iA[i]].getUAB() + components[iA[j]].getUAB()) / 2.;
1236 double volABij = sqrt(components[iA[i]].getVolA() * components[iA[j]].getVolA())
1237 * pow(sqrt(s_ij[idxi] * s_ij[idxj]) / (0.5 * (s_ij[idxi] + s_ij[idxj])), 3);
1238 delta_ij[idxa] = ghs[iA[i] * ncomp + iA[j]] * (exp(eABij / _T) - 1) * pow(s_ij[iA[i] * ncomp + iA[j]], 3) * volABij;
1239 }
1240 idxa += 1;
1241 }
1242 XA[i] = (-1 + sqrt(1 + 8 * den * delta_ij[i * num_sites + i])) / (4 * den * delta_ij[i * num_sites + i]);
1243 if (!std::isfinite(XA[i])) {
1244 XA[i] = 0.02;
1245 }
1246 }
1247
1248 vector<double> ddelta_dx(static_cast<std::size_t>(num_sites) * num_sites * ncomp, 0);
1249 int idx_ddelta = 0;
1250 for (size_t k = 0; k < ncomp; k++) {
1251 std::size_t idxi = 0UL; // index for the ii-th compound
1252 std::size_t idxj = 0UL; // index for the jj-th compound
1253 idxa = 0;
1254 for (auto i = 0UL; i < static_cast<std::size_t>(num_sites); i++) {
1255 idxi = iA[i] * ncomp + iA[i];
1256 for (auto j = 0UL; j < static_cast<std::size_t>(num_sites); j++) {
1257 idxj = iA[j] * ncomp + iA[j];
1258 if (assoc_matrix[idxa] != 0) {
1259 double eABij = (components[iA[i]].getUAB() + components[iA[j]].getUAB()) / 2.;
1260 double volABij = sqrt(components[iA[i]].getVolA() * components[iA[j]].getVolA())
1261 * pow(sqrt(s_ij[idxi] * s_ij[idxj]) / (0.5 * (s_ij[idxi] + s_ij[idxj])), 3);
1262 double dghsd_dx =
1263 PI / 6. * components[k].getM()
1264 * (pow(d[k], 3) / (1 - zeta[3]) / (1 - zeta[3])
1265 + 3 * d[iA[i]] * d[iA[j]] / (d[iA[i]] + d[iA[j]])
1266 * (d[k] * d[k] / (1 - zeta[3]) / (1 - zeta[3]) + 2 * pow(d[k], 3) * zeta[2] / pow(1 - zeta[3], 3))
1267 + 2 * pow((d[iA[i]] * d[iA[j]] / (d[iA[i]] + d[iA[j]])), 2)
1268 * (2 * d[k] * d[k] * zeta[2] / pow(1 - zeta[3], 3) + 3 * (pow(d[k], 3) * zeta[2] * zeta[2] / pow(1 - zeta[3], 4))));
1269 ddelta_dx[idx_ddelta] = dghsd_dx * (exp(eABij / _T) - 1) * pow(s_ij[iA[i] * ncomp + iA[j]], 3) * volABij;
1270 }
1271 idx_ddelta += 1;
1272 idxa += 1;
1273 }
1274 }
1275 }
1276
1277 int ctr = 0;
1278 double dif = 1000.;
1279 vector<double> XA_old = XA;
1280 while ((ctr < 100) && (dif > 1e-15)) {
1281 ctr += 1;
1282 XA = XA_find(XA_old, delta_ij, den, x_assoc);
1283 dif = 0.;
1284 for (int i = 0; i < num_sites; i++) {
1285 dif += abs(XA[i] - XA_old[i]);
1286 }
1287 for (int i = 0; i < num_sites; i++) {
1288 XA_old[i] = (XA[i] + XA_old[i]) / 2.0;
1289 }
1290 }
1291
1292 vector<double> dXA_dx(num_sites * ncomp, 0);
1293 dXA_dx = dXAdx_find(assoc_num, delta_ij, den, XA, ddelta_dx, x_assoc);
1294
1295 int ij = 0;
1296 for (size_t i = 0; i < ncomp; i++) {
1297 for (int j = 0; j < num_sites; j++) {
1298 mu_assoc[i] += mole_fractions[iA[j]] * den * dXA_dx[ij] * (1 / XA[j] - 0.5);
1299 ij += 1;
1300 }
1301 }
1302
1303 for (int i = 0; i < num_sites; i++) {
1304 mu_assoc[iA[i]] += log(XA[i]) - 0.5 * XA[i] + 0.5;
1305 }
1306 }
1307
1308 // Ion term ---------------------------------------------------------------
1309 vector<double> mu_ion(ncomp, 0);
1310 if (ion_term) {
1311 vector<double> q(ncomp);
1312 for (size_t i = 0; i < ncomp; i++) {
1313 q[i] = components[i].getZ() * E_CHRG;
1314 }
1315
1316 summ = 0.;
1317 for (size_t i = 0; i < ncomp; i++) {
1318 summ += components[i].getZ() * components[i].getZ() * mole_fractions[i];
1319 }
1320 double kappa =
1321 sqrt(den * E_CHRG * E_CHRG / kb / _T / (dielc * perm_vac) * summ); // the inverse Debye screening length. Equation 4 in Held et al. 2008.
1322
1323 if (kappa != 0) {
1324 vector<double> chi(ncomp);
1325 vector<double> sigma_k(ncomp);
1326 double summ1 = 0.;
1327 double summ2 = 0.;
1328 for (size_t i = 0; i < ncomp; i++) {
1329 chi[i] = 3 / pow(kappa * d[i], 3) * (1.5 + log(1 + kappa * d[i]) - 2 * (1 + kappa * d[i]) + 0.5 * pow(1 + kappa * d[i], 2));
1330 sigma_k[i] = -2 * chi[i] + 3 / (1 + kappa * d[i]);
1331 summ1 += q[i] * q[i] * mole_fractions[i] * sigma_k[i];
1332 summ2 += mole_fractions[i] * q[i] * q[i];
1333 }
1334
1335 for (size_t i = 0; i < ncomp; i++) {
1336 mu_ion[i] = -q[i] * q[i] * kappa / 24. / PI / kb / _T / (dielc * perm_vac) * (2 * chi[i] + summ1 / summ2);
1337 }
1338 }
1339 }
1340
1341 CoolPropDbl Z = calc_compressibility_factor();
1342
1343 vector<double> mu(ncomp, 0);
1344 vector<CoolPropDbl> fugcoef(ncomp, 0);
1345 for (size_t i = 0; i < ncomp; i++) {
1346 mu[i] = mu_hc[i] + mu_disp[i] + mu_polar[i] + mu_assoc[i] + mu_ion[i];
1347 fugcoef[i] = exp(mu[i] - log(Z)); // the fugacity coefficients
1348 }
1349
1350 return fugcoef;
1351}
1352
1353CoolPropDbl PCSAFTBackend::calc_gibbsmolar_residual() {
1354 CoolPropDbl ares = calc_alphar();
1355 CoolPropDbl Z = calc_compressibility_factor();
1356
1357 // Residual Gibbs energy on the (T, V) basis CoolProp uses everywhere
1358 // else (matches HEOS::calc_gibbsmolar_residual = R*T*(alphar + delta*
1359 // dalphar_ddelta) and satisfies g_res = h_res - T*s_res). Gross &
1360 // Sadowski's A.50 includes an additional -ln(Z)*RT term that converts
1361 // to (T, P) basis (the residual chemical potential at fixed P);
1362 // dropping it removes the inconsistency reported in #1943.
1363 CoolPropDbl gres = (ares + (Z - 1)) * kb * N_AV * _T;
1364 return gres;
1365}
1366
1367CoolPropDbl PCSAFTBackend::calc_compressibility_factor() {
1368 auto ncomp = N; // number of components
1369 vector<double> d(ncomp);
1370 for (auto i = 0UL; i < ncomp; i++) {
1371 d[i] = components[i].getSigma() * (1 - 0.12 * exp(-3 * components[i].getU() / _T));
1372 }
1373 if (ion_term) {
1374 for (auto i = 0UL; i < ncomp; i++) {
1375 if (components[i].getZ() != 0) {
1376 d[i] =
1377 components[i].getSigma() * (1 - 0.12); // for ions the diameter is assumed to be temperature independent (see Held et al. 2014)
1378 }
1379 }
1380 }
1381
1382 double den = _rhomolar * N_AV / 1.0e30;
1383
1384 vector<double> zeta(4, 0);
1385 double summ = NAN;
1386 for (int i = 0; i < 4; i++) {
1387 summ = 0;
1388 for (size_t j = 0; j < ncomp; j++) {
1389 summ += mole_fractions[j] * components[j].getM() * pow(d[j], i);
1390 }
1391 zeta[i] = PI / 6 * den * summ;
1392 }
1393
1394 double eta = zeta[3];
1395 double m_avg = 0;
1396 for (size_t i = 0; i < ncomp; i++) {
1397 m_avg += mole_fractions[i] * components[i].getM();
1398 }
1399
1400 vector<double> ghs(ncomp * ncomp, 0);
1401 vector<double> denghs(ncomp * ncomp, 0);
1402 vector<double> e_ij(ncomp * ncomp, 0);
1403 vector<double> s_ij(ncomp * ncomp, 0);
1404 double m2es3 = 0.;
1405 double m2e2s3 = 0.;
1406 int idx = -1;
1407 for (size_t i = 0; i < ncomp; i++) {
1408 for (size_t j = 0; j < ncomp; j++) {
1409 idx += 1;
1410 s_ij[idx] = (components[i].getSigma() + components[j].getSigma()) / 2.;
1411 if (ion_term) {
1412 if (components[i].getZ() * components[j].getZ()
1413 <= 0) { // for two cations or two anions e_ij is kept at zero to avoid dispersion between like ions (see Held et al. 2014)
1414 if (k_ij.empty()) {
1415 e_ij[idx] = sqrt(components[i].getU() * components[j].getU());
1416 } else {
1417 e_ij[idx] = sqrt(components[i].getU() * components[j].getU()) * (1 - (k_ij[idx] + k_ijT[idx] * _T));
1418 }
1419 }
1420 } else {
1421 if (k_ij.empty()) {
1422 e_ij[idx] = sqrt(components[i].getU() * components[j].getU());
1423 } else {
1424 e_ij[idx] = sqrt(components[i].getU() * components[j].getU()) * (1 - (k_ij[idx] + k_ijT[idx] * _T));
1425 }
1426 }
1427 m2es3 = m2es3 + mole_fractions[i] * mole_fractions[j] * components[i].getM() * components[j].getM() * e_ij[idx] / _T * pow(s_ij[idx], 3);
1428 m2e2s3 =
1429 m2e2s3
1430 + mole_fractions[i] * mole_fractions[j] * components[i].getM() * components[j].getM() * pow(e_ij[idx] / _T, 2) * pow(s_ij[idx], 3);
1431 ghs[idx] = 1 / (1 - zeta[3]) + (d[i] * d[j] / (d[i] + d[j])) * 3 * zeta[2] / (1 - zeta[3]) / (1 - zeta[3])
1432 + pow(d[i] * d[j] / (d[i] + d[j]), 2) * 2 * zeta[2] * zeta[2] / pow(1 - zeta[3], 3);
1433 denghs[idx] =
1434 zeta[3] / (1 - zeta[3]) / (1 - zeta[3])
1435 + (d[i] * d[j] / (d[i] + d[j])) * (3 * zeta[2] / (1 - zeta[3]) / (1 - zeta[3]) + 6 * zeta[2] * zeta[3] / pow(1 - zeta[3], 3))
1436 + pow(d[i] * d[j] / (d[i] + d[j]), 2)
1437 * (4 * zeta[2] * zeta[2] / pow(1 - zeta[3], 3) + 6 * zeta[2] * zeta[2] * zeta[3] / pow(1 - zeta[3], 4));
1438 }
1439 }
1440
1441 double Zhs = zeta[3] / (1 - zeta[3]) + 3. * zeta[1] * zeta[2] / zeta[0] / (1. - zeta[3]) / (1. - zeta[3])
1442 + (3. * pow(zeta[2], 3.) - zeta[3] * pow(zeta[2], 3.)) / zeta[0] / pow(1. - zeta[3], 3.);
1443
1444 static double a0[7] = {0.9105631445, 0.6361281449, 2.6861347891, -26.547362491, 97.759208784, -159.59154087, 91.297774084};
1445 static double a1[7] = {-0.3084016918, 0.1860531159, -2.5030047259, 21.419793629, -65.255885330, 83.318680481, -33.746922930};
1446 static double a2[7] = {-0.0906148351, 0.4527842806, 0.5962700728, -1.7241829131, -4.1302112531, 13.776631870, -8.6728470368};
1447 static double b0[7] = {0.7240946941, 2.2382791861, -4.0025849485, -21.003576815, 26.855641363, 206.55133841, -355.60235612};
1448 static double b1[7] = {-0.5755498075, 0.6995095521, 3.8925673390, -17.215471648, 192.67226447, -161.82646165, -165.20769346};
1449 static double b2[7] = {0.0976883116, -0.2557574982, -9.1558561530, 20.642075974, -38.804430052, 93.626774077, -29.666905585};
1450
1451 vector<double> a(7, 0);
1452 vector<double> b(7, 0);
1453 for (int i = 0; i < 7; i++) {
1454 a[i] = a0[i] + (m_avg - 1.) / m_avg * a1[i] + (m_avg - 1.) / m_avg * (m_avg - 2.) / m_avg * a2[i];
1455 b[i] = b0[i] + (m_avg - 1.) / m_avg * b1[i] + (m_avg - 1.) / m_avg * (m_avg - 2.) / m_avg * b2[i];
1456 }
1457
1458 double detI1_det = 0.0;
1459 double detI2_det = 0.0;
1460 double I2 = 0.0;
1461 for (int i = 0; i < 7; i++) {
1462 detI1_det += a[i] * (i + 1) * pow(eta, i);
1463 detI2_det += b[i] * (i + 1) * pow(eta, i);
1464 I2 += b[i] * pow(eta, i);
1465 }
1466 double C1 = 1.
1467 / (1. + m_avg * (8 * eta - 2 * eta * eta) / pow(1 - eta, 4)
1468 + (1 - m_avg) * (20 * eta - 27 * eta * eta + 12 * pow(eta, 3) - 2 * pow(eta, 4)) / pow((1 - eta) * (2 - eta), 2.0));
1469 double C2 = -1. * C1 * C1
1470 * (m_avg * (-4 * eta * eta + 20 * eta + 8) / pow(1 - eta, 5)
1471 + (1 - m_avg) * (2 * pow(eta, 3) + 12 * eta * eta - 48 * eta + 40) / pow((1 - eta) * (2 - eta), 3.0));
1472
1473 summ = 0.0;
1474 for (size_t i = 0; i < ncomp; i++) {
1475 summ += mole_fractions[i] * (components[i].getM() - 1) / ghs[i * ncomp + i] * denghs[i * ncomp + i];
1476 }
1477
1478 double Zid = 1.0;
1479 double Zhc = m_avg * Zhs - summ;
1480 double Zdisp = -2 * PI * den * detI1_det * m2es3 - PI * den * m_avg * (C1 * detI2_det + C2 * eta * I2) * m2e2s3;
1481
1482 // Dipole term (Gross and Vrabec term) --------------------------------------
1483 double Zpolar = 0;
1484 if (polar_term) {
1485 double A2 = 0.;
1486 double A3 = 0.;
1487 double dA2_det = 0.;
1488 double dA3_det = 0.;
1489 vector<double> adip(5, 0);
1490 vector<double> bdip(5, 0);
1491 vector<double> cdip(5, 0);
1492 vector<double> dipmSQ(ncomp, 0);
1493 double J2 = NAN, detJ2_det = NAN, J3 = NAN, detJ3_det = NAN;
1494
1495 static double a0dip[5] = {0.3043504, -0.1358588, 1.4493329, 0.3556977, -2.0653308};
1496 static double a1dip[5] = {0.9534641, -1.8396383, 2.0131180, -7.3724958, 8.2374135};
1497 static double a2dip[5] = {-1.1610080, 4.5258607, 0.9751222, -12.281038, 5.9397575};
1498 static double b0dip[5] = {0.2187939, -1.1896431, 1.1626889, 0, 0};
1499 static double b1dip[5] = {-0.5873164, 1.2489132, -0.5085280, 0, 0};
1500 static double b2dip[5] = {3.4869576, -14.915974, 15.372022, 0, 0};
1501 static double c0dip[5] = {-0.0646774, 0.1975882, -0.8087562, 0.6902849, 0};
1502 static double c1dip[5] = {-0.9520876, 2.9924258, -2.3802636, -0.2701261, 0};
1503 static double c2dip[5] = {-0.6260979, 1.2924686, 1.6542783, -3.4396744, 0};
1504
1505 const static double conv = 7242.702976750923; // conversion factor, see the note below Table 2 in Gross and Vrabec 2006
1506
1507 for (size_t i = 0; i < ncomp; i++) {
1508 dipmSQ[i] = pow(components[i].getDipm(), 2.) / (components[i].getM() * components[i].getU() * pow(components[i].getSigma(), 3.)) * conv;
1509 }
1510
1511 double m_ij = NAN;
1512 for (size_t i = 0; i < ncomp; i++) {
1513 for (size_t j = 0; j < ncomp; j++) {
1514 m_ij = sqrt(components[i].getM() * components[j].getM());
1515 if (m_ij > 2) {
1516 m_ij = 2;
1517 }
1518 J2 = 0.;
1519 detJ2_det = 0.;
1520 for (int l = 0; l < 5; l++) {
1521 adip[l] = a0dip[l] + (m_ij - 1) / m_ij * a1dip[l] + (m_ij - 1) / m_ij * (m_ij - 2) / m_ij * a2dip[l];
1522 bdip[l] = b0dip[l] + (m_ij - 1) / m_ij * b1dip[l] + (m_ij - 1) / m_ij * (m_ij - 2) / m_ij * b2dip[l];
1523 J2 += (adip[l] + bdip[l] * e_ij[i * ncomp + j] / _T)
1524 * pow(eta, l); // i*ncomp+j needs to be used for e_ij because it is formatted as a 1D vector
1525 detJ2_det += (adip[l] + bdip[l] * e_ij[i * ncomp + j] / _T) * (l + 1) * pow(eta, l);
1526 }
1527 A2 += mole_fractions[i] * mole_fractions[j] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * pow(s_ij[i * ncomp + i], 3)
1528 * pow(s_ij[j * ncomp + j], 3) / pow(s_ij[i * ncomp + j], 3) * components[i].getDipnum() * components[j].getDipnum() * dipmSQ[i]
1529 * dipmSQ[j] * J2;
1530 dA2_det += mole_fractions[i] * mole_fractions[j] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T * pow(s_ij[i * ncomp + i], 3)
1531 * pow(s_ij[j * ncomp + j], 3) / pow(s_ij[i * ncomp + j], 3) * components[i].getDipnum() * components[j].getDipnum()
1532 * dipmSQ[i] * dipmSQ[j] * detJ2_det;
1533 }
1534 }
1535
1536 double m_ijk = NAN;
1537 for (size_t i = 0; i < ncomp; i++) {
1538 for (size_t j = 0; j < ncomp; j++) {
1539 for (size_t k = 0; k < ncomp; k++) {
1540 m_ijk = pow((components[i].getM() * components[j].getM() * components[k].getM()), 1 / 3.);
1541 if (m_ijk > 2) {
1542 m_ijk = 2;
1543 }
1544 J3 = 0.;
1545 detJ3_det = 0.;
1546 for (int l = 0; l < 5; l++) {
1547 cdip[l] = c0dip[l] + (m_ijk - 1) / m_ijk * c1dip[l] + (m_ijk - 1) / m_ijk * (m_ijk - 2) / m_ijk * c2dip[l];
1548 J3 += cdip[l] * pow(eta, l);
1549 detJ3_det += cdip[l] * (l + 2) * pow(eta, (l + 1));
1550 }
1551 A3 += mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T
1552 * e_ij[k * ncomp + k] / _T * pow(s_ij[i * ncomp + i], 3) * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3)
1553 / s_ij[i * ncomp + j] / s_ij[i * ncomp + k] / s_ij[j * ncomp + k] * components[i].getDipnum() * components[j].getDipnum()
1554 * components[k].getDipnum() * dipmSQ[i] * dipmSQ[j] * dipmSQ[k] * J3;
1555 dA3_det += mole_fractions[i] * mole_fractions[j] * mole_fractions[k] * e_ij[i * ncomp + i] / _T * e_ij[j * ncomp + j] / _T
1556 * e_ij[k * ncomp + k] / _T * pow(s_ij[i * ncomp + i], 3) * pow(s_ij[j * ncomp + j], 3) * pow(s_ij[k * ncomp + k], 3)
1557 / s_ij[i * ncomp + j] / s_ij[i * ncomp + k] / s_ij[j * ncomp + k] * components[i].getDipnum()
1558 * components[j].getDipnum() * components[k].getDipnum() * dipmSQ[i] * dipmSQ[j] * dipmSQ[k] * detJ3_det;
1559 }
1560 }
1561 }
1562
1563 A2 = -PI * den * A2;
1564 A3 = -4 / 3. * PI * PI * den * den * A3;
1565 dA2_det = -PI * den / eta * dA2_det;
1566 dA3_det = -4 / 3. * PI * PI * den / eta * den / eta * dA3_det;
1567
1568 if (A2 != 0) { // when the mole fraction of the polar compounds is 0 then A2 = 0 and division by 0 occurs
1569 Zpolar = eta * ((dA2_det * (1 - A3 / A2) + (dA3_det * A2 - A3 * dA2_det) / A2) / (1 - A3 / A2) / (1 - A3 / A2));
1570 }
1571 }
1572
1573 // Association term -------------------------------------------------------
1574 double Zassoc = 0;
1575 if (assoc_term) {
1576 int num_sites = 0;
1577 vector<int> iA; //indices of associating compounds
1578 for (auto it = assoc_num.begin(); it != assoc_num.end(); ++it) {
1579 num_sites += *it;
1580 for (int i = 0; i < *it; i++) {
1581 iA.push_back(static_cast<int>(it - assoc_num.begin()));
1582 }
1583 }
1584
1585 vector<double> x_assoc(num_sites); // mole fractions of only the associating compounds
1586 for (int i = 0; i < num_sites; i++) {
1587 x_assoc[i] = mole_fractions[iA[i]];
1588 }
1589
1590 // these indices are necessary because we are only using 1D vectors
1591 vector<double> XA(num_sites, 0);
1592 vector<double> delta_ij(static_cast<std::size_t>(num_sites) * num_sites, 0);
1593 std::size_t idxa = 0UL;
1594 std::size_t idxi = 0UL; // index for the ii-th compound
1595 std::size_t idxj = 0UL; // index for the jj-th compound
1596 for (auto i = 0UL; i < static_cast<std::size_t>(num_sites); i++) {
1597 idxi = iA[i] * ncomp + iA[i];
1598 for (auto j = 0UL; j < static_cast<std::size_t>(num_sites); j++) {
1599 idxj = iA[j] * ncomp + iA[j];
1600 if (assoc_matrix[idxa] != 0) {
1601 double eABij = (components[iA[i]].getUAB() + components[iA[j]].getUAB()) / 2.;
1602 double volABij = sqrt(components[iA[i]].getVolA() * components[iA[j]].getVolA())
1603 * pow(sqrt(s_ij[idxi] * s_ij[idxj]) / (0.5 * (s_ij[idxi] + s_ij[idxj])), 3);
1604 delta_ij[idxa] = ghs[iA[i] * ncomp + iA[j]] * (exp(eABij / _T) - 1) * pow(s_ij[iA[i] * ncomp + iA[j]], 3) * volABij;
1605 }
1606 idxa += 1;
1607 }
1608 XA[i] = (-1 + sqrt(1 + 8 * den * delta_ij[i * num_sites + i])) / (4 * den * delta_ij[i * num_sites + i]);
1609 if (!std::isfinite(XA[i])) {
1610 XA[i] = 0.02;
1611 }
1612 }
1613
1614 vector<double> ddelta_dx(static_cast<std::size_t>(num_sites) * num_sites * ncomp, 0);
1615 int idx_ddelta = 0;
1616 for (size_t k = 0; k < ncomp; k++) {
1617 std::size_t idxi = 0UL; // index for the ii-th compound
1618 std::size_t idxj = 0UL; // index for the jj-th compound
1619 idxa = 0;
1620 for (auto i = 0UL; i < static_cast<std::size_t>(num_sites); i++) {
1621 idxi = iA[i] * ncomp + iA[i];
1622 for (auto j = 0UL; j < static_cast<std::size_t>(num_sites); j++) {
1623 idxj = iA[j] * ncomp + iA[j];
1624 if (assoc_matrix[idxa] != 0) {
1625 double eABij = (components[iA[i]].getUAB() + components[iA[j]].getUAB()) / 2.;
1626 double volABij = sqrt(components[iA[i]].getVolA() * components[iA[j]].getVolA())
1627 * pow(sqrt(s_ij[idxi] * s_ij[idxj]) / (0.5 * (s_ij[idxi] + s_ij[idxj])), 3);
1628 double dghsd_dx =
1629 PI / 6. * components[k].getM()
1630 * (pow(d[k], 3) / (1 - zeta[3]) / (1 - zeta[3])
1631 + 3 * d[iA[i]] * d[iA[j]] / (d[iA[i]] + d[iA[j]])
1632 * (d[k] * d[k] / (1 - zeta[3]) / (1 - zeta[3]) + 2 * pow(d[k], 3) * zeta[2] / pow(1 - zeta[3], 3))
1633 + 2 * pow((d[iA[i]] * d[iA[j]] / (d[iA[i]] + d[iA[j]])), 2)
1634 * (2 * d[k] * d[k] * zeta[2] / pow(1 - zeta[3], 3) + 3 * (pow(d[k], 3) * zeta[2] * zeta[2] / pow(1 - zeta[3], 4))));
1635 ddelta_dx[idx_ddelta] = dghsd_dx * (exp(eABij / _T) - 1) * pow(s_ij[iA[i] * ncomp + iA[j]], 3) * volABij;
1636 }
1637 idx_ddelta += 1;
1638 idxa += 1;
1639 }
1640 }
1641 }
1642
1643 int ctr = 0;
1644 double dif = 1000.;
1645 vector<double> XA_old = XA;
1646 while ((ctr < 100) && (dif > 1e-14)) {
1647 ctr += 1;
1648 XA = XA_find(XA_old, delta_ij, den, x_assoc);
1649 dif = 0.;
1650 for (int i = 0; i < num_sites; i++) {
1651 dif += abs(XA[i] - XA_old[i]);
1652 }
1653 for (int i = 0; i < num_sites; i++) {
1654 XA_old[i] = (XA[i] + XA_old[i]) / 2.0;
1655 }
1656 }
1657
1658 vector<double> dXA_dx(num_sites * ncomp, 0);
1659 dXA_dx = dXAdx_find(assoc_num, delta_ij, den, XA, ddelta_dx, x_assoc);
1660
1661 summ = 0.;
1662 int ij = 0;
1663 for (size_t i = 0; i < ncomp; i++) {
1664 for (int j = 0; j < num_sites; j++) {
1665 summ += mole_fractions[i] * den * mole_fractions[iA[j]] * (1 / XA[j] - 0.5) * dXA_dx[ij];
1666 ij += 1;
1667 }
1668 }
1669
1670 Zassoc = summ;
1671 }
1672
1673 // Ion term ---------------------------------------------------------------
1674 double Zion = 0;
1675 if (ion_term) {
1676 vector<double> q(ncomp);
1677 for (size_t i = 0; i < ncomp; i++) {
1678 q[i] = components[i].getZ() * E_CHRG;
1679 }
1680
1681 summ = 0.;
1682 for (size_t i = 0; i < ncomp; i++) {
1683 summ += pow(components[i].getZ(), 2.) * mole_fractions[i];
1684 }
1685
1686 double kappa =
1687 sqrt(den * E_CHRG * E_CHRG / kb / _T / (dielc * perm_vac) * summ); // the inverse Debye screening length. Equation 4 in Held et al. 2008.
1688
1689 if (kappa != 0) {
1690 double chi = NAN, sigma_k = NAN;
1691 summ = 0.;
1692 for (size_t i = 0; i < ncomp; i++) {
1693 chi = 3 / pow(kappa * d[i], 3) * (1.5 + log(1 + kappa * d[i]) - 2 * (1 + kappa * d[i]) + 0.5 * pow(1 + kappa * d[i], 2));
1694 sigma_k = -2 * chi + 3 / (1 + kappa * d[i]);
1695 summ += q[i] * q[i] * mole_fractions[i] * sigma_k;
1696 }
1697 Zion = -1 * kappa / 24. / PI / kb / _T / (dielc * perm_vac) * summ;
1698 }
1699 }
1700
1701 double Z = Zid + Zhc + Zdisp + Zpolar + Zassoc + Zion;
1702 return Z;
1703}
1704
1705void PCSAFTBackend::post_update(bool optional_checks) {
1706 // Check the values that must always be set
1707 if (!ValidNumber(_p)) {
1708 throw ValueError("p is not a valid number");
1709 }
1710 if (_T < 0) {
1711 throw ValueError("T is less than zero");
1712 }
1713 if (!ValidNumber(_T)) {
1714 throw ValueError("T is not a valid number");
1715 }
1716 if (_rhomolar < 0) {
1717 throw ValueError("rhomolar is less than zero");
1718 }
1719 if (!ValidNumber(_rhomolar)) {
1720 throw ValueError("rhomolar is not a valid number");
1721 }
1722
1723 if (optional_checks) {
1724 if (!ValidNumber(_Q)) {
1725 throw ValueError("Q is not a valid number");
1726 }
1727 if (_phase == iphase_unknown) {
1728 throw ValueError("_phase is unknown");
1729 }
1730 }
1731}
1732
1733void PCSAFTBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) {
1734 if (get_debug_level() > 10) {
1735 std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
1736 get_input_pair_short_desc(input_pair).c_str(), value1, value2)
1737 << '\n';
1738 }
1739
1740 // Mass-quality input pair on a true mixture: solve iteratively for Qmolar
1741 // before delegating to the molar-pair flash. Pure / pseudo-pure (size==1)
1742 // goes through mass_to_molar_inputs in the existing flow.
1743 if (CoolProp::is_Qmass_pair(input_pair) && mole_fractions.size() > 1) {
1744 update_Qmass_pair(input_pair, value1, value2);
1745 return;
1746 }
1747
1748 // Converting input to CoolPropDbl
1749 CoolPropDbl ld_value1 = value1, ld_value2 = value2;
1750 value1 = ld_value1;
1751 value2 = ld_value2;
1752
1753 // Clear the state
1754 clear();
1755
1756 if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) {
1757 throw ValueError("Mole fractions must be set");
1758 }
1759
1760 if (SatL->mole_fractions.empty()) {
1761 SatL->set_mole_fractions(mole_fractions);
1762 }
1763 if (SatV->mole_fractions.empty()) {
1764 SatV->set_mole_fractions(mole_fractions);
1765 double summ = 0;
1766 for (size_t i = 0; i < N; i++) {
1767 if (SatV->components[i].getZ() != 0) { // we make the assumption that ions do not appear in the vapor phase
1768 SatV->mole_fractions[i] = 0;
1769 } else {
1770 summ += SatV->mole_fractions[i];
1771 }
1772 }
1773 for (size_t i = 0; i < N; i++) {
1774 SatV->mole_fractions[i] = SatV->mole_fractions[i] / summ;
1775 }
1776 }
1777
1778 // If the inputs are in mass units, convert them to molar units
1779 mass_to_molar_inputs(input_pair, value1, value2);
1780
1781 switch (input_pair) {
1782 case PT_INPUTS:
1783 _p = value1;
1784 _T = value2;
1785 if (water_present) {
1786 components[water_idx].calc_water_sigma(_T);
1787 dielc = dielc_water(
1788 _T); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
1789 }
1790
1791 if (imposed_phase_index != iphase_not_imposed) {
1792 // Use the imposed phase index
1793 _phase = imposed_phase_index;
1794 } else {
1795 _phase = calc_phase_internal(input_pair);
1796 }
1797 _rhomolar = solver_rho_Tp(value2 /*T*/, value1 /*p*/, _phase /*phase*/);
1798 break;
1799 case QT_INPUTS:
1800 _Q = value1;
1801 _T = value2;
1802 SatL->_Q = value1;
1803 SatV->_Q = value1;
1804 SatL->_T = value2;
1805 SatV->_T = value2;
1806 _phase = iphase_twophase;
1807 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1808 if (water_present) {
1809 components[water_idx].calc_water_sigma(_T);
1810 SatL->components[water_idx].calc_water_sigma(_T);
1811 SatV->components[water_idx].calc_water_sigma(_T);
1812 dielc = dielc_water(
1813 _T); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
1814 SatL->dielc = dielc_water(_T);
1815 SatV->dielc = dielc_water(_T);
1816 }
1817 flash_QT(*this);
1818 break;
1819 case PQ_INPUTS:
1820 _p = value1;
1821 _Q = value2;
1822 SatL->_p = value1;
1823 SatV->_p = value1;
1824 SatL->_Q = value2;
1825 SatV->_Q = value2;
1826 _phase = iphase_twophase;
1827 if ((_Q < 0) || (_Q > 1)) {
1828 throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1829 }
1830 flash_PQ(*this);
1831 break;
1832 case DmolarT_INPUTS:
1833 _rhomolar = value1;
1834 _T = value2;
1835 SatL->_rhomolar = value1;
1836 SatV->_rhomolar = value1;
1837 SatL->_T = value2;
1838 SatV->_T = value2;
1839 if (water_present) {
1840 components[water_idx].calc_water_sigma(_T);
1841 SatL->components[water_idx].calc_water_sigma(_T);
1842 SatV->components[water_idx].calc_water_sigma(_T);
1843 dielc = dielc_water(
1844 _T); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
1845 SatL->dielc = dielc_water(_T);
1846 SatV->dielc = dielc_water(_T);
1847 }
1848 _p = update_DmolarT(_rhomolar);
1849
1850 if (imposed_phase_index != iphase_not_imposed) {
1851 // Use the imposed phase index
1852 _phase = imposed_phase_index;
1853 } else {
1854 _phase =
1855 calc_phase_internal(input_pair); // if in the two phase region, the pressure is updated by this function to equal the bubble point
1856 }
1857 break;
1858 case SmolarT_INPUTS:
1859 case DmolarP_INPUTS:
1860 case DmolarHmolar_INPUTS:
1861 case DmolarSmolar_INPUTS:
1862 case DmolarUmolar_INPUTS:
1863 case HmolarP_INPUTS:
1864 case PSmolar_INPUTS:
1865 case PUmolar_INPUTS:
1866 case HmolarSmolar_INPUTS:
1867 case QSmolar_INPUTS:
1868 case HmolarQ_INPUTS:
1869 case DmolarQ_INPUTS:
1870 default:
1871 throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
1872 }
1873
1874 // set Q, if not already set
1875 if (!ValidNumber(_Q)) {
1876 if (_phase == iphase_gas) {
1877 _Q = 1;
1878 } else if (_phase == iphase_liquid) {
1879 _Q = 0;
1880 }
1881 }
1882
1883 post_update();
1884}
1885
1886phases PCSAFTBackend::calc_phase_internal(CoolProp::input_pairs input_pair) {
1887 phases phase = iphase_unknown;
1888
1889 double p_input = NAN, rho_input = NAN;
1890 double p_bub = NAN, p_dew = NAN, p_equil = NAN;
1891 switch (input_pair) {
1892 case PT_INPUTS:
1893 p_input = _p;
1894 rho_input = _rhomolar;
1895 // first try to estimate without a full flash calculation
1896 _Q = 0;
1897 SatL->_Q = _Q;
1898 SatV->_Q = _Q;
1899 SatL->_T = _T;
1900 SatV->_T = _T;
1901 p_equil = estimate_flash_p(*this);
1902 if (p_input > 1.6 * p_equil) {
1903 phase = iphase_liquid;
1904 } else if (p_input < 0.5 * p_equil) {
1905 phase = iphase_gas;
1906 } else {
1907 // if the pressure is too close to the estimated bubble point, then do a full flash calculation to determine the phase
1908 _Q = 0;
1909 SatL->_Q = _Q;
1910 SatV->_Q = _Q;
1911 SatL->_T = _T;
1912 SatV->_T = _T;
1913 try {
1914 flash_QT(*this);
1915 } catch (const SolutionError& /* ex */) {
1916 phase = iphase_supercritical;
1917 break;
1918 }
1919 p_bub = _p;
1920 _p = p_input;
1921 _rhomolar = rho_input;
1922 if (_p > p_bub) {
1923 phase = iphase_liquid;
1924 } else if (_p == p_bub) {
1925 phase = iphase_twophase;
1926 } else {
1927 _Q = 1;
1928 SatL->_Q = _Q;
1929 SatV->_Q = _Q;
1930 flash_QT(*this);
1931 p_dew = _p;
1932 _p = p_input;
1933 _rhomolar = rho_input;
1934 if (_p < p_dew) {
1935 phase = iphase_gas;
1936 } else if ((_p <= p_bub) && (_p >= p_dew)) {
1937 phase = iphase_twophase;
1938 } else {
1939 phase = iphase_unknown;
1940 }
1941 }
1942 }
1943 break;
1944 case DmolarT_INPUTS: {
1945 double rho_bub = NAN, rho_dew = NAN;
1946 p_input = _p;
1947 rho_input = _rhomolar;
1948
1949 _Q = 0;
1950 SatL->_Q = _Q;
1951 SatV->_Q = _Q;
1952 SatL->_T = _T;
1953 SatV->_T = _T;
1954 try {
1955 flash_QT(*this);
1956 } catch (const SolutionError& /* ex */) {
1957 phase = iphase_supercritical;
1958 break;
1959 }
1960 rho_bub = _rhomolar;
1961 p_bub = _p;
1962 _p = p_input;
1963 _rhomolar = rho_input;
1964 if (_rhomolar > rho_bub) {
1965 phase = iphase_liquid;
1966 } else if (_rhomolar == rho_bub) {
1967 phase = iphase_twophase;
1968 _p = p_bub;
1969 _Q = 1 - (_rhomolar - SatV->_rhomolar) / (SatL->_rhomolar - SatV->_rhomolar);
1970 } else {
1971 _Q = 1;
1972 SatL->_Q = _Q;
1973 SatV->_Q = _Q;
1974 flash_QT(*this);
1975 rho_dew = _rhomolar;
1976 _p = p_input;
1977 _rhomolar = rho_input;
1978 if (_rhomolar < rho_dew) {
1979 phase = iphase_gas;
1980 } else if ((_rhomolar <= rho_bub) && (_rhomolar >= rho_dew)) {
1981 phase = iphase_twophase;
1982 _p = p_bub;
1983 _Q = 1 - (_rhomolar - SatV->_rhomolar) / (SatL->_rhomolar - SatV->_rhomolar);
1984 }
1985 }
1986 break;
1987 }
1988 default:
1989 throw ValueError(
1990 format("Phase determination for this pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
1991 }
1992
1993 return phase;
1994}
1995
1996void PCSAFTBackend::flash_QT(PCSAFTBackend& PCSAFT) {
1997 bool solution_found = false;
1998 double p_guess = 0;
1999 double p = 0;
2000 try {
2001 p_guess = estimate_flash_p(PCSAFT);
2002 p = outerTQ(p_guess, PCSAFT);
2003 solution_found = true;
2004 } catch (const SolutionError& /* ex */) { // NOLINT(bugprone-empty-catch)
2005 // Initial-guess attempt failed; fall through to the pressure sweep below.
2006 } catch (const ValueError& /* ex */) { // NOLINT(bugprone-empty-catch)
2007 // Same as above — fall through to the pressure sweep.
2008 }
2009
2010 // if solution hasn't been found, try cycling through a range of pressures
2011 if (!solution_found) {
2012 double p_lbound = -6; // here we're using log10 of the pressure
2013 double p_ubound = 9;
2014 double p_step = 0.1;
2015 p_guess = p_lbound;
2016 while (p_guess < p_ubound && !solution_found) {
2017 try {
2018 p = outerTQ(pow(10, p_guess), PCSAFT);
2019 solution_found = true;
2020 } catch (const SolutionError& /* ex */) {
2021 p_guess += p_step;
2022 } catch (const ValueError& /* ex */) {
2023 p_guess += p_step;
2024 }
2025 }
2026 }
2027
2028 if (!solution_found) {
2029 throw SolutionError("solution could not be found for TQ flash");
2030 }
2031
2032 // Load the outputs
2033 PCSAFT._p = p;
2034 PCSAFT._rhomolar = 1 / (PCSAFT._Q / PCSAFT.SatV->_rhomolar + (1 - PCSAFT._Q) / PCSAFT.SatL->_rhomolar);
2035 PCSAFT._phase = iphase_twophase;
2036}
2037
2038void PCSAFTBackend::flash_PQ(PCSAFTBackend& PCSAFT) {
2039 bool solution_found = false;
2040 double t_guess = 0;
2041 double t = 0;
2042 try {
2043 t_guess = estimate_flash_t(PCSAFT);
2044 t = outerPQ(t_guess, PCSAFT);
2045 solution_found = true;
2046 } catch (const SolutionError& /* ex */) { // NOLINT(bugprone-empty-catch)
2047 // Initial-guess attempt failed; fall through to the temperature sweep below.
2048 } catch (const ValueError& /* ex */) { // NOLINT(bugprone-empty-catch)
2049 // Same as above — fall through to the temperature sweep.
2050 }
2051
2052 // if solution hasn't been found, try calling the flash function directly with a range of initial temperatures
2053 if (!solution_found) {
2054 double t_lbound = 1;
2055 double t_ubound = 800;
2056 double t_step = 10;
2057 if (PCSAFT.ion_term) {
2058 t_lbound = 264;
2059 t_ubound = 350;
2060 }
2061 t_guess = t_ubound;
2062 while (t_guess > t_lbound && !solution_found) {
2063 try {
2064 t = outerPQ(t_guess, PCSAFT);
2065 solution_found = true;
2066 } catch (const SolutionError& /* ex */) {
2067 t_guess -= t_step;
2068 } catch (const ValueError& /* ex */) {
2069 t_guess -= t_step;
2070 }
2071 }
2072 }
2073
2074 if (!solution_found) {
2075 throw SolutionError("solution could not be found for PQ flash");
2076 }
2077
2078 // Load the outputs
2079 PCSAFT._T = t;
2080 PCSAFT._rhomolar = 1 / (PCSAFT._Q / PCSAFT.SatV->_rhomolar + (1 - PCSAFT._Q) / PCSAFT.SatL->_rhomolar);
2081 PCSAFT._phase = iphase_twophase;
2082}
2083
2084double PCSAFTBackend::outerPQ(double t_guess, PCSAFTBackend& PCSAFT) {
2085 // Based on the algorithm proposed in H. A. J. Watson, M. Vikse, T. Gundersen, and P. I. Barton, “Reliable Flash Calculations: Part 1. Nonsmooth Inside-Out Algorithms,” Ind. Eng. Chem. Res., vol. 56, no. 4, pp. 960–973, Feb. 2017, doi: 10.1021/acs.iecr.6b03956.
2086 auto ncomp = N; // number of components
2087 double TOL = 1e-8;
2088 int MAXITER = 200;
2089
2090 // Define the residual to be driven to zero
2091 class SolverInnerResid : public FuncWrapper1D
2092 {
2093 public:
2094 PCSAFTBackend& PCSAFT;
2095 CoolPropDbl kb0;
2096 vector<CoolPropDbl> u;
2097
2098 SolverInnerResid(PCSAFTBackend& PCSAFT, CoolPropDbl kb0, vector<CoolPropDbl> u) : PCSAFT(PCSAFT), kb0(kb0), u(std::move(u)) {}
2099 CoolPropDbl call(CoolPropDbl R) override {
2100 auto ncomp = PCSAFT.components.size();
2101 double error = 0;
2102
2103 vector<double> pp(ncomp, 0);
2104 double L = 0;
2105 for (auto i = 0U; i < ncomp; i++) {
2106 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2107 pp[i] = PCSAFT.mole_fractions[i] / (1 - R + kb0 * R * exp(u[i]));
2108 L += pp[i];
2109 } else {
2110 L += PCSAFT.mole_fractions[i];
2111 }
2112 }
2113 L = (1 - R) * L;
2114
2115 error = pow((L + PCSAFT._Q - 1), 2.);
2116 return error;
2117 };
2118 };
2119
2120 double x_ions = 0.; // overall mole fraction of ions in the system
2121 for (size_t i = 0; i < ncomp; i++) {
2122 if (PCSAFT.ion_term && PCSAFT.components[i].getZ() != 0) {
2123 x_ions += PCSAFT.mole_fractions[i];
2124 }
2125 }
2126
2127 // initialize variables
2128 vector<double> k(ncomp, 0), u(ncomp, 0), kprime(ncomp, 0), uprime(ncomp, 0);
2129 double Tref = t_guess - 1;
2130 double Tprime = t_guess + 1;
2131 double t = t_guess;
2132
2133 PCSAFT.SatL->_T = t; // _T must be updated because the density calculation depends on it
2134 PCSAFT.SatV->_T = t;
2135
2136 // calculate sigma for water, if it is present
2137 if (PCSAFT.water_present) {
2138 PCSAFT.components[water_idx].calc_water_sigma(t);
2139 PCSAFT.SatL->components[water_idx].calc_water_sigma(t);
2140 PCSAFT.SatV->components[water_idx].calc_water_sigma(t);
2141 PCSAFT.dielc =
2142 dielc_water(t); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
2143 PCSAFT.SatL->dielc = dielc_water(t);
2144 PCSAFT.SatV->dielc = dielc_water(t);
2145 }
2146
2147 // calculate initial guess for compositions based on fugacity coefficients and Raoult's Law.
2148 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(t, PCSAFT.SatL->_p, iphase_liquid);
2149 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(t, PCSAFT.SatV->_p, iphase_gas);
2150 if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2151 throw SolutionError("liquid and vapor densities are the same.");
2152 }
2153 vector<double> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2154 vector<double> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2155
2156 double xv_sum = 0;
2157 double xl_sum = 0;
2158 for (size_t i = 0; i < ncomp; i++) {
2159 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) { // this if statement sets k to 0 for ionic components
2160 k[i] = fugcoef_l[i] / fugcoef_v[i];
2161 } else {
2162 k[i] = 0;
2163 }
2164 PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2165 xl_sum += PCSAFT.SatL->mole_fractions[i];
2166 PCSAFT.SatV->mole_fractions[i] = k[i] * PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2167 xv_sum += PCSAFT.SatV->mole_fractions[i];
2168 }
2169
2170 if (xv_sum != 1) {
2171 for (size_t i = 0; i < ncomp; i++) {
2172 PCSAFT.SatV->mole_fractions[i] = PCSAFT.SatV->mole_fractions[i] / xv_sum;
2173 }
2174 }
2175
2176 if (xl_sum != 1) {
2177 for (size_t i = 0; i < ncomp; i++) {
2178 PCSAFT.SatL->mole_fractions[i] = PCSAFT.SatL->mole_fractions[i] / xl_sum;
2179 }
2180 }
2181
2182 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(t, PCSAFT.SatL->_p, iphase_liquid);
2183 fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2184 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(t, PCSAFT.SatV->_p, iphase_gas);
2185 fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2186 for (size_t i = 0; i < ncomp; i++) {
2187 k[i] = fugcoef_l[i] / fugcoef_v[i];
2188 }
2189
2190 PCSAFT.SatL->_T = Tprime; // _T must be updated because the density calculation depends on it
2191 PCSAFT.SatV->_T = Tprime;
2192
2193 if (PCSAFT.water_present) {
2194 PCSAFT.components[water_idx].calc_water_sigma(Tprime);
2195 PCSAFT.SatL->components[water_idx].calc_water_sigma(Tprime);
2196 PCSAFT.SatV->components[water_idx].calc_water_sigma(Tprime);
2197 PCSAFT.dielc = dielc_water(
2198 Tprime); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
2199 PCSAFT.SatL->dielc = dielc_water(Tprime);
2200 PCSAFT.SatV->dielc = dielc_water(Tprime);
2201 }
2202 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(Tprime, PCSAFT.SatL->_p, iphase_liquid);
2203 fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2204 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(Tprime, PCSAFT.SatV->_p, iphase_gas);
2205 fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2206 for (size_t i = 0; i < ncomp; i++) {
2207 kprime[i] = fugcoef_l[i] / fugcoef_v[i];
2208 }
2209
2210 vector<double> t_weight(ncomp);
2211 double t_sum = 0;
2212 for (size_t i = 0; i < ncomp; i++) {
2213 double dlnk_dt = (kprime[i] - k[i]) / (Tprime - t);
2214 t_weight[i] = PCSAFT.SatV->mole_fractions[i] * dlnk_dt / (1 + PCSAFT._Q * (k[i] - 1));
2215 t_sum += t_weight[i];
2216 }
2217
2218 double kb = 0;
2219 for (size_t i = 0; i < ncomp; i++) {
2220 double wi = t_weight[i] / t_sum;
2221 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2222 kb += wi * std::log(k[i]);
2223 }
2224 }
2225 kb = std::exp(kb);
2226
2227 t_sum = 0;
2228 for (size_t i = 0; i < ncomp; i++) {
2229 double dlnk_dt = (kprime[i] - k[i]) / (Tprime - t);
2230 t_weight[i] = PCSAFT.SatV->mole_fractions[i] * dlnk_dt / (1 + PCSAFT._Q * (kprime[i] - 1));
2231 t_sum += t_weight[i];
2232 }
2233
2234 double kbprime = 0;
2235 for (size_t i = 0; i < ncomp; i++) {
2236 double wi = t_weight[i] / t_sum;
2237 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2238 kbprime += wi * std::log(kprime[i]);
2239 }
2240 }
2241 kbprime = std::exp(kbprime);
2242 double kb0 = kbprime;
2243
2244 for (size_t i = 0; i < ncomp; i++) {
2245 u[i] = std::log(k[i] / kb);
2246 uprime[i] = std::log(kprime[i] / kbprime);
2247 }
2248
2249 double B = std::log(kbprime / kb) / (1 / Tprime - 1 / t);
2250 double A = std::log(kb) - B * (1 / t - 1 / Tref);
2251
2252 // solve
2253 SolverInnerResid resid(*this, kb0, u);
2254
2255 vector<double> pp(ncomp, 0);
2256 double maxdif = 1e10 * TOL;
2257 int itr = 0;
2258 double Rmin = 0, Rmax = 1;
2259 while (maxdif > TOL && itr < MAXITER) {
2260 // save previous values for calculating the difference at the end of the iteration
2261 vector<double> u_old = u;
2262 double A_old = A;
2263
2264 resid.u = u;
2265 double R0 = kb * PCSAFT._Q / (kb * PCSAFT._Q + kb0 * (1 - PCSAFT._Q));
2266 double R = R0;
2267 if (resid.call(R) > TOL) {
2268 R = BoundedSecant(resid, R0, Rmin, Rmax, DBL_EPSILON, TOL, MAXITER);
2269 }
2270
2271 double pp_sum = 0;
2272 double eupp_sum = 0;
2273 for (size_t i = 0; i < ncomp; i++) {
2274 pp[i] = PCSAFT.mole_fractions[i] / (1 - R + kb0 * R * std::exp(u[i]));
2275 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2276 pp_sum += pp[i];
2277 eupp_sum += std::exp(u[i]) * pp[i];
2278 }
2279 }
2280 kb = pp_sum / eupp_sum;
2281
2282 t = 1 / (1 / Tref + (std::log(kb) - A) / B);
2283 for (size_t i = 0; i < ncomp; i++) {
2284 if (x_ions == 0) {
2285 PCSAFT.SatL->mole_fractions[i] = pp[i] / pp_sum;
2286 PCSAFT.SatV->mole_fractions[i] = std::exp(u[i]) * pp[i] / eupp_sum;
2287 } else if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2288 PCSAFT.SatL->mole_fractions[i] = pp[i] / pp_sum * (1 - x_ions / (1 - PCSAFT._Q));
2289 PCSAFT.SatV->mole_fractions[i] = std::exp(u[i]) * pp[i] / eupp_sum;
2290 } else {
2291 PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 - PCSAFT._Q);
2292 PCSAFT.SatV->mole_fractions[i] = 0;
2293 }
2294 }
2295
2296 PCSAFT.SatL->_T = t; // _T must be updated because the density calculation depends on it
2297 PCSAFT.SatV->_T = t;
2298
2299 if (PCSAFT.water_present) {
2300 PCSAFT.components[water_idx].calc_water_sigma(t);
2301 PCSAFT.SatL->components[water_idx].calc_water_sigma(t);
2302 PCSAFT.SatV->components[water_idx].calc_water_sigma(t);
2303 PCSAFT.dielc = dielc_water(
2304 t); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
2305 PCSAFT.SatL->dielc = dielc_water(t);
2306 PCSAFT.SatV->dielc = dielc_water(t);
2307 }
2308 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(t, PCSAFT._p, iphase_liquid);
2309 vector<CoolPropDbl> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2310 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(t, PCSAFT._p, iphase_gas);
2311 vector<CoolPropDbl> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2312 for (size_t i = 0; i < ncomp; i++) {
2313 k[i] = fugcoef_l[i] / fugcoef_v[i];
2314 u[i] = std::log(k[i] / kb);
2315 }
2316
2317 if (itr == 0) {
2318 B = std::log(kbprime / kb) / (1 / Tprime - 1 / t);
2319 if (B > 0) {
2320 throw SolutionError("B > 0 in outerPQ");
2321 }
2322 }
2323 A = std::log(kb) - B * (1 / t - 1 / Tref);
2324
2325 maxdif = std::abs(A - A_old);
2326 for (size_t i = 0; i < ncomp; i++) {
2327 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2328 double dif = std::abs(u[i] - u_old[i]);
2329 if (dif > maxdif) {
2330 maxdif = dif;
2331 }
2332 }
2333 }
2334
2335 itr += 1;
2336 }
2337
2338 if (!std::isfinite(t) || maxdif > 1e-3 || t < 0) {
2339 throw SolutionError("outerPQ did not converge to a solution");
2340 }
2341
2342 return t;
2343}
2344
2345double PCSAFTBackend::outerTQ(double p_guess, PCSAFTBackend& PCSAFT) {
2346 // Based on the algorithm proposed in H. A. J. Watson, M. Vikse, T. Gundersen, and P. I. Barton, “Reliable Flash Calculations: Part 1. Nonsmooth Inside-Out Algorithms,” Ind. Eng. Chem. Res., vol. 56, no. 4, pp. 960–973, Feb. 2017, doi: 10.1021/acs.iecr.6b03956.
2347 auto ncomp = N; // number of components
2348 double TOL = 1e-8;
2349 int MAXITER = 200;
2350
2351 // Define the residual to be driven to zero
2352 class SolverInnerResid : public FuncWrapper1D
2353 {
2354 public:
2355 PCSAFTBackend& PCSAFT;
2356 CoolPropDbl kb0;
2357 vector<CoolPropDbl> u;
2358
2359 SolverInnerResid(PCSAFTBackend& PCSAFT, CoolPropDbl kb0, vector<CoolPropDbl> u) : PCSAFT(PCSAFT), kb0(kb0), u(std::move(u)) {}
2360 CoolPropDbl call(CoolPropDbl R) override {
2361 auto ncomp = PCSAFT.components.size();
2362 double error = 0;
2363
2364 vector<double> pp(ncomp, 0);
2365 double L = 0;
2366
2367 for (size_t i = 0; i < ncomp; i++) {
2368 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2369 pp[i] = PCSAFT.mole_fractions[i] / (1 - R + kb0 * R * exp(u[i]));
2370 L += pp[i];
2371 } else {
2372 L += PCSAFT.mole_fractions[i];
2373 }
2374 }
2375 L = (1 - R) * L;
2376
2377 error = pow((L + PCSAFT._Q - 1), 2.);
2378 return error;
2379 };
2380 };
2381
2382 double x_ions = 0.; // overall mole fraction of ions in the system
2383 for (size_t i = 0; i < ncomp; i++) {
2384 if (PCSAFT.ion_term && PCSAFT.components[i].getZ() != 0) {
2385 x_ions += PCSAFT.mole_fractions[i];
2386 }
2387 }
2388
2389 // initialize variables
2390 vector<double> k(ncomp, 0), u(ncomp, 0), kprime(ncomp, 0), uprime(ncomp, 0);
2391 double Pref = p_guess - 0.01 * p_guess;
2392 double Pprime = p_guess + 0.01 * p_guess;
2393 if (p_guess > 1e6) { // when close to the critical pressure then we need to have Pprime be less than p_guess
2394 Pprime = p_guess - 0.005 * p_guess;
2395 }
2396 double p = p_guess;
2397
2398 // calculate initial guess for compositions based on fugacity coefficients and Raoult's Law.
2399 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, p, iphase_liquid);
2400 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, p, iphase_gas);
2401 if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2402 throw SolutionError("liquid and vapor densities are the same.");
2403 }
2404 vector<double> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2405 vector<double> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2406
2407 double xv_sum = 0;
2408 double xl_sum = 0;
2409 for (size_t i = 0; i < ncomp; i++) {
2410 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) { // this if statement sets k to 0 for ionic components
2411 k[i] = fugcoef_l[i] / fugcoef_v[i];
2412 } else {
2413 k[i] = 0;
2414 }
2415 PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2416 xl_sum += PCSAFT.SatL->mole_fractions[i];
2417 PCSAFT.SatV->mole_fractions[i] = k[i] * PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2418 xv_sum += PCSAFT.SatV->mole_fractions[i];
2419 }
2420
2421 if (xv_sum != 1) {
2422 for (size_t i = 0; i < ncomp; i++) {
2423 PCSAFT.SatV->mole_fractions[i] = PCSAFT.SatV->mole_fractions[i] / xv_sum;
2424 }
2425 }
2426
2427 if (xl_sum != 1) {
2428 for (size_t i = 0; i < ncomp; i++) {
2429 PCSAFT.SatL->mole_fractions[i] = PCSAFT.SatL->mole_fractions[i] / xl_sum;
2430 }
2431 }
2432
2433 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, p, iphase_liquid);
2434 fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2435 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, p, iphase_gas);
2436 fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2437 for (size_t i = 0; i < ncomp; i++) {
2438 k[i] = fugcoef_l[i] / fugcoef_v[i];
2439 u[i] = std::log(k[i] / kb);
2440 }
2441
2442 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, Pprime, iphase_liquid);
2443 fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2444 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, Pprime, iphase_gas);
2445 fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2446 for (size_t i = 0; i < ncomp; i++) {
2447 kprime[i] = fugcoef_l[i] / fugcoef_v[i];
2448 }
2449
2450 vector<double> t_weight(ncomp);
2451 double t_sum = 0;
2452 for (size_t i = 0; i < ncomp; i++) {
2453 double dlnk_dt = (kprime[i] - k[i]) / (Pprime - p);
2454 t_weight[i] = PCSAFT.SatV->mole_fractions[i] * dlnk_dt / (1 + PCSAFT._Q * (k[i] - 1));
2455 t_sum += t_weight[i];
2456 }
2457
2458 double kb = 0;
2459 for (size_t i = 0; i < ncomp; i++) {
2460 double wi = t_weight[i] / t_sum;
2461 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2462 kb += wi * std::log(k[i]);
2463 }
2464 }
2465 kb = std::exp(kb);
2466
2467 t_sum = 0;
2468 for (size_t i = 0; i < ncomp; i++) {
2469 double dlnk_dt = (kprime[i] - k[i]) / (Pprime - p);
2470 t_weight[i] = PCSAFT.SatV->mole_fractions[i] * dlnk_dt / (1 + PCSAFT._Q * (kprime[i] - 1));
2471 t_sum += t_weight[i];
2472 }
2473
2474 double kbprime = 0;
2475 for (size_t i = 0; i < ncomp; i++) {
2476 double wi = t_weight[i] / t_sum;
2477 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2478 kbprime += wi * std::log(kprime[i]);
2479 }
2480 }
2481 kbprime = std::exp(kbprime);
2482 double kb0 = kbprime;
2483
2484 for (size_t i = 0; i < ncomp; i++) {
2485 u[i] = std::log(k[i] / kb);
2486 uprime[i] = std::log(kprime[i] / kbprime);
2487 }
2488
2489 double B = std::log(kbprime / kb) / (1 / Pprime - 1 / p);
2490 double A = std::log(kb) - B * (1 / p - 1 / Pref);
2491
2492 if (B < 0) {
2493 throw SolutionError("B < 0 in outerTQ");
2494 }
2495
2496 // solve
2497 SolverInnerResid resid(*this, kb0, u);
2498
2499 vector<double> pp(ncomp, 0);
2500 double maxdif = 1e10 * TOL;
2501 int itr = 0;
2502 double Rmin = 0, Rmax = 1;
2503 while (maxdif > TOL && itr < MAXITER) {
2504 // save previous values for calculating the difference at the end of the iteration
2505 vector<double> u_old = u;
2506 double A_old = A;
2507
2508 double R0 = kb * PCSAFT._Q / (kb * PCSAFT._Q + kb0 * (1 - PCSAFT._Q));
2509 resid.u = u;
2510 double R = R0;
2511 if (resid.call(R) > TOL) {
2512 R = BoundedSecant(resid, R0, Rmin, Rmax, DBL_EPSILON, TOL, MAXITER);
2513 }
2514
2515 double pp_sum = 0;
2516 double eupp_sum = 0;
2517 for (size_t i = 0; i < ncomp; i++) {
2518 pp[i] = PCSAFT.mole_fractions[i] / (1 - R + kb0 * R * std::exp(u[i]));
2519 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2520 pp_sum += pp[i];
2521 eupp_sum += std::exp(u[i]) * pp[i];
2522 }
2523 }
2524 kb = pp_sum / eupp_sum;
2525
2526 p = 1 / (1 / Pref + (std::log(kb) - A) / B);
2527 for (size_t i = 0; i < ncomp; i++) {
2528 if (x_ions == 0) {
2529 PCSAFT.SatL->mole_fractions[i] = pp[i] / pp_sum;
2530 PCSAFT.SatV->mole_fractions[i] = std::exp(u[i]) * pp[i] / eupp_sum;
2531 } else if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2532 PCSAFT.SatL->mole_fractions[i] = pp[i] / pp_sum * (1 - x_ions / (1 - PCSAFT._Q));
2533 PCSAFT.SatV->mole_fractions[i] = std::exp(u[i]) * pp[i] / eupp_sum;
2534 } else {
2535 PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 - PCSAFT._Q);
2536 PCSAFT.SatV->mole_fractions[i] = 0;
2537 }
2538 }
2539
2540 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, p, iphase_liquid);
2541 vector<CoolPropDbl> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2542 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, p, iphase_gas);
2543 vector<CoolPropDbl> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2544 for (size_t i = 0; i < ncomp; i++) {
2545 k[i] = fugcoef_l[i] / fugcoef_v[i];
2546 u[i] = std::log(k[i] / kb);
2547 }
2548
2549 if (itr == 0) {
2550 B = std::log(kbprime / kb) / (1 / Pprime - 1 / p);
2551 }
2552 A = std::log(kb) - B * (1 / p - 1 / Pref);
2553
2554 maxdif = std::abs(A - A_old);
2555 for (size_t i = 0; i < ncomp; i++) {
2556 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2557 double dif = std::abs(u[i] - u_old[i]);
2558 if (dif > maxdif) {
2559 maxdif = dif;
2560 } else if (!std::isfinite(dif)) {
2561 maxdif = dif;
2562 }
2563 }
2564 }
2565 itr += 1;
2566 }
2567
2568 if (!std::isfinite(p) || !std::isfinite(maxdif) || maxdif > 0.1 || p < 0) {
2569 throw SolutionError("outerTQ did not converge to a solution");
2570 }
2571
2572 return p;
2573}
2574
2575double PCSAFTBackend::estimate_flash_t(PCSAFTBackend& PCSAFT) {
2579 double t_guess = _HUGE;
2580 auto ncomp = N; // number of components
2581
2582 double x_ions = 0.; // overall mole fraction of ions in the system
2583 for (size_t i = 0; i < ncomp; i++) {
2584 if (PCSAFT.ion_term && PCSAFT.components[i].getZ() != 0) {
2585 x_ions += PCSAFT.mole_fractions[i];
2586 }
2587 }
2588
2589 bool guess_found = false;
2590 double t_step = 30;
2591 double t_start = 571;
2592 double t_lbound = 1;
2593 if (PCSAFT.ion_term) {
2594 t_step = 15;
2595 t_start = 350;
2596 t_lbound = 264;
2597 }
2598 while (!guess_found && t_start > t_lbound) {
2599 // initialize variables
2600 double Tprime = t_start - 50;
2601 double t = t_start;
2602
2603 PCSAFT.SatL->_T = t; // _T must be updated because the density calculation depends on it
2604 PCSAFT.SatV->_T = t;
2605
2606 // calculate sigma for water, if it is present
2607 if (PCSAFT.water_present) {
2608 PCSAFT.components[water_idx].calc_water_sigma(t);
2609 PCSAFT.SatL->components[water_idx].calc_water_sigma(t);
2610 PCSAFT.SatV->components[water_idx].calc_water_sigma(t);
2611 PCSAFT.dielc = dielc_water(
2612 t); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
2613 PCSAFT.SatL->dielc = dielc_water(t);
2614 PCSAFT.SatV->dielc = dielc_water(t);
2615 }
2616
2617 try {
2618 double p1 = estimate_flash_p(PCSAFT);
2619 PCSAFT.SatL->_T = Tprime;
2620 PCSAFT.SatV->_T = Tprime;
2621 double p2 = estimate_flash_p(PCSAFT);
2622 PCSAFT.SatL->_T = t; // reset to initial value
2623 PCSAFT.SatV->_T = t;
2624
2625 double slope = (std::log10(p1) - std::log10(p2)) / (1 / t - 1 / Tprime);
2626 double intercept = std::log10(p1) - slope * (1 / t);
2627 t_guess = slope / (std::log10(PCSAFT._p) - intercept);
2628 guess_found = true;
2629 } catch (const SolutionError& /* ex */) {
2630 t_start -= t_step;
2631 }
2632 }
2633
2634 if (!guess_found) {
2635 throw SolutionError("an estimate for the VLE temperature could not be found");
2636 }
2637
2638 return t_guess;
2639}
2640
2641double PCSAFTBackend::estimate_flash_p(PCSAFTBackend& PCSAFT) {
2645 double p_guess = _HUGE;
2646 auto ncomp = N; // number of components
2647
2648 double x_ions = 0.; // overall mole fraction of ions in the system
2649 for (auto i = 0U; i < ncomp; i++) {
2650 if (PCSAFT.ion_term && PCSAFT.components[i].getZ() != 0) {
2651 x_ions += PCSAFT.mole_fractions[i];
2652 }
2653 }
2654
2655 bool guess_found = false;
2656 double p_start = 10000;
2657 while (!guess_found && p_start < 1e7) {
2658 // initialize variables
2659 vector<double> k(ncomp, 0), u(ncomp, 0), kprime(ncomp, 0), uprime(ncomp, 0);
2660 double Pprime = 0.99 * p_start;
2661 double p = p_start;
2662
2663 // calculate initial guess for compositions based on fugacity coefficients and Raoult's Law.
2664 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, p, iphase_liquid);
2665 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, p, iphase_gas);
2666 if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2667 p_start = p_start + 2e5;
2668 continue;
2669 }
2670 vector<double> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2671 vector<double> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2672
2673 double xv_sum = 0;
2674 double xl_sum = 0;
2675 for (size_t i = 0; i < ncomp; i++) {
2676 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2677 k[i] = fugcoef_l[i] / fugcoef_v[i];
2678 } else {
2679 k[i] = 0; // set k to 0 for ionic components
2680 }
2681 PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2682 xl_sum += PCSAFT.SatL->mole_fractions[i];
2683 PCSAFT.SatV->mole_fractions[i] = k[i] * PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2684 xv_sum += PCSAFT.SatV->mole_fractions[i];
2685 }
2686
2687 if (xv_sum != 1) {
2688 for (size_t i = 0; i < ncomp; i++) {
2689 PCSAFT.SatV->mole_fractions[i] = PCSAFT.SatV->mole_fractions[i] / xv_sum;
2690 }
2691 }
2692
2693 if (xl_sum != 1) {
2694 for (size_t i = 0; i < ncomp; i++) {
2695 PCSAFT.SatL->mole_fractions[i] = PCSAFT.SatL->mole_fractions[i] / xl_sum;
2696 }
2697 }
2698
2699 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT.SatL->_T, p, iphase_liquid);
2700 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT.SatV->_T, p, iphase_gas);
2701 if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2702 p_start = p_start + 2e5;
2703 continue;
2704 }
2705 fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2706 fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2707 double numer = 0;
2708 double denom = 0;
2709 for (size_t i = 0; i < ncomp; i++) {
2710 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2711 numer += PCSAFT.SatL->mole_fractions[i] * fugcoef_l[i];
2712 denom += PCSAFT.SatV->mole_fractions[i] * fugcoef_v[i];
2713 }
2714 }
2715 double ratio = numer / denom;
2716
2717 PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT.SatL->_T, Pprime, iphase_liquid);
2718 PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT.SatV->_T, Pprime, iphase_gas);
2719 if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2720 p_start = p_start + 2e5;
2721 continue;
2722 }
2723 fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2724 fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2725 numer = 0;
2726 denom = 0;
2727 for (size_t i = 0; i < ncomp; i++) {
2728 if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2729 numer += PCSAFT.SatL->mole_fractions[i] * fugcoef_l[i];
2730 denom += PCSAFT.SatV->mole_fractions[i] * fugcoef_v[i];
2731 }
2732 }
2733 double ratio_prime = numer / denom;
2734
2735 double slope = (std::log10(ratio) - std::log10(ratio_prime)) / (std::log10(p) - std::log10(Pprime));
2736 double intercept = std::log10(ratio) - slope * std::log10(p);
2737 p_guess = std::pow(10, -intercept / slope);
2738
2739 guess_found = true;
2740 }
2741
2742 if (!guess_found) {
2743 throw SolutionError("an estimate for the VLE pressure could not be found");
2744 }
2745
2746 return p_guess;
2747}
2748
2749CoolPropDbl PCSAFTBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, phases phase) {
2750 // Define the residual to be driven to zero
2751 class SolverRhoResid : public FuncWrapper1D
2752 {
2753 public:
2754 PCSAFTBackend& PCSAFT;
2755 CoolPropDbl T, p;
2756
2757 SolverRhoResid(PCSAFTBackend& PCSAFT, CoolPropDbl T, CoolPropDbl p) : PCSAFT(PCSAFT), T(T), p(p) {}
2758 CoolPropDbl call(CoolPropDbl rhomolar) override {
2759 CoolPropDbl peos = PCSAFT.update_DmolarT(rhomolar);
2760 double cost = (peos - p) / p;
2761 if (ValidNumber(cost)) {
2762 return cost;
2763 } else {
2764 return 1.0e20;
2765 }
2766 };
2767 };
2768
2769 SolverRhoResid resid(*this, T, p);
2770
2771 // split into grid and find bounds for each root
2772 vector<double> x_lo, x_hi;
2773 int num_pts = 20;
2774 double limit_lower = -8; // first use a log scale for the low density region
2775 double limit_upper = -1;
2776 double rho_guess = 1e-13;
2777 double rho_guess_prev = rho_guess;
2778 double err_prev = (update_DmolarT(reduced_to_molar(rho_guess, T)) - p) / p;
2779 for (int i = 0; i < num_pts; i++) {
2780 rho_guess = pow(10, (limit_upper - limit_lower) / (double)num_pts * i + limit_lower);
2781 double err = (update_DmolarT(reduced_to_molar(rho_guess, T)) - p) / p;
2782 if (err * err_prev < 0) {
2783 x_lo.push_back(rho_guess_prev);
2784 x_hi.push_back(rho_guess);
2785 }
2786 err_prev = err;
2787 rho_guess_prev = rho_guess;
2788 }
2789
2790 limit_lower = 0.1; // for the high density region the log scale is not needed
2791 limit_upper = 0.7405;
2792 for (int i = 0; i < num_pts; i++) {
2793 rho_guess = (limit_upper - limit_lower) / (double)num_pts * i + limit_lower;
2794 double err = (update_DmolarT(reduced_to_molar(rho_guess, T)) - p) / p;
2795 if (err * err_prev < 0) {
2796 x_lo.push_back(rho_guess_prev);
2797 x_hi.push_back(rho_guess);
2798 }
2799 err_prev = err;
2800 rho_guess_prev = rho_guess;
2801 }
2802
2803 // solve for appropriate root(s)
2804 double rho = _HUGE;
2805 double x_lo_molar = 1e-8, x_hi_molar = 1e7;
2806
2807 if (x_lo.size() == 1) {
2808 // rho_guess = reduced_to_molar((x_lo[0] + x_hi[0]) / 2., T);
2809 x_lo_molar = reduced_to_molar(x_lo[0], T);
2810 x_hi_molar = reduced_to_molar(x_hi[0], T);
2811 rho = Brent(resid, x_lo_molar, x_hi_molar, DBL_EPSILON, 1e-8, 200);
2812 } else if (x_lo.size() <= 3 && !x_lo.empty()) {
2813 if ((phase == iphase_liquid) || (phase == iphase_supercritical_liquid)) {
2814 // rho_guess = reduced_to_molar((x_lo.back() + x_hi.back()) / 2., T);
2815 x_lo_molar = reduced_to_molar(x_lo.back(), T);
2816 x_hi_molar = reduced_to_molar(x_hi.back(), T);
2817 rho = Brent(resid, x_lo_molar, x_hi_molar, DBL_EPSILON, 1e-8, 200);
2818 } else if ((phase == iphase_gas) || (phase == iphase_supercritical_gas) || (phase == iphase_supercritical)) {
2819 // rho_guess = reduced_to_molar((x_lo[0] + x_hi[0]) / 40., T); // starting with a lower guess often provides better results
2820 x_lo_molar = reduced_to_molar(x_lo[0], T);
2821 x_hi_molar = reduced_to_molar(x_hi[0], T);
2822 rho = Brent(resid, x_lo_molar, x_hi_molar, DBL_EPSILON, 1e-8, 200);
2823 }
2824 } else if (x_lo.size() > 3) {
2825 // if multiple roots to check, then find the one with the minimum gibbs energy. Reference: Privat R, Gani R, Jaubert JN. Are safe results obtained when the PC-SAFT equation of state is applied to ordinary pure chemicals?. Fluid Phase Equilibria. 2010 Aug 15;295(1):76-92.
2826 double g_min = 1e60;
2827 for (size_t i = 0; i < x_lo.size(); i++) {
2828 x_lo_molar = reduced_to_molar(x_lo[i], T);
2829 x_hi_molar = reduced_to_molar(x_hi[i], T);
2830 double rho_i = Brent(resid, x_lo_molar, x_hi_molar, DBL_EPSILON, 1e-8, 200);
2831 double rho_original = this->_rhomolar;
2832 this->_rhomolar = rho_i;
2833 // Phase selection between roots at the SAME (T, P) requires the
2834 // residual Gibbs energy on the (T, P) basis — i.e. with the
2835 // -ln(Z)*RT correction. calc_gibbsmolar_residual itself is on
2836 // the (T, ρ) basis to match HEOS / satisfy g_res = h_res - T*s_res
2837 // (#1943); add the conversion term locally so phase selection
2838 // still picks the correct (lowest-g(T, P)) root.
2839 double Z_i = calc_compressibility_factor();
2840 double g_i = calc_gibbsmolar_residual() - log(Z_i) * kb * N_AV * _T;
2841 this->_rhomolar = rho_original;
2842 if (g_i < g_min) {
2843 g_min = g_i;
2844 rho = rho_i;
2845 }
2846 }
2847 } else {
2848 int num_pts = 25;
2849 double err_min = 1e40;
2850 double rho_min = NAN;
2851 for (int i = 0; i < num_pts; i++) {
2852 double rho_guess = (0.7405 - 1e-8) / (double)num_pts * i + 1e-8;
2853 double err = (update_DmolarT(reduced_to_molar(rho_guess, T)) - p) / p;
2854 if (abs(err) < err_min) {
2855 err_min = abs(err);
2856 rho_min = reduced_to_molar(rho_guess, T);
2857 }
2858 }
2859 rho = rho_min;
2860 }
2861
2862 return rho;
2863}
2864
2865CoolPropDbl PCSAFTBackend::reduced_to_molar(CoolPropDbl nu, CoolPropDbl T) {
2866 vector<CoolPropDbl> d(N);
2867 CoolPropDbl summ = 0.;
2868 for (int i = 0; i < N; i++) {
2869 d[i] = components[i].getSigma() * (1 - 0.12 * exp(-3 * components[i].getU() / T));
2870 summ += mole_fractions[i] * components[i].getM() * pow(d[i], 3.);
2871 }
2872 return 6 / PI * nu / summ * 1.0e30 / N_AV;
2873}
2874
2875CoolPropDbl PCSAFTBackend::calc_molar_mass() {
2876 double summer = 0;
2877 for (unsigned int i = 0; i < N; ++i) {
2878 summer += mole_fractions[i] * components[i].molar_mass();
2879 }
2880 return summer;
2881}
2882
2883AbstractState::PhaseMolarMasses PCSAFTBackend::calc_phase_molar_masses() {
2884 if (mole_fractions.size() == 1) {
2885 const double mm = molar_mass();
2886 return {mm, mm};
2887 }
2888 double MM_l = 0;
2889 double MM_v = 0;
2890 for (std::size_t i = 0; i < N; ++i) {
2891 const double mm_i = components[i].molar_mass();
2892 MM_l += static_cast<double>(SatL->mole_fractions[i]) * mm_i;
2893 MM_v += static_cast<double>(SatV->mole_fractions[i]) * mm_i;
2894 }
2895 return {MM_l, MM_v};
2896}
2897
2898vector<double> PCSAFTBackend::XA_find(vector<double> XA_guess, vector<double> delta_ij, double den, vector<double> x) {
2900 auto num_sites = XA_guess.size();
2901 vector<double> XA = XA_guess;
2902
2903 int idxij = -1; // index for delta_ij
2904 for (auto i = 0U; i < num_sites; i++) {
2905 double summ = 0.;
2906 for (int j = 0; j < num_sites; j++) {
2907 idxij += 1;
2908 summ += den * x[j] * XA_guess[j] * delta_ij[idxij];
2909 }
2910 XA[i] = 1. / (1. + summ);
2911 }
2912
2913 return XA;
2914}
2915
2916vector<double> PCSAFTBackend::dXAdt_find(vector<double> delta_ij, double den, vector<double> XA, vector<double> ddelta_dt, vector<double> x) {
2918 auto num_sites = XA.size();
2919 Eigen::MatrixXd B = Eigen::MatrixXd::Zero(num_sites, 1);
2920 Eigen::MatrixXd A = Eigen::MatrixXd::Zero(num_sites, num_sites);
2921
2922 double summ = NAN;
2923 int ij = 0;
2924 for (auto i = 0U; i < num_sites; i++) {
2925 summ = 0;
2926 for (int j = 0; j < num_sites; j++) {
2927 B(i) -= x[j] * XA[j] * ddelta_dt[ij];
2928 A(i, j) = x[j] * delta_ij[ij];
2929 summ += x[j] * XA[j] * delta_ij[ij];
2930 ij += 1;
2931 }
2932 A(i, i) = pow(1 + den * summ, 2.) / den;
2933 }
2934
2935 Eigen::MatrixXd solution = A.lu().solve(B); //Solves linear system of equations
2936 vector<double> dXA_dt(num_sites);
2937 for (int i = 0; i < num_sites; i++) {
2938 dXA_dt[i] = solution(i);
2939 }
2940 return dXA_dt;
2941}
2942
2943vector<double> PCSAFTBackend::dXAdx_find(vector<int> assoc_num, vector<double> delta_ij, double den, vector<double> XA, vector<double> ddelta_dx,
2944 vector<double> x) {
2947 auto num_sites = XA.size();
2948 auto ncomp = assoc_num.size();
2949 Eigen::MatrixXd B(num_sites * ncomp, 1);
2950 Eigen::MatrixXd A = Eigen::MatrixXd::Zero(num_sites * ncomp, num_sites * ncomp);
2951
2952 double sum1 = NAN, sum2 = NAN;
2953 int idx1 = 0;
2954 int ij = 0;
2955 for (auto i = 0U; i < ncomp; i++) {
2956 for (auto j = 0U; j < num_sites; j++) {
2957 sum1 = 0;
2958 for (int k = 0; k < num_sites; k++) {
2959 sum1 = sum1 + den * x[k] * (XA[k] * ddelta_dx[i * num_sites * num_sites + j * num_sites + k]);
2960 A(ij, i * num_sites + k) = XA[j] * XA[j] * den * x[k] * delta_ij[j * num_sites + k];
2961 }
2962
2963 sum2 = 0;
2964 for (int l = 0; l < assoc_num[i]; l++) {
2965 sum2 = sum2 + XA[idx1 + l] * delta_ij[idx1 * num_sites + l * num_sites + j];
2966 }
2967
2968 A(ij, ij) = A(ij, ij) + 1;
2969 B(ij) = -1 * XA[j] * XA[j] * (sum1 + sum2);
2970 ij += 1;
2971 }
2972 idx1 += assoc_num[i];
2973 }
2974
2975 Eigen::MatrixXd solution = A.lu().solve(B); //Solves linear system of equations
2976 vector<double> dXA_dx(num_sites * ncomp);
2977 for (int i = 0; i < num_sites * ncomp; i++) {
2978 dXA_dx[i] = solution(i);
2979 }
2980 return dXA_dx;
2981}
2982
2983void PCSAFTBackend::set_assoc_matrix() {
2984 vector<int>
2985 charge; // whether the association site has a partial positive charge (i.e. hydrogen), negative charge, or elements of both (e.g. for acids modelled as type 1)
2986
2987 for (int i = 0; i < N; i++) {
2988 vector<std::string> assoc_scheme = components[i].getAssocScheme();
2989 auto num_sites = 0;
2990 auto num = assoc_scheme.size();
2991 for (auto j = 0U; j < num; j++) {
2992 switch (get_scheme_index(assoc_scheme[j])) {
2993 case i1: {
2994 charge.push_back(0);
2995 num_sites += 1;
2996 break;
2997 }
2998 case i2a: {
2999 vector<int> tmp{0, 0};
3000 charge.insert(charge.end(), tmp.begin(), tmp.end());
3001 num_sites += 2;
3002 break;
3003 }
3004 case i2b: {
3005 vector<int> tmp{-1, 1};
3006 charge.insert(charge.end(), tmp.begin(), tmp.end());
3007 num_sites += 2;
3008 break;
3009 }
3010 case i3a: {
3011 vector<int> tmp{0, 0, 0};
3012 charge.insert(charge.end(), tmp.begin(), tmp.end());
3013 num_sites += 3;
3014 break;
3015 }
3016 case i3b: {
3017 vector<int> tmp{-1, -1, 1};
3018 charge.insert(charge.end(), tmp.begin(), tmp.end());
3019 num_sites += 3;
3020 break;
3021 }
3022 case i4a: {
3023 vector<int> tmp{0, 0, 0, 0};
3024 charge.insert(charge.end(), tmp.begin(), tmp.end());
3025 num_sites += 4;
3026 break;
3027 }
3028 case i4b: {
3029 vector<int> tmp{1, 1, 1, -1};
3030 charge.insert(charge.end(), tmp.begin(), tmp.end());
3031 num_sites += 4;
3032 break;
3033 }
3034 case i4c: {
3035 vector<int> tmp{-1, -1, 1, 1};
3036 charge.insert(charge.end(), tmp.begin(), tmp.end());
3037 num_sites += 4;
3038 break;
3039 }
3040 default:
3041 throw ValueError(format("%s is not a valid association type.", assoc_scheme[j]));
3042 }
3043 }
3044
3045 assoc_num.push_back(num_sites);
3046 }
3047
3048 for (auto i1 = charge.begin(); i1 != charge.end(); i1++) {
3049 for (int& i2 : charge) {
3050 if (*i1 == 0 || i2 == 0) {
3051 assoc_matrix.push_back(1);
3052 } else if (*i1 == 1 && i2 == -1) {
3053 assoc_matrix.push_back(1);
3054 } else if (*i1 == -1 && i2 == 1) {
3055 assoc_matrix.push_back(1);
3056 } else {
3057 assoc_matrix.push_back(0);
3058 }
3059 }
3060 }
3061}
3062
3063double PCSAFTBackend::dielc_water(double t) {
3079 double dielc = NAN;
3080 if (t < 263.15) {
3081 throw ValueError("The current function for the dielectric constant for water is only valid for temperatures above 263.15 K.");
3082 } else if (t <= 368.15) {
3083 dielc = 7.6555618295E-04 * _T * _T - 8.1783881423E-01 * _T + 2.5419616803E+02;
3084 } else if (t <= 443.15) {
3085 dielc = 0.0005003272124 * _T * _T - 0.6285556029 * _T + 220.4467027;
3086 } else {
3087 throw ValueError("The current function for the dielectric constant for water is only valid for temperatures less than 443.15 K.");
3088 }
3089 return dielc;
3090}
3091
3092} /* namespace CoolProp */
Generated by doxygen 1.9.4