Loading [MathJax]/extensions/tex2jax.js
Cantera  3.2.0a1
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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 return cp_mole();
56}
57
59{
60 // Get the nondimensional standard state enthalpies
61 getEnthalpy_RT(hbar);
62 // dimensionalize it.
63 double T = temperature();
64 for (size_t k = 0; k < m_kk; k++) {
65 hbar[k] *= GasConstant * T;
66 }
67
68 // Update the activity coefficients, This also update the internally stored
69 // molalities.
72 for (size_t k = 0; k < m_kk; k++) {
73 hbar[k] -= GasConstant * T * T * dlnActCoeffdT_Scaled_[k];
74 }
75}
76
78{
79 getCp_R(cpbar);
80 // dimensionalize it.
81 for (size_t k = 0; k < m_kk; k++) {
82 cpbar[k] *= GasConstant;
83 }
84}
85
87{
88 // Get the nondimensional standard state entropies
89 getEntropy_R(sbar);
90 double T = temperature();
91
92 // Update the activity coefficients, This also update the internally stored
93 // molalities.
96
97 for (size_t k = 0; k < m_kk; k++) {
98 double xx = std::max(moleFractions_[k], SmallNumber);
99 sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
100 }
101 // dimensionalize it.
102 for (size_t k = 0; k < m_kk; k++) {
103 sbar[k] *= GasConstant;
104 }
105}
106
108{
109 // Get the standard state values in m^3 kmol-1
110 getStandardVolumes(vbar);
111 for (size_t iK = 0; iK < m_kk; iK++) {
112 vbar[iK] += 0.0;
113 }
114}
115
117{
118 if (m_input.hasKey("interactions")) {
119 for (const auto& item : m_input["interactions"].asVector<AnyMap>()) {
120 auto& species = item["species"].asVector<string>(2);
121 vector<double> h_excess = item.convertVector("excess-enthalpy", "J/kmol");
122 vector<double> s_excess = item.convertVector("excess-entropy", "J/kmol/K");
124 h_excess.data(), h_excess.size(),
125 s_excess.data(), s_excess.size());
126 }
127 }
128 initLengths();
130}
131
133{
135 vector<AnyMap> interactions;
136 for (size_t n = 0; n < m_pSpecies_A_ij.size(); n++) {
137 AnyMap interaction;
138 interaction["species"] = vector<string>{
140 vector<double> h = m_HE_m_ij[n];
141 vector<double> s = m_SE_m_ij[n];
142 while (h.size() > 1 && h.back() == 0) {
143 h.pop_back();
144 }
145 while (s.size() > 1 && s.back() == 0) {
146 s.pop_back();
147 }
148 interaction["excess-enthalpy"].setQuantity(std::move(h), "J/kmol");
149 interaction["excess-entropy"].setQuantity(std::move(s), "J/kmol/K");
150 interactions.push_back(std::move(interaction));
151 }
152 phaseNode["interactions"] = std::move(interactions);
153}
154
156{
158}
159
161{
162 double T = temperature();
163 lnActCoeff_Scaled_.assign(m_kk, 0.0);
164
165 // Scaling: I moved the division of RT higher so that we are always dealing
166 // with G/RT dimensionless terms within the routine. There is a severe
167 // problem with roundoff error in these calculations. The dimensionless
168 // terms help.
169 for (size_t i = 0; i < m_HE_m_ij.size(); i++) {
170 size_t iA = m_pSpecies_A_ij[i];
171 size_t iB = m_pSpecies_B_ij[i];
172 double XA = moleFractions_[iA];
173 double XB = moleFractions_[iB];
174 double deltaX = XA - XB;
175 const vector<double>& he_vec = m_HE_m_ij[i];
176 const vector<double>& se_vec = m_SE_m_ij[i];
177 double poly = 1.0;
178 double polyMm1 = 1.0;
179 double sum = 0.0;
180 double sumMm1 = 0.0;
181 double sum2 = 0.0;
182 for (size_t m = 0; m < he_vec.size(); m++) {
183 double A_ge = (he_vec[m] - T * se_vec[m]) / (GasConstant * T);
184 sum += A_ge * poly;
185 sum2 += A_ge * (m + 1) * poly;
186 poly *= deltaX;
187 if (m >= 1) {
188 sumMm1 += (A_ge * polyMm1 * m);
189 polyMm1 *= deltaX;
190 }
191 }
192 double oneMXA = 1.0 - XA;
193 double oneMXB = 1.0 - XB;
194 for (size_t k = 0; k < m_kk; k++) {
195 if (iA == k) {
196 lnActCoeff_Scaled_[k] += (oneMXA * XB * sum) + (XA * XB * sumMm1 * (oneMXA + XB));
197 } else if (iB == k) {
198 lnActCoeff_Scaled_[k] += (oneMXB * XA * sum) + (XA * XB * sumMm1 * (-oneMXB - XA));
199 } else {
200 lnActCoeff_Scaled_[k] += -(XA * XB * sum2);
201 }
202 }
203 // Debug against formula in literature
204 }
205}
206
208{
209 double T = temperature();
210 dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
211
212 for (size_t i = 0; i < m_HE_m_ij.size(); i++) {
213 size_t iA = m_pSpecies_A_ij[i];
214 size_t iB = m_pSpecies_B_ij[i];
215 double XA = moleFractions_[iA];
216 double XB = moleFractions_[iB];
217 double deltaX = XA - XB;
218 double poly = 1.0;
219 double sum = 0.0;
220 const vector<double>& he_vec = m_HE_m_ij[i];
221 double sumMm1 = 0.0;
222 double polyMm1 = 1.0;
223 double sum2 = 0.0;
224 for (size_t m = 0; m < he_vec.size(); m++) {
225 double h_e = - he_vec[m] / (GasConstant * T * T);
226 sum += h_e * poly;
227 sum2 += h_e * (m + 1) * poly;
228 poly *= deltaX;
229 if (m >= 1) {
230 sumMm1 += (h_e * polyMm1 * m);
231 polyMm1 *= deltaX;
232 }
233 }
234 double oneMXA = 1.0 - XA;
235 double oneMXB = 1.0 - XB;
236 for (size_t k = 0; k < m_kk; k++) {
237 if (iA == k) {
238 dlnActCoeffdT_Scaled_[k] += (oneMXA * XB * sum) + (XA * XB * sumMm1 * (oneMXA + XB));
239 } else if (iB == k) {
240 dlnActCoeffdT_Scaled_[k] += (oneMXB * XA * sum) + (XA * XB * sumMm1 * (-oneMXB - XA));
241 } else {
242 dlnActCoeffdT_Scaled_[k] += -(XA * XB * sum2);
243 }
244 }
245 }
246
247 for (size_t k = 0; k < m_kk; k++) {
249 }
250}
251
252void RedlichKisterVPSSTP::getdlnActCoeffdT(double* dlnActCoeffdT) const
253{
255 for (size_t k = 0; k < m_kk; k++) {
256 dlnActCoeffdT[k] = dlnActCoeffdT_Scaled_[k];
257 }
258}
259
260void RedlichKisterVPSSTP::getd2lnActCoeffdT2(double* d2lnActCoeffdT2) const
261{
263 for (size_t k = 0; k < m_kk; k++) {
264 d2lnActCoeffdT2[k] = d2lnActCoeffdT2_Scaled_[k];
265 }
266}
267
269{
270 double T = temperature();
271 dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
272
273 for (size_t i = 0; i < m_HE_m_ij.size(); i++) {
274 size_t iA = m_pSpecies_A_ij[i];
275 size_t iB = m_pSpecies_B_ij[i];
276 double XA = moleFractions_[iA];
277 double XB = moleFractions_[iB];
278 double deltaX = XA - XB;
279 double poly = 1.0;
280 double sum = 0.0;
281 const vector<double>& he_vec = m_HE_m_ij[i];
282 const vector<double>& se_vec = m_SE_m_ij[i];
283 double sumMm1 = 0.0;
284 double polyMm1 = 1.0;
285 double polyMm2 = 1.0;
286 double sumMm2 = 0.0;
287 for (size_t m = 0; m < he_vec.size(); m++) {
288 double A_ge = (he_vec[m] - T * se_vec[m]) / (GasConstant * T);;
289 sum += A_ge * poly;
290 poly *= deltaX;
291 if (m >= 1) {
292 sumMm1 += (A_ge * polyMm1 * m);
293 polyMm1 *= deltaX;
294 }
295 if (m >= 2) {
296 sumMm2 += (A_ge * polyMm2 * m * (m - 1.0));
297 polyMm2 *= deltaX;
298 }
299 }
300
301 for (size_t k = 0; k < m_kk; k++) {
302 if (iA == k) {
304 XA * (- (1-XA+XB) * sum + 2*(1.0 - XA) * XB * sumMm1
305 + sumMm1 * (XB * (1 - 2*XA + XB) - XA * (1 - XA + 2*XB))
306 + 2 * XA * XB * sumMm2 * (1.0 - XA + XB));
307 } else if (iB == k) {
309 XB * (- (1-XB+XA) * sum - 2*(1.0 - XB) * XA * sumMm1
310 + sumMm1 * (XA * (2*XB - XA - 1) - XB * (-2*XA + XB - 1))
311 - 2 * XA * XB * sumMm2 * (-XA - 1 + XB));
312 }
313 }
314 }
315}
316
318{
319 double T = temperature();
321
322 for (size_t i = 0; i < m_HE_m_ij.size(); i++) {
323 size_t iA = m_pSpecies_A_ij[i];
324 size_t iB = m_pSpecies_B_ij[i];
325 double XA = moleFractions_[iA];
326 double XB = moleFractions_[iB];
327 double deltaX = XA - XB;
328 double poly = 1.0;
329 double sum = 0.0;
330 const vector<double>& he_vec = m_HE_m_ij[i];
331 const vector<double>& se_vec = m_SE_m_ij[i];
332 double sumMm1 = 0.0;
333 double polyMm1 = 1.0;
334 double polyMm2 = 1.0;
335 double sum2 = 0.0;
336 double sum2Mm1 = 0.0;
337 double sumMm2 = 0.0;
338 for (size_t m = 0; m < he_vec.size(); m++) {
339 double A_ge = he_vec[m] - T * se_vec[m];
340 sum += A_ge * poly;
341 sum2 += A_ge * (m + 1) * poly;
342 poly *= deltaX;
343 if (m >= 1) {
344 sumMm1 += (A_ge * polyMm1 * m);
345 sum2Mm1 += (A_ge * polyMm1 * m * (1.0 + m));
346 polyMm1 *= deltaX;
347 }
348 if (m >= 2) {
349 sumMm2 += (A_ge * polyMm2 * m * (m - 1.0));
350 polyMm2 *= deltaX;
351 }
352 }
353
354 for (size_t k = 0; k < m_kk; k++) {
355 if (iA == k) {
356 dlnActCoeff_dX_(k, iA) += (- XB * sum + (1.0 - XA) * XB * sumMm1
357 + XB * sumMm1 * (1.0 - 2.0 * XA + XB)
358 + XA * XB * sumMm2 * (1.0 - XA + XB));
359
360 dlnActCoeff_dX_(k, iB) += ((1.0 - XA) * sum - (1.0 - XA) * XB * sumMm1
361 + XA * sumMm1 * (1.0 + 2.0 * XB - XA)
362 - XA * XB * sumMm2 * (1.0 - XA + XB));
363 } else if (iB == k) {
364 dlnActCoeff_dX_(k, iA) += ((1.0 - XB) * sum + (1.0 - XA) * XB * sumMm1
365 + XB * sumMm1 * (1.0 - 2.0 * XA + XB)
366 + XA * XB * sumMm2 * (1.0 - XA + XB));
367
368 dlnActCoeff_dX_(k, iB) += (- XA * sum - (1.0 - XB) * XA * sumMm1
369 + XA * sumMm1 * (XB - XA - (1.0 - XB))
370 - XA * XB * sumMm2 * (-XA - (1.0 - XB)));
371 } else {
372 dlnActCoeff_dX_(k, iA) += (- XB * sum2 - XA * XB * sum2Mm1);
373 dlnActCoeff_dX_(k, iB) += (- XA * sum2 + XA * XB * sum2Mm1);
374 }
375 }
376 }
377}
378
379void RedlichKisterVPSSTP::getdlnActCoeffds(const double dTds, const double* const dXds,
380 double* dlnActCoeffds) const
381{
384 for (size_t k = 0; k < m_kk; k++) {
385 dlnActCoeffds[k] = dlnActCoeffdT_Scaled_[k] * dTds;
386 for (size_t j = 0; j < m_kk; j++) {
387 dlnActCoeffds[k] += dlnActCoeff_dX_(k, j) * dXds[j];
388 }
389 }
390}
391
392void RedlichKisterVPSSTP::getdlnActCoeffdlnN_diag(double* dlnActCoeffdlnN_diag) const
393{
395 for (size_t j = 0; j < m_kk; j++) {
396 dlnActCoeffdlnN_diag[j] = dlnActCoeff_dX_(j, j);
397 for (size_t k = 0; k < m_kk; k++) {
398 dlnActCoeffdlnN_diag[k] -= dlnActCoeff_dX_(j, k) * moleFractions_[k];
399 }
400 }
401}
402
403void RedlichKisterVPSSTP::getdlnActCoeffdlnX_diag(double* dlnActCoeffdlnX_diag) const
404{
406 for (size_t k = 0; k < m_kk; k++) {
407 dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
408 }
409}
410
411void RedlichKisterVPSSTP::getdlnActCoeffdlnN(const size_t ld, double* dlnActCoeffdlnN)
412{
414 double* data = & dlnActCoeffdlnN_(0,0);
415 for (size_t k = 0; k < m_kk; k++) {
416 for (size_t m = 0; m < m_kk; m++) {
417 dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m];
418 }
419 }
420}
421
423 const string& speciesA, const string& speciesB,
424 const double* excess_enthalpy, size_t n_enthalpy,
425 const double* excess_entropy, size_t n_entropy)
426{
427 size_t kA = speciesIndex(speciesA);
428 size_t kB = speciesIndex(speciesB);
429 if (kA == npos) {
430 throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
431 "Species '{}' not present in phase", speciesA);
432 } else if (kB == npos) {
433 throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
434 "Species '{}' not present in phase", speciesB);
435 }
436 if (charge(kA) != 0) {
437 throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
438 "Species '{}' should be neutral", speciesA);
439 } else if (charge(kB) != 0) {
440 throw CanteraError("RedlichKisterVPSSTP::addBinaryInteraction",
441 "Species '{}' should be neutral", speciesB);
442 }
443
444 m_pSpecies_A_ij.push_back(kA);
445 m_pSpecies_B_ij.push_back(kB);
446 m_HE_m_ij.emplace_back(excess_enthalpy, excess_enthalpy + n_enthalpy);
447 m_SE_m_ij.emplace_back(excess_entropy, excess_entropy + n_entropy);
448 size_t N = max(n_enthalpy, n_entropy);
449 m_HE_m_ij.back().resize(N, 0.0);
450 m_SE_m_ij.back().resize(N, 0.0);
451 dlnActCoeff_dX_.resize(N, N, 0.0);
452}
453
454}
(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:432
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1477
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:855
double temperature() const
Temperature (K).
Definition Phase.h:563
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:874
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:539
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.
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.
void getdlnActCoeffdT(double *dlnActCoeffdT) const override
Get the array of temperature derivatives of the log activity coefficients.
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 double cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
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:595
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.