Cantera  2.0
GibbsExcessVPSSTP.cpp
Go to the documentation of this file.
1 /**
2  * @file GibbsExcessVPSSTP.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::GibbsExcessVPSSTP GibbsExcessVPSSTP\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 <iomanip>
23 
24 using namespace std;
25 
26 namespace Cantera
27 {
28 
29 /*
30  * Default constructor.
31  *
32  */
33 GibbsExcessVPSSTP::GibbsExcessVPSSTP() :
35  moleFractions_(0),
36  lnActCoeff_Scaled_(0),
37  dlnActCoeffdT_Scaled_(0),
38  d2lnActCoeffdT2_Scaled_(0),
39  dlnActCoeffdlnN_diag_(0),
40  dlnActCoeffdlnX_diag_(0),
41  dlnActCoeffdlnN_(0,0),
42  m_pp(0)
43 {
44 }
45 
46 /*
47  * Copy Constructor:
48  *
49  * Note this stuff will not work until the underlying phase
50  * has a working copy constructor
51  */
54  moleFractions_(0),
55  lnActCoeff_Scaled_(0),
56  dlnActCoeffdT_Scaled_(0),
57  d2lnActCoeffdT2_Scaled_(0),
58  dlnActCoeffdlnN_diag_(0),
59  dlnActCoeffdlnX_diag_(0),
60  dlnActCoeffdlnN_(0,0),
61  m_pp(0)
62 {
64 }
65 
66 /*
67  * operator=()
68  *
69  * Note this stuff will not work until the underlying phase
70  * has a working assignment operator
71  */
74 {
75  if (&b == this) {
76  return *this;
77  }
78 
80 
88  m_pp = b.m_pp;
89 
90  return *this;
91 }
92 
93 /*
94  *
95  * ~GibbsExcessVPSSTP(): (virtual)
96  *
97  * Destructor: does nothing:
98  *
99  */
101 {
102 }
103 
104 /*
105  * This routine duplicates the current object and returns
106  * a pointer to ThermoPhase.
107  */
110 {
111  GibbsExcessVPSSTP* mtp = new GibbsExcessVPSSTP(*this);
112  return (ThermoPhase*) mtp;
113 }
114 
115 /*
116  * -------------- Utilities -------------------------------
117  */
118 
119 void GibbsExcessVPSSTP::setMassFractions(const doublereal* const y)
120 {
123 }
124 
125 void GibbsExcessVPSSTP::setMassFractions_NoNorm(const doublereal* const y)
126 {
129 }
130 
131 void GibbsExcessVPSSTP::setMoleFractions(const doublereal* const x)
132 {
135 }
136 
137 void GibbsExcessVPSSTP::setMoleFractions_NoNorm(const doublereal* const x)
138 {
141 }
142 
143 
144 void GibbsExcessVPSSTP::setConcentrations(const doublereal* const c)
145 {
148 }
149 
150 
151 // Equation of state type flag.
152 /*
153  * The ThermoPhase base class returns
154  * zero. Subclasses should define this to return a unique
155  * non-zero value. Known constants defined for this purpose are
156  * listed in mix_defs.h. The GibbsExcessVPSSTP class also returns
157  * zero, as it is a non-complete class.
158  */
160 {
161  return 0;
162 }
163 
164 
165 
166 /*
167  * ------------ Molar Thermodynamic Properties ----------------------
168  */
169 
170 /*
171  *
172  * ------------ Mechanical Properties ------------------------------
173  *
174  */
175 
176 /*
177  * Set the pressure at constant temperature. Units: Pa.
178  * This method sets a constant within the object.
179  * The mass density is not a function of pressure.
180  */
182 {
183  setState_TP(temperature(), p);
184 }
185 
187 {
188  doublereal* vbar = NULL;
189  vbar = new doublereal[m_kk];
190  // double *vbar = &m_pp[0];
192 
193  doublereal vtotal = 0.0;
194  for (size_t i = 0; i < m_kk; i++) {
195  vtotal += vbar[i] * moleFractions_[i];
196  }
197  doublereal dd = meanMolecularWeight() / vtotal;
198  Phase::setDensity(dd);
199  delete [] vbar;
200 }
201 
202 void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p)
203 {
205  /*
206  * Store the current pressure
207  */
208  m_Pcurrent = p;
209  /*
210  * update the standard state thermo
211  * -> This involves calling the water function and setting the pressure
212  */
214 
215  /*
216  * Calculate the partial molar volumes, and then the density of the fluid
217  */
218  calcDensity();
219 }
220 
221 
222 
223 /*
224  * - Activities, Standard States, Activity Concentrations -----------
225  */
227 {
228  getActivities(c);
229 }
230 
231 
233 {
234  return 1.0;
235 }
236 
237 doublereal GibbsExcessVPSSTP::logStandardConc(size_t k) const
238 {
239  return 0.0;
240 }
241 
242 void GibbsExcessVPSSTP::getActivities(doublereal* ac) const
243 {
246  for (size_t k = 0; k < m_kk; k++) {
247  ac[k] *= moleFractions_[k];
248  }
249 }
250 
251 void GibbsExcessVPSSTP::getActivityCoefficients(doublereal* const ac) const
252 {
253 
255 
256  // Protect against roundoff when taking exponentials
257  for (size_t k = 0; k < m_kk; k++) {
258  if (ac[k] > 700.) {
259  ac[k] = exp(700.0);
260  } else if (ac[k] < -700.) {
261  ac[k] = exp(-700.0);
262  } else {
263  ac[k] = exp(ac[k]);
264  }
265  }
266 }
267 //====================================================================================================================
268 
270 {
271  getChemPotentials(mu);
272  double ve = Faraday * electricPotential();
273  for (size_t k = 0; k < m_kk; k++) {
274  mu[k] += ve*charge(k);
275  }
276 }
277 
278 /*
279  * ------------ Partial Molar Properties of the Solution ------------
280  */
281 
282 // Return an array of partial molar volumes for the
283 // species in the mixture. Units: m^3/kmol.
284 /*
285  * Frequently, for this class of thermodynamics representations,
286  * the excess Volume due to mixing is zero. Here, we set it as
287  * a default. It may be overridden in derived classes.
288  *
289  * @param vbar Output vector of species partial molar volumes.
290  * Length = m_kk. units are m^3/kmol.
291  */
292 void GibbsExcessVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
293 {
294  /*
295  * Get the standard state values in m^3 kmol-1
296  */
297  getStandardVolumes(vbar);
298 }
299 
300 
301 
302 doublereal GibbsExcessVPSSTP::err(std::string msg) const
303 {
304  throw CanteraError("GibbsExcessVPSSTP","Base class method "
305  +msg+" called. Equation of state type: "+int2str(eosType()));
306  return 0;
307 }
308 
309 
310 
311 double GibbsExcessVPSSTP::checkMFSum(const doublereal* const x) const
312 {
313  doublereal norm = accumulate(x, x + m_kk, 0.0);
314  if (fabs(norm - 1.0) > 1.0E-9) {
315  throw CanteraError("GibbsExcessVPSSTP::checkMFSum",
316  "(MF sum - 1) exceeded tolerance of 1.0E-9:" + fp2str(norm));
317  }
318  return norm;
319 }
320 
321 /*
322  * Returns the units of the standard and general concentrations
323  * Note they have the same units, as their divisor is
324  * defined to be equal to the activity of the kth species
325  * in the solution, which is unitless.
326  *
327  * This routine is used in print out applications where the
328  * units are needed. Usually, MKS units are assumed throughout
329  * the program and in the XML input files.
330  *
331  * On return uA contains the powers of the units (MKS assumed)
332  * of the standard concentrations and generalized concentrations
333  * for the kth species.
334  *
335  * uA[0] = kmol units - default = 1
336  * uA[1] = m units - default = -nDim(), the number of spatial
337  * dimensions in the Phase class.
338  * uA[2] = kg units - default = 0;
339  * uA[3] = Pa(pressure) units - default = 0;
340  * uA[4] = Temperature units - default = 0;
341  * uA[5] = time units - default = 0
342  */
343 void GibbsExcessVPSSTP::getUnitsStandardConc(double* uA, int k, int sizeUA) const
344 {
345  for (int i = 0; i < sizeUA; i++) {
346  if (i == 0) {
347  uA[0] = 0.0;
348  }
349  if (i == 1) {
350  uA[1] = 0.0;
351  }
352  if (i == 2) {
353  uA[2] = 0.0;
354  }
355  if (i == 3) {
356  uA[3] = 0.0;
357  }
358  if (i == 4) {
359  uA[4] = 0.0;
360  }
361  if (i == 5) {
362  uA[5] = 0.0;
363  }
364  }
365 }
366 
367 
368 /*
369  * @internal Initialize. This method is provided to allow
370  * subclasses to perform any initialization required after all
371  * species have been added. For example, it might be used to
372  * resize internal work arrays that must have an entry for
373  * each species. The base class implementation does nothing,
374  * and subclasses that do not require initialization do not
375  * need to overload this method. When importing a CTML phase
376  * description, this method is called just prior to returning
377  * from function importPhase.
378  *
379  * @see importCTML.cpp
380  */
382 {
383  initLengths();
386 }
387 
388 
389 // Initialize lengths of local variables after all species have
390 // been identified.
392 {
393  m_kk = nSpecies();
394  moleFractions_.resize(m_kk);
395  lnActCoeff_Scaled_.resize(m_kk);
396  dlnActCoeffdT_Scaled_.resize(m_kk);
398  dlnActCoeffdlnX_diag_.resize(m_kk);
399  dlnActCoeffdlnN_diag_.resize(m_kk);
401  m_pp.resize(m_kk);
402 }
403 
404 
405 }
406