diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index e094fb92b..bde3f730f 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -11,7 +11,7 @@ #include "cantera/thermo/mix_defs.h" #include "Kinetics.h" #include "cantera/kinetics/RxnMolChange.h" - +#include "Reaction.h" #include "cantera/base/utilities.h" #include "RateCoeffMgr.h" @@ -210,6 +210,7 @@ public: virtual void init(); virtual void addReaction(ReactionData& r); + virtual void addReaction(shared_ptr r); virtual void finalize(); virtual bool ready() const; //! @} @@ -389,8 +390,12 @@ public: int phaseStability(const size_t iphase) const; virtual void determineFwdOrdersBV(ReactionData& rdata, vector_fp& fwdFullorders); + virtual void determineFwdOrdersBV(ElectrochemicalReaction& r, vector_fp& fwdFullorders); protected: + void addElementaryReaction(InterfaceReaction& rdata); + void addGlobalReaction(InterfaceReaction& r); + //! Temporary work vector of length m_kk vector_fp m_grt; diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 0f461265f..890de5b0f 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -209,6 +209,8 @@ public: //! Forward value of the apparent Electrochemical transfer coefficient doublereal beta; + + bool exchange_current_density_formulation; }; } diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index c3d8b5774..febb6a0d0 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -896,6 +896,122 @@ void InterfaceKinetics::addReaction(ReactionData& r) } } +void InterfaceKinetics::addReaction(shared_ptr r_base) +{ + InterfaceReaction& r = dynamic_cast(*r_base); + // Create a SurfaceArrhenius rate calculator and set the coverage dependencies + SurfaceArrhenius rate(r.rate); + + // Turn on the global flag indicating surface coverage dependence + if (!r.coverage_deps.empty()) { + m_has_coverage_dependence = true; + } + + for (map::const_iterator iter = r.coverage_deps.begin(); + iter != r.coverage_deps.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + rate.addCoverageDependence(k, iter->second.a, iter->second.m, iter->second.E); + } + + m_rates.install(m_ii, rate); + + // Store activation energy + m_E.push_back(rate.activationEnergy_R()); + + ElectrochemicalReaction* re = dynamic_cast(&r); + if (re) { + m_has_electrochem_rxns = true; + m_beta.push_back(re->beta); + m_ctrxn.push_back(m_ii); + if (re->exchange_current_density_formulation) { + m_has_exchange_current_density_formulation = true; + m_ctrxn_ecdf.push_back(1); + } else { + m_ctrxn_ecdf.push_back(0); + } + m_ctrxn_resistivity_.push_back(re->film_resistivity); + + if (r.reaction_type == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN || + r.reaction_type == BUTLERVOLMER_RXN || + r.reaction_type == SURFACEAFFINITY_RXN || + r.reaction_type == GLOBAL_RXN) { + // Specify alternative forms of the electrochemical reaction + if (r.reaction_type == BUTLERVOLMER_RXN) { + m_ctrxn_BVform.push_back(1); + } else if (r.reaction_type == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN) { + m_ctrxn_BVform.push_back(2); + } else { + // set the default to be the normal forward / reverse calculation method + m_ctrxn_BVform.push_back(0); + } + if (!r.orders.empty()) { + vector_fp orders(nTotalSpecies(), 0.0); + for (Composition::const_iterator iter = r.orders.begin(); + iter != r.orders.end(); + ++iter) { + orders[kineticsSpeciesIndex(iter->first)] = iter->second; + } + RxnOrders* ro = new RxnOrders(); + ro->fill(orders); + m_ctrxn_ROPOrdersList_.push_back(ro); + m_ctrxn_FwdOrdersList_.push_back(0); + + // Fill in the Fwd Orders dependence here for B-V reactions + if (r.reaction_type == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN || + r.reaction_type == BUTLERVOLMER_RXN) { + vector_fp fwdFullorders(m_kk, 0.0); + determineFwdOrdersBV(*re, fwdFullorders); + RxnOrders* ro = new RxnOrders(); + ro->fill(fwdFullorders); + m_ctrxn_FwdOrdersList_[m_ii] = ro; + } + } else { + m_ctrxn_ROPOrdersList_.push_back(0); + m_ctrxn_FwdOrdersList_.push_back(0); + } + + } else { + m_ctrxn_BVform.push_back(0); + m_ctrxn_ROPOrdersList_.push_back(0); + m_ctrxn_FwdOrdersList_.push_back(0); + if (re->film_resistivity > 0.0) { + throw CanteraError("InterfaceKinetics::addReaction()", + "film resistivity set for elementary reaction"); + } + } + } + + if (r.reversible) { + m_revindex.push_back(nReactions()); + m_nrev++; + } else { + m_irrev.push_back(nReactions()); + m_nirrev++; + } + Kinetics::addReaction(r_base); + + m_rxnPhaseIsReactant.push_back(std::vector(nPhases(), false)); + m_rxnPhaseIsProduct.push_back(std::vector(nPhases(), false)); + + size_t i = m_ii - 1; + for (Composition::const_iterator iter = r.reactants.begin(); + iter != r.reactants.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + size_t p = speciesPhaseIndex(k); + m_rxnPhaseIsReactant[i][p] = true; + } + for (Composition::const_iterator iter = r.products.begin(); + iter != r.products.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + size_t p = speciesPhaseIndex(k); + m_rxnPhaseIsProduct[i][p] = true; + } +} + + void InterfaceKinetics::setIOFlag(int ioFlag) { m_ioFlag = ioFlag; @@ -1102,6 +1218,49 @@ void InterfaceKinetics::determineFwdOrdersBV(ReactionData& rdata, std::vector& fwdFullOrders) +{ + // Start out with the full ROP orders vector. + // This vector will have the BV exchange current density orders in it. + fwdFullOrders.assign(nTotalSpecies(), 0.0); + for (Composition::const_iterator iter = r.orders.begin(); + iter != r.orders.end(); + ++iter) { + fwdFullOrders[kineticsSpeciesIndex(iter->first)] = iter->second; + } + + // forward and reverse beta values + double betaf = r.beta; + double betar = 1.0 - betaf; + + // Loop over the reactants doing away with the BV terms. + // This should leave the reactant terms only, even if they are non-mass action. + for (Composition::const_iterator iter = r.reactants.begin(); + iter != r.reactants.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + fwdFullOrders[k] += betaf * iter->second; + // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet) + if (abs(fwdFullOrders[k]) < 0.00001) { + fwdFullOrders[k] = 0.0; + } + } + + // Loop over the products doing away with the BV terms. + // This should leave the reactant terms only, even if they are non-mass action. + for (Composition::const_iterator iter = r.products.begin(); + iter != r.products.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + fwdFullOrders[k] -= betaf * iter->second; + // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet) + if (abs(fwdFullOrders[k]) < 0.00001) { + fwdFullOrders[k] = 0.0; + } + } +} + + void EdgeKinetics::finalize() { // Note we can't call the Interface::finalize() routine because we need to check for a dimension of 1 below. diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 3772c1765..1c9e7cc00 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -109,6 +109,7 @@ ElectrochemicalReaction::ElectrochemicalReaction(const Composition& reactants_, , equilibrium_constant_power(1.0) , affinity_power(1.0) , beta(0.0) + , exchange_current_density_formulation(false) { } diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index dcd1fd792..ee8610e5f 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -1,7 +1,9 @@ #include "gtest/gtest.h" #include "cantera/kinetics/importKinetics.h" #include "cantera/thermo/IdealGasPhase.h" +#include "cantera/thermo/SurfPhase.h" #include "cantera/kinetics/GasKinetics.h" +#include "cantera/kinetics/InterfaceKinetics.h" #include "cantera/base/Array.h" using namespace Cantera; @@ -158,3 +160,68 @@ TEST_F(KineticsFromScratch, add_chebyshev_reaction) kin.finalize(); check_rates(4); } + +class InterfaceKineticsFromScratch : public testing::Test +{ +public: + InterfaceKineticsFromScratch() + : gas("../data/sofc-test.xml", "gas") + , gas_ref("../data/sofc-test.xml", "gas") + , surf("../data/sofc-test.xml", "metal_surface") + , surf_ref("../data/sofc-test.xml", "metal_surface") + { + std::vector th; + th.push_back(&surf_ref); + th.push_back(&gas_ref); + importKinetics(surf_ref.xml(), th, &kin_ref); + kin.addPhase(surf); + kin.addPhase(gas); + kin.init(); + } + + IdealGasPhase gas; + IdealGasPhase gas_ref; + SurfPhase surf; + SurfPhase surf_ref; + InterfaceKinetics kin; + InterfaceKinetics kin_ref; + + //! iRef is the index of the corresponding reaction in the reference mech + void check_rates(int iRef) { + ASSERT_EQ((size_t) 1, kin.nReactions()); + + std::string X = "H2:0.2 O2:0.5 H2O:0.1 N2:0.2"; + std::string Xs = "H(m):0.1 O(m):0.2 OH(m):0.3 (m):0.4"; + gas.setState_TPX(1200, 5*OneAtm, X); + gas_ref.setState_TPX(1200, 5*OneAtm, X); + surf.setState_TP(1200, 5*OneAtm); + surf_ref.setState_TP(1200, 5*OneAtm); + surf.setCoveragesByName(Xs); + surf_ref.setCoveragesByName(Xs); + + vector_fp k(1), k_ref(kin_ref.nReactions()); + + kin.getFwdRateConstants(&k[0]); + kin_ref.getFwdRateConstants(&k_ref[0]); + EXPECT_FLOAT_EQ(k_ref[iRef], k[0]); + + kin.getRevRateConstants(&k[0]); + kin_ref.getRevRateConstants(&k_ref[0]); + EXPECT_FLOAT_EQ(k_ref[iRef], k[0]); + } +}; + +TEST_F(InterfaceKineticsFromScratch, add_surface_reaction) +{ + // Reaction 3 on the metal surface + // surface_reaction( "H(m) + O(m) <=> OH(m) + (m)", + // [5.00000E+22, 0, 100.0], id = 'metal-rxn4') + Composition reac = parseCompString("H(m):1 O(m):1"); + Composition prod = parseCompString("OH(m):1 (m):1"); + Arrhenius rate(5e21, 0, 100.0e6 / GasConstant); // kJ/mol -> J/kmol + + shared_ptrR(new InterfaceReaction(reac, prod, rate)); + kin.addReaction(R); + kin.finalize(); + check_rates(3); +}