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 };