Loading [MathJax]/extensions/tex2jax.js
Cantera  3.2.0a1
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
MargulesVPSSTP.cpp
Go to the documentation of this file.
1/**
2 * @file MargulesVPSSTP.cpp
3 * Definitions for ThermoPhase object for phases which
4 * employ excess Gibbs free energy formulations related to Margules
5 * expansions (see @ref thermoprops
6 * and class @link Cantera::MargulesVPSSTP MargulesVPSSTP@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
16namespace Cantera
17{
18MargulesVPSSTP::MargulesVPSSTP(const string& inputFile, const string& id_)
19{
20 initThermoFile(inputFile, id_);
21}
22
23// -- Activities, Standard States, Activity Concentrations -----------
24
26{
27 // Update the activity coefficients
29
30 // take the exp of the internally stored coefficients.
31 for (size_t k = 0; k < m_kk; k++) {
32 lnac[k] = lnActCoeff_Scaled_[k];
33 }
34}
35
36// ------------ Partial Molar Properties of the Solution ------------
37
39{
40 // First get the standard chemical potentials in molar form. This requires
41 // updates of standard state as a function of T and P
43
44 // Update the activity coefficients
46 for (size_t k = 0; k < m_kk; k++) {
47 double xx = std::max(moleFractions_[k], SmallNumber);
48 mu[k] += RT() * (log(xx) + lnActCoeff_Scaled_[k]);
49 }
50}
51
53{
54 return cp_mole() - GasConstant;
55}
56
58{
59 // Get the nondimensional standard state enthalpies
60 getEnthalpy_RT(hbar);
61
62 // dimensionalize it.
63 for (size_t k = 0; k < m_kk; k++) {
64 hbar[k] *= RT();
65 }
66
67 // Update the activity coefficients, This also update the internally stored
68 // molalities.
71 for (size_t k = 0; k < m_kk; k++) {
72 hbar[k] -= RT() * temperature() * dlnActCoeffdT_Scaled_[k];
73 }
74}
75
76void MargulesVPSSTP::getPartialMolarCp(double* cpbar) const
77{
78 // Get the nondimensional standard state entropies
79 getCp_R(cpbar);
80 double T = temperature();
81
82 // Update the activity coefficients, This also update the internally stored
83 // molalities.
86
87 for (size_t k = 0; k < m_kk; k++) {
88 cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
89 }
90 // dimensionalize it.
91 for (size_t k = 0; k < m_kk; k++) {
92 cpbar[k] *= GasConstant;
93 }
94}
95
97{
98 // Get the nondimensional standard state entropies
99 getEntropy_R(sbar);
100 double T = temperature();
101
102 // Update the activity coefficients, This also update the internally stored
103 // molalities.
106
107 for (size_t k = 0; k < m_kk; k++) {
108 double xx = std::max(moleFractions_[k], SmallNumber);
109 sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
110 }
111
112 // dimensionalize it.
113 for (size_t k = 0; k < m_kk; k++) {
114 sbar[k] *= GasConstant;
115 }
116}
117
119{
120 double T = temperature();
121
122 // Get the standard state values in m^3 kmol-1
123 getStandardVolumes(vbar);
124
125 for (size_t i = 0; i < numBinaryInteractions_; i++) {
126 size_t iA = m_pSpecies_A_ij[i];
127 size_t iB = m_pSpecies_B_ij[i];
128 double XA = moleFractions_[iA];
129 double XB = moleFractions_[iB];
130 double g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
131 double g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
132 const double temp1 = g0 + g1 * XB;
133 const double all = -1.0*XA*XB*temp1 - XA*XB*XB*g1;
134
135 for (size_t iK = 0; iK < m_kk; iK++) {
136 vbar[iK] += all;
137 }
138 vbar[iA] += XB * temp1;
139 vbar[iB] += XA * temp1 + XA*XB*g1;
140 }
141}
142
144{
145 initLengths();
146 if (m_input.hasKey("interactions")) {
147 for (auto& item : m_input["interactions"].asVector<AnyMap>()) {
148 auto& species = item["species"].asVector<string>(2);
149 vector<double> h(2), s(2), vh(2), vs(2);
150 if (item.hasKey("excess-enthalpy")) {
151 h = item.convertVector("excess-enthalpy", "J/kmol", 2);
152 }
153 if (item.hasKey("excess-entropy")) {
154 s = item.convertVector("excess-entropy", "J/kmol/K", 2);
155 }
156 if (item.hasKey("excess-volume-enthalpy")) {
157 vh = item.convertVector("excess-volume-enthalpy", "m^3/kmol", 2);
158 }
159 if (item.hasKey("excess-volume-entropy")) {
160 vs = item.convertVector("excess-volume-entropy", "m^3/kmol/K", 2);
161 }
163 h[0], h[1], s[0], s[1], vh[0], vh[1], vs[0], vs[1]);
164 }
165 }
167}
168
170{
172 vector<AnyMap> interactions;
173 for (size_t n = 0; n < m_pSpecies_A_ij.size(); n++) {
174 AnyMap interaction;
175 interaction["species"] = vector<string>{
177 if (m_HE_b_ij[n] != 0 || m_HE_c_ij[n] != 0) {
178 interaction["excess-enthalpy"].setQuantity(
179 {m_HE_b_ij[n], m_HE_c_ij[n]}, "J/kmol");
180 }
181 if (m_SE_b_ij[n] != 0 || m_SE_c_ij[n] != 0) {
182 interaction["excess-entropy"].setQuantity(
183 {m_SE_b_ij[n], m_SE_c_ij[n]}, "J/kmol/K");
184 }
185 if (m_VHE_b_ij[n] != 0 || m_VHE_c_ij[n] != 0) {
186 interaction["excess-volume-enthalpy"].setQuantity(
187 {m_VHE_b_ij[n], m_VHE_c_ij[n]}, "m^3/kmol");
188 }
189 if (m_VSE_b_ij[n] != 0 || m_VSE_c_ij[n] != 0) {
190 interaction["excess-volume-entropy"].setQuantity(
191 {m_VSE_b_ij[n], m_VSE_c_ij[n]}, "m^3/kmol/K");
192 }
193 interactions.push_back(std::move(interaction));
194 }
195 phaseNode["interactions"] = std::move(interactions);
196}
197
199{
201}
202
203void MargulesVPSSTP::addBinaryInteraction(const string& speciesA,
204 const string& speciesB, double h0, double h1, double s0, double s1,
205 double vh0, double vh1, double vs0, double vs1)
206{
207 size_t kA = speciesIndex(speciesA);
208 size_t kB = speciesIndex(speciesB);
209 // The interaction is silently ignored if either species is not defined in
210 // the current phase.
211 if (kA == npos || kB == npos) {
212 return;
213 }
214 m_pSpecies_A_ij.push_back(kA);
215 m_pSpecies_B_ij.push_back(kB);
216
217 m_HE_b_ij.push_back(h0);
218 m_HE_c_ij.push_back(h1);
219 m_SE_b_ij.push_back(s0);
220 m_SE_c_ij.push_back(s1);
221 m_VHE_b_ij.push_back(vh0);
222 m_VHE_c_ij.push_back(vh1);
223 m_VSE_b_ij.push_back(vs0);
224 m_VSE_c_ij.push_back(vs1);
226}
227
228
230{
231 double T = temperature();
232 lnActCoeff_Scaled_.assign(m_kk, 0.0);
233 for (size_t i = 0; i < numBinaryInteractions_; i++) {
234 size_t iA = m_pSpecies_A_ij[i];
235 size_t iB = m_pSpecies_B_ij[i];
236 double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT();
237 double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT();
238 double XA = moleFractions_[iA];
239 double XB = moleFractions_[iB];
240 const double XAXB = XA * XB;
241 const double g0g1XB = (g0 + g1 * XB);
242 const double all = -1.0 * XAXB * g0g1XB - XAXB * XB * g1;
243 for (size_t iK = 0; iK < m_kk; iK++) {
244 lnActCoeff_Scaled_[iK] += all;
245 }
246 lnActCoeff_Scaled_[iA] += XB * g0g1XB;
247 lnActCoeff_Scaled_[iB] += XA * g0g1XB + XAXB * g1;
248 }
249}
250
252{
253 double invT = 1.0 / temperature();
254 double invRTT = 1.0 / GasConstant*invT*invT;
255 dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
256 d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0);
257 for (size_t i = 0; i < numBinaryInteractions_; i++) {
258 size_t iA = m_pSpecies_A_ij[i];
259 size_t iB = m_pSpecies_B_ij[i];
260 double XA = moleFractions_[iA];
261 double XB = moleFractions_[iB];
262 double g0 = -m_HE_b_ij[i] * invRTT;
263 double g1 = -m_HE_c_ij[i] * invRTT;
264 const double XAXB = XA * XB;
265 const double g0g1XB = (g0 + g1 * XB);
266 const double all = -1.0 * XAXB * g0g1XB - XAXB * XB * g1;
267 const double mult = 2.0 * invT;
268 const double dT2all = mult * all;
269 for (size_t iK = 0; iK < m_kk; iK++) {
270 dlnActCoeffdT_Scaled_[iK] += all;
271 d2lnActCoeffdT2_Scaled_[iK] -= dT2all;
272 }
273 dlnActCoeffdT_Scaled_[iA] += XB * g0g1XB;
274 dlnActCoeffdT_Scaled_[iB] += XA * g0g1XB + XAXB * g1;
275 d2lnActCoeffdT2_Scaled_[iA] -= mult * XB * g0g1XB;
276 d2lnActCoeffdT2_Scaled_[iB] -= mult * (XA * g0g1XB + XAXB * g1);
277 }
278}
279
280void MargulesVPSSTP::getdlnActCoeffdT(double* dlnActCoeffdT) const
281{
283 for (size_t k = 0; k < m_kk; k++) {
284 dlnActCoeffdT[k] = dlnActCoeffdT_Scaled_[k];
285 }
286}
287
288void MargulesVPSSTP::getd2lnActCoeffdT2(double* d2lnActCoeffdT2) const
289{
291 for (size_t k = 0; k < m_kk; k++) {
292 d2lnActCoeffdT2[k] = d2lnActCoeffdT2_Scaled_[k];
293 }
294}
295
296void MargulesVPSSTP::getdlnActCoeffds(const double dTds, const double* const dXds,
297 double* dlnActCoeffds) const
298{
299 double T = temperature();
301 for (size_t iK = 0; iK < m_kk; iK++) {
302 dlnActCoeffds[iK] = 0.0;
303 }
304
305 for (size_t i = 0; i < numBinaryInteractions_; i++) {
306 size_t iA = m_pSpecies_A_ij[i];
307 size_t iB = m_pSpecies_B_ij[i];
308 double XA = moleFractions_[iA];
309 double XB = moleFractions_[iB];
310 double dXA = dXds[iA];
311 double dXB = dXds[iB];
312 double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT();
313 double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT();
314 const double g02g1XB = g0 + 2*g1*XB;
315 const double g2XAdXB = 2*g1*XA*dXB;
316 const double all = (-XB * dXA - XA *dXB) * g02g1XB - XB *g2XAdXB;
317 for (size_t iK = 0; iK < m_kk; iK++) {
318 dlnActCoeffds[iK] += all + dlnActCoeffdT_Scaled_[iK]*dTds;
319 }
320 dlnActCoeffds[iA] += dXB * g02g1XB;
321 dlnActCoeffds[iB] += dXA * g02g1XB + g2XAdXB;
322 }
323}
324
326{
327 double T = temperature();
328 dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
329
330 for (size_t iK = 0; iK < m_kk; iK++) {
331 double XK = moleFractions_[iK];
332
333 for (size_t i = 0; i < numBinaryInteractions_; i++) {
334 size_t iA = m_pSpecies_A_ij[i];
335 size_t iB = m_pSpecies_B_ij[i];
336 size_t delAK = 0;
337 size_t delBK = 0;
338
339 if (iA==iK) {
340 delAK = 1;
341 } else if (iB==iK) {
342 delBK = 1;
343 }
344
345 double XA = moleFractions_[iA];
346 double XB = moleFractions_[iB];
347
348 double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT();
349 double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT();
350
351 dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB)));
352 }
354 }
355}
356
358{
359 double T = temperature();
361
362 // Loop over the activity coefficient gamma_k
363 for (size_t iK = 0; iK < m_kk; iK++) {
364 for (size_t iM = 0; iM < m_kk; iM++) {
365 double XM = moleFractions_[iM];
366 for (size_t i = 0; i < numBinaryInteractions_; i++) {
367 size_t iA = m_pSpecies_A_ij[i];
368 size_t iB = m_pSpecies_B_ij[i];
369 double delAK = 0.0;
370 double delBK = 0.0;
371 double delAM = 0.0;
372 double delBM = 0.0;
373 if (iA==iK) {
374 delAK = 1.0;
375 } else if (iB==iK) {
376 delBK = 1.0;
377 }
378 if (iA==iM) {
379 delAM = 1.0;
380 } else if (iB==iM) {
381 delBM = 1.0;
382 }
383
384 double XA = moleFractions_[iA];
385 double XB = moleFractions_[iB];
386 double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT();
387 double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT();
388 dlnActCoeffdlnN_(iK,iM) += g0*((delAM-XA)*(delBK-XB)+(delAK-XA)*(delBM-XB));
389 dlnActCoeffdlnN_(iK,iM) += 2*g1*((delAM-XA)*(delBK-XB)*XB+(delAK-XA)*(delBM-XB)*XB+(delBM-XB)*(delBK-XB)*XA);
390 }
391 dlnActCoeffdlnN_(iK,iM) = XM*dlnActCoeffdlnN_(iK,iM);
392 }
393 }
394}
395
397{
398 double T = temperature();
399 dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
400
401 for (size_t i = 0; i < numBinaryInteractions_; i++) {
402 size_t iA = m_pSpecies_A_ij[i];
403 size_t iB = m_pSpecies_B_ij[i];
404
405 double XA = moleFractions_[iA];
406 double XB = moleFractions_[iB];
407
408 double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT();
409 double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT();
410
411 dlnActCoeffdlnX_diag_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB);
412 dlnActCoeffdlnX_diag_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB);
413 }
414}
415
416void MargulesVPSSTP::getdlnActCoeffdlnN_diag(double* dlnActCoeffdlnN_diag) const
417{
419 for (size_t k = 0; k < m_kk; k++) {
420 dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k];
421 }
422}
423
424void MargulesVPSSTP::getdlnActCoeffdlnX_diag(double* dlnActCoeffdlnX_diag) const
425{
427 for (size_t k = 0; k < m_kk; k++) {
428 dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
429 }
430}
431
432void MargulesVPSSTP::getdlnActCoeffdlnN(const size_t ld, double* dlnActCoeffdlnN)
433{
435 double* data = & dlnActCoeffdlnN_(0,0);
436 for (size_t k = 0; k < m_kk; k++) {
437 for (size_t m = 0; m < m_kk; m++) {
438 dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m];
439 }
440 }
441}
442
443}
(see Thermodynamic Properties and class MargulesVPSSTP).
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
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 ...
vector< double > dlnActCoeffdlnN_diag_
Storage for the current derivative values of the gradients with respect to logarithm of the mole frac...
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 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.
vector< double > m_VHE_c_ij
Enthalpy term for the ternary mole fraction interaction of the excess Gibbs free energy expression.
vector< double > m_SE_b_ij
Entropy term for the binary mole fraction interaction of the excess Gibbs free energy expression.
size_t numBinaryInteractions_
number of binary interaction expressions
vector< double > m_SE_c_ij
Entropy term for the ternary mole fraction interaction of the excess Gibbs free energy expression.
void s_update_dlnActCoeff_dlnN_diag() const
Update the derivative of the log of the activity coefficients wrt log(moles) - diagonal only.
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
vector< double > m_VSE_c_ij
Entropy term for the ternary mole fraction interaction of the excess Gibbs free energy expression.
void s_update_dlnActCoeff_dlnN() const
Update the derivative of the log of the activity coefficients wrt log(moles_m)
vector< double > m_HE_b_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.
vector< double > m_VSE_b_ij
Entropy term for the binary mole fraction interaction of the excess Gibbs free energy expression.
void getPartialMolarCp(double *cpbar) const override
Returns an array of partial molar entropies 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
Update the derivative of the log of the activity coefficients wrt log(mole fraction)
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...
vector< double > m_VHE_b_ij
Enthalpy term for the binary mole fraction interaction of the excess Gibbs free energy expression.
void addBinaryInteraction(const string &speciesA, const string &speciesB, double h0, double h1, double s0, double s1, double vh0, double vh1, double vs0, double vs1)
Add a binary species interaction with the specified parameters.
MargulesVPSSTP(const string &inputFile="", const string &id="")
Construct a MargulesVPSSTP object from an input file.
vector< double > m_HE_c_ij
Enthalpy term for the ternary mole fraction interaction of the excess Gibbs free energy expression.
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 ...
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
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.