diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 890de5b0f..a961b0f3e 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -168,7 +168,7 @@ class InterfaceReaction : public Reaction { public: InterfaceReaction(const Composition& reactants, const Composition& products, - const Arrhenius& rate); + const Arrhenius& rate, bool isStick=false); //! Adjustments to the Arrhenius rate expression dependent on surface //! species coverages. Three coverage parameters (a, E, m) are used for each @@ -179,6 +179,10 @@ public: //! The rate coefficient, without taking into account the coverage //! dependencies. Arrhenius rate; + + // Set to true if `rate` is a parameterization of the sticking coefficient + // rather than the forward rate constant + bool is_sticking_coefficient; }; diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 909594b5b..ed717c44c 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -18,8 +18,6 @@ namespace Cantera { class Array2D; -class Plog; -class SurfaceArrhenius; //! Arrhenius reaction rate type depends only on temperature /** @@ -96,7 +94,19 @@ public: s << ");" << std::endl; } - //! @deprecated. To be removed after Cantera 2.2 + //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending + //! on the reaction order) + double preExponentialFactor() const { + return m_A; + } + + //! Return the temperature exponent *b* + double temperatureExponent() const { + return m_b; + } + + //! Return the activation energy divided by the gas constant (i.e. the + //! activation temperature) [K] doublereal activationEnergy_R() const { return m_E; } @@ -108,9 +118,6 @@ public: protected: doublereal m_logA, m_b, m_E, m_A; - - friend class Plog; - friend class SurfaceArrhenius; }; @@ -137,7 +144,7 @@ public: } SurfaceArrhenius(); - explicit SurfaceArrhenius(const Arrhenius& rate); + explicit SurfaceArrhenius(double A, double b, double Ta); explicit SurfaceArrhenius(const ReactionData& rdata); void addCoverageDependence(size_t k, doublereal a, @@ -187,7 +194,6 @@ public: (m_E + m_ecov)*recipT + m_mcov); } - //! @deprecated. To be removed after Cantera 2.2 doublereal activationEnergy_R() const { return m_E + m_ecov; } diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index febb6a0d0..5cca5e941 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -900,7 +900,54 @@ 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); + double A_rate = r.rate.preExponentialFactor(); + double b_rate = r.rate.temperatureExponent(); + + if (r.is_sticking_coefficient) { + // Identify the interface phase + size_t iInterface = npos; + size_t min_dim = 4; + for (size_t i = 0; i < nPhases(); i++) { + if (thermo(i).nDim() < min_dim) { + iInterface = i; + min_dim = thermo(i).nDim(); + } + } + + b_rate += 0.5; + + // Identify the sticking species, and adjust the A-factor + bool foundStick = false; + for (Composition::const_iterator iter = r.reactants.begin(); + iter != r.reactants.end(); + ++iter) { + size_t iPhase = speciesPhaseIndex(kineticsSpeciesIndex(iter->first)); + const ThermoPhase& p = thermo(iPhase); + size_t k = p.speciesIndex(iter->first); + if (iPhase == iInterface) { + // Interface species. Convert from coverages used in the + // sticking probability expression to the concentration units + // used in the mass action rate expression + double order = getValue(r.orders, iter->first, iter->second); + A_rate /= pow(p.standardConcentration(k), order); + } else { + // Non-interface species. There should be exactly one of these + if (foundStick) { + throw CanteraError("InterfaceKinetics::addReaction", + "Multiple non-interface species found" + "in sticking reaction: '" + r.equation() + "'"); + } + foundStick = true; + A_rate *= sqrt(GasConstant/(2*Pi*p.molecularWeight(k))); + } + } + if (!foundStick) { + throw CanteraError("InterfaceKinetics::addReaction", + "No non-interface species found" + "in sticking reaction: '" + r.equation() + "'"); + } + } + SurfaceArrhenius rate(A_rate, b_rate, r.rate.activationEnergy_R()); // Turn on the global flag indicating surface coverage dependence if (!r.coverage_deps.empty()) { diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 1c9e7cc00..ad39bb533 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -95,9 +95,11 @@ ChebyshevReaction::ChebyshevReaction(const Composition& reactants_, InterfaceReaction::InterfaceReaction(const Composition& reactants_, const Composition& products_, - const Arrhenius& rate_) + const Arrhenius& rate_, + bool isStick) : Reaction(INTERFACE_RXN, reactants_, products_) , rate(rate_) + , is_sticking_coefficient(isStick) { } diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 7f5d7c432..4ca7d7bf8 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -50,11 +50,11 @@ SurfaceArrhenius::SurfaceArrhenius() { } -SurfaceArrhenius::SurfaceArrhenius(const Arrhenius& rate) - : m_logA(rate.m_logA) - , m_b(rate.m_b) - , m_E(rate.m_E) - , m_A(rate.m_A) +SurfaceArrhenius::SurfaceArrhenius(double A, double b, double Ta) + : m_logA(std::log(A)) + , m_b(b) + , m_E(Ta) + , m_A(A) , m_acov(0.0) , m_ecov(0.0) , m_mcov(0.0) @@ -223,9 +223,9 @@ Plog::Plog(const std::multimap& rates) maxRates_ = std::max(rateCount, maxRates_); j++; - A_.push_back(iter->second.m_A); - n_.push_back(iter->second.m_b); - Ea_.push_back(iter->second.m_E); + A_.push_back(iter->second.preExponentialFactor()); + n_.push_back(iter->second.temperatureExponent()); + Ea_.push_back(iter->second.activationEnergy_R()); } // For pressures with only one Arrhenius expression, it is more diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index ee8610e5f..3c35cd8ae 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -225,3 +225,18 @@ TEST_F(InterfaceKineticsFromScratch, add_surface_reaction) kin.finalize(); check_rates(3); } + +TEST_F(InterfaceKineticsFromScratch, add_sticking_reaction) +{ + // Reaction 0 on the metal surface + // surface_reaction( "H2 + (m) + (m) <=> H(m) + H(m)", + // stick(0.1, 0, 0), id = 'metal-rxn1') + Composition reac = parseCompString("H2:1 (m):2"); + Composition prod = parseCompString("H(m):2"); + Arrhenius rate(0.1, 0, 0.0); + + shared_ptrR(new InterfaceReaction(reac, prod, rate, true)); + kin.addReaction(R); + kin.finalize(); + check_rates(0); +}