diff --git a/Cantera/src/base/xml.cpp b/Cantera/src/base/xml.cpp index 2ea75093c..1484c1b6f 100644 --- a/Cantera/src/base/xml.cpp +++ b/Cantera/src/base/xml.cpp @@ -52,9 +52,8 @@ namespace Cantera { */ XML_Error(int line=0) : m_line(line), - m_msg(0) + m_msg("Error in XML file") { - m_msg = "Error in XML file"; if (line > 0) { m_msg += " at line " + int2str(line+1); } diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 664544779..5683e640c 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -1571,6 +1571,7 @@ namespace VCSnonideal { for (eT = 0; eT < nebase; eT++) { ename = tPhase->elementName(eT); m_elementNames[e] = ename; + m_elementType[e] = tPhase->elementType(eT); e++; } diff --git a/Cantera/src/equil/vcs_defs.h b/Cantera/src/equil/vcs_defs.h index 534d5ca54..2c50ffaa1 100644 --- a/Cantera/src/equil/vcs_defs.h +++ b/Cantera/src/equil/vcs_defs.h @@ -314,6 +314,11 @@ namespace VCSnonideal { * constraint to one category. * @{ */ + + + //! An element constraint that is current turned off +#define VCS_ELEM_TYPE_TURNEDOFF -1 + //! Normal element constraint consisting of positive coefficients for the //! formula matrix. /*! @@ -335,11 +340,37 @@ namespace VCSnonideal { */ #define VCS_ELEM_TYPE_CHARGENEUTRALITY 2 + //! Constraint associated with maintaing a fixed lattice stoichiometry int eh + //! solids + /*! + * The constraint may have positive or negative values. The lattice 0 species will + * have negative values while higher lattices will have positive values + */ +#define VCS_ELEM_TYPE_LATTICERATIO 3 + + //! Constraint associated with maintaining frozen kinetic equilibria in + //! some functional groups within molecules + /*! + * We seek here to say that some functional groups or ionic states should be + * treated as if they are separate elements given the time scale of the problem. + * This will be abs positive constraint. We have not implemented any examples yet. + * A requirement will be that we must be able to add and subtract these contraints. + */ +#define VCS_ELEM_TYPE_KINETICFROZEN 4 + + //! Constraint associated with the maintenance of a surface phase + /*! + * We don't have any examples of this yet either. However, surfaces only exist + * because they are interfaces between bulk layers. If we want to treat surfaces + * within thermodynamic systems we must come up with a way to constrain their total + * number. + */ +#define VCS_ELEM_TYPE_SURFACECONSTRAINT 5 //! Other constraint equations /*! * currently there are none */ -#define VCS_ELEM_TYPE_OTHERCONSTRAINT 3 +#define VCS_ELEM_TYPE_OTHERCONSTRAINT 6 //@} /*! diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 75444a566..251b1496d 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -518,9 +518,18 @@ namespace VCSnonideal { * FormulaMatrix[] -> Copy the formula matrix over */ for (i = 0; i < nspecies; i++) { + bool nonzero = false; for (j = 0; j < nelements; j++) { + if (pub->FormulaMatrix[j][i] != 0.0) { + nonzero = true; + } m_formulaMatrix[j][i] = pub->FormulaMatrix[j][i]; } + if (!nonzero) { + plogf("vcs_prob_specifyFully:: species %d %s has a zero formula matrix!\n", i, + pub->SpName[i].c_str()); + return VCS_PUB_BAD; + } } /* @@ -573,17 +582,31 @@ namespace VCSnonideal { /* * Formulate the Goal Element Abundance Vector */ + double sum; if (pub->gai.size() != 0) { - for (i = 0; i < nelements; i++) m_elemAbundancesGoal[i] = pub->gai[i]; + for (i = 0; i < nelements; i++) { + m_elemAbundancesGoal[i] = pub->gai[i]; + if (pub->m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO) { + if (m_elemAbundancesGoal[i] < 1.0E-10) { + m_elemAbundancesGoal[i] = 0.0; + } + } + } } else { if (m_doEstimateEquil == 0) { for (j = 0; j < nelements; j++) { m_elemAbundancesGoal[j] = 0.0; for (kspec = 0; kspec < nspecies; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + sum += m_molNumSpecies_old[kspec]; m_elemAbundancesGoal[j] += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec]; } } + if (pub->m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO) { + if (m_elemAbundancesGoal[j] < 1.0E-10 * sum) { + m_elemAbundancesGoal[j] = 0.0; + } + } } } else { plogf("%sElement Abundances, m_elemAbundancesGoal[], not specified\n", ser); diff --git a/Cantera/src/numerics/ResidEval.h b/Cantera/src/numerics/ResidEval.h index da65c4cfa..3e25f8129 100644 --- a/Cantera/src/numerics/ResidEval.h +++ b/Cantera/src/numerics/ResidEval.h @@ -31,7 +31,7 @@ namespace Cantera { * Virtual base class for DAE residual function evaluators. * Classes derived from ResidEval evaluate the residual function * \f[ - \vec{F}(t,\vec{y}, \vec{y^\prime}) + * \vec{F}(t,\vec{y}, \vec{y^\prime}) * \f] * The DAE solver attempts to find a solution y(t) such that F = 0. * @ingroup DAE_Group @@ -59,13 +59,28 @@ namespace Cantera { return c_NONE; } + //! Initialization function + virtual void initSizes() + { + int neq = nEquations(); + m_alg.resize(neq, 0); + } + /** * Specify that solution component k is purely algebraic - * that is, the derivative of this component does not appear * in the residual function. */ - virtual void setAlgebraic(const int k) { m_alg[k] = 1; } - virtual bool isAlgebraic(const int k) {return (m_alg[k] == 1); } + virtual void setAlgebraic(const int k) { + if ((int) m_alg.size() < (k+1)) { + initSizes(); + } + m_alg[k] = 1; + } + + virtual bool isAlgebraic(const int k) { + return (m_alg[k] == 1); + } /** @@ -106,6 +121,7 @@ namespace Cantera { */ virtual int getInitialConditions(const doublereal t0, doublereal * const y, doublereal * const ydot) { + initSizes(); throw CanteraError("ResidEval::GetInitialConditions()", "base class called"); return 1; } @@ -143,7 +159,12 @@ namespace Cantera { protected: - std::map m_alg; + //! Mapping vector that stores whether a degree of freedom is a DAE or not + /*! + * The first index is the equation number. The second index is 1 if it is a DAE, + * and zero if it is not. + */ + std::vector m_alg; std::map m_constrain; private: diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h index 08a23fc35..a1799253e 100644 --- a/Cantera/src/numerics/ResidJacEval.h +++ b/Cantera/src/numerics/ResidJacEval.h @@ -60,7 +60,6 @@ namespace Cantera { public: - //!Default constructor /*! * @param atol Initial value of the global tolerance (defaults to 1.0E-13) @@ -208,7 +207,6 @@ namespace Cantera { const doublereal * const y, const doublereal * const ydot); - //! Return a vector of delta y's for calculation of the numerical Jacobian /*! * There is a default algorithm provided. @@ -233,7 +231,6 @@ namespace Cantera { doublereal * const delta_y, const doublereal * const solnWeights = 0); - //! Returns a vector of column scale factors that can be used to column scale Jacobians. /*! * Default to yScales[] = 1.0 @@ -322,7 +319,6 @@ namespace Cantera { SquareMatrix &J, doublereal * const resid); - protected: //! constant value of atol @@ -335,5 +331,3 @@ namespace Cantera { #endif - - diff --git a/Cantera/src/thermo/Constituents.cpp b/Cantera/src/thermo/Constituents.cpp index 2d42d9564..ff7446d80 100644 --- a/Cantera/src/thermo/Constituents.cpp +++ b/Cantera/src/thermo/Constituents.cpp @@ -118,6 +118,9 @@ namespace Cantera { return m_Elements->atomicNumber(m); } + int Constituents::elementType(int m) const{ + return m_Elements->elementType(m); + } /* * Add an element to the set. @@ -155,9 +158,9 @@ namespace Cantera { */ void Constituents:: addUniqueElement(const std::string& symbol, doublereal weight, - int atomicNumber, doublereal entropy298) + int atomicNumber, doublereal entropy298, int elem_type) { - m_Elements->addUniqueElement(symbol, weight, atomicNumber, entropy298); + m_Elements->addUniqueElement(symbol, weight, atomicNumber, entropy298, elem_type); } void Constituents:: @@ -291,10 +294,10 @@ namespace Cantera { m_speciesNames.push_back(name); m_speciesCharge.push_back(charge); m_speciesSize.push_back(size); - int m_mm = m_Elements->nElements(); + int ne = m_Elements->nElements(); // Create a changeable copy of the element composition. We now change the charge potentially - vector_fp compNew(m_mm); - for (int m = 0; m < m_mm; m++) { + vector_fp compNew(ne); + for (int m = 0; m < ne; m++) { compNew[m] = comp[m]; } double wt = 0.0; @@ -313,30 +316,17 @@ namespace Cantera { } } } else { - m_Elements->m_elementsFrozen = false; - addUniqueElement("E", 0.000545, 0, 0.0); - m_Elements->m_elementsFrozen = true; - m_mm = m_Elements->nElements(); - if (m_kk > 0) { - vector_fp old(m_speciesComp); - m_speciesComp.resize(m_kk*m_mm, 0.0); - for (int k = 0; k < m_kk; k++) { - int m_old = m_mm - 1; - for (int m = 0; m < m_old; m++) { - m_speciesComp[k * m_mm + m] = old[k * (m_old) + m]; - } - m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0; - } - } + addUniqueElementAfterFreeze("E", 0.000545, 0, 0.0, CT_ELEM_TYPE_ELECTRONCHARGE); + ne = m_Elements->nElements(); eindex = m_Elements->elementIndex("E"); - compNew.resize(m_mm); - compNew[m_mm-1] = - charge; + compNew.resize(ne); + compNew[ne - 1] = - charge; //comp[eindex] = -charge; // throw CanteraError("Constituents::addSpecies", // "Element List doesn't include E, yet this species has charge:" + name); } } - for (int m = 0; m < m_mm; m++) { + for (int m = 0; m < ne; m++) { m_speciesComp.push_back(compNew[m]); wt += compNew[m] * aw[m]; } @@ -470,6 +460,8 @@ namespace Cantera { return m_speciesComp[m_mm * k + m]; } + + //==================================================================================================================== /* * * getAtoms() @@ -485,6 +477,39 @@ namespace Cantera { } } + + + //==================================================================================================================== + int Constituents::addUniqueElementAfterFreeze(const std::string& symbol, doublereal weight, int atomicNumber, + doublereal entropy298, int elem_type) + { + int ii = elementIndex(symbol); + if (ii != -1) { + return ii; + } + // Check to see that the element isn't really in the list + m_Elements->m_elementsFrozen = false; + addUniqueElement(symbol, weight, atomicNumber, entropy298, elem_type); + m_Elements->m_elementsFrozen = true; + int m_mm = m_Elements->nElements(); + ii = elementIndex(symbol); + if (ii != m_mm-1) { + throw CanteraError("Constituents::addElementAfterFreeze()", "confused"); + } + if (m_kk > 0) { + vector_fp old(m_speciesComp); + m_speciesComp.resize(m_kk*m_mm, 0.0); + for (int k = 0; k < m_kk; k++) { + int m_old = m_mm - 1; + for (int m = 0; m < m_old; m++) { + m_speciesComp[k * m_mm + m] = old[k * (m_old) + m]; + } + m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0; + } + } + return ii; + } + //==================================================================================================================== /* * This copy constructor just calls the assignment operator * for this class. diff --git a/Cantera/src/thermo/Constituents.h b/Cantera/src/thermo/Constituents.h index be134749f..3a31522c3 100644 --- a/Cantera/src/thermo/Constituents.h +++ b/Cantera/src/thermo/Constituents.h @@ -110,7 +110,6 @@ namespace Cantera { /// exception, ElementRangeError, is thrown. std::string elementName(int m) const; - /// Index of element named 'name'. /// The index is an integer /// assigned to each element in the order it was added, @@ -141,6 +140,8 @@ namespace Cantera { */ int atomicNumber(int m) const; + int elementType(int m) const; + /// Return a read-only reference to the vector of element names. const std::vector& elementNames() const; @@ -193,7 +194,7 @@ namespace Cantera { */ void addUniqueElement(const std::string& symbol, doublereal weight, int atomicNumber = 0, - doublereal entropy298 = ENTROPY298_UNKNOWN); + doublereal entropy298 = ENTROPY298_UNKNOWN, int elem_type = CT_ELEM_TYPE_ABSPOS); //! Adde an element, checking for uniqueness /*! @@ -216,6 +217,23 @@ namespace Cantera { /// True if freezeElements has been called. bool elementsFrozen(); + //! Add an element after the elements have been frozen, checking for uniqueness + /*! + * The uniqueness is checked by comparing the string symbol. If + * not unique, nothing is done. + * + * @param symbol String symbol of the element + * @param weight Atomic weight of the element (kg kmol-1). + * @param atomicNumber Atomic number of the element (unitless) + * @param entropy298 Entropy of the element at 298 K and 1 bar + * in its most stable form. The default is + * the value ENTROPY298_UNKNOWN, which is + * interpreted as an unknown, and if used + * will cause Cantera to throw an error. + */ + int addUniqueElementAfterFreeze(const std::string& symbol, doublereal weight, int atomicNumber, + doublereal entropy298 = ENTROPY298_UNKNOWN, int elem_type = CT_ELEM_TYPE_ABSPOS); + //@} /// Returns the number of species in the phase diff --git a/Cantera/src/thermo/DebyeHuckel.h b/Cantera/src/thermo/DebyeHuckel.h index 16e3b5dc8..f2ce3d375 100644 --- a/Cantera/src/thermo/DebyeHuckel.h +++ b/Cantera/src/thermo/DebyeHuckel.h @@ -1255,13 +1255,6 @@ namespace Cantera { * -------------- Utilities ------------------------------- */ - - /** - * Return a reference to the species thermodynamic property - * manager. @todo This method will fail if no species thermo - * manager has been installed. - */ - SpeciesThermo& speciesThermo() { return *m_spthermo; } //! Initialize the object's internal lengths after species are set /** diff --git a/Cantera/src/thermo/Elements.cpp b/Cantera/src/thermo/Elements.cpp index 2c9e61bdc..7641a2365 100644 --- a/Cantera/src/thermo/Elements.cpp +++ b/Cantera/src/thermo/Elements.cpp @@ -222,6 +222,7 @@ namespace Cantera { Elements::Elements() : m_mm(0), m_elementsFrozen(false), + m_elem_type(0), numSubscribers(0) { } @@ -256,7 +257,7 @@ namespace Cantera { m_atomicNumbers = right.m_atomicNumbers; m_elementNames = right.m_elementNames; m_entropy298 = right.m_entropy298; - + m_elem_type = right.m_elem_type; numSubscribers = 0; return *this; @@ -335,7 +336,43 @@ namespace Cantera { AssertTrace(m >= 0 && m < m_mm); return (m_entropy298[m]); } - + //==================================================================================================================== + //! Return the element constraint type + /*! + * Possible types include: + * + * CT_ELEM_TYPE_TURNEDOFF -1 + * CT_ELEM_TYPE_ABSPOS 0 + * CT_ELEM_TYPE_ELECTRONCHARGE 1 + * CT_ELEM_TYPE_CHARGENEUTRALITY 2 + * CT_ELEM_TYPE_LATTICERATIO 3 + * CT_ELEM_TYPE_KINETICFROZEN 4 + * CT_ELEM_TYPE_SURFACECONSTRAINT 5 + * CT_ELEM_TYPE_OTHERCONSTRAINT 6 + * + * The default is CT_ELEM_TYPE_ABSPOS + */ + int Elements::elementType(int m) const + { + return m_elem_type[m]; + } + //==================================================================================================================== + // Change the element type of the mth constraint + /* + * Reassigns an element type + * + * @param m Element index + * @param elem_type New elem type to be assigned + * + * @return Returns the old element type + */ + int Elements::changeElementType(int m, int elem_type) + { + int old = m_elem_type[m]; + m_elem_type[m] = elem_type; + return old; + } + //==================================================================================================================== /* * * Add an element to the current set of elements in the current object. @@ -367,16 +404,22 @@ namespace Cantera { #ifdef USE_DGG_CODE m_definedElements[symbol] = nElements() + 1; #endif + if (symbol == "E") { + m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE); + } else { + m_elem_type.push_back(CT_ELEM_TYPE_ABSPOS); + } + m_mm++; } - + //=========================================================================================================== void Elements:: addElement(const XML_Node& e) { doublereal weight = atof(e["atomicWt"].c_str()); string symbol = e["name"]; addElement(symbol, weight); } - + //=========================================================================================================== /* * addUniqueElement(): * @@ -393,7 +436,7 @@ namespace Cantera { #ifdef USE_DGG_CODE void Elements:: addUniqueElement(const std::string& symbol, doublereal weight, int atomicNumber, - doublereal entropy298) + doublereal entropy298, int elem_type) { if (m_elementsFrozen) throw ElementsFrozen("addElement"); @@ -413,13 +456,17 @@ namespace Cantera { m_elementNames.push_back(symbol); m_atomicNumbers.push_back(atomicNumber); m_entropy298.push_back(entropy298); + if (symbol == "E") { + m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE); + } else { + m_elem_type.push_back(elem_type); + } m_mm++; } else { if (m_atomicWeights[i] != weight) { throw CanteraError("AddUniqueElement", - "Duplicate Elements (" + symbol + - ") have different weights"); + "Duplicate Elements (" + symbol + ") have different weights"); } } } @@ -427,7 +474,8 @@ namespace Cantera { #else void Elements:: addUniqueElement(const std::string& symbol, - doublereal weight, int atomicNumber, doublereal entropy298) + doublereal weight, int atomicNumber, doublereal entropy298, + int elem_type) { if (weight == -12345.0) { weight = LookupWtElements(symbol); @@ -458,12 +506,16 @@ namespace Cantera { m_elementNames.push_back(symbol); m_atomicNumbers.push_back(atomicNumber); m_entropy298.push_back(entropy298); + if (symbol == "E") { + m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE); + } else { + m_elem_type.push_back(elem_type); + } m_mm++; } else { if (m_atomicWeights[i] != weight) { throw CanteraError("AddUniqueElement", - "Duplicate Elements (" + symbol + - ") have different weights"); + "Duplicate Elements (" + symbol + ") have different weights"); } } } @@ -506,6 +558,8 @@ namespace Cantera { m_mm = 0; m_atomicWeights.resize(0); m_elementNames.resize(0); + m_entropy298.resize(0); + m_elem_type.resize(0); m_elementsFrozen = false; } diff --git a/Cantera/src/thermo/Elements.h b/Cantera/src/thermo/Elements.h index fa0e5dc2b..8f738acad 100644 --- a/Cantera/src/thermo/Elements.h +++ b/Cantera/src/thermo/Elements.h @@ -26,6 +26,74 @@ namespace Cantera { class XML_Node; class ElementRangeError; + /*! + * @name Types of Element Constraint Equations + * + * There may be several different types of element constraints handled + * by the equilibrium program and by Cantera in other contexts. + * These defines are used to assign each constraint to one category. + * @{ + */ + + //! An element constraint that is current turned off +#define CT_ELEM_TYPE_TURNEDOFF -1 + + //! Normal element constraint consisting of positive coefficients for the + //! formula matrix. + /*! + * All species have positive coefficients within the formula matrix. + * With this constraint, we may employ various strategies to handle + * small values of the element number successfully. + */ +#define CT_ELEM_TYPE_ABSPOS 0 + + //! This refers to conservation of electrons + /*! + * Electrons may have positive or negative values in the Formula matrix. + */ +#define CT_ELEM_TYPE_ELECTRONCHARGE 1 + + //! This refers to a charge neutrality of a single phase + /*! + * Charge neutrality may have positive or negative values in the Formula matrix. + */ +#define CT_ELEM_TYPE_CHARGENEUTRALITY 2 + + //! Constraint associated with maintaing a fixed lattice stoichiometry int eh + //! solids + /*! + * The constraint may have positive or negative values. The lattice 0 species will + * have negative values while higher lattices will have positive values + */ +#define CT_ELEM_TYPE_LATTICERATIO 3 + + //! Constraint associated with maintaining frozen kinetic equilibria in + //! some functional groups within molecules + /*! + * We seek here to say that some functional groups or ionic states should be + * treated as if they are separate elements given the time scale of the problem. + * This will be abs positive constraint. We have not implemented any examples yet. + * A requirement will be that we must be able to add and subtract these contraints. + */ +#define CT_ELEM_TYPE_KINETICFROZEN 4 + + //! Constraint associated with the maintenance of a surface phase + /*! + * We don't have any examples of this yet either. However, surfaces only exist + * because they are interfaces between bulk layers. If we want to treat surfaces + * within thermodynamic systems we must come up with a way to constrain their total + * number. + */ +#define CT_ELEM_TYPE_SURFACECONSTRAINT 5 + + //! Other constraint equations + /*! + * currently there are none + */ +#define CT_ELEM_TYPE_OTHERCONSTRAINT 6 + //@} + + //! Positive number indicating we don't know the gibbs free energy //! of the element in its most stable state at 298.15 K and 1 bar. //#define GIBSSFE298_UNKNOWN 123456789. @@ -120,6 +188,37 @@ namespace Cantera { */ doublereal entropyElement298(int m) const; + //! Return the element constraint type + /*! + * Possible types include: + * + * CT_ELEM_TYPE_ABSPOS 0 + * CT_ELEM_TYPE_ELECTRONCHARGE 1 + * CT_ELEM_TYPE_CHARGENEUTRALITY 2 + * CT_ELEM_TYPE_LATTICERATIO 3 + * CT_ELEM_TYPE_KINETICFROZEN 4 + * CT_ELEM_TYPE_SURFACECONSTRAINT 5 + * CT_ELEM_TYPE_OTHERCONSTRAINT 6 + * + * The default is CT_ELEM_TYPE_ABSPOS + * + * @param m Element index + * + * @return Returns the element type + */ + int elementType(int m) const; + + //! Change the element type of the mth constraint + /*! + * Reassigns an element type + * + * @param m Element index + * @param elem_type New elem type to be assigned + * + * @return Returns the old element type + */ + int changeElementType(int m, int elem_type); + /// vector of element atomic weights const vector_fp& atomicWeights() const { return m_atomicWeights; } @@ -199,7 +298,7 @@ namespace Cantera { */ void addUniqueElement(const std::string& symbol, doublereal weight = -12345.0, int atomicNumber = 0, - doublereal entropy298 = ENTROPY298_UNKNOWN); + doublereal entropy298 = ENTROPY298_UNKNOWN, int elem_type = CT_ELEM_TYPE_ABSPOS); //! Add an element to the current set of elements in the current object. /*! @@ -258,7 +357,7 @@ namespace Cantera { * If this is true, then no elements may be added to the * object. */ - bool m_elementsFrozen; + bool m_elementsFrozen; /** * Vector of element atomic weights: @@ -285,6 +384,9 @@ namespace Cantera { */ vector_fp m_entropy298; + //! Vector of element types + vector_int m_elem_type; + /** * Number of Constituents Objects that use this object * diff --git a/Cantera/src/thermo/FixedChemPotSSTP.cpp b/Cantera/src/thermo/FixedChemPotSSTP.cpp index 0f845b1f6..42a99d96e 100644 --- a/Cantera/src/thermo/FixedChemPotSSTP.cpp +++ b/Cantera/src/thermo/FixedChemPotSSTP.cpp @@ -144,6 +144,8 @@ namespace Cantera { ss.addChild("h", sval); ss.addChild("s", "0.0"); saveSpeciesData(0, s); + delete s; + s = 0; } //==================================================================================================================== diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index 17339386b..2526b9ff9 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -2053,14 +2053,6 @@ namespace Cantera { * -------------- Utilities ------------------------------- */ - /** - * Return a reference to the species thermodynamic property - * manager. - * - * @todo This method will fail if no species thermo - * manager has been installed. - */ - SpeciesThermo& speciesThermo() { return *m_spthermo; } //! Initialization of a HMWSoln phase using an xml file /*! diff --git a/Cantera/src/thermo/IdealMolalSoln.h b/Cantera/src/thermo/IdealMolalSoln.h index 737fe349a..151eb3083 100644 --- a/Cantera/src/thermo/IdealMolalSoln.h +++ b/Cantera/src/thermo/IdealMolalSoln.h @@ -777,14 +777,6 @@ namespace Cantera { /* * -------------- Utilities ------------------------------- */ - - /*! - * Return a reference to the species thermodynamic property - * manager. @todo This method will fail if no species thermo - * manager has been installed. - */ - SpeciesThermo& speciesThermo() { return *m_spthermo; } - //! Initialization routine for an IdealMolalSoln phase. /*! diff --git a/Cantera/src/thermo/IdealSolidSolnPhase.cpp b/Cantera/src/thermo/IdealSolidSolnPhase.cpp index 6fd4107e6..8841712eb 100644 --- a/Cantera/src/thermo/IdealSolidSolnPhase.cpp +++ b/Cantera/src/thermo/IdealSolidSolnPhase.cpp @@ -117,7 +117,7 @@ namespace Cantera { IdealSolidSolnPhase *ii = new IdealSolidSolnPhase(*this); return (ThermoPhase*) ii; } -//==================================================================================================================== + //==================================================================================================================== /** * Equation of state flag. Returns the value cIdealGas, defined * in mix_defs.h. @@ -679,8 +679,7 @@ namespace Cantera { * property manager. They are polynomial functions of temperature. * @see SpeciesThermo */ - void IdealSolidSolnPhase:: - getPartialMolarEnthalpies(doublereal* hbar) const { + void IdealSolidSolnPhase::getPartialMolarEnthalpies(doublereal* hbar) const { const array_fp& _h = enthalpy_RT_ref(); doublereal rt = GasConstant * temperature(); scale(_h.begin(), _h.end(), hbar, rt); @@ -892,8 +891,7 @@ namespace Cantera { * units = m^3 / kmol */ void IdealSolidSolnPhase::getStandardVolumes(doublereal *vol) const { - copy(m_speciesMolarVolume.begin(), - m_speciesMolarVolume.end(), vol); + copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), vol); } @@ -1227,8 +1225,7 @@ namespace Cantera { "Unknown standardConc model: " + formStringa); } } else { - throw CanteraError(subname.c_str(), - "Unspecified standardConc model"); + throw CanteraError(subname.c_str(), "Unspecified standardConc model"); } /* diff --git a/Cantera/src/thermo/IdealSolidSolnPhase.h b/Cantera/src/thermo/IdealSolidSolnPhase.h index a9f87b29c..279463915 100644 --- a/Cantera/src/thermo/IdealSolidSolnPhase.h +++ b/Cantera/src/thermo/IdealSolidSolnPhase.h @@ -579,9 +579,9 @@ namespace Cantera { /// @name Partial Molar Properties of the Solution ----------------------------- //@{ - /** - * Returns an array of partial molar enthalpies for the species - * in the mixture. + + //! 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 * pure species enthalpies @@ -1057,8 +1057,9 @@ namespace Cantera { */ doublereal m_Pcurrent; + //! Vector of molar volumes for each species in the solution /** - * Species molar volume \f$ m^3 kmol^-1 \f$ + * Species molar volumes \f$ m^3 kmol^-1 \f$ */ array_fp m_speciesMolarVolume; diff --git a/Cantera/src/thermo/LatticePhase.cpp b/Cantera/src/thermo/LatticePhase.cpp index f86a26823..a0c900464 100644 --- a/Cantera/src/thermo/LatticePhase.cpp +++ b/Cantera/src/thermo/LatticePhase.cpp @@ -28,11 +28,20 @@ #include #include +using namespace std; + namespace Cantera { // Base Empty constructor LatticePhase::LatticePhase() : - m_tlast(0.0) + m_mm(0), + m_tmin(0.0), + m_tmax(0.0), + m_Pref(OneAtm), + m_Pcurrent(OneAtm), + m_tlast(0.0), + m_speciesMolarVolume(0), + m_site_density(0.0) { } @@ -41,7 +50,14 @@ namespace Cantera { * @param right Object to be copied */ LatticePhase::LatticePhase(const LatticePhase &right) : - m_tlast(0.0) + m_mm(0), + m_tmin(0.0), + m_tmax(0.0), + m_Pref(OneAtm), + m_Pcurrent(OneAtm), + m_tlast(0.0), + m_speciesMolarVolume(0), + m_site_density(0.0) { *this = operator=(right); } @@ -56,15 +72,16 @@ namespace Cantera { m_mm = right.m_mm; m_tmin = right.m_tmin; m_tmax = right.m_tmax; - m_p0 = right.m_p0; + m_Pref = right.m_Pref; + m_Pcurrent = right.m_Pcurrent; m_tlast = right.m_tlast; m_h0_RT = right.m_h0_RT; m_cp0_R = right.m_cp0_R; m_g0_RT = right.m_g0_RT; m_s0_R = right.m_s0_R; - m_press = right.m_press; m_vacancy = right.m_vacancy; - m_molar_density = right.m_molar_density; + m_speciesMolarVolume = right.m_speciesMolarVolume; + m_site_density = right.m_site_density; } return *this; } @@ -195,99 +212,176 @@ namespace Cantera { return GasConstant * (mean_X(&entropy_R_ref()[0]) - sum_xlogx()); } - + //==================================================================================================================== doublereal LatticePhase::gibbs_mole() const { return enthalpy_mole() - temperature() * entropy_mole(); } - + //==================================================================================================================== doublereal LatticePhase::cp_mole() const { return GasConstant * mean_X(&cp_R_ref()[0]); } - + //==================================================================================================================== doublereal LatticePhase::cv_mole() const { return cp_mole(); } - - - - void LatticePhase::setPressure(doublereal p) { - m_press = p; - setMolarDensity(m_molar_density); + //==================================================================================================================== + doublereal LatticePhase::calcDensity() { + setMolarDensity(m_site_density); + doublereal mw = meanMolecularWeight(); + doublereal dens = mw * m_site_density; + /* + * Calculate the molarVolume of the solution (m**3 kmol-1) + */ + // const doublereal * const dtmp = moleFractdivMMW(); + // doublereal invDens = dot(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), dtmp); + /* + * Set the density in the parent State object directly, + * by calling the State::setDensity() function. + */ + // doublereal dens = 1.0/invDens; + // State::setDensity(dens); + return dens; } - + //==================================================================================================================== + void LatticePhase::setPressure(doublereal p) { + m_Pcurrent = p; + calcDensity(); + } + //==================================================================================================================== + void LatticePhase::setMoleFractions(const doublereal * const x) { + State::setMoleFractions(x); + calcDensity(); + } + //==================================================================================================================== + void LatticePhase::setMoleFractions_NoNorm(const doublereal * const x) { + State::setMoleFractions(x); + calcDensity(); + } + //==================================================================================================================== + void LatticePhase::setMassFractions(const doublereal * const y) { + State::setMassFractions(y); + calcDensity(); + } + //==================================================================================================================== + void LatticePhase::setMassFractions_NoNorm(const doublereal * const y) { + State::setMassFractions_NoNorm(y); + calcDensity(); + } + //==================================================================================================================== + void LatticePhase::setConcentrations(const doublereal * const c) { + State::setConcentrations(c); + calcDensity(); + } + //==================================================================================================================== void LatticePhase::getActivityConcentrations(doublereal* c) const { getMoleFractions(c); } - + //==================================================================================================================== void LatticePhase::getActivityCoefficients(doublereal* ac) const { for (int k = 0; k < m_kk; k++) { ac[k] = 1.0; } } - + //==================================================================================================================== doublereal LatticePhase::standardConcentration(int k) const { return 1.0; } - + //==================================================================================================================== doublereal LatticePhase::logStandardConc(int k) const { return 0.0; } - + //==================================================================================================================== void LatticePhase::getChemPotentials(doublereal* mu) const { - doublereal vdp = ((pressure() - m_spthermo->refPressure())/ - molarDensity()); + doublereal delta_p = m_Pcurrent - m_Pref; doublereal xx; - doublereal rt = temperature() * GasConstant; + doublereal RT = temperature() * GasConstant; const array_fp& g_RT = gibbs_RT_ref(); for (int k = 0; k < m_kk; k++) { xx = fmaxx(SmallNumber, moleFraction(k)); - mu[k] = rt*(g_RT[k] + log(xx)) + vdp; + mu[k] = RT * (g_RT[k] + log(xx)) + + delta_p * m_speciesMolarVolume[k]; + } + + } + //==================================================================================================================== + void LatticePhase::getPartialMolarEnthalpies(doublereal* hbar) const { + const array_fp& _h = enthalpy_RT_ref(); + doublereal rt = GasConstant * temperature(); + scale(_h.begin(), _h.end(), hbar, rt); + } + //==================================================================================================================== + void LatticePhase::getPartialMolarEntropies(doublereal* sbar) const { + const array_fp& _s = entropy_R_ref(); + doublereal r = GasConstant; + doublereal xx; + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(SmallNumber, moleFraction(k)); + sbar[k] = r * (_s[k] - log(xx)); } } - + //==================================================================================================================== + void LatticePhase::getPartialMolarCp(doublereal* cpbar) const { + getCp_R(cpbar); + for (int k = 0; k < m_kk; k++) { + cpbar[k] *= GasConstant; + } + } + //==================================================================================================================== void LatticePhase::getPartialMolarVolumes(doublereal* vbar) const { getStandardVolumes(vbar); } - + //==================================================================================================================== void LatticePhase::getStandardChemPotentials(doublereal* mu0) const { const array_fp& gibbsrt = gibbs_RT_ref(); scale(gibbsrt.begin(), gibbsrt.end(), mu0, _RT()); } - + //==================================================================================================================== void LatticePhase::getPureGibbs(doublereal* gpure) const { const array_fp& gibbsrt = gibbs_RT_ref(); - scale(gibbsrt.begin(), gibbsrt.end(), gpure, _RT()); - } - - void LatticePhase::getEnthalpy_RT(doublereal* hrt) const { - const array_fp& _h = enthalpy_RT_ref(); - std::copy(_h.begin(), _h.end(), hrt); - doublereal tmp = (pressure() - m_p0) / (molarDensity() * GasConstant * temperature()); + doublereal delta_p = (m_Pcurrent - m_Pref); + double RT = GasConstant * temperature(); for (int k = 0; k < m_kk; k++) { - hrt[k] += tmp; + gpure[k] = RT * gibbsrt[k] + delta_p * m_speciesMolarVolume[k]; } } - + //==================================================================================================================== + void LatticePhase::getEnthalpy_RT(doublereal* hrt) const { + const array_fp& _h = enthalpy_RT_ref(); + doublereal delta_prt = ((m_Pcurrent - m_Pref) / (GasConstant * temperature())); + for (int k = 0; k < m_kk; k++) { + hrt[k] = _h[k] + delta_prt * m_speciesMolarVolume[k]; + } + } + //==================================================================================================================== void LatticePhase::getEntropy_R(doublereal* sr) const { const array_fp& _s = entropy_R_ref(); std::copy(_s.begin(), _s.end(), sr); } - + //==================================================================================================================== void LatticePhase::getGibbs_RT(doublereal* grt) const { const array_fp& gibbsrt = gibbs_RT_ref(); - std::copy(gibbsrt.begin(), gibbsrt.end(), grt); + doublereal RT = _RT(); + doublereal delta_prt = (m_Pcurrent - m_Pref)/ RT; + for (int k = 0; k < m_kk; k++) { + grt[k] = gibbsrt[k] + delta_prt * m_speciesMolarVolume[k]; + } } + //==================================================================================================================== + void LatticePhase::getGibbs_ref(doublereal *g) const { + getGibbs_RT_ref(g); + for (int k = 0; k < m_kk; k++) { + g[k] *= GasConstant * temperature(); + } + } + //=================================================================================================================== void LatticePhase::getCp_R(doublereal* cpr) const { const array_fp& _cpr = cp_R_ref(); std::copy(_cpr.begin(), _cpr.end(), cpr); } - + //=================================================================================================================== void LatticePhase::getStandardVolumes(doublereal* vbar) const { - doublereal vv = 1.0/m_molar_density; - for (int k = 0; k < m_kk; k++) { - vbar[k] = vv; - } + copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), vbar); } //======================================================================================================= // Returns the vector of nondimensional Enthalpies of the reference state at the current temperature @@ -310,6 +404,13 @@ namespace Cantera { _updateThermo(); return m_g0_RT; } + //==================================================================================================================== + void LatticePhase::getGibbs_RT_ref(doublereal *grt) const { + _updateThermo(); + for (int k = 0; k < m_kk; k++) { + grt[k] = m_g0_RT[k]; + } + } //======================================================================================================= // Returns a reference to the dimensionless reference state Entropy vector. /* @@ -330,7 +431,21 @@ namespace Cantera { _updateThermo(); return m_cp0_R; } - //======================================================================================================= + //==================================================================================================================== + // Initialize the ThermoPhase object after all species have been set up + /* + * @internal Initialize. + * + * This method performs any initialization required after all + * species have been added. For example, it is used to + * resize internal work arrays that must have an entry for + * each species. + * This method is called from ThermoPhase::initThermoXML(), + * which is called from importPhase(), + * just prior to returning from the function, importPhase(). + * + * @see importCTML.cpp + */ void LatticePhase::initThermo() { m_kk = nSpecies(); m_mm = nElements(); @@ -338,14 +453,62 @@ namespace Cantera { doublereal tmax = m_spthermo->maxTemp(); if (tmin > 0.0) m_tmin = tmin; if (tmax > 0.0) m_tmax = tmax; - m_p0 = refPressure(); + m_Pref = refPressure(); int leng = m_kk; m_h0_RT.resize(leng); m_g0_RT.resize(leng); m_cp0_R.resize(leng); m_s0_R.resize(leng); - setMolarDensity(m_molar_density); + m_speciesMolarVolume.resize(leng, 0.0); + + ThermoPhase::initThermo(); + } + //==================================================================================================================== + void LatticePhase::initThermoXML(XML_Node& phaseNode, std::string id) { + std::string subname = "LatticePhase::initThermoXML"; + /* + * Check on the thermo field. Must have: + * + */ + if (phaseNode.hasChild("thermo")) { + XML_Node& thNode = phaseNode.child("thermo"); + std::string mStringa = thNode.attrib("model"); + std::string mString = lowercase(mStringa); + if (mString != "lattice") { + throw CanteraError(subname.c_str(), + "Unknown thermo model: " + mStringa); + } + } else { + throw CanteraError(subname.c_str(), + "Unspecified thermo model"); + } + /* + * Now go get the molar volumes. use the default if not found + */ + XML_Node& speciesList = phaseNode.child("speciesArray"); + XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"], &phaseNode.root()); + const std::vector &sss = speciesNames(); + + for (int k = 0; k < m_kk; k++) { + m_speciesMolarVolume[k] = m_site_density; + XML_Node* s = speciesDB->findByAttr("name", sss[k]); + if (!s) { + throw CanteraError(" LatticePhase::initThermoXML", "database problems"); + } + XML_Node *ss = s->findByName("standardState"); + if (ss) { + if (ss->findByName("molarVolume")) { + m_speciesMolarVolume[k] = getFloat(*ss, "molarVolume", "toSI"); + } + } + } + + /* + * Call the base initThermo, which handles setting the initial + * state. + */ + ThermoPhase::initThermoXML(phaseNode, id); } //===================================================================================================== // Update the species reference state thermodynamic functions @@ -355,10 +518,6 @@ namespace Cantera { */ void LatticePhase::_updateThermo() const { doublereal tnow = temperature(); - if (fabs(molarDensity() - m_molar_density)/m_molar_density > 0.0001) { - throw CanteraError("_updateThermo","molar density changed from " - +fp2str(m_molar_density)+" to "+fp2str(molarDensity())); - } if (m_tlast != tnow) { m_spthermo->update(tnow, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]); m_tlast = tnow; @@ -370,8 +529,8 @@ namespace Cantera { } //===================================================================================================== void LatticePhase::setParameters(int n, doublereal* const c) { - m_molar_density = c[0]; - setMolarDensity(m_molar_density); + m_site_density = c[0]; + setMolarDensity(m_site_density); } //===================================================================================================== void LatticePhase::getParameters(int &n, doublereal * const c) const { @@ -382,7 +541,7 @@ namespace Cantera { //===================================================================================================== void LatticePhase::setParametersFromXML(const XML_Node& eosdata) { eosdata._require("model", "Lattice"); - m_molar_density = getFloat(eosdata, "site_density", "toSI"); + m_site_density = getFloat(eosdata, "site_density", "toSI"); m_vacancy = getChildValue(eosdata, "vacancy_species"); } //===================================================================================================== diff --git a/Cantera/src/thermo/LatticePhase.h b/Cantera/src/thermo/LatticePhase.h index b2622955a..eabde4c25 100644 --- a/Cantera/src/thermo/LatticePhase.h +++ b/Cantera/src/thermo/LatticePhase.h @@ -453,7 +453,7 @@ namespace Cantera { * independent value of the pressure. */ virtual doublereal pressure() const { - return m_press; + return m_Pcurrent; } //! Set the internally storred pressure (Pa) at constant @@ -465,7 +465,67 @@ namespace Cantera { * @param p Input Pressure (Pa) */ virtual void setPressure(doublereal p); - + + //! Calculate the density of the mixture using the partial + //! molar volumes and mole fractions as input + /*! + * The formula for this is + * + * \f[ + * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} + * \f] + * + * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are + * the molecular weights, and \f$V_k\f$ are the pure species + * molar volumes. + * + * Note, the basis behind this formula is that in an ideal + * solution the partial molar volumes are equal to the pure + * species molar volumes. We have additionally specified + * in this class that the pure species molar volumes are + * independent of temperature and pressure. + * + * NOTE: This is a non-virtual function, which is not a + * member of the ThermoPhase base class. + */ + doublereal calcDensity(); + + //! Set the mole fractions + /*! + * @param x Input vector of mole fractions. + * Length: m_kk. + */ + virtual void setMoleFractions(const doublereal * const x); + + //! Set the mole fractions, but don't normalize them to one. + /*! + * @param x Input vector of mole fractions. + * Length: m_kk. + */ + virtual void setMoleFractions_NoNorm(const doublereal * const x); + + //! Set the mass fractions, and normalize them to one. + /*! + * @param y Input vector of mass fractions. + * Length: m_kk. + */ + virtual void setMassFractions(const doublereal * const y); + + //! Set the mass fractions, but don't normalize them to one + /*! + * @param y Input vector of mass fractions. + * Length: m_kk. + */ + virtual void setMassFractions_NoNorm(const doublereal * const y); + + //! Set the concentration, + /*! + * @param c Input vector of concentrations. + * Length: m_kk. + */ + virtual void setConcentrations(const doublereal * const c); + + //@} /// @name Activities, Standard States, and Activity Concentrations /** @@ -547,6 +607,71 @@ namespace Cantera { */ virtual void getChemPotentials(doublereal* mu) const; + + //@} + /// @name Partial Molar Properties of the Solution ----------------------------- + //@{ + + + /** + * 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 + * pure species enthalpies + * \f[ + * \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k + * \f] + * The reference-state pure-species enthalpies, \f$ \hat h^{ref}_k(T) \f$, + * at the reference pressure,\f$ P_{ref} \f$, + * are computed by the species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + * + * @param hbar Output vector containing partial molar enthalpies. + * Length: m_kk. + */ + virtual void getPartialMolarEnthalpies(doublereal* hbar) const; + + /** + * Returns an array of partial molar entropies of the species in the + * solution. Units: J/kmol/K. + * For this phase, the partial molar entropies are equal to the + * pure species entropies plus the ideal solution contribution. + * \f[ + * \bar s_k(T,P) = \hat s^0_k(T) - R log(X_k) + * \f] + * The reference-state pure-species entropies,\f$ \hat s^{ref}_k(T) \f$, + * at the reference pressure, \f$ P_{ref} \f$, are computed by the + * species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + * + * @param sbar Output vector containing partial molar entropies. + * Length: m_kk. + */ + virtual void getPartialMolarEntropies(doublereal* sbar) const; + + /** + * Returns an array of partial molar Heat Capacities at constant + * pressure of the species in the + * solution. Units: J/kmol/K. + * For this phase, the partial molar heat capacities are equal + * to the standard state heat capacities. + * + * @param cpbar Output vector of partial heat capacities. Length: m_kk. + */ + virtual void getPartialMolarCp(doublereal* cpbar) const; + + //! Return an array of partial molar volumes for the + //! species in the mixture. Units: m^3/kmol. + /*! + * @param vbar Output vector of speciar partial molar volumes. + * Length = m_kk. units are m^3/kmol. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + + //! Get the array of chemical potentials at unit activity for the //! species standard states at the current T and P of the solution. /*! @@ -568,14 +693,7 @@ namespace Cantera { */ virtual void getPureGibbs(doublereal* gpure) const; - //! Return an array of partial molar volumes for the - //! species in the mixture. Units: m^3/kmol. - /*! - * @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 //@{ @@ -702,6 +820,25 @@ namespace Cantera { */ const array_fp& gibbs_RT_ref() const; + //! Returns the vector of nondimensional + //! Gibbs Free Energies of the reference state at the current temperature + //! of the solution and the reference pressure for the species. + /*! + * @param grt Output vector containing the nondimensional reference state + * Gibbs Free energies. Length: m_kk. + */ + virtual void getGibbs_RT_ref(doublereal *grt) const; + + //! Returns the vector of the gibbs function of the reference state at the current temperature + //! of the solution and the reference pressure for the species. + /*! + * units = J/kmol + * + * @param g Output vector containing the reference state + * Gibbs Free energies. Length: m_kk. Units: J/kmol. + */ + virtual void getGibbs_ref(doublereal *g) const; + //! Returns a reference to the dimensionless reference state Entropy vector. /*! * This function is part of the layer that checks/recalculates the reference @@ -736,6 +873,33 @@ namespace Cantera { */ virtual void initThermo(); + + //! Import and initialize a ThermoPhase object using an XML tree. + /*! + * Here we read extra information about the XML description + * of a phase. Regular information about elements and species + * and their reference state thermodynamic information + * have already been read at this point. + * For example, we do not need to call this function for + * ideal gas equations of state. + * This function is called from importPhase() + * after the elements and the + * species are initialized with default ideal solution + * level data. + * + * @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. + */ + virtual void initThermoXML(XML_Node& phaseNode, std::string id); + //! Set the equation of state parameters from the argument list /*! * @internal @@ -792,6 +956,7 @@ namespace Cantera { protected: + //! Number of elements int m_mm; @@ -810,7 +975,17 @@ namespace Cantera { doublereal m_tmax; //! Reference state pressure - doublereal m_p0; + doublereal m_Pref; + + + //! The current pressure + /*! + * Since the density isn't a function of pressure, but only of the + * mole fractions, we need to independently specify the pressure. + * The density variable which is inherited as part of the State class, + * m_dens, is always kept current whenever T, P, or X[] change. + */ + doublereal m_Pcurrent; //! Current value of the temperature (Kelvin) mutable doublereal m_tlast; @@ -827,8 +1002,6 @@ namespace Cantera { //! Temporary storage for the reference state entropies at the current temperature mutable array_fp m_s0_R; - //! Current value of the pressure (Pa) - doublereal m_press; //! String name for the species which represents a vacency //! in the lattice @@ -837,13 +1010,21 @@ namespace Cantera { */ std::string m_vacancy; - //! Molar density of the lattice solid + //! Vector of molar volumes for each species in the solution + /** + * Species molar volumes \f$ m^3 kmol^-1 \f$ + */ + array_fp m_speciesMolarVolume; + + //! Site Density of the lattice solid /*! - * Currently, this does not change as a function of T, P or composition + * Currently, this is imposed as a function of T, P or composition * * units are kmol m-3 */ - doublereal m_molar_density; + doublereal m_site_density; + + // doublereal m_molar_lattice_volume; private: diff --git a/Cantera/src/thermo/LatticeSolidPhase.cpp b/Cantera/src/thermo/LatticeSolidPhase.cpp index 7e1aac142..6b0d04385 100644 --- a/Cantera/src/thermo/LatticeSolidPhase.cpp +++ b/Cantera/src/thermo/LatticeSolidPhase.cpp @@ -23,8 +23,16 @@ #include "LatticePhase.h" #include "SpeciesThermo.h" #include "ThermoFactory.h" +#include "SpeciesThermoFactory.h" +#include "GeneralSpeciesThermo.h" #include +#ifndef MIN +# define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) +#endif +#ifndef MAX +# define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif using namespace std; //====================================================================================================================== @@ -39,7 +47,9 @@ namespace Cantera { m_molar_density(0.0), m_nlattice(0), m_lattice(0), - m_x(0) + m_x(0), + theta_(0), + tmpV_(0) { } //==================================================================================================================== @@ -54,7 +64,9 @@ namespace Cantera { m_molar_density(0.0), m_nlattice(0), m_lattice(0), - m_x(0) + m_x(0), + theta_(0), + tmpV_(0) { *this = operator=(right); } @@ -74,12 +86,15 @@ namespace Cantera { m_nlattice = right.m_nlattice; deepStdVectorPointerCopy(right.m_lattice, m_lattice); m_x = right.m_x; + theta_ = right.theta_; + tmpV_ = right.tmpV_; } return *this; } //==================================================================================================================== // Destructor LatticeSolidPhase::~LatticeSolidPhase() { + // We own the sublattices. So we have to delete the sublattices for (int n = 0; n < m_nlattice; n++) { delete m_lattice[n]; m_lattice[n] = 0; @@ -98,59 +113,121 @@ namespace Cantera { LatticeSolidPhase *igp = new LatticeSolidPhase(*this); return (ThermoPhase *) igp; } + + //==================================================================================================================== + // Minimum temperature for which the thermodynamic data for the species + // or phase are valid. + /* + * If no argument is supplied, the + * value returned will be the lowest temperature at which the + * data for \e all species are valid. Otherwise, the value + * will be only for species \a k. This function is a wrapper + * that calls the species thermo minTemp function. + * + * @param k index of the species. Default is -1, which will return the max of the min value + * over all species. + */ + doublereal LatticeSolidPhase::minTemp(int k) const { + if (k >= 0) { + for (int n = 0; n < m_nlattice; n++) { + if (lkstart_[n+1] < k) { + double ml = (m_lattice[n])->minTemp(k-lkstart_[n]); + return ml; + } + } + } + doublereal mm = 1.0E300; + for (int n = 0; n < m_nlattice; n++) { + double ml = (m_lattice[n])->minTemp(-1); + mm = MIN(mm, ml); + } + return mm; + } + //==================================================================================================================== + // Maximum temperature for which the thermodynamic data for the species + // or phase are valid. + /* + * If no argument is supplied, the + * value returned will be the lowest temperature at which the + * data for \e all species are valid. Otherwise, the value + * will be only for species \a k. This function is a wrapper + * that calls the species thermo minTemp function. + * + * @param k index of the species. Default is -1, which will return the max of the min value + * over all species. + */ + doublereal LatticeSolidPhase::maxTemp(int k) const { + if (k >= 0) { + for (int n = 0; n < m_nlattice; n++) { + if (lkstart_[n+1] < k) { + double ml = (m_lattice[n])->maxTemp(k - lkstart_[n]); + return ml; + } + } + } + doublereal mm = -1.0E300; + for (int n = 0; n < m_nlattice; n++) { + double ml = (m_lattice[n])->maxTemp(-1); + mm = MAX(mm, ml); + } + return mm; + } + //==================================================================================================================== + /* + * Returns the reference pressure in Pa. This function is a wrapper + * that calls the species thermo refPressure function. + */ + doublereal LatticeSolidPhase::refPressure() const { + return m_lattice[0]->refPressure(); + } //==================================================================================================================== doublereal LatticeSolidPhase:: enthalpy_mole() const { _updateThermo(); - doublereal ndens, sum = 0.0; + doublereal sum = 0.0; int n; for (n = 0; n < m_nlattice; n++) { - ndens = m_lattice[n]->molarDensity(); - sum += ndens * m_lattice[n]->enthalpy_mole(); + sum += theta_[n] * m_lattice[n]->enthalpy_mole(); } - return sum/molarDensity(); + return sum; } //==================================================================================================================== doublereal LatticeSolidPhase::intEnergy_mole() const { _updateThermo(); - doublereal ndens, sum = 0.0; + doublereal sum = 0.0; int n; for (n = 0; n < m_nlattice; n++) { - ndens = m_lattice[n]->molarDensity(); - sum += ndens * m_lattice[n]->intEnergy_mole(); + sum += theta_[n] * m_lattice[n]->intEnergy_mole(); } - return sum/molarDensity(); + return sum; } //==================================================================================================================== doublereal LatticeSolidPhase::entropy_mole() const { _updateThermo(); - doublereal ndens, sum = 0.0; + doublereal sum = 0.0; int n; for (n = 0; n < m_nlattice; n++) { - ndens = m_lattice[n]->molarDensity(); - sum += ndens * m_lattice[n]->entropy_mole(); + sum += theta_[n] * m_lattice[n]->entropy_mole(); } - return sum/molarDensity(); + return sum; } //==================================================================================================================== doublereal LatticeSolidPhase::gibbs_mole() const { _updateThermo(); - doublereal ndens, sum = 0.0; + doublereal sum = 0.0; for (int n = 0; n < m_nlattice; n++) { - ndens = m_lattice[n]->molarDensity(); - sum += ndens * m_lattice[n]->gibbs_mole(); + sum += theta_[n] * m_lattice[n]->gibbs_mole(); } - return sum/molarDensity(); + return sum; } //==================================================================================================================== doublereal LatticeSolidPhase::cp_mole() const { _updateThermo(); doublereal sum = 0.0; for (int n = 0; n < m_nlattice; n++) { - doublereal ndens = m_lattice[n]->molarDensity(); - sum += ndens * m_lattice[n]->cp_mole(); + sum += theta_[n] * m_lattice[n]->cp_mole(); } - return sum/molarDensity(); + return sum; } //==================================================================================================================== void LatticeSolidPhase::getActivityConcentrations(doublereal* c) const { @@ -176,12 +253,49 @@ namespace Cantera { return 0.0; } + //==================================================================================================================== + // Set the pressure at constant temperature. Units: Pa. + /* + * + * @param p Pressure (units - Pa) + */ + void LatticeSolidPhase::setPressure(doublereal p) { + m_press = p; + for (int n = 0; n < m_nlattice; n++) { + m_lattice[n]->setPressure(m_press); + } + calcDensity(); + } + //==================================================================================================================== + // Calculate the density of the solid mixture + /* + * The formula for this is + * + * \f[ + * \rho = \sum_n{ \rho_n \theta_n } + * \f] + * + * where \f$ \rho_n \f$ is the density of the nth sublattice + * + * Note this is a nonvirtual function. + */ + doublereal LatticeSolidPhase::calcDensity() { + double sum = 0.0; + for (int n = 0; n < m_nlattice; n++) { + sum += theta_[n] * m_lattice[n]->density(); + } + State::setDensity(sum); + return sum; + } //==================================================================================================================== // Set the mole fractions to the specified values, and then // normalize them so that they sum to 1.0 for each of the subphases /* + * On input, the mole fraction vector is assumed to sum to one for each of the sublattices. The sublattices + * are updated with this mole fraction vector. The mole fractions are also storred within this object, after + * they are normalized to one by dividing by the number of sublattices. * - * @param x Input vector of mole fractions. There is no restriction + * @param x Input vector of mole fractions. There is no restriction * on the sum of the mole fraction vector. Internally, * this object will pass portions of this vector to the sublattices which assume that the portions * individually sum to one. @@ -189,19 +303,16 @@ namespace Cantera { */ void LatticeSolidPhase::setMoleFractions(const doublereal* const x) { int nsp, strt = 0; - doublereal sum = 0.0; for (int n = 0; n < m_nlattice; n++) { nsp = m_lattice[n]->nSpecies(); - m_lattice[n]->setMoleFractions(x+strt); - for (int k = 0; k < nsp; k++) { - sum += x[strt + k]; - } + m_lattice[n]->setMoleFractions(x + strt); strt += nsp; } for (int k = 0; k < strt; k++) { - m_x[k] = x[k] / sum; + m_x[k] = x[k] / m_nlattice; } State::setMoleFractions(DATA_PTR(m_x)); + calcDensity(); } //==================================================================================================================== // Get the species mole fraction vector. @@ -213,6 +324,7 @@ namespace Cantera { */ void LatticeSolidPhase::getMoleFractions(doublereal* const x) const { int nsp, strt = 0; + // the ifdef block should be the way we calculate this.!!!!! State::getMoleFractions(x); doublereal sum; for (int n = 0; n < m_nlattice; n++) { @@ -241,28 +353,103 @@ namespace Cantera { } } //==================================================================================================================== + // 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. + * + * This returns the underlying lattice chemical potentials + * + * @param mu Output vector of species chemical + * potentials. Length: m_kk. Units: J/kmol + */ void LatticeSolidPhase::getChemPotentials(doublereal* mu) const { _updateThermo(); int strt = 0; for (int n = 0; n < m_nlattice; n++) { - doublereal dratio = m_lattice[n]->molarDensity()/molarDensity(); + int nlsp = m_lattice[n]->nSpecies(); m_lattice[n]->getChemPotentials(mu+strt); - scale(mu + strt, mu + strt + m_lattice[n]->nSpecies(), mu + strt, dratio); - strt += m_lattice[n]->nSpecies(); + strt += nlsp; } } //==================================================================================================================== + void LatticeSolidPhase::getPartialMolarEnthalpies(doublereal* hbar) const { + _updateThermo(); + int strt = 0; + for (int n = 0; n < m_nlattice; n++) { + int nlsp = m_lattice[n]->nSpecies(); + m_lattice[n]->getPartialMolarEnthalpies(hbar + strt); + strt += nlsp; + } + } + //==================================================================================================================== + void LatticeSolidPhase::getPartialMolarEntropies(doublereal* sbar) const { + _updateThermo(); + int strt = 0; + for (int n = 0; n < m_nlattice; n++) { + int nlsp = m_lattice[n]->nSpecies(); + m_lattice[n]->getPartialMolarEntropies(sbar + strt); + strt += nlsp; + } + } + //==================================================================================================================== + void LatticeSolidPhase::getPartialMolarCp(doublereal* cpbar) const { + _updateThermo(); + int strt = 0; + for (int n = 0; n < m_nlattice; n++) { + int nlsp = m_lattice[n]->nSpecies(); + m_lattice[n]->getPartialMolarCp(cpbar + strt); + strt += nlsp; + } + } + //==================================================================================================================== + void LatticeSolidPhase::getPartialMolarVolumes(doublereal* vbar) const { + _updateThermo(); + int strt = 0; + for (int n = 0; n < m_nlattice; n++) { + int nlsp = m_lattice[n]->nSpecies(); + m_lattice[n]->getPartialMolarVolumes(vbar + strt); + strt += nlsp; + } + } + //==================================================================================================================== + // Get the array of standard state chemical potentials at unit activity for the species + // at their standard states at the current T and P of the solution. + /* + * These are the standard state chemical potentials \f$ \mu^0_k(T,P) + * \f$. The values are evaluated at the current + * temperature and pressure of the solution. + * + * This returns the underlying lattice standard chemical potentials, as the units are kmol-1 of + * the sublattice species. + * + * @param mu0 Output vector of chemical potentials. + * Length: m_kk. Units: J/kmol + */ void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const { _updateThermo(); int strt = 0; for (int n = 0; n < m_nlattice; n++) { - doublereal dratio = m_lattice[n]->molarDensity()/molarDensity(); m_lattice[n]->getStandardChemPotentials(mu0+strt); - scale(mu0 + strt, mu0 + strt + m_lattice[n]->nSpecies(), mu0 + strt, dratio); strt += m_lattice[n]->nSpecies(); } } //==================================================================================================================== + void LatticeSolidPhase::getGibbs_RT_ref(doublereal *grt) const { + _updateThermo(); + for (int n = 0; n < m_nlattice; n++) { + m_lattice[n]->getGibbs_RT_ref(grt + lkstart_[n]); + } + } + //==================================================================================================================== + void LatticeSolidPhase::getGibbs_ref(doublereal *g) const { + getGibbs_RT_ref(g); + for (int k = 0; k < m_kk; k++) { + g[k] *= GasConstant * temperature(); + } + } + //==================================================================================================================== // Add in species from Slave phases /* * This hook is used for cSS_CONVENTION_SLAVE phases @@ -272,84 +459,127 @@ namespace Cantera { void LatticeSolidPhase::installSlavePhases(Cantera::XML_Node* phaseNode) { int m, k; + int kk = 0; + int kstart = 0; + SpeciesThermoFactory* spFactory = SpeciesThermoFactory::factory(); + SpeciesThermo * spthermo_ptr = new GeneralSpeciesThermo(); + setSpeciesThermo(spthermo_ptr); + m_speciesData.clear(); + + XML_Node& eosdata = phaseNode->child("thermo"); + XML_Node& la = eosdata.child("LatticeArray"); + std::vector lattices; + la.getChildren("phase",lattices); for (int n = 0; n < m_nlattice; n++) { LatticePhase *lp = m_lattice[n]; + XML_Node* phaseNode_ptr = lattices[n]; int nsp = lp->nSpecies(); vector constArr(lp->nElements()); + const vector_fp& aws = lp->atomicWeights(); + for (int es = 0; es < lp->nElements(); es++) { + string esName = lp->elementName(es); + double wt = aws[es]; + int an = lp->atomicNumber(es); + int e298 = lp->entropyElement298(es); + int et = lp->elementType(es); + addUniqueElementAfterFreeze(esName, wt, an, e298, et); + } + const std::vector & spNode = lp->speciesData(); + kstart = kk; + + for (k = 0; k < nsp; k++) { std::string sname = lp->speciesName(k); std::map comp; lp->getAtoms(k, DATA_PTR(constArr)); + int nel = nElements(); + vector_fp ecomp(nel, 0.0); for (m = 0; m < lp->nElements(); m++) { if (constArr[m] != 0.0) { - std::string ename = lp->elementName(m); - comp[ename] = constArr[m]; - } - } - int nel = nElements(); - vector_fp ecomp(nel, 0.0); - for (m = 0; m < nel; m++) { - double anum = comp[elementName(m)]; - if (anum != 0.0) { - ecomp[m] = anum; + std::string oldEname = lp->elementName(m); + int newIndex = elementIndex(oldEname); + if (newIndex < 0) { + throw CanteraError("LatticeSolidPhase::installSlavePhases", "confused"); + } + ecomp[newIndex] = constArr[m]; } } double chrg = lp->charge(k); double sz = lp->size(k); addUniqueSpecies(sname, &ecomp[0], chrg, sz); + spFactory->installThermoForSpecies(kk, *(spNode[k]), this, *m_spthermo, phaseNode_ptr); + + m_speciesData.push_back(new XML_Node(*(spNode[k]))); + kk++; + } + /* + * Add in the lattice stoichiometry constraint + */ + if (n > 0) { + string econ = "LC_"; + econ += int2str(n); + econ += "_" + id(); + int m = addUniqueElementAfterFreeze(econ, 0.0, 0, 0.0, CT_ELEM_TYPE_LATTICERATIO); + m_mm = nElements(); + LatticePhase *lp0 = m_lattice[0]; + int nsp0 = lp0->nSpecies(); + for (k = 0; k < nsp0; k++) { + m_speciesComp[k * m_mm + m] = -theta_[0]; + } + for (k = 0; k < nsp; k++) { + int ks = kstart + k; + m_speciesComp[ks * m_mm + m] = theta_[n]; + } } } } + //==================================================================================================================== + // Initialize the ThermoPhase object after all species have been set up + /* + * @internal Initialize. + * + * This method is provided to allow subclasses to perform any initialization required after all + * species have been added. For example, it might be used to + * resize internal work arrays that must have an entry for + * each species. The base class implementation does nothing, + * and subclasses that do not require initialization do not + * need to overload this method. When importing a CTML phase + * description, this method is called from ThermoPhase::initThermoXML(), + * which is called from importPhase(), just prior to returning from function importPhase(). + * + * @see importCTML.cpp + */ void LatticeSolidPhase::initThermo() { m_kk = nSpecies(); m_mm = nElements(); - m_x.resize(m_kk); + initLengths(); int nsp, k, loc = 0; - doublereal ndens; - m_molar_density = 0.0; for (int n = 0; n < m_nlattice; n++) { nsp = m_lattice[n]->nSpecies(); - ndens = m_lattice[n]->molarDensity(); + lkstart_[n] = loc; + nspLattice_[n] = nsp; for (k = 0; k < nsp; k++) { - m_x[loc] = ndens * m_lattice[n]->moleFraction(k); + m_x[loc] =m_lattice[n]->moleFraction(k) / (double) m_nlattice; loc++; } - m_molar_density += ndens; + lkstart_[n+1] = loc; } setMoleFractions(DATA_PTR(m_x)); - - // const vector& spnames = speciesNames(); - // int n, k, kl, namesize; - // int nl = m_sitedens.size(); - // string s; - // m_lattice.resize(m_kk,-1); - // vector_fp conc(m_kk, 0.0); - - // compositionMap xx; - // for (n = 0; n < nl; n++) { - // for (k = 0; k < m_kk; k++) { - // xx[speciesName(k)] = -1.0; - // } - // parseCompString(m_sp[n], xx); - // for (k = 0; k < m_kk; k++) { - // if (xx[speciesName(k)] != -1.0) { - // conc[k] = m_sitedens[n]*xx[speciesName(k)]; - // m_lattice[k] = n; - // } - // } - - // } - // for (k = 0; k < m_kk; k++) { - // if (m_lattice[k] == -1) { - // throw CanteraError("LatticeSolidPhase::" - // "setParametersFromXML","Species "+speciesName(k) - // +" not a member of any lattice."); - // } - // } - // setMoleFractions(DATA_PTR(conc)); + ThermoPhase::initThermo(); + } + //==================================================================================================================== + // Initialize vectors that depend on the number of species and sublattices + /* + * + */ + void LatticeSolidPhase::initLengths() { + theta_.resize(m_nlattice,0); + nspLattice_.resize(m_nlattice); + lkstart_.resize(m_nlattice+1); + m_x.resize(m_kk, 0.0); + tmpV_.resize(m_kk, 0.0); } - //==================================================================================================================== void LatticeSolidPhase::_updateThermo() const { doublereal tnow = temperature(); @@ -386,6 +616,14 @@ namespace Cantera { setMoleFractions(DATA_PTR(m_x)); } //==================================================================================================================== + + + //==================================================================================================================== + // Set the parameters from the XML file + /*! + * Currently, this is the spot that we read in all of the sublattice phases. + * The SetParametersFromXML() call is carried out at + */ void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata) { eosdata._require("model","LatticeSolid"); XML_Node& la = eosdata.child("LatticeArray"); @@ -398,9 +636,78 @@ namespace Cantera { XML_Node& i = *lattices[n]; m_lattice.push_back((LatticePhase*)newPhase(i)); } + std::vector pnam; + std::vector pval; + XML_Node& ls = eosdata.child("LatticeStoichiometry"); + int np = getPairs(ls, pnam, pval); + theta_.resize(nl); + for (int i = 0; i < np; i++) { + double val = fpValueCheck(pval[i]); + bool found = false; + for (int j = 0; j < nl; j++) { + ThermoPhase &tp = *(m_lattice[j]); + string idj = tp.id(); + if (idj == pnam[i]) { + theta_[j] = val; + found = true; + break; + } + } + if (!found) { + throw CanteraError("", "not found"); + } + } + + } + //==================================================================================================================== + // Return a changeable reference to the calculation manager + // for species reference-state thermodynamic properties + /* + * + * @param k Speices id. The default is -1, meaning return the default + * + * @internal + */ + SpeciesThermo& LatticeSolidPhase::speciesThermo(int k) { + return *m_spthermo; + /* + int kk; + if (k >= 0) { + for (int n = 0; n < m_nlattice; n++) { + if (lkstart_[n+1] < k) { + kk = k - lkstart_[n]; + return m_lattice[n]->speciesThermo(kk); + } + } + } + return m_lattice[0]->speciesThermo(-1); + */ } //==================================================================================================================== +#ifdef H298MODIFY_CAPABILITY + + //! Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1) + /*! + * The 298K heat of formation is defined as the enthalpy change to create the standard state + * of the species from its constituent elements in their standard states at 298 K and 1 bar. + * + * @param k Species k + * @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar + */ + void LatticeSolidPhase::modifyOneHf298SS(const int k, const doublereal Hf298New) { + for (int n = 0; n < m_nlattice; n++) { + if (lkstart_[n+1] < k) { + int kk = k-lkstart_[n]; + SpeciesThermo& l_spthermo = m_lattice[n]->speciesThermo(); + l_spthermo.modifyOneHf298(kk, Hf298New); + } + } + m_tlast += 0.0001234; + _updateThermo(); + } +#endif + //==================================================================================================================== doublereal LatticeSolidPhase::err(std::string msg) const { throw CanteraError("LatticeSolidPhase","Unimplemented " + msg); diff --git a/Cantera/src/thermo/LatticeSolidPhase.h b/Cantera/src/thermo/LatticeSolidPhase.h index f7acfd2e0..961c42316 100644 --- a/Cantera/src/thermo/LatticeSolidPhase.h +++ b/Cantera/src/thermo/LatticeSolidPhase.h @@ -34,17 +34,49 @@ namespace Cantera { - //! A phase that is comprised of an additive combination of other lattice phases + //! A phase that is comprised of a fixed additive combination of other lattice phases /*! - * This is the main way Cantera describes semiconductors and other solid phases. - * This %ThermoPhase object calculates its properties as a sum over other LatticePhase objects. Each of the %LatticePhase - * objects is a ThermoPhase object by itself. + * This is the main way %Cantera describes semiconductors and other solid phases. + * This %ThermoPhase object calculates its properties as a sum over other %LatticePhase objects. Each of the %LatticePhase + * objects is a %ThermoPhase object by itself. + * + * The results from this LatticeSolidPhase model reduces to the LatticePhase model when there is one + * lattice phase and the molar densities of the sublattice and the molar density within the LatticeSolidPhase + * have the same values. * + * + * The mole fraction vector is redefined witin the the LatticeSolidPhase object. Each of the mole + * fractions sum to one on each of the sublattices. The routine getMoleFraction() and setMoleFraction() + * have been redefined to use this convention. + * + *
+ *

Specification of Species Standard %State Properties

+ *
+ * + * The standard state properties are calculated in the normal way for each of the sublattices. The normal way + * here means that a thermodynamic polynomial in temperature is developed. Also, a constant volume approximation + * for the pressure dependence is assumed. All of these properties are on a Joules per kmol of sublattice + * constituent basis. + * + *
+ *

Specification of Solution Thermodynamic Properties

+ *
+ * The sum over the %LatticePhase objects is carried out by weighting each %LatticePhase object - * value with the molarDensity of the LatticePhase. Then the resulting quantity is divided by + * value with the molar density (kmol m-3) of its %LatticePhase. Then the resulting quantity is divided by * the molar density of the total compound. The LatticeSolidPhase object therefore only contains a - * listing of the number of Lattice Phases - * that comprises the solid and it contains a value for the molar density of the entire mixture. + * listing of the number of %LatticePhase object + * that comprises the solid, and it contains a value for the molar density of the entire mixture. + * This is the same thing as saying that + * + * \f[ + * L_i = L^{solid} \theta_i + * \f] + * + * \f$ L_i \f$ is the molar volume of the ith lattice. \f$ L^{solid} \f$ is the molar volume of the entire + * solid. \f$ \theta_i \f$ is a fixed weighting factor for the ith lattice representing the lattice + * stoichiometric coefficient. For this object the \f$ \theta_i \f$ values are fixed. + * * * Let's take FeS2 as an example, which may be thought of as a combination of two lattices: Fe and S lattice. * The Fe sublattice has a molar density of 1 gmol cm-3. The S sublattice has a molar density of 2 gmol cm-3. @@ -53,14 +85,34 @@ namespace Cantera { * associated with the sublattices. The Fe sublattice will have a weight of 1.0 associated with it. The * S sublattice will have a weight of 2.0 associated with it. * - * Currently, the molar density is set to a constant. + * + *
+ *

Specification of Solution Density Properties

+ *
* - * The results from this LatticeSolidPhase model reduces to the LatticePhase model when there is one - * lattice phase and the molar densities of the sublattice and the molar density within the LatticeSolidPhase - * have the same values. + * Currently, molar density is not a constant within the object, even though the species molar volumes are a + * constant. The basic idea is that a swelling of one of the sublattices will result in a swelling of + * of all of the lattices. Therefore, the molar volumes of the individual lattices are not independent of + * one another. * - * The mole fraction vector has been redefined within the LatticeSolidPhase object. The mole fractions sum - * to one within each of the individual lattice phases. The routine getMoleFraction() and setMoleFraction() + * The molar volume of the Lattice solid is calculated from the following formula + * + * \f[ + * V = \sum_i{ \theta_i V_i^{lattice}} + * \f] + * + * where \f$ V_i^{lattice} \f$ is the molar volume of the ith sublattice. This is calculated from the + * following standard formula. + * + * + * \f[ + * V_i = \sum_k{ \X_k V_k} + * \f] + * + * where k is a species in the ith sublattice. + * + * The mole fraction vector is redefined witin the the LatticeSolidPhase object. Each of the mole + * fractions sum to one on each of the sublattices. The routine getMoleFraction() and setMoleFraction() * have been redefined to use this convention. * * (This object is still under construction) @@ -104,6 +156,39 @@ namespace Cantera { */ virtual int eosType() const { return cLatticeSolid; } + //! Minimum temperature for which the thermodynamic data for the species + //! or phase are valid. + /*! + * If no argument is supplied, the + * value returned will be the lowest temperature at which the + * data for \e all species are valid. Otherwise, the value + * will be only for species \a k. This function is a wrapper + * that calls the species thermo minTemp function. + * + * @param k index of the species. Default is -1, which will return the max of the min value + * over all species. + */ + virtual doublereal minTemp(int k = -1) const; + + //! Maximum temperature for which the thermodynamic data for the species + //! are valid. + /*! + * If no argument is supplied, the + * value returned will be the highest temperature at which the + * data for \e all species are valid. Otherwise, the value + * will be only for species \a k. This function is a wrapper + * that calls the species thermo maxTemp function. + * + * @param k index of the species. Default is -1, which will return the min of the max value + * over all species. + */ + virtual doublereal maxTemp(int k = -1) const; + + + //! Returns the reference pressure in Pa. This function is a wrapper + //! that calls the species thermo refPressure function. + virtual doublereal refPressure() const ; + //! This method returns the convention used in specification //! of the standard state, of which there are currently two, //! temperature based, and variable pressure based. @@ -116,12 +201,11 @@ namespace Cantera { //! Return the Molar Enthalpy. Units: J/kmol. /*! - * The molar enthalpy is determined by the following formula, where \f$ C_n \f$ is the - * lattice molar density of the nth lattice, and \f$ C_T \f$ is the molar density - * of the solid compound. + * The molar enthalpy is determined by the following formula, where \f$ \theta_n \f$ is the + * lattice stoichiometric coefficient of the nth lattice * * \f[ - * \tilde h(T,P) = \frac{\sum_n C_n \tilde h_n(T,P) }{C_T}, + * \tilde h(T,P) = {\sum_n \theta_n \tilde h_n(T,P) } * \f] * * \f$ \tilde h_n(T,P) \f$ is the enthalpy of the nth lattice. @@ -133,12 +217,11 @@ namespace Cantera { //! Return the Molar Internal Energy. Units: J/kmol. /*! - * The molar internal energy is determined by the following formula, where \f$ C_n \f$ is the - * lattice molar density of the nth lattice, and \f$ C_T \f$ is the molar density - * of the solid compound. + * The molar enthalpy is determined by the following formula, where \f$ \theta_n \f$ is the + * lattice stoichiometric coefficient of the nth lattice * * \f[ - * \tilde u(T,P) = \frac{\sum_n C_n \tilde u_n(T,P) }{C_T}, + * \tilde u(T,P) = {\sum_n \theta_n \tilde u_n(T,P) } * \f] * * \f$ \tilde u_n(T,P) \f$ is the internal energy of the nth lattice. @@ -148,13 +231,12 @@ namespace Cantera { virtual doublereal intEnergy_mole() const; //! Return the Molar Entropy. Units: J/kmol/K. - /*! - * The molar entropy is determined by the following formula, where \f$ C_n \f$ is the - * lattice molar density of the nth lattice, and \f$ C_T \f$ is the molar density - * of the solid compound. + /*! + * The molar enthalpy is determined by the following formula, where \f$ \theta_n \f$ is the + * lattice stoichiometric coefficient of the nth lattice * * \f[ - * \tilde s(T,P) = \frac{\sum_n C_n \tilde s_n(T,P) }{C_T}, + * \tilde s(T,P) = \sum_n \theta_n \tilde s_n(T,P) * \f] * * \f$ \tilde s_n(T,P) \f$ is the molar entropy of the nth lattice. @@ -163,14 +245,13 @@ namespace Cantera { */ virtual doublereal entropy_mole() const; - //! Return the Molar Enthalpy. Units: J/kmol. + //! Return the Molar Gibbs energy. Units: J/kmol. /*! - * The molar enthalpy is determined by the following formula, where \f$ C_n \f$ is the - * lattice molar density of the nth lattice, and \f$ C_T \f$ is the molar density - * of the solid compound. + * The molar gibbs free energy is determined by the following formula, where \f$ \theta_n \f$ is the + * lattice stoichiometric coefficient of the nth lattice * * \f[ - * \tilde h(T,P) = \frac{\sum_n C_n \tilde h_n(T,P) }{C_T}, + * \tilde h(T,P) = {\sum_n \theta_n \tilde h_n(T,P) } * \f] * * \f$ \tilde h_n(T,P) \f$ is the enthalpy of the nth lattice. @@ -226,16 +307,30 @@ namespace Cantera { * * @param p Pressure (units - Pa) */ - virtual void setPressure(doublereal p) { - m_press = p; - setMolarDensity(m_molar_density); - } + virtual void setPressure(doublereal p); + + //! Calculate the density of the solid mixture + /*! + * The formula for this is + * + * \f[ + * \rho = \sum_n{ \rho_n \theta_n } + * \f] + * + * where \f$ \rho_n \f$ is the density of the nth sublattice + * + * Note this is a nonvirtual function. + */ + doublereal calcDensity(); //! Set the mole fractions to the specified values, and then //! normalize them so that they sum to 1.0 for each of the subphases - /*! + /* + * On input, the mole fraction vector is assumed to sum to one for each of the sublattices. The sublattices + * are updated with this mole fraction vector. The mole fractions are also storred within this object, after + * they are normalized to one by dividing by the number of sublattices. * - * @param x Input vector of mole fractions. There is no restriction + * @param x Input vector of mole fractions. There is no restriction * on the sum of the mole fraction vector. Internally, * this object will pass portions of this vector to the sublattices which assume that the portions * individually sum to one. @@ -358,20 +453,88 @@ namespace Cantera { * species in solution at the current temperature, pressure * and mole fraction of the solution. * + * This returns the underlying lattice chemical potentials, as the units are kmol-1 of + * the sublattice species. + * * @param mu Output vector of species chemical * potentials. Length: m_kk. Units: J/kmol */ virtual void getChemPotentials(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 + * pure species enthalpies + * \f[ + * \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k + * \f] + * The reference-state pure-species enthalpies, \f$ \hat h^{ref}_k(T) \f$, + * at the reference pressure,\f$ P_{ref} \f$, + * are computed by the species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + * + * @param hbar Output vector containing partial molar enthalpies. + * Length: m_kk. + */ + virtual void getPartialMolarEnthalpies(doublereal* hbar) const; + + /** + * Returns an array of partial molar entropies of the species in the + * solution. Units: J/kmol/K. + * For this phase, the partial molar entropies are equal to the + * pure species entropies plus the ideal solution contribution. + * \f[ + * \bar s_k(T,P) = \hat s^0_k(T) - R log(X_k) + * \f] + * The reference-state pure-species entropies,\f$ \hat s^{ref}_k(T) \f$, + * at the reference pressure, \f$ P_{ref} \f$, are computed by the + * species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + * + * @param sbar Output vector containing partial molar entropies. + * Length: m_kk. + */ + virtual void getPartialMolarEntropies(doublereal* sbar) const; + + /** + * Returns an array of partial molar Heat Capacities at constant + * pressure of the species in the + * solution. Units: J/kmol/K. + * For this phase, the partial molar heat capacities are equal + * to the standard state heat capacities. + * + * @param cpbar Output vector of partial heat capacities. Length: m_kk. + */ + virtual void getPartialMolarCp(doublereal* cpbar) const; + + /** + * returns an array of partial molar volumes of the species + * in the solution. Units: m^3 kmol-1. + * + * For this solution, thepartial molar volumes are equal to the + * constant species molar volumes. + * + * @param vbar Output vector of partial molar volumes. Length: m_kk. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + //! Get the array of standard state chemical potentials at unit activity for the species //! at their standard states at the current T and P of the solution. /*! * These are the standard state chemical potentials \f$ \mu^0_k(T,P) * \f$. The values are evaluated at the current - * temperature and pressure of the solution + * temperature and pressure of the solution. + * + * + * This returns the underlying lattice standard chemical potentials, as the units are kmol-1 of + * the sublattice species. * * @param mu0 Output vector of chemical potentials. - * Length: m_kk. + * Length: m_kk. Units: J/kmol */ virtual void getStandardChemPotentials(doublereal* mu0) const; @@ -400,6 +563,36 @@ namespace Cantera { * @param k index of the species (defaults to zero) */ virtual doublereal logStandardConc(int k=0) const; + //@} + /// @name Thermodynamic Values for the Species Reference States -------------------- + //@{ + + /** + * Returns the vector of nondimensional + * enthalpies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * + * This function fills in its one entry in hrt[] by calling + * the underlying species thermo function for the + * dimensionless gibbs free energy, calculated from the + * dimensionless enthalpy and entropy. + */ + virtual void getGibbs_RT_ref(doublereal *grt) const; + + + /** + * Returns the vector of the + * gibbs function of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * units = J/kmol + * + * This function fills in its one entry in g[] by calling + * the underlying species thermo functions for the + * gibbs free energy, calculated from enthalpy and the + * entropy, and the multiplying by RT. + */ + virtual void getGibbs_ref(doublereal *g) const; + //! Initialize the ThermoPhase object after all species have been set up /*! @@ -420,6 +613,9 @@ namespace Cantera { */ virtual void initThermo(); + //! Initialize vectors that depend on the number of species and sublattices + void initLengths(); + //! Add in species from Slave phases /*! * This hook is used for cSS_CONVENTION_SLAVE phases @@ -454,6 +650,16 @@ namespace Cantera { */ void setLatticeMoleFractionsByName(int n, std::string x); + //! Return a changeable reference to the calculation manager + //! for species reference-state thermodynamic properties + /*! + * This routine returns the calculation manager for the sublattice + * + * @param k Speices id. The default is -1, meaning return the default + * + * @internal + */ + virtual SpeciesThermo& speciesThermo(int k = -1); #ifdef H298MODIFY_CAPABILITY @@ -465,10 +671,7 @@ namespace Cantera { * @param k Species k * @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar */ - virtual void modifyOneHf298SS(const int k, const doublereal Hf298New) { - m_spthermo->modifyOneHf298(k, Hf298New); - m_tlast += 0.0001234; - } + virtual void modifyOneHf298SS(const int k, const doublereal Hf298New); #endif private: @@ -507,6 +710,16 @@ namespace Cantera { */ mutable vector_fp m_x; + //! Lattice stoichiometric coefficients + std::vector theta_; + + //! Temporary vector + mutable vector_fp tmpV_; + + std::vector nspLattice_; + + std::vector lkstart_; + private: //! Update the reference thermodynamic functions diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index ea91a0cd7..d6d2dc451 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -525,7 +525,9 @@ namespace Cantera { /*************************************************************** * Add the elements. ***************************************************************/ - th->addElementsFromXML(phase); + if (ssConvention != cSS_CONVENTION_SLAVE) { + th->addElementsFromXML(phase); + } /*************************************************************** * Add the species. @@ -538,11 +540,13 @@ namespace Cantera { vector sparrays; phase.getChildren("speciesArray", sparrays); int jsp, nspa = static_cast(sparrays.size()); - if (nspa == 0) { - throw CanteraError("importPhase", - "phase, " + th->id() + ", has zero \"speciesArray\" XML nodes.\n" + if (ssConvention != cSS_CONVENTION_SLAVE) { + if (nspa == 0) { + throw CanteraError("importPhase", + "phase, " + th->id() + ", has zero \"speciesArray\" XML nodes.\n" + " There must be at least one speciesArray nodes " "with one or more species"); + } } vector dbases; vector_int sprule(nspa,0); @@ -608,7 +612,7 @@ namespace Cantera { // If the phase has a species thermo manager already installed, // delete it since we are adding new species. - delete &th->speciesThermo(); + //delete &th->speciesThermo(); // Decide whether the the phase has a variable pressure ss or not SpeciesThermo* spth = 0; diff --git a/Cantera/src/thermo/ThermoPhase.cpp b/Cantera/src/thermo/ThermoPhase.cpp index 42ec63f7d..0d8a11618 100644 --- a/Cantera/src/thermo/ThermoPhase.cpp +++ b/Cantera/src/thermo/ThermoPhase.cpp @@ -822,9 +822,47 @@ namespace Cantera { if (i == 4) uA[4] = 0.0; if (i == 5) uA[5] = 0.0; } + } + //================================================================================================================= + // Install a species thermodynamic property manager. + /* + * The species thermodynamic property manager + * computes properties of the pure species for use in + * constructing solution properties. It is meant for internal + * use, and some classes derived from ThermoPhase may not use + * any species thermodynamic property manager. This method is + * called by function importPhase() in importCTML.cpp. + * + * @param spthermo input pointer to the species thermodynamic property + * manager. + * + * @internal + */ + void ThermoPhase::setSpeciesThermo(SpeciesThermo* spthermo) { + if (m_spthermo) { + if (m_spthermo != spthermo) { + delete m_spthermo; + } + } + m_spthermo = spthermo; + } + //================================================================================================================= + // Return a changeable reference to the calculation manager + // for species reference-state thermodynamic properties + /* + * + * @param k Speices id. The default is -1, meaning return the default + * + * @internal + */ + SpeciesThermo& ThermoPhase::speciesThermo(int k) { + if (!m_spthermo) { + throw CanteraError("ThermoPhase::speciesThermo()", + "species reference state thermo manager was not set"); + } + return *m_spthermo; } //================================================================================================================= - /* * initThermoFile(): * diff --git a/Cantera/src/thermo/ThermoPhase.h b/Cantera/src/thermo/ThermoPhase.h index 8e35524d0..a27d999b8 100644 --- a/Cantera/src/thermo/ThermoPhase.h +++ b/Cantera/src/thermo/ThermoPhase.h @@ -758,7 +758,7 @@ namespace Cantera { * Returns the reference pressure in Pa. This function is a wrapper * that calls the species thermo refPressure function. */ - doublereal refPressure() const { + virtual doublereal refPressure() const { return m_spthermo->refPressure(); } @@ -775,7 +775,7 @@ namespace Cantera { * @param k index of the species. Default is -1, which will return the max of the min value * over all species. */ - doublereal minTemp(int k = -1) const { + virtual doublereal minTemp(int k = -1) const { return m_spthermo->minTemp(k); } @@ -844,7 +844,7 @@ namespace Cantera { * @param k index of the species. Default is -1, which will return the min of the max value * over all species. */ - doublereal maxTemp(int k = -1) const { + virtual doublereal maxTemp(int k = -1) const { return m_spthermo->maxTemp(k); } @@ -1883,19 +1883,17 @@ namespace Cantera { * * @internal */ - void setSpeciesThermo(SpeciesThermo* spthermo) - { m_spthermo = spthermo; } + void setSpeciesThermo(SpeciesThermo* spthermo); //! Return a changeable reference to the calculation manager //! for species reference-state thermodynamic properties /*! * - * @todo This method will fail if no species thermo - * manager has been installed. + * @param k Speices id. The default is -1, meaning return the default * * @internal */ - SpeciesThermo& speciesThermo() { return *m_spthermo; } + virtual SpeciesThermo& speciesThermo(int k = -1); /** * @internal diff --git a/Cantera/src/thermo/mix_defs.h b/Cantera/src/thermo/mix_defs.h index 02ab93eb9..b104d12dc 100644 --- a/Cantera/src/thermo/mix_defs.h +++ b/Cantera/src/thermo/mix_defs.h @@ -47,7 +47,7 @@ namespace Cantera { const int cMetalSHEelectrons = 9; // SHE electrode electrons const int cLatticeSolid = 20; // LatticeSolidPhase.h - const int cLattice = 21; + const int cLattice = 21; //LatticePhase.h // pure fluids with liquid/vapor eqs of state const int cPureFluid = 10;