[Kinetics] InterfaceReaction supports sticking coefficients

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;
 };
 
 

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;
     }

src/kinetics/InterfaceKinetics.cpp

@@ -900,7 +900,54 @@ void InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base)
 {
     InterfaceReaction& r = dynamic_cast<InterfaceReaction&>(*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()) {

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)
 {
 }
 

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<double, Arrhenius>& 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

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_ptr<InterfaceReaction>R(new InterfaceReaction(reac, prod, rate, true));
+    kin.addReaction(R);
+    kin.finalize();
+    check_rates(0);
+}