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

Add RedlichKisterVPSSTP::s_update_dlnActCoeff_dlnX_diag()

Browse source 
Fixes bug in RedlichKisterVPSSTP::getlnActCoeffdlnX_diag() which would always return 0

Resolves #293
This commit is contained in:
pwcnorthrop 2015-10-07 14:53:37 -05:00 • committed by Ray Speth
parent db769a47c7
commit 125d2e0e23
2 changed files with 61 additions and 1 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
include/cantera/thermo/RedlichKisterVPSSTP.h
View file

@ -432,6 +432,16 @@ private:
*/
void s_update_dlnActCoeff_dX_() const;
//! Internal routine that calculates the total derivative of the activity
//! coefficients with respect to the log of the mole fractions.
/*!
* This function will be called to update the internally stored vector of
* the total derivatives (i.e. not assuming other mole fractions are
* constant) of the natural logarithm of the activity coefficients with
* respect to the log of the mole fraction.
*/
void s_update_dlnActCoeff_dlnX_diag() const;
public:
//! Utility routine that calculates a literature expression
/*!

52
src/thermo/RedlichKisterVPSSTP.cpp
View file

@ -383,6 +383,56 @@ void RedlichKisterVPSSTP::getd2lnActCoeffdT2(doublereal* d2lnActCoeffdT2) const
}
}
void RedlichKisterVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
{
double T = temperature();
dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
double deltaX = XA - XB;
size_t N = m_N_ij[i];
double poly = 1.0;
double sum = 0.0;
vector_fp& he_vec = m_HE_m_ij[i];
vector_fp& se_vec = m_SE_m_ij[i];
double sumMm1 = 0.0;
double polyMm1 = 1.0;
double polyMm2 = 1.0;
double sumMm2 = 0.0;
for (size_t m = 0; m < N; m++) {
double A_ge = (he_vec[m] - T * se_vec[m]) / (GasConstant * T);;
sum += A_ge * poly;
poly *= deltaX;
if (m >= 1) {
sumMm1 += (A_ge * polyMm1 * m);
polyMm1 *= deltaX;
}
if (m >= 2) {
sumMm2 += (A_ge * polyMm2 * m * (m - 1.0));
polyMm2 *= deltaX;
}
}
for (size_t k = 0; k < m_kk; k++) {
if (iA == k) {
dlnActCoeffdlnX_diag_[k] +=
XA * (- (1-XA+XB) * sum + 2*(1.0 - XA) * XB * sumMm1
+ sumMm1 * (XB * (1 - 2*XA + XB) - XA * (1 - XA + 2*XB))
+ 2 * XA * XB * sumMm2 * (1.0 - XA + XB));
} else if (iB == k) {
dlnActCoeffdlnX_diag_[k] +=
XB * (- (1-XB+XA) * sum - 2*(1.0 - XB) * XA * sumMm1
+ sumMm1 * (XA * (2*XB - XA - 1) - XB * (-2*XA + XB - 1))
- 2 * XA * XB * sumMm2 * (-XA - 1 + XB));
}
}
}
}
void RedlichKisterVPSSTP::s_update_dlnActCoeff_dX_() const
{
doublereal T = temperature();
@ -472,7 +522,7 @@ void RedlichKisterVPSSTP::getdlnActCoeffdlnN_diag(doublereal* dlnActCoeffdlnN_di
void RedlichKisterVPSSTP::getdlnActCoeffdlnX_diag(doublereal* dlnActCoeffdlnX_diag) const
{
s_update_dlnActCoeff_dX_();
s_update_dlnActCoeff_dlnX_diag();
for (size_t k = 0; k < m_kk; k++) {
dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
}