From 6ba38684de768987793dfd8d1fdbb1ab2fb6b858 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 18 Dec 2008 15:46:27 +0000 Subject: [PATCH] Added pH scaling to HMWSoln, and checked it against EQ3 result. This is preliminary. --- Cantera/src/thermo/HMWSoln.cpp | 44 +++++++++++++- Cantera/src/thermo/HMWSoln.h | 48 +++++++++------ Cantera/src/thermo/MolalityVPSSTP.cpp | 85 +++++++++++++++++++++++++-- Cantera/src/thermo/MolalityVPSSTP.h | 79 ++++++++++++++++++++++++- 4 files changed, 227 insertions(+), 29 deletions(-) diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp index 75f611cba..724f8a9d6 100644 --- a/Cantera/src/thermo/HMWSoln.cpp +++ b/Cantera/src/thermo/HMWSoln.cpp @@ -917,10 +917,14 @@ namespace Cantera { double xmolSolvent = moleFraction(m_indexSolvent); ac[m_indexSolvent] = exp(m_lnActCoeffMolal[m_indexSolvent]) * xmolSolvent; + /* + * Apply the pH scale + */ + applyphScale(ac); } /* - * getMolalityActivityCoefficients() (virtual, const) + * getUnscaledMolalityActivityCoefficients() (virtual, const) * * Get the array of non-dimensional Molality based * activity coefficients at @@ -931,7 +935,7 @@ namespace Cantera { * Note, most of the work is done in an internal private routine */ void HMWSoln:: - getMolalityActivityCoefficients(doublereal* acMolality) const { + getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const { updateStandardStateThermo(); A_Debye_TP(-1.0, -1.0); s_update_lnMolalityActCoeff(); @@ -5521,7 +5525,7 @@ namespace Cantera { } - /***********************************************************************************************/ + /**********************************************************************************************/ /* * Calculate the lambda interactions. @@ -5827,6 +5831,40 @@ namespace Cantera { } } + + // Apply the current phScale to a set of activity Coefficients or activities + /* + * See the Eq3/6 Manual for a thorough discussion. + * + * @param acMolality input/Output vector containing the molality based + * activity coefficients. length: m_kk. + */ + void HMWSoln::applyphScale(doublereal *acMolality) const { + if (m_pHScalingType == PHSCALE_PITZER) return; + if (m_pHScalingType != PHSCALE_NBS) { + throw CanteraError("", "shoudln't be here"); + } + + /* + * Find the ionic strength + */ + doublereal Is = m_IionicMolality; + doublereal sqrtIs = sqrt(Is); + + /* + * Find the Debye Huckel coefficient + */ + doublereal A = m_A_Debye; + doublereal lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs); + doublereal lnGammaCLMs1 = m_lnActCoeffMolal[m_indexCLM]; + doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1); + + for (int k = 1; k < m_kk; k++) { + acMolality[k] *= exp(m_speciesCharge[k] * afac); + } + } + + int HMWSoln::debugPrinting() { #ifdef DEBUG_MODE return m_debugCalc; diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index ff0c918ec..e138d5ec7 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -1647,19 +1647,6 @@ namespace Cantera { */ virtual void getActivities(doublereal* ac) const; - //! Get the array of non-dimensional molality-based - //! activity coefficients at - //! the current solution temperature, pressure, and solution concentration. - /*! - * note solvent is on molar scale. The solvent molar - * based activity coefficient is returned. - * - * @param acMolality Vector of Molality-based activity coefficients - * Length: m_kk - */ - virtual void - getMolalityActivityCoefficients(doublereal* acMolality) const; - //@} /// @name Partial Molar Properties of the Solution ----------------- //@{ @@ -2035,8 +2022,8 @@ namespace Cantera { */ virtual void initThermoXML(XML_Node& phaseNode, std::string id); - /** - * Report the molar volume of species k + //! Report the molar volume of species k + /*! * * units - \f$ m^3 kmol^-1 \f$ * @@ -2064,10 +2051,11 @@ namespace Cantera { virtual double A_Debye_TP(double temperature = -1.0, double pressure = -1.0) const; - /** - * Value of the derivative of the Debye Huckel constant with - * respect to temperature as a function of temperature - * and pressure. + + //! Value of the derivative of the Debye Huckel constant with + //! respect to temperature as a function of temperature + //! and pressure. + /*! * * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) * @@ -2197,6 +2185,28 @@ namespace Cantera { */ void printCoeffs () const; + //! Get the array of unscaled non-dimensional molality based + //! activity coefficients at the current solution temperature, + //! pressure, and solution concentration. + /*! + * See Denbigh p. 278 for a thorough discussion. This class must be overwritten in + * classes which derive from %MolalityVPSSTP. This function takes over from the + * molar-based activity coefficient calculation, getActivityCoefficients(), in + * derived classes. + * + * @param acMolality Output vector containing the molality based activity coefficients. + * length: m_kk. + */ + void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const; + + //! Apply the current phScale to a set of activity Coefficients or activities + /*! + * See the Eq3/6 Manual for a thorough discussion. + * + * @param acMolality input/Output vector containing the molality based + * activity coefficients. length: m_kk. + */ + void applyphScale(doublereal *acMolality) const; //@} diff --git a/Cantera/src/thermo/MolalityVPSSTP.cpp b/Cantera/src/thermo/MolalityVPSSTP.cpp index 25ea98838..86d3c6d37 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.cpp +++ b/Cantera/src/thermo/MolalityVPSSTP.cpp @@ -118,6 +118,42 @@ namespace Cantera { * -------------- Utilities ------------------------------- */ + // Equation of state type flag. + /* + * The ThermoPhase base class returns + * zero. Subclasses should define this to return a unique + * non-zero value. Known constants defined for this purpose are + * listed in mix_defs.h. The MolalityVPSSTP class also returns + * zero, as it is a non-complete class. + */ + int MolalityVPSSTP::eosType() const { + return 0; + } + + // Set the pH scale, which determines the scale for single-ion activity + // coefficients. + /* + * Single ion activity coefficients are not unique in terms of the + * representing actual measureable quantities. + */ + void MolalityVPSSTP::setpHScale(const int pHscaleType) { + m_pHScalingType = pHscaleType; + if (pHscaleType != PHSCALE_PITZER && pHscaleType != PHSCALE_NBS) { + throw CanteraError("MolalityVPSSTP::setpHScale", + "Unknown scale type: " + int2str(pHscaleType)); + } + } + + // Reports the pH scale, which determines the scale for single-ion activity + // coefficients. + /* + * Single ion activity coefficients are not unique in terms of the + * representing actual measureable quantities. + */ + int MolalityVPSSTP::pHScale() const { + return m_pHScalingType; + } + /* * setSolvent(): * Utilities for Solvent ID and Molality @@ -441,11 +477,25 @@ namespace Cantera { } } - + // Get the array of non-dimensional molality based + // activity coefficients at the current solution temperature, + // pressure, and solution concentration. + /* + * See Denbigh p. 278 for a thorough discussion. This class must be overwritten in + * classes which derive from %MolalityVPSSTP. This function takes over from the + * molar-based activity coefficient calculation, getActivityCoefficients(), in + * derived classes. + * + * Note these activity coefficients have the current pH scale applied to them. + * + * @param acMolality Output vector containing the molality based activity coefficients. + * length: m_kk. + */ void MolalityVPSSTP::getMolalityActivityCoefficients(doublereal *acMolality) const { - err("getMolalityActivityCoefficients"); + getUnscaledMolalityActivityCoefficients(acMolality); + applyphScale(acMolality); } - + /* * osmotic coefficient: * @@ -608,7 +658,34 @@ namespace Cantera { m_indexCLM = findCLMIndex(); } - // Returns the index of the Cl- species. + // Get the array of unscaled non-dimensional molality based + // activity coefficients at the current solution temperature, + // pressure, and solution concentration. + /* + * See Denbigh p. 278 for a thorough discussion. This class must be overwritten in + * classes which derive from %MolalityVPSSTP. This function takes over from the + * molar-based activity coefficient calculation, getActivityCoefficients(), in + * derived classes. + * + * @param acMolality Output vector containing the molality based activity coefficients. + * length: m_kk. + */ + void MolalityVPSSTP::getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const { + err("getUnscaledMolalityActivityCoefficients"); + } + + // Apply the current phScale to a set of activity Coefficients or activities + /* + * See the Eq3/6 Manual for a thorough discussion. + * + * @param acMolality input/Output vector containing the molality based + * activity coefficients. length: m_kk. + */ + void MolalityVPSSTP::applyphScale(doublereal *acMolality) const { + err("applyphScale"); + } + + // Returns the index of the Cl- species. /* * The Cl- species is special in the sense that it's single ion * molalality-based activity coefficient is used in the specification diff --git a/Cantera/src/thermo/MolalityVPSSTP.h b/Cantera/src/thermo/MolalityVPSSTP.h index d36e74bb9..11fbc9206 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.h +++ b/Cantera/src/thermo/MolalityVPSSTP.h @@ -154,6 +154,27 @@ namespace Cantera { * * All objects that derive from this are assumed to have molality based standard states. * + * Molality based activity coefficients are scaled according to the current + * pH scale. See the Eq3/6 manual for details. + * + * Activity coefficients for species k may be altered between scales s1 to s2 + * using the following formula + * + * \f[ + * ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) + * + \frac{z_k}{z_j} \left( ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right) + * \f] + * + * where j is any one species. For the NBS scale, j is equal to the Cl- species + * and + * + * \f[ + * ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}} + * \f] + * + * The Pitzer scale doesn't actually change anything. The pitzer scale is defined + * as the raw unscaled activity coefficients produced by the underlying objects. + * * @todo Make two solvent minimum fractions. One would be for calculation of the non-ideal * factors. The other one would be for purposes of stoichiometry evaluation. the * stoichiometry evaluation one would be a 1E-13 limit. Anything less would create @@ -219,7 +240,7 @@ namespace Cantera { * listed in mix_defs.h. The MolalityVPSSTP class also returns * zero, as it is a non-complete class. */ - virtual int eosType() const { return 0; } + virtual int eosType() const; //! Set the pH scale, which determines the scale for single-ion activity @@ -230,6 +251,14 @@ namespace Cantera { */ void setpHScale(const int pHscaleType); + //! Reports the pH scale, which determines the scale for single-ion activity + //! coefficients. + /*! + * Single ion activity coefficients are not unique in terms of the + * representing actual measureable quantities. + */ + int pHScale() const; + /** * @} * @name Molar Thermodynamic Properties @@ -526,7 +555,7 @@ namespace Cantera { */ void getActivityCoefficients(doublereal* ac) const; - //! Get the array of non-dimensional molality based + //! Get the array of non-dimensional molality based //! activity coefficients at the current solution temperature, //! pressure, and solution concentration. /*! @@ -535,10 +564,30 @@ namespace Cantera { * molar-based activity coefficient calculation, getActivityCoefficients(), in * derived classes. * + * These molality based activity coefficients are scaled according to the current + * pH scale. See the Eq3/6 manual for details. + * + * Activity coefficients for species k may be altered between scales s1 to s2 + * using the following formula + * + * \f[ + * ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) + * + \frac{z_k}{z_j} \left( ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right) + * \f] + * + * where j is any one species. For the NBS scale, j is equal to the Cl- species + * and + * + * \f[ + * ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}} + * \f] + * * @param acMolality Output vector containing the molality based activity coefficients. * length: m_kk. */ virtual void getMolalityActivityCoefficients(doublereal *acMolality) const; + + //! Calculate the osmotic coefficient /*! @@ -739,6 +788,31 @@ namespace Cantera { */ virtual std::string report(bool show_thermo = true) const; + protected: + + //! Get the array of unscaled non-dimensional molality based + //! activity coefficients at the current solution temperature, + //! pressure, and solution concentration. + /*! + * See Denbigh p. 278 for a thorough discussion. This class must be overwritten in + * classes which derive from %MolalityVPSSTP. This function takes over from the + * molar-based activity coefficient calculation, getActivityCoefficients(), in + * derived classes. + * + * @param acMolality Output vector containing the molality based activity coefficients. + * length: m_kk. + */ + virtual void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const; + + //! Apply the current phScale to a set of activity Coefficients or activities + /*! + * See the Eq3/6 Manual for a thorough discussion. + * + * @param acMolality input/Output vector containing the molality based + * activity coefficients. length: m_kk. + */ + virtual void applyphScale(doublereal *acMolality) const; + private: //! Returns the index of the Cl- species. /*! @@ -825,7 +899,6 @@ namespace Cantera { }; - //! Scale to be used for the output of single-ion activity coefficients //! is that used by Pitzer. /*!