research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Added the zeta interaction term.

Browse source
This commit is contained in:
Harry Moffat 2009-01-09 04:17:31 +00:00
parent 1bd9dc858a
commit 477d0080bc
4 changed files with 387 additions and 9 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

277
Cantera/src/thermo/HMWSoln.cpp
View file

@ -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];

10
Cantera/src/thermo/HMWSoln.h
View file

@ -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_();

108
Cantera/src/thermo/HMWSoln_input.cpp
View file

@ -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);
}
}
}

1
Cantera/src/thermo/VPSSMgr_Water_HKFT.cpp
View file

@ -15,7 +15,6 @@
*/
/*
*
* $Date$
* $Revision$
*/