diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp new file mode 100644 index 000000000..96e8b79dc --- /dev/null +++ b/Cantera/src/thermo/HMWSoln.cpp @@ -0,0 +1,5265 @@ +/** + * @file HMWSoln.cpp + * + * Member functions of Pitzer activity coefficient implementation. + */ +/* + * 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 "HMWSoln.h" +#include "importCTML.h" +#include "WaterProps.h" +#include "WaterPDSS.h" + +namespace Cantera { + + /** + * Default constructor + */ + HMWSoln::HMWSoln() : + MolalityVPSSTP(), + m_formPitzer(PITZERFORM_BASE), + m_formPitzerTemp(PITZER_TEMP_CONSTANT), + m_formGC(2), + m_Pcurrent(OneAtm), + m_IionicMolality(0.0), + m_maxIionicStrength(100.0), + m_TempPitzerRef(298.15), + m_IionicMolalityStoich(0.0), + m_form_A_Debye(A_DEBYE_WATER), + m_A_Debye(1.172576), // units = sqrt(kg/gmol) + m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m + m_waterSS(0), + m_densWaterSS(1000.), + m_waterProps(0), + m_debugCalc(0) + { + for (int i = 0; i < 17; i++) { + elambda[i] = 0.0; + elambda1[i] = 0.0; + } + } + /** + * Working constructors + * + * The two constructors below are the normal way + * the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + */ + HMWSoln::HMWSoln(string inputFile, string id) : + MolalityVPSSTP(), + m_formPitzer(PITZERFORM_BASE), + m_formPitzerTemp(PITZER_TEMP_CONSTANT), + m_formGC(2), + m_Pcurrent(OneAtm), + m_IionicMolality(0.0), + m_maxIionicStrength(100.0), + m_TempPitzerRef(298.15), + m_IionicMolalityStoich(0.0), + m_form_A_Debye(A_DEBYE_WATER), + m_A_Debye(1.172576), // units = sqrt(kg/gmol) + m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m + m_waterSS(0), + m_densWaterSS(1000.), + m_waterProps(0), + m_debugCalc(0) + { + for (int i = 0; i < 17; i++) { + elambda[i] = 0.0; + elambda1[i] = 0.0; + } + constructPhaseFile(inputFile, id); + } + + HMWSoln::HMWSoln(XML_Node& phaseRoot, string id) : + MolalityVPSSTP(), + m_formPitzer(PITZERFORM_BASE), + m_formPitzerTemp(PITZER_TEMP_CONSTANT), + m_formGC(2), + m_Pcurrent(OneAtm), + m_IionicMolality(0.0), + m_maxIionicStrength(100.0), + m_TempPitzerRef(298.15), + m_IionicMolalityStoich(0.0), + m_form_A_Debye(A_DEBYE_WATER), + m_A_Debye(1.172576), // units = sqrt(kg/gmol) + m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m + m_waterSS(0), + m_densWaterSS(1000.), + m_waterProps(0), + m_debugCalc(0) + { + for (int i = 0; i < 17; i++) { + elambda[i] = 0.0; + elambda1[i] = 0.0; + } + constructPhaseXML(phaseRoot, id); + } + + /** + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + HMWSoln::HMWSoln(const HMWSoln &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 + */ + HMWSoln& HMWSoln:: + operator=(const HMWSoln &b) { + if (&b != this) { + MolalityVPSSTP::operator=(b); + m_formPitzer = b.m_formPitzer; + m_formPitzerTemp = b.m_formPitzerTemp; + m_formGC = b.m_formGC; + m_Pcurrent = b.m_Pcurrent; + m_Aionic = b.m_Aionic; + m_IionicMolality = b.m_IionicMolality; + m_maxIionicStrength = b.m_maxIionicStrength; + m_TempPitzerRef = b.m_TempPitzerRef; + m_IionicMolalityStoich= b.m_IionicMolalityStoich; + m_form_A_Debye = b.m_form_A_Debye; + m_A_Debye = b.m_A_Debye; + m_B_Debye = b.m_B_Debye; + if (!m_waterSS) { + m_waterSS = new WaterPDSS(this, 0); + } + m_waterSS = b.m_waterSS; + m_densWaterSS = b.m_densWaterSS; + if (!m_waterProps) { + m_waterProps = new WaterProps(*b.m_waterProps); + } else { + m_waterProps = b.m_waterProps; + } + m_waterSS = b.m_waterSS; + 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_Beta0MX_ij = b.m_Beta0MX_ij; + m_Beta0MX_ij_L = b.m_Beta0MX_ij_L; + m_Beta0MX_ij_LL = b.m_Beta0MX_ij_LL; + m_Beta0MX_ij_P = b.m_Beta0MX_ij_P; + m_Beta0MX_ij_coeff = b.m_Beta0MX_ij_coeff; + m_Beta1MX_ij = b.m_Beta1MX_ij; + m_Beta1MX_ij_L = b.m_Beta1MX_ij_L; + m_Beta1MX_ij_LL = b.m_Beta1MX_ij_LL; + m_Beta1MX_ij_P = b.m_Beta1MX_ij_P; + m_Beta1MX_ij_coeff = b.m_Beta1MX_ij_coeff; + m_Beta2MX_ij = b.m_Beta2MX_ij; + m_Beta2MX_ij_L = b.m_Beta2MX_ij_L; + m_Beta2MX_ij_LL = b.m_Beta2MX_ij_LL; + m_Beta2MX_ij_P = b.m_Beta2MX_ij_P; + m_Alpha1MX_ij = b.m_Alpha1MX_ij; + m_CphiMX_ij = b.m_CphiMX_ij; + m_CphiMX_ij_L = b.m_CphiMX_ij_L; + m_CphiMX_ij_LL = b.m_CphiMX_ij_LL; + m_CphiMX_ij_P = b.m_CphiMX_ij_P; + m_CphiMX_ij_coeff = b.m_CphiMX_ij_coeff; + m_Theta_ij = b.m_Theta_ij; + m_Theta_ij_L = b.m_Theta_ij_L; + m_Theta_ij_LL = b.m_Theta_ij_LL; + m_Theta_ij_P = b.m_Theta_ij_P; + m_Psi_ijk = b.m_Psi_ijk; + m_Psi_ijk_L = b.m_Psi_ijk_L; + m_Psi_ijk_LL = b.m_Psi_ijk_LL; + m_Psi_ijk_P = b.m_Psi_ijk_P; + m_Lambda_ij = b.m_Lambda_ij; + m_Lambda_ij_L = b.m_Lambda_ij_L; + m_Lambda_ij_LL = b.m_Lambda_ij_LL; + m_Lambda_ij_P = b.m_Lambda_ij_P; + m_lnActCoeffMolal = b.m_lnActCoeffMolal; + m_dlnActCoeffMolaldT = b.m_dlnActCoeffMolaldT; + m_d2lnActCoeffMolaldT2= b.m_d2lnActCoeffMolaldT2; + m_dlnActCoeffMolaldP = b.m_dlnActCoeffMolaldP; + + m_gfunc_IJ = b.m_gfunc_IJ; + m_hfunc_IJ = b.m_hfunc_IJ; + m_BMX_IJ = b.m_BMX_IJ; + m_BMX_IJ_L = b.m_BMX_IJ_L; + m_BMX_IJ_LL = b.m_BMX_IJ_LL; + m_BMX_IJ_P = b.m_BMX_IJ_P; + m_BprimeMX_IJ = b.m_BprimeMX_IJ; + m_BprimeMX_IJ_L = b.m_BprimeMX_IJ_L; + m_BprimeMX_IJ_LL = b.m_BprimeMX_IJ_LL; + m_BprimeMX_IJ_P = b.m_BprimeMX_IJ_P; + m_BphiMX_IJ = b.m_BphiMX_IJ; + m_BphiMX_IJ_L = b.m_BphiMX_IJ_L; + m_BphiMX_IJ_LL = b.m_BphiMX_IJ_LL; + m_BphiMX_IJ_P = b.m_BphiMX_IJ_P; + m_Phi_IJ = b.m_Phi_IJ; + m_Phi_IJ_L = b.m_Phi_IJ_L; + m_Phi_IJ_LL = b.m_Phi_IJ_LL; + m_Phi_IJ_P = b.m_Phi_IJ_P; + m_Phiprime_IJ = b.m_Phiprime_IJ; + m_PhiPhi_IJ = b.m_PhiPhi_IJ; + m_PhiPhi_IJ_L = b.m_PhiPhi_IJ_L; + m_PhiPhi_IJ_LL = b.m_PhiPhi_IJ_LL; + m_PhiPhi_IJ_P = b.m_PhiPhi_IJ_P; + m_CMX_IJ = b.m_CMX_IJ; + m_CMX_IJ_L = b.m_CMX_IJ_L; + m_CMX_IJ_LL = b.m_CMX_IJ_LL; + m_CMX_IJ_P = b.m_CMX_IJ_P; + m_gamma = b.m_gamma; + + m_CounterIJ = b.m_CounterIJ; + m_debugCalc = b.m_debugCalc; + } + return *this; + } + + + /** + * Test matrix for this object + * + * + * test problems: + * 1 = NaCl problem - 5 species - + * the thermo is read in from an XML file + * + * speci molality charge + * Cl- 6.0954 6.0997E+00 -1 + * H+ 1.0000E-08 2.1628E-09 1 + * Na+ 6.0954E+00 6.0997E+00 1 + * OH- 7.5982E-07 1.3977E-06 -1 + * HMW_params____beta0MX__beta1MX__beta2MX__CphiMX_____alphaMX__thetaij + * 10 + * 1 2 0.1775 0.2945 0.0 0.00080 2.0 0.0 + * 1 3 0.0765 0.2664 0.0 0.00127 2.0 0.0 + * 1 4 0.0 0.0 0.0 0.0 0.0 -0.050 + * 2 3 0.0 0.0 0.0 0.0 0.0 0.036 + * 2 4 0.0 0.0 0.0 0.0 0.0 0.0 + * 3 4 0.0864 0.253 0.0 0.0044 2.0 0.0 + * Triplet_interaction_parameters_psiaa'_or_psicc' + * 2 + * 1 2 3 -0.004 + * 1 3 4 -0.006 + */ + HMWSoln::HMWSoln(int testProb) : + MolalityVPSSTP(), + m_formPitzer(PITZERFORM_BASE), + m_formPitzerTemp(PITZER_TEMP_CONSTANT), + m_formGC(2), + m_Pcurrent(OneAtm), + m_IionicMolality(0.0), + m_maxIionicStrength(30.0), + m_TempPitzerRef(298.15), + m_IionicMolalityStoich(0.0), + m_form_A_Debye(A_DEBYE_WATER), + m_A_Debye(1.172576), // units = sqrt(kg/gmol) + m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m + m_waterSS(0), + m_densWaterSS(1000.), + m_waterProps(0), + m_debugCalc(0) + { + if (testProb != 1) { + printf("unknown test problem\n"); + exit(-1); + } + + constructPhaseFile("HMW_NaCl.xml", ""); + + int i = speciesIndex("Cl-"); + int j = speciesIndex("H+"); + int n = i * m_kk + j; + int ct = m_CounterIJ[n]; + m_Beta0MX_ij[ct] = 0.1775; + m_Beta1MX_ij[ct] = 0.2945; + m_CphiMX_ij[ct] = 0.0008; + m_Alpha1MX_ij[ct]= 2.000; + + + i = speciesIndex("Cl-"); + j = speciesIndex("Na+"); + n = i * m_kk + j; + ct = m_CounterIJ[n]; + m_Beta0MX_ij[ct] = 0.0765; + m_Beta1MX_ij[ct] = 0.2664; + m_CphiMX_ij[ct] = 0.00127; + m_Alpha1MX_ij[ct]= 2.000; + + + i = speciesIndex("Cl-"); + j = speciesIndex("OH-"); + n = i * m_kk + j; + ct = m_CounterIJ[n]; + m_Theta_ij[ct] = -0.05; + + i = speciesIndex("H+"); + j = speciesIndex("Na+"); + n = i * m_kk + j; + ct = m_CounterIJ[n]; + m_Theta_ij[ct] = 0.036; + + i = speciesIndex("Na+"); + j = speciesIndex("OH-"); + n = i * m_kk + j; + ct = m_CounterIJ[n]; + m_Beta0MX_ij[ct] = 0.0864; + m_Beta1MX_ij[ct] = 0.253; + m_CphiMX_ij[ct] = 0.0044; + m_Alpha1MX_ij[ct]= 2.000; + + i = speciesIndex("Cl-"); + j = speciesIndex("H+"); + int k = speciesIndex("Na+"); + double param = -0.004; + n = i * m_kk *m_kk + j * m_kk + k ; + m_Psi_ijk[n] = param; + n = i * m_kk *m_kk + k * m_kk + j ; + m_Psi_ijk[n] = param; + n = j * m_kk *m_kk + i * m_kk + k ; + m_Psi_ijk[n] = param; + n = j * m_kk *m_kk + k * m_kk + i ; + m_Psi_ijk[n] = param; + n = k * m_kk *m_kk + j * m_kk + i ; + m_Psi_ijk[n] = param; + n = k * m_kk *m_kk + i * m_kk + j ; + m_Psi_ijk[n] = param; + + i = speciesIndex("Cl-"); + j = speciesIndex("Na+"); + k = speciesIndex("OH-"); + param = -0.006; + n = i * m_kk *m_kk + j * m_kk + k ; + m_Psi_ijk[n] = param; + n = i * m_kk *m_kk + k * m_kk + j ; + m_Psi_ijk[n] = param; + n = j * m_kk *m_kk + i * m_kk + k ; + m_Psi_ijk[n] = param; + n = j * m_kk *m_kk + k * m_kk + i ; + m_Psi_ijk[n] = param; + n = k * m_kk *m_kk + j * m_kk + i ; + m_Psi_ijk[n] = param; + n = k * m_kk *m_kk + i * m_kk + j ; + m_Psi_ijk[n] = param; + + printCoeffs(); + } + + /** + * ~HMWSoln(): (virtual) + * + * Destructor: does nothing: + */ + HMWSoln::~HMWSoln() { + delete m_waterProps; + delete m_waterSS; + } + + /** + * duplMyselfAsThermoPhase(): + * + * This routine operates at the ThermoPhase level to + * duplicate the current object. It uses the copy constructor + * defined above. + */ + ThermoPhase* HMWSoln::duplMyselfAsThermoPhase() { + HMWSoln* mtp = new HMWSoln(*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 HMWSoln::eosType() const { + int res; + switch (m_formGC) { + case 0: + res = cHMWSoln0; + break; + case 1: + res = cHMWSoln1; + break; + case 2: + res = cHMWSoln2; + 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 HMWSoln::enthalpy_mole() const { + getPartialMolarEnthalpies(DATA_PTR(m_tmpV)); + getMoleFractions(DATA_PTR(m_pp)); + double val = mean_X(DATA_PTR(m_tmpV)); +#ifdef DEBUG_HKM + double val0 = 0.0; + for (int k = 0; k < m_kk; k++) { + val0 += m_tmpV[k] * m_pp[k]; + } + //if (val != val0) { + // printf("ERROR\n"); + //} +#endif + return val; + } + + doublereal HMWSoln::relative_enthalpy() const { + getPartialMolarEnthalpies(DATA_PTR(m_tmpV)); + double hbar = mean_X(DATA_PTR(m_tmpV)); + getEnthalpy_RT(DATA_PTR(m_gamma)); + double RT = GasConstant * temperature(); + for (int k = 0; k < m_kk; k++) { + m_gamma[k] *= RT; + } + double h0bar = mean_X(DATA_PTR(m_gamma)); + return (hbar - h0bar); + } + + + + doublereal HMWSoln::relative_molal_enthalpy() const { + double L = relative_enthalpy(); + getMoleFractions(DATA_PTR(m_tmpV)); + double xanion = 0.0; + int kcation = -1; + double xcation = 0.0; + int kanion = -1; + const double *charge = DATA_PTR(m_speciesCharge); + for (int k = 0; k < m_kk; k++) { + if (charge[k] > 0.0) { + if (m_tmpV[k] > xanion) { + xanion = m_tmpV[k]; + kanion = k; + } + } else if (charge[k] < 0.0) { + if (m_tmpV[k] > xcation) { + xcation = m_tmpV[k]; + kcation = k; + } + } + } + if (kcation < 0 || kanion < 0) { + return L; + } + double xuse = xcation; + int kuse = kcation; + double factor = 1; + if (xanion < xcation) { + xuse = xanion; + kuse = kanion; + if (charge[kcation] != 1.0) { + factor = charge[kcation]; + } + } else { + if (charge[kanion] != 1.0) { + factor = charge[kanion]; + } + } + xuse = xuse / factor; + L = L / xuse; + return L; + } + + /** + * Molar internal energy of the solution. Units: J/kmol. + */ + doublereal HMWSoln::intEnergy_mole() const { + double hh = enthalpy_mole(); + double pres = pressure(); + double molarV = 1.0/molarDensity(); + double uu = hh - pres * molarV; + return uu; + } + + doublereal HMWSoln::entropy_mole() const { + getPartialMolarEntropies(DATA_PTR(m_tmpV)); + return mean_X(DATA_PTR(m_tmpV)); + } + + /// Molar Gibbs function. Units: J/kmol. + doublereal HMWSoln::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 HMWSoln::cp_mole() const { + getPartialMolarCp(DATA_PTR(m_tmpV)); + double val = mean_X(DATA_PTR(m_tmpV)); + return val; + } + + /// Molar heat capacity at constant volume. Units: J/kmol/K. + doublereal HMWSoln::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 HMWSoln::pressure() const { + return m_Pcurrent; + } + + /** + * 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. + */ + void HMWSoln::setPressure(doublereal p) { +#ifdef DEBUG_HKM + //printf("setPressure: %g\n", p); +#endif + double temp = temperature(); + /* + * Call the water SS and set it's internal state + */ + m_waterSS->setTempPressure(temp, p); + + /* + * Store the internal density of the water SS. + * Note, we would have to do this for all other + * species if they had pressure dependent properties. + */ + m_densWaterSS = m_waterSS->density(); + /* + * Store the current pressure + */ + m_Pcurrent = p; + /* + * Calculate all of the other standard volumes + * -> note these are constant for now + */ + /* + * Get the partial molar volumes of all of the + * species. -> note this is a lookup for + * water, here since it was done above. + */ + double *vbar = &m_pp[0]; + getPartialMolarVolumes(vbar); + + /* + * Get mole fractions of all species. + */ + double *x = &m_tmpV[0]; + getMoleFractions(x); + + /* + * Calculate the solution molar volume and the + * solution density. + */ + doublereal vtotal = 0.0; + for (int i = 0; i < m_kk; i++) { + vtotal += vbar[i] * x[i]; + } + doublereal dd = meanMolecularWeight() / vtotal; + + /* + * Now, update the State class with the results. This + * store the denisty. + */ + State::setDensity(dd); + + } + + /** + * 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 HMWSoln::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 HMWSoln::thermalExpansionCoeff() const { + throw CanteraError("HMWSoln::thermalExpansionCoeff", + "unimplemented"); + return 0.0; + } + + /** + * Overwritten setDensity() function is necessary because the + * density is not an indendent variable. + * + * This function will now throw an error condition + * + * Note, in general, setting the phase density is now a nonlinear + * calculation. P and T are the fundamental variables. This + * routine should be revamped to do the nonlinear problem + * + * @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 HMWSoln::setDensity(doublereal rho) { + double dens_old = density(); + + if (rho != dens_old) { + throw CanteraError("HMWSoln::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 an overwritten function from the State.h + * class + */ + void HMWSoln::setMolarDensity(doublereal rho) { + throw CanteraError("HMWSoln::setMolarDensity", + "Density is not an independent variable"); + } + + /** + * Overwritten setTemperature(double) from State.h. This + * function sets the temperature, and makes sure that + * the value propagates to underlying objects. + */ + void HMWSoln::setTemperature(double temp) { + m_waterSS->setTemperature(temp); + State::setTemperature(temp); + } + + // + // ------- 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 HMWSoln::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 HMWSoln::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 HMWSoln::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 HMWSoln::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 HMWSoln::getActivities(doublereal* ac) const { + /* + * Update the molality array, m_molalities() + * This requires an update due to mole fractions + */ + s_update_lnMolalityActCoeff(); + /* + * Now calculate the array of activities. + */ + 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 HMWSoln:: + getMolalityActivityCoefficients(doublereal* acMolality) const { + + A_Debye_TP(-1.0, -1.0); + s_update_lnMolalityActCoeff(); + 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 HMWSoln::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. + */ + s_update_lnMolalityActCoeff(); + /* + * + */ + 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 HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const { + getEnthalpy_RT(hbar); + + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnMolalityActCoeff(); + s_update_dlnMolalityActCoeff_dT(); + double T = temperature(); + double RT = GasConstant * T; + for (int k = 0; k < m_kk; k++) { + hbar[k] *= RT; + } + double RTT = RT * 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 + * + * Combining this with the expression H = G + TS yields: + * + * \f[ + * \bar s_k(T,P) = \hat s^0_k(T) - R log(M0 * molality[k] ac[k]) + * - R T^2 d log(ac[k]) / dT + * \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 HMWSoln:: + 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 stored molalities. + */ + s_update_lnMolalityActCoeff(); + + 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. + */ + s_update_dlnMolalityActCoeff_dT(); + 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 + * species standard state molar volumes. However, extensions + * to this will be implemented in the future. + * + * The general relation is + * + * vbar_i = d(chemPot_i)/dP at const T, n + * + * = V0_i + d(Gex)/dP)_T,M + * + * = V0_i + RT d(lnActCoeffi)dP _T,M + * + * So, if the activity coefficients depended on pressure this + * function would be nontrivial. + */ + void HMWSoln::getPartialMolarVolumes(doublereal* vbar) const { + /* + * Get the standard state values + */ + getStandardVolumes(vbar); + /* + * Update the derivatives wrt the activity coefficients. + */ + s_update_lnMolalityActCoeff(); + s_Pitzer_dlnMolalityActCoeff_dP(); + double T = temperature(); + double RT = GasConstant * T; + for (int k = 0; k < m_kk; k++) { + vbar[k] += RT * m_dlnActCoeffMolaldP[k]; + } + } + + /* + * 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 HMWSoln::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. + */ + s_update_lnMolalityActCoeff(); + s_update_dlnMolalityActCoeff_dT(); + s_update_d2lnMolalityActCoeff_dT2(); + 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 HMWSoln::getStandardChemPotentials(doublereal* mu) const { + getGibbs_ref(mu); + doublereal pref; + doublereal delta_p; + for (int k = 1; k < m_kk; k++) { + pref = m_spthermo->refPressure(k); + delta_p = m_Pcurrent - pref; + mu[k] += delta_p * m_speciesSize[k]; + } + mu[0] = m_waterSS->gibbs_mole(); + } + + /** + * 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 HMWSoln::getGibbs_RT(doublereal* grt) const { + getStandardChemPotentials(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 HMWSoln::getPureGibbs(doublereal* gpure) const { + getStandardChemPotentials(gpure); + } + + /** + * + * 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 HMWSoln:: + getEnthalpy_RT(doublereal* hrt) const { + getEnthalpy_RT_ref(hrt); + doublereal pref; + doublereal delta_p; + double RT = _RT(); + for (int k = 1; k < m_kk; k++) { + pref = m_spthermo->refPressure(k); + delta_p = m_Pcurrent - pref; + hrt[k] += delta_p/ RT * m_speciesSize[k]; + } + hrt[0] = m_waterSS->enthalpy_mole(); + hrt[0] /= RT; + } + + /** + * 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 HMWSoln:: + getEntropy_R(doublereal* sr) const { + getEntropy_R_ref(sr); + sr[0] = m_waterSS->entropy_mole(); + sr[0] /= GasConstant; + } + + /** + * 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 HMWSoln::getCp_R(doublereal* cpr) const { + getCp_R_ref(cpr); + cpr[0] = m_waterSS->cp_mole(); + cpr[0] /= GasConstant; + } + + /** + * Get the molar volumes of each species in their standard + * states at the current + * T and P of the solution. + * units = m^3 / kmol + * + * The water calculation is done separately. + */ + void HMWSoln::getStandardVolumes(doublereal *vol) const { + copy(m_speciesSize.begin(), + m_speciesSize.end(), vol); + double dd = m_waterSS->density(); + vol[0] = molecularWeight(0)/dd; + } + + /* + * ------ Thermodynamic Values for the Species Reference States --- + */ + + // -> This is handled by VPStandardStatesTP + + /* + * -------------- Utilities ------------------------------- + */ + + /** + * @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 HMWSoln::setParameters(int n, doublereal* c) { + } + void HMWSoln::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 HMWSoln::setParametersFromXML(const XML_Node& eosdata) { + } + + /** + * Get the saturation pressure for a given temperature. + * Note the limitations of this function. Stability considerations + * concernting multiphase equilibrium are ignored in this + * calculation. Therefore, the call is made directly to the SS of + * water underneath. The object is put back into its original + * state at the end of the call. + */ + doublereal HMWSoln::satPressure(doublereal t) const { + double p_old = pressure(); + double t_old = temperature(); + double pres = m_waterSS->satPressure(t); + /* + * Set the underlying object back to its original state. + */ + m_waterSS->setState_TP(t_old, p_old); + return pres; + } + + /** + * Report the molar volume of species k + * + * units - \f$ m^3 kmol^-1 \f$ + */ + double HMWSoln::speciesMolarVolume(int k) const { + double vol = m_speciesSize[k]; + if (k == 0) { + double dd = m_waterSS->density(); + vol = molecularWeight(0)/dd; + } + return vol; + } + + + /** + * A_Debye_TP() (virtual) + * + * Returns the A_Debye parameter as a function of temperature + * and pressure. This function also sets the internal value + * of the parameter within the object, if it is changeable. + * + * The default is to assume that it is constant, given + * in the initialization process and storred in the + * member double, m_A_Debye + * + * A_Debye = (1/(8 Pi)) sqrt(2 Na dw /1000) + * (e e/(epsilon R T))^3/2 + * + * where epsilon = e_rel * e_naught + * + * Note, this is si units. Frequently, gaussian units are + * used in Pitzer's papers where D is used, D = epsilon/(4 Pi) + * units = A_Debye has units of sqrt(gmol kg-1). + */ + double HMWSoln::A_Debye_TP(double tempArg, double presArg) const { + double T = temperature(); + double A; + if (tempArg != -1.0) { + T = tempArg; + } + double P = pressure(); + if (presArg != -1.0) { + P = presArg; + } + + switch (m_form_A_Debye) { + case A_DEBYE_CONST: + A = m_A_Debye; + break; + case A_DEBYE_WATER: + A = m_waterProps->ADebye(T, P, 0); + //A = WaterProps::ADebye(T, P, 0); + m_A_Debye = A; + break; + default: + printf("shouldn't be here\n"); + exit(-1); + } + return A; + } + + /** + * dA_DebyedT_TP() (virtual) + * + * Returns the derivative of the A_Debye parameter with + * respect to temperature as a function of temperature + * and pressure. + * + * units = A_Debye has units of sqrt(gmol kg-1). + * Temp has units of Kelvin. + */ + double HMWSoln::dA_DebyedT_TP(double tempArg, double presArg) const { + double T = temperature(); + if (tempArg != -1.0) { + T = tempArg; + } + double P = pressure(); + if (presArg != -1.0) { + P = presArg; + } + double dAdT; + switch (m_form_A_Debye) { + case A_DEBYE_CONST: + dAdT = 0.0; + break; + case A_DEBYE_WATER: + dAdT = m_waterProps->ADebye(T, P, 1); + //dAdT = WaterProps::ADebye(T, P, 1); + break; + default: + printf("shouldn't be here\n"); + exit(-1); + } + return dAdT; + } + + /** + * dA_DebyedP_TP() (virtual) + * + * Returns the derivative of the A_Debye parameter with + * respect to pressure, as a function of temperature + * and pressure. + * + * units = A_Debye has units of sqrt(gmol kg-1). + * Pressure has units of pascals. + */ + double HMWSoln::dA_DebyedP_TP(double tempArg, double presArg) const { + double T = temperature(); + if (tempArg != -1.0) { + T = tempArg; + } + double P = pressure(); + if (presArg != -1.0) { + P = presArg; + } + double dAdP; + switch (m_form_A_Debye) { + case A_DEBYE_CONST: + dAdP = 0.0; + break; + case A_DEBYE_WATER: + dAdP = m_waterProps->ADebye(T, P, 3); + break; + default: + printf("shouldn't be here\n"); + exit(-1); + } + return dAdP; + } + + + /** + * Calculate the DH Parameter used for the Enthalpy calcalations + * + * ADebye_L = 4 R T**2 d(Aphi) / dT + * + * where Aphi = A_Debye/3 + * + * units -> J / (kmolK) * sqrt( kg/gmol) + * + */ + double HMWSoln::ADebye_L(double tempArg, double presArg) const { + double dAdT = dA_DebyedT_TP(); + double dAphidT = dAdT /3.0; + double T = temperature(); + if (tempArg != -1.0) { + T = tempArg; + } + double retn = dAphidT * (4.0 * GasConstant * T * T); + return retn; + } + + /** + * Calculate the DH Parameter used for the Volume calcalations + * + * ADebye_V = - 4 R T d(Aphi) / dP + * + * where Aphi = A_Debye/3 + * + * units -> J / (kmolK) * sqrt( kg/gmol) + * + */ + double HMWSoln::ADebye_V(double tempArg, double presArg) const { + double dAdP = dA_DebyedP_TP(); + double dAphidP = dAdP /3.0; + double T = temperature(); + if (tempArg != -1.0) { + T = tempArg; + } + double retn = - dAphidP * (4.0 * GasConstant * T); + return retn; + } + + /** + * Return Pitzer's definition of A_J. This is basically the + * temperature derivative of A_L, and the second derivative + * of Aphi + * It's the DH parameter used in heat capacity calculations + * + * A_J = 2 A_L/T + 4 * R * T * T * d2(A_phi)/dT2 + * + * Units = sqrt(kg/gmol) (R) + * + * where + * ADebye_L = 4 R T**2 d(Aphi) / dT + * + * where Aphi = A_Debye/3 + * + * units -> J / (kmolK) * sqrt( kg/gmol) + * + */ + double HMWSoln::ADebye_J(double tempArg, double presArg) const { + double T = temperature(); + if (tempArg != -1.0) { + T = tempArg; + } + double A_L = ADebye_L(T, presArg); + double d2 = d2A_DebyedT2_TP(T, presArg); + double d2Aphi = d2 / 3.0; + double retn = 2.0 * A_L / T + 4.0 * GasConstant * T * T *d2Aphi; + return retn; + } + + /** + * 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. + * + * units = A_Debye has units of sqrt(gmol kg-1). + * Temp has units of Kelvin. + */ + double HMWSoln::d2A_DebyedT2_TP(double tempArg, double presArg) const { + double T = temperature(); + if (tempArg != -1.0) { + T = tempArg; + } + double P = pressure(); + if (presArg != -1.0) { + P = presArg; + } + double d2AdT2; + switch (m_form_A_Debye) { + case A_DEBYE_CONST: + d2AdT2 = 0.0; + break; + case A_DEBYE_WATER: + d2AdT2 = m_waterProps->ADebye(T, P, 2); + break; + default: + printf("shouldn't be here\n"); + exit(-1); + } + return d2AdT2; + } + + /* + * ----------- Critical State Properties -------------------------- + */ + + /* + * ---------- Other Property Functions + */ + double HMWSoln::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 HMWSoln::err(string msg) const { + throw CanteraError("HMWSoln", + "Unfinished func called: " + msg ); + return 0.0; + } + + + /** + * initLengths(): + * + * This internal function adjusts the lengths of arrays based on + * the number of species. This is done before these arrays are + * populated with parameter values. + */ + void HMWSoln::initLengths() { + m_kk = nSpecies(); + MolalityVPSSTP::initThermo(); + + /* + * Resize lengths equal to the number of species in + * the phase. + */ + int leng = m_kk; + m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral); + m_speciesSize.resize(leng); + m_Aionic.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); + + + int maxCounterIJlen = 1 + (leng-1) * (leng-2) / 2; + + /* + * Figure out the size of the temperature coefficient + * arrays + */ + int TCoeffLength = 1; + if (m_formPitzerTemp == PITZER_TEMP_LINEAR) { + TCoeffLength = 2; + } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) { + TCoeffLength = 5; + } + + m_Beta0MX_ij.resize(maxCounterIJlen, 0.0); + m_Beta0MX_ij_L.resize(maxCounterIJlen, 0.0); + m_Beta0MX_ij_LL.resize(maxCounterIJlen, 0.0); + m_Beta0MX_ij_P.resize(maxCounterIJlen, 0.0); + m_Beta0MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0); + + m_Beta1MX_ij.resize(maxCounterIJlen, 0.0); + m_Beta1MX_ij_L.resize(maxCounterIJlen, 0.0); + m_Beta1MX_ij_LL.resize(maxCounterIJlen, 0.0); + m_Beta1MX_ij_P.resize(maxCounterIJlen, 0.0); + m_Beta1MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0); + + m_Beta2MX_ij.resize(maxCounterIJlen, 0.0); + m_Beta2MX_ij_L.resize(maxCounterIJlen, 0.0); + m_Beta2MX_ij_LL.resize(maxCounterIJlen, 0.0); + m_Beta2MX_ij_P.resize(maxCounterIJlen, 0.0); + + m_CphiMX_ij.resize(maxCounterIJlen, 0.0); + m_CphiMX_ij_L.resize(maxCounterIJlen, 0.0); + m_CphiMX_ij_LL.resize(maxCounterIJlen, 0.0); + m_CphiMX_ij_P.resize(maxCounterIJlen, 0.0); + m_CphiMX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0); + + m_Alpha1MX_ij.resize(maxCounterIJlen, 0.0); + m_Theta_ij.resize(maxCounterIJlen, 0.0); + m_Theta_ij_L.resize(maxCounterIJlen, 0.0); + m_Theta_ij_LL.resize(maxCounterIJlen, 0.0); + m_Theta_ij_P.resize(maxCounterIJlen, 0.0); + + m_Psi_ijk.resize(m_kk*m_kk*m_kk, 0.0); + m_Psi_ijk_L.resize(m_kk*m_kk*m_kk, 0.0); + m_Psi_ijk_LL.resize(m_kk*m_kk*m_kk, 0.0); + m_Psi_ijk_P.resize(m_kk*m_kk*m_kk, 0.0); + + m_Lambda_ij.resize(leng, leng, 0.0); + m_Lambda_ij_L.resize(leng, leng, 0.0); + m_Lambda_ij_LL.resize(leng, leng, 0.0); + m_Lambda_ij_P.resize(leng, leng, 0.0); + + m_lnActCoeffMolal.resize(leng, 0.0); + m_dlnActCoeffMolaldT.resize(leng, 0.0); + m_d2lnActCoeffMolaldT2.resize(leng, 0.0); + m_dlnActCoeffMolaldP.resize(leng, 0.0); + + m_CounterIJ.resize(m_kk*m_kk, 0); + + m_gfunc_IJ.resize(maxCounterIJlen, 0.0); + m_hfunc_IJ.resize(maxCounterIJlen, 0.0); + m_BMX_IJ.resize(maxCounterIJlen, 0.0); + m_BMX_IJ_L.resize(maxCounterIJlen, 0.0); + m_BMX_IJ_LL.resize(maxCounterIJlen, 0.0); + m_BMX_IJ_P.resize(maxCounterIJlen, 0.0); + m_BprimeMX_IJ.resize(maxCounterIJlen, 0.0); + m_BprimeMX_IJ_L.resize(maxCounterIJlen, 0.0); + m_BprimeMX_IJ_LL.resize(maxCounterIJlen, 0.0); + m_BprimeMX_IJ_P.resize(maxCounterIJlen, 0.0); + m_BphiMX_IJ.resize(maxCounterIJlen, 0.0); + m_BphiMX_IJ_L.resize(maxCounterIJlen, 0.0); + m_BphiMX_IJ_LL.resize(maxCounterIJlen, 0.0); + m_BphiMX_IJ_P.resize(maxCounterIJlen, 0.0); + m_Phi_IJ.resize(maxCounterIJlen, 0.0); + m_Phi_IJ_L.resize(maxCounterIJlen, 0.0); + m_Phi_IJ_LL.resize(maxCounterIJlen, 0.0); + m_Phi_IJ_P.resize(maxCounterIJlen, 0.0); + m_Phiprime_IJ.resize(maxCounterIJlen, 0.0); + m_PhiPhi_IJ.resize(maxCounterIJlen, 0.0); + m_PhiPhi_IJ_L.resize(maxCounterIJlen, 0.0); + m_PhiPhi_IJ_LL.resize(maxCounterIJlen, 0.0); + m_PhiPhi_IJ_P.resize(maxCounterIJlen, 0.0); + m_CMX_IJ.resize(maxCounterIJlen, 0.0); + m_CMX_IJ_L.resize(maxCounterIJlen, 0.0); + m_CMX_IJ_LL.resize(maxCounterIJlen, 0.0); + m_CMX_IJ_P.resize(maxCounterIJlen, 0.0); + + m_gamma.resize(leng, 0.0); + + counterIJ_setup(); + } + + + + /** + * Calcuate the natural log of the molality-based + * activity coefficients. + * + */ + void HMWSoln::s_update_lnMolalityActCoeff() const { + + /* + * Calculate the molalities. Currently, the molalities + * may not be current with respect to the contents of the + * State objects' data. + */ + calcMolalities(); + /* + * Calculate the stoichiometric ionic charge. This isn't used in the + * Pitzer formulation. + */ + m_IionicMolalityStoich = 0.0; + for (int k = 0; k < m_kk; k++) { + double z_k = m_speciesCharge[k]; + double zs_k1 = m_speciesCharge_Stoich[k]; + if (z_k == zs_k1) { + m_IionicMolalityStoich += m_molalities[k] * z_k * z_k; + } else { + double 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; + } + + /* + * Update the temperature dependence of the pitzer coefficients + * and their derivatives + */ + s_updatePitzerCoeffWRTemp(); + + /* + * Now do the main calculation. + */ + s_updatePitzerSublnMolalityActCoeff(); + } + + /* + * Set up a counter variable for keeping track of symmetric binary + * interactactions amongst the solute species. + * + * n = m_kk*i + j + * m_Counter[n] = counter + */ + void HMWSoln::counterIJ_setup(void) const { + int n, nc, i, j; + m_CounterIJ.resize(m_kk * m_kk); + int counter = 0; + for (i = 0; i < m_kk; i++) { + n = i; + nc = m_kk * i; + m_CounterIJ[n] = 0; + m_CounterIJ[nc] = 0; + } + for (i = 1; i < (m_kk - 1); i++) { + n = m_kk * i + i; + m_CounterIJ[n] = 0; + for (j = (i+1); j < m_kk; j++) { + n = m_kk * j + i; + nc = m_kk * i + j; + counter++; + m_CounterIJ[n] = counter; + m_CounterIJ[nc] = counter; + } + } + } + + /** + * Calculates the Pitzer coefficients' dependence on the + * temperature. It will also calculate the temperature + * derivatives of the coefficients, as they are important + * in the calculation of the latent heats and the + * heat capacities of the mixtures. + * + * @param doDerivs If >= 1, then the routine will calculate + * the first derivative. If >= 2, the + * routine will calculate the first and second + * temperature derivative. + * default = 2 + */ + void HMWSoln::s_updatePitzerCoeffWRTemp(int doDerivs) const { + + int i, j, n, counterIJ; + const double *beta0MX_coeff; + const double *beta1MX_coeff; + const double *CphiMX_coeff; + double T = temperature(); + double Tr = m_TempPitzerRef; + double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0; + if (m_formPitzerTemp == PITZER_TEMP_LINEAR) { + tlin = T - Tr; + } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) { + tlin = T - Tr; + tquad = T * T - Tr * Tr; + tln = log(T/ Tr); + tinv = 1.0/T - 1.0/Tr; + } + + for (i = 1; i < (m_kk - 1); i++) { + for (j = (i+1); j < m_kk; j++) { + + + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + + beta0MX_coeff = m_Beta0MX_ij_coeff.ptrColumn(counterIJ); + beta1MX_coeff = m_Beta1MX_ij_coeff.ptrColumn(counterIJ); + CphiMX_coeff = m_CphiMX_ij_coeff.ptrColumn(counterIJ); + + switch (m_formPitzerTemp) { + case PITZER_TEMP_CONSTANT: + break; + case PITZER_TEMP_LINEAR: + m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0] + + beta0MX_coeff[1]*tlin; + m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1]; + m_Beta0MX_ij_LL[counterIJ] = 0.0; + m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0] + + beta1MX_coeff[1]*tlin; + m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1]; + m_Beta1MX_ij_LL[counterIJ] = 0.0; + m_CphiMX_ij [counterIJ] = CphiMX_coeff[0] + + CphiMX_coeff[1]*tlin; + m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1]; + m_CphiMX_ij_LL[counterIJ] = 0.0; + break; + + case PITZER_TEMP_COMPLEX1: + m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0] + + beta0MX_coeff[1]*tlin + + beta0MX_coeff[2]*tquad + + beta0MX_coeff[3]*tinv + + beta0MX_coeff[4]*tln; + + m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0] + + beta1MX_coeff[1]*tlin + + beta1MX_coeff[2]*tquad + + beta1MX_coeff[3]*tinv + + beta1MX_coeff[4]*tln; + + m_CphiMX_ij[counterIJ] = CphiMX_coeff[0] + + CphiMX_coeff[1]*tlin + + CphiMX_coeff[2]*tquad + + CphiMX_coeff[3]*tinv + + CphiMX_coeff[4]*tln; + + m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1] + + beta0MX_coeff[2]*2.0*T + - beta0MX_coeff[3]/(T*T) + + beta0MX_coeff[4]/T; + + m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1] + + beta1MX_coeff[2]*2.0*T + - beta1MX_coeff[3]/(T*T) + + beta1MX_coeff[4]/T; + + + m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1] + + CphiMX_coeff[2]*2.0*T + - CphiMX_coeff[3]/(T*T) + + CphiMX_coeff[4]/T; + + doDerivs = 2; + if (doDerivs > 1) { + m_Beta0MX_ij_LL[counterIJ] = + + beta0MX_coeff[2]*2.0 + + 2.0*beta0MX_coeff[3]/(T*T*T) + - beta0MX_coeff[4]/(T*T); + + m_Beta1MX_ij_LL[counterIJ] = + + beta1MX_coeff[2]*2.0 + + 2.0*beta1MX_coeff[3]/(T*T*T) + - beta1MX_coeff[4]/(T*T); + + m_CphiMX_ij_LL[counterIJ] = + + CphiMX_coeff[2]*2.0 + + 2.0*CphiMX_coeff[3]/(T*T*T) + - CphiMX_coeff[4]/(T*T); + } + +#ifdef DEBUG_HKM + /* + * Turn terms off for debugging + */ + //m_Beta0MX_ij_L[counterIJ] = 0; + //m_Beta0MX_ij_LL[counterIJ] = 0; + //m_Beta1MX_ij_L[counterIJ] = 0; + //m_Beta1MX_ij_LL[counterIJ] = 0; + //m_CphiMX_ij_L[counterIJ] = 0; + //m_CphiMX_ij_LL[counterIJ] = 0; +#endif + break; + } + + + + } + } + + } + /** + * Calculate the Pitzer portion of the activity coefficients. + * + * This is the main routine in the whole module. It calculates the + * molality based activity coefficients for the solutes, and + * the activity of water. + */ + void HMWSoln:: + s_updatePitzerSublnMolalityActCoeff() const { + + /* + * HKM -> Assumption is made that the solvent is + * species 0. + */ + if (m_indexSolvent != 0) { + printf("Wrong index solvent value!\n"); + exit(-1); + } + +#ifdef DEBUG_HKM + int printE = 0; + if (temperature() == 323.15) { + printE = 0; + } +#endif + double wateract; + string sni, snj, snk; + + /* + * This is the molality of the species in solution. + */ + const double *molality = DATA_PTR(m_molalities); + /* + * These are the charges of the species accessed from Constituents.h + */ + const double *charge = DATA_PTR(m_speciesCharge); + + /* + * These are data inputs about the Pitzer correlation. They come + * from the input file for the Pitzer model. + */ + const double *beta0MX = DATA_PTR(m_Beta0MX_ij); + const double *beta1MX = DATA_PTR(m_Beta1MX_ij); + const double *beta2MX = DATA_PTR(m_Beta2MX_ij); + const double *CphiMX = DATA_PTR(m_CphiMX_ij); + const double *thetaij = DATA_PTR(m_Theta_ij); + const double *alphaMX = DATA_PTR(m_Alpha1MX_ij); + + const double *psi_ijk = DATA_PTR(m_Psi_ijk); + //n = k + j * m_kk + i * m_kk * m_kk; + + + double *gamma = DATA_PTR(m_gamma); + /* + * Local variables defined by Coltrin + */ + double etheta[5][5], etheta_prime[5][5], sqrtIs; + /* + * Molality based ionic strength of the solution + */ + double Is = 0.0; + /* + * Molarcharge of the solution: In Pitzer's notation, + * this is his variable called "Z". + */ + double molarcharge = 0.0; + /* + * molalitysum is the sum of the molalities over all solutes, + * even those with zero charge. + */ + double molalitysum = 0.0; + + double *g = DATA_PTR(m_gfunc_IJ); + double *hfunc = DATA_PTR(m_hfunc_IJ); + double *BMX = DATA_PTR(m_BMX_IJ); + double *BprimeMX = DATA_PTR(m_BprimeMX_IJ); + double *BphiMX = DATA_PTR(m_BphiMX_IJ); + double *Phi = DATA_PTR(m_Phi_IJ); + double *Phiprime = DATA_PTR(m_Phiprime_IJ); + double *Phiphi = DATA_PTR(m_PhiPhi_IJ); + double *CMX = DATA_PTR(m_CMX_IJ); + + + double x, g12rooti, gprime12rooti; + double Aphi, F, zsqF; + double sum1, sum2, sum3, sum4, sum5, term1; + double sum_m_phi_minus_1, osmotic_coef, lnwateract; + + int z1, z2; + int n, i, j, k, m, counterIJ, counterIJ2; + +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf("\n Debugging information from hmw_act \n"); + } +#endif + /* + * Make sure the counter variables are setup + */ + counterIJ_setup(); + + /* + * ---------- Calculate common sums over solutes --------------------- + */ + for (n = 1; n < m_kk; n++) { + // ionic strength + Is += charge[n] * charge[n] * molality[n]; + // total molar charge + molarcharge += fabs(charge[n]) * molality[n]; + molalitysum += molality[n]; + } + Is *= 0.5; + if (Is > m_maxIionicStrength) { + Is = m_maxIionicStrength; + } + /* + * Store the ionic molality in the object for reference. + */ + m_IionicMolality = Is; + sqrtIs = sqrt(Is); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 1: \n"); + printf(" ionic strenth = %14.7le \n total molar " + "charge = %14.7le \n", Is, molarcharge); + } +#endif + + /* + * The following call to calc_lambdas() calculates all 16 elements + * of the elambda and elambda1 arrays, given the value of the + * ionic strength (Is) + */ + calc_lambdas(Is); + + /* + * ----- Step 2: Find the coefficients E-theta and ------------------- + * E-thetaprime for all combinations of positive + * unlike charges up to 4 + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 2: \n"); + } +#endif + for (z1 = 1; z1 <=4; z1++) { + for (z2 =1; z2 <=4; z2++) { + calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", + z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); + } +#endif + } + } + +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 3: \n"); + printf(" Species Species g(x) " + " hfunc(x) \n"); + } +#endif + + /* + * + * calculate g(x) and hfunc(x) for each cation-anion pair MX + * In the original literature, hfunc, was called gprime. However, + * it's not the derivative of g(x), so I renamed it. + */ + for (i = 1; i < (m_kk - 1); i++) { + for (j = (i+1); j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * Only loop over oppositely charge species + */ + if (charge[i]*charge[j] < 0) { + /* + * x is a reduced function variable + */ + x = sqrtIs * alphaMX[counterIJ]; + if (x > 1.0E-100) { + g[counterIJ] = 2.0*(1.0-(1.0 + x) * exp(-x)) / (x*x); + hfunc[counterIJ] = -2.0* + (1.0-(1.0 + x + 0.5*x*x) * exp(-x)) / (x*x); + } + else { + g[counterIJ] = 0.0; + hfunc[counterIJ] = 0.0; + } + } + else { + g[counterIJ] = 0.0; + hfunc[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(), + g[counterIJ], hfunc[counterIJ]); + } +#endif + } + } + + /* + * --------- SUBSECTION TO CALCULATE BMX, BprimeMX, BphiMX ---------- + * --------- Agrees with Pitzer, Eq. (49), (51), (55) + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 4: \n"); + printf(" Species Species BMX " + "BprimeMX BphiMX \n"); + } +#endif + x = 12.0 * sqrtIs; + if (x > 1.0E-100) { + g12rooti = 2.0*(1.0-(1.0 + x) * exp(-x)) / (x*x); + gprime12rooti = -2.0*(1.0-(1.0 + x + 0.5*x*x) * exp(-x)) / (x*x); + } else { + g12rooti = 0.0; + gprime12rooti = 0.0; + } + + for (i = 1; i < m_kk - 1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + +#ifdef DEBUG_HKM + if (printE) { + if (counterIJ == 2) { + printf("%s %s\n", speciesName(i).c_str(), + speciesName(j).c_str()); + printf("beta0MX[%d] = %g\n", counterIJ, beta0MX[counterIJ]); + printf("beta1MX[%d] = %g\n", counterIJ, beta1MX[counterIJ]); + } + } +#endif + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0.0) { + BMX[counterIJ] = beta0MX[counterIJ] + + beta1MX[counterIJ] * g[counterIJ] + + beta2MX[counterIJ] * g12rooti; +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf("%d %g: %g %g %g\n", + counterIJ, BMX[counterIJ], beta0MX[counterIJ], + beta1MX[counterIJ], g[counterIJ]); + } +#endif + if (Is > 1.0E-150) { + BprimeMX[counterIJ] = (beta1MX[counterIJ] * hfunc[counterIJ]/Is + + beta2MX[counterIJ] * gprime12rooti/Is); + } else { + BprimeMX[counterIJ] = 0.0; + } + BphiMX[counterIJ] = BMX[counterIJ] + Is*BprimeMX[counterIJ]; + } + else { + BMX[counterIJ] = 0.0; + BprimeMX[counterIJ] = 0.0; + BphiMX[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %11.7f %11.7f %11.7f \n", + sni.c_str(), snj.c_str(), + BMX[counterIJ], BprimeMX[counterIJ], BphiMX[counterIJ] ); + } +#endif + } + } + + /* + * --------- SUBSECTION TO CALCULATE CMX ---------- + * --------- Agrees with Pitzer, Eq. (53). + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 5: \n"); + printf(" Species Species CMX \n"); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0.0) { + CMX[counterIJ] = CphiMX[counterIJ]/ + (2.0* sqrt(fabs(charge[i]*charge[j]))); + } + else { + CMX[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (printE) { + if (counterIJ == 2) { + printf("%s %s\n", speciesName(i).c_str(), + speciesName(j).c_str()); + printf("CphiMX[%d] = %g\n", counterIJ, CphiMX[counterIJ]); + } + } +#endif +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(), + CMX[counterIJ]); + } +#endif + } + } + + /* + * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- + * --------- Agrees with Pitzer, Eq. 72, 73, 74 + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 6: \n"); + printf(" Species Species Phi_ij " + " Phiprime_ij Phi^phi_ij \n"); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] > 0) { + z1 = (int) fabs(charge[i]); + z2 = (int) fabs(charge[j]); + Phi[counterIJ] = thetaij[counterIJ] + etheta[z1][z2]; + Phiprime[counterIJ] = etheta_prime[z1][z2]; + Phiphi[counterIJ] = Phi[counterIJ] + Is * Phiprime[counterIJ]; + } + else { + Phi[counterIJ] = 0.0; + Phiprime[counterIJ] = 0.0; + Phiphi[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", + sni.c_str(), snj.c_str(), + Phi[counterIJ], Phiprime[counterIJ], Phiphi[counterIJ] ); + } +#endif + } + } + + /* + * ------------- SUBSECTION FOR CALCULATION OF F ---------------------- + * ------------ Agrees with Pitzer Eqn. (65) -------------------------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 7: \n"); + } +#endif + // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol)) + // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data + // ( A * ln(10) = 1.17593) + // Aphi = A_Debye_Huckel * 2.30258509 / 3.0; + Aphi = m_A_Debye / 3.0; + F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); +#ifdef DEBUG_HKM + if (printE) { + printf("Aphi = %20.13g\n", Aphi); + } +#endif +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" initial value of F = %10.6f \n", F ); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0) { + F = F + molality[i]*molality[j] * BprimeMX[counterIJ]; + } + /* + * Both species have a non-zero charge, and they + * have the same sign + */ + if (charge[i]*charge[j] > 0) { + F = F + molality[i]*molality[j] * Phiprime[counterIJ]; + } +#ifdef DEBUG_HKM + if (m_debugCalc) printf(" F = %10.6f \n", F ); +#endif + } + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 8: \n"); + } +#endif + + for (i = 1; i < m_kk; i++) { + + /* + * -------- SUBSECTION FOR CALCULATING THE ACTCOEFF FOR CATIONS ----- + * -------- -> equations agree with my notes, Eqn. (118). + * -> Equations agree with Pitzer, eqn.(63) + */ + if (charge[i] > 0 ) { + // species i is the cation (positive) to calc the actcoeff + zsqF = charge[i]*charge[i]*F; + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + for (j = 1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + + if (charge[j] < 0.0) { + // sum over all anions + sum1 = sum1 + molality[j]* + (2.0*BMX[counterIJ]+molarcharge*CMX[counterIJ]); + if (j < m_kk-1) { + /* + * This term is the ternary interaction involving the + * non-duplicate sum over double anions, j, k, with + * respect to the cation, i. + */ + for (k = j+1; k < m_kk; k++) { + // an inner sum over all anions + if (charge[k] < 0.0) { + n = k + j * m_kk + i * m_kk * m_kk; + sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n]; + } + } + } + } + + + if (charge[j] > 0.0) { + // sum over all cations + if (j != i) sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]); + for (k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + // two inner sums over anions + + n = k + j * m_kk + i * m_kk * m_kk; + sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n]; + /* + * Find the counterIJ for the j,k interaction + */ + n = m_kk*j + k; + counterIJ2 = m_CounterIJ[n]; + sum4 = sum4 + (fabs(charge[i])* + molality[j]*molality[k]*CMX[counterIJ2]); + } + } + } + + /* + * Handle neutral j species + */ + if (charge[j] == 0) { + sum5 = sum5 + molality[j]*2.0*m_Lambda_ij(j,i); + } + } + /* + * Add all of the contributions up to yield the log of the + * solute activity coefficients (molality scale) + */ + m_lnActCoeffMolal[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5; + gamma[i] = exp(m_lnActCoeffMolal[i]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n", + sni.c_str(), m_lnActCoeffMolal[i], gamma[i]); + printf(" %12g %12g %12g %12g %12g %12g\n", + zsqF, sum1, sum2, sum3, sum4, sum5); + } +#endif + } + + /* + * -------- SUBSECTION FOR CALCULATING THE ACTCOEFF FOR ANIONS ------ + * -------- -> equations agree with my notes, Eqn. (119). + * -> Equations agree with Pitzer, eqn.(64) + */ + if (charge[i] < 0 ) { + // species i is an anion (negative) + zsqF = charge[i]*charge[i]*F; + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + for (j = 1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + + /* + * For Anions, do the cation interactions. + */ + if (charge[j] > 0) { + sum1 = sum1 + molality[j]* + (2.0*BMX[counterIJ]+molarcharge*CMX[counterIJ]); + if (j < m_kk-1) { + for (k = j+1; k < m_kk; k++) { + // an inner sum over all cations + if (charge[k] > 0) { + n = k + j * m_kk + i * m_kk * m_kk; + sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n]; + } + } + } + } + + /* + * For Anions, do the other anion interactions. + */ + if (charge[j] < 0.0) { + // sum over all anions + if (j != i) { + sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]); + } + for (k = 1; k < m_kk; k++) { + if (charge[k] > 0.0) { + // two inner sums over cations + n = k + j * m_kk + i * m_kk * m_kk; + sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n]; + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*j + k; + counterIJ2 = m_CounterIJ[n]; + sum4 = sum4 + + (fabs(charge[i])* + molality[j]*molality[k]*CMX[counterIJ2]); + } + } + } + + /* + * for Anions, do the neutral species interaction + */ + if (charge[j] == 0.0) { + sum5 = sum5 + molality[j]*2.0*m_Lambda_ij(j,i); + } + } + m_lnActCoeffMolal[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5; + gamma[i] = exp(m_lnActCoeffMolal[i]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n", + sni.c_str(), m_lnActCoeffMolal[i], gamma[i]); + printf(" %12g %12g %12g %12g %12g %12g\n", + zsqF, sum1, sum2, sum3, sum4, sum5); + } +#endif + } + /* + * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF ------- + * ------ -> equations agree with my notes, + * -> Equations agree with Pitzer, + */ + if (charge[i] == 0.0 ) { + sum1 = 0.0; + for (j = 1; j < m_kk; j++) { + sum1 = sum1 + molality[j]*2.0*m_Lambda_ij(i,j); + } + m_lnActCoeffMolal[i] = sum1; + gamma[i] = exp(m_lnActCoeffMolal[i]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n", + sni.c_str(), m_lnActCoeffMolal[i], gamma[i]); + } +#endif + } + + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 9: \n"); + } +#endif + /* + * -------- SUBSECTION FOR CALCULATING THE OSMOTIC COEFF --------- + * -------- -> equations agree with my notes, Eqn. (117). + * -> Equations agree with Pitzer, eqn.(62) + */ + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + double sum6 = 0.0; + /* + * term1 is the DH term in the osmotic coefficient expression + * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer + * implementations. + * Is = Ionic strength on the molality scale (units of (gmol/kg)) + * Aphi = A_Debye / 3 (units of sqrt(kg/gmol)) + */ + term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is)); + + for (j = 1; j < m_kk; j++) { + /* + * Loop Over Cations + */ + if (charge[j] > 0.0) { + for (k = 1; k < m_kk; k++){ + if (charge[k] < 0.0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + + sum1 = sum1 + molality[j]*molality[k]* + (BphiMX[counterIJ] + molarcharge*CMX[counterIJ]); + } + } + + for (k = j+1; k < m_kk; k++) { + if (j == (m_kk-1)) { + // we should never reach this step + printf("logic error 1 in Step 9 of hmw_act"); + exit(1); + } + if (charge[k] > 0.0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + * between 2 cations. + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + sum2 = sum2 + molality[j]*molality[k]*Phiphi[counterIJ]; + for (m = 1; m < m_kk; m++) { + if (charge[m] < 0.0) { + // species m is an anion + n = m + k * m_kk + j * m_kk * m_kk; + sum2 = sum2 + + molality[j]*molality[k]*molality[m]*psi_ijk[n]; + } + } + } + } + } + + /* + * Loop Over Anions + */ + if (charge[j] < 0) { + for (k = j+1; k < m_kk; k++) { + if (j == m_kk-1) { + // we should never reach this step + printf("logic error 2 in Step 9 of hmw_act"); + exit(1); + } + if (charge[k] < 0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + * between two anions + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + + sum3 = sum3 + molality[j]*molality[k]*Phiphi[counterIJ]; + for (m = 1; m < m_kk; m++) { + if (charge[m] > 0.0) { + n = m + k * m_kk + j * m_kk * m_kk; + sum3 = sum3 + + molality[j]*molality[k]*molality[m]*psi_ijk[n]; + } + } + } + } + } + + /* + * Loop Over Neutral Species + */ + if (charge[j] == 0) { + for (k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + sum4 = sum4 + molality[j]*molality[k]*m_Lambda_ij(j,k); + } + if (charge[k] > 0.0) { + sum5 = sum5 + molality[j]*molality[k]*m_Lambda_ij(j,k); + } + if (charge[k] == 0.0) { + if (k > j) { + sum6 = sum6 + molality[j]*molality[k]*m_Lambda_ij(j,k); + } else if (k == j) { + sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_ij(j,k); + } + } + } + } + } + sum_m_phi_minus_1 = 2.0 * + (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6); + /* + * Calculate the osmotic coefficient from + * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i) + */ + if (molalitysum > 1.0E-150) { + osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysum); + } else { + osmotic_coef = 1.0; + } +#ifdef DEBUG_HKM + if (printE) { + + printf("OsmCoef - 1 = %20.13g\n", osmotic_coef - 1.0); + } +#endif +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f " + "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n", + term1, sum1, sum2, sum3, sum4, sum5); + printf(" sum_m_phi_minus_1=%10.6f osmotic_coef=%10.6f\n", + sum_m_phi_minus_1, osmotic_coef); + } + + if (m_debugCalc) { + printf(" Step 10: \n"); + } +#endif + lnwateract = -(m_weightSolvent/1000.0) * molalitysum * osmotic_coef; + wateract = exp(lnwateract); + + /* + * In Cantera, we define the activity coefficient of the solvent as + * + * act_0 = actcoeff_0 * Xmol_0 + * + * We have just computed act_0. However, this routine returns + * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). + */ + double xmolSolvent = moleFraction(m_indexSolvent); + m_lnActCoeffMolal[0] = lnwateract - log(xmolSolvent); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Weight of Solvent = %16.7g\n", m_weightSolvent); + printf(" molalitySum = %16.7g\n", molalitysum); + printf(" ln_a_water=%10.6f a_water=%10.6f\n\n", + lnwateract, wateract); + } +#endif + } + + /** + * s_update_dlnMolalityActCoeff_dT() (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 HMWSoln::s_update_dlnMolalityActCoeff_dT() const { + + for (int k = 0; k < m_kk; k++) { + m_dlnActCoeffMolaldT[k] = 0.0; + } + s_Pitzer_dlnMolalityActCoeff_dT(); + } + + /*************************************************************************************/ + + /** + * Calculate the Pitzer portion of the temperature + * derivative of the log activity coefficients. + * This is an internal routine. + * + * It may be assumed that the + * Pitzer activity coefficient routine is called immediately + * preceding the calling of this routine. Therefore, some + * quantities do not need to be recalculated in this routine. + * + */ + void HMWSoln::s_Pitzer_dlnMolalityActCoeff_dT() const { + + /* + * HKM -> Assumption is made that the solvent is + * species 0. + */ +#ifdef DEBUG_HKM + m_debugCalc = 0; +#endif + if (m_indexSolvent != 0) { + printf("Wrong index solvent value!\n"); + exit(-1); + } + + double d_wateract_dT; + string sni, snj, snk; + + const double *molality = DATA_PTR(m_molalities); + const double *charge = DATA_PTR(m_speciesCharge); + const double *beta0MX_L = DATA_PTR(m_Beta0MX_ij_L); + const double *beta1MX_L = DATA_PTR(m_Beta1MX_ij_L); + const double *beta2MX_L = DATA_PTR(m_Beta2MX_ij_L); + const double *CphiMX_L = DATA_PTR(m_CphiMX_ij_L); + const double *thetaij_L = DATA_PTR(m_Theta_ij_L); + const double *alphaMX = DATA_PTR(m_Alpha1MX_ij); + const double *psi_ijk_L = DATA_PTR(m_Psi_ijk_L); + double *gamma = DATA_PTR(m_gamma); + /* + * Local variables defined by Coltrin + */ + double etheta[5][5], etheta_prime[5][5], sqrtIs; + /* + * Molality based ionic strength of the solution + */ + double Is = 0.0; + /* + * Molarcharge of the solution: In Pitzer's notation, + * this is his variable called "Z". + */ + double molarcharge = 0.0; + /* + * molalitysum is the sum of the molalities over all solutes, + * even those with zero charge. + */ + double molalitysum = 0.0; + + double *g = DATA_PTR(m_gfunc_IJ); + double *hfunc = DATA_PTR(m_hfunc_IJ); + double *BMX_L = DATA_PTR(m_BMX_IJ_L); + double *BprimeMX_L= DATA_PTR(m_BprimeMX_IJ_L); + double *BphiMX_L = DATA_PTR(m_BphiMX_IJ_L); + double *Phi_L = DATA_PTR(m_Phi_IJ_L); + double *Phiprime = DATA_PTR(m_Phiprime_IJ); + double *Phiphi_L = DATA_PTR(m_PhiPhi_IJ_L); + double *CMX_L = DATA_PTR(m_CMX_IJ_L); + + double x, g12rooti, gprime12rooti; + double Aphi, dFdT, zsqdFdT; + double sum1, sum2, sum3, sum4, sum5, term1; + double sum_m_phi_minus_1, d_osmotic_coef_dT, d_lnwateract_dT; + + int z1, z2; + int n, i, j, k, m, counterIJ, counterIJ2; + +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf("\n Debugging information from " + "s_Pitzer_dlnMolalityActCoeff_dT()\n"); + } +#endif + /* + * Make sure the counter variables are setup + */ + counterIJ_setup(); + + /* + * ---------- Calculate common sums over solutes --------------------- + */ + for (n = 1; n < m_kk; n++) { + // ionic strength + Is += charge[n] * charge[n] * molality[n]; + // total molar charge + molarcharge += fabs(charge[n]) * molality[n]; + molalitysum += molality[n]; + } + Is *= 0.5; + if (Is > m_maxIionicStrength) { + Is = m_maxIionicStrength; + } + /* + * Store the ionic molality in the object for reference. + */ + m_IionicMolality = Is; + sqrtIs = sqrt(Is); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 1: \n"); + printf(" ionic strenth = %14.7le \n total molar " + "charge = %14.7le \n", Is, molarcharge); + } +#endif + + /* + * The following call to calc_lambdas() calculates all 16 elements + * of the elambda and elambda1 arrays, given the value of the + * ionic strength (Is) + */ + calc_lambdas(Is); + + /* + * ----- Step 2: Find the coefficients E-theta and ------------------- + * E-thetaprime for all combinations of positive + * unlike charges up to 4 + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 2: \n"); + } +#endif + for (z1 = 1; z1 <=4; z1++) { + for (z2 =1; z2 <=4; z2++) { + calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", + z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); + } +#endif + } + } + +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 3: \n"); + printf(" Species Species g(x) " + " hfunc(x) \n"); + } +#endif + + /* + * + * calculate g(x) and hfunc(x) for each cation-anion pair MX + * In the original literature, hfunc, was called gprime. However, + * it's not the derivative of g(x), so I renamed it. + */ + for (i = 1; i < (m_kk - 1); i++) { + for (j = (i+1); j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * Only loop over oppositely charge species + */ + if (charge[i]*charge[j] < 0) { + /* + * x is a reduced function variable + */ + x = sqrtIs * alphaMX[counterIJ]; + if (x > 1.0E-100) { + g[counterIJ] = 2.0*(1.0-(1.0 + x) * exp(-x)) / (x*x); + hfunc[counterIJ] = -2.0* + (1.0-(1.0 + x + 0.5*x*x) * exp(-x)) / (x*x); + } + else { + g[counterIJ] = 0.0; + hfunc[counterIJ] = 0.0; + } + } + else { + g[counterIJ] = 0.0; + hfunc[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(), + g[counterIJ], hfunc[counterIJ]); + } +#endif + } + } + + /* + * ------- SUBSECTION TO CALCULATE BMX_L, BprimeMX_L, BphiMX_L ---------- + * ------- These are now temperature derivatives of the + * previously calculated quantities. + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 4: \n"); + printf(" Species Species BMX " + "BprimeMX BphiMX \n"); + } +#endif + x = 12.0 * sqrtIs; + if (x > 1.0E-100) { + g12rooti = 2.0*(1.0-(1.0 + x) * exp(-x)) / (x*x); + gprime12rooti = -2.0*(1.0-(1.0 + x + 0.5*x*x) * exp(-x)) / (x*x); + } else { + g12rooti = 0.0; + gprime12rooti = 0.0; + } + + for (i = 1; i < m_kk - 1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0.0) { + BMX_L[counterIJ] = beta0MX_L[counterIJ] + + beta1MX_L[counterIJ] * g[counterIJ] + + beta2MX_L[counterIJ] * g12rooti; +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf("%d %g: %g %g %g\n", + counterIJ, BMX_L[counterIJ], beta0MX_L[counterIJ], + beta1MX_L[counterIJ], g[counterIJ]); + } +#endif + if (Is > 1.0E-150) { + BprimeMX_L[counterIJ] = (beta1MX_L[counterIJ] * hfunc[counterIJ]/Is + + beta2MX_L[counterIJ] * gprime12rooti/Is); + } else { + BprimeMX_L[counterIJ] = 0.0; + } + BphiMX_L[counterIJ] = BMX_L[counterIJ] + Is*BprimeMX_L[counterIJ]; + } + else { + BMX_L[counterIJ] = 0.0; + BprimeMX_L[counterIJ] = 0.0; + BphiMX_L[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %11.7f %11.7f %11.7f \n", + sni.c_str(), snj.c_str(), + BMX_L[counterIJ], BprimeMX_L[counterIJ], BphiMX_L[counterIJ]); + } +#endif + } + } + + /* + * --------- SUBSECTION TO CALCULATE CMX_L ---------- + * --------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 5: \n"); + printf(" Species Species CMX \n"); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0.0) { + CMX_L[counterIJ] = CphiMX_L[counterIJ]/ + (2.0* sqrt(fabs(charge[i]*charge[j]))); + } + else { + CMX_L[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(), + CMX_L[counterIJ]); + } +#endif + } + } + + /* + * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- + * -------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 6: \n"); + printf(" Species Species Phi_ij " + " Phiprime_ij Phi^phi_ij \n"); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] > 0) { + z1 = (int) fabs(charge[i]); + z2 = (int) fabs(charge[j]); + //Phi[counterIJ] = thetaij_L[counterIJ] + etheta[z1][z2]; + Phi_L[counterIJ] = thetaij_L[counterIJ]; + //Phiprime[counterIJ] = etheta_prime[z1][z2]; + Phiprime[counterIJ] = 0.0; + Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ]; + } + else { + Phi_L[counterIJ] = 0.0; + Phiprime[counterIJ] = 0.0; + Phiphi_L[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", + sni.c_str(), snj.c_str(), + Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ] ); + } +#endif + } + } + + /* + * ----------- SUBSECTION FOR CALCULATION OF dFdT --------------------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 7: \n"); + } +#endif + // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol)) + // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data + // ( A * ln(10) = 1.17593) + // Aphi = A_Debye_Huckel * 2.30258509 / 3.0; + Aphi = m_A_Debye / 3.0; + + double dA_DebyedT = dA_DebyedT_TP(); + double dAphidT = dA_DebyedT /3.0; +#ifdef DEBUG_HKM + //dAphidT = 0.0; +#endif + //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + // + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); + //dAphidT = Al / (4.0 * GasConstant * T * T); + dFdT = -dAphidT * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" initial value of dFdT = %10.6f \n", dFdT ); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0) { + dFdT = dFdT + molality[i]*molality[j] * BprimeMX_L[counterIJ]; + } + /* + * Both species have a non-zero charge, and they + * have the same sign, e.g., both positive or both negative. + */ + if (charge[i]*charge[j] > 0) { + dFdT = dFdT + molality[i]*molality[j] * Phiprime[counterIJ]; + } +#ifdef DEBUG_HKM + if (m_debugCalc) printf(" dFdT = %10.6f \n", dFdT); +#endif + } + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 8: \n"); + } +#endif + + for (i = 1; i < m_kk; i++) { + + /* + * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS ----- + * -- + */ + if (charge[i] > 0 ) { + // species i is the cation (positive) to calc the actcoeff + zsqdFdT = charge[i]*charge[i]*dFdT; + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + for (j = 1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + + if (charge[j] < 0.0) { + // sum over all anions + sum1 = sum1 + molality[j]* + (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]); + if (j < m_kk-1) { + /* + * This term is the ternary interaction involving the + * non-duplicate sum over double anions, j, k, with + * respect to the cation, i. + */ + for (k = j+1; k < m_kk; k++) { + // an inner sum over all anions + if (charge[k] < 0.0) { + n = k + j * m_kk + i * m_kk * m_kk; + sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n]; + } + } + } + } + + + if (charge[j] > 0.0) { + // sum over all cations + if (j != i) { + sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]); + } + for (k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + // two inner sums over anions + + n = k + j * m_kk + i * m_kk * m_kk; + sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n]; + /* + * Find the counterIJ for the j,k interaction + */ + n = m_kk*j + k; + counterIJ2 = m_CounterIJ[n]; + sum4 = sum4 + (fabs(charge[i])* + molality[j]*molality[k]*CMX_L[counterIJ2]); + } + } + } + + /* + * Handle neutral j species + */ + if (charge[j] == 0) { + sum5 = sum5 + molality[j]*2.0*m_Lambda_ij_L(j,i); + } + } + /* + * Add all of the contributions up to yield the log of the + * solute activity coefficients (molality scale) + */ + m_dlnActCoeffMolaldT[i] = + zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5; + gamma[i] = exp(m_dlnActCoeffMolaldT[i]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n", + sni.c_str(), m_dlnActCoeffMolaldT[i], gamma[i]); + printf(" %12g %12g %12g %12g %12g %12g\n", + zsqdFdT, sum1, sum2, sum3, sum4, sum5); + } +#endif + } + + /* + * ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR ANIONS ------ + * + */ + if (charge[i] < 0 ) { + // species i is an anion (negative) + zsqdFdT = charge[i]*charge[i]*dFdT; + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + for (j = 1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + + /* + * For Anions, do the cation interactions. + */ + if (charge[j] > 0) { + sum1 = sum1 + molality[j]* + (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]); + if (j < m_kk-1) { + for (k = j+1; k < m_kk; k++) { + // an inner sum over all cations + if (charge[k] > 0) { + n = k + j * m_kk + i * m_kk * m_kk; + sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n]; + } + } + } + } + + /* + * For Anions, do the other anion interactions. + */ + if (charge[j] < 0.0) { + // sum over all anions + if (j != i) { + sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]); + } + for (k = 1; k < m_kk; k++) { + if (charge[k] > 0.0) { + // two inner sums over cations + n = k + j * m_kk + i * m_kk * m_kk; + sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n]; + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*j + k; + counterIJ2 = m_CounterIJ[n]; + sum4 = sum4 + + (fabs(charge[i])* + molality[j]*molality[k]*CMX_L[counterIJ2]); + } + } + } + + /* + * for Anions, do the neutral species interaction + */ + if (charge[j] == 0.0) { + sum5 = sum5 + molality[j]*2.0*m_Lambda_ij_L(j,i); + } + } + m_dlnActCoeffMolaldT[i] = + zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5; + gamma[i] = exp(m_dlnActCoeffMolaldT[i]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n", + sni.c_str(), m_dlnActCoeffMolaldT[i], gamma[i]); + printf(" %12g %12g %12g %12g %12g %12g\n", + zsqdFdT, sum1, sum2, sum3, sum4, sum5); + } +#endif + } + /* + * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF ------- + * ------ -> equations agree with my notes, + * -> Equations agree with Pitzer, + */ + if (charge[i] == 0.0 ) { + sum1 = 0.0; + for (j = 1; j < m_kk; j++) { + sum1 = sum1 + molality[j]*2.0*m_Lambda_ij_L(i,j); + } + m_dlnActCoeffMolaldT[i] = sum1; + gamma[i] = exp(m_dlnActCoeffMolaldT[i]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n", + sni.c_str(), m_dlnActCoeffMolaldT[i], gamma[i]); + } +#endif + } + + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 9: \n"); + } +#endif + /* + * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT --------- + * + */ + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + double sum6 = 0.0; + /* + * term1 is the temperature derivative of the + * DH term in the osmotic coefficient expression + * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer + * implementations. + * Is = Ionic strength on the molality scale (units of (gmol/kg)) + * Aphi = A_Debye / 3 (units of sqrt(kg/gmol)) + */ + term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is)); + + for (j = 1; j < m_kk; j++) { + /* + * Loop Over Cations + */ + if (charge[j] > 0.0) { + for (k = 1; k < m_kk; k++){ + if (charge[k] < 0.0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + + sum1 = sum1 + molality[j]*molality[k]* + (BphiMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]); + } + } + + for (k = j+1; k < m_kk; k++) { + if (j == (m_kk-1)) { + // we should never reach this step + printf("logic error 1 in Step 9 of hmw_act"); + exit(1); + } + if (charge[k] > 0.0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + * between 2 cations. + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + sum2 = sum2 + molality[j]*molality[k]*Phiphi_L[counterIJ]; + for (m = 1; m < m_kk; m++) { + if (charge[m] < 0.0) { + // species m is an anion + n = m + k * m_kk + j * m_kk * m_kk; + sum2 = sum2 + + molality[j]*molality[k]*molality[m]*psi_ijk_L[n]; + } + } + } + } + } + + /* + * Loop Over Anions + */ + if (charge[j] < 0) { + for (k = j+1; k < m_kk; k++) { + if (j == m_kk-1) { + // we should never reach this step + printf("logic error 2 in Step 9 of hmw_act"); + exit(1); + } + if (charge[k] < 0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + * between two anions + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + + sum3 = sum3 + molality[j]*molality[k]*Phiphi_L[counterIJ]; + for (m = 1; m < m_kk; m++) { + if (charge[m] > 0.0) { + n = m + k * m_kk + j * m_kk * m_kk; + sum3 = sum3 + + molality[j]*molality[k]*molality[m]*psi_ijk_L[n]; + } + } + } + } + } + + /* + * Loop Over Neutral Species + */ + if (charge[j] == 0) { + for (k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + sum4 = sum4 + molality[j]*molality[k]*m_Lambda_ij_L(j,k); + } + if (charge[k] > 0.0) { + sum5 = sum5 + molality[j]*molality[k]*m_Lambda_ij_L(j,k); + } + if (charge[k] == 0.0) { + if (k > j) { + sum6 = sum6 + molality[j]*molality[k]*m_Lambda_ij_L(j,k); + } else if (k == j) { + sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_ij_L(j,k); + } + } + } + } + } + sum_m_phi_minus_1 = 2.0 * + (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6); + /* + * Calculate the osmotic coefficient from + * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i) + */ + if (molalitysum > 1.0E-150) { + d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum); + } else { + d_osmotic_coef_dT = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f " + "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n", + term1, sum1, sum2, sum3, sum4, sum5); + printf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dT =%10.6f\n", + sum_m_phi_minus_1, d_osmotic_coef_dT); + } + + if (m_debugCalc) { + printf(" Step 10: \n"); + } +#endif + d_lnwateract_dT = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT; + d_wateract_dT = exp(d_lnwateract_dT); + + /* + * In Cantera, we define the activity coefficient of the solvent as + * + * act_0 = actcoeff_0 * Xmol_0 + * + * We have just computed act_0. However, this routine returns + * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). + */ + //double xmolSolvent = moleFraction(m_indexSolvent); + m_dlnActCoeffMolaldT[0] = d_lnwateract_dT; +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" d_ln_a_water_dT = %10.6f d_a_water_dT=%10.6f\n\n", + d_lnwateract_dT, d_wateract_dT); + } +#endif + } + + /*************************************************************************************/ + + + /** + * s_update_d2lnMolalityActCoeff_dT2() (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. + * This is an internal routine + * + * We assume that the activity coefficients and first temperature + * derivatives of the activity coefficients are current. + * + * It may be assumed that the + * Pitzer activity coefficient and first deriv routine are called immediately + * preceding the calling of this routine. Therefore, some + * quantities do not need to be recalculated in this routine. + * + * solvent activity coefficient is on the molality + * scale. It's derivatives are too. + */ + void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const { + + /* + * HKM -> Assumption is made that the solvent is + * species 0. + */ +#ifdef DEBUG_HKM + m_debugCalc = 0; +#endif + if (m_indexSolvent != 0) { + printf("Wrong index solvent value!\n"); + exit(-1); + } + + double d2_wateract_dT2; + string sni, snj, snk; + + const double *molality = DATA_PTR(m_molalities); + const double *charge = DATA_PTR(m_speciesCharge); + const double *beta0MX_LL= DATA_PTR(m_Beta0MX_ij_LL); + const double *beta1MX_LL= DATA_PTR(m_Beta1MX_ij_LL); + const double *beta2MX_LL= DATA_PTR(m_Beta2MX_ij_LL); + const double *CphiMX_LL = DATA_PTR(m_CphiMX_ij_LL); + const double *thetaij_LL= DATA_PTR(m_Theta_ij_LL); + const double *alphaMX = DATA_PTR(m_Alpha1MX_ij); + const double *psi_ijk_LL= DATA_PTR(m_Psi_ijk_LL); + + /* + * Local variables defined by Coltrin + */ + double etheta[5][5], etheta_prime[5][5], sqrtIs; + /* + * Molality based ionic strength of the solution + */ + double Is = 0.0; + /* + * Molarcharge of the solution: In Pitzer's notation, + * this is his variable called "Z". + */ + double molarcharge = 0.0; + /* + * molalitysum is the sum of the molalities over all solutes, + * even those with zero charge. + */ + double molalitysum = 0.0; + + double *g = DATA_PTR(m_gfunc_IJ); + double *hfunc = DATA_PTR(m_hfunc_IJ); + double *BMX_LL = DATA_PTR(m_BMX_IJ_LL); + double *BprimeMX_LL=DATA_PTR(m_BprimeMX_IJ_LL); + double *BphiMX_LL= DATA_PTR(m_BphiMX_IJ_LL); + double *Phi_LL = DATA_PTR(m_Phi_IJ_LL); + double *Phiprime = DATA_PTR(m_Phiprime_IJ); + double *Phiphi_LL= DATA_PTR(m_PhiPhi_IJ_LL); + double *CMX_LL = DATA_PTR(m_CMX_IJ_LL); + + + double x, g12rooti, gprime12rooti; + double d2FdT2, zsqd2FdT2; + double sum1, sum2, sum3, sum4, sum5, term1; + double sum_m_phi_minus_1, d2_osmotic_coef_dT2, d2_lnwateract_dT2; + + int z1, z2; + int n, i, j, k, m, counterIJ, counterIJ2; + +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf("\n Debugging information from " + "s_Pitzer_d2lnMolalityActCoeff_dT2()\n"); + } +#endif + /* + * Make sure the counter variables are setup + */ + counterIJ_setup(); + + + /* + * ---------- Calculate common sums over solutes --------------------- + */ + for (n = 1; n < m_kk; n++) { + // ionic strength + Is += charge[n] * charge[n] * molality[n]; + // total molar charge + molarcharge += fabs(charge[n]) * molality[n]; + molalitysum += molality[n]; + } + Is *= 0.5; + if (Is > m_maxIionicStrength) { + Is = m_maxIionicStrength; + } + /* + * Store the ionic molality in the object for reference. + */ + m_IionicMolality = Is; + sqrtIs = sqrt(Is); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 1: \n"); + printf(" ionic strenth = %14.7le \n total molar " + "charge = %14.7le \n", Is, molarcharge); + } +#endif + + /* + * The following call to calc_lambdas() calculates all 16 elements + * of the elambda and elambda1 arrays, given the value of the + * ionic strength (Is) + */ + calc_lambdas(Is); + + /* + * ----- Step 2: Find the coefficients E-theta and ------------------- + * E-thetaprime for all combinations of positive + * unlike charges up to 4 + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 2: \n"); + } +#endif + for (z1 = 1; z1 <=4; z1++) { + for (z2 =1; z2 <=4; z2++) { + calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", + z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); + } +#endif + } + } + +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 3: \n"); + printf(" Species Species g(x) " + " hfunc(x) \n"); + } +#endif + + /* + * + * calculate g(x) and hfunc(x) for each cation-anion pair MX + * In the original literature, hfunc, was called gprime. However, + * it's not the derivative of g(x), so I renamed it. + */ + for (i = 1; i < (m_kk - 1); i++) { + for (j = (i+1); j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * Only loop over oppositely charge species + */ + if (charge[i]*charge[j] < 0) { + /* + * x is a reduced function variable + */ + x = sqrtIs * alphaMX[counterIJ]; + if (x > 1.0E-100) { + g[counterIJ] = 2.0*(1.0-(1.0 + x) * exp(-x)) / (x*x); + hfunc[counterIJ] = -2.0* + (1.0-(1.0 + x + 0.5*x*x) * exp(-x)) / (x*x); + } + else { + g[counterIJ] = 0.0; + hfunc[counterIJ] = 0.0; + } + } + else { + g[counterIJ] = 0.0; + hfunc[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(), + g[counterIJ], hfunc[counterIJ]); + } +#endif + } + } + /* + * ------- SUBSECTION TO CALCULATE BMX_L, BprimeMX_LL, BphiMX_L ---------- + * ------- These are now temperature derivatives of the + * previously calculated quantities. + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 4: \n"); + printf(" Species Species BMX " + "BprimeMX BphiMX \n"); + } +#endif + x = 12.0 * sqrtIs; + if (x > 1.0E-100) { + g12rooti = 2.0*(1.0-(1.0 + x) * exp(-x)) / (x*x); + gprime12rooti = -2.0*(1.0-(1.0 + x + 0.5*x*x) * exp(-x)) / (x*x); + } else { + g12rooti = 0.0; + gprime12rooti = 0.0; + } + + for (i = 1; i < m_kk - 1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0.0) { + BMX_LL[counterIJ] = beta0MX_LL[counterIJ] + + beta1MX_LL[counterIJ] * g[counterIJ] + + beta2MX_LL[counterIJ] * g12rooti; +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf("%d %g: %g %g %g\n", + counterIJ, BMX_LL[counterIJ], beta0MX_LL[counterIJ], + beta1MX_LL[counterIJ], g[counterIJ]); + } +#endif + if (Is > 1.0E-150) { + BprimeMX_LL[counterIJ] = (beta1MX_LL[counterIJ] * hfunc[counterIJ]/Is + + beta2MX_LL[counterIJ] * gprime12rooti/Is); + } else { + BprimeMX_LL[counterIJ] = 0.0; + } + BphiMX_LL[counterIJ] = BMX_LL[counterIJ] + Is*BprimeMX_LL[counterIJ]; + } + else { + BMX_LL[counterIJ] = 0.0; + BprimeMX_LL[counterIJ] = 0.0; + BphiMX_LL[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %11.7f %11.7f %11.7f \n", + sni.c_str(), snj.c_str(), + BMX_LL[counterIJ], BprimeMX_LL[counterIJ], BphiMX_LL[counterIJ]); + } +#endif + } + } + + /* + * --------- SUBSECTION TO CALCULATE CMX_LL ---------- + * --------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 5: \n"); + printf(" Species Species CMX \n"); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0.0) { + CMX_LL[counterIJ] = CphiMX_LL[counterIJ]/ + (2.0* sqrt(fabs(charge[i]*charge[j]))); + } else { + CMX_LL[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(), + CMX_LL[counterIJ]); + } +#endif + } + } + + /* + * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- + * -------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 6: \n"); + printf(" Species Species Phi_ij " + " Phiprime_ij Phi^phi_ij \n"); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] > 0) { + z1 = (int) fabs(charge[i]); + z2 = (int) fabs(charge[j]); + //Phi[counterIJ] = thetaij[counterIJ] + etheta[z1][z2]; + //Phi_L[counterIJ] = thetaij_L[counterIJ]; + Phi_LL[counterIJ] = thetaij_LL[counterIJ]; + //Phiprime[counterIJ] = etheta_prime[z1][z2]; + Phiprime[counterIJ] = 0.0; + //Phiphi[counterIJ] = Phi[counterIJ] + Is * Phiprime[counterIJ]; + //Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ]; + Phiphi_LL[counterIJ] = Phi_LL[counterIJ]; + } + else { + Phi_LL[counterIJ] = 0.0; + Phiprime[counterIJ] = 0.0; + Phiphi_LL[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + //printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", + // sni.c_str(), snj.c_str(), + // Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ] ); + printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", + sni.c_str(), snj.c_str(), + Phi_LL[counterIJ], Phiprime[counterIJ], Phiphi_LL[counterIJ] ); + } +#endif + } + } + + /* + * ----------- SUBSECTION FOR CALCULATION OF d2FdT2 --------------------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 7: \n"); + } +#endif + // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol)) + // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data + // ( A * ln(10) = 1.17593) + // Aphi = A_Debye_Huckel * 2.30258509 / 3.0; + // Aphi = m_A_Debye / 3.0; + + //double dA_DebyedT = dA_DebyedT_TP(); + //double dAphidT = dA_DebyedT /3.0; + double d2AphidT2 = d2A_DebyedT2_TP() / 3.0; +#ifdef DEBUG_HKM + //d2AphidT2 = 0.0; +#endif + //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + // + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); + //dAphidT = Al / (4.0 * GasConstant * T * T); + //dFdT = -dAphidT * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + // + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); + d2FdT2 = -d2AphidT2 * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" initial value of d2FdT2 = %10.6f \n", d2FdT2 ); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0) { + d2FdT2 = d2FdT2 + molality[i]*molality[j] * BprimeMX_LL[counterIJ]; + } + /* + * Both species have a non-zero charge, and they + * have the same sign, e.g., both positive or both negative. + */ + if (charge[i]*charge[j] > 0) { + d2FdT2 = d2FdT2 + molality[i]*molality[j] * Phiprime[counterIJ]; + } +#ifdef DEBUG_HKM + if (m_debugCalc) printf(" d2FdT2 = %10.6f \n", d2FdT2); +#endif + } + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 8: \n"); + } +#endif + + for (i = 1; i < m_kk; i++) { + + /* + * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS ----- + * -- + */ + if (charge[i] > 0 ) { + // species i is the cation (positive) to calc the actcoeff + zsqd2FdT2 = charge[i]*charge[i]*d2FdT2; + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + for (j = 1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + + if (charge[j] < 0.0) { + // sum over all anions + sum1 = sum1 + molality[j]* + (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]); + if (j < m_kk-1) { + /* + * This term is the ternary interaction involving the + * non-duplicate sum over double anions, j, k, with + * respect to the cation, i. + */ + for (k = j+1; k < m_kk; k++) { + // an inner sum over all anions + if (charge[k] < 0.0) { + n = k + j * m_kk + i * m_kk * m_kk; + sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n]; + } + } + } + } + + + if (charge[j] > 0.0) { + // sum over all cations + if (j != i) { + sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]); + } + for (k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + // two inner sums over anions + + n = k + j * m_kk + i * m_kk * m_kk; + sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n]; + /* + * Find the counterIJ for the j,k interaction + */ + n = m_kk*j + k; + counterIJ2 = m_CounterIJ[n]; + sum4 = sum4 + (fabs(charge[i])* + molality[j]*molality[k]*CMX_LL[counterIJ2]); + } + } + } + + /* + * Handle neutral j species + */ + if (charge[j] == 0) { + sum5 = sum5 + molality[j]*2.0*m_Lambda_ij_LL(j,i); + } + } + /* + * Add all of the contributions up to yield the log of the + * solute activity coefficients (molality scale) + */ + m_d2lnActCoeffMolaldT2[i] = + zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5; +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s d2lngammadT2[i]=%10.6f \n", + sni.c_str(), m_d2lnActCoeffMolaldT2[i]); + printf(" %12g %12g %12g %12g %12g %12g\n", + zsqd2FdT2, sum1, sum2, sum3, sum4, sum5); + } +#endif + } + + + /* + * ------ SUBSECTION FOR CALCULATING THE d2ACTCOEFFdT2 FOR ANIONS ------ + * + */ + if (charge[i] < 0 ) { + // species i is an anion (negative) + zsqd2FdT2 = charge[i]*charge[i]*d2FdT2; + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + for (j = 1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + + /* + * For Anions, do the cation interactions. + */ + if (charge[j] > 0) { + sum1 = sum1 + molality[j]* + (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]); + if (j < m_kk-1) { + for (k = j+1; k < m_kk; k++) { + // an inner sum over all cations + if (charge[k] > 0) { + n = k + j * m_kk + i * m_kk * m_kk; + sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n]; + } + } + } + } + + /* + * For Anions, do the other anion interactions. + */ + if (charge[j] < 0.0) { + // sum over all anions + if (j != i) { + sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]); + } + for (k = 1; k < m_kk; k++) { + if (charge[k] > 0.0) { + // two inner sums over cations + n = k + j * m_kk + i * m_kk * m_kk; + sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n]; + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*j + k; + counterIJ2 = m_CounterIJ[n]; + sum4 = sum4 + + (fabs(charge[i])* + molality[j]*molality[k]*CMX_LL[counterIJ2]); + } + } + } + + /* + * for Anions, do the neutral species interaction + */ + if (charge[j] == 0.0) { + sum5 = sum5 + molality[j]*2.0*m_Lambda_ij_LL(j,i); + } + } + m_d2lnActCoeffMolaldT2[i] = + zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5; +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s d2lngammadT2[i]=%10.6f\n", + sni.c_str(), m_d2lnActCoeffMolaldT2[i]); + printf(" %12g %12g %12g %12g %12g %12g\n", + zsqd2FdT2, sum1, sum2, sum3, sum4, sum5); + } +#endif + } + /* + * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF ------- + * ------ -> equations agree with my notes, + * -> Equations agree with Pitzer, + */ + if (charge[i] == 0.0 ) { + sum1 = 0.0; + for (j = 1; j < m_kk; j++) { + sum1 = sum1 + molality[j]*2.0*m_Lambda_ij_LL(i,j); + } + m_d2lnActCoeffMolaldT2[i] = sum1; +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s d2lngammadT2[i]=%10.6f \n", + sni.c_str(), m_d2lnActCoeffMolaldT2[i]); + } +#endif + } + + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 9: \n"); + } +#endif + + /* + * ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 --------- + * + */ + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + double sum6 = 0.0; + /* + * term1 is the temperature derivative of the + * DH term in the osmotic coefficient expression + * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer + * implementations. + * Is = Ionic strength on the molality scale (units of (gmol/kg)) + * Aphi = A_Debye / 3 (units of sqrt(kg/gmol)) + */ + term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is)); + + for (j = 1; j < m_kk; j++) { + /* + * Loop Over Cations + */ + if (charge[j] > 0.0) { + for (k = 1; k < m_kk; k++){ + if (charge[k] < 0.0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + + sum1 = sum1 + molality[j]*molality[k]* + (BphiMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]); + } + } + + for (k = j+1; k < m_kk; k++) { + if (j == (m_kk-1)) { + // we should never reach this step + printf("logic error 1 in Step 9 of hmw_act"); + exit(1); + } + if (charge[k] > 0.0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + * between 2 cations. + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + sum2 = sum2 + molality[j]*molality[k]*Phiphi_LL[counterIJ]; + for (m = 1; m < m_kk; m++) { + if (charge[m] < 0.0) { + // species m is an anion + n = m + k * m_kk + j * m_kk * m_kk; + sum2 = sum2 + + molality[j]*molality[k]*molality[m]*psi_ijk_LL[n]; + } + } + } + } + } + + /* + * Loop Over Anions + */ + if (charge[j] < 0) { + for (k = j+1; k < m_kk; k++) { + if (j == m_kk-1) { + // we should never reach this step + printf("logic error 2 in Step 9 of hmw_act"); + exit(1); + } + if (charge[k] < 0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + * between two anions + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + + sum3 = sum3 + molality[j]*molality[k]*Phiphi_LL[counterIJ]; + for (m = 1; m < m_kk; m++) { + if (charge[m] > 0.0) { + n = m + k * m_kk + j * m_kk * m_kk; + sum3 = sum3 + + molality[j]*molality[k]*molality[m]*psi_ijk_LL[n]; + } + } + } + } + } + + /* + * Loop Over Neutral Species + */ + if (charge[j] == 0) { + for (k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + sum4 = sum4 + molality[j]*molality[k]*m_Lambda_ij_LL(j,k); + } + if (charge[k] > 0.0) { + sum5 = sum5 + molality[j]*molality[k]*m_Lambda_ij_LL(j,k); + } + if (charge[k] == 0.0) { + if (k > j) { + sum6 = sum6 + molality[j]*molality[k]*m_Lambda_ij_LL(j,k); + } else if (k == j) { + sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_ij_LL(j,k); + } + } + } + } + } + sum_m_phi_minus_1 = 2.0 * + (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6); + /* + * Calculate the osmotic coefficient from + * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i) + */ + if (molalitysum > 1.0E-150) { + d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum); + } else { + d2_osmotic_coef_dT2 = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f " + "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n", + term1, sum1, sum2, sum3, sum4, sum5); + printf(" sum_m_phi_minus_1=%10.6f d2_osmotic_coef_dT2=%10.6f\n", + sum_m_phi_minus_1, d2_osmotic_coef_dT2); + } + + if (m_debugCalc) { + printf(" Step 10: \n"); + } +#endif + d2_lnwateract_dT2 = -(m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2; + + /* + * In Cantera, we define the activity coefficient of the solvent as + * + * act_0 = actcoeff_0 * Xmol_0 + * + * We have just computed act_0. However, this routine returns + * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). + */ + m_d2lnActCoeffMolaldT2[0] = d2_lnwateract_dT2; + +#ifdef DEBUG_HKM + if (m_debugCalc) { + d2_wateract_dT2 = exp(d2_lnwateract_dT2); + printf(" d2_ln_a_water_dT2 = %10.6f d2_a_water_dT2=%10.6f\n\n", + d2_lnwateract_dT2, d2_wateract_dT2); + } +#endif + } + + /***********************************************************************************************/ + + /** + * s_Pitzer_dlnMolalityActCoeff_dP() (private, const ) + * + * Using internally stored values, this function calculates + * the pressure 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 HMWSoln::s_Pitzer_dlnMolalityActCoeff_dP() const { + + for (int k = 0; k < m_kk; k++) { + m_dlnActCoeffMolaldP[k] = 0.0; + } + s_update_dlnMolalityActCoeff_dP(); + } + + /** + * s_update_dlnMolalityActCoeff_dP() (private, const ) + * + * Using internally stored values, this function calculates + * the pressure derivative of the logarithm of the + * activity coefficient for all species in the mechanism. + * This is an internal routine + * + * We assume that the activity coefficients are current. + * + * It may be assumed that the + * Pitzer activity coefficient and first deriv routine are called immediately + * preceding the calling of this routine. Therefore, some + * quantities do not need to be recalculated in this routine. + * + * solvent activity coefficient is on the molality + * scale. It's derivatives are too. + */ + void HMWSoln::s_update_dlnMolalityActCoeff_dP() const { + + + /* + * HKM -> Assumption is made that the solvent is + * species 0. + */ +#ifdef DEBUG_HKM + m_debugCalc = 0; +#endif + if (m_indexSolvent != 0) { + printf("Wrong index solvent value!\n"); + exit(-1); + } + + double d_wateract_dP; + string sni, snj, snk; + + const double *molality = DATA_PTR(m_molalities); + const double *charge = DATA_PTR(m_speciesCharge); + const double *beta0MX_P = DATA_PTR(m_Beta0MX_ij_P); + const double *beta1MX_P = DATA_PTR(m_Beta1MX_ij_P); + const double *beta2MX_P = DATA_PTR(m_Beta2MX_ij_P); + const double *CphiMX_P = DATA_PTR(m_CphiMX_ij_P); + const double *thetaij_P = DATA_PTR(m_Theta_ij_P); + const double *alphaMX = DATA_PTR(m_Alpha1MX_ij); + const double *psi_ijk_P = DATA_PTR(m_Psi_ijk_P); + + /* + * Local variables defined by Coltrin + */ + double etheta[5][5], etheta_prime[5][5], sqrtIs; + /* + * Molality based ionic strength of the solution + */ + double Is = 0.0; + /* + * Molarcharge of the solution: In Pitzer's notation, + * this is his variable called "Z". + */ + double molarcharge = 0.0; + /* + * molalitysum is the sum of the molalities over all solutes, + * even those with zero charge. + */ + double molalitysum = 0.0; + + double *g = DATA_PTR(m_gfunc_IJ); + double *hfunc = DATA_PTR(m_hfunc_IJ); + double *BMX_P = DATA_PTR(m_BMX_IJ_P); + double *BprimeMX_P= DATA_PTR(m_BprimeMX_IJ_P); + double *BphiMX_P = DATA_PTR(m_BphiMX_IJ_P); + double *Phi_P = DATA_PTR(m_Phi_IJ_P); + double *Phiprime = DATA_PTR(m_Phiprime_IJ); + double *Phiphi_P = DATA_PTR(m_PhiPhi_IJ_P); + double *CMX_P = DATA_PTR(m_CMX_IJ_P); + + double x, g12rooti, gprime12rooti; + double Aphi, dFdP, zsqdFdP; + double sum1, sum2, sum3, sum4, sum5, term1; + double sum_m_phi_minus_1, d_osmotic_coef_dP, d_lnwateract_dP; + + int z1, z2; + int n, i, j, k, m, counterIJ, counterIJ2; + + double currTemp = temperature(); + double currPres = pressure(); + +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf("\n Debugging information from " + "s_Pitzer_dlnMolalityActCoeff_dP()\n"); + } +#endif + /* + * Make sure the counter variables are setup + */ + counterIJ_setup(); + + /* + * ---------- Calculate common sums over solutes --------------------- + */ + for (n = 1; n < m_kk; n++) { + // ionic strength + Is += charge[n] * charge[n] * molality[n]; + // total molar charge + molarcharge += fabs(charge[n]) * molality[n]; + molalitysum += molality[n]; + } + Is *= 0.5; + if (Is > m_maxIionicStrength) { + Is = m_maxIionicStrength; + } + /* + * Store the ionic molality in the object for reference. + */ + m_IionicMolality = Is; + sqrtIs = sqrt(Is); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 1: \n"); + printf(" ionic strenth = %14.7le \n total molar " + "charge = %14.7le \n", Is, molarcharge); + } +#endif + + /* + * The following call to calc_lambdas() calculates all 16 elements + * of the elambda and elambda1 arrays, given the value of the + * ionic strength (Is) + */ + calc_lambdas(Is); + + + /* + * ----- Step 2: Find the coefficients E-theta and ------------------- + * E-thetaprime for all combinations of positive + * unlike charges up to 4 + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 2: \n"); + } +#endif + for (z1 = 1; z1 <=4; z1++) { + for (z2 =1; z2 <=4; z2++) { + calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", + z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); + } +#endif + } + } + +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 3: \n"); + printf(" Species Species g(x) " + " hfunc(x) \n"); + } +#endif + + + /* + * + * calculate g(x) and hfunc(x) for each cation-anion pair MX + * In the original literature, hfunc, was called gprime. However, + * it's not the derivative of g(x), so I renamed it. + */ + for (i = 1; i < (m_kk - 1); i++) { + for (j = (i+1); j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * Only loop over oppositely charge species + */ + if (charge[i]*charge[j] < 0) { + /* + * x is a reduced function variable + */ + x = sqrtIs * alphaMX[counterIJ]; + if (x > 1.0E-100) { + g[counterIJ] = 2.0*(1.0-(1.0 + x) * exp(-x)) / (x*x); + hfunc[counterIJ] = -2.0* + (1.0-(1.0 + x + 0.5*x*x) * exp(-x)) / (x*x); + } + else { + g[counterIJ] = 0.0; + hfunc[counterIJ] = 0.0; + } + } + else { + g[counterIJ] = 0.0; + hfunc[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(), + g[counterIJ], hfunc[counterIJ]); + } +#endif + } + } + + + /* + * ------- SUBSECTION TO CALCULATE BMX_L, BprimeMX_L, BphiMX_L ---------- + * ------- These are now temperature derivatives of the + * previously calculated quantities. + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 4: \n"); + printf(" Species Species BMX " + "BprimeMX BphiMX \n"); + } +#endif + x = 12.0 * sqrtIs; + if (x > 1.0E-100) { + g12rooti = 2.0*(1.0-(1.0 + x) * exp(-x)) / (x*x); + gprime12rooti = -2.0*(1.0-(1.0 + x + 0.5*x*x) * exp(-x)) / (x*x); + } else { + g12rooti = 0.0; + gprime12rooti = 0.0; + } + + for (i = 1; i < m_kk - 1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0.0) { + BMX_P[counterIJ] = beta0MX_P[counterIJ] + + beta1MX_P[counterIJ] * g[counterIJ] + + beta2MX_P[counterIJ] * g12rooti; +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf("%d %g: %g %g %g\n", + counterIJ, BMX_P[counterIJ], beta0MX_P[counterIJ], + beta1MX_P[counterIJ], g[counterIJ]); + } +#endif + if (Is > 1.0E-150) { + BprimeMX_P[counterIJ] = (beta1MX_P[counterIJ] * hfunc[counterIJ]/Is + + beta2MX_P[counterIJ] * gprime12rooti/Is); + } else { + BprimeMX_P[counterIJ] = 0.0; + } + BphiMX_P[counterIJ] = BMX_P[counterIJ] + Is*BprimeMX_P[counterIJ]; + } + else { + BMX_P[counterIJ] = 0.0; + BprimeMX_P[counterIJ] = 0.0; + BphiMX_P[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %11.7f %11.7f %11.7f \n", + sni.c_str(), snj.c_str(), + BMX_P[counterIJ], BprimeMX_P[counterIJ], BphiMX_P[counterIJ]); + } +#endif + } + } + + + /* + * --------- SUBSECTION TO CALCULATE CMX_L ---------- + * --------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 5: \n"); + printf(" Species Species CMX \n"); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0.0) { + CMX_P[counterIJ] = CphiMX_P[counterIJ]/ + (2.0* sqrt(fabs(charge[i]*charge[j]))); + } + else { + CMX_P[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(), + CMX_P[counterIJ]); + } +#endif + } + } + + /* + * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- + * -------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 6: \n"); + printf(" Species Species Phi_ij " + " Phiprime_ij Phi^phi_ij \n"); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] > 0) { + z1 = (int) fabs(charge[i]); + z2 = (int) fabs(charge[j]); + //Phi[counterIJ] = thetaij_L[counterIJ] + etheta[z1][z2]; + Phi_P[counterIJ] = thetaij_P[counterIJ]; + //Phiprime[counterIJ] = etheta_prime[z1][z2]; + Phiprime[counterIJ] = 0.0; + Phiphi_P[counterIJ] = Phi_P[counterIJ] + Is * Phiprime[counterIJ]; + } + else { + Phi_P[counterIJ] = 0.0; + Phiprime[counterIJ] = 0.0; + Phiphi_P[counterIJ] = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + snj = speciesName(j); + printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", + sni.c_str(), snj.c_str(), + Phi_P[counterIJ], Phiprime[counterIJ], Phiphi_P[counterIJ] ); + } +#endif + } + } + + /* + * ----------- SUBSECTION FOR CALCULATION OF dFdT --------------------- + */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 7: \n"); + } +#endif + // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol)) + // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data + // ( A * ln(10) = 1.17593) + // Aphi = A_Debye_Huckel * 2.30258509 / 3.0; + Aphi = m_A_Debye / 3.0; + + double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres); + double dAphidP = dA_DebyedP /3.0; +#ifdef DEBUG_HKM + //dAphidT = 0.0; +#endif + //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + // + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); + //dAphidT = Al / (4.0 * GasConstant * T * T); + dFdP = -dAphidP * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" initial value of dFdP = %10.6f \n", dFdP ); + } +#endif + for (i = 1; i < m_kk-1; i++) { + for (j = i+1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + /* + * both species have a non-zero charge, and one is positive + * and the other is negative + */ + if (charge[i]*charge[j] < 0) { + dFdP = dFdP + molality[i]*molality[j] * BprimeMX_P[counterIJ]; + } + /* + * Both species have a non-zero charge, and they + * have the same sign, e.g., both positive or both negative. + */ + if (charge[i]*charge[j] > 0) { + dFdP = dFdP + molality[i]*molality[j] * Phiprime[counterIJ]; + } +#ifdef DEBUG_HKM + if (m_debugCalc) printf(" dFdP = %10.6f \n", dFdP); +#endif + } + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 8: \n"); + } +#endif + + + for (i = 1; i < m_kk; i++) { + + /* + * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS ----- + * -- + */ + if (charge[i] > 0 ) { + // species i is the cation (positive) to calc the actcoeff + zsqdFdP = charge[i]*charge[i]*dFdP; + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + for (j = 1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + + if (charge[j] < 0.0) { + // sum over all anions + sum1 = sum1 + molality[j]* + (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]); + if (j < m_kk-1) { + /* + * This term is the ternary interaction involving the + * non-duplicate sum over double anions, j, k, with + * respect to the cation, i. + */ + for (k = j+1; k < m_kk; k++) { + // an inner sum over all anions + if (charge[k] < 0.0) { + n = k + j * m_kk + i * m_kk * m_kk; + sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n]; + } + } + } + } + + + + if (charge[j] > 0.0) { + // sum over all cations + if (j != i) { + sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]); + } + for (k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + // two inner sums over anions + + n = k + j * m_kk + i * m_kk * m_kk; + sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n]; + /* + * Find the counterIJ for the j,k interaction + */ + n = m_kk*j + k; + counterIJ2 = m_CounterIJ[n]; + sum4 = sum4 + (fabs(charge[i])* + molality[j]*molality[k]*CMX_P[counterIJ2]); + } + } + } + + /* + * Handle neutral j species + */ + if (charge[j] == 0) { + sum5 = sum5 + molality[j]*2.0*m_Lambda_ij_L(j,i); + } + } + + /* + * Add all of the contributions up to yield the log of the + * solute activity coefficients (molality scale) + */ + m_dlnActCoeffMolaldP[i] = + zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5; + +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s lngamma[i]=%10.6f \n", + sni.c_str(), m_dlnActCoeffMolaldP[i]); + printf(" %12g %12g %12g %12g %12g %12g\n", + zsqdFdP, sum1, sum2, sum3, sum4, sum5); + } +#endif + } + + /* + * ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR ANIONS ------ + * + */ + if (charge[i] < 0 ) { + // species i is an anion (negative) + zsqdFdP = charge[i]*charge[i]*dFdP; + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + for (j = 1; j < m_kk; j++) { + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*i + j; + counterIJ = m_CounterIJ[n]; + + /* + * For Anions, do the cation interactions. + */ + if (charge[j] > 0) { + sum1 = sum1 + molality[j]* + (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]); + if (j < m_kk-1) { + for (k = j+1; k < m_kk; k++) { + // an inner sum over all cations + if (charge[k] > 0) { + n = k + j * m_kk + i * m_kk * m_kk; + sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n]; + } + } + } + } + + /* + * For Anions, do the other anion interactions. + */ + if (charge[j] < 0.0) { + // sum over all anions + if (j != i) { + sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]); + } + for (k = 1; k < m_kk; k++) { + if (charge[k] > 0.0) { + // two inner sums over cations + n = k + j * m_kk + i * m_kk * m_kk; + sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n]; + /* + * Find the counterIJ for the symmetric binary interaction + */ + n = m_kk*j + k; + counterIJ2 = m_CounterIJ[n]; + sum4 = sum4 + + (fabs(charge[i])* + molality[j]*molality[k]*CMX_P[counterIJ2]); + } + } + } + + /* + * for Anions, do the neutral species interaction + */ + if (charge[j] == 0.0) { + sum5 = sum5 + molality[j]*2.0*m_Lambda_ij_L(j,i); + } + } + m_dlnActCoeffMolaldP[i] = + zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5; +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s lndactcoeffmolaldP[i]=%10.6f \n", + sni.c_str(), m_dlnActCoeffMolaldP[i]); + printf(" %12g %12g %12g %12g %12g %12g\n", + zsqdFdP, sum1, sum2, sum3, sum4, sum5); + } +#endif + } + + + /* + * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF ------- + * ------ -> equations agree with my notes, + * -> Equations agree with Pitzer, + */ + if (charge[i] == 0.0 ) { + sum1 = 0.0; + for (j = 1; j < m_kk; j++) { + sum1 = sum1 + molality[j]*2.0*m_Lambda_ij_L(i,j); + } + m_dlnActCoeffMolaldP[i] = sum1; +#ifdef DEBUG_HKM + if (m_debugCalc) { + sni = speciesName(i); + printf(" %-16s dlnActCoeffMolaldP[i]=%10.6f \n", + sni.c_str(), m_dlnActCoeffMolaldP[i]); + } +#endif + } + + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Step 9: \n"); + } +#endif + + /* + * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT --------- + * + */ + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + sum5 = 0.0; + double sum6 = 0.0; + /* + * term1 is the temperature derivative of the + * DH term in the osmotic coefficient expression + * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer + * implementations. + * Is = Ionic strength on the molality scale (units of (gmol/kg)) + * Aphi = A_Debye / 3 (units of sqrt(kg/gmol)) + */ + term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is)); + + for (j = 1; j < m_kk; j++) { + /* + * Loop Over Cations + */ + if (charge[j] > 0.0) { + for (k = 1; k < m_kk; k++){ + if (charge[k] < 0.0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + + sum1 = sum1 + molality[j]*molality[k]* + (BphiMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]); + } + } + + for (k = j+1; k < m_kk; k++) { + if (j == (m_kk-1)) { + // we should never reach this step + printf("logic error 1 in Step 9 of hmw_act"); + exit(1); + } + if (charge[k] > 0.0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + * between 2 cations. + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + sum2 = sum2 + molality[j]*molality[k]*Phiphi_P[counterIJ]; + for (m = 1; m < m_kk; m++) { + if (charge[m] < 0.0) { + // species m is an anion + n = m + k * m_kk + j * m_kk * m_kk; + sum2 = sum2 + + molality[j]*molality[k]*molality[m]*psi_ijk_P[n]; + } + } + } + } + } + + + /* + * Loop Over Anions + */ + if (charge[j] < 0) { + for (k = j+1; k < m_kk; k++) { + if (j == m_kk-1) { + // we should never reach this step + printf("logic error 2 in Step 9 of hmw_act"); + exit(1); + } + if (charge[k] < 0) { + /* + * Find the counterIJ for the symmetric j,k binary interaction + * between two anions + */ + n = m_kk*j + k; + counterIJ = m_CounterIJ[n]; + + sum3 = sum3 + molality[j]*molality[k]*Phiphi_P[counterIJ]; + for (m = 1; m < m_kk; m++) { + if (charge[m] > 0.0) { + n = m + k * m_kk + j * m_kk * m_kk; + sum3 = sum3 + + molality[j]*molality[k]*molality[m]*psi_ijk_P[n]; + } + } + } + } + } + + /* + * Loop Over Neutral Species + */ + if (charge[j] == 0) { + for (k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + sum4 = sum4 + molality[j]*molality[k]*m_Lambda_ij_P(j,k); + } + if (charge[k] > 0.0) { + sum5 = sum5 + molality[j]*molality[k]*m_Lambda_ij_P(j,k); + } + if (charge[k] == 0.0) { + if (k > j) { + sum6 = sum6 + molality[j]*molality[k]*m_Lambda_ij_P(j,k); + } else if (k == j) { + sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_ij_P(j,k); + } + } + } + } + } + sum_m_phi_minus_1 = 2.0 * + (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6); + + + /* + * Calculate the osmotic coefficient from + * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i) + */ + if (molalitysum > 1.0E-150) { + d_osmotic_coef_dP = 0.0 + (sum_m_phi_minus_1 / molalitysum); + } else { + d_osmotic_coef_dP = 0.0; + } +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f " + "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n", + term1, sum1, sum2, sum3, sum4, sum5); + printf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dP =%10.6f\n", + sum_m_phi_minus_1, d_osmotic_coef_dP); + } + + if (m_debugCalc) { + printf(" Step 10: \n"); + } +#endif + d_lnwateract_dP = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP; + d_wateract_dP = exp(d_lnwateract_dP); + + /* + * In Cantera, we define the activity coefficient of the solvent as + * + * act_0 = actcoeff_0 * Xmol_0 + * + * We have just computed act_0. However, this routine returns + * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). + */ + //double xmolSolvent = moleFraction(m_indexSolvent); + m_dlnActCoeffMolaldP[0] = d_lnwateract_dP; +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" d_ln_a_water_dP = %10.6f d_a_water_dP=%10.6f\n\n", + d_lnwateract_dP, d_wateract_dP); + } +#endif + + + + } + + /***********************************************************************************************/ + + /* + * Calculate the lambda interactions. + * + * Calculate E-lambda terms for charge combinations of like sign, + * using method of Pitzer (1975). + * + * This code snipet is included from Bethke, Appendix 2. + */ + void HMWSoln::calc_lambdas(double is) const { + double aphi, dj, jfunc, jprime, t, x, zprod; + int i, ij, j; + /* + * Coefficients c1-c4 are used to approximate + * the integral function "J"; + * aphi is the Debye-Huckel constant at 25 C + */ + + double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528; + + aphi = 0.392; /* Value at 25 C */ +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" Is = %g\n", is); + } +#endif + if (is < 1.0E-150) { + for (i = 0; i < 17; i++) { + elambda[i] = 0.0; + elambda1[i] = 0.0; + } + return; + } + /* + * Calculate E-lambda terms for charge combinations of like sign, + * using method of Pitzer (1975). Charges up to 4 are calculated. + */ + + for (i=1; i<=4; i++) { + for (j=i; j<=4; j++) { + ij = i*j; + /* + * calculate the product of the charges + */ + zprod = (double)ij; + /* + * calculate Xmn (A1) from Harvie, Weare (1980). + */ + x = 6.0* zprod * aphi * sqrt(is); /* eqn 23 */ + + jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4))); /* eqn 47 */ + + t = c3 * c4 * pow(x,c4); + dj = c1* pow(x,(-c2-1.0)) * (c2+t) * exp(-c3*pow(x,c4)); + jprime = (jfunc/x)*(1.0 + jfunc*dj); + + elambda[ij] = zprod*jfunc / (4.0*is); /* eqn 14 */ + elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is)) + - elambda[ij])/is; +#ifdef DEBUG_HKM + if (m_debugCalc) { + printf(" ij = %d, elambda = %g, elambda1 = %g\n", + ij, elambda[ij], elambda1[ij]); + } +#endif + } + } + } + + /* + * Calculate the etheta interaction. + * This interaction accounts for the mixing effects of like-signed + * ions with different charges. There is fairly extensive literature + * on this effect. See the notes. + * This interaction will be nonzero for species with the same charge. + * + * This code snipet is included from Bethke, Appendix 2. + */ + void HMWSoln::calc_thetas(int z1, int z2, + double *etheta, double *etheta_prime) const { + int i, j; + double f1, f2; + + /* + * Calculate E-theta(i) and E-theta'(I) using method of + * Pitzer (1987) + */ + i = abs(z1); + j = abs(z2); + +#ifdef DEBUG_HKM + if (i > 4 || j > 4) { + printf("we shouldn't be here\n"); + exit(-1); + } +#endif + + if ((i == 0) || (j == 0)) { + printf("ERROR calc_thetas called with one species being neutral\n"); + exit(-1); + } + + /* + * Check to see if the charges are of opposite sign. If they are of + * opposite sign then their etheta interaction is zero. + */ + if (z1*z2 < 0) { + *etheta = 0.0; + *etheta_prime = 0.0; + } + /* + * Actually calculate the interaction. + */ + else { + f1 = (double)i / (2.0 * j); + f2 = (double)j / (2.0 * i); + *etheta = elambda[i*j] - f1*elambda[j*j] - f2*elambda[i*i]; + *etheta_prime = elambda1[i*j] - f1*elambda1[j*j] - f2*elambda1[i*i]; + } + } + + /** + * This routine prints out the input pitzer coefficients for the + * current mechanism + */ + void HMWSoln::printCoeffs() const { + int i, j, k; + string sni, snj; + calcMolalities(); + const double *charge = DATA_PTR(m_speciesCharge); + double *molality = DATA_PTR(m_molalities); + double *moleF = DATA_PTR(m_tmpV); + /* + * Update the coefficients wrt Temperature + * Calculate the derivatives as well + */ + s_updatePitzerCoeffWRTemp(2); + getMoleFractions(moleF); + + printf("Index Name MoleF Molality Charge\n"); + for (k = 0; k < m_kk; k++) { + sni = speciesName(k); + printf("%2d %-16s %14.7le %14.7le %5.1f \n", + k, sni.c_str(), moleF[k], molality[k], charge[k]); + } + + printf("\n Species Species beta0MX " + "beta1MX beta2MX CphiMX alphaMX thetaij \n"); + for (i = 1; i < m_kk - 1; i++) { + sni = speciesName(i); + for (j = i+1; j < m_kk; j++) { + snj = speciesName(j); + int n = i * m_kk + j; + int ct = m_CounterIJ[n]; + printf(" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n", + sni.c_str(), snj.c_str(), + m_Beta0MX_ij[ct], m_Beta1MX_ij[ct], + m_Beta2MX_ij[ct], m_CphiMX_ij[ct], + m_Alpha1MX_ij[ct], m_Theta_ij[ct] ); + + + } + } + + printf("\n Species Species Species " + "psi \n"); + for (i = 1; i < m_kk; i++) { + sni = speciesName(i); + for (j = 1; j < m_kk; j++) { + snj = speciesName(j); + for (k = 1; k < m_kk; k++) { + string snk = speciesName(k); + int n = k + j * m_kk + i * m_kk * m_kk; + if (m_Psi_ijk[n] != 0.0) { + printf(" %-16s %-16s %-16s %9.5f \n", + sni.c_str(), snj.c_str(), + snk.c_str(), m_Psi_ijk[n]); + } + } + } + } + } + + + /*****************************************************************************/ +} +/*****************************************************************************/ diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h new file mode 100644 index 000000000..d21e38bc7 --- /dev/null +++ b/Cantera/src/thermo/HMWSoln.h @@ -0,0 +1,1492 @@ +/** + * @file HMWSoln.h + * + * Header file for Pitzer activity coefficient implementation + */ +/* + * 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_HMWSOLN_H +#define CT_HMWSOLN_H + +#include "MolalityVPSSTP.h" +#include "electrolytes.h" + +namespace Cantera { + + /** + * @defgroup thermoprops Thermodynamic Properties + * + * These classes are used to compute thermodynamic properties. + */ + + /** + * HMWSoln.h + * + * Major Parameters: + * The form of the Pitzer expression refers to the + * form of the Gibbs free energy expression. The temperature + * dependence of the Pitzer coefficients are handled by + * another parameter. + * + * m_formPitzer = Form of the Pitzer expression + * + * PITZERFORM_BASE = 0 + * + * Only one form is supported atm. This parameter is included for + * future expansion. + * + */ +#define PITZERFORM_BASE 0 + + + /* + * Formulations for the temperature dependence of the Pitzer + * coefficients. Note, the temperature dependence of the + * Gibbs free energy also depends on the temperature dependence + * of the standard state and the temperature dependence of the + * Debye-Huckel constant, which includes the dielectric constant + * and the density. Therefore, this expression defines only part + * of the temperature dependence for the mixture thermodynamic + * functions. + * + * PITZER_TEMP_CONSTANT + * All coefficients are considered constant wrt temperature + * PITZER_TEMP_LINEAR + * All coefficients are assumed to have a linear dependence + * wrt to temperature. + * PITZER_TEMP_COMPLEX1 + * All coefficnets are assumed to have a complex functional + * based dependence wrt temperature; See: + * (Silvester, Pitzer, J. Phys. Chem. 81, 19 1822 (1977)). + * + * beta0 = q0 + q3(1/T - 1/Tr) + q4(ln(T/Tr)) + + * q1(T - Tr) + q2(T**2 - Tr**2) + */ +#define PITZER_TEMP_CONSTANT 0 +#define PITZER_TEMP_LINEAR 1 +#define PITZER_TEMP_COMPLEX1 2 + + /* + * Acceptable ways to calculate the value of A_Debye + */ +#define A_DEBYE_CONST 0 +#define A_DEBYE_WATER 1 + + class WaterProps; + class WaterPDSS; + + /** + * Definition of the HMWSoln object + */ + class HMWSoln : public MolalityVPSSTP { + + public: + + /// Constructors + HMWSoln(); + + HMWSoln(const HMWSoln &); + HMWSoln& operator=(const HMWSoln&); + + HMWSoln(string inputFile, string id = ""); + HMWSoln(XML_Node& phaseRef, string id = ""); + + /** + * This is a special constructor, used to replicate test problems + * during the initial verification of the object + */ + HMWSoln(int testProb); + + /// Destructor. + virtual ~HMWSoln(); + + + 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; + + /** + * Excess molar enthalpy of the solution from + * the mixing process. Units: J/ kmol. + * + * Note this is kmol of the total solution. + */ + virtual doublereal relative_enthalpy() const; + + /** + * Excess molar enthalpy of the solution from + * the mixing process on a molality basis. + * Units: J/ (kmol add salt). + * + * Note this is kmol of the guessed at salt composition + */ + virtual doublereal relative_molal_enthalpy() 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); + + /** + * Calculate the density of the mixture using the partial + * molar volumes and mole fractions as input + * + * The formula for this is + * + * \f[ + * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} + * \f] + * + * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are + * the molecular weights, and \f$V_k\f$ are the pure species + * molar volumes. + * + * Note, the basis behind this formula is that in an ideal + * solution the partial molar volumes are equal to the pure + * species molar volumes. We have additionally specified + * in this class that the pure species molar volumes are + * independent of temperature and pressure. + * + * NOTE: This is a non-virtual function, which is not a + * member of the ThermoPhase base class. + */ + void calcDensity(); + + /** + * Overwritten setDensity() function is necessary because the + * density is not an indendent variable. + * + * This function will now throw an error condition + * + * @internal May have to adjust the strategy here to make + * the eos for these materials slightly compressible, in order + * to create a condition where the density is a function of + * the pressure. + * + * This function will now throw an error condition. + * + * NOTE: This is an overwritten function from the State.h + * class + */ + void setDensity(doublereal rho); + + /** + * Overwritten setMolarDensity() function is necessary because the + * density is not an indendent variable. + * + * This function will now throw an error condition. + * + * NOTE: This is an overwritten function from the State.h + * class + */ + void setMolarDensity(doublereal rho); + + /** + * Overwritten setTemperature(double) from State.h. This + * function sets the temperature, and makes sure that + * the value propagates to underlying objects. + */ + virtual void setTemperature(doublereal temp); + + /** + * 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; + + 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; } + + /* + * constructPhaseFile() (virtual from HMWSoln) + * + * Import, construct, and initialize a HMWSoln phase + * specification from an XML tree into the current object. + * + * This routine is a precursor to constructPhaseXML(XML_Node*) + * routine, which does most of the work. + */ + virtual void constructPhaseFile(string inputFile, string id); + + /* + * constructPhaseXML (virtual from HMWSoln) + * + * This is the main routine for constructing the phase. + * + * Most of the work is carried out by the cantera base + * routine, importPhase(). That routine imports all of the + * species and element data, including the standard states + * of the species. + * + * Then, In this routine, we read the information + * particular to the specification of the activity + * coefficient model for the Pitzer parameterization. + */ + virtual void constructPhaseXML(XML_Node& phaseNode, string id); + + /** + * @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(); + + /* + * initThermoXML() (virtual from ThermoPhase) + * + * This gets called from importPhase(). It processes the XML file + * after the species are set up. This is the main routine for + * reading in activity coefficient parameters. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + virtual void initThermoXML(XML_Node& phaseNode, 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; + + + /** + * Value of the Debye Huckel constant as a function of temperature + * and pressure. + * + * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) + * + * Units = sqrt(kg/gmol) + */ + virtual double A_Debye_TP(double temperature = -1.0, + double pressure = -1.0) const; + + /** + * Value of the derivative of the Debye Huckel constant with + * respect to temperature as a function of temperature + * and pressure. + * + * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) + * + * Units = sqrt(kg/gmol) + */ + virtual double dA_DebyedT_TP(double temperature = -1.0, + double pressure = -1.0) const; + + /** + * Value of the derivative of the Debye Huckel constant with + * respect to pressure, as a function of temperature + * and pressure. + * + * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) + * + * Units = sqrt(kg/gmol) + */ + virtual double dA_DebyedP_TP(double temperature = -1.0, + double pressure = -1.0) const; + + /** + * Return Pitzer's definition of A_L. This is basically the + * derivative of the A_phi multiplied by 4 R T**2 + * + * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) + * dA_phidT = d(A_Debye)/dT / 3.0 + * A_L = dA_phidT * (4 * R * T * T) + * + * Units = sqrt(kg/gmol) (RT) + * + */ + double ADebye_L(double temperature = -1.0, + double pressure = -1.0) const; + + + /** + * Return Pitzer's definition of A_J. This is basically the + * temperature derivative of A_L, and the second derivative + * of A_phi + * + * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) + * dA_phidT = d(A_Debye)/dT / 3.0 + * A_J = 2 A_L/T + 4 * R * T * T * d2(A_phi)/dT2 + * + * Units = sqrt(kg/gmol) (R) + */ + double ADebye_J(double temperature = -1.0, + double pressure = -1.0) const; + /** + * Return Pitzer's definition of A_V. This is the + * derivative wrt pressure of A_phi multiplied by - 4 R T + * + * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) + * dA_phidT = d(A_Debye)/dP / 3.0 + * A_V = - dA_phidP * (4 * R * T) + * + * Units = sqrt(kg/gmol) (RT) / Pascal + * + */ + double ADebye_V(double temperature = -1.0, + double pressure = -1.0) const; + /** + * Value of the 2nd derivative of the Debye Huckel constant with + * respect to temperature as a function of temperature + * and pressure. + * + * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) + * + * Units = sqrt(kg/gmol) + */ + 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; + + /** + * + * formPitzer(): + * + * Returns the form of the Pitzer parameterization used + */ + int formPitzer() const { return m_formPitzer; } + + /** + * Print out all of the input coefficients. + */ + void printCoeffs () const; + + + //@} + + protected: + + /** + * This is the form of the Pitzer 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: + * + * PITZERFORM_BASE = 0 (only one supported atm) + * + */ + int m_formPitzer; + + /** + * This is the form of the temperature dependence of Pitzer + * parameterization used in the model. + * + * PITZER_TEMP_CONSTANT 0 + * PITZER_TEMP_LINEAR 1 + * PITZER_TEMP_COMPLEX1 2 + */ + int m_formPitzerTemp; + + /** + * 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. This is now the independent variable + * as it must be for multicomponent solutions. + */ + 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 + * Associated Salts, if present in the mechanism, + * don't contribute to the value of the ionic strength + * in this version of the Ionic strength. + */ + mutable double m_IionicMolality; + + /** + * Maximum value of the ionic strength allowed in the + * calculation of the activity coefficients. + */ + double m_maxIionicStrength; + + /** + * Reference Temperature for the Pitzer formulations. + */ + double m_TempPitzerRef; + + protected: + /** + * Stoichiometric ionic strength on the molality scale. + * This differs from m_IionicMolality in the sense that + * associated salts are treated as unassociated salts, + * when calculating the Ionic strength by this method. + */ + mutable double m_IionicMolalityStoich; + + public: + /** + * Form of the constant outside the Debye-Huckel term + * called A. It's normally a function of temperature + * and pressure. However, it can be set from the + * input file in order to aid in numerical comparisons. + * Acceptable forms: + * + * A_DEBYE_CONST 0 + * A_DEBYE_WATER 1 + * + * The A_DEBYE_WATER form may be used for water solvents + * with needs to cover varying temperatures and pressures. + * Note, the dielectric constant of water is a relatively + * strong function of T, and its variability must be + * accounted for, + */ + mutable int m_form_A_Debye; + + protected: + /** + * A_Debye -> this expression appears on the top of the + * ln actCoeff term in the general Debye-Huckel + * expression + * It depends on temperature. And, therefore, + * most be recalculated whenever T or P changes. + * + * A_Debye = (F e B_Debye) / (8 Pi epsilon R T) + * + * where B_Debye = F / sqrt(epsilon R T/2) + * (dw/1000)^(1/2) + * + * A_Debye = (1/ (8 Pi)) (2 Pi * Na * dw/1000)^(1/2) + * (e * e / (epsilon * kb * T))^(3/2) + * + * Units = sqrt(kg/gmol) + * + * Nominal value = 1.172576 sqrt(kg/gmol) + * based on: + * epsilon/epsilon_0 = 78.54 + * (water at 25C) + * epsilon_0 = 8.854187817E-12 C2 N-1 m-2 + * e = 1.60217653 E-19 C + * 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 + * dw = C_0 * M_0 (density of water) (kg/m3) + * = 1.0E3 at 25C + */ + 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; + + /** + * Water standard state -> derived from the + * equation of state for water. + */ + WaterPDSS *m_waterSS; + double m_densWaterSS; + + /** + * Water property calculator + */ + WaterProps *m_waterProps; + + /** + * 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 Pitzer/HMW formulation. + * Beta0_ij[i][j] is the value of the Beta0 coefficient + * for the ij salt. It will be nonzero iff i and j are + * both charged and have opposite sign. The array is also + * symmetric. + * counterIJ where counterIJ = m_counterIJ[i][j] + * is used to access this array. + */ + mutable vector_fp m_Beta0MX_ij; + mutable vector_fp m_Beta0MX_ij_L; + mutable vector_fp m_Beta0MX_ij_LL; + mutable vector_fp m_Beta0MX_ij_P; + mutable Array2D m_Beta0MX_ij_coeff; + + /** + * Array of 2D data used in the Pitzer/HMW formulation. + * Beta1_ij[i][j] is the value of the Beta1 coefficient + * for the ij salt. It will be nonzero iff i and j are + * both charged and have opposite sign. The array is also + * symmetric. + * counterIJ where counterIJ = m_counterIJ[i][j] + * is used to access this array. + */ + mutable vector_fp m_Beta1MX_ij; + mutable vector_fp m_Beta1MX_ij_L; + mutable vector_fp m_Beta1MX_ij_LL; + mutable vector_fp m_Beta1MX_ij_P; + mutable Array2D m_Beta1MX_ij_coeff; + + /** + * Array of 2D data used in the Pitzer/HMW formulation. + * Beta2_ij[i][j] is the value of the Beta2 coefficient + * for the ij salt. It will be nonzero iff i and j are + * both charged and have opposite sign, and i and j + * both have charges of 2 or more. The array is also + * symmetric. + * counterIJ where counterIJ = m_counterIJ[i][j] + * is used to access this array. + */ + vector_fp m_Beta2MX_ij; + vector_fp m_Beta2MX_ij_L; + vector_fp m_Beta2MX_ij_LL; + vector_fp m_Beta2MX_ij_P; + + /** + * Array of 2D data used in the Pitzer/HMW formulation. + * Alpha1MX_ij[i][j] is the value of the alpha1 coefficient + * for the ij interaction. It will be nonzero iff i and j are + * both charged and have opposite sign, and i and j + * both have charges of 2 or more. The array is also + * symmetric. + * counterIJ where counterIJ = m_counterIJ[i][j] + * is used to access this array. + */ + vector_fp m_Alpha1MX_ij; + + /** + * Array of 2D data used in the Pitzer/HMW formulation. + * CphiMX_ij[i][j] is the value of the Cphi coefficient + * for the ij interaction. It will be nonzero iff i and j are + * both charged and have opposite sign, and i and j + * both have charges of 2 or more. The array is also + * symmetric. + * counterIJ where counterIJ = m_counterIJ[i][j] + * is used to access this array. + */ + mutable vector_fp m_CphiMX_ij; + mutable vector_fp m_CphiMX_ij_L; + mutable vector_fp m_CphiMX_ij_LL; + mutable vector_fp m_CphiMX_ij_P; + mutable Array2D m_CphiMX_ij_coeff; + + /** + * Array of 2D data used in the Pitzer/HMW formulation. + * Theta_ij[i][j] is the value of the theta coefficient + * for the ij interaction. It will be nonzero for charged + * ions with the same sign. It is symmetric. + * counterIJ where counterIJ = m_counterIJ[i][j] + * is used to access this array. + * + * HKM Recent Pitzer papers have used a functional form + * for Theta_ij, which depends on the ionic strength. + */ + vector_fp m_Theta_ij; + vector_fp m_Theta_ij_L; + vector_fp m_Theta_ij_LL; + vector_fp m_Theta_ij_P; + + /** + * Array of 3D data sed in the Pitzer/HMW formulation. + * Psi_ijk[n] is the value of the psi coefficient for the + * ijk interaction where + * + * n = k + j * m_kk + i * m_kk * m_kk; + * + * It is potentially nonzero everywhere. + * The first two coordinates are symmetric wrt cations, + * and the last two coordinates are symmetric wrt anions. + */ + vector_fp m_Psi_ijk; + vector_fp m_Psi_ijk_L; + vector_fp m_Psi_ijk_LL; + vector_fp m_Psi_ijk_P; + + /* + * Array of 2D data used in the Pitzer/HMW formulation. + * Lambda_ij[i][j] represents the lambda coefficient for the + * ij interaction. This is a general interaction representing + * neutral species. The neutral species occupy the first + * index, i.e., i. The charged species occupy the j coordinate. + * neutral, neutral interactions are also included here. + */ + Array2D m_Lambda_ij; + Array2D m_Lambda_ij_L; + Array2D m_Lambda_ij_LL; + Array2D m_Lambda_ij_P; + + /** + * Logarithm of the activity coefficients on the molality + * scale. + * mutable because we change this if the composition + * or temperature or pressure changes. + */ + mutable vector_fp m_lnActCoeffMolal; + mutable vector_fp m_dlnActCoeffMolaldT; + mutable vector_fp m_d2lnActCoeffMolaldT2; + mutable vector_fp m_dlnActCoeffMolaldP; + + /* + * -------- Temporary Variables Used in the Activity Coeff Calc + */ + + /* + * Set up a counter variable for keeping track of symmetric binary + * interactions amongst the solute species. + * + * n = m_kk*i + j + * m_CounterIJ[n] = counterIJ + */ + mutable array_int m_CounterIJ; + + /** + * This is elambda, MEC + */ + mutable double elambda[17]; + + /** + * This is elambda1, MEC + */ + mutable double elambda1[17]; + + /** + * Various temporary arrays used in the calculation of + * the Pitzer activity coefficents. + * The subscript, L, denotes the same quantity's derivative + * wrt temperature + */ + mutable vector_fp m_gfunc_IJ; + mutable vector_fp m_hfunc_IJ; + mutable vector_fp m_BMX_IJ; + mutable vector_fp m_BMX_IJ_L; + mutable vector_fp m_BMX_IJ_LL; + mutable vector_fp m_BMX_IJ_P; + mutable vector_fp m_BprimeMX_IJ; + mutable vector_fp m_BprimeMX_IJ_L; + mutable vector_fp m_BprimeMX_IJ_LL; + mutable vector_fp m_BprimeMX_IJ_P; + mutable vector_fp m_BphiMX_IJ; + mutable vector_fp m_BphiMX_IJ_L; + mutable vector_fp m_BphiMX_IJ_LL; + mutable vector_fp m_BphiMX_IJ_P; + mutable vector_fp m_Phi_IJ; + mutable vector_fp m_Phi_IJ_L; + mutable vector_fp m_Phi_IJ_LL; + mutable vector_fp m_Phi_IJ_P; + mutable vector_fp m_Phiprime_IJ; + mutable vector_fp m_PhiPhi_IJ; + mutable vector_fp m_PhiPhi_IJ_L; + mutable vector_fp m_PhiPhi_IJ_LL; + mutable vector_fp m_PhiPhi_IJ_P; + mutable vector_fp m_CMX_IJ; + mutable vector_fp m_CMX_IJ_L; + mutable vector_fp m_CMX_IJ_LL; + mutable vector_fp m_CMX_IJ_P; + + mutable vector_fp m_gamma; + + 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 s_updateDHlnMolalityActCoeff() const; + void s_update_lnMolalityActCoeff() const; + public: + void s_Pitzer_dlnMolalityActCoeff_dT() const; + void s_Pitzer_dlnMolalityActCoeff_dP() const; + private: + /** + * This function calculates the temperature derivative of the + * natural logarithm of the molality activity coefficients. + */ + void s_update_dlnMolalityActCoeff_dT() const; + + /** + * This function calcultes the temperature second derivative + * of the natural logarithm of the molality activity + * coefficients. + */ + void s_update_d2lnMolalityActCoeff_dT2() const; + /** + * This function calculates the pressure derivative of the + * natural logarithm of the molality activity coefficients. + */ + void s_update_dlnMolalityActCoeff_dP() const; + + /** + * This function calculates the temperature derivatives + * of the Pitzer coefficients + */ + void s_updatePitzerCoeffWRTemp(int doDerivs = 2) const; + + /** + * This function does the main pitzer coefficient + * calculation + */ + void s_updatePitzerSublnMolalityActCoeff() const; + /* + * Calculate the lambda interactions. + * + * + * Calculate E-lambda terms for charge combinations of like sign, + * using method of Pitzer (1975). + */ + void calc_lambdas(double is) const; + + /** + * Calculate etheta and etheta_prime + * + * This interaction will be nonzero for species with the + * same charge. this routine is not to be called for + * neutral species; it core dumps or error exits. + * + * MEC implementation routine. + * + * @param z1 charge of the first molecule + * @param z2 charge of the second molecule + * @param etheta return pointer containing etheta + * @param etheta_prime Return pointer containing etheta_prime. + * + * This routine uses the internal variables, + * elambda[] and elambda1[]. + * + * There is no prohibition against calling + * + */ + void calc_thetas(int z1, int z2, + double *etheta, double *etheta_prime) const; + + + void counterIJ_setup(void) const; + void readXMLBinarySalt(XML_Node &BinSalt); + void readXMLThetaAnion(XML_Node &BinSalt); + void readXMLThetaCation(XML_Node &BinSalt); + void readXMLPsiCommonAnion(XML_Node &BinSalt); + void readXMLPsiCommonCation(XML_Node &BinSalt); + void readXMLLambdaNeutral(XML_Node &BinSalt); + + + public: + /* + * Turn on copious debug printing + */ + mutable int m_debugCalc; + + }; + +} + +#endif + diff --git a/Cantera/src/thermo/HMWSoln_input.cpp b/Cantera/src/thermo/HMWSoln_input.cpp new file mode 100644 index 000000000..e9d33ec98 --- /dev/null +++ b/Cantera/src/thermo/HMWSoln_input.cpp @@ -0,0 +1,1130 @@ +/** + * @file HMWSoln_input.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$ + */ + +#include "HMWSoln.h" +#include "importCTML.h" +#include "WaterProps.h" +#include "WaterPDSS.h" + +namespace Cantera { + + /** + * 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; + } + + /** + * Process an XML node called "SimpleSaltParameters. + * This node contains all of the parameters necessary to describe + * the Pitzer model for that particular binary salt. + * This function reads the XML file and writes the coefficients + * it finds to an internal data structures. + */ + void HMWSoln::readXMLBinarySalt(XML_Node &BinSalt) { + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + int nParamsFound, i; + vector_fp vParams; + string iName = BinSalt.attrib("cation"); + if (iName == "") { + throw CanteraError("HMWSoln::readXMLBinarySalt", "no cation attrib"); + } + string jName = BinSalt.attrib("anion"); + if (jName == "") { + throw CanteraError("HMWSoln::readXMLBinarySalt", "no anion attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species + */ + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + string ispName = speciesName(iSpecies); + if (charge[iSpecies] <= 0) { + throw CanteraError("HMWSoln::readXMLBinarySalt", "cation charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + string jspName = speciesName(jSpecies); + if (charge[jSpecies] >= 0) { + throw CanteraError("HMWSoln::readXMLBinarySalt", "anion charge problem"); + } + + int n = iSpecies * m_kk + jSpecies; + int counter = m_CounterIJ[n]; + int num = BinSalt.nChildren(); + for (int iChild = 0; iChild < num; iChild++) { + XML_Node &xmlChild = BinSalt.child(iChild); + stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + /* + * Process the binary salt child elements + */ + if (nodeName == "beta0") { + /* + * Get the string containing all of the values + */ + getFloatArray(xmlChild, vParams, false, "", "beta0"); + nParamsFound = vParams.size(); + if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) { + if (nParamsFound != 1) { + throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_Beta0MX_ij[counter] = vParams[0]; + m_Beta0MX_ij_coeff(0,counter) = m_Beta0MX_ij[counter]; + } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) { + if (nParamsFound != 2) { + throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_Beta0MX_ij_coeff(0,counter) = vParams[0]; + m_Beta0MX_ij_coeff(1,counter) = vParams[1]; + m_Beta0MX_ij[counter] = vParams[0]; + } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) { + if (nParamsFound != 5) { + throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + for (i = 0; i < nParamsFound; i++) { + m_Beta0MX_ij_coeff(i, counter) = vParams[i]; + } + m_Beta0MX_ij[counter] = vParams[0]; + } + } + if (nodeName == "beta1") { + + /* + * Get the string containing all of the values + */ + getFloatArray(xmlChild, vParams, false, "", "beta1"); + nParamsFound = vParams.size(); + if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) { + if (nParamsFound != 1) { + throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_Beta1MX_ij[counter] = vParams[0]; + m_Beta1MX_ij_coeff(0,counter) = m_Beta1MX_ij[counter]; + } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) { + if (nParamsFound != 2) { + throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_Beta1MX_ij_coeff(0,counter) = vParams[0]; + m_Beta1MX_ij_coeff(1,counter) = vParams[1]; + m_Beta1MX_ij[counter] = vParams[0]; + } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) { + if (nParamsFound < 3) { + throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + for (i = 0; i < nParamsFound; i++) { + m_Beta1MX_ij_coeff(i, counter) = vParams[i]; + } + m_Beta1MX_ij[counter] = vParams[0]; + } + + } + if (nodeName == "beta2") { + stemp = xmlChild.value(); + m_Beta2MX_ij[counter] = atofCheck(stemp.c_str()); + } + if (nodeName == "cphi") { + /* + * Get the string containing all of the values + */ + getFloatArray(xmlChild, vParams, false, "", "Cphi"); + nParamsFound = vParams.size(); + if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) { + if (nParamsFound != 1) { + throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_CphiMX_ij[counter] = vParams[0]; + m_CphiMX_ij_coeff(0,counter) = m_CphiMX_ij[counter]; + } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) { + if (nParamsFound != 2) { + throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_CphiMX_ij_coeff(0,counter) = vParams[0]; + m_CphiMX_ij_coeff(1,counter) = vParams[1]; + m_CphiMX_ij[counter] = vParams[0]; + } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) { + if (nParamsFound != 5) { + throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + for (i = 0; i < nParamsFound; i++) { + m_CphiMX_ij_coeff(i, counter) = vParams[i]; + } + m_CphiMX_ij[counter] = vParams[0]; + } + } + + if (nodeName == "alpha1") { + stemp = xmlChild.value(); + m_Alpha1MX_ij[counter] = atofCheck(stemp.c_str()); + } + } + } + + /** + * Process an XML node called "ThetaAnion". + * This node contains all of the parameters necessary to describe + * the binary interactions between two anions. + */ + void HMWSoln::readXMLThetaAnion(XML_Node &BinSalt) { + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + string iName = BinSalt.attrib("anion1"); + if (iName == "") { + throw CanteraError("HMWSoln::readXMLThetaAnion", "no anion1 attrib"); + } + string jName = BinSalt.attrib("anion2"); + if (jName == "") { + throw CanteraError("HMWSoln::readXMLThetaAnion", "no anion2 attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species + */ + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + if (charge[iSpecies] >= 0) { + throw CanteraError("HMWSoln::readXMLThetaAnion", "anion1 charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + if (charge[jSpecies] >= 0) { + throw CanteraError("HMWSoln::readXMLThetaAnion", "anion2 charge problem"); + } + + int n = iSpecies * m_kk + jSpecies; + int counter = m_CounterIJ[n]; + int num = BinSalt.nChildren(); + for (int i = 0; i < num; i++) { + XML_Node &xmlChild = BinSalt.child(i); + stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + if (nodeName == "theta") { + stemp = xmlChild.value(); + double old = m_Theta_ij[counter]; + m_Theta_ij[counter] = atofCheck(stemp.c_str()); + if (old != 0.0) { + if (old != m_Theta_ij[counter]) { + throw CanteraError("HMWSoln::readXMLThetaAnion", "conflicting values"); + } + } + } + } + } + + /** + * Process an XML node called "ThetaCation". + * This node contains all of the parameters necessary to describe + * the binary interactions between two cation. + */ + void HMWSoln::readXMLThetaCation(XML_Node &BinSalt) { + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + string iName = BinSalt.attrib("cation1"); + if (iName == "") { + throw CanteraError("HMWSoln::readXMLThetaCation", "no cation1 attrib"); + } + string jName = BinSalt.attrib("cation2"); + if (jName == "") { + throw CanteraError("HMWSoln::readXMLThetaCation", "no cation2 attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species + */ + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + if (charge[iSpecies] <= 0) { + throw CanteraError("HMWSoln::readXMLThetaCation", "cation1 charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + if (charge[jSpecies] <= 0) { + throw CanteraError("HMWSoln::readXMLThetaCation", "cation2 charge problem"); + } + + int n = iSpecies * m_kk + jSpecies; + int counter = m_CounterIJ[n]; + int num = BinSalt.nChildren(); + for (int i = 0; i < num; i++) { + XML_Node &xmlChild = BinSalt.child(i); + stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + if (nodeName == "theta") { + stemp = xmlChild.value(); + double old = m_Theta_ij[counter]; + m_Theta_ij[counter] = atofCheck(stemp.c_str()); + if (old != 0.0) { + if (old != m_Theta_ij[counter]) { + throw CanteraError("HMWSoln::readXMLThetaCation", "conflicting values"); + } + } + } + } + } + + /** + * Process an XML node called "readXMLPsiCommonCation". + * This node contains all of the parameters necessary to describe + * the binary interactions between two anions and one common cation. + */ + void HMWSoln::readXMLPsiCommonCation(XML_Node &BinSalt) { + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + string kName = BinSalt.attrib("cation"); + if (kName == "") { + throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no cation attrib"); + } + string iName = BinSalt.attrib("anion1"); + if (iName == "") { + throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no anion1 attrib"); + } + string jName = BinSalt.attrib("anion2"); + if (jName == "") { + throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no anion2 attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species + */ + int kSpecies = speciesIndex(kName); + if (kSpecies < 0) { + return; + } + if (charge[kSpecies] <= 0) { + throw CanteraError("HMWSoln::readXMLPsiCommonCation", + "cation charge problem"); + } + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + if (charge[iSpecies] >= 0) { + throw CanteraError("HMWSoln::readXMLPsiCommonCation", + "anion1 charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + if (charge[jSpecies] >= 0) { + throw CanteraError("HMWSoln::readXMLPsiCommonCation", + "anion2 charge problem"); + } + + int n = iSpecies * m_kk + jSpecies; + int counter = m_CounterIJ[n]; + int num = BinSalt.nChildren(); + for (int i = 0; i < num; i++) { + XML_Node &xmlChild = BinSalt.child(i); + stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + if (nodeName == "theta") { + stemp = xmlChild.value(); + double old = m_Theta_ij[counter]; + m_Theta_ij[counter] = atofCheck(stemp.c_str()); + if (old != 0.0) { + if (old != m_Theta_ij[counter]) { + throw CanteraError("HMWSoln::readXMLPsiCommonCation", + "conflicting values"); + } + } + } + if (nodeName == "psi") { + stemp = xmlChild.value(); + double param = atofCheck(stemp.c_str()); + n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ; + m_Psi_ijk[n] = param; + n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies ; + m_Psi_ijk[n] = param; + n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies ; + m_Psi_ijk[n] = param; + n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies ; + m_Psi_ijk[n] = param; + n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies ; + m_Psi_ijk[n] = param; + n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies ; + m_Psi_ijk[n] = param; + } + } + } + + + + /** + * Process an XML node called "PsiCommonAnion". + * This node contains all of the parameters necessary to describe + * the binary interactions between two cations and one common anion. + */ + void HMWSoln::readXMLPsiCommonAnion(XML_Node &BinSalt) { + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + string kName = BinSalt.attrib("anion"); + if (kName == "") { + throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no anion attrib"); + } + string iName = BinSalt.attrib("cation1"); + if (iName == "") { + throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no cation1 attrib"); + } + string jName = BinSalt.attrib("cation2"); + if (jName == "") { + throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no cation2 attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species + */ + int kSpecies = speciesIndex(kName); + if (kSpecies < 0) { + return; + } + if (charge[kSpecies] >= 0) { + throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "anion charge problem"); + } + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + if (charge[iSpecies] <= 0) { + throw CanteraError("HMWSoln::readXMLPsiCommonAnion", + "cation1 charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + if (charge[jSpecies] <= 0) { + throw CanteraError("HMWSoln::readXMLPsiCommonAnion", + "cation2 charge problem"); + } + + int n = iSpecies * m_kk + jSpecies; + int counter = m_CounterIJ[n]; + int num = BinSalt.nChildren(); + for (int i = 0; i < num; i++) { + XML_Node &xmlChild = BinSalt.child(i); + stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + if (nodeName == "theta") { + stemp = xmlChild.value(); + double old = m_Theta_ij[counter]; + m_Theta_ij[counter] = atofCheck(stemp.c_str()); + if (old != 0.0) { + if (old != m_Theta_ij[counter]) { + throw CanteraError("HMWSoln::readXMLPsiCommonAnion", + "conflicting values"); + } + } + } + if (nodeName == "psi") { + stemp = xmlChild.value(); + double param = atofCheck(stemp.c_str()); + n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ; + m_Psi_ijk[n] = param; + n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies ; + m_Psi_ijk[n] = param; + n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies ; + m_Psi_ijk[n] = param; + n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies ; + m_Psi_ijk[n] = param; + n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies ; + m_Psi_ijk[n] = param; + n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies ; + m_Psi_ijk[n] = param; + } + } + } + + /** + * Process an XML node called "LambdaNeutral". + * This node contains all of the parameters necessary to describe + * the binary interactions between one neutral species and + * any other species (neutral or otherwise) in the mechanism. + */ + void HMWSoln::readXMLLambdaNeutral(XML_Node &BinSalt) { + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + string iName = BinSalt.attrib("neutral"); + if (iName == "") { + throw CanteraError("HMWSoln::readXMLLambdaNeutral", "no neutral attrib"); + } + string jName = BinSalt.attrib("speciesj"); + if (jName == "") { + throw CanteraError("HMWSoln::readXMLLambdaNeutral", "no speciesj attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species + */ + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + if (charge[iSpecies] != 0) { + throw CanteraError("HMWSoln::readXMLLambdaNeutral", "neutral charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + + int num = BinSalt.nChildren(); + for (int i = 0; i < num; i++) { + XML_Node &xmlChild = BinSalt.child(i); + stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + if (nodeName == "lambda") { + stemp = xmlChild.value(); + double old = m_Lambda_ij(iSpecies,jSpecies); + m_Lambda_ij(iSpecies,jSpecies) = atofCheck(stemp.c_str()); + if (old != 0.0) { + if (old != m_Lambda_ij(iSpecies,jSpecies)) { + throw CanteraError("HMWSoln::readXMLLambdaNeutral", "conflicting values"); + } + } + } + } + } + + /** + * Initialization routine for a HMWSoln phase. + * + * This is a virtual routine. This routine will call initThermo() + * for the parent class as well. + */ + void HMWSoln::initThermo() { + MolalityVPSSTP::initThermo(); + initLengths(); + } + + /** + * Import, construct, and initialize a HMWSoln phase + * specification from an XML tree into the current object. + * + * This routine is a precursor to constructPhaseXML(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void HMWSoln::constructPhaseFile(string inputFile, string id) { + + if (inputFile.size() == 0) { + throw CanteraError("HMWSoln:constructPhaseFile", + "input file is null"); + } + string path = findInputFile(inputFile); + ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("HMWSoln:constructPhaseFile","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + XML_Node &phaseNode_XML = xml(); + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("HMWSoln:constructPhaseFile", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + fxml_phase->copy(&phaseNode_XML); + constructPhaseXML(*fxml_phase, id); + delete fxml; + } + + /** + * Import, construct, and initialize a HMWSoln phase + * specification from an XML tree into the current object. + * + * Most of the work is carried out by the cantera base + * routine, importPhase(). That routine imports all of the + * species and element data, including the standard states + * of the species. + * + * Then, In this routine, we read the information + * particular to the specification of the activity + * coefficient model for the Pitzer parameterization. + * + * We also read information about the molar volumes of the + * standard states if present in the XML file. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void HMWSoln::constructPhaseXML(XML_Node& phaseNode, string id) { + string stemp; + if (id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("HMWSoln::constructPhaseXML", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("HMWSoln::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; + stemp = scNode.attrib("model"); + string formString = lowercase(stemp); + 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("HMWSoln::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("HMWSoln::constructPhaseXML", + "badly formed solvent XML node"); + } + solventName = nameSolventa[0]; + } + + /* + * Determine the form of the Pitzer model, + * We will use this information to size arrays below. + */ + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& scNode = thermoNode.child("activityCoefficients"); + m_formPitzer = m_formPitzer; + stemp = scNode.attrib("model"); + string formString = lowercase(stemp); + if (formString != "") { + if (formString == "pitzer" || formString == "default") { + m_formPitzer = PITZERFORM_BASE; + } else if (formString == "base") { + m_formPitzer = PITZERFORM_BASE; + } else { + throw CanteraError("HMWSoln::constructPhaseXML", + "Unknown Pitzer ActivityCoeff model: " + + formString); + } + } + + /* + * Determine the form of the temperature dependence + * of the Pitzer activity coefficient model. + */ + stemp = scNode.attrib("TempModel"); + formString = lowercase(stemp); + if (formString != "") { + if (formString == "constant" || formString == "default") { + m_formPitzerTemp = PITZER_TEMP_CONSTANT; + } else if (formString == "linear") { + m_formPitzerTemp = PITZER_TEMP_LINEAR; + } else if (formString == "complex" || formString == "complex1") { + m_formPitzerTemp = PITZER_TEMP_COMPLEX1; + } else { + throw CanteraError("HMWSoln::constructPhaseXML", + "Unknown Pitzer ActivityCoeff Temp model: " + + formString); + } + } + + /* + * Determine the reference temperature + * of the Pitzer activity coefficient model's temperature + * dependence formulation: defaults to 25C + */ + stemp = scNode.attrib("TempReference"); + formString = lowercase(stemp); + if (formString != "") { + m_TempPitzerRef = atofCheck(formString.c_str()); + } else { + m_TempPitzerRef = 273.15 + 25; + } + + } + + /* + * 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("HMWSoln::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 HMWSoln:: + 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"); + + /* + * 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("HMWSoln::initThermoXML", + "badly formed solvent XML node"); + } + solventName = nameSolventa[0]; + } + + /* + * 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. + */ + for (k = 0; k < m_kk; k++) { + string sname = speciesName(k); + if (solventName == sname) { + setSolvent(k); + if (k != 0) { + throw CanteraError("HMWSoln::initThermoXML", + "Solvent must be species 0 atm"); + } + m_indexSolvent = k; + break; + } + } + if (m_indexSolvent == -1) { + cout << "HMWSoln::initThermo: Solvent Name not found" + << endl; + throw CanteraError("HMWSoln::initThermoXML", + "Solvent name not found"); + } + if (m_indexSolvent != 0) { + throw CanteraError("HMWSoln::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 + } + + /* + * Initialize the water standard state model + */ + if (m_waterSS) delete m_waterSS; + m_waterSS = new WaterPDSS(this, 0); + + /* + * Initialize the water property calculator. It will share + * the internal eos water calculator. + */ + m_waterProps = new WaterProps(m_waterSS); + + /* + * 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")) { + XML_Node &ADebye = acNode.child("A_Debye"); + m_form_A_Debye = A_DEBYE_CONST; + stemp = "model"; + if (ADebye.hasAttrib(stemp)) { + string atemp = ADebye.attrib(stemp); + stemp = lowercase(atemp); + if (stemp == "water") { + m_form_A_Debye = A_DEBYE_WATER; + } + } + if (m_form_A_Debye == A_DEBYE_CONST) { + 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 the Maximum Ionic Strength + */ + if (acNode.hasChild("maxIonicStrength")) { + m_maxIionicStrength = getFloat(acNode, "maxIonicStrength"); +#ifdef DEBUG_HKM_NOT + cout << "m_maxIionicStrength = " + < 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( + // "HMWSoln::initThermo error", + // "no species match was found" + // ); + } else { + double val = fpValue(_b->second); + m_speciesCharge_Stoich[kk] = val; + } + } + } + } + + /* + * Loop through the children getting multiple instances of + * parameters + */ + if (acNodePtr) { + int n = acNodePtr->nChildren(); + for (int i = 0; i < n; i++) { + XML_Node &xmlACChild = acNodePtr->child(i); + stemp = xmlACChild.name(); + string nodeName = lowercase(stemp); + /* + * Process a binary salt field, or any of the other XML fields + * that make up the Pitzer Database. Entries will be ignored + * if any of the species in the entry isn't in the solution. + */ + if (nodeName == "binarysaltparameters") { + readXMLBinarySalt(xmlACChild); + } else if (nodeName == "thetaanion") { + readXMLThetaAnion(xmlACChild); + } else if (nodeName == "thetacation") { + readXMLThetaCation(xmlACChild); + } else if (nodeName == "psicommonanion") { + readXMLPsiCommonAnion(xmlACChild); + } else if (nodeName == "psicommoncation") { + readXMLPsiCommonCation(xmlACChild); + } else if (nodeName == "lambdaneutral") { + readXMLLambdaNeutral(xmlACChild); + } + } + } + + + } + + /* + * 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("HMWSoln::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("HMWSoln::initThermoXML", + "Bad electrolyte type: " + est); + } + } + } + } + } + + + /* + * Lastly set the state + */ + if (phaseNode.hasChild("state")) { + XML_Node& stateNode = phaseNode.child("state"); + setStateFromXML(stateNode); + } + + } +} diff --git a/Cantera/src/thermo/IdealMolalSoln.cpp b/Cantera/src/thermo/IdealMolalSoln.cpp index b7f7316bf..58f31d694 100644 --- a/Cantera/src/thermo/IdealMolalSoln.cpp +++ b/Cantera/src/thermo/IdealMolalSoln.cpp @@ -779,7 +779,7 @@ namespace Cantera { */ /** - * Initialization routine for an IdealMolalSoln phase: + * Initialization routine for an IdealMolalSoln phase. * * This is a virtual routine. This routine will call initThermo() * for the parent class as well. diff --git a/Cantera/src/thermo/IdealMolalSoln.h b/Cantera/src/thermo/IdealMolalSoln.h index 16621544d..44f5ae2eb 100644 --- a/Cantera/src/thermo/IdealMolalSoln.h +++ b/Cantera/src/thermo/IdealMolalSoln.h @@ -8,7 +8,11 @@ * solution assumes that all molality-based activity * coefficients are equal to one. */ - +/* + * 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. + */ /* * $Author$ * $Date$ diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index 05ed803e2..261fccedc 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -26,13 +26,14 @@ ELECTRO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o \ MolalityVPSSTP.o VPStandardStateTP.o \ IdealSolidSolnPhase.o IdealMolalSoln.o \ WaterPropsIAPWSphi.o WaterPropsIAPWS.o WaterProps.o \ - PDSS.o WaterPDSS.o WaterTP.o + PDSS.o WaterPDSS.o WaterTP.o \ + HMWSoln.o HMWSoln_input.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 + PDSS.h WaterPDSS.h WaterTP.h HMWSoln.h electrolytes.h endif ifeq ($(do_issp),1) ISSP_OBJ = IdealSolidSolnPhase.o diff --git a/Cantera/src/thermo/MolalityVPSSTP.h b/Cantera/src/thermo/MolalityVPSSTP.h index afd86b1dd..9a8b08e89 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.h +++ b/Cantera/src/thermo/MolalityVPSSTP.h @@ -8,14 +8,12 @@ * calculating liquid electrolyte thermodynamics. */ /* - * Copywrite (2005) Sandia Corporation. Under the terms of + * 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. */ /* - * $Author$ - * $Date$ - * $Revision$ + * $Id$ */ #ifndef CT_MOLALITYVPSSTP_H