Cantera  2.1.2
MolarityIonicVPSSTP.cpp
Go to the documentation of this file.
1 /**
2  * @file MolarityIonicVPSSTP.cpp
3  * Definitions for intermediate ThermoPhase object for phases which
4  * employ excess gibbs free energy formulations
5  * (see \ref thermoprops
6  * and class \link Cantera::MolarityIonicVPSSTP MolarityIonicVPSSTP\endlink).
7  *
8  * Header file for a derived class of ThermoPhase that handles
9  * variable pressure standard state methods for calculating
10  * thermodynamic properties that are further based upon expressions
11  * for the excess gibbs free energy expressed as a function of
12  * the mole fractions.
13  */
14 /*
15  * Copyright (2009) Sandia Corporation. Under the terms of
16  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
17  * U.S. Government retains certain rights in this software.
18  */
22 
23 #include <cstdio>
24 
25 using namespace std;
26 
27 namespace Cantera
28 {
29 
30 MolarityIonicVPSSTP::MolarityIonicVPSSTP() :
32  PBType_(PBTYPE_PASSTHROUGH),
33  numPBSpecies_(m_kk),
34  indexSpecialSpecies_(npos),
35  numCationSpecies_(0),
36  numAnionSpecies_(0),
37  numPassThroughSpecies_(0),
38  neutralPBindexStart(0)
39 {
40 }
41 
42 MolarityIonicVPSSTP::MolarityIonicVPSSTP(const std::string& inputFile,
43  const std::string& id_) :
45  PBType_(PBTYPE_PASSTHROUGH),
46  numPBSpecies_(m_kk),
47  indexSpecialSpecies_(npos),
48  numCationSpecies_(0),
49  numAnionSpecies_(0),
50  numPassThroughSpecies_(0),
51  neutralPBindexStart(0)
52 {
53  initThermoFile(inputFile, id_);
54 }
55 
57  const std::string& id_) :
59  PBType_(PBTYPE_PASSTHROUGH),
60  numPBSpecies_(m_kk),
61  indexSpecialSpecies_(npos),
62  numCationSpecies_(0),
63  numAnionSpecies_(0),
64  numPassThroughSpecies_(0),
65  neutralPBindexStart(0)
66 {
67  importPhase(*findXMLPhase(&phaseRoot, id_), this);
68 }
69 
72  PBType_(PBTYPE_PASSTHROUGH),
73  numPBSpecies_(m_kk),
74  indexSpecialSpecies_(npos),
75  numCationSpecies_(0),
76  numAnionSpecies_(0),
77  numPassThroughSpecies_(0),
78  neutralPBindexStart(0)
79 {
80  *this = operator=(b);
81 }
82 
85 {
86  if (&b != this) {
88  }
89 
90  PBType_ = b.PBType_;
93  PBMoleFractions_ = b.PBMoleFractions_;
96  anionList_ = b.anionList_;
97  numAnionSpecies_ = b.numAnionSpecies_;
98  passThroughList_ = b.passThroughList_;
99  numPassThroughSpecies_ = b.numPassThroughSpecies_;
100  neutralPBindexStart = b.neutralPBindexStart;
101  moleFractionsTmp_ = b.moleFractionsTmp_;
102 
103  return *this;
104 }
105 
108 {
109  return new MolarityIonicVPSSTP(*this);
110 }
111 
112 /*
113  * -------------- Utilities -------------------------------
114  */
115 
117 {
118  return 0;
119 }
120 
121 /*
122  * - Activities, Standard States, Activity Concentrations -----------
123  */
124 
126 {
127  /*
128  * Update the activity coefficients
129  */
131 
132  /*
133  * take the exp of the internally stored coefficients.
134  */
135  for (size_t k = 0; k < m_kk; k++) {
136  lnac[k] = lnActCoeff_Scaled_[k];
137  }
138 }
139 
140 void MolarityIonicVPSSTP::getChemPotentials(doublereal* mu) const
141 {
142  doublereal xx;
143  /*
144  * First get the standard chemical potentials in
145  * molar form.
146  * -> this requires updates of standard state as a function
147  * of T and P
148  */
150  /*
151  * Update the activity coefficients
152  */
154  /*
155  *
156  */
157  doublereal RT = GasConstant * temperature();
158  for (size_t k = 0; k < m_kk; k++) {
159  xx = std::max(moleFractions_[k], SmallNumber);
160  mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]);
161  }
162 }
163 
165 {
166  getChemPotentials(mu);
167  double ve = Faraday * electricPotential();
168  for (size_t k = 0; k < m_kk; k++) {
169  mu[k] += ve*charge(k);
170  }
171 }
172 
174 {
175  /*
176  * Get the nondimensional standard state enthalpies
177  */
178  getEnthalpy_RT(hbar);
179  /*
180  * dimensionalize it.
181  */
182  double T = temperature();
183  double RT = GasConstant * T;
184  for (size_t k = 0; k < m_kk; k++) {
185  hbar[k] *= RT;
186  }
187  /*
188  * Update the activity coefficients, This also update the
189  * internally stored molalities.
190  */
193  double RTT = RT * T;
194  for (size_t k = 0; k < m_kk; k++) {
195  hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k];
196  }
197 }
198 
199 void MolarityIonicVPSSTP::getPartialMolarCp(doublereal* cpbar) const
200 {
201  /*
202  * Get the nondimensional standard state entropies
203  */
204  getCp_R(cpbar);
205  double T = temperature();
206  /*
207  * Update the activity coefficients, This also update the
208  * internally stored molalities.
209  */
212 
213  for (size_t k = 0; k < m_kk; k++) {
214  cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
215  }
216  /*
217  * dimensionalize it.
218  */
219  for (size_t k = 0; k < m_kk; k++) {
220  cpbar[k] *= GasConstant;
221  }
222 }
223 
225 {
226  double xx;
227  /*
228  * Get the nondimensional standard state entropies
229  */
230  getEntropy_R(sbar);
231  double T = temperature();
232  /*
233  * Update the activity coefficients, This also update the
234  * internally stored molalities.
235  */
238 
239  for (size_t k = 0; k < m_kk; k++) {
240  xx = std::max(moleFractions_[k], SmallNumber);
241  sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
242  }
243  /*
244  * dimensionalize it.
245  */
246  for (size_t k = 0; k < m_kk; k++) {
247  sbar[k] *= GasConstant;
248  }
249 }
250 
251 void MolarityIonicVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
252 {
253  /*
254  * Get the standard state values in m^3 kmol-1
255  */
256  getStandardVolumes(vbar);
257  for (size_t iK = 0; iK < m_kk; iK++) {
258  vbar[iK] += 0.0;
259  }
260 }
261 
263 {
264  size_t k;
265  size_t kCat;
266  size_t kMax;
267  doublereal sumCat;
268  doublereal sumAnion;
269  doublereal chP, chM;
270  doublereal sum = 0.0;
271  doublereal sumMax;
272  switch (PBType_) {
273  case PBTYPE_PASSTHROUGH:
274  for (k = 0; k < m_kk; k++) {
275  PBMoleFractions_[k] = moleFractions_[k];
276  }
277  break;
278  case PBTYPE_SINGLEANION:
279  sumCat = 0.0;
280  sumAnion = 0.0;
281  for (k = 0; k < m_kk; k++) {
282  moleFractionsTmp_[k] = moleFractions_[k];
283  }
284  kMax = npos;
285  sumMax = 0.0;
286  for (k = 0; k < cationList_.size(); k++) {
287  kCat = cationList_[k];
288  chP = m_speciesCharge[kCat];
289  if (moleFractions_[kCat] > sumMax) {
290  kMax = k;
291  sumMax = moleFractions_[kCat];
292  }
293  sumCat += chP * moleFractions_[kCat];
294  }
295  k = anionList_[0];
296  chM = m_speciesCharge[k];
297  sumAnion = moleFractions_[k] * chM;
298  sum = sumCat - sumAnion;
299  if (fabs(sum) > 1.0E-16) {
300  moleFractionsTmp_[cationList_[kMax]] -= sum / m_speciesCharge[kMax];
301  sum = 0.0;
302  for (k = 0; k < numCationSpecies_; k++) {
303  sum += moleFractionsTmp_[k];
304  }
305  for (k = 0; k < numCationSpecies_; k++) {
306  moleFractionsTmp_[k]/= sum;
307  }
308  }
309 
310  for (k = 0; k < numCationSpecies_; k++) {
311  PBMoleFractions_[k] = moleFractionsTmp_[cationList_[k]];
312  }
313  for (k = 0; k < numPassThroughSpecies_; k++) {
314  PBMoleFractions_[neutralPBindexStart + k] = moleFractions_[passThroughList_[k]];
315  }
316 
317  sum = std::max(0.0, PBMoleFractions_[0]);
318  for (k = 1; k < numPBSpecies_; k++) {
319  sum += PBMoleFractions_[k];
320 
321  }
322  for (k = 0; k < numPBSpecies_; k++) {
323  PBMoleFractions_[k] /= sum;
324  }
325 
326  break;
327  case PBTYPE_SINGLECATION:
328  throw CanteraError("eosType", "Unknown type");
329 
330  break;
331 
332  case PBTYPE_MULTICATIONANION:
333  throw CanteraError("eosType", "Unknown type");
334 
335  break;
336  default:
337  throw CanteraError("eosType", "Unknown type");
338  break;
339 
340  }
341 }
342 
344 {
345  for (size_t k = 0; k < m_kk; k++) {
346  lnActCoeff_Scaled_[k] = 0.0;
347  }
348 }
349 
351 {
352 
353 
354 }
355 
357 {
358 
359 }
360 
361 /*
362  * ------------ Partial Molar Properties of the Solution ------------
363  */
364 
365 doublereal MolarityIonicVPSSTP::err(const std::string& msg) const
366 {
367  throw CanteraError("MolarityIonicVPSSTP","Base class method "
368  +msg+" called. Equation of state type: "+int2str(eosType()));
369  return 0;
370 }
371 
373 {
375  initLengths();
376  /*
377  * Go find the list of cations and anions
378  */
379  double ch;
380  numCationSpecies_ = 0;
381  cationList_.clear();
382  anionList_.clear();
383  passThroughList_.clear();
384  for (size_t k = 0; k < m_kk; k++) {
385  ch = m_speciesCharge[k];
386  if (ch > 0.0) {
387  cationList_.push_back(k);
389  } else if (ch < 0.0) {
390  anionList_.push_back(k);
391  numAnionSpecies_++;
392  } else {
393  passThroughList_.push_back(k);
394  numPassThroughSpecies_++;
395  }
396  }
397  numPBSpecies_ = numCationSpecies_ + numAnionSpecies_ - 1;
398  neutralPBindexStart = numPBSpecies_;
399  PBType_ = PBTYPE_MULTICATIONANION;
400  if (numAnionSpecies_ == 1) {
401  PBType_ = PBTYPE_SINGLEANION;
402  } else if (numCationSpecies_ == 1) {
403  PBType_ = PBTYPE_SINGLECATION;
404  }
405  if (numAnionSpecies_ == 0 && numCationSpecies_ == 0) {
406  PBType_ = PBTYPE_PASSTHROUGH;
407  }
408 }
409 
411 {
412  m_kk = nSpecies();
413  moleFractionsTmp_.resize(m_kk);
414 }
415 
416 void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
417 {
418  std::string subname = "MolarityIonicVPSSTP::initThermoXML";
419  std::string stemp;
420 
421  if ((int) id.size() > 0) {
422  string idp = phaseNode.id();
423  if (idp != id) {
424  throw CanteraError(subname, "phasenode and Id are incompatible");
425  }
426  }
427 
428  /*
429  * Check on the thermo field. Must have one of:
430  * <thermo model="MolarityIonicVPSS" />
431  * <thermo model="MolarityIonicVPSSTP" />
432  */
433  if (!phaseNode.hasChild("thermo")) {
434  throw CanteraError(subname, "no thermo XML node");
435  }
436  XML_Node& thermoNode = phaseNode.child("thermo");
437  std::string mStringa = thermoNode.attrib("model");
438  std::string mString = lowercase(mStringa);
439  if (mString != "molarityionicvpss" && mString != "molarityionicvpsstp") {
440  throw CanteraError(subname.c_str(),
441  "Unknown thermo model: " + mStringa + " - This object only knows \"MolarityIonicVPSSTP\" ");
442  }
443 
444  /*
445  * Go get all of the coefficients and factors in the
446  * activityCoefficients XML block
447  */
448  XML_Node* acNodePtr = 0;
449  if (thermoNode.hasChild("activityCoefficients")) {
450  XML_Node& acNode = thermoNode.child("activityCoefficients");
451  acNodePtr = &acNode;
452  mStringa = acNode.attrib("model");
453  mString = lowercase(mStringa);
454  // if (mString != "redlich-kister") {
455  // throw CanteraError(subname.c_str(),
456  // "Unknown activity coefficient model: " + mStringa);
457  //}
458  size_t n = acNodePtr->nChildren();
459  for (size_t i = 0; i < n; i++) {
460  XML_Node& xmlACChild = acNodePtr->child(i);
461  stemp = xmlACChild.name();
462  std::string nodeName = lowercase(stemp);
463  /*
464  * Process a binary interaction
465  */
466  if (nodeName == "binaryneutralspeciesparameters") {
467  readXMLBinarySpecies(xmlACChild);
468  }
469  }
470  }
471 
472 
473  /*
474  * Go down the chain
475  */
476  GibbsExcessVPSSTP::initThermoXML(phaseNode, id);
477 }
478 
480 {
481  std::string xname = xmLBinarySpecies.name();
482 
483 }
484 
485 std::string MolarityIonicVPSSTP::report(bool show_thermo) const
486 {
487  char p[800];
488  string s = "";
489  try {
490  if (name() != "") {
491  sprintf(p, " \n %s:\n", name().c_str());
492  s += p;
493  }
494  sprintf(p, " \n temperature %12.6g K\n", temperature());
495  s += p;
496  sprintf(p, " pressure %12.6g Pa\n", pressure());
497  s += p;
498  sprintf(p, " density %12.6g kg/m^3\n", density());
499  s += p;
500  sprintf(p, " mean mol. weight %12.6g amu\n", meanMolecularWeight());
501  s += p;
502 
503  doublereal phi = electricPotential();
504  sprintf(p, " potential %12.6g V\n", phi);
505  s += p;
506 
507  size_t kk = nSpecies();
508  vector_fp x(kk);
509  vector_fp molal(kk);
510  vector_fp mu(kk);
511  vector_fp muss(kk);
512  vector_fp acMolal(kk);
513  vector_fp actMolal(kk);
514  getMoleFractions(&x[0]);
515 
516  getChemPotentials(&mu[0]);
517  getStandardChemPotentials(&muss[0]);
518  getActivities(&actMolal[0]);
519 
520 
521  if (show_thermo) {
522  sprintf(p, " \n");
523  s += p;
524  sprintf(p, " 1 kg 1 kmol\n");
525  s += p;
526  sprintf(p, " ----------- ------------\n");
527  s += p;
528  sprintf(p, " enthalpy %12.6g %12.4g J\n",
530  s += p;
531  sprintf(p, " internal energy %12.6g %12.4g J\n",
533  s += p;
534  sprintf(p, " entropy %12.6g %12.4g J/K\n",
536  s += p;
537  sprintf(p, " Gibbs function %12.6g %12.4g J\n",
538  gibbs_mass(), gibbs_mole());
539  s += p;
540  sprintf(p, " heat capacity c_p %12.6g %12.4g J/K\n",
541  cp_mass(), cp_mole());
542  s += p;
543  try {
544  sprintf(p, " heat capacity c_v %12.6g %12.4g J/K\n",
545  cv_mass(), cv_mole());
546  s += p;
547  } catch (CanteraError& e) {
548  e.save();
549  sprintf(p, " heat capacity c_v <not implemented> \n");
550  s += p;
551  }
552  }
553 
554  } catch (CanteraError& e) {
555  e.save();
556  }
557  return s;
558 }
559 
560 }
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:534
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
doublereal electricPotential() const
Returns the electric potential of this phase (V).
Definition: ThermoPhase.h:391
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
Definition: ThermoPhase.h:279
XML_Node * findXMLPhase(XML_Node *root, const std::string &idtarget)
Search an XML_Node tree for a named phase XML_Node.
Definition: xml.cpp:1104
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:534
GibbsExcessVPSSTP & operator=(const GibbsExcessVPSSTP &b)
Assignment operator.
size_t numCationSpecies_
Number of cations in the mixture.
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of each species in their standard states at the current T and P of the solution...
virtual void calcPseudoBinaryMoleFractions() const
Calculate pseudo binary mole fractions.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
virtual void getActivities(doublereal *ac) const
Get the array of non-dimensional activities (molality based for this class and classes that derive fr...
vector_fp m_speciesCharge
Vector of species charges. length m_kk.
Definition: Phase.h:731
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
virtual void getPartialMolarCp(doublereal *cpbar) const
Returns an array of partial molar entropies for the species in the mixture.
virtual void getLnActivityCoefficients(doublereal *lnac) const
Get the array of non-dimensional molar-based ln activity coefficients at the current solution tempera...
void readXMLBinarySpecies(XML_Node &xmlBinarySpecies)
Process an XML node called "binaryNeutralSpeciesParameters".
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
doublereal size(size_t k) const
This routine returns the size of species k.
Definition: Phase.h:398
std::vector< size_t > cationList_
Vector of cation indices in the mixture.
std::string lowercase(const std::string &s)
Cast a copy of a string to lower case.
Definition: stringUtils.cpp:58
std::vector< doublereal > d2lnActCoeffdT2_Scaled_
Storage for the current derivative values of the gradients with respect to temperature of the log of ...
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:519
std::string name() const
Return the name of the phase.
Definition: Phase.cpp:129
Header for intermediate ThermoPhase object for phases which employ gibbs excess free energy based for...
MolarityIonicVPSSTP & operator=(const MolarityIonicVPSSTP &b)
Assignment operator.
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
Definition: ThermoPhase.h:274
XML_Node & child(const size_t n) const
Return a changeable reference to the n'th child of the current node.
Definition: xml.cpp:584
doublereal intEnergy_mass() const
Specific internal energy.
Definition: ThermoPhase.h:940
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:101
virtual std::string report(bool show_thermo=true) const
returns a summary of the state of the phase as a string
bool importPhase(XML_Node &phase, ThermoPhase *th, SpeciesThermoFactory *spfactory)
Import a phase information into an empty thermophase object.
virtual doublereal intEnergy_mole() const
Molar internal energy. Units: J/kmol.
Definition: ThermoPhase.h:269
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
Definition: ThermoPhase.h:264
void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object.
doublereal gibbs_mass() const
Specific Gibbs function.
Definition: ThermoPhase.h:954
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
void s_update_dlnActCoeff_dT() const
Update the derivative of the log of the activity coefficients wrt T.
doublereal entropy_mass() const
Specific entropy.
Definition: ThermoPhase.h:947
doublereal pressure() const
Returns the current pressure of the phase.
std::string name() const
Returns the name of the XML node.
Definition: xml.h:390
doublereal cp_mass() const
Specific heat at constant pressure.
Definition: ThermoPhase.h:961
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
doublereal err(const std::string &msg) const
Error function.
doublereal cv_mass() const
Specific heat at constant volume.
Definition: ThermoPhase.h:968
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:574
void getElectrochemPotentials(doublereal *mu) const
Get the species electrochemical potentials.
void initLengths()
Initialize lengths of local variables after all species have been identified.
virtual int eosType() const
Equation of state type flag.
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
Definition: ThermoPhase.h:289
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:252
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Enthalpy functions for the standard state species at the current T an...
doublereal temperature() const
Temperature (K).
Definition: Phase.h:528
std::string id() const
Return the id attribute, if present.
Definition: xml.cpp:467
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:139
std::vector< doublereal > lnActCoeff_Scaled_
Storage for the current values of the activity coefficients of the species.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:165
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies for the species in the mixture.
doublereal enthalpy_mass() const
Specific enthalpy.
Definition: ThermoPhase.h:933
std::vector< doublereal > moleFractions_
Storage for the current values of the mole fractions of the species.
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:588
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:66
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Initialize a ThermoPhase object, potentially reading activity coefficient information from an XML dat...
Contains declarations for string manipulation functions within Cantera.
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
size_t numPBSpecies_
Number of pseudo binary species.
virtual ThermoPhase * duplMyselfAsThermoPhase() const
Duplication routine for objects which inherit from ThermoPhase.
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
Definition: ThermoPhase.h:284
size_t m_kk
Number of species in the phase.
Definition: Phase.h:716
virtual void initThermoFile(const std::string &inputFile, const std::string &id)
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity.
size_t nChildren(bool discardComments=false) const
Return the number of children.
Definition: xml.cpp:594
void s_update_dlnActCoeff_dX_() const
Internal routine that calculates the derivative of the activity coefficients wrt the mole fractions...
size_t indexSpecialSpecies_
index of special species
std::vector< doublereal > dlnActCoeffdT_Scaled_
Storage for the current derivative values of the gradients with respect to temperature of the log of ...
void s_update_lnActCoeff() const
Update the activity coefficients.
doublereal 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:504
void save()
Function to put this error onto Cantera's error stack.