diff --git a/Cantera/src/equil/MultiPhase.cpp b/Cantera/src/equil/MultiPhase.cpp index 96e49176c..53c17ef45 100644 --- a/Cantera/src/equil/MultiPhase.cpp +++ b/Cantera/src/equil/MultiPhase.cpp @@ -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(); diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index 4cdc7580e..a7f87e4ec 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -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) diff --git a/Cantera/src/thermo/RedlichKisterVPSSTP.cpp b/Cantera/src/thermo/RedlichKisterVPSSTP.cpp new file mode 100644 index 000000000..a97ee19d8 --- /dev/null +++ b/Cantera/src/thermo/RedlichKisterVPSSTP.cpp @@ -0,0 +1,1136 @@ +/** + * @file RedlichKisterVPSSTP.cpp + * Definitions for ThermoPhase object for phases which + * employ excess gibbs free energy formulations related to RedlichKister + * expansions (see \ref thermoprops + * and class \link Cantera::RedlichKisterVPSSTP RedlichKisterVPSSTP\endlink). + * + */ +/* + * Copywrite (2009) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Date: 2011-04-14 12:24:13 -0600 (Thu, 14 Apr 2011) $ + * $Revision: 713 $ + */ + + +#include "RedlichKisterVPSSTP.h" +#include "ThermoFactory.h" +#include + +using namespace std; + +namespace Cantera { + + static const double xxSmall = 1.0E-150; + //==================================================================================================================== + /* + * Default constructor. + * + */ + RedlichKisterVPSSTP::RedlichKisterVPSSTP() : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + m_pSpecies_A_ij(0), + m_pSpecies_B_ij(0), + m_N_ij(0), + m_HE_m_ij(0), + m_SE_m_ij(0), + formRedlichKister_(0), + formTempModel_(0), + dlnActCoeff_dX_() + { + } + //==================================================================================================================== + /* + * 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. + + */ + RedlichKisterVPSSTP::RedlichKisterVPSSTP(std::string inputFile, std::string id) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + m_pSpecies_A_ij(0), + m_pSpecies_B_ij(0), + m_N_ij(0), + m_HE_m_ij(0), + m_SE_m_ij(0), + formRedlichKister_(0), + formTempModel_(0), + dlnActCoeff_dX_() + { + constructPhaseFile(inputFile, id); + } + //==================================================================================================================== + RedlichKisterVPSSTP::RedlichKisterVPSSTP(XML_Node& phaseRoot, std::string id) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + m_pSpecies_A_ij(0), + m_pSpecies_B_ij(0), + m_N_ij(0), + m_HE_m_ij(0), + m_SE_m_ij(0), + formRedlichKister_(0), + formTempModel_(0), + dlnActCoeff_dX_() + { + constructPhaseXML(phaseRoot, id); + } + //==================================================================================================================== + // Special constructor for a hard-coded problem + /* + * + * LiKCl treating the PseudoBinary layer as passthrough. + * -> test to predict the eutectic and liquidus correctly. + * + */ + RedlichKisterVPSSTP::RedlichKisterVPSSTP(int testProb) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + m_pSpecies_A_ij(0), + m_pSpecies_B_ij(0), + m_N_ij(0), + m_HE_m_ij(0), + m_SE_m_ij(0), + formRedlichKister_(0), + formTempModel_(0), + dlnActCoeff_dX_() + { + constructPhaseFile("LiKCl_liquid.xml", ""); + numBinaryInteractions_ = 1; + + m_HE_m_ij.resize(0); + m_SE_m_ij.resize(0); + + vector_fp he(2); + he[0] = 0.0; + he[1] = 0.0; + vector_fp se(2); + se[0] = 0.0; + se[1] = 0.0; + + m_HE_m_ij.push_back(he); + m_SE_m_ij.push_back(se); + m_N_ij.push_back(1); + m_pSpecies_A_ij.resize(1); + m_pSpecies_B_ij.resize(1); + + int iLiLi = speciesIndex("LiLi"); + if (iLiLi < 0) { + throw CanteraError("RedlichKisterVPSSTP test1 constructor", + "Unable to find LiLi"); + } + m_pSpecies_A_ij[0] = iLiLi; + + + int iVLi = speciesIndex("VLi"); + if (iVLi < 0) { + throw CanteraError("RedlichKisterVPSSTP test1 constructor", + "Unable to find VLi"); + } + m_pSpecies_B_ij[0] = iVLi; + + + } + //==================================================================================================================== + /* + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + RedlichKisterVPSSTP::RedlichKisterVPSSTP(const RedlichKisterVPSSTP &b) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + m_pSpecies_A_ij(0), + m_pSpecies_B_ij(0), + m_N_ij(0), + m_HE_m_ij(0), + m_SE_m_ij(0), + formRedlichKister_(0), + formTempModel_(0), + dlnActCoeff_dX_() + { + RedlichKisterVPSSTP::operator=(b); + } + //==================================================================================================================== + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + RedlichKisterVPSSTP& RedlichKisterVPSSTP:: + operator=(const RedlichKisterVPSSTP &b) { + if (&b == this) { + return *this; + } + + GibbsExcessVPSSTP::operator=(b); + + numBinaryInteractions_ = b.numBinaryInteractions_ ; + m_pSpecies_A_ij = b.m_pSpecies_A_ij; + m_pSpecies_B_ij = b.m_pSpecies_B_ij; + m_N_ij = b.m_N_ij; + m_HE_m_ij = b.m_HE_m_ij; + m_SE_m_ij = b.m_SE_m_ij; + formRedlichKister_ = b.formRedlichKister_; + formTempModel_ = b.formTempModel_; + dlnActCoeff_dX_ = b.dlnActCoeff_dX_; + + return *this; + } + //==================================================================================================================== + /* + * + * ~RedlichKisterVPSSTP(): (virtual) + * + * Destructor: does nothing: + * + */ + RedlichKisterVPSSTP::~RedlichKisterVPSSTP() { + } + //==================================================================================================================== + /* + * This routine duplicates the current object and returns + * a pointer to ThermoPhase. + */ + ThermoPhase* + RedlichKisterVPSSTP::duplMyselfAsThermoPhase() const { + RedlichKisterVPSSTP* mtp = new RedlichKisterVPSSTP(*this); + return (ThermoPhase *) mtp; + } + + //==================================================================================================================== + // 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 RedlichKisterVPSSTP class also returns + * zero, as it is a non-complete class. + */ + int RedlichKisterVPSSTP::eosType() const { + return 0; + } + //==================================================================================================================== + /* + * Import, construct, and initialize a phase + * specification from an XML tree into the current object. + * + * This routine is a precursor to constructPhaseXML(XML_Node*) + * routine, which does most of the work. + * + * @param infile 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 RedlichKisterVPSSTP::constructPhaseFile(std::string inputFile, std::string id) { + + if ((int) inputFile.size() == 0) { + throw CanteraError("RedlichKisterVPSSTP:constructPhaseFile", + "input file is null"); + } + string path = findInputFile(inputFile); + std::ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("RedlichKisterVPSSTP:constructPhaseFile","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + XML_Node &phaseNode_XML = xml(); + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("RedlichKisterVPSSTP:constructPhaseFile", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + fxml_phase->copy(&phaseNode_XML); + constructPhaseXML(*fxml_phase, id); + delete fxml; + } + //==================================================================================================================== + /* + * Import, construct, and initialize a HMWSoln phase + * specification from an XML tree into the current object. + * + * Most of the work is carried out by the cantera base + * routine, importPhase(). That routine imports all of the + * species and element data, including the standard states + * of the species. + * + * Then, In this routine, we read the information + * particular to the specification of the activity + * coefficient model for the Pitzer parameterization. + * + * We also read information about the molar volumes of the + * standard states if present in the XML file. + * + * @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 RedlichKisterVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id) { + string stemp; + if ((int) id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("RedlichKisterVPSSTP::constructPhaseXML", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("RedlichKisterVPSSTP::constructPhaseXML", + "no thermo XML node"); + } + XML_Node& thermoNode = phaseNode.child("thermo"); + + /* + * Make sure that the thermo model is RedlichKister + */ + stemp = thermoNode.attrib("model"); + string formString = lowercase(stemp); + if (formString != "redlich-kister") { + throw CanteraError("RedlichKisterVPSSTP::constructPhaseXML", + "model name isn't Redlich-Kister: " + formString); + + } + + /* + * Call the Cantera importPhase() function. This will import + * all of the species into the phase. This will also handle + * all of the solvent and solute standard states + */ + bool m_ok = importPhase(phaseNode, this); + if (!m_ok) { + throw CanteraError("RedlichKisterVPSSTP::constructPhaseXML","importPhase failed "); + } + + } + //==================================================================================================================== + + + + /* + * ------------ Molar Thermodynamic Properties ---------------------- + */ + + + /* + * - Activities, Standard States, Activity Concentrations ----------- + */ + + // 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. + * + * Here we define the activity concentrations as equal + * to the activities, because the standard concentration is 1. + * + * @param c Output array of generalized concentrations. The + * units depend upon the implementation of the + * reaction rate expressions within the phase. + */ + void RedlichKisterVPSSTP::getActivityConcentrations(doublereal* c) const { + getActivities(c); + } + //==================================================================================================================== + doublereal RedlichKisterVPSSTP::standardConcentration(int k) const { + //err("standardConcentration"); + //return -1.0; + return 1.0; + } + //==================================================================================================================== + doublereal RedlichKisterVPSSTP::logStandardConc(int k) const { + //err("logStandardConc"); + //return -1.0; + return 0.0; + } + //==================================================================================================================== + // 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. + */ + void RedlichKisterVPSSTP::getActivityCoefficients(doublereal* ac) const { + /* + * Update the activity coefficients + */ + s_update_lnActCoeff(); + + /* + * take the exp of the internally storred coefficients. + */ + for (int k = 0; k < m_kk; k++) { + ac[k] = exp(lnActCoeff_Scaled_[k]); + } + } + //==================================================================================================================== + /* + * ------------ Partial Molar Properties of the Solution ------------ + */ + //==================================================================================================================== + void RedlichKisterVPSSTP::getElectrochemPotentials(doublereal* mu) const { + getChemPotentials(mu); + double ve = Faraday * electricPotential(); + for (int k = 0; k < m_kk; k++) { + mu[k] += ve*charge(k); + } + } + //==================================================================================================================== + void RedlichKisterVPSSTP::getChemPotentials(doublereal* mu) const { + doublereal xx; + /* + * First get the standard chemical potentials in + * molar form. + * -> this requires updates of standard state as a function + * of T and P + */ + getStandardChemPotentials(mu); + /* + * Update the activity coefficients + */ + s_update_lnActCoeff(); + /* + * + */ + doublereal RT = GasConstant * temperature(); + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(moleFractions_[k], xxSmall); + mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]); + } + } + //==================================================================================================================== + //Molar enthalpy. Units: J/kmol. + doublereal RedlichKisterVPSSTP::enthalpy_mole() const { + int kk = nSpecies(); + double hbar[kk], h = 0; + getPartialMolarEnthalpies(hbar); + for (int i = 0; i < kk; i++){ + h += moleFractions_[i]*hbar[i]; + } + return h; + } + //==================================================================================================================== + /// Molar entropy. Units: J/kmol. + doublereal RedlichKisterVPSSTP::entropy_mole() const { + int kk = nSpecies(); + double sbar[kk], s = 0; + getPartialMolarEntropies(sbar); + for (int i = 0; i < kk; i++){ + s += moleFractions_[i]*sbar[i]; + } + return s; + } + //==================================================================================================================== + /// Molar heat capacity at constant pressure. Units: J/kmol/K. + doublereal RedlichKisterVPSSTP::cp_mole() const { + int kk = nSpecies(); + double cpbar[kk], cp = 0; + getPartialMolarCp(cpbar); + for (int i = 0; i < kk; i++){ + cp += moleFractions_[i]*cpbar[i]; + } + return cp; + } + //==================================================================================================================== + /// Molar heat capacity at constant volume. Units: J/kmol/K. + doublereal RedlichKisterVPSSTP::cv_mole() const { + return cp_mole() - GasConstant; + } + //==================================================================================================================== + // 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] + * + */ + void RedlichKisterVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const { + /* + * Get the nondimensional standard state enthalpies + */ + getEnthalpy_RT(hbar); + /* + * dimensionalize it. + */ + double T = temperature(); + double RT = GasConstant * T; + for (int k = 0; k < m_kk; k++) { + hbar[k] *= RT; + } + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnActCoeff(); + s_update_dlnActCoeff_dT(); + double RTT = RT * T; + for (int k = 0; k < m_kk; k++) { + hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k]; + } + } + //==================================================================================================================== + // Returns an array of partial molar heat capacities 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} + * \f] + * + */ + void RedlichKisterVPSSTP::getPartialMolarCp(doublereal* cpbar) const { + /* + * Get the nondimensional standard state entropies + */ + getCp_R(cpbar); + double T = temperature(); + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnActCoeff(); + s_update_dlnActCoeff_dT(); + + for (int k = 0; k < m_kk; k++) { + cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k]; + } + /* + * dimensionalize it. + */ + for (int k = 0; k < m_kk; k++) { + cpbar[k] *= GasConstant; + } + } + //==================================================================================================================== + // 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} + * \f] + * + */ + void RedlichKisterVPSSTP::getPartialMolarEntropies(doublereal* sbar) const { + double xx; + /* + * Get the nondimensional standard state entropies + */ + getEntropy_R(sbar); + double T = temperature(); + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnActCoeff(); + s_update_dlnActCoeff_dT(); + + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(moleFractions_[k], xxSmall); + sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k]; + } + /* + * dimensionalize it. + */ + for (int k = 0; k < m_kk; k++) { + sbar[k] *= GasConstant; + } + } + + /* + * ------------ Partial Molar Properties of the Solution ------------ + */ + //==================================================================================================================== + // 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. + */ + void RedlichKisterVPSSTP::getPartialMolarVolumes(doublereal* vbar) const { + int iK; + /* + * Get the standard state values in m^3 kmol-1 + */ + getStandardVolumes(vbar); + for ( iK = 0; iK < m_kk; iK++ ){ + + vbar[iK] += 0.0; + } + } + //==================================================================================================================== + doublereal RedlichKisterVPSSTP::err(std::string msg) const { + throw CanteraError("RedlichKisterVPSSTP","Base class method " + +msg+" called. Equation of state type: "+int2str(eosType())); + return 0; + } + //==================================================================================================================== + /* + * @internal Initialize. This method is provided to allow + * subclasses to perform any initialization required after all + * species have been added. For example, it might be used to + * resize internal work arrays that must have an entry for + * each species. The base class implementation does nothing, + * and subclasses that do not require initialization do not + * need to overload this method. When importing a CTML phase + * description, this method is called just prior to returning + * from function importPhase. + * + * @see importCTML.cpp + */ + void RedlichKisterVPSSTP::initThermo() { + initLengths(); + GibbsExcessVPSSTP::initThermo(); + } + //==================================================================================================================== + // Initialize lengths of local variables after all species have + // been identified. + void RedlichKisterVPSSTP::initLengths() { + m_kk = nSpecies(); + dlnActCoeffdlnN_.resize(m_kk, m_kk); + } + //==================================================================================================================== + /* + * initThermoXML() (virtual from ThermoPhase) + * Import and initialize a ThermoPhase object + * + * @param phaseNode This object must be the phase node of a complete XML tree + * description of the phase, including all of the species data. In other words while "phase" must + * point to an XML phase object, it must have sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done to see if phaseNode is pointing to the phase + * with the correct id. + */ + void RedlichKisterVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { + std::string subname = "RedlichKisterVPSSTP::initThermoXML"; + std::string stemp; + + /* + * Check on the thermo field. Must have: + * + */ + + XML_Node& thermoNode = phaseNode.child("thermo"); + std::string mStringa = thermoNode.attrib("model"); + std::string mString = lowercase(mStringa); + if (mString != "redlich-kister") { + throw CanteraError(subname.c_str(), + "Unknown thermo model: " + mStringa + " - This object only knows \"Redlich-Kister\" "); + } + + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + XML_Node *acNodePtr = 0; + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& acNode = thermoNode.child("activityCoefficients"); + acNodePtr = &acNode; + std::string mStringa = acNode.attrib("model"); + std::string mString = lowercase(mStringa); + if (mString != "redlich-kister") { + throw CanteraError(subname.c_str(), + "Unknown activity coefficient model: " + mStringa); + } + int n = acNodePtr->nChildren(); + for (int i = 0; i < n; i++) { + XML_Node &xmlACChild = acNodePtr->child(i); + stemp = xmlACChild.name(); + std::string nodeName = lowercase(stemp); + /* + * Process a binary salt field, or any of the other XML fields + * that make up the Pitzer Database. Entries will be ignored + * if any of the species in the entry isn't in the solution. + */ + if (nodeName == "binaryneutralspeciesparameters") { + readXMLBinarySpecies(xmlACChild); + } + } + } + /* + * Go down the chain + */ + GibbsExcessVPSSTP::initThermoXML(phaseNode, id); + } + //=================================================================================================================== + // Update the activity coefficients + /* + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + * + */ + void RedlichKisterVPSSTP::s_update_lnActCoeff() const { + int iA, iB, m, k; + doublereal XA, XB; + doublereal T = temperature(); + doublereal RT = GasConstant * T; + + fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + + for (int i = 0; i < numBinaryInteractions_; i++) { + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + doublereal deltaX = XA - XB; + int N = m_N_ij[i]; + vector_fp &he_vec = m_HE_m_ij[i]; + vector_fp &se_vec = m_SE_m_ij[i]; + doublereal poly = 1.0; + doublereal polyMm1 = 1.0; + doublereal sum = 0.0; + doublereal sumMm1 = 0.0; + doublereal sum2 = 0.0; + for (m = 0; m < N; m++) { + doublereal A_ge = he_vec[m] - T * se_vec[m]; + sum += A_ge * poly; + sum2 += A_ge * (m + 1) * poly; + poly *= deltaX; + if (m >= 1) { + sumMm1 += (A_ge * polyMm1 * m); + polyMm1 *= deltaX; + } + } + doublereal oneMXA = 1.0 - XA; + doublereal oneMXB = 1.0 - XB; + for (k = 0; k < m_kk; k++) { + if (iA == k) { + lnActCoeff_Scaled_[k] += (oneMXA * XB * sum) + (XA * XB * sumMm1 * (oneMXA + XB)); + } else if (iB == k) { + lnActCoeff_Scaled_[k] += (oneMXB * XA * sum) + (XA * XB * sumMm1 * (-oneMXB - XA)); + } else { + lnActCoeff_Scaled_[k] += -(XA * XB * sum2); + } + } + // Debug against formula in literature +#ifdef DEBUG_MODE_NOT + double lnA = 0.0; + double lnB = 0.0; + double polyk = 1.0; + double fac = 2.0 * XA - 1.0; + for (m = 0; m < N; m++) { + doublereal A_ge = he_vec[m] - T * se_vec[m]; + lnA += A_ge * oneMXA * oneMXA * polyk * (1.0 + 2.0 * XA * m / fac); + lnB += A_ge * XA * XA * polyk * (1.0 - 2.0 * oneMXA * m / fac); + polyk *= fac; + } + // This gives the same result as above + // printf("RT lnActCoeff_Scaled_[iA] = %15.8E , lnA = %15.8E\n", lnActCoeff_Scaled_[iA], lnA); + // printf("RT lnActCoeff_Scaled_[iB] = %15.8E , lnB = %15.8E\n", lnActCoeff_Scaled_[iB], lnB); + +#endif + + } + for (k = 0; k < m_kk; k++) { + lnActCoeff_Scaled_[k] /= RT; + } + } + //=================================================================================================================== + // Update the derivative of the log of the activity coefficients wrt T + /* + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + * + + */ + void RedlichKisterVPSSTP::s_update_dlnActCoeff_dT() const { + int iA, iB, m, k; + doublereal XA, XB; + // doublereal T = temperature(); + + fvo_zero_dbl_1(dlnActCoeffdT_Scaled_, m_kk); + fvo_zero_dbl_1(d2lnActCoeffdT2_Scaled_, m_kk); + + for (int i = 0; i < numBinaryInteractions_; i++) { + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + doublereal deltaX = XA - XB; + int N = m_N_ij[i]; + doublereal poly = 1.0; + doublereal sum = 0.0; + + vector_fp &se_vec = m_SE_m_ij[i]; + doublereal sumMm1 = 0.0; + doublereal polyMm1 = 1.0; + doublereal sum2 = 0.0; + for (m = 0; m < N; m++) { + doublereal A_ge = - se_vec[m]; + sum += A_ge * poly; + sum2 += A_ge * (m + 1) * poly; + poly *= deltaX; + if (m >= 1) { + sumMm1 += (A_ge * polyMm1 * m); + polyMm1 *= deltaX; + } + } + doublereal oneMXA = 1.0 - XA; + doublereal oneMXB = 1.0 - XB; + for (k = 0; k < m_kk; k++) { + if (iA == k) { + dlnActCoeffdT_Scaled_[k] += (oneMXA * XB * sum) + (XA * XB * sumMm1 * (oneMXA + XB)); + } else if (iB == k) { + dlnActCoeffdT_Scaled_[k] += (oneMXB * XA * sum) + (XA * XB * sumMm1 * (-oneMXB - XA)); + } else { + dlnActCoeffdT_Scaled_[k] += -(XA * XB * sum2); + } + } + } + } + //==================================================================================================================== + void RedlichKisterVPSSTP::getdlnActCoeffdT(doublereal *dlnActCoeffdT) const { + s_update_dlnActCoeff_dT(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdT[k] = dlnActCoeffdT_Scaled_[k]; + } + } + //==================================================================================================================== + void RedlichKisterVPSSTP::getd2lnActCoeffdT2(doublereal *d2lnActCoeffdT2) const { + s_update_dlnActCoeff_dT(); + for (int k = 0; k < m_kk; k++) { + d2lnActCoeffdT2[k] = d2lnActCoeffdT2_Scaled_[k]; + } + } + //==================================================================================================================== + void RedlichKisterVPSSTP::s_update_dlnActCoeff_dX_() const { + + + int iA, iB, m, k; + doublereal XA, XB; + doublereal T = temperature(); + + dlnActCoeff_dX_.zero(); + + for (int i = 0; i < numBinaryInteractions_; i++) { + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + doublereal deltaX = XA - XB; + int N = m_N_ij[i]; + doublereal poly = 1.0; + doublereal sum = 0.0; + vector_fp &he_vec = m_HE_m_ij[i]; + vector_fp &se_vec = m_SE_m_ij[i]; + doublereal sumMm1 = 0.0; + doublereal polyMm1 = 1.0; + doublereal polyMm2 = 1.0; + doublereal sum2 = 0.0; + doublereal sum2Mm1 = 0.0; + doublereal sumMm2 = 0.0; + for (m = 0; m < N; m++) { + doublereal A_ge = he_vec[m] - T * se_vec[m]; + sum += A_ge * poly; + sum2 += A_ge * (m + 1) * poly; + poly *= deltaX; + if (m >= 1) { + sumMm1 += (A_ge * polyMm1 * m); + sum2Mm1 += (A_ge * polyMm1 * m * (1.0 + m)); + polyMm1 *= deltaX; + } + if (m >= 2) { + sumMm2 += (A_ge * polyMm2 * m * (m - 1.0)); + polyMm2 *= deltaX; + } + } + + for (k = 0; k < m_kk; k++) { + if (iA == k) { + + dlnActCoeff_dX_(k, iA) += (- XB * sum + (1.0 - XA) * XB * sumMm1 + + XB * sumMm1 * (1.0 - 2.0 * XA + XB) + + XA * XB * sumMm2 * (1.0 - XA + XB)); + + dlnActCoeff_dX_(k, iB) += ((1.0 - XA) * sum - (1.0 - XA) * XB * sumMm1 + + XA * sumMm1 * (1.0 + 2.0 * XB - XA) + - XA * XB * sumMm2 * (1.0 - XA + XB)); + + } else if (iB == k) { + + dlnActCoeff_dX_(k, iA) += ((1.0 - XB) * sum + (1.0 - XA) * XB * sumMm1 + + XB * sumMm1 * (1.0 - 2.0 * XA + XB) + + XA * XB * sumMm2 * (1.0 - XA + XB)); + + dlnActCoeff_dX_(k, iB) += (- XA * sum - (1.0 - XB) * XA * sumMm1 + + XA * sumMm1 * (XB - XA - (1.0 - XB)) + - XA * XB * sumMm2 * (-XA - (1.0 - XB))); + } else { + + dlnActCoeff_dX_(k, iA) += ( - XB * sum2 - XA * XB * sum2Mm1); + + dlnActCoeff_dX_(k, iB) += ( - XA * sum2 + XA * XB * sum2Mm1); + + } + } + } + } + //==================================================================================================================== + // 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. + */ + void RedlichKisterVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal * const dXds, + doublereal *dlnActCoeffds) const { + s_update_dlnActCoeff_dT(); + s_update_dlnActCoeff_dX_(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffds[k] = dlnActCoeffdT_Scaled_[k] * dTds; + for (int l = 0; l < m_kk; l++) { + dlnActCoeffds[k] += dlnActCoeff_dX_(k, l) * dXds[l]; + } + } + } + + //==================================================================================================================== + void RedlichKisterVPSSTP::getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const { + s_update_dlnActCoeff_dX_(); + for (int l = 0; l < m_kk; l++) { + dlnActCoeffdlnN_diag[l] = dlnActCoeff_dX_(l, l); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdlnN_diag[k] -= dlnActCoeff_dX_(l, k) * moleFractions_[k]; + } + } + } + //==================================================================================================================== + void RedlichKisterVPSSTP::getdlnActCoeffdlnX_diag(doublereal *dlnActCoeffdlnX_diag) const { + s_update_dlnActCoeff_dX_(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k]; + } + } + //==================================================================================================================== + void RedlichKisterVPSSTP::getdlnActCoeffdlnN(const int ld, doublereal *dlnActCoeffdlnN) { + s_update_dlnActCoeff_dX_(); + double *data = & dlnActCoeffdlnN_(0,0); + for (int k = 0; k < m_kk; k++) { + for (int m = 0; m < m_kk; m++) { + dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m]; + } + } + } + //==================================================================================================================== + void RedlichKisterVPSSTP::resizeNumInteractions(const int num) { + numBinaryInteractions_ = num; + m_pSpecies_A_ij.resize(num, -1); + m_pSpecies_B_ij.resize(num, -1); + m_N_ij.resize(num, -1); + m_HE_m_ij.resize(num); + m_SE_m_ij.resize(num); + dlnActCoeff_dX_.resize(num, num, 0.0); + } + //==================================================================================================================== + // Process an XML node called "binaryNeutralSpeciesParameters" + /* + * This node contains all of the parameters necessary to describe the RedlichKister Interaction for + * a single binary interaction. This function reads the XML file and writes the coefficients + * it finds to an internal data structures. + */ + void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node &xmLBinarySpecies) { + std::string xname = xmLBinarySpecies.name(); + if (xname != "binaryNeutralSpeciesParameters") { + throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", + "Incorrect name for processing this routine: " + xname); + } + double *charge = DATA_PTR(m_speciesCharge); + std::string stemp; + int nParamsFound = 0; + int Npoly = 0; + vector_fp hParams, sParams, vParams; + std::string iName = xmLBinarySpecies.attrib("speciesA"); + if (iName == "") { + throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", "no speciesA attrib"); + } + std::string jName = xmLBinarySpecies.attrib("speciesB"); + if (jName == "") { + throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", "no speciesB attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species. This means that the interaction doesn't occur for the current + * implementation of the phase. + */ + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + string ispName = speciesName(iSpecies); + if (charge[iSpecies] != 0) { + throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", "speciesA charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + std::string jspName = speciesName(jSpecies); + if (charge[jSpecies] != 0) { + throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", "speciesB charge problem"); + } + /* + * Ok we have found a valid interaction + */ + numBinaryInteractions_++; + int iSpot = numBinaryInteractions_ - 1; + m_pSpecies_A_ij.resize(numBinaryInteractions_); + m_pSpecies_B_ij.resize(numBinaryInteractions_); + m_pSpecies_A_ij[iSpot] = iSpecies; + m_pSpecies_B_ij[iSpot] = jSpecies; + + int num = xmLBinarySpecies.nChildren(); + for (int iChild = 0; iChild < num; iChild++) { + XML_Node &xmlChild = xmLBinarySpecies.child(iChild); + stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + /* + * Process the binary species interaction child elements + */ + if (nodeName == "excessenthalpy") { + /* + * Get the string containing all of the values + */ + ctml::getFloatArray(xmlChild, hParams, true, "toSI", "excessEnthalpy"); + nParamsFound = hParams.size(); + if (nParamsFound > Npoly) { + Npoly = nParamsFound; + } + + } + + if (nodeName == "excessentropy") { + /* + * Get the string containing all of the values + */ + ctml::getFloatArray(xmlChild, sParams, true, "toSI", "excessEntropy"); + nParamsFound = sParams.size(); + if (nParamsFound > Npoly) { + Npoly = nParamsFound; + } + } + } + hParams.resize(Npoly, 0.0); + sParams.resize(Npoly, 0.0); + m_HE_m_ij.push_back(hParams); + m_SE_m_ij.push_back(sParams); + m_N_ij.push_back(Npoly); + resizeNumInteractions(numBinaryInteractions_); + } + //==================================================================================================================== +#ifdef DEBUG_MODE + void RedlichKisterVPSSTP::Vint(double &VintOut, double &voltsOut) { + int iA, iB, m; + doublereal XA, XB; + doublereal T = temperature(); + doublereal RT = GasConstant * T; + double Volts = 0.0; + + fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + + for (int i = 0; i < numBinaryInteractions_; i++) { + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + if (XA <= 1.0E-14) { + XA = 1.0E-14; + } + if (XA >= (1.0 - 1.0E-14)) { + XA = 1.0 - 1.0E-14; + } + + int N = m_N_ij[i]; + vector_fp &he_vec = m_HE_m_ij[i]; + vector_fp &se_vec = m_SE_m_ij[i]; + double fac = 2.0 * XA - 1.0; + if (fabs(fac) < 1.0E-13) { + fac = 1.0E-13; + } + double polykp1 = fac; + double poly1mk = fac; + + for (m = 0; m < N; m++) { + doublereal A_ge = he_vec[m]; + Volts += A_ge * ( polykp1 - (2.0 * XA * m * (1.0-XA) ) / poly1mk ); + polykp1 *= fac; + poly1mk /= fac; + } + } + Volts /= Faraday; + + double termp = RT * log((1.0 - XA)/XA) / Faraday; + + VintOut = Volts; + voltsOut = Volts + termp; + } + #endif + //==================================================================================================================== +} + diff --git a/Cantera/src/thermo/RedlichKisterVPSSTP.h b/Cantera/src/thermo/RedlichKisterVPSSTP.h new file mode 100644 index 000000000..d11e9e2bd --- /dev/null +++ b/Cantera/src/thermo/RedlichKisterVPSSTP.h @@ -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. + * + * + *
+ *

Specification of Species Standard %State Properties

+ *
+ * + * 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. + * + * + *
+ *

Specification of Solution Thermodynamic Properties

+ *
+ * + * The molar excess Gibbs free energy is given by the following formula which is a sum over interactions i. + * Each of the interactions are binary interactions involving two of the species in the phase, denoted, Ai + * and Bi. + * 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 k 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 k 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 k 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 k 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 k 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] + * + *
+ *

%Application within %Kinetics Managers

+ *
+ * + * \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 k is independent of k 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. + * + * + *
+ *

Instantiation of the Class

+ *
+ * + * + * 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 (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 + * + *
+ *

XML Example

+ *
+ * An example of an XML Element named phase setting up a IdealGasPhase + * object named silane is given below. + * + * + * @verbatim + + + Si H He + + H2 H HE SIH4 SI SIH SIH2 SIH3 H3SISIH SI2H6 + H2SISIH2 SI3H8 SI2 SI3 + + + + + + + @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 mth + * species with respect to the number of moles of the kth 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 + + + + + diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index ff9480d97..0cf6338fc 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -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; } diff --git a/Cantera/src/thermo/ThermoFactory.h b/Cantera/src/thermo/ThermoFactory.h index 01ec78d12..b4af0abef 100644 --- a/Cantera/src/thermo/ThermoFactory.h +++ b/Cantera/src/thermo/ThermoFactory.h @@ -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 diff --git a/Cantera/src/thermo/mix_defs.h b/Cantera/src/thermo/mix_defs.h index 1aef185ea..a74fb4ccf 100644 --- a/Cantera/src/thermo/mix_defs.h +++ b/Cantera/src/thermo/mix_defs.h @@ -79,6 +79,8 @@ namespace Cantera { const int cMargulesVPSSTP = 301; + const int cRedlichKisterVPSSTP = 303; + const int cMolarityIonicVPSSTP = 401; const int cMixedSolventElectrolyte = 402; diff --git a/docs/Cantera.cfg.in b/docs/Cantera.cfg.in index f3f3a4dab..5b260f5d0 100755 --- a/docs/Cantera.cfg.in +++ b/docs/Cantera.cfg.in @@ -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 \