Added pH scaling to HMWSoln, and checked it against EQ3 result. This

is preliminary.
This commit is contained in:
Harry Moffat 2008-12-18 15:46:27 +00:00
parent 0687ef6eda
commit 6ba38684de
4 changed files with 227 additions and 29 deletions

View file

@ -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;

View file

@ -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;
//@}

View file

@ -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

View file

@ -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.
/*!