research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Added a few more files to handle liquid electrochemistry thermo.

Browse source
This commit is contained in:
Harry Moffat 2005-11-14 18:49:56 +00:00
parent bcd9750364
commit 5a9ad810eb
6 changed files with 1703 additions and 3 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
Cantera/src/thermo/Makefile.in
View file

@ -18,8 +18,10 @@ do_ranlib = @DO_RANLIB@
CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
# Extended Cantera Thermodynamics Object Files
CATHERMO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o
CATHERMO_H = SingleSpeciesTP.h StoichSubstanceSSTP.h
CATHERMO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o \
MolalityVPSSTP.o VPStandardStateTP.o
CATHERMO_H = SingleSpeciesTP.h StoichSubstanceSSTP.h \
MolalityVPSSTP.h VPStandardStateTP.h
CXX_INCLUDES = -I.. @CXX_INCLUDES@
LIB = @buildlib@/libcaThermo.a

505
Cantera/src/thermo/MolalityVPSSTP.cpp Normal file
View file

@ -0,0 +1,505 @@
/**
*
* @file MolalityVPSSTP.cpp
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
/*
* $Author$
* $Date$
* $Revision$
*/
#ifndef MAX
#define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
#endif
#include "MolalityVPSSTP.h"
namespace Cantera {
/*
* Default constructor.
*
* This doesn't do much more than initialize constants with
* default values for water at 25C.
*/
MolalityVPSSTP::MolalityVPSSTP() :
VPStandardStateTP(),
m_indexSolvent(0),
m_weightSolvent(18.0),
m_xmolSolventMIN(0.01),
m_Mnaught(18.0E-3)
{
}
/**
* Copy Constructor:
*
* Note this stuff will not work until the underlying phase
* has a working copy constructor
*/
MolalityVPSSTP::MolalityVPSSTP(const MolalityVPSSTP &b) :
VPStandardStateTP(),
m_indexSolvent(b.m_indexSolvent),
m_xmolSolventMIN(b.m_xmolSolventMIN),
m_Mnaught(b.m_Mnaught),
m_molalities(b.m_molalities)
{
throw CanteraError("MolalityVPSSTP::operator=()",
"Not Implemented Fully");
*this = operator=(b);
}
/*
* operator=()
*
* Note this stuff will not work until the underlying phase
* has a working assignment operator
*/
MolalityVPSSTP& MolalityVPSSTP::
operator=(const MolalityVPSSTP &b) {
if (&b != this) {
VPStandardStateTP::operator=(b);
m_indexSolvent = b.m_indexSolvent;
m_weightSolvent = b.m_weightSolvent;
m_xmolSolventMIN = b.m_xmolSolventMIN;
m_Mnaught = b.m_Mnaught;
m_molalities = b.m_molalities;
}
throw CanteraError("MolalityVPSSTP::operator=()",
"Not Implemented Fully");
return *this;
}
/**
*
* ~MolalityVPSSTP(): (virtual)
*
* Destructor: does nothing:
*
*/
MolalityVPSSTP::~MolalityVPSSTP() {
}
/**
* This routine duplicates the current object and returns
* a pointer to ThermoPhase.
*/
ThermoPhase*
MolalityVPSSTP::duplMyselfAsThermoPhase() {
MolalityVPSSTP* mtp = new MolalityVPSSTP(*this);
return (ThermoPhase *) mtp;
}
/*
* -------------- Utilities -------------------------------
*/
/*
* setSolvent():
* Utilities for Solvent ID and Molality
* Here we also calculate and store the molecular weight
* of the solvent and the m_Mnaught parameter.
*/
void MolalityVPSSTP::setSolvent(int k) {
if (k < 0 || k >= m_kk) {
throw CanteraError("MolalityVPSSTP::setSolute ", "trouble");
}
m_indexSolvent = k;
m_weightSolvent = molecularWeight(k);
m_Mnaught = m_weightSolvent / 1000.;
}
/*
* return the solvent id index number.
*/
int MolalityVPSSTP::solventIndex() const {
return m_indexSolvent;
}
/*
* Sets the minimum mole fraction in the molality formulation
*/
void MolalityVPSSTP::
setMoleFSolventMin(doublereal xmolSolventMIN) {
if (xmolSolventMIN <= 0.0) {
throw CanteraError("MolalityVPSSTP::setSolute ", "trouble");
} else if (xmolSolventMIN > 0.9) {
throw CanteraError("MolalityVPSSTP::setSolute ", "trouble");
}
m_xmolSolventMIN = xmolSolventMIN;
}
/*
* Returns the minimum mole fraction in the molality
* formulation.
*/
doublereal MolalityVPSSTP::moleFSolventMin() const {
return m_xmolSolventMIN;
}
/**
* getMolalities():
* We calculate the vector of molalities of the species
* in the phase
* m_i = (n_i) / (1000 * M_o * n_o_p)
*
* where M_o is the molecular weight of the solvent
* n_o is the mole fraction of the solvent
* n_i is the mole fraction of the solute.
* n_o_p = max (n_o_min, n_o)
* n_o_min = minimum mole fraction of solvent allowed
* in the denominator.
*/
void MolalityVPSSTP::getMolalities(doublereal * const molal) const {
getMoleFractions(molal);
double xmolSolvent = molal[m_indexSolvent];
if (xmolSolvent < m_xmolSolventMIN) {
xmolSolvent = m_xmolSolventMIN;
}
double denomInv = 1.0/
(m_Mnaught * xmolSolvent);
for (int k = 0; k < m_kk; k++) {
molal[k] *= denomInv;
}
for (int k = 0; k < m_kk; k++) {
m_molalities[k] = molal[k];
}
}
/**
* setMolalities():
* We are supplied with the molalities of all of the
* solute species. We then calculate the mole fractions of all
* species and update the ThermoPhase object.
*
* m_i = (n_i) / (W_o/1000 * n_o_p)
*
* where M_o is the molecular weight of the solvent
* n_o is the mole fraction of the solvent
* n_i is the mole fraction of the solute.
* n_o_p = max (n_o_min, n_o)
* n_o_min = minimum mole fraction of solvent allowed
* in the denominator.
*/
void MolalityVPSSTP::setMolalities(const doublereal * const molal) {
double Lsum = 1.0 / m_Mnaught;
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
m_molalities[k] = molal[k];
Lsum += molal[k];
}
}
double tmp = 1.0 / Lsum;
m_molalities[m_indexSolvent] = tmp / m_Mnaught;
double sum = m_molalities[m_indexSolvent];
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
m_molalities[k] = tmp * molal[k];
sum += m_molalities[k];
}
}
if (sum != 1.0) {
tmp = 1.0 / sum;
for (int k = 0; k < m_kk; k++) {
m_molalities[k] *= tmp;
}
}
setMoleFractions(m_molalities.begin());
/*
* Essentially we don't trust the input: We calculate
* the molalities from the mole fractions that we
* just obtained.
*/
getMolalities(m_molalities.begin());
}
/*
* setMolalitiesByName()
*
* This routine sets the molalities by name
* HKM -> Might need to be more complicated here, setting
* neutrals so that the existing mole fractions are
* preserved.
*/
void MolalityVPSSTP::setMolalitiesByName(compositionMap& mMap) {
int kk = nSpecies();
doublereal x;
/*
* Get a vector of mole fractions
*/
vector_fp mf(kk, 0.0);
getMoleFractions(mf.begin());
double xmolS = mf[m_indexSolvent];
double xmolSmin = max(xmolS, m_xmolSolventMIN);
compositionMap::iterator p;
for (int k = 0; k < kk; k++) {
p = mMap.find(speciesName(k));
if (p != mMap.end()) {
x = mMap[speciesName(k)];
if (x > 0.0) {
mf[k] = x * m_Mnaught * xmolSmin;
}
}
}
/*
* check charge neutrality
*/
int largePos = -1;
double cPos = 0.0;
int largeNeg = -1;
double cNeg = 0.0;
double sum = 0.0;
for (int k = 0; k < kk; k++) {
double ch = charge(k);
if (mf[k] > 0.0) {
if (ch > 0.0) {
if (ch * mf[k] > cPos) {
largePos = k;
cPos = ch * mf[k];
}
}
if (ch < 0.0) {
if (fabs(ch) * mf[k] > cNeg) {
largeNeg = k;
cNeg = fabs(ch) * mf[k];
}
}
}
sum += mf[k] * ch;
}
if (sum != 0.0) {
if (sum > 0.0) {
if (cPos > sum) {
mf[largePos] -= sum / charge(largePos);
} else {
throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
"unbalanced charges");
}
} else {
if (cNeg > (-sum)) {
mf[largeNeg] -= (-sum) / fabs(charge(largeNeg));
} else {
throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
"unbalanced charges");
}
}
}
sum = 0.0;
for (int k = 0; k < kk; k++) {
sum += mf[k];
}
sum = 1.0/sum;
for (int k = 0; k < kk; k++) {
mf[k] *= sum;
}
setMoleFractions(mf.begin());
/*
* After we formally set the mole fractions, we
* calculate the molalities again and store it in
* this object.
*/
getMolalities(m_molalities.begin());
}
/*
* setMolalitiesByNames()
*
* Set the molalities of the solutes by name
*/
void MolalityVPSSTP::setMolalitiesByName(const string& x) {
compositionMap xx;
int kk = nSpecies();
for (int k = 0; k < kk; k++) {
xx[speciesName(k)] = -1.0;
}
parseCompString(x, xx);
setMolalitiesByName(xx);
}
/*
* Update the internal array that contains the molalities of the
* species.
*/
void MolalityVPSSTP::updateMolalities() const {
getMolalities(m_molalities.begin());
}
/*
* ------------ Molar Thermodynamic Properties ----------------------
*/
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
/**
* This method returns the activity convention.
* Currently, there are two activity conventions
* Molar-based activities
* Unit activity of species at either a hypothetical pure
* solution of the species or at a hypothetical
* pure ideal solution at infinite dilution
* cAC_CONVENTION_MOLAR 0
* - default
*
* Molality based activities
* (unit activity of solutes at a hypothetical 1 molal
* solution referenced to infinite dilution at all
* pressures and temperatures).
* (solvent is still on molar basis).
* cAC_CONVENTION_MOLALITY 1
*
* We set the convention to molality here.
*/
int MolalityVPSSTP::activityConvention() const {
return cAC_CONVENTION_MOLALITY;
}
/**
* Get the array of non-dimensional activity coefficients at
* the current solution temperature, pressure, and
* solution concentration.
* These are mole fraction based activity coefficients. In this
* object, their calculation is based on translating the values
* of Molality based activity coefficients.
* See Denbigh p. 278 for a thorough discussion.
*
* Note, the solvent is treated differently. getMolalityActivityCoeff()
* returns the molar based solvent activity coefficient already.
* Therefore, we do not have to divide by x_s here.
*/
void MolalityVPSSTP::getActivityCoefficients(doublereal* ac) const {
getMolalityActivityCoefficients(ac);
double xmolSolvent = moleFraction(m_indexSolvent);
if (xmolSolvent < m_xmolSolventMIN) {
xmolSolvent = m_xmolSolventMIN;
}
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
ac[k] /= xmolSolvent;
}
}
}
/**
* osmotic coefficient:
*
* Calculate the osmotic coefficient of the solvent. Note there
* are lots of definitions of the osmotic coefficient floating
* around. We use the one defined in the Pitzer paper:
*
* Definition:
* - sum(m_i) * M0 * oc = ln(activity_solvent)
*/
doublereal MolalityVPSSTP::osmoticCoefficient() const {
vector_fp act(m_kk);
getActivities(act.begin());
double sum = 0;
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
sum += MAX(m_molalities[k], 0.0);
}
}
double oc = 1.0;
double lac = log(act[m_indexSolvent]);
if (sum > 1.0E-200) {
oc = - lac / (m_Mnaught * sum);
}
return oc;
}
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
doublereal MolalityVPSSTP::err(string msg) const {
throw CanteraError("MolalityVPSSTP","Base class method "
+msg+" called. Equation of state type: "+int2str(eosType()));
return 0;
}
/**
* Returns the units of the standard and general concentrations
* Note they have the same units, as their divisor is
* defined to be equal to the activity of the kth species
* in the solution, which is unitless.
*
* This routine is used in print out applications where the
* units are needed. Usually, MKS units are assumed throughout
* the program and in the XML input files.
*
* On return uA contains the powers of the units (MKS assumed)
* of the standard concentrations and generalized concentrations
* for the kth species.
*
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
*/
void MolalityVPSSTP::getUnitsStandardConc(double *uA, int k, int sizeUA) {
for (int i = 0; i < sizeUA; i++) {
if (i == 0) uA[0] = 1.0;
if (i == 1) uA[1] = -nDim();
if (i == 2) uA[2] = 0.0;
if (i == 3) uA[3] = 0.0;
if (i == 4) uA[4] = 0.0;
if (i == 5) uA[5] = 0.0;
}
}
/*
* Set the thermodynamic state.
*/
void MolalityVPSSTP::setStateFromXML(const XML_Node& state) {
VPStandardStateTP::setStateFromXML(state);
string comp = getString(state,"soluteMolalities");
if (comp != "") {
setMolalitiesByName(comp);
}
if (state.hasChild("pressure")) {
double p = getFloat(state, "pressure", "pressure");
setPressure(p);
}
}
/**
* @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 MolalityVPSSTP::initThermo() {
VPStandardStateTP::initThermo();
m_molalities.resize(m_kk);
}
}

455
Cantera/src/thermo/MolalityVPSSTP.h Normal file
View file

@ -0,0 +1,455 @@
/**
* @file MolalityVPSSTP.h
*
* 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 (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
/*
* $Author$
* $Date$
* $Revision$
*/
#ifndef CT_MOLALITYVPSSTP_H
#define CT_MOLALITYVPSSTP_H
#include "VPStandardStateTP.h"
namespace Cantera {
/**
* @ingroup thermoprops
*/
/**
* MolalityVPSSTP is 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.
*/
class MolalityVPSSTP : public VPStandardStateTP {
public:
/// Constructors
MolalityVPSSTP();
MolalityVPSSTP(const MolalityVPSSTP &);
/// Assignment operator
MolalityVPSSTP& operator=(const MolalityVPSSTP&);
/// Destructor.
virtual ~MolalityVPSSTP();
/**
* 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();
/**
*
* @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 { return 0; }
/**
* @}
* @name Molar Thermodynamic Properties
* @{
*/
/**
* @}
* @name Utilities for Solvent ID and Molality
* @{
*/
/**
* This routine sets the index number of the solvent for
* the phase.
*
* Note, having a solvent
* is a precursor to many things having to do with molality.
*
* @param k the solvent index number
*/
void setSolvent(int k);
/**
* Sets the minimum mole fraction in the molality formulation.
* Note the molality formulation is singular in the limit that
* the solvent mole fraction goes to zero. Numerically, how
* this limit is treated and resolved is an ongoing issue within
* Cantera.
*/
void setMoleFSolventMin(doublereal xmolSolventMIN);
/**
* Returns the solvent index.
*/
int solventIndex() const;
/**
* Returns the minimum mole fraction in the molality
* formulation.
*/
doublereal moleFSolventMin() const;
/**
* getMolalities()
* This function will return the molalities of the
* species.
*
*/
void getMolalities(doublereal * const molal) const;
void setMolalities(const doublereal * const molal);
void setMolalitiesByName(compositionMap& xMap);
void setMolalitiesByName(const string &);
void updateMolalities() const;
/**
* @}
* @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 the activity convention.
* Currently, there are two activity conventions
* Molar-based activities
* Unit activity of species at either a hypothetical pure
* solution of the species or at a hypothetical
* pure ideal solution at infinite dilution
* cAC_CONVENTION_MOLAR 0
* - default
*
* Molality based acvtivities
* (unit activity of solutes at a hypothetical 1 molal
* solution referenced to infinite dilution at all
* pressures and temperatures).
* cAC_CONVENTION_MOLALITY 1
*
* We set the convention to molality here.
*/
int activityConvention() const;
/**
* This method returns an array of generalized concentrations
* \f$ C_k\f$ that are defined such that
* \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$
* is a standard concentration
* defined below. These generalized concentrations are used
* by kinetics manager classes to compute the forward and
* reverse rates of elementary reactions.
*
* @param c Array of generalized concentrations. The
* units depend upon the implementation of the
* reaction rate expressions within the phase.
*/
virtual void getActivityConcentrations(doublereal* c) const {
err("getActivityConcentrations");
}
/**
* 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.
*/
virtual doublereal standardConcentration(int k=0) const {
err("standardConcentration");
return -1.0;
}
/**
* Returns the natural logarithm of the standard
* concentration of the kth species
*/
virtual doublereal logStandardConc(int k=0) const {
err("logStandardConc");
return -1.0;
}
/**
* Returns the units of the standard and generalized
* concentrations Note they have the same units, as their
* ratio is defined to be equal to the activity of the kth
* species in the solution, which is unitless.
*
* This routine is used in print out applications where the
* units are needed. Usually, MKS units are assumed throughout
* the program and in the XML input files.
*
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
*/
virtual void getUnitsStandardConc(double *uA, int k = 0,
int sizeUA = 6);
/**
* Get the array of non-dimensional activities (molality
* based for this class and classes that derive from it) at
* the current solution temperature, pressure, and
* solution concentration.
*/
virtual void getActivities(doublereal* ac) const {
err("getActivities");
}
/**
* Get the array of non-dimensional activity coefficients at
* the current solution temperature, pressure, and
* solution concentration.
* These are mole fraction based activity coefficients. In this
* object, their calculation is based on translating the values
* of Molality based activity coefficients.
* See Denbigh p. 278 for a thorough discussion
*/
void getActivityCoefficients(doublereal* ac) const;
/**
* Get the array of non-dimensional molality based
* activity coefficients at the current solution temperature,
* pressure, and solution concentration.
* See Denbigh p. 278 for a thorough discussion
*/
virtual void getMolalityActivityCoefficients(doublereal *acMolality)
const {
err("getMolalityActivityCoefficients");
}
/**
* Calculate the osmotic coefficient
* units = dimensionless
*/
virtual double osmoticCoefficient() const;
//@}
/// @name Partial Molar Properties of the Solution
//@{
/**
* Get the species electrochemical potentials.
* These are partial molar quantities.
* This method adds a term \f$ Fz_k \phi_k \f$ to the
* to each chemical potential.
*
* Units: J/kmol
*/
void getElectrochemPotentials(doublereal* mu) const {
getChemPotentials(mu);
double ve = Faraday * electricPotential();
for (int k = 0; k < m_kk; k++) {
mu[k] += ve*charge(k);
}
}
//@}
/// @name Properties of the Standard State of the Species in the Solution
//@{
//@}
/// @name Thermodynamic Values for the Species Reference States
//@{
///////////////////////////////////////////////////////
//
// The methods below are not virtual, and should not
// be overloaded.
//
//////////////////////////////////////////////////////
/**
* @name Specific Properties
* @{
*/
/**
* @name Setting the State
*
* These methods set all or part of the thermodynamic
* state.
* @{
*/
//@}
/**
* @name Chemical Equilibrium
* Routines that implement the Chemical equilibrium capability
* for a single phase, based on the element-potential method.
* @{
*/
/**
* This method is used by the ChemEquil element-potential
* based equilibrium solver.
* It sets the state such that the chemical potentials of the
* species within the current phase satisfy
* \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m}
* \left(\frac{\lambda_m} {\hat R T}\right) \f] where
* \f$ \lambda_m \f$ is the element potential of element m. The
* temperature is unchanged. Any phase (ideal or not) that
* implements this method can be equilibrated by ChemEquil.
*/
virtual void setToEquilState(const doublereal* lambda_RT) {
err("setToEquilState");
}
// called by function 'equilibrate' in ChemEquil.h to transfer
// the element potentials to this object
void setElementPotentials(const vector_fp& lambda) {
m_lambda = lambda;
}
void getElementPotentials(doublereal* lambda) {
copy(m_lambda.begin(), m_lambda.end(), lambda);
}
//@}
/**
* Set equation of state parameter values from XML
* entries. This method is called by function importPhase in
* file importCTML.cpp when processing a phase definition in
* an input file. It should be overloaded in subclasses to set
* any parameters that are specific to that particular phase
* model.
*
* The MolalityVPSSTP object defines a new method for setting
* the concentrations of a phase. The new method is defined by a
* block called "soluteMolalities". If this block
* is found, the concentrations within that phase are
* set to the "name":"molalities pairs found within that
* XML block. The solvent concentration is then set
* to everything else.
*
* @param eosdata An XML_Node object corresponding to
* the "thermo" entry for this phase in the input file.
*
*/
virtual void setStateFromXML(const XML_Node& state);
/// The following methods are used in the process of constructing
/// the phase and setting its parameters from a specification in an
/// input file. They are not normally used in application programs.
/// To see how they are used, see files importCTML.cpp and
/// ThermoFactory.cpp.
/**
* @internal Initialize. This method is provided to allow
* subclasses to perform any initialization required after all
* species have been added. For example, it might be used to
* resize internal work arrays that must have an entry for
* each species. The base class implementation does nothing,
* and subclasses that do not require initialization do not
* need to overload this method. When importing a CTML phase
* description, this method is called just prior to returning
* from function importPhase.
*
* @see importCTML.cpp
*/
virtual void initThermo();
protected:
int m_indexSolvent;
doublereal m_weightSolvent;
/*
* In any molality implementation, it makes sense to have
* a minimum solvent mole fraction requirement, since the
* implementation becomes singular in the xmolSolvent=0
* limit. The default is to set it to 0.01.
* We then modify the molality definition to ensure that
* molal_solvent = 0 when xmol_solvent = 0.
*/
doublereal m_xmolSolventMIN;
/*
* This is the multiplication factor that goes inside
* log expressions involving the molalities of species.
* Its equal to Wt_0 / 1000.
* where Wt_0 = weight of solvent (kg/kmol)
*/
doublereal m_Mnaught;
mutable vector_fp m_molalities;
private:
doublereal err(string msg) const;
};
}
#endif

1
Cantera/src/thermo/StoichSubstanceSSTP.h
View file

@ -108,7 +108,6 @@ namespace Cantera {
*/
virtual doublereal thermalExpansionCoeff() const ;
//@}
/**
* @}

300
Cantera/src/thermo/VPStandardStateTP.cpp Normal file
View file

@ -0,0 +1,300 @@
/**
*
* @file VPStandardStateTP.cpp
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
/*
* $Author$
* $Date$
* $Revision$
*/
// turn off warnings under Windows
#ifdef WIN32
#pragma warning(disable:4786)
#pragma warning(disable:4503)
#endif
#include "VPStandardStateTP.h"
namespace Cantera {
/*
* Default constructor
*/
VPStandardStateTP::VPStandardStateTP() :
ThermoPhase(),
m_tlast(-1.0)
{
}
/*
* Copy Constructor:
*
* Note this stuff will not work until the underlying phase
* has a working copy constructor.
*
* The copy constructor just calls the assignment operator
* to do the heavy lifting.
*/
VPStandardStateTP::VPStandardStateTP(const VPStandardStateTP &b) :
ThermoPhase(),
m_tlast(-1.0)
{
*this = b;
}
/*
* operator=()
*
* Note this stuff will not work until the underlying phase
* has a working assignment operator
*/
VPStandardStateTP& VPStandardStateTP::
operator=(const VPStandardStateTP &b) {
if (&b != this) {
/*
* Mostly, this is a passthrough to the underlying
* assignment operator for the ThermoPhae parent object.
*/
ThermoPhase::operator=(b);
/*
* However, we have to handle data that we own.
*/
m_tlast = b.m_tlast;
m_h0_RT = b.m_h0_RT;
m_cp0_R = b.m_cp0_R;
m_g0_RT = b.m_g0_RT;
m_s0_R = b.m_s0_R;
}
return *this;
}
/*
* ~VPStandardStateTP(): (virtual)
*
* This destructor does nothing. All of the owned objects
* handle themselves.
*/
VPStandardStateTP::~VPStandardStateTP() {
}
/*
* Duplication function.
* This calls the copy constructor for this object.
*/
ThermoPhase* VPStandardStateTP::duplMyselfAsThermoPhase() {
VPStandardStateTP* vptp = new VPStandardStateTP(*this);
return (ThermoPhase *) vptp;
}
/*
* -------------- Utilities -------------------------------
*/
/*
* ------------Molar Thermodynamic Properties -------------------------
*/
doublereal VPStandardStateTP::err(string msg) const {
throw CanteraError("VPStandardStateTP","Base class method "
+msg+" called. Equation of state type: "+int2str(eosType()));
return 0;
}
/**
* Returns the units of the standard and general concentrations
* Note they have the same units, as their divisor is
* defined to be equal to the activity of the kth species
* in the solution, which is unitless.
*
* This routine is used in print out applications where the
* units are needed. Usually, MKS units are assumed throughout
* the program and in the XML input files.
*
* On return uA contains the powers of the units (MKS assumed)
* of the standard concentrations and generalized concentrations
* for the kth species.
*
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
*/
void VPStandardStateTP::
getUnitsStandardConc(double *uA, int k, int sizeUA) {
for (int i = 0; i < sizeUA; i++) {
if (i == 0) uA[0] = 1.0;
if (i == 1) uA[1] = -nDim();
if (i == 2) uA[2] = 0.0;
if (i == 3) uA[3] = 0.0;
if (i == 4) uA[4] = 0.0;
if (i == 5) uA[5] = 0.0;
}
}
/*
* ---- Partial Molar Properties of the Solution -----------------
*/
/**
* Get the array of non-dimensional species chemical potentials
* These are partial molar Gibbs free energies.
* \f$ \mu_k / \hat R T \f$.
* Units: unitless
*
* We close the loop on this function, here, calling
* getChemPotentials() and then dividing by RT.
*/
void VPStandardStateTP::getChemPotentials_RT(doublereal* muRT) const{
getChemPotentials(muRT);
doublereal invRT = 1.0 / _RT();
for (int k = 0; k < m_kk; k++) {
muRT[k] *= invRT;
}
}
/*
* ----- Thermodynamic Values for the Species Reference States ----
*/
/**
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
*/
void VPStandardStateTP::getEnthalpy_RT_ref(doublereal *hrt) const {
/*
* Call the function that makes sure the local copy of
* the species reference thermo functions are up to date
* for the current temperature.
*/
_updateRefStateThermo();
/*
* Copy the enthalpy function into return vector.
*/
copy(m_h0_RT.begin(), m_h0_RT.end(), hrt);
}
/**
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
*/
void VPStandardStateTP::getGibbs_RT_ref(doublereal *grt) const {
/*
* Call the function that makes sure the local copy of
* the species reference thermo functions are up to date
* for the current temperature.
*/
_updateRefStateThermo();
/*
* Copy the gibbs function into return vector.
*/
copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
}
/**
* Returns the vector of the
* gibbs function of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* units = J/kmol
*
* This is filled in here so that derived classes don't have to
* take care of it.
*/
void VPStandardStateTP::getGibbs_ref(doublereal *g) const {
getGibbs_RT_ref(g);
double RT = _RT();
for (int k = 0; k < m_kk; k++) {
g[k] *= RT;
}
}
/**
* Returns the vector of nondimensional
* entropies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
*/
void VPStandardStateTP::getEntropy_R_ref(doublereal *er) const {
/*
* Call the function that makes sure the local copy of
* the species reference thermo functions are up to date
* for the current temperature.
*/
_updateRefStateThermo();
/*
* Copy the gibbs function into return vector.
*/
copy(m_s0_R.begin(), m_s0_R.end(), er);
}
/**
* Returns the vector of nondimensional
* constant pressure heat capacities of the reference state
* at the current temperature of the solution
* and reference pressure for the species.
*/
void VPStandardStateTP::getCp_R_ref(doublereal *cpr) const {
/*
* Call the function that makes sure the local copy of
* the species reference thermo functions are up to date
* for the current temperature.
*/
_updateRefStateThermo();
/*
* Copy the gibbs function into return vector.
*/
copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
}
/**
* Perform initializations after all species have been
* added.
*/
void VPStandardStateTP::initThermo() {
ThermoPhase::initThermo();
m_kk = nSpecies();
int leng = m_kk;
m_h0_RT.resize(leng);
m_g0_RT.resize(leng);
m_cp0_R.resize(leng);
m_s0_R.resize(leng);
}
/**
* void _updateRefStateThermo() (private, const)
*
* This function gets called for every call to functions in this
* class. It checks to see whether the temperature has changed and
* thus the reference thermodynamics functions for all of the species
* must be recalculated.
* If the temperature has changed, the species thermo manager is called
* to recalculate G, Cp, H, and S at the current temperature.
*/
void VPStandardStateTP::_updateRefStateThermo() const {
doublereal tnow = temperature();
if (m_tlast != tnow) {
m_spthermo->update(tnow, m_cp0_R.begin(), m_h0_RT.begin(),
m_s0_R.begin());
m_tlast = tnow;
for (int k = 0; k < m_kk; k++) {
m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
}
}
}
}

439
Cantera/src/thermo/VPStandardStateTP.h Normal file
View file

@ -0,0 +1,439 @@
/**
* @file VPStandardStateTP.h
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties. These include most of the
* methods for calculating liquid electrolyte thermodynamics.
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
/*
* $Author$
* $Date$
* $Revision$
*/
#ifndef CT_VPSTANDARDSTATETP_H
#define CT_VPSTANDARDSTATETP_H
#include "ThermoPhase.h"
namespace Cantera {
class XML_Node;
/**
* @ingroup thermoprops
*
* This is a filter class for ThermoPhase that implements
* a variable pressure standard state for ThermoPhase objects.
*
* In addition support for the molality unit scale is provided.
*
* Currently, it really is just a shell. The ThermoPhase object
* itself is based around the general concepts of
* VPStandardStateTP. Therefore, there really isn't much going
* on here.
* However, this may change. The ThermoPhase object itself
* could change. Additionally, this object may revolve around
* the molality unit scale in the near future. We will have to see
* how things fare.
*/
class VPStandardStateTP : public ThermoPhase {
public:
/// Constructor.
VPStandardStateTP();
/// Copy Constructor.
VPStandardStateTP(const VPStandardStateTP &);
/// Assignment operator
VPStandardStateTP& operator=(const VPStandardStateTP &);
/// Destructor.
virtual ~VPStandardStateTP();
/*
* Duplication routine
*/
virtual ThermoPhase *duplMyselfAsThermoPhase();
/**
*
* @name Utilities
* @{
*/
/**
* Equation of state type flag. The base class returns
* zero. Subclasses should define this to return a unique
* non-zero value. Constants defined for this purpose are
* listed in mix_defs.h.
*/
virtual int eosType() const { return 0; }
/**
* @}
* @name Molar Thermodynamic Properties of the Solution
* @{
*/
/*
* These are handled by inherited objects. At this level,
* this pass-through routine doesn't add anything to the
* ThermoPhase description.
*/
/**
* @}
* @name Mechanical Properties
* @{
*/
/*
* These are handled by inherited objects. At this level,
* this pass-through routine doesn't add anything to the
* ThermoPhase description.
*/
/**
* @}
* @name Electric Potential
*
* The phase may be at some non-zero electrical
* potential. These methods set or get the value of the
* electric potential.
* @{
*/
/*
* These are handled by inherited objects. At this level,
* this pass-through routine doesn't add anything to the
* ThermoPhase description.
*/
/**
* @}
* @name Activities 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)\f$ is
* the chemical potential at unit activity, which depends only
* on temperature.
* @{
*/
/**
* Returns the units of the standard and generalized
* concentrations Note they have the same units, as their
* ratio is defined to be equal to the activity of the kth
* species in the solution, which is unitless.
*
* This routine is used in print out applications where the
* units are needed. Usually, MKS units are assumed throughout
* the program and in the XML input files.
*
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
*/
virtual void getUnitsStandardConc(double *uA, int k = 0,
int sizeUA = 6);
//@}
/// @name Partial Molar Properties of the Solution
//@{
/**
* Get the array of non-dimensional species chemical potentials
* These are partial molar Gibbs free energies.
* \f$ \mu_k / \hat R T \f$.
* Units: unitless
*
* We close the loop on this function, here, calling
* getChemPotentials() and then dividing by RT.
*/
virtual void getChemPotentials_RT(doublereal* mu) const;
//@}
/// @name Properties of the Standard State of the Species in the Solution
//@{
/*
* These are handled by inherited objects. At this level,
* this pass-through routine doesn't add anything to the
* ThermoPhase description.
*
* However, we assume these methods exist for inherited objects.
* Therefore, we will bring the error routines up to this object
*/
/**
* Get the array of chemical potentials at unit activity.
* These
* are the standard state chemical potentials \f$ \mu^0_k(T,P)
* \f$.. The values are evaluated at the current
* temperature and pressure.
*/
virtual void getStandardChemPotentials(doublereal* mu) const {
err("getStandardChemPotentials");
}
/**
* Get the nondimensional Enthalpy functions for the species
* at their standard states at the current
* <I>T</I> and <I>P</I> of the solution.
*/
virtual void getEnthalpy_RT(doublereal* hrt) const {
err("getEnthalpy_RT");
}
/**
* Get the array of nondimensional Enthalpy functions for the
* standard state species
* at the current <I>T</I> and <I>P</I> of the solution.
*/
virtual void getEntropy_R(doublereal* sr) const {
err("getEntropy_R");
}
/**
* Get the nondimensional Gibbs functions for the species
* at their standard states of solution at the current T and P
* of the solution.
*/
virtual void getGibbs_RT(doublereal* grt) const {
err("getGibbs_RT");
}
/**
* Get the nondimensional Gibbs functions for the standard
* state of the species at the current T and P.
*/
virtual void getPureGibbs(doublereal* gpure) const {
err("getPureGibbs");
}
/**
* Returns the vector of nondimensional
* internal Energies of the standard state at the current temperature
* and pressure of the solution for each species.
*/
virtual void getIntEnergy_RT(doublereal *urt) const {
err("getIntEnergy_RT");
}
/**
* Get the nondimensional Heat Capacities at constant
* pressure for the standard state of the species
* at the current T and P.
*/
virtual void getCp_R(doublereal* cpr) const {
err("getCp_R");
}
/**
* Get the molar volumes of each species in their standard
* states at the current
* <I>T</I> and <I>P</I> of the solution.
* units = m^3 / kmol
*/
virtual void getStandardVolumes(doublereal *vol) const {
err("getStandardVolumes");
}
//@}
/// @name Thermodynamic Values for the Species Reference States --------------------
//@{
/**
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
*/
virtual void getEnthalpy_RT_ref(doublereal *hrt) const;
/**
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
*/
virtual void getGibbs_RT_ref(doublereal *grt) const;
/**
* Returns the vector of the
* gibbs function of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* units = J/kmol
*/
virtual void getGibbs_ref(doublereal *g) const;
/**
* Returns the vector of nondimensional
* entropies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
*/
virtual void getEntropy_R_ref(doublereal *er) const;
/**
* Returns the vector of nondimensional
* constant pressure heat capacities of the reference state
* at the current temperature of the solution
* and reference pressure for the species.
*/
virtual void getCp_R_ref(doublereal *cprt) const;
///////////////////////////////////////////////////////
//
// The methods below are not virtual, and should not
// be overloaded.
//
//////////////////////////////////////////////////////
/**
* @name Specific Properties
* @{
*/
/**
* @name Setting the State
*
* These methods set all or part of the thermodynamic
* state.
* @{
*/
//@}
/**
* @name Chemical Equilibrium
* Chemical equilibrium.
* @{
*/
//@}
/**
* Set equation of state parameter values from XML
* entries. This method is called by function importPhase in
* file importCTML.cpp when processing a phase definition in
* an input file. It should be overloaded in subclasses to set
* any parameters that are specific to that particular phase
* model.
*
* @param eosdata An XML_Node object corresponding to
* the "thermo" entry for this phase in the input file.
*/
virtual void setParametersFromXML(const XML_Node& eosdata) {}
//---------------------------------------------------------
/// @name Critical state properties.
/// These methods are only implemented by some subclasses.
//@{
//@}
/// @name Saturation properties.
/// These methods are only implemented by subclasses that
/// implement full liquid-vapor equations of state.
///
//@}
/// The following methods are used in the process of constructing
/// the phase and setting its parameters from a specification in an
/// input file. They are not normally used in application programs.
/// To see how they are used, see files importCTML.cpp and
/// ThermoFactory.cpp.
/**
* @internal Initialize. This method is provided to allow
* subclasses to perform any initialization required after all
* species have been added. For example, it might be used to
* resize internal work arrays that must have an entry for
* each species. The base class implementation does nothing,
* and subclasses that do not require initialization do not
* need to overload this method. When importing a CTML phase
* description, this method is called just prior to returning
* from function importPhase.
*
* @see importCTML.cpp
*/
virtual void initThermo();
protected:
/*
* The last temperature at which the reference thermodynamic
* properties were calculated at.
*/
mutable doublereal m_tlast;
/**
* Vector containing the species reference enthalpies at T = m_tlast
*/
mutable vector_fp m_h0_RT;
/**
* Vector containing the species reference constant pressure
* heat capacities at T = m_tlast
*/
mutable vector_fp m_cp0_R;
/**
* Vector containing the species reference Gibbs functions
* at T = m_tlast
*/
mutable vector_fp m_g0_RT;
/**
* Vector containing the species reference entropies
* at T = m_tlast
*/
mutable vector_fp m_s0_R;
private:
/**
* VPStandardStateTP has its own err routine
*
*/
doublereal err(string msg) const;
/**
* This function gets called for every call to functions in this
* class. It checks to see whether the temperature has changed and
* thus the reference thermodynamics functions for all of the species
* must be recalculated.
* If the temperature has changed, the species thermo manager is called
* to recalculate G, Cp, H, and S at the current temperature.
*/
void _updateRefStateThermo() const;
};
}
#endif