diff --git a/Cantera/src/thermo/DebyeHuckel.cpp b/Cantera/src/thermo/DebyeHuckel.cpp index ca4985702..e2a050846 100644 --- a/Cantera/src/thermo/DebyeHuckel.cpp +++ b/Cantera/src/thermo/DebyeHuckel.cpp @@ -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 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 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 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); } } } diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index c27f3b891..5504faa23 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -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 diff --git a/Cantera/src/thermo/HMWSoln_input.cpp b/Cantera/src/thermo/HMWSoln_input.cpp index 068fb9968..b50d7b6c5 100644 --- a/Cantera/src/thermo/HMWSoln_input.cpp +++ b/Cantera/src/thermo/HMWSoln_input.cpp @@ -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 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 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 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 diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index 3826c3398..15871d082 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -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) diff --git a/Cantera/src/thermo/MineralEQ3.cpp b/Cantera/src/thermo/MineralEQ3.cpp new file mode 100644 index 000000000..cdaef71b2 --- /dev/null +++ b/Cantera/src/thermo/MineralEQ3.cpp @@ -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 + +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 + * T and P 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 T and P 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: + * + * + * + * 3.52 + * + * + */ + 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; + } + +} + + diff --git a/Cantera/src/thermo/MineralEQ3.h b/Cantera/src/thermo/MineralEQ3.h new file mode 100644 index 000000000..3c8dbf64f --- /dev/null +++ b/Cantera/src/thermo/MineralEQ3.h @@ -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. + * + * + * Specification of Species Standard %State Properties + * + * 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. + * + * + * Specification of Solution Thermodynamic Properties + * + * All solution properties are obtained from the standard state + * species functions, since there is only one species in the phase. + * + * Application within %Kinetics Managers + * + * 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. + * + * Instanteation of the Class + * + * 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 + * + * XML Example + * + * 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 + + + + Na Cl + + NaCl(S) + + 2.165 + + + + + + + + + + Na:1 Cl:1 + + + + 50.72389, 6.672267, -2.517167, + 10.15934, -0.200675, -427.2115, + 130.3973 + + + + 2.165 + + @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 T and P 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 T and P 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 T and P 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 T and P 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 T and P 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 T and P 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 + + + 3.52 + + @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 diff --git a/Cantera/src/thermo/PDSS_HKFT.h b/Cantera/src/thermo/PDSS_HKFT.h index 3e71049c2..33c4a8ed9 100644 --- a/Cantera/src/thermo/PDSS_HKFT.h +++ b/Cantera/src/thermo/PDSS_HKFT.h @@ -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 diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index 2628a5e07..fd79bf219 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -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(); diff --git a/Cantera/src/thermo/ThermoPhase.cpp b/Cantera/src/thermo/ThermoPhase.cpp index a3998db62..61a04230b 100644 --- a/Cantera/src/thermo/ThermoPhase.cpp +++ b/Cantera/src/thermo/ThermoPhase.cpp @@ -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 & 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. */ diff --git a/Cantera/src/thermo/ThermoPhase.h b/Cantera/src/thermo/ThermoPhase.h index 40d0025ae..a6cbf0b25 100755 --- a/Cantera/src/thermo/ThermoPhase.h +++ b/Cantera/src/thermo/ThermoPhase.h @@ -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 & 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 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