MolalityVPSSTP.h

Go to the documentation of this file.
00001 /**
00002  *  @file MolalityVPSSTP.h
00003  *   Header for intermediate ThermoPhase object for phases which
00004  *   employ molality based activity coefficient formulations
00005  *  (see \ref thermoprops 
00006  * and class \link Cantera::MolalityVPSSTP MolalityVPSSTP\endlink).
00007  *
00008  * Header file for a derived class of ThermoPhase that handles
00009  * variable pressure standard state methods for calculating
00010  * thermodynamic properties that are further based upon activities
00011  * based on the molality scale.  These include most of the methods for
00012  * calculating liquid electrolyte thermodynamics.
00013  */
00014 /*
00015  * Copywrite (2006) Sandia Corporation. Under the terms of 
00016  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
00017  * U.S. Government retains certain rights in this software.
00018  */
00019 /*
00020  *  $Id: MolalityVPSSTP.h 279 2009-12-05 19:08:43Z hkmoffa $
00021  */
00022 
00023 #ifndef CT_MOLALITYVPSSTP_H
00024 #define CT_MOLALITYVPSSTP_H
00025 
00026 #include "VPStandardStateTP.h"
00027 
00028 namespace Cantera {
00029 
00030   /**
00031    * @ingroup thermoprops
00032    */
00033 
00034   /*!
00035    * MolalityVPSSTP is a derived class of ThermoPhase that handles
00036    * variable pressure standard state methods for calculating
00037    * thermodynamic properties that are further based on
00038    * molality-scaled activities.
00039    * This category incorporates most of the methods
00040    * for calculating liquid electrolyte thermodynamics that have been
00041    * developed since the 1970's.
00042    *
00043    * This class adds additional functions onto the %ThermoPhase interface
00044    * that handle molality based standard states. The %ThermoPhase
00045    * class includes a member function, ThermoPhase::activityConvention() 
00046    * that indicates which convention the activities are based on. The
00047    * default is to assume activities are based on the molar convention. 
00048    * However, classes which derive from the MolalityVPSSTP class return
00049    * <b>cAC_CONVENTION_MOLALITY</b> from this member function.
00050    *
00051    * The molality of a solute, \f$ m_i \f$, is defined as
00052    *
00053    * \f[
00054    *     m_i = \frac{n_i}{\tilde{M}_o n_o} 
00055    * \f]
00056    * where
00057    * \f[
00058    *     \tilde{M}_o = \frac{M_o}{1000} 
00059    * \f]
00060    *    
00061    * where \f$ M_o \f$ is the molecular weight of the solvent. The molality
00062    * has units of gmol kg<SUP>-1</SUP>. For the solute, the molality may be
00063    * considered as the amount of gmol's of solute per kg of solvent, a natural
00064    * experimental quantity.
00065    *
00066    * The formulas for calculating mole fractions if given the molalities of
00067    * the solutes is stated below. First calculate \f$ L^{sum} \f$, an intermediate
00068    * quantity.
00069    *
00070    *   \f[
00071    *       L^{sum} = \frac{1}{\tilde{M}_o X_o} = \frac{1}{\tilde{M}_o} + \sum_{i\ne o} m_i
00072    *  \f]
00073    *  Then,
00074    *  \f[
00075    *       X_o =   \frac{1}{\tilde{M}_o L^{sum}}  
00076    *  \f]
00077    *  \f[
00078    *       X_i =   \frac{m_i}{L^{sum}}  
00079    *  \f]
00080    * where \f$ X_o \f$ is the mole fraction of solvent, and \f$ X_o \f$ is the
00081    * mole fraction of solute <I>i</I>. Thus, the molality scale and the mole fraction
00082    * scale offer a one-to-one mapping between each other, except in the limit
00083    * of a zero solvent mole fraction.
00084    *
00085    * The standard states for thermodynamic objects that derive from <b>MolalityVPSSTP</b>
00086    * are on the unit molality basis. Chemical potentials
00087    * of the solutes,  \f$ \mu_k \f$, and the solvent, \f$ \mu_o \f$, which are based 
00088    * on the molality form, have the following general format:
00089    *
00090    * \f[
00091    *    \mu_k = \mu^{\triangle}_k(T,P) + R T ln(\gamma_k^{\triangle} \frac{m_k}{m^\triangle}) 
00092    * \f]
00093    * \f[
00094    *    \mu_o = \mu^o_o(T,P) + RT ln(a_o) 
00095    * \f]
00096    *
00097    * where \f$ \gamma_k^{\triangle} \f$ is the molality based activity coefficient for species
00098    * \f$k\f$. 
00099    *
00100    * The chemical potential of the solvent is thus expressed in a different format
00101    * than the chemical potential of the solutes. Additionally, the activity of the
00102    * solvent, \f$ a_o \f$, is further reexpressed in terms of an osmotic coefficient,
00103    * \f$ \phi \f$.
00104    *   \f[
00105    *       \phi = \frac{- ln(a_o)}{\tilde{M}_o \sum_{i \ne o} m_i}
00106    *   \f]
00107    * 
00108    *  MolalityVPSSTP::osmoticCoefficient() returns the value of \f$ \phi \f$. 
00109    *  Note there  are a few of definitions of the osmotic coefficient floating
00110    *  around. We use the one defined in 
00111    *  (Activity Coefficients in Electrolyte Solutions, K. S. Pitzer
00112    *   CRC Press, Boca Raton, 1991, p. 85, Eqn. 28). This definition is most clearly
00113    *  related to theoretical calculation.
00114    *
00115    * The molar-based activity coefficients \f$ \gamma_k \f$ may be calculated 
00116    * from the molality-based
00117    * activity coefficients, \f$ \gamma_k^\triangle \f$ by the following 
00118    * formula.
00119    * \f[
00120    *     \gamma_k = \frac{\gamma_k^\triangle}{X_o}
00121    * \f]
00122    * For purposes of establishing a convention, the molar activity coefficient of the
00123    * solvent is set equal to the molality-based activity coefficient of the
00124    * solvent: 
00125    * \f[
00126    *     \gamma_o = \gamma_o^\triangle
00127    * \f]
00128    *
00129    *  The molality-based and molarity-based standard states may be related to one
00130    *  another by the following formula.
00131    *
00132    * \f[
00133    *   \mu_k^\triangle(T,P) = \mu_k^o(T,P) + R T \ln(\tilde{M}_o m^\triangle)
00134    * \f]
00135    *
00136    *  An important convention is followed in all routines that derive from <b>%MolalityVPSSTP</b>.
00137    *  Standard state thermodynamic functions and reference state thermodynamic functions
00138    *  return the molality-based quantities. Also all functions which return 
00139    *  activities return the molality-based activities. The reason for this convention
00140    *  has been discussed in supporting memos. However, it's important because the
00141    *  term in the equation above is non-trivial. For example it's equal
00142    *  to 2.38 kcal gmol<SUP>-1</SUP> for water at 298 K. 
00143    *
00144    *
00145    * In order to prevent a singularity, this class includes the concept of a minimum
00146    * value for the solvent mole fraction. All calculations involving the formulation
00147    * of activity coefficients and other non-ideal solution behavior adhere to 
00148    * this concept of a minimul value for the solvent mole fraction. This makes sense
00149    * because these solution behavior were all designed and measured far away from 
00150    * the zero solvent singularity condition and are not applicable in that limit. 
00151    *   
00152    *
00153    * This objects add a layer that supports molality. It inherits from VPStandardStateTP.
00154    *
00155    *  All objects that derive from this are assumed to have molality based standard states.
00156    *
00157    *  Molality based activity coefficients are scaled according to the current
00158    *  pH scale. See the Eq3/6 manual for details.
00159    *
00160    *  Activity coefficients for species k may be altered between scales s1 to s2
00161    *  using the following formula
00162    *
00163    *   \f[
00164    *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) 
00165    *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
00166    *   \f]
00167    *
00168    *  where j is any one species. For the NBS scale, j is equal to the Cl- species
00169    *  and 
00170    *
00171    *  \f[
00172    *       ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}}
00173    *  \f]
00174    *
00175    *  The Pitzer scale doesn't actually change anything. The pitzer scale is defined
00176    *  as the raw unscaled activity coefficients produced by the underlying objects.
00177    *
00178    *  <H3> SetState Strategy  </H3>
00179    *
00180    *   The MolalityVPSSTP object does not have a setState strategy concerning the
00181    *   molalities. It does not keep track of whether the molalities have changed.
00182    *   It's strictly an interfacial layer that writes the current mole fractions to the
00183    *   State object. When molalities are needed it recalculates the molalities from 
00184    *   the State object's mole fraction vector.
00185    *
00186    *
00187    * @todo Make two solvent minimum fractions. One would be for calculation of the non-ideal
00188    *       factors. The other one would be for purposes of stoichiometry evaluation. the
00189    *       stoichiometry evaluation one would be a 1E-13 limit. Anything less would create
00190    *       problems with roundoff error.
00191    *
00192    */
00193   class MolalityVPSSTP : public VPStandardStateTP  {
00194 
00195   public:
00196         
00197     /// Constructors 
00198     /*!
00199      * This doesn't do much more than initialize constants with
00200      * default values for water at 25C. Water molecular weight 
00201      * comes from the default elements.xml file. It actually
00202      * differs slightly from the IAPWS95 value of 18.015268. However,
00203      * density conservation and therefore element conservation
00204      * is the more important principle to follow.
00205      */
00206     MolalityVPSSTP();
00207 
00208     //! Copy constructor
00209     /*!
00210      *  Note this stuff will not work until the underlying phase
00211      *  has a working copy constructor
00212      *
00213      * @param b class to be copied
00214      */
00215     MolalityVPSSTP(const MolalityVPSSTP &b);
00216 
00217     /// Assignment operator
00218     /*!
00219      *  Note this stuff will not work until the underlying phase
00220      *  has a working assignment operator
00221      *
00222      * @param b class to be copied.
00223      */
00224     MolalityVPSSTP& operator=(const     MolalityVPSSTP&b);
00225 
00226     /// Destructor. 
00227     virtual ~MolalityVPSSTP();
00228 
00229     //! Duplication routine for objects which inherit from  ThermoPhase.
00230     /*!
00231      *  This virtual routine can be used to duplicate thermophase objects
00232      *  inherited from ThermoPhase even if the application only has
00233      *  a pointer to ThermoPhase to work with.
00234      */
00235     virtual ThermoPhase *duplMyselfAsThermoPhase() const;
00236     
00237     /**
00238      *   
00239      * @name  Utilities  
00240      * @{
00241      */
00242 
00243    
00244     //! Equation of state type flag.
00245     /*!
00246      * The ThermoPhase base class returns
00247      * zero. Subclasses should define this to return a unique
00248      * non-zero value. Known constants defined for this purpose are
00249      * listed in mix_defs.h. The MolalityVPSSTP class also returns
00250      * zero, as it is a non-complete class.
00251      */
00252     virtual int eosType() const;
00253 
00254     //! Set the pH scale, which determines the scale for single-ion activity 
00255     //! coefficients.
00256     /*!
00257      *  Single ion activity coefficients are not unique in terms of the
00258      *  representing actual measureable quantities.
00259      *
00260      * @param pHscaleType  Integer representing the pHscale 
00261      */
00262     void setpHScale(const int pHscaleType);
00263 
00264     //! Reports the pH scale, which determines the scale for single-ion activity 
00265     //! coefficients.
00266     /*!
00267      *  Single ion activity coefficients are not unique in terms of the
00268      *  representing actual measureable quantities. 
00269      *
00270      * @return Return the pHscale type
00271      */
00272     int pHScale() const;
00273 
00274     /**
00275      * @} 
00276      * @name  Molar Thermodynamic Properties 
00277      * @{
00278      */
00279 
00280 
00281     /**
00282      * @}
00283      * @name Utilities for Solvent ID and Molality
00284      * @{
00285      */
00286 
00287     /**
00288      * This routine sets the index number of the solvent for
00289      * the phase.
00290      *  
00291      *  Note, having a solvent 
00292      *   is a precursor to many things having to do with molality.
00293      *
00294      * @param k the solvent index number
00295      */
00296     void setSolvent(int k);
00297 
00298     /**
00299      * Sets the minimum mole fraction in the molality formulation.
00300      * Note the molality formulation is singular in the limit that
00301      * the solvent mole fraction goes to zero. Numerically, how
00302      * this limit is treated and resolved is an ongoing issue within
00303      * Cantera.
00304      *
00305      * @param xmolSolventMIN  Input double containing the minimum mole fraction
00306      */
00307     void setMoleFSolventMin(doublereal xmolSolventMIN);
00308 
00309     //! Returns the solvent index.
00310     int solventIndex() const;
00311 
00312     /**
00313      * Returns the minimum mole fraction in the molality 
00314      * formulation.
00315      */
00316     doublereal moleFSolventMin() const;
00317 
00318     //! Calculates the molality of all species and stores the result internally.
00319     /*!
00320      *   We calculate the vector of molalities of the species
00321      *   in the phase and store the result internally:
00322      *   \f[
00323      *     m_i = \frac{X_i}{1000 * M_o * X_{o,p}}
00324      *   \f]
00325      *    where 
00326      *    - \f$ M_o \f$ is the molecular weight of the solvent
00327      *    - \f$ X_o \f$ is the mole fraction of the solvent
00328      *    - \f$ X_i \f$ is the mole fraction of the solute.
00329      *    - \f$ X_{o,p} = max (X_{o}^{min}, X_o) \f$
00330      *    - \f$ X_{o}^{min} \f$ = minimum mole fraction of solvent allowed
00331      *              in the denominator.
00332      */
00333     void calcMolalities() const;
00334 
00335     //!  This function will return the molalities of the species.
00336     /*!
00337      *   We calculate the vector of molalities of the species
00338      *   in the phase
00339      * \f[
00340      *     m_i = \frac{X_i}{1000 * M_o * X_{o,p}}
00341      * \f]
00342      *    where 
00343      *    - \f$ M_o \f$ is the molecular weight of the solvent
00344      *    - \f$ X_o \f$ is the mole fraction of the solvent
00345      *    - \f$ X_i \f$ is the mole fraction of the solute.
00346      *    - \f$ X_{o,p} = \max (X_{o}^{min}, X_o) \f$
00347      *    - \f$ X_{o}^{min} \f$ = minimum mole fraction of solvent allowed
00348      *              in the denominator.
00349      *
00350      * @param molal       Output vector of molalities. Length: m_kk.
00351      */
00352     void getMolalities(doublereal * const molal) const;
00353 
00354     //! Set the molalities of the solutes in a phase
00355     /*!
00356      * Note, the entry for the solvent is not used.
00357      *   We are supplied with the molalities of all of the
00358      *   solute species. We then calculate the mole fractions of all
00359      *   species and update the %ThermoPhase object.
00360      *  \f[
00361      *     m_i = \frac{X_i}{M_o/1000 * X_{o,p}}
00362      *  \f]
00363      *    where 
00364      *    -  \f$M_o\f$ is the molecular weight of the solvent
00365      *    -  \f$X_o\f$ is the mole fraction of the solvent
00366      *    -  \f$X_i\f$ is the mole fraction of the solute.
00367      *    -  \f$X_{o,p} = \max(X_o^{min}, X_o)\f$
00368      *    -  \f$X_o^{min}\f$ = minimum mole fraction of solvent allowed
00369      *                     in the denominator.
00370      *
00371      * The formulas for calculating mole fractions are
00372      *  \f[
00373      *       L^{sum} = \frac{1}{\tilde{M}_o X_o} = \frac{1}{\tilde{M}_o} + \sum_{i\ne o} m_i
00374      *  \f]
00375      *  Then,
00376      *  \f[
00377      *       X_o =   \frac{1}{\tilde{M}_o L^{sum}}  
00378      *  \f]
00379      *  \f[
00380      *       X_i =   \frac{m_i}{L^{sum}}  
00381      *  \f]
00382      *  It is currently an error if the solvent mole fraction is attempted to be set 
00383      *  to a value lower than  \f$X_o^{min}\f$.
00384      *
00385      * @param molal   Input vector of molalities. Length: m_kk.
00386      */
00387     void setMolalities(const doublereal * const molal);
00388 
00389     //! Set the molalities of a phase
00390     /*!
00391      * Set the molalities of the solutes in a phase. Note, the entry for the
00392      * solvent is not used.
00393      *
00394      * @param xMap  Composition Map containing the molalities.
00395      */
00396     void setMolalitiesByName(compositionMap& xMap);
00397 
00398     //! Set the molalities of a phase
00399     /*!
00400      * Set the molalities of the solutes in a phase. Note, the entry for the
00401      * solvent is not used.
00402      *
00403      * @param name  String containing the information for a composition map.
00404      */
00405     void setMolalitiesByName(const std::string& name);
00406 
00407     /**
00408      * @}
00409      * @name Mechanical Properties
00410      * @{
00411      */
00412 
00413     /**
00414      * @} 
00415      * @name Potential Energy
00416      * 
00417      * Species may have an additional potential energy due to the
00418      * presence of external gravitation or electric fields. These
00419      * methods allow specifying a potential energy for individual
00420      * species.
00421      * @{
00422      */
00423 
00424     /**
00425      * @}
00426      * @name Activities, Standard States, and Activity Concentrations
00427      *
00428      * The activity \f$a_k\f$ of a species in solution is
00429      * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
00430      * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
00431      * the chemical potential at unit activity, which depends only
00432      * on temperature and pressure.
00433      * @{
00434      */
00435 
00436     /**
00437      * This method returns the activity convention.
00438      * Currently, there are two activity conventions
00439      *  Molar-based activities
00440      *       Unit activity of species at either a hypothetical pure
00441      *       solution of the species or at a hypothetical
00442      *       pure ideal solution at infinite dilution
00443      *   cAC_CONVENTION_MOLAR 0
00444      *      - default
00445      *  
00446      *  Molality based acvtivities
00447      *       (unit activity of solutes at a hypothetical 1 molal
00448      *        solution referenced to infinite dilution at all
00449      *        pressures and temperatures).
00450      *   cAC_CONVENTION_MOLALITY 1
00451      *
00452      *  We set the convention to molality here.
00453      */
00454     int activityConvention() const;
00455 
00456     /**
00457      * This method returns an array of generalized concentrations
00458      * \f$ C_k\f$ that are defined such that 
00459      * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ 
00460      * is a standard concentration
00461      * defined below.  These generalized concentrations are used
00462      * by kinetics manager classes to compute the forward and
00463      * reverse rates of elementary reactions. 
00464      *
00465      * @param c Array of generalized concentrations. The 
00466      *           units depend upon the implementation of the
00467      *           reaction rate expressions within the phase.
00468      */
00469     virtual void getActivityConcentrations(doublereal* c) const;
00470 
00471     /**
00472      * The standard concentration \f$ C^0_k \f$ used to normalize
00473      * the generalized concentration. In many cases, this quantity
00474      * will be the same for all species in a phase - for example,
00475      * for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this
00476      * reason, this method returns a single value, instead of an
00477      * array.  However, for phases in which the standard
00478      * concentration is species-specific (e.g. surface species of
00479      * different sizes), this method may be called with an
00480      * optional parameter indicating the species.
00481      *
00482      * @param k species index. Defaults to zero.
00483      */
00484     virtual doublereal standardConcentration(int k=0) const;
00485 
00486     /**
00487      * Returns the natural logarithm of the standard 
00488      * concentration of the kth species
00489      *
00490      * @param k  species index
00491      */
00492     virtual doublereal logStandardConc(int k=0) const;
00493 
00494     /**
00495      * Returns the units of the standard and generalized
00496      * concentrations Note they have the same units, as their
00497      * ratio is defined to be equal to the activity of the kth
00498      * species in the solution, which is unitless.
00499      *
00500      * This routine is used in print out applications where the
00501      * units are needed. Usually, MKS units are assumed throughout
00502      * the program and in the XML input files.
00503      *
00504      * @param uA Output vector containing the units
00505      *  uA[0] = kmol units - default  = 1
00506      *  uA[1] = m    units - default  = -nDim(), the number of spatial
00507      *                                dimensions in the Phase class.
00508      *  uA[2] = kg   units - default  = 0;
00509      *  uA[3] = Pa(pressure) units - default = 0;
00510      *  uA[4] = Temperature units - default = 0;
00511      *  uA[5] = time units - default = 0
00512      * @param k species index. Defaults to 0.
00513      * @param sizeUA output int containing the size of the vector.
00514      *        Currently, this is equal to 6.
00515      */
00516     virtual void getUnitsStandardConc(double *uA, int k = 0,
00517                                       int sizeUA = 6) const;
00518 
00519     
00520     //! Get the array of non-dimensional activities (molality
00521     //! based for this class and classes that derive from it) at
00522     //! the current solution temperature, pressure, and solution concentration.
00523     /*!
00524      * All standard state properties for molality-based phases are
00525      * evaluated consistent with the molality scale. Therefore, this function
00526      * must return molality-based activities.
00527      *
00528      * \f[
00529      *  a_i^\triangle = \gamma_k^{\triangle} \frac{m_k}{m^\triangle}
00530      * \f]
00531      *
00532      * This function must be implemented in derived classes.
00533      *
00534      * @param ac     Output vector of molality-based activities. Length: m_kk.
00535      */
00536     virtual void getActivities(doublereal* ac) const;
00537 
00538     //! Get the array of non-dimensional activity coefficients at
00539     //! the current solution temperature, pressure, and solution concentration.
00540     /*!
00541      * These are mole-fraction based activity coefficients. In this
00542      * object, their calculation is based on translating the values
00543      * of the molality-based activity coefficients.
00544      *  See Denbigh p. 278 for a thorough discussion.
00545      *
00546      * The molar-based activity coefficients \f$ \gamma_k \f$ may be calculated from the 
00547      * molality-based
00548      * activity coefficients, \f$ \gamma_k^\triangle \f$ by the following 
00549      * formula.
00550      * \f[
00551      *     \gamma_k = \frac{\gamma_k^\triangle}{X_o}
00552      * \f]
00553      *
00554      * For purposes of establishing a convention, the molar activity coefficient of the
00555      * solvent is set equal to the molality-based activity coefficient of the
00556      * solvent:
00557      * 
00558      * \f[
00559      *     \gamma_o = \gamma_o^\triangle
00560      * \f]
00561      *
00562      * Derived classes don't need to overload this function. This function is
00563      * handled at this level.
00564      *
00565      * @param ac  Output vector containing the mole-fraction based activity coefficients.
00566      *            length: m_kk.
00567      */
00568     void getActivityCoefficients(doublereal* ac) const;
00569 
00570     //!  Get the array of non-dimensional molality based 
00571     //!  activity coefficients at the current solution temperature, 
00572     //!  pressure, and  solution concentration.
00573     /*!
00574      *  See Denbigh p. 278 for a thorough discussion. This class must be overwritten in
00575      *  classes which derive from %MolalityVPSSTP. This function takes over from the
00576      *  molar-based activity coefficient calculation, getActivityCoefficients(), in
00577      *  derived classes.
00578      *
00579      *  These molality based activity coefficients are scaled according to the current
00580      *  pH scale. See the Eq3/6 manual for details.
00581      *
00582      *  Activity coefficients for species k may be altered between scales s1 to s2
00583      *  using the following formula
00584      *
00585      *   \f[
00586      *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) 
00587      *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
00588      *   \f]
00589      *
00590      *  where j is any one species. For the NBS scale, j is equal to the Cl- species
00591      *  and 
00592      *
00593      *  \f[
00594      *       ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}}
00595      *  \f]
00596      *
00597      * @param acMolality Output vector containing the molality based activity coefficients.
00598      *                   length: m_kk.
00599      */
00600     virtual void getMolalityActivityCoefficients(doublereal *acMolality) const;
00601 
00602   
00603    
00604     //! Calculate the osmotic coefficient
00605     /*!
00606      *   \f[
00607      *       \phi = \frac{- ln(a_o)}{\tilde{M}_o \sum_{i \ne o} m_i}
00608      *   \f]
00609      * 
00610      *  Note there  are a few of definitions of the osmotic coefficient floating
00611      *  around. We use the one defined in 
00612      *  (Activity Coefficients in Electrolyte Solutions, K. S. Pitzer
00613      *   CRC Press, Boca Raton, 1991, p. 85, Eqn. 28). This definition is most clearly
00614      *  related to theoretical calculation.
00615      *
00616      *  units = dimensionless
00617      */
00618     virtual double osmoticCoefficient() const;
00619 
00620     //@}
00621     /// @name  Partial Molar Properties of the Solution 
00622     //@{
00623 
00624 
00625     /**
00626      * Get the species electrochemical potentials. 
00627      * These are partial molar quantities.
00628      * This method adds a term \f$ Fz_k \phi_k \f$ to the 
00629      * to each chemical potential.
00630      *
00631      * Units: J/kmol
00632      *
00633      * @param mu     output vector containing the species electrochemical potentials.
00634      *               Length: m_kk.
00635      */
00636     void getElectrochemPotentials(doublereal* mu) const;
00637 
00638  
00639     //@}
00640     /// @name  Properties of the Standard State of the Species in the Solution
00641     //@{
00642 
00643      
00644 
00645     //@}
00646     /// @name Thermodynamic Values for the Species Reference States
00647     //@{
00648 
00649 
00650     ///////////////////////////////////////////////////////
00651     //
00652     //  The methods below are not virtual, and should not
00653     //  be overloaded.
00654     //
00655     //////////////////////////////////////////////////////
00656 
00657     /**
00658      * @name Specific Properties
00659      * @{
00660      */
00661 
00662 
00663     /**
00664      * @name Setting the State
00665      *
00666      * These methods set all or part of the thermodynamic
00667      * state.
00668      * @{
00669      */
00670 
00671     //@}
00672 
00673     /**
00674      * @name Chemical Equilibrium
00675      * Routines that implement the Chemical equilibrium capability
00676      * for a single phase, based on the element-potential method.
00677      * @{
00678      */
00679 
00680     /**
00681      * This method is used by the ChemEquil element-potential
00682      * based equilibrium solver.
00683      * It sets the state such that the chemical potentials of the 
00684      * species within the current phase satisfy
00685      * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m}
00686      * \left(\frac{\lambda_m} {\hat R T}\right) \f] where 
00687      * \f$ \lambda_m \f$ is the element potential of element m. The
00688      * temperature is unchanged.  Any phase (ideal or not) that
00689      * implements this method can be equilibrated by ChemEquil.
00690      *
00691      * @param lambda_RT Input vector containing the dimensionless 
00692      *                  element potentials.
00693      */
00694     virtual void setToEquilState(const doublereal* lambda_RT);
00695 
00696 
00697     //@}
00698 
00699 
00700     //! Set equation of state parameter values from XML entries.
00701     /*!
00702      * This method is called by function importPhase() in
00703      * file importCTML.cpp when processing a phase definition in
00704      * an input file. It should be overloaded in subclasses to set
00705      * any parameters that are specific to that particular phase
00706      * model. 
00707      *
00708      * The MolalityVPSSTP object defines a new method for setting
00709      * the concentrations of a phase. The new method is defined by a
00710      * block called "soluteMolalities". If this block
00711      * is found, the concentrations within that phase are
00712      * set to the "name":"molalities pairs found within that
00713      * XML block. The solvent concentration is then set
00714      * to everything else.
00715      *
00716      * The function first calls the overloaded function ,
00717      * VPStandardStateTP::setStateFromXML(), to pick up the parent class
00718      * behavior.
00719      *
00720      * usage: Overloaded functions should call this function
00721      *        before carrying out their own behavior.
00722      *   
00723      * @param state An XML_Node object corresponding to
00724      *              the "state" entry for this phase in the input file.
00725      */
00726     virtual void setStateFromXML(const XML_Node& state);
00727 
00728     /// The following methods are used in the process of constructing
00729     /// the phase and setting its parameters from a specification in an 
00730     /// input file. They are not normally used in application programs.
00731     /// To see how they are used, see files importCTML.cpp and 
00732     /// ThermoFactory.cpp.
00733 
00734 
00735     /*!
00736      * @internal Initialize. This method is provided to allow
00737      * subclasses to perform any initialization required after all
00738      * species have been added. For example, it might be used to
00739      * resize internal work arrays that must have an entry for
00740      * each species.  The base class implementation does nothing,
00741      * and subclasses that do not require initialization do not
00742      * need to overload this method.  When importing a CTML phase
00743      * description, this method is called just prior to returning
00744      * from function importPhase.
00745      *
00746      * @see importCTML.cpp
00747      */
00748     virtual void initThermo();
00749 
00750 
00751     /**
00752      *   Import and initialize a ThermoPhase object
00753      *
00754      * @param phaseNode This object must be the phase node of a
00755      *             complete XML tree
00756      *             description of the phase, including all of the
00757      *             species data. In other words while "phase" must
00758      *             point to an XML phase object, it must have
00759      *             sibling nodes "speciesData" that describe
00760      *             the species in the phase.
00761      * @param id   ID of the phase. If nonnull, a check is done
00762      *             to see if phaseNode is pointing to the phase
00763      *             with the correct id. 
00764      */
00765     void initThermoXML(XML_Node& phaseNode, std::string id);
00766 
00767     
00768     //! Set the temperature (K), pressure (Pa), and molalities
00769     //!(gmol kg-1) of the solutes
00770     /*!
00771      * @param t           Temperature (K)
00772      * @param p           Pressure (Pa)
00773      * @param molalities  Input vector of molalities of the solutes.
00774      *                    Length: m_kk.
00775      */
00776     void setState_TPM(doublereal t, doublereal p, 
00777                       const doublereal * const molalities);
00778 
00779     //! Set the temperature (K), pressure (Pa), and molalities.
00780     /*!
00781      * @param t           Temperature (K)
00782      * @param p           Pressure (Pa)
00783      * @param m           compositionMap containing the molalities
00784      */
00785     void setState_TPM(doublereal t, doublereal p, compositionMap& m);
00786 
00787     //! Set the temperature (K), pressure (Pa), and molalities. 
00788     /*!
00789      * @param t           Temperature (K)
00790      * @param p           Pressure (Pa)
00791      * @param m           String which gets translated into a composition
00792      *                    map for the molalities of the solutes.
00793      */
00794     void setState_TPM(doublereal t, doublereal p, const std::string& m);
00795 
00796     //! returns a summary of the state of the phase as a string
00797     /*!
00798      * @param show_thermo If true, extra information is printed out
00799      *                    about the thermodynamic state of the system.
00800      */
00801     virtual std::string report(bool show_thermo = true) const;
00802 
00803   protected:
00804 
00805     //! Get the array of unscaled non-dimensional molality based 
00806     //!  activity coefficients at the current solution temperature, 
00807     //!  pressure, and  solution concentration.
00808     /*!
00809      *  See Denbigh p. 278 for a thorough discussion. This class must be overwritten in
00810      *  classes which derive from %MolalityVPSSTP. This function takes over from the
00811      *  molar-based activity coefficient calculation, getActivityCoefficients(), in
00812      *  derived classes.
00813      *
00814      * @param acMolality Output vector containing the molality based activity coefficients.
00815      *                   length: m_kk.
00816      */
00817     virtual void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const;
00818 
00819     //! Apply the current phScale to a set of activity Coefficients or activities
00820     /*!
00821      *  See the Eq3/6 Manual for a thorough discussion.
00822      *
00823      * @param acMolality input/Output vector containing the molality based 
00824      *                   activity coefficients. length: m_kk.
00825      */
00826     virtual void applyphScale(doublereal *acMolality) const;
00827 
00828   private:
00829     //! Returns the index of the Cl- species.
00830     /*!
00831      *  The Cl- species is special in the sense that it's single ion
00832      *  molalality-based activity coefficient is used in the specification
00833      *  of the pH scale for single ions. Therefore, we need to know
00834      *  what species index is Cl-. If the species isn't in the species
00835      *  list then this routine returns -1, and we can't use the NBS
00836      *  pH scale. 
00837      *    
00838      *  Right now we use a restrictive interpretation. The species
00839      *  must be named "Cl-". It must consist of exactly one Cl and one E 
00840      *  atom. 
00841      */
00842     virtual int findCLMIndex() const;
00843 
00844     //! Initialize lengths of local variables after all species have
00845     //! been identified.
00846     void initLengths();
00847             
00848   protected:
00849 
00850     //! Index of the solvent
00851     /*!
00852      * Currently the index of the solvent is hard-coded to the value 0
00853      */
00854     int        m_indexSolvent;
00855 
00856     //! Scaling to be used for output of single-ion species activity
00857     //! coefficients.
00858     /*!
00859      *   Index of the species to be used in the single-ion scaling
00860      *   law. This is the indentity of the Cl- species for the PHSCALE_NBS
00861      *   scaling.
00862      *   Either PHSCALE_PITZER or PHSCALE_NBS
00863      */
00864     int        m_pHScalingType;
00865 
00866     //! Index of the phScale species
00867     /*!
00868      *   Index of the species to be used in the single-ion scaling
00869      *   law. This is the indentity of the Cl- species for the PHSCALE_NBS
00870      *   scaling
00871      */
00872     int        m_indexCLM;
00873 
00874     //! Molecular weight of the Solvent
00875     doublereal m_weightSolvent;
00876 
00877     /*!
00878      * In any molality implementation, it makes sense to have
00879      * a minimum solvent mole fraction requirement, since the 
00880      * implementation becomes singular in the xmolSolvent=0
00881      * limit. The default is to set it to 0.01.
00882      * We then modify the molality definition to ensure that
00883      * molal_solvent = 0 when xmol_solvent = 0. 
00884      */
00885     doublereal m_xmolSolventMIN;
00886 
00887     //! This is the multiplication factor that goes inside 
00888     //! log expressions involving the molalities of species.
00889     /*!
00890      * Its equal to Wt_0 / 1000,
00891      *     where Wt_0 = weight of solvent (kg/kmol)
00892      */
00893     doublereal m_Mnaught;
00894 
00895     //! Current value of the molalities of the species in the phase.
00896     /*!
00897      * Note this vector is a mutable quantity.
00898      * units are (kg/kmol)
00899      */
00900     mutable vector_fp  m_molalities;
00901 
00902   private:
00903     //! Error function
00904     /*!
00905      *  Print an error string and exit
00906      *
00907      * @param msg  Message to be printed
00908      */
00909     doublereal err(std::string msg) const;
00910 
00911   };
00912 
00913 
00914   //! Scale to be used for the output of single-ion activity coefficients
00915   //! is that used by Pitzer.
00916   /*!
00917    *  This is the internal scale used within the code. One property is that
00918    *  the activity coefficients for the cation and anion of a single salt
00919    *  will be equal. This scale is the one presumed by the formulation of the
00920    *  single-ion activity coefficients described in this report.
00921    *
00922    *  Activity coefficients for species k may be altered between scales s1 to s2
00923    *  using the following formula
00924    *
00925    *   \f[
00926    *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) 
00927    *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
00928    *   \f]
00929    *
00930    *  where j is any one species.
00931    *
00932    *
00933    */
00934   const int PHSCALE_PITZER = 0;
00935 
00936   //! Scale to be used for evaluation of single-ion activity coefficients
00937   //! is that used by the NBS standard for evaluation of the pH variable.
00938   /*!
00939    *  This is not the internal scale used within the code.
00940    *
00941    *  Activity coefficients for species k may be altered between scales s1 to s2
00942    *  using the following formula
00943    *
00944    *   \f[
00945    *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) 
00946    *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
00947    *   \f]
00948    *
00949    *  where j is any one species. For the NBS scale, j is equal to the Cl- species
00950    *  and 
00951    *
00952    *  \f[
00953    *       ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}}
00954    *  \f]
00955    *
00956    *  This is the NBS pH scale, which is used in all conventional pH 
00957    *  measurements. and is based on the Bates-Guggenheim quations.
00958    *
00959    */
00960   const int PHSCALE_NBS    = 1;
00961 
00962 }
00963         
00964 #endif
00965 
00966 
00967 
00968 
00969 
Generated by  doxygen 1.6.3