Cantera  2.0
PseudoBinaryVPSSTP.cpp
Go to the documentation of this file.
1 /**
2  * @file PseudoBinaryVPSSTP.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::PseudoBinaryVPSSTP PseudoBinaryVPSSTP\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  */
21 
22 #include <cmath>
23 #include <cstdio>
24 
25 using namespace std;
26 
27 namespace Cantera
28 {
29 
30 /*
31  * Default constructor.
32  *
33  */
34 PseudoBinaryVPSSTP::PseudoBinaryVPSSTP() :
36  PBType_(PBTYPE_PASSTHROUGH),
37  numPBSpecies_(m_kk),
38  indexSpecialSpecies_(npos),
39  numCationSpecies_(0),
40  numAnionSpecies_(0),
41  numPassThroughSpecies_(0),
42  neutralPBindexStart(0),
43  cationPhase_(0),
44  anionPhase_(0)
45 {
46 }
47 
48 /*
49  * Copy Constructor:
50  *
51  * Note this stuff will not work until the underlying phase
52  * has a working copy constructor
53  */
56  PBType_(PBTYPE_PASSTHROUGH),
57  numPBSpecies_(m_kk),
58  indexSpecialSpecies_(npos),
59  numCationSpecies_(0),
60  numAnionSpecies_(0),
61  numPassThroughSpecies_(0),
62  neutralPBindexStart(0),
63  cationPhase_(0),
64  anionPhase_(0)
65 {
66  *this = operator=(b);
67 }
68 
69 /*
70  * operator=()
71  *
72  * Note this stuff will not work until the underlying phase
73  * has a working assignment operator
74  */
77 {
78  if (&b != this) {
80  }
81 
82  PBType_ = b.PBType_;
85  PBMoleFractions_ = b.PBMoleFractions_;
86  cationList_ = b.cationList_;
87  numCationSpecies_ = b.numCationSpecies_;
88  anionList_ = b.anionList_;
89  numAnionSpecies_ = b.numAnionSpecies_;
90  passThroughList_ = b.passThroughList_;
91  numPassThroughSpecies_ = b.numPassThroughSpecies_;
92  neutralPBindexStart = b.neutralPBindexStart;
93  cationPhase_ = b.cationPhase_;
94  anionPhase_ = b.anionPhase_;
95  moleFractionsTmp_ = b.moleFractionsTmp_;
96 
97  return *this;
98 }
99 
100 /**
101  *
102  * ~PseudoBinaryVPSSTP(): (virtual)
103  *
104  * Destructor: does nothing:
105  *
106  */
108 {
109 }
110 
111 /*
112  * This routine duplicates the current object and returns
113  * a pointer to ThermoPhase.
114  */
117 {
118  PseudoBinaryVPSSTP* mtp = new PseudoBinaryVPSSTP(*this);
119  return (ThermoPhase*) mtp;
120 }
121 
122 /*
123  * -------------- Utilities -------------------------------
124  */
125 
126 
127 // Equation of state type flag.
128 /*
129  * The ThermoPhase base class returns
130  * zero. Subclasses should define this to return a unique
131  * non-zero value. Known constants defined for this purpose are
132  * listed in mix_defs.h. The PseudoBinaryVPSSTP class also returns
133  * zero, as it is a non-complete class.
134  */
136 {
137  return 0;
138 }
139 
140 
141 
142 /*
143  * ------------ Molar Thermodynamic Properties ----------------------
144  */
145 
146 
147 /*
148  * - Activities, Standard States, Activity Concentrations -----------
149  */
150 
151 
153 {
154  err("standardConcentration");
155  return -1.0;
156 }
157 
158 doublereal PseudoBinaryVPSSTP::logStandardConc(size_t k) const
159 {
160  err("logStandardConc");
161  return -1.0;
162 }
163 
164 
165 
167 {
168  getChemPotentials(mu);
169  double ve = Faraday * electricPotential();
170  for (size_t k = 0; k < m_kk; k++) {
171  mu[k] += ve*charge(k);
172  }
173 }
174 
176 {
177  size_t k;
178  doublereal sumCat;
179  doublereal sumAnion;
180  doublereal sum = 0.0;
181  switch (PBType_) {
182  case PBTYPE_PASSTHROUGH:
183  for (k = 0; k < m_kk; k++) {
184  PBMoleFractions_[k] = moleFractions_[k];
185  }
186  break;
187  case PBTYPE_SINGLEANION:
188  sumCat = 0.0;
189  sumAnion = 0.0;
190  for (k = 0; k < m_kk; k++) {
191  moleFractionsTmp_[k] = moleFractions_[k];
192  }
193  for (k = 0; k < cationList_.size(); k++) {
194  sumCat += moleFractions_[cationList_[k]];
195  }
196  sumAnion = moleFractions_[anionList_[k]];
197  PBMoleFractions_[0] = sumCat -sumAnion;
198  moleFractionsTmp_[indexSpecialSpecies_] -= PBMoleFractions_[0];
199 
200 
201  for (k = 0; k < numCationSpecies_; k++) {
202  PBMoleFractions_[1+k] = moleFractionsTmp_[cationList_[k]];
203  }
204 
205  for (k = 0; k < numPassThroughSpecies_; k++) {
206  PBMoleFractions_[neutralPBindexStart + k] =
207  moleFractions_[cationList_[k]];
208  }
209 
210  sum = std::max(0.0, PBMoleFractions_[0]);
211  for (k = 1; k < numPBSpecies_; k++) {
212  sum += PBMoleFractions_[k];
213  }
214  for (k = 0; k < numPBSpecies_; k++) {
215  PBMoleFractions_[k] /= sum;
216  }
217 
218  break;
219  case PBTYPE_SINGLECATION:
220  throw CanteraError("eosType", "Unknown type");
221 
222  break;
223 
224  case PBTYPE_MULTICATIONANION:
225  throw CanteraError("eosType", "Unknown type");
226 
227  break;
228  default:
229  throw CanteraError("eosType", "Unknown type");
230  break;
231 
232  }
233 }
234 
235 /*
236  * ------------ Partial Molar Properties of the Solution ------------
237  */
238 
239 
240 doublereal PseudoBinaryVPSSTP::err(std::string msg) const
241 {
242  throw CanteraError("PseudoBinaryVPSSTP","Base class method "
243  +msg+" called. Equation of state type: "+int2str(eosType()));
244  return 0;
245 }
246 
247 
248 /*
249  * @internal Initialize. This method is provided to allow
250  * subclasses to perform any initialization required after all
251  * species have been added. For example, it might be used to
252  * resize internal work arrays that must have an entry for
253  * each species. The base class implementation does nothing,
254  * and subclasses that do not require initialization do not
255  * need to overload this method. When importing a CTML phase
256  * description, this method is called just prior to returning
257  * from function importPhase.
258  *
259  * @see importCTML.cpp
260  */
262 {
263  initLengths();
265 }
266 
267 
268 // Initialize lengths of local variables after all species have
269 // been identified.
271 {
272  m_kk = nSpecies();
273  moleFractions_.resize(m_kk);
274 }
275 
276 /*
277  * initThermoXML() (virtual from ThermoPhase)
278  * Import and initialize a ThermoPhase object
279  *
280  * @param phaseNode This object must be the phase node of a
281  * complete XML tree
282  * description of the phase, including all of the
283  * species data. In other words while "phase" must
284  * point to an XML phase object, it must have
285  * sibling nodes "speciesData" that describe
286  * the species in the phase.
287  * @param id ID of the phase. If nonnull, a check is done
288  * to see if phaseNode is pointing to the phase
289  * with the correct id.
290  */
291 void PseudoBinaryVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id)
292 {
293 
294 
295  GibbsExcessVPSSTP::initThermoXML(phaseNode, id);
296 }
297 
298 /**
299  * Format a summary of the mixture state for output.
300  */
301 std::string PseudoBinaryVPSSTP::report(bool show_thermo) const
302 {
303 
304 
305  char p[800];
306  string s = "";
307  try {
308  if (name() != "") {
309  sprintf(p, " \n %s:\n", name().c_str());
310  s += p;
311  }
312  sprintf(p, " \n temperature %12.6g K\n", temperature());
313  s += p;
314  sprintf(p, " pressure %12.6g Pa\n", pressure());
315  s += p;
316  sprintf(p, " density %12.6g kg/m^3\n", density());
317  s += p;
318  sprintf(p, " mean mol. weight %12.6g amu\n", meanMolecularWeight());
319  s += p;
320 
321  doublereal phi = electricPotential();
322  sprintf(p, " potential %12.6g V\n", phi);
323  s += p;
324 
325  size_t kk = nSpecies();
326  vector_fp x(kk);
327  vector_fp molal(kk);
328  vector_fp mu(kk);
329  vector_fp muss(kk);
330  vector_fp acMolal(kk);
331  vector_fp actMolal(kk);
332  getMoleFractions(&x[0]);
333 
334  getChemPotentials(&mu[0]);
335  getStandardChemPotentials(&muss[0]);
336  getActivities(&actMolal[0]);
337 
338 
339  if (show_thermo) {
340  sprintf(p, " \n");
341  s += p;
342  sprintf(p, " 1 kg 1 kmol\n");
343  s += p;
344  sprintf(p, " ----------- ------------\n");
345  s += p;
346  sprintf(p, " enthalpy %12.6g %12.4g J\n",
348  s += p;
349  sprintf(p, " internal energy %12.6g %12.4g J\n",
351  s += p;
352  sprintf(p, " entropy %12.6g %12.4g J/K\n",
354  s += p;
355  sprintf(p, " Gibbs function %12.6g %12.4g J\n",
356  gibbs_mass(), gibbs_mole());
357  s += p;
358  sprintf(p, " heat capacity c_p %12.6g %12.4g J/K\n",
359  cp_mass(), cp_mole());
360  s += p;
361  try {
362  sprintf(p, " heat capacity c_v %12.6g %12.4g J/K\n",
363  cv_mass(), cv_mole());
364  s += p;
365  } catch (CanteraError& err) {
366  err.save();
367  sprintf(p, " heat capacity c_v <not implemented> \n");
368  s += p;
369  }
370  }
371 
372  } catch (CanteraError& err) {
373  err.save();
374  }
375  return s;
376 }
377 
378 
379 }
380