Add Blowers Masel reaction to GasKinetics and InterfaceKinetics

AUTHORS

@@ -52,4 +52,5 @@ Laurien Vandewalle (@lavdwall)
 Bryan Weber (@bryanwweber), University of Connecticut
 Armin Wehrfritz (@awehrfritz)
 Richard West (@rwest), Northeastern University
+Chao Xu (@12Chao), Northeastern University
 Thorsten Zirwes (@g3bk47), Karlsruhe Institute of Technology

README.rst

@@ -89,7 +89,7 @@ possible.
 Development Site
 ================
 
-The current development version is 2.6.0a1. The current stable version is
+The current development version is 2.6.0a2. The current stable version is
 2.5.1. The `latest Cantera source code <https://github.com/Cantera/cantera>`_,
 the `issue tracker <https://github.com/Cantera/cantera/issues>`_ for bugs and
 enhancement requests, `downloads of Cantera releases and binary installers

SConstruct

@@ -708,7 +708,7 @@ for arg in ARGUMENTS:
         print('Encountered unexpected command line argument: %r' % arg)
         sys.exit(1)
 
-env['cantera_version'] = "2.6.0a1"
+env["cantera_version"] = "2.6.0a2"
 # For use where pre-release tags are not permitted (MSI, sonames)
 env['cantera_pure_version'] = re.match(r'(\d+\.\d+\.\d+)', env['cantera_version']).group(0)
 env['cantera_short_version'] = re.match(r'(\d+\.\d+)', env['cantera_version']).group(0)

doc/doxygen/Doxyfile

@@ -38,7 +38,7 @@ PROJECT_NAME           = Cantera
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = 2.6.0a1
+PROJECT_NUMBER         = 2.6.0a2
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a

doc/sphinx/cython/kinetics.rst

@@ -58,11 +58,21 @@ ChebyshevReaction
 .. autoclass:: ChebyshevReaction(reactants='', products='')
    :no-undoc-members:
 
+BlowersMaselReaction
+^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: BlowersMaselReaction(reactants='', products='')
+   :no-undoc-members:
+
 InterfaceReaction
 ^^^^^^^^^^^^^^^^^
 .. autoclass:: InterfaceReaction(reactants='', products='')
    :no-undoc-members:
 
+BlowersMaselInterfaceReaction
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: BlowersMaselInterfaceReaction(reactants='', products='')
+   :no-undoc-members:
+
 Auxilliary Reaction Data
 ------------------------
 
@@ -85,6 +95,10 @@ SriFalloff
 .. autoclass:: SriFalloff(coeffs=(), init=True)
    :no-undoc-members:
 
+BlowersMasel
+^^^^^^^^^^^^
+.. autoclass:: BlowersMasel(A, b, E0, w)
+
 Reaction Path Analysis
 ----------------------
 

doc/sphinx/yaml/reactions.rst

@@ -26,12 +26,17 @@ The fields common to all ``reaction`` entries are:
     - :ref:`chemically-activated <sec-yaml-chemically-activated>`
     - :ref:`pressure-dependent-Arrhenius <sec-yaml-pressure-dependent-Arrhenius>`
     - :ref:`Chebyshev <sec-yaml-Chebyshev>`
+    - :ref:`Blowers-Masel <sec-yaml-Blowers-Masel>`
+    - :ref:`surface-Blowers-Masel <sec-yaml-surface-Blowers-Masel>`
+
+    Reactions without a specified ``type`` on surfaces or edges are
+    automatically treated as :ref:`interface <sec-yaml-interface-reaction>`
+    reactions, and reactions that involve charge transfer between phases are
+    automatically treated as :ref:`electrochemical <sec-yaml-electrochemical-reaction>`
+    reactions. Reactions on surfaces or edges specifying ``type`` as
+    ``Blowers-Masel`` are treated as
+    :ref:`surface-Blowers-Masel <sec-yaml-surface-Blowers-Masel>`.
 
-    Reactions on surfaces or edges are automatically treated as
-    :ref:`interface <sec-yaml-interface-reaction>` reactions, and reactions that
-    involve charge transfer between phases are automatically treated as
-    :ref:`electrochemical <sec-yaml-electrochemical-reaction>` reactions, without the
-    need to specify the ``type``.
 
 ``duplicate``
     Boolean indicating whether the reaction is a known duplicate of another
@@ -255,6 +260,32 @@ Example::
            [-2.26210e-01,  1.69190e-01,  4.85810e-03, -2.38030e-03],
            [-1.43220e-01,  7.71110e-02,  1.27080e-02, -6.41540e-04]]
 
+.. _sec-yaml-Blowers-Masel:
+
+``Blowers-Masel``
+-----------------
+
+A reaction with parameters to calculate rate constant based on Blowers Masel
+approximation as `described here <https://cantera.org/science/reactions.html#sec-blowers-masel>`__.
+
+Additional fields are:
+
+``rate-constant``
+    A list of values containing the pre-exponential factor :math:`A`, the
+    temperature exponent :math:`b`, the intrinsic activation energy :math:`E_{a0}`,
+    and the average of the bond dissociation energy of the bond breaking and that
+    being formed in the reaction :math:`w`.
+
+``negative-A``
+    A boolean indicating whether a negative value for the pre-exponential factor
+    is allowed. The default is ``false``.
+
+Example::
+
+    equation: O + H2 <=> H + OH
+    type: Blowers-Masel
+    rate-constant: {A: 3.87e+04 cm^2/mol/s, b: 2.7, Ea0: 6260.0 cal/mol, w: 1e9 cal/mol}
+
 
 .. _sec-yaml-interface-reaction:
 
@@ -263,7 +294,7 @@ Example::
 
 A reaction occuring on a surface between two bulk phases, or along an edge
 at the intersection of two surfaces, as
-`described here <https://cantera.org/science/reactions.html#surface-reactions>`__.
+`described here <https://cantera.org/science/reactions.html#sec-surface>`__.
 
 Includes the fields of an :ref:`sec-yaml-elementary` reaction plus:
 
@@ -321,3 +352,23 @@ Example::
     equation: LiC6 <=> Li+(e) + C6
     rate-constant: [5.74, 0.0, 0.0]
     beta: 0.4
+
+.. _sec-yaml-surface-Blowers-Masel:
+
+``surface-Blowers-Masel``
+-------------------------
+
+A reaction occurring on a surface between two bulk phases, or along an edge
+at the intersection of two surfaces, which the rate constant can be calculated
+by Blowers Masel Approximation with Arrhenius expression as
+`described here <https://cantera.org/science/reactions.html#surface-blowers-masel-reactions>`__.
+
+Includes the fields of a :ref:`sec-yaml-Blowers-Masel` reaction and
+the fields of an :ref:`sec-yaml-interface-reaction` reaction.
+
+Example::
+
+    equation: 2 H(s) => H2 + 2 Pt(s)
+    type: Blowers-Masel
+    rate-constant: {A: 3.7e21 cm^2/mol/s, b: 0, Ea0: 67400 J/mol, w: 1000000 J/mol}
+    coverage-dependencies: {H(s): {a: 0, m: 0, E: -6000 J/mol}}

include/cantera/kinetics/GasKinetics.h

@@ -88,6 +88,7 @@ class GasKinetics : public BulkKinetics
 
     Rate1<Plog> m_plog_rates;
     Rate1<ChebyshevRate> m_cheb_rates;
+    Rate1<BlowersMasel> m_blowersmasel_rates;
 
     //! @name Reaction rate data
     //!@{
@@ -109,11 +110,13 @@ class GasKinetics : public BulkKinetics
     void addFalloffReaction(FalloffReaction& r);
     void addPlogReaction(PlogReaction& r);
     void addChebyshevReaction(ChebyshevReaction& r);
+    void addBlowersMaselReaction(BlowersMaselReaction& r);
 
     void modifyThreeBodyReaction(size_t i, ThreeBodyReaction& r);
     void modifyFalloffReaction(size_t i, FalloffReaction& r);
     void modifyPlogReaction(size_t i, PlogReaction& r);
     void modifyChebyshevReaction(size_t i, ChebyshevReaction& r);
+    void modifyBlowersMaselReaction(size_t i, BlowersMaselReaction& r);
 
     //! Update the equilibrium constants in molar units.
     void updateKc();

include/cantera/kinetics/InterfaceKinetics.h

@@ -19,6 +19,7 @@ namespace Cantera
 class SurfPhase;
 class ImplicitSurfChem;
 class InterfaceReaction;
+class BlowersMaselInterfaceReaction;
 
 //! A kinetics manager for heterogeneous reaction mechanisms. The reactions are
 //! assumed to occur at a 2D interface between two 3D phases.
@@ -384,6 +385,17 @@ class InterfaceKinetics : public Kinetics
     SurfaceArrhenius buildSurfaceArrhenius(size_t i, InterfaceReaction& r,
                                            bool replace);
 
+    //! Build a BMSurfaceArrhenius object from a Reaction, taking into account
+    //! the possible sticking coefficient form and coverage dependencies
+    //! @param i  Reaction number
+    //! @param r  Reaction object containing rate coefficient parameters
+    //! @param replace  True if replacing an existing reaction
+    //! @todo This function duplicated most of the code from buildSurfaceArrhenius
+    //! to return a slightly different reaction rate class and could be refactored in the
+    //! future.
+    BMSurfaceArrhenius buildBMSurfaceArrhenius(size_t i, BlowersMaselInterfaceReaction& r,
+                                               bool replace);
+
     //! Temporary work vector of length m_kk
     vector_fp m_grt;
 
@@ -401,6 +413,13 @@ class InterfaceKinetics : public Kinetics
      */
     Rate1<SurfaceArrhenius> m_rates;
 
+    //! Templated class containing the vector of surface Blowers Masel reactions for this interface
+    /*!
+     *  The templated class is described in RateCoeffMgr.h
+     *  The class BMSurfaceArrhenius is described in RxnRates.h
+     */
+    Rate1<BMSurfaceArrhenius> m_blowers_masel_rates;
+
     bool m_redo_rates;
 
     //! Vector of irreversible reaction numbers

include/cantera/kinetics/Kinetics.h

@@ -945,6 +945,9 @@ class Kinetics
     //! Net rate-of-progress for each reaction
     vector_fp m_ropnet;
 
+    //! The enthalpy change for each reaction to calculate Blowers-Masel rates
+    vector_fp m_dH;
+
     //! @see skipUndeclaredSpecies()
     bool m_skipUndeclaredSpecies;
 

include/cantera/kinetics/RateCoeffMgr.h

@@ -69,6 +69,13 @@ class Rate1
         }
     }
 
+    void updateBlowersMasel(double T, double logT, double* values, double* deltaH) {
+        double recipT = 1.0/T;
+        for (size_t i=0; i != m_rates.size(); i++) {
+            values[m_rxn[i]] = m_rates[i].updateRC(logT, recipT, deltaH[m_rxn[i]]);
+        }
+    }
+
     size_t nReactions() const {
         return m_rates.size();
     }

include/cantera/kinetics/Reaction.h

@@ -374,6 +374,7 @@ class InterfaceReaction : public ElementaryReaction
     InterfaceReaction();
     InterfaceReaction(const Composition& reactants, const Composition& products,
                       const Arrhenius& rate, bool isStick=false);
+    virtual void calculateRateCoeffUnits(const Kinetics& kin);
     virtual void getParameters(AnyMap& reactionNode) const;
 
     virtual std::string type() const {
@@ -416,6 +417,60 @@ class ElectrochemicalReaction : public InterfaceReaction
     bool exchange_current_density_formulation;
 };
 
+//! A reaction with rate parameters for Blowers-Masel approximation
+class BlowersMaselReaction: public Reaction
+{
+public:
+    BlowersMaselReaction();
+    BlowersMaselReaction(const Composition& reactants,
+                         const Composition& products, const BlowersMasel& rate);
+    virtual void getParameters(AnyMap& reactionNode) const;
+    virtual void validate();
+
+    virtual std::string type() const {
+        return "Blowers-Masel";
+    }
+
+    BlowersMasel rate;
+
+    bool allow_negative_pre_exponential_factor;
+};
+
+//! A reaction occurring on an interface (i.e. a SurfPhase or an EdgePhase)
+//! with the rate calculated with Blowers-Masel approximation.
+class BlowersMaselInterfaceReaction : public BlowersMaselReaction
+{
+public:
+    BlowersMaselInterfaceReaction();
+    BlowersMaselInterfaceReaction(const Composition& reactants, const Composition& products,
+                      const BlowersMasel& rate, bool isStick=false);
+    virtual void getParameters(AnyMap& reactionNode) const;
+    virtual void calculateRateCoeffUnits(const Kinetics& kin);
+
+    virtual std::string type() const {
+        return "surface-Blowers-Masel";
+    }
+    //! Adjustments to the Arrhenius rate expression dependent on surface
+    //! species coverages. Three coverage parameters (a, E, m) are used for each
+    //! species on which the rate depends. See SurfaceArrhenius for details on
+    //! the parameterization.
+    std::map<std::string, CoverageDependency> coverage_deps;
+
+    //! Set to true if `rate` is a parameterization of the sticking coefficient
+    //! rather than the forward rate constant
+    bool is_sticking_coefficient;
+
+    //! Set to true if `rate` is a sticking coefficient which should be
+    //! translated into a rate coefficient using the correction factor developed
+    //! by Motz & Wise for reactions with high (near-unity) sticking
+    //! coefficients. Defaults to 'false'.
+    bool use_motz_wise_correction;
+
+    //! For reactions with multiple non-surface species, the sticking species
+    //! needs to be explicitly identified.
+    std::string sticking_species;
+};
+
 //! Create Reaction objects for all `<reaction>` nodes in an XML document.
 //!
 //! The `<reaction>` nodes are assumed to be children of the `<reactionData>`
@@ -488,6 +543,12 @@ void setupElectrochemicalReaction(ElectrochemicalReaction&,
 void setupElectrochemicalReaction(ElectrochemicalReaction&,
                                   const AnyMap&, const Kinetics&);
 
+//! @internal May be changed without notice in future versions
+void setupBlowersMaselReaction(BlowersMaselReaction&,
+                               const AnyMap&, const Kinetics&);
+//! @internal May be changed without notice in future versions
+void setupBlowersMaselInterfaceReaction(BlowersMaselInterfaceReaction&,
+                                        const AnyMap&, const Kinetics&);
 }
 
 #endif