[Thermo] Fix enthalpy/entropy calculations in RedlichKisterVPSSTP

Derivatives of the activity coefficients were being calculated incorrectly,
leading to bad values for enthalpies, entropies, and specific heat
capacities.

Partially resolves #1320

include/cantera/thermo/RedlichKisterVPSSTP.h

@@ -51,7 +51,9 @@ namespace Cantera
  *    G^E_{i} =   n X_{Ai} X_{Bi} \sum_m \left( A^{i}_m {\left( X_{Ai} -  X_{Bi} \right)}^m \right)
  * \f]
  *
- * and where we can break down the Gibbs free energy contributions into enthalpy and entropy contributions
+ * where n is the total moles in the solution and where we can break down the Gibbs free
+ * energy contributions into enthalpy and entropy contributions by defining
+ * \f$ A^i_m = H^i_m - T S^i_m \f$ :
  *
  * \f[
  *    H^E_i = n X_{Ai} X_{Bi} \sum_m \left( H^{i}_m {\left( X_{Ai} -  X_{Bi} \right)}^m \right)
@@ -61,8 +63,8 @@ namespace Cantera
  *    S^E_i = n X_{Ai} X_{Bi} \sum_m \left( S^{i}_m {\left( X_{Ai} -  X_{Bi} \right)}^m \right)
  * \f]
  *
- * where n is the total moles in the solution. The activity of a species defined
- * in the phase is given by an excess Gibbs free energy formulation.
+ * The activity of a species defined in the phase is given by an excess Gibbs free
+ * energy formulation:
  *
  * \f[
  *      a_k = \gamma_k  X_k
@@ -80,6 +82,20 @@ namespace Cantera
  *                           + \sum_i \delta_{Ai,k} X_{Ai} X_{Bi} \sum_m \left(  A^{i}_0 +  A^{i}_m {\left( X_{Ai} -  X_{Bi} \right)}^{m-1} (1 - X_{Ai} + X_{Bi}) \right)
  * \f]
  *
+ * Evaluating thermodynamic properties requires the following derivatives of
+ * \f$ \ln(\gamma_k) \f$:
+ *
+ * \f[
+ *    \frac{d \ln( \gamma_k )}{dT} = - \frac{1}{RT^2} \left( \sum_i \delta_{Ai,k} (1 - X_{Ai}) X_{Bi} \sum_m \left( H^{i}_m {\left( X_{Ai} -  X_{Bi} \right)}^m \right)
+ *        + \sum_i \delta_{Ai,k} X_{Ai} X_{Bi} \sum_m \left(  H^{i}_0 +  H^{i}_m {\left( X_{Ai} -  X_{Bi} \right)}^{m-1} (1 - X_{Ai} + X_{Bi}) \right) \right)
+ * \f]
+ *
+ * and
+ *
+ * \f[
+ *    \frac{d^2 \ln( \gamma_k )}{dT^2} = -\frac{2}{T} \frac{d \ln( \gamma_k )}{dT}
+ * \f]
+ *
  * This object inherits from the class VPStandardStateTP. Therefore, the
  * specification and calculation of all standard state and reference state
  * values are handled at that level. Various functional forms for the standard

src/thermo/RedlichKisterVPSSTP.cpp

@@ -254,8 +254,8 @@ void RedlichKisterVPSSTP::s_update_lnActCoeff() const
 
 void RedlichKisterVPSSTP::s_update_dlnActCoeff_dT() const
 {
+    doublereal T = temperature();
     dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
-    d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0);
 
     for (size_t i = 0; i < numBinaryInteractions_; i++) {
         size_t iA = m_pSpecies_A_ij[i];
@@ -266,17 +266,17 @@ void RedlichKisterVPSSTP::s_update_dlnActCoeff_dT() const
         size_t N = m_N_ij[i];
         doublereal poly = 1.0;
         doublereal sum = 0.0;
-        vector_fp& se_vec = m_SE_m_ij[i];
+        vector_fp& he_vec = m_HE_m_ij[i];
         doublereal sumMm1 = 0.0;
         doublereal polyMm1 = 1.0;
         doublereal sum2 = 0.0;
         for (size_t m = 0; m < N; m++) {
-            doublereal A_ge = - se_vec[m];
-            sum += A_ge * poly;
-            sum2 += A_ge * (m + 1) * poly;
+            double h_e = - he_vec[m] / (GasConstant * T * T);
+            sum += h_e * poly;
+            sum2 += h_e * (m + 1) * poly;
             poly *= deltaX;
             if (m >= 1) {
-                sumMm1 += (A_ge * polyMm1 * m);
+                sumMm1 += (h_e * polyMm1 * m);
                 polyMm1 *= deltaX;
             }
         }
@@ -292,6 +292,10 @@ void RedlichKisterVPSSTP::s_update_dlnActCoeff_dT() const
             }
         }
     }
+
+    for (size_t k = 0; k < m_kk; k++) {
+        d2lnActCoeffdT2_Scaled_[k] = -2 / T * dlnActCoeffdT_Scaled_[k];
+    }
 }
 
 void RedlichKisterVPSSTP::getdlnActCoeffdT(doublereal* dlnActCoeffdT) const

test/data/consistency-cases.yaml

@@ -322,7 +322,7 @@ Redlich-Kister:
     file: thermo-models.yaml
     phase: Redlich-Kister-LiC6
     known-failures:
-      c[vp]_eq_.+: "Implementation of cv is incorrect. See GitHub Issue #1320"
+      cv_eq_.+: "Implementation of cv is incorrect. See GitHub Issue #1320"
   states:
   - {T: 300, P: 1 atm, X: {Li(C6): 0.85, V(C6): 0.15}}
   - {T: 440, P: 1 atm, X: {Li(C6): 0.75, V(C6): 0.25}}