diff --git a/Cantera/src/thermo/DebyeHuckel.cpp b/Cantera/src/thermo/DebyeHuckel.cpp new file mode 100644 index 000000000..f655b2aba --- /dev/null +++ b/Cantera/src/thermo/DebyeHuckel.cpp @@ -0,0 +1,2217 @@ +/** + * @file DebyeHuckel.cpp + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#ifndef MAX +#define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif + +#include "DebyeHuckel.h" +#include "importCTML.h" + + +namespace Cantera { + + /** + * Default constructor + */ + DebyeHuckel::DebyeHuckel() : + MolalityVPSSTP(), + m_formDH(DHFORM_DILUTE_LIMIT), + m_formGC(2), + m_Pcurrent(OneAtm), + m_IionicMolality(0.0), + m_maxIionicStrength(30.0), + m_useHelgesonFixedForm(false), + m_IionicMolalityStoich(0.0), + m_A_Debye(1.172576), // units = sqrt(kg/gmol) + m_B_Debye(3.28640E9) // units = sqrt(kg/gmol) / m + { + m_npActCoeff.resize(3); + m_npActCoeff[0] = 0.1127; + m_npActCoeff[1] = -0.01049; + m_npActCoeff[2] = 1.545E-3; + } + /** + * Working constructors + * + * The two constructors below are the normal way + * the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + */ + DebyeHuckel::DebyeHuckel(string inputFile, string id) : + MolalityVPSSTP(), + m_formDH(DHFORM_DILUTE_LIMIT), + m_formGC(2), + m_Pcurrent(OneAtm), + m_IionicMolality(0.0), + m_maxIionicStrength(30.0), + m_useHelgesonFixedForm(false), + m_IionicMolalityStoich(0.0), + m_A_Debye(1.172576), // units = sqrt(kg/gmol) + m_B_Debye(3.28640E9) // units = sqrt(kg/gmol) / m + { + m_npActCoeff.resize(3); + m_npActCoeff[0] = 0.1127; + m_npActCoeff[1] = -0.01049; + m_npActCoeff[2] = 1.545E-3; + constructPhaseFile(inputFile, id); + } + + DebyeHuckel::DebyeHuckel(XML_Node& phaseRoot, string id) : + MolalityVPSSTP(), + m_formDH(DHFORM_DILUTE_LIMIT), + m_formGC(2), + m_Pcurrent(OneAtm), + m_IionicMolality(0.0), + m_maxIionicStrength(3.0), + m_useHelgesonFixedForm(false), + m_IionicMolalityStoich(0.0), + m_A_Debye(1.172576), // units = sqrt(kg/gmol) + m_B_Debye(3.28640E9) // units = sqrt(kg/gmol) / m + { + m_npActCoeff.resize(3); + m_npActCoeff[0] = 0.1127; + m_npActCoeff[1] = -0.01049; + m_npActCoeff[2] = 1.545E-3; + constructPhaseXML(phaseRoot, id); + } + + /** + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + DebyeHuckel::DebyeHuckel(const DebyeHuckel &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 + */ + DebyeHuckel& DebyeHuckel:: + operator=(const DebyeHuckel &b) { + if (&b != this) { + MolalityVPSSTP::operator=(b); + m_formDH = b.m_formDH; + m_formGC = b.m_formGC; + m_Pcurrent = b.m_Pcurrent; + m_Aionic = b.m_Aionic; + m_npActCoeff = b.m_npActCoeff; + m_IionicMolality = b.m_IionicMolality; + m_maxIionicStrength = b.m_maxIionicStrength; + m_useHelgesonFixedForm= b.m_useHelgesonFixedForm; + m_IionicMolalityStoich= b.m_IionicMolalityStoich; + m_A_Debye = b.m_A_Debye; + m_B_Debye = b.m_B_Debye; + m_B_Dot = b.m_B_Dot; + m_expg0_RT = b.m_expg0_RT; + m_pe = b.m_pe; + m_pp = b.m_pp; + m_tmpV = b.m_tmpV; + m_speciesCharge_Stoich= b.m_speciesCharge_Stoich; + m_Beta_ij = b.m_Beta_ij; + m_lnActCoeffMolal = b.m_lnActCoeffMolal; + m_d2lnActCoeffMolaldT2= b.m_d2lnActCoeffMolaldT2; + } + return *this; + } + + + /** + * ~DebyeHuckel(): (virtual) + * + * Destructor: does nothing: + */ + DebyeHuckel::~DebyeHuckel() { + } + + /** + * duplMyselfAsThermoPhase(): + * + * This routine operates at the ThermoPhase level to + * duplicate the current object. It uses the copy constructor + * defined above. + */ + ThermoPhase* DebyeHuckel::duplMyselfAsThermoPhase() { + DebyeHuckel* mtp = new DebyeHuckel(*this); + return (ThermoPhase *) mtp; + } + + /** + * 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. + */ + int DebyeHuckel::eosType() const { + int res; + switch (m_formGC) { + case 0: + res = cDebyeHuckel0; + break; + case 1: + res = cDebyeHuckel1; + break; + case 2: + res = cDebyeHuckel2; + break; + default: + throw CanteraError("eosType", "Unknown type"); + break; + } + return res; + } + + // + // -------- Molar Thermodynamic Properties of the Solution --------------- + // + /** + * Molar enthalpy of the solution. Units: J/kmol. + */ + doublereal DebyeHuckel::enthalpy_mole() const { + return err("not implemented"); + } + + /** + * Molar internal energy of the solution. Units: J/kmol. + */ + doublereal DebyeHuckel::intEnergy_mole() const { + getPartialMolarEnthalpies(DATA_PTR(m_tmpV)); + return mean_X(DATA_PTR(m_tmpV)); + } + + doublereal DebyeHuckel::entropy_mole() const { + getPartialMolarEntropies(DATA_PTR(m_tmpV)); + return mean_X(DATA_PTR(m_tmpV)); + } + + /// Molar Gibbs function. Units: J/kmol. + doublereal DebyeHuckel::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 DebyeHuckel::cp_mole() const { + getPartialMolarCp(DATA_PTR(m_tmpV)); + return mean_X(DATA_PTR(m_tmpV)); + } + + /// Molar heat capacity at constant volume. Units: J/kmol/K. + doublereal DebyeHuckel::cv_mole() const { + //getPartialMolarCv(m_tmpV.begin()); + //return mean_X(m_tmpV.begin()); + err("not implemented"); + return 0.0; + } + + // + // ------- Mechanical Equation of State Properties ------------------------ + // + + /** + * Pressure. Units: Pa. + * For this incompressible system, we return the internally storred + * independent value of the pressure. + */ + doublereal DebyeHuckel::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 DebyeHuckel::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 DebyeHuckel::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 DebyeHuckel::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 DebyeHuckel::setMolarDensity(doublereal conc) { + double concI = molarDensity(); + if (conc != concI) { + throw CanteraError("Idea;MolalSoln::setMolarDensity", + "molarDensity/density 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 DebyeHuckel::getActivityConcentrations(doublereal* c) const { + double c_solvent = standardConcentration(); + getActivities(c); + for (int k = 0; k < m_kk; k++) { + c[k] *= c_solvent; + } + } + + /** + * 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 DebyeHuckel::standardConcentration(int k) const { + double mvSolvent = m_speciesSize[m_indexSolvent]; + return 1.0 / mvSolvent; + } + + /** + * Returns the natural logarithm of the standard + * concentration of the kth species + */ + doublereal DebyeHuckel::logStandardConc(int k) const { + double c_solvent = standardConcentration(k); + return log(c_solvent); + } + + /** + * 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 DebyeHuckel::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; + } + } + + + /** + * Get the array of non-dimensional activities at + * the current solution temperature, pressure, and + * solution concentration. + * (note solvent activity coefficient is on the molar scale). + * + */ + void DebyeHuckel::getActivities(doublereal* ac) const { + /* + * Update the molality array, m_molalities() + * This requires an update due to mole fractions + */ + _updatelnMolalityActCoeff(); + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]); + } + } + double xmolSolvent = moleFraction(m_indexSolvent); + ac[m_indexSolvent] = + exp(m_lnActCoeffMolal[m_indexSolvent]) * xmolSolvent; + } + + /** + * getMolalityActivityCoefficients() (virtual, const) + * + * Get the array of non-dimensional Molality based + * activity coefficients at + * the current solution temperature, pressure, and + * solution concentration. + * (note solvent activity coefficient is on the molar scale). + * + * Note, most of the work is done in an internal private routine + */ + void DebyeHuckel:: + getMolalityActivityCoefficients(doublereal* acMolality) const { + _updatelnMolalityActCoeff(); + copy(m_lnActCoeffMolal.begin(), m_lnActCoeffMolal.end(), acMolality); + for (int k = 0; k < m_kk; k++) { + acMolality[k] = exp(acMolality[k]); + } + } + + // + // ------ 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) + * \f] + * + * \f[ + * \mu_solvent = \mu^{o}_solvent(T,P) + + * R T ((X_solvent - 1.0) / X_solvent) + * \f] + */ + void DebyeHuckel::getChemPotentials(doublereal* mu) const{ + double xx; + const double xxSmall = 1.0E-150; + /* + * First get the standard chemical potentials in + * molar form. + * -> this requires updates of standard state as a function + * of T and P + */ + getStandardChemPotentials(mu); + /* + * Update the activity coefficients + * This also updates the internal molality array. + */ + _updatelnMolalityActCoeff(); + /* + * + */ + doublereal RT = GasConstant * temperature(); + double xmolSolvent = moleFraction(m_indexSolvent); + for (int k = 0; k < m_kk; k++) { + if (m_indexSolvent != k) { + xx = MAX(m_molalities[k], xxSmall); + mu[k] += RT * (log(xx) + m_lnActCoeffMolal[k]); + } + } + xx = MAX(xmolSolvent, xxSmall); + mu[m_indexSolvent] += + RT * (log(xx) + m_lnActCoeffMolal[m_indexSolvent]); + } + + + /** + * Returns an array of partial molar enthalpies for the species + * in the mixture. + * Units (J/kmol) + * + * We calculate this quantity partially from the relation and + * partially by calling the standard state enthalpy function. + * + * hbar_i = - T**2 * d(chemPot_i/T)/dT + * + * We calculate + */ + void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const { + getEnthalpy_RT(hbar); + + /* + * Check to see whether activity coefficients are temperature + * dependent. If they are, then calculate the their temperature + * derivatives and add them into the result. + */ + double dAdT = dA_DebyedT_TP(); + if (dAdT != 0.0) { + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + _updatelnMolalityActCoeff(); + _updatedlnMolalityActCoeffdT(); + double T = temperature(); + double RTT = GasConstant * T * T; + for (int k = 0; k < m_kk; k++) { + hbar[k] -= RTT * m_dlnActCoeffMolaldT[k]; + } + } + } + + /** + * + * 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 + */ + void DebyeHuckel:: + getPartialMolarEntropies(doublereal* sbar) const { + int k; + /* + * Get the standard state entropies at the temperature + * and pressure of the solution. + */ + getEntropy_R(sbar); + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + _updatelnMolalityActCoeff(); + + doublereal R = GasConstant; + doublereal mm; + /* + * First we will add in the obvious dependence on the T + * term out front of the log activity term + */ + for (k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + mm = fmaxx(SmallNumber, m_molalities[k]); + sbar[k] -= R * (log(mm) + m_lnActCoeffMolal[k]); + } + } + double xmolSolvent = moleFraction(m_indexSolvent); + mm = fmaxx(SmallNumber, xmolSolvent); + sbar[m_indexSolvent] -= + R *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]); + /* + * Check to see whether activity coefficients are temperature + * dependent. If they are, then calculate the their temperature + * derivatives and add them into the result. + */ + double dAdT = dA_DebyedT_TP(); + if (dAdT != 0.0) { + _updatedlnMolalityActCoeffdT(); + double RT = R * temperature(); + for (k = 0; k < m_kk; k++) { + sbar[k] -= RT * m_dlnActCoeffMolaldT[k]; + } + } + } + + /** + * getPartialMolarVolumes() (virtual, const) + * + * 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. + * + * The general relation is + * + * vbar_i = d(chemPot_i)/dP at const T, n + * + * So, if the activity coefficients depended on pressure this + * function would be nontrivial. + */ + void DebyeHuckel::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 DebyeHuckel::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; + } + + /* + * Check to see whether activity coefficients are temperature + * dependent. If they are, then calculate the their temperature + * derivatives and add them into the result. + */ + double dAdT = dA_DebyedT_TP(); + if (dAdT != 0.0) { + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + _updatelnMolalityActCoeff(); + _updatedlnMolalityActCoeffdT(); + _updated2lnMolalityActCoeffdT2(); + double T = temperature(); + double RT = GasConstant * T; + double RTT = RT * T; + for (int k = 0; k < m_kk; k++) { + cpbar[k] -= (2.0 * RT * m_dlnActCoeffMolaldT[k] + + RTT * m_d2lnActCoeffMolaldT2[k]); + } + } + } + + + /* + * -------- 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 DebyeHuckel::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_speciesSize[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 DebyeHuckel::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 DebyeHuckel::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_speciesSize[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 DebyeHuckel:: + 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_speciesSize[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 DebyeHuckel:: + 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 DebyeHuckel::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 DebyeHuckel::getStandardVolumes(doublereal *vol) const { + copy(m_speciesSize.begin(), + m_speciesSize.end(), vol); + } + + + /* + * ------ Thermodynamic Values for the Species Reference States --- + */ + + // -> This is handled by VPStandardStatesTP + + /* + * -------------- Utilities ------------------------------- + */ + + /** + * Initialization routine for a DebyeHuckel phase. + * + * This is a virtual routine. This routine will call initThermo() + * for the parent class as well. + */ + void DebyeHuckel::initThermo() { + MolalityVPSSTP::initThermo(); + initLengths(); + } + + /** + * constructPhaseFile + * + * Initialization of a Debye-Huckel 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 DebyeHuckel::constructPhaseFile(string inputFile, string id) { + + if (inputFile.size() == 0) { + throw CanteraError("DebyeHuckel::initThermo", + "input file is null"); + } + string path = findInputFile(inputFile); + ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("DebyeHuckel::initThermo","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("DebyeHuckel::initThermo", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + fxml_phase->copy(&phaseNode_XML); + constructPhaseXML(*fxml_phase, id); + delete fxml; + } + + /** + * interp_est() (static) + * + * utility function to assign an integer value from a string + * for the ElectrolyteSpeciesType field. + */ + static int interp_est(string estString) { + const char *cc = estString.c_str(); + if (!strcasecmp(cc, "solvent")) { + return cEST_solvent; + } else if (!strcasecmp(cc, "chargedspecies")) { + return cEST_chargedSpecies; + } else if (!strcasecmp(cc, "weakAcidAssociated")) { + return cEST_weakAcidAssociated; + } else if (!strcasecmp(cc, "strongAcidAssociated")) { + return cEST_strongAcidAssociated; + } else if (!strcasecmp(cc, "polarNeutral")) { + return cEST_polarNeutral; + } else if (!strcasecmp(cc, "nonpolarNeutral")) { + return cEST_nonpolarNeutral; + } + int retn, rval; + if ((retn = sscanf(cc, "%d", &rval)) != 1) { + return -1; + } + return rval; + } + + /** + * Import and initialize a DebyeHuckel 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 DebyeHuckel::constructPhaseXML(XML_Node& phaseNode, string id) { + + if (id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("DebyeHuckel::constructPhaseXML", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("DebyeHuckel::constructPhaseXML", + "no thermo XML node"); + } + XML_Node& thermoNode = phaseNode.child("thermo"); + + /* + * Possibly 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; + printf("exit standardConc = unity not done\n"); + exit(-1); + } else if (formString == "molar_volume") { + m_formGC = 1; + printf("exit standardConc = molar_volume not done\n"); + exit(-1); + } else if (formString == "solvent_volume") { + m_formGC = 2; + } else { + throw CanteraError("DebyeHuckel::constructPhaseXML", + "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("DebyeHuckel::constructPhaseXML", + "badly formed solvent XML node"); + } + solventName = nameSolventa[0]; + } + + /* + * Determine the form of the Debye-Huckel model, + * m_formDH. We will use this information to size arrays below. + */ + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& scNode = thermoNode.child("activityCoefficients"); + m_formDH = DHFORM_DILUTE_LIMIT; + string formString = scNode.attrib("model"); + if (formString != "") { + if (formString == "Dilute_limit") { + m_formDH = DHFORM_DILUTE_LIMIT; + } else if (formString == "Bdot_with_variable_a") { + m_formDH = DHFORM_BDOT_AK ; + } else if (formString == "Bdot_with_common_a") { + m_formDH = DHFORM_BDOT_ACOMMON; + } else if (formString == "Beta_ij") { + m_formDH = DHFORM_BETAIJ; + } else if (formString == "Pitzer_with_Beta_ij") { + m_formDH = DHFORM_PITZER_BETAIJ; + } else { + throw CanteraError("DebyeHuckel::constructPhaseXML", + "Unknown standardConc model: " + formString); + } + } + } else { + /* + * If there is no XML node named "activityCoefficients", assume + * that we are doing the extreme dilute limit assumption + */ + m_formDH = DHFORM_DILUTE_LIMIT; + } + + /* + * Call the Cantera importPhase() function. This will import + * all of the species into the phase. This will also handle + * all of the solvent and solute standard states + */ + bool m_ok = importPhase(phaseNode, this); + if (!m_ok) { + throw CanteraError("DebyeHuckel::constructPhaseXML", + "importPhase failed "); + } + + } + + /** + * Process the XML file after species are set up. + * + * This gets called from importPhase(). It processes the XML file + * after the species are set up. This is the main routine for + * reading in activity coefficient parameters. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void DebyeHuckel:: + initThermoXML(XML_Node& phaseNode, string id) { + int k; + string stemp; + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("HMWSoln::initThermoXML", + "no thermo XML node"); + } + XML_Node& thermoNode = phaseNode.child("thermo"); + + /* + * Initialize all of the lengths of arrays in the object + * now that we know what species are in the phase. + */ + initLengths(); + + /* + * Reconcile the solvent name and index. + */ + /* + * 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("DebyeHuckel::initThermoXML", + "badly formed solvent XML node"); + } + solventName = nameSolventa[0]; + } + for (k = 0; k < m_kk; k++) { + string sname = speciesName(k); + if (solventName == sname) { + m_indexSolvent = k; + break; + } + } + if (m_indexSolvent == -1) { + cout << "DebyeHuckel::initThermoXML: Solvent Name not found" + << endl; + throw CanteraError("DebyeHuckel::initThermoXML", + "Solvent name not found"); + } + if (m_indexSolvent != 0) { + throw CanteraError("DebyeHuckel::initThermoXML", + "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 (k = 0; k < m_kk; k++) { + XML_Node* s = speciesDB->findByAttr("name", sss[k]); + XML_Node *ss = s->findByName("standardState"); + m_speciesSize[k] = getFloat(*ss, "molarVolume", "-"); +#ifdef DEBUG_HKM_NOT + cout << "species " << sss[k] << " has volume " << + m_speciesSize[k] << endl; +#endif + } + + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + XML_Node *acNodePtr = 0; + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& acNode = thermoNode.child("activityCoefficients"); + acNodePtr = &acNode; + /* + * Look for parameters for A_Debye + */ + if (acNode.hasChild("A_Debye")) { + m_A_Debye = getFloat(acNode, "A_Debye"); +#ifdef DEBUG_HKM_NOT + cout << "A_Debye = " << m_A_Debye << endl; +#endif + } + + /* + * Look for parameters for B_Debye + */ + if (acNode.hasChild("B_Debye")) { + m_B_Debye = getFloat(acNode, "B_Debye"); +#ifdef DEBUG_HKM_NOT + cout << "B_Debye = " << m_B_Debye << endl; +#endif + } + + /* + * Look for parameters for B_dot + */ + if (acNode.hasChild("B_dot")) { + if (m_formDH == DHFORM_BETAIJ || + m_formDH == DHFORM_DILUTE_LIMIT || + m_formDH == DHFORM_PITZER_BETAIJ) { + throw CanteraError("DebyeHuckel:init", + "B_dot entry in the wrong DH form"); + } + double bdot_common = getFloat(acNode, "B_dot"); +#ifdef DEBUG_HKM_NOT + cout << "B_dot = " << bdot_common << endl; +#endif + /* + * Set B_dot parameters for charged species + */ + for (int k = 0; k < m_kk; k++) { + double z_k = charge(k); + if (fabs (z_k) > 0.0001) { + m_B_Dot[k] = bdot_common; + } else { + m_B_Dot[k] = 0.0; + } + } + } + + /* + * Look for Parameters for the Maximum Ionic Strength + */ + if (acNode.hasChild("maxIonicStrength")) { + m_maxIionicStrength = getFloat(acNode, "maxIonicStrength"); +#ifdef DEBUG_HKM_NOT + cout << "m_maxIionicStrength = " + < That's how the activity coefficient is + * parameterized. In this case only do we allow the + * code to read in these parameters. + */ + if (m_formDH == DHFORM_BDOT_AK) { + /* + * Define a string-string map, and interpret the + * value of the xml element as binary pairs separated + * by colons, e.g.: + * Na+:3.0 + * Cl-:4.0 + * H+:9.0 + * OH-:3.5 + * Read them into the map. + */ + map m; + getMap(irNode, m); + /* + * Iterate over the map pairs, interpreting the + * first string as a species in the current phase. + * If no match is made, silently ignore the + * lack of agreement (HKM -> may be changed in the + * future). + */ + map::const_iterator _b = m.begin(); + for (; _b != m.end(); ++_b) { + int kk = speciesIndex(_b->first); + if (kk < 0) { + //throw CanteraError( + // "DebyeHuckel::initThermoXML error", + // "no species match was found" + // ); + } else { + m_Aionic[kk] = fpValue(_b->second) * Afactor; + } + } + } + } + /* + * Get the matrix of coefficients for the Beta + * binary interaction parameters. We assume here that + * this matrix is symmetric, so that we only have to + * input 1/2 of the values. + */ + if (acNode.hasChild("DHBetaMatrix")) { + if (m_formDH == DHFORM_BETAIJ || + m_formDH == DHFORM_PITZER_BETAIJ) { + XML_Node& irNode = acNode.child("DHBetaMatrix"); + const vector& sn = speciesNames(); + getMatrixValues(irNode, sn, sn, m_Beta_ij, true, true); + } else { + throw CanteraError("DebyeHuckel::initThermoXML:", + "DHBetaMatrix found for wrong type"); + } + } + + /* + * Fill in parameters for the calculation of the + * stoichiometric Ionic Strength + * + * The default is that stoich charge is the same as the + * regular charge. + */ + m_speciesCharge_Stoich.resize(m_kk, 0.0); + for (k = 0; k < m_kk; k++) { + m_speciesCharge_Stoich[k] = m_speciesCharge[k]; + } + /* + * First look at the species database. + * -> Look for the subelement "stoichIsMods" + * in each of the species SS databases. + */ + const XML_Node *phaseSpecies = speciesData(); + if (phaseSpecies) { + string kname, jname; + vector xspecies; + phaseSpecies->getChildren("species",xspecies); + int jj = xspecies.size(); + for (k = 0; k < m_kk; k++) { + int jmap = -1; + kname = speciesName(k); + for (int j = 0; j < jj; j++) { + const XML_Node& sp = *xspecies[j]; + jname = sp["name"]; + if (jname == kname) { + jmap = j; + break; + } + } + if (jmap > -1) { + const XML_Node& sp = *xspecies[jmap]; + if (sp.hasChild("stoichIsMods")) { + double val = getFloat(sp, "stoichIsMods"); + m_speciesCharge_Stoich[k] = val; + } + } + } + } + /* + * Now look at the activity coefficient database + */ + if (acNodePtr) { + if (acNodePtr->hasChild("stoichIsMods")) { + XML_Node& sIsNode = acNodePtr->child("stoichIsMods"); + + map msIs; + getMap(sIsNode, msIs); + map::const_iterator _b = msIs.begin(); + for (; _b != msIs.end(); ++_b) { + int kk = speciesIndex(_b->first); + if (kk < 0) { + //throw CanteraError( + // "DebyeHuckel::initThermoXML error", + // "no species match was found" + // ); + } else { + double val = fpValue(_b->second); + m_speciesCharge_Stoich[kk] = val; + } + } + } + } + } + + /* + * Fill in the vector specifying the electrolyte species + * type + * + * First fill in default values. Everthing is either + * a charge species, a nonpolar neutral, or the solvent. + */ + for (k = 0; k < m_kk; k++) { + if (fabs(m_speciesCharge[k]) > 0.0001) { + m_electrolyteSpeciesType[k] = cEST_chargedSpecies; + if (fabs(m_speciesCharge_Stoich[k] - m_speciesCharge[k]) + > 0.0001) { + m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated; + } + } else if (fabs(m_speciesCharge_Stoich[k]) > 0.0001) { + m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated; + } else { + m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral; + } + } + m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent; + /* + * First look at the species database. + * -> Look for the subelement "stoichIsMods" + * in each of the species SS databases. + */ + const XML_Node *phaseSpecies = speciesData(); + const XML_Node *spPtr = 0; + if (phaseSpecies) { + string kname; + for (k = 0; k < m_kk; k++) { + kname = speciesName(k); + spPtr = speciesXML_Node(kname, phaseSpecies); + if (!spPtr) { + if (spPtr->hasChild("electrolyteSpeciesType")) { + string est = getString(*spPtr, "electrolyteSpeciesType"); + if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) { + throw CanteraError("DebyeHuckel:initThermoXML", + "Bad electrolyte type: " + est); + } + } + } + } + } + /* + * Then look at the phase thermo specification + */ + if (acNodePtr) { + if (acNodePtr->hasChild("electrolyteSpeciesType")) { + XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType"); + map msEST; + getMap(ESTNode, msEST); + map::const_iterator _b = msEST.begin(); + for (; _b != msEST.end(); ++_b) { + int kk = speciesIndex(_b->first); + if (kk < 0) { + } else { + string est = _b->second; + if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) { + throw CanteraError("DebyeHuckel:initThermoXML", + "Bad electrolyte type: " + est); + } + } + } + } + } + + /* + * Lastly 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 DebyeHuckel::setParameters(int n, doublereal* c) { + } + void DebyeHuckel::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 (initThermoXML) + * It just didn't seem to fit. + */ + void DebyeHuckel::setParametersFromXML(const XML_Node& eosdata) { + } + + /** + * Report the molar volume of species k + * + * units - \f$ m^3 kmol^-1 \f$ + */ + double DebyeHuckel::speciesMolarVolume(int k) const { + return m_speciesSize[k]; + } + + + /** + * A_Debye_TP() (virtual) + * + * Returns the A_Debye parameter as a function of temperature + * and pressure. + * + * The default is to assume that it is constant, given + * in the initialization process and storred in the + * member double, m_A_Debye + */ + double DebyeHuckel::A_Debye_TP(double tempArg, double presArg) const { + double T = temperature(); + if (tempArg != -1.0) { + T = tempArg; + } + return m_A_Debye; + } + + /** + * dA_DebyedT_TP() (virtual) + * + * Returns the derivative of the A_Debye parameter with + * respect to temperature as a function of temperature + * and pressure. + * + * The default is to assume that it is equal to zero + * -> note, placeholder until a better formalism is + * put in place. + */ + double DebyeHuckel::dA_DebyedT_TP(double tempArg, double presArg) const { + double T = temperature(); + if (tempArg != -1.0) { + T = tempArg; + } + return 0.0; + } + + /** + * d2A_DebyedT2_TP() (virtual) + * + * Returns the 2nd derivative of the A_Debye parameter with + * respect to temperature as a function of temperature + * and pressure. + * + * The default is to assume that it is equal to zero + * -> note, placeholder until a better formalism is + * put in place. + */ + double DebyeHuckel::d2A_DebyedT2_TP(double tempArg, double presArg) const { + double T = temperature(); + if (tempArg != -1.0) { + T = tempArg; + } + return 0.0; + } + + /* + * ----------- Critical State Properties -------------------------- + */ + + /* + * ---------- Other Property Functions + */ + double DebyeHuckel::AionicRadius(int k) const { + return m_Aionic[k]; + } + + /* + * ------------ Private and Restricted Functions ------------------ + */ + + /** + * Bail out of functions with an error exit if they are not + * implemented. + */ + doublereal DebyeHuckel::err(string msg) const { + throw CanteraError("DebyeHuckel", + "Unfinished func called: " + msg ); + return 0.0; + } + + + /** + * initLengths(): + * + * This internal function adjusts the lengths of arrays based on + * the number of species + */ + void DebyeHuckel::initLengths() { + m_kk = nSpecies(); + MolalityVPSSTP::initThermo(); + + /* + * Obtain the limits of the temperature from the species + * thermo handler's limits. + */ + int leng = m_kk; + m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral); + m_speciesSize.resize(leng); + m_Aionic.resize(leng, 0.0); + m_lnActCoeffMolal.resize(leng, 0.0); + m_dlnActCoeffMolaldT.resize(leng, 0.0); + m_d2lnActCoeffMolaldT2.resize(leng, 0.0); + m_B_Dot.resize(leng, 0.0); + m_expg0_RT.resize(leng, 0.0); + m_pe.resize(leng, 0.0); + m_pp.resize(leng, 0.0); + m_tmpV.resize(leng, 0.0); + if (m_formDH == DHFORM_BETAIJ || + m_formDH == DHFORM_PITZER_BETAIJ) { + m_Beta_ij.resize(leng, leng, 0.0); + } + } + + /** + * nonpolarActCoeff() (private) + * + * Static function that implements the non-polar species + * salt-out modifications. + * Returns the calculated activity coefficients. + */ + double DebyeHuckel::_nonpolarActCoeff(double IionicMolality) const { + double I2 = IionicMolality * IionicMolality; + double l10actCoeff = + m_npActCoeff[0] * IionicMolality + + m_npActCoeff[1] * I2 + + m_npActCoeff[2] * I2 * IionicMolality; + return pow(10.0 , l10actCoeff); + } + + + /** + * _osmoticCoeffHelgesonFixedForm() + * + * Formula for the osmotic coefficient that occurs in + * the GWB. It is originally from Helgeson for a variable + * NaCl brine. It's to be used with extreme caution. + */ + double DebyeHuckel:: + _osmoticCoeffHelgesonFixedForm() const { + const double a0 = 1.454; + const double b0 = 0.02236; + const double c0 = 9.380E-3; + const double d0 = -5.362E-4; + double Is = m_IionicMolalityStoich; + if (Is <= 0.0) { + return 0.0; + } + double Is2 = Is * Is; + double bhat = 1.0 + a0 * sqrt(Is); + double func = bhat - 2.0 * log(bhat) - 1.0/bhat; + double v1 = m_A_Debye / (a0 * a0 * a0 * Is) * func; + double oc = 1.0 - v1 + b0 * Is / 2.0 + 2.0 * c0 * Is2 / 3.0 + + 3.0 * d0 * Is2 * Is / 4.0; + return oc; + } + + + /** + * _activityWaterHelgesonFixedForm() + * + * Formula for the log of the activity of the water + * solvent that occurs in + * the GWB. It is originally from Helgeson for a variable + * NaCl brine. It's to be used with extreme caution. + */ + double DebyeHuckel:: + _lnactivityWaterHelgesonFixedForm() const { + /* + * Update the internally storred vector of molalities + */ + calcMolalities(); + double oc = _osmoticCoeffHelgesonFixedForm(); + double sum = 0.0; + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + sum += MAX(m_molalities[k], 0.0); + } + } + if (sum > 2.0 * m_maxIionicStrength) { + sum = 2.0 * m_maxIionicStrength; + }; + double lac = - m_Mnaught * sum * oc; + return lac; + } + + /** + * _updatelnMolalityActCoeff(): + * + * Using internally stored values, this function calculates + * the activity coefficients for all species. + * + * The ln(activity_solvent) is first calculated for the + * solvent. Then the molar based activity coefficient + * is calculated and returned. + * + * ( Note this is the main routine for implementing the + * activity coefficient formulation.) + */ + void DebyeHuckel::_updatelnMolalityActCoeff() const { + double z_k, zs_k1, zs_k2; + /* + * Update the internally storred vector of molalities + */ + calcMolalities(); + /* + * Calculate the apparent (real) ionic strength. + * + * Note this is not the stoichiometric ionic strengh, + * where reactions of ions forming neutral salts + * are ignorred in calculating the ionic strength. + */ + m_IionicMolality = 0.0; + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + m_IionicMolality += m_molalities[k] * z_k * z_k; + } + m_IionicMolality /= 2.0; + + if (m_IionicMolality > m_maxIionicStrength) { + m_IionicMolality = m_maxIionicStrength; + } + + /* + * Calculate the stoichiometric ionic charge + */ + m_IionicMolalityStoich = 0.0; + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + zs_k1 = m_speciesCharge_Stoich[k]; + if (z_k == zs_k1) { + m_IionicMolalityStoich += m_molalities[k] * z_k * z_k; + } else { + zs_k2 = z_k - zs_k1; + m_IionicMolalityStoich + += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2); + } + } + m_IionicMolalityStoich /= 2.0; + + if (m_IionicMolalityStoich > m_maxIionicStrength) { + m_IionicMolalityStoich = m_maxIionicStrength; + } + + /** + * Possibly update the storred value of the + * Debye-Huckel parameter A_Debye + * This parameter appears on the top of the activity + * coefficient expression. + * It depends on T (and P), as it depends explicity + * on the temperature. Also, the dielectric constant + * is usually a fairly strong function of T, also. + */ + m_A_Debye = A_Debye_TP(); + + /* + * Calculate a safe value for the mole fraction + * of the solvent + */ + double xmolSolvent = moleFraction(m_indexSolvent); + xmolSolvent = MAX(8.689E-3, xmolSolvent); + + int est; + double ac_nonPolar = 1.0; + double numTmp = m_A_Debye * sqrt(m_IionicMolality); + double denomTmp = m_B_Debye * sqrt(m_IionicMolality); + double coeff; + double lnActivitySolvent = 0.0; + double tmp; + double tmpLn; + double y, yp1, sigma; + switch (m_formDH) { + case DHFORM_DILUTE_LIMIT: + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + m_lnActCoeffMolal[k] = - z_k * z_k * numTmp; + } + lnActivitySolvent = + (xmolSolvent - 1.0)/xmolSolvent + + 2.0 / 3.0 * m_A_Debye * m_Mnaught * + m_IionicMolality * sqrt(m_IionicMolality); + break; + + case DHFORM_BDOT_AK: + ac_nonPolar = _nonpolarActCoeff(m_IionicMolality); + for (int k = 0; k < m_kk; k++) { + est = m_electrolyteSpeciesType[k]; + if (est == cEST_nonpolarNeutral) { + m_lnActCoeffMolal[k] = log(ac_nonPolar); + } else { + z_k = m_speciesCharge[k]; + m_lnActCoeffMolal[k] = + - z_k * z_k * numTmp / (1.0 + denomTmp * m_Aionic[k]) + + log(10.0) * m_B_Dot[k] * m_IionicMolality; + } + } + + lnActivitySolvent = (xmolSolvent - 1.0)/xmolSolvent; + coeff = 2.0 / 3.0 * m_A_Debye * m_Mnaught + * sqrt(m_IionicMolality); + tmp = 0.0; + if (denomTmp > 0.0) { + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent || m_Aionic[k] != 0.0) { + y = denomTmp * m_Aionic[k]; + yp1 = y + 1.0; + sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1)); + z_k = m_speciesCharge[k]; + tmp += m_molalities[k] * z_k * z_k * sigma / 2.0; + } + } + } + lnActivitySolvent += coeff * tmp; + tmp = 0.0; + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + if ((k != m_indexSolvent) && (z_k != 0.0)) { + tmp += m_B_Dot[k] * m_molalities[k]; + } + } + lnActivitySolvent -= + m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0; + + /* + * Special section to implement the Helgeson fixed form + * for the water brine activity coefficient. + */ + if (m_useHelgesonFixedForm) { + lnActivitySolvent = _lnactivityWaterHelgesonFixedForm(); + } + break; + + case DHFORM_BDOT_ACOMMON: + denomTmp *= m_Aionic[0]; + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + m_lnActCoeffMolal[k] = + - z_k * z_k * numTmp / (1.0 + denomTmp) + + log(10.0) * m_B_Dot[k] * m_IionicMolality; + } + if (denomTmp > 0.0) { + y = denomTmp; + yp1 = y + 1.0; + sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1)); + } else { + sigma = 0.0; + } + lnActivitySolvent = + (xmolSolvent - 1.0)/xmolSolvent + + 2.0 /3.0 * m_A_Debye * m_Mnaught * + m_IionicMolality * sqrt(m_IionicMolality) * sigma; + tmp = 0.0; + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + if ((k != m_indexSolvent) && (z_k != 0.0)) { + tmp += m_B_Dot[k] * m_molalities[k]; + } + } + lnActivitySolvent -= + m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0; + + break; + + case DHFORM_BETAIJ: + denomTmp = m_B_Debye * m_Aionic[0]; + denomTmp *= sqrt(m_IionicMolality); + lnActivitySolvent = + (xmolSolvent - 1.0)/xmolSolvent; + + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + z_k = m_speciesCharge[k]; + m_lnActCoeffMolal[k] = + - z_k * z_k * numTmp / (1.0 + denomTmp); + for (int j = 0; j < m_kk; j++) { + double beta = m_Beta_ij.value(k, j); +#ifdef DEBUG_HKM_NOT + if (beta != 0.0) { + printf("b: k = %d, j = %d, betakj = %g\n", + k, j, beta); + } +#endif + m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * beta; + } + } + } + if (denomTmp > 0.0) { + y = denomTmp; + yp1 = y + 1.0; + sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 -2.0*log(yp1)); + } else { + sigma = 0.0; + } + lnActivitySolvent = + (xmolSolvent - 1.0)/xmolSolvent + + 2.0 /3.0 * m_A_Debye * m_Mnaught * + m_IionicMolality * sqrt(m_IionicMolality) * sigma; + tmp = 0.0; + for (int k = 0; k < m_kk; k++) { + for (int j = 0; j < m_kk; j++) { + tmp += + m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j]; + } + } + lnActivitySolvent -= m_Mnaught * tmp; + break; + + case DHFORM_PITZER_BETAIJ: + denomTmp = m_B_Debye * sqrt(m_IionicMolality); + denomTmp *= m_Aionic[0]; + numTmp = m_A_Debye * sqrt(m_IionicMolality); + tmpLn = log(1.0 + denomTmp); + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + z_k = m_speciesCharge[k]; + m_lnActCoeffMolal[k] = + - z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp); + m_lnActCoeffMolal[k] += + - 2.0 * z_k * z_k * m_A_Debye * tmpLn / + (3.0 * m_B_Debye * m_Aionic[0]); + for (int j = 0; j < m_kk; j++) { + m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * + m_Beta_ij.value(k, j); + } + } + } + sigma = 1.0 / (1.0 + denomTmp); + lnActivitySolvent = + (xmolSolvent - 1.0)/xmolSolvent + + 2.0 /3.0 * m_A_Debye * m_Mnaught * + m_IionicMolality * sqrt(m_IionicMolality) * sigma; + tmp = 0.0; + for (int k = 0; k < m_kk; k++) { + for (int j = 0; j < m_kk; j++) { + tmp += + m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j]; + } + } + lnActivitySolvent -= m_Mnaught * tmp; + break; + + default: + printf("ERROR\n"); + exit(-1); + } + /* + * Above, we calculated the ln(activitySolvent). Translate that + * into the molar-based activity coefficient by dividing by + * the solvent mole fraction. Solvents are not on the molality + * scale. + */ + xmolSolvent = moleFraction(m_indexSolvent); + m_lnActCoeffMolal[m_indexSolvent] = + lnActivitySolvent - log(xmolSolvent); + } + + /** + * _updatedMolalityActCoeffdT() (private, const ) + * + * Using internally stored values, this function calculates + * the temperature derivative of the logarithm of the + * activity coefficient for all species in the mechanism. + * + * We assume that the activity coefficients are current. + * + * solvent activity coefficient is on the molality + * scale. It's derivative is too. + */ + void DebyeHuckel::_updatedlnMolalityActCoeffdT() const { + double z_k, coeff, tmp, y, yp1, sigma, tmpLn; + int k; + double dAdT = dA_DebyedT_TP(); + if (dAdT == 0.0) { + for (k = 0; k < m_kk; k++) { + m_dlnActCoeffMolaldT[k] = 0.0; + } + return; + } + /* + * Calculate a safe value for the mole fraction + * of the solvent + */ + double xmolSolvent = moleFraction(m_indexSolvent); + xmolSolvent = MAX(8.689E-3, xmolSolvent); + + + double sqrtI = sqrt(m_IionicMolality); + double numdAdTTmp = dAdT * sqrtI; + double denomTmp = m_B_Debye * sqrtI; + + switch (m_formDH) { + case DHFORM_DILUTE_LIMIT: + for (int k = 0; k < m_kk; k++) { + m_dlnActCoeffMolaldT[k] = + m_lnActCoeffMolal[k] * dAdT / m_A_Debye; + } + break; + + case DHFORM_BDOT_AK: + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + m_dlnActCoeffMolaldT[k] = + - z_k * z_k * numdAdTTmp / (1.0 + denomTmp * m_Aionic[k]); + } + + m_dlnActCoeffMolaldT[m_indexSolvent] = 0.0; + + coeff = 2.0 / 3.0 * dAdT * m_Mnaught * sqrtI; + tmp = 0.0; + if (denomTmp > 0.0) { + for (int k = 0; k < m_kk; k++) { + y = denomTmp * m_Aionic[k]; + yp1 = y + 1.0; + sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1)); + z_k = m_speciesCharge[k]; + tmp += m_molalities[k] * z_k * z_k * sigma / 2.0; + } + } + m_dlnActCoeffMolaldT[m_indexSolvent] += coeff * tmp; + break; + + case DHFORM_BDOT_ACOMMON: + denomTmp *= m_Aionic[0]; + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + m_dlnActCoeffMolaldT[k] = + - z_k * z_k * numdAdTTmp / (1.0 + denomTmp); + } + if (denomTmp > 0.0) { + y = denomTmp; + yp1 = y + 1.0; + sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1)); + } else { + sigma = 0.0; + } + m_dlnActCoeffMolaldT[m_indexSolvent] = + 2.0 /3.0 * dAdT * m_Mnaught * m_IionicMolality * sqrtI * sigma; + break; + + case DHFORM_BETAIJ: + denomTmp *= m_Aionic[0]; + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + z_k = m_speciesCharge[k]; + m_dlnActCoeffMolaldT[k] = + - z_k * z_k * numdAdTTmp / (1.0 + denomTmp); + } + } + if (denomTmp > 0.0) { + y = denomTmp; + yp1 = y + 1.0; + sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1)); + } else { + sigma = 0.0; + } + m_dlnActCoeffMolaldT[m_indexSolvent] = + (xmolSolvent - 1.0)/xmolSolvent + + 2.0 /3.0 * dAdT * m_Mnaught * + m_IionicMolality * sqrtI * sigma; + break; + + case DHFORM_PITZER_BETAIJ: + denomTmp *= m_Aionic[0]; + tmpLn = log(1.0 + denomTmp); + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + z_k = m_speciesCharge[k]; + m_dlnActCoeffMolaldT[k] = + - z_k * z_k * numdAdTTmp / (1.0 + denomTmp) + - 2.0 * z_k * z_k * dAdT * tmpLn + / (m_B_Debye * m_Aionic[0]); + m_dlnActCoeffMolaldT[k] /= 3.0; + } + } + + sigma = 1.0 / ( 1.0 + denomTmp); + m_dlnActCoeffMolaldT[m_indexSolvent] = + (xmolSolvent - 1.0)/xmolSolvent + + 2.0 /3.0 * dAdT * m_Mnaught * + m_IionicMolality * sqrtI * sigma; + break; + + default: + printf("ERROR\n"); + exit(-1); + break; + } + } + + /** + * _updated2lnMolalityActCoeffdT2() (private, const ) + * + * Using internally stored values, this function calculates + * the temperature 2nd derivative of the logarithm of the + * activity coefficient + * for all species in the mechanism. + * + * We assume that the activity coefficients are current. + * + * solvent activity coefficient is on the molality + * scale. It's derivatives are too. + */ + + void DebyeHuckel::_updated2lnMolalityActCoeffdT2() const { + double z_k, coeff, tmp, y, yp1, sigma, tmpLn; + int k; + double dAdT = dA_DebyedT_TP(); + if (dAdT == 0.0) { + for (k = 0; k < m_kk; k++) { + m_dlnActCoeffMolaldT[k] = 0.0; + } + return; + } + double d2AdT2 = d2A_DebyedT2_TP(); + + /* + * Calculate a safe value for the mole fraction + * of the solvent + */ + double xmolSolvent = moleFraction(m_indexSolvent); + xmolSolvent = MAX(8.689E-3, xmolSolvent); + + + double sqrtI = sqrt(m_IionicMolality); + double numd2AdT2Tmp = d2AdT2 * sqrtI; + double denomTmp = m_B_Debye * sqrtI; + + switch (m_formDH) { + case DHFORM_DILUTE_LIMIT: + for (int k = 0; k < m_kk; k++) { + m_d2lnActCoeffMolaldT2[k] = + m_lnActCoeffMolal[k] * d2AdT2 / m_A_Debye; + } + break; + + case DHFORM_BDOT_AK: + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + m_d2lnActCoeffMolaldT2[k] = + - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp * m_Aionic[k]); + } + + m_d2lnActCoeffMolaldT2[m_indexSolvent] = 0.0; + + coeff = 2.0 / 3.0 * d2AdT2 * m_Mnaught * sqrtI; + tmp = 0.0; + if (denomTmp > 0.0) { + for (int k = 0; k < m_kk; k++) { + y = denomTmp * m_Aionic[k]; + yp1 = y + 1.0; + sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1)); + z_k = m_speciesCharge[k]; + tmp += m_molalities[k] * z_k * z_k * sigma / 2.0; + } + } + m_d2lnActCoeffMolaldT2[m_indexSolvent] += coeff * tmp; + break; + + case DHFORM_BDOT_ACOMMON: + denomTmp *= m_Aionic[0]; + for (int k = 0; k < m_kk; k++) { + z_k = m_speciesCharge[k]; + m_d2lnActCoeffMolaldT2[k] = + - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp); + } + if (denomTmp > 0.0) { + y = denomTmp; + yp1 = y + 1.0; + sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1)); + } else { + sigma = 0.0; + } + m_d2lnActCoeffMolaldT2[m_indexSolvent] = + 2.0 /3.0 * d2AdT2 * m_Mnaught * + m_IionicMolality * sqrtI * sigma; + break; + + case DHFORM_BETAIJ: + denomTmp *= m_Aionic[0]; + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + z_k = m_speciesCharge[k]; + m_d2lnActCoeffMolaldT2[k] = + - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp); + } + } + if (denomTmp > 0.0) { + y = denomTmp; + yp1 = y + 1.0; + sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 -2.0*log(yp1)); + } else { + sigma = 0.0; + } + m_d2lnActCoeffMolaldT2[m_indexSolvent] = + (xmolSolvent - 1.0)/xmolSolvent + + 2.0 /3.0 * d2AdT2 * m_Mnaught * + m_IionicMolality * sqrtI * sigma; + break; + + case DHFORM_PITZER_BETAIJ: + denomTmp *= m_Aionic[0]; + tmpLn = log(1.0 + denomTmp); + for (int k = 0; k < m_kk; k++) { + if (k != m_indexSolvent) { + z_k = m_speciesCharge[k]; + m_dlnActCoeffMolaldT[k] = + - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp) + - 2.0 * z_k * z_k * d2AdT2 * tmpLn + / (m_B_Debye * m_Aionic[0]); + m_dlnActCoeffMolaldT[k] /= 3.0; + } + } + + sigma = 1.0 / ( 1.0 + denomTmp); + m_dlnActCoeffMolaldT[m_indexSolvent] = + (xmolSolvent - 1.0)/xmolSolvent + + 2.0 /3.0 * d2AdT2 * m_Mnaught * + m_IionicMolality * sqrtI * sigma; + break; + + default: + printf("ERROR\n"); + exit(-1); + break; + } + } + + +} + + diff --git a/Cantera/src/thermo/DebyeHuckel.h b/Cantera/src/thermo/DebyeHuckel.h new file mode 100644 index 000000000..b633bcfe9 --- /dev/null +++ b/Cantera/src/thermo/DebyeHuckel.h @@ -0,0 +1,1130 @@ +/** + * @file DebyeHuckel.h + * + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id$ + */ + +#ifndef CT_DEBYEHUCKEL_H +#define CT_DEBYEHUCKEL_H + +#include "MolalityVPSSTP.h" +#include "electrolytes.h" +#include "Array.h" + +namespace Cantera { + + /** + * @defgroup thermoprops Thermodynamic Properties + * + * These classes are used to compute thermodynamic properties. + */ + + /** + * DebyeHuckel.h + * + * Major Parameters: + * + * m_formDH = Form of the Debye-Huckel expression + * + * DHFORM_DILUTE_LIMIT = 0 + * + * This form assumes a dilute limit to DH, and is mainly + * for informational purposes: + * + * ln(gamma_k)/RT = -z_k**2 * alpha * sqrt(I) + * + * where I = 1/2 sum_k( molality_k * z_k**2) + * + * DHFORM_BDOT_AK = 1 + * + * This form assumes Bethke's format for the DH coefficient + * + * ln(gamma_k)/RT = -z_k**2 * alpha * sqrt(I) / (1 + B * a_k * sqrt(I)) + * + bdot_k * I + * + * (note, this particular form where a_k can differ in + * multielectrolyte + * solutions has problems wrt a gibbs-duhem analysis. However + * we include it here because there is a lot of data fit to it) + * + * DHFORM_BDOT_AUNIFORM = 2 + * + * This form assumes Bethke's format for the DH coefficient + * + * ln(gamma_k)/RT = -z_k**2 * alpha * sqrt(I) / (1 + B * a * sqrt(I)) + * + bdot_k * I + * + * The value of a is determined at the beginning of the + * calculation, and not changed. + * + * DHFORM_BETAIJ = 3 + * + * This form assumes a linear expansion in a virial coefficient form + * It is used extensively in Newmann's book, and is the beginning of + * more complex treatments for stronger electrolytes, like Pitzer + * and HMW treatments. + * + * ln(gamma_k)/RT = -z_k**2 * alpha * sqrt(I) / (1 + B * a * sqrt(I)) + * + 2* sum_j (beta_jk m_j) + * + * DHFORM_PITZER_BETAIJ = 4 + * + * This form assumes an activity coefficient formulation consistent + * with a truncated form of Pitzer's formulation. + * + * ln(gamma_k)/RT = -z_k**2 * alpha * sqrt(I) / (1 + B * a * sqrt(I)) + * -2 * z_k**2 * alpha * ln(1 + B * a * sqrt(I)) / (B * a) + * + 2 * sum_j (beta_jk m_j) + * + */ +#define DHFORM_DILUTE_LIMIT 0 +#define DHFORM_BDOT_AK 1 +#define DHFORM_BDOT_ACOMMON 2 +#define DHFORM_BETAIJ 3 +#define DHFORM_PITZER_BETAIJ 4 + + + /** + * Definition of the DebyeHuckel object + */ + class DebyeHuckel : public MolalityVPSSTP { + + public: + + /// Constructors + DebyeHuckel(); + DebyeHuckel(const DebyeHuckel &); + DebyeHuckel& operator=(const DebyeHuckel&); + + DebyeHuckel(string inputFile, string id = ""); + DebyeHuckel(XML_Node& phaseRef, string id = ""); + + /// Destructor. + virtual ~DebyeHuckel(); + + + 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; + + /** + * @} + * @name Molar Thermodynamic Properties of the Solution -------------- + * @{ + */ + + /// Molar enthalpy. Units: J/kmol. + /** + * Molar enthalpy of the solution. Units: J/kmol. + * (HKM -> Bump up to Parent object) + */ + virtual doublereal enthalpy_mole() const; + + /// Molar internal energy. Units: J/kmol. + /** + * Molar internal energy of the solution. Units: J/kmol. + * (HKM -> Bump up to Parent object) + */ + 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 + * + * (HKM -> Bump up to Parent object) + */ + virtual doublereal entropy_mole() const; + + /// Molar Gibbs function. Units: J/kmol. + /* + * (HKM -> Bump up to Parent object) + */ + virtual doublereal gibbs_mole() const; + + /// Molar heat capacity at constant pressure. Units: J/kmol/K. + /* + * (HKM -> Bump up to Parent object) + */ + virtual doublereal cp_mole() const; + + /// Molar heat capacity at constant volume. Units: J/kmol/K. + /* + * (HKM -> Bump up to Parent object) + */ + 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 a virtual function overwritten from the State.h + * class + */ + virtual void setMolarDensity(doublereal conc); + + /** + * 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, 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 the pressure. + * Activity is assumed to be molality-based here. + * @{ + */ + + /** + * 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 molality-based 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; + + 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$. + * Activity is molality based in this object. + * 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 + * standard states + * 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; } + + + /** + * @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(); + + /* + * Initialization of a DebyeHuckel 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. + */ + virtual void constructPhaseFile(string infile, string id=""); + + /* + * Import and initialize a DebyeHuckel 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. + */ + virtual void constructPhaseXML(XML_Node& phaseNode, string 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; + + + /** + * + */ + virtual double A_Debye_TP(double temperature = -1.0, + double pressure = -1.0) const; + + virtual double dA_DebyedT_TP(double temperature = -1.0, + double pressure = -1.0) const; + + virtual double d2A_DebyedT2_TP(double temperature = -1.0, + double pressure = -1.0) const; + + /* + * AionicRadius() + * + * Reports the ionic radius of the kth species + */ + double AionicRadius(int k = 0) const; + + /** + * + * formDH(): + * + * Returns the form of the Debye-Huckel parameterization used + */ + int formDH() const { return m_formDH; } + + Array2D& get_Beta_ij() { return m_Beta_ij; } + + private: + /* Static function that implements the non-polar species + * salt-out modifications. + * Returns the calculated activity coefficients. + */ + double _nonpolarActCoeff(double IionicMolality) const; + + /** + * Formula for the osmotic coefficient that occurs in + * the GWB. It is originally from Helgeson for a variable + * NaCl brine. It's to be used with extreme caution. + */ + double _osmoticCoeffHelgesonFixedForm() const; + double _lnactivityWaterHelgesonFixedForm() const; + + //@} + + + protected: + + /** + * This is the form of the Debye-Huckel parameterization + * used in this model. + * The options are described at the top of this document, + * and in the general documentation. + * The list is repeated here: + * + * DHFORM_DILUTE_LIMIT = 0 (default) + * DHFORM_BDOT_AK = 1 + * DHFORM_BDOT_AUNIFORM = 2 + * DHFORM_BETAIJ = 3 + * DHFORM_PITZER_BETAIJ = 4 + */ + int m_formDH; + + /** + * Format for the generalized concentration: + * + * 0 = unity + * 1 = molar_volume + * 2 = solvent_volume (default) + * + * The generalized concentrations can have three different forms + * depending on the value of the member attribute m_formGC, which + * is supplied in the constructor. + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + * The value and form of the generalized concentration will affect + * reaction rate constants involving species in this phase. + * + * (HKM Note: Using option #1 may lead to spurious results and + * has been included only with warnings. The reason is that it + * molar volumes of electrolytes may often be negative. The + * molar volume of H+ is defined to be zero too. Either options + * 0 or 2 are the appropriate choice. Option 0 leads to + * bulk reaction rate constants which have units of s-1. + * Option 2 leads to bulk reaction rate constants for + * bimolecular rxns which have units of m-3 kmol-1 s-1.) + */ + int m_formGC; + + /** + * Current pressure in Pascal + */ + double m_Pcurrent; + + + vector_int m_electrolyteSpeciesType; + + /** + * Species molar volumes \f$ m^3 kmol^-1 \f$ + * -> m_speciesSize in Constituents.h + */ + //array_fp m_speciesMolarVolume; + + /** + * a_k = Size of the ionic species in the DH formulation + * units = meters + */ + array_fp m_Aionic; + + /** + * Current value of the ionic strength on the molality scale + */ + mutable double m_IionicMolality; + + /** + * Maximum value of the ionic strength allowed in the + * calculation of the activity coefficients. + */ + double m_maxIionicStrength; + + /** + * If true, then the fixed for of Helgeson's activity + * for water is used instead of the rigoruous form + * obtained from Gibbs-Duhem relation. This should be + * used with caution, and is really only included as a + * validation exercise. + */ + public: + bool m_useHelgesonFixedForm; + protected: + /** + * Stoichiometric ionic strength on the molality scale + */ + mutable double m_IionicMolalityStoich; + + /** + * A_Debye -> this expression appears on the top of the + * ln actCoeff term in the general Debye-Huckel + * expression + * It depends on temperature + * + * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) + * + * Units = sqrt(kg/gmol) + * + * Nominal value = 1.172576 sqrt(kg/gmol) + * based on: + * epsilon/epsilon_0 = 78.54 + * (water at 25C) + * epsilon_0 = 8.854187817E12 C2 N-1 m-2 + * e = 8.314472E3 kg m2 s-2 kmol-1 K-1 + * F = 9.6485309E7 C kmol-1 + * R = 8.314472E3 kg m2 s-2 kmol-1 K-1 + * T = 298.15 K + * B_Debye = 3.28640E9 sqrt(kg/gmol)/m + */ + mutable double m_A_Debye; + + /** + * B_Debye -> this expression appears on the bottom of the + * ln actCoeff term in the general Debye-Huckel + * expression + * It depends on temperature + * + * B_Bebye = F / sqrt( epsilon R T / 2 ) + * + * Units = sqrt(kg/gmol) / m + * + * Nominal value = 3.28640E9 sqrt(kg/gmol) / m + * based on: + * epsilon/epsilon_0 = 78.54 + * (water at 25C) + * epsilon_0 = 8.854187817E12 C2 N-1 m-2 + * e = 8.314472E3 kg m2 s-2 kmol-1 K-1 + * F = 9.6485309E7 C kmol-1 + * R = 8.314472E3 kg m2 s-2 kmol-1 K-1 + * T = 298.15 K + */ + double m_B_Debye; + + /** + * B_Dot -> This expression is an extension of the + * Debye-Huckel expression used in some formulations + * to extend DH to higher molalities. + * B_dot is specific to the major ionic pair. + */ + array_fp m_B_Dot; + + /** + * m_npActCoeff -> These are coefficients to describe + * the increase in activity coeff for non-polar molecules + * due to the electrolyte becoming stronger (the so-called + * salt-out effect) + */ + array_fp m_npActCoeff; + + /** + * 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; + + /** + * Stoichiometric species charge -> This is for calculations + * of the ionic strength which ignore ion-ion pairing into + * neutral molecules. The Stoichiometric species charge is the + * charge of one of the ion that would occur if the species broke + * into two charged ion pairs. + * NaCl -> m_speciesCharge_Stoich = -1; + * HSO4- -> H+ + SO42- = -2 + * -> The other charge is calculated. + * For species that aren't ion pairs, its equal to the + * m_speciesCharge[] value. + */ + vector_fp m_speciesCharge_Stoich; + + /** + * Array of 2D data used in the DHFORM_BETAIJ formulation + * Beta_ij.value(i,j) is the coefficient of the jth species + * for the specification of the chemical potential of the ith + * species. + */ + Array2D m_Beta_ij; + + /** + * Logarithm of the activity coefficients on the molality + * scale. + * mutable because we change this if the composition + * or temperature or pressure changes. + */ + mutable array_fp m_lnActCoeffMolal; + mutable array_fp m_dlnActCoeffMolaldT; + mutable array_fp m_d2lnActCoeffMolaldT2; + + private: + doublereal err(string msg) const; + + + void initLengths(); + + /* + * This function will be called to update the internally storred + * natural logarithm of the molality activity coefficients + */ + void _updatelnMolalityActCoeff() const; + + void _updatedlnMolalityActCoeffdT() const; + void _updated2lnMolalityActCoeffdT2() const; + }; + +} + +#endif + + + + + diff --git a/Cantera/src/thermo/HMWSoln_input.cpp b/Cantera/src/thermo/HMWSoln_input.cpp index e9d33ec98..b1c6d6528 100644 --- a/Cantera/src/thermo/HMWSoln_input.cpp +++ b/Cantera/src/thermo/HMWSoln_input.cpp @@ -757,8 +757,7 @@ namespace Cantera { if (!m_ok) { throw CanteraError("HMWSoln::constructPhaseXML","importPhase failed "); } - - + } /** @@ -1117,7 +1116,6 @@ namespace Cantera { } } - /* * Lastly set the state */ diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index 261fccedc..f03a5c80f 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -27,13 +27,14 @@ ELECTRO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o \ IdealSolidSolnPhase.o IdealMolalSoln.o \ WaterPropsIAPWSphi.o WaterPropsIAPWS.o WaterProps.o \ PDSS.o WaterPDSS.o WaterTP.o \ - HMWSoln.o HMWSoln_input.o + HMWSoln.o HMWSoln_input.o DebyeHuckel.o ELECTRO_H = SingleSpeciesTP.h StoichSubstanceSSTP.h \ MolalityVPSSTP.h VPStandardStateTP.h \ IdealSolidSolnPhase.h IdealMolalSoln.h \ WaterPropsIAPWSphi.h WaterPropsIAPWS.h WaterProps.h \ - PDSS.h WaterPDSS.h WaterTP.h HMWSoln.h electrolytes.h + PDSS.h WaterPDSS.h WaterTP.h HMWSoln.h electrolytes.h \ + DebyeHuckel.h endif ifeq ($(do_issp),1) ISSP_OBJ = IdealSolidSolnPhase.o