Implement and test MaskellSolidSolnPhase::getActivityCoefficients().

include/cantera/thermo/MaskellSolidSolnPhase.h

@@ -349,6 +349,7 @@ protected:
 
     //! Index of the species whose mole fraction defines the extent of reduction r
     int product_species_index;
+    int reactant_species_index;
 
 private:
     // Functions to calculate some of the pieces of the mixing terms.

src/thermo/MaskellSolidSolnPhase.cpp

@@ -29,7 +29,8 @@ MaskellSolidSolnPhase::MaskellSolidSolnPhase() :
     m_g0_RT(2),
     m_s0_R(2),
     h_mixing(0.0),
-    product_species_index(0)
+    product_species_index(0),
+    reactant_species_index(1)
 {
 }
 //=====================================================================================================
@@ -42,7 +43,8 @@ MaskellSolidSolnPhase::MaskellSolidSolnPhase(const MaskellSolidSolnPhase& b) :
     m_g0_RT(2),
     m_s0_R(2),
     h_mixing(0.0),
-    product_species_index(0)
+    product_species_index(0),
+    reactant_species_index(1)
 {
     *this = b;
 }
@@ -163,8 +165,17 @@ void MaskellSolidSolnPhase::setMolarDensity(const doublereal n)
 void MaskellSolidSolnPhase::
 getActivityCoefficients(doublereal* ac) const
 {
- // throw CanteraError("MaskellSolidSolnPhase::getActivityCoefficients()",
+    _updateThermo();
- //                    "needs to be implemented");
+    const doublereal r = moleFraction(product_species_index);
+    const doublereal pval = p(r);
+    const doublereal fmval = fm(r);
+    const doublereal rfm = r * fmval;
+    const doublereal RT = GasConstant * temperature();
+    const doublereal A = (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval))  /
+                         (std::pow(1 - r - rfm, 1 + pval) * (1 - r));
+    const doublereal B = pval * h_mixing / RT;
+    ac[product_species_index] = 1 / (A * r * r) * std::exp(-B);
+    ac[reactant_species_index] = A * r / (1 - r) * std::exp(B);
 }
 
 void MaskellSolidSolnPhase::
@@ -176,12 +187,12 @@ getChemPotentials(doublereal* mu) const
     const doublereal fmval = fm(r);
     const doublereal rfm = r * fmval;
     const doublereal RT = GasConstant * temperature();
-    const doublereal muDelta = -pval * h_mixing - GasConstant * temperature()
+    const doublereal DgbarDr = pval * h_mixing +
-                       * std::log( (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval) * r)  /
+                               GasConstant * temperature() *
-                                   (std::pow(1 - r - rfm, 1 + pval) * (1 - r)) );
+                               std::log( (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval) * r)  /
-    const int sign = (product_species_index == 0) ? 1 : -1;
+                               (std::pow(1 - r - rfm, 1 + pval) * (1 - r)) );
-    mu[0] = RT * m_g0_RT[0] + sign * muDelta;
+    mu[product_species_index] = RT * m_g0_RT[product_species_index] - DgbarDr;
-    mu[1] = RT * m_g0_RT[1] - sign * muDelta;
+    mu[reactant_species_index] = RT * m_g0_RT[reactant_species_index] + DgbarDr;
 }
 
 void MaskellSolidSolnPhase::
@@ -280,12 +291,17 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string
             XML_Node& scNode = thNode.child("product_species");
             std::string product_species_name = scNode.value();
             product_species_index = speciesIndex(product_species_name);
-            if( product_species_index == npos )
+            if( product_species_index == static_cast<int>(npos) )
             {
               throw CanteraError(subname.c_str(),
                                  "Species " + product_species_name + " not found.");
             }
-            std::cout << "parsed product_species_index = " << product_species_index << std::endl;
+            if( product_species_index == 0 )
+            {
+              reactant_species_index = 1;
+            } else {
+              reactant_species_index = 0;
+            }
         }
 
     } else {

test/thermo/MaskellSolidSolnPhase_Test.cpp

@@ -63,6 +63,7 @@ TEST_F(MaskellSolidSolnPhase_Test, chem_potentials)
     set_r(0.5);
 
     MaskellSolidSolnPhase * maskell_phase = dynamic_cast<MaskellSolidSolnPhase *>(test_phase);
+    ASSERT_TRUE(maskell_phase != NULL);
 
     maskell_phase->set_h_mix(0.);
     const double expected_result_0[9] = {1.2338461168724738e7, 8.011774549216799e6, 4.990989640314685e6, 2.415973128783114e6, 0., -2.415973128783114e6, -4.99098964031469e6, -8.0117745492168e6, -1.2338461168724738e7};
@@ -77,6 +78,33 @@ TEST_F(MaskellSolidSolnPhase_Test, chem_potentials)
     check_chemPotentials(expected_result_minus_5000);
 }
 
+TEST_F(MaskellSolidSolnPhase_Test, activityCoeffs)
+{
+    const std::string valid_file("../data/MaskellSolidSolnPhase_valid.xml");
+    initializeTestPhaseWithXML(valid_file);
+    test_phase->setState_TP(298., 1.);
+    set_r(0.5);
+
+    MaskellSolidSolnPhase * maskell_phase = dynamic_cast<MaskellSolidSolnPhase *>(test_phase);
+    ASSERT_TRUE(maskell_phase != NULL);
+
+    // Test that mu0 + RT log(activityCoeff * MoleFrac) == mu
+    const double RT = GasConstant * 298.;
+    std::vector<double> mu0(2);
+    std::vector<double> activityCoeffs(2);
+    std::vector<double> chemPotentials(2);
+    for(int i=0; i < 9; ++i)
+    {
+        const double r = 0.1 * (i+1);
+        set_r(r);
+        test_phase->getChemPotentials(&chemPotentials[0]);
+        test_phase->getActivityCoefficients(&activityCoeffs[0]);
+        test_phase->getStandardChemPotentials(&mu0[0]);
+        EXPECT_NEAR(chemPotentials[0], mu0[0] + RT*std::log(activityCoeffs[0] * r), 1.e-6);
+        EXPECT_NEAR(chemPotentials[1], mu0[1] + RT*std::log(activityCoeffs[1] * (1-r)), 1.e-6);
+    }
+}
+
 TEST_F(MaskellSolidSolnPhase_Test, partialMolarVolumes)
 {
     const std::string valid_file("../data/MaskellSolidSolnPhase_valid.xml");