From 49b413e61e2ce21cd575bf18391c129eedbe349e Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 3 Mar 2009 21:08:31 +0000 Subject: [PATCH] Added an initial attempt at a Margules excess gibbs free energy model. It's just hard coded to match the Eutectic of LiKCl at the moment. --- Cantera/src/thermo/GibbsExcessVPSSTP.cpp | 83 ++- Cantera/src/thermo/GibbsExcessVPSSTP.h | 71 ++- Cantera/src/thermo/Makefile.in | 4 +- Cantera/src/thermo/MargulesVPSSTP.cpp | 586 ++++++++++++++++++++++ Cantera/src/thermo/MargulesVPSSTP.h | 424 ++++++++++++++++ Cantera/src/thermo/PseudoBinaryVPSSTP.cpp | 368 ++++++++++++++ Cantera/src/thermo/PseudoBinaryVPSSTP.h | 392 +++++++++++++++ Cantera/src/thermo/VPSSMgrFactory.cpp | 65 ++- Cantera/src/thermo/VPSSMgr_ConstVol.cpp | 4 +- Cantera/src/thermo/VPStandardStateTP.h | 13 +- 10 files changed, 1975 insertions(+), 35 deletions(-) create mode 100644 Cantera/src/thermo/MargulesVPSSTP.cpp create mode 100644 Cantera/src/thermo/MargulesVPSSTP.h create mode 100644 Cantera/src/thermo/PseudoBinaryVPSSTP.cpp create mode 100644 Cantera/src/thermo/PseudoBinaryVPSSTP.h diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp index 9a947fc67..131572b10 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp @@ -59,7 +59,11 @@ namespace Cantera { if (&b != this) { VPStandardStateTP::operator=(b); } + moleFractions_ = b.moleFractions_; + lnActCoeff_Scaled_ = b.lnActCoeff_Scaled_; + m_pp = b.m_pp; + return *this; } @@ -144,6 +148,52 @@ namespace Cantera { * ------------ Molar Thermodynamic Properties ---------------------- */ + /* + * + * ------------ Mechanical Properties ------------------------------ + * + */ + + /* + * Set the pressure at constant temperature. Units: Pa. + * This method sets a constant within the object. + * The mass density is not a function of pressure. + */ + void GibbsExcessVPSSTP::setPressure(doublereal p) { +#ifdef DEBUG_MODE + //printf("setPressure: %g\n", p); +#endif + /* + * Store the current pressure + */ + m_Pcurrent = p; + /* + * update the standard state thermo + * -> This involves calling the water function and setting the pressure + */ + updateStandardStateThermo(); + + /* + * Calculate all of the other standard volumes + * -> note these are constant for now + */ + calcDensity(); + } + + void GibbsExcessVPSSTP::calcDensity() { + double *vbar = &m_pp[0]; + getPartialMolarVolumes(vbar); + + doublereal vtotal = 0.0; + for (int i = 0; i < m_kk; i++) { + vtotal += vbar[i] * moleFractions_[i]; + } + doublereal dd = meanMolecularWeight() / vtotal; + State::setDensity(dd); + } + + + /* * - Activities, Standard States, Activity Concentrations ----------- @@ -181,6 +231,24 @@ namespace Cantera { * ------------ Partial Molar Properties of the Solution ------------ */ + // Return an array of partial molar volumes for the + // species in the mixture. Units: m^3/kmol. + /* + * Frequently, for this class of thermodynamics representations, + * the excess Volume due to mixing is zero. Here, we set it as + * a default. It may be overriden in derived classes. + * + * @param vbar Output vector of speciar partial molar volumes. + * Length = m_kk. units are m^3/kmol. + */ + void GibbsExcessVPSSTP::getPartialMolarVolumes(doublereal* vbar) const { + /* + * Get the standard state values in m^3 kmol-1 + */ + getStandardVolumes(vbar); + } + + doublereal GibbsExcessVPSSTP::err(std::string msg) const { throw CanteraError("GibbsExcessVPSSTP","Base class method " @@ -188,6 +256,16 @@ namespace Cantera { return 0; } + //@} + /// @name Properties of the Standard State of the Species in the Solution + //@{ + + + + //@} + /// @name Thermodynamic Values for the Species Reference States + //@{ + double GibbsExcessVPSSTP::checkMFSum(const doublereal * const x) const { doublereal norm = accumulate(x, x + m_kk, 0.0); if (fabs(norm - 1.0) > 1.0E-9) { @@ -256,6 +334,8 @@ namespace Cantera { void GibbsExcessVPSSTP::initLengths() { m_kk = nSpecies(); moleFractions_.resize(m_kk); + lnActCoeff_Scaled_.resize(m_kk); + m_pp.resize(m_kk); } /* @@ -275,8 +355,7 @@ namespace Cantera { */ void GibbsExcessVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { - initLengths(); - + VPStandardStateTP::initThermoXML(phaseNode, id); } diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.h b/Cantera/src/thermo/GibbsExcessVPSSTP.h index 3309352f6..f37deb271 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.h +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.h @@ -107,8 +107,6 @@ namespace Cantera { /// Assignment operator /*! - * Note this stuff will not work until the underlying phase - * has a working assignment operator * * @param b class to be copied. */ @@ -151,13 +149,7 @@ namespace Cantera { */ - /** - * @} - * @name Utilities for Solvent ID and Molality - * @{ - */ - - + /** @@ -166,6 +158,46 @@ namespace Cantera { * @{ */ + //! Set the internally storred pressure (Pa) at constant + //! temperature and composition + /*! + * This method sets the pressure within the object. + * The water model is a completely compressible model. + * Also, the dielectric constant is pressure dependent. + * + * @param p input Pressure (Pa) + * + * @todo Implement a variable pressure capability + */ + virtual void setPressure(doublereal p); + + private: + /** + * 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. + */ + void calcDensity(); + + public: /** * @} * @name Potential Energy @@ -279,7 +311,18 @@ namespace Cantera { */ void getElectrochemPotentials(doublereal* mu) const; - + //! Return an array of partial molar volumes for the + //! species in the mixture. Units: m^3/kmol. + /*! + * Frequently, for this class of thermodynamics representations, + * the excess Volume due to mixing is zero. Here, we set it as + * a default. It may be overriden in derived classes. + * + * @param vbar Output vector of speciar partial molar volumes. + * Length = m_kk. units are m^3/kmol. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + //@} /// @name Properties of the Standard State of the Species in the Solution //@{ @@ -458,11 +501,17 @@ namespace Cantera { double checkMFSum(const doublereal * const x) const; - private: + protected: //! Storage for the current values of the mole fractions of the species mutable std::vector moleFractions_; + //! Storage for the current values of the activity coefficients of the + //! species, divided by RT + mutable std::vector lnActCoeff_Scaled_; + + mutable std::vector m_pp; + }; diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index 1e618236a..7be968840 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -79,9 +79,9 @@ ELECTRO_H = MolalityVPSSTP.h VPStandardStateTP.h \ endif ifeq ($(do_issp),1) ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o \ - GibbsExcessVPSSTP.o + GibbsExcessVPSSTP.o PseudoBinaryVPSSTP.o MargulesVPSSTP.o ISSP_H = IdealSolidSolnPhase.h StoichSubstanceSSTP.h SingleSpeciesTP.h MineralEQ3.h \ - GibbsExcessVPSSTP.h + GibbsExcessVPSSTP.h PseudoBinaryVPSSTP.h MargulesVPSSTP.h endif CATHERMO_OBJ = $(THERMO_OBJ) $(ELECTRO_OBJ) $(ISSP_OBJ) diff --git a/Cantera/src/thermo/MargulesVPSSTP.cpp b/Cantera/src/thermo/MargulesVPSSTP.cpp new file mode 100644 index 000000000..e338eda48 --- /dev/null +++ b/Cantera/src/thermo/MargulesVPSSTP.cpp @@ -0,0 +1,586 @@ +/** + * @file MargulesVPSSTP.cpp + * Definitions for ThermoPhase object for phases which + * employ excess gibbs free energy formulations related to Margules + * expansions (see \ref thermoprops + * and class \link Cantera::MargulesVPSSTP MargulesVPSSTP\endlink). + * + */ +/* + * 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 "MargulesVPSSTP.h" +#include "ThermoFactory.h" + + +using namespace std; + +namespace Cantera { + + static const double xxSmall = 1.0E-150; + /* + * Default constructor. + * + */ + MargulesVPSSTP::MargulesVPSSTP() : + PseudoBinaryVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + } + + /* + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + MargulesVPSSTP::MargulesVPSSTP(const MargulesVPSSTP &b) : + PseudoBinaryVPSSTP() + { + *this = operator=(b); + } + + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + MargulesVPSSTP& MargulesVPSSTP:: + operator=(const MargulesVPSSTP &b) { + if (&b != this) { + PseudoBinaryVPSSTP::operator=(b); + } + + numBinaryInteractions_ = b.numBinaryInteractions_ ; + m_HE_b_ij = b.m_HE_b_ij; + m_HE_c_ij = b.m_HE_c_ij; + m_HE_d_ij = b.m_HE_d_ij; + m_SE_b_ij = b.m_SE_b_ij; + m_SE_c_ij = b.m_SE_c_ij; + m_SE_d_ij = b.m_SE_d_ij; + m_pSpecies_A_ij = b.m_pSpecies_A_ij; + m_pSpecies_B_ij = b.m_pSpecies_B_ij; + formMargules_ = b.formMargules_; + formTempModel_ = b.formTempModel_; + + return *this; + } + + /** + * + * ~MargulesVPSSTP(): (virtual) + * + * Destructor: does nothing: + * + */ + MargulesVPSSTP::~MargulesVPSSTP() { + } + + /* + * This routine duplicates the current object and returns + * a pointer to ThermoPhase. + */ + ThermoPhase* + MargulesVPSSTP::duplMyselfAsThermoPhase() const { + MargulesVPSSTP* mtp = new MargulesVPSSTP(*this); + return (ThermoPhase *) mtp; + } + + // Special constructor for a hard-coded problem + /* + * + * LiKCl treating the PseudoBinary layer as passthrough. + * -> test to predict the eutectic and liquidus correctly. + * + */ + MargulesVPSSTP::MargulesVPSSTP(int testProb) : + PseudoBinaryVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + + + constructPhaseFile("LiKCl_liquid.xml", ""); + + + numBinaryInteractions_ = 1; + + m_HE_b_ij.resize(1); + m_HE_c_ij.resize(1); + m_HE_d_ij.resize(1); + + m_SE_b_ij.resize(1); + m_SE_c_ij.resize(1); + m_SE_d_ij.resize(1); + + m_pSpecies_A_ij.resize(1); + m_pSpecies_B_ij.resize(1); + + + m_HE_b_ij[0] = -17570E3; + m_HE_c_ij[0] = -377.0E3; + m_HE_d_ij[0] = 0.0; + + m_SE_b_ij[0] = -7.627E3; + m_SE_c_ij[0] = 4.958E3; + m_SE_d_ij[0] = 0.0; + + int iLiCl = speciesIndex("LiCl(L)"); + if (iLiCl < 0) { + throw CanteraError("MargulesVPSSTP test1 constructor", + "Unable to find LiCl(L)"); + } + m_pSpecies_B_ij[0] = iLiCl; + + + int iKCl = speciesIndex("KCl(L)"); + if (iKCl < 0) { + throw CanteraError("MargulesVPSSTP test1 constructor", + "Unable to find KCl(L)"); + } + m_pSpecies_A_ij[0] = iKCl; + } + + + /* + * -------------- 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 MargulesVPSSTP class also returns + * zero, as it is a non-complete class. + */ + int MargulesVPSSTP::eosType() const { + return 0; + } + + /* + * Import, construct, and initialize a phase + * specification from an XML tree into the current object. + * + * This routine is a precursor to constructPhaseXML(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void MargulesVPSSTP::constructPhaseFile(std::string inputFile, std::string id) { + + if (inputFile.size() == 0) { + throw CanteraError("MargulesVPSSTP:constructPhaseFile", + "input file is null"); + } + string path = findInputFile(inputFile); + std::ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("MargulesVPSSTP:constructPhaseFile","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + XML_Node &phaseNode_XML = xml(); + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("MargulesVPSSTP:constructPhaseFile", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + fxml_phase->copy(&phaseNode_XML); + constructPhaseXML(*fxml_phase, id); + delete fxml; + } + + /* + * Import, construct, and initialize a HMWSoln phase + * specification from an XML tree into the current object. + * + * Most of the work is carried out by the cantera base + * routine, importPhase(). That routine imports all of the + * species and element data, including the standard states + * of the species. + * + * Then, In this routine, we read the information + * particular to the specification of the activity + * coefficient model for the Pitzer parameterization. + * + * We also read information about the molar volumes of the + * standard states if present in the XML file. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void MargulesVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id) { + string stemp; + if (id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("MargulesVPSSTP::constructPhaseXML", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("MargulesVPSSTP::constructPhaseXML", + "no thermo XML node"); + } + XML_Node& thermoNode = phaseNode.child("thermo"); + + /* + * Possibly change the form of the standard concentrations + */ + + /* + * Get the Name of the Solvent: + * solventName + */ + string solventName = ""; + if (thermoNode.hasChild("solvent")) { + XML_Node& scNode = thermoNode.child("solvent"); + vector nameSolventa; + getStringArray(scNode, nameSolventa); + int nsp = static_cast(nameSolventa.size()); + if (nsp != 1) { + throw CanteraError("MargulesVPSSTP::constructPhaseXML", + "badly formed solvent XML node"); + } + solventName = nameSolventa[0]; + } + + /* + * Determine the form of the Pitzer model, + * We will use this information to size arrays below. + */ + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& scNode = thermoNode.child("activityCoefficients"); + + stemp = scNode.attrib("model"); + string formString = lowercase(stemp); + if (formString != "") { + if (formString == "margules" || formString == "default") { + formMargules_ = 0; + + } else { + throw CanteraError("MargulesVPSSTP::constructPhaseXML", + "Unknown ActivityCoeff model: " + + formString); + } + } + + /* + * Determine the form of the temperature dependence + * of the Pitzer activity coefficient model. + */ + stemp = scNode.attrib("TempModel"); + formString = lowercase(stemp); + if (formString != "") { + if (formString == "constant" || formString == "default") { + formTempModel_ = 0; + } else { + throw CanteraError("MargulesVPSSTP::constructPhaseXML", + "Unknown Pitzer ActivityCoeff Temp model: " + + formString); + } + } + + + } + + /* + * Call the Cantera importPhase() function. This will import + * all of the species into the phase. This will also handle + * all of the solvent and solute standard states + */ + bool m_ok = importPhase(phaseNode, this); + if (!m_ok) { + throw CanteraError("MargulesVPSSTP::constructPhaseXML","importPhase failed "); + } + + } + + + /* + * ------------ Molar Thermodynamic Properties ---------------------- + */ + + + /* + * - Activities, Standard States, Activity Concentrations ----------- + */ + + + doublereal MargulesVPSSTP::standardConcentration(int k) const { + err("standardConcentration"); + return -1.0; + } + + doublereal MargulesVPSSTP::logStandardConc(int k) const { + err("logStandardConc"); + return -1.0; + } + + // Get the array of non-dimensional molar-based activity coefficients at + // the current solution temperature, pressure, and solution concentration. + /* + * @param ac Output vector of activity coefficients. Length: m_kk. + */ + void MargulesVPSSTP::getActivityCoefficients(doublereal* ac) const { + /* + * Update the activity coefficients + */ + s_update_lnActCoeff(); + + /* + * take the exp of the internally storred coefficients. + */ + for (int k = 0; k < m_kk; k++) { + ac[k] = exp(lnActCoeff_Scaled_[k]); + } + } + + + void MargulesVPSSTP::getElectrochemPotentials(doublereal* mu) const { + getChemPotentials(mu); + double ve = Faraday * electricPotential(); + for (int k = 0; k < m_kk; k++) { + mu[k] += ve*charge(k); + } + } + + + void MargulesVPSSTP::getChemPotentials(doublereal* mu) const { + doublereal xx; + /* + * First get the standard chemical potentials in + * molar form. + * -> this requires updates of standard state as a function + * of T and P + */ + getStandardChemPotentials(mu); + /* + * Update the activity coefficients + */ + s_update_lnActCoeff(); + /* + * + */ + doublereal RT = GasConstant * temperature(); + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(moleFractions_[k], xxSmall); + mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]); + } + } + + + + + + /* + * ------------ Partial Molar Properties of the Solution ------------ + */ + + + doublereal MargulesVPSSTP::err(std::string msg) const { + throw CanteraError("MargulesVPSSTP","Base class method " + +msg+" called. Equation of state type: "+int2str(eosType())); + return 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 MargulesVPSSTP::initThermo() { + initLengths(); + PseudoBinaryVPSSTP::initThermo(); + } + + + // Initialize lengths of local variables after all species have + // been identified. + void MargulesVPSSTP::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 MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { + + PseudoBinaryVPSSTP::initThermoXML(phaseNode, id); + } + + // Update the activity coefficients + /* + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + * + * he = XAXB(B + C(XA - XB) + d ( X_A X_B) + */ + void MargulesVPSSTP::s_update_lnActCoeff() const { + int iA, iB; + double XA, XB, g0 , g1; + double T = temperature(); + + for (int i = 0; i < m_kk; i++) { + lnActCoeff_Scaled_[i] = 0.0; + } + double RT = GasConstant * temperature(); + for (int i = 0; i < numBinaryInteractions_; i++) { + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT ; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + + lnActCoeff_Scaled_[iA] += XB * XB * (g0 + g1 * (XB - XA)); + lnActCoeff_Scaled_[iB] += XA * XA * g0 + XA * XB * g1 * (2 * XA); + + } + + + } + + /** + * Format a summary of the mixture state for output. + */ + std::string MargulesVPSSTP::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/MargulesVPSSTP.h b/Cantera/src/thermo/MargulesVPSSTP.h new file mode 100644 index 000000000..f9d9cb55e --- /dev/null +++ b/Cantera/src/thermo/MargulesVPSSTP.h @@ -0,0 +1,424 @@ +/** + * @file Margules.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_MARGULESVPSSTP_H +#define CT_MARGULESVPSSTP_H + +#include "PseudoBinaryVPSSTP.h" + +namespace Cantera { + + /** + * @ingroup thermoprops + */ + + /*! + * MargulesVPSSTP is a derived class of PseudoBinaryVPSSTP. + * + * + * + */ + class MargulesVPSSTP : public PseudoBinaryVPSSTP { + + public: + + //! Constructor + /*! + * 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. + */ + MargulesVPSSTP(); + + //! Special constructor for a hard-coded problem + /*! + * + * LiKCl treating the PseudoBinary layer as passthrough. + * -> test to predict the eutectic and liquidus correctly. + * + */ + MargulesVPSSTP(int testProb); + + //! Copy constructor + /*! + * Note this stuff will not work until the underlying phase + * has a working copy constructor + * + * @param b class to be copied + */ + MargulesVPSSTP(const MargulesVPSSTP&b); + + //! Assignment operator + /*! + * + * @param b class to be copied. + */ + MargulesVPSSTP& operator=(const MargulesVPSSTP &b); + + //! Destructor + virtual ~MargulesVPSSTP(); + + //! 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; + + //! Initialization of a phase using an xml file + /*! + * This routine is a precursor to + * routine, which does most of the work. + * + * @param inputFile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void constructPhaseFile(std::string inputFile, std::string id); + + //! Import and initialize a phase + //! specification in an XML tree into the current object. + /*! + * Here we read an XML description of the phase. + * We import descriptions of the elements that make up the + * species in a phase. + * We import information about the species, including their + * reference state thermodynamic polynomials. We then freeze + * the state of the species. + * + * Then, we read the species molar volumes from the xml + * tree to finish the initialization. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void constructPhaseXML(XML_Node& phaseNode, std::string id); + + /** + * @} + * @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; + + //! Get the array of non-dimensional molar-based activity coefficients at + //! the current solution temperature, pressure, and solution concentration. + /*! + * @param ac Output vector of activity coefficients. Length: m_kk. + */ + virtual void getActivityCoefficients(doublereal* ac) const; + + + + + //@} + /// @name Partial Molar Properties of the Solution + //@{ + + //! Get the species chemical potentials. Units: J/kmol. + /*! + * This function returns a vector of chemical potentials of the + * species in solution at the current temperature, pressure + * and mole fraction of the solution. + * + * @param mu Output vector of species chemical + * potentials. Length: m_kk. Units: J/kmol + */ + virtual void getChemPotentials(doublereal* mu) const; + + + //! Get the species electrochemical potentials. + /*! + * 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(); + + //! Update the activity coefficients + /*! + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + */ + void s_update_lnActCoeff() const; + + + private: + //! Error function + /*! + * Print an error string and exit + * + * @param msg Message to be printed + */ + doublereal err(std::string msg) const; + + protected: + + + //! number of binary interaction expressions + + int numBinaryInteractions_; + + mutable vector_fp m_HE_b_ij; + + mutable vector_fp m_HE_c_ij; + + mutable vector_fp m_HE_d_ij; + + + mutable vector_fp m_SE_b_ij; + + mutable vector_fp m_SE_c_ij; + + mutable vector_fp m_SE_d_ij; + + vector_int m_pSpecies_A_ij; + vector_int m_pSpecies_B_ij; + + + int formMargules_; + int formTempModel_; + + private: + + + }; + + + +} + +#endif + + + + + diff --git a/Cantera/src/thermo/PseudoBinaryVPSSTP.cpp b/Cantera/src/thermo/PseudoBinaryVPSSTP.cpp new file mode 100644 index 000000000..c369facfc --- /dev/null +++ b/Cantera/src/thermo/PseudoBinaryVPSSTP.cpp @@ -0,0 +1,368 @@ +/** + * @file PseudoBinaryVPSSTP.cpp + * Definitions for intermediate ThermoPhase object for phases which + * employ excess gibbs free energy formulations + * (see \ref thermoprops + * and class \link Cantera::PseudoBinaryVPSSTP PseudoBinaryVPSSTP\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 "PseudoBinaryVPSSTP.h" +using namespace std; + +namespace Cantera { + + /* + * Default constructor. + * + */ + PseudoBinaryVPSSTP::PseudoBinaryVPSSTP() : + GibbsExcessVPSSTP(), + PBType_(PBTYPE_PASSTHROUGH), + numPBSpecies_(m_kk), + indexSpecialSpecies_(-1), + numCationSpecies_(0), + numAnionSpecies_(0), + numPassThroughSpecies_(0), + neutralPBindexStart(0), + cationPhase_(0), + anionPhase_(0) + { + } + + /* + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + PseudoBinaryVPSSTP::PseudoBinaryVPSSTP(const PseudoBinaryVPSSTP &b) : + GibbsExcessVPSSTP(), + PBType_(PBTYPE_PASSTHROUGH), + numPBSpecies_(m_kk), + indexSpecialSpecies_(-1), + numCationSpecies_(0), + numAnionSpecies_(0), + numPassThroughSpecies_(0), + neutralPBindexStart(0), + cationPhase_(0), + anionPhase_(0) + { + *this = operator=(b); + } + + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + PseudoBinaryVPSSTP& PseudoBinaryVPSSTP:: + operator=(const PseudoBinaryVPSSTP &b) { + if (&b != this) { + GibbsExcessVPSSTP::operator=(b); + } + + PBType_ = b.PBType_; + numPBSpecies_ = b.numPBSpecies_; + indexSpecialSpecies_ = b.indexSpecialSpecies_; + PBMoleFractions_ = b.PBMoleFractions_; + cationList_ = b.cationList_; + numCationSpecies_ = b.numCationSpecies_; + anionList_ = b.anionList_; + numAnionSpecies_ = b.numAnionSpecies_; + passThroughList_ = b.passThroughList_; + numPassThroughSpecies_ = b.numPassThroughSpecies_; + neutralPBindexStart = b.neutralPBindexStart; + cationPhase_ = b.cationPhase_; + anionPhase_ = b.anionPhase_; + moleFractionsTmp_ = b.moleFractionsTmp_; + + return *this; + } + + /** + * + * ~PseudoBinaryVPSSTP(): (virtual) + * + * Destructor: does nothing: + * + */ + PseudoBinaryVPSSTP::~PseudoBinaryVPSSTP() { + } + + /* + * This routine duplicates the current object and returns + * a pointer to ThermoPhase. + */ + ThermoPhase* + PseudoBinaryVPSSTP::duplMyselfAsThermoPhase() const { + PseudoBinaryVPSSTP* mtp = new PseudoBinaryVPSSTP(*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 PseudoBinaryVPSSTP class also returns + * zero, as it is a non-complete class. + */ + int PseudoBinaryVPSSTP::eosType() const { + return 0; + } + + + + /* + * ------------ Molar Thermodynamic Properties ---------------------- + */ + + + /* + * - Activities, Standard States, Activity Concentrations ----------- + */ + + + doublereal PseudoBinaryVPSSTP::standardConcentration(int k) const { + err("standardConcentration"); + return -1.0; + } + + doublereal PseudoBinaryVPSSTP::logStandardConc(int k) const { + err("logStandardConc"); + return -1.0; + } + + + + + void PseudoBinaryVPSSTP::getElectrochemPotentials(doublereal* mu) const { + getChemPotentials(mu); + double ve = Faraday * electricPotential(); + for (int k = 0; k < m_kk; k++) { + mu[k] += ve*charge(k); + } + } + + void PseudoBinaryVPSSTP::calcPseudoBinaryMoleFractions() const { + int k; + doublereal sumCat; + doublereal sumAnion; + doublereal sum = 0.0; + switch (PBType_) { + case PBTYPE_PASSTHROUGH: + for (k = 0; k < m_kk; k++) { + PBMoleFractions_[k] = moleFractions_[k]; + } + break; + case PBTYPE_SINGLEANION: + sumCat = 0.0; + sumAnion = 0.0; + for (k = 0; k < m_kk; k++) { + moleFractionsTmp_[k] = moleFractions_[k]; + } + for (k = 0; k < (int) cationList_.size(); k++) { + sumCat += moleFractions_[cationList_[k]]; + } + sumAnion = moleFractions_[anionList_[k]]; + PBMoleFractions_[0] = sumCat -sumAnion; + moleFractionsTmp_[indexSpecialSpecies_] -= PBMoleFractions_[0]; + + + for (k = 0; k < numCationSpecies_; k++) { + PBMoleFractions_[1+k] = moleFractionsTmp_[cationList_[k]]; + } + + for (k = 0; k < numPassThroughSpecies_; k++) { + PBMoleFractions_[neutralPBindexStart + k] = + moleFractions_[cationList_[k]]; + } + + sum = fmax(0.0, PBMoleFractions_[0]); + for (k = 1; k < numPBSpecies_; k++) { + sum += PBMoleFractions_[k]; + } + for (k = 0; k < numPBSpecies_; k++) { + PBMoleFractions_[k] /= sum; + } + + break; + case PBTYPE_SINGLECATION: + throw CanteraError("eosType", "Unknown type"); + + break; + + case PBTYPE_MULTICATIONANION: + throw CanteraError("eosType", "Unknown type"); + + break; + default: + throw CanteraError("eosType", "Unknown type"); + break; + + } + } + + /* + * ------------ Partial Molar Properties of the Solution ------------ + */ + + + doublereal PseudoBinaryVPSSTP::err(std::string msg) const { + throw CanteraError("PseudoBinaryVPSSTP","Base class method " + +msg+" called. Equation of state type: "+int2str(eosType())); + return 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 PseudoBinaryVPSSTP::initThermo() { + initLengths(); + GibbsExcessVPSSTP::initThermo(); + } + + + // Initialize lengths of local variables after all species have + // been identified. + void PseudoBinaryVPSSTP::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 PseudoBinaryVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { + + + GibbsExcessVPSSTP::initThermoXML(phaseNode, id); + } + + /** + * Format a summary of the mixture state for output. + */ + std::string PseudoBinaryVPSSTP::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/PseudoBinaryVPSSTP.h b/Cantera/src/thermo/PseudoBinaryVPSSTP.h new file mode 100644 index 000000000..3cc9a9d5b --- /dev/null +++ b/Cantera/src/thermo/PseudoBinaryVPSSTP.h @@ -0,0 +1,392 @@ +/** + * @file PseudoBinaryVPSSTP.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_PSEUDOBINARYVPSSTP_H +#define CT_PSEUDOBINARYVPSSTP_H + +#include "GibbsExcessVPSSTP.h" + +namespace Cantera { + + /** + * @ingroup thermoprops + */ + + /*! + * PseudoBinaryVPSSTP is a derived class of ThermoPhase + * GibbsExcessVPSSTP 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. + * + * This layer will massage the mole fraction vector to implement + * cation and anion based mole numbers in an optional manner + * + * The way that it collects the cation and anion based mole numbers + * is via holding two extra ThermoPhase objects. These + * can include standard states for salts. + * + * + */ + class PseudoBinaryVPSSTP : public GibbsExcessVPSSTP { + + 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. + */ + PseudoBinaryVPSSTP(); + + //! Copy constructor + /*! + * Note this stuff will not work until the underlying phase + * has a working copy constructor + * + * @param b class to be copied + */ + PseudoBinaryVPSSTP(const PseudoBinaryVPSSTP&b); + + /// Assignment operator + /*! + * + * @param b class to be copied. + */ + PseudoBinaryVPSSTP& operator=(const PseudoBinaryVPSSTP&b); + + /// Destructor. + virtual ~PseudoBinaryVPSSTP(); + + //! 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; + + + + + //@} + /// @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. + * @{ + */ + + //! Calculate pseudo binary mole fractions + /*! + * + */ + virtual void calcPseudoBinaryMoleFractions() const; + + + //@} + + /** + * @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; + + protected: + + int PBType_; + + //! Number of pseudo binary species + int numPBSpecies_; + + //! index of special species + int indexSpecialSpecies_; + + mutable std::vector PBMoleFractions_; + + std::vector cationList_; + int numCationSpecies_; + + std::vectoranionList_; + int numAnionSpecies_; + + std::vector passThroughList_; + int numPassThroughSpecies_; + int neutralPBindexStart; + + ThermoPhase *cationPhase_; + + ThermoPhase *anionPhase_; + + mutable std::vector moleFractionsTmp_; + + private: + + + }; + +#define PBTYPE_PASSTHROUGH 0 +#define PBTYPE_SINGLEANION 1 +#define PBTYPE_SINGLECATION 2 +#define PBTYPE_MULTICATIONANION 3 + + + +} + +#endif + + + + + diff --git a/Cantera/src/thermo/VPSSMgrFactory.cpp b/Cantera/src/thermo/VPSSMgrFactory.cpp index 9460a3225..a5a006e50 100644 --- a/Cantera/src/thermo/VPSSMgrFactory.cpp +++ b/Cantera/src/thermo/VPSSMgrFactory.cpp @@ -68,17 +68,23 @@ namespace Cantera { * @param spData_node Species Data XML node. This node contains a list * of species XML nodes underneath it. * - * @todo Make sure that spDadta_node is species Data XML node by checking its name is speciesData + * @todo Make sure that spDadta_node is species Data XML node by checking + * its name is speciesData */ static void getVPSSMgrTypes(std::vector & spDataNodeList, - int &has_nasa, - int& has_shomate, - int& has_simple, + int &has_nasa_idealGas, + int &has_nasa_constVol, + int& has_shomate_idealGas, + int& has_shomate_constVol, + int& has_simple_idealGas, + int& has_simple_constVol, int &has_water, int &has_tpx, int &has_hptx, int &has_other) { + XML_Node *ss_ptr = 0; + string ssModel = "idealGas"; size_t ns = spDataNodeList.size(); for (size_t n = 0; n < ns; n++) { bool ifound = false; @@ -98,21 +104,49 @@ namespace Cantera { if (!ifound) { if (spNode->hasChild("thermo")) { const XML_Node& th = spNode->child("thermo"); + if (spNode->hasChild("standardState")) { + ss_ptr = &(spNode->child("standardState")); + ssModel = ss_ptr->attrib("model"); + } if (th.hasChild("NASA")) { - has_nasa++; + if (ssModel == "idealGas") { + has_nasa_idealGas++; + } else if (ssModel == "constant_incompressible" || + ssModel == "constantVolume") { + has_nasa_constVol++; + } else { + throw UnknownVPSSMgrModel("getVPSSMgrTypes:", + spNode->attrib("name")); + } ifound = true; } if (th.hasChild("Shomate")) { - has_shomate++; + if (ssModel == "idealGas") { + has_shomate_idealGas++; + } else if (ssModel == "constant_incompressible" || + ssModel == "constantVolume") { + has_shomate_constVol++; + } else { + throw UnknownVPSSMgrModel("getVPSSMgrTypes:", + spNode->attrib("name")); + } ifound = true; } if (th.hasChild("const_cp")){ - has_simple = 1; + if (ssModel == "idealGas") { + has_simple_idealGas++; + } else if (ssModel == "constant_incompressible" || + ssModel == "constantVolume") { + has_simple_constVol++; + } else { + throw UnknownVPSSMgrModel("getVPSSMgrTypes:", + spNode->attrib("name")); + } ifound = true; } if (th.hasChild("poly")) { if (th.child("poly")["order"] == "1") { - has_simple = 1; + has_simple_constVol = 1; ifound = true; } else throw CanteraError("newSpeciesThermo", "poly with order > 1 not yet supported"); @@ -252,12 +286,14 @@ namespace Cantera { } - int inasa = 0, ishomate = 0, isimple = 0, iwater = 0, itpx = 0, iother = 0; + int inasaIG = 0, inasaCV = 0, ishomateIG = 0, ishomateCV = 0, + isimpleIG = 0, isimpleCV = 0, + iwater = 0, itpx = 0, iother = 0; int ihptx = 0; try { - getVPSSMgrTypes(spDataNodeList, inasa, ishomate, isimple, iwater, - itpx, ihptx, iother); + getVPSSMgrTypes(spDataNodeList, inasaIG, inasaCV, ishomateIG, ishomateCV, + isimpleIG, isimpleCV, iwater, itpx, ihptx, iother); } catch (UnknownSpeciesThermoModel) { iother = 1; popError(); @@ -270,6 +306,13 @@ namespace Cantera { vpss = new VPSSMgr_Water_HKFT(vp_ptr, spth); } } + if (vpss == 0) { + if (inasaCV || ishomateCV || isimpleCV) { + if (!inasaIG && !ishomateIG && !isimpleIG && !itpx && !ihptx && !iother) { + vpss = new VPSSMgr_ConstVol(vp_ptr, spth); + } + } + } if (vpss == 0) { vpss = new VPSSMgr_General(vp_ptr, spth); } diff --git a/Cantera/src/thermo/VPSSMgr_ConstVol.cpp b/Cantera/src/thermo/VPSSMgr_ConstVol.cpp index 4f725e081..bba01572a 100644 --- a/Cantera/src/thermo/VPSSMgr_ConstVol.cpp +++ b/Cantera/src/thermo/VPSSMgr_ConstVol.cpp @@ -112,7 +112,7 @@ namespace Cantera { "no standardState Node for species " + s->name()); } std::string model = (*ss)["model"]; - if (model != "constant_incompressible") { + if (model != "constant_incompressible" && model != "constantVolume") { throw CanteraError("VPSSMgr_ConstVol::initThermoXML", "standardState model for species isn't constant_incompressible: " + s->name()); } @@ -135,7 +135,7 @@ namespace Cantera { "no standardState Node for species " + speciesNode.name()); } std::string model = (*ss)["model"]; - if (model != "constant_incompressible") { + if (model != "constant_incompressible" && model != "constantVolume") { throw CanteraError("VPSSMgr_ConstVol::initThermoXML", "standardState model for species isn't " "constant_incompressible: " + speciesNode.name()); diff --git a/Cantera/src/thermo/VPStandardStateTP.h b/Cantera/src/thermo/VPStandardStateTP.h index 8d5f872f2..5c381b8ec 100644 --- a/Cantera/src/thermo/VPStandardStateTP.h +++ b/Cantera/src/thermo/VPStandardStateTP.h @@ -1,12 +1,12 @@ /** * @file VPStandardStateTP.h - * Header file for a derived class of ThermoPhase that handles - * variable pressure standard state methods for calculating - * thermodynamic properties (see \ref thermoprops and - * class \link Cantera::VPStandardStateTP VPStandardStateTP\endlink). + * Header file for a derived class of ThermoPhase that handles + * variable pressure standard state methods for calculating + * thermodynamic properties (see \ref thermoprops and + * class \link Cantera::VPStandardStateTP VPStandardStateTP\endlink). * - * These include most of the - * methods for calculating liquid electrolyte thermodynamics. + * These include most of the + * methods for calculating liquid electrolyte thermodynamics. */ /* * Copywrite (2005) Sandia Corporation. Under the terms of @@ -14,7 +14,6 @@ * U.S. Government retains certain rights in this software. */ /* - * $Author$ * $Date$ * $Revision$ */