Cantera  3.1.0a1
RedlichKisterVPSSTP.cpp
Go to the documentation of this file.
1 /**
2  * @file RedlichKisterVPSSTP.cpp
3  * Definitions for ThermoPhase object for phases which
4  * employ excess Gibbs free energy formulations related to RedlichKister
5  * expansions (see @ref thermoprops
6  * and class @link Cantera::RedlichKisterVPSSTP RedlichKisterVPSSTP@endlink).
7  */
8 
9 // This file is part of Cantera. See License.txt in the top-level directory or
10 // at https://cantera.org/license.txt for license and copyright information.
11 
12 #include "cantera/thermo/RedlichKisterVPSSTP.h"
13 #include "cantera/thermo/ThermoFactory.h"
14 #include "cantera/base/stringUtils.h"
15 
16 using namespace std;
17 
18 namespace Cantera
19 {
20 RedlichKisterVPSSTP::RedlichKisterVPSSTP(const string& inputFile, const string& id_)
21 {
22  initThermoFile(inputFile, id_);
23 }
24 
25 // - Activities, Standard States, Activity Concentrations -----------
26 
27 void RedlichKisterVPSSTP::getLnActivityCoefficients(double* lnac) const
28 {
29  // Update the activity coefficients
30  s_update_lnActCoeff();
31 
32  for (size_t k = 0; k < m_kk; k++) {
33  lnac[k] = lnActCoeff_Scaled_[k];
34  }
35 }
36 
37 // ------------ Partial Molar Properties of the Solution ------------
38 
39 void RedlichKisterVPSSTP::getChemPotentials(double* mu) const
40 {
41  // First get the standard chemical potentials in molar form. This requires
42  // updates of standard state as a function of T and P
43  getStandardChemPotentials(mu);
44  // Update the activity coefficients
45  s_update_lnActCoeff();
46 
47  for (size_t k = 0; k < m_kk; k++) {
48  double xx = std::max(moleFractions_[k], SmallNumber);
49  mu[k] += RT() * (log(xx) + lnActCoeff_Scaled_[k]);
50  }
51 }
52 
53 double RedlichKisterVPSSTP::enthalpy_mole() const
54 {
55  double h = 0;
56  vector<double> hbar(m_kk);
57  getPartialMolarEnthalpies(&hbar[0]);
58  for (size_t i = 0; i < m_kk; i++) {
59  h += moleFractions_[i]*hbar[i];
60  }
61  return h;
62 }
63 
64 double RedlichKisterVPSSTP::entropy_mole() const
65 {
66  double s = 0;
67  vector<double> sbar(m_kk);
68  getPartialMolarEntropies(&sbar[0]);
69  for (size_t i = 0; i < m_kk; i++) {
70  s += moleFractions_[i]*sbar[i];
71  }
72  return s;
73 }
74 
75 double RedlichKisterVPSSTP::cp_mole() const
76 {
77  double cp = 0;
78  vector<double> cpbar(m_kk);
79  getPartialMolarCp(&cpbar[0]);
80  for (size_t i = 0; i < m_kk; i++) {
81  cp += moleFractions_[i]*cpbar[i];
82  }
83  return cp;
84 }
85 
86 double RedlichKisterVPSSTP::cv_mole() const
87 {
88  return cp_mole();
89 }
90 
91 void RedlichKisterVPSSTP::getPartialMolarEnthalpies(double* hbar) const
92 {
93  // Get the nondimensional standard state enthalpies
94  getEnthalpy_RT(hbar);
95  // dimensionalize it.
96  double T = temperature();
97  for (size_t k = 0; k < m_kk; k++) {
98  hbar[k] *= GasConstant * T;
99  }
100 
101  // Update the activity coefficients, This also update the internally stored
102  // molalities.
103  s_update_lnActCoeff();
104  s_update_dlnActCoeff_dT();
105  for (size_t k = 0; k < m_kk; k++) {
106  hbar[k] -= GasConstant * T * T * dlnActCoeffdT_Scaled_[k];
107  }
108 }
109 
110 void RedlichKisterVPSSTP::getPartialMolarCp(double* cpbar) const
111 {
112  getCp_R(cpbar);
113  // dimensionalize it.
114  for (size_t k = 0; k < m_kk; k++) {
115  cpbar[k] *= GasConstant;
116  }
117 }
118 
119 void RedlichKisterVPSSTP::getPartialMolarEntropies(double* sbar) const
120 {
121  // Get the nondimensional standard state entropies
122  getEntropy_R(sbar);
123  double T = temperature();
124 
125  // Update the activity coefficients, This also update the internally stored
126  // molalities.
127  s_update_lnActCoeff();
128  s_update_dlnActCoeff_dT();
129 
130  for (size_t k = 0; k < m_kk; k++) {
131  double xx = std::max(moleFractions_[k], SmallNumber);
132  sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
133  }
134  // dimensionalize it.
135  for (size_t k = 0; k < m_kk; k++) {
136  sbar[k] *= GasConstant;
137  }
138 }
139 
140 void RedlichKisterVPSSTP::getPartialMolarVolumes(double* vbar) const
141 {
142  // Get the standard state values in m^3 kmol-1
143  getStandardVolumes(vbar);
144  for (size_t iK = 0; iK < m_kk; iK++) {
145  vbar[iK] += 0.0;
146  }
147 }
148 
149 void RedlichKisterVPSSTP::initThermo()
150 {
151  if (m_input.hasKey("interactions")) {
152  for (const auto& item : m_input["interactions"].asVector<AnyMap>()) {
153  auto& species = item["species"].asVector<string>(2);
154  vector<double> h_excess = item.convertVector("excess-enthalpy", "J/kmol");
155  vector<double> s_excess = item.convertVector("excess-entropy", "J/kmol/K");
156  addBinaryInteraction(species[0], species[1],
157  h_excess.data(), h_excess.size(),
158  s_excess.data(), s_excess.size());
159  }
160  }
161  initLengths();
162  GibbsExcessVPSSTP::initThermo();
163 }
164 
165 void RedlichKisterVPSSTP::getParameters(AnyMap& phaseNode) const
166 {
167  GibbsExcessVPSSTP::getParameters(phaseNode);
168  vector<AnyMap> interactions;
169  for (size_t n = 0; n < m_pSpecies_A_ij.size(); n++) {
170  AnyMap interaction;
171  interaction["species"] = vector<string>{
172  speciesName(m_pSpecies_A_ij[n]), speciesName(m_pSpecies_B_ij[n])};
173  vector<double> h = m_HE_m_ij[n];
174  vector<double> s = m_SE_m_ij[n];
175  while (h.size() > 1 && h.back() == 0) {
176  h.pop_back();
177  }
178  while (s.size() > 1 && s.back() == 0) {
179  s.pop_back();
180  }
181  interaction["excess-enthalpy"].setQuantity(std::move(h), "J/kmol");
182  interaction["excess-entropy"].setQuantity(std::move(s), "J/kmol/K");
183  interactions.push_back(std::move(interaction));
184  }
185  phaseNode["interactions"] = std::move(interactions);
186 }
187 
188 void RedlichKisterVPSSTP::initLengths()
189 {
190  dlnActCoeffdlnN_.resize(m_kk, m_kk);
191 }
192 
193 void RedlichKisterVPSSTP::s_update_lnActCoeff() const
194 {
195  double T = temperature();
196  lnActCoeff_Scaled_.assign(m_kk, 0.0);
197 
198  // Scaling: I moved the division of RT higher so that we are always dealing
199  // with G/RT dimensionless terms within the routine. There is a severe
200  // problem with roundoff error in these calculations. The dimensionless
201  // terms help.
202  for (size_t i = 0; i < m_HE_m_ij.size(); i++) {
203  size_t iA = m_pSpecies_A_ij[i];
204  size_t iB = m_pSpecies_B_ij[i];
205  double XA = moleFractions_[iA];
206  double XB = moleFractions_[iB];
207  double deltaX = XA - XB;
208  const vector<double>& he_vec = m_HE_m_ij[i];
209  const vector<double>& se_vec = m_SE_m_ij[i];
210  double poly = 1.0;
211  double polyMm1 = 1.0;
212  double sum = 0.0;
213  double sumMm1 = 0.0;
214  double sum2 = 0.0;
215  for (size_t m = 0; m < he_vec.size(); m++) {
216  double A_ge = (he_vec[m] - T * se_vec[m]) / (GasConstant * T);
217  sum += A_ge * poly;
218  sum2 += A_ge * (m + 1) * poly;
219  poly *= deltaX;
220  if (m >= 1) {
221  sumMm1 += (A_ge * polyMm1 * m);
222  polyMm1 *= deltaX;
223  }
224  }
225  double oneMXA = 1.0 - XA;
226  double oneMXB = 1.0 - XB;
227  for (size_t k = 0; k < m_kk; k++) {
228  if (iA == k) {
229  lnActCoeff_Scaled_[k] += (oneMXA * XB * sum) + (XA * XB * sumMm1 * (oneMXA + XB));
230  } else if (iB == k) {
231  lnActCoeff_Scaled_[k] += (oneMXB * XA * sum) + (XA * XB * sumMm1 * (-oneMXB - XA));
232  } else {
233  lnActCoeff_Scaled_[k] += -(XA * XB * sum2);
234  }
235  }
236  // Debug against formula in literature
237  }
238 }
239 
240 void RedlichKisterVPSSTP::s_update_dlnActCoeff_dT() const
241 {
242  double T = temperature();
243  dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
244 
245  for (size_t i = 0; i < m_HE_m_ij.size(); i++) {
246  size_t iA = m_pSpecies_A_ij[i];
247  size_t iB = m_pSpecies_B_ij[i];
248  double XA = moleFractions_[iA];
249  double XB = moleFractions_[iB];
250  double deltaX = XA - XB;
251  double poly = 1.0;
252  double sum = 0.0;
253  const vector<double>& he_vec = m_HE_m_ij[i];
254  double sumMm1 = 0.0;
255  double polyMm1 = 1.0;
256  double sum2 = 0.0;
257  for (size_t m = 0; m < he_vec.size(); m++) {
258  double h_e = - he_vec[m] / (GasConstant * T * T);
259  sum += h_e * poly;
260  sum2 += h_e * (m + 1) * poly;
261  poly *= deltaX;
262  if (m >= 1) {
263  sumMm1 += (h_e * polyMm1 * m);
264  polyMm1 *= deltaX;
265  }
266  }
267  double oneMXA = 1.0 - XA;
268  double oneMXB = 1.0 - XB;
269  for (size_t k = 0; k < m_kk; k++) {
270  if (iA == k) {
271  dlnActCoeffdT_Scaled_[k] += (oneMXA * XB * sum) + (XA * XB * sumMm1 * (oneMXA + XB));
272  } else if (iB == k) {
273  dlnActCoeffdT_Scaled_[k] += (oneMXB * XA * sum) + (XA * XB * sumMm1 * (-oneMXB - XA));
274  } else {
275  dlnActCoeffdT_Scaled_[k] += -(XA * XB * sum2);
276  }
277  }
278  }
279 
280  for (size_t k = 0; k < m_kk; k++) {
281  d2lnActCoeffdT2_Scaled_[k] = -2 / T * dlnActCoeffdT_Scaled_[k];
282  }
283 }
284 
285 void RedlichKisterVPSSTP::getdlnActCoeffdT(double* dlnActCoeffdT) const
286 {
287  s_update_dlnActCoeff_dT();
288  for (size_t k = 0; k < m_kk; k++) {
289  dlnActCoeffdT[k] = dlnActCoeffdT_Scaled_[k];
290  }
291 }
292 
293 void RedlichKisterVPSSTP::getd2lnActCoeffdT2(double* d2lnActCoeffdT2) const
294 {
295  s_update_dlnActCoeff_dT();
296  for (size_t k = 0; k < m_kk; k++) {
297  d2lnActCoeffdT2[k] = d2lnActCoeffdT2_Scaled_[k];
298  }
299 }
300 
301 void RedlichKisterVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
302 {
303  double T = temperature();
304  dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
305 
306  for (size_t i = 0; i < m_HE_m_ij.size(); i++) {
307  size_t iA = m_pSpecies_A_ij[i];
308  size_t iB = m_pSpecies_B_ij[i];
309  double XA = moleFractions_[iA];
310  double XB = moleFractions_[iB];
311  double deltaX = XA - XB;
312  double poly = 1.0;
313  double sum = 0.0;
314  const vector<double>& he_vec = m_HE_m_ij[i];
315  const vector<double>& se_vec = m_SE_m_ij[i];
316  double sumMm1 = 0.0;
317  double polyMm1 = 1.0;
318  double polyMm2 = 1.0;
319  double sumMm2 = 0.0;
320  for (size_t m = 0; m < he_vec.size(); m++) {
321  double A_ge = (he_vec[m] - T * se_vec[m]) / (GasConstant * T);;
322  sum += A_ge * poly;
323  poly *= deltaX;
324  if (m >= 1) {
325  sumMm1 += (A_ge * polyMm1 * m);
326  polyMm1 *= deltaX;
327  }
328  if (m >= 2) {
329  sumMm2 += (A_ge * polyMm2 * m * (m - 1.0));
330  polyMm2 *= deltaX;
331  }
332  }
333 
334  for (size_t k = 0; k < m_kk; k++) {
335  if (iA == k) {
336  dlnActCoeffdlnX_diag_[k] +=
337  XA * (- (1-XA+XB) * sum + 2*(1.0 - XA) * XB * sumMm1
338  + sumMm1 * (XB * (1 - 2*XA + XB) - XA * (1 - XA + 2*XB))
339  + 2 * XA * XB * sumMm2 * (1.0 - XA + XB));
340  } else if (iB == k) {
341  dlnActCoeffdlnX_diag_[k] +=
342  XB * (- (1-XB+XA) * sum - 2*(1.0 - XB) * XA * sumMm1
343  + sumMm1 * (XA * (2*XB - XA - 1) - XB * (-2*XA + XB - 1))
344  - 2 * XA * XB * sumMm2 * (-XA - 1 + XB));
345  }
346  }
347  }
348 }
349 
350 void RedlichKisterVPSSTP::s_update_dlnActCoeff_dX_() const
351 {
352  double T = temperature();
353  dlnActCoeff_dX_.zero();
354 
355  for (size_t i = 0; i < m_HE_m_ij.size(); i++) {
356  size_t iA = m_pSpecies_A_ij[i];
357  size_t iB = m_pSpecies_B_ij[i];
358  double XA = moleFractions_[iA];
359  double XB = moleFractions_[iB];
360  double deltaX = XA - XB;
361  double poly = 1.0;
362  double sum = 0.0;
363  const vector<double>& he_vec = m_HE_m_ij[i];
364  const vector<double>& se_vec = m_SE_m_ij[i];
365  double sumMm1 = 0.0;
366  double polyMm1 = 1.0;
367  double polyMm2 = 1.0;
368  double sum2 = 0.0;
369  double sum2Mm1 = 0.0;
370  double sumMm2 = 0.0;
371  for (size_t m = 0; m < he_vec.size(); m++) {
372  double A_ge = he_vec[m] - T * se_vec[m];
373  sum += A_ge * poly;
374  sum2 += A_ge * (m + 1) * poly;
375  poly *= deltaX;
376  if (m >= 1) {
377  sumMm1 += (A_ge * polyMm1 * m);
378  sum2Mm1 += (A_ge * polyMm1 * m * (1.0 + m));
379  polyMm1 *= deltaX;
380  }
381  if (m >= 2) {
382  sumMm2 += (A_ge * polyMm2 * m * (m - 1.0));
383  polyMm2 *= deltaX;
384  }
385  }
386 
387  for (size_t k = 0; k < m_kk; k++) {
388  if (iA == k) {
389  dlnActCoeff_dX_(k, iA) += (- XB * sum + (1.0 - XA) * XB * sumMm1
390  + XB * sumMm1 * (1.0 - 2.0 * XA + XB)
391  + XA * XB * sumMm2 * (1.0 - XA + XB));
392 
393  dlnActCoeff_dX_(k, iB) += ((1.0 - XA) * sum - (1.0 - XA) * XB * sumMm1
394  + XA * sumMm1 * (1.0 + 2.0 * XB - XA)
395  - XA * XB * sumMm2 * (1.0 - XA + XB));
396  } else if (iB == k) {
397  dlnActCoeff_dX_(k, iA) += ((1.0 - XB) * sum + (1.0 - XA) * XB * sumMm1
398  + XB * sumMm1 * (1.0 - 2.0 * XA + XB)
399  + XA * XB * sumMm2 * (1.0 - XA + XB));
400 
401  dlnActCoeff_dX_(k, iB) += (- XA * sum - (1.0 - XB) * XA * sumMm1
402  + XA * sumMm1 * (XB - XA - (1.0 - XB))
403  - XA * XB * sumMm2 * (-XA - (1.0 - XB)));
404  } else {
405  dlnActCoeff_dX_(k, iA) += (- XB * sum2 - XA * XB * sum2Mm1);
406  dlnActCoeff_dX_(k, iB) += (- XA * sum2 + XA * XB * sum2Mm1);
407  }
408  }
409  }
410 }
411 
412 void RedlichKisterVPSSTP::getdlnActCoeffds(const double dTds, const double* const dXds,
413  double* dlnActCoeffds) const
414 {
415  s_update_dlnActCoeff_dT();
416  s_update_dlnActCoeff_dX_();
417  for (size_t k = 0; k < m_kk; k++) {
418  dlnActCoeffds[k] = dlnActCoeffdT_Scaled_[k] * dTds;
419  for (size_t j = 0; j < m_kk; j++) {
420  dlnActCoeffds[k] += dlnActCoeff_dX_(k, j) * dXds[j];
421  }
422  }
423 }
424 
425 void RedlichKisterVPSSTP::getdlnActCoeffdlnN_diag(double* dlnActCoeffdlnN_diag) const
426 {
427  s_update_dlnActCoeff_dX_();
428  for (size_t j = 0; j < m_kk; j++) {
429  dlnActCoeffdlnN_diag[j] = dlnActCoeff_dX_(j, j);
430  for (size_t k = 0; k < m_kk; k++) {
431  dlnActCoeffdlnN_diag[k] -= dlnActCoeff_dX_(j, k) * moleFractions_[k];
432  }
433  }
434 }
435 
436 void RedlichKisterVPSSTP::getdlnActCoeffdlnX_diag(double* dlnActCoeffdlnX_diag) const
437 {
438  s_update_dlnActCoeff_dlnX_diag();
439  for (size_t k = 0; k < m_kk; k++) {
440  dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
441  }
442 }
443 
444 void RedlichKisterVPSSTP::getdlnActCoeffdlnN(const size_t ld, double* dlnActCoeffdlnN)
445 {
446  s_update_dlnActCoeff_dX_();
447  double* data = & dlnActCoeffdlnN_(0,0);
448  for (size_t k = 0; k < m_kk; k++) {
449  for (size_t m = 0; m < m_kk; m++) {
450  dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m];
451  }
452  }
453 }
454 
455 void RedlichKisterVPSSTP::addBinaryInteraction(
456  const string& speciesA, const string& speciesB,
457  const double* excess_enthalpy, size_t n_enthalpy,
458  const double* excess_entropy, size_t n_entropy)
459 {
460  size_t kA = speciesIndex(speciesA);
461  size_t kB = speciesIndex(speciesB);
462  if (kA == npos) {
463  throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
464  "Species '{}' not present in phase", speciesA);
465  } else if (kB == npos) {
466  throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
467  "Species '{}' not present in phase", speciesB);
468  }
469  if (charge(kA) != 0) {
470  throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
471  "Species '{}' should be neutral", speciesA);
472  } else if (charge(kB) != 0) {
473  throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
474  "Species '{}' should be neutral", speciesB);
475  }
476 
477  m_pSpecies_A_ij.push_back(kA);
478  m_pSpecies_B_ij.push_back(kB);
479  m_HE_m_ij.emplace_back(excess_enthalpy, excess_enthalpy + n_enthalpy);
480  m_SE_m_ij.emplace_back(excess_entropy, excess_entropy + n_entropy);
481  size_t N = max(n_enthalpy, n_entropy);
482  m_HE_m_ij.back().resize(N, 0.0);
483  m_SE_m_ij.back().resize(N, 0.0);
484  dlnActCoeff_dX_.resize(N, N, 0.0);
485 }
486 
487 }
RedlichKisterVPSSTP.h
(see Thermodynamic Properties and class RedlichKisterVPSSTP).
ThermoFactory.h
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
Cantera::AnyMap
A map of string keys to values whose type can vary at runtime.
Definition: AnyMap.h:427
Cantera::CanteraError
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
Cantera::GasConstant
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:120
Cantera
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
Cantera::npos
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:180
Cantera::SmallNumber
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:158
stringUtils.h
Contains declarations for string manipulation functions within Cantera.