Loading [MathJax]/extensions/TeX/AMSsymbols.js
CoolProp  6.7.0
An open-source fluid property and humid air property database
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
11 using std::vector;
12 
13 /*
14 References
15 ----------
16 * J. Gross and G. Sadowski, “Perturbed-Chain SAFT:  An Equation of State
17 Based on a Perturbation Theory for Chain Molecules,” Ind. Eng. Chem.
18 Res., vol. 40, no. 4, pp. 1244–1260, Feb. 2001.
19 * M. Kleiner and G. Sadowski, “Modeling of Polar Systems Using PCP-SAFT: 
20 An Approach to Account for Induced-Association Interactions,” J. Phys.
21 Chem. C, vol. 111, no. 43, pp. 15544–15553, Nov. 2007.
22 * Gross Joachim and Vrabec Jadran, “An equation‐of‐state contribution
23 for polar components: Dipolar molecules,” AIChE J., vol. 52, no. 3,
24 pp. 1194–1204, Feb. 2006.
25 * A. J. de Villiers, C. E. Schwarz, and A. J. Burger, “Improving
26 vapour–liquid-equilibria predictions for mixtures with non-associating polar
27 components using sPC-SAFT extended with two dipolar terms,” Fluid Phase
28 Equilibria, vol. 305, no. 2, pp. 174–184, Jun. 2011.
29 * S. H. Huang and M. Radosz, “Equation of state for small, large,
30 polydisperse, and associating molecules,” Ind. Eng. Chem. Res., vol. 29,
31 no. 11, pp. 2284–2294, Nov. 1990.
32 * S. H. Huang and M. Radosz, “Equation of state for small, large,
33 polydisperse, and associating molecules: extension to fluid mixtures,”
34 Ind. 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,
36 polydisperse, and associating molecules: extension to fluid mixtures.
37 [Erratum to document cited in CA115(8):79950j],” Ind. Eng. Chem. Res.,
38 vol. 32, no. 4, pp. 762–762, Apr. 1993.
39 * J. Gross and G. Sadowski, “Application of the Perturbed-Chain SAFT
40 Equation of State to Associating Systems,” Ind. Eng. Chem. Res., vol.
41 41, no. 22, pp. 5510–5515, Oct. 2002.
42 * L. F. Cameretti, G. Sadowski, and J. M. Mollerup, “Modeling of Aqueous
43 Electrolyte Solutions with Perturbed-Chain Statistical Associated Fluid
44 Theory,” 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
46 Electrolyte Solutions with Perturbed-Chain Statistical Association Fluid
47 Theory,” Ind. Eng. Chem. Res., vol. 44, no. 23, pp. 8944–8944, Nov. 2005.
48 * C. Held, L. F. Cameretti, and G. Sadowski, “Modeling aqueous
49 electrolyte solutions: Part 1. Fully dissociated electrolytes,” Fluid
50 Phase Equilibria, vol. 270, no. 1, pp. 87–96, Aug. 2008.
51 * C. Held, T. Reschke, S. Mohammad, A. Luza, and G. Sadowski, “ePC-SAFT
52 revised,” Chem. Eng. Res. Des., vol. 92, no. 12, pp. 2884–2897, Dec. 2014.
53 */
54 
55 namespace CoolProp {
56 
57 PCSAFTBackend::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 
124 PCSAFTBackend::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 
190 void 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 
201 void 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 
220 PCSAFTBackend* 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 
226 void 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 
233 CoolPropDbl PCSAFTBackend::update_DmolarT(CoolPropDbl rho) {
234  _rhomolar = rho;
235  return this->calc_pressure();
236 }
237 
238 CoolPropDbl 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 
246 CoolPropDbl PCSAFTBackend::calc_alphar(void) {
247  int ncomp = N; // number of components
248  vector<double> d(ncomp);
249  for (int i = 0; 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 (int i = 0; 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(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  int idxa = 0;
441  int idxi = 0; // index for the ii-th compound
442  int idxj = 0; // index for the jj-th compound
443  for (int 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 
516 CoolPropDbl PCSAFTBackend::calc_dadt(void) {
517  int ncomp = N; // number of components
518  vector<double> d(ncomp), dd_dt(ncomp);
519  for (int i = 0; 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 (int i = 0; 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(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  int idxa = 0;
758  int idxi = 0; // index for the ii-th compound
759  int idxj = 0; // index for the jj-th compound
760  for (int i = 0; i < num_sites; i++) {
761  idxi = iA[i]*ncomp+iA[i];
762  for (int 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 
845 CoolPropDbl 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 
853 CoolPropDbl 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 
861 vector<CoolPropDbl> PCSAFTBackend::calc_fugacity_coefficients(void) {
862  int ncomp = N; // number of components
863  vector<double> d(ncomp);
864  for (int i = 0; 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 (int i = 0; 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(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  int idxa = 0;
1215  int idxi = 0; // index for the ii-th compound
1216  int idxj = 0; // index for the jj-th compound
1217  for (int i = 0; i < num_sites; i++) {
1218  idxi = iA[i]*ncomp+iA[i];
1219  for (int j = 0; j < 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  int idxi = 0; // index for the ii-th compound
1239  int idxj = 0; // index for the jj-th compound
1240  idxa = 0;
1241  for (int i = 0; i < num_sites; i++) {
1242  idxi = iA[i]*ncomp+iA[i];
1243  for (int j = 0; j < 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 
1340 CoolPropDbl 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 
1348 CoolPropDbl PCSAFTBackend::calc_compressibility_factor(void) {
1349  int ncomp = N; // number of components
1350  vector<double> d(ncomp);
1351  for (int i = 0; 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 (int i = 0; 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(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  int idxa = 0;
1569  int idxi = 0; // index for the ii-th compound
1570  int idxj = 0; // index for the jj-th compound
1571  for (int i = 0; i < num_sites; i++) {
1572  idxi = iA[i]*ncomp+iA[i];
1573  for (int j = 0; j < 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  int idxi = 0; // index for the ii-th compound
1593  int idxj = 0; // index for the jj-th compound
1594  idxa = 0;
1595  for (int i = 0; i < num_sites; i++) {
1596  idxi = iA[i]*ncomp+iA[i];
1597  for (int j = 0; j < 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 
1680 void 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 
1708 void 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 
1850 phases 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  }
1878  catch (const SolutionError& ex) {
1879  phase = iphase_supercritical;
1880  break;
1881  }
1882  p_bub = _p;
1883  _p = p_input; _rhomolar = rho_input;
1884  if (_p > p_bub) {
1885  phase = iphase_liquid;
1886  }
1887  else if (_p == p_bub) {
1888  phase = iphase_twophase;
1889  }
1890  else {
1891  _Q = 1;
1892  SatL->_Q = _Q; SatV->_Q = _Q;
1893  flash_QT(*this);
1894  p_dew = _p;
1895  _p = p_input; _rhomolar = rho_input;
1896  if (_p < p_dew) {
1897  phase = iphase_gas;
1898  }
1899  else if ((_p <= p_bub) && (_p >= p_dew)) {
1900  phase = iphase_twophase;
1901  }
1902  else{
1903  phase = iphase_unknown;
1904  }
1905  }
1906  }
1907  break;
1908  case DmolarT_INPUTS:
1909  double rho_bub, rho_dew;
1910  p_input = _p; rho_input = _rhomolar;
1911 
1912  _Q = 0;
1913  SatL->_Q = _Q;
1914  SatV->_Q = _Q;
1915  SatL->_T = _T;
1916  SatV->_T = _T;
1917  try {
1918  flash_QT(*this);
1919  } catch (const SolutionError& ex) {
1920  phase = iphase_supercritical;
1921  break;
1922  }
1923  rho_bub = _rhomolar;
1924  p_bub = _p;
1925  _p = p_input;
1926  _rhomolar = rho_input;
1927  if (_rhomolar > rho_bub) {
1928  phase = iphase_liquid;
1929  } else if (_rhomolar == rho_bub) {
1930  phase = iphase_twophase;
1931  _p = p_bub;
1932  _Q = 1 - (_rhomolar - SatV->_rhomolar) / (SatL->_rhomolar - SatV->_rhomolar);
1933  } else {
1934  _Q = 1;
1935  SatL->_Q = _Q;
1936  SatV->_Q = _Q;
1937  flash_QT(*this);
1938  rho_dew = _rhomolar;
1939  _p = p_input;
1940  _rhomolar = rho_input;
1941  if (_rhomolar < rho_dew) {
1942  phase = iphase_gas;
1943  } else if ((_rhomolar <= rho_bub) && (_rhomolar >= rho_dew)) {
1944  phase = iphase_twophase;
1945  _p = p_bub;
1946  _Q = 1 - (_rhomolar - SatV->_rhomolar) / (SatL->_rhomolar - SatV->_rhomolar);
1947  }
1948  }
1949  break;
1950  default:
1951  throw ValueError(
1952  format("Phase determination for this pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
1953  }
1954 
1955  return phase;
1956 }
1957 
1958 
1959 void PCSAFTBackend::flash_QT(PCSAFTBackend &PCSAFT) {
1960  bool solution_found = false;
1961  double p_guess = 0;
1962  double p = 0;
1963  try {
1964  p_guess = estimate_flash_p(PCSAFT);
1965  p = outerTQ(p_guess, PCSAFT);
1966  solution_found = true;
1967  }
1968  catch (const SolutionError& ex) {}
1969  catch (const ValueError& ex) {}
1970 
1971  // if solution hasn't been found, try cycling through a range of pressures
1972  if (!solution_found) {
1973  double p_lbound = -6; // here we're using log10 of the pressure
1974  double p_ubound = 9;
1975  double p_step = 0.1;
1976  p_guess = p_lbound;
1977  while (p_guess < p_ubound && !solution_found) {
1978  try {
1979  p = outerTQ(pow(10, p_guess), PCSAFT);
1980  solution_found = true;
1981  } catch (const SolutionError& ex) {
1982  p_guess += p_step;
1983  } catch (const ValueError& ex) {
1984  p_guess += p_step;
1985  }
1986  }
1987  }
1988 
1989  if (!solution_found) {
1990  throw SolutionError("solution could not be found for TQ flash");
1991  }
1992 
1993  // Load the outputs
1994  PCSAFT._p = p;
1995  PCSAFT._rhomolar = 1/(PCSAFT._Q/PCSAFT.SatV->_rhomolar + (1 - PCSAFT._Q)/PCSAFT.SatL->_rhomolar);
1996  PCSAFT._phase = iphase_twophase;
1997 }
1998 
1999 
2000 void PCSAFTBackend::flash_PQ(PCSAFTBackend &PCSAFT) {
2001  bool solution_found = false;
2002  double t_guess = 0;
2003  double t = 0;
2004  try {
2005  t_guess = estimate_flash_t(PCSAFT);
2006  t = outerPQ(t_guess, PCSAFT);
2007  solution_found = true;
2008  }
2009  catch (const SolutionError& ex) {}
2010  catch (const ValueError& ex) {}
2011 
2012  // if solution hasn't been found, try calling the flash function directly with a range of initial temperatures
2013  if (!solution_found) {
2014  double t_lbound = 1;
2015  double t_ubound = 800;
2016  double t_step = 10;
2017  if (PCSAFT.ion_term) {
2018  t_lbound = 264;
2019  t_ubound = 350;
2020  }
2021  t_guess = t_ubound;
2022  while (t_guess > t_lbound && !solution_found) {
2023  try {
2024  t = outerPQ(t_guess, PCSAFT);
2025  solution_found = true;
2026  } catch (const SolutionError& ex) {
2027  t_guess -= t_step;
2028  } catch (const ValueError& ex) {
2029  t_guess -= t_step;
2030  }
2031  }
2032  }
2033 
2034  if (!solution_found) {
2035  throw SolutionError("solution could not be found for PQ flash");
2036  }
2037 
2038  // Load the outputs
2039  PCSAFT._T = t;
2040  PCSAFT._rhomolar = 1/(PCSAFT._Q/PCSAFT.SatV->_rhomolar + (1 - PCSAFT._Q)/PCSAFT.SatL->_rhomolar);
2041  PCSAFT._phase = iphase_twophase;
2042 }
2043 
2044 
2045 double PCSAFTBackend::outerPQ(double t_guess, PCSAFTBackend &PCSAFT) {
2046  // 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.
2047  int ncomp = N; // number of components
2048  double TOL = 1e-8;
2049  double MAXITER = 200;
2050 
2051  // Define the residual to be driven to zero
2052  class SolverInnerResid : public FuncWrapper1D
2053  {
2054  public:
2055  PCSAFTBackend &PCSAFT;
2056  CoolPropDbl kb0;
2057  vector<CoolPropDbl> u;
2058 
2059  SolverInnerResid(PCSAFTBackend &PCSAFT, CoolPropDbl kb0, vector<CoolPropDbl> u)
2060  : PCSAFT(PCSAFT), kb0(kb0), u(u){}
2061  CoolPropDbl call(CoolPropDbl R){
2062  int ncomp = PCSAFT.components.size();
2063  double error = 0;
2064 
2065  vector<double> pp(ncomp, 0);
2066  double L = 0;
2067  for (int i = 0; i < ncomp; i++) {
2068  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2069  pp[i] = PCSAFT.mole_fractions[i] / (1 - R + kb0 * R * exp(u[i]));
2070  L += pp[i];
2071  } else {
2072  L += PCSAFT.mole_fractions[i];
2073  }
2074  }
2075  L = (1 - R) * L;
2076 
2077  error = pow((L + PCSAFT._Q - 1), 2.);
2078  return error;
2079  };
2080  };
2081 
2082  double x_ions = 0.; // overall mole fraction of ions in the system
2083  for (int i = 0; i < ncomp; i++) {
2084  if (PCSAFT.ion_term && PCSAFT.components[i].getZ() != 0) {
2085  x_ions += PCSAFT.mole_fractions[i];
2086  }
2087  }
2088 
2089  // initialize variables
2090  vector<double> k(ncomp, 0), u(ncomp, 0), kprime(ncomp, 0), uprime(ncomp, 0);
2091  double Tref = t_guess - 1;
2092  double Tprime = t_guess + 1;
2093  double t = t_guess;
2094 
2095  PCSAFT.SatL->_T = t; // _T must be updated because the density calculation depends on it
2096  PCSAFT.SatV->_T = t;
2097 
2098  // calculate sigma for water, if it is present
2099  if (PCSAFT.water_present) {
2100  PCSAFT.components[water_idx].calc_water_sigma(t);
2101  PCSAFT.SatL->components[water_idx].calc_water_sigma(t);
2102  PCSAFT.SatV->components[water_idx].calc_water_sigma(t);
2103  PCSAFT.dielc = dielc_water(t); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
2104  PCSAFT.SatL->dielc = dielc_water(t);
2105  PCSAFT.SatV->dielc = dielc_water(t);
2106  }
2107 
2108  // calculate initial guess for compositions based on fugacity coefficients and Raoult's Law.
2109  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(t, PCSAFT.SatL->_p, iphase_liquid);
2110  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(t, PCSAFT.SatV->_p, iphase_gas);
2111  if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2112  throw SolutionError("liquid and vapor densities are the same.");
2113  }
2114  vector<double> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2115  vector<double> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2116 
2117  double xv_sum = 0;
2118  double xl_sum = 0;
2119  for (int i = 0; i < ncomp; i++) {
2120  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) { // this if statement sets k to 0 for ionic components
2121  k[i] = fugcoef_l[i] / fugcoef_v[i];
2122  } else {
2123  k[i] = 0;
2124  }
2125  PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2126  xl_sum += PCSAFT.SatL->mole_fractions[i];
2127  PCSAFT.SatV->mole_fractions[i] = k[i] * PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2128  xv_sum += PCSAFT.SatV->mole_fractions[i];
2129  }
2130 
2131  if (xv_sum != 1) {
2132  for (int i = 0; i < ncomp; i++) {
2133  PCSAFT.SatV->mole_fractions[i] = PCSAFT.SatV->mole_fractions[i] / xv_sum;
2134  }
2135  }
2136 
2137  if (xl_sum != 1) {
2138  for (int i = 0; i < ncomp; i++) {
2139  PCSAFT.SatL->mole_fractions[i] = PCSAFT.SatL->mole_fractions[i] / xl_sum;
2140  }
2141  }
2142 
2143  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(t, PCSAFT.SatL->_p, iphase_liquid);
2144  fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2145  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(t, PCSAFT.SatV->_p, iphase_gas);
2146  fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2147  for (int i = 0; i < ncomp; i++) {
2148  k[i] = fugcoef_l[i] / fugcoef_v[i];
2149  }
2150 
2151  PCSAFT.SatL->_T = Tprime; // _T must be updated because the density calculation depends on it
2152  PCSAFT.SatV->_T = Tprime;
2153 
2154  if (PCSAFT.water_present) {
2155  PCSAFT.components[water_idx].calc_water_sigma(Tprime);
2156  PCSAFT.SatL->components[water_idx].calc_water_sigma(Tprime);
2157  PCSAFT.SatV->components[water_idx].calc_water_sigma(Tprime);
2158  PCSAFT.dielc = dielc_water(Tprime); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
2159  PCSAFT.SatL->dielc = dielc_water(Tprime);
2160  PCSAFT.SatV->dielc = dielc_water(Tprime);
2161  }
2162  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(Tprime, PCSAFT.SatL->_p, iphase_liquid);
2163  fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2164  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(Tprime, PCSAFT.SatV->_p, iphase_gas);
2165  fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2166  for (int i = 0; i < ncomp; i++) {
2167  kprime[i] = fugcoef_l[i] / fugcoef_v[i];
2168  }
2169 
2170  vector<double> t_weight(ncomp);
2171  double t_sum = 0;
2172  for (int i = 0; i < ncomp; i++) {
2173  double dlnk_dt = (kprime[i] - k[i]) / (Tprime - t);
2174  t_weight[i] = PCSAFT.SatV->mole_fractions[i] * dlnk_dt / (1 + PCSAFT._Q * (k[i] - 1));
2175  t_sum += t_weight[i];
2176  }
2177 
2178  double kb = 0;
2179  for (int i = 0; i < ncomp; i++) {
2180  double wi = t_weight[i] / t_sum;
2181  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2182  kb += wi * std::log(k[i]);
2183  }
2184  }
2185  kb = std::exp(kb);
2186 
2187  t_sum = 0;
2188  for (int i = 0; i < ncomp; i++) {
2189  double dlnk_dt = (kprime[i] - k[i]) / (Tprime - t);
2190  t_weight[i] = PCSAFT.SatV->mole_fractions[i] * dlnk_dt / (1 + PCSAFT._Q * (kprime[i] - 1));
2191  t_sum += t_weight[i];
2192  }
2193 
2194  double kbprime = 0;
2195  for (int i = 0; i < ncomp; i++) {
2196  double wi = t_weight[i] / t_sum;
2197  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2198  kbprime += wi * std::log(kprime[i]);
2199  }
2200  }
2201  kbprime = std::exp(kbprime);
2202  double kb0 = kbprime;
2203 
2204  for (int i = 0; i < ncomp; i++) {
2205  u[i] = std::log(k[i] / kb);
2206  uprime[i] = std::log(kprime[i] / kbprime);
2207  }
2208 
2209  double B = std::log(kbprime / kb) / (1/Tprime - 1/t);
2210  double A = std::log(kb) - B * (1/t - 1/Tref);
2211 
2212  // solve
2213  SolverInnerResid resid(*this, kb0, u);
2214 
2215  vector<double> pp(ncomp, 0);
2216  double maxdif = 1e10 * TOL;
2217  int itr = 0;
2218  double Rmin = 0, Rmax = 1;
2219  while (maxdif > TOL && itr < MAXITER) {
2220  // save previous values for calculating the difference at the end of the iteration
2221  vector<double> u_old = u;
2222  double A_old = A;
2223 
2224  resid.u = u;
2225  double R0 = kb * PCSAFT._Q / (kb * PCSAFT._Q + kb0 * (1 - PCSAFT._Q));
2226  double R = R0;
2227  if (resid.call(R) > TOL) {
2228  R = BoundedSecant(resid, R0, Rmin, Rmax, DBL_EPSILON, TOL, MAXITER);
2229  }
2230 
2231  double pp_sum = 0;
2232  double eupp_sum = 0;
2233  for (int i = 0; i < ncomp; i++) {
2234  pp[i] = PCSAFT.mole_fractions[i] / (1 - R + kb0 * R * std::exp(u[i]));
2235  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2236  pp_sum += pp[i];
2237  eupp_sum += std::exp(u[i]) * pp[i];
2238  }
2239  }
2240  kb = pp_sum / eupp_sum;
2241 
2242  t = 1 / (1 / Tref + (std::log(kb) - A) / B);
2243  for (int i = 0; i < ncomp; i++) {
2244  if (x_ions == 0) {
2245  PCSAFT.SatL->mole_fractions[i] = pp[i] / pp_sum;
2246  PCSAFT.SatV->mole_fractions[i] = std::exp(u[i]) * pp[i] / eupp_sum;
2247  }
2248  else if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2249  PCSAFT.SatL->mole_fractions[i] = pp[i] / pp_sum * (1 - x_ions / (1 - PCSAFT._Q));
2250  PCSAFT.SatV->mole_fractions[i] = std::exp(u[i]) * pp[i] / eupp_sum;
2251  }
2252  else {
2253  PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 - PCSAFT._Q);
2254  PCSAFT.SatV->mole_fractions[i] = 0;
2255  }
2256  }
2257 
2258  PCSAFT.SatL->_T = t; // _T must be updated because the density calculation depends on it
2259  PCSAFT.SatV->_T = t;
2260 
2261  if (PCSAFT.water_present) {
2262  PCSAFT.components[water_idx].calc_water_sigma(t);
2263  PCSAFT.SatL->components[water_idx].calc_water_sigma(t);
2264  PCSAFT.SatV->components[water_idx].calc_water_sigma(t);
2265  PCSAFT.dielc = dielc_water(t); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
2266  PCSAFT.SatL->dielc = dielc_water(t);
2267  PCSAFT.SatV->dielc = dielc_water(t);
2268  }
2269  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(t, PCSAFT._p, iphase_liquid);
2270  vector<CoolPropDbl> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2271  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(t, PCSAFT._p, iphase_gas);
2272  vector<CoolPropDbl> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2273  for (int i = 0; i < ncomp; i++) {
2274  k[i] = fugcoef_l[i] / fugcoef_v[i];
2275  u[i] = std::log(k[i] / kb);
2276  }
2277 
2278  if (itr == 0) {
2279  B = std::log(kbprime / kb) / (1/Tprime - 1/t);
2280  if (B > 0) {
2281  throw SolutionError("B > 0 in outerPQ");
2282  }
2283  }
2284  A = std::log(kb) - B * (1/t - 1/Tref);
2285 
2286  maxdif = std::abs(A - A_old);
2287  for (int i = 0; i < ncomp; i++) {
2288  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2289  double dif = std::abs(u[i] - u_old[i]);
2290  if (dif > maxdif) {
2291  maxdif = dif;
2292  }
2293  }
2294  }
2295 
2296  itr += 1;
2297  }
2298 
2299  if (!std::isfinite(t) || maxdif > 1e-3 || t < 0) {
2300  throw SolutionError("outerPQ did not converge to a solution");
2301  }
2302 
2303  return t;
2304 }
2305 
2306 
2307 double PCSAFTBackend::outerTQ(double p_guess, PCSAFTBackend &PCSAFT) {
2308  // 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.
2309  int ncomp = N; // number of components
2310  double TOL = 1e-8;
2311  double MAXITER = 200;
2312 
2313  // Define the residual to be driven to zero
2314  class SolverInnerResid : public FuncWrapper1D
2315  {
2316  public:
2317  PCSAFTBackend &PCSAFT;
2318  CoolPropDbl kb0;
2319  vector<CoolPropDbl> u;
2320 
2321  SolverInnerResid(PCSAFTBackend &PCSAFT, CoolPropDbl kb0, vector<CoolPropDbl> u)
2322  : PCSAFT(PCSAFT), kb0(kb0), u(u){}
2323  CoolPropDbl call(CoolPropDbl R){
2324  int ncomp = PCSAFT.components.size();
2325  double error = 0;
2326 
2327  vector<double> pp(ncomp, 0);
2328  double L = 0;
2329 
2330  for (int i = 0; i < ncomp; i++) {
2331  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2332  pp[i] = PCSAFT.mole_fractions[i] / (1 - R + kb0 * R * exp(u[i]));
2333  L += pp[i];
2334  } else {
2335  L += PCSAFT.mole_fractions[i];
2336  }
2337  }
2338  L = (1 - R) * L;
2339 
2340  error = pow((L + PCSAFT._Q - 1), 2.);
2341  return error;
2342  };
2343  };
2344 
2345  double x_ions = 0.; // overall mole fraction of ions in the system
2346  for (int i = 0; i < ncomp; i++) {
2347  if (PCSAFT.ion_term && PCSAFT.components[i].getZ() != 0) {
2348  x_ions += PCSAFT.mole_fractions[i];
2349  }
2350  }
2351 
2352  // initialize variables
2353  vector<double> k(ncomp, 0), u(ncomp, 0), kprime(ncomp, 0), uprime(ncomp, 0);
2354  double Pref = p_guess - 0.01 * p_guess;
2355  double Pprime = p_guess + 0.01 * p_guess;
2356  if (p_guess > 1e6) { // when close to the critical pressure then we need to have Pprime be less than p_guess
2357  Pprime = p_guess - 0.005 * p_guess;
2358  }
2359  double p = p_guess;
2360 
2361  // calculate initial guess for compositions based on fugacity coefficients and Raoult's Law.
2362  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, p, iphase_liquid);
2363  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, p, iphase_gas);
2364  if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2365  throw SolutionError("liquid and vapor densities are the same.");
2366  }
2367  vector<double> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2368  vector<double> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2369 
2370  double xv_sum = 0;
2371  double xl_sum = 0;
2372  for (int i = 0; i < ncomp; i++) {
2373  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) { // this if statement sets k to 0 for ionic components
2374  k[i] = fugcoef_l[i] / fugcoef_v[i];
2375  } else {
2376  k[i] = 0;
2377  }
2378  PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2379  xl_sum += PCSAFT.SatL->mole_fractions[i];
2380  PCSAFT.SatV->mole_fractions[i] = k[i] * PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2381  xv_sum += PCSAFT.SatV->mole_fractions[i];
2382  }
2383 
2384  if (xv_sum != 1) {
2385  for (int i = 0; i < ncomp; i++) {
2386  PCSAFT.SatV->mole_fractions[i] = PCSAFT.SatV->mole_fractions[i] / xv_sum;
2387  }
2388  }
2389 
2390  if (xl_sum != 1) {
2391  for (int i = 0; i < ncomp; i++) {
2392  PCSAFT.SatL->mole_fractions[i] = PCSAFT.SatL->mole_fractions[i] / xl_sum;
2393  }
2394  }
2395 
2396  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, p, iphase_liquid);
2397  fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2398  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, p, iphase_gas);
2399  fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2400  for (int i = 0; i < ncomp; i++) {
2401  k[i] = fugcoef_l[i] / fugcoef_v[i];
2402  u[i] = std::log(k[i] / kb);
2403  }
2404 
2405  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, Pprime, iphase_liquid);
2406  fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2407  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, Pprime, iphase_gas);
2408  fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2409  for (int i = 0; i < ncomp; i++) {
2410  kprime[i] = fugcoef_l[i] / fugcoef_v[i];
2411  }
2412 
2413  vector<double> t_weight(ncomp);
2414  double t_sum = 0;
2415  for (int i = 0; i < ncomp; i++) {
2416  double dlnk_dt = (kprime[i] - k[i]) / (Pprime - p);
2417  t_weight[i] = PCSAFT.SatV->mole_fractions[i] * dlnk_dt / (1 + PCSAFT._Q * (k[i] - 1));
2418  t_sum += t_weight[i];
2419  }
2420 
2421  double kb = 0;
2422  for (int i = 0; i < ncomp; i++) {
2423  double wi = t_weight[i] / t_sum;
2424  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2425  kb += wi * std::log(k[i]);
2426  }
2427  }
2428  kb = std::exp(kb);
2429 
2430  t_sum = 0;
2431  for (int i = 0; i < ncomp; i++) {
2432  double dlnk_dt = (kprime[i] - k[i]) / (Pprime - p);
2433  t_weight[i] = PCSAFT.SatV->mole_fractions[i] * dlnk_dt / (1 + PCSAFT._Q * (kprime[i] - 1));
2434  t_sum += t_weight[i];
2435  }
2436 
2437  double kbprime = 0;
2438  for (int i = 0; i < ncomp; i++) {
2439  double wi = t_weight[i] / t_sum;
2440  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2441  kbprime += wi * std::log(kprime[i]);
2442  }
2443  }
2444  kbprime = std::exp(kbprime);
2445  double kb0 = kbprime;
2446 
2447  for (int i = 0; i < ncomp; i++) {
2448  u[i] = std::log(k[i] / kb);
2449  uprime[i] = std::log(kprime[i] / kbprime);
2450  }
2451 
2452  double B = std::log(kbprime / kb) / (1/Pprime - 1/p);
2453  double A = std::log(kb) - B * (1/p - 1/Pref);
2454 
2455  if (B < 0) {
2456  throw SolutionError("B < 0 in outerTQ");
2457  }
2458 
2459  // solve
2460  SolverInnerResid resid(*this, kb0, u);
2461 
2462  vector<double> pp(ncomp, 0);
2463  double maxdif = 1e10 * TOL;
2464  int itr = 0;
2465  double Rmin = 0, Rmax = 1;
2466  while (maxdif > TOL && itr < MAXITER) {
2467  // save previous values for calculating the difference at the end of the iteration
2468  vector<double> u_old = u;
2469  double A_old = A;
2470 
2471  double R0 = kb * PCSAFT._Q / (kb * PCSAFT._Q + kb0 * (1 - PCSAFT._Q));
2472  resid.u = u;
2473  double R = R0;
2474  if (resid.call(R) > TOL) {
2475  R = BoundedSecant(resid, R0, Rmin, Rmax, DBL_EPSILON, TOL, MAXITER);
2476  }
2477 
2478  double pp_sum = 0;
2479  double eupp_sum = 0;
2480  for (int i = 0; i < ncomp; i++) {
2481  pp[i] = PCSAFT.mole_fractions[i] / (1 - R + kb0 * R * std::exp(u[i]));
2482  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2483  pp_sum += pp[i];
2484  eupp_sum += std::exp(u[i]) * pp[i];
2485  }
2486  }
2487  kb = pp_sum / eupp_sum;
2488 
2489  p = 1 / (1 / Pref + (std::log(kb) - A) / B);
2490  for (int i = 0; i < ncomp; i++) {
2491  if (x_ions == 0) {
2492  PCSAFT.SatL->mole_fractions[i] = pp[i] / pp_sum;
2493  PCSAFT.SatV->mole_fractions[i] = std::exp(u[i]) * pp[i] / eupp_sum;
2494  }
2495  else if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2496  PCSAFT.SatL->mole_fractions[i] = pp[i] / pp_sum * (1 - x_ions/(1 - PCSAFT._Q));
2497  PCSAFT.SatV->mole_fractions[i] = std::exp(u[i]) * pp[i] / eupp_sum;
2498  }
2499  else {
2500  PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 - PCSAFT._Q);
2501  PCSAFT.SatV->mole_fractions[i] = 0;
2502  }
2503  }
2504 
2505  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, p, iphase_liquid);
2506  vector<CoolPropDbl> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2507  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, p, iphase_gas);
2508  vector<CoolPropDbl> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2509  for (int i = 0; i < ncomp; i++) {
2510  k[i] = fugcoef_l[i] / fugcoef_v[i];
2511  u[i] = std::log(k[i] / kb);
2512  }
2513 
2514  if (itr == 0) {
2515  B = std::log(kbprime / kb) / (1/Pprime - 1/p);
2516  }
2517  A = std::log(kb) - B * (1/p - 1/Pref);
2518 
2519  maxdif = std::abs(A - A_old);
2520  for (int i = 0; i < ncomp; i++) {
2521  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2522  double dif = std::abs(u[i] - u_old[i]);
2523  if (dif > maxdif) {
2524  maxdif = dif;
2525  } else if (!std::isfinite(dif)) {
2526  maxdif = dif;
2527  }
2528  }
2529  }
2530  itr += 1;
2531  }
2532 
2533  if (!std::isfinite(p) || !std::isfinite(maxdif) || maxdif > 0.1 || p < 0) {
2534  throw SolutionError("outerTQ did not converge to a solution");
2535  }
2536 
2537  return p;
2538 }
2539 
2540 double PCSAFTBackend::estimate_flash_t(PCSAFTBackend &PCSAFT) {
2544  double t_guess = _HUGE;
2545  int ncomp = N; // number of components
2546 
2547  double x_ions = 0.; // overall mole fraction of ions in the system
2548  for (int i = 0; i < ncomp; i++) {
2549  if (PCSAFT.ion_term && PCSAFT.components[i].getZ() != 0) {
2550  x_ions += PCSAFT.mole_fractions[i];
2551  }
2552  }
2553 
2554  bool guess_found = false;
2555  double t_step = 30;
2556  double t_start = 571;
2557  double t_lbound = 1;
2558  if (PCSAFT.ion_term) {
2559  t_step = 15;
2560  t_start = 350;
2561  t_lbound = 264;
2562  }
2563  while (!guess_found && t_start > t_lbound) {
2564  // initialize variables
2565  double Tprime = t_start - 50;
2566  double t = t_start;
2567 
2568  PCSAFT.SatL->_T = t; // _T must be updated because the density calculation depends on it
2569  PCSAFT.SatV->_T = t;
2570 
2571  // calculate sigma for water, if it is present
2572  if (PCSAFT.water_present) {
2573  PCSAFT.components[water_idx].calc_water_sigma(t);
2574  PCSAFT.SatL->components[water_idx].calc_water_sigma(t);
2575  PCSAFT.SatV->components[water_idx].calc_water_sigma(t);
2576  PCSAFT.dielc = dielc_water(t); // Right now only aqueous mixtures are supported. Other solvents could be modeled by replacing the dielc_water function.
2577  PCSAFT.SatL->dielc = dielc_water(t);
2578  PCSAFT.SatV->dielc = dielc_water(t);
2579  }
2580 
2581  try {
2582  double p1 = estimate_flash_p(PCSAFT);
2583  PCSAFT.SatL->_T = Tprime;
2584  PCSAFT.SatV->_T = Tprime;
2585  double p2 = estimate_flash_p(PCSAFT);
2586  PCSAFT.SatL->_T = t; // reset to initial value
2587  PCSAFT.SatV->_T = t;
2588 
2589  double slope = (std::log10(p1) - std::log10(p2)) / (1/t - 1/Tprime);
2590  double intercept = std::log10(p1) - slope * (1/t);
2591  t_guess = slope / (std::log10(PCSAFT._p) - intercept);
2592  guess_found = true;
2593  } catch (const SolutionError& ex) {
2594  t_start -= t_step;
2595  }
2596  }
2597 
2598  if (!guess_found) {
2599  throw SolutionError("an estimate for the VLE temperature could not be found");
2600  }
2601 
2602  return t_guess;
2603 }
2604 
2605 
2606 double PCSAFTBackend::estimate_flash_p(PCSAFTBackend &PCSAFT) {
2610  double p_guess = _HUGE;
2611  int ncomp = N; // number of components
2612 
2613  double x_ions = 0.; // overall mole fraction of ions in the system
2614  for (int i = 0; i < ncomp; i++) {
2615  if (PCSAFT.ion_term && PCSAFT.components[i].getZ() != 0) {
2616  x_ions += PCSAFT.mole_fractions[i];
2617  }
2618  }
2619 
2620  bool guess_found = false;
2621  double p_start = 10000;
2622  while (!guess_found && p_start < 1e7) {
2623  // initialize variables
2624  vector<double> k(ncomp, 0), u(ncomp, 0), kprime(ncomp, 0), uprime(ncomp, 0);
2625  double Pprime = 0.99 * p_start;
2626  double p = p_start;
2627 
2628  // calculate initial guess for compositions based on fugacity coefficients and Raoult's Law.
2629  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT._T, p, iphase_liquid);
2630  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT._T, p, iphase_gas);
2631  if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2632  p_start = p_start + 2e5;
2633  continue;
2634  }
2635  vector<double> fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2636  vector<double> fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2637 
2638 
2639  double xv_sum = 0;
2640  double xl_sum = 0;
2641  for (int i = 0; i < ncomp; i++) {
2642  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2643  k[i] = fugcoef_l[i] / fugcoef_v[i];
2644  } else {
2645  k[i] = 0; // set k to 0 for ionic components
2646  }
2647  PCSAFT.SatL->mole_fractions[i] = PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2648  xl_sum += PCSAFT.SatL->mole_fractions[i];
2649  PCSAFT.SatV->mole_fractions[i] = k[i] * PCSAFT.mole_fractions[i] / (1 + PCSAFT._Q * (k[i] - 1));
2650  xv_sum += PCSAFT.SatV->mole_fractions[i];
2651  }
2652 
2653  if (xv_sum != 1) {
2654  for (int i = 0; i < ncomp; i++) {
2655  PCSAFT.SatV->mole_fractions[i] = PCSAFT.SatV->mole_fractions[i] / xv_sum;
2656  }
2657  }
2658 
2659  if (xl_sum != 1) {
2660  for (int i = 0; i < ncomp; i++) {
2661  PCSAFT.SatL->mole_fractions[i] = PCSAFT.SatL->mole_fractions[i] / xl_sum;
2662  }
2663  }
2664 
2665  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT.SatL->_T, p, iphase_liquid);
2666  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT.SatV->_T, p, iphase_gas);
2667  if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2668  p_start = p_start + 2e5;
2669  continue;
2670  }
2671  fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2672  fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2673  double numer = 0;
2674  double denom = 0;
2675  for (int i = 0; i < ncomp; i++) {
2676  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2677  numer += PCSAFT.SatL->mole_fractions[i] * fugcoef_l[i];
2678  denom += PCSAFT.SatV->mole_fractions[i] * fugcoef_v[i];
2679  }
2680  }
2681  double ratio = numer / denom;
2682 
2683  PCSAFT.SatL->_rhomolar = PCSAFT.SatL->solver_rho_Tp(PCSAFT.SatL->_T, Pprime, iphase_liquid);
2684  PCSAFT.SatV->_rhomolar = PCSAFT.SatV->solver_rho_Tp(PCSAFT.SatV->_T, Pprime, iphase_gas);
2685  if ((PCSAFT.SatL->_rhomolar - PCSAFT.SatV->_rhomolar) < 1e-4) {
2686  p_start = p_start + 2e5;
2687  continue;
2688  }
2689  fugcoef_l = PCSAFT.SatL->calc_fugacity_coefficients();
2690  fugcoef_v = PCSAFT.SatV->calc_fugacity_coefficients();
2691  numer = 0;
2692  denom = 0;
2693  for (int i = 0; i < ncomp; i++) {
2694  if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
2695  numer += PCSAFT.SatL->mole_fractions[i] * fugcoef_l[i];
2696  denom += PCSAFT.SatV->mole_fractions[i] * fugcoef_v[i];
2697  }
2698  }
2699  double ratio_prime = numer / denom;
2700 
2701  double slope = (std::log10(ratio) - std::log10(ratio_prime)) / (std::log10(p) - std::log10(Pprime));
2702  double intercept = std::log10(ratio) - slope * std::log10(p);
2703  p_guess = std::pow(10, -intercept / slope);
2704 
2705  guess_found = true;
2706  }
2707 
2708  if (!guess_found) {
2709  throw SolutionError("an estimate for the VLE pressure could not be found");
2710  }
2711 
2712  return p_guess;
2713 }
2714 
2715 CoolPropDbl PCSAFTBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, phases phase) {
2716  // Define the residual to be driven to zero
2717  class SolverRhoResid : public FuncWrapper1D
2718  {
2719  public:
2720  PCSAFTBackend& PCSAFT;
2721  CoolPropDbl T, p;
2722 
2723  SolverRhoResid(PCSAFTBackend& PCSAFT, CoolPropDbl T, CoolPropDbl p) : PCSAFT(PCSAFT), T(T), p(p) {}
2724  CoolPropDbl call(CoolPropDbl rhomolar) {
2725  CoolPropDbl peos = PCSAFT.update_DmolarT(rhomolar);
2726  double cost = (peos - p) / p;
2727  if (ValidNumber(cost)) {
2728  return cost;
2729  } else {
2730  return 1.0e20;
2731  }
2732  };
2733  };
2734 
2735  SolverRhoResid resid(*this, T, p);
2736 
2737  // split into grid and find bounds for each root
2738  vector<double> x_lo, x_hi;
2739  int num_pts = 20;
2740  double limit_lower = -8; // first use a log scale for the low density region
2741  double limit_upper = -1;
2742  double rho_guess = 1e-13;
2743  double rho_guess_prev = rho_guess;
2744  double err_prev = (update_DmolarT(reduced_to_molar(rho_guess, T)) - p) / p;
2745  for (int i = 0; i < num_pts; i++) {
2746  rho_guess = pow(10, (limit_upper - limit_lower) / (double)num_pts * i + limit_lower);
2747  double err = (update_DmolarT(reduced_to_molar(rho_guess, T)) - p) / p;
2748  if (err * err_prev < 0) {
2749  x_lo.push_back(rho_guess_prev);
2750  x_hi.push_back(rho_guess);
2751  }
2752  err_prev = err;
2753  rho_guess_prev = rho_guess;
2754  }
2755 
2756  limit_lower = 0.1; // for the high density region the log scale is not needed
2757  limit_upper = 0.7405;
2758  for (int i = 0; i < num_pts; i++) {
2759  rho_guess = (limit_upper - limit_lower) / (double)num_pts * i + limit_lower;
2760  double err = (update_DmolarT(reduced_to_molar(rho_guess, T)) - p) / p;
2761  if (err * err_prev < 0) {
2762  x_lo.push_back(rho_guess_prev);
2763  x_hi.push_back(rho_guess);
2764  }
2765  err_prev = err;
2766  rho_guess_prev = rho_guess;
2767  }
2768 
2769  // solve for appropriate root(s)
2770  double rho = _HUGE;
2771  double x_lo_molar = 1e-8, x_hi_molar = 1e7;
2772 
2773  if (x_lo.size() == 1) {
2774  rho_guess = reduced_to_molar((x_lo[0] + x_hi[0]) / 2., T);
2775  x_lo_molar = reduced_to_molar(x_lo[0], T);
2776  x_hi_molar = reduced_to_molar(x_hi[0], T);
2777  rho = Brent(resid, x_lo_molar, x_hi_molar, DBL_EPSILON, 1e-8, 200);
2778  } else if (x_lo.size() <= 3 && !x_lo.empty()) {
2779  if ((phase == iphase_liquid) || (phase == iphase_supercritical_liquid)) {
2780  rho_guess = reduced_to_molar((x_lo.back() + x_hi.back()) / 2., T);
2781  x_lo_molar = reduced_to_molar(x_lo.back(), T);
2782  x_hi_molar = reduced_to_molar(x_hi.back(), T);
2783  rho = Brent(resid, x_lo_molar, x_hi_molar, DBL_EPSILON, 1e-8, 200);
2784  } else if ((phase == iphase_gas) || (phase == iphase_supercritical_gas) || (phase == iphase_supercritical)) {
2785  rho_guess = reduced_to_molar((x_lo[0] + x_hi[0]) / 40., T); // starting with a lower guess often provides better results
2786  x_lo_molar = reduced_to_molar(x_lo[0], T);
2787  x_hi_molar = reduced_to_molar(x_hi[0], T);
2788  rho = Brent(resid, x_lo_molar, x_hi_molar, DBL_EPSILON, 1e-8, 200);
2789  }
2790  } else if (x_lo.size() > 3) {
2791  // 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.
2792  double g_min = 1e60;
2793  for (int i = 0; i < x_lo.size(); i++) {
2794  rho_guess = reduced_to_molar((x_lo[i] + x_hi[i]) / 2., T);
2795  x_lo_molar = reduced_to_molar(x_lo[i], T);
2796  x_hi_molar = reduced_to_molar(x_hi[i], T);
2797  double rho_i = Brent(resid, x_lo_molar, x_hi_molar, DBL_EPSILON, 1e-8, 200);
2798  double rho_original = this->_rhomolar;
2799  this->_rhomolar = rho_i;
2800  double g_i = calc_gibbsmolar_residual();
2801  this->_rhomolar = rho_original;
2802  if (g_i < g_min) {
2803  g_min = g_i;
2804  rho = rho_i;
2805  }
2806  }
2807  } else {
2808  int num_pts = 25;
2809  double err_min = 1e40;
2810  double rho_min;
2811  for (int i = 0; i < num_pts; i++) {
2812  double rho_guess = (0.7405 - 1e-8) / (double)num_pts * i + 1e-8;
2813  double err = (update_DmolarT(reduced_to_molar(rho_guess, T)) - p) / p;
2814  if (abs(err) < err_min) {
2815  err_min = abs(err);
2816  rho_min = reduced_to_molar(rho_guess, T);
2817  }
2818  }
2819  rho = rho_min;
2820  }
2821 
2822  return rho;
2823 }
2824 
2825 CoolPropDbl PCSAFTBackend::reduced_to_molar(CoolPropDbl nu, CoolPropDbl T) {
2826  vector<CoolPropDbl> d(N);
2827  CoolPropDbl summ = 0.;
2828  for (int i = 0; i < N; i++) {
2829  d[i] = components[i].getSigma() * (1 - 0.12 * exp(-3 * components[i].getU() / T));
2830  summ += mole_fractions[i] * components[i].getM() * pow(d[i], 3.);
2831  }
2832  return 6 / PI * nu / summ * 1.0e30 / N_AV;
2833 }
2834 
2835 CoolPropDbl PCSAFTBackend::calc_molar_mass(void) {
2836  double summer = 0;
2837  for (unsigned int i = 0; i < N; ++i) {
2838  summer += mole_fractions[i] * components[i].molar_mass();
2839  }
2840  return summer;
2841 }
2842 
2843 
2844 vector<double> PCSAFTBackend::XA_find(vector<double> XA_guess, vector<double> delta_ij, double den,
2845  vector<double> x) {
2847  int num_sites = XA_guess.size();
2848  vector<double> XA = XA_guess;
2849 
2850  int idxij = -1; // index for delta_ij
2851  for (int i = 0; i < num_sites; i++) {
2852  double summ = 0.;
2853  for (int j = 0; j < num_sites; j++) {
2854  idxij += 1;
2855  summ += den*x[j]*XA_guess[j]*delta_ij[idxij];
2856  }
2857  XA[i] = 1./(1.+summ);
2858  }
2859 
2860  return XA;
2861 }
2862 
2863 
2864 vector<double> PCSAFTBackend::dXAdt_find(vector<double> delta_ij, double den,
2865  vector<double> XA, vector<double> ddelta_dt, vector<double> x) {
2867  int num_sites = XA.size();
2868  Eigen::MatrixXd B = Eigen::MatrixXd::Zero(num_sites, 1);
2869  Eigen::MatrixXd A = Eigen::MatrixXd::Zero(num_sites, num_sites);
2870 
2871  double summ;
2872  int ij = 0;
2873  for (int i = 0; i < num_sites; i++) {
2874  summ = 0;
2875  for (int j = 0; j < num_sites; j++) {
2876  B(i) -= x[j]*XA[j]*ddelta_dt[ij];
2877  A(i,j) = x[j]*delta_ij[ij];
2878  summ += x[j]*XA[j]*delta_ij[ij];
2879  ij += 1;
2880  }
2881  A(i,i) = pow(1+den*summ, 2.)/den;
2882  }
2883 
2884  Eigen::MatrixXd solution = A.lu().solve(B); //Solves linear system of equations
2885  vector<double> dXA_dt(num_sites);
2886  for (int i = 0; i < num_sites; i++) {
2887  dXA_dt[i] = solution(i);
2888  }
2889  return dXA_dt;
2890 }
2891 
2892 
2893 vector<double> PCSAFTBackend::dXAdx_find(vector<int> assoc_num, vector<double> delta_ij,
2894  double den, vector<double> XA, vector<double> ddelta_dx, vector<double> x) {
2897  int num_sites = XA.size();
2898  int ncomp = assoc_num.size();
2899  Eigen::MatrixXd B(num_sites*ncomp, 1);
2900  Eigen::MatrixXd A = Eigen::MatrixXd::Zero(num_sites*ncomp, num_sites*ncomp);
2901 
2902  double sum1, sum2;
2903  int idx1 = 0;
2904  int ij = 0;
2905  for (int i = 0; i < ncomp; i++) {
2906  for (int j = 0; j < num_sites; j++) {
2907  sum1 = 0;
2908  for (int k = 0; k < num_sites; k++) {
2909  sum1 = sum1 + den*x[k]*(XA[k]*ddelta_dx[i*num_sites*num_sites + j*num_sites + k]);
2910  A(ij,i*num_sites+k) = XA[j]*XA[j]*den*x[k]*delta_ij[j*num_sites+k];
2911  }
2912 
2913  sum2 = 0;
2914  for (int l = 0; l < assoc_num[i]; l++) {
2915  sum2 = sum2 + XA[idx1+l]*delta_ij[idx1*num_sites+l*num_sites+j];
2916  }
2917 
2918  A(ij,ij) = A(ij,ij) + 1;
2919  B(ij) = -1*XA[j]*XA[j]*(sum1 + sum2);
2920  ij += 1;
2921  }
2922  idx1 += assoc_num[i];
2923  }
2924 
2925  Eigen::MatrixXd solution = A.lu().solve(B); //Solves linear system of equations
2926  vector<double> dXA_dx(num_sites*ncomp);
2927  for (int i = 0; i < num_sites*ncomp; i++) {
2928  dXA_dx[i] = solution(i);
2929  }
2930  return dXA_dx;
2931 }
2932 
2933 
2934 void PCSAFTBackend::set_assoc_matrix(){
2935  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)
2936 
2937  for (int i = 0; i < N; i++){
2938  vector<std::string> assoc_scheme = components[i].getAssocScheme();
2939  int num_sites = 0;
2940  int num = assoc_scheme.size();
2941  for (int j = 0; j < num; j++) {
2942  switch(get_scheme_index(assoc_scheme[j])) {
2943  case i1: {
2944  charge.push_back(0);
2945  num_sites += 1;
2946  break;
2947  }
2948  case i2a: {
2949  vector<int> tmp{0, 0};
2950  charge.insert(charge.end(), tmp.begin(), tmp.end());
2951  num_sites += 2;
2952  break;
2953  }
2954  case i2b: {
2955  vector<int> tmp{-1, 1};
2956  charge.insert(charge.end(), tmp.begin(), tmp.end());
2957  num_sites += 2;
2958  break;
2959  }
2960  case i3a: {
2961  vector<int> tmp{0, 0, 0};
2962  charge.insert(charge.end(), tmp.begin(), tmp.end());
2963  num_sites += 3;
2964  break;
2965  }
2966  case i3b: {
2967  vector<int> tmp{-1, -1, 1};
2968  charge.insert(charge.end(), tmp.begin(), tmp.end());
2969  num_sites += 3;
2970  break;
2971  }
2972  case i4a: {
2973  vector<int> tmp{0, 0, 0, 0};
2974  charge.insert(charge.end(), tmp.begin(), tmp.end());
2975  num_sites += 4;
2976  break;
2977  }
2978  case i4b: {
2979  vector<int> tmp{1, 1, 1, -1};
2980  charge.insert(charge.end(), tmp.begin(), tmp.end());
2981  num_sites += 4;
2982  break;
2983  }
2984  case i4c: {
2985  vector<int> tmp{-1, -1, 1, 1};
2986  charge.insert(charge.end(), tmp.begin(), tmp.end());
2987  num_sites += 4;
2988  break;
2989  }
2990  default:
2991  throw ValueError(format("%s is not a valid association type.", assoc_scheme[j]));
2992  }
2993  }
2994 
2995  assoc_num.push_back(num_sites);
2996  }
2997 
2998  for (std::vector<int>::iterator i1 = charge.begin(); i1 != charge.end(); i1++) {
2999  for (std::vector<int>::iterator i2 = charge.begin(); i2 != charge.end(); i2++) {
3000  if (*i1 == 0 || *i2 == 0) {
3001  assoc_matrix.push_back(1);
3002  }
3003  else if (*i1 == 1 && *i2 == -1) {
3004  assoc_matrix.push_back(1);
3005  }
3006  else if (*i1 == -1 && *i2 == 1) {
3007  assoc_matrix.push_back(1);
3008  }
3009  else {
3010  assoc_matrix.push_back(0);
3011  }
3012  }
3013  }
3014 }
3015 
3016 
3017 double PCSAFTBackend::dielc_water(double t) {
3033  double dielc;
3034  if (t < 263.15) {
3035  throw ValueError("The current function for the dielectric constant for water is only valid for temperatures above 263.15 K.");
3036  } else if (t <= 368.15) {
3037  dielc = 7.6555618295E-04 * _T * _T - 8.1783881423E-01 * _T + 2.5419616803E+02;
3038  } else if (t <= 443.15) {
3039  dielc = 0.0005003272124 * _T * _T - 0.6285556029 * _T + 220.4467027;
3040  } else {
3041  throw ValueError("The current function for the dielectric constant for water is only valid for temperatures less than 443.15 K.");
3042  }
3043  return dielc;
3044 }
3045 
3046 } /* namespace CoolProp */
Generated on Sat Feb 22 2025 01:20:23 for CoolProp by doxygen 1.9.1