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