diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp index 09d399d53..dd8374a8b 100644 --- a/Cantera/src/thermo/HMWSoln.cpp +++ b/Cantera/src/thermo/HMWSoln.cpp @@ -2762,6 +2762,27 @@ namespace Cantera { } } #endif + /* + * Zeta interaction term + */ + for (int k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + int izeta = j; + int jzeta = i; + n = izeta * m_kk * m_kk + jzeta * m_kk + k; + double zeta = psi_ijk[n]; + if (zeta != 0.0) { + sum5 = sum5 + molality[j]*molality[k]*zeta; +#ifdef DEBUG_MODE + if (m_debugCalc) { + snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Zeta term on %-16s m_n m_a zeta_nMa = %10.5f\n", snj.c_str(), + molality[j]*molality[k]*psi_ijk[n]); + } +#endif + } + } + } } } /* @@ -2913,6 +2934,28 @@ namespace Cantera { } } #endif + /* + * Zeta interaction term + */ + for (int k = 1; k < m_kk; k++) { + if (charge[k] > 0.0) { + int izeta = j; + int jzeta = k; + int kzeta = i; + n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta; + double zeta = psi_ijk[n]; + if (zeta != 0.0) { + sum5 = sum5 + molality[j]*molality[k]*zeta; +#ifdef DEBUG_MODE + if (m_debugCalc) { + snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Zeta term on %-16s m_n m_c zeta_ncX = %10.5f\n", snj.c_str(), + molality[j]*molality[k]*psi_ijk[n]); + } +#endif + } + } + } } } m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5; @@ -2930,18 +2973,64 @@ namespace Cantera { * ------ -> equations agree with my notes, * -> Equations agree with Pitzer, */ - if (charge[i] == 0.0 ) { + if (charge[i] == 0.0) { +#ifdef DEBUG_MODE + if (m_debugCalc) { + sni = speciesName(i); + printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str()); + } +#endif sum1 = 0.0; + sum3 = 0.0; for (j = 1; j < m_kk; j++) { sum1 = sum1 + molality[j]*2.0*m_Lambda_nj(i,j); +#ifdef DEBUG_MODE + if (m_debugCalc) { + if (m_Lambda_nj(i,j) != 0.0) { + snj = speciesName(j) + ":"; + printf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj.c_str(), + molality[j]*2.0*m_Lambda_nj(i,j)); + } + } +#endif + + /* + * Zeta term -> we piggyback on the psi term + */ + if (charge[j] > 0.0) { + for (k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + n = k + j * m_kk + i * m_kk * m_kk; + sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n]; +#ifdef DEBUG_MODE + if (m_debugCalc) { + if (psi_ijk[n] != 0.0) { + snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), + molality[j]*molality[k]*psi_ijk[n]); + } + } +#endif + } + } + } } sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]; - m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2; +#ifdef DEBUG_MODE + if (m_debugCalc) { + if (m_Mu_nnn[i] != 0.0) { + printf(" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n", + 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]); + } + } +#endif + + m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3; gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]); #ifdef DEBUG_MODE if (m_debugCalc) { sni = speciesName(i); - printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n", + printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n", sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]); } #endif @@ -3066,6 +3155,19 @@ namespace Cantera { sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj(j,k); } } + if (charge[k] < 0.0) { + int izeta = j; + for (m = 1; m < m_kk; m++) { + if (charge[m] > 0.0) { + int jzeta = m; + n = k + jzeta * m_kk + izeta * m_kk * m_kk; + double zeta = psi_ijk[n]; + if (zeta != 0.0) { + sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta; + } + } + } + } } sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn[j]; } @@ -3640,6 +3742,20 @@ namespace Cantera { if (charge[j] == 0) { sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i); } + /* + * Zeta interaction term + */ + for (int k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + int izeta = j; + int jzeta = i; + n = izeta * m_kk * m_kk + jzeta * m_kk + k; + double zeta_L = psi_ijk_L[n]; + if (zeta_L != 0.0) { + sum5 = sum5 + molality[j]*molality[k]*zeta_L; + } + } + } } /* * Add all of the contributions up to yield the log of the @@ -3725,6 +3841,18 @@ namespace Cantera { */ if (charge[j] == 0.0) { sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i); + for (int k = 1; k < m_kk; k++) { + if (charge[k] > 0.0) { + int izeta = j; + int jzeta = k; + int kzeta = i; + n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta; + double zeta_L = psi_ijk_L[n]; + if (zeta_L != 0.0) { + sum5 = sum5 + molality[j]*molality[k]*zeta_L; + } + } + } } } m_dlnActCoeffMolaldT_Unscaled[i] = @@ -3747,11 +3875,23 @@ namespace Cantera { */ if (charge[i] == 0.0 ) { sum1 = 0.0; + sum3 = 0.0; for (j = 1; j < m_kk; j++) { sum1 = sum1 + molality[j]*2.0*m_Lambda_nj_L(i,j); + /* + * Zeta term -> we piggyback on the psi term + */ + if (charge[j] > 0.0) { + for (k = 1; k < m_kk; k++) { + 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]; + } + } + } } sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_L[i]; - m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2; + m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2 + sum3; d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]); #ifdef DEBUG_MODE if (m_debugCalc) { @@ -3881,6 +4021,19 @@ namespace Cantera { sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_L(j,k); } } + if (charge[k] < 0.0) { + int izeta = j; + for (m = 1; m < m_kk; m++) { + if (charge[m] > 0.0) { + int jzeta = m; + n = k + jzeta * m_kk + izeta * m_kk * m_kk; + double zeta_L = psi_ijk_L[n]; + if (zeta_L != 0.0) { + sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L; + } + } + } + } } sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn_L[j]; } @@ -4456,6 +4609,20 @@ namespace Cantera { */ if (charge[j] == 0) { sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_LL(j,i); + /* + * Zeta interaction term + */ + for (int k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + int izeta = j; + int jzeta = i; + n = izeta * m_kk * m_kk + jzeta * m_kk + k; + double zeta_LL = psi_ijk_LL[n]; + if (zeta_LL != 0.0) { + sum5 = sum5 + molality[j]*molality[k]*zeta_LL; + } + } + } } } /* @@ -4542,6 +4709,21 @@ namespace Cantera { */ if (charge[j] == 0.0) { sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_LL(j,i); + /* + * Zeta interaction term + */ + for (int k = 1; k < m_kk; k++) { + if (charge[k] > 0.0) { + int izeta = j; + int jzeta = k; + int kzeta = i; + n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta; + double zeta_LL = psi_ijk_LL[n]; + if (zeta_LL != 0.0) { + sum5 = sum5 + molality[j]*molality[k]*zeta_LL; + } + } + } } } m_d2lnActCoeffMolaldT2_Unscaled[i] = @@ -4563,11 +4745,23 @@ namespace Cantera { */ if (charge[i] == 0.0 ) { sum1 = 0.0; + sum3 = 0.0; for (j = 1; j < m_kk; j++) { sum1 = sum1 + molality[j]*2.0*m_Lambda_nj_LL(i,j); + /* + * Zeta term -> we piggyback on the psi term + */ + if (charge[j] > 0.0) { + for (k = 1; k < m_kk; k++) { + 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]; + } + } + } } sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i]; - m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2; + m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3; #ifdef DEBUG_MODE if (m_debugCalc) { sni = speciesName(i); @@ -4697,6 +4891,19 @@ namespace Cantera { sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_LL(j,k); } } + if (charge[k] < 0.0) { + int izeta = j; + for (m = 1; m < m_kk; m++) { + if (charge[m] > 0.0) { + int jzeta = m; + n = k + jzeta * m_kk + izeta * m_kk * m_kk; + double zeta_LL = psi_ijk_LL[n]; + if (zeta_LL != 0.0) { + sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL; + } + } + } + } } sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_LL[j]; @@ -5269,7 +5476,21 @@ namespace Cantera { * for Anions, do the neutral species interaction */ if (charge[j] == 0) { - sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i); + sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_P(j,i); + /* + * Zeta interaction term + */ + for (int k = 1; k < m_kk; k++) { + if (charge[k] < 0.0) { + int izeta = j; + int jzeta = i; + n = izeta * m_kk * m_kk + jzeta * m_kk + k; + double zeta_P = psi_ijk_P[n]; + if (zeta_P != 0.0) { + sum5 = sum5 + molality[j]*molality[k]*zeta_P; + } + } + } } } @@ -5357,7 +5578,22 @@ namespace Cantera { * for Anions, do the neutral species interaction */ if (charge[j] == 0.0) { - sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i); + sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_P(j,i); + /* + * Zeta interaction term + */ + for (int k = 1; k < m_kk; k++) { + if (charge[k] > 0.0) { + int izeta = j; + int jzeta = k; + int kzeta = i; + n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta; + double zeta_P = psi_ijk_P[n]; + if (zeta_P != 0.0) { + sum5 = sum5 + molality[j]*molality[k]*zeta_P; + } + } + } } } m_dlnActCoeffMolaldP_Unscaled[i] = @@ -5378,11 +5614,23 @@ namespace Cantera { */ if (charge[i] == 0.0) { sum1 = 0.0; + sum3 = 0.0; for (j = 1; j < m_kk; j++) { sum1 += molality[j]*2.0*m_Lambda_nj_P(i,j); + /* + * Zeta term -> we piggyback on the psi term + */ + if (charge[j] > 0.0) { + for (k = 1; k < m_kk; k++) { + 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]; + } + } + } } sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i]; - m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2; + m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3; #ifdef DEBUG_MODE if (m_debugCalc) { sni = speciesName(i); @@ -5513,6 +5761,19 @@ namespace Cantera { sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_P(j,k); } } + if (charge[k] < 0.0) { + int izeta = j; + for (m = 1; m < m_kk; m++) { + if (charge[m] > 0.0) { + int jzeta = m; + n = k + jzeta * m_kk + izeta * m_kk * m_kk; + double zeta_P = psi_ijk_P[n]; + if (zeta_P != 0.0) { + sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P; + } + } + } + } } sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_P[j]; diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index 71fd47399..b2041e10f 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -3403,6 +3403,16 @@ namespace Cantera { */ void readXMLMunnnNeutral(XML_Node &BinSalt); + //! Process an XML node called "zetaCation" + /*! + * This node contains all of the parameters necessary to describe + * the ternary interactions between one neutral, one cation, and one anion. + * + * @param BinSalt reference to the XML_Node named psiCommonCation + * containing the + * neutral - cation - anion interaction + */ + void readXMLZetaCation(XML_Node &BinSalt); //! Precalculate the IMS Cutoff parameters for typeCutoff = 2 void calcIMSCutoffParams_(); diff --git a/Cantera/src/thermo/HMWSoln_input.cpp b/Cantera/src/thermo/HMWSoln_input.cpp index 2f96f6f04..80dc8aa22 100644 --- a/Cantera/src/thermo/HMWSoln_input.cpp +++ b/Cantera/src/thermo/HMWSoln_input.cpp @@ -731,6 +731,8 @@ namespace Cantera { } } } + + /** * Process an XML node called "LambdaNeutral". @@ -895,7 +897,111 @@ namespace Cantera { } } + /* + * Process an XML node called "readXMLZetaCation". + * This node contains all of the parameters necessary to describe + * the ternary interactions between a neutral, a cation and an anion + */ + void HMWSoln::readXMLZetaCation(XML_Node &BinSalt) { + string xname = BinSalt.name(); + if (xname != "zetaCation") { + throw CanteraError("HMWSoln::readXMLZetaCation", + "Incorrect name for processing this routine: " + xname); + } + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + vector_fp vParams; + int nParamsFound = 0; + + string iName = BinSalt.attrib("neutral"); + if (iName == "") { + throw CanteraError("HMWSoln::readXMLZetaCation", "no neutral attrib"); + } + + string jName = BinSalt.attrib("cation1"); + if (jName == "") { + throw CanteraError("HMWSoln::readXMLZetaCation", "no cation1 attrib"); + } + string kName = BinSalt.attrib("anion1"); + if (kName == "") { + throw CanteraError("HMWSoln::readXMLZetaCation", "no anion1 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.0) { + throw CanteraError("HMWSoln::readXMLZetaCation", "neutral charge problem"); + } + + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + if (charge[jSpecies] <= 0.0) { + throw CanteraError("HMWSoln::readXLZetaCation", "cation1 charge problem"); + } + + int kSpecies = speciesIndex(kName); + if (kSpecies < 0) { + return; + } + if (charge[kSpecies] >= 0.0) { + throw CanteraError("HMWSoln::readXMLZetaCation", "anion1 charge problem"); + } + + 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 == "zeta") { + getFloatArray(xmlChild, vParams, false, "", "zeta"); + nParamsFound = vParams.size(); + int n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ; + + if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) { + if (nParamsFound != 1) { + throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for " + + iName + "::" + jName + "::" + kName, + "wrong number of params found"); + } + m_Psi_ijk_coeff(0,n) = vParams[0]; + m_Psi_ijk[n] = vParams[0]; + } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) { + if (nParamsFound != 2) { + throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for " + + iName + "::" + jName + "::" + kName, + "wrong number of params found"); + } + m_Psi_ijk_coeff(0,n) = vParams[0]; + m_Psi_ijk_coeff(1,n) = vParams[1]; + m_Psi_ijk[n] = vParams[0]; + } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) { + if (nParamsFound == 1) { + vParams.resize(5, 0.0); + nParamsFound = 5; + } else if (nParamsFound != 5) { + throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for " + + iName + "::" + jName + "::" + kName, + "wrong number of params found"); + } + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + } + + // There are no duplicate entries + } + } + } + /* @@ -1435,6 +1541,8 @@ namespace Cantera { readXMLPsiCommonCation(xmlACChild); } else if (nodeName == "lambdaneutral") { readXMLLambdaNeutral(xmlACChild); + } else if (nodeName == "zetacation") { + readXMLZetaCation(xmlACChild); } } } diff --git a/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp b/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp index 7e2211a35..9de92b01b 100644 --- a/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp +++ b/Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp @@ -15,7 +15,6 @@ */ /* - * * $Date$ * $Revision$ */