From 5f01cc30bbb45da43fca49765e73bf3542c3daf5 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 1 Nov 2011 21:09:58 +0000 Subject: [PATCH] Worked on MolarityIonicVPSSTP, filling it in a little bit. Changed GibbsExcessVPSSTP so that getLnActivityCoefficient() is central. The straight activity coefficients may not be representable within machine limits. --- Cantera/src/thermo/GibbsExcessVPSSTP.cpp | 22 + Cantera/src/thermo/GibbsExcessVPSSTP.h | 7 + Cantera/src/thermo/MargulesVPSSTP.cpp | 10 +- Cantera/src/thermo/MargulesVPSSTP.h | 8 +- Cantera/src/thermo/MolarityIonicVPSSTP.cpp | 453 ++++++++++++++++++++- Cantera/src/thermo/MolarityIonicVPSSTP.h | 222 +++++++++- Cantera/src/thermo/RedlichKisterVPSSTP.cpp | 20 +- Cantera/src/thermo/RedlichKisterVPSSTP.h | 9 +- 8 files changed, 713 insertions(+), 38 deletions(-) diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp index 339b0edd9..8f357a9e5 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp @@ -233,7 +233,29 @@ namespace Cantera { } } + //==================================================================================================================== + // Get the array of non-dimensional molar-based activity coefficients at + // the current solution temperature, pressure, and solution concentration. + /* + * @param ac Output vector of activity coefficients. Length: m_kk. + */ + void GibbsExcessVPSSTP::getActivityCoefficients(doublereal * const ac) const { + + getLnActivityCoefficients(ac); + // Protect against roundoff when taking exponentials + for (int k = 0; k < m_kk; k++) { + if (ac[k] > 700.) { + ac[k] = exp(700.); + } else if (ac[k] < -700.) { + ac[k] = exp(-700); + } else { + ac[k] = exp(ac[k]); + } + } + } + //==================================================================================================================== + void GibbsExcessVPSSTP::getElectrochemPotentials(doublereal* mu) const { getChemPotentials(mu); double ve = Faraday * electricPotential(); diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.h b/Cantera/src/thermo/GibbsExcessVPSSTP.h index 521f7f344..8de94d8ff 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.h +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.h @@ -289,6 +289,13 @@ namespace Cantera { */ virtual void getActivities(doublereal* ac) const; + //! Get the array of non-dimensional molar-based ln activity coefficients at + //! the current solution temperature, pressure, and solution concentration. + /*! + * @param lnac Output vector of ln activity coefficients. Length: m_kk. + */ + virtual void getActivityCoefficients(doublereal * const ac) const; + //! Get the array of temperature derivatives of the log activity coefficients /*! diff --git a/Cantera/src/thermo/MargulesVPSSTP.cpp b/Cantera/src/thermo/MargulesVPSSTP.cpp index 15fc3edb2..ebe7c89aa 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.cpp +++ b/Cantera/src/thermo/MargulesVPSSTP.cpp @@ -376,12 +376,12 @@ namespace Cantera { return 0.0; } //==================================================================================================================== - // Get the array of non-dimensional molar-based activity coefficients at + // Get the array of non-dimensional molar-based ln activity coefficients at // the current solution temperature, pressure, and solution concentration. /* - * @param ac Output vector of activity coefficients. Length: m_kk. + * @param lnac Output vector of ln activity coefficients. Length: m_kk. */ - void MargulesVPSSTP::getActivityCoefficients(doublereal* ac) const { + void MargulesVPSSTP::getLnActivityCoefficients(doublereal* lnac) const { /* * Update the activity coefficients */ @@ -391,10 +391,10 @@ namespace Cantera { * take the exp of the internally storred coefficients. */ for (int k = 0; k < m_kk; k++) { - ac[k] = exp(lnActCoeff_Scaled_[k]); + lnac[k] = lnActCoeff_Scaled_[k]; } } - + //==================================================================================================================== /* * ------------ Partial Molar Properties of the Solution ------------ */ diff --git a/Cantera/src/thermo/MargulesVPSSTP.h b/Cantera/src/thermo/MargulesVPSSTP.h index ebeb9d5e2..5d675fb14 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.h +++ b/Cantera/src/thermo/MargulesVPSSTP.h @@ -536,15 +536,13 @@ namespace Cantera { */ virtual doublereal logStandardConc(int k=0) const; - //! Get the array of non-dimensional molar-based activity coefficients at + //! Get the array of non-dimensional molar-based ln activity coefficients at //! the current solution temperature, pressure, and solution concentration. /*! - * @param ac Output vector of activity coefficients. Length: m_kk. + * @param ac Output vector of ln activity coefficients. Length: m_kk. */ - virtual void getActivityCoefficients(doublereal* ac) const; + virtual void getLnActivityCoefficients(doublereal* lnac) const; - - //@} /// @name Partial Molar Properties of the Solution diff --git a/Cantera/src/thermo/MolarityIonicVPSSTP.cpp b/Cantera/src/thermo/MolarityIonicVPSSTP.cpp index e657aea3d..1e7831ba9 100644 --- a/Cantera/src/thermo/MolarityIonicVPSSTP.cpp +++ b/Cantera/src/thermo/MolarityIonicVPSSTP.cpp @@ -23,12 +23,13 @@ #include "MolarityIonicVPSSTP.h" - +#include "ThermoFactory.h" #include using namespace std; namespace Cantera { + static const double xxSmall = 1.0E-150; //==================================================================================================================== /* * Default constructor. @@ -44,6 +45,40 @@ namespace Cantera { numPassThroughSpecies_(0), neutralPBindexStart(0) { + } + //==================================================================================================================== + /* + * Working constructors + * + * The two constructors below are the normal way + * the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + */ + MolarityIonicVPSSTP::MolarityIonicVPSSTP(std::string inputFile, std::string id) : + GibbsExcessVPSSTP(), + PBType_(PBTYPE_PASSTHROUGH), + numPBSpecies_(m_kk), + indexSpecialSpecies_(-1), + numCationSpecies_(0), + numAnionSpecies_(0), + numPassThroughSpecies_(0), + neutralPBindexStart(0) + { + constructPhaseFile(inputFile, id); + } + //==================================================================================================================== + MolarityIonicVPSSTP::MolarityIonicVPSSTP(XML_Node& phaseRoot, std::string id) : + GibbsExcessVPSSTP(), + PBType_(PBTYPE_PASSTHROUGH), + numPBSpecies_(m_kk), + indexSpecialSpecies_(-1), + numCationSpecies_(0), + numAnionSpecies_(0), + numPassThroughSpecies_(0), + neutralPBindexStart(0) + { + constructPhaseXML(phaseRoot, id); } //==================================================================================================================== /* @@ -130,8 +165,120 @@ namespace Cantera { return 0; } - + //==================================================================================================================== + /* + * Import, construct, and initialize a phase + * specification from an XML tree into the current object. + * + * This routine is a precursor to constructPhaseXML(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void MolarityIonicVPSSTP::constructPhaseFile(std::string inputFile, std::string id) { + if ((int) inputFile.size() == 0) { + throw CanteraError("MolarityIonicVPSSTP:constructPhaseFile", + "input file is null"); + } + string path = findInputFile(inputFile); + std::ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("MolarityIonicVPSSTP:constructPhaseFile","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + XML_Node &phaseNode_XML = xml(); + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("MolarityIonicVPSSTP:constructPhaseFile", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + fxml_phase->copy(&phaseNode_XML); + constructPhaseXML(*fxml_phase, id); + delete fxml; + } + //==================================================================================================================== + /* + * Import, construct, and initialize a HMWSoln phase + * specification from an XML tree into the current object. + * + * Most of the work is carried out by the cantera base + * routine, importPhase(). That routine imports all of the + * species and element data, including the standard states + * of the species. + * + * Then, In this routine, we read the information + * particular to the specification of the activity + * coefficient model for the Pitzer parameterization. + * + * We also read information about the molar volumes of the + * standard states if present in the XML file. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void MolarityIonicVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id) { + string stemp; + if ((int) id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML", + "no thermo XML node"); + } + XML_Node& thermoNode = phaseNode.child("thermo"); + + /* + * Make sure that the thermo model is MolarityIonic + */ + stemp = thermoNode.attrib("model"); + string formString = lowercase(stemp); + if (formString != "molarityionicvpss") { + throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML", + "model name isn't MolarityIonicVPSS: " + formString); + + } + + /* + * Call the Cantera importPhase() function. This will import + * all of the species into the phase. This will also handle + * all of the solvent and solute standard states + */ + bool m_ok = importPhase(phaseNode, this); + if (!m_ok) { + throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML","importPhase failed "); + } + + } + //==================================================================================================================== /* * ------------ Molar Thermodynamic Properties ---------------------- */ @@ -140,16 +287,78 @@ namespace Cantera { /* * - Activities, Standard States, Activity Concentrations ----------- */ - + // This method returns an array of generalized concentrations + /* + * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / + * C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration + * defined below and \f$ a_k \f$ are activities used in the + * thermodynamic functions. These activity (or generalized) + * concentrations are used + * by kinetics manager classes to compute the forward and + * reverse rates of elementary reactions. Note that they may + * or may not have units of concentration --- they might be + * partial pressures, mole fractions, or surface coverages, + * for example. + * + * Here we define the activity concentrations as equal + * to the activities, because the standard concentration is 1. + * + * @param c Output array of generalized concentrations. The + * units depend upon the implementation of the + * reaction rate expressions within the phase. + */ + void MolarityIonicVPSSTP::getActivityConcentrations(doublereal* c) const { + getActivities(c); + } //==================================================================================================================== doublereal MolarityIonicVPSSTP::standardConcentration(int k) const { - err("standardConcentration"); return -1.0; } //==================================================================================================================== doublereal MolarityIonicVPSSTP::logStandardConc(int k) const { - err("logStandardConc"); - return -1.0; + return 0.0; + } + //==================================================================================================================== + // Get the array of non-dimensional molar-based activity coefficients at + // the current solution temperature, pressure, and solution concentration. + /* + * @param ac Output vector of activity coefficients. Length: m_kk. + */ + void MolarityIonicVPSSTP::getLnActivityCoefficients(doublereal* lnac) const { + /* + * Update the activity coefficients + */ + s_update_lnActCoeff(); + + /* + * take the exp of the internally storred coefficients. + */ + for (int k = 0; k < m_kk; k++) { + lnac[k] = lnActCoeff_Scaled_[k]; + } + } + //==================================================================================================================== + void MolarityIonicVPSSTP::getChemPotentials(doublereal* mu) const { + doublereal xx; + /* + * First get the standard chemical potentials in + * molar form. + * -> this requires updates of standard state as a function + * of T and P + */ + getStandardChemPotentials(mu); + /* + * Update the activity coefficients + */ + s_update_lnActCoeff(); + /* + * + */ + doublereal RT = GasConstant * temperature(); + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(moleFractions_[k], xxSmall); + mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]); + } } //==================================================================================================================== @@ -160,6 +369,145 @@ namespace Cantera { mu[k] += ve*charge(k); } } + + //==================================================================================================================== + // Returns an array of partial molar enthalpies for the species + // in the mixture. + /* + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * molality-based activity coefficent wrt temperature + * + * \f[ + * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * \f] + * + */ + void MolarityIonicVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const { + /* + * Get the nondimensional standard state enthalpies + */ + getEnthalpy_RT(hbar); + /* + * dimensionalize it. + */ + double T = temperature(); + double RT = GasConstant * T; + for (int k = 0; k < m_kk; k++) { + hbar[k] *= RT; + } + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnActCoeff(); + s_update_dlnActCoeff_dT(); + double RTT = RT * T; + for (int k = 0; k < m_kk; k++) { + hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k]; + } + } + //==================================================================================================================== + // Returns an array of partial molar heat capacities for the species + // in the mixture. + /* + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * activity coefficent wrt temperature + * + * \f[ + * ??????????? \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * \f] + * + */ + void MolarityIonicVPSSTP::getPartialMolarCp(doublereal* cpbar) const { + /* + * Get the nondimensional standard state entropies + */ + getCp_R(cpbar); + double T = temperature(); + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnActCoeff(); + s_update_dlnActCoeff_dT(); + + for (int k = 0; k < m_kk; k++) { + cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k]; + } + /* + * dimensionalize it. + */ + for (int k = 0; k < m_kk; k++) { + cpbar[k] *= GasConstant; + } + } + //==================================================================================================================== + // Returns an array of partial molar entropies for the species + // in the mixture. + /* + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * activity coefficent wrt temperature + * + * \f[ + * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * \f] + * + */ + void MolarityIonicVPSSTP::getPartialMolarEntropies(doublereal* sbar) const { + double xx; + /* + * Get the nondimensional standard state entropies + */ + getEntropy_R(sbar); + double T = temperature(); + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnActCoeff(); + s_update_dlnActCoeff_dT(); + + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(moleFractions_[k], xxSmall); + sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k]; + } + /* + * dimensionalize it. + */ + for (int k = 0; k < m_kk; k++) { + sbar[k] *= GasConstant; + } + } + // Return an array of partial molar volumes for the + // species in the mixture. Units: m^3/kmol. + /* + * Frequently, for this class of thermodynamics representations, + * the excess Volume due to mixing is zero. Here, we set it as + * a default. It may be overriden in derived classes. + * + * @param vbar Output vector of speciar partial molar volumes. + * Length = m_kk. units are m^3/kmol. + */ + void MolarityIonicVPSSTP::getPartialMolarVolumes(doublereal* vbar) const { + int iK; + /* + * Get the standard state values in m^3 kmol-1 + */ + getStandardVolumes(vbar); + for ( iK = 0; iK < m_kk; iK++ ){ + + vbar[iK] += 0.0; + } + } //==================================================================================================================== void MolarityIonicVPSSTP::calcPseudoBinaryMoleFractions() const { int k; @@ -241,6 +589,38 @@ namespace Cantera { } } //==================================================================================================================== + + // Update the activity coefficients + /* + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + * + */ + void MolarityIonicVPSSTP::s_update_lnActCoeff() const { + int k; + for (k = 0; k < m_kk; k++) { + lnActCoeff_Scaled_[k] = 0.0; + } + } + //==================================================================================================================== + void MolarityIonicVPSSTP::s_update_dlnActCoeff_dT() const { + + + } + //==================================================================================================================== + // Internal routine that calculates the derivative of the activity coefficients wrt + // the mole fractions. + /* + * This routine calculates the the derivative of the activity coefficients wrt to mole fraction + * with all other mole fractions held constant. This is strictly not permitted. However, if the + * resulting matrix is multiplied by a permissible deltaX vector then everything is ok. + * + * This is the natural way to handle concentration derivatives in this routine. + */ + void MolarityIonicVPSSTP::s_update_dlnActCoeff_dX_() const { + + } + //==================================================================================================================== /* * ------------ Partial Molar Properties of the Solution ------------ */ @@ -323,11 +703,70 @@ namespace Cantera { * with the correct id. */ void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { + std::string subname = "MolarityIonicVPSSTP::initThermoXML"; + std::string stemp; + /* + * Check on the thermo field. Must have: + * + */ + + XML_Node& thermoNode = phaseNode.child("thermo"); + std::string mStringa = thermoNode.attrib("model"); + std::string mString = lowercase(mStringa); + if (mString != "molarityionicvpss") { + throw CanteraError(subname.c_str(), + "Unknown thermo model: " + mStringa + " - This object only knows \"MolarityIonicVPSS\" "); + } + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + XML_Node *acNodePtr = 0; + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& acNode = thermoNode.child("activityCoefficients"); + acNodePtr = &acNode; + std::string mStringa = acNode.attrib("model"); + std::string mString = lowercase(mStringa); + // if (mString != "redlich-kister") { + // throw CanteraError(subname.c_str(), + // "Unknown activity coefficient model: " + mStringa); + //} + int n = acNodePtr->nChildren(); + for (int i = 0; i < n; i++) { + XML_Node &xmlACChild = acNodePtr->child(i); + stemp = xmlACChild.name(); + std::string nodeName = lowercase(stemp); + /* + * Process a binary interaction + */ + if (nodeName == "binaryneutralspeciesparameters") { + readXMLBinarySpecies(xmlACChild); + } + } + } - + + /* + * Go down the chain + */ GibbsExcessVPSSTP::initThermoXML(phaseNode, id); } //==================================================================================================================== + // Process an XML node called "binaryNeutralSpeciesParameters" + /* + * This node contains all of the parameters necessary to describe + * a single binary interaction. This function reads the XML file and writes the coefficients + * it finds to an internal data structures. + */ + void MolarityIonicVPSSTP::readXMLBinarySpecies(XML_Node &xmLBinarySpecies) { + std::string xname = xmLBinarySpecies.name(); + + } + //==================================================================================================================== /* * Format a summary of the mixture state for output. */ diff --git a/Cantera/src/thermo/MolarityIonicVPSSTP.h b/Cantera/src/thermo/MolarityIonicVPSSTP.h index dd6c339e8..e75531a8e 100644 --- a/Cantera/src/thermo/MolarityIonicVPSSTP.h +++ b/Cantera/src/thermo/MolarityIonicVPSSTP.h @@ -38,8 +38,7 @@ namespace Cantera { * thermodynamic properties that are further based on * expressing the Excess Gibbs free energy as a function of * the mole fractions (or pseudo mole fractions) of the consitituents. - * This category is the workhorse for describing ionic systems which - * are not on the molality scale. + * This category is the workhorse for describing ionic systems which are not on the molality scale. * * This class adds additional functions onto the %ThermoPhase interface * that handles the calculation of the excess Gibbs free energy. The %ThermoPhase @@ -56,7 +55,7 @@ namespace Cantera { * it is expected that there exists a charge balance at all times. * One of the ions must be a "special ion" in the sense that its' thermodynamic * functions are set to zero, and the thermo functions of all other - * ions are based on a valuation relative to the special ion. + * ions are based on a valuation relative to that special ion. * */ class MolarityIonicVPSSTP : public GibbsExcessVPSSTP { @@ -74,6 +73,31 @@ namespace Cantera { */ MolarityIonicVPSSTP(); + //! Construct and initialize a MolarityIonicVPSSTP ThermoPhase object + //! directly from an xml input file + /*! + * Working constructors + * + * The two constructors below are the normal way the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the XML database to get the info for the phase. + * + * @param inputFile Name of the input file containing the phase XML data + * to set up the object + * @param id ID of the phase in the input file. Defaults to the + * empty string. + */ + MolarityIonicVPSSTP(std::string inputFile, std::string id = ""); + + //! Construct and initialize a MolarityIonicVPSSTP ThermoPhase object + //! directly from an XML database + /*! + * @param phaseRef XML phase node containing the description of the phase + * @param id id attribute containing the name of the phase. + * (default is the empty string) + */ + MolarityIonicVPSSTP(XML_Node& phaseRef, std::string id = ""); + + //! Copy constructor /*! * Note this stuff will not work until the underlying phase @@ -118,7 +142,46 @@ namespace Cantera { */ virtual int eosType() const; - + //! Initialization of a phase using an xml file + /*! + * This routine is a precursor to + * routine, which does most of the work. + * + * @param inputFile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void constructPhaseFile(std::string inputFile, std::string id); + + //! Import and initialize a phase + //! specification in an XML tree into the current object. + /*! + * Here we read an XML description of the phase. + * We import descriptions of the elements that make up the + * species in a phase. + * We import information about the species, including their + * reference state thermodynamic polynomials. We then freeze + * the state of the species. + * + * Then, we read the species molar volumes from the xml + * tree to finish the initialization. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void constructPhaseXML(XML_Node& phaseNode, std::string id); /** * @} @@ -165,7 +228,24 @@ namespace Cantera { * @{ */ - + //! This method returns an array of generalized concentrations + /*! + * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / + * C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration + * defined below and \f$ a_k \f$ are activities used in the + * thermodynamic functions. These activity (or generalized) + * concentrations are used + * by kinetics manager classes to compute the forward and + * reverse rates of elementary reactions. Note that they may + * or may not have units of concentration --- they might be + * partial pressures, mole fractions, or surface coverages, + * for example. + * + * @param c Output array of generalized concentrations. The + * units depend upon the implementation of the + * reaction rate expressions within the phase. + */ + virtual void getActivityConcentrations(doublereal* c) const; /** @@ -191,6 +271,15 @@ namespace Cantera { */ virtual doublereal logStandardConc(int k=0) const; + //! Get the array of non-dimensional molar-based ln activity coefficients at + //! the current solution temperature, pressure, and solution concentration. + /*! + * @param lnac Output vector of ln activity coefficients. Length: m_kk. + */ + virtual void getLnActivityCoefficients(doublereal* ac) const; + + + @@ -198,6 +287,16 @@ namespace Cantera { /// @name Partial Molar Properties of the Solution //@{ + //! Get the species chemical potentials. Units: J/kmol. + /*! + * This function returns a vector of chemical potentials of the + * species in solution at the current temperature, pressure + * and mole fraction of the solution. + * + * @param mu Output vector of species chemical + * potentials. Length: m_kk. Units: J/kmol + */ + virtual void getChemPotentials(doublereal* mu) const; /** * Get the species electrochemical potentials. @@ -212,7 +311,79 @@ namespace Cantera { */ void getElectrochemPotentials(doublereal* mu) const; + //! Returns an array of partial molar enthalpies for the species + //! in the mixture. + /*! + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * molality-based activity coefficent wrt temperature + * + * \f[ + * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * \f] + * + * @param hbar Vector of returned partial molar enthalpies + * (length m_kk, units = J/kmol) + */ + virtual void getPartialMolarEnthalpies(doublereal* hbar) const; + + //! Returns an array of partial molar entropies for the species + //! in the mixture. + /*! + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * activity coefficent wrt temperature + * + * \f[ + * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * - R \ln( \gamma_k X_k) + * - R T \frac{d \ln(\gamma_k) }{dT} + * \f] + * + * @param sbar Vector of returned partial molar entropies + * (length m_kk, units = J/kmol/K) + */ + virtual void getPartialMolarEntropies(doublereal* sbar) const; + + //! Returns an array of partial molar entropies for the species + //! in the mixture. + /*! + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * activity coefficent wrt temperature + * + * \f[ + * ??????????????? + * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * - R \ln( \gamma_k X_k) + * - R T \frac{d \ln(\gamma_k) }{dT} + * ??????????????? + * \f] + * + * @param cpbar Vector of returned partial molar heat capacities + * (length m_kk, units = J/kmol/K) + */ + virtual void getPartialMolarCp(doublereal* cpbar) const; + //! Return an array of partial molar volumes for the + //! species in the mixture. Units: m^3/kmol. + /*! + * Frequently, for this class of thermodynamics representations, + * the excess Volume due to mixing is zero. Here, we set it as + * a default. It may be overriden in derived classes. + * + * @param vbar Output vector of speciar partial molar volumes. + * Length = m_kk. units are m^3/kmol. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + + //@} /// @name Properties of the Standard State of the Species in the Solution //@{ @@ -318,11 +489,48 @@ namespace Cantera { private: - //! Initialize lengths of local variables after all species have - //! been identified. + //! Initialize lengths of local variables after all species have been identified. void initLengths(); + + //! Process an XML node called "binaryNeutralSpeciesParameters" + /*! + * This node contains all of the parameters necessary to describe + * the Redlich-Kister model for a particular binary interaction. + * This function reads the XML file and writes the coefficients + * it finds to an internal data structures. + * + * @param xmlBinarySpecies Reference to the XML_Node named "binaryNeutralSpeciesParameters" + * containing the binary interaction + */ + void readXMLBinarySpecies(XML_Node &xmlBinarySpecies); + //! Update the activity coefficients + /*! + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + */ + void s_update_lnActCoeff() const; + + //! Update the derivative of the log of the activity coefficients wrt T + /*! + * This function will be called to update the internally storred + * derivative of the natural logarithm of the activity coefficients + * wrt temperature. + */ + void s_update_dlnActCoeff_dT() const; + + //! Internal routine that calculates the derivative of the activity coefficients wrt + //! the mole fractions. + /*! + * This routine calculates the the derivative of the activity coefficients wrt to mole fraction + * with all other mole fractions held constant. This is strictly not permitted. However, if the + * resulting matrix is multiplied by a permissible deltaX vector then everything is ok. + * + * This is the natural way to handle concentration derivatives in this routine. + */ + void s_update_dlnActCoeff_dX_() const; + private: //! Error function diff --git a/Cantera/src/thermo/RedlichKisterVPSSTP.cpp b/Cantera/src/thermo/RedlichKisterVPSSTP.cpp index a97ee19d8..b6dc57f3d 100644 --- a/Cantera/src/thermo/RedlichKisterVPSSTP.cpp +++ b/Cantera/src/thermo/RedlichKisterVPSSTP.cpp @@ -387,7 +387,7 @@ namespace Cantera { /* * @param ac Output vector of activity coefficients. Length: m_kk. */ - void RedlichKisterVPSSTP::getActivityCoefficients(doublereal* ac) const { + void RedlichKisterVPSSTP::getLnActivityCoefficients(doublereal* lnac) const { /* * Update the activity coefficients */ @@ -397,7 +397,7 @@ namespace Cantera { * take the exp of the internally storred coefficients. */ for (int k = 0; k < m_kk; k++) { - ac[k] = exp(lnActCoeff_Scaled_[k]); + lnac[k] = lnActCoeff_Scaled_[k]; } } //==================================================================================================================== @@ -728,6 +728,12 @@ namespace Cantera { doublereal RT = GasConstant * T; fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + + /* + * Scaling: I moved the division of RT higher so that we are always dealing with G/RT dimensionless terms + * within the routine. There is a severe problem with roundoff error in these calculations. The + * dimensionless terms help. + */ for (int i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; @@ -744,7 +750,7 @@ namespace Cantera { doublereal sumMm1 = 0.0; doublereal sum2 = 0.0; for (m = 0; m < N; m++) { - doublereal A_ge = he_vec[m] - T * se_vec[m]; + doublereal A_ge = (he_vec[m] - T * se_vec[m]) / RT; sum += A_ge * poly; sum2 += A_ge * (m + 1) * poly; poly *= deltaX; @@ -771,7 +777,7 @@ namespace Cantera { double polyk = 1.0; double fac = 2.0 * XA - 1.0; for (m = 0; m < N; m++) { - doublereal A_ge = he_vec[m] - T * se_vec[m]; + doublereal A_ge = (he_vec[m] - T * se_vec[m]) / RT; lnA += A_ge * oneMXA * oneMXA * polyk * (1.0 + 2.0 * XA * m / fac); lnB += A_ge * XA * XA * polyk * (1.0 - 2.0 * oneMXA * m / fac); polyk *= fac; @@ -783,9 +789,7 @@ namespace Cantera { #endif } - for (k = 0; k < m_kk; k++) { - lnActCoeff_Scaled_[k] /= RT; - } + } //=================================================================================================================== // Update the derivative of the log of the activity coefficients wrt T @@ -1117,7 +1121,7 @@ namespace Cantera { double poly1mk = fac; for (m = 0; m < N; m++) { - doublereal A_ge = he_vec[m]; + doublereal A_ge = he_vec[m] - T * se_vec[m]; Volts += A_ge * ( polykp1 - (2.0 * XA * m * (1.0-XA) ) / poly1mk ); polykp1 *= fac; poly1mk /= fac; diff --git a/Cantera/src/thermo/RedlichKisterVPSSTP.h b/Cantera/src/thermo/RedlichKisterVPSSTP.h index d11e9e2bd..d4d927d62 100644 --- a/Cantera/src/thermo/RedlichKisterVPSSTP.h +++ b/Cantera/src/thermo/RedlichKisterVPSSTP.h @@ -533,15 +533,12 @@ namespace Cantera { */ virtual doublereal logStandardConc(int k=0) const; - //! Get the array of non-dimensional molar-based activity coefficients at + //! Get the array of non-dimensional molar-based ln activity coefficients at //! the current solution temperature, pressure, and solution concentration. /*! - * @param ac Output vector of activity coefficients. Length: m_kk. + * @param lnac Output vector of ln activity coefficients. Length: m_kk. */ - virtual void getActivityCoefficients(doublereal* ac) const; - - - + virtual void getLnActivityCoefficients(doublereal* ac) const; //@} /// @name Partial Molar Properties of the Solution