Cantera  2.0
GibbsExcessVPSSTP.h
Go to the documentation of this file.
1 /**
2  * @file GibbsExcessVPSSTP.h
3  * Header for intermediate ThermoPhase object for phases which
4  * employ gibbs excess free energy based 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 activities
11  * based on the molality scale. These include most of the methods for
12  * calculating liquid electrolyte thermodynamics.
13  */
14 /*
15  * Copyright (2006) 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  */
19 #ifndef CT_GIBBSEXCESSVPSSTP_H
20 #define CT_GIBBSEXCESSVPSSTP_H
21 
22 #include "VPStandardStateTP.h"
23 
24 namespace Cantera
25 {
26 
27 /**
28  * @ingroup thermoprops
29  */
30 
31 /*!
32  * GibbsExcessVPSSTP is a derived class of ThermoPhase that handles
33  * variable pressure standard state methods for calculating
34  * thermodynamic properties that are further based on
35  * expressing the Excess Gibbs free energy as a function of
36  * the mole fractions (or pseudo mole fractions) of constituents.
37  * This category is the workhorse for describing molten salts,
38  * solid-phase mixtures of semiconductors, and mixtures of miscible
39  * and semi-miscible compounds.
40  *
41  * It includes
42  * . regular solutions
43  * . Margules expansions
44  * . NTRL equation
45  * . Wilson's equation
46  * . UNIQUAC equation of state.
47  *
48  * This class adds additional functions onto the %ThermoPhase interface
49  * that handles the calculation of the excess Gibbs free energy. The %ThermoPhase
50  * class includes a member function, ThermoPhase::activityConvention()
51  * that indicates which convention the activities are based on. The
52  * default is to assume activities are based on the molar convention.
53  * That default is used here.
54  *
55  * All of the Excess Gibbs free energy formulations in this area employ
56  * symmetrical formulations.
57  *
58  *
59  * Chemical potentials
60  * of species k, \f$ \mu_o \f$, has the following general format:
61  *
62  * \f[
63  * \mu_k = \mu^o_k(T,P) + R T ln( \gamma_k X_k )
64  * \f]
65  *
66  *
67  * where \f$ \gamma_k^{\triangle} \f$ is a molar based activity coefficient for species
68  * \f$k\f$.
69  *
70  * GibbsExcessVPSSTP contains an internal vector with the current mole
71  * fraction vector. That's one of its primary usages. In order to keep the mole fraction
72  * vector constant, all of the setState functions are redesigned at this layer.
73  *
74  *
75  * <H3>
76  * Activity Concentrations: Relationship of %ThermoPhase to %Kinetics Expressions
77  * </H3>
78  *
79  * As explained in a similar discussion in the ThermoPhase class, the actual units used
80  * in kinetics expressions must be specified in the ThermoPhase class for the corresponding
81  * species. These units vary with the field of study. %Cantera uses the concept of
82  * activity concentrations to represent this. Activity concentrations are used directly
83  * in the expressions for kinetics. Standard concentrations are used as the multiplicative
84  * constant that takes the activity of a species and turns it into an activity concentration.
85  * Standard concentrations must not depend on the concentration of the species in the phase.
86  *
87  * Here we set a standard for the specification of the standard concentrations for this class
88  * and all child classes underneath it. We specify here that the standard concentration is
89  * equal to 1 for all species. Therefore, the activities appear directly in kinetics expressions
90  * involving species in underlying %GibbsExcessVPSSTP phases.
91  *
92  * <H3>
93  * SetState Strategy
94  * </H3>
95  *
96  * All setState functions that set the internal state of the ThermoPhase object are
97  * overloaded at this level, so that a current mole fraction vector is maintained within
98  * the object.
99  *
100  *
101  */
103 {
104 
105 public:
106 
107  /// Constructors
108  /*!
109  * This doesn't do much more than initialize constants with
110  * default values for water at 25C. Water molecular weight
111  * comes from the default elements.xml file. It actually
112  * differs slightly from the IAPWS95 value of 18.015268. However,
113  * density conservation and therefore element conservation
114  * is the more important principle to follow.
115  */
117 
118  //! Copy constructor
119  /*!
120  * Note this stuff will not work until the underlying phase
121  * has a working copy constructor
122  *
123  * @param b class to be copied
124  */
126 
127  /// Assignment operator
128  /*!
129  *
130  * @param b class to be copied.
131  */
133 
134  /// Destructor.
135  virtual ~GibbsExcessVPSSTP();
136 
137  //! Duplication routine for objects which inherit from ThermoPhase.
138  /*!
139  * This virtual routine can be used to duplicate thermophase objects
140  * inherited from ThermoPhase even if the application only has
141  * a pointer to ThermoPhase to work with.
142  */
143  virtual ThermoPhase* duplMyselfAsThermoPhase() const;
144 
145  /**
146  *
147  * @name Utilities
148  * @{
149  */
150 
151 
152  //! Equation of state type flag.
153  /*!
154  * The ThermoPhase base class returns
155  * zero. Subclasses should define this to return a unique
156  * non-zero value. Known constants defined for this purpose are
157  * listed in mix_defs.h. The MolalityVPSSTP class also returns
158  * zero, as it is a non-complete class.
159  */
160  virtual int eosType() const;
161 
162 
163 
164  /**
165  * @}
166  * @name Molar Thermodynamic Properties
167  * @{
168  */
169 
170 
171 
172 
173 
174  /**
175  * @}
176  * @name Mechanical Properties
177  * @{
178  */
179 
180  //! Set the internally stored pressure (Pa) at constant
181  //! temperature and composition
182  /*!
183  * This method sets the pressure within the object.
184  * The water model is a completely compressible model.
185  * Also, the dielectric constant is pressure dependent.
186  *
187  * @param p input Pressure (Pa)
188  *
189  * @todo Implement a variable pressure capability
190  */
191  virtual void setPressure(doublereal p);
192 
193 protected:
194 
195  /**
196  * Calculate the density of the mixture using the partial
197  * molar volumes and mole fractions as input
198  *
199  * The formula for this is
200  *
201  * \f[
202  * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}}
203  * \f]
204  *
205  * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are
206  * the molecular weights, and \f$V_k\f$ are the pure species
207  * molar volumes.
208  *
209  * Note, the basis behind this formula is that in an ideal
210  * solution the partial molar volumes are equal to the pure
211  * species molar volumes. We have additionally specified
212  * in this class that the pure species molar volumes are
213  * independent of temperature and pressure.
214  *
215  * NOTE: This is a non-virtual function, which is not a
216  * member of the ThermoPhase base class.
217  */
218  void calcDensity();
219 
220 public:
221  /**
222  * @}
223  * @name Potential Energy
224  *
225  * Species may have an additional potential energy due to the
226  * presence of external gravitation or electric fields. These
227  * methods allow specifying a potential energy for individual
228  * species.
229  * @{
230  */
231 
232  /**
233  * @}
234  * @name Activities, Standard States, and Activity Concentrations
235  *
236  * The activity \f$a_k\f$ of a species in solution is
237  * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
238  * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
239  * the chemical potential at unit activity, which depends only
240  * on temperature and pressure.
241  * @{
242  */
243 
244  //! This method returns an array of generalized concentrations
245  /*!
246  * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k /
247  * C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration
248  * defined below and \f$ a_k \f$ are activities used in the
249  * thermodynamic functions. These activity (or generalized)
250  * concentrations are used by kinetics manager classes to compute the forward and
251  * reverse rates of elementary reactions. Note that they may
252  * or may not have units of concentration --- they might be
253  * partial pressures, mole fractions, or surface coverages,
254  * for example.
255  *
256  * @param c Output array of generalized concentrations. The
257  * units depend upon the implementation of the
258  * reaction rate expressions within the phase.
259  */
260  virtual void getActivityConcentrations(doublereal* c) const;
261 
262 
263 
264  /**
265  * The standard concentration \f$ C^0_k \f$ used to normalize
266  * the generalized concentration. In many cases, this quantity
267  * will be the same for all species in a phase - for example,
268  * for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this
269  * reason, this method returns a single value, instead of an
270  * array. However, for phases in which the standard
271  * concentration is species-specific (e.g. surface species of
272  * different sizes), this method may be called with an
273  * optional parameter indicating the species.
274  *
275  * The standard concentration for defaulted to 1. In other words
276  * the activity concentration is assumed to be 1.
277  *
278  * @param k species index. Defaults to zero.
279  */
280  virtual doublereal standardConcentration(size_t k=0) const;
281 
282  /**
283  * Returns the natural logarithm of the standard
284  * concentration of the kth species
285  *
286  * @param k species index
287  */
288  virtual doublereal logStandardConc(size_t k=0) const;
289 
290  /**
291  * Returns the units of the standard and generalized
292  * concentrations Note they have the same units, as their
293  * ratio is defined to be equal to the activity of the kth
294  * species in the solution, which is unitless.
295  *
296  * This routine is used in print out applications where the
297  * units are needed. Usually, MKS units are assumed throughout
298  * the program and in the XML input files.
299  *
300  * @param uA Output vector containing the units
301  * uA[0] = kmol units - default = 1
302  * uA[1] = m units - default = -nDim(), the number of spatial
303  * dimensions in the Phase class.
304  * uA[2] = kg units - default = 0;
305  * uA[3] = Pa(pressure) units - default = 0;
306  * uA[4] = Temperature units - default = 0;
307  * uA[5] = time units - default = 0
308  * @param k species index. Defaults to 0.
309  * @param sizeUA output int containing the size of the vector.
310  * Currently, this is equal to 6.
311  */
312  virtual void getUnitsStandardConc(double* uA, int k = 0,
313  int sizeUA = 6) const;
314 
315 
316  //! Get the array of non-dimensional activities (molality
317  //! based for this class and classes that derive from it) at
318  //! the current solution temperature, pressure, and solution concentration.
319  /*!
320  * \f[
321  * a_i^\triangle = \gamma_k^{\triangle} \frac{m_k}{m^\triangle}
322  * \f]
323  *
324  * This function must be implemented in derived classes.
325  *
326  * @param ac Output vector of molality-based activities. Length: m_kk.
327  */
328  virtual void getActivities(doublereal* ac) const;
329 
330  //! Get the array of non-dimensional molar-based activity coefficients at
331  //! the current solution temperature, pressure, and solution concentration.
332  /*!
333  * @param ac Output vector of activity coefficients. Length: m_kk.
334  */
335  virtual void getActivityCoefficients(doublereal* ac) const;
336 
337  //! Get the array of temperature derivatives of the log activity coefficients
338  /*!
339  * This function is a virtual class, but it first appears in GibbsExcessVPSSTP
340  * class and derived classes from GibbsExcessVPSSTP.
341  *
342  * units = 1/Kelvin
343  *
344  * @param dlnActCoeffdT Output vector of temperature derivatives of the
345  * log Activity Coefficients. length = m_kk
346  */
347  virtual void getdlnActCoeffdT(doublereal* dlnActCoeffdT) const {
348  err("getdlnActCoeffdT");
349  }
350 
351  //! Get the array of derivatives of the log activity coefficients with respect to the log of the species mole numbers
352  /*!
353  * Implementations should take the derivative of the logarithm of the activity coefficient with respect to a
354  * species log mole number (with all other species mole numbers held constant). The default treatment in the
355  * %ThermoPhase object is to set this vector to zero.
356  *
357  * units = 1 / kmol
358  *
359  * dlnActCoeffdlnN[ ld * k + m] will contain the derivative of log act_coeff for the <I>m</I><SUP>th</SUP>
360  * species with respect to the number of moles of the <I>k</I><SUP>th</SUP> species.
361  *
362  * \f[
363  * \frac{d \ln(\gamma_m) }{d \ln( n_k ) }\Bigg|_{n_i}
364  * \f]
365  *
366  * @param ld Number of rows in the matrix
367  * @param dlnActCoeffdlnN Output vector of derivatives of the
368  * log Activity Coefficients. length = m_kk * m_kk
369  */
370  virtual void getdlnActCoeffdlnN(const size_t ld, doublereal* const dlnActCoeffdlnN) {
371  err(" getdlnActCoeffdlnN: nonzero and nonimplemented");
372  }
373 
374  //! Get the array of log concentration-like derivatives of the
375  //! log activity coefficients
376  /*!
377  * This function is a virtual method. For ideal mixtures
378  * (unity activity coefficients), this can return zero.
379  * Implementations should take the derivative of the
380  * logarithm of the activity coefficient with respect to the
381  * logarithm of the concentration-like variable (i.e. number of moles in
382  * in a unit volume. ) that represents the standard state.
383  * This quantity is to be used in conjunction with derivatives of
384  * that concentration-like variable when the derivative of the chemical
385  * potential is taken.
386  *
387  * units = dimensionless
388  *
389  * @param dlnActCoeffdlnX Output vector of derivatives of the
390  * log Activity Coefficients. length = m_kk
391  */
392  virtual void getdlnActCoeffdlnX(doublereal* dlnActCoeffdlnX) const {
393  err("getdlnActCoeffdlnX");
394  }
395 
396 
397  //@}
398  /// @name Partial Molar Properties of the Solution
399  //@{
400 
401 
402  /**
403  * Get the species electrochemical potentials.
404  * These are partial molar quantities.
405  * This method adds a term \f$ Fz_k \phi_k \f$ to the
406  * to each chemical potential.
407  *
408  * Units: J/kmol
409  *
410  * @param mu output vector containing the species electrochemical potentials.
411  * Length: m_kk.
412  */
413  void getElectrochemPotentials(doublereal* mu) const;
414 
415  //! Return an array of partial molar volumes for the
416  //! species in the mixture. Units: m^3/kmol.
417  /*!
418  * Frequently, for this class of thermodynamics representations,
419  * the excess Volume due to mixing is zero. Here, we set it as
420  * a default. It may be overridden in derived classes.
421  *
422  * @param vbar Output vector of species partial molar volumes.
423  * Length = m_kk. units are m^3/kmol.
424  */
425  virtual void getPartialMolarVolumes(doublereal* vbar) const;
426 
427  //@}
428  /// @name Properties of the Standard State of the Species in the Solution
429  //@{
430 
431 
432 
433  //@}
434  /// @name Thermodynamic Values for the Species Reference States
435  //@{
436 
437 
438  ///////////////////////////////////////////////////////
439  //
440  // The methods below are not virtual, and should not
441  // be overloaded.
442  //
443  //////////////////////////////////////////////////////
444 
445  /**
446  * @name Specific Properties
447  * @{
448  */
449 
450 
451  /**
452  * @name Setting the State
453  *
454  * These methods set all or part of the thermodynamic
455  * state.
456  * @{
457  */
458 
459 
460  //! Set the temperature (K) and pressure (Pa)
461  /*!
462  * Set the temperature and pressure.
463  *
464  * @param t Temperature (K)
465  * @param p Pressure (Pa)
466  */
467  virtual void setState_TP(doublereal t, doublereal p);
468 
469  //@}
470 
471  /**
472  * @name Chemical Equilibrium
473  * Routines that implement the Chemical equilibrium capability
474  * for a single phase, based on the element-potential method.
475  * @{
476  */
477 
478  /**
479  * Set the mass fractions to the specified values, and then
480  * normalize them so that they sum to 1.0.
481  * @param y Array of unnormalized mass fraction values (input).
482  * Must have a length greater than or equal to the number of
483  * species.
484  *
485  * @param y Input vector of mass fractions.
486  * Length is m_kk.
487  */
488  virtual void setMassFractions(const doublereal* const y);
489 
490  /**
491  * Set the mass fractions to the specified values without
492  * normalizing. This is useful when the normalization
493  * condition is being handled by some other means, for example
494  * by a constraint equation as part of a larger set of
495  * equations.
496  *
497  * @param y Input vector of mass fractions.
498  * Length is m_kk.
499  */
500  virtual void setMassFractions_NoNorm(const doublereal* const y);
501 
502 
503  /**
504  * Set the mole fractions to the specified values, and then
505  * normalize them so that they sum to 1.0.
506  * @param x Array of unnormalized mole fraction values (input).
507  * Must have a length greater than or equal to the number of
508  * species.
509  *
510  * @param x Input vector of mole fractions.
511  * Length is m_kk.
512  */
513  virtual void setMoleFractions(const doublereal* const x);
514 
515  /**
516  * Set the mole fractions to the specified values without
517  * normalizing. This is useful when the normalization
518  * condition is being handled by some other means, for example
519  * by a constraint equation as part of a larger set of
520  * equations.
521  *
522  * @param x Input vector of mole fractions.
523  * Length is m_kk.
524  */
525  virtual void setMoleFractions_NoNorm(const doublereal* const x);
526 
527  /**
528  * Set the concentrations to the specified values within the
529  * phase.
530  *
531  * @param c The input vector to this routine is in dimensional
532  * units. For volumetric phases c[k] is the
533  * concentration of the kth species in kmol/m3.
534  * For surface phases, c[k] is the concentration
535  * in kmol/m2. The length of the vector is the number
536  * of species in the phase.
537  */
538  virtual void setConcentrations(const doublereal* const c);
539 
540  //@}
541 
542 
543 
544  /// The following methods are used in the process of constructing
545  /// the phase and setting its parameters from a specification in an
546  /// input file. They are not normally used in application programs.
547  /// To see how they are used, see files importCTML.cpp and
548  /// ThermoFactory.cpp.
549 
550 
551  /*!
552  * @internal Initialize. This method is provided to allow
553  * subclasses to perform any initialization required after all
554  * species have been added. For example, it might be used to
555  * resize internal work arrays that must have an entry for
556  * each species. The base class implementation does nothing,
557  * and subclasses that do not require initialization do not
558  * need to overload this method. When importing a CTML phase
559  * description, this method is called just prior to returning
560  * from function importPhase.
561  *
562  * @see importCTML.cpp
563  */
564  virtual void initThermo();
565 
566 
567 private:
568 
569  //! Initialize lengths of local variables after all species have
570  //! been identified.
571  void initLengths();
572 
573  //! Error function
574  /*!
575  * Print an error string and exit
576  *
577  * @param msg Message to be printed
578  */
579  doublereal err(std::string msg) const;
580 
581 protected:
582 
583  //! utility routine to check mole fraction sum
584  /*!
585  * @param x vector of mole fractions.
586  */
587  double checkMFSum(const doublereal* const x) const;
588 
589 protected:
590 
591  // HKM get rid of _Scaled_ prefix
592 
593  //! Storage for the current values of the mole fractions of the species
594  /*!
595  * This vector is kept up-to-date when the setState functions are called.
596  * Therefore, it may be considered to be an independent variable.
597  *
598  * Note in order to do this, the setState functions are redefined to always
599  * keep this vector current.
600  */
601  mutable std::vector<doublereal> moleFractions_;
602 
603  //! Storage for the current values of the activity coefficients of the
604  //! species
605  mutable std::vector<doublereal> lnActCoeff_Scaled_;
606 
607  //! Storage for the current derivative values of the
608  //! gradients with respect to temperature of the
609  //! log of the activity coefficients of the species
610  mutable std::vector<doublereal> dlnActCoeffdT_Scaled_;
611 
612  //! Storage for the current derivative values of the
613  //! gradients with respect to temperature of the
614  //! log of the activity coefficients of the species
615  mutable std::vector<doublereal> d2lnActCoeffdT2_Scaled_;
616 
617  //! Storage for the current derivative values of the
618  //! gradients with respect to logarithm of the mole fraction of the
619  //! log of the activity coefficients of the species @deprecated
620  mutable std::vector<doublereal> dlnActCoeffdlnN_diag_;
621 
622  //! Storage for the current derivative values of the
623  //! gradients with respect to logarithm of the mole fraction of the
624  //! log of the activity coefficients of the species @deprecated
625  mutable std::vector<doublereal> dlnActCoeffdlnX_diag_;
626 
627  //! Storage for the current derivative values of the gradients with respect to logarithm of the species mole number of the
628  //! log of the activity coefficients of the species
629  /*!
630  * dlnActCoeffdlnN_(k, m) is the derivative of ln(gamma_k) wrt ln mole number of species m
631  */
633 
634  //! Temporary storage space that is fair game
635  mutable std::vector<doublereal> m_pp;
636 
637 };
638 
639 
640 }
641 
642 #endif
643 
644 
645 
646 
647