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

Added the RedlichKisterVPSSTP thermo model to the system.

Browse source 
Fixed an error in MultiPhase that Carlos found.
This commit is contained in:
Harry Moffat 2011-10-31 22:05:07 +00:00
parent 2560431f07
commit 5d1e55596d
8 changed files with 2142 additions and 10 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
Cantera/src/equil/MultiPhase.cpp
View file

@ -293,6 +293,9 @@ namespace Cantera {
}
//====================================================================================================================
int MultiPhase::speciesIndex(std::string speciesName, std::string phaseName) {
if (!m_init) {
init();
}
int p = phaseIndex(phaseName);
if (p < 0) {
throw CanteraError("MultiPhase::speciesIndex", "phase not found: " + phaseName);
@ -382,7 +385,7 @@ namespace Cantera {
else
return false;
}
//====================================================================================================================
/// The Gibbs free energy of the mixture (J).
doublereal MultiPhase::gibbs() const {
index_t i;
@ -455,6 +458,9 @@ namespace Cantera {
/// Set the mole fractions of phase \a n to the values in
/// array \a x.
void MultiPhase::setPhaseMoleFractions(const index_t n, const doublereal* const x) {
if (!m_init) {
init();
}
phase_t* p = m_phase[n];
p->setState_TPX(m_temp, m_press, x);
int nsp = p->nSpecies();

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

@ -81,11 +81,11 @@ ifeq ($(do_issp),1)
ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o \
GibbsExcessVPSSTP.o MolarityIonicVPSSTP.o MargulesVPSSTP.o \
IonsFromNeutralVPSSTP.o PDSS_IonsFromNeutral.o FixedChemPotSSTP.o \
MixedSolventElectrolyte.o
MixedSolventElectrolyte.o RedlichKisterVPSSTP.o
ISSP_H = IdealSolidSolnPhase.h StoichSubstanceSSTP.h SingleSpeciesTP.h MineralEQ3.h \
GibbsExcessVPSSTP.h MolarityIonicVPSSTP.h MargulesVPSSTP.h \
IonsFromNeutralVPSSTP.h PDSS_IonsFromNeutral.h FixedChemPotSSTP.h \
MixedSolventElectrolyte.h
MixedSolventElectrolyte.h RedlichKisterVPSSTP.h
endif
CATHERMO_OBJ = $(THERMO_OBJ) $(ELECTRO_OBJ) $(ISSP_OBJ)

1136
Cantera/src/thermo/RedlichKisterVPSSTP.cpp Normal file
View file

File diff suppressed because it is too large Load diff

983
Cantera/src/thermo/RedlichKisterVPSSTP.h Normal file
View file

@ -0,0 +1,983 @@
/**
* @file RedlichKisterVPSSTP.h
* Header for intermediate ThermoPhase object for phases which
* employ gibbs excess free energy based formulations
* (see \ref thermoprops
* and class \link Cantera::RedlichKisterVPSSTP RedlichKisterVPSSTP\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: MargulesVPSSTP.h 782 2011-10-19 20:45:06Z hkmoffa $
*/
#ifndef CT_REDLICHKISTERVPSSTP_H
#define CT_REDLICHKISTERVPSSTP_H
#include "GibbsExcessVPSSTP.h"
#include "Array.h"
namespace Cantera {
/**
* @ingroup thermoprops
*/
//! RedlichKisterVPSSTP is a derived class of GibbsExcessVPSSTP that employs
//! the Redlich-Kister approximation for the excess gibbs free energy
/*!
*
* %RedlichKisterVPSSTP derives from class GibbsExcessVPSSTP which is derived
* from VPStandardStateTP, and overloads the virtual methods defined there with ones that
* use expressions appropriate for the Redlich Kister 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.
*
*
* <HR>
* <H2> Specification of Species Standard %State Properties </H2>
* <HR>
*
* All species are defined to have standard states that depend upon both
* the temperature and the pressure. The Redlich-Kister 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.
*
*
* <HR>
* <H2> Specification of Solution Thermodynamic Properties </H2>
* <HR>
*
* The molar excess Gibbs free energy is given by the following formula which is a sum over interactions <I>i</I>.
* Each of the interactions are binary interactions involving two of the species in the phase, denoted, <I>Ai</I>
* and <I>Bi</I>.
* This is the generalization of the Redlich-Kister formulation for a phase that has more than 2 species.
*
* \f[
* G^E = \sum_{i} G^E_{i}
* \f]
*
* where
*
* \f[
* G^E_{i} = n X_{Ai} X_{Bi} \sum_m \left( A^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
* \f]
*
* and where we can break down the gibbs free energy contributions into enthalpy and entropy contributions
*
* \f[
* H^E_i = n X_{Ai} X_{Bi} \sum_m \left( H^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
* \f]
*
* \f[
* S^E_i = n X_{Ai} X_{Bi} \sum_m \left( S^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
* \f]
*
* where n is the total moles in the solution.
*
* 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[
* R T \ln( \gamma_k )= \frac{d(n G^E)}{d(n_k)}\Bigg|_{n_i}
* \f]
*
* Taking the derivatives results in the following expression
* \f[
* R T \ln( \gamma_k )= \sum_i \delta_{Ai,k} (1 - X_{Ai}) X_{Bi} \sum_m \left( A^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
* + \sum_i \delta_{Ai,k} X_{Ai} X_{Bi} \sum_m \left( A^{i}_0 + A^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^{m-1} (1 - X_{Ai} + X_{Bi}) \right)
* \f]
* where
*
*
* This object inherits from the class VPStandardStateTP. Therefore, the specification and
* calculation of all standard state and reference state values are handled at that level. Various functional
* forms for the standard state are permissible.
* The chemical potential for species <I>k</I> is equal to
*
* \f[
* \mu_k(T,P) = \mu^o_k(T, P) + R T \ln(\gamma_k X_k)
* \f]
*
* The partial molar entropy for species <I>k</I> is given by the following relation,
*
* \f[
* \tilde{s}_k(T,P) = s^o_k(T,P) - R \ln( \gamma_k X_k )
* - R T \frac{d \ln(\gamma_k) }{dT}
* \f]
*
* The partial molar enthalpy for species <I>k</I> is given by
*
* \f[
* \tilde{h}_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
* \f]
*
* The partial molar volume for species <I>k</I> is
*
* \f[
* \tilde V_k(T,P) = V^o_k(T,P) + R T \frac{d \ln(\gamma_k) }{dP}
* \f]
*
* The partial molar Heat Capacity for species <I>k</I> is
*
* \f[
* \tilde{C}_{p,k}(T,P) = C^o_{p,k}(T,P) - 2 R T \frac{d \ln( \gamma_k )}{dT}
* - R T^2 \frac{d^2 \ln(\gamma_k) }{{dT}^2}
* \f]
*
* <HR>
* <H2> %Application within %Kinetics Managers </H2>
* <HR>
*
* \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 <I>k</I> is independent of <I>k</I> 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.
*
*
* <HR>
* <H2> Instantiation of the Class </H2>
* <HR>
*
*
* 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 <IdealGasPhase *>(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
*
* <HR>
* <H2> XML Example </H2>
* <HR>
* An example of an XML Element named phase setting up a IdealGasPhase
* object named silane is given below.
*
*
* @verbatim
<!-- phase silane -->
<phase dim="3" id="silane">
<elementArray datasrc="elements.xml"> Si H He </elementArray>
<speciesArray datasrc="#species_data">
H2 H HE SIH4 SI SIH SIH2 SIH3 H3SISIH SI2H6
H2SISIH2 SI3H8 SI2 SI3
</speciesArray>
<reactionArray datasrc="#reaction_data"/>
<thermo model="IdealGas"/>
<kinetics model="GasKinetics"/>
<transport model="None"/>
</phase>
@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 RedlichKisterVPSSTP : public GibbsExcessVPSSTP {
public:
//! Constructor
/*!
* This doesn't do much more than initialize constants with
* default values.
*/
RedlichKisterVPSSTP();
//! Construct and initialize a RedlichKisterVPSSTP 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.
*/
RedlichKisterVPSSTP(std::string inputFile, std::string id = "");
//! Construct and initialize a RedlichKisterVPSSTP 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)
*/
RedlichKisterVPSSTP(XML_Node& phaseRef, std::string id = "");
//! Special constructor for a hard-coded problem
/*!
*
* @param testProb Hard-coded value. Only the value of 1 is
* used. It's for
* a LiKCl system
* -> test to predict the eutectic and liquidus correctly.
*/
RedlichKisterVPSSTP(int testProb);
//! Copy constructor
/*!
* Note this stuff will not work until the underlying phase
* has a working copy constructor
*
* @param b class to be copied
*/
RedlichKisterVPSSTP(const RedlichKisterVPSSTP& b);
//! Assignment operator
/*!
*
* @param b class to be copied.
*/
RedlichKisterVPSSTP& operator=(const RedlichKisterVPSSTP &b);
//! Destructor
virtual ~RedlichKisterVPSSTP();
//! Duplication routine for objects which inherit from ThermoPhase.
/*!
* This virtual routine can be used to duplicate thermophase objects
* inherited from ThermoPhase even if the application only has
* a pointer to ThermoPhase to work with.
*/
virtual ThermoPhase *duplMyselfAsThermoPhase() const;
/**
*
* @name Utilities
* @{
*/
//! Equation of state type flag.
/*!
* The ThermoPhase base class returns
* zero. Subclasses should define this to return a unique
* non-zero value. Known constants defined for this purpose are
* listed in mix_defs.h. The MolalityVPSSTP class also returns
* zero, as it is a non-complete class.
*/
virtual int eosType() const;
//! Initialization of a phase using an xml file
/*!
* This routine is a precursor to
* routine, which does most of the work.
*
* @param inputFile XML file containing the description of the
* phase
*
* @param id Optional parameter identifying the name of the
* phase. If none is given, the first XML
* phase element will be used.
*/
void constructPhaseFile(std::string inputFile, std::string id);
//! Import and initialize a phase
//! specification in an XML tree into the current object.
/*!
* Here we read an XML description of the phase.
* We import descriptions of the elements that make up the
* species in a phase.
* We import information about the species, including their
* reference state thermodynamic polynomials. We then freeze
* the state of the species.
*
* Then, we read the species molar volumes from the xml
* tree to finish the initialization.
*
* @param phaseNode This object must be the phase node of a
* complete XML tree
* description of the phase, including all of the
* species data. In other words while "phase" must
* point to an XML phase object, it must have
* sibling nodes "speciesData" that describe
* the species in the phase.
*
* @param id ID of the phase. If nonnull, a check is done
* to see if phaseNode is pointing to the phase
* with the correct id.
*/
void constructPhaseXML(XML_Node& phaseNode, std::string id);
/**
* @}
* @name Molar Thermodynamic Properties
* @{
*/
/**
* @}
* @name Utilities for Solvent ID and Molality
* @{
*/
/**
* @}
* @name Mechanical Properties
* @{
*/
/**
* @}
* @name Potential Energy
*
* Species may have an additional potential energy due to the
* presence of external gravitation or electric fields. These
* methods allow specifying a potential energy for individual
* species.
* @{
*/
/**
* @}
* @name Activities, Standard States, and Activity Concentrations
*
* The activity \f$a_k\f$ of a species in solution is
* related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
* + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
* the chemical potential at unit activity, which depends only
* on temperature and pressure.
* @{
*/
//! 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;
/**
* The standard concentration \f$ C^0_k \f$ used to normalize
* the generalized concentration. In many cases, this quantity
* will be the same for all species in a phase - for example,
* for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this
* reason, this method returns a single value, instead of an
* array. However, for phases in which the standard
* concentration is species-specific (e.g. surface species of
* different sizes), this method may be called with an
* optional parameter indicating the species.
*
* @param k species index. Defaults to zero.
*/
virtual doublereal standardConcentration(int k=0) const;
/**
* Returns the natural logarithm of the standard
* concentration of the kth species
*
* @param k species index
*/
virtual doublereal logStandardConc(int k=0) const;
//! Get the array of non-dimensional molar-based activity coefficients at
//! the current solution temperature, pressure, and solution concentration.
/*!
* @param ac Output vector of activity coefficients. Length: m_kk.
*/
virtual void getActivityCoefficients(doublereal* ac) const;
//@}
/// @name Partial Molar Properties of the Solution
//@{
//! Get the species chemical potentials. Units: J/kmol.
/*!
* This function returns a vector of chemical potentials of the
* species in solution at the current temperature, pressure
* and mole fraction of the solution.
*
* @param mu Output vector of species chemical
* potentials. Length: m_kk. Units: J/kmol
*/
virtual void getChemPotentials(doublereal* mu) const;
/// Molar enthalpy. Units: J/kmol.
virtual doublereal enthalpy_mole() const;
/// Molar entropy. Units: J/kmol.
virtual doublereal entropy_mole() const;
/// Molar heat capacity at constant pressure. Units: J/kmol/K.
virtual doublereal cp_mole() const;
/// Molar heat capacity at constant volume. Units: J/kmol/K.
virtual doublereal cv_mole() const;
//! Returns an array of partial molar enthalpies for the species
//! in the mixture.
/*!
* Units (J/kmol)
*
* For this phase, the partial molar enthalpies are equal to the
* standard state enthalpies modified by the derivative of the
* molality-based activity coefficent wrt temperature
*
* \f[
* \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
* \f]
*
* @param hbar Vector of returned partial molar enthalpies
* (length m_kk, units = J/kmol)
*/
virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
//! Returns an array of partial molar entropies for the species
//! in the mixture.
/*!
* Units (J/kmol)
*
* For this phase, the partial molar enthalpies are equal to the
* standard state enthalpies modified by the derivative of the
* activity coefficent wrt temperature
*
* \f[
* \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
* - R \ln( \gamma_k X_k)
* - R T \frac{d \ln(\gamma_k) }{dT}
* \f]
*
* @param sbar Vector of returned partial molar entropies
* (length m_kk, units = J/kmol/K)
*/
virtual void getPartialMolarEntropies(doublereal* sbar) const;
//! Returns an array of partial molar entropies for the species
//! in the mixture.
/*!
* Units (J/kmol)
*
* For this phase, the partial molar enthalpies are equal to the
* standard state enthalpies modified by the derivative of the
* activity coefficent wrt temperature
*
* \f[
* ???????????????
* \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
* - R \ln( \gamma_k X_k)
* - R T \frac{d \ln(\gamma_k) }{dT}
* ???????????????
* \f]
*
* @param cpbar Vector of returned partial molar heat capacities
* (length m_kk, units = J/kmol/K)
*/
virtual void getPartialMolarCp(doublereal* cpbar) const;
//! Return an array of partial molar volumes for the
//! species in the mixture. Units: m^3/kmol.
/*!
* Frequently, for this class of thermodynamics representations,
* the excess Volume due to mixing is zero. Here, we set it as
* a default. It may be overriden in derived classes.
*
* @param vbar Output vector of speciar partial molar volumes.
* Length = m_kk. units are m^3/kmol.
*/
virtual void getPartialMolarVolumes(doublereal* vbar) const;
//! 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., units = J/kmol
*/
void getElectrochemPotentials(doublereal* mu) const;
//! Get the array of temperature second derivatives of the log activity coefficients
/*!
* This function is a virtual class, but it first appears in GibbsExcessVPSSTP
* class and derived classes from GibbsExcessVPSSTP.
*
* units = 1/Kelvin
*
* @param d2lnActCoeffdT2 Output vector of temperature 2nd derivatives of the
* log Activity Coefficients. length = m_kk
*
*/
virtual void getd2lnActCoeffdT2(doublereal *d2lnActCoeffdT2) const;
//! Get the array of temperature derivatives of the log activity coefficients
/*!
* This function is a virtual class, but it first appears in GibbsExcessVPSSTP
* class and derived classes from GibbsExcessVPSSTP.
*
* units = 1/Kelvin
*
* @param dlnActCoeffdT Output vector of temperature derivatives of the
* log Activity Coefficients. length = m_kk
*
*/
virtual void getdlnActCoeffdT(doublereal *dlnActCoeffdT) 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);
/**
* @}
* @name Derivatives of Thermodynamic Variables needed for Applications
* @{
*/
//! Get the change in activity coefficients w.r.t. change in state (temp, mole fraction, etc.) along
//! a line in parameter space or along a line in physical space
/*!
*
* @param dTds Input of temperature change along the path
* @param dXds Input vector of changes in mole fraction along the path. length = m_kk
* Along the path length it must be the case that the mole fractions sum to one.
* @param dlnActCoeffds Output vector of the directional derivatives of the
* log Activity Coefficients along the path. length = m_kk
* units are 1/units(s). if s is a physical coordinate then the units are 1/m.
*/
virtual void getdlnActCoeffds(const doublereal dTds, const doublereal * const dXds, doublereal *dlnActCoeffds) const;
//! Get the array of log concentration-like derivatives of the
//! log activity coefficients - diagonal component
/*!
* This function is a virtual method. For ideal mixtures
* (unity activity coefficients), this can return zero.
* Implementations should take the derivative of the
* logarithm of the activity coefficient with respect to the
* logarithm of the mole fraction.
*
* units = dimensionless
*
* @param dlnActCoeffdlnX_diag Output vector of the diagonal component of the log(mole fraction)
* derivatives of the log Activity Coefficients.
* length = m_kk
*/
virtual void getdlnActCoeffdlnX_diag(doublereal *dlnActCoeffdlnX_diag) const;
//! Get the array of derivatives of the log activity coefficients wrt mole numbers - diagonal only
/*!
* This function is a virtual method. For ideal mixtures
* (unity activity coefficients), this can return zero.
* Implementations should take the derivative of the
* logarithm of the activity coefficient with respect to the
* logarithm of the concentration-like variable (i.e. mole fraction,
* molality, etc.) that represents the standard state.
*
* units = dimensionless
*
* @param dlnActCoeffdlnN_diag Output vector of the diagonal entries for the log(mole fraction)
* derivatives of the log Activity Coefficients.
* length = m_kk
*/
virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const;
//! Get the array of derivatives of the ln activity coefficients with respect to the ln species mole numbers
/*!
* Implementations should take the derivative of the logarithm of the activity coefficient with respect to a
* log of a species mole number (with all other species mole numbers held constant)
*
* units = 1 / kmol
*
* dlnActCoeffdlnN[ ld * k + m] will contain the derivative of log act_coeff for the <I>m</I><SUP>th</SUP>
* species with respect to the number of moles of the <I>k</I><SUP>th</SUP> species.
*
* \f[
* \frac{d \ln(\gamma_m) }{d \ln( n_k ) }\Bigg|_{n_i}
* \f]
*
* @param ld Number of rows in the matrix
* @param dlnActCoeffdlnN Output vector of derivatives of the
* log Activity Coefficients. length = m_kk * m_kk
*/
virtual void getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN) ;
//@}
private:
//! Process an XML node called "binaryNeutralSpeciesParameters"
/*!
* This node contains all of the parameters necessary to describe
* the Redlich-Kister model for a particular binary interaction.
* This function reads the XML file and writes the coefficients
* it finds to an internal data structures.
*
* @param xmlBinarySpecies Reference to the XML_Node named "binaryNeutralSpeciesParameters"
* containing the binary interaction
*/
void readXMLBinarySpecies(XML_Node &xmlBinarySpecies);
//! Resize internal arrays within the object that depend upon the number
//! of binary Redlich-Kister interaction terms
/*!
* @param num Number of binary Redlich-Kister interaction terms
*/
void resizeNumInteractions(const int num);
//! 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;
//! Update the derivative of the log of the activity coefficients wrt T
/*!
* This function will be called to update the internally storred
* derivative of the natural logarithm of the activity coefficients
* wrt temperature.
*/
void s_update_dlnActCoeff_dT() const;
//! Internal routine that calculates the derivative of the activity coefficients wrt
//! the mole fractions.
/*!
* This routine calculates the the derivative of the activity coefficients wrt to mole fraction
* with all other mole fractions held constant. This is strictly not permitted. However, if the
* resulting matrix is multiplied by a permissible deltaX vector then everything is ok.
*
* This is the natural way to handle concentration derivatives in this routine.
*/
void s_update_dlnActCoeff_dX_() const;
#ifdef DEBUG_MODE
public:
//! Utility routine that calculates a literature expression
/*!
* @param VintOut Output contribution to the voltage corresponding to nonideal term
* @param voltsOut Output contribution to the voltage corresponding to nonideal term and mf term
*/
void Vint(double &VintOut, double &voltsOut) ;
#endif
private:
//! Error function
/*!
* Print an error string and exit
*
* @param msg Message to be printed
*/
doublereal err(std::string msg) const;
protected:
//! number of binary interaction expressions
int numBinaryInteractions_;
//! vector of species indices representing species A in the interaction
/*!
* Each Redlich-Kister excess Gibbs free energy term involves two species, A and B.
* This vector identifies species A.
*/
vector_int m_pSpecies_A_ij;
//! vector of species indices representing species B in the interaction
/*!
* Each Redlich-Kisterexcess Gibbs free energy term involves two species, A and B.
* This vector identifies species B.
*/
vector_int m_pSpecies_B_ij;
//! Vector of the length of the polynomial for the interaction.
vector_int m_N_ij;
//! Enthalpy term for the binary mole fraction interaction of the
//! excess gibbs free energy expression
mutable std::vector< vector_fp> m_HE_m_ij;
//! Entropy term for the binary mole fraction interaction of the
//! excess gibbs free energy expression
mutable std::vector< vector_fp> m_SE_m_ij;
//! form of the RedlichKister interaction expression
/*!
* Currently there is only one form.
*/
int formRedlichKister_;
//! form of the temperatuer dependence of the Redlich-Kister interaction expression
/*!
* Currently there is only one form -> constant wrt temperature.
*/
int formTempModel_;
//! Two dimensional array of derivatives of activity coefficients wrt mole fractions
mutable Array2D dlnActCoeff_dX_;
};
}
#endif

12
Cantera/src/thermo/ThermoFactory.cpp
View file

@ -28,6 +28,7 @@
#ifdef WITH_IDEAL_SOLUTIONS
#include "IdealSolidSolnPhase.h"
#include "MargulesVPSSTP.h"
#include "RedlichKisterVPSSTP.h"
#include "IonsFromNeutralVPSSTP.h"
#include "PhaseCombo_Interaction.h"
#endif
@ -99,7 +100,7 @@ namespace Cantera {
/*!
* @deprecated This entire structure could be replaced with a std::map
*/
static int ntypes = 22;
static int ntypes = 23;
//! Define the string name of the %ThermoPhase types that are handled by this factory routine
static string _types[] = {"IdealGas", "Incompressible",
@ -109,7 +110,7 @@ namespace Cantera {
"IdealMolalSolution", "IdealGasVPSS",
"MineralEQ3", "MetalSHEelectrons", "Margules", "PhaseCombo_Interaction",
"IonsFromNeutralMolecule", "FixedChemPot", "MolarityIonicVPSSTP",
"MixedSolventElectrolyte"
"MixedSolventElectrolyte", "Redlich-Kister"
};
//! Define the integer id of the %ThermoPhase types that are handled by this factory routine
@ -120,7 +121,7 @@ namespace Cantera {
cIdealMolalSoln, cVPSS_IdealGas,
cMineralEQ3, cMetalSHEelectrons,
cMargulesVPSSTP, cPhaseCombo_Interaction, cIonsFromNeutral, cFixedChemPot,
cMolarityIonicVPSSTP, cMixedSolventElectrolyte
cMolarityIonicVPSSTP, cMixedSolventElectrolyte, cRedlichKisterVPSSTP
};
/*
@ -162,6 +163,10 @@ namespace Cantera {
th = new MargulesVPSSTP();
break;
case cRedlichKisterVPSSTP:
th = new RedlichKisterVPSSTP();
break;
case cPhaseCombo_Interaction:
th = new PhaseCombo_Interaction();
break;
@ -811,7 +816,6 @@ namespace Cantera {
factory->installThermoForSpecies(k, s, &th, *spthermo_ptr, phaseNode_ptr);
}
return true;
}

5
Cantera/src/thermo/ThermoFactory.h
View file

@ -255,8 +255,7 @@ namespace Cantera {
*
* @ingroup thermoprops
*/
bool importPhase(XML_Node& phase, ThermoPhase* th,
SpeciesThermoFactory* spfactory = 0);
bool importPhase(XML_Node& phase, ThermoPhase* th, SpeciesThermoFactory* spfactory = 0);
//! Install a species into a ThermoPhase object, which defines
//! the phase thermodynamics and speciation.
@ -305,7 +304,7 @@ namespace Cantera {
VPSSMgr *vpss_ptr = 0,
SpeciesThermoFactory* factory = 0);
//!Search an XML tree for species data.
//! Search an XML tree for species data.
/*!
* This utility routine will search the XML tree for the species
* named by the string, kname. It will return the XML_Node

2
Cantera/src/thermo/mix_defs.h
View file

@ -79,6 +79,8 @@ namespace Cantera {
const int cMargulesVPSSTP = 301;
const int cRedlichKisterVPSSTP = 303;
const int cMolarityIonicVPSSTP = 401;
const int cMixedSolventElectrolyte = 402;

2
docs/Cantera.cfg.in
View file

@ -707,6 +707,8 @@ FILE_PATTERNS = Kinetics.h Kinetics.cpp \
GibbsExcessVPSSTP.cpp \
MargulesVPSSTP.h \
MargulesVPSSTP.cpp \
RedlichKisterVPSSTP.h \
RedlichKisterVPSSTP.cpp \
IonsFromNeutralVPSSTP.h \
IonsFromNeutralVPSSTP.cpp \
VPSSMgr.h \