Started coding up general excess gibbs free energy thermo formulations

based on mole fraction expansions.
This commit is contained in:
Harry Moffat 2009-02-23 20:41:03 +00:00
parent fe281b1a0e
commit 07bcac001b
4 changed files with 729 additions and 3 deletions

View file

@ -147,7 +147,6 @@ namespace Cantera {
/// Return a read-only reference to the vector of atomic weights.
const vector_fp& atomicWeights() const;
/// Number of elements.
int nElements() const;

View file

@ -0,0 +1,315 @@
/**
* @file GibbsExcessVPSSTP.cpp
* Definitions for intermediate ThermoPhase object for phases which
* employ excess gibbs free energy formulations
* (see \ref thermoprops
* and class \link Cantera::GibbsExcessVPSSTP GibbsExcessVPSSTP\endlink).
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based upon expressions
* for the excess gibbs free energy expressed as a function of
* the mole fractions
*/
/*
* Copywrite (2009) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
/*
* $Date$
* $Revision$
*/
#include "GibbsExcessVPSSTP.h"
using namespace std;
namespace Cantera {
/*
* Default constructor.
*
*/
GibbsExcessVPSSTP::GibbsExcessVPSSTP() :
VPStandardStateTP()
{
}
/*
* Copy Constructor:
*
* Note this stuff will not work until the underlying phase
* has a working copy constructor
*/
GibbsExcessVPSSTP::GibbsExcessVPSSTP(const GibbsExcessVPSSTP &b) :
VPStandardStateTP()
{
*this = operator=(b);
}
/*
* operator=()
*
* Note this stuff will not work until the underlying phase
* has a working assignment operator
*/
GibbsExcessVPSSTP& GibbsExcessVPSSTP::
operator=(const GibbsExcessVPSSTP &b) {
if (&b != this) {
VPStandardStateTP::operator=(b);
}
return *this;
}
/**
*
* ~GibbsExcessVPSSTP(): (virtual)
*
* Destructor: does nothing:
*
*/
GibbsExcessVPSSTP::~GibbsExcessVPSSTP() {
}
/*
* This routine duplicates the current object and returns
* a pointer to ThermoPhase.
*/
ThermoPhase*
GibbsExcessVPSSTP::duplMyselfAsThermoPhase() const {
GibbsExcessVPSSTP* mtp = new GibbsExcessVPSSTP(*this);
return (ThermoPhase *) mtp;
}
/*
* -------------- Utilities -------------------------------
*/
// Equation of state type flag.
/*
* The ThermoPhase base class returns
* zero. Subclasses should define this to return a unique
* non-zero value. Known constants defined for this purpose are
* listed in mix_defs.h. The GibbsExcessVPSSTP class also returns
* zero, as it is a non-complete class.
*/
int GibbsExcessVPSSTP::eosType() const {
return 0;
}
/*
* ------------ Molar Thermodynamic Properties ----------------------
*/
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
doublereal GibbsExcessVPSSTP::standardConcentration(int k) const {
err("standardConcentration");
return -1.0;
}
doublereal GibbsExcessVPSSTP::logStandardConc(int k) const {
err("logStandardConc");
return -1.0;
}
void GibbsExcessVPSSTP::getActivities(doublereal* ac) const {
getActivityCoefficients(ac);
getMoleFractions(DATA_PTR(moleFractions_));
for (int k = 0; k < m_kk; k++) {
ac[k] *= moleFractions_[k];
}
}
void GibbsExcessVPSSTP::getElectrochemPotentials(doublereal* mu) const {
getChemPotentials(mu);
double ve = Faraday * electricPotential();
for (int k = 0; k < m_kk; k++) {
mu[k] += ve*charge(k);
}
}
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
doublereal GibbsExcessVPSSTP::err(std::string msg) const {
throw CanteraError("GibbsExcessVPSSTP","Base class method "
+msg+" called. Equation of state type: "+int2str(eosType()));
return 0;
}
/*
* Returns the units of the standard and general concentrations
* Note they have the same units, as their divisor is
* defined to be equal to the activity of the kth species
* in the solution, which is unitless.
*
* This routine is used in print out applications where the
* units are needed. Usually, MKS units are assumed throughout
* the program and in the XML input files.
*
* On return uA contains the powers of the units (MKS assumed)
* of the standard concentrations and generalized concentrations
* for the kth species.
*
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
*/
void GibbsExcessVPSSTP::getUnitsStandardConc(double *uA, int k, int sizeUA) const {
for (int i = 0; i < sizeUA; i++) {
if (i == 0) uA[0] = 1.0;
if (i == 1) uA[1] = -nDim();
if (i == 2) uA[2] = 0.0;
if (i == 3) uA[3] = 0.0;
if (i == 4) uA[4] = 0.0;
if (i == 5) uA[5] = 0.0;
}
}
/*
* @internal Initialize. This method is provided to allow
* subclasses to perform any initialization required after all
* species have been added. For example, it might be used to
* resize internal work arrays that must have an entry for
* each species. The base class implementation does nothing,
* and subclasses that do not require initialization do not
* need to overload this method. When importing a CTML phase
* description, this method is called just prior to returning
* from function importPhase.
*
* @see importCTML.cpp
*/
void GibbsExcessVPSSTP::initThermo() {
initLengths();
VPStandardStateTP::initThermo();
}
// Initialize lengths of local variables after all species have
// been identified.
void GibbsExcessVPSSTP::initLengths() {
m_kk = nSpecies();
moleFractions_.resize(m_kk);
}
/*
* initThermoXML() (virtual from ThermoPhase)
* Import and initialize a ThermoPhase object
*
* @param phaseNode This object must be the phase node of a
* complete XML tree
* description of the phase, including all of the
* species data. In other words while "phase" must
* point to an XML phase object, it must have
* sibling nodes "speciesData" that describe
* the species in the phase.
* @param id ID of the phase. If nonnull, a check is done
* to see if phaseNode is pointing to the phase
* with the correct id.
*/
void GibbsExcessVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) {
initLengths();
VPStandardStateTP::initThermoXML(phaseNode, id);
}
/**
* Format a summary of the mixture state for output.
*/
std::string GibbsExcessVPSSTP::report(bool show_thermo) const {
char p[800];
string s = "";
try {
if (name() != "") {
sprintf(p, " \n %s:\n", name().c_str());
s += p;
}
sprintf(p, " \n temperature %12.6g K\n", temperature());
s += p;
sprintf(p, " pressure %12.6g Pa\n", pressure());
s += p;
sprintf(p, " density %12.6g kg/m^3\n", density());
s += p;
sprintf(p, " mean mol. weight %12.6g amu\n", meanMolecularWeight());
s += p;
doublereal phi = electricPotential();
sprintf(p, " potential %12.6g V\n", phi);
s += p;
int kk = nSpecies();
array_fp x(kk);
array_fp molal(kk);
array_fp mu(kk);
array_fp muss(kk);
array_fp acMolal(kk);
array_fp actMolal(kk);
getMoleFractions(&x[0]);
getChemPotentials(&mu[0]);
getStandardChemPotentials(&muss[0]);
getActivities(&actMolal[0]);
if (show_thermo) {
sprintf(p, " \n");
s += p;
sprintf(p, " 1 kg 1 kmol\n");
s += p;
sprintf(p, " ----------- ------------\n");
s += p;
sprintf(p, " enthalpy %12.6g %12.4g J\n",
enthalpy_mass(), enthalpy_mole());
s += p;
sprintf(p, " internal energy %12.6g %12.4g J\n",
intEnergy_mass(), intEnergy_mole());
s += p;
sprintf(p, " entropy %12.6g %12.4g J/K\n",
entropy_mass(), entropy_mole());
s += p;
sprintf(p, " Gibbs function %12.6g %12.4g J\n",
gibbs_mass(), gibbs_mole());
s += p;
sprintf(p, " heat capacity c_p %12.6g %12.4g J/K\n",
cp_mass(), cp_mole());
s += p;
try {
sprintf(p, " heat capacity c_v %12.6g %12.4g J/K\n",
cv_mass(), cv_mole());
s += p;
}
catch(CanteraError) {
sprintf(p, " heat capacity c_v <not implemented> \n");
s += p;
}
}
} catch (CanteraError) {
;
}
return s;
}
}

View file

@ -0,0 +1,410 @@
/**
* @file gibbsExcessVPSSTP.h
* Header for intermediate ThermoPhase object for phases which
* employ gibbs excess free energy based formulations
* (see \ref thermoprops
* and class \link Cantera::gibbsExcessVPSSTP gibbsExcessVPSSTP\endlink).
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based upon activities
* based on the molality scale. These include most of the methods for
* calculating liquid electrolyte thermodynamics.
*/
/*
* Copywrite (2006) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
/*
* $Id$
*/
#ifndef CT_GIBBSEXCESSVPSSTP_H
#define CT_GIBBSEXCESSVPSSTP_H
#include "VPStandardStateTP.h"
namespace Cantera {
/**
* @ingroup thermoprops
*/
/*!
* GibbsExcessVPSSTP is a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based on
* expressing the Excess Gibbs free energy as a function of
* the mole fractions (or pseudo mole fractions) of consitituents.
* This category is the workhorse for describing molten salts,
* solid-phase mixtures of semiconductors, and mixtures of miscible
* and semi-miscible compounds.
*
* It includes
* . regular solutions
* . Margueles expansions
* . NTRL equation
* . Wilson's equation
* . UNIQUAC equation of state.
*
* This class adds additional functions onto the %ThermoPhase interface
* that handles the calculation of the excess Gibbs free energy. The %ThermoPhase
* class includes a member function, ThermoPhase::activityConvention()
* that indicates which convention the activities are based on. The
* default is to assume activities are based on the molar convention.
* That default is used here.
*
* All of the Excess Gibbs free energy formulations in this area employ
* symmetrical formulations.
*
*
* Chemical potentials
* of species k, \f$ \mu_o \f$, has the following general format:
*
* \f[
* \mu_k = \mu^o_k(T,P) + R T ln( \gamma_k X_k )
* \f]
*
*
* where \f$ \gamma_k^{\triangle} \f$ is a molar based activity coefficient for species
* \f$k\f$.
*
* GibbsExcessVPSSTP contains an internal vector with the current mole
* fraction vector. That's one of its primary usages.
*
* <H3> SetState Strategy </H3>
*
* The gibbsExcessVPSSTP object does not have a setState strategy.
* It's strictly an interfacial layer that writes the current mole fractions to the
* State object.
*
*
*/
class GibbsExcessVPSSTP : public VPStandardStateTP {
public:
/// Constructors
/*!
* This doesn't do much more than initialize constants with
* default values for water at 25C. Water molecular weight
* comes from the default elements.xml file. It actually
* differs slightly from the IAPWS95 value of 18.015268. However,
* density conservation and therefore element conservation
* is the more important principle to follow.
*/
GibbsExcessVPSSTP();
//! Copy constructor
/*!
* Note this stuff will not work until the underlying phase
* has a working copy constructor
*
* @param b class to be copied
*/
GibbsExcessVPSSTP(const GibbsExcessVPSSTP &b);
/// Assignment operator
/*!
* Note this stuff will not work until the underlying phase
* has a working assignment operator
*
* @param b class to be copied.
*/
GibbsExcessVPSSTP& operator=(const GibbsExcessVPSSTP &b);
/// Destructor.
virtual ~GibbsExcessVPSSTP();
//! Duplication routine for objects which inherit from ThermoPhase.
/*!
* This virtual routine can be used to duplicate thermophase objects
* inherited from ThermoPhase even if the application only has
* a pointer to ThermoPhase to work with.
*/
virtual ThermoPhase *duplMyselfAsThermoPhase() const;
/**
*
* @name Utilities
* @{
*/
//! Equation of state type flag.
/*!
* The ThermoPhase base class returns
* zero. Subclasses should define this to return a unique
* non-zero value. Known constants defined for this purpose are
* listed in mix_defs.h. The MolalityVPSSTP class also returns
* zero, as it is a non-complete class.
*/
virtual int eosType() const;
/**
* @}
* @name Molar Thermodynamic Properties
* @{
*/
/**
* @}
* @name Utilities for Solvent ID and Molality
* @{
*/
/**
* @}
* @name Mechanical Properties
* @{
*/
/**
* @}
* @name Potential Energy
*
* Species may have an additional potential energy due to the
* presence of external gravitation or electric fields. These
* methods allow specifying a potential energy for individual
* species.
* @{
*/
/**
* @}
* @name Activities, Standard States, and Activity Concentrations
*
* The activity \f$a_k\f$ of a species in solution is
* related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
* + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
* the chemical potential at unit activity, which depends only
* on temperature and pressure.
* @{
*/
/**
* The standard concentration \f$ C^0_k \f$ used to normalize
* the generalized concentration. In many cases, this quantity
* will be the same for all species in a phase - for example,
* for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this
* reason, this method returns a single value, instead of an
* array. However, for phases in which the standard
* concentration is species-specific (e.g. surface species of
* different sizes), this method may be called with an
* optional parameter indicating the species.
*
* @param k species index. Defaults to zero.
*/
virtual doublereal standardConcentration(int k=0) const;
/**
* Returns the natural logarithm of the standard
* concentration of the kth species
*
* @param k species index
*/
virtual doublereal logStandardConc(int k=0) const;
/**
* Returns the units of the standard and generalized
* concentrations Note they have the same units, as their
* ratio is defined to be equal to the activity of the kth
* species in the solution, which is unitless.
*
* This routine is used in print out applications where the
* units are needed. Usually, MKS units are assumed throughout
* the program and in the XML input files.
*
* @param uA Output vector containing the units
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
* @param k species index. Defaults to 0.
* @param sizeUA output int containing the size of the vector.
* Currently, this is equal to 6.
*/
virtual void getUnitsStandardConc(double *uA, int k = 0,
int sizeUA = 6) const;
//! Get the array of non-dimensional activities (molality
//! based for this class and classes that derive from it) at
//! the current solution temperature, pressure, and solution concentration.
/*!
* All standard state properties for molality-based phases are
* evaluated consistent with the molality scale. Therefore, this function
* must return molality-based activities.
*
* \f[
* a_i^\triangle = \gamma_k^{\triangle} \frac{m_k}{m^\triangle}
* \f]
*
* This function must be implemented in derived classes.
*
* @param ac Output vector of molality-based activities. Length: m_kk.
*/
virtual void getActivities(doublereal* ac) const;
//@}
/// @name Partial Molar Properties of the Solution
//@{
/**
* Get the species electrochemical potentials.
* These are partial molar quantities.
* This method adds a term \f$ Fz_k \phi_k \f$ to the
* to each chemical potential.
*
* Units: J/kmol
*
* @param mu output vector containing the species electrochemical potentials.
* Length: m_kk.
*/
void getElectrochemPotentials(doublereal* mu) const;
//@}
/// @name Properties of the Standard State of the Species in the Solution
//@{
//@}
/// @name Thermodynamic Values for the Species Reference States
//@{
///////////////////////////////////////////////////////
//
// The methods below are not virtual, and should not
// be overloaded.
//
//////////////////////////////////////////////////////
/**
* @name Specific Properties
* @{
*/
/**
* @name Setting the State
*
* These methods set all or part of the thermodynamic
* state.
* @{
*/
//@}
/**
* @name Chemical Equilibrium
* Routines that implement the Chemical equilibrium capability
* for a single phase, based on the element-potential method.
* @{
*/
//@}
/// The following methods are used in the process of constructing
/// the phase and setting its parameters from a specification in an
/// input file. They are not normally used in application programs.
/// To see how they are used, see files importCTML.cpp and
/// ThermoFactory.cpp.
/*!
* @internal Initialize. This method is provided to allow
* subclasses to perform any initialization required after all
* species have been added. For example, it might be used to
* resize internal work arrays that must have an entry for
* each species. The base class implementation does nothing,
* and subclasses that do not require initialization do not
* need to overload this method. When importing a CTML phase
* description, this method is called just prior to returning
* from function importPhase.
*
* @see importCTML.cpp
*/
virtual void initThermo();
/**
* Import and initialize a ThermoPhase object
*
* @param phaseNode This object must be the phase node of a
* complete XML tree
* description of the phase, including all of the
* species data. In other words while "phase" must
* point to an XML phase object, it must have
* sibling nodes "speciesData" that describe
* the species in the phase.
* @param id ID of the phase. If nonnull, a check is done
* to see if phaseNode is pointing to the phase
* with the correct id.
*/
void initThermoXML(XML_Node& phaseNode, std::string id);
//! returns a summary of the state of the phase as a string
/*!
* @param show_thermo If true, extra information is printed out
* about the thermodynamic state of the system.
*/
virtual std::string report(bool show_thermo = true) const;
private:
//! Initialize lengths of local variables after all species have
//! been identified.
void initLengths();
private:
//! Error function
/*!
* Print an error string and exit
*
* @param msg Message to be printed
*/
doublereal err(std::string msg) const;
private:
//! Storage for the current values of the mole fractions of the species
mutable std::vector<doublereal> moleFractions_;
};
}
#endif

View file

@ -78,8 +78,10 @@ ELECTRO_H = MolalityVPSSTP.h VPStandardStateTP.h \
VPSSMgr_Water_ConstVol.h
endif
ifeq ($(do_issp),1)
ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o
ISSP_H = IdealSolidSolnPhase.h StoichSubstanceSSTP.h SingleSpeciesTP.h MineralEQ3.h
ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o \
GibbsExcessVPSSTP.o
ISSP_H = IdealSolidSolnPhase.h StoichSubstanceSSTP.h SingleSpeciesTP.h MineralEQ3.h \
GibbsExcessVPSSTP.h
endif
CATHERMO_OBJ = $(THERMO_OBJ) $(ELECTRO_OBJ) $(ISSP_OBJ)