From 06276fba6daa02687c749d59e06a0856a6711e3f Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 2 Sep 2009 17:17:52 +0000 Subject: [PATCH] Added a new thermo object in. IonsFromNeutralVPSSTP : This is initially used for molten salts. Molten salt thermo is treated via the neutral molecule molten salt expressions. However, ions migrate by themselves. therefore, you must treat ions within the transport equations. This thermophase object treats the ions as species, and is made to be consistent with the neutral molecules thermophase object that is used to express the non-ideal thermodynamics of the molten salt solution. Currently, MargulesVPSSTP is used for this as the LiKCl system is parameterized using this excess Gibbs free energy formulation. MargulesVPSSTP Various upgrades and fillers. It now reads and is initiated from the input file. setState_TP(), setTemperature(), and setPressure() within VPStandardStateTP() formulations There was some circular logic, previously. Now setTemperature() and setPressure() are subserviant to setState_TP(). Basically, it doesn't make sense to calculate the standard state without first getting both the temperature and pressure (not just one and then the other) before you calculate the standard state for VPStandardStateTP derivatives. Think about the water EOS, which is a real fluid. You don't want to be trying to calculate the standard state within the enveloppe where neither gas nor liquid is stable. PseudoBinaryVPSSTP Deprecated. --- Cantera/src/thermo/DebyeHuckel.cpp | 16 +- Cantera/src/thermo/DebyeHuckel.h | 9 + Cantera/src/thermo/GibbsExcessVPSSTP.cpp | 36 +- Cantera/src/thermo/GibbsExcessVPSSTP.h | 12 +- Cantera/src/thermo/HMWSoln.cpp | 54 +- Cantera/src/thermo/HMWSoln.h | 10 + Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp | 1139 ++++++++++++++++++ Cantera/src/thermo/IonsFromNeutralVPSSTP.h | 695 +++++++++++ Cantera/src/thermo/Makefile.in | 10 +- Cantera/src/thermo/MargulesVPSSTP.cpp | 285 +++-- Cantera/src/thermo/MargulesVPSSTP.h | 305 ++++- Cantera/src/thermo/PDSS_IonsFromNeutral.cpp | 484 ++++++++ Cantera/src/thermo/PDSS_IonsFromNeutral.h | 462 +++++++ Cantera/src/thermo/PseudoBinaryVPSSTP.cpp | 1 - Cantera/src/thermo/ThermoFactory.cpp | 12 +- Cantera/src/thermo/ThermoPhase.cpp | 1 - Cantera/src/thermo/VPSSMgr_ConstVol.cpp | 38 + Cantera/src/thermo/VPSSMgr_ConstVol.h | 23 + Cantera/src/thermo/VPSSMgr_General.cpp | 17 +- Cantera/src/thermo/VPStandardStateTP.cpp | 16 +- Cantera/src/thermo/VPStandardStateTP.h | 22 +- Cantera/src/thermo/mix_defs.h | 5 +- 22 files changed, 3501 insertions(+), 151 deletions(-) create mode 100644 Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp create mode 100644 Cantera/src/thermo/IonsFromNeutralVPSSTP.h create mode 100644 Cantera/src/thermo/PDSS_IonsFromNeutral.cpp create mode 100644 Cantera/src/thermo/PDSS_IonsFromNeutral.h diff --git a/Cantera/src/thermo/DebyeHuckel.cpp b/Cantera/src/thermo/DebyeHuckel.cpp index 8d2be1136..f5b5aeb3d 100644 --- a/Cantera/src/thermo/DebyeHuckel.cpp +++ b/Cantera/src/thermo/DebyeHuckel.cpp @@ -312,16 +312,13 @@ namespace Cantera { return m_Pcurrent; } - /* - * 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 DebyeHuckel::setPressure(doublereal p) { + setState_TP(temperature(), p); + } -#ifdef DEBUG_MODE - //printf("setPressure: %g\n", p); -#endif + void DebyeHuckel::setState_TP(doublereal t, doublereal p) { + + State::setTemperature(t); /* * Store the current pressure */ @@ -458,8 +455,7 @@ namespace Cantera { * the value propagates to underlying objects. */ void DebyeHuckel::setTemperature(const doublereal temp) { - _updateStandardStateThermo(); - State::setTemperature(temp); + setState_TP(temp, m_Pcurrent); } diff --git a/Cantera/src/thermo/DebyeHuckel.h b/Cantera/src/thermo/DebyeHuckel.h index 3b8285a2a..6c9f2f2ee 100644 --- a/Cantera/src/thermo/DebyeHuckel.h +++ b/Cantera/src/thermo/DebyeHuckel.h @@ -846,6 +846,15 @@ namespace Cantera { */ virtual void setTemperature(const doublereal temp); + //! Set the temperature (K) and pressure (Pa) + /*! + * Set the temperature and pressure. + * + * @param t Temperature (K) + * @param p Pressure (Pa) + */ + virtual void setState_TP(doublereal t, doublereal p); + /** * The isothermal compressibility. Units: 1/Pa. * The isothermal compressibility is defined as diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp index ad23a911d..d58daaac0 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp @@ -160,24 +160,7 @@ namespace Cantera { * 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(); + setState_TP(temperature(), p); } void GibbsExcessVPSSTP::calcDensity() { @@ -192,6 +175,23 @@ namespace Cantera { State::setDensity(dd); } + void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p) { + State::setTemperature(t); + /* + * Store the current pressure + */ + m_Pcurrent = p; + /* + * update the standard state thermo + * -> This involves calling the water function and setting the pressure + */ + updateStandardStateThermo(); + + /* + * Calculate the partial molar volumes, and then the density of the fluid + */ + calcDensity(); + } diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.h b/Cantera/src/thermo/GibbsExcessVPSSTP.h index d405d3fba..e8a33cce1 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.h +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.h @@ -171,7 +171,7 @@ namespace Cantera { */ virtual void setPressure(doublereal p); - private: + protected: /** * Calculate the density of the mixture using the partial * molar volumes and mole fractions as input @@ -356,6 +356,14 @@ namespace Cantera { */ + //! Set the temperature (K) and pressure (Pa) + /*! + * Set the temperature and pressure. + * + * @param t Temperature (K) + * @param p Pressure (Pa) + */ + virtual void setState_TP(doublereal t, doublereal p); //@} @@ -402,6 +410,7 @@ namespace Cantera { * Length is m_kk. */ virtual void setMoleFractions(const doublereal* const x); + /** * Set the mole fractions to the specified values without * normalizing. This is useful when the normalization @@ -510,6 +519,7 @@ namespace Cantera { //! species, divided by RT mutable std::vector lnActCoeff_Scaled_; + //! Temporary storage space that is fair game mutable std::vector m_pp; }; diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp index d97522f64..c221ab0cf 100644 --- a/Cantera/src/thermo/HMWSoln.cpp +++ b/Cantera/src/thermo/HMWSoln.cpp @@ -763,30 +763,7 @@ namespace Cantera { * The mass density is not a function of pressure. */ void HMWSoln::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(); - - /* - * Store the internal density of the water SS. - * Note, we would have to do this for all other - * species if they had pressure dependent properties. - */ - m_densWaterSS = m_waterSS->density(); - /* - * Calculate all of the other standard volumes - * -> note these are constant for now - */ - calcDensity(); + setState_TP(temperature(), p); } void HMWSoln::calcDensity() { @@ -889,9 +866,36 @@ namespace Cantera { * the value propagates to underlying objects. */ void HMWSoln::setTemperature(const doublereal temp) { + setState_TP(temp, m_Pcurrent); + } + + /* + * Overwritten setTemperature(double) from State.h. This + * function sets the temperature, and makes sure that + * the value propagates to underlying objects. + */ + void HMWSoln::setState_TP(doublereal temp, doublereal pres) { State::setTemperature(temp); - //m_waterSS->setTemperature(temp); + /* + * Store the current pressure + */ + m_Pcurrent = pres; + + /* + * update the standard state thermo + * -> This involves calling the water function and setting the pressure + */ updateStandardStateThermo(); + /* + * Store the internal density of the water SS. + * Note, we would have to do this for all other + * species if they had pressure dependent properties. + */ + m_densWaterSS = m_waterSS->density(); + /* + * Calculate all of the other standard volumes + * -> note these are constant for now + */ calcDensity(); } diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index e6b82a2f2..ff9b5c0a5 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -1533,6 +1533,16 @@ namespace Cantera { */ virtual void setTemperature(const doublereal temp); + //! Set the temperature (K) and pressure (Pa) + /*! + * Set the temperature and pressure. + * + * @param t Temperature (K) + * @param p Pressure (Pa) + */ + virtual void setState_TP(doublereal t, doublereal p); + + /** * The isothermal compressibility. Units: 1/Pa. * The isothermal compressibility is defined as diff --git a/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp new file mode 100644 index 000000000..365ab3354 --- /dev/null +++ b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -0,0 +1,1139 @@ +/** + * @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: 2009/03/27 01:08:55 $ + * $Revision: 1.2 $ + */ + + +#include "IonsFromNeutralVPSSTP.h" +#include "ThermoFactory.h" + +#include "PDSS_IonsFromNeutral.h" +#include "mix_defs.h" + +#include + +using namespace std; + +#ifndef MIN +# define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) +#endif + +namespace Cantera { + + /* + * Default constructor. + * + */ + IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP() : + GibbsExcessVPSSTP(), + ionSolnType_(cIonSolnType_SINGLEANION), + numNeutralMoleculeSpecies_(0), + indexSpecialSpecies_(-1), + indexSecondSpecialSpecies_(-1), + numCationSpecies_(0), + numAnionSpecies_(0), + numPassThroughSpecies_(0), + neutralMoleculePhase_(0), + cationPhase_(0), + anionPhase_(0) + { + } + /* + * Working constructors + * + * The two constructors below are the normal way + * the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + + */ + IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(std::string inputFile, std::string id, + ThermoPhase *neutralPhase) : + GibbsExcessVPSSTP(), + ionSolnType_(cIonSolnType_SINGLEANION), + numNeutralMoleculeSpecies_(0), + indexSpecialSpecies_(-1), + indexSecondSpecialSpecies_(-1), + numCationSpecies_(0), + numAnionSpecies_(0), + numPassThroughSpecies_(0), + neutralMoleculePhase_(neutralPhase), + cationPhase_(0), + anionPhase_(0) + { + constructPhaseFile(inputFile, id); + } + + IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(XML_Node& phaseRoot, std::string id, + ThermoPhase *neutralPhase) : + GibbsExcessVPSSTP(), + ionSolnType_(cIonSolnType_SINGLEANION), + numNeutralMoleculeSpecies_(0), + indexSpecialSpecies_(-1), + indexSecondSpecialSpecies_(-1), + numCationSpecies_(0), + numAnionSpecies_(0), + numPassThroughSpecies_(0), + neutralMoleculePhase_(neutralPhase), + cationPhase_(0), + anionPhase_(0) + { + constructPhaseXML(phaseRoot, id); + } + + + + /* + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const IonsFromNeutralVPSSTP &b) : + GibbsExcessVPSSTP(), + ionSolnType_(cIonSolnType_SINGLEANION), + numNeutralMoleculeSpecies_(0), + indexSpecialSpecies_(-1), + indexSecondSpecialSpecies_(-1), + numCationSpecies_(0), + numAnionSpecies_(0), + numPassThroughSpecies_(0), + neutralMoleculePhase_(0), + cationPhase_(0), + anionPhase_(0) + { + *this = operator=(b); + } + + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + IonsFromNeutralVPSSTP& IonsFromNeutralVPSSTP:: + operator=(const IonsFromNeutralVPSSTP &b) { + if (&b != this) { + GibbsExcessVPSSTP::operator=(b); + } + + ionSolnType_ = b.ionSolnType_; + numNeutralMoleculeSpecies_ = b.numNeutralMoleculeSpecies_; + indexSpecialSpecies_ = b.indexSpecialSpecies_; + indexSecondSpecialSpecies_ = b.indexSecondSpecialSpecies_; + NeutralMolecMoleFractions_ = b.NeutralMolecMoleFractions_; + cationList_ = b.cationList_; + numCationSpecies_ = b.numCationSpecies_; + anionList_ = b.anionList_; + numAnionSpecies_ = b.numAnionSpecies_; + passThroughList_ = b.passThroughList_; + numPassThroughSpecies_ = b.numPassThroughSpecies_; + /* + * This is a shallow copy. We need to figure this out + */ + neutralMoleculePhase_ = b.neutralMoleculePhase_; + if (neutralMoleculePhase_) { + exit(-1); + } + + cationPhase_ = b.cationPhase_; + anionPhase_ = b.anionPhase_; + moleFractionsTmp_ = b.moleFractionsTmp_; + + return *this; + } + + /** + * + * ~IonsFromNeutralVPSSTP(): (virtual) + * + * Destructor: does nothing: + * + */ + IonsFromNeutralVPSSTP::~IonsFromNeutralVPSSTP() { + } + + /* + * This routine duplicates the current object and returns + * a pointer to ThermoPhase. + */ + ThermoPhase* + IonsFromNeutralVPSSTP::duplMyselfAsThermoPhase() const { + IonsFromNeutralVPSSTP* mtp = new IonsFromNeutralVPSSTP(*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 IonsFromNeutralVPSSTP class also returns + * zero, as it is a non-complete class. + */ + int IonsFromNeutralVPSSTP::eosType() const { + return cIonFromNeutral; + } + + + + /* + * ------------ Molar Thermodynamic Properties ---------------------- + */ + /* + * Molar enthalpy of the solution. Units: J/kmol. + */ + doublereal IonsFromNeutralVPSSTP::enthalpy_mole() const { + getPartialMolarEnthalpies(DATA_PTR(m_pp)); + return mean_X(DATA_PTR(m_pp)); + } + + /** + * Molar internal energy of the solution. Units: J/kmol. + * + * This is calculated from the soln enthalpy and then + * subtracting pV. + */ + doublereal IonsFromNeutralVPSSTP::intEnergy_mole() const { + double hh = enthalpy_mole(); + double pres = pressure(); + double molarV = 1.0/molarDensity(); + double uu = hh - pres * molarV; + return uu; + } + + /** + * Molar soln entropy at constant pressure. Units: J/kmol/K. + * + * This is calculated from the partial molar entropies. + */ + doublereal IonsFromNeutralVPSSTP::entropy_mole() const { + getPartialMolarEntropies(DATA_PTR(m_pp)); + return mean_X(DATA_PTR(m_pp)); + } + + /// Molar Gibbs function. Units: J/kmol. + doublereal IonsFromNeutralVPSSTP::gibbs_mole() const { + getChemPotentials(DATA_PTR(m_pp)); + return mean_X(DATA_PTR(m_pp)); + } + /** Molar heat capacity at constant pressure. Units: J/kmol/K. + * + * Returns the solution heat capacition at constant pressure. + * This is calculated from the partial molar heat capacities. + */ + doublereal IonsFromNeutralVPSSTP::cp_mole() const { + getPartialMolarCp(DATA_PTR(m_pp)); + double val = mean_X(DATA_PTR(m_pp)); + return val; + } + + /// Molar heat capacity at constant volume. Units: J/kmol/K. + doublereal IonsFromNeutralVPSSTP::cv_mole() const { + // Need to revisit this, as it is wrong + getPartialMolarCp(DATA_PTR(m_pp)); + return mean_X(DATA_PTR(m_pp)); + //err("not implemented"); + //return 0.0; + } + + /* + * - Activities, Standard States, Activity Concentrations ----------- + */ + + // This method returns an array of generalized concentrations + /* + * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / + * C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration + * defined below and \f$ a_k \f$ are activities used in the + * thermodynamic functions. These activity (or generalized) + * concentrations are used + * by kinetics manager classes to compute the forward and + * reverse rates of elementary reactions. Note that they may + * or may not have units of concentration --- they might be + * partial pressures, mole fractions, or surface coverages, + * for example. + * + * @param c Output array of generalized concentrations. The + * units depend upon the implementation of the + * reaction rate expressions within the phase. + */ + void IonsFromNeutralVPSSTP::getActivityConcentrations(doublereal* c) const { + err(""); + } + + // 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 IonsFromNeutralVPSSTP::getActivityCoefficients(doublereal* ac) const { + + // This stuff has moved to the setState routines + // calcNeutralMoleculeMoleFractions(); + // neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_)); + // neutralMoleculePhase_->getStandardChemPotentials(DATA_PTR(muNeutralMolecule_)); + + /* + * 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]); + } + } + + + doublereal IonsFromNeutralVPSSTP::standardConcentration(int k) const { + err("standardConcentration"); + return -1.0; + } + + doublereal IonsFromNeutralVPSSTP::logStandardConc(int k) const { + err("logStandardConc"); + return -1.0; + } + + /* + * --------- 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 + */ + void + IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const { + int k, icat, jNeut; + doublereal xx, fact2; + /* + * Transfer the mole fractions to the slave neutral molecule + * phase + * Note we may move this in the future. + */ + //calcNeutralMoleculeMoleFractions(); + //neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_)); + + /* + * Get the standard chemical potentials of netural molecules + */ + neutralMoleculePhase_->getStandardChemPotentials(DATA_PTR(muNeutralMolecule_)); + + doublereal RT_ = GasConstant * temperature(); + + switch (ionSolnType_) { + case cIonSolnType_PASSTHROUGH: + neutralMoleculePhase_->getChemPotentials(mu); + break; + case cIonSolnType_SINGLEANION: + neutralMoleculePhase_->getActivityCoefficients(DATA_PTR(gammaNeutralMolecule_)); + + fact2 = 2.0 * RT_ * log(2.0); + + // Do the cation list + for (k = 0; k < (int) cationList_.size(); k++) { + //! Get the id for the next cation + icat = cationList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + xx = fmaxx(SmallNumber, moleFractions_[icat]); + mu[icat] = muNeutralMolecule_[jNeut] + fact2 + RT_ * log(gammaNeutralMolecule_[jNeut] * xx); + } + + // Do the anion list + icat = anionList_[0]; + jNeut = fm_invert_ionForNeutral[icat]; + xx = fmaxx(SmallNumber, moleFractions_[icat]); + mu[icat] = RT_ * log(xx); + + // Do the list of neutral molecules + for (k = 0; k < numPassThroughSpecies_; k++) { + icat = passThroughList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + xx = fmaxx(SmallNumber, moleFractions_[icat]); + mu[icat] = muNeutralMolecule_[jNeut] + RT_ * log( gammaNeutralMolecule_[jNeut] * xx); + } + break; + + case cIonSolnType_SINGLECATION: + throw CanteraError("eosType", "Unknown type"); + break; + case cIonSolnType_MULTICATIONANION: + throw CanteraError("eosType", "Unknown type"); + break; + default: + throw CanteraError("eosType", "Unknown type"); + break; + } + } + + // This is temporary. We will get rid of this + void IonsFromNeutralVPSSTP::setTemperature(doublereal t) { + double p = pressure(); + IonsFromNeutralVPSSTP::setState_TP(t, p); + } + + // This is temporary. We will get rid of this + void IonsFromNeutralVPSSTP::setPressure(doublereal p) { + double t = temperature(); + IonsFromNeutralVPSSTP::setState_TP(t, p); + } + + // Set the temperature (K) and pressure (Pa) + /* + * Setting the pressure may involve the solution of a nonlinear equation. + * + * @param t Temperature (K) + * @param p Pressure (Pa) + */ + void IonsFromNeutralVPSSTP::setState_TP(doublereal t, doublereal p) { + /* + * This is a two phase process. First, we calculate the standard states + * within the neutral molecule phase. + */ + neutralMoleculePhase_->setState_TP(t, p); + VPStandardStateTP::setState_TP(t,p); + + /* + * Calculate the partial molar volumes, and then the density of the fluid + */ + calcDensity(); + } + + // Calculate ion mole fractions from neutral molecule + // mole fractions. + /* + * @param mf Dump the mole fractions into this vector. + */ + void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal * const mf) const { + int k; + doublereal fmij; + /* + * Download the neutral mole fraction vector into the + * vector, NeutralMolecMoleFractions_[] + */ + neutralMoleculePhase_->getMoleFractions(DATA_PTR(NeutralMolecMoleFractions_)); + + // Zero the mole fractions + fbo_zero_dbl_1(mf, m_kk); + + /* + * Use the formula matrix to calculate the relative mole numbers. + */ + for (int jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { + for (k = 0; k < m_kk; k++) { + fmij = fm_neutralMolec_ions_[k + jNeut * m_kk]; + mf[k] += fmij * NeutralMolecMoleFractions_[jNeut]; + } + } + + /* + * Normalize the new mole fractions + */ + doublereal sum = 0.0; + for (k = 0; k < m_kk; k++) { + sum += mf[k]; + } + for (k = 0; k < m_kk; k++) { + mf[k] /= sum; + } + + } + + // Calculate neutral molecule mole fractions + /* + * This routine calculates the neutral molecule mole + * fraction given the vector of ion mole fractions, + * i.e., the mole fractions from this ThermoPhase. + * Note, this routine basically assumes that there + * is charge neutrality. If there isn't, then it wouldn't + * make much sense. + * + * for the case of cIonSolnType_SINGLEANION, some slough + * in the charge neutrality is allowed. The cation number + * is followed, while the difference in charge neutrality + * is dumped into the anion mole number to fix the imbalance. + */ + void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const { + int k, icat, jNeut; + doublereal sumCat; + doublereal sumAnion; + doublereal fmij; + doublereal sum = 0.0; + + //! Zero the vector we are trying to find. + for (k = 0; k < numNeutralMoleculeSpecies_; k++) { + NeutralMolecMoleFractions_[k] = 0.0; + } +#ifdef DEBUG_MODE + sum = -1.0; + for (k = 0; k < m_kk; k++) { + sum += moleFractions_[k]; + } + if (fabs(sum) > 1.0E-11) { + throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", + "molefracts don't sum to one: " + fp2str(sum)); + } +#endif + + // bool fmSimple = true; + + switch (ionSolnType_) { + + case cIonSolnType_PASSTHROUGH: + + for (k = 0; k < m_kk; k++) { + NeutralMolecMoleFractions_[k] = moleFractions_[k]; + } + break; + + case cIonSolnType_SINGLEANION: + + sumCat = 0.0; + sumAnion = 0.0; + for (k = 0; k < numNeutralMoleculeSpecies_; k++) { + NeutralMolecMoleFractions_[k] = 0.0; + } + + for (k = 0; k < (int) cationList_.size(); k++) { + //! Get the id for the next cation + icat = cationList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + if (jNeut >= 0) { + fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk]; + AssertTrace(fmij != 0.0); + NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij; + } + } + + for (k = 0; k < numPassThroughSpecies_; k++) { + icat = passThroughList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk]; + NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij; + } + +#ifdef DEBUG_MODE + for (k = 0; k < m_kk; k++) { + moleFractionsTmp_[k] = moleFractions_[k]; + } + for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { + for (k = 0; k < m_kk; k++) { + fmij = fm_neutralMolec_ions_[k + jNeut * m_kk]; + moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut]; + } + } + for (k = 0; k < m_kk; k++) { + if (fabs(moleFractionsTmp_[k]) > 1.0E-13) { + //! Check to see if we have in fact found the inverse. + if (anionList_[0] != k) { + throw CanteraError("", "neutral molecule calc error"); + } else { + //! For the single anion case, we will allow some slippage + if (fabs(moleFractionsTmp_[k]) > 1.0E-5) { + throw CanteraError("", "neutral molecule calc error - anion"); + } + } + } + } +#endif + + // Normalize the Neutral Molecule mole fractions + sum = 0.0; + for (k = 0; k < numNeutralMoleculeSpecies_; k++) { + sum += NeutralMolecMoleFractions_[k]; + } + for (k = 0; k < numNeutralMoleculeSpecies_; k++) { + NeutralMolecMoleFractions_[k] /= sum; + } + + break; + + case cIonSolnType_SINGLECATION: + + throw CanteraError("eosType", "Unknown type"); + + break; + + case cIonSolnType_MULTICATIONANION: + + throw CanteraError("eosType", "Unknown type"); + break; + + default: + + throw CanteraError("eosType", "Unknown type"); + break; + + } + } + + void IonsFromNeutralVPSSTP::setMassFractions(const doublereal* const y) { + GibbsExcessVPSSTP::setMassFractions(y); + calcNeutralMoleculeMoleFractions(); + neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_)); + } + + void IonsFromNeutralVPSSTP::setMassFractions_NoNorm(const doublereal* const y) { + GibbsExcessVPSSTP::setMassFractions_NoNorm(y); + calcNeutralMoleculeMoleFractions(); + neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_)); + } + + void IonsFromNeutralVPSSTP::setMoleFractions(const doublereal* const y) { + GibbsExcessVPSSTP::setMoleFractions(y); + calcNeutralMoleculeMoleFractions(); + neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_)); + } + + void IonsFromNeutralVPSSTP::setMoleFractions_NoNorm(const doublereal* const y) { + GibbsExcessVPSSTP::setMoleFractions_NoNorm(y); + calcNeutralMoleculeMoleFractions(); + neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_)); + } + + + void IonsFromNeutralVPSSTP::setConcentrations(const doublereal* const y) { + GibbsExcessVPSSTP::setConcentrations(y); + calcNeutralMoleculeMoleFractions(); + neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_)); + } + + /* + * ------------ Partial Molar Properties of the Solution ------------ + */ + + + doublereal IonsFromNeutralVPSSTP::err(std::string msg) const { + throw CanteraError("IonsFromNeutralVPSSTP","Base class method " + +msg+" called. Equation of state type: "+int2str(eosType())); + 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 IonsFromNeutralVPSSTP::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 IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id) { + string stemp; + if (id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML", + "no thermo XML node"); + } + XML_Node& thermoNode = phaseNode.child("thermo"); + + + + /* + * Make sure that the thermo model is IonsFromNeutralMolecule + */ + stemp = thermoNode.attrib("model"); + string formString = lowercase(stemp); + if (formString != "ionsfromneutralmolecule") { + throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML", + "model name isn't IonsFromNeutralMolecule: " + formString); + } + + /* + * Find the Neutral Molecule Phase + */ + if (!thermoNode.hasChild("neutralMoleculePhase")) { + throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML", + "no neutralMoleculePhase XML node"); + } + XML_Node& neutralMoleculeNode = thermoNode.child("neutralMoleculePhase"); + + string nsource = neutralMoleculeNode["datasrc"]; + XML_Node *neut_ptr = get_XML_Node(nsource, 0); + if (!neut_ptr) { + throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML", + "neut_ptr = 0"); + } + + /* + * Create the neutralMolecule ThermoPhase if we haven't already + */ + if (!neutralMoleculePhase_) { + neutralMoleculePhase_ = newPhase(*neut_ptr); + } + + /* + * 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("IonsFromNeutralVPSSTP::constructPhaseXML", + "importPhase failed "); + } + + } + + + /* + * @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 IonsFromNeutralVPSSTP::initThermo() { + initLengths(); + GibbsExcessVPSSTP::initThermo(); + } + + + // Initialize lengths of local variables after all species have + // been identified. + void IonsFromNeutralVPSSTP::initLengths() { + m_kk = nSpecies(); + numNeutralMoleculeSpecies_ = neutralMoleculePhase_->nSpecies(); + moleFractions_.resize(m_kk); + fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk); + fm_invert_ionForNeutral.resize(numNeutralMoleculeSpecies_); + NeutralMolecMoleFractions_.resize(numNeutralMoleculeSpecies_); + cationList_.resize(m_kk); + anionList_.resize(m_kk); + passThroughList_.resize(m_kk); + moleFractionsTmp_.resize(m_kk); + muNeutralMolecule_.resize(numNeutralMoleculeSpecies_); + gammaNeutralMolecule_.resize(numNeutralMoleculeSpecies_); + + } + + static double factorOverlap(const std::vector& elnamesVN , + const std::vector& elemVectorN, + const int nElementsN, + const std::vector& elnamesVI , + const std::vector& elemVectorI, + const int nElementsI) + { + double fMax = 1.0E100; + for (int mi = 0; mi < nElementsI; mi++) { + if (elnamesVI[mi] != "E") { + if (elemVectorI[mi] > 1.0E-13) { + double eiNum = elemVectorI[mi]; + for (int mn = 0; mn < nElementsN; mn++) { + if (elnamesVI[mi] == elnamesVN[mn]) { + if (elemVectorN[mn] <= 1.0E-13) { + return 0.0; + } + fMax = MIN(fMax, elemVectorN[mn]/eiNum); + } + } + } + } + } + return fMax; + } + + /* + * 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 IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { + int k; + /* + * variables that need to be populated + * + * cationList_ + * numCationSpecies_; + */ + + numCationSpecies_ = 0; + cationList_.clear(); + for (k = 0; k < m_kk; k++) { + if (charge(k) > 0) { + cationList_.push_back(k); + numCationSpecies_++; + } + } + + numAnionSpecies_ = 0; + anionList_.clear(); + for (k = 0; k < m_kk; k++) { + if (charge(k) < 0) { + anionList_.push_back(k); + numAnionSpecies_++; + } + } + + numPassThroughSpecies_= 0; + passThroughList_.clear(); + for (k = 0; k < m_kk; k++) { + if (charge(k) == 0) { + passThroughList_.push_back(k); + numPassThroughSpecies_++; + } + } + + PDSS_IonsFromNeutral *speciesSS = 0; + indexSpecialSpecies_ = -1; + for (k = 0; k < m_kk; k++) { + speciesSS = dynamic_cast(providePDSS(k)); + if (!speciesSS) { + throw CanteraError("initThermoXML", "Dynamic cast failed"); + } + if (speciesSS->specialSpecies_ == 1) { + indexSpecialSpecies_ = k; + } + if (speciesSS->specialSpecies_ == 2) { + indexSecondSpecialSpecies_ = k; + } + } + + + int nElementsN = neutralMoleculePhase_->nElements(); + const std::vector& elnamesVN = neutralMoleculePhase_->elementNames(); + std::vector elemVectorN(nElementsN); + std::vector elemVectorN_orig(nElementsN); + + int nElementsI = nElements(); + const std::vector& elnamesVI = elementNames(); + vector elemVectorI(nElementsI); + + vector fm_tmp(m_kk); + for (int jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { + fm_invert_ionForNeutral[jNeut] = -1; + } + for (int jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { + for (int m = 0; m < nElementsN; m++) { + elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m); + } + elemVectorN_orig = elemVectorN; + fvo_zero_dbl_1(fm_tmp, m_kk); + + for (int m = 0; m < nElementsI; m++) { + elemVectorI[m] = nAtoms(indexSpecialSpecies_, m); + } + double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN, + elnamesVI ,elemVectorI, nElementsI); + if (fac > 0.0) { + for (int m = 0; m < nElementsN; m++) { + std::string mName = elnamesVN[m]; + for (int mi = 0; mi < nElementsI; mi++) { + std::string eName = elnamesVI[mi]; + if (mName == eName) { + elemVectorN[m] -= fac * elemVectorI[mi]; + } + + } + } + } + fm_neutralMolec_ions_[indexSpecialSpecies_ + jNeut * m_kk ] += fac; + + + for (k = 0; k < m_kk; k++) { + for (int m = 0; m < nElementsI; m++) { + elemVectorI[m] = nAtoms(k, m); + } + double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN, + elnamesVI ,elemVectorI, nElementsI); + if (fac > 0.0) { + for (int m = 0; m < nElementsN; m++) { + std::string mName = elnamesVN[m]; + for (int mi = 0; mi < nElementsI; mi++) { + std::string eName = elnamesVI[mi]; + if (mName == eName) { + elemVectorN[m] -= fac * elemVectorI[mi]; + } + + } + } + bool notTaken = true; + for (int iNeut = 0; iNeut < jNeut; iNeut++) { + if (fm_invert_ionForNeutral[iNeut] == k) { + notTaken = false; + } + } + if (notTaken) { + fm_invert_ionForNeutral[jNeut] = k; + } + } + fm_neutralMolec_ions_[k + jNeut * m_kk] += fac; + } + + // Ok check the work + for (int m = 0; m < nElementsN; m++) { + if (fabs(elemVectorN[m]) > 1.0E-13) { + throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML", + "Simple formula matrix generation failed"); + } + } + + + } + /* + * This includes the setStateFromXML calls + */ + GibbsExcessVPSSTP::initThermoXML(phaseNode, id); + + /* + * There is one extra step here. We assure ourselves that we + * have charge conservation. + */ + } + + // Update the activity coefficients + /* + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + * + * he = X_A X_B(B + C(X_A - X_B)) + */ + void IonsFromNeutralVPSSTP::s_update_lnActCoeff() const { + int k, icat, jNeut; + doublereal fmij; + /* + * Get the activity coefficiens of the neutral molecules + */ + neutralMoleculePhase_->getActivityCoefficients(DATA_PTR(gammaNeutralMolecule_)); + + switch (ionSolnType_) { + case cIonSolnType_PASSTHROUGH: + break; + case cIonSolnType_SINGLEANION: + + // Do the cation list + for (k = 0; k < (int) cationList_.size(); k++) { + //! Get the id for the next cation + icat = cationList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk]; + lnActCoeff_Scaled_[icat] = fmij * log(gammaNeutralMolecule_[jNeut]); + } + + // Do the anion list + icat = anionList_[0]; + jNeut = fm_invert_ionForNeutral[icat]; + lnActCoeff_Scaled_[icat]= 0.0; + + // Do the list of neutral molecules + for (k = 0; k < numPassThroughSpecies_; k++) { + icat = passThroughList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + lnActCoeff_Scaled_[icat] = log(gammaNeutralMolecule_[jNeut]); + } + break; + + case cIonSolnType_SINGLECATION: + throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type"); + break; + case cIonSolnType_MULTICATIONANION: + throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type"); + break; + default: + throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type"); + break; + } + + } + + /** + * Format a summary of the mixture state for output. + */ + std::string IonsFromNeutralVPSSTP::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/IonsFromNeutralVPSSTP.h b/Cantera/src/thermo/IonsFromNeutralVPSSTP.h new file mode 100644 index 000000000..33c92858e --- /dev/null +++ b/Cantera/src/thermo/IonsFromNeutralVPSSTP.h @@ -0,0 +1,695 @@ +/** + * @file IonsFromNeutralVPSSTP.h + * Header for intermediate ThermoPhase object for phases which + * consist of ions whose thermodynamics is calculated from neutral molecule thermodynamics. + * (see \ref thermoprops + * and class \link Cantera::IonsFromNeutralVPSSTP IonsFromNeutralVPSSTP\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: PseudoBinaryVPSSTP.h,v 1.1 2009/03/03 21:08:31 hkmoffa Exp $ + */ + +#ifndef CT_IONSFROMNEUTRALVPSSTP_H +#define CT_IONSFROMNEUTRALVPSSTP_H + +#include "GibbsExcessVPSSTP.h" + +namespace Cantera { + + //! enums for molten salt ion solution types + /*! + * Types identify how complicated the solution is. If there + * is just mixing on one of the sublattices but not the other, + * then the math is considerably simpler. + */ + enum IonSolnType_enumType { + cIonSolnType_PASSTHROUGH = 2000 , + cIonSolnType_SINGLEANION , + cIonSolnType_SINGLECATION , + cIonSolnType_MULTICATIONANION + }; + + /*! + * 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 IonsFromNeutralVPSSTP : public GibbsExcessVPSSTP { + + public: + + /// Constructors + /*! + * + */ + IonsFromNeutralVPSSTP(); + + //! Construct and initialize an HMWSoln ThermoPhase object + //! directly from an asci input file + /*! + * Working constructors + * + * The two constructors below are the normal way + * the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + * + * @param inputFile Name of the input file containing the phase XML data + * to set up the object + * @param id ID of the phase in the input file. Defaults to the + * empty string. + */ + IonsFromNeutralVPSSTP(std::string inputFile, std::string id = "", + ThermoPhase *neutralPhase = 0); + + + //! Construct and initialize an HMWSoln ThermoPhase object + //! directly from an XML database + /*! + * @param phaseRef XML phase node containing the description of the phase + * @param id id attribute containing the name of the phase. + * (default is the empty string) + */ + IonsFromNeutralVPSSTP(XML_Node& phaseRoot, std::string id = "", + ThermoPhase *neutralPhase = 0); + + + //! Copy constructor + /*! + * Note this stuff will not work until the underlying phase + * has a working copy constructor + * + * @param b class to be copied + */ + IonsFromNeutralVPSSTP(const IonsFromNeutralVPSSTP &b); + + /// Assignment operator + /*! + * + * @param b class to be copied. + */ + IonsFromNeutralVPSSTP & operator=(const IonsFromNeutralVPSSTP &b); + + /// Destructor. + virtual ~IonsFromNeutralVPSSTP(); + + //! 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 + * @{ + */ + + //! Return the Molar enthalpy. Units: J/kmol. + /*! + * This is calculated from the partial molar enthalpies of the species + */ + virtual doublereal enthalpy_mole() const; + + /** + * Molar internal energy. J/kmol. + * * + * This is calculated from the soln enthalpy and then + * subtracting pV. + */ + virtual doublereal intEnergy_mole() const; + + /** + * Molar entropy. Units: J/kmol/K. + * + * + */ + virtual doublereal entropy_mole() const; + + /** + * Molar Gibbs free Energy for an ideal gas. + * Units = J/kmol. + */ + virtual doublereal gibbs_mole() const; + + /** + * Molar heat capacity at constant pressure. Units: J/kmol/K. + * For an ideal gas mixture, + * + */ + virtual doublereal cp_mole() const; + + /** + * Molar heat capacity at constant volume. Units: J/kmol/K. + * + */ + virtual doublereal cv_mole() const; + + + /** + * @} + * @name Utilities + * @{ + */ + + + + + /** + * @} + * @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. + * @{ + */ + + //! This method returns an array of generalized concentrations + /*! + * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / + * C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration + * defined below and \f$ a_k \f$ are activities used in the + * thermodynamic functions. These activity (or generalized) + * concentrations are used + * by kinetics manager classes to compute the forward and + * reverse rates of elementary reactions. Note that they may + * or may not have units of concentration --- they might be + * partial pressures, mole fractions, or surface coverages, + * for example. + * + * @param c Output array of generalized concentrations. The + * units depend upon the implementation of the + * reaction rate expressions within the phase. + */ + virtual void getActivityConcentrations(doublereal* c) const; + + //! 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; + + + /** + * 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 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; + + + + //@} + /// @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. + * @{ + */ + + virtual void setTemperature(doublereal t); + virtual void setPressure(doublereal p); + + //! Set the temperature (K) and pressure (Pa) + /*! + * Setting the pressure may involve the solution of a nonlinear equation. + * + * @param t Temperature (K) + * @param p Pressure (Pa) + */ + virtual void setState_TP(doublereal t, doublereal p); + + + //! Calculate ion mole fractions from neutral molecule + //! mole fractions. + /*! + * @param mf Dump the mole fractions into this vector. + */ + virtual void calcIonMoleFractions(doublereal *const mf) const; + + //! Calculate neutral molecule mole fractions + /*! + * This routine calculates the neutral molecule mole + * fraction given the vector of ion mole fractions, + * i.e., the mole fractions from this ThermoPhase. + * Note, this routine basically assumes that there + * is charge neutrality. If there isn't, then it wouldn't + * make much sense. + * + * for the case of cIonSolnType_SINGLEANION, some slough + * in the charge neutrality is allowed. The cation number + * is followed, while the difference in charge neutrality + * is dumped into the anion mole number to fix the imbalance. + */ + virtual void calcNeutralMoleculeMoleFractions() const; + + /** + * Set the mass fractions to the specified values, and then + * normalize them so that they sum to 1.0. + * @param y Array of unnormalized mass fraction values (input). + * Must have a length greater than or equal to the number of + * species. + * + * @param y Input vector of mass fractions. + * Length is m_kk. + */ + virtual void setMassFractions(const doublereal* const y); + + /** + * Set the mass fractions to the specified values without + * normalizing. This is useful when the normalization + * condition is being handled by some other means, for example + * by a constraint equation as part of a larger set of + * equations. + * + * @param y Input vector of mass fractions. + * Length is m_kk. + */ + virtual void setMassFractions_NoNorm(const doublereal* const y); + + /** + * Set the mole fractions to the specified values, and then + * normalize them so that they sum to 1.0. + * @param x Array of unnormalized mole fraction values (input). + * Must have a length greater than or equal to the number of + * species. + * + * @param x Input vector of mole fractions. + * Length is m_kk. + */ + virtual void setMoleFractions(const doublereal* const x); + + /** + * Set the mole fractions to the specified values without + * normalizing. This is useful when the normalization + * condition is being handled by some other means, for example + * by a constraint equation as part of a larger set of + * equations. + * + * @param x Input vector of mole fractions. + * Length is m_kk. + */ + virtual void setMoleFractions_NoNorm(const doublereal* const x); + + /** + * Set the concentrations to the specified values within the + * phase. + * + * @param c The input vector to this routine is in dimensional + * units. For volumetric phases c[k] is the + * concentration of the kth species in kmol/m3. + * For surface phases, c[k] is the concentration + * in kmol/m2. The length of the vector is the number + * of species in the phase. + */ + virtual void setConcentrations(const doublereal* const c); + + //@} + + /** + * @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. + + //! Initialization of a HMWSoln phase using an xml file + /*! + * This routine is a precursor to initThermo(XML_Node*) + * 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 HMWSoln 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); + + + /*! + * @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: + + //! Ion solution type + /*! + * There is either mixing on the anion, cation, or both lattices. + * There is also a passthrough option + * + * Defaults to cIonSolnType_SINGLEANION, so that LiKCl can be hardwired + */ + IonSolnType_enumType ionSolnType_; + + //! Number of neutral molecule species + /*! + * This is equal to the number of species in the + * neutralMoleculePhase_ ThermoPhase. + */ + int numNeutralMoleculeSpecies_; + + //! Index of special species + int indexSpecialSpecies_; + + //! Index of special species + int indexSecondSpecialSpecies_; + + //! Formula Matrix for composition of neutral molecules + //! in terms of the molecules in this ThermoPhase + /*! + * fm_neutralMolec_ions[ i + jNeut * NumNeut ] + * + * This is the number of ions of type i in the neutral + * molecule jNeut. + */ + std::vector fm_neutralMolec_ions_; + + //! Mapping between ion species and neutral molecule for quick invert. + /*! + * For the case of fm_invert_simple_ = true, we assume that there + * is a quick way to invert the formula matrix so that we can + * quickly calculate the neutral molecule mole fraction + * given the ion mole fraction vector. + * + * We assume that for a selected set of ion species, that that + * ion is only in the neutral molecule, jNeut. + * + * therefore, + * + * NeutralMolecMoleFractions_[jNeut] += moleFractions_[i_ion] / fmij; + * + * where fmij is the number of ions in neutral molecule jNeut. + * + * Thus, we formulate the neutral molecule mole fraction NeutralMolecMoleFractions_[] + * vector from this association. We further assume that there are + * no other associations. If fm_invert_simple_ is not true, + * then we need to do a formal inversion which takes a great + * deal of time and is not currently implemented. + */ + std::vector fm_invert_ionForNeutral; + + //! Mole fractions using the Neutral Molecule Mole fraction basis + mutable std::vector NeutralMolecMoleFractions_; + + + //! List of the species in this ThermoPhase which are cation species + std::vector cationList_; + + //! Number of cation species + int numCationSpecies_; + + //! List of the species in this ThermoPhase which are anion species + std::vectoranionList_; + + //! Number of anion species + int numAnionSpecies_; + + //! List of the species in this ThermoPhase which are passed + //! through to the neutralMoleculePhase ThermoPhase. + /*! + * These have neutral charges. + */ + std::vector passThroughList_; + + //! Number of the species in this ThermoPhase which are passed + //! through to the neutralMoleculePhase ThermoPhase + int numPassThroughSpecies_; + + public: + ThermoPhase *neutralMoleculePhase_; + protected: + //! ThermoPhase for the cation lattice + /*! + * Currently this is unimplemented and may be deleted + */ + ThermoPhase *cationPhase_; + + //! ThermoPhase for the anion lattice + /*! + * Currently this is unimplemented and may be deleted + */ + ThermoPhase *anionPhase_; + + //! Temporary mole fraction vector + mutable std::vector moleFractionsTmp_; + + mutable std::vector muNeutralMolecule_; + mutable std::vector gammaNeutralMolecule_; + + private: + + + }; + + + + + +} + +#endif diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index eb261de8f..3a31d01d4 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -66,22 +66,24 @@ ELECTRO_OBJ = MolalityVPSSTP.o VPStandardStateTP.o \ WaterPropsIAPWSphi.o WaterPropsIAPWS.o WaterProps.o \ PDSS.o PDSS_Water.o PDSS_HKFT.o \ HMWSoln.o HMWSoln_input.o DebyeHuckel.o \ - WaterSSTP.o \ + WaterSSTP.o \ VPSSMgr_Water_ConstVol.o VPSSMgr_Water_HKFT.o ELECTRO_H = MolalityVPSSTP.h VPStandardStateTP.h \ IdealMolalSoln.h \ WaterPropsIAPWSphi.h WaterPropsIAPWS.h WaterProps.h \ PDSS.h PDSS_Water.h PDSS_HKFT.h \ - HMWSoln.h electrolytes.h \ + HMWSoln.h electrolytes.h \ DebyeHuckel.h WaterSSTP.h VPSSMgr_Water_HKFT.h \ VPSSMgr_Water_ConstVol.h endif ifeq ($(do_issp),1) ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o \ - GibbsExcessVPSSTP.o PseudoBinaryVPSSTP.o MargulesVPSSTP.o + GibbsExcessVPSSTP.o PseudoBinaryVPSSTP.o MargulesVPSSTP.o \ + IonsFromNeutralVPSSTP.o PDSS_IonsFromNeutral.o ISSP_H = IdealSolidSolnPhase.h StoichSubstanceSSTP.h SingleSpeciesTP.h MineralEQ3.h \ - GibbsExcessVPSSTP.h PseudoBinaryVPSSTP.h MargulesVPSSTP.h + GibbsExcessVPSSTP.h PseudoBinaryVPSSTP.h MargulesVPSSTP.h \ + IonsFromNeutralVPSSTP.h PDSS_IonsFromNeutral.h endif CATHERMO_OBJ = $(THERMO_OBJ) $(ELECTRO_OBJ) $(ISSP_OBJ) diff --git a/Cantera/src/thermo/MargulesVPSSTP.cpp b/Cantera/src/thermo/MargulesVPSSTP.cpp index b1b03e855..37ca850da 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.cpp +++ b/Cantera/src/thermo/MargulesVPSSTP.cpp @@ -31,13 +31,41 @@ namespace Cantera { * */ MargulesVPSSTP::MargulesVPSSTP() : - PseudoBinaryVPSSTP(), + GibbsExcessVPSSTP(), numBinaryInteractions_(0), formMargules_(0), formTempModel_(0) { } + /* + * Working constructors + * + * The two constructors below are the normal way + * the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + + */ + MargulesVPSSTP::MargulesVPSSTP(std::string inputFile, std::string id) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + constructPhaseFile(inputFile, id); + } + + MargulesVPSSTP::MargulesVPSSTP(XML_Node& phaseRoot, std::string id) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + constructPhaseXML(phaseRoot, id); + } + + /* * Copy Constructor: * @@ -45,7 +73,7 @@ namespace Cantera { * has a working copy constructor */ MargulesVPSSTP::MargulesVPSSTP(const MargulesVPSSTP &b) : - PseudoBinaryVPSSTP() + GibbsExcessVPSSTP() { *this = operator=(b); } @@ -59,7 +87,7 @@ namespace Cantera { MargulesVPSSTP& MargulesVPSSTP:: operator=(const MargulesVPSSTP &b) { if (&b != this) { - PseudoBinaryVPSSTP::operator=(b); + GibbsExcessVPSSTP::operator=(b); } numBinaryInteractions_ = b.numBinaryInteractions_ ; @@ -105,7 +133,7 @@ namespace Cantera { * */ MargulesVPSSTP::MargulesVPSSTP(int testProb) : - PseudoBinaryVPSSTP(), + GibbsExcessVPSSTP(), numBinaryInteractions_(0), formMargules_(0), formTempModel_(0) @@ -262,63 +290,14 @@ namespace Cantera { 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); - } - } - - + * Make sure that the thermo model is Margules + */ + stemp = thermoNode.attrib("model"); + string formString = lowercase(stemp); + if (formString != "margules") { + throw CanteraError("MargulesVPSSTP::constructPhaseXML", + "model name isn't Margules: " + formString); + } /* @@ -406,7 +385,7 @@ namespace Cantera { } } - + @@ -437,7 +416,7 @@ namespace Cantera { */ void MargulesVPSSTP::initThermo() { initLengths(); - PseudoBinaryVPSSTP::initThermo(); + GibbsExcessVPSSTP::initThermo(); } @@ -445,7 +424,7 @@ namespace Cantera { // been identified. void MargulesVPSSTP::initLengths() { m_kk = nSpecies(); - moleFractions_.resize(m_kk); + } /* @@ -464,8 +443,64 @@ namespace Cantera { * with the correct id. */ void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { + string subname = "MargulesVPSSTP::initThermoXML"; + string stemp; + + /* + * Check on the thermo field. Must have: + * + */ + + XML_Node& thermoNode = phaseNode.child("thermo"); + string mStringa = thermoNode.attrib("model"); + string mString = lowercase(mStringa); + if (mString != "margules") { + throw CanteraError(subname.c_str(), + "Unknown thermo model: " + mStringa); + } + + + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + XML_Node *acNodePtr = 0; + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& acNode = thermoNode.child("activityCoefficients"); + acNodePtr = &acNode; + string mStringa = thermoNode.attrib("model"); + string mString = lowercase(mStringa); + if (mString != "margules") { + throw CanteraError(subname.c_str(), + "Unknown activity coefficient model: " + mStringa); + } + int n = acNodePtr->nChildren(); + for (int i = 0; i < n; i++) { + XML_Node &xmlACChild = acNodePtr->child(i); + stemp = xmlACChild.name(); + string nodeName = lowercase(stemp); + /* + * Process a binary salt field, or any of the other XML fields + * that make up the Pitzer Database. Entries will be ignored + * if any of the species in the entry isn't in the solution. + */ + if (nodeName == "binaryneutralspeciesparameters") { + readXMLBinarySpecies(xmlACChild); + + } + } + } + + /* + * Go down the chain + */ + GibbsExcessVPSSTP::initThermoXML(phaseNode, id); + - PseudoBinaryVPSSTP::initThermoXML(phaseNode, id); } // Update the activity coefficients @@ -473,16 +508,15 @@ namespace Cantera { * 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) + * he = X_A X_B(B + C(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; - } + fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + double RT = GasConstant * temperature(); for (int i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; @@ -496,13 +530,122 @@ namespace Cantera { lnActCoeff_Scaled_[iA] += XB * XB * (g0 + g1 * (XB - XA)); lnActCoeff_Scaled_[iB] += XA * XA * g0 + XA * XB * g1 * (2 * XA); - } + } + void MargulesVPSSTP::resizeNumInteractions(const int num) { + numBinaryInteractions_ = num; + m_HE_b_ij.resize(num, 0.0); + m_HE_c_ij.resize(num, 0.0); + m_HE_d_ij.resize(num, 0.0); + m_SE_b_ij.resize(num, 0.0); + m_SE_c_ij.resize(num, 0.0); + m_SE_d_ij.resize(num, 0.0); + + m_pSpecies_A_ij.resize(num, -1); + m_pSpecies_B_ij.resize(num, -1); + } - /** + + /* + * Process an XML node called "binaryNeutralSpeciesParameters" + * This node contains all of the parameters necessary to describe + * the Margules Interaction for a single binary interaction + * This function reads the XML file and writes the coefficients + * it finds to an internal data structures. + */ + void MargulesVPSSTP::readXMLBinarySpecies(XML_Node &xmLBinarySpecies) { + string xname = xmLBinarySpecies.name(); + if (xname != "binaryNeutralSpeciesParameters") { + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", + "Incorrect name for processing this routine: " + xname); + } + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + int nParamsFound; + vector_fp vParams; + string iName = xmLBinarySpecies.attrib("speciesA"); + if (iName == "") { + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "no speciesA attrib"); + } + string jName = xmLBinarySpecies.attrib("speciesB"); + if (jName == "") { + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "no speciesB attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species + */ + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + string ispName = speciesName(iSpecies); + if (charge[iSpecies] != 0) { + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "speciesA charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + string jspName = speciesName(jSpecies); + if (charge[jSpecies] != 0) { + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "speciesB charge problem"); + } + + resizeNumInteractions(numBinaryInteractions_ + 1); + int iSpot = numBinaryInteractions_ - 1; + m_pSpecies_A_ij[iSpot] = iSpecies; + m_pSpecies_B_ij[iSpot] = jSpecies; + + int num = xmLBinarySpecies.nChildren(); + for (int iChild = 0; iChild < num; iChild++) { + XML_Node &xmlChild = xmLBinarySpecies.child(iChild); + stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + /* + * Process the binary species interaction child elements + */ + if (nodeName == "excessenthalpy") { + /* + * Get the string containing all of the values + */ + getFloatArray(xmlChild, vParams, true, "", "excessEnthalpy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEnthalpy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_HE_b_ij[iSpot] = vParams[0]; + m_HE_c_ij[iSpot] = vParams[1]; + } + + if (nodeName == "excessentropy") { + /* + * Get the string containing all of the values + */ + getFloatArray(xmlChild, vParams, true, "", "excessEntropy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEntropy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_SE_b_ij[iSpot] = vParams[0]; + m_SE_c_ij[iSpot] = vParams[1]; + } + + + } + + } + + /** * Format a summary of the mixture state for output. */ std::string MargulesVPSSTP::report(bool show_thermo) const { diff --git a/Cantera/src/thermo/MargulesVPSSTP.h b/Cantera/src/thermo/MargulesVPSSTP.h index 471aa03c4..cf5f25cea 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.h +++ b/Cantera/src/thermo/MargulesVPSSTP.h @@ -3,7 +3,7 @@ * 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). + * and class \link Cantera::MargulesVPSSTP MargulesVPSSTP\endlink). * * Header file for a derived class of ThermoPhase that handles * variable pressure standard state methods for calculating @@ -31,13 +31,268 @@ namespace Cantera { * @ingroup thermoprops */ + + //! MargulesVPSSTP is a derived class of GibbsExcessVPSSTP that employs + //! the Margules approximation for the excess gibbs free energy /*! - * MargulesVPSSTP is a derived class of PseudoBinaryVPSSTP. + * + * %MargulesVPSSTP derives from class GibbsExcessVPSSTP which is derived + * from VPStandardStateTP, + * and overloads the virtual methods defined there with ones that + * use expressions appropriate for the Margules Excess gibbs free energy + * approximation. + * + * The independent unknowns are pressure, temperature, and mass fraction. + * + * Several concepts are introduced. The first concept is there are temporary + * variables for holding the species standard state values + * of Cp, H, S, G, and V at the + * last temperature and pressure called. These functions are not recalculated + * if a new call is made using the previous temperature and pressure. Currently, + * these variables and the calculation method are handled by the VPSSMgr class, + * for which VPStandardStateTP owns a pointer to. + * + * To support the above functionality, pressure and temperature variables, + * m_plast_ss and m_tlast_ss, are kept which store the last pressure and temperature + * used in the evaluation of standard state properties. + * + * This class is usually used for nearly incompressible phases. For those phases, it + * makes sense to change the equation of state independent variable from + * density to pressure. The variable m_Pcurrent contains the current value of the + * pressure within the phase. * * + *
+ *

Specification of Species Standard %State Properties

+ *
* - */ - class MargulesVPSSTP : public PseudoBinaryVPSSTP { + * All species are defined to have standard states that depend upon both + * the temperature and the pressure. The Margules approximation assumes + * symmetric standard states, where all of the standard state assume + * that the species are in pure component states at the temperatue + * and pressure of the solution. I don't think it prevents, however, + * some species from being dilute in the solution. + * + * + *
+ *

Specification of Solution Thermodynamic Properties

+ *
+ * + * The excess Gibbs free energy + * + * \f[ + * G^E = \sum_i \left( H_{Ei} - T S_{Ei} \right) + * \f] + * \f[ + * H^E_i = X_{Ai} X_{Bi} \left( h_{o,i} + h_{1,i} X_{Bi} \right) + * \f] + * \f[ + * S^E_i = X_{Ai} X_{Bi} \left( s_{o,i} + s_{1,i} X_{Bi} \right) + * \f] + * + * + * The activity of a species defined in the phase is given by an excess + * Gibbs free energy formulation. + * + * \f[ + * a_k = \gamma_k X_k + * \f] + * + * where \f$ X_k \f$ is the mole fraction of species k. + * The chemical potential for species k is equal to + * + * \f[ + * \mu_k(T,P) = \mu^o_k(T, P) + R T \log(\gamma_k X_k) + * \f] + * + * In terms of the reference state, the above can be rewritten + * + * + * \f[ + * \mu_k(T,P) = \mu^{ref}_k(T, P) + R T \log(\frac{P X_k}{P_{ref}}) + * \f] + * + * The partial molar entropy for species k is given by the following relation, + * + * \f[ + * \tilde{s}_k(T,P) = s^o_k(T,P) - R \log(X_k) = s^{ref}_k(T) - R \log(\frac{P X_k}{P_{ref}}) + * \f] + * + * The partial molar enthalpy for species k is + * + * \f[ + * \tilde{h}_k(T,P) = h^o_k(T,P) = h^{ref}_k(T) + * \f] + * + * The partial molar Internal Energy for species k is + * + * \f[ + * \tilde{u}_k(T,P) = u^o_k(T,P) = u^{ref}_k(T) + * \f] + * + * The partial molar Heat Capacity for species k is + * + * \f[ + * \tilde{Cp}_k(T,P) = Cp^o_k(T,P) = Cp^{ref}_k(T) + * \f] + * + *
+ *

%Application within %Kinetics Managers

+ *
+ * + * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / + * C^s_k, \f$ where \f$ C^s_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. + * The activity concentration,\f$ C^a_k \f$,is given by the following expression. + * + * \f[ + * C^a_k = C^s_k X_k = \frac{P}{R T} X_k + * \f] + * + * The standard concentration for species k is independent of k and equal to + * + * \f[ + * C^s_k = C^s = \frac{P}{R T} + * \f] + * + * For example, a bulk-phase binary gas reaction between species j and k, producing + * a new gas species l would have the + * following equation for its rate of progress variable, \f$ R^1 \f$, which has + * units of kmol m-3 s-1. + * + * \f[ + * R^1 = k^1 C_j^a C_k^a = k^1 (C^s a_j) (C^s a_k) + * \f] + * where + * \f[ + * C_j^a = C^s a_j \mbox{\quad and \quad} C_k^a = C^s a_k + * \f] + * + * + * \f$ C_j^a \f$ is the activity concentration of species j, and + * \f$ C_k^a \f$ is the activity concentration of species k. \f$ C^s \f$ + * is the standard concentration. \f$ a_j \f$ is + * the activity of species j which is equal to the mole fraction of j. + * + * The reverse rate constant can then be obtained from the law of microscopic reversibility + * and the equilibrium expression for the system. + * + * \f[ + * \frac{a_j a_k}{ a_l} = K_a^{o,1} = \exp(\frac{\mu^o_l - \mu^o_j - \mu^o_k}{R T} ) + * \f] + * + * \f$ K_a^{o,1} \f$ is the dimensionless form of the equilibrium constant, associated with + * the pressure dependent standard states \f$ \mu^o_l(T,P) \f$ and their associated activities, + * \f$ a_l \f$, repeated here: + * + * \f[ + * \mu_l(T,P) = \mu^o_l(T, P) + R T \log(a_l) + * \f] + * + * We can switch over to expressing the equilibrium constant in terms of the reference + * state chemical potentials + * + * \f[ + * K_a^{o,1} = \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) * \frac{P_{ref}}{P} + * \f] + * + * The concentration equilibrium constant, \f$ K_c \f$, may be obtained by changing over + * to activity concentrations. When this is done: + * + * \f[ + * \frac{C^a_j C^a_k}{ C^a_l} = C^o K_a^{o,1} = K_c^1 = + * \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) * \frac{P_{ref}}{RT} + * \f] + * + * %Kinetics managers will calculate the concentration equilibrium constant, \f$ K_c \f$, + * using the second and third part of the above expression as a definition for the concentration + * equilibrium constant. + * + * For completeness, the pressure equilibrium constant may be obtained as well + * + * \f[ + * \frac{P_j P_k}{ P_l P_{ref}} = K_p^1 = \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) + * \f] + * + * \f$ K_p \f$ is the simplest form of the equilibrium constant for ideal gases. However, it isn't + * necessarily the simplest form of the equilibrium constant for other types of phases; \f$ K_c \f$ is + * used instead because it is completely general. + * + * The reverse rate of progress may be written down as + * \f[ + * R^{-1} = k^{-1} C_l^a = k^{-1} (C^o a_l) + * \f] + * + * where we can use the concept of microscopic reversibility to + * write the reverse rate constant in terms of the + * forward reate constant and the concentration equilibrium + * constant, \f$ K_c \f$. + * + * \f[ + * k^{-1} = k^1 K^1_c + * \f] + * + * \f$k^{-1} \f$ has units of s-1. + * + * + *
+ *

Instantiation of the Class

+ *
+ * + * + * The constructor for this phase is located in the default ThermoFactory + * for %Cantera. A new %IdealGasPhase may be created by the following code + * snippet: + * + * @code + * XML_Node *xc = get_XML_File("silane.xml"); + * XML_Node * const xs = xc->findNameID("phase", "silane"); + * ThermoPhase *silane_tp = newPhase(*xs); + * IdealGasPhase *silaneGas = dynamic_cast (silane_tp); + * @endcode + * + * or by the following constructor: + * + * @code + * XML_Node *xc = get_XML_File("silane.xml"); + * XML_Node * const xs = xc->findNameID("phase", "silane"); + * IdealGasPhase *silaneGas = new IdealGasPhase(*xs); + * @endcode + * + *
+ *

XML Example

+ *
+ * An example of an XML Element named phase setting up a IdealGasPhase + * object named silane is given below. + * + * + * @verbatim + + + Si H He + + H2 H HE SIH4 SI SIH SIH2 SIH3 H3SISIH SI2H6 + H2SISIH2 SI3H8 SI2 SI3 + + + + + + + @endverbatim + * + * The model attribute "IdealGas" of the thermo XML element identifies the phase as + * being of the type handled by the IdealGasPhase object. + * + * @ingroup thermoprops + * + + */ + class MargulesVPSSTP : public GibbsExcessVPSSTP { public: @@ -52,6 +307,33 @@ namespace Cantera { */ MargulesVPSSTP(); + //! Construct and initialize a MargulesVPSSTP ThermoPhase object + //! directly from an xml input file + /*! + * Working constructors + * + * The two constructors below are the normal way + * the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + * + * @param inputFile Name of the input file containing the phase XML data + * to set up the object + * @param id ID of the phase in the input file. Defaults to the + * empty string. + */ + MargulesVPSSTP(std::string inputFile, std::string id = ""); + + //! Construct and initialize a MargulesVPSSTP ThermoPhase object + //! directly from an XML database + /*! + * @param phaseRef XML phase node containing the description of the phase + * @param id id attribute containing the name of the phase. + * (default is the empty string) + */ + MargulesVPSSTP(XML_Node& phaseRef, std::string id = ""); + + //! Special constructor for a hard-coded problem /*! * @@ -68,7 +350,7 @@ namespace Cantera { * * @param b class to be copied */ - MargulesVPSSTP(const MargulesVPSSTP&b); + MargulesVPSSTP(const MargulesVPSSTP&b); //! Assignment operator /*! @@ -357,6 +639,19 @@ namespace Cantera { private: + //! Process an XML node called "binaryNeutralSpeciesParameters" + /*! + * This node contains all of the parameters necessary to describe + * the Margules model for a particular binary interaction. + * This function reads the XML file and writes the coefficients + * it finds to an internal data structures. + * + * @param BinSalt reference to the XML_Node named "binaryNeutralSpeciesParameters" + * containing the binary interaction + */ + void readXMLBinarySpecies(XML_Node &xmLBinarySpecies); + + void resizeNumInteractions(const int num); //! Initialize lengths of local variables after all species have diff --git a/Cantera/src/thermo/PDSS_IonsFromNeutral.cpp b/Cantera/src/thermo/PDSS_IonsFromNeutral.cpp new file mode 100644 index 000000000..1c598b04a --- /dev/null +++ b/Cantera/src/thermo/PDSS_IonsFromNeutral.cpp @@ -0,0 +1,484 @@ +/** + * @file PDSS_IonsFromNeutral.cpp + * Implementation of a pressure dependent standard state + * virtual function. + */ +/* + * 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: PDSS_IonsFromNeutral.cpp,v 1.8 2009/01/04 06:34:20 hkmoffa Exp $ + */ + +#include "ct_defs.h" +#include "xml.h" +#include "ctml.h" +#include "PDSS_IonsFromNeutral.h" +#include "ThermoFactory.h" +#include "IonsFromNeutralVPSSTP.h" + +#include "VPStandardStateTP.h" + +using namespace std; + +namespace Cantera { + /** + * Basic list of constructors and duplicators + */ + + PDSS_IonsFromNeutral::PDSS_IonsFromNeutral(VPStandardStateTP *tp, int spindex) : + PDSS(tp, spindex), + neutralMoleculePhase_(0), + numMult_(0), + add2RTln2_(true), + specialSpecies_(0) + { + m_pdssType = cPDSS_IONSFROMNEUTRAL; + } + + + PDSS_IonsFromNeutral::PDSS_IonsFromNeutral(VPStandardStateTP *tp, int spindex, + std::string inputFile, std::string id) : + PDSS(tp, spindex), + neutralMoleculePhase_(0), + numMult_(0), + add2RTln2_(true), + specialSpecies_(0) + { + m_pdssType = cPDSS_IONSFROMNEUTRAL; + constructPDSSFile(tp, spindex, inputFile, id); + } + + + + PDSS_IonsFromNeutral::PDSS_IonsFromNeutral(VPStandardStateTP *tp, int spindex, const XML_Node& speciesNode, + const XML_Node& phaseRoot, bool spInstalled) : + PDSS(tp, spindex), + neutralMoleculePhase_(0), + numMult_(0), + add2RTln2_(true), + specialSpecies_(0) + { + if (!spInstalled) { + throw CanteraError("PDSS_IonsFromNeutral", "sp installing not done yet"); + } + m_pdssType = cPDSS_IONSFROMNEUTRAL; + std::string id = ""; + constructPDSSXML(tp, spindex, speciesNode, phaseRoot, id); + } + + + + PDSS_IonsFromNeutral::PDSS_IonsFromNeutral(const PDSS_IonsFromNeutral &b) : + PDSS(b) + { + /* + * Use the assignment operator to do the brunt + * of the work for the copy construtor. + */ + *this = b; + } + + /** + * Assignment operator + */ + PDSS_IonsFromNeutral& PDSS_IonsFromNeutral::operator=(const PDSS_IonsFromNeutral&b) { + if (&b == this) return *this; + PDSS::operator=(b); + + m_tmin = b.m_tmin; + m_tmax = b.m_tmax; + neutralMoleculePhase_ = b.neutralMoleculePhase_; + numMult_ = b.numMult_; + idNeutralMoleculeVec = b.idNeutralMoleculeVec; + factorVec = b.factorVec; + add2RTln2_ = b.add2RTln2_; + specialSpecies_ = b.specialSpecies_; + + return *this; + } + + PDSS_IonsFromNeutral::~PDSS_IonsFromNeutral() { + } + + //! Duplicator + PDSS* PDSS_IonsFromNeutral::duplMyselfAsPDSS() const { + PDSS_IonsFromNeutral * idg = new PDSS_IonsFromNeutral(*this); + return (PDSS *) idg; + } + + /** + * constructPDSSXML: + * + * Initialization of a PDSS_IonsFromNeutral object using an + * xml file. + + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void PDSS_IonsFromNeutral::constructPDSSXML(VPStandardStateTP *tp, int spindex, + const XML_Node& speciesNode, + const XML_Node& phaseNode, std::string id) { + const XML_Node *tn = speciesNode.findByName("thermo"); + if (!tn) { + throw CanteraError("PDSS_IonsFromNeutral::constructPDSSXML", + "no thermo Node for species " + speciesNode.name()); + } + std::string model = lowercase((*tn)["model"]); + if (model != "ionfromneutral") { + throw CanteraError("PDSS_IonsFromNeutral::constructPDSSXML", + "thermo model for species isn't IonsFromNeutral: " + + speciesNode.name()); + } + const XML_Node *nsm = tn->findByName("neutralSpeciesMultipliers"); + if (!nsm) { + throw CanteraError("PDSS_IonsFromNeutral::constructPDSSXML", + "no Thermo::neutralSpeciesMultipliers Node for species " + speciesNode.name()); + } + + IonsFromNeutralVPSSTP *ionPhase = dynamic_cast(tp); + neutralMoleculePhase_ = ionPhase->neutralMoleculePhase_; + + std::vector key; + std::vector val; + + /* + * + */ + numMult_ = ctml::getPairs(*nsm, key, val); + idNeutralMoleculeVec.resize(numMult_); + factorVec.resize(numMult_); + tmpNM.resize(neutralMoleculePhase_->nSpecies()); + + for (int i = 0; i < numMult_; i++) { + idNeutralMoleculeVec[i] = neutralMoleculePhase_->speciesIndex(key[i]); + factorVec[i] = fpValueCheck(val[i]); + } + specialSpecies_ = 0; + const XML_Node *ss = tn->findByName("specialSpecies"); + if (ss) { + specialSpecies_ = 1; + } + const XML_Node *sss = tn->findByName("secondSpecialSpecies"); + if (sss) { + specialSpecies_ = 2; + } + add2RTln2_ = true; + if (specialSpecies_ == 1) { + add2RTln2_ = false; + } + + } + + + void PDSS_IonsFromNeutral::constructPDSSFile(VPStandardStateTP *tp, int spindex, + std::string inputFile, std::string id) { + + if (inputFile.size() == 0) { + throw CanteraError("PDSS_IonsFromNeutral::constructPDSSFile", + "input file is null"); + } + std::string path = findInputFile(inputFile); + ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("PDSS_IonsFromNeutral::constructPDSSFile","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("PDSS_IonsFromNeutral::constructPDSSFile", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + + XML_Node& speciesList = fxml_phase->child("speciesArray"); + XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"], + &(fxml_phase->root())); + const vector&sss = tp->speciesNames(); + + const XML_Node* s = speciesDB->findByAttr("name", sss[spindex]); + + constructPDSSXML(tp, spindex, *s, *fxml_phase, id); + delete fxml; + } + + void PDSS_IonsFromNeutral::initThermoXML(const XML_Node& phaseNode, std::string &id) { + PDSS::initThermoXML(phaseNode, id); + } + + void PDSS_IonsFromNeutral::initThermo() { + PDSS::initThermo(); + SpeciesThermo &sp = m_tp->speciesThermo(); + m_p0 = sp.refPressure(m_spindex); + m_minTemp = m_spthermo->minTemp(m_spindex); + m_maxTemp = m_spthermo->maxTemp(m_spindex); + } + + /** + * Return the molar enthalpy in units of J kmol-1 + */ + doublereal + PDSS_IonsFromNeutral::enthalpy_mole() const { + doublereal val = enthalpy_RT(); + doublereal RT = GasConstant * m_temp; + return (val * RT); + } + + doublereal + PDSS_IonsFromNeutral::enthalpy_RT() const { + neutralMoleculePhase_->getEnthalpy_RT(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + return val; + } + + + /** + * Calculate the internal energy in mks units of + * J kmol-1 + */ + doublereal + PDSS_IonsFromNeutral::intEnergy_mole() const { + doublereal val = m_h0_RT_ptr[m_spindex] - 1.0; + doublereal RT = GasConstant * m_temp; + return (val * RT); + } + + /** + * Calculate the entropy in mks units of + * J kmol-1 K-1 + */ + doublereal + PDSS_IonsFromNeutral::entropy_mole() const { + doublereal val = entropy_R(); + return (val * GasConstant); + } + + doublereal + PDSS_IonsFromNeutral::entropy_R() const { + neutralMoleculePhase_->getEntropy_R(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + if (add2RTln2_) { + val -= 2.0 * log(2.0); + } + return val; + } + + /** + * Calculate the Gibbs free energy in mks units of + * J kmol-1 K-1. + */ + doublereal + PDSS_IonsFromNeutral::gibbs_mole() const { + doublereal val = gibbs_RT(); + doublereal RT = GasConstant * m_temp; + return (val * RT); + } + + doublereal + PDSS_IonsFromNeutral::gibbs_RT() const { + neutralMoleculePhase_->getGibbs_RT(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + if (add2RTln2_) { + val += 2.0 * log(2.0); + } + return val; + } + + /** + * Calculate the constant pressure heat capacity + * in mks units of J kmol-1 K-1 + */ + doublereal + PDSS_IonsFromNeutral::cp_mole() const { + doublereal val = cp_R(); + return (val * GasConstant); + } + + doublereal + PDSS_IonsFromNeutral::cp_R() const { + neutralMoleculePhase_->getCp_R(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + return val; + } + + doublereal + PDSS_IonsFromNeutral::molarVolume() const { + neutralMoleculePhase_->getStandardVolumes(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + return val; + } + + + doublereal + PDSS_IonsFromNeutral::density() const { + return (m_pres * m_mw / (GasConstant * m_temp)); + } + + /* + * Calculate the constant volume heat capacity + * in mks units of J kmol-1 K-1 + */ + doublereal + PDSS_IonsFromNeutral::cv_mole() const { + throw CanteraError("PDSS_IonsFromNeutral::cv_mole()", "unimplemented"); + return 0.0; + } + + + doublereal + PDSS_IonsFromNeutral::gibbs_RT_ref() const { + neutralMoleculePhase_->getGibbs_RT_ref(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + if (add2RTln2_) { + val += 2.0 * log(2.0); + } + return val; + } + + doublereal PDSS_IonsFromNeutral::enthalpy_RT_ref() const { + neutralMoleculePhase_->getEnthalpy_RT_ref(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + return val; + } + + doublereal PDSS_IonsFromNeutral::entropy_R_ref() const { + neutralMoleculePhase_->getEntropy_R_ref(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + if (add2RTln2_) { + val -= 2.0 * log(2.0); + } + return val; + } + + doublereal PDSS_IonsFromNeutral::cp_R_ref() const { + neutralMoleculePhase_->getCp_R_ref(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + return val; + } + + doublereal PDSS_IonsFromNeutral::molarVolume_ref() const { + neutralMoleculePhase_->getStandardVolumes_ref(DATA_PTR(tmpNM)); + doublereal val = 0.0; + for (int i = 0; i < numMult_; i++) { + int jNeut = idNeutralMoleculeVec[i]; + val += factorVec[i] * tmpNM[jNeut]; + } + return val; + } + + /* + * Calculate the pressure (Pascals), given the temperature and density + * Temperature: kelvin + * rho: density in kg m-3 + */ + doublereal PDSS_IonsFromNeutral::pressure() const { + return m_pres; + } + + void PDSS_IonsFromNeutral::setPressure(doublereal p) { + m_pres = p; + neutralMoleculePhase_->setPressure(p); + } + + + /// critical temperature + doublereal PDSS_IonsFromNeutral::critTemperature() const { + throw CanteraError("PDSS_IonsFromNeutral::critTemperature()", "unimplemented"); + return (0.0); + } + + /// critical pressure + doublereal PDSS_IonsFromNeutral::critPressure() const { + throw CanteraError("PDSS_IonsFromNeutral::critPressure()", "unimplemented"); + return (0.0); + } + + /// critical density + doublereal PDSS_IonsFromNeutral::critDensity() const { + throw CanteraError("PDSS_IonsFromNeutral::critDensity()", "unimplemented"); + return (0.0); + } + + + /* + * Return the temperature + * + * Obtain the temperature from the owning VPStandardStateTP object + * if you can. + */ + doublereal PDSS_IonsFromNeutral::temperature() const { + m_temp = m_vpssmgr_ptr->temperature(); + return m_temp; + } + + void PDSS_IonsFromNeutral::setTemperature(doublereal temp) { + m_temp = temp; + neutralMoleculePhase_->setTemperature(temp); + } + + + void PDSS_IonsFromNeutral::setState_TP(doublereal temp, doublereal pres) { + m_pres = pres; + m_temp = temp; + neutralMoleculePhase_->setState_TP(temp, pres); + } + + void PDSS_IonsFromNeutral::setState_TR(doublereal temp, doublereal rho) { + neutralMoleculePhase_->setState_TR(temp, rho); + } + + /// saturation pressure + doublereal PDSS_IonsFromNeutral::satPressure(doublereal t){ + throw CanteraError("PDSS_IonsFromNeutral::satPressure()", "unimplemented"); + /*NOTREACHED*/ + return (0.0); + } + + +} diff --git a/Cantera/src/thermo/PDSS_IonsFromNeutral.h b/Cantera/src/thermo/PDSS_IonsFromNeutral.h new file mode 100644 index 000000000..8c545355c --- /dev/null +++ b/Cantera/src/thermo/PDSS_IonsFromNeutral.h @@ -0,0 +1,462 @@ +/** + * @file PDSS_IonsFromNeutral.h + * Declarations for the class PDSS_IonsFromNeutral ( + * which handles calculations for a single ion in a fluid, whose properties + * are calculated from another neutral molecule. + * (see \ref pdssthermo and class \link Cantera::PDSS_IonsFromNeutral PDSS_IonsFromNeutral\endlink). + */ +/* + * 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: PDSS_IdealGas.h,v 1.4 2008/10/13 21:01:48 hkmoffa Exp $ + */ + +#ifndef CT_PDSS_IONSFROMNEUTRAL_H +#define CT_PDSS_IONSFROMNEUTRAL_H + +#include "PDSS.h" + + +namespace Cantera { + class XML_Node; + class VPStandardStateTP; + class ThermoPhase; + + + //! Derived class for pressure dependent standard states of an ideal gas species + /*! + * This class is for a single Ideal Gas species. + * + * @ingroup pdssthermo + */ + class PDSS_IonsFromNeutral : public PDSS { + + public: + + /** + * @name Constructors + * @{ + */ + + //! Constructor + /*! + * @param tp Pointer to the ThermoPhase object pertaining to the phase + * @param spindex Species index of the species in the phase + */ + PDSS_IonsFromNeutral(VPStandardStateTP *tp, int spindex); + + //! Copy Constructur + /*! + * @param b Object to be copied + */ + PDSS_IonsFromNeutral(const PDSS_IonsFromNeutral& b); + + //! Assignment operator + /*! + * @param b Object to be copeid + */ + PDSS_IonsFromNeutral& operator=(const PDSS_IonsFromNeutral& b); + + //! Constructor that initializes the object by examining the input file + //! of the ThermoPhase object + /*! + * This function calls the constructPDSSFile member function. + * + * @param tp Pointer to the ThermoPhase object pertaining to the phase + * @param spindex Species index of the species in the phase + * @param inputFile String name of the input file + * @param id String name of the phase in the input file. The default + * is the empty string, in which case the first phase in the + * file is used. + */ + PDSS_IonsFromNeutral(VPStandardStateTP *tp, int spindex, + std::string inputFile, std::string id = ""); + + + //! Constructor that initializes the object by examining the input file + //! of the ThermoPhase object + /*! + * This function calls the constructPDSSXML member function. + * + * @param vptp_ptr Pointer to the ThermoPhase object pertaining to the phase + * @param spindex Species index of the species in the phase + * @param speciesNode Reference to the species XML tree. + * @param phaseRef Reference to the XML tree containing the phase information. + * @param spInstalled Boolean indicating whether the species is installed yet + * or not. + */ + PDSS_IonsFromNeutral(VPStandardStateTP *vptp_ptr, int spindex, const XML_Node& speciesNode, + const XML_Node& phaseRef, bool spInstalled); + + + //! Destructor + virtual ~PDSS_IonsFromNeutral(); + + //! Duplicator + virtual PDSS *duplMyselfAsPDSS() const; + + /** + * @} + * @name Utilities + * @{ + */ + virtual int pdssType() const { return -1; } + + + /** + * @} + * @name Molar Thermodynamic Properties of the Species Standard State + * in the Solution + * @{ + */ + + //! Return the molar enthalpy in units of J kmol-1 + /*! + * Returns the species standard state enthalpy in J kmol-1 at the + * current temperature and pressure. + * + * @return returns the species standard state enthalpy in J kmol-1 + */ + virtual doublereal enthalpy_mole() const; + + //! Return the standard state molar enthalpy divided by RT + /*! + * Returns the species standard state enthalpy divided by RT at the + * current temperature and pressure. + * + * @return returns the species standard state enthalpy in unitless form + */ + virtual doublereal enthalpy_RT() const; + + //! Return the molar internal Energy in units of J kmol-1 + /*! + * Returns the species standard state internal Energy in J kmol-1 at the + * current temperature and pressure. + * + * @return returns the species standard state internal Energy in J kmol-1 + */ + virtual doublereal intEnergy_mole() const; + + //! Return the molar entropy in units of J kmol-1 K-1 + /*! + * Returns the species standard state entropy in J kmol-1 K-1 at the + * current temperature and pressure. + * + * @return returns the species standard state entropy in J kmol-1 K-1 + */ + virtual doublereal entropy_mole() const; + + //! Return the standard state entropy divided by RT + /*! + * Returns the species standard state entropy divided by RT at the + * current temperature and pressure. + * + * @return returns the species standard state entropy divided by RT + */ + virtual doublereal entropy_R() const; + + //! Return the molar gibbs free energy in units of J kmol-1 + /*! + * Returns the species standard state gibbs free energy in J kmol-1 at the + * current temperature and pressure. + * + * @return returns the species standard state gibbs free energy in J kmol-1 + */ + virtual doublereal gibbs_mole() const; + + //! Return the molar gibbs free energy divided by RT + /*! + * Returns the species standard state gibbs free energy divided by RT at the + * current temperature and pressure. + * + * @return returns the species standard state gibbs free energy divided by RT + */ + virtual doublereal gibbs_RT() const; + + //! Return the molar const pressure heat capacity in units of J kmol-1 K-1 + /*! + * Returns the species standard state Cp in J kmol-1 K-1 at the + * current temperature and pressure. + * + * @return returns the species standard state Cp in J kmol-1 K-1 + */ + virtual doublereal cp_mole() const; + + //! Return the molar const pressure heat capacity divided by RT + /*! + * Returns the species standard state Cp divided by RT at the + * current temperature and pressure. + * + * @return returns the species standard state Cp divided by RT + */ + virtual doublereal cp_R() const; + + //! Return the molar const volume heat capacity in units of J kmol-1 K-1 + /*! + * Returns the species standard state Cv in J kmol-1 K-1 at the + * current temperature and pressure. + * + * @return returns the species standard state Cv in J kmol-1 K-1 + */ + virtual doublereal cv_mole() const; + + //! Return the molar volume at standard state + /*! + * Returns the species standard state molar volume at the + * current temperature and pressure + * + * @return returns the standard state molar volume divided by R + * units are m**3 kmol-1. + */ + virtual doublereal molarVolume() const; + + //! Return the standard state density at standard state + /*! + * Returns the species standard state density at the + * current temperature and pressure + * + * @return returns the standard state density + * units are kg m-3 + */ + virtual doublereal density() const; + + /** + * @} + * @name Properties of the Reference State of the Species + * in the Solution + * @{ + */ + + //! Return the molar gibbs free energy divided by RT at reference pressure + /*! + * Returns the species reference state gibbs free energy divided by RT at the + * current temperature. + * + * @return returns the reference state gibbs free energy divided by RT + */ + virtual doublereal gibbs_RT_ref() const; + + //! Return the molar enthalpy divided by RT at reference pressure + /*! + * Returns the species reference state enthalpy divided by RT at the + * current temperature. + * + * @return returns the reference state enthalpy divided by RT + */ + virtual doublereal enthalpy_RT_ref() const; + + //! Return the molar entropy divided by R at reference pressure + /*! + * Returns the species reference state entropy divided by R at the + * current temperature. + * + * @return returns the reference state entropy divided by R + */ + virtual doublereal entropy_R_ref() const; + + //! Return the molar heat capacity divided by R at reference pressure + /*! + * Returns the species reference state heat capacity divided by R at the + * current temperature. + * + * @return returns the reference state heat capacity divided by R + */ + virtual doublereal cp_R_ref() const; + + //! Return the molar volume at reference pressure + /*! + * Returns the species reference state molar volume at the + * current temperature. + * + * @return returns the reference state molar volume divided by R + * units are m**3 kmol-1. + */ + virtual doublereal molarVolume_ref() const; + + /* + * Get the difference in the standard state thermodynamic properties + * between the reference pressure, po, and the current pressure. + */ + + /** + * @} + * @name Mechanical Equation of State Properties + * @{ + */ + + //! Returns the pressure (Pa) + virtual doublereal pressure() const; + + //! Sets the pressure in the object + /*! + * Currently, this sets the pressure in the PDSS object. + * It is indeterminant what happens to the owning VPStandardStateTP + * object and to the VPSSMgr object. + * + * @param pres Pressure to be set (Pascal) + */ + virtual void setPressure(doublereal pres); + + //! Set the internal temperature + /*! + * @param temp Temperature (Kelvin) + */ + virtual void setTemperature(doublereal temp); + + //! Return the current storred temperature + doublereal temperature() const; + + //! Set the internal temperature and pressure + /*! + * @param temp Temperature (Kelvin) + * @param pres pressure (Pascals) + */ + virtual void setState_TP(doublereal temp, doublereal pres); + + //! Set the internal temperature and density + /*! + * @param temp Temperature (Kelvin) + * @param rho Density (Pascals) + */ + virtual void setState_TR(doublereal temp, doublereal rho); + + /** + * @} + * @name Miscellaneous properties of the standard state + * @{ + */ + + /// critical temperature + virtual doublereal critTemperature() const; + + /// critical pressure + virtual doublereal critPressure() const; + + /// critical density + virtual doublereal critDensity() const; + + /// saturation pressure + /*! + * @param t Temperature (Kelvin) + */ + virtual doublereal satPressure(doublereal t); + + /** + * @} + * @name Initialization of the Object + * @{ + */ + + //! Initialization of a PDSS object using an + //! input XML file. + /*! + * + * This routine is a precursor to constructPDSSXML(XML_Node*) + * routine, which does most of the work. + * + * @param vptp_ptr Pointer to the Variable pressure %ThermoPhase object + * This object must have already been malloced. + * + * @param spindex Species index within the phase + * + * @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 constructPDSSFile(VPStandardStateTP *vptp_ptr, int spindex, + std::string inputFile, std::string id); + + //!Initialization of a PDSS object using an xml tree + /*! + * This routine is a driver for the initialization of the + * object. + * + * basic logic: + * initThermo() (cascade) + * getStuff from species Part of XML file + * initThermoXML(phaseNode) (cascade) + * + * @param vptp_ptr Pointer to the Variable pressure %ThermoPhase object + * This object must have already been malloced. + * + * @param spindex Species index within the phase + * + * @param phaseNode Reference to the phase Information for the phase + * that owns this species. + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void constructPDSSXML(VPStandardStateTP *vptp_ptr, int spindex, + const XML_Node& speciesNode, + const XML_Node& phaseNode, std::string id); + + //! Initialization routine for the PDSS object based on the phaseNode + /*! + * This is a cascading call, where each level should call the + * the parent level. + * + * @param phaseNode Reference to the phase Information for the phase + * that owns this species. + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + virtual void initThermoXML(const XML_Node& phaseNode, std::string& id); + + //! Initialization routine for all of the shallow pointers + /*! + * This is a cascading call, where each level should call the + * the parent level. + * + * The initThermo() routines get called before the initThermoXML() routines + * from the constructPDSSXML() routine. + * + * + * Calls initPtrs(); + */ + virtual void initThermo(); + + //@} + + + + protected: + + //! Maximum temperature the standard states are good for + doublereal m_tmin; + + //! Minimum temperature the standard states are good for + doublereal m_tmax; + + + //! Pointer to the Neutral Molecule thermophase object + ThermoPhase *neutralMoleculePhase_; + + public: + int numMult_; + + std::vector idNeutralMoleculeVec; + + std::vector factorVec; + + bool add2RTln2_; + + mutable std::vector tmpNM; + + int specialSpecies_; + }; +} + +#endif + + + diff --git a/Cantera/src/thermo/PseudoBinaryVPSSTP.cpp b/Cantera/src/thermo/PseudoBinaryVPSSTP.cpp index 63a416f8c..12ee677af 100644 --- a/Cantera/src/thermo/PseudoBinaryVPSSTP.cpp +++ b/Cantera/src/thermo/PseudoBinaryVPSSTP.cpp @@ -160,7 +160,6 @@ namespace Cantera { - void PseudoBinaryVPSSTP::getElectrochemPotentials(doublereal* mu) const { getChemPotentials(mu); double ve = Faraday * electricPotential(); diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index a5fa73a5c..288c6c414 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -27,6 +27,7 @@ #ifdef WITH_IDEAL_SOLUTIONS #include "IdealSolidSolnPhase.h" +#include "MargulesVPSSTP.h" #endif #ifdef WITH_PURE_FLUIDS @@ -85,13 +86,13 @@ namespace Cantera { boost::mutex ThermoFactory::thermo_mutex; #endif - static int ntypes = 16; + static int ntypes = 17; static string _types[] = {"IdealGas", "Incompressible", "Surface", "Edge", "Metal", "StoichSubstance", "PureFluid", "LatticeSolid", "Lattice", "HMW", "IdealSolidSolution", "DebyeHuckel", "IdealMolalSolution", "IdealGasVPSS", - "MineralEQ3", "electrodeElectron" + "MineralEQ3", "electrodeElectron", "Margules" }; static int _itypes[] = {cIdealGas, cIncompressible, @@ -99,7 +100,8 @@ namespace Cantera { cPureFluid, cLatticeSolid, cLattice, cHMW, cIdealSolidSolnPhase, cDebyeHuckel, cIdealMolalSoln, cVPSS_IdealGas, - cMineralEQ3, cElectrodeElectron + cMineralEQ3, cElectrodeElectron, + cMargulesVPSSTP }; /* @@ -136,6 +138,10 @@ namespace Cantera { case cIdealSolidSolnPhase: th = new IdealSolidSolnPhase(); break; + + case cMargulesVPSSTP: + th = new MargulesVPSSTP(); + break; #endif #ifdef WITH_METAL diff --git a/Cantera/src/thermo/ThermoPhase.cpp b/Cantera/src/thermo/ThermoPhase.cpp index 856d66979..215237cca 100644 --- a/Cantera/src/thermo/ThermoPhase.cpp +++ b/Cantera/src/thermo/ThermoPhase.cpp @@ -910,7 +910,6 @@ namespace Cantera { * Set the thermodynamic state. */ void ThermoPhase::setStateFromXML(const XML_Node& state) { - string comp = getChildValue(state,"moleFractions"); if (comp != "") setMoleFractionsByName(comp); diff --git a/Cantera/src/thermo/VPSSMgr_ConstVol.cpp b/Cantera/src/thermo/VPSSMgr_ConstVol.cpp index f67c35942..bf05a8640 100644 --- a/Cantera/src/thermo/VPSSMgr_ConstVol.cpp +++ b/Cantera/src/thermo/VPSSMgr_ConstVol.cpp @@ -86,6 +86,44 @@ namespace Cantera { // m_Vss[k] constant } } + + /* + * Returns the vector of nondimensional + * Gibbs free energies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * + * @param grt Output vector contains the nondimensional Gibbs free energies + * of the reference state of the species + * length = m_kk, units = dimensionless. + */ + void VPSSMgr_ConstVol::getGibbs_RT_ref(doublereal *grt) const { + if (m_useTmpRefStateStorage) { + std::copy(m_g0_RT.begin(), m_g0_RT.end(), grt); + doublereal _rt = GasConstant * m_tlast; + scale(grt, grt + m_kk, grt, _rt); + } else { + throw CanteraError("VPSSMgr_ConstVol::getGibbs_RT_ref", + "unimplemented without m_useTmpRefStateStorage"); + } + } + + + // Get the molar volumes of the species reference states at the current + // T and P_ref of the solution. + /* + * units = m^3 / kmol + * + * @param vol Output vector containing the standard state volumes. + * Length: m_kk. + */ + void VPSSMgr_ConstVol::getStandardVolumes_ref(doublereal *vol) const { + if (m_useTmpStandardStateStorage) { + std::copy(m_Vss.begin(), m_Vss.end(), vol); + } else { + throw CanteraError("VPSSMgr_ConstVol::getStandardVolumes_ref", + "unimplemented without m_useTmpRefStateStorage"); + } + } void VPSSMgr_ConstVol::initThermo() { VPSSMgr::initThermo(); diff --git a/Cantera/src/thermo/VPSSMgr_ConstVol.h b/Cantera/src/thermo/VPSSMgr_ConstVol.h index 1e5835f10..14a867530 100644 --- a/Cantera/src/thermo/VPSSMgr_ConstVol.h +++ b/Cantera/src/thermo/VPSSMgr_ConstVol.h @@ -126,6 +126,29 @@ namespace Cantera { * _updateRefStateThermo() is defined in the parent object. */ //@{ + + /*! + * Returns the vector of nondimensional + * Gibbs free energies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * + * @param grt Output vector contains the nondimensional Gibbs free energies + * of the reference state of the species + * length = m_kk, units = dimensionless. + */ + virtual void getGibbs_RT_ref(doublereal *grt) const ; + + + //! Get the molar volumes of the species reference states at the current + //! T and P_ref of the solution. + /*! + * units = m^3 / kmol + * + * @param vol Output vector containing the standard state volumes. + * Length: m_kk. + */ + virtual void getStandardVolumes_ref(doublereal *vol) const ; + //@} //! @name Initialization Methods - For Internal use diff --git a/Cantera/src/thermo/VPSSMgr_General.cpp b/Cantera/src/thermo/VPSSMgr_General.cpp index be66ec1fa..043310b8a 100644 --- a/Cantera/src/thermo/VPSSMgr_General.cpp +++ b/Cantera/src/thermo/VPSSMgr_General.cpp @@ -31,6 +31,7 @@ #include "PDSS_Water.h" #include "PDSS_ConstVol.h" #include "PDSS_HKFT.h" +#include "PDSS_IonsFromNeutral.h" #include "GeneralSpeciesThermo.h" using namespace std; @@ -175,9 +176,23 @@ namespace Cantera { "failed dynamic cast"); } genSpthermo->installPDSShandler(k, kPDSS, this); + + } else if (model == "IonFromNeutral") { + if (!genSpthermo) { + throw CanteraError("VPSSMgr_General::returnPDSS_ptr", + "failed dynamic cast"); + } + doST = false; + kPDSS = new PDSS_IonsFromNeutral(m_vptp_ptr, k, speciesNode, *phaseNode_ptr, true); + if (!kPDSS) { + throw CanteraError("VPSSMgr_General::returnPDSS_ptr", + "new PDSS_IonsFromNeutral failed"); + } + genSpthermo->installPDSShandler(k, kPDSS, this); + } else { throw CanteraError("VPSSMgr_General::returnPDSS_ptr", - "unknown"); + "unknown standard state formulation: " + model); } return kPDSS; } diff --git a/Cantera/src/thermo/VPStandardStateTP.cpp b/Cantera/src/thermo/VPStandardStateTP.cpp index b4b523ee0..5bd0a5a7e 100644 --- a/Cantera/src/thermo/VPStandardStateTP.cpp +++ b/Cantera/src/thermo/VPStandardStateTP.cpp @@ -338,8 +338,13 @@ namespace Cantera { void VPStandardStateTP::setTemperature(doublereal t) { - State::setTemperature(t); - //updateStandardStateThermo(); + setState_TP(t, m_Pcurrent); + updateStandardStateThermo(); + } + + void VPStandardStateTP::setPressure(doublereal p) { + setState_TP(temperature(), p); + updateStandardStateThermo(); } @@ -359,12 +364,13 @@ namespace Cantera { /* * Now, we still need to do the calculations for general ThermoPhase objects. * So, we switch back to a virtual function call, setTemperature, and - * setPressure to recalculate stuff at the higher level. At this point, + * setPressure to recalculate stuff for child ThermoPhase objects of + * the VPStandardStateTP object. At this point, * we haven't touched m_tlast or m_plast, so some calculations may still * need to be done at the ThermoPhase object level. */ - setTemperature(t); - setPressure(pres); + //setTemperature(t); + //setPressure(pres); } diff --git a/Cantera/src/thermo/VPStandardStateTP.h b/Cantera/src/thermo/VPStandardStateTP.h index 6b81889bb..72f7cfebd 100644 --- a/Cantera/src/thermo/VPStandardStateTP.h +++ b/Cantera/src/thermo/VPStandardStateTP.h @@ -45,7 +45,7 @@ namespace Cantera { * for which VPStandardStateTP owns a pointer to. * * To support the above functionality, pressure and temperature variables, - * m_plast_ss and m_tlast_ss, are kept which store the last pressure and temperature + * m_Plast_ss and m_Tlast_ss, are kept which store the last pressure and temperature * used in the evaluation of standard state properties. * * This class is usually used for nearly incompressible phases. For those phases, it @@ -257,18 +257,30 @@ namespace Cantera { //! Set the temperature of the phase /*! - * Currently this just passes down to State::setTemperature() - * without doing anything. Calculations are changing temperatures are triggered - * later. + * Currently this passes down to setState_TP(). It does not + * make sense to calculate the standard state without first + * setting T and P. * * @param T Temperature (kelvin) */ virtual void setTemperature(const doublereal T); + + //! Set the internally storred pressure (Pa) at constant + //! temperature and composition + /*! + * Currently this passes down to setState_TP(). It does not + * make sense to calculate the standard state without first + * setting T and P. + * + * @param p input Pressure (Pa) + */ + virtual void setPressure(const doublereal p); + //! Set the temperature and pressure at the same time /*! - * Note this function currently triggers a reevalulation of the standard + * Note this function triggers a reevalulation of the standard * state quantities. * * @param T temperature (kelvin) diff --git a/Cantera/src/thermo/mix_defs.h b/Cantera/src/thermo/mix_defs.h index fef148800..cd5f6f3df 100755 --- a/Cantera/src/thermo/mix_defs.h +++ b/Cantera/src/thermo/mix_defs.h @@ -70,7 +70,9 @@ namespace Cantera { const int cIdealSolnGasVPSS = 500; const int cIdealSolnGasVPSS_iscv = 501; + const int cMargulesVPSSTP = 301; + const int cIonFromNeutral = 2000; //! Variable Pressure Standard State ThermoPhase objects const int cVPSS_IdealGas = 1001; @@ -87,7 +89,8 @@ namespace Cantera { cPDSS_CONSTVOL, cPDSS_MOLAL_CONSTVOL, cPDSS_WATER, - cPDSS_MOLAL_HKFT + cPDSS_MOLAL_HKFT, + cPDSS_IONSFROMNEUTRAL };