From 07bcac001b0e8b82d775d91025fd3bfeade2c676 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Mon, 23 Feb 2009 20:41:03 +0000 Subject: [PATCH] Started coding up general excess gibbs free energy thermo formulations based on mole fraction expansions. --- Cantera/src/thermo/Constituents.h | 1 - Cantera/src/thermo/GibbsExcessVPSSTP.cpp | 315 +++++++++++++++++ Cantera/src/thermo/GibbsExcessVPSSTP.h | 410 +++++++++++++++++++++++ Cantera/src/thermo/Makefile.in | 6 +- 4 files changed, 729 insertions(+), 3 deletions(-) create mode 100644 Cantera/src/thermo/GibbsExcessVPSSTP.cpp create mode 100644 Cantera/src/thermo/GibbsExcessVPSSTP.h diff --git a/Cantera/src/thermo/Constituents.h b/Cantera/src/thermo/Constituents.h index d5d5fe58d..a38622906 100755 --- a/Cantera/src/thermo/Constituents.h +++ b/Cantera/src/thermo/Constituents.h @@ -147,7 +147,6 @@ namespace Cantera { /// Return a read-only reference to the vector of atomic weights. const vector_fp& atomicWeights() const; - /// Number of elements. int nElements() const; diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp new file mode 100644 index 000000000..d2ee007f1 --- /dev/null +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp @@ -0,0 +1,315 @@ +/** + * @file GibbsExcessVPSSTP.cpp + * Definitions for intermediate ThermoPhase object for phases which + * employ excess gibbs free energy formulations + * (see \ref thermoprops + * and class \link Cantera::GibbsExcessVPSSTP GibbsExcessVPSSTP\endlink). + * + * Header file for a derived class of ThermoPhase that handles + * variable pressure standard state methods for calculating + * thermodynamic properties that are further based upon expressions + * for the excess gibbs free energy expressed as a function of + * the mole fractions + */ +/* + * Copywrite (2009) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Date$ + * $Revision$ + */ + + +#include "GibbsExcessVPSSTP.h" +using namespace std; + +namespace Cantera { + + /* + * Default constructor. + * + */ + GibbsExcessVPSSTP::GibbsExcessVPSSTP() : + VPStandardStateTP() + { + } + + /* + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + GibbsExcessVPSSTP::GibbsExcessVPSSTP(const GibbsExcessVPSSTP &b) : + VPStandardStateTP() + { + *this = operator=(b); + } + + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + GibbsExcessVPSSTP& GibbsExcessVPSSTP:: + operator=(const GibbsExcessVPSSTP &b) { + if (&b != this) { + VPStandardStateTP::operator=(b); + } + return *this; + } + + /** + * + * ~GibbsExcessVPSSTP(): (virtual) + * + * Destructor: does nothing: + * + */ + GibbsExcessVPSSTP::~GibbsExcessVPSSTP() { + } + + /* + * This routine duplicates the current object and returns + * a pointer to ThermoPhase. + */ + ThermoPhase* + GibbsExcessVPSSTP::duplMyselfAsThermoPhase() const { + GibbsExcessVPSSTP* mtp = new GibbsExcessVPSSTP(*this); + return (ThermoPhase *) mtp; + } + + /* + * -------------- Utilities ------------------------------- + */ + + // Equation of state type flag. + /* + * The ThermoPhase base class returns + * zero. Subclasses should define this to return a unique + * non-zero value. Known constants defined for this purpose are + * listed in mix_defs.h. The GibbsExcessVPSSTP class also returns + * zero, as it is a non-complete class. + */ + int GibbsExcessVPSSTP::eosType() const { + return 0; + } + + + + /* + * ------------ Molar Thermodynamic Properties ---------------------- + */ + + + /* + * - Activities, Standard States, Activity Concentrations ----------- + */ + + + doublereal GibbsExcessVPSSTP::standardConcentration(int k) const { + err("standardConcentration"); + return -1.0; + } + + doublereal GibbsExcessVPSSTP::logStandardConc(int k) const { + err("logStandardConc"); + return -1.0; + } + + void GibbsExcessVPSSTP::getActivities(doublereal* ac) const { + getActivityCoefficients(ac); + getMoleFractions(DATA_PTR(moleFractions_)); + for (int k = 0; k < m_kk; k++) { + ac[k] *= moleFractions_[k]; + } + } + + + void GibbsExcessVPSSTP::getElectrochemPotentials(doublereal* mu) const { + getChemPotentials(mu); + double ve = Faraday * electricPotential(); + for (int k = 0; k < m_kk; k++) { + mu[k] += ve*charge(k); + } + } + + /* + * ------------ Partial Molar Properties of the Solution ------------ + */ + + + doublereal GibbsExcessVPSSTP::err(std::string msg) const { + throw CanteraError("GibbsExcessVPSSTP","Base class method " + +msg+" called. Equation of state type: "+int2str(eosType())); + return 0; + } + + /* + * Returns the units of the standard and general concentrations + * Note they have the same units, as their divisor 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. + * + * On return uA contains the powers of the units (MKS assumed) + * of the standard concentrations and generalized concentrations + * for the kth species. + * + * 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 GibbsExcessVPSSTP::getUnitsStandardConc(double *uA, int k, int sizeUA) const { + for (int i = 0; i < sizeUA; i++) { + if (i == 0) uA[0] = 1.0; + if (i == 1) uA[1] = -nDim(); + if (i == 2) uA[2] = 0.0; + if (i == 3) uA[3] = 0.0; + if (i == 4) uA[4] = 0.0; + if (i == 5) uA[5] = 0.0; + } + } + + + /* + * @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 GibbsExcessVPSSTP::initThermo() { + initLengths(); + VPStandardStateTP::initThermo(); + + } + + + // Initialize lengths of local variables after all species have + // been identified. + void GibbsExcessVPSSTP::initLengths() { + m_kk = nSpecies(); + moleFractions_.resize(m_kk); + } + + /* + * initThermoXML() (virtual from ThermoPhase) + * Import and initialize a ThermoPhase object + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void GibbsExcessVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { + + initLengths(); + + VPStandardStateTP::initThermoXML(phaseNode, id); + } + + /** + * Format a summary of the mixture state for output. + */ + std::string GibbsExcessVPSSTP::report(bool show_thermo) const { + + + char p[800]; + string s = ""; + try { + if (name() != "") { + sprintf(p, " \n %s:\n", name().c_str()); + s += p; + } + sprintf(p, " \n temperature %12.6g K\n", temperature()); + s += p; + sprintf(p, " pressure %12.6g Pa\n", pressure()); + s += p; + sprintf(p, " density %12.6g kg/m^3\n", density()); + s += p; + sprintf(p, " mean mol. weight %12.6g amu\n", meanMolecularWeight()); + s += p; + + doublereal phi = electricPotential(); + sprintf(p, " potential %12.6g V\n", phi); + s += p; + + int kk = nSpecies(); + array_fp x(kk); + array_fp molal(kk); + array_fp mu(kk); + array_fp muss(kk); + array_fp acMolal(kk); + array_fp actMolal(kk); + getMoleFractions(&x[0]); + + getChemPotentials(&mu[0]); + getStandardChemPotentials(&muss[0]); + getActivities(&actMolal[0]); + + + if (show_thermo) { + sprintf(p, " \n"); + s += p; + sprintf(p, " 1 kg 1 kmol\n"); + s += p; + sprintf(p, " ----------- ------------\n"); + s += p; + sprintf(p, " enthalpy %12.6g %12.4g J\n", + enthalpy_mass(), enthalpy_mole()); + s += p; + sprintf(p, " internal energy %12.6g %12.4g J\n", + intEnergy_mass(), intEnergy_mole()); + s += p; + sprintf(p, " entropy %12.6g %12.4g J/K\n", + entropy_mass(), entropy_mole()); + s += p; + sprintf(p, " Gibbs function %12.6g %12.4g J\n", + gibbs_mass(), gibbs_mole()); + s += p; + sprintf(p, " heat capacity c_p %12.6g %12.4g J/K\n", + cp_mass(), cp_mole()); + s += p; + try { + sprintf(p, " heat capacity c_v %12.6g %12.4g J/K\n", + cv_mass(), cv_mole()); + s += p; + } + catch(CanteraError) { + sprintf(p, " heat capacity c_v \n"); + s += p; + } + } + + } catch (CanteraError) { + ; + } + return s; + } + + +} + diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.h b/Cantera/src/thermo/GibbsExcessVPSSTP.h new file mode 100644 index 000000000..57db5781f --- /dev/null +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.h @@ -0,0 +1,410 @@ +/** + * @file gibbsExcessVPSSTP.h + * Header for intermediate ThermoPhase object for phases which + * employ gibbs excess free energy based formulations + * (see \ref thermoprops + * and class \link Cantera::gibbsExcessVPSSTP gibbsExcessVPSSTP\endlink). + * + * Header file for a derived class of ThermoPhase that handles + * variable pressure standard state methods for calculating + * thermodynamic properties that are further based upon activities + * based on the molality scale. These include most of the methods for + * calculating liquid electrolyte thermodynamics. + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#ifndef CT_GIBBSEXCESSVPSSTP_H +#define CT_GIBBSEXCESSVPSSTP_H + +#include "VPStandardStateTP.h" + +namespace Cantera { + + /** + * @ingroup thermoprops + */ + + /*! + * GibbsExcessVPSSTP is a derived class of ThermoPhase that handles + * variable pressure standard state methods for calculating + * thermodynamic properties that are further based on + * expressing the Excess Gibbs free energy as a function of + * the mole fractions (or pseudo mole fractions) of consitituents. + * This category is the workhorse for describing molten salts, + * solid-phase mixtures of semiconductors, and mixtures of miscible + * and semi-miscible compounds. + * + * It includes + * . regular solutions + * . Margueles expansions + * . NTRL equation + * . Wilson's equation + * . UNIQUAC equation of state. + * + * This class adds additional functions onto the %ThermoPhase interface + * that handles the calculation of the excess Gibbs free energy. The %ThermoPhase + * class includes a member function, ThermoPhase::activityConvention() + * that indicates which convention the activities are based on. The + * default is to assume activities are based on the molar convention. + * That default is used here. + * + * All of the Excess Gibbs free energy formulations in this area employ + * symmetrical formulations. + * + * + * Chemical potentials + * of species k, \f$ \mu_o \f$, has the following general format: + * + * \f[ + * \mu_k = \mu^o_k(T,P) + R T ln( \gamma_k X_k ) + * \f] + * + * + * where \f$ \gamma_k^{\triangle} \f$ is a molar based activity coefficient for species + * \f$k\f$. + * + * GibbsExcessVPSSTP contains an internal vector with the current mole + * fraction vector. That's one of its primary usages. + * + *

SetState Strategy

+ * + * The gibbsExcessVPSSTP object does not have a setState strategy. + * It's strictly an interfacial layer that writes the current mole fractions to the + * State object. + * + * + */ + class GibbsExcessVPSSTP : public VPStandardStateTP { + + public: + + /// Constructors + /*! + * This doesn't do much more than initialize constants with + * default values for water at 25C. Water molecular weight + * comes from the default elements.xml file. It actually + * differs slightly from the IAPWS95 value of 18.015268. However, + * density conservation and therefore element conservation + * is the more important principle to follow. + */ + GibbsExcessVPSSTP(); + + //! Copy constructor + /*! + * Note this stuff will not work until the underlying phase + * has a working copy constructor + * + * @param b class to be copied + */ + GibbsExcessVPSSTP(const GibbsExcessVPSSTP &b); + + /// Assignment operator + /*! + * Note this stuff will not work until the underlying phase + * has a working assignment operator + * + * @param b class to be copied. + */ + GibbsExcessVPSSTP& operator=(const GibbsExcessVPSSTP &b); + + /// Destructor. + virtual ~GibbsExcessVPSSTP(); + + //! Duplication routine for objects which inherit from ThermoPhase. + /*! + * This virtual routine can be used to duplicate thermophase objects + * inherited from ThermoPhase even if the application only has + * a pointer to ThermoPhase to work with. + */ + virtual ThermoPhase *duplMyselfAsThermoPhase() const; + + /** + * + * @name Utilities + * @{ + */ + + + //! Equation of state type flag. + /*! + * The ThermoPhase base class returns + * zero. Subclasses should define this to return a unique + * non-zero value. Known constants defined for this purpose are + * listed in mix_defs.h. The MolalityVPSSTP class also returns + * zero, as it is a non-complete class. + */ + virtual int eosType() const; + + + + /** + * @} + * @name Molar Thermodynamic Properties + * @{ + */ + + + /** + * @} + * @name Utilities for Solvent ID and Molality + * @{ + */ + + + + + /** + * @} + * @name Mechanical Properties + * @{ + */ + + /** + * @} + * @name Potential Energy + * + * Species may have an additional potential energy due to the + * presence of external gravitation or electric fields. These + * methods allow specifying a potential energy for individual + * species. + * @{ + */ + + /** + * @} + * @name Activities, Standard States, and Activity Concentrations + * + * The activity \f$a_k\f$ of a species in solution is + * related to the chemical potential by \f[ \mu_k = \mu_k^0(T) + * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is + * the chemical potential at unit activity, which depends only + * on temperature and pressure. + * @{ + */ + + + + + /** + * The standard concentration \f$ C^0_k \f$ used to normalize + * the generalized concentration. In many cases, this quantity + * will be the same for all species in a phase - for example, + * for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this + * reason, this method returns a single value, instead of an + * array. However, for phases in which the standard + * concentration is species-specific (e.g. surface species of + * different sizes), this method may be called with an + * optional parameter indicating the species. + * + * @param k species index. Defaults to zero. + */ + virtual doublereal standardConcentration(int k=0) const; + + /** + * Returns the natural logarithm of the standard + * concentration of the kth species + * + * @param k species index + */ + virtual doublereal logStandardConc(int k=0) 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. + * + * @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(double *uA, int k = 0, + int sizeUA = 6) const; + + + //! Get the array of non-dimensional activities (molality + //! based for this class and classes that derive from it) at + //! the current solution temperature, pressure, and solution concentration. + /*! + * All standard state properties for molality-based phases are + * evaluated consistent with the molality scale. Therefore, this function + * must return molality-based activities. + * + * \f[ + * a_i^\triangle = \gamma_k^{\triangle} \frac{m_k}{m^\triangle} + * \f] + * + * This function must be implemented in derived classes. + * + * @param ac Output vector of molality-based activities. Length: m_kk. + */ + virtual void getActivities(doublereal* ac) const; + + + + //@} + /// @name Partial Molar Properties of the Solution + //@{ + + + /** + * Get the species electrochemical potentials. + * These are partial molar quantities. + * This method adds a term \f$ Fz_k \phi_k \f$ to the + * to each chemical potential. + * + * Units: J/kmol + * + * @param mu output vector containing the species electrochemical potentials. + * Length: m_kk. + */ + void getElectrochemPotentials(doublereal* mu) const; + + + //@} + /// @name Properties of the Standard State of the Species in the Solution + //@{ + + + + //@} + /// @name Thermodynamic Values for the Species Reference States + //@{ + + + /////////////////////////////////////////////////////// + // + // The methods below are not virtual, and should not + // be overloaded. + // + ////////////////////////////////////////////////////// + + /** + * @name Specific Properties + * @{ + */ + + + /** + * @name Setting the State + * + * These methods set all or part of the thermodynamic + * state. + * @{ + */ + + //@} + + /** + * @name Chemical Equilibrium + * Routines that implement the Chemical equilibrium capability + * for a single phase, based on the element-potential method. + * @{ + */ + + + //@} + + + + /// The following methods are used in the process of constructing + /// the phase and setting its parameters from a specification in an + /// input file. They are not normally used in application programs. + /// To see how they are used, see files importCTML.cpp and + /// ThermoFactory.cpp. + + + /*! + * @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(); + + + /** + * Import and initialize a ThermoPhase object + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void initThermoXML(XML_Node& phaseNode, std::string id); + + + //! returns a summary of the state of the phase as a string + /*! + * @param show_thermo If true, extra information is printed out + * about the thermodynamic state of the system. + */ + virtual std::string report(bool show_thermo = true) const; + + + private: + + + //! Initialize lengths of local variables after all species have + //! been identified. + void initLengths(); + + + + private: + //! Error function + /*! + * Print an error string and exit + * + * @param msg Message to be printed + */ + doublereal err(std::string msg) const; + + private: + + //! Storage for the current values of the mole fractions of the species + mutable std::vector moleFractions_; + + }; + + +} + +#endif + + + + + diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index 993e57594..1e618236a 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -78,8 +78,10 @@ ELECTRO_H = MolalityVPSSTP.h VPStandardStateTP.h \ VPSSMgr_Water_ConstVol.h endif ifeq ($(do_issp),1) -ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o -ISSP_H = IdealSolidSolnPhase.h StoichSubstanceSSTP.h SingleSpeciesTP.h MineralEQ3.h +ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o \ + GibbsExcessVPSSTP.o +ISSP_H = IdealSolidSolnPhase.h StoichSubstanceSSTP.h SingleSpeciesTP.h MineralEQ3.h \ + GibbsExcessVPSSTP.h endif CATHERMO_OBJ = $(THERMO_OBJ) $(ELECTRO_OBJ) $(ISSP_OBJ)