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

[Thermo] Solvent is always the first species

Browse source 
No useful capabilities are provided by allowing the solvent species to vary, and
there are many places where the solvent was already implicitly assumed to be the
first species.
This commit is contained in:
Ray Speth 2017-08-02 22:39:40 -04:00
parent e3afaf5e61
commit e33fe6904d
6 changed files with 160 additions and 356 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

9
include/cantera/thermo/MolalityVPSSTP.h
View file

@ -225,10 +225,14 @@ public:
* molality.
*
* @param k the solvent index number
* @deprecated The solvent is always the first species in the phase. To be
* removed after Cantera 2.4.
*/
void setSolvent(size_t k);
//! Returns the solvent index.
//! @deprecated The solvent is always the first species in the phase. To be
//! removed after Cantera 2.4.
size_t solventIndex() const;
/**
@ -563,11 +567,6 @@ private:
virtual size_t findCLMIndex() const;
protected:
//! Index of the solvent. Currently the index of the solvent is hard-coded
//! to the value 0
size_t m_indexSolvent;
//! Scaling to be used for output of single-ion species activity
//! coefficients.
/*!

252
src/thermo/DebyeHuckel.cpp
View file

@ -146,7 +146,7 @@ void DebyeHuckel::getActivityConcentrations(doublereal* c) const
doublereal DebyeHuckel::standardConcentration(size_t k) const
{
double mvSolvent = m_speciesSize[m_indexSolvent];
double mvSolvent = m_speciesSize[0];
return 1.0 / mvSolvent;
}
@ -157,14 +157,11 @@ void DebyeHuckel::getActivities(doublereal* ac) const
// Update the molality array, m_molalities(). This requires an update due to
// mole fractions
s_update_lnMolalityActCoeff();
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
}
for (size_t k = 1; k < m_kk; k++) {
ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
}
double xmolSolvent = moleFraction(m_indexSolvent);
ac[m_indexSolvent] =
exp(m_lnActCoeffMolal[m_indexSolvent]) * xmolSolvent;
double xmolSolvent = moleFraction(0);
ac[0] = exp(m_lnActCoeffMolal[0]) * xmolSolvent;
}
void DebyeHuckel::getMolalityActivityCoefficients(doublereal* acMolality) const
@ -191,16 +188,13 @@ void DebyeHuckel::getChemPotentials(doublereal* mu) const
// Update the activity coefficients. This also updates the internal molality
// array.
s_update_lnMolalityActCoeff();
double xmolSolvent = moleFraction(m_indexSolvent);
for (size_t k = 0; k < m_kk; k++) {
if (m_indexSolvent != k) {
xx = std::max(m_molalities[k], SmallNumber);
mu[k] += RT() * (log(xx) + m_lnActCoeffMolal[k]);
}
double xmolSolvent = moleFraction(0);
for (size_t k = 1; k < m_kk; k++) {
xx = std::max(m_molalities[k], SmallNumber);
mu[k] += RT() * (log(xx) + m_lnActCoeffMolal[k]);
}
xx = std::max(xmolSolvent, SmallNumber);
mu[m_indexSolvent] +=
RT() * (log(xx) + m_lnActCoeffMolal[m_indexSolvent]);
mu[0] += RT() * (log(xx) + m_lnActCoeffMolal[0]);
}
void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const
@ -246,15 +240,13 @@ void DebyeHuckel::getPartialMolarEntropies(doublereal* sbar) const
// First we will add in the obvious dependence on the T term out front of
// the log activity term
doublereal mm;
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
mm = std::max(SmallNumber, m_molalities[k]);
sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal[k]);
}
for (size_t k = 1; k < m_kk; k++) {
mm = std::max(SmallNumber, m_molalities[k]);
sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal[k]);
}
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
mm = std::max(SmallNumber, xmolSolvent);
sbar[m_indexSolvent] -= GasConstant *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
sbar[0] -= GasConstant *(log(mm) + m_lnActCoeffMolal[0]);
// Check to see whether activity coefficients are temperature dependent. If
// they are, then calculate the their temperature derivatives and add them
@ -428,40 +420,6 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
setDebyeHuckelModel("Dilute_limit");
}
// Reconcile the solvent name and index.
// Get the Name of the Solvent:
// <solvent> solventName </solvent>
std::string solventName = "";
if (thermoNode.hasChild("solvent")) {
XML_Node& scNode = thermoNode.child("solvent");
vector<std::string> nameSolventa;
getStringArray(scNode, nameSolventa);
if (nameSolventa.size() != 1) {
throw CanteraError("DebyeHuckel::initThermoXML",
"badly formed solvent XML node");
}
solventName = nameSolventa[0];
}
for (size_t k = 0; k < m_kk; k++) {
std::string sname = speciesName(k);
if (solventName == sname) {
m_indexSolvent = k;
break;
}
}
if (m_indexSolvent == npos) {
cout << "DebyeHuckel::initThermoXML: Solvent Name not found"
<< endl;
throw CanteraError("DebyeHuckel::initThermoXML",
"Solvent name not found");
}
if (m_indexSolvent != 0) {
throw CanteraError("DebyeHuckel::initThermoXML",
"Solvent " + solventName +
" should be first species");
}
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
XML_Node* acNodePtr = 0;
@ -822,10 +780,8 @@ double DebyeHuckel::_lnactivityWaterHelgesonFixedForm() const
calcMolalities();
double oc = _osmoticCoeffHelgesonFixedForm();
double sum = 0.0;
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
sum += std::max(m_molalities[k], 0.0);
}
for (size_t k = 1; k < m_kk; k++) {
sum += std::max(m_molalities[k], 0.0);
}
if (sum > 2.0 * m_maxIionicStrength) {
sum = 2.0 * m_maxIionicStrength;
@ -876,7 +832,7 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
m_A_Debye = A_Debye_TP();
// Calculate a safe value for the mole fraction of the solvent
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
int est;
@ -920,7 +876,7 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
tmp = 0.0;
if (denomTmp > 0.0) {
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent || m_Aionic[k] != 0.0) {
if (k != 0 || m_Aionic[k] != 0.0) {
y = denomTmp * m_Aionic[k];
yp1 = y + 1.0;
sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
@ -931,9 +887,9 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
}
lnActivitySolvent += coeff * tmp;
tmp = 0.0;
for (size_t k = 0; k < m_kk; k++) {
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
if ((k != m_indexSolvent) && (z_k != 0.0)) {
if (z_k != 0.0) {
tmp += m_B_Dot[k] * m_molalities[k];
}
}
@ -967,9 +923,9 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
2.0 /3.0 * m_A_Debye * m_Mnaught *
m_IionicMolality * sqrt(m_IionicMolality) * sigma;
tmp = 0.0;
for (size_t k = 0; k < m_kk; k++) {
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
if ((k != m_indexSolvent) && (z_k != 0.0)) {
if (z_k != 0.0) {
tmp += m_B_Dot[k] * m_molalities[k];
}
}
@ -983,15 +939,13 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
lnActivitySolvent =
(xmolSolvent - 1.0)/xmolSolvent;
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_lnActCoeffMolal[k] =
- z_k * z_k * numTmp / (1.0 + denomTmp);
for (size_t j = 0; j < m_kk; j++) {
double beta = m_Beta_ij.value(k, j);
m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * beta;
}
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
m_lnActCoeffMolal[k] =
- z_k * z_k * numTmp / (1.0 + denomTmp);
for (size_t j = 0; j < m_kk; j++) {
double beta = m_Beta_ij.value(k, j);
m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * beta;
}
}
if (denomTmp > 0.0) {
@ -1020,18 +974,16 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
denomTmp *= m_Aionic[0];
numTmp = m_A_Debye * sqrt(m_IionicMolality);
tmpLn = log(1.0 + denomTmp);
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_lnActCoeffMolal[k] =
- z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp);
m_lnActCoeffMolal[k] +=
- 2.0 * z_k * z_k * m_A_Debye * tmpLn /
(3.0 * m_B_Debye * m_Aionic[0]);
for (size_t j = 0; j < m_kk; j++) {
m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] *
m_Beta_ij.value(k, j);
}
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
m_lnActCoeffMolal[k] =
- z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp);
m_lnActCoeffMolal[k] +=
- 2.0 * z_k * z_k * m_A_Debye * tmpLn /
(3.0 * m_B_Debye * m_Aionic[0]);
for (size_t j = 0; j < m_kk; j++) {
m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] *
m_Beta_ij.value(k, j);
}
}
sigma = 1.0 / (1.0 + denomTmp);
@ -1056,9 +1008,8 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
// Above, we calculated the ln(activitySolvent). Translate that into the
// molar-based activity coefficient by dividing by the solvent mole
// fraction. Solvents are not on the molality scale.
xmolSolvent = moleFraction(m_indexSolvent);
m_lnActCoeffMolal[m_indexSolvent] =
lnActivitySolvent - log(xmolSolvent);
xmolSolvent = moleFraction(0);
m_lnActCoeffMolal[0] = lnActivitySolvent - log(xmolSolvent);
}
void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
@ -1074,7 +1025,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
}
// Calculate a safe value for the mole fraction of the solvent
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);
double numdAdTTmp = dAdT * sqrtI;
@ -1089,7 +1040,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
}
d_lnActivitySolvent_dT = 2.0 / 3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrt(m_IionicMolality);
m_dlnActCoeffMolaldT[m_indexSolvent] = d_lnActivitySolvent_dT;
m_dlnActCoeffMolaldT[0] = d_lnActivitySolvent_dT;
break;
case DHFORM_BDOT_AK:
@ -1099,7 +1050,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
- z_k * z_k * numdAdTTmp / (1.0 + denomTmp * m_Aionic[k]);
}
m_dlnActCoeffMolaldT[m_indexSolvent] = 0.0;
m_dlnActCoeffMolaldT[0] = 0.0;
coeff = 2.0 / 3.0 * dAdT * m_Mnaught * sqrtI;
tmp = 0.0;
if (denomTmp > 0.0) {
@ -1111,7 +1062,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
}
}
m_dlnActCoeffMolaldT[m_indexSolvent] += coeff * tmp;
m_dlnActCoeffMolaldT[0] += coeff * tmp;
break;
case DHFORM_BDOT_ACOMMON:
@ -1128,19 +1079,15 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
} else {
sigma = 0.0;
}
m_dlnActCoeffMolaldT[m_indexSolvent] =
2.0 /3.0 * dAdT * m_Mnaught *
m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_BETAIJ:
denomTmp *= m_Aionic[0];
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldT[k] =
- z_k * z_k * numdAdTTmp / (1.0 + denomTmp);
}
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldT[k] = -z_k*z_k * numdAdTTmp / (1.0 + denomTmp);
}
if (denomTmp > 0.0) {
y = denomTmp;
@ -1149,28 +1096,23 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
} else {
sigma = 0.0;
}
m_dlnActCoeffMolaldT[m_indexSolvent] =
2.0 /3.0 * dAdT * m_Mnaught *
m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_PITZER_BETAIJ:
denomTmp *= m_Aionic[0];
tmpLn = log(1.0 + denomTmp);
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldT[k] =
- z_k * z_k * numdAdTTmp / (1.0 + denomTmp)
- 2.0 * z_k * z_k * dAdT * tmpLn
/ (m_B_Debye * m_Aionic[0]);
m_dlnActCoeffMolaldT[k] /= 3.0;
}
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldT[k] =
- z_k * z_k * numdAdTTmp / (1.0 + denomTmp)
- 2.0 * z_k * z_k * dAdT * tmpLn / (m_B_Debye * m_Aionic[0]);
m_dlnActCoeffMolaldT[k] /= 3.0;
}
sigma = 1.0 / (1.0 + denomTmp);
m_dlnActCoeffMolaldT[m_indexSolvent] =
2.0 /3.0 * dAdT * m_Mnaught *
m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
@ -1193,7 +1135,7 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
}
// Calculate a safe value for the mole fraction of the solvent
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);
double numd2AdT2Tmp = d2AdT2 * sqrtI;
@ -1214,7 +1156,7 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
- z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp * m_Aionic[k]);
}
m_d2lnActCoeffMolaldT2[m_indexSolvent] = 0.0;
m_d2lnActCoeffMolaldT2[0] = 0.0;
coeff = 2.0 / 3.0 * d2AdT2 * m_Mnaught * sqrtI;
tmp = 0.0;
if (denomTmp > 0.0) {
@ -1226,7 +1168,7 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
}
}
m_d2lnActCoeffMolaldT2[m_indexSolvent] += coeff * tmp;
m_d2lnActCoeffMolaldT2[0] += coeff * tmp;
break;
case DHFORM_BDOT_ACOMMON:
@ -1243,19 +1185,15 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
} else {
sigma = 0.0;
}
m_d2lnActCoeffMolaldT2[m_indexSolvent] =
2.0 /3.0 * d2AdT2 * m_Mnaught *
m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_BETAIJ:
denomTmp *= m_Aionic[0];
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_d2lnActCoeffMolaldT2[k] =
- z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp);
}
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
m_d2lnActCoeffMolaldT2[k] = -z_k*z_k * numd2AdT2Tmp / (1.0 + denomTmp);
}
if (denomTmp > 0.0) {
y = denomTmp;
@ -1264,28 +1202,23 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
} else {
sigma = 0.0;
}
m_d2lnActCoeffMolaldT2[m_indexSolvent] =
2.0 /3.0 * d2AdT2 * m_Mnaught *
m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_PITZER_BETAIJ:
denomTmp *= m_Aionic[0];
tmpLn = log(1.0 + denomTmp);
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_d2lnActCoeffMolaldT2[k] =
- z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp)
- 2.0 * z_k * z_k * d2AdT2 * tmpLn
/ (m_B_Debye * m_Aionic[0]);
m_d2lnActCoeffMolaldT2[k] /= 3.0;
}
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
m_d2lnActCoeffMolaldT2[k] =
- z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp)
- 2.0 * z_k * z_k * d2AdT2 * tmpLn / (m_B_Debye * m_Aionic[0]);
m_d2lnActCoeffMolaldT2[k] /= 3.0;
}
sigma = 1.0 / (1.0 + denomTmp);
m_d2lnActCoeffMolaldT2[m_indexSolvent] =
2.0 /3.0 * d2AdT2 * m_Mnaught *
m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
@ -1308,7 +1241,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
}
// Calculate a safe value for the mole fraction of the solvent
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);
double numdAdPTmp = dAdP * sqrtI;
@ -1334,7 +1267,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
}
}
m_dlnActCoeffMolaldP[m_indexSolvent] = 0.0;
m_dlnActCoeffMolaldP[0] = 0.0;
coeff = 2.0 / 3.0 * dAdP * m_Mnaught * sqrtI;
tmp = 0.0;
if (denomTmp > 0.0) {
@ -1346,7 +1279,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
}
}
m_dlnActCoeffMolaldP[m_indexSolvent] += coeff * tmp;
m_dlnActCoeffMolaldP[0] += coeff * tmp;
break;
case DHFORM_BDOT_ACOMMON:
@ -1363,19 +1296,16 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
} else {
sigma = 0.0;
}
m_dlnActCoeffMolaldP[m_indexSolvent] =
m_dlnActCoeffMolaldP[0] =
2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_BETAIJ:
denomTmp *= m_Aionic[0];
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldP[k] =
- z_k * z_k * numdAdPTmp / (1.0 + denomTmp);
}
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldP[k] = - z_k*z_k * numdAdPTmp / (1.0 + denomTmp);
}
if (denomTmp > 0.0) {
y = denomTmp;
@ -1384,28 +1314,24 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
} else {
sigma = 0.0;
}
m_dlnActCoeffMolaldP[m_indexSolvent] =
2.0 /3.0 * dAdP * m_Mnaught *
m_dlnActCoeffMolaldP[0] = 2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_PITZER_BETAIJ:
denomTmp *= m_Aionic[0];
tmpLn = log(1.0 + denomTmp);
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldP[k] =
- z_k * z_k * numdAdPTmp / (1.0 + denomTmp)
- 2.0 * z_k * z_k * dAdP * tmpLn
/ (m_B_Debye * m_Aionic[0]);
m_dlnActCoeffMolaldP[k] /= 3.0;
}
for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldP[k] =
- z_k * z_k * numdAdPTmp / (1.0 + denomTmp)
- 2.0 * z_k * z_k * dAdP * tmpLn
/ (m_B_Debye * m_Aionic[0]);
m_dlnActCoeffMolaldP[k] /= 3.0;
}
sigma = 1.0 / (1.0 + denomTmp);
m_dlnActCoeffMolaldP[m_indexSolvent] =
2.0 /3.0 * dAdP * m_Mnaught *
m_dlnActCoeffMolaldP[0] = 2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;

82
src/thermo/HMWSoln.cpp
View file

@ -307,7 +307,7 @@ void HMWSoln::getActivityConcentrations(doublereal* c) const
doublereal HMWSoln::standardConcentration(size_t k) const
{
getStandardVolumes(m_tmpV.data());
double mvSolvent = m_tmpV[m_indexSolvent];
double mvSolvent = m_tmpV[0];
if (k > 0) {
return m_Mnaught / mvSolvent;
}
@ -323,14 +323,11 @@ void HMWSoln::getActivities(doublereal* ac) const
s_update_lnMolalityActCoeff();
// Now calculate the array of activities.
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
}
for (size_t k = 1; k < m_kk; k++) {
ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
}
double xmolSolvent = moleFraction(m_indexSolvent);
ac[m_indexSolvent] =
exp(m_lnActCoeffMolal_Scaled[m_indexSolvent]) * xmolSolvent;
double xmolSolvent = moleFraction(0);
ac[0] = exp(m_lnActCoeffMolal_Scaled[0]) * xmolSolvent;
}
void HMWSoln::getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const
@ -357,16 +354,13 @@ void HMWSoln::getChemPotentials(doublereal* mu) const
// Update the activity coefficients. This also updates the internal molality
// array.
s_update_lnMolalityActCoeff();
double xmolSolvent = moleFraction(m_indexSolvent);
for (size_t k = 0; k < m_kk; k++) {
if (m_indexSolvent != k) {
xx = std::max(m_molalities[k], SmallNumber);
mu[k] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[k]);
}
double xmolSolvent = moleFraction(0);
for (size_t k = 1; k < m_kk; k++) {
xx = std::max(m_molalities[k], SmallNumber);
mu[k] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[k]);
}
xx = std::max(xmolSolvent, SmallNumber);
mu[m_indexSolvent] +=
RT() * (log(xx) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
mu[0] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[0]);
}
void HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const
@ -406,15 +400,13 @@ void HMWSoln::getPartialMolarEntropies(doublereal* sbar) const
// First we will add in the obvious dependence on the T term out front of
// the log activity term
doublereal mm;
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
mm = std::max(SmallNumber, m_molalities[k]);
sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
}
for (size_t k = 1; k < m_kk; k++) {
mm = std::max(SmallNumber, m_molalities[k]);
sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
}
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
mm = std::max(SmallNumber, xmolSolvent);
sbar[m_indexSolvent] -= GasConstant *(log(mm) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
sbar[0] -= GasConstant *(log(mm) + m_lnActCoeffMolal_Scaled[0]);
// Check to see whether activity coefficients are temperature dependent. If
// they are, then calculate the their temperature derivatives and add them
@ -781,7 +773,7 @@ void HMWSoln::s_update_lnMolalityActCoeff() const
// Now do the main calculation.
s_updatePitzer_lnMolalityActCoeff();
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
double lnxs = log(xx);
@ -929,7 +921,7 @@ void HMWSoln::calcMolalitiesCropped() const
if (cropMethod == 1) {
double* molF = m_gamma_tmp.data();
getMoleFractions(molF);
double xmolSolvent = molF[m_indexSolvent];
double xmolSolvent = molF[0];
if (xmolSolvent >= MC_X_o_cutoff_) {
return;
}
@ -1231,12 +1223,6 @@ void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const
void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
{
// HKM -> Assumption is made that the solvent is species 0.
if (m_indexSolvent != 0) {
throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
"Wrong index solvent value!");
}
// Use the CROPPED molality of the species in solution.
const vector_fp& molality = m_molalitiesCropped;
@ -1913,7 +1899,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
//
// We have just computed act_0. However, this routine returns
// ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
if (m_debugCalc) {
@ -1959,12 +1945,6 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
// immediately preceding the calling of this routine. Therefore, some
// quantities do not need to be recalculated in this routine.
// HKM -> Assumption is made that the solvent is species 0.
if (m_indexSolvent != 0) {
throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
"Wrong index solvent value!");
}
const vector_fp& molality = m_molalitiesCropped;
double* d_gamma_dT_Unscaled = m_gamma_tmp.data();
@ -2569,12 +2549,6 @@ void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const
void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
{
// HKM -> Assumption is made that the solvent is species 0.
if (m_indexSolvent != 0) {
throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
"Wrong index solvent value!");
}
const double* molality = m_molalitiesCropped.data();
// Local variables defined by Coltrin
@ -3171,12 +3145,6 @@ void HMWSoln::s_update_dlnMolalityActCoeff_dP() const
void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
{
// HKM -> Assumption is made that the solvent is species 0.
if (m_indexSolvent != 0) {
throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
"Wrong index solvent value!");
}
const double* molality = m_molalitiesCropped.data();
// Local variables defined by Coltrin
@ -3830,27 +3798,27 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
// Calculate the molalities. Currently, the molalities may not be current
// with respect to the contents of the State objects' data.
calcMolalities();
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
if (IMS_typeCutoff_ == 0) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (IMS_typeCutoff_ == 1) {
if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
double tmp = log(xx * IMS_gamma_k_min_);
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
IMS_lnActCoeffMolal_[0] = log(IMS_gamma_o_min_);
return;
} else {
// If we are in the middle region, calculate the connecting polynomials
@ -3887,7 +3855,7 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
IMS_lnActCoeffMolal_[0] = lngammao;
}
} else if (IMS_typeCutoff_ == 2) {
// Exponentials - trial 2
@ -3895,7 +3863,7 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else {
double xoverc = xmolSolvent/IMS_cCut_;
@ -3921,7 +3889,7 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
IMS_lnActCoeffMolal_[0] = lngammao;
}
}
return;

41
src/thermo/HMWSoln_input.cpp
View file

@ -993,49 +993,10 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
// Get the Name of the Solvent:
// <solvent> solventName </solvent>
string solventName = "";
if (thermoNode.hasChild("solvent")) {
XML_Node& scNode = thermoNode.child("solvent");
vector<string> nameSolventa;
getStringArray(scNode, nameSolventa);
if (nameSolventa.size() != 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 (size_t 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 == npos) {
std::cout << "HMWSoln::initThermo: Solvent Name not found"
<< std::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 specification of the standard states for species in the
// solution. This includes the molar volumes data blocks for incompressible
// species.
@ -1239,7 +1200,7 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral;
}
}
m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent;
m_electrolyteSpeciesType[0] = cEST_solvent;
// First look at the species database. Look for the subelement
// "stoichIsMods" in each of the species SS databases.

100
src/thermo/IdealMolalSoln.cpp
View file

@ -180,10 +180,10 @@ doublereal IdealMolalSoln::standardConcentration(size_t k) const
case 0:
break;
case 1:
return c0 = 1.0 /m_speciesMolarVolume[m_indexSolvent];
return c0 = 1.0 /m_speciesMolarVolume[0];
break;
case 2:
c0 = 1.0 / m_speciesMolarVolume[m_indexSolvent];
c0 = 1.0 / m_speciesMolarVolume[0];
break;
}
return c0;
@ -200,12 +200,11 @@ void IdealMolalSoln::getActivities(doublereal* ac) const
for (size_t k = 0; k < m_kk; k++) {
ac[k] = m_molalities[k];
}
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
// Limit the activity coefficient to be finite as the solvent mole
// fraction goes to zero.
xmolSolvent = std::max(m_xmolSolventMIN, xmolSolvent);
ac[m_indexSolvent] =
exp((xmolSolvent - 1.0)/xmolSolvent);
ac[0] = exp((xmolSolvent - 1.0)/xmolSolvent);
} else {
s_updateIMS_lnMolalityActCoeff();
@ -214,9 +213,8 @@ void IdealMolalSoln::getActivities(doublereal* ac) const
for (size_t k = 1; k < m_kk; k++) {
ac[k] = m_molalities[k] * exp(IMS_lnActCoeffMolal_[k]);
}
double xmolSolvent = moleFraction(m_indexSolvent);
ac[m_indexSolvent] =
exp(IMS_lnActCoeffMolal_[m_indexSolvent]) * xmolSolvent;
double xmolSolvent = moleFraction(0);
ac[0] = exp(IMS_lnActCoeffMolal_[0]) * xmolSolvent;
}
}
@ -226,12 +224,11 @@ void IdealMolalSoln::getMolalityActivityCoefficients(doublereal* acMolality) con
for (size_t k = 0; k < m_kk; k++) {
acMolality[k] = 1.0;
}
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
// Limit the activity coefficient to be finite as the solvent mole
// fraction goes to zero.
xmolSolvent = std::max(m_xmolSolventMIN, xmolSolvent);
acMolality[m_indexSolvent] =
exp((xmolSolvent - 1.0)/xmolSolvent) / xmolSolvent;
acMolality[0] = exp((xmolSolvent - 1.0)/xmolSolvent) / xmolSolvent;
} else {
s_updateIMS_lnMolalityActCoeff();
std::copy(IMS_lnActCoeffMolal_.begin(), IMS_lnActCoeffMolal_.end(), acMolality);
@ -245,9 +242,6 @@ void IdealMolalSoln::getMolalityActivityCoefficients(doublereal* acMolality) con
void IdealMolalSoln::getChemPotentials(doublereal* mu) const
{
// Assertion is made for speed
AssertThrow(m_indexSolvent == 0, "solvent not the first species");
// First get the standard chemical potentials. This requires updates of
// standard state as a function of T and P These are defined at unit
// molality.
@ -258,7 +252,7 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
calcMolalities();
// get the solvent mole fraction
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
if (IMS_typeCutoff_ == 0 || xmolSolvent > 3.* IMS_X_o_cutoff_/2.0) {
for (size_t k = 1; k < m_kk; k++) {
@ -269,8 +263,7 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
// Do the solvent
// -> see my notes
double xx = std::max(xmolSolvent, SmallNumber);
mu[m_indexSolvent] +=
(RT() * (xmolSolvent - 1.0) / xx);
mu[0] += (RT() * (xmolSolvent - 1.0) / xx);
} else {
// Update the activity coefficients. This also updates the internal
// molality array.
@ -281,8 +274,7 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
mu[k] += RT() * (log(xx) + IMS_lnActCoeffMolal_[k]);
}
double xx = std::max(xmolSolvent, SmallNumber);
mu[m_indexSolvent] +=
RT() * (log(xx) + IMS_lnActCoeffMolal_[m_indexSolvent]);
mu[0] += RT() * (log(xx) + IMS_lnActCoeffMolal_[0]);
}
}
@ -299,14 +291,12 @@ void IdealMolalSoln::getPartialMolarEntropies(doublereal* sbar) const
getEntropy_R(sbar);
calcMolalities();
if (IMS_typeCutoff_ == 0) {
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
doublereal mm = std::max(SmallNumber, m_molalities[k]);
sbar[k] -= GasConstant * log(mm);
}
for (size_t k = 1; k < m_kk; k++) {
doublereal mm = std::max(SmallNumber, m_molalities[k]);
sbar[k] -= GasConstant * log(mm);
}
double xmolSolvent = moleFraction(m_indexSolvent);
sbar[m_indexSolvent] -= (GasConstant * (xmolSolvent - 1.0) / xmolSolvent);
double xmolSolvent = moleFraction(0);
sbar[0] -= (GasConstant * (xmolSolvent - 1.0) / xmolSolvent);
} else {
// Update the activity coefficients, This also update the internally
// stored molalities.
@ -315,15 +305,13 @@ void IdealMolalSoln::getPartialMolarEntropies(doublereal* sbar) const
// First we will add in the obvious dependence on the T term out front
// of the log activity term
doublereal mm;
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
mm = std::max(SmallNumber, m_molalities[k]);
sbar[k] -= GasConstant * (log(mm) + IMS_lnActCoeffMolal_[k]);
}
for (size_t k = 1; k < m_kk; k++) {
mm = std::max(SmallNumber, m_molalities[k]);
sbar[k] -= GasConstant * (log(mm) + IMS_lnActCoeffMolal_[k]);
}
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
mm = std::max(SmallNumber, xmolSolvent);
sbar[m_indexSolvent] -= GasConstant *(log(mm) + IMS_lnActCoeffMolal_[m_indexSolvent]);
sbar[0] -= GasConstant *(log(mm) + IMS_lnActCoeffMolal_[0]);
}
}
@ -377,19 +365,6 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
setStandardConcentrationModel(scNode["model"]);
}
// Get the Name of the Solvent:
// <solvent> solventName </solvent>
std::string solventName = "";
if (thermoNode.hasChild("solvent")) {
std::vector<std::string> nameSolventa;
getStringArray(thermoNode.child("solvent"), nameSolventa);
if (nameSolventa.size() != 1) {
throw CanteraError("IdealMolalSoln::initThermoXML",
"badly formed solvent XML node");
}
solventName = nameSolventa[0];
}
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
std::string modelString = acNode.attrib("model");
@ -425,25 +400,6 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
setCutoffModel("none");
}
}
// Reconcile the solvent name and index.
for (size_t k = 0; k < m_kk; k++) {
if (solventName == speciesName(k)) {
m_indexSolvent = k;
break;
}
}
if (m_indexSolvent == npos) {
std::cout << "IdealMolalSoln::initThermo: Solvent Name not found"
<< std::endl;
throw CanteraError("IdealMolalSoln::initThermo",
"Solvent name not found");
}
if (m_indexSolvent != 0) {
throw CanteraError("IdealMolalSoln::initThermo",
"Solvent " + solventName +
" should be first species");
}
}
void IdealMolalSoln::initThermo()
@ -494,28 +450,28 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
// with respect to the contents of the State objects' data.
calcMolalities();
double xmolSolvent = moleFraction(m_indexSolvent);
double xmolSolvent = moleFraction(0);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
if (IMS_typeCutoff_ == 0) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (IMS_typeCutoff_ == 1) {
if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
double tmp = log(xx * IMS_gamma_k_min_);
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
IMS_lnActCoeffMolal_[0] = log(IMS_gamma_o_min_);
return;
} else {
// If we are in the middle region, calculate the connecting polynomials
@ -552,7 +508,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
IMS_lnActCoeffMolal_[0] = lngammao;
}
} else if (IMS_typeCutoff_ == 2) {
// Exponentials - trial 2
@ -560,7 +516,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else {
double xoverc = xmolSolvent/IMS_cCut_;
@ -584,7 +540,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
IMS_lnActCoeffMolal_[0] = lngammao;
}
}
}

32
src/thermo/MolalityVPSSTP.cpp
View file

@ -22,7 +22,6 @@ namespace Cantera
{
MolalityVPSSTP::MolalityVPSSTP() :
m_indexSolvent(0),
m_pHScalingType(PHSCALE_PITZER),
m_indexCLM(npos),
m_weightSolvent(18.01528),
@ -53,20 +52,15 @@ int MolalityVPSSTP::pHScale() const
void MolalityVPSSTP::setSolvent(size_t k)
{
if (k >= m_kk) {
throw CanteraError("MolalityVPSSTP::setSolute ",
"bad value");
}
m_indexSolvent = k;
AssertThrowMsg(m_indexSolvent==0, "MolalityVPSSTP::setSolvent",
"Molality-based methods limit solvent id to being 0");
m_weightSolvent = molecularWeight(k);
m_Mnaught = m_weightSolvent / 1000.;
warn_deprecated("MolalityVPSSTP::setSolvent", "Solvent is always the first"
" species. To be removed after Cantera 2.4.");
}
size_t MolalityVPSSTP::solventIndex() const
{
return m_indexSolvent;
warn_deprecated("MolalityVPSSTP::solventIndex", "Solvent is always the"
" first species. To be removed after Cantera 2.4.");
return 0;
}
void MolalityVPSSTP::setMoleFSolventMin(doublereal xmolSolventMIN)
@ -87,7 +81,7 @@ doublereal MolalityVPSSTP::moleFSolventMin() const
void MolalityVPSSTP::calcMolalities() const
{
getMoleFractions(m_molalities.data());
double xmolSolvent = std::max(m_molalities[m_indexSolvent], m_xmolSolventMIN);
double xmolSolvent = std::max(m_molalities[0], m_xmolSolventMIN);
double denomInv = 1.0/ (m_Mnaught * xmolSolvent);
for (size_t k = 0; k < m_kk; k++) {
m_molalities[k] *= denomInv;
@ -110,8 +104,8 @@ void MolalityVPSSTP::setMolalities(const doublereal* const molal)
Lsum += molal[k];
}
double tmp = 1.0 / Lsum;
m_molalities[m_indexSolvent] = tmp / m_Mnaught;
double sum = m_molalities[m_indexSolvent];
m_molalities[0] = tmp / m_Mnaught;
double sum = m_molalities[0];
for (size_t k = 1; k < m_kk; k++) {
m_molalities[k] = tmp * molal[k];
sum += m_molalities[k];
@ -137,7 +131,7 @@ void MolalityVPSSTP::setMolalitiesByName(const compositionMap& mMap)
// Get a vector of mole fractions
vector_fp mf(m_kk, 0.0);
getMoleFractions(mf.data());
double xmolSmin = std::max(mf[m_indexSolvent], m_xmolSolventMIN);
double xmolSmin = std::max(mf[0], m_xmolSolventMIN);
for (size_t k = 0; k < m_kk; k++) {
double mol_k = getValue(mMap, speciesName(k), 0.0);
if (mol_k > 0) {
@ -228,8 +222,7 @@ void MolalityVPSSTP::getActivities(doublereal* ac) const
void MolalityVPSSTP::getActivityCoefficients(doublereal* ac) const
{
getMolalityActivityCoefficients(ac);
AssertThrow(m_indexSolvent==0, "MolalityVPSSTP::getActivityCoefficients");
double xmolSolvent = std::max(moleFraction(m_indexSolvent), m_xmolSolventMIN);
double xmolSolvent = std::max(moleFraction(0), m_xmolSolventMIN);
for (size_t k = 1; k < m_kk; k++) {
ac[k] /= xmolSolvent;
}
@ -254,7 +247,7 @@ doublereal MolalityVPSSTP::osmoticCoefficient() const
}
double oc = 1.0;
if (sum > 1.0E-200) {
oc = - log(act[m_indexSolvent]) / (m_Mnaught * sum);
oc = - log(act[0]) / (m_Mnaught * sum);
}
return oc;
}
@ -371,7 +364,8 @@ bool MolalityVPSSTP::addSpecies(shared_ptr<Species> spec)
if (added) {
if (m_kk == 1) {
// The solvent defaults to species 0
setSolvent(0);
m_weightSolvent = molecularWeight(0);
m_Mnaught = m_weightSolvent / 1000.;
}
m_molalities.push_back(0.0);
}