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