diff --git a/include/cantera/thermo/MaskellSolidSolnPhase.h b/include/cantera/thermo/MaskellSolidSolnPhase.h
index 02478ba36..405012748 100644
--- a/include/cantera/thermo/MaskellSolidSolnPhase.h
+++ b/include/cantera/thermo/MaskellSolidSolnPhase.h
@@ -34,26 +34,8 @@ namespace Cantera
class MaskellSolidSolnPhase : public VPStandardStateTP
{
public:
- /**
- * The generalized concentrations can have three different forms
- * depending on the value of the member attribute #m_formGC, which
- * is supplied in the constructor or read from the xml data file.
- *
- * @param formCG This parameter initializes the #m_formGC variable.
- */
MaskellSolidSolnPhase();
- //! Construct and initialize an MaskellSolidSolnPhase ThermoPhase object
- //! directly from an XML database
- /*!
- * @param root XML tree containing a description of the phase.
- * The tree must be positioned at the XML element
- * named phase with id, "id", on input to this routine.
- * @param id The name of this phase. This is used to look up
- * the phase in the XML datafile.
- */
- //MaskellSolidSolnPhase(XML_Node& root, const std::string& id="");
-
//! Copy Constructor
MaskellSolidSolnPhase(const MaskellSolidSolnPhase&);
@@ -67,6 +49,61 @@ public:
*/
virtual ThermoPhase* duplMyselfAsThermoPhase() const;
+ /**
+ * This method returns the array of generalized
+ * concentrations. The generalized concentrations are used
+ * in the evaluation of the rates of progress for reactions
+ * involving species in this phase. The generalized
+ * concentration divided by the standard concentration is also
+ * equal to the activity of species.
+ *
+ * For this implementation the activity is defined to be the
+ * mole fraction of the species. The generalized concentration
+ * is defined to be equal to the mole fraction divided by
+ * the partial molar volume. The generalized concentrations
+ * for species in this phase therefore have units of
+ * kmol m-3. Rate constants must reflect this fact.
+ *
+ * On a general note, the following must be true.
+ * For an ideal solution, the generalized concentration must consist
+ * of the mole fraction multiplied by a constant. The constant may be
+ * fairly arbitrarily chosen, with differences adsorbed into the
+ * reaction rate expression. 1/V_N, 1/V_k, or 1 are equally good,
+ * as long as the standard concentration is adjusted accordingly.
+ * However, it must be a constant (and not the concentration, btw,
+ * which is a function of the mole fractions) in order for the
+ * ideal solution properties to hold at the same time having the
+ * standard concentration to be independent of the mole fractions.
+ *
+ * In this implementation the form of the generalized concentrations
+ * depend upon the member attribute, #m_formGC.
+ *
+ * HKM Note: We have absorbed the pressure dependence of the pure species
+ * state into the thermodynamics functions. Therefore the
+ * standard state on which the activities are based depend
+ * on both temperature and pressure. If we hadn't, it would have
+ * appeared in this function in a very awkward exp[] format.
+ *
+ * @param c Pointer to array of doubles of length m_kk, which on exit
+ * will contain the generalized concentrations.
+ */
+ virtual void getActivityConcentrations(doublereal* c) const;
+
+ //! Return the standard concentration for the kth species
+ /*!
+ * The standard concentration \f$ C^0_k \f$ used to normalize the
+ * generalized concentration. In many cases, this quantity will be the
+ * same for all species in a phase. However, for this case, we will return
+ * a distinct concentration for each species. This is the inverse of the
+ * species molar volume. Units for the standard concentration are kmol
+ * m-3.
+ *
+ * @param k Species number: this is a require parameter,
+ * a change from the ThermoPhase base class, where it was
+ * an optional parameter.
+ */
+ virtual doublereal standardConcentration(size_t k) const;
+
//! @name Molar Thermodynamic Properties of the Solution
//! @{
/**
@@ -310,6 +347,9 @@ protected:
//! Value of the enthalpy change on mixing due to protons changing from type B to type A configurations.
doublereal h_mixing;
+ //! Index of the species whose mole fraction defines the extent of reduction r
+ int product_species_index;
+
private:
// Functions to calculate some of the pieces of the mixing terms.
doublereal s() const;
diff --git a/src/thermo/MaskellSolidSolnPhase.cpp b/src/thermo/MaskellSolidSolnPhase.cpp
index 9430d31bb..f72101d69 100644
--- a/src/thermo/MaskellSolidSolnPhase.cpp
+++ b/src/thermo/MaskellSolidSolnPhase.cpp
@@ -28,7 +28,8 @@ MaskellSolidSolnPhase::MaskellSolidSolnPhase() :
m_cp0_R(2),
m_g0_RT(2),
m_s0_R(2),
- h_mixing(0.0)
+ h_mixing(0.0),
+ product_species_index(0)
{
}
//=====================================================================================================
@@ -40,7 +41,8 @@ MaskellSolidSolnPhase::MaskellSolidSolnPhase(const MaskellSolidSolnPhase& b) :
m_cp0_R(2),
m_g0_RT(2),
m_s0_R(2),
- h_mixing(0.0)
+ h_mixing(0.0),
+ product_species_index(0)
{
*this = b;
}
@@ -59,6 +61,26 @@ ThermoPhase* MaskellSolidSolnPhase::duplMyselfAsThermoPhase() const
return new MaskellSolidSolnPhase(*this);
}
//=====================================================================================================
+void MaskellSolidSolnPhase::
+getActivityConcentrations(doublereal* c) const
+{
+ std::vector pmv(m_kk);
+ getPartialMolarVolumes(&pmv[0]);
+ const doublereal* const dtmp = moleFractdivMMW();
+ const double mmw = meanMolecularWeight();
+ for (size_t k = 0; k < m_kk; k++) {
+ c[k] = dtmp[k] * mmw / pmv[k];
+ }
+}
+
+doublereal MaskellSolidSolnPhase::standardConcentration(size_t k) const
+{
+ std::vector pmv(m_kk);
+ getPartialMolarVolumes(&pmv[0]);
+ doublereal result = 1.0 / pmv[k];
+ return result;
+}
+
/********************************************************************
* Molar Thermodynamic Properties of the Solution
********************************************************************/
@@ -68,21 +90,21 @@ enthalpy_mole() const
{
_updateThermo();
const doublereal h0 = GasConstant * temperature() * mean_X(&m_h0_RT[0]);
- const doublereal r = moleFraction(0);
+ const doublereal r = moleFraction(product_species_index);
const doublereal fmval = fm(r);
return h0 + r * fmval * h_mixing;
}
//=====================================================================================================
doublereal xlogx(doublereal x)
{
- return x * std::log(x);
+ return x * std::log(x);
}
doublereal MaskellSolidSolnPhase::entropy_mole() const
{
_updateThermo();
const doublereal s0 = GasConstant * mean_X(&m_s0_R[0]);
- const doublereal r = moleFraction(0);
+ const doublereal r = moleFraction(product_species_index);
const doublereal fmval = fm(r);
const doublereal rfm = r * fmval;
return s0 + GasConstant * (xlogx(1-rfm) - xlogx(rfm) - xlogx(1-r-rfm) - xlogx((1-fmval)*r) - xlogx(1-r) - xlogx(r));
@@ -149,7 +171,7 @@ void MaskellSolidSolnPhase::
getChemPotentials(doublereal* mu) const
{
_updateThermo();
- const doublereal r = moleFraction(0);
+ const doublereal r = moleFraction(product_species_index);
const doublereal pval = p(r);
const doublereal fmval = fm(r);
const doublereal rfm = r * fmval;
@@ -157,8 +179,9 @@ getChemPotentials(doublereal* mu) const
const doublereal muDelta = -pval * h_mixing - GasConstant * temperature()
* std::log( (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval) * r) /
(std::pow(1 - r - rfm, 1 + pval) * (1 - r)) );
- mu[0] = RT * m_g0_RT[0] + muDelta;
- mu[1] = RT * m_g0_RT[1] - muDelta;
+ const int sign = (product_species_index == 0) ? 1 : -1;
+ mu[0] = RT * m_g0_RT[0] + sign * muDelta;
+ mu[1] = RT * m_g0_RT[1] - sign * muDelta;
}
void MaskellSolidSolnPhase::
@@ -252,6 +275,19 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string
throw CanteraError(subname.c_str(),
"Mixing enthalpy parameter not specified.");
}
+
+ if (thNode.hasChild("product_species")) {
+ XML_Node& scNode = thNode.child("product_species");
+ std::string product_species_name = scNode.value();
+ product_species_index = speciesIndex(product_species_name);
+ if( product_species_index == npos )
+ {
+ throw CanteraError(subname.c_str(),
+ "Species " + product_species_name + " not found.");
+ }
+ std::cout << "parsed product_species_index = " << product_species_index << std::endl;
+ }
+
} else {
throw CanteraError(subname.c_str(),
"Unspecified thermo model");
diff --git a/test/data/MaskellSolidSolnPhase_valid.xml b/test/data/MaskellSolidSolnPhase_valid.xml
index dfe4ee689..4b07a596f 100644
--- a/test/data/MaskellSolidSolnPhase_valid.xml
+++ b/test/data/MaskellSolidSolnPhase_valid.xml
@@ -12,7 +12,7 @@
-1000.
- 1.0
+ H(s)