Cantera  3.1.0a1
Loading...
Searching...
No Matches
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
15
16using namespace std;
17
18namespace Cantera
19{
20RedlichKisterVPSSTP::RedlichKisterVPSSTP(const string& inputFile, const string& id_)
21{
22 initThermoFile(inputFile, id_);
23}
24
25// - Activities, Standard States, Activity Concentrations -----------
26
28{
29 // Update the activity coefficients
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
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
44 // Update the activity coefficients
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
54{
55 double h = 0;
56 vector<double> hbar(m_kk);
58 for (size_t i = 0; i < m_kk; i++) {
59 h += moleFractions_[i]*hbar[i];
60 }
61 return h;
62}
63
65{
66 double s = 0;
67 vector<double> sbar(m_kk);
69 for (size_t i = 0; i < m_kk; i++) {
70 s += moleFractions_[i]*sbar[i];
71 }
72 return s;
73}
74
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
87{
88 return cp_mole();
89}
90
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.
105 for (size_t k = 0; k < m_kk; k++) {
106 hbar[k] -= GasConstant * T * T * dlnActCoeffdT_Scaled_[k];
107 }
108}
109
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
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.
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
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
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");
157 h_excess.data(), h_excess.size(),
158 s_excess.data(), s_excess.size());
159 }
160 }
161 initLengths();
163}
164
166{
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>{
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
189{
191}
192
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
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++) {
282 }
283}
284
285void RedlichKisterVPSSTP::getdlnActCoeffdT(double* dlnActCoeffdT) const
286{
288 for (size_t k = 0; k < m_kk; k++) {
289 dlnActCoeffdT[k] = dlnActCoeffdT_Scaled_[k];
290 }
291}
292
293void RedlichKisterVPSSTP::getd2lnActCoeffdT2(double* d2lnActCoeffdT2) const
294{
296 for (size_t k = 0; k < m_kk; k++) {
297 d2lnActCoeffdT2[k] = d2lnActCoeffdT2_Scaled_[k];
298 }
299}
300
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) {
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) {
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
351{
352 double T = temperature();
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
412void RedlichKisterVPSSTP::getdlnActCoeffds(const double dTds, const double* const dXds,
413 double* dlnActCoeffds) const
414{
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
425void RedlichKisterVPSSTP::getdlnActCoeffdlnN_diag(double* dlnActCoeffdlnN_diag) const
426{
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
436void RedlichKisterVPSSTP::getdlnActCoeffdlnX_diag(double* dlnActCoeffdlnX_diag) const
437{
439 for (size_t k = 0; k < m_kk; k++) {
440 dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
441 }
442}
443
444void RedlichKisterVPSSTP::getdlnActCoeffdlnN(const size_t ld, double* dlnActCoeffdlnN)
445{
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
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}
(see Thermodynamic Properties and class RedlichKisterVPSSTP).
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
void zero()
Set all of the entries to zero.
Definition Array.h:127
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition Array.cpp:47
Base class for exceptions thrown by Cantera classes.
vector< double > d2lnActCoeffdT2_Scaled_
Storage for the current derivative values of the gradients with respect to temperature of the log of ...
Array2D dlnActCoeffdlnN_
Storage for the current derivative values of the gradients with respect to logarithm of the species m...
vector< double > lnActCoeff_Scaled_
Storage for the current values of the activity coefficients of the species.
vector< double > dlnActCoeffdlnX_diag_
Storage for the current derivative values of the gradients with respect to logarithm of the mole frac...
vector< double > moleFractions_
Storage for the current values of the mole fractions of the species.
vector< double > dlnActCoeffdT_Scaled_
Storage for the current derivative values of the gradients with respect to temperature of the log of ...
size_t m_kk
Number of species in the phase.
Definition Phase.h:842
double temperature() const
Temperature (K).
Definition Phase.h:562
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:129
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition Phase.cpp:856
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition Phase.h:538
void getdlnActCoeffds(const double dTds, const double *const dXds, double *dlnActCoeffds) const override
Get the change in activity coefficients wrt changes in state (temp, mole fraction,...
void s_update_dlnActCoeff_dX_() const
Internal routine that calculates the derivative of the activity coefficients wrt the mole fractions.
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
Array2D dlnActCoeff_dX_
Two dimensional array of derivatives of activity coefficients wrt mole fractions.
vector< vector< double > > m_SE_m_ij
Entropy term for the binary mole fraction interaction of the excess Gibbs free energy expression.
void getParameters(AnyMap &phaseNode) const override
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
void getPartialMolarVolumes(double *vbar) const override
Return an array of partial molar volumes for the species in the mixture.
vector< size_t > m_pSpecies_A_ij
vector of species indices representing species A in the interaction
double cv_mole() const override
Molar heat capacity at constant volume. Units: J/kmol/K.
void s_update_dlnActCoeff_dT() const
Update the derivative of the log of the activity coefficients wrt T.
vector< size_t > m_pSpecies_B_ij
vector of species indices representing species B in the interaction
void addBinaryInteraction(const string &speciesA, const string &speciesB, const double *excess_enthalpy, size_t n_enthalpy, const double *excess_entropy, size_t n_entropy)
Add a binary species interaction with the specified parameters.
vector< vector< double > > m_HE_m_ij
Enthalpy term for the binary mole fraction interaction of the excess Gibbs free energy expression.
double entropy_mole() const override
Molar entropy. Units: J/kmol/K.
void getdlnActCoeffdT(double *dlnActCoeffdT) const override
Get the array of temperature derivatives of the log activity coefficients.
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
void getPartialMolarCp(double *cpbar) const override
Returns an array of partial molar heat capacities for the species in the mixture.
void initLengths()
Initialize lengths of local variables after all species have been identified.
void getLnActivityCoefficients(double *lnac) const override
Get the array of non-dimensional molar-based ln activity coefficients at the current solution tempera...
void s_update_dlnActCoeff_dlnX_diag() const
Internal routine that calculates the total derivative of the activity coefficients with respect to th...
void s_update_lnActCoeff() const
Update the activity coefficients.
void getdlnActCoeffdlnX_diag(double *dlnActCoeffdlnX_diag) const override
Get the array of ln mole fraction derivatives of the log activity coefficients - diagonal component o...
RedlichKisterVPSSTP(const string &inputFile="", const string &id="")
Construct a RedlichKisterVPSSTP object from an input file.
void getdlnActCoeffdlnN_diag(double *dlnActCoeffdlnN_diag) const override
Get the array of log species mole number derivatives of the log activity coefficients.
void getd2lnActCoeffdT2(double *d2lnActCoeffdT2) const
Get the array of temperature second derivatives of the log activity coefficients.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies for the species in the mixture.
void getdlnActCoeffdlnN(const size_t ld, double *const dlnActCoeffdlnN) override
Get the array of derivatives of the log activity coefficients with respect to the log of the species ...
virtual void getParameters(AnyMap &phaseNode) const
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
double RT() const
Return the Gas Constant multiplied by the current temperature.
virtual void initThermo()
Initialize the ThermoPhase object after all species have been set up.
void initThermoFile(const string &inputFile, const string &id)
Initialize a ThermoPhase object using an input file.
AnyMap m_input
Data supplied via setParameters.
void getEntropy_R(double *sr) const override
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
void getStandardChemPotentials(double *mu) const override
Get the array of chemical potentials at unit activity for the species at their standard states at the...
void getCp_R(double *cpr) const override
Get the nondimensional Heat Capacities at constant pressure for the species standard states at the cu...
void getEnthalpy_RT(double *hrt) const override
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
void getStandardVolumes(double *vol) const override
Get the molar volumes of the species standard states at the current T and P of the solution.
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition ct_defs.h:120
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
const double SmallNumber
smallest number to compare to zero.
Definition ct_defs.h:158
Contains declarations for string manipulation functions within Cantera.