MaskellSolidPhaseSoln thermo updates.

Fix sign error in chemical potentials and update
activity coefficients accordingly.

Correct implementation of getActivityConcentrations()

Additional unit tests.
This commit is contained in:
Victor Brunini 2014-03-07 17:14:13 +00:00
parent 0981391ac4
commit 1565bad198
3 changed files with 49 additions and 42 deletions

View file

@ -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<SUP>-3</SUP>. 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.
*/

View file

@ -68,10 +68,10 @@ getActivityConcentrations(doublereal* c) const
{
std::vector<doublereal> 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::

View file

@ -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<MaskellSolidSolnPhase *>(test_phase) != NULL);
std::vector<double> 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<MaskellSolidSolnPhase *>(test_phase) != NULL);
std::vector<double> 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<MaskellSolidSolnPhase *>(test_phase) != NULL);
// Check to make sure activityConcentration_i == standardConcentration_i * gamma_i * X_i
std::vector<double> standardConcs(2);
std::vector<double> activityCoeffs(2);
std::vector<double> 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);
}
}
};