From 0981391ac469c8b1ea996f50aaaec307cddfc6bb Mon Sep 17 00:00:00 2001 From: Victor Brunini Date: Fri, 7 Mar 2014 00:53:33 +0000 Subject: [PATCH] Implement and test MaskellSolidSolnPhase::getActivityCoefficients(). --- .../cantera/thermo/MaskellSolidSolnPhase.h | 1 + src/thermo/MaskellSolidSolnPhase.cpp | 40 +++++++++++++------ test/thermo/MaskellSolidSolnPhase_Test.cpp | 28 +++++++++++++ 3 files changed, 57 insertions(+), 12 deletions(-) diff --git a/include/cantera/thermo/MaskellSolidSolnPhase.h b/include/cantera/thermo/MaskellSolidSolnPhase.h index 405012748..cdeecb6f9 100644 --- a/include/cantera/thermo/MaskellSolidSolnPhase.h +++ b/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. diff --git a/src/thermo/MaskellSolidSolnPhase.cpp b/src/thermo/MaskellSolidSolnPhase.cpp index f72101d69..b471a0355 100644 --- a/src/thermo/MaskellSolidSolnPhase.cpp +++ b/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()", - // "needs to be implemented"); + _updateThermo(); + 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() - * std::log( (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval) * r) / - (std::pow(1 - r - rfm, 1 + pval) * (1 - r)) ); - const int sign = (product_species_index == 0) ? 1 : -1; - mu[0] = RT * m_g0_RT[0] + sign * muDelta; - mu[1] = RT * m_g0_RT[1] - sign * muDelta; + const doublereal DgbarDr = pval * h_mixing + + GasConstant * temperature() * + std::log( (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval) * r) / + (std::pow(1 - r - rfm, 1 + pval) * (1 - r)) ); + mu[product_species_index] = RT * m_g0_RT[product_species_index] - DgbarDr; + 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(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 { diff --git a/test/thermo/MaskellSolidSolnPhase_Test.cpp b/test/thermo/MaskellSolidSolnPhase_Test.cpp index 3b0daa9e5..018f94b9d 100644 --- a/test/thermo/MaskellSolidSolnPhase_Test.cpp +++ b/test/thermo/MaskellSolidSolnPhase_Test.cpp @@ -63,6 +63,7 @@ TEST_F(MaskellSolidSolnPhase_Test, chem_potentials) set_r(0.5); MaskellSolidSolnPhase * maskell_phase = dynamic_cast(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(test_phase); + ASSERT_TRUE(maskell_phase != NULL); + + // Test that mu0 + RT log(activityCoeff * MoleFrac) == mu + const double RT = GasConstant * 298.; + std::vector mu0(2); + std::vector activityCoeffs(2); + std::vector 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");