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