From 7abb459bbe2128f08db0c55c4856392a7a252465 Mon Sep 17 00:00:00 2001 From: Victor Brunini Date: Thu, 6 Mar 2014 22:30:19 +0000 Subject: [PATCH] Implement MaskellSolidSolnPhase::standardConcentrations() and getActivityConcentrations() --- .../cantera/thermo/MaskellSolidSolnPhase.h | 76 ++++++++++++++----- src/thermo/MaskellSolidSolnPhase.cpp | 52 +++++++++++-- test/data/MaskellSolidSolnPhase_valid.xml | 2 +- 3 files changed, 103 insertions(+), 27 deletions(-) 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)