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