diff --git a/Cantera/src/thermo/IdealMolalSoln.cpp b/Cantera/src/thermo/IdealMolalSoln.cpp new file mode 100644 index 000000000..2e28991ec --- /dev/null +++ b/Cantera/src/thermo/IdealMolalSoln.cpp @@ -0,0 +1,1026 @@ +/** + * + * @file IdealMolalSoln.cpp + */ +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Author$ + * $Date$ + * $Revision$ + */ + + +#ifndef MAX +#define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif + +#include "IdealMolalSoln.h" +#include "importCTML.h" + +namespace Cantera { + + /** + * Default constructor + */ + IdealMolalSoln::IdealMolalSoln() : + MolalityVPSSTP(), + m_Pcurrent(OneAtm), + m_formGC(2) + { + } + + + /** + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + IdealMolalSoln::IdealMolalSoln(const IdealMolalSoln &b) : + MolalityVPSSTP(b) + { + /* + * Use the assignment operator to do the brunt + * of the work for the copy construtor. + */ + *this = b; + } + + /** + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + IdealMolalSoln& IdealMolalSoln:: + operator=(const IdealMolalSoln &b) { + if (&b != this) { + MolalityVPSSTP::operator=(b); + m_speciesMolarVolume = b.m_speciesMolarVolume; + m_Pcurrent = b.m_Pcurrent; + m_formGC = b.m_formGC; + m_expg0_RT = b.m_expg0_RT; + m_pe = b.m_pe; + m_pp = b.m_pp; + m_tmpV = b.m_tmpV; + } + return *this; + } + + IdealMolalSoln::IdealMolalSoln(string inputFile, string id) : + MolalityVPSSTP() + { + constructPhaseFile(inputFile, id); + } + + IdealMolalSoln::IdealMolalSoln(XML_Node& root, string id) : + MolalityVPSSTP() + { + constructPhaseXML(root, id); + } + + /** + * + * ~IdealMolalSoln(): (virtual) + * + * Destructor: does nothing: + * + */ + IdealMolalSoln::~IdealMolalSoln() { + } + + /** + * + */ + ThermoPhase* IdealMolalSoln::duplMyselfAsThermoPhase() { + IdealMolalSoln* mtp = new IdealMolalSoln(*this); + return (ThermoPhase *) mtp; + } + + // + // -------- Molar Thermodynamic Properties of the Solution --------------- + // + /** + * Molar enthalpy of the solution. Units: J/kmol. + */ + doublereal IdealMolalSoln::enthalpy_mole() const { + getPartialMolarEnthalpies(DATA_PTR(m_tmpV)); + getMoleFractions(DATA_PTR(m_pp)); + double val = mean_X(DATA_PTR(m_tmpV)); + return val; + } + + /** + * Molar internal energy of the solution. Units: J/kmol. + */ + doublereal IdealMolalSoln::intEnergy_mole() const { + getPartialMolarEnthalpies(DATA_PTR(m_tmpV)); + return mean_X(DATA_PTR(m_tmpV)); + } + + doublereal IdealMolalSoln::entropy_mole() const { + getPartialMolarEntropies(DATA_PTR(m_tmpV)); + return mean_X(DATA_PTR(m_tmpV)); + } + + /// Molar Gibbs function. Units: J/kmol. + doublereal IdealMolalSoln::gibbs_mole() const { + getChemPotentials(DATA_PTR(m_tmpV)); + return mean_X(DATA_PTR(m_tmpV)); + } + + /// Molar heat capacity at constant pressure. Units: J/kmol/K. + doublereal IdealMolalSoln::cp_mole() const { + getPartialMolarCp(DATA_PTR(m_tmpV)); + double val = mean_X(DATA_PTR(m_tmpV)); + return val; + } + + + /// Molar heat capacity at constant volume. Units: J/kmol/K. + doublereal IdealMolalSoln::cv_mole() const { + return err("not implemented"); + } + + // + // ------- Mechanical Equation of State Properties ------------------------ + // + + /** + * Pressure. Units: Pa. + * For this incompressible system, we return the internally storred + * independent value of the pressure. + */ + doublereal IdealMolalSoln::pressure() const { + return m_Pcurrent; + } + + /** + * The isothermal compressibility. Units: 1/Pa. + * The isothermal compressibility is defined as + * \f[ + * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T + * \f] + * + * It's equal to zero for this model, since the molar volume + * doesn't change with pressure or temperature. + */ + doublereal IdealMolalSoln::isothermalCompressibility() const { + return 0.0; + } + + /** + * The thermal expansion coefficient. Units: 1/K. + * The thermal expansion coefficient is defined as + * + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + * + * It's equal to zero for this model, since the molar volume + * doesn't change with pressure or temperature. + */ + doublereal IdealMolalSoln::thermalExpansionCoeff() const { + return 0.0; + } + + /** + * Overwritten setDensity() function is necessary because the + * density is not an indendent variable. + * + * This function will now throw an error condition + * + * @internal May have to adjust the strategy here to make + * the eos for these materials slightly compressible, in order + * to create a condition where the density is a function of + * the pressure. + * + * This function will now throw an error condition. + * + * NOTE: This is an overwritten function from the State.h + * class + */ + void IdealMolalSoln::setDensity(doublereal rho) { + double dens = density(); + if (rho != dens) { + throw CanteraError("Idea;MolalSoln::setDensity", + "Density is not an independent variable"); + } + } + + /** + * Overwritten setMolarDensity() function is necessary because the + * density is not an indendent variable. + * + * This function will now throw an error condition. + * + * NOTE: This is a virtual function, overwritten function from the State.h + * class + */ + void IdealMolalSoln::setMolarDensity(doublereal conc) { + double concI = State::molarDensity(); + if (conc != concI) { + throw CanteraError("IdealMolalSoln::setMolarDensity", + "molarDensity/denisty is not an independent variable"); + } + } + + // + // ------- Activities and Activity Concentrations + // + + /** + * This method returns an array of generalized concentrations + * \f$ C_k\f$ that are defined such that + * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ + * is a standard concentration + * defined below. These generalized concentrations are used + * by kinetics manager classes to compute the forward and + * reverse rates of elementary reactions. + * + * @param c Array of generalized concentrations. The + * units depend upon the implementation of the + * reaction rate expressions within the phase. + */ + void IdealMolalSoln::getActivityConcentrations(doublereal* c) const { + if (m_formGC != 1) { + double c_solvent = standardConcentration(); + getActivities(c); + for (int k = 0; k < m_kk; k++) { + c[k] *= c_solvent; + } + } else { + getActivities(c); + for (int k = 0; k < m_kk; k++) { + double c0 = standardConcentration(k); + c[k] *= c0; + } + } + } + + /** + * 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. + * + * For the time being we will use the concentration of pure + * solvent for the the standard concentration of all species. + * This has the effect of making reaction rates + * based on the molality of species proportional to the + * molality of the species. + */ + doublereal IdealMolalSoln::standardConcentration(int k) const { + double c0, mvSolvent; + switch (m_formGC) { + case 0: + c0 = 1.0; + break; + case 1: + c0 = 1.0 /m_speciesMolarVolume[m_indexSolvent]; + break; + case 2: + mvSolvent = m_speciesMolarVolume[m_indexSolvent]; + c0 = 1.0 / mvSolvent; + break; + } + return c0; + } + + /** + * Returns the natural logarithm of the standard + * concentration of the kth species + */ + doublereal IdealMolalSoln::logStandardConc(int k) const { + double c0 = standardConcentration(k); + return log(c0); + } + + /** + * Returns the units of the standard and general concentrations + * Note they have the same units, as their divisor is + * defined to be equal to the activity of the kth species + * in the solution, which is unitless. + * + * This routine is used in print out applications where the + * units are needed. Usually, MKS units are assumed throughout + * the program and in the XML input files. + * + * On return uA contains the powers of the units (MKS assumed) + * of the standard concentrations and generalized concentrations + * for the kth species. + * + * uA[0] = kmol units - default = 1 + * uA[1] = m units - default = -nDim(), the number of spatial + * dimensions in the Phase class. + * uA[2] = kg units - default = 0; + * uA[3] = Pa(pressure) units - default = 0; + * uA[4] = Temperature units - default = 0; + * uA[5] = time units - default = 0 + */ + void IdealMolalSoln::getUnitsStandardConc(double *uA, int k, int sizeUA) { + int eos = eosType(); + if (eos == 0) { + for (int i = 0; i < sizeUA; i++) { + uA[i] = 0.0; + } + } else { + for (int i = 0; i < sizeUA; i++) { + if (i == 0) uA[0] = 1.0; + if (i == 1) uA[1] = -nDim(); + if (i == 2) uA[2] = 0.0; + if (i == 3) uA[3] = 0.0; + if (i == 4) uA[4] = 0.0; + if (i == 5) uA[5] = 0.0; + } + } + } + + /** + * Get the array of non-dimensional molality-based + * activities at the current solution temperature, + * pressure, and solution concentration. + * + * The max against 8.689E-3 is to limit the activity + * coefficient to be greater than 1.0E-50. + */ + void IdealMolalSoln::getActivities(doublereal* ac) const { + /* + * Update the molality array, m_molalities() + * This requires an update due to mole fractions + */ + calcMolalities(); + for (int k = 0; k < m_kk; k++) { + ac[k] = m_molalities[k]; + } + double xmolSolvent = moleFraction(m_indexSolvent); + xmolSolvent = MAX(8.689E-3, xmolSolvent); + ac[m_indexSolvent] = + exp((xmolSolvent - 1.0)/xmolSolvent); + } + + /** + * Get the array of non-dimensional Molality based + * activity coefficients at + * the current solution temperature, pressure, and + * solution concentration. + * See Denbigh + * (note solvent activity coefficient is on the molar scale). + * + * The max against 5.0E-3 (1/200) is to limit the activity + * coefficient to be greater than 1.0E-50. + */ + void IdealMolalSoln:: + getMolalityActivityCoefficients(doublereal* acMolality) const { + for (int k = 0; k < m_kk; k++) { + acMolality[k] = 1.0; + } + double xmolSolvent = moleFraction(m_indexSolvent); + xmolSolvent = MAX(8.689E-3, xmolSolvent); + acMolality[m_indexSolvent] = + exp((xmolSolvent - 1.0)/xmolSolvent) / xmolSolvent; + } + + // + // ------ 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. + * + * \f[ + * \mu_k = \mu^{o}_k(T,P) + R T ln(m_k / m_units) + * \f] + * + * where m_units is equal to 1 gmol kg-1 + * + * \f[ + * \mu_solvent = \mu^{o}_solvent(T,P) + + * R T ((X_solvent - 1.0) / X_solvent) + * \f] + */ + void IdealMolalSoln::getChemPotentials(doublereal* mu) const{ + double xx; + const double xxSmall = 1.0E-150; + /* + * First get the standard chemical potentials + * -> this requires updates of standard state as a function + * of T and P + * These are defined at unit molality. + */ + getStandardChemPotentials(mu); + /* + * Update the molality array, m_molalities() + * This requires an update due to mole fractions + */ + calcMolalities(); + /* + * + */ + doublereal RT = GasConstant * temperature(); + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + xx = MAX(m_molalities[k], xxSmall); + mu[k] += RT * log(xx); + } + } + /* + * Do the solvent + * -> see my notes + */ + double xmolSolvent = moleFraction(m_indexSolvent); + xx = MAX(xmolSolvent, xxSmall); + mu[m_indexSolvent] += + (RT * (xmolSolvent - 1.0) / xx); + } + + /** + * 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 + * SS species enthalpies + * \f[ + * \bar h_k(T,P) = \hat h^{0}_k(T,P) + * \f] + * + * note hbar = ubar + T d(ubar/dT) + * see note about partial molar entropies. + */ + void IdealMolalSoln::getPartialMolarEnthalpies(doublereal* hbar) const { + getEnthalpy_RT(hbar); + doublereal RT = _RT(); + for (int k = 0; k < m_kk; k++) { + hbar[k] *= RT; + } + } + + /** + * + * getPartialMolarEntropies() (virtual, const) + * + * Returns an array of partial molar entropies of the species in the + * solution. Units: J/kmol. + * + * Maxwell's equations provide an insight in how to calculate this + * (p.215 Smith and Van Ness) + * + * d(chemPot_i)/dT = -sbar_i + * + * + * For this phase, the partial molar entropies are equal to the + * SS species entropies plus the ideal solution contribution. + * + * \f[ + * \bar s_k(T,P) = \hat s^0_k(T) - R log( molality[k] / m_units) + * \f] + * \f[ + * \bar s_solvent(T,P) = \hat s^0_solvent(T) + * - R ((xmolSolvent - 1.0) / xmolSolvent) + * \f] + * + * The reference-state pure-species entropies,\f$ \hat s^0_k(T) \f$, + * at the reference pressure, \f$ P_{ref} \f$, are computed by the + * species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + */ + void IdealMolalSoln:: + getPartialMolarEntropies(doublereal* sbar) const { + getEntropy_R(sbar); + doublereal R = GasConstant; + doublereal mm; + calcMolalities(); + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + mm = fmaxx(SmallNumber, m_molalities[k]); + sbar[k] -= R * log(mm); + } + } + double xmolSolvent = moleFraction(m_indexSolvent); + sbar[m_indexSolvent] -= (R * (xmolSolvent - 1.0) / xmolSolvent); + } + + /** + * returns an array of partial molar volumes of the species + * in the solution. Units: m^3 kmol-1. + * + * For this solution, the partial molar volumes are equal to the + * constant species molar volumes. + */ + void IdealMolalSoln::getPartialMolarVolumes(doublereal* vbar) const { + getStandardVolumes(vbar); + } + + /* + * Partial molar heat capacity of the solution: + * The kth partial molar heat capacity is equal to + * the temperature derivative of the partial molar + * enthalpy of the kth species in the solution at constant + * P and composition (p. 220 Smith and Van Ness). + * + * Cp = -T d2(chemPot_i)/dT2 + */ + void IdealMolalSoln::getPartialMolarCp(doublereal* cpbar) const { + /* + * Get the nondimensional gibbs standard state of the + * species at the T and P of the solution. + */ + getCp_R(cpbar); + + for (int k = 0; k < m_kk; k++) { + cpbar[k] *= GasConstant; + } + } + + /* + * -------- Properties of the Standard State of the Species + * in the Solution ------------------ + */ + + /** + * getStandardChemPotentials() (virtual, const) + * + * + * Get the standard state chemical potentials of the species. + * This is the array of chemical potentials at unit activity + * (Mole fraction scale) + * \f$ \mu^0_k(T,P) \f$. + * We define these here as the chemical potentials of the pure + * species at the temperature and pressure of the solution. + * This function is used in the evaluation of the + * equilibrium constant Kc. Therefore, Kc will also depend + * on T and P. This is the norm for liquid and solid systems. + * + * units = J / kmol + */ + void IdealMolalSoln::getStandardChemPotentials(doublereal* mu) const { + getGibbs_ref(mu); + doublereal pref; + doublereal delta_p; + for (int k = 0; k < m_kk; k++) { + pref = m_spthermo->refPressure(k); + delta_p = m_Pcurrent - pref; + mu[k] += delta_p * m_speciesMolarVolume[k]; + } + } + + /** + * Get the nondimensional gibbs function for the species + * standard states at the current T and P of the solution. + * + * \f[ + * \mu^0_k(T,P) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure + * species k at the reference pressure, \f$P_{ref}\f$. + * + * @param grt Vector of length m_kk, which on return sr[k] + * will contain the nondimensional + * standard state gibbs function for species k. + */ + void IdealMolalSoln::getGibbs_RT(doublereal* grt) const { + getPureGibbs(grt); + doublereal invRT = 1.0 / _RT(); + for (int k = 0; k < m_kk; k++) { + grt[k] *= invRT; + } + } + + /** + * + * getPureGibbs() + * + * Get the Gibbs functions for the pure species + * at the current T and P of the solution. + * We assume an incompressible constant partial molar + * volume here: + * \f[ + * \mu^0_k(T,p) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k + * \f] + * where \f$V_k\f$ is the molar volume of pure species k<\I>. + * \f$ u^{ref}_k(T)\f$ is the chemical potential of pure + * species k<\I> at the reference pressure, \f$P_{ref}\f$. + */ + void IdealMolalSoln::getPureGibbs(doublereal* gpure) const { + getGibbs_ref(gpure); + doublereal pref; + doublereal delta_p; + for (int k = 0; k < m_kk; k++) { + pref = m_spthermo->refPressure(k); + delta_p = m_Pcurrent - pref; + gpure[k] += delta_p * m_speciesMolarVolume[k]; + } + } + + /** + * + * getEnthalpy_RT() (virtual, const) + * + * Get the array of nondimensional Enthalpy functions for the ss + * species at the current T and P of the solution. + * We assume an incompressible constant partial molar + * volume here: + * \f[ + * h^0_k(T,P) = h^{ref}_k(T) + (P - P_{ref}) * V_k + * \f] + * where \f$V_k\f$ is the molar volume of SS species k<\I>. + * \f$ h^{ref}_k(T)\f$ is the enthalpy of the SS + * species k<\I> at the reference pressure, \f$P_{ref}\f$. + */ + void IdealMolalSoln:: + getEnthalpy_RT(doublereal* hrt) const { + getEnthalpy_RT_ref(hrt); + doublereal pref; + doublereal delta_p; + double RT = _RT(); + for (int k = 0; k < m_kk; k++) { + pref = m_spthermo->refPressure(k); + delta_p = m_Pcurrent - pref; + hrt[k] += delta_p/ RT * m_speciesMolarVolume[k]; + } + } + + /** + * getEntropy_R() (virtual, const) + * + * Get the nondimensional Entropies for the species + * standard states at the current T and P of the solution. + * + * Note, this is equal to the reference state entropies + * due to the zero volume expansivity: + * i.e., (dS/dp)_T = (dV/dT)_P = 0.0 + * + * @param sr Vector of length m_kk, which on return sr[k] + * will contain the nondimensional + * standard state entropy of species k. + */ + void IdealMolalSoln:: + getEntropy_R(doublereal* sr) const { + getEntropy_R_ref(sr); + } + + /** + * Get the nondimensional heat capacity at constant pressure + * function for the species + * standard states at the current T and P of the solution. + * \f[ + * Cp^0_k(T,P) = Cp^{ref}_k(T) + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ Cp^{ref}_k(T)\f$ is the constant pressure heat capacity + * of species k at the reference pressure, \f$p_{ref}\f$. + * + * @param cpr Vector of length m_kk, which on return cpr[k] + * will contain the nondimensional + * constant pressure heat capacity for species k. + */ + void IdealMolalSoln::getCp_R(doublereal* cpr) const { + getCp_R_ref(cpr); + } + + /** + * Get the molar volumes of each species in their standard + * states at the current + * T and P of the solution. + * units = m^3 / kmol + */ + void IdealMolalSoln::getStandardVolumes(doublereal *vol) const { + copy(m_speciesMolarVolume.begin(), + m_speciesMolarVolume.end(), vol); + } + + + /* + * ------ Thermodynamic Values for the Species Reference States --- + */ + + // -> This is handled by VPStandardStatesTP + + /* + * -------------- Utilities ------------------------------- + */ + + /** + * Initialization of an IdealSolidSolnPhase phase: + * Note this function is pretty much useless because it doesn't + * get the xml tree passed to it. Suggest a change. + */ + void IdealMolalSoln::initThermo() { + initLengths(); + MolalityVPSSTP::initThermo(); + } + + /** + * Initialization of an IdealSolidSolnPhase phase using an + * xml file + * + * This routine is a precursor to initThermo(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 IdealMolalSoln::constructPhaseFile(string inputFile, string id) { + + if (inputFile.size() == 0) { + throw CanteraError("IdealSolidSolnPhase::constructPhaseFile", + "input file is null"); + } + string path = findInputFile(inputFile); + ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("IdealMolalSoln::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("IdealMolalSoln::constructPhaseFile", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + fxml_phase->copy(&phaseNode_XML); + constructPhaseXML(*fxml_phase, id); + delete fxml; + } + + /** + * Import and initialize an IdealSolidSolnPhase 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 IdealMolalSoln::constructPhaseXML(XML_Node& phaseNode, string id) { + if (id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("IdealMolalSoln::constructPhaseXML", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("IdealMolalSoln::constructPhaseXML", + "no thermo XML node"); + } + + /* + * Call the Cantera importPhase() function. This will import + * all of the species into the phase. Then, it will call + * initThermoXML() below. + */ + bool m_ok = importPhase(phaseNode, this); + if (!m_ok) { + throw CanteraError("constructPhaseXML","importPhase failed "); + } + } + + /** + * Import and initialize an IdealSolidSolnPhase 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 IdealMolalSoln::initThermoXML(XML_Node& phaseNode, string id) { + + /* + * Initialize the whole thermo object, using a virtual function. + */ + initThermo(); + + if (id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("initThermo", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("IdealMolalSoln::initThermo", + "no thermo XML node"); + } + XML_Node& thermoNode = phaseNode.child("thermo"); + + /* + * Possible change the form of the standard concentrations + */ + if (thermoNode.hasChild("standardConc")) { + XML_Node& scNode = thermoNode.child("standardConc"); + m_formGC = 2; + string formString = scNode.attrib("model"); + if (formString != "") { + if (formString == "unity") { + m_formGC = 0; + } else if (formString == "molar_volume") { + m_formGC = 1; + } else if (formString == "solvent_volume") { + m_formGC = 2; + } else { + throw CanteraError("IdealMolalSoln::initThermo", + "Unknown standardConc model: " + formString); + } + } + } + + /* + * Get the Name of the Solvent: + * solventName + */ + string solventName = ""; + if (thermoNode.hasChild("solvent")) { + XML_Node& scNode = thermoNode.child("solvent"); + vector nameSolventa; + getStringArray(scNode, nameSolventa); + int nsp = static_cast(nameSolventa.size()); + if (nsp != 1) { + throw CanteraError("IdealMolalSoln::initThermoXML", + "badly formed solvent XML node"); + } + solventName = nameSolventa[0]; + } + + + /* + * Reconcile the solvent name and index. + */ + for (int k = 0; k < m_kk; k++) { + string sname = speciesName(k); + if (solventName == sname) { + m_indexSolvent = k; + break; + } + } + if (m_indexSolvent == -1) { + cout << "IdealMolalSoln::initThermo: Solvent Name not found" + << endl; + throw CanteraError("IdealMolalSoln::initThermo", + "Solvent name not found"); + } + if (m_indexSolvent != 0) { + throw CanteraError("IdealMolalSoln::initThermo", + "Solvent " + solventName + + " should be first species"); + } + + /* + * Now go get the molar volumes + */ + XML_Node& speciesList = phaseNode.child("speciesArray"); + XML_Node* speciesDB = + get_XML_NameID("speciesData", speciesList["datasrc"], + &phaseNode.root()); + const vector&sss = speciesNames(); + + for (int k = 0; k < m_kk; k++) { + XML_Node* s = speciesDB->findByAttr("name", sss[k]); + XML_Node *ss = s->findByName("standardState"); + m_speciesMolarVolume[k] = getFloat(*ss, "molarVolume", "-"); +#ifdef DEBUG_HKM + cout << "species " << sss[k] << " has volume " << + m_speciesMolarVolume[k] << endl; +#endif + } + + /* + * Set the state + */ + if (phaseNode.hasChild("state")) { + XML_Node& stateNode = phaseNode.child("state"); + setStateFromXML(stateNode); + } + + } + + /** + * @internal + * Set equation of state parameters. The number and meaning of + * these depends on the subclass. + * @param n number of parameters + * @param c array of \i n coefficients + * + */ + void IdealMolalSoln::setParameters(int n, doublereal* c) { + } + void IdealMolalSoln::getParameters(int &n, doublereal * const c) { + } + /** + * Set equation of state parameter values from XML + * entries. This method is called by function importPhase in + * file importCTML.cpp when processing a phase definition in + * an input file. It should be overloaded in subclasses to set + * any parameters that are specific to that particular phase + * model. + * + * @param eosdata An XML_Node object corresponding to + * the "thermo" entry for this phase in the input file. + * + * HKM -> Right now, the parameters are set elsewhere (initThermo) + * It just didn't seem to fit. + */ + void IdealMolalSoln::setParametersFromXML(const XML_Node& eosdata) { + } + + /* + * ----------- Critical State Properties -------------------------- + */ + + /* + * ------------ Private and Restricted Functions ------------------ + */ + + /** + * Bail out of functions with an error exit if they are not + * implemented. + */ + doublereal IdealMolalSoln::err(string msg) const { + throw CanteraError("IdealMolalSoln", + "Unfinished func called: " + msg ); + return 0.0; + } + + /** + * This internal function adjusts the lengths of arrays + */ + void IdealMolalSoln::initLengths() { + m_kk = nSpecies(); + /* + * Obtain the limits of the temperature from the species + * thermo handler's limits. + */ + int leng = m_kk; + m_expg0_RT.resize(leng); + m_pe.resize(leng, 0.0); + m_pp.resize(leng); + m_speciesMolarVolume.resize(leng); + m_tmpV.resize(leng); + } + +} + diff --git a/Cantera/src/thermo/IdealMolalSoln.h b/Cantera/src/thermo/IdealMolalSoln.h new file mode 100644 index 000000000..5b6222889 --- /dev/null +++ b/Cantera/src/thermo/IdealMolalSoln.h @@ -0,0 +1,841 @@ +/** + * @file IdealMolalSoln.h + * + * 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 on the molality scale. The Ideal molal + * solution assumes that all molality-based activity + * coefficients are equal to one. + */ + +/* + * $Author$ + * $Date$ + * $Revision$ + */ + +#ifndef CT_IDEALMOLALSOLN_H +#define CT_IDEALMOLALSOLN_H + +#include "MolalityVPSSTP.h" + +namespace Cantera { + + /** + * @defgroup thermoprops Thermodynamic Properties + * + * These classes are used to compute thermodynamic properties. + */ + + class IdealMolalSoln : public MolalityVPSSTP { + + public: + + /// Constructors + IdealMolalSoln(); + IdealMolalSoln(const IdealMolalSoln &); + IdealMolalSoln& operator=(const IdealMolalSoln&); + + IdealMolalSoln(string inputFile, string id = ""); + IdealMolalSoln(XML_Node& phaseRef, string id = ""); + + /// Destructor. + virtual ~IdealMolalSoln(); + + + ThermoPhase *duplMyselfAsThermoPhase(); + + /** + * + * @name Utilities + * @{ + */ + + /** + * Equation of state type flag. The base class returns + * zero. Subclasses should define this to return a unique + * non-zero value. Constants defined for this purpose are + * listed in mix_defs.h. + */ + virtual int eosType() const { return 0; } + + + /** + * @} + * @name Molar Thermodynamic Properties of the Solution --------------- + * @{ + */ + + /// Molar enthalpy. Units: J/kmol. + /** + * Molar enthalpy of the solution. Units: J/kmol. + */ + virtual doublereal enthalpy_mole() const; + + /// Molar internal energy. Units: J/kmol. + /** + * Molar internal energy of the solution. Units: J/kmol. + */ + virtual doublereal intEnergy_mole() const; + + /// Molar entropy. Units: J/kmol/K. + /** + * Molar entropy of the solution. Units: J/kmol/K. + * For an ideal, constant partial molar volume solution mixture with + * pure species phases which exhibit zero volume expansivity: + * \f[ + * \hat s(T, P, X_k) = \sum_k X_k \hat s^0_k(T) + * - \hat R \sum_k X_k log(X_k) + * \f] + * The reference-state pure-species entropies + * \f$ \hat s^0_k(T,p_{ref}) \f$ are computed by the + * species thermodynamic + * property manager. The pure species entropies are independent of + * temperature since the volume expansivities are equal to zero. + * @see SpeciesThermo + */ + virtual doublereal entropy_mole() const; + + /// Molar Gibbs function. Units: J/kmol. + + virtual doublereal gibbs_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; + + //@} + /** @name Mechanical Equation of State Properties ------------------------- + //@{ + * + * In this equation of state implementation, the density is a + * function only of the mole fractions. Therefore, it can't be + * an independent variable. Instead, the pressure is used as the + * independent variable. Functions which try to set the thermodynamic + * state by calling setDensity() may cause an exception to be + * thrown. + */ + + + /** + * Pressure. Units: Pa. + * For this incompressible system, we return the internally storred + * independent value of the pressure. + */ + virtual doublereal pressure() const; + + /** + * Set the pressure at constant temperature. Units: Pa. + * This method sets a constant within the object. + * The mass density is not a function of pressure. + */ + virtual void setPressure(doublereal p) { + m_Pcurrent = p; + } + + + /** + * Calculate the density of the mixture using the partial + * molar volumes and mole fractions as input + * + * The formula for this is + * + * \f[ + * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} + * \f] + * + * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are + * the molecular weights, and \f$V_k\f$ are the pure species + * molar volumes. + * + * Note, the basis behind this formula is that in an ideal + * solution the partial molar volumes are equal to the pure + * species molar volumes. We have additionally specified + * in this class that the pure species molar volumes are + * independent of temperature and pressure. + * + * NOTE: This is a non-virtual function, which is not a + * member of the ThermoPhase base class. + */ + void calcDensity(); + + /** + * Overwritten setDensity() function is necessary because the + * density is not an indendent variable. + * + * This function will now throw an error condition + * + * @internal May have to adjust the strategy here to make + * the eos for these materials slightly compressible, in order + * to create a condition where the density is a function of + * the pressure. + * + * This function will now throw an error condition. + * + * NOTE: This is an overwritten function from the State.h + * class + */ + void setDensity(doublereal rho); + + /** + * Overwritten setMolarDensity() function is necessary because the + * density is not an indendent variable. + * + * This function will now throw an error condition. + * + * NOTE: This is an overwritten function from the State.h + * class + */ + void setMolarDensity(doublereal rho); + + /** + * The isothermal compressibility. Units: 1/Pa. + * The isothermal compressibility is defined as + * \f[ + * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T + * \f] + */ + virtual doublereal isothermalCompressibility() const; + + /** + * The thermal expansion coefficient. Units: 1/K. + * The thermal expansion coefficient is defined as + * + * \f[ + * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P + * \f] + */ + virtual doublereal thermalExpansionCoeff() const; + + /** + * @} + * @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. + * @{ + */ + + /** + * Set the potential energy of species k to pe. + * Units: J/kmol. + * This function must be reimplemented in inherited classes + * of ThermoPhase. + */ + virtual void setPotentialEnergy(int k, doublereal pe) { + err("setPotentialEnergy"); + } + + /** + * Get the potential energy of species k. + * Units: J/kmol. + * This function must be reimplemented in inherited classes + * of ThermoPhase. + */ + virtual doublereal potentialEnergy(int k) const { + return err("potentialEnergy"); + } + + /** + * Set the electric potential of this phase (V). + * This is used by classes InterfaceKinetics and EdgeKinetics to + * compute the rates of charge-transfer reactions, and in computing + * the electrochemical potentials of the species. + */ + void setElectricPotential(doublereal v) { + m_phi = v; + } + + /// The electric potential of this phase (V). + doublereal electricPotential() const { return m_phi; } + + + /** + * @} + * @name Activities 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)\f$ is + * the chemical potential at unit activity, which depends only + * on temperature and the pressure. + * @{ + */ + + /** + * This method returns an array of generalized concentrations + * \f$ C_k\f$ that are defined such that + * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ + * is a standard concentration + * defined below. These generalized concentrations are used + * by kinetics manager classes to compute the forward and + * reverse rates of elementary reactions. + * + * @param c 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. + */ + virtual doublereal standardConcentration(int k=0) const; + + /** + * Returns the natural logarithm of the standard + * concentration of the kth species + */ + virtual doublereal logStandardConc(int k=0) const; + + /** + * Returns the units of the standard and generalized + * concentrations Note they have the same units, as their + * ratio is defined to be equal to the activity of the kth + * species in the solution, which is unitless. + * + * This routine is used in print out applications where the + * units are needed. Usually, MKS units are assumed throughout + * the program and in the XML input files. + * + * uA[0] = kmol units - default = 1 + * uA[1] = m units - default = -nDim(), the number of spatial + * dimensions in the Phase class. + * uA[2] = kg units - default = 0; + * uA[3] = Pa(pressure) units - default = 0; + * uA[4] = Temperature units - default = 0; + * uA[5] = time units - default = 0 + */ + virtual void getUnitsStandardConc(double *uA, int k = 0, + int sizeUA = 6); + + /** + * Get the array of non-dimensional activities at + * the current solution temperature, pressure, and + * solution concentration. + * + * (note solvent is on molar scale) + */ + virtual void getActivities(doublereal* ac) const; + + /** + * Get the array of non-dimensional molality-based + * activity coefficients at the current solution temperature, + * pressure, and solution concentration. + * + * + * (note solvent is on molar scale. The solvent molar + * based activity coefficient is returned). + */ + virtual void + getMolalityActivityCoefficients(doublereal* acMolality) 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. + * \f[ + * \mu_k = \mu^{ref}_k(T) + V_k * (p - p_o) + R T ln(X_k) + * \f] + * or another way to phrase this is + * \f[ + * \mu_k = \mu^o_k(T,p) + R T ln(X_k) + * \f] + * where \f$ \mu^o_k(T,p) = \mu^{ref}_k(T) + V_k * (p - p_o)\f$ + */ + virtual void getChemPotentials(doublereal* mu) 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 + */ + void getElectrochemPotentials(doublereal* mu) const { + getChemPotentials(mu); + double ve = Faraday * electricPotential(); + for (int k = 0; k < m_kk; k++) { + mu[k] += ve*charge(k); + } + } + + /** + * 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 + * pure species enthalpies + * \f[ + * \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k + * \f] + * The reference-state pure-species enthalpies, \f$ \hat h^{ref}_k(T) \f$, + * at the reference pressure,\f$ P_{ref} \f$, + * are computed by the species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + */ + virtual void getPartialMolarEnthalpies(doublereal* hbar) const; + + /** + * getPartialMolarEntropies() (virtual, const) + * + * Returns an array of partial molar entropies of the species in the + * solution. Units: J/kmol. + * + * Maxwell's equations provide an insight in how to calculate this + * (p.215 Smith and Van Ness) + * + * d(chemPot_i)/dT = -sbar_i + * + * + * For this phase, the partial molar entropies are equal to the + * SS species entropies plus the ideal solution contribution.following + * contribution: + * \f[ + * \bar s_k(T,P) = \hat s^0_k(T) - R log(M0 * molality[k]) + * \f] + * \f[ + * \bar s_solvent(T,P) = \hat s^0_solvent(T) + * - R ((xmolSolvent - 1.0) / xmolSolvent) + * \f] + * + * The reference-state pure-species entropies,\f$ \hat s^0_k(T) \f$, + * at the reference pressure, \f$ P_{ref} \f$, are computed by the + * species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + */ + virtual void getPartialMolarEntropies(doublereal* sbar) const; + + /** + * returns an array of partial molar volumes of the species + * in the solution. Units: m^3 kmol-1. + * + * For this solution, thepartial molar volumes are equal to the + * constant species molar volumes. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + + /* + * Partial molar heat capacity of the solution: + * The kth partial molar heat capacity is equal to + * the temperature derivative of the partial molar + * enthalpy of the kth species in the solution at constant + * P and composition (p. 220 Smith and Van Ness). + * + * Cp = -T d2(chemPot_i)/dT2 + */ + virtual void getPartialMolarCp(doublereal* cpbar) const; + + //@} + /// @name Properties of the Standard State of the Species + // in the Solution -- + //@{ + + /** + * Get the standard state chemical potentials of the species. + * This is the array of chemical potentials at unit activity + * \f$ \mu^0_k(T,P) \f$. + * We define these here as the chemical potentials of the pure + * species at the temperature and pressure of the solution. + * This function is used in the evaluation of the + * equilibrium constant Kc. Therefore, Kc will also depend + * on T and P. This is the norm for liquid and solid systems. + * + * units = J / kmol + */ + virtual void getStandardChemPotentials(doublereal* mu) const; + + /** + * Get the nondimensional gibbs function for the species + * standard states at the current T and P of the solution. + * + * \f[ + * \mu^0_k(T,P) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure + * species k at the reference pressure, \f$P_{ref}\f$. + * + * @param grt Vector of length m_kk, which on return sr[k] + * will contain the nondimensional + * standard state gibbs function for species k. + */ + virtual void getGibbs_RT(doublereal* grt) const; + + /** + * Get the nondimensional Gibbs functions for the standard + * state of the species at the current T and P. + */ + virtual void getPureGibbs(doublereal* gpure) const; + + /** + * + * getEnthalpy_RT() (virtual, const) + * + * Get the array of nondimensional Enthalpy functions for the ss + * species at the current T and P of the solution. + * We assume an incompressible constant partial molar + * volume here: + * \f[ + * h^0_k(T,P) = h^{ref}_k(T) + (P - P_{ref}) * V_k + * \f] + * where \f$V_k\f$ is the molar volume of SS species k<\I>. + * \f$ h^{ref}_k(T)\f$ is the enthalpy of the SS + * species k<\I> at the reference pressure, \f$P_{ref}\f$. + */ + virtual void getEnthalpy_RT(doublereal* hrt) const; + + /** + * Get the nondimensional Entropies for the species + * standard states at the current T and P of the solution. + * + * Note, this is equal to the reference state entropies + * due to the zero volume expansivity: + * i.e., (dS/dp)_T = (dV/dT)_P = 0.0 + * + * @param sr Vector of length m_kk, which on return sr[k] + * will contain the nondimensional + * standard state entropy of species k. + */ + virtual void getEntropy_R(doublereal* sr) const; + + /** + * Get the nondimensional heat capacity at constant pressure + * function for the species + * standard states at the current T and P of the solution. + * \f[ + * Cp^0_k(T,P) = Cp^{ref}_k(T) + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ Cp^{ref}_k(T)\f$ is the constant pressure heat capacity + * of species k at the reference pressure, \f$p_{ref}\f$. + * + * @param cpr Vector of length m_kk, which on return cpr[k] + * will contain the nondimensional + * constant pressure heat capacity for species k. + */ + virtual void getCp_R(doublereal* cpr) const; + + /** + * Get the molar volumes of each species in their standard + * states at the current + * T and P of the solution. + * units = m^3 / kmol + */ + virtual void getStandardVolumes(doublereal *vol) const; + + //@} + /// @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 + * Chemical equilibrium. + * @{ + */ + + /** + * This method is used by the ChemEquil equilibrium solver. + * It sets the state such that the chemical potentials satisfy + * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m} + * \left(\frac{\lambda_m} {\hat R T}\right) \f] where + * \f$ \lambda_m \f$ is the element potential of element m. The + * temperature is unchanged. Any phase (ideal or not) that + * implements this method can be equilibrated by ChemEquil. + */ + virtual void setToEquilState(const doublereal* lambda_RT) { + err("setToEquilState"); + } + + // called by function 'equilibrate' in ChemEquil.h to transfer + // the element potentials to this object + void setElementPotentials(const vector_fp& lambda) { + m_lambda = lambda; + } + + void getElementPotentials(doublereal* lambda) { + copy(m_lambda.begin(), m_lambda.end(), lambda); + } + + //@} + + + /** + * @internal + * Set equation of state parameters. The number and meaning of + * these depends on the subclass. + * @param n number of parameters + * @param c array of \i n coefficients + * + */ + virtual void setParameters(int n, doublereal* c); + virtual void getParameters(int &n, doublereal * const c); + + /** + * Set equation of state parameter values from XML + * entries. This method is called by function importPhase in + * file importCTML.cpp when processing a phase definition in + * an input file. It should be overloaded in subclasses to set + * any parameters that are specific to that particular phase + * model. + * + * @param eosdata An XML_Node object corresponding to + * the "thermo" entry for this phase in the input file. + */ + virtual void setParametersFromXML(const XML_Node& eosdata); + + //--------------------------------------------------------- + /// @name Critical state properties. + /// These methods are only implemented by some subclasses. + + //@{ + + /// Critical temperature (K). + virtual doublereal critTemperature() const { + err("critTemperature"); return -1.0; + } + + /// Critical pressure (Pa). + virtual doublereal critPressure() const { + err("critPressure"); return -1.0; + } + + /// Critical density (kg/m3). + virtual doublereal critDensity() const { + err("critDensity"); return -1.0; + } + + //@} + + /// @name Saturation properties. + /// These methods are only implemented by subclasses that + /// implement full liquid-vapor equations of state. + /// + virtual doublereal satTemperature(doublereal p) const { + err("satTemperature"); return -1.0; + } + + virtual doublereal satPressure(doublereal t) const { + err("satPressure"); return -1.0; + } + + virtual doublereal vaporFraction() const { + err("vaprFraction"); return -1.0; + } + + virtual void setState_Tsat(doublereal t, doublereal x) { + err("setState_sat"); + } + + virtual void setState_Psat(doublereal p, doublereal x) { + err("setState_sat"); + } + + //@} + + + /* + * -------------- Utilities ------------------------------- + */ + + /** + * @internal Install a species thermodynamic property + * manager. The species thermodynamic property manager + * computes properties of the pure species for use in + * constructing solution properties. It is meant for internal + * use, and some classes derived from ThermoPhase may not use + * any species thermodynamic property manager. + */ + void setSpeciesThermo(SpeciesThermo* spthermo) + { m_spthermo = spthermo; } + + /** + * Return a reference to the species thermodynamic property + * manager. @todo This method will fail if no species thermo + * manager has been installed. + */ + SpeciesThermo& speciesThermo() { return *m_spthermo; } + + + + + /** + * initThermo() (virtual from ThermoPhase) + * + * This internal routine is responsible for setting up + * the internal storage. + */ + virtual void initThermo(); + + /** + * constructPhaseFile (virtual from here) + * + * Initialization of an IdealSolidSolnPhase phase using an + * xml file identified by its path + * + * 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. + */ + virtual void constructPhaseFile(string infile, string id=""); + + /** + * constructPhaseXML (virtual from here) + * + * This is the main routine for constructing the phase. + * It processes the XML file, and then it calls importPhase(). + * Then, initThermoXML() is called after importPhase(). + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + virtual void constructPhaseXML(XML_Node& phaseNode, string id); + + /** + * initThermoXML (virtual from ThermoPhase) + * + * + * This routine is called from importPhase() to finish + * up the initialization of the thermo object. It reads in the + * species molar volumes. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + virtual void initThermoXML(XML_Node& phaseNode, string id=""); + + + + /** + * Report the molar volume of species k + * + * units - \f$ m^3 kmol^-1 \f$ + */ + double speciesMolarVolume(int k) const; + + /** + * Fill in a return vector containing the species molar volumes + * units - \f$ m^3 kmol^-1 \f$ + */ + void getSpeciesMolarVolumes(double *smv) const; + //@} + + + + + protected: + /** + * Species molar volume \f$ m^3 kmol^-1 \f$ + */ + array_fp m_speciesMolarVolume; + /* + * Current pressure in Pascal + */ + double m_Pcurrent; + int m_formGC; + + /** + * Vector containing the species reference exp(-G/RT) functions + * at T = m_tlast + */ + mutable vector_fp m_expg0_RT; + + /** + * Vector of potential energies for the species. + */ + mutable vector_fp m_pe; + + /** + * Temporary array used in equilibrium calculations + */ + mutable vector_fp m_pp; + + /** + * vector of size m_kk, used as a temporary holding area. + */ + mutable vector_fp m_tmpV; + + private: + doublereal err(string msg) const; + + + void initLengths(); + }; + +} + +#endif + + + + + diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index eb93ebbf4..c986363f9 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -14,16 +14,30 @@ INCDIR = ../../../build/include/cantera/kernel/thermo INSTALL_TSC = ../../../bin/install_tsc do_ranlib = @DO_RANLIB@ +do_electro = @COMPILE_ELECTROLYTES@ +do_issp = @COMPILE_IDEAL_SOLUTIONS@ CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) # Extended Cantera Thermodynamics Object Files -CATHERMO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o \ - MolalityVPSSTP.o VPStandardStateTP.o \ - IdealSolidSolnPhase.o -CATHERMO_H = SingleSpeciesTP.h StoichSubstanceSSTP.h \ + +ifeq ($(do_electro),1) +ELECTRO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o \ + MolalityVPSSTP.o VPStandardStateTP.o \ + IdealSolidSolnPhase.o IdealMolalSoln.o +ELECTRO_H = SingleSpeciesTP.h StoichSubstanceSSTP.h \ MolalityVPSSTP.h VPStandardStateTP.h \ - IdealSolidSolnPhase.h + IdealSolidSolnPhase.h IdealMolalSoln.h +endif +ifeq ($(do_issp),1) +ISSP_OBJ = IdealSolidSolnPhase.o +ISSP_H = IdealSolidSolnPhase.h +endif + +CATHERMO_OBJ = $(ELECTRO_OBJ) $(ISSP_OBJ) + +CATHERMO_H = $(ELECTRO_H) $(ISSP_H) + CXX_INCLUDES = -I.. @CXX_INCLUDES@ LIB = @buildlib@/libcaThermo.a @@ -68,7 +82,7 @@ depends: @MAKE@ .depends .depends: $(DEPENDS) - cat *.d > .depends + cat $(DEPENDS) > .depends TAGS: etags *.h *.cpp diff --git a/Cantera/src/thermo/MolalityVPSSTP.cpp b/Cantera/src/thermo/MolalityVPSSTP.cpp index 4ccac77bf..5d428dd3e 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.cpp +++ b/Cantera/src/thermo/MolalityVPSSTP.cpp @@ -22,521 +22,536 @@ namespace Cantera { - /* - * Default 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. - */ - MolalityVPSSTP::MolalityVPSSTP() : - VPStandardStateTP(), - m_indexSolvent(0), - m_weightSolvent(18.01528), - m_xmolSolventMIN(0.01), - m_Mnaught(18.01528E-3) - { - } + /* + * Default 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. + */ + MolalityVPSSTP::MolalityVPSSTP() : + VPStandardStateTP(), + m_indexSolvent(0), + m_weightSolvent(18.01528), + m_xmolSolventMIN(0.01), + m_Mnaught(18.01528E-3) + { + } - /** - * Copy Constructor: - * - * Note this stuff will not work until the underlying phase - * has a working copy constructor - */ - MolalityVPSSTP::MolalityVPSSTP(const MolalityVPSSTP &b) : - VPStandardStateTP(), - m_indexSolvent(b.m_indexSolvent), - m_xmolSolventMIN(b.m_xmolSolventMIN), - m_Mnaught(b.m_Mnaught), - m_molalities(b.m_molalities) - { - throw CanteraError("MolalityVPSSTP::operator=()", - "Not Implemented Fully"); - *this = operator=(b); - } + /** + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + MolalityVPSSTP::MolalityVPSSTP(const MolalityVPSSTP &b) : + VPStandardStateTP(), + m_indexSolvent(b.m_indexSolvent), + m_xmolSolventMIN(b.m_xmolSolventMIN), + m_Mnaught(b.m_Mnaught), + m_molalities(b.m_molalities) + { + throw CanteraError("MolalityVPSSTP::operator=()", + "Not Implemented Fully"); + *this = operator=(b); + } - /* - * operator=() - * - * Note this stuff will not work until the underlying phase - * has a working assignment operator - */ - MolalityVPSSTP& MolalityVPSSTP:: - operator=(const MolalityVPSSTP &b) { - if (&b != this) { - VPStandardStateTP::operator=(b); - m_indexSolvent = b.m_indexSolvent; - m_weightSolvent = b.m_weightSolvent; - m_xmolSolventMIN = b.m_xmolSolventMIN; - m_Mnaught = b.m_Mnaught; - m_molalities = b.m_molalities; - } - throw CanteraError("MolalityVPSSTP::operator=()", - "Not Implemented Fully"); - return *this; + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + MolalityVPSSTP& MolalityVPSSTP:: + operator=(const MolalityVPSSTP &b) { + if (&b != this) { + VPStandardStateTP::operator=(b); + m_indexSolvent = b.m_indexSolvent; + m_weightSolvent = b.m_weightSolvent; + m_xmolSolventMIN = b.m_xmolSolventMIN; + m_Mnaught = b.m_Mnaught; + m_molalities = b.m_molalities; } + throw CanteraError("MolalityVPSSTP::operator=()", + "Not Implemented Fully"); + return *this; + } - /** - * - * ~MolalityVPSSTP(): (virtual) - * - * Destructor: does nothing: - * - */ - MolalityVPSSTP::~MolalityVPSSTP() { - } - - /** - * This routine duplicates the current object and returns - * a pointer to ThermoPhase. - */ - ThermoPhase* - MolalityVPSSTP::duplMyselfAsThermoPhase() { - MolalityVPSSTP* mtp = new MolalityVPSSTP(*this); - return (ThermoPhase *) mtp; - } - - /* - * -------------- Utilities ------------------------------- - */ - - /* - * setSolvent(): - * Utilities for Solvent ID and Molality - * Here we also calculate and store the molecular weight - * of the solvent and the m_Mnaught parameter. - * @param k index of the solvent. - */ - void MolalityVPSSTP::setSolvent(int k) { - if (k < 0 || k >= m_kk) { - throw CanteraError("MolalityVPSSTP::setSolute ", - "bad value"); - } - m_indexSolvent = k; - m_weightSolvent = molecularWeight(k); - m_Mnaught = m_weightSolvent / 1000.; + /** + * + * ~MolalityVPSSTP(): (virtual) + * + * Destructor: does nothing: + * + */ + MolalityVPSSTP::~MolalityVPSSTP() { + } + + /** + * This routine duplicates the current object and returns + * a pointer to ThermoPhase. + */ + ThermoPhase* + MolalityVPSSTP::duplMyselfAsThermoPhase() { + MolalityVPSSTP* mtp = new MolalityVPSSTP(*this); + return (ThermoPhase *) mtp; + } + + /* + * -------------- Utilities ------------------------------- + */ + + /* + * setSolvent(): + * Utilities for Solvent ID and Molality + * Here we also calculate and store the molecular weight + * of the solvent and the m_Mnaught parameter. + * @param k index of the solvent. + */ + void MolalityVPSSTP::setSolvent(int k) { + if (k < 0 || k >= m_kk) { + throw CanteraError("MolalityVPSSTP::setSolute ", + "bad value"); } + m_indexSolvent = k; + m_weightSolvent = molecularWeight(k); + m_Mnaught = m_weightSolvent / 1000.; + } - /* - * return the solvent id index number. - */ - int MolalityVPSSTP::solventIndex() const { - return m_indexSolvent; - } + /* + * return the solvent id index number. + */ + int MolalityVPSSTP::solventIndex() const { + return m_indexSolvent; + } - /** - * Sets the minimum mole fraction in the molality formulation. The - * minimum mole fraction must be in the range 0 to 0.9. - */ - void MolalityVPSSTP:: - setMoleFSolventMin(doublereal xmolSolventMIN) { - if (xmolSolventMIN <= 0.0) { - throw CanteraError("MolalityVPSSTP::setSolute ", "trouble"); - } else if (xmolSolventMIN > 0.9) { - throw CanteraError("MolalityVPSSTP::setSolute ", "trouble"); - } - m_xmolSolventMIN = xmolSolventMIN; + /** + * Sets the minimum mole fraction in the molality formulation. The + * minimum mole fraction must be in the range 0 to 0.9. + */ + void MolalityVPSSTP:: + setMoleFSolventMin(doublereal xmolSolventMIN) { + if (xmolSolventMIN <= 0.0) { + throw CanteraError("MolalityVPSSTP::setSolute ", "trouble"); + } else if (xmolSolventMIN > 0.9) { + throw CanteraError("MolalityVPSSTP::setSolute ", "trouble"); } + m_xmolSolventMIN = xmolSolventMIN; + } - /** - * Returns the minimum mole fraction in the molality formulation. - */ - doublereal MolalityVPSSTP::moleFSolventMin() const { - return m_xmolSolventMIN; + /** + * Returns the minimum mole fraction in the molality formulation. + */ + doublereal MolalityVPSSTP::moleFSolventMin() const { + return m_xmolSolventMIN; + } + + /** + * calcMolalities(): + * We calculate the vector of molalities of the species + * in the phase and store the result internally: + * \f[ + * m_i = (n_i) / (1000 * M_o * n_{o,p}) + * \f] + * where + * - \f$ M_o \f$ is the molecular weight of the solvent + * - \f$ n_o \f$ is the mole fraction of the solvent + * - \f$ n_i \f$ is the mole fraction of the solute. + * - \f$ n_{o,p} = max (n_{o, min}, n_o) \f$ + * - \f$ n_{o,min} \f$ = minimum mole fraction of solvent allowed + * in the denominator. + */ + void MolalityVPSSTP::calcMolalities() const { + getMoleFractions(DATA_PTR(m_molalities)); + double xmolSolvent = m_molalities[m_indexSolvent]; + if (xmolSolvent < m_xmolSolventMIN) { + xmolSolvent = m_xmolSolventMIN; } - - /** - * calcMolalities(): - * We calculate the vector of molalities of the species - * in the phase and store the result internally: - * \f[ - * m_i = (n_i) / (1000 * M_o * n_{o,p}) - * \f] - * where - * - \f$ M_o \f$ is the molecular weight of the solvent - * - \f$ n_o \f$ is the mole fraction of the solvent - * - \f$ n_i \f$ is the mole fraction of the solute. - * - \f$ n_{o,p} = max (n_{o, min}, n_o) \f$ - * - \f$ n_{o,min} \f$ = minimum mole fraction of solvent allowed - * in the denominator. - */ - void MolalityVPSSTP::calcMolalities() const { - getMoleFractions(DATA_PTR(m_molalities)); - double xmolSolvent = m_molalities[m_indexSolvent]; - if (xmolSolvent < m_xmolSolventMIN) { - xmolSolvent = m_xmolSolventMIN; - } - double denomInv = 1.0/ (m_Mnaught * xmolSolvent); - for (int k = 0; k < m_kk; k++) { - m_molalities[k] *= denomInv; - } + double denomInv = 1.0/ (m_Mnaught * xmolSolvent); + for (int k = 0; k < m_kk; k++) { + m_molalities[k] *= denomInv; } + } - /** - * getMolalities(): - * We calculate the vector of molalities of the species - * in the phase - * \f[ - * m_i = (n_i) / (1000 * M_o * n_{o,p}) - * \f] - * where - * - \f$ M_o \f$ is the molecular weight of the solvent - * - \f$ n_o \f$ is the mole fraction of the solvent - * - \f$ n_i \f$ is the mole fraction of the solute. - * - \f$ n_{o,p} = max (n_{o, min}, n_o) \f$ - * - \f$ n_{o,min} \f$ = minimum mole fraction of solvent allowed - * in the denominator. - */ - void MolalityVPSSTP::getMolalities(doublereal * const molal) const { - calcMolalities(); - for (int k = 0; k < m_kk; k++) { - molal[k] = m_molalities[k]; + /** + * getMolalities(): + * We calculate the vector of molalities of the species + * in the phase + * \f[ + * m_i = (n_i) / (1000 * M_o * n_{o,p}) + * \f] + * where + * - \f$ M_o \f$ is the molecular weight of the solvent + * - \f$ n_o \f$ is the mole fraction of the solvent + * - \f$ n_i \f$ is the mole fraction of the solute. + * - \f$ n_{o,p} = max (n_{o, min}, n_o) \f$ + * - \f$ n_{o,min} \f$ = minimum mole fraction of solvent allowed + * in the denominator. + */ + void MolalityVPSSTP::getMolalities(doublereal * const molal) const { + calcMolalities(); + for (int k = 0; k < m_kk; k++) { + molal[k] = m_molalities[k]; + } + } + + /** + * setMolalities(): + * We are supplied with the molalities of all of the + * solute species. We then calculate the mole fractions of all + * species and update the ThermoPhase object. + * + * m_i = (n_i) / (W_o/1000 * n_o_p) + * + * where M_o is the molecular weight of the solvent + * n_o is the mole fraction of the solvent + * n_i is the mole fraction of the solute. + * n_o_p = max (n_o_min, n_o) + * n_o_min = minimum mole fraction of solvent allowed + * in the denominator. + */ + void MolalityVPSSTP::setMolalities(const doublereal * const molal) { + + double Lsum = 1.0 / m_Mnaught; + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + m_molalities[k] = molal[k]; + Lsum += molal[k]; } } - - /** - * setMolalities(): - * We are supplied with the molalities of all of the - * solute species. We then calculate the mole fractions of all - * species and update the ThermoPhase object. - * - * m_i = (n_i) / (W_o/1000 * n_o_p) - * - * where M_o is the molecular weight of the solvent - * n_o is the mole fraction of the solvent - * n_i is the mole fraction of the solute. - * n_o_p = max (n_o_min, n_o) - * n_o_min = minimum mole fraction of solvent allowed - * in the denominator. - */ - void MolalityVPSSTP::setMolalities(const doublereal * const molal) { - - double Lsum = 1.0 / m_Mnaught; - for (int k = 0; k < m_kk; k++) { - if (k != m_indexSolvent) { - m_molalities[k] = molal[k]; - Lsum += molal[k]; - } - } - double tmp = 1.0 / Lsum; - m_molalities[m_indexSolvent] = tmp / m_Mnaught; - double sum = m_molalities[m_indexSolvent]; - for (int k = 0; k < m_kk; k++) { - if (k != m_indexSolvent) { - m_molalities[k] = tmp * molal[k]; - sum += m_molalities[k]; - } - } - if (sum != 1.0) { - tmp = 1.0 / sum; - for (int k = 0; k < m_kk; k++) { - m_molalities[k] *= tmp; - } - } - setMoleFractions(DATA_PTR(m_molalities)); - /* - * Essentially we don't trust the input: We calculate - * the molalities from the mole fractions that we - * just obtained. - */ - calcMolalities(); + double tmp = 1.0 / Lsum; + m_molalities[m_indexSolvent] = tmp / m_Mnaught; + double sum = m_molalities[m_indexSolvent]; + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + m_molalities[k] = tmp * molal[k]; + sum += m_molalities[k]; + } } - + if (sum != 1.0) { + tmp = 1.0 / sum; + for (int k = 0; k < m_kk; k++) { + m_molalities[k] *= tmp; + } + } + setMoleFractions(DATA_PTR(m_molalities)); /* - * setMolalitiesByName() - * - * This routine sets the molalities by name - * HKM -> Might need to be more complicated here, setting - * neutrals so that the existing mole fractions are - * preserved. + * Essentially we don't trust the input: We calculate + * the molalities from the mole fractions that we + * just obtained. */ - void MolalityVPSSTP::setMolalitiesByName(compositionMap& mMap) { - int kk = nSpecies(); - doublereal x; - /* - * Get a vector of mole fractions - */ - vector_fp mf(kk, 0.0); - getMoleFractions(DATA_PTR(mf)); - double xmolS = mf[m_indexSolvent]; - double xmolSmin = max(xmolS, m_xmolSolventMIN); - compositionMap::iterator p; - for (int k = 0; k < kk; k++) { - p = mMap.find(speciesName(k)); - if (p != mMap.end()) { - x = mMap[speciesName(k)]; - if (x > 0.0) { - mf[k] = x * m_Mnaught * xmolSmin; - } + calcMolalities(); + } + + /* + * setMolalitiesByName() + * + * This routine sets the molalities by name + * HKM -> Might need to be more complicated here, setting + * neutrals so that the existing mole fractions are + * preserved. + */ + void MolalityVPSSTP::setMolalitiesByName(compositionMap& mMap) { + int kk = nSpecies(); + doublereal x; + /* + * Get a vector of mole fractions + */ + vector_fp mf(kk, 0.0); + getMoleFractions(DATA_PTR(mf)); + double xmolS = mf[m_indexSolvent]; + double xmolSmin = max(xmolS, m_xmolSolventMIN); + compositionMap::iterator p; + for (int k = 0; k < kk; k++) { + p = mMap.find(speciesName(k)); + if (p != mMap.end()) { + x = mMap[speciesName(k)]; + if (x > 0.0) { + mf[k] = x * m_Mnaught * xmolSmin; + } + } + } + /* + * check charge neutrality + */ + int largePos = -1; + double cPos = 0.0; + int largeNeg = -1; + double cNeg = 0.0; + double sum = 0.0; + for (int k = 0; k < kk; k++) { + double ch = charge(k); + if (mf[k] > 0.0) { + if (ch > 0.0) { + if (ch * mf[k] > cPos) { + largePos = k; + cPos = ch * mf[k]; } } - /* - * check charge neutrality - */ - int largePos = -1; - double cPos = 0.0; - int largeNeg = -1; - double cNeg = 0.0; - double sum = 0.0; - for (int k = 0; k < kk; k++) { - double ch = charge(k); - if (mf[k] > 0.0) { - if (ch > 0.0) { - if (ch * mf[k] > cPos) { - largePos = k; - cPos = ch * mf[k]; - } - } - if (ch < 0.0) { - if (fabs(ch) * mf[k] > cNeg) { - largeNeg = k; - cNeg = fabs(ch) * mf[k]; - } - } - } - sum += mf[k] * ch; + if (ch < 0.0) { + if (fabs(ch) * mf[k] > cNeg) { + largeNeg = k; + cNeg = fabs(ch) * mf[k]; + } } - if (sum != 0.0) { - if (sum > 0.0) { - if (cPos > sum) { - mf[largePos] -= sum / charge(largePos); - } else { - throw CanteraError("MolalityVPSSTP:setMolalitiesbyName", - "unbalanced charges"); - } - } else { - if (cNeg > (-sum)) { - mf[largeNeg] -= (-sum) / fabs(charge(largeNeg)); - } else { - throw CanteraError("MolalityVPSSTP:setMolalitiesbyName", - "unbalanced charges"); - } - } + } + sum += mf[k] * ch; + } + if (sum != 0.0) { + if (sum > 0.0) { + if (cPos > sum) { + mf[largePos] -= sum / charge(largePos); + } else { + throw CanteraError("MolalityVPSSTP:setMolalitiesbyName", + "unbalanced charges"); + } + } else { + if (cNeg > (-sum)) { + mf[largeNeg] -= (-sum) / fabs(charge(largeNeg)); + } else { + throw CanteraError("MolalityVPSSTP:setMolalitiesbyName", + "unbalanced charges"); + } + } - } - sum = 0.0; - for (int k = 0; k < kk; k++) { - sum += mf[k]; - } - sum = 1.0/sum; - for (int k = 0; k < kk; k++) { - mf[k] *= sum; - } - setMoleFractions(DATA_PTR(mf)); - /* - * After we formally set the mole fractions, we - * calculate the molalities again and store it in - * this object. - */ - calcMolalities(); } - + sum = 0.0; + for (int k = 0; k < kk; k++) { + sum += mf[k]; + } + sum = 1.0/sum; + for (int k = 0; k < kk; k++) { + mf[k] *= sum; + } + setMoleFractions(DATA_PTR(mf)); /* - * setMolalitiesByNames() - * - * Set the molalities of the solutes by name + * After we formally set the mole fractions, we + * calculate the molalities again and store it in + * this object. */ - void MolalityVPSSTP::setMolalitiesByName(const string& x) { - compositionMap xx; - int kk = nSpecies(); - for (int k = 0; k < kk; k++) { - xx[speciesName(k)] = -1.0; - } - parseCompString(x, xx); - setMolalitiesByName(xx); + calcMolalities(); + } + + /* + * setMolalitiesByNames() + * + * Set the molalities of the solutes by name + */ + void MolalityVPSSTP::setMolalitiesByName(const string& x) { + compositionMap xx; + int kk = nSpecies(); + for (int k = 0; k < kk; k++) { + xx[speciesName(k)] = -1.0; } + parseCompString(x, xx); + setMolalitiesByName(xx); + } - /* - * ------------ Molar Thermodynamic Properties ---------------------- - */ + /* + * ------------ Molar Thermodynamic Properties ---------------------- + */ - /* - * - Activities, Standard States, Activity Concentrations ----------- - */ + /* + * - Activities, Standard States, Activity Concentrations ----------- + */ - /** - * This method returns the activity convention. - * Currently, there are two activity conventions - * Molar-based activities - * Unit activity of species at either a hypothetical pure - * solution of the species or at a hypothetical - * pure ideal solution at infinite dilution - * cAC_CONVENTION_MOLAR 0 - * - default - * - * Molality based activities - * (unit activity of solutes at a hypothetical 1 molal - * solution referenced to infinite dilution at all - * pressures and temperatures). - * (solvent is still on molar basis). - * cAC_CONVENTION_MOLALITY 1 - * - * We set the convention to molality here. - */ - int MolalityVPSSTP::activityConvention() const { - return cAC_CONVENTION_MOLALITY; + /** + * This method returns the activity convention. + * Currently, there are two activity conventions + * Molar-based activities + * Unit activity of species at either a hypothetical pure + * solution of the species or at a hypothetical + * pure ideal solution at infinite dilution + * cAC_CONVENTION_MOLAR 0 + * - default + * + * Molality based activities + * (unit activity of solutes at a hypothetical 1 molal + * solution referenced to infinite dilution at all + * pressures and temperatures). + * (solvent is still on molar basis). + * cAC_CONVENTION_MOLALITY 1 + * + * We set the convention to molality here. + */ + int MolalityVPSSTP::activityConvention() const { + return cAC_CONVENTION_MOLALITY; + } + + /** + * Get the array of non-dimensional activity coefficients at + * the current solution temperature, pressure, and + * solution concentration. + * These are mole fraction based activity coefficients. In this + * object, their calculation is based on translating the values + * of Molality based activity coefficients. + * See Denbigh p. 278 for a thorough discussion. + * + * Note, the solvent is treated differently. getMolalityActivityCoeff() + * returns the molar based solvent activity coefficient already. + * Therefore, we do not have to divide by x_s here. + */ + void MolalityVPSSTP::getActivityCoefficients(doublereal* ac) const { + getMolalityActivityCoefficients(ac); + double xmolSolvent = moleFraction(m_indexSolvent); + if (xmolSolvent < m_xmolSolventMIN) { + xmolSolvent = m_xmolSolventMIN; } - - /** - * Get the array of non-dimensional activity coefficients at - * the current solution temperature, pressure, and - * solution concentration. - * These are mole fraction based activity coefficients. In this - * object, their calculation is based on translating the values - * of Molality based activity coefficients. - * See Denbigh p. 278 for a thorough discussion. - * - * Note, the solvent is treated differently. getMolalityActivityCoeff() - * returns the molar based solvent activity coefficient already. - * Therefore, we do not have to divide by x_s here. - */ - void MolalityVPSSTP::getActivityCoefficients(doublereal* ac) const { - getMolalityActivityCoefficients(ac); - double xmolSolvent = moleFraction(m_indexSolvent); - if (xmolSolvent < m_xmolSolventMIN) { - xmolSolvent = m_xmolSolventMIN; - } - for (int k = 0; k < m_kk; k++) { - if (k != m_indexSolvent) { - ac[k] /= xmolSolvent; - } - } + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + ac[k] /= xmolSolvent; + } } + } - /** - * osmotic coefficient: - * - * Calculate the osmotic coefficient of the solvent. Note there - * are lots of definitions of the osmotic coefficient floating - * around. We use the one defined in the Pitzer's book: - * (Activity Coeff in Electrolyte Solutions, K. S. Pitzer - * CRC Press, Boca Raton, 1991, p. 85, Eqn. 28). - * - * Definition: - * - sum(m_i) * Mnaught * oc = ln(activity_solvent) - */ - doublereal MolalityVPSSTP::osmoticCoefficient() const { - /* - * First, we calculate the activities all over again - */ - vector_fp act(m_kk); - getActivities(DATA_PTR(act)); - /* - * Then, we calculate the sum of the solvent molalities - */ - double sum = 0; - for (int k = 0; k < m_kk; k++) { - if (k != m_indexSolvent) { - sum += MAX(m_molalities[k], 0.0); - } - } - double oc = 1.0; - double lac = log(act[m_indexSolvent]); - if (sum > 1.0E-200) { - oc = - lac / (m_Mnaught * sum); - } - return oc; - } - + /** + * osmotic coefficient: + * + * Calculate the osmotic coefficient of the solvent. Note there + * are lots of definitions of the osmotic coefficient floating + * around. We use the one defined in the Pitzer's book: + * (Activity Coeff in Electrolyte Solutions, K. S. Pitzer + * CRC Press, Boca Raton, 1991, p. 85, Eqn. 28). + * + * Definition: + * - sum(m_i) * Mnaught * oc = ln(activity_solvent) + */ + doublereal MolalityVPSSTP::osmoticCoefficient() const { /* - * ------------ Partial Molar Properties of the Solution ------------ + * First, we calculate the activities all over again */ - - - doublereal MolalityVPSSTP::err(string msg) const { - throw CanteraError("MolalityVPSSTP","Base class method " - +msg+" called. Equation of state type: "+int2str(eosType())); - return 0; - } - - /** - * Returns the units of the standard and general concentrations - * Note they have the same units, as their divisor is - * defined to be equal to the activity of the kth species - * in the solution, which is unitless. - * - * This routine is used in print out applications where the - * units are needed. Usually, MKS units are assumed throughout - * the program and in the XML input files. - * - * On return uA contains the powers of the units (MKS assumed) - * of the standard concentrations and generalized concentrations - * for the kth species. - * - * uA[0] = kmol units - default = 1 - * uA[1] = m units - default = -nDim(), the number of spatial - * dimensions in the Phase class. - * uA[2] = kg units - default = 0; - * uA[3] = Pa(pressure) units - default = 0; - * uA[4] = Temperature units - default = 0; - * uA[5] = time units - default = 0 + vector_fp act(m_kk); + getActivities(DATA_PTR(act)); + /* + * Then, we calculate the sum of the solvent molalities */ - void MolalityVPSSTP::getUnitsStandardConc(double *uA, int k, int sizeUA) { - for (int i = 0; i < sizeUA; i++) { - if (i == 0) uA[0] = 1.0; - if (i == 1) uA[1] = -nDim(); - if (i == 2) uA[2] = 0.0; - if (i == 3) uA[3] = 0.0; - if (i == 4) uA[4] = 0.0; - if (i == 5) uA[5] = 0.0; - } + double sum = 0; + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + sum += MAX(m_molalities[k], 0.0); + } } + double oc = 1.0; + double lac = log(act[m_indexSolvent]); + if (sum > 1.0E-200) { + oc = - lac / (m_Mnaught * sum); + } + return oc; + } + + /* + * ------------ Partial Molar Properties of the Solution ------------ + */ + + + doublereal MolalityVPSSTP::err(string msg) const { + throw CanteraError("MolalityVPSSTP","Base class method " + +msg+" called. Equation of state type: "+int2str(eosType())); + return 0; + } + + /** + * Returns the units of the standard and general concentrations + * Note they have the same units, as their divisor is + * defined to be equal to the activity of the kth species + * in the solution, which is unitless. + * + * This routine is used in print out applications where the + * units are needed. Usually, MKS units are assumed throughout + * the program and in the XML input files. + * + * On return uA contains the powers of the units (MKS assumed) + * of the standard concentrations and generalized concentrations + * for the kth species. + * + * uA[0] = kmol units - default = 1 + * uA[1] = m units - default = -nDim(), the number of spatial + * dimensions in the Phase class. + * uA[2] = kg units - default = 0; + * uA[3] = Pa(pressure) units - default = 0; + * uA[4] = Temperature units - default = 0; + * uA[5] = time units - default = 0 + */ + void MolalityVPSSTP::getUnitsStandardConc(double *uA, int k, int sizeUA) { + for (int i = 0; i < sizeUA; i++) { + if (i == 0) uA[0] = 1.0; + if (i == 1) uA[1] = -nDim(); + if (i == 2) uA[2] = 0.0; + if (i == 3) uA[3] = 0.0; + if (i == 4) uA[4] = 0.0; + if (i == 5) uA[5] = 0.0; + } + } + /* + * Set the thermodynamic state. + */ + void MolalityVPSSTP::setStateFromXML(const XML_Node& state) { + VPStandardStateTP::setStateFromXML(state); + string comp = getString(state,"soluteMolalities"); + if (comp != "") { + setMolalitiesByName(comp); + } + if (state.hasChild("pressure")) { + double p = getFloat(state, "pressure", "pressure"); + setPressure(p); + } + } + + /** + * Set the temperature (K), pressure (Pa), and molalities + * (gmol kg-1) of the solutes + */ + void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p, + const doublereal * const molalities) { + setMolalities(molalities); + setTemperature(t); + setPressure(p); + } + + /** Set the temperature (K), pressure (Pa), and molalities. */ + void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p, compositionMap& m) { + setMolalitiesByName(m); + setTemperature(t); + setPressure(p); + } + + /** Set the temperature (K), pressure (Pa), and molality. */ + void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p, const string& m) { + setMolalitiesByName(m); + setTemperature(t); + setPressure(p); + } + + + /** + * @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 MolalityVPSSTP::initThermo() { + initLengths(); + VPStandardStateTP::initThermo(); + /* - * Set the thermodynamic state. + * The solvent defaults to species 0 */ - void MolalityVPSSTP::setStateFromXML(const XML_Node& state) { - VPStandardStateTP::setStateFromXML(state); - string comp = getString(state,"soluteMolalities"); - if (comp != "") { - setMolalitiesByName(comp); - } - if (state.hasChild("pressure")) { - double p = getFloat(state, "pressure", "pressure"); - setPressure(p); - } - } - - /** - * Set the temperature (K), pressure (Pa), and molalities - * (gmol kg-1) of the solutes - */ - void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p, - const doublereal * const molalities) { - setMolalities(molalities); - setTemperature(t); - setPressure(p); - } - - /** - * @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 MolalityVPSSTP::initThermo() { - initLengths(); - VPStandardStateTP::initThermo(); - - /* - * The solvent defaults to species 0 - */ - setSolvent(0); - } + setSolvent(0); + } void MolalityVPSSTP::initLengths() { int m_kk = nSpecies(); diff --git a/Cantera/src/thermo/MolalityVPSSTP.h b/Cantera/src/thermo/MolalityVPSSTP.h index dcdc2ad57..20cfa5db0 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.h +++ b/Cantera/src/thermo/MolalityVPSSTP.h @@ -25,457 +25,462 @@ namespace Cantera { - /** - * @ingroup thermoprops - */ + /** + * @ingroup thermoprops + */ - /** - * MolalityVPSSTP is 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. - */ - class MolalityVPSSTP : public VPStandardStateTP { + /** + * MolalityVPSSTP is 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. + */ + class MolalityVPSSTP : public VPStandardStateTP { - public: + public: - /// Constructors - MolalityVPSSTP(); - MolalityVPSSTP(const MolalityVPSSTP &); - /// Assignment operator - MolalityVPSSTP& operator=(const MolalityVPSSTP&); + /// Constructors + MolalityVPSSTP(); + MolalityVPSSTP(const MolalityVPSSTP &); + /// Assignment operator + MolalityVPSSTP& operator=(const MolalityVPSSTP&); - /// Destructor. - virtual ~MolalityVPSSTP(); + /// Destructor. + virtual ~MolalityVPSSTP(); - /** - * 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(); + /** + * 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(); - /** - * - * @name Utilities - * @{ - */ + /** + * + * @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 { return 0; } + /** + * 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 { return 0; } - /** - * @} - * @name Molar Thermodynamic Properties - * @{ - */ + /** + * @} + * @name Molar Thermodynamic Properties + * @{ + */ - /** - * @} - * @name Utilities for Solvent ID and Molality - * @{ - */ + /** + * @} + * @name Utilities for Solvent ID and Molality + * @{ + */ - /** - * This routine sets the index number of the solvent for - * the phase. - * - * Note, having a solvent - * is a precursor to many things having to do with molality. - * - * @param k the solvent index number - */ - void setSolvent(int k); + /** + * This routine sets the index number of the solvent for + * the phase. + * + * Note, having a solvent + * is a precursor to many things having to do with molality. + * + * @param k the solvent index number + */ + void setSolvent(int k); - /** - * Sets the minimum mole fraction in the molality formulation. - * Note the molality formulation is singular in the limit that - * the solvent mole fraction goes to zero. Numerically, how - * this limit is treated and resolved is an ongoing issue within - * Cantera. - */ - void setMoleFSolventMin(doublereal xmolSolventMIN); + /** + * Sets the minimum mole fraction in the molality formulation. + * Note the molality formulation is singular in the limit that + * the solvent mole fraction goes to zero. Numerically, how + * this limit is treated and resolved is an ongoing issue within + * Cantera. + */ + void setMoleFSolventMin(doublereal xmolSolventMIN); - /** - * Returns the solvent index. - */ - int solventIndex() const; + /** + * Returns the solvent index. + */ + int solventIndex() const; - /** - * Returns the minimum mole fraction in the molality - * formulation. - */ - doublereal moleFSolventMin() const; + /** + * Returns the minimum mole fraction in the molality + * formulation. + */ + doublereal moleFSolventMin() const; - /** - * Calculates the molality of all species and - * stores the result internally. - */ - void calcMolalities() const; + /** + * Calculates the molality of all species and + * stores the result internally. + */ + void calcMolalities() const; - /** - * getMolalities() - * This function will return the molalities of the - * species. - * - */ - void getMolalities(doublereal * const molal) const; + /** + * getMolalities() + * This function will return the molalities of the + * species. + * + */ + void getMolalities(doublereal * const molal) const; - void setMolalities(const doublereal * const molal); - void setMolalitiesByName(compositionMap& xMap); - void setMolalitiesByName(const string &); + void setMolalities(const doublereal * const molal); + void setMolalitiesByName(compositionMap& xMap); + void setMolalitiesByName(const string &); - /** - * @} - * @name Mechanical Properties - * @{ - */ + /** + * @} + * @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 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. - * @{ - */ + /** + * @} + * @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 the activity convention. - * Currently, there are two activity conventions - * Molar-based activities - * Unit activity of species at either a hypothetical pure - * solution of the species or at a hypothetical - * pure ideal solution at infinite dilution - * cAC_CONVENTION_MOLAR 0 - * - default - * - * Molality based acvtivities - * (unit activity of solutes at a hypothetical 1 molal - * solution referenced to infinite dilution at all - * pressures and temperatures). - * cAC_CONVENTION_MOLALITY 1 - * - * We set the convention to molality here. - */ - int activityConvention() const; + /** + * This method returns the activity convention. + * Currently, there are two activity conventions + * Molar-based activities + * Unit activity of species at either a hypothetical pure + * solution of the species or at a hypothetical + * pure ideal solution at infinite dilution + * cAC_CONVENTION_MOLAR 0 + * - default + * + * Molality based acvtivities + * (unit activity of solutes at a hypothetical 1 molal + * solution referenced to infinite dilution at all + * pressures and temperatures). + * cAC_CONVENTION_MOLALITY 1 + * + * We set the convention to molality here. + */ + int activityConvention() const; - /** - * This method returns an array of generalized concentrations - * \f$ C_k\f$ that are defined such that - * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ - * is a standard concentration - * defined below. These generalized concentrations are used - * by kinetics manager classes to compute the forward and - * reverse rates of elementary reactions. - * - * @param c Array of generalized concentrations. The - * units depend upon the implementation of the - * reaction rate expressions within the phase. - */ - virtual void getActivityConcentrations(doublereal* c) const { - err("getActivityConcentrations"); - } + /** + * This method returns an array of generalized concentrations + * \f$ C_k\f$ that are defined such that + * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ + * is a standard concentration + * defined below. These generalized concentrations are used + * by kinetics manager classes to compute the forward and + * reverse rates of elementary reactions. + * + * @param c Array of generalized concentrations. The + * units depend upon the implementation of the + * reaction rate expressions within the phase. + */ + virtual void getActivityConcentrations(doublereal* c) const { + err("getActivityConcentrations"); + } - /** - * 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. - */ - virtual doublereal standardConcentration(int k=0) const { - err("standardConcentration"); - return -1.0; - } + /** + * 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. + */ + virtual doublereal standardConcentration(int k=0) const { + err("standardConcentration"); + return -1.0; + } - /** - * Returns the natural logarithm of the standard - * concentration of the kth species - */ - virtual doublereal logStandardConc(int k=0) const { - err("logStandardConc"); - return -1.0; - } + /** + * Returns the natural logarithm of the standard + * concentration of the kth species + */ + virtual doublereal logStandardConc(int k=0) const { + err("logStandardConc"); + return -1.0; + } - /** - * Returns the units of the standard and generalized - * concentrations Note they have the same units, as their - * ratio is defined to be equal to the activity of the kth - * species in the solution, which is unitless. - * - * This routine is used in print out applications where the - * units are needed. Usually, MKS units are assumed throughout - * the program and in the XML input files. - * - * uA[0] = kmol units - default = 1 - * uA[1] = m units - default = -nDim(), the number of spatial - * dimensions in the Phase class. - * uA[2] = kg units - default = 0; - * uA[3] = Pa(pressure) units - default = 0; - * uA[4] = Temperature units - default = 0; - * uA[5] = time units - default = 0 - */ - virtual void getUnitsStandardConc(double *uA, int k = 0, - int sizeUA = 6); + /** + * Returns the units of the standard and generalized + * concentrations Note they have the same units, as their + * ratio is defined to be equal to the activity of the kth + * species in the solution, which is unitless. + * + * This routine is used in print out applications where the + * units are needed. Usually, MKS units are assumed throughout + * the program and in the XML input files. + * + * uA[0] = kmol units - default = 1 + * uA[1] = m units - default = -nDim(), the number of spatial + * dimensions in the Phase class. + * uA[2] = kg units - default = 0; + * uA[3] = Pa(pressure) units - default = 0; + * uA[4] = Temperature units - default = 0; + * uA[5] = time units - default = 0 + */ + virtual void getUnitsStandardConc(double *uA, int k = 0, + int sizeUA = 6); - /** - * Get the array of non-dimensional activities (molality - * based for this class and classes that derive from it) at - * the current solution temperature, pressure, and - * solution concentration. - */ - virtual void getActivities(doublereal* ac) const { - err("getActivities"); - } + /** + * Get the array of non-dimensional activities (molality + * based for this class and classes that derive from it) at + * the current solution temperature, pressure, and + * solution concentration. + */ + virtual void getActivities(doublereal* ac) const { + err("getActivities"); + } - /** - * Get the array of non-dimensional activity coefficients at - * the current solution temperature, pressure, and - * solution concentration. - * These are mole fraction based activity coefficients. In this - * object, their calculation is based on translating the values - * of Molality based activity coefficients. - * See Denbigh p. 278 for a thorough discussion - */ - void getActivityCoefficients(doublereal* ac) const; + /** + * Get the array of non-dimensional activity coefficients at + * the current solution temperature, pressure, and + * solution concentration. + * These are mole fraction based activity coefficients. In this + * object, their calculation is based on translating the values + * of Molality based activity coefficients. + * See Denbigh p. 278 for a thorough discussion + */ + void getActivityCoefficients(doublereal* ac) const; - /** - * Get the array of non-dimensional molality based - * activity coefficients at the current solution temperature, - * pressure, and solution concentration. - * See Denbigh p. 278 for a thorough discussion - */ - virtual void getMolalityActivityCoefficients(doublereal *acMolality) - const { - err("getMolalityActivityCoefficients"); - } + /** + * Get the array of non-dimensional molality based + * activity coefficients at the current solution temperature, + * pressure, and solution concentration. + * See Denbigh p. 278 for a thorough discussion + */ + virtual void getMolalityActivityCoefficients(doublereal *acMolality) + const { + err("getMolalityActivityCoefficients"); + } - /** - * Calculate the osmotic coefficient - * units = dimensionless - */ - virtual double osmoticCoefficient() const; + /** + * Calculate the osmotic coefficient + * units = dimensionless + */ + virtual double osmoticCoefficient() const; - //@} - /// @name Partial Molar Properties of the Solution - //@{ + //@} + /// @name Partial Molar Properties of the Solution + //@{ - /** - * 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 - */ - void getElectrochemPotentials(doublereal* mu) const { - getChemPotentials(mu); - double ve = Faraday * electricPotential(); - for (int k = 0; k < m_kk; k++) { - mu[k] += ve*charge(k); - } - } + /** + * 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 + */ + void getElectrochemPotentials(doublereal* mu) const { + getChemPotentials(mu); + double ve = Faraday * electricPotential(); + for (int k = 0; k < m_kk; k++) { + mu[k] += ve*charge(k); + } + } - //@} - /// @name Properties of the Standard State of the Species in the Solution - //@{ + //@} + /// @name Properties of the Standard State of the Species in the Solution + //@{ - //@} - /// @name Thermodynamic Values for the Species Reference States - //@{ + //@} + /// @name Thermodynamic Values for the Species Reference States + //@{ - /////////////////////////////////////////////////////// - // - // The methods below are not virtual, and should not - // be overloaded. - // - ////////////////////////////////////////////////////// + /////////////////////////////////////////////////////// + // + // The methods below are not virtual, and should not + // be overloaded. + // + ////////////////////////////////////////////////////// - /** - * @name Specific Properties - * @{ - */ + /** + * @name Specific Properties + * @{ + */ - /** - * @name Setting the State - * - * These methods set all or part of the thermodynamic - * state. - * @{ - */ + /** + * @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. - * @{ - */ + /** + * @name Chemical Equilibrium + * Routines that implement the Chemical equilibrium capability + * for a single phase, based on the element-potential method. + * @{ + */ - /** - * This method is used by the ChemEquil element-potential - * based equilibrium solver. - * It sets the state such that the chemical potentials of the - * species within the current phase satisfy - * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m} - * \left(\frac{\lambda_m} {\hat R T}\right) \f] where - * \f$ \lambda_m \f$ is the element potential of element m. The - * temperature is unchanged. Any phase (ideal or not) that - * implements this method can be equilibrated by ChemEquil. - */ - virtual void setToEquilState(const doublereal* lambda_RT) { - err("setToEquilState"); - } + /** + * This method is used by the ChemEquil element-potential + * based equilibrium solver. + * It sets the state such that the chemical potentials of the + * species within the current phase satisfy + * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m} + * \left(\frac{\lambda_m} {\hat R T}\right) \f] where + * \f$ \lambda_m \f$ is the element potential of element m. The + * temperature is unchanged. Any phase (ideal or not) that + * implements this method can be equilibrated by ChemEquil. + */ + virtual void setToEquilState(const doublereal* lambda_RT) { + err("setToEquilState"); + } - // called by function 'equilibrate' in ChemEquil.h to transfer - // the element potentials to this object - void setElementPotentials(const vector_fp& lambda) { - m_lambda = lambda; - } + // called by function 'equilibrate' in ChemEquil.h to transfer + // the element potentials to this object + void setElementPotentials(const vector_fp& lambda) { + m_lambda = lambda; + } - void getElementPotentials(doublereal* lambda) { - copy(m_lambda.begin(), m_lambda.end(), lambda); - } + void getElementPotentials(doublereal* lambda) { + copy(m_lambda.begin(), m_lambda.end(), lambda); + } - //@} + //@} - /** - * Set equation of state parameter values from XML - * entries. This method is called by function importPhase in - * file importCTML.cpp when processing a phase definition in - * an input file. It should be overloaded in subclasses to set - * any parameters that are specific to that particular phase - * model. - * - * The MolalityVPSSTP object defines a new method for setting - * the concentrations of a phase. The new method is defined by a - * block called "soluteMolalities". If this block - * is found, the concentrations within that phase are - * set to the "name":"molalities pairs found within that - * XML block. The solvent concentration is then set - * to everything else. - * - * @param eosdata An XML_Node object corresponding to - * the "thermo" entry for this phase in the input file. - * - */ - virtual void setStateFromXML(const XML_Node& state); + /** + * Set equation of state parameter values from XML + * entries. This method is called by function importPhase in + * file importCTML.cpp when processing a phase definition in + * an input file. It should be overloaded in subclasses to set + * any parameters that are specific to that particular phase + * model. + * + * The MolalityVPSSTP object defines a new method for setting + * the concentrations of a phase. The new method is defined by a + * block called "soluteMolalities". If this block + * is found, the concentrations within that phase are + * set to the "name":"molalities pairs found within that + * XML block. The solvent concentration is then set + * to everything else. + * + * @param eosdata An XML_Node object corresponding to + * the "thermo" entry for this phase in the input file. + * + */ + virtual void setStateFromXML(const XML_Node& state); - /// 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. + /// 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(); + /** + * @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, string id); + /** + * 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, string id); - /** - * Set the temperature (K), pressure (Pa), and molalities - * (gmol kg-1) of the solutes - */ - void setState_TPM(doublereal t, doublereal p, - const doublereal * const molalities); + /** + * Set the temperature (K), pressure (Pa), and molalities + * (gmol kg-1) of the solutes + */ + void setState_TPM(doublereal t, doublereal p, + const doublereal * const molalities); + /** Set the temperature (K), pressure (Pa), and molalities. */ + void setState_TPM(doublereal t, doublereal p, compositionMap& m); - private: - void initLengths(); + /** Set the temperature (K), pressure (Pa), and molalities. */ + void setState_TPM(doublereal t, doublereal p, const string& m); + + private: + void initLengths(); - protected: + protected: - int m_indexSolvent; - doublereal m_weightSolvent; - /* - * In any molality implementation, it makes sense to have - * a minimum solvent mole fraction requirement, since the - * implementation becomes singular in the xmolSolvent=0 - * limit. The default is to set it to 0.01. - * We then modify the molality definition to ensure that - * molal_solvent = 0 when xmol_solvent = 0. - */ - doublereal m_xmolSolventMIN; - /* - * This is the multiplication factor that goes inside - * log expressions involving the molalities of species. - * Its equal to Wt_0 / 1000. - * where Wt_0 = weight of solvent (kg/kmol) - */ - doublereal m_Mnaught; + int m_indexSolvent; + doublereal m_weightSolvent; + /* + * In any molality implementation, it makes sense to have + * a minimum solvent mole fraction requirement, since the + * implementation becomes singular in the xmolSolvent=0 + * limit. The default is to set it to 0.01. + * We then modify the molality definition to ensure that + * molal_solvent = 0 when xmol_solvent = 0. + */ + doublereal m_xmolSolventMIN; + /* + * This is the multiplication factor that goes inside + * log expressions involving the molalities of species. + * Its equal to Wt_0 / 1000. + * where Wt_0 = weight of solvent (kg/kmol) + */ + doublereal m_Mnaught; - mutable vector_fp m_molalities; - private: - doublereal err(string msg) const; + mutable vector_fp m_molalities; + private: + doublereal err(string msg) const; - }; + }; }