From 1565bad198afd5363b7fe06cc38f40790622892d Mon Sep 17 00:00:00 2001 From: Victor Brunini Date: Fri, 7 Mar 2014 17:14:13 +0000 Subject: [PATCH] MaskellSolidPhaseSoln thermo updates. Fix sign error in chemical potentials and update activity coefficients accordingly. Correct implementation of getActivityConcentrations() Additional unit tests. --- .../cantera/thermo/MaskellSolidSolnPhase.h | 27 ----------- src/thermo/MaskellSolidSolnPhase.cpp | 16 +++---- test/thermo/MaskellSolidSolnPhase_Test.cpp | 48 ++++++++++++++++--- 3 files changed, 49 insertions(+), 42 deletions(-) diff --git a/include/cantera/thermo/MaskellSolidSolnPhase.h b/include/cantera/thermo/MaskellSolidSolnPhase.h index cdeecb6f9..ee0ddf7f9 100644 --- a/include/cantera/thermo/MaskellSolidSolnPhase.h +++ b/include/cantera/thermo/MaskellSolidSolnPhase.h @@ -57,33 +57,6 @@ public: * concentration divided by the standard concentration is also * equal to the activity of species. * - * For this implementation the activity is defined to be the - * mole fraction of the species. The generalized concentration - * is defined to be equal to the mole fraction divided by - * the partial molar volume. The generalized concentrations - * for species in this phase therefore have units of - * kmol m-3. Rate constants must reflect this fact. - * - * On a general note, the following must be true. - * For an ideal solution, the generalized concentration must consist - * of the mole fraction multiplied by a constant. The constant may be - * fairly arbitrarily chosen, with differences adsorbed into the - * reaction rate expression. 1/V_N, 1/V_k, or 1 are equally good, - * as long as the standard concentration is adjusted accordingly. - * However, it must be a constant (and not the concentration, btw, - * which is a function of the mole fractions) in order for the - * ideal solution properties to hold at the same time having the - * standard concentration to be independent of the mole fractions. - * - * In this implementation the form of the generalized concentrations - * depend upon the member attribute, #m_formGC. - * - * HKM Note: We have absorbed the pressure dependence of the pure species - * state into the thermodynamics functions. Therefore the - * standard state on which the activities are based depend - * on both temperature and pressure. If we hadn't, it would have - * appeared in this function in a very awkward exp[] format. - * * @param c Pointer to array of doubles of length m_kk, which on exit * will contain the generalized concentrations. */ diff --git a/src/thermo/MaskellSolidSolnPhase.cpp b/src/thermo/MaskellSolidSolnPhase.cpp index b471a0355..46fdc12a3 100644 --- a/src/thermo/MaskellSolidSolnPhase.cpp +++ b/src/thermo/MaskellSolidSolnPhase.cpp @@ -68,10 +68,10 @@ getActivityConcentrations(doublereal* c) const { std::vector pmv(m_kk); getPartialMolarVolumes(&pmv[0]); - const doublereal* const dtmp = moleFractdivMMW(); - const double mmw = meanMolecularWeight(); - for (size_t k = 0; k < m_kk; k++) { - c[k] = dtmp[k] * mmw / pmv[k]; + getActivityCoefficients(c); + for(unsigned sp=0; sp < m_kk; ++sp) + { + c[sp] *= moleFraction(sp) / pmv[sp]; } } @@ -174,8 +174,8 @@ getActivityCoefficients(doublereal* ac) const 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); + ac[product_species_index] = A * std::exp(B); + ac[reactant_species_index] = 1 / (A * r * (1-r) ) * std::exp(-B); } void MaskellSolidSolnPhase:: @@ -191,8 +191,8 @@ getChemPotentials(doublereal* mu) const 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; + 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:: diff --git a/test/thermo/MaskellSolidSolnPhase_Test.cpp b/test/thermo/MaskellSolidSolnPhase_Test.cpp index 018f94b9d..be2cd2684 100644 --- a/test/thermo/MaskellSolidSolnPhase_Test.cpp +++ b/test/thermo/MaskellSolidSolnPhase_Test.cpp @@ -37,8 +37,8 @@ public: const double r = 0.1 * (i+1); set_r(r); test_phase->getChemPotentials(&chemPotentials[0]); - EXPECT_NEAR(expected_result[i], chemPotentials[0], 1.e-6); - EXPECT_NEAR(1000.-expected_result[i], chemPotentials[1], 1.e-6); + EXPECT_NEAR(-expected_result[i], chemPotentials[0], 1.e-6); + EXPECT_NEAR(1000.+expected_result[i], chemPotentials[1], 1.e-6); } } }; @@ -78,6 +78,18 @@ TEST_F(MaskellSolidSolnPhase_Test, chem_potentials) check_chemPotentials(expected_result_minus_5000); } +TEST_F(MaskellSolidSolnPhase_Test, partialMolarVolumes) +{ + const std::string valid_file("../data/MaskellSolidSolnPhase_valid.xml"); + initializeTestPhaseWithXML(valid_file); + ASSERT_TRUE(dynamic_cast(test_phase) != NULL); + + std::vector pmv(2); + test_phase->getPartialMolarVolumes(&pmv[0]); + EXPECT_EQ(0.005, pmv[0]); + EXPECT_EQ(0.01, pmv[1]); +} + TEST_F(MaskellSolidSolnPhase_Test, activityCoeffs) { const std::string valid_file("../data/MaskellSolidSolnPhase_valid.xml"); @@ -105,16 +117,38 @@ TEST_F(MaskellSolidSolnPhase_Test, activityCoeffs) } } -TEST_F(MaskellSolidSolnPhase_Test, partialMolarVolumes) +TEST_F(MaskellSolidSolnPhase_Test, standardConcentrations) { const std::string valid_file("../data/MaskellSolidSolnPhase_valid.xml"); initializeTestPhaseWithXML(valid_file); ASSERT_TRUE(dynamic_cast(test_phase) != NULL); - std::vector pmv(2); - test_phase->getPartialMolarVolumes(&pmv[0]); - EXPECT_EQ(0.005, pmv[0]); - EXPECT_EQ(0.01, pmv[1]); + EXPECT_DOUBLE_EQ(1.0 / 0.005, test_phase->standardConcentration(0)); + EXPECT_DOUBLE_EQ(1.0 / 0.01, test_phase->standardConcentration(1)); +} + +TEST_F(MaskellSolidSolnPhase_Test, activityConcentrations) +{ + const std::string valid_file("../data/MaskellSolidSolnPhase_valid.xml"); + initializeTestPhaseWithXML(valid_file); + ASSERT_TRUE(dynamic_cast(test_phase) != NULL); + + // Check to make sure activityConcentration_i == standardConcentration_i * gamma_i * X_i + std::vector standardConcs(2); + std::vector activityCoeffs(2); + std::vector activityConcentrations(2); + for(int i=0; i < 9; ++i) + { + const double r = 0.1 * (i+1); + set_r(r); + test_phase->getActivityCoefficients(&activityCoeffs[0]); + standardConcs[0] = test_phase->standardConcentration(0); + standardConcs[1] = test_phase->standardConcentration(1); + test_phase->getActivityConcentrations(&activityConcentrations[0]); + + EXPECT_NEAR(standardConcs[0] * r * activityCoeffs[0], activityConcentrations[0], 1.e-6); + EXPECT_NEAR(standardConcs[1] * (1-r) * activityCoeffs[1], activityConcentrations[1], 1.e-6); + } } };