diff --git a/include/cantera/kinetics/importKinetics.h b/include/cantera/kinetics/importKinetics.h index 5793550c3..f5bb7777e 100644 --- a/include/cantera/kinetics/importKinetics.h +++ b/include/cantera/kinetics/importKinetics.h @@ -110,6 +110,26 @@ bool importKinetics(const XML_Node& phase, std::vector th, */ bool buildSolutionFromXML(XML_Node& root, const std::string& id, const std::string& nm, ThermoPhase* th, Kinetics* kin); + +//! Check to ensure that all electrochemical reactions are specified correctly +/*! + * This function ensures the user has correctly specified all electrochemical + * reactions. The routine counts the amount of charge (i.e. number of electron + * elements specified for each species in each phase) for both reactants and + * products. If net charge transfer phases during a reaction, the reaction is + * electrochemical. If not already specified as such, the function defines the + * reaction as electrochemical, corrects the reaction attributes, and sets + * beta = 0.5. + * + * @param p This is an XML node containing a description of the owning + * phase for the kinetics object. + * @param kin This is a pointer to a kinetics manager class. + * @param r This is the reaction node that is being evaluated + * @return The function always returns true. +*/ +bool checkElectrochemReaction(const XML_Node& p, Kinetics& kin, const XML_Node& r); + + } #endif diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 470e4de3e..0eac0d4aa 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -255,7 +255,7 @@ InterfaceReaction::InterfaceReaction(const Composition& reactants_, ElectrochemicalReaction::ElectrochemicalReaction() : film_resistivity(0.0) - , beta(0.0) + , beta(0.5) , exchange_current_density_formulation(false) { } @@ -265,12 +265,11 @@ ElectrochemicalReaction::ElectrochemicalReaction(const Composition& reactants_, const Arrhenius& rate_) : InterfaceReaction(reactants_, products_, rate_) , film_resistivity(0.0) - , beta(0.0) + , beta(0.5) , exchange_current_density_formulation(false) { } - Arrhenius readArrhenius(const XML_Node& arrhenius_node) { return Arrhenius(getFloat(arrhenius_node, "A", "toSI"), diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index 9aeb4f837..19c8ac45c 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -79,6 +79,7 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin, // if no 'include' directive, then include all reactions if (incl.empty()) { for (size_t i = 0; i < allrxns.size(); i++) { + checkElectrochemReaction(p,kin,*allrxns[i]); kin.addReaction(newReaction(*allrxns[i])); ++itot; } @@ -110,6 +111,7 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin, // do a lexical min max and operation. This sometimes // has surprising results. if ((rxid >= imin) && (rxid <= imax)) { + checkElectrochemReaction(p,kin,*r); kin.addReaction(newReaction(*r)); ++itot; } @@ -146,7 +148,7 @@ bool importKinetics(const XML_Node& phase, std::vector th, } } - // if other phases are involved in the reaction mechanism, they must be + // If other phases are involved in the reaction mechanism, they must be // listed in a 'phaseArray' child element. Homogeneous mechanisms do not // need to include a phaseArray element. vector phase_ids; @@ -212,4 +214,87 @@ bool buildSolutionFromXML(XML_Node& root, const std::string& id, return true; } +bool checkElectrochemReaction(const XML_Node& p, Kinetics& kin, const XML_Node& r) +{ + // If other phases are involved in the reaction mechanism, they must be + // listed in a 'phaseArray' child element. Homogeneous mechanisms do not + // need to include a phaseArray element. + vector phase_ids; + if (p.hasChild("phaseArray")) { + const XML_Node& pa = p.child("phaseArray"); + getStringArray(pa, phase_ids); + } + phase_ids.push_back(p["id"]); + + // Get reaction product and reactant information + Composition reactants = parseCompString(r.child("reactants").value()); + Composition products = parseCompString(r.child("products").value()); + + + // If the reaction has undeclared species don't perform electrochemical check + for (const auto& sp : reactants) { + if (kin.kineticsSpeciesIndex(sp.first) == npos) { + return true; + } + } + + for (const auto& sp : products) { + if (kin.kineticsSpeciesIndex(sp.first) == npos) { + return true; + } + } + + // Initialize the electron counter for each phase + std::vector e_counter(phase_ids.size(), 0.0); + + // Find the amount of electrons in the products for each phase + for (const auto& sp : products) { + const ThermoPhase& ph = kin.speciesPhase(sp.first); + size_t k = ph.speciesIndex(sp.first); + double stoich = sp.second; + for (size_t m = 0; m < phase_ids.size(); m++) { + if (phase_ids[m] == ph.id()) { + e_counter[m] += stoich * ph.charge(k); + break; + } + } + } + + // Subtract the amount of electrons in the reactants for each phase + for (const auto& sp : reactants) { + const ThermoPhase& ph = kin.speciesPhase(sp.first); + size_t k = ph.speciesIndex(sp.first); + double stoich = sp.second; + for (size_t m = 0; m < phase_ids.size(); m++) { + if (phase_ids[m] == ph.id()) { + e_counter[m] -= stoich * ph.charge(k); + break; + } + } + } + + // If the electrons change phases then the reaction is electrochemical + bool echemical = false; + for(size_t m = 0; m < phase_ids.size(); m++) { + if (fabs(e_counter[m]) > 1e-4) { + echemical = true; + break; + } + } + + // If the reaction is electrochemical, ensure the reaction is identified as + // electrochemical. If not already specified beta is assumed to be 0.5 + std::string type = ba::to_lower_copy(r["type"]); + if (!r.child("rateCoeff").hasChild("electrochem")) { + if ((type != "butlervolmer_noactivitycoeffs" && + type != "butlervolmer" && + type != "surfaceaffinity") && + echemical) { + XML_Node& f = r.child("rateCoeff").addChild("electrochem",""); + f.addAttribute("beta",0.5); + } + } + return true; +} + } diff --git a/test/data/sofc-test.xml b/test/data/sofc-test.xml index 1c5570d3f..ea7bbaf2a 100644 --- a/test/data/sofc-test.xml +++ b/test/data/sofc-test.xml @@ -412,7 +412,6 @@ H(m) + O''(ox) [=] (m) + electron + OH'(ox) - 5.000000E+10 0.0