diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index b359d5579..690877d38 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -67,7 +67,7 @@ ELECTRO_OBJ = MolalityVPSSTP.o VPStandardStateTP.o \ PDSS.o PDSS_Water.o PDSS_HKFT.o \ HMWSoln.o HMWSoln_input.o DebyeHuckel.o \ WaterSSTP.o MetalSHEelectrons.o \ - VPSSMgr_Water_ConstVol.o VPSSMgr_Water_HKFT.o + VPSSMgr_Water_ConstVol.o VPSSMgr_Water_HKFT.o PhaseCombo_Interaction.o ELECTRO_H = MolalityVPSSTP.h VPStandardStateTP.h \ IdealMolalSoln.h \ @@ -75,7 +75,7 @@ ELECTRO_H = MolalityVPSSTP.h VPStandardStateTP.h \ PDSS.h PDSS_Water.h PDSS_HKFT.h \ HMWSoln.h electrolytes.h \ DebyeHuckel.h WaterSSTP.h MetalSHEelectrons.h VPSSMgr_Water_HKFT.h \ - VPSSMgr_Water_ConstVol.h + VPSSMgr_Water_ConstVol.h PhaseCombo_Interaction.h endif ifeq ($(do_issp),1) ISSP_OBJ = IdealSolidSolnPhase.o StoichSubstanceSSTP.o SingleSpeciesTP.o MineralEQ3.o \ diff --git a/Cantera/src/thermo/MargulesVPSSTP.cpp b/Cantera/src/thermo/MargulesVPSSTP.cpp index f039c561d..b3e85eda7 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.cpp +++ b/Cantera/src/thermo/MargulesVPSSTP.cpp @@ -715,7 +715,7 @@ namespace Cantera { if (thermoNode.hasChild("activityCoefficients")) { XML_Node& acNode = thermoNode.child("activityCoefficients"); acNodePtr = &acNode; - string mStringa = thermoNode.attrib("model"); + string mStringa = acNode.attrib("model"); string mString = lowercase(mStringa); if (mString != "margules") { throw CanteraError(subname.c_str(), diff --git a/Cantera/src/thermo/PhaseCombo_Interaction.cpp b/Cantera/src/thermo/PhaseCombo_Interaction.cpp new file mode 100644 index 000000000..f3ac98ab1 --- /dev/null +++ b/Cantera/src/thermo/PhaseCombo_Interaction.cpp @@ -0,0 +1,1276 @@ +/** + * @file + * + */ +/* + * 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: 2010-11-12 14:37:41 -0700 (Fri, 12 Nov 2010) $ + * $Revision: 641 $ + */ + + +#include "PhaseCombo_Interaction.h" +#include "ThermoFactory.h" +#include + +using namespace std; + +namespace Cantera { + + static const double xxSmall = 1.0E-150; + //==================================================================================================================== + /* + * Default constructor. + * + * HKM - Checked for Transition + */ + PhaseCombo_Interaction::PhaseCombo_Interaction() : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + } + //==================================================================================================================== + /* + * 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. + * + * HKM - Checked for Transition + */ + PhaseCombo_Interaction::PhaseCombo_Interaction(std::string inputFile, std::string id) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + constructPhaseFile(inputFile, id); + } + //==================================================================================================================== + // + /* + * + * HKM - Checked for Transition + */ + PhaseCombo_Interaction::PhaseCombo_Interaction(XML_Node& phaseRoot, std::string id) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + constructPhaseXML(phaseRoot, id); + } + + //==================================================================================================================== + /* + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + * + * HKM - Checked for Transition + */ + PhaseCombo_Interaction::PhaseCombo_Interaction(const PhaseCombo_Interaction &b) : + GibbsExcessVPSSTP() + { + PhaseCombo_Interaction::operator=(b); + } + //==================================================================================================================== + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + * + * HKM - Checked for Transition + */ + PhaseCombo_Interaction& PhaseCombo_Interaction:: + operator=(const PhaseCombo_Interaction &b) { + if (&b == this) { + return *this; + } + + GibbsExcessVPSSTP::operator=(b); + + numBinaryInteractions_ = b.numBinaryInteractions_ ; + m_HE_b_ij = b.m_HE_b_ij; + m_HE_c_ij = b.m_HE_c_ij; + m_HE_d_ij = b.m_HE_d_ij; + m_SE_b_ij = b.m_SE_b_ij; + m_SE_c_ij = b.m_SE_c_ij; + m_SE_d_ij = b.m_SE_d_ij; + m_VHE_b_ij = b.m_VHE_b_ij; + m_VHE_c_ij = b.m_VHE_c_ij; + m_VHE_d_ij = b.m_VHE_d_ij; + m_VSE_b_ij = b.m_VSE_b_ij; + m_VSE_c_ij = b.m_VSE_c_ij; + m_VSE_d_ij = b.m_VSE_d_ij; + m_pSpecies_A_ij = b.m_pSpecies_A_ij; + m_pSpecies_B_ij = b.m_pSpecies_B_ij; + formMargules_ = b.formMargules_; + formTempModel_ = b.formTempModel_; + + return *this; + } + //==================================================================================================================== + /** + * + * ~PhaseCombo_Interaction(): (virtual) + * + * Destructor: does nothing: + * + * HKM - Checked for Transition + */ + PhaseCombo_Interaction::~PhaseCombo_Interaction() { + } + //==================================================================================================================== + /* + * This routine duplicates the current object and returnsa pointer to ThermoPhase. + * + * HKM - Checked for Transition + */ + ThermoPhase* + PhaseCombo_Interaction::duplMyselfAsThermoPhase() const { + PhaseCombo_Interaction* mtp = new PhaseCombo_Interaction(*this); + return (ThermoPhase *) mtp; + } + //==================================================================================================================== + // Special constructor for a hard-coded problem + /* + * + * LiKCl treating the PseudoBinary layer as passthrough. + * -> test to predict the eutectic and liquidus correctly. + * + */ + PhaseCombo_Interaction::PhaseCombo_Interaction(int testProb) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + + + constructPhaseFile("PhaseCombo_Interaction.xml", ""); + + + numBinaryInteractions_ = 1; + + m_HE_b_ij.resize(1); + m_HE_c_ij.resize(1); + m_HE_d_ij.resize(1); + + m_SE_b_ij.resize(1); + m_SE_c_ij.resize(1); + m_SE_d_ij.resize(1); + + m_VHE_b_ij.resize(1); + m_VHE_c_ij.resize(1); + m_VHE_d_ij.resize(1); + + m_VSE_b_ij.resize(1); + m_VSE_c_ij.resize(1); + m_VSE_d_ij.resize(1); + + m_pSpecies_A_ij.resize(1); + m_pSpecies_B_ij.resize(1); + + + + m_HE_b_ij[0] = -17570E3; + m_HE_c_ij[0] = -377.0E3; + m_HE_d_ij[0] = 0.0; + + m_SE_b_ij[0] = -7.627E3; + m_SE_c_ij[0] = 4.958E3; + m_SE_d_ij[0] = 0.0; + + + int iLiT = speciesIndex("LiTFe1S2(S)"); + if (iLiT < 0) { + throw CanteraError("PhaseCombo_Interaction test1 constructor", + "Unable to find LiTFe1S2(S)"); + } + m_pSpecies_A_ij[0] = iLiT; + + + int iLi2 = speciesIndex("Li2Fe1S2(S)"); + if (iLi2 < 0) { + throw CanteraError("PhaseCombo_Interaction test1 constructor", + "Unable to find Li2Fe1S2(S)"); + } + m_pSpecies_B_ij[0] = iLi2; + throw CanteraError("", "unimplemented"); + } + //==================================================================================================================== + + /* + * -------------- 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 PhaseCombo_Interaction class also returns + * zero, as it is a non-complete class. + */ + int PhaseCombo_Interaction::eosType() const { + return cPhaseCombo_Interaction; + } + //==================================================================================================================== + /* + * 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. + * + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::constructPhaseFile(std::string inputFile, std::string id) { + + if ((int) inputFile.size() == 0) { + throw CanteraError("PhaseCombo_Interaction:constructPhaseFile", + "input file is null"); + } + string path = findInputFile(inputFile); + std::ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("PhaseCombo_Interaction: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("PhaseCombo_Interaction: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. + * + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::constructPhaseXML(XML_Node& phaseNode, std::string id) { + string stemp; + if ((int) id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("PhaseCombo_Interaction::constructPhaseXML", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("PhaseCombo_Interaction::constructPhaseXML", + "no thermo XML node"); + } + XML_Node& thermoNode = phaseNode.child("thermo"); + + /* + * Make sure that the thermo model is PhaseCombo_Interaction + */ + stemp = thermoNode.attrib("model"); + string formString = lowercase(stemp); + if (formString != "phasecombo_interaction") { + throw CanteraError("PhaseCombo_Interaction::constructPhaseXML", + "model name isn't PhaseCombo_Interaction: " + formString); + } + + /* + * Call the Cantera importPhase() function. This will import + * all of the species into the phase. This will also handle + * all of the species standard states + */ + bool m_ok = importPhase(phaseNode, this); + if (!m_ok) { + throw CanteraError("PhaseCombo_Interaction::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 PhaseCombo_Interaction::getActivityConcentrations(doublereal* c) const { + getActivities(c); + } + //==================================================================================================================== + doublereal PhaseCombo_Interaction::standardConcentration(int k) const { + //err("standardConcentration"); + //return -1.0; + return 1.0; + } + //==================================================================================================================== + doublereal PhaseCombo_Interaction::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 PhaseCombo_Interaction::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 PhaseCombo_Interaction::getElectrochemPotentials(doublereal* mu) const { + getChemPotentials(mu); + double ve = Faraday * electricPotential(); + for (int k = 0; k < m_kk; k++) { + mu[k] += ve*charge(k); + } + } + + //==================================================================================================================== + void PhaseCombo_Interaction::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 PhaseCombo_Interaction::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 PhaseCombo_Interaction::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 PhaseCombo_Interaction::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 PhaseCombo_Interaction::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 PhaseCombo_Interaction::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 PhaseCombo_Interaction::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 PhaseCombo_Interaction::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 PhaseCombo_Interaction::getPartialMolarVolumes(doublereal* vbar) const { + + int iA, iB, iK, delAK, delBK; + double XA, XB, XK, g0 , g1; + double T = temperature(); + + /* + * Get the standard state values in m^3 kmol-1 + */ + getStandardVolumes(vbar); + + for ( iK = 0; iK < m_kk; iK++ ){ + delAK = 0; + delBK = 0; + XK = moleFractions_[iK]; + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]); + g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]); + + vbar[iK] += XA*XB*(g0+g1*XB)+((delAK-XA)*XB+XA*(delBK-XB))*(g0+g1*XB)+XA*XB*(delBK-XB)*g1; + } + } + } + //==================================================================================================================== + doublereal PhaseCombo_Interaction::err(std::string msg) const { + throw CanteraError("PhaseCombo_Interaction","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 PhaseCombo_Interaction::initThermo() { + initLengths(); + GibbsExcessVPSSTP::initThermo(); + } + + //==================================================================================================================== + // Initialize lengths of local variables after all species have + // been identified. + void PhaseCombo_Interaction::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 PhaseCombo_Interaction::initThermoXML(XML_Node& phaseNode, std::string id) { + string subname = "PhaseCombo_Interaction::initThermoXML"; + string stemp; + + /* + * Check on the thermo field. Must have: + * + */ + + XML_Node& thermoNode = phaseNode.child("thermo"); + string mStringa = thermoNode.attrib("model"); + string mString = lowercase(mStringa); + if (mString != "phasecombo_interaction") { + throw CanteraError(subname.c_str(), "Unknown thermo model: " + mStringa); + } + + + /* + * 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; + string mStringa = acNode.attrib("model"); + string mString = lowercase(mStringa); + if (mString != "margules") { + 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(); + 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 + * + * he = X_A X_B(B + C X_B) + * + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::s_update_lnActCoeff() const { + int iA, iB, iK, delAK, delBK; + doublereal XA, XB, g0 , g1; + doublereal xx; + doublereal T = temperature(); + doublereal RT = GasConstant*T; + fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + + for (iK = 0; iK < m_kk; iK++) { + /* + * We never sample the end of the mole fraction domains + */ + xx = fmaxx(moleFractions_[iK], xxSmall); + /* + * First wipe out the ideal solution mixing term + */ + lnActCoeff_Scaled_[iK] = - log(xx); + + /* + * Then add in the Margules interaction terms. that's it! + */ + for (int i = 0; i < numBinaryInteractions_; i++) { + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + delAK = 0; + delBK = 0; + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + lnActCoeff_Scaled_[iK] += (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1; + } + } + } + //=================================================================================================================== + // 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 + * + * he = X_A X_B(B + C X_B) + * + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::s_update_dlnActCoeff_dT() const { + int iA, iB, iK, delAK, delBK; + doublereal XA, XB, g0, g1; + doublereal T = temperature(); + doublereal RTT = GasConstant*T*T; + fvo_zero_dbl_1(dlnActCoeffdT_Scaled_, m_kk); + fvo_zero_dbl_1(d2lnActCoeffdT2_Scaled_, m_kk); + for (iK = 0; iK < m_kk; iK++) { + for (int i = 0; i < numBinaryInteractions_; i++) { + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + delAK = 0; + delBK = 0; + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + g0 = -m_HE_b_ij[i] / RTT; + g1 = -m_HE_c_ij[i] / RTT; + double temp = (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1; + dlnActCoeffdT_Scaled_[iK] += temp; + d2lnActCoeffdT2_Scaled_[iK] -= 2.0 * temp / T; + } + } + } + //==================================================================================================================== + // + /* + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::getdlnActCoeffdT(doublereal *dlnActCoeffdT) const { + s_update_dlnActCoeff_dT(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdT[k] = dlnActCoeffdT_Scaled_[k]; + } + } + //==================================================================================================================== + // + /* + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::getd2lnActCoeffdT2(doublereal *d2lnActCoeffdT2) const { + s_update_dlnActCoeff_dT(); + for (int k = 0; k < m_kk; k++) { + d2lnActCoeffdT2[k] = d2lnActCoeffdT2_Scaled_[k]; + } + } + //==================================================================================================================== + + // 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. + * + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::getdlnActCoeffds(const doublereal dTds, const doublereal * const dXds, + doublereal *dlnActCoeffds) const { + + + int iA, iB, iK, delAK, delBK; + doublereal XA, XB, XK, g0 , g1, dXA, dXB; + doublereal T = temperature(); + doublereal RT = GasConstant*T; + doublereal xx; + + //fvo_zero_dbl_1(dlnActCoeff, m_kk); + s_update_dlnActCoeff_dT(); + + for (iK = 0; iK < m_kk; iK++) { + + XK = moleFractions_[iK]; + + /* + * We never sample the end of the mole fraction domains + */ + xx = fmaxx(moleFractions_[iK], xxSmall); + /* + * First wipe out the ideal solution mixing term + */ + if (xx > xxSmall) { + dlnActCoeffds[iK] += - 1.0 / xx; + } + + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + delAK = 0; + delBK = 0; + + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + dXA = dXds[iA]; + dXB = dXds[iB]; + + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + + dlnActCoeffds[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB + + dlnActCoeffdT_Scaled_[iK]*dTds; + } + } + } + //==================================================================================================================== + // Update the derivative of the log of the activity coefficients wrt the log of the corresponding species number density + /* + * This function will be called to update the internally stored gradients of the + * logarithm of the activity coefficients. These are used in the determination + * of the diffusion coefficients. + * + * he = X_A X_B(B + C X_B) + * + * This function only carries out the diagonal calculation + * + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN_diag() const { + int iA, iB, iK, delAK, delBK; + doublereal XA, XB, XK, g0 , g1; + doublereal T = temperature(); + doublereal RT = GasConstant*T; + doublereal xx; + + fvo_zero_dbl_1(dlnActCoeffdlnN_diag_, m_kk); + + for (iK = 0; iK < m_kk; iK++) { + + XK = moleFractions_[iK]; + /* + * We never sample the end of the mole fraction domains + */ + xx = fmaxx(moleFractions_[iK], xxSmall); + /* + * First wipe out the ideal solution mixing term + */ + // lnActCoeff_Scaled_[iK] = - log(xx); + if (xx > xxSmall) { + dlnActCoeffdlnN_diag_[iK] = - 1.0 + xx; + } + + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + delAK = 0; + delBK = 0; + + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + + dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB))); + } + dlnActCoeffdlnN_diag_[iK] = XK*dlnActCoeffdlnN_diag_[iK]; + } + + } + //==================================================================================================================== + // Update the derivative of the log of the activity coefficients wrt ln N_k + /* + * This function will be called to update the internally storred gradients of the + * logarithm of the activity coefficients. These are used in the determination + * of the diffusion coefficients. + * + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN() const { + int iA, iB; + doublereal delAK, delBK; + double XA, XB, g0 , g1, XK, XM; + double xx , delKM; + double T = temperature(); + double RT = GasConstant*T; + + doublereal delAM, delBM; + + dlnActCoeffdlnN_.zero(); + + /* + * Loop over the activity coefficient gamma_k + */ + for (int iK = 0; iK < m_kk; iK++) { + XK = moleFractions_[iK]; + /* + * We never sample the end of the mole fraction domains + */ + xx = fmaxx(moleFractions_[iK], xxSmall); + + for (int iM = 0; iM < m_kk; iM++) { + XM = moleFractions_[iM]; + + if (xx > xxSmall) { + delKM = 0.0; + if (iK == iM) delKM = 1.0; + // this gets multiplied by XM at the bottom + dlnActCoeffdlnN_(iK,iM) += - delKM/XM + 1.0; + } + + + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + delAK = 0.0; + delBK = 0.0; + delAM = 0.0; + delBM = 0.0; + if (iA==iK) delAK = 1.0; + else if (iB==iK) delBK = 1.0; + if (iA==iM) delAM = 1.0; + else if (iB==iM) delBM = 1.0; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + + dlnActCoeffdlnN_(iK,iM) += g0*((delAM-XA)*(delBK-XB)+(delAK-XA)*(delBM-XB)); + dlnActCoeffdlnN_(iK,iM) += 2*g1*((delAM-XA)*(delBK-XB)*XB+(delAK-XA)*(delBM-XB)*XB+(delBM-XB)*(delBK-XB)*XA); + + } + dlnActCoeffdlnN_(iK,iM) = XM * dlnActCoeffdlnN_(iK,iM); + } + } + } + //==================================================================================================================== + void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnX_diag() const { + + int iA, iB; + doublereal XA, XB, g0 , g1; + doublereal T = temperature(); + + fvo_zero_dbl_1(dlnActCoeffdlnX_diag_, m_kk); + + doublereal RT = GasConstant * T; + + + 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]; + + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + + dlnActCoeffdlnX_diag_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB); + dlnActCoeffdlnX_diag_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB); + } + throw CanteraError("", "unimplemented"); + } + + //==================================================================================================================== + // + /* + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const { + s_update_dlnActCoeff_dlnN_diag(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k]; + } + } + //==================================================================================================================== + // + /* + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::getdlnActCoeffdlnX_diag(doublereal *dlnActCoeffdlnX_diag) const { + s_update_dlnActCoeff_dlnX_diag(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k]; + } + } + //==================================================================================================================== + // + /* + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::getdlnActCoeffdlnN(const int ld, doublereal *dlnActCoeffdlnN) const { + s_update_dlnActCoeff_dlnN(); + 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]; + } + } + } + //==================================================================================================================== + // + /* + * HKM - Checked for Transition + */ + void PhaseCombo_Interaction::resizeNumInteractions(const int num) { + numBinaryInteractions_ = num; + m_HE_b_ij.resize(num, 0.0); + m_HE_c_ij.resize(num, 0.0); + m_HE_d_ij.resize(num, 0.0); + m_SE_b_ij.resize(num, 0.0); + m_SE_c_ij.resize(num, 0.0); + m_SE_d_ij.resize(num, 0.0); + m_VHE_b_ij.resize(num, 0.0); + m_VHE_c_ij.resize(num, 0.0); + m_VHE_d_ij.resize(num, 0.0); + m_VSE_b_ij.resize(num, 0.0); + m_VSE_c_ij.resize(num, 0.0); + m_VSE_d_ij.resize(num, 0.0); + + m_pSpecies_A_ij.resize(num, -1); + m_pSpecies_B_ij.resize(num, -1); + throw CanteraError("", "unimplemented"); + } + //==================================================================================================================== + + /* + * Process an XML node called "binaryNeutralSpeciesParameters" + * This node contains all of the parameters necessary to describe + * the Margules Interaction for a single binary interaction + * This function reads the XML file and writes the coefficients + * it finds to an internal data structures. + */ + void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node &xmLBinarySpecies) { + string xname = xmLBinarySpecies.name(); + if (xname != "binaryNeutralSpeciesParameters") { + throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies", + "Incorrect name for processing this routine: " + xname); + } + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + int nParamsFound; + vector_fp vParams; + string iName = xmLBinarySpecies.attrib("speciesA"); + if (iName == "") { + throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies", "no speciesA attrib"); + } + string jName = xmLBinarySpecies.attrib("speciesB"); + if (jName == "") { + throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies", "no speciesB attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species + */ + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + string ispName = speciesName(iSpecies); + if (charge[iSpecies] != 0) { + throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies", "speciesA charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + string jspName = speciesName(jSpecies); + if (charge[jSpecies] != 0) { + throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies", "speciesB charge problem"); + } + + resizeNumInteractions(numBinaryInteractions_ + 1); + int iSpot = numBinaryInteractions_ - 1; + 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 + */ + getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessEnthalpy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_HE_b_ij[iSpot] = vParams[0]; + m_HE_c_ij[iSpot] = vParams[1]; + } + + if (nodeName == "excessentropy") { + /* + * Get the string containing all of the values + */ + getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessEntropy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_SE_b_ij[iSpot] = vParams[0]; + m_SE_c_ij[iSpot] = vParams[1]; + } + + if (nodeName == "excessvolume_enthalpy") { + /* + * Get the string containing all of the values + */ + getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessVolume_Enthalpy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_VHE_b_ij[iSpot] = vParams[0]; + m_VHE_c_ij[iSpot] = vParams[1]; + } + + if (nodeName == "excessvolume_entropy") { + /* + * Get the string containing all of the values + */ + getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessVolume_Entropy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_VSE_b_ij[iSpot] = vParams[0]; + m_VSE_c_ij[iSpot] = vParams[1]; + } + + + } + } + //==================================================================================================================== +} +//====================================================================================================================== diff --git a/Cantera/src/thermo/PhaseCombo_Interaction.h b/Cantera/src/thermo/PhaseCombo_Interaction.h new file mode 100644 index 000000000..c5dacb318 --- /dev/null +++ b/Cantera/src/thermo/PhaseCombo_Interaction.h @@ -0,0 +1,705 @@ +/** + * @file + */ +/* + * 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 641 2010-11-12 21:37:41Z hkmoffa $ + */ + +#ifndef CT_LIFES_X_ONEPHASE_VPSSTP_H +#define CT_LIFES_X_ONEPHASE_VPSSTP_H + +#include "PseudoBinaryVPSSTP.h" +#include "GibbsExcessVPSSTP.h" + +namespace Cantera { + + /** + * @ingroup thermoprops + */ + + + //! MargulesVPSSTP is a derived class of GibbsExcessVPSSTP that employs + //! the Margules approximation for the excess gibbs free energy + + class PhaseCombo_Interaction : public GibbsExcessVPSSTP { + + public: + + //! Constructor + /*! + * This doesn't do much more than initialize constants with + * default values for water at 25C. Water molecular weight + * comes from the default elements.xml file. It actually + * differs slightly from the IAPWS95 value of 18.015268. However, + * density conservation and therefore element conservation + * is the more important principle to follow. + */ + PhaseCombo_Interaction(); + + //! Construct and initialize a PhaseCombo_Interaction 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. + */ + PhaseCombo_Interaction(std::string inputFile, std::string id = ""); + + //! Construct and initialize a PhaseCombo_Interaction 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) + */ + PhaseCombo_Interaction(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. + */ + PhaseCombo_Interaction(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 + */ + PhaseCombo_Interaction(const PhaseCombo_Interaction& b); + + //! Assignment operator + /*! + * + * @param b class to be copied. + */ + PhaseCombo_Interaction& operator=(const PhaseCombo_Interaction &b); + + //! Destructor + virtual ~PhaseCombo_Interaction(); + + //! 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 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 log 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) const; + + //@} + + private: + + //! Process an XML node called "binaryNeutralSpeciesParameters" + /*! + * This node contains all of the parameters necessary to describe + * the Margules 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 Margules interaction terms + /*! + * @param num Number of binary Margules 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; + + //! Update the derivative of the log of the activity coefficients + //! wrt log(mole fraction) + /*! + * This function will be called to update the internally storred + * derivative of the natural logarithm of the activity coefficients + * wrt logarithm of the mole fractions. + */ + void s_update_dlnActCoeff_dlnX_diag() const; + + //! Update the derivative of the log of the activity coefficients + //! wrt log(moles) - diagonal only + /*! + * This function will be called to update the internally storred diagonal entries for the + * derivative of the natural logarithm of the activity coefficients + * wrt logarithm of the moles. + */ + void s_update_dlnActCoeff_dlnN_diag() const; + + //! Update the derivative of the log of the activity coefficients wrt log(moles_m) + /*! + * This function will be called to update the internally storred + * derivative of the natural logarithm of the activity coefficients + * wrt logarithm of the mole number of species + */ + void s_update_dlnActCoeff_dlnN() const; + + + 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_; + + //! Enthalpy term for the binary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_HE_b_ij; + + //! Enthalpy term for the ternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_HE_c_ij; + + //! Enthalpy term for the quaternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_HE_d_ij; + + //! Entropy term for the binary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_SE_b_ij; + + //! Entropy term for the ternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_SE_c_ij; + + //! Entropy term for the quaternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_SE_d_ij; + + //! Enthalpy term for the binary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VHE_b_ij; + + //! Enthalpy term for the ternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VHE_c_ij; + + //! Enthalpy term for the quaternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VHE_d_ij; + + //! Entropy term for the binary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VSE_b_ij; + + //! Entropy term for the ternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VSE_c_ij; + + //! Entropy term for the quaternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VSE_d_ij; + + + + //! vector of species indices representing species A in the interaction + /*! + * Each Margules 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 Margules excess Gibbs free energy term involves two species, A and B. + * This vector identifies species B. + */ + vector_int m_pSpecies_B_ij; + + //! form of the Margules interaction expression + /*! + * Currently there is only one form. + */ + int formMargules_; + + //! form of the temperatuer dependence of the Margules interaction expression + /*! + * Currently there is only one form -> constant wrt temperature. + */ + int formTempModel_; + + + }; + + + +} + +#endif + + + + + diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index d6d2dc451..18454b841 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -29,6 +29,7 @@ #include "IdealSolidSolnPhase.h" #include "MargulesVPSSTP.h" #include "IonsFromNeutralVPSSTP.h" +#include "PhaseCombo_Interaction.h" #endif #ifdef WITH_PURE_FLUIDS @@ -99,7 +100,7 @@ namespace Cantera { "PureFluid", "LatticeSolid", "Lattice", "HMW", "IdealSolidSolution", "DebyeHuckel", "IdealMolalSolution", "IdealGasVPSS", - "MineralEQ3", "MetalSHEelectrons", "Margules", + "MineralEQ3", "MetalSHEelectrons", "Margules", "PhaseCombo_Interaction", "IonsFromNeutralMolecule", "FixedChemPot" }; @@ -110,7 +111,7 @@ namespace Cantera { cHMW, cIdealSolidSolnPhase, cDebyeHuckel, cIdealMolalSoln, cVPSS_IdealGas, cMineralEQ3, cMetalSHEelectrons, - cMargulesVPSSTP, cIonsFromNeutral, cFixedChemPot + cMargulesVPSSTP, cPhaseCombo_Interaction, cIonsFromNeutral, cFixedChemPot }; /* @@ -152,6 +153,10 @@ namespace Cantera { th = new MargulesVPSSTP(); break; + case cPhaseCombo_Interaction: + th = new PhaseCombo_Interaction(); + break; + case cIonsFromNeutral: th = new IonsFromNeutralVPSSTP(); break; diff --git a/Cantera/src/thermo/mix_defs.h b/Cantera/src/thermo/mix_defs.h index b104d12dc..b81981de7 100644 --- a/Cantera/src/thermo/mix_defs.h +++ b/Cantera/src/thermo/mix_defs.h @@ -75,6 +75,7 @@ namespace Cantera { const int cMargulesVPSSTP = 301; + const int cPhaseCombo_Interaction = 305; const int cIonsFromNeutral = 2000;