From a96f86540381cb693dd91a22e725d29303ca1dcd Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 27 Oct 2010 17:00:46 +0000 Subject: [PATCH] Added FixedChemPotSSTP object. This is a utility object that will provide a bath species at a certain chemical potential to the equilibrium solver. --- Cantera/src/thermo/FixedChemPotSSTP.cpp | 492 ++++++++++++++++++ Cantera/src/thermo/FixedChemPotSSTP.h | 630 ++++++++++++++++++++++++ Cantera/src/thermo/Makefile.in | 4 +- Cantera/src/thermo/SingleSpeciesTP.cpp | 14 +- Cantera/src/thermo/ThermoFactory.cpp | 14 +- Cantera/src/thermo/mix_defs.h | 3 + 6 files changed, 1148 insertions(+), 9 deletions(-) create mode 100644 Cantera/src/thermo/FixedChemPotSSTP.cpp create mode 100644 Cantera/src/thermo/FixedChemPotSSTP.h diff --git a/Cantera/src/thermo/FixedChemPotSSTP.cpp b/Cantera/src/thermo/FixedChemPotSSTP.cpp new file mode 100644 index 000000000..38aabe132 --- /dev/null +++ b/Cantera/src/thermo/FixedChemPotSSTP.cpp @@ -0,0 +1,492 @@ +/** + * @file FixedChemPotSSTP.cpp + * Definition file for the FixedChemPotSSTP class, which represents a fixed-composition + * incompressible substance with a constant chemical potential (see \ref thermoprops and + * class \link Cantera::FixedChemPotSSTP FixedChemPotSSTP\endlink) + */ + +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + * + */ + +/* + * $Id: FixedChemPotSSTP.cpp 255 2009-11-09 23:36:49Z hkmoffa $ + */ + +#include "ct_defs.h" +#include "mix_defs.h" +#include "FixedChemPotSSTP.h" +#include "SpeciesThermo.h" +#include "ThermoFactory.h" + +#include + +namespace Cantera { + //==================================================================================================================== + /* + * ---- Constructors ------- + */ + //==================================================================================================================== + /* + * Default Constructor for the FixedChemPotSSTP class + */ + FixedChemPotSSTP::FixedChemPotSSTP() : + SingleSpeciesTP(), + chemPot_(0.0) + { + } + //==================================================================================================================== + // Create and initialize a FixedChemPotSSTP ThermoPhase object + // from an asci input file + /* + * @param infile name of the input file + * @param id name of the phase id in the file. + * If this is blank, the first phase in the file is used. + */ + FixedChemPotSSTP::FixedChemPotSSTP(std::string infile, std::string id) : + SingleSpeciesTP(), + chemPot_(0.0) + { + XML_Node* root = get_XML_File(infile); + if (id == "-") id = ""; + XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, root); + if (!xphase) { + throw CanteraError("FixedChemPotSSTP::FixedChemPotSSTP", + "Couldn't find phase name in file:" + id); + } + // Check the model name to ensure we have compatibility + const XML_Node& th = xphase->child("thermo"); + std::string model = th["model"]; + if (model != "StoichSubstance" && model != "StoichSubstanceSSTP" && model != "FixedChemPot") { + throw CanteraError("FixedChemPotSSTP::FixedChemPotSSTP", + "thermo model attribute must be FixedChemPot or StoichSubstance"); + } + importPhase(*xphase, this); + } + //==================================================================================================================== + // Full Constructor. + /* + * @param phaseRef XML node pointing to a FixedChemPotSSTP description + * @param id Id of the phase. + */ + FixedChemPotSSTP::FixedChemPotSSTP(XML_Node& xmlphase, std::string id) : + SingleSpeciesTP(), + chemPot_(0.0) + { + if (id != "") { + std::string idxml = xmlphase["id"]; + if (id != idxml) { + throw CanteraError("FixedChemPotSSTP::FixedChemPotSSTP", + "id's don't match"); + } + } + const XML_Node& th = xmlphase.child("thermo"); + std::string model = th["model"]; + if (model != "StoichSubstance" && model != "StoichSubstanceSSTP" && model != "FixedChemPotSSTP") { + throw CanteraError("FixedChemPotSSTP::FixedChemPotSSTP", + "thermo model attribute must be StoichSubstance or FixedChemPot"); + } + importPhase(xmlphase, this); + + if (model == "StoichSubstance" || model == "StoichSubstanceSSTP") { + _updateThermo(); + chemPot_ = (m_h0_RT[0] - m_s0_R[0]) * GasConstant * temperature(); + } + } + //==================================================================================================================== + // Copy constructor + /* + * @param right Object to be copied + */ + FixedChemPotSSTP::FixedChemPotSSTP(const FixedChemPotSSTP &right) : + SingleSpeciesTP() + { + *this = operator=(right); + } + //==================================================================================================================== + // Assignment operator + /* + * @param right Object to be copied + */ + FixedChemPotSSTP & + FixedChemPotSSTP::operator=(const FixedChemPotSSTP & right) { + if (&right != this) { + SingleSpeciesTP::operator=(right); + + chemPot_ = right.chemPot_; + } + return *this; + } + //==================================================================================================================== + /* + * Destructor for the routine (virtual) + * + */ + FixedChemPotSSTP::~FixedChemPotSSTP() + { + } + //==================================================================================================================== + // Duplication function + /* + * This virtual function is used to create a duplicate of the + * current phase. It's used to duplicate the phase when given + * a ThermoPhase pointer to the phase. + * + * @return It returns a ThermoPhase pointer. + */ + ThermoPhase *FixedChemPotSSTP::duplMyselfAsThermoPhase() const { + FixedChemPotSSTP *stp = new FixedChemPotSSTP(*this); + return (ThermoPhase *) stp; + } + //==================================================================================================================== + + /* + * ---- Utilities ----- + */ + + /* + * Equation of state flag. Returns the value cStoichSubstance, + * defined in mix_defs.h. + */ + int FixedChemPotSSTP::eosType() const { + return cFixedChemPot; + } + + /* + * ---- Molar Thermodynamic properties of the solution ---- + */ + + /* + * ----- Mechanical Equation of State ------ + */ + //==================================================================================================================== + /* + * Pressure. Units: Pa. + * For an incompressible substance, the density is independent + * of pressure. This method simply returns the stored + * pressure value. + */ + doublereal FixedChemPotSSTP::pressure() const { + return m_press; + } + //==================================================================================================================== + /* + * Set the pressure at constant temperature. Units: Pa. + * For an incompressible substance, the density is + * independent of pressure. Therefore, this method only + * stores the specified pressure value. It does not + * modify the density. + */ + void FixedChemPotSSTP::setPressure(doublereal p) { + m_press = p; + } + //==================================================================================================================== + /* + * The isothermal compressibility. Units: 1/Pa. + * The isothermal compressibility is defined as + * \f[ + * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T + * \f] + * + * It's equal to zero for this model, since the molar volume + * doesn't change with pressure or temperature. + */ + doublereal FixedChemPotSSTP::isothermalCompressibility() const { + return 0.0; + } + //==================================================================================================================== + /* + * The thermal expansion coefficient. Units: 1/K. + * The thermal expansion coefficient is defined as + * + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + * + * It's equal to zero for this model, since the molar volume + * doesn't change with pressure or temperature. + */ + doublereal FixedChemPotSSTP::thermalExpansionCoeff() const { + return 0.0; + } + //==================================================================================================================== + /* + * ---- Chemical Potentials and Activities ---- + */ + //==================================================================================================================== + /* + * This method returns the array of generalized + * concentrations. For a stoichiometric substance, there is + * only one species, and the generalized concentration is 1.0. + */ + void FixedChemPotSSTP:: + getActivityConcentrations(doublereal* c) const { + c[0] = 1.0; + } + //==================================================================================================================== + /* + * The standard concentration. This is defined as the concentration + * by which the generalized concentration is normalized to produce + * the activity. + */ + doublereal FixedChemPotSSTP::standardConcentration(int k) const { + return 1.0; + } + //==================================================================================================================== + /* + * Returns the natural logarithm of the standard + * concentration of the kth species + */ + doublereal FixedChemPotSSTP::logStandardConc(int k) const { + return 0.0; + } + //==================================================================================================================== + /* + * Returns the units of the standard and generalized + * concentrations Note they have the same units, as their + * ratio is defined to be equal to the activity of the kth + * species in the solution, which is unitless. + * + * This routine is used in print out applications where the + * units are needed. Usually, MKS units are assumed throughout + * the program and in the XML input files. + * + * uA[0] = kmol units - default = 1 + * uA[1] = m units - default = -nDim(), the number of spatial + * dimensions in the Phase class. + * uA[2] = kg units - default = 0; + * uA[3] = Pa(pressure) units - default = 0; + * uA[4] = Temperature units - default = 0; + * uA[5] = time units - default = 0 + */ + void FixedChemPotSSTP:: + getUnitsStandardConc(doublereal *uA, int k, int sizeUA) const { + for (int i = 0; i < 6; i++) { + uA[i] = 0; + } + } + //==================================================================================================================== + /* + * ---- Partial Molar Properties of the Solution ---- + */ + void FixedChemPotSSTP::getPartialMolarVolumes(doublereal* vbar) const { + vbar[0] = 0.0; + } + //==================================================================================================================== + /* + * ---- Properties of the Standard State of the Species in the Solution + * ---- + */ + //==================================================================================================================== + /* + * Get the array of chemical potentials at unit activity + * \f$ \mu^0_k \f$. + * + * For a stoichiometric substance, there is no activity term in + * the chemical potential expression, and therefore the + * standard chemical potential and the chemical potential + * are both equal to the molar Gibbs function. + */ + void FixedChemPotSSTP:: + getStandardChemPotentials(doublereal* mu0) const { + mu0[0] = chemPot_; + } + //==================================================================================================================== + /* + * Get the nondimensional Enthalpy functions for the species + * at their standard states at the current + * T and P of the solution. + * Molar enthalpy. Units: J/kmol. For an incompressible, + * stoichiometric substance, the internal energy is + * independent of pressure, and therefore the molar enthalpy + * is \f[ \hat h(T, P) = \hat u(T) + P \hat v \f], where the + * molar specific volume is constant. + */ + void FixedChemPotSSTP::getEnthalpy_RT(doublereal* hrt) const { + double rt = _RT(); + hrt[0] = chemPot_ / rt; + } + //==================================================================================================================== + /* + * Get the array of nondimensional Entropy functions for the + * standard state species + * at the current T and P of the solution. + */ + void FixedChemPotSSTP::getEntropy_R(doublereal* sr) const { + sr[0] = 0.0; + } + //==================================================================================================================== + /* + * Get the nondimensional Gibbs functions for the species + * at their standard states of solution at the current T and P + * of the solution + */ + void FixedChemPotSSTP::getGibbs_RT(doublereal* grt) const { + double rt = _RT(); + grt[0] = chemPot_ / rt; + } + //==================================================================================================================== + /* + * Get the nondimensional Gibbs functions for the standard + * state of the species at the current T and P. + */ + void FixedChemPotSSTP::getCp_R(doublereal* cpr) const { + cpr[0] = 0.0; + } + //==================================================================================================================== + /* + * Molar internal energy (J/kmol). + * For an incompressible, + * stoichiometric substance, the molar internal energy is + * independent of pressure. Since the thermodynamic properties + * are specified by giving the standard-state enthalpy, the + * term \f$ P_0 \hat v\f$ is subtracted from the specified molar + * enthalpy to compute the molar internal energy. + */ + void FixedChemPotSSTP::getIntEnergy_RT(doublereal* urt) const { + urt[0] = chemPot_; + } + //==================================================================================================================== + // Get the molar volumes of each species in their standard + // states at the current T and P of the solution. + /* + * units = m^3 / kmol + * + * We set this to zero + * + * @param vbar On output this contains the standard volume of the species + * and phase (m^3/kmol). Vector of length 1 + */ + void FixedChemPotSSTP::getStandardVolumes(doublereal* vbar) const { + vbar[0] = 0.0; + } + //==================================================================================================================== + /* + * ---- Thermodynamic Values for the Species Reference States ---- + */ + //==================================================================================================================== + void FixedChemPotSSTP::getIntEnergy_RT_ref(doublereal* urt) const { + urt[0] = chemPot_; + } + //==================================================================================================================== + void FixedChemPotSSTP::getEnthalpy_RT_ref(doublereal* hrt) const { + double rt = _RT(); + hrt[0] = chemPot_ / rt; + } + //==================================================================================================================== + void FixedChemPotSSTP::getEntropy_R_ref(doublereal* sr) const { + sr[0] = 0.0; + } + //==================================================================================================================== + void FixedChemPotSSTP::getGibbs_RT_ref(doublereal* grt) const { + double rt = _RT(); + grt[0] = chemPot_ / rt; + } + //==================================================================================================================== + void FixedChemPotSSTP::getGibbs_ref(doublereal* g) const { + g[0] = chemPot_; + } + //==================================================================================================================== + void FixedChemPotSSTP::getCp_R_ref(doublereal* cpr) const { + cpr[0] = 0.0; + } + //==================================================================================================================== + /* + * ---- Saturation Properties + */ + //==================================================================================================================== + /* + * ---- Initialization and Internal functions + */ + //==================================================================================================================== + /* + * @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 just prior to returning + * from function importPhase. + * + * @see importCTML.cpp + */ + void FixedChemPotSSTP::initThermo() { + /* + * Call the base class thermo initializer + */ + SingleSpeciesTP::initThermo(); + } + //==================================================================================================================== + + void FixedChemPotSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("FixedChemPotSSTP::initThermoXML", "no thermo XML node"); + } + XML_Node &tnode = phaseNode.child("thermo"); + std::string model = tnode["model"]; + if (model != "StoichSubstance" && model != "FixedChemPot" && model != "StoichSubstanceSSTP") { + throw CanteraError("FixedChemPotSSTP::initThermoXML()", + "thermo model attribute must be FixedChemPot or StoichSubstance or StoichSubstanceSSTP"); + } + if (model == "FixedChemPot") { + double val = getFloatDefaultUnits(tnode, "chemicalPotential", "J/kmol"); + chemPot_ = val; + } + SingleSpeciesTP::initThermoXML(phaseNode, id); + + + } + //==================================================================================================================== + /* + * setParameters: + * + * Generic routine that is used to set the parameters used + * by this model. + * C[0] = density of phase [ kg/m3 ] + */ + void FixedChemPotSSTP::setParameters(int n, doublereal * const c) { + chemPot_ = c[0]; + } + //==================================================================================================================== + /* + * getParameters: + * + * Generic routine that is used to get the parameters used + * by this model. + * n = 1 + * C[0] = density of phase [ kg/m3 ] + */ + void FixedChemPotSSTP::getParameters(int &n, doublereal * const c) const { + n = 1; + c[0] = chemPot_; + } + //==================================================================================================================== + void FixedChemPotSSTP::setParametersFromXML(const XML_Node& eosdata) { + std::string model = eosdata["model"]; + if (model != "StoichSubstance" && model != "FixedChemPot" && model != "StoichSubstanceSSTP") { + throw CanteraError("FixedChemPotSSTP::setParametersFromXML", + "thermo model attribute must be FixedChemPot or StoichSubstance or StoichSubstanceSSTP"); + } + if (model == "FixedChemPotSSTP") { + doublereal val = getFloatDefaultUnits(eosdata, "chemicalPotential", "J/kmol"); + chemPot_ = val; + } + } + //==================================================================================================================== + // Function to set the chemical potential directly + /* + * @param chemPot Value of the chemical potential (units J/kmol) + */ + void FixedChemPotSSTP::setChemicalPotential(doublereal chemPot) { + chemPot_ = chemPot; + } + //==================================================================================================================== +} diff --git a/Cantera/src/thermo/FixedChemPotSSTP.h b/Cantera/src/thermo/FixedChemPotSSTP.h new file mode 100644 index 000000000..1df7ec613 --- /dev/null +++ b/Cantera/src/thermo/FixedChemPotSSTP.h @@ -0,0 +1,630 @@ +/** + * @file FixedChemPotSSTP.h + * Header file for the FixedChemPotSSTP class, which represents a fixed-composition + * incompressible substance with a constant chemical potential (see \ref thermoprops and + * class \link Cantera::FixedChemPotSSTP FixedChemPotSSTP\endlink) + */ + +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +/* + * $Date: 2009-11-09 16:36:49 -0700 (Mon, 09 Nov 2009) $ + * $Revision: 255 $ + */ + +#ifndef CT_FIXEDCHEMPOTSSTP_H +#define CT_FIXEDCHEMPOTSSTP_H + +#include "mix_defs.h" +#include "SingleSpeciesTP.h" +#include "SpeciesThermo.h" + +namespace Cantera { + + //! Class %FixedChemPotSSTP represents a stoichiometric (fixed + //! composition) incompressible substance. + /*! + * This class internally changes the independent degree of freedom from + * density to pressure. This is necessary because the phase is + * incompressible. It uses a zero volume approximation. + * + * + * Specification of Species Standard %State Properties + * + * This class inherits from SingleSpeciesTP. + * It uses a single value for the chemical potential which is assumed to be constant + * with respect to temperature and pressure. + * + * The reference state thermodynamics is inherited from SingleSpeciesTP. However, + * it's only used to set the initial chemical potential to the value + * of the chemical potential at the starting conditions. Thereafter, + * it is ignored. + * + * For a zero volume material, the internal energy and the enthalpy are + * equal to the chemical potential. The entropy, the heat capacity, and the molar volume + * are equal to zero. + * + * + * Specification of Solution Thermodynamic Properties + * + * All solution properties are obtained from the standard state + * species functions, since there is only one species in the phase. + * + * Application within %Kinetics Managers + * + * The standard concentration is equal to 1.0. This means that the + * kinetics operator works on an (activities basis). Since this + * is a stoichiometric substance, this means that the concentration + * of this phase drops out of kinetics expressions. + * + * An example of a reaction using this is a sticking coefficient + * reaction of a substance in an ideal gas phase on a surface with a bulk phase + * species in this phase. In this case, the rate of progress for this + * reaction, \f$ R_s \f$, may be expressed via the following equation: + * \f[ + * R_s = k_s C_{gas} + * \f] + * where the units for \f$ R_s \f$ are kmol m-2 s-1. \f$ C_{gas} \f$ has units + * of kmol m-3. Therefore, the kinetic rate constant, \f$ k_s \f$, has + * units of m s-1. Nowhere does the concentration of the bulk phase + * appear in the rate constant expression, since it's a stoichiometric + * phase, and the activity is always equal to 1.0. + * + * Instanteation of the Class + * + * The constructor for this phase is located in the default ThermoFactory + * for %Cantera. A new %FixedChemPotSSTP may be created by a standalone xml file + * which is given below. + * + * It may also be created by the following code snippets. The code + * includes the special member function setChemicalPotential( chempot), which + * sets the chemical potential to a specific value in J / kmol. + * + * @code + * sprintf(file_ID,"%s#Li(Fixed)", iFile); + * XML_Node *xm = get_XML_NameID("phase", file_ID, 0); + * FixedChemPotSSTP *LiFixed = new FixedChemPotSSTP(*xm); + // Set the chemical potential to -2.3E7 J/kmol + * LiFixed->setChemicalPotential(-2.3E7.) + * @endcode + * + * or by the following call to importPhase(): + * + * @code + * sprintf(file_ID,"%s#NaCl(S)", iFile); + * XML_Node *xm = get_XML_NameID("phase", file_ID, 0); + * FixedChemPotSSTP solid; + * importPhase(*xm, &solid); + * @endcode + * + * XML Example + * + * The phase model name for this is called StoichSubstance. It must be supplied + * as the model attribute of the thermo XML element entry. + * Within the phase XML block, + * the density of the phase must be specified. An example of an XML file + * this phase is given below. + * + * @verbatim + + + + + + + + Li + + + LiFixed + + + -2.3E7 + + + + + + + + + + Li:1 + + + + 50.72389, 6.672267, -2.517167, + 10.15934, -0.200675, -427.2115, + 130.3973 + + + + + + + @endverbatim + * + * The model attribute, "FixedChemPot", on the thermo element + * identifies the phase as being a FixedChemPotSSTP object. + * + * @ingroup thermoprops + */ + class FixedChemPotSSTP : public SingleSpeciesTP { + + public: + + //! Default constructor for the FixedChemPotSSTP class + FixedChemPotSSTP(); + + //! Construct and initialize a FixedChemPotSSTP ThermoPhase object + //! directly from an asci input file + /*! + * @param infile name of the input file + * @param id name of the phase id in the file. + * If this is blank, the first phase in the file is used. + */ + FixedChemPotSSTP(std::string infile, std::string id = ""); + + //! Construct and initialize a FixedChemPotSSTP ThermoPhase object + //! directly from an XML database + /*! + * @param phaseRef XML node pointing to a FixedChemPotSSTP description + * @param id Id of the phase. + */ + FixedChemPotSSTP(XML_Node& phaseRef, std::string id = ""); + + //! Copy constructor + /*! + * @param right Object to be copied + */ + FixedChemPotSSTP(const FixedChemPotSSTP &right); + + //! Assignment operator + /*! + * @param right Object to be copied + */ + FixedChemPotSSTP & operator=(const FixedChemPotSSTP & right); + + //! Destructor for the routine (virtual) + virtual ~FixedChemPotSSTP(); + + //! Duplication function + /*! + * This virtual function is used to create a duplicate of the + * current phase. It's used to duplicate the phase when given + * a ThermoPhase pointer to the phase. + * + * @return It returns a ThermoPhase pointer. + */ + ThermoPhase *duplMyselfAsThermoPhase() const; + + /** + * + * @name Utilities + * @{ + */ + + /** + * Equation of state flag. + * + * Returns the value cStoichSubstance, defined in mix_defs.h. + */ + virtual int eosType() const; + + /** + * @} + * @name Molar Thermodynamic Properties of the Solution + * @{ + */ + + /** + * @} + * @name Mechanical Equation of State + * @{ + */ + + + //! Report the Pressure. Units: Pa. + /*! + * For an incompressible substance, the density is independent + * of pressure. This method simply returns the storred + * pressure value. + */ + virtual doublereal pressure() const; + + //! Set the pressure at constant temperature. Units: Pa. + /*! + * For an incompressible substance, the density is + * independent of pressure. Therefore, this method only + * stores the specified pressure value. It does not + * modify the density. + * + * @param p Pressure (units - Pa) + */ + virtual void setPressure(doublereal p); + + //! Returns the isothermal compressibility. Units: 1/Pa. + /*! + * The isothermal compressibility is defined as + * \f[ + * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T + * \f] + */ + virtual doublereal isothermalCompressibility() const; + + //! Return the volumetric thermal expansion coefficient. Units: 1/K. + /*! + * The thermal expansion coefficient is defined as + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + */ + virtual doublereal thermalExpansionCoeff() const ; + + /** + * @} + * @name Activities, Standard States, and Activity Concentrations + * + * This section is largely handled by parent classes, since there + * is only one species. Therefore, the activity is equal to one. + * @{ + */ + + //! This method returns an array of generalized concentrations + /*! + * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / + * C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration + * defined below and \f$ a_k \f$ are activities used in the + * thermodynamic functions. These activity (or generalized) + * concentrations are used + * by kinetics manager classes to compute the forward and + * reverse rates of elementary reactions. + * + * For a stoichiomeetric substance, there is + * only one species, and the generalized concentration is 1.0. + * + * @param c Output array of generalized concentrations. The + * units depend upon the implementation of the + * reaction rate expressions within the phase. + */ + virtual void getActivityConcentrations(doublereal* c) const; + + //! Return the standard concentration for the kth species + /*! + * The standard concentration \f$ C^0_k \f$ used to normalize + * the activity (i.e., generalized) concentration. + * This phase assumes that the kinetics operator works on an + * dimensionless basis. Thus, the standard concentration is + * equal to 1.0. + * + * @param k Optional parameter indicating the species. The default + * is to assume this refers to species 0. + * @return + * Returns The standard Concentration as 1.0 + */ + virtual doublereal standardConcentration(int k=0) const; + + //! Natural logarithm of the standard concentration of the kth species. + /*! + * @param k index of the species (defaults to zero) + */ + virtual doublereal logStandardConc(int k=0) const; + + //! Get the array of chemical potentials at unit activity for the species + //! at their standard states at the current T and P of the solution. + /*! + * For a stoichiometric substance, there is no activity term in + * the chemical potential expression, and therefore the + * standard chemical potential and the chemical potential + * are both equal to the molar Gibbs function. + * + * 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 + * + * @param mu0 Output vector of chemical potentials. + * Length: m_kk. + */ + virtual void getStandardChemPotentials(doublereal* mu0) const; + + //! Returns the units of the standard and generalized concentrations. + /*! + * Note they have the same units, as their + * ratio is defined to be equal to the activity of the kth + * species in the solution, which is unitless. + * + * This routine is used in print out applications where the + * units are needed. Usually, MKS units are assumed throughout + * the program and in the XML input files. + * + * The base %ThermoPhase class assigns thedefault quantities + * of (kmol/m3) for all species. + * Inherited classes are responsible for overriding the default + * values if necessary. + * + * @param uA Output vector containing the units + * uA[0] = kmol units - default = 1 + * uA[1] = m units - default = -nDim(), the number of spatial + * dimensions in the Phase class. + * uA[2] = kg units - default = 0; + * uA[3] = Pa(pressure) units - default = 0; + * uA[4] = Temperature units - default = 0; + * uA[5] = time units - default = 0 + * @param k species index. Defaults to 0. + * @param sizeUA output int containing the size of the vector. + * Currently, this is equal to 6. + */ + virtual void getUnitsStandardConc(doublereal *uA, int k = 0, + int sizeUA = 6) const; + + //@} + /// @name Partial Molar Properties of the Solution + /// + /// These properties are handled by the parent class, + /// SingleSpeciesTP + //@{ + + //! Get the species partial molar volumes. Units: m^3/kmol. + /*! + * This is the phase molar volume. \f$ V(T,P) = V_o(T,P) \f$. + * + * set to zero. + * + * @param vbar On return, contains the molar volume of the single species + * and the phase. Units are m^3 / kmol. Length = 1 + */ + void getPartialMolarVolumes(doublereal* vbar) const; + + //@} + /// @name Properties of the Standard State of the Species in the Solution + //@{ + + //! Get the nondimensional Enthalpy functions for the species + //! at their standard states at the current T and P of the solution. + /*! + * @param hrt Output vector of nondimensional standard state enthalpies. + * Length: m_kk. + */ + virtual void getEnthalpy_RT(doublereal* hrt) const; + + //! Get the array of nondimensional Entropy functions for the + //! standard state species at the current T and P of the solution. + /*! + * @param sr Output vector of nondimensional standard state entropies. + * Length: m_kk. + */ + virtual void getEntropy_R(doublereal* sr) const; + + //! Get the nondimensional Gibbs functions for the species + //! in their standard states at the current T and P of the solution. + /*! + * @param grt Output vector of nondimensional standard state gibbs free energies + * Length: m_kk. + */ + virtual void getGibbs_RT(doublereal* grt) const; + + //! Get the nondimensional Heat Capacities at constant + //! pressure for the species standard states + //! at the current T and P of the solution + /*! + * @param cpr Output vector of nondimensional standard state heat capacities + * Length: m_kk. + */ + virtual void getCp_R(doublereal* cpr) const; + + //! Returns the vector of nondimensional Internal Energies of the standard + //! state species at the current T and P of the solution + /*! + * For an incompressible, + * stoichiometric substance, the molar internal energy is + * independent of pressure. Since the thermodynamic properties + * are specified by giving the standard-state enthalpy, the + * term \f$ P_{ref} \hat v\f$ is subtracted from the specified reference molar + * enthalpy to compute the standard state molar internal energy. + * + * @param urt output vector of nondimensional standard state + * internal energies of the species. Length: m_kk. + */ + virtual void getIntEnergy_RT(doublereal* urt) const; + + //! Get the molar volumes of each species in their standard + //! states at the current T and P of the solution. + /* + * units = m^3 / kmol + * + * We set this to zero + * + * @param vbar On output this contains the standard volume of the species + * and phase (m^3/kmol). Vector of length 1 + */ + virtual void getStandardVolumes(doublereal* vbar) const; + + //@} + /// @name Thermodynamic Values for the Species Reference States + //@{ + + //! Returns the vector of nondimensional + //! internal Energies of the reference state at the current temperature + //! of the solution and the reference pressure for each species. + /*! + * @param urt Output vector of nondimensional reference state + * internal energies of the species. + * Length: m_kk + */ + virtual void getIntEnergy_RT_ref(doublereal *urt) const; + + //@} + /// @name Thermodynamic Values for the Species Reference State + /// + + /*! + * 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 is resolved in this class. It is assumed that the m_spthermo species thermo + * pointer is populated and yields the reference state. + * + * @param hrt Output vector containing the nondimensional reference state enthalpies + * Length: m_kk. + */ + virtual void getEnthalpy_RT_ref(doublereal *hrt) const; + + /*! + * 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 is resolved in this class. It is assumed that the m_spthermo species thermo + * pointer is populated and yields the reference state. + * + * @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 + * + * This function is resolved in this class. It is assumed that the m_spthermo species thermo + * pointer is populated and yields the reference state. + * + * @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 the vector of nondimensional + * entropies of the reference state at the current temperature + * of the solution and the reference pressure for each species. + * + * This function is resolved in this class. It is assumed that the m_spthermo species thermo + * pointer is populated and yields the reference state. + * + * @param er Output vector containing the nondimensional reference state + * entropies. Length: m_kk. + */ + virtual void getEntropy_R_ref(doublereal *er) const; + + /*! + * Returns the vector of nondimensional + * constant pressure heat capacities of the reference state + * at the current temperature of the solution + * and reference pressure for each species. + * + * This function is resolved in this class. It is assumed that the m_spthermo species thermo + * pointer is populated and yields the reference state. + * + * @param cprt Output vector of nondimensional reference state + * heat capacities at constant pressure for the species. + * Length: m_kk + */ + virtual void getCp_R_ref(doublereal *cprt) const; + + + + + /* + * ---- Critical State Properties + */ + + + /* + * ---- Saturation Properties + */ + + /* + * @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 just prior to returning + * from function importPhase. + * + * @see importCTML.cpp + */ + virtual void initThermo(); + + + virtual void initThermoXML(XML_Node& phaseNode, std::string id); + + //! Set the equation of state parameters + /*! + * @internal + * The number and meaning of these depends on the subclass. + * + * @param n number of parameters + * @param c array of \a n coefficients + * c[0] = density of phase [ kg/m3 ] + */ + virtual void setParameters(int n, doublereal * const c); + + //! Get the equation of state parameters in a vector + /*! + * @internal + * + * @param n number of parameters + * @param c array of \a n coefficients + * + * For this phase: + * - n = 1 + * - c[0] = density of phase [ kg/m3 ] + */ + virtual void getParameters(int &n, doublereal * const c) const; + + //! Set equation of state parameter values from XML entries. + /*! + * This method is called by function importPhase() in + * file importCTML.cpp when processing a phase definition in + * an input file. It should be overloaded in subclasses to set + * any parameters that are specific to that particular phase + * model. Note, this method is called before the phase is + * initialzed with elements and/or species. + * + * For this phase, the chemical potential is set + * + * @param eosdata An XML_Node object corresponding to + * the "thermo" entry for this phase in the input file. + * + * eosdata points to the thermo block, and looks like this: + * + * @verbatim + + + -2.7E7 + + @endverbatim + * + */ + virtual void setParametersFromXML(const XML_Node& eosdata); + + + //! Function to set the chemical potential directly + /*! + * @param chemPot Value of the chemical potential (units J/kmol) + */ + void setChemicalPotential(doublereal chemPot); + + protected: + + + double chemPot_; + + }; + + +} + +#endif diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index 73669357e..b359d5579 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -80,10 +80,10 @@ endif ifeq ($(do_issp),1) ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o \ GibbsExcessVPSSTP.o PseudoBinaryVPSSTP.o MargulesVPSSTP.o \ - IonsFromNeutralVPSSTP.o PDSS_IonsFromNeutral.o + IonsFromNeutralVPSSTP.o PDSS_IonsFromNeutral.o FixedChemPotSSTP.o ISSP_H = IdealSolidSolnPhase.h StoichSubstanceSSTP.h SingleSpeciesTP.h MineralEQ3.h \ GibbsExcessVPSSTP.h PseudoBinaryVPSSTP.h MargulesVPSSTP.h \ - IonsFromNeutralVPSSTP.h PDSS_IonsFromNeutral.h + IonsFromNeutralVPSSTP.h PDSS_IonsFromNeutral.h FixedChemPotSSTP.h endif CATHERMO_OBJ = $(THERMO_OBJ) $(ELECTRO_OBJ) $(ISSP_OBJ) diff --git a/Cantera/src/thermo/SingleSpeciesTP.cpp b/Cantera/src/thermo/SingleSpeciesTP.cpp index 8e8372c47..7860bb859 100644 --- a/Cantera/src/thermo/SingleSpeciesTP.cpp +++ b/Cantera/src/thermo/SingleSpeciesTP.cpp @@ -346,7 +346,7 @@ namespace Cantera { * This member function is resolved here. A single species phase obtains its * thermo from the standard state function. * - * @param cpbar On return, Contains the molar volume of the single species + * @param vbar On return, Contains the molar volume of the single species * and the phase. Units are m^3 / kmol. Length = 1 */ void SingleSpeciesTP::getPartialMolarVolumes(doublereal* vbar) const { @@ -514,7 +514,11 @@ namespace Cantera { void SingleSpeciesTP::setState_UV(doublereal u, doublereal v, doublereal tol) { doublereal dt; - setDensity(1.0/v); + if (v == 0.0) { + setDensity(1.0E100); + } else { + setDensity(1.0/v); + } for (int n = 0; n < 50; n++) { dt = (u - intEnergy_mass())/cv_mass(); if (dt > 100.0) dt = 100.0; @@ -548,7 +552,11 @@ namespace Cantera { void SingleSpeciesTP::setState_SV(doublereal s, doublereal v, doublereal tol) { doublereal dt; - setDensity(1.0/v); + if (v == 0.0) { + setDensity(1.0E100); + } else { + setDensity(1.0/v); + } for (int n = 0; n < 50; n++) { dt = (s - entropy_mass())*temperature()/cv_mass(); if (dt > 100.0) dt = 100.0; diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index cf04e5e00..ea91a0cd7 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -51,7 +51,6 @@ #ifdef WITH_STOICH_SUBSTANCE #ifdef USE_SSTP #include "StoichSubstanceSSTP.h" - #else #include "StoichSubstance.h" #endif @@ -60,6 +59,7 @@ #ifdef WITH_STOICH_SUBSTANCE #include "MineralEQ3.h" #include "MetalSHEelectrons.h" +#include "FixedChemPotSSTP.h" #endif //#include "importCTML.h" @@ -91,7 +91,7 @@ namespace Cantera { /*! * @deprecated This entire structure could be replaced with a std::map */ - static int ntypes = 18; + static int ntypes = 19; //! Define the string name of the %ThermoPhase types that are handled by this factory routine static string _types[] = {"IdealGas", "Incompressible", @@ -100,7 +100,7 @@ namespace Cantera { "HMW", "IdealSolidSolution", "DebyeHuckel", "IdealMolalSolution", "IdealGasVPSS", "MineralEQ3", "MetalSHEelectrons", "Margules", - "IonsFromNeutralMolecule" + "IonsFromNeutralMolecule", "FixedChemPot" }; //! Define the integer id of the %ThermoPhase types that are handled by this factory routine @@ -110,7 +110,7 @@ namespace Cantera { cHMW, cIdealSolidSolnPhase, cDebyeHuckel, cIdealMolalSoln, cVPSS_IdealGas, cMineralEQ3, cMetalSHEelectrons, - cMargulesVPSSTP, cIonsFromNeutral + cMargulesVPSSTP, cIonsFromNeutral, cFixedChemPot }; /* @@ -173,6 +173,12 @@ namespace Cantera { break; #endif +#ifdef WITH_STOICH_SUBSTANCE + case cFixedChemPot: + th = new FixedChemPotSSTP; + break; +#endif + #ifdef WITH_STOICH_SUBSTANCE case cMineralEQ3: th = new MineralEQ3(); diff --git a/Cantera/src/thermo/mix_defs.h b/Cantera/src/thermo/mix_defs.h index 50482d169..02ab93eb9 100644 --- a/Cantera/src/thermo/mix_defs.h +++ b/Cantera/src/thermo/mix_defs.h @@ -55,6 +55,9 @@ namespace Cantera { /// An edge between two 2D surfaces const int cEdge = 6; + //! Stoichiometric compound with a constant chemical potential + const int cFixedChemPot = 70; + /// Constant partial molar volume solution IdealSolidSolnPhase.h const int cIdealSolidSolnPhase = 5009;