[Kinetics] Validate balance of surface sites for interface reactions

The number of surface sites should be the same for the reactants and products.

Fixes #412
This commit is contained in:
Ray Speth 2017-01-20 17:21:42 -05:00
parent 886d7b7cdc
commit ac5337130a
2 changed files with 35 additions and 0 deletions

View file

@ -604,6 +604,31 @@ void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS)
bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base)
{
if (!m_surf) {
init();
}
// Check that the number of surface sites is balanced
double reac_sites = 0.0;
double prod_sites = 0.0;
for (const auto& reactant : r_base->reactants) {
size_t k = m_surf->speciesIndex(reactant.first);
if (k != npos) {
reac_sites += reactant.second * m_surf->size(k);
}
}
for (const auto& product : r_base->products) {
size_t k = m_surf->speciesIndex(product.first);
if (k != npos) {
prod_sites += product.second * m_surf->size(k);
}
}
if (fabs(reac_sites - prod_sites) > 1e-5 * (reac_sites + prod_sites)) {
throw CanteraError("InterfaceKinetics::addReaction", "Number of surface"
" sites not balanced in reaction {}.\nReactant sites: {}\n"
"Product sites: {}", r_base->equation(), reac_sites, prod_sites);
}
size_t i = nReactions();
bool added = Kinetics::addReaction(r_base);
if (!added) {

View file

@ -379,6 +379,16 @@ TEST_F(InterfaceKineticsFromScratch, add_sticking_reaction)
check_rates(0);
}
TEST_F(InterfaceKineticsFromScratch, unbalanced_sites)
{
Composition reac = parseCompString("H(m):1 O(m):1");
Composition prod = parseCompString("OH(m):1");
Arrhenius rate(5e21, 0, 100.0e6 / GasConstant);
auto R = make_shared<InterfaceReaction>(reac, prod, rate);
ASSERT_THROW(kin.addReaction(R), CanteraError);
}
class KineticsAddSpecies : public testing::Test
{
public: