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

Changed m_speciesData to be a vector of pointers. This is how it is

Browse source 
used anyway.
This commit is contained in:
Harry Moffat 2008-12-24 18:18:23 +00:00
parent 12e553f234
commit 4d4f43e679
10 changed files with 1316 additions and 111 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

65
Cantera/src/thermo/DebyeHuckel.cpp
View file

@ -1480,32 +1480,29 @@ namespace Cantera {
* -> Look for the subelement "stoichIsMods"
* in each of the species SS databases.
*/
const XML_Node *phaseSpecies = speciesData();
if (phaseSpecies) {
std::string kname, jname;
vector<XML_Node*> xspecies;
phaseSpecies->getChildren("species",xspecies);
int jj = xspecies.size();
for (k = 0; k < m_kk; k++) {
int jmap = -1;
kname = speciesName(k);
for (int j = 0; j < jj; j++) {
const XML_Node& sp = *xspecies[j];
jname = sp["name"];
if (jname == kname) {
jmap = j;
break;
}
std::vector<const XML_Node *> xspecies= speciesData();
std::string kname, jname;
int jj = xspecies.size();
for (k = 0; k < m_kk; k++) {
int jmap = -1;
kname = speciesName(k);
for (int j = 0; j < jj; j++) {
const XML_Node& sp = *xspecies[j];
jname = sp["name"];
if (jname == kname) {
jmap = j;
break;
}
if (jmap > -1) {
const XML_Node& sp = *xspecies[jmap];
if (sp.hasChild("stoichIsMods")) {
double val = getFloat(sp, "stoichIsMods");
m_speciesCharge_Stoich[k] = val;
}
}
if (jmap > -1) {
const XML_Node& sp = *xspecies[jmap];
if (sp.hasChild("stoichIsMods")) {
double val = getFloat(sp, "stoichIsMods");
m_speciesCharge_Stoich[k] = val;
}
}
}
/*
* Now look at the activity coefficient database
*/
@ -1558,20 +1555,18 @@ namespace Cantera {
* -> Look for the subelement "stoichIsMods"
* in each of the species SS databases.
*/
const XML_Node *phaseSpecies = speciesData();
std::vector<const XML_Node *> xspecies= speciesData();
const XML_Node *spPtr = 0;
if (phaseSpecies) {
std::string kname;
for (k = 0; k < m_kk; k++) {
kname = speciesName(k);
spPtr = speciesXML_Node(kname, phaseSpecies);
if (!spPtr) {
if (spPtr->hasChild("electrolyteSpeciesType")) {
std::string est = getString(*spPtr, "electrolyteSpeciesType");
if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
throw CanteraError("DebyeHuckel:initThermoXML",
"Bad electrolyte type: " + est);
}
std::string kname;
for (k = 0; k < m_kk; k++) {
kname = speciesName(k);
spPtr = xspecies[k];
if (!spPtr) {
if (spPtr->hasChild("electrolyteSpeciesType")) {
std::string est = getString(*spPtr, "electrolyteSpeciesType");
if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
throw CanteraError("DebyeHuckel:initThermoXML",
"Bad electrolyte type: " + est);
}
}
}

9
Cantera/src/thermo/HMWSoln.h
View file

@ -1987,7 +1987,9 @@ namespace Cantera {
*/
void constructPhaseXML(XML_Node& phaseNode, std::string id);
/**
//! Internal initialization required after all species have
//! been added
/*!
* @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
@ -2002,8 +2004,9 @@ namespace Cantera {
*/
virtual void initThermo();
/*
* initThermoXML() (virtual from ThermoPhase)
//! Initialize the phase parameters from an XML file.
/*!
* initThermoXML() (virtual from ThermoPhase)
*
* This gets called from importPhase(). It processes the XML file
* after the species are set up. This is the main routine for

68
Cantera/src/thermo/HMWSoln_input.cpp
View file

@ -1359,32 +1359,30 @@ namespace Cantera {
* -> Look for the subelement "stoichIsMods"
* in each of the species SS databases.
*/
const XML_Node *phaseSpecies = speciesData();
if (phaseSpecies) {
string kname, jname;
vector<XML_Node*> xspecies;
phaseSpecies->getChildren("species", xspecies);
int jj = xspecies.size();
for (k = 0; k < m_kk; k++) {
int jmap = -1;
kname = speciesName(k);
for (int j = 0; j < jj; j++) {
const XML_Node& sp = *xspecies[j];
jname = sp["name"];
if (jname == kname) {
jmap = j;
break;
}
std::vector<const XML_Node *> xspecies = speciesData();
string kname, jname;
int jj = xspecies.size();
for (k = 0; k < m_kk; k++) {
int jmap = -1;
kname = speciesName(k);
for (int j = 0; j < jj; j++) {
const XML_Node& sp = *xspecies[j];
jname = sp["name"];
if (jname == kname) {
jmap = j;
break;
}
if (jmap > -1) {
const XML_Node& sp = *xspecies[jmap];
if (sp.hasChild("stoichIsMods")) {
double val = getFloat(sp, "stoichIsMods");
m_speciesCharge_Stoich[k] = val;
}
}
if (jmap > -1) {
const XML_Node& sp = *xspecies[jmap];
if (sp.hasChild("stoichIsMods")) {
double val = getFloat(sp, "stoichIsMods");
m_speciesCharge_Stoich[k] = val;
}
}
}
/*
* Now look at the activity coefficient database
*/
@ -1470,23 +1468,21 @@ namespace Cantera {
* -> Look for the subelement "stoichIsMods"
* in each of the species SS databases.
*/
const XML_Node *phaseSpecies = speciesData();
std::vector<const XML_Node *> xspecies = speciesData();
const XML_Node *spPtr = 0;
if (phaseSpecies) {
string kname;
for (k = 0; k < m_kk; k++) {
kname = speciesName(k);
spPtr = speciesXML_Node(kname, phaseSpecies);
if (!spPtr) {
if (spPtr->hasChild("electrolyteSpeciesType")) {
string est = getString(*spPtr, "electrolyteSpeciesType");
if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
throw CanteraError("HMWSoln::initThermoXML",
"Bad electrolyte type: " + est);
}
string kname;
for (k = 0; k < m_kk; k++) {
kname = speciesName(k);
spPtr = xspecies[k];
if (!spPtr) {
if (spPtr->hasChild("electrolyteSpeciesType")) {
string est = getString(*spPtr, "electrolyteSpeciesType");
if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
throw CanteraError("HMWSoln::initThermoXML",
"Bad electrolyte type: " + est);
}
}
}
}
}
/*
* Then look at the phase thermo specification

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

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

557
Cantera/src/thermo/MineralEQ3.cpp Normal file
View file

@ -0,0 +1,557 @@
/**
* @file MineralEQ3.cpp
* Definition file for the MineralEQ3 class, which represents a fixed-composition
* incompressible substance (see \ref thermoprops and
* class \link Cantera::MineralEQ3 MineralEQ3\endlink)
*/
/*
* 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.
*
* Copyright 2001 California Institute of Technology
*/
/*
* $Id$
*/
#include "ct_defs.h"
#include "mix_defs.h"
#include "MineralEQ3.h"
#include "SpeciesThermo.h"
#include "ThermoFactory.h"
#include "MineralEQ3.h"
#include <string>
using namespace std;
namespace Cantera {
/*
* ---- Constructors -------
*/
/*
* Default Constructor for the MineralEQ3 class
*/
MineralEQ3::MineralEQ3():
StoichSubstanceSSTP ()
{
}
// Create and initialize a MineralEQ3 ThermoPhase object
// from an asci input file
/*
* @param infile name of the input file
* @param id name of the phase id in the file.
* If this is blank, the first phase in the file is used.
*/
MineralEQ3::MineralEQ3(std::string infile, std::string id) :
StoichSubstanceSSTP()
{
XML_Node* root = get_XML_File(infile);
if (id == "-") id = "";
XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, root);
if (!xphase) {
throw CanteraError("MineralEQ3::MineralEQ3",
"Couldn't find phase name in file:" + id);
}
// Check the model name to ensure we have compatibility
const XML_Node& th = xphase->child("thermo");
std::string model = th["model"];
if (model != "StoichSubstance" && model != "MineralEQ3") {
throw CanteraError("MineralEQ3::MineralEQ3",
"thermo model attribute must be StoichSubstance");
}
importPhase(*xphase, this);
}
// Full Constructor.
/*
* @param phaseRef XML node pointing to a MineralEQ3 description
* @param id Id of the phase.
*/
MineralEQ3::MineralEQ3(XML_Node& xmlphase, std::string id) :
StoichSubstanceSSTP()
{
if (id != "") {
std::string idxml = xmlphase["id"];
if (id != idxml) {
throw CanteraError("MineralEQ3::MineralEQ3",
"id's don't match");
}
}
const XML_Node& th = xmlphase.child("thermo");
std::string model = th["model"];
if (model != "StoichSubstance" && model != "MineralEQ3") {
throw CanteraError("MineralEQ3::MineralEQ3",
"thermo model attribute must be StoichSubstance");
}
importPhase(xmlphase, this);
}
//! Copy constructor
/*!
* @param right Object to be copied
*/
MineralEQ3::MineralEQ3(const MineralEQ3 &right) :
StoichSubstanceSSTP()
{
*this = operator=(right);
}
//! Assignment operator
/*!
* @param right Object to be copied
*/
MineralEQ3 &
MineralEQ3::operator=(const MineralEQ3 & right) {
if (&right == this) {
return *this;
}
StoichSubstanceSSTP::operator=(right);
m_Mu0_tr_pr = right.m_Mu0_tr_pr;
m_Entrop_tr_pr = right.m_Entrop_tr_pr;
m_deltaG_formation_tr_pr = right.m_deltaG_formation_tr_pr;
m_deltaH_formation_tr_pr = right.m_deltaH_formation_tr_pr;
m_V0_tr_pr = right.m_V0_tr_pr;
m_a = right.m_a;
m_b = right.m_b;
m_c = right.m_c;
return *this;
}
/*
* Destructor for the routine (virtual)
*
*/
MineralEQ3::~MineralEQ3()
{
}
// Duplication function
/*
* This virtual function is used to create a duplicate of the
* current phase. It's used to duplicate the phase when given
* a ThermoPhase pointer to the phase.
*
* @return It returns a ThermoPhase pointer.
*/
ThermoPhase *MineralEQ3::duplMyselfAsThermoPhase() const {
MineralEQ3 *stp = new MineralEQ3(*this);
return (ThermoPhase *) stp;
}
/*
* ---- Utilities -----
*/
/*
* Equation of state flag. Returns the value cStoichSubstance,
* defined in mix_defs.h.
*/
int MineralEQ3::eosType() const {
return cStoichSubstance;
}
/*
* ---- Molar Thermodynamic properties of the solution ----
*/
/**
* ----- Mechanical Equation of State ------
*/
/*
* Pressure. Units: Pa.
* For an incompressible substance, the density is independent
* of pressure. This method simply returns the stored
* pressure value.
*/
doublereal MineralEQ3::pressure() const {
return m_press;
}
/*
* Set the pressure at constant temperature. Units: Pa.
* For an incompressible substance, the density is
* independent of pressure. Therefore, this method only
* stores the specified pressure value. It does not
* modify the density.
*/
void MineralEQ3::setPressure(doublereal p) {
m_press = p;
}
/*
* The isothermal compressibility. Units: 1/Pa.
* The isothermal compressibility is defined as
* \f[
* \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
* \f]
*
* It's equal to zero for this model, since the molar volume
* doesn't change with pressure or temperature.
*/
doublereal MineralEQ3::isothermalCompressibility() const {
return 0.0;
}
/*
* The thermal expansion coefficient. Units: 1/K.
* The thermal expansion coefficient is defined as
*
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*
* It's equal to zero for this model, since the molar volume
* doesn't change with pressure or temperature.
*/
doublereal MineralEQ3::thermalExpansionCoeff() const {
return 0.0;
}
/*
* ---- Chemical Potentials and Activities ----
*/
/*
* This method returns the array of generalized
* concentrations. For a stoichiomeetric substance, there is
* only one species, and the generalized concentration is 1.0.
*/
void MineralEQ3::
getActivityConcentrations(doublereal* c) const {
c[0] = 1.0;
}
/*
* The standard concentration. This is defined as the concentration
* by which the generalized concentration is normalized to produce
* the activity.
*/
doublereal MineralEQ3::standardConcentration(int k) const {
return 1.0;
}
/*
* Returns the natural logarithm of the standard
* concentration of the kth species
*/
doublereal MineralEQ3::logStandardConc(int k) const {
return 0.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
*/
void MineralEQ3::
getUnitsStandardConc(doublereal *uA, int k, int sizeUA) const {
for (int i = 0; i < 6; i++) {
uA[i] = 0;
}
}
/*
* ---- Partial Molar Properties of the Solution ----
*/
/*
* ---- Properties of the Standard State of the Species in the Solution
* ----
*/
/*
* Get the array of chemical potentials at unit activity
* \f$ \mu^0_k \f$.
*
* For a stoichiometric substance, there is no activity term in
* the chemical potential expression, and therefore the
* standard chemical potential and the chemical potential
* are both equal to the molar Gibbs function.
*/
void MineralEQ3::
getStandardChemPotentials(doublereal* mu0) const {
getGibbs_RT(mu0);
mu0[0] *= GasConstant * temperature();
}
/*
* 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.
* Molar enthalpy. Units: J/kmol. For an incompressible,
* stoichiometric substance, the internal energy is
* independent of pressure, and therefore the molar enthalpy
* is \f[ \hat h(T, P) = \hat u(T) + P \hat v \f], where the
* molar specific volume is constant.
*/
void MineralEQ3::getEnthalpy_RT(doublereal* hrt) const {
getEnthalpy_RT_ref(hrt);
doublereal RT = GasConstant * temperature();
doublereal presCorrect = (m_press - m_p0) / molarDensity();
hrt[0] += presCorrect / RT;
}
/*
* Get the array of nondimensional Entropy functions for the
* standard state species
* at the current <I>T</I> and <I>P</I> of the solution.
*/
void MineralEQ3::getEntropy_R(doublereal* sr) const {
getEntropy_R_ref(sr);
}
/*
* Get the nondimensional Gibbs functions for the species
* at their standard states of solution at the current T and P
* of the solution
*/
void MineralEQ3::getGibbs_RT(doublereal* grt) const {
getEnthalpy_RT(grt);
grt[0] -= m_s0_R[0];
}
/*
* Get the nondimensional Gibbs functions for the standard
* state of the species at the current T and P.
*/
void MineralEQ3::getCp_R(doublereal* cpr) const {
_updateThermo();
cpr[0] = m_cp0_R[0];
}
/*
* Molar internal energy (J/kmol).
* For an incompressible,
* stoichiometric substance, the molar internal energy is
* independent of pressure. Since the thermodynamic properties
* are specified by giving the standard-state enthalpy, the
* term \f$ P_0 \hat v\f$ is subtracted from the specified molar
* enthalpy to compute the molar internal energy.
*/
void MineralEQ3::getIntEnergy_RT(doublereal* urt) const {
_updateThermo();
doublereal RT = GasConstant * temperature();
doublereal PV = m_p0 / molarDensity();
urt[0] = m_h0_RT[0] - PV / RT;
}
/*
* ---- Thermodynamic Values for the Species Reference States ----
*/
/*
* Molar internal energy or the reference state at the current
* temperature, T (J/kmol).
* For an incompressible,
* stoichiometric substance, the molar internal energy is
* independent of pressure. Since the thermodynamic properties
* are specified by giving the standard-state enthalpy, the
* term \f$ P_0 \hat v\f$ is subtracted from the specified molar
* enthalpy to compute the molar internal energy.
*
* Note, this is equal to the standard state internal energy
* evaluated at the reference pressure.
*/
void MineralEQ3::getIntEnergy_RT_ref(doublereal* urt) const {
_updateThermo();
doublereal RT = GasConstant * temperature();
doublereal PV = m_p0 / molarDensity();
urt[0] = m_h0_RT[0] - PV / RT;
}
/*
* ---- Saturation Properties
*/
/*
* ---- Initialization and Internal functions
*/
/**
* @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 MineralEQ3::initThermo() {
/*
* Call the base class thermo initializer
*/
StoichSubstanceSSTP::initThermo();
}
/**
* setParameters:
*
* Generic routine that is used to set the parameters used
* by this model.
* C[0] = density of phase [ kg/m3 ]
*/
void MineralEQ3::setParameters(int n, doublereal * const c) {
doublereal rho = c[0];
setDensity(rho);
}
/**
* getParameters:
*
* Generic routine that is used to get the parameters used
* by this model.
* n = 1
* C[0] = density of phase [ kg/m3 ]
*/
void MineralEQ3::getParameters(int &n, doublereal * const c) const {
doublereal rho = density();
n = 1;
c[0] = rho;
}
// Initialize the phase parameters from an XML file.
/*
* initThermoXML() (virtual from ThermoPhase)
*
* This gets called from importPhase(). It processes the XML file
* after the species are set up. This is the main routine for
* reading in activity coefficient parameters.
*
* @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 MineralEQ3::initThermoXML(XML_Node& phaseNode, std::string id) {
/*
* Find the Thermo XML node
*/
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("HMWSoln::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
}
/*
* Reads an xml data block for the parameters needed by this
* routine. eosdata is a reference to the xml thermo block, and looks
* like this:
*
* <phase id="stoichsolid" >
* <thermo model="StoichSubstance">
* <density units="g/cm3">3.52</density>
* </thermo>
* </phase>
*/
void MineralEQ3::setParametersFromXML(const XML_Node& eosdata) {
std::string model = eosdata["model"];
if (model != "StoichSubstance" && model != "MineralEQ3") {
throw CanteraError("MineralEQ3::MineralEQ3",
"thermo model attribute must be StoichSubstance");
}
doublereal rho = getFloat(eosdata, "density", "toSI");
setDensity(rho);
}
doublereal MineralEQ3::LookupGe(const std::string& elemName) {
#ifdef OLDWAY
int num = sizeof(geDataTable) / sizeof(struct GeData);
string s3 = elemName.substr(0,3);
for (int i = 0; i < num; i++) {
//if (!std::strncmp(elemName.c_str(), aWTable[i].name, 3)) {
if (s3 == geDataTable[i].name) {
return (geDataTable[i].GeValue);
}
}
throw CanteraError("LookupGe", "element " + s + " not found");
return -1.0;
#else
int iE = elementIndex(elemName);
if (iE < 0) {
throw CanteraError("PDSS_HKFT::LookupGe", "element " + elemName + " not found");
}
doublereal geValue = entropyElement298(iE);
if (geValue == ENTROPY298_UNKNOWN) {
throw CanteraError("PDSS_HKFT::LookupGe",
"element " + elemName + " doesn not have a supplied entropy298");
}
geValue *= (-298.15);
return geValue;
#endif
}
void MineralEQ3::convertDGFormation() {
/*
* Ok let's get the element compositions and conversion factors.
*/
int ne = nElements();
doublereal na;
doublereal ge;
string ename;
doublereal totalSum = 0.0;
for (int m = 0; m < ne; m++) {
na = nAtoms(0, m);
if (na > 0.0) {
ename = elementName(m);
ge = LookupGe(ename);
totalSum += na * ge;
}
}
// Add in the charge
// if (m_charge_j != 0.0) {
// ename = "H";
// ge = LookupGe(ename);
// totalSum -= m_charge_j * ge;
//}
// Ok, now do the calculation. Convert to joules kmol-1
doublereal dg = m_deltaG_formation_tr_pr * 4.184 * 1.0E3;
//! Store the result into an internal variable.
m_Mu0_tr_pr = dg + totalSum;
}
}

645
Cantera/src/thermo/MineralEQ3.h Normal file
View file

@ -0,0 +1,645 @@
/**
* @file StoichSubstanceSSTP.h
* Header file for the StoichSubstanceSSTP class, which represents a fixed-composition
* incompressible substance (see \ref thermoprops and
* class \link Cantera::StoichSubstanceSSTP StoichSubstanceSSTP\endlink)
*/
/*
* 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_MINERALEQ3_H
#define CT_MINERALEQ3_H
#include "mix_defs.h"
#include "SingleSpeciesTP.h"
#include "SpeciesThermo.h"
#include "StoichSubstanceSSTP.h"
namespace Cantera {
//! Class %StoichSubstanceSSTP represents a stoichiometric (fixed
//! composition) incompressible substance.
/*!
* This class internally changes the independent degree of freedom from
* density to pressure. This is necessary because the phase is
* incompressible. It uses a constant volume approximation.
*
*
* <b> Specification of Species Standard %State Properties </b>
*
* This class inherits from SingleSpeciesTP.
* It is assumed that the reference state thermodynamics may be
* obtained by a pointer to a populated species thermodynamic property
* manager class (see ThermoPhase::m_spthermo). How to relate pressure
* changes to the reference state thermodynamics is resolved at this level.
*
* For an incompressible,
* stoichiometric substance, the molar internal energy is
* independent of pressure. Since the thermodynamic properties
* are specified by giving the standard-state enthalpy, the
* term \f$ P_0 \hat v\f$ is subtracted from the specified molar
* enthalpy to compute the molar internal energy. The entropy is
* assumed to be independent of the pressure.
*
* The enthalpy function is given by the following relation.
*
* \f[
* \raggedright h^o_k(T,P) =
* h^{ref}_k(T) + \tilde v \left( P - P_{ref} \right)
* \f]
*
* For an incompressible,
* stoichiometric substance, the molar internal energy is
* independent of pressure. Since the thermodynamic properties
* are specified by giving the standard-state enthalpy, the
* term \f$ P_{ref} \tilde v\f$ is subtracted from the specified reference molar
* enthalpy to compute the molar internal energy.
*
* \f[
* u^o_k(T,P) = h^{ref}_k(T) - P_{ref} \tilde v
* \f]
*
* The standard state heat capacity and entropy are independent
* of pressure. The standard state gibbs free energy is obtained
* from the enthalpy and entropy functions.
*
*
* <b> Specification of Solution Thermodynamic Properties </b>
*
* All solution properties are obtained from the standard state
* species functions, since there is only one species in the phase.
*
* <b> Application within %Kinetics Managers </b>
*
* The standard concentration is equal to 1.0. This means that the
* kinetics operator works on an (activities basis). Since this
* is a stoichiometric substance, this means that the concentration
* of this phase drops out of kinetics expressions.
*
* An example of a reaction using this is a sticking coefficient
* reaction of a substance in an ideal gas phase on a surface with a bulk phase
* species in this phase. In this case, the rate of progress for this
* reaction, \f$ R_s \f$, may be expressed via the following equation:
* \f[
* R_s = k_s C_{gas}
* \f]
* where the units for \f$ R_s \f$ are kmol m-2 s-1. \f$ C_{gas} \f$ has units
* of kmol m-3. Therefore, the kinetic rate constant, \f$ k_s \f$, has
* units of m s-1. Nowhere does the concentration of the bulk phase
* appear in the rate constant expression, since it's a stoichiometric
* phase and the activity is always equal to 1.0.
*
* <b> Instanteation of the Class </b>
*
* The constructor for this phase is NOT located in the default ThermoFactory
* for %Cantera. However, a new %StoichSubstanceSSTP may be created by
* the following code snippets:
*
* @code
* sprintf(file_ID,"%s#NaCl(S)", iFile);
* XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
* StoichSubstanceSSTP *solid = new StoichSubstanceSSTP(*xm);
* @endcode
*
* or by the following call to importPhase():
*
* @code
* sprintf(file_ID,"%s#NaCl(S)", iFile);
* XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
* StoichSubstanceSSTP solid;
* importPhase(*xm, &solid);
* @endcode
*
* <b> XML Example </b>
*
* The phase model name for this is called StoichSubstance. It must be supplied
* as the model attribute of the thermo XML element entry.
* Within the phase XML block,
* the density of the phase must be specified. An example of an XML file
* this phase is given below.
*
* @verbatim
<!-- phase NaCl(S) -->
<phase dim="3" id="NaCl(S)">
<elementArray datasrc="elements.xml">
Na Cl
</elementArray>
<speciesArray datasrc="#species_NaCl(S)"> NaCl(S) </speciesArray>
<thermo model="StoichSubstanceSSTP">
<density units="g/cm3">2.165</density>
</thermo>
<transport model="None"/>
<kinetics model="none"/>
</phase>
<!-- species definitions -->
<speciesData id="species_NaCl(S)">
<!-- species NaCl(S) -->
<species name="NaCl(S)">
<atomArray> Na:1 Cl:1 </atomArray>
<thermo>
<Shomate Pref="1 bar" Tmax="1075.0" Tmin="250.0">
<floatArray size="7">
50.72389, 6.672267, -2.517167,
10.15934, -0.200675, -427.2115,
130.3973
</floatArray>
</Shomate>
</thermo>
<density units="g/cm3">2.165</density>
</species>
</speciesData> @endverbatim
*
* The model attribute, "StoichSubstanceSSTP", on the thermo element identifies the phase as being
* a StoichSubstanceSSTP object.
*
* @ingroup thermoprops
*/
class MineralEQ3 : public StoichSubstanceSSTP {
public:
//! Default constructor for the StoichSubstanceSSTP class
MineralEQ3();
//! Construct and initialize a StoichSubstanceSSTP ThermoPhase object
//! directly from an asci input file
/*!
* @param infile name of the input file
* @param id name of the phase id in the file.
* If this is blank, the first phase in the file is used.
*/
MineralEQ3(std::string infile, std::string id = "");
//! Construct and initialize a StoichSubstanceSSTP ThermoPhase object
//! directly from an XML database
/*!
* @param phaseRef XML node pointing to a StoichSubstanceSSTP description
* @param id Id of the phase.
*/
MineralEQ3(XML_Node& phaseRef, std::string id = "");
//! Copy constructor
/*!
* @param right Object to be copied
*/
MineralEQ3(const MineralEQ3 &right);
//! Assignment operator
/*!
* @param right Object to be copied
*/
MineralEQ3 & operator=(const MineralEQ3 & right);
//! Destructor for the routine (virtual)
virtual ~MineralEQ3();
//! Duplication function
/*!
* This virtual function is used to create a duplicate of the
* current phase. It's used to duplicate the phase when given
* a ThermoPhase pointer to the phase.
*
* @return It returns a ThermoPhase pointer.
*/
ThermoPhase *duplMyselfAsThermoPhase() const;
/**
*
* @name Utilities
* @{
*/
/**
* Equation of state flag.
*
* Returns the value cStoichSubstance, defined in mix_defs.h.
*/
virtual int eosType() const;
/**
* @}
* @name Molar Thermodynamic Properties of the Solution
* @{
*/
/**
* @}
* @name Mechanical Equation of State
* @{
*/
//! Report the Pressure. Units: Pa.
/*!
* For an incompressible substance, the density is independent
* of pressure. This method simply returns the storred
* pressure value.
*/
virtual doublereal pressure() const;
//! Set the pressure at constant temperature. Units: Pa.
/*!
* For an incompressible substance, the density is
* independent of pressure. Therefore, this method only
* stores the specified pressure value. It does not
* modify the density.
*
* @param p Pressure (units - Pa)
*/
virtual void setPressure(doublereal p);
//! Returns the isothermal compressibility. Units: 1/Pa.
/*!
* The isothermal compressibility is defined as
* \f[
* \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
* \f]
*/
virtual doublereal isothermalCompressibility() const;
//! Return the volumetric thermal expansion coefficient. Units: 1/K.
/*!
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
virtual doublereal thermalExpansionCoeff() const ;
/**
* @}
* @name Activities, Standard States, and Activity Concentrations
*
* This section is largely handled by parent classes, since there
* is only one species. Therefore, the activity is equal to one.
* @{
*/
//! 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.
*
* For a stoichiomeetric substance, there is
* only one species, and the generalized concentration is 1.0.
*
* @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;
//! Return the standard concentration for the kth species
/*!
* The standard concentration \f$ C^0_k \f$ used to normalize
* the activity (i.e., generalized) concentration.
* This phase assumes that the kinetics operator works on an
* dimensionless basis. Thus, the standard concentration is
* equal to 1.0.
*
* @param k Optional parameter indicating the species. The default
* is to assume this refers to species 0.
* @return
* Returns The standard Concentration as 1.0
*/
virtual doublereal standardConcentration(int k=0) const;
//! Natural logarithm of the standard concentration of the kth species.
/*!
* @param k index of the species (defaults to zero)
*/
virtual doublereal logStandardConc(int k=0) const;
//! Get the array of chemical potentials at unit activity for the species
//! at their standard states at the current <I>T</I> and <I>P</I> of the solution.
/*!
* For a stoichiometric substance, there is no activity term in
* the chemical potential expression, and therefore the
* standard chemical potential and the chemical potential
* are both equal to the molar Gibbs function.
*
* These are the standard state chemical potentials \f$ \mu^0_k(T,P)
* \f$. The values are evaluated at the current
* temperature and pressure of the solution
*
* @param mu0 Output vector of chemical potentials.
* Length: m_kk.
*/
virtual void getStandardChemPotentials(doublereal* mu0) 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.
*
* The base %ThermoPhase class assigns thedefault quantities
* of (kmol/m3) for all species.
* Inherited classes are responsible for overriding the default
* values if necessary.
*
* @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(doublereal *uA, int k = 0,
int sizeUA = 6) const;
//@}
/// @name Partial Molar Properties of the Solution
///
/// These properties are handled by the parent class,
/// SingleSpeciesTP
//@{
//@}
/// @name Properties of the Standard State of the Species in the Solution
//@{
//! 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.
/*!
* @param hrt Output vector of nondimensional standard state enthalpies.
* Length: m_kk.
*/
virtual void getEnthalpy_RT(doublereal* hrt) const;
//! Get the array of nondimensional Entropy functions for the
//! standard state species at the current <I>T</I> and <I>P</I> of the solution.
/*!
* @param sr Output vector of nondimensional standard state entropies.
* Length: m_kk.
*/
virtual void getEntropy_R(doublereal* sr) const;
//! Get the nondimensional Gibbs functions for the species
//! in their standard states at the current <I>T</I> and <I>P</I> of the solution.
/*!
* @param grt Output vector of nondimensional standard state gibbs free energies
* Length: m_kk.
*/
virtual void getGibbs_RT(doublereal* grt) const;
//! Get the nondimensional Heat Capacities at constant
//! pressure for the species standard states
//! at the current <I>T</I> and <I>P</I> of the solution
/*!
* @param cpr Output vector of nondimensional standard state heat capacities
* Length: m_kk.
*/
virtual void getCp_R(doublereal* cpr) const;
//! Returns the vector of nondimensional Internal Energies of the standard
//! state species at the current <I>T</I> and <I>P</I> of the solution
/*!
* For an incompressible,
* stoichiometric substance, the molar internal energy is
* independent of pressure. Since the thermodynamic properties
* are specified by giving the standard-state enthalpy, the
* term \f$ P_{ref} \hat v\f$ is subtracted from the specified reference molar
* enthalpy to compute the standard state molar internal energy.
*
* @param urt output vector of nondimensional standard state
* internal energies of the species. Length: m_kk.
*/
virtual void getIntEnergy_RT(doublereal* urt) const;
//@}
/// @name Thermodynamic Values for the Species Reference States
//@{
//! Returns the vector of nondimensional
//! internal Energies of the reference state at the current temperature
//! of the solution and the reference pressure for each species.
/*!
* @param urt Output vector of nondimensional reference state
* internal energies of the species.
* Length: m_kk
*/
virtual void getIntEnergy_RT_ref(doublereal *urt) const;
/*
* ---- Critical State Properties
*/
/*
* ---- Saturation Properties
*/
//! Initialization of a HMWSoln phase using an xml file
/*!
* This routine is a precursor to initThermo(XML_Node*)
* routine, which does most of the work.
*
* @param inputFile XML file containing the description of the
* phase
*
* @param id Optional parameter identifying the name of the
* phase. If none is given, the first XML
* phase element will be used.
*/
void constructPhaseFile(std::string inputFile, std::string id);
//! Import and initialize a HMWSoln phase
//! specification in an XML tree into the current object.
/*!
* Here we read an XML description of the phase.
* We import descriptions of the elements that make up the
* species in a phase.
* We import information about the species, including their
* reference state thermodynamic polynomials. We then freeze
* the state of the species.
*
* Then, we read the species molar volumes from the xml
* tree to finish the initialization.
*
* @param phaseNode This object must be the phase node of a
* complete XML tree
* description of the phase, including all of the
* species data. In other words while "phase" must
* point to an XML phase object, it must have
* sibling nodes "speciesData" that describe
* the species in the phase.
*
* @param id ID of the phase. If nonnull, a check is done
* to see if phaseNode is pointing to the phase
* with the correct id.
*/
void constructPhaseXML(XML_Node& phaseNode, std::string id);
//! Internal initialization required after all species have
//! been added
/*!
* @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();
//! Initialize the phase parameters from an XML file.
/*!
* initThermoXML() (virtual from ThermoPhase)
*
* This gets called from importPhase(). It processes the XML file
* after the species are set up. This is the main routine for
* reading in activity coefficient parameters.
*
* @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.
*/
virtual void initThermoXML(XML_Node& phaseNode, std::string id);
//! Set the equation of state parameters
/*!
* @internal
* The number and meaning of these depends on the subclass.
*
* @param n number of parameters
* @param c array of \a n coefficients
* c[0] = density of phase [ kg/m3 ]
*/
virtual void setParameters(int n, doublereal * const c);
//! Get the equation of state parameters in a vector
/*!
* @internal
*
* @param n number of parameters
* @param c array of \a n coefficients
*
* For this phase:
* - n = 1
* - c[0] = density of phase [ kg/m3 ]
*/
virtual void getParameters(int &n, doublereal * const c) const;
//! 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. Note, this method is called before the phase is
* initialzed with elements and/or species.
*
* For this phase, the density of the phase is specified in this block.
*
* @param eosdata An XML_Node object corresponding to
* the "thermo" entry for this phase in the input file.
*
* eosdata points to the thermo block, and looks like this:
*
* @verbatim
<phase id="stoichsolid" >
<thermo model="StoichSubstance">
<density units="g/cm3">3.52</density>
</thermo>
</phase> @endverbatim
*
*/
virtual void setParametersFromXML(const XML_Node& eosdata);
doublereal LookupGe(const std::string& elemName);
void convertDGFormation();
protected:
//! Value of the Absolute Gibbs Free Energy NIST scale at T_r and P_r
/*!
* This is the NIST scale value of Gibbs free energy at T_r = 298.15
* and P_r = 1 atm.
*
* J kmol-1
*/
doublereal m_Mu0_tr_pr;
//! Input value of S_j at Tr and Pr (cal gmol-1 K-1)
/*!
* Tr = 298.15 Pr = 1 atm
*/
doublereal m_Entrop_tr_pr;
//! Input Value of deltaG of Formation at Tr and Pr (cal gmol-1)
/*!
* Tr = 298.15 Pr = 1 atm
*
* This is the delta G for the formation reaction of the
* ion from elements in their stable state at Tr, Pr.
*/
doublereal m_deltaG_formation_tr_pr;
//! Input Value of deltaH of Formation at Tr and Pr (cal gmol-1)
/*!
* Tr = 298.15 Pr = 1 atm
*
* This is the delta H for the formation reaction of the
* ion from elements in their stable state at Tr, Pr.
*/
doublereal m_deltaH_formation_tr_pr;
//! Input Value of the molar volume at T_r and P_r
/*!
* cm^3 / gmol
*/
doublereal m_V0_tr_pr;
//! a coefficient (cal gmol-1 K-1)
doublereal m_a;
//! b coefficient (cal gmol-1 K-2) x 10^3
doublereal m_b;
//! c coefficient (cal K gmol-1 K) x 10^-5
doublereal m_c;
};
}
#endif

24
Cantera/src/thermo/PDSS_HKFT.h
View file

@ -579,7 +579,7 @@ namespace Cantera {
doublereal m_r_e_j;
//! Value of deltaG of Formation at Tr and Pr (cal gmol-1)
//! Input value of deltaG of Formation at Tr and Pr (cal gmol-1)
/*!
* Tr = 298.15 Pr = 1 atm
*
@ -588,7 +588,7 @@ namespace Cantera {
*/
doublereal m_deltaG_formation_tr_pr;
//! Value of deltaH of Formation at Tr and Pr (cal gmol-1)
//! Input value of deltaH of Formation at Tr and Pr (cal gmol-1)
/*!
* Tr = 298.15 Pr = 1 atm
*
@ -597,40 +597,40 @@ namespace Cantera {
*/
doublereal m_deltaH_formation_tr_pr;
//! Value of the Absolute Gibbs Free Energy NIST scale at tr and pr
//! Value of the Absolute Gibbs Free Energy NIST scale at T_r and P_r
/*!
* this is the NIST scale value of Gibbs free energy at T_r = 298.15
* This is the NIST scale value of Gibbs free energy at T_r = 298.15
* and P_r = 1 atm.
*
* J kmol-1
*/
doublereal m_Mu0_tr_pr;
//! Value of S_j at Tr and Pr (cal gmol-1 K-1)
//! Input value of S_j at Tr and Pr (cal gmol-1 K-1)
/*!
* Tr = 298.15 Pr = 1 atm
*/
doublereal m_Entrop_tr_pr;
//! a1 coefficient (cal gmol-1 bar-1)
//! Input a1 coefficient (cal gmol-1 bar-1)
doublereal m_a1;
//! a2 coefficient (cal gmol-1)
//! Input a2 coefficient (cal gmol-1)
doublereal m_a2;
//! a3 coefficient (cal K gmol-1 bar-1)
//! Input a3 coefficient (cal K gmol-1 bar-1)
doublereal m_a3;
//! a4 coefficient (cal K gmol-1)
//! Input a4 coefficient (cal K gmol-1)
doublereal m_a4;
//! c1 coefficient (cal gmol-1 K-1)
//! Input c1 coefficient (cal gmol-1 K-1)
doublereal m_c1;
//! c2 coefficient (cal K gmol-1)
//! Input c2 coefficient (cal K gmol-1)
doublereal m_c2;
//! omega_pr_tr coefficient(cal gmol-1)
//! Input omega_pr_tr coefficient(cal gmol-1)
doublereal m_omega_pr_tr;
//! y = dZdT = 1/(esp*esp) desp/dT at 298.15 and 1 bar

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

@ -541,6 +541,7 @@ namespace Cantera {
bool ok = installSpecies(k, *s, *th, spth, spRuleList[i],
&phase, vp_spth, spfactory);
if (ok) {
th->saveSpeciesData(k, s);
++k;
}
}
@ -548,8 +549,6 @@ namespace Cantera {
// done adding species.
th->freezeSpecies();
th->saveSpeciesData(db);
// Perform any required subclass-specific initialization.
th->initThermo();

19
Cantera/src/thermo/ThermoPhase.cpp
View file

@ -881,7 +881,24 @@ namespace Cantera {
"Number of species is less than or equal to zero");
}
}
void ThermoPhase::saveSpeciesData(const int k, const XML_Node* const data) {
if ((int) m_speciesData.size() < (k + 1)) {
m_speciesData.resize(k+1, 0);
}
m_speciesData[k] = data;
}
//! Return a pointer to the XML tree containing the species
/// data for this phase.
const std::vector<const XML_Node *> & ThermoPhase::speciesData() const {
if ((int) m_speciesData.size() != m_kk) {
throw CanteraError("ThermoPhase::speciesData",
"m_speciesData is the wrong size");
}
return m_speciesData;
}
/*
* Set the thermodynamic state.
*/

33
Cantera/src/thermo/ThermoPhase.h
View file

@ -1780,21 +1780,11 @@ namespace Cantera {
* @param data Pointer to the XML_Node data containing
* information about the species in the phase.
*/
void saveSpeciesData(const XML_Node* data) {
m_speciesData = data;
}
void saveSpeciesData(const int k, const XML_Node* data);
/// Return a pointer to the XML tree containing the species
/// data for this phase.
const XML_Node* speciesData() {
if (!m_speciesData) {
throw CanteraError("ThermoPhase::speciesData",
"m_speciesData is NULL");
}
return m_speciesData;
}
//! Return a pointer to the vector of XML nodes containing the species
//! data for this phase.
const std::vector<const XML_Node *> & speciesData() const;
//! Install a species thermodynamic property manager.
/*!
@ -2012,11 +2002,13 @@ namespace Cantera {
*/
SpeciesThermo* m_spthermo;
/// Pointer to the XML tree containing the species
/// data for this phase. This is used to access data needed to
/// construct the transport manager and other properties
/// later in the initialization process.
const XML_Node* m_speciesData;
//! Vector of pointers to the species databases.
/*!
* This is used to access data needed to
* construct the transport manager and other properties
* later in the initialization process.
*/
std::vector<const XML_Node *> m_speciesData;
//! Index number of the phase
/*!
@ -2038,7 +2030,8 @@ namespace Cantera {
/// -> length equal to number of elements
vector_fp m_lambdaRRT;
//! Boolean indicating whether there is a valid set of saved element potentials for this phase
//! Boolean indicating whether there is a valid set of saved element potentials
//! for this phase
bool m_hasElementPotentials;
//! Boolean indicating whether a charge neutrality condition is a necessity