diff --git a/AUTHORS b/AUTHORS index 599256f18d4..bb379deb596 100644 --- a/AUTHORS +++ b/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 diff --git a/README.rst b/README.rst index 37c98561277..4c7b81d1532 100644 --- a/README.rst +++ b/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 `_, the `issue tracker `_ for bugs and enhancement requests, `downloads of Cantera releases and binary installers diff --git a/SConstruct b/SConstruct index 7e9ce0457f5..21f3b530c5f 100644 --- a/SConstruct +++ b/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) diff --git a/doc/doxygen/Doxyfile b/doc/doxygen/Doxyfile index 5509fe699b4..e0a1cf47409 100644 --- a/doc/doxygen/Doxyfile +++ b/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 diff --git a/doc/sphinx/cython/kinetics.rst b/doc/sphinx/cython/kinetics.rst index fea8f723d2b..c31cc0fb496 100644 --- a/doc/sphinx/cython/kinetics.rst +++ b/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 ---------------------- diff --git a/doc/sphinx/yaml/reactions.rst b/doc/sphinx/yaml/reactions.rst index a0a9b396960..c3b5ee82d68 100644 --- a/doc/sphinx/yaml/reactions.rst +++ b/doc/sphinx/yaml/reactions.rst @@ -26,12 +26,17 @@ The fields common to all ``reaction`` entries are: - :ref:`chemically-activated ` - :ref:`pressure-dependent-Arrhenius ` - :ref:`Chebyshev ` + - :ref:`Blowers-Masel ` + - :ref:`surface-Blowers-Masel ` + + Reactions without a specified ``type`` on surfaces or edges are + automatically treated as :ref:`interface ` + reactions, and reactions that involve charge transfer between phases are + automatically treated as :ref:`electrochemical ` + reactions. Reactions on surfaces or edges specifying ``type`` as + ``Blowers-Masel`` are treated as + :ref:`surface-Blowers-Masel `. - Reactions on surfaces or edges are automatically treated as - :ref:`interface ` reactions, and reactions that - involve charge transfer between phases are automatically treated as - :ref:`electrochemical ` 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 `__. + +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 `__. +`described here `__. 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 `__. + +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}} diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index 1973323a7c4..7ef097ef32b 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -88,6 +88,7 @@ class GasKinetics : public BulkKinetics Rate1 m_plog_rates; Rate1 m_cheb_rates; + Rate1 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(); diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index fec3d801b59..fe74ed39480 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/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 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 m_blowers_masel_rates; + bool m_redo_rates; //! Vector of irreversible reaction numbers diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index a355f8a687c..620cb89c792 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/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; diff --git a/include/cantera/kinetics/RateCoeffMgr.h b/include/cantera/kinetics/RateCoeffMgr.h index 4f89f464a74..bc3ccbf4a07 100644 --- a/include/cantera/kinetics/RateCoeffMgr.h +++ b/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(); } diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 9598c3e1ff9..c709b5b2b3e 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/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 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 `` nodes in an XML document. //! //! The `` nodes are assumed to be children of the `` @@ -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 diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index e2a90951e94..5f526ebf220 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -95,6 +95,115 @@ class Arrhenius doublereal m_logA, m_b, m_E, m_A; }; +//! Blowers Masel reaction rate type depends on the enthalpy of reaction +/** + * The Blowers Masel approximation is written by Paul Blowers, + * Rich Masel (DOI: https://doi.org/10.1002/aic.690461015) to + * adjust the activation energy based on enthalpy change of a reaction: + * + * \f{eqnarray*}{ + * E_a &=& 0\; \text{if }\Delta H < -4E_0 \\ + * E_a &=& \Delta H\; \text{if }\Delta H > 4E_0 \\ + * E_a &=& \frac{(w + \Delta H / 2)(V_P - 2w + + * \Delta H)^2}{(V_P^2 - 4w^2 + (\Delta H)^2)}\; \text{Otherwise} + * \f} + * where + * \f[ + * V_P = \frac{2w (w + E_0)}{w - E_0}, + * \f] + * \f$ w \f$ is the average bond dissociation energy of the bond breaking + * and that being formed in the reaction. Since the expression is + * very insensitive to \f$ w \f$ for \f$ w >= 2 E_0 \f$, \f$ w \f$ + * can be approximated to an arbitrary high value like 1000 kJ/mol. + * + * After the activation energy is determined by Blowers-Masel approximation, + * it can be plugged into Arrhenius function to calculate the rate constant. + * \f[ + * k_f = A T^b \exp (-E_a/RT) + * \f] + */ +class BlowersMasel +{ +public: + //! Default constructor. + BlowersMasel(); + + /// Constructor. + /// @param A pre-exponential. The unit system is + /// (kmol, m, s). The actual units depend on the reaction + /// order and the dimensionality (surface or bulk). + /// @param b Temperature exponent. Non-dimensional. + /// @param E0 Intrinsic activation energy in temperature units. Kelvin. + /// @param w average bond dissociation energy of the bond being formed and + /// broken in the reaction, in temperature units. Kelvin. + + BlowersMasel(double A, double b, double E0, double w); + + void getParameters(AnyMap& rateNode, const Units& rate_units) const; + + //! Update concentration-dependent parts of the rate coefficient. + /*! + * For this class, there are no concentration-dependent parts, so this + * method does nothing. + */ + void update_C(const doublereal* c) { + } + + /** + * Update the value the rate constant. + * + * This function returns the actual value of the rate constant. It can be + * safely called for negative values of the pre-exponential factor. + */ + doublereal updateRC(doublereal logT, doublereal recipT, doublereal deltaH) { + double m_E = activationEnergy_R(deltaH); + return m_A * std::exp(m_b * logT - m_E * recipT); + } + + //! 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 actual activation energy (a function of the delta H of reaction) + //! divided by the gas constant (i.e. the activation temperature) [K] + doublereal activationEnergy_R(doublereal deltaH) { + double Ea; // will be in temperature units (Kelvin) + double deltaH_R = deltaH / GasConstant; // deltaH in temperature units (Kelvin) + if (deltaH_R < -4 * m_E0) { + Ea = 0; + } else if (deltaH_R > 4 * m_E0) { + Ea = deltaH_R; + } else { + // m_w is in Kelvin + // vp is in Kelvin + double vp = 2 * m_w * ((m_w + m_E0) / (m_w - m_E0)); + double vp_2w_dH = (vp - 2 * m_w + deltaH_R); // (Vp - 2 w + dH) + Ea = (m_w + deltaH_R / 2) * (vp_2w_dH * vp_2w_dH) / + (vp * vp - 4 * m_w * m_w + deltaH_R * deltaH_R); // in Kelvin + } + return Ea; + } + + //! Return the intrinsic activation energy divided by the gas constant (i.e. the + //! intrinsic activation temperature) [K] + doublereal activationEnergy_R0() const { + return m_E0; + } + //! Return the bond dissociation energy *w* divided by the gas constant[K] + doublereal bondEnergy() const { + return m_w; + } + +protected: + doublereal m_logA, m_b, m_A, m_w, m_E0; +}; /** * An Arrhenius rate with coverage-dependent terms. @@ -440,6 +549,144 @@ class ChebyshevRate vector_fp dotProd_; //!< dot product of chebCoeffs with the reduced pressure polynomial }; +/** + * A Blowers Masel rate with coverage-dependent terms. + * + * The rate expression is given by [Kee, R. J., Coltrin, M. E., & Glarborg, P. + * (2005). Chemically reacting flow: theory and practice. John Wiley & Sons. + * Eq 11.113]: + * \f[ + * k_f = A T^b \exp \left( + * \ln 10 \sum a_k \theta_k + * - \frac{1}{RT} \left( E_a + \sum E_k\theta_k \right) + * + \sum m_k \ln \theta_k + * \right) + * \f] + * or, equivalently, and as implemented in Cantera, + * \f[ + * k_f = A T^b \exp \left( - \frac{E_a}{RT} \right) + * \prod_k 10^{a_k \theta_k} \theta_k^{m_k} + * \exp \left( \frac{- E_k \theta_k}{RT} \right) + * \f] + * where the parameters \f$ (a_k, E_k, m_k) \f$ describe the dependency on the + * surface coverage of species \f$ k, \theta_k \f$. + * They are coupled with Blowers Masel approximation written by [Paul Blowers, + * Rich Masel(2004). Engineering Approximations For Activation Energies in Hydrogen + * Transfer Reactions. Eq (10)-(11)], which can be used to change the activation + * based on enthalpy of the reaction: + * + * \f{eqnarray*}{ + * E_a &=& 0\; \text{if }\Delta H < -4E_0 \\ + * E_a &=& \Delta H\; \text{if }\Delta H > 4E_0 \\ + * E_a &=& \frac{(w + \Delta H / 2)(V_P - 2w + + * \Delta H)^2}{(V_P^2 - 4w^2 + (\Delta H)^2)}\; \text{Otherwise} + * \f} + * where + * \f[ + * V_P = 2w (w + E_0) / (w - E_0), + * \f] + * \f$ w \f$ is the average bond dissociation energy of the bond breaking + * and that being formed in the reaction. Since the expression is + * very insensitive to \f$ w \f$ for \f$ w >= 2 E_0 \f$, \f$ w \f$ + * can be approximated to an arbitrary high value like 1000 kJ/mol. + */ +class BMSurfaceArrhenius +{ + +public: + BMSurfaceArrhenius(); + explicit BMSurfaceArrhenius(double A, double b, double Ta, double w); + + //! Add a coverage dependency for species *k*, with exponential dependence + //! *a*, power-law exponent *m*, and activation energy dependence *e*, + //! where *e* is in Kelvin, i.e. energy divided by the molar gas constant. + void addCoverageDependence(size_t k, double a, + double m, double e); + + void update_C(const doublereal* theta) { + m_acov = 0.0; + m_ecov = 0.0; + m_mcov = 0.0; + size_t k; + doublereal th; + for (size_t n = 0; n < m_ac.size(); n++) { + k = m_sp[n]; + m_acov += m_ac[n] * theta[k]; + m_ecov += m_ec[n] * theta[k]; + } + for (size_t n = 0; n < m_mc.size(); n++) { + k = m_msp[n]; + th = std::max(theta[k], Tiny); + m_mcov += m_mc[n]*std::log(th); + } + } + + /** + * Update the value of the rate constant. + * + * This function returns the actual value of the rate constant. It can be + * safely called for negative values of the pre-exponential factor. + */ + doublereal updateRC(doublereal logT, doublereal recipT, doublereal deltaH) { + double Ea = activationEnergy_R(deltaH); + return m_A * std::exp(std::log(10.0)*m_acov + m_b*logT - + (Ea)*recipT + m_mcov); + } + + //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending + //! on the reaction order) accounting coverage dependence. + /*! + * Returns reaction pre-exponent accounting for both *a* and *m*. + */ + doublereal preExponentialFactor() const { + return m_A * std::exp(std::log(10.0)*m_acov + m_mcov); + } + + //! Return effective temperature exponent + doublereal temperatureExponent() const { + return m_b; + } + + //! Return the activation energy divided by the gas constant (i.e. the + //! activation temperature) [K] based on the input enthalpy + //! accounting coverage dependence. + doublereal activationEnergy_R(doublereal deltaH) { + double deltaH_R = deltaH / GasConstant; // deltaH in temperature units (Kelvin) + double E; + if (deltaH_R < -4 * m_E0) { + E = 0; + } else if (deltaH_R > 4 * m_E0) { + E = deltaH_R; + } else { + // m_w is in Kelvin + // vp is in Kelvin + double vp = 2 * m_w * ((m_w + m_E0) / (m_w - m_E0)); + double vp_2w_dH = (vp - 2 * m_w + deltaH_R); // (Vp - 2 w + dH) + E = (m_w + deltaH_R / 2) * (vp_2w_dH * vp_2w_dH) / + (vp * vp - 4 * m_w * m_w + deltaH_R * deltaH_R); // in Kelvin + } + return E + m_ecov; + } + + //! Return the intrinsic activation energy divided by the gas constant (i.e. the + //! activation temperature) [K], accounting coverage dependence. + doublereal activationEnergy_R0() const { + return m_E0 + m_ecov; + } + + //! Return the bond energy *w* divided by the gas constant[K] + doublereal bondEnergy() const { + return m_w; + } + +protected: + doublereal m_b, m_A, m_E0, m_w; + doublereal m_acov, m_ecov, m_mcov; + std::vector m_sp, m_msp; + vector_fp m_ac, m_ec, m_mc; +}; + + } #endif diff --git a/include/cantera/kinetics/reaction_defs.h b/include/cantera/kinetics/reaction_defs.h index 7a1aedee97a..7e3050b6427 100644 --- a/include/cantera/kinetics/reaction_defs.h +++ b/include/cantera/kinetics/reaction_defs.h @@ -72,6 +72,11 @@ const int CUSTOMPY_RXN = 99; */ const int CHEMACT_RXN = 8; +/** + * A reaction to calculate the reaction rate with Blowers Masel approximation + */ +const int BLOWERSMASEL_RXN = 9; + /** * A reaction occurring on a surface. * NOTE: This is a bit ambiguous, and will be taken out in the future @@ -83,6 +88,11 @@ const int SURFACE_RXN = 20; //! A reaction occurring on an interface, e.g a surface or edge. const int INTERFACE_RXN = 20; +//! A reaction occurring on an interface, e.g a surface or edge. +//! The rate of the reaction can be calculated and updated by Blowers- +//! Masel approximation +const int BMINTERFACE_RXN =21; + /** * A reaction occurring at a one-dimensional interface between two surface phases. * NOTE: This is a bit ambiguous, and will be taken out in the future diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index cb9dce619fd..c2ddaaf4dd3 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -451,6 +451,21 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": CxxTestReaction() cbool allow_negative_pre_exponential_factor + cdef cppclass CxxBlowersMasel "Cantera::BlowersMasel": + CxxBlowersMasel() + CxxBlowersMasel(double, double, double, double) + double updateRC(double, double, double) + double preExponentialFactor() + double temperatureExponent() + double activationEnergy_R(double) + double activationEnergy_R0() + double bondEnergy() + + cdef cppclass CxxBlowersMaselReaction "Cantera::BlowersMaselReaction"(CxxReaction): + CxxBlowersMaselReaction() + CxxBlowersMasel rate + cbool allow_negative_pre_exponential_factor + cdef cppclass CxxCoverageDependency "Cantera::CoverageDependency": CxxCoverageDependency(double, double, double) double a @@ -463,6 +478,12 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cbool use_motz_wise_correction string sticking_species + cdef cppclass CxxBlowersMaselInterfaceReaction "Cantera::BlowersMaselInterfaceReaction" (CxxBlowersMaselReaction): + stdmap[string, CxxCoverageDependency] coverage_deps + cbool is_sticking_coefficient + cbool use_motz_wise_correction + string sticking_species + cdef extern from "cantera/kinetics/FalloffFactory.h" namespace "Cantera": cdef shared_ptr[CxxFalloff] CxxNewFalloff "Cantera::newFalloff" (string, vector[double]) except +translate_exception @@ -1101,6 +1122,10 @@ cdef class Arrhenius: cdef CxxArrhenius* rate cdef Reaction reaction # parent reaction, to prevent garbage collection +cdef class BlowersMasel: + cdef CxxBlowersMasel* rate + cdef Reaction reaction + cdef class Falloff: cdef shared_ptr[CxxFalloff] _falloff cdef CxxFalloff* falloff diff --git a/interfaces/cython/cantera/ck2yaml.py b/interfaces/cython/cantera/ck2yaml.py index 2808dfa3f6a..8173c0ac884 100644 --- a/interfaces/cython/cantera/ck2yaml.py +++ b/interfaces/cython/cantera/ck2yaml.py @@ -1898,7 +1898,7 @@ def write_yaml(self, name='gas', out_name='mech.yaml'): metadata = BlockMap([ ('generator', 'ck2yaml'), ('input-files', FlowList(files)), - ('cantera-version', '2.6.0a1'), + ('cantera-version', '2.6.0a2'), ('date', formatdate(localtime=True)), ]) if desc.strip(): diff --git a/interfaces/cython/cantera/cti2yaml.py b/interfaces/cython/cantera/cti2yaml.py index 447f90c71db..426fdb10657 100644 --- a/interfaces/cython/cantera/cti2yaml.py +++ b/interfaces/cython/cantera/cti2yaml.py @@ -1622,7 +1622,7 @@ def convert(filename=None, output_name=None, text=None): # information regarding conversion metadata = BlockMap([ ('generator', 'cti2yaml'), - ('cantera-version', '2.6.0a1'), + ('cantera-version', '2.6.0a2'), ('date', formatdate(localtime=True)), ]) if filename is not None: diff --git a/interfaces/cython/cantera/ctml2yaml.py b/interfaces/cython/cantera/ctml2yaml.py index 89991fb8689..955df98c807 100644 --- a/interfaces/cython/cantera/ctml2yaml.py +++ b/interfaces/cython/cantera/ctml2yaml.py @@ -2616,7 +2616,7 @@ def convert( metadata = BlockMap( { "generator": "ctml2yaml", - "cantera-version": "2.6.0a1", + "cantera-version": "2.6.0a2", "date": formatdate(localtime=True), } ) diff --git a/interfaces/cython/cantera/examples/kinetics/blowers_masel.py b/interfaces/cython/cantera/examples/kinetics/blowers_masel.py new file mode 100644 index 00000000000..90db4174c69 --- /dev/null +++ b/interfaces/cython/cantera/examples/kinetics/blowers_masel.py @@ -0,0 +1,114 @@ +""" +A simple example to demonstrate the difference between Blowers-Masel +reaction and elementary reaction. + +The first two reactions have the same reaction equations with Arrhenius and +Blowers-Masel rate parameters, respectively. The Blowers-Masel parameters are the same +as the Arrhenius parameters with an additional value, bond energy. + +First we show that the forward rate constants of the first 2 different reactions are +different because of the different rate expression, then we print the forward rate +constants for reaction 2 and reaction 3 to show that even two reactions that have the +same Blowers-Masel parameters can have different forward rate constants. + +The first plot generated shows how the rate constant changes with respect to temperature +for elementary and Blower-Masel reactions. The second plot shows the activation energy +change of a Blowers-Masel reaction with respect to the delta enthalpy of the reaction. + +Requires: cantera >= 2.6.0, matplotlib >= 2.0 +""" + +import cantera as ct +import numpy as np +import matplotlib.pyplot as plt + +#Create an elementary reaction O+H2<=>H+OH +r1 = ct.ElementaryReaction({'O':1, 'H2':1}, {'H':1, 'OH':1}) +r1.rate = ct.Arrhenius(3.87e1, 2.7, 6260*1000*4.184) + +#Create a Blowers-Masel reaction O+H2<=>H+OH +r2 = ct.BlowersMaselReaction({'O':1, 'H2':1}, {'H':1, 'OH':1}) +r2.rate = ct.BlowersMasel(3.87e1, 2.7, 6260*1000*4.184, 1e9) + +#Create a Blowers-Masel reaction with same parameters with r2 +#reaction equation is H+CH4<=>CH3+H2 +r3 = ct.BlowersMaselReaction({'H':1, 'CH4':1}, {'CH3':1, 'H2':1}) +r3.rate = ct.BlowersMasel(3.87e1, 2.7, 6260*1000*4.184, 1e9) + +gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=ct.Solution('gri30.yaml').species(), reactions=[r1, r2, r3]) + +gas.TP = 300, ct.one_atm + +r1_rc = gas.forward_rate_constants[0] +r2_rc = gas.forward_rate_constants[1] +r3_rc = gas.forward_rate_constants[2] + +print("The first and second reactions have same reaction equation," + " but they have different reaction types, so the forward rate" + " constant of the first reaction is {0:.3f} kmol/(m^3.s)," + " the forward rate constant of the second reaction is {1:.3f} kmol/(m^3.s).".format(r1_rc, r2_rc)) + +print("The rate parameters of second and the third reactions are same," + " but the forward rate constant of second reaction is {0:.3f} kmol/(m^3.s)," + " the forward rate constant of the third reaction is" + " {1:.3f} kmol/(m^3.s).".format(r2_rc, r3_rc)) + +# Comparing the reaction forward rate constant change of +# Blowers-Masel reaction and elementary reaction with +# respect to the temperature. +r1_kf = [] +r2_kf = [] +T_range = np.arange(300, 3500, 100) +for temp in T_range: + gas.TP = temp, ct.one_atm + r1_kf.append(gas.forward_rate_constants[0]) + r2_kf.append(gas.forward_rate_constants[1]) +plt.plot(T_range, r1_kf, label='Reaction 1 (Elementary)') +plt.plot(T_range, r2_kf, label='Reaction 2 (Blowers-Masel)') +plt.xlabel("Temperature(K)") +plt.ylabel(r"Forward Rate Constant (kmol/(m$^3\cdot$ s))") +plt.title("Comparison of $k_f$ vs. Temperature For Reaction 1 and 2",y=1.1) +plt.legend() +plt.savefig("kf_to_T.png") + +# This is the function to change the enthalpy of a species +# so that the enthalpy change of reactions involving this +#species can be changed +def change_species_enthalpy(gas, species_name, dH): + """ + Find the species by name and change it's enthalpy by dH (in J/kmol) + """ + index = gas.species_index(species_name) + + species = gas.species(index) + # dx is in fact (delta H / R). Note that R in cantera is 8314.462 J/kmol/K + dx = dH / ct.gas_constant + perturbed_coeffs = species.thermo.coeffs.copy() + perturbed_coeffs[6] += dx + perturbed_coeffs[13] += dx + + species.thermo = ct.NasaPoly2(species.thermo.min_temp, species.thermo.max_temp, + species.thermo.reference_pressure, perturbed_coeffs) + + gas.modify_species(index, species) + +# Plot the activation energy change of reaction 2 with respect to the +# enthalpy change +E0 = gas.reaction(1).rate.intrinsic_activation_energy +upper_limit_enthalpy = 5 * E0 +lower_limit_enthalpy = -5 * E0 + +Ea_list = [] +deltaH_list = np.linspace(lower_limit_enthalpy, upper_limit_enthalpy, 100) +for deltaH in deltaH_list: + change_species_enthalpy(gas, 'H', deltaH - gas.delta_enthalpy[1]) + Ea_list.append(gas.reaction(1).rate.activation_energy(gas.delta_enthalpy[1])) + +plt.figure() +plt.plot(deltaH_list, Ea_list) +plt.xlabel("Enthalpy Change (J/kmol)") +plt.ylabel("Activation Energy Change (J/kmol)") +plt.title(r"$E_a$ vs. $\Delta H$ For O+H2<=>H+OH", y=1.1) +plt.savefig("Ea_to_H.png") +plt.show() diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index eb0135f1c68..3db13faf693 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -910,6 +910,103 @@ cdef class ChebyshevReaction(Reaction): r.rate.update_C(&logP) return r.rate.updateRC(logT, recipT) +cdef class BlowersMasel: + """ + A reaction rate coefficient which depends on temperature and enthalpy change + of the reaction follows Blowers-Masel approximation and modified Arrhenius form + described in `Arrhenius`. The functions are used by reactions defined using + `BlowersMaselReaction` and `BlowersMaselInterfaceReaction`. + """ + def __cinit__(self, A=0, b=0, E0=0, w=0, init=True): + if init: + self.rate = new CxxBlowersMasel(A, b, E0 / gas_constant, w / gas_constant) + self.reaction = None + + def __dealloc__(self): + if self.reaction is None: + del self.rate + + property pre_exponential_factor: + """ + The pre-exponential factor *A* in units of m, kmol, and s raised to + powers depending on the reaction order. + """ + def __get__(self): + return self.rate.preExponentialFactor() + + property temperature_exponent: + """ + The temperature exponent *b*. + """ + def __get__(self): + return self.rate.temperatureExponent() + + def activation_energy(self, float dH): + """ + The activation energy *E* [J/kmol] + + :param dH: The enthalpy of reaction [J/kmol] at the current temperature + """ + return self.rate.activationEnergy_R(dH) * gas_constant + + property bond_energy: + """ + Average bond dissociation energy of the bond being formed and broken + in the reaction *E* [J/kmol]. + """ + def __get__(self): + return self.rate.bondEnergy() * gas_constant + + property intrinsic_activation_energy: + """ + The intrinsic activation energy *E0* [J/kmol]. + """ + def __get__(self): + return self.rate.activationEnergy_R0() * gas_constant + + def __repr__(self): + return 'BlowersMasel(A={:g}, b={:g}, E0={:g}, w={:g})'.format( + self.pre_exponential_factor, self.temperature_exponent, + self.intrinsic_activation_energy, self.bond_energy) + + def __call__(self, float T, float dH): + cdef double logT = np.log(T) + cdef double recipT = 1/T + return self.rate.updateRC(logT, recipT, dH) + +cdef wrapBlowersMasel(CxxBlowersMasel* rate, Reaction reaction): + r = BlowersMasel(init=False) + r.rate = rate + r.reaction = reaction + return r + +cdef class BlowersMaselReaction(Reaction): + """ + A reaction which follows mass-action kinetics with Blowers Masel + reaction rate. + """ + reaction_type = "Blowers-Masel" + + property rate: + """ Get/Set the `Arrhenius` rate coefficient for this reaction. """ + def __get__(self): + cdef CxxBlowersMaselReaction* r = self.reaction + return wrapBlowersMasel(&(r.rate), self) + def __set__(self, BlowersMasel rate): + cdef CxxBlowersMaselReaction* r = self.reaction + r.rate = deref(rate.rate) + + property allow_negative_pre_exponential_factor: + """ + Get/Set whether the rate coefficient is allowed to have a negative + pre-exponential factor. + """ + def __get__(self): + cdef CxxBlowersMaselReaction* r = self.reaction + return r.allow_negative_pre_exponential_factor + def __set__(self, allow): + cdef CxxBlowersMaselReaction* r = self.reaction + r.allow_negative_pre_exponential_factor = allow cdef class CustomReaction(Reaction): """ @@ -1076,3 +1173,75 @@ cdef class InterfaceReaction(ElementaryReaction): def __set__(self, species): cdef CxxInterfaceReaction* r = self.reaction r.sticking_species = stringify(species) + +cdef class BlowersMaselInterfaceReaction(BlowersMaselReaction): + """ + A reaction occurring on an `Interface` (i.e. a surface or an edge) + with the rate parameterization of `BlowersMasel`. + """ + reaction_type = "surface-Blowers-Masel" + + property coverage_deps: + """ + Get/Set a dict containing adjustments to the Arrhenius rate expression + dependent on surface species coverages. The keys of the dict are species + names, and the values are tuples specifying the three coverage + parameters ``(a, m, E)`` which are the modifiers for the pre-exponential + factor [m, kmol, s units], the temperature exponent [nondimensional], + and the activation energy [J/kmol], respectively. + """ + def __get__(self): + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction + deps = {} + cdef pair[string,CxxCoverageDependency] item + for item in r.coverage_deps: + deps[pystr(item.first)] = (item.second.a, item.second.m, + item.second.E * gas_constant) + return deps + def __set__(self, deps): + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction + r.coverage_deps.clear() + cdef str species + for species, D in deps.items(): + r.coverage_deps[stringify(species)] = CxxCoverageDependency( + D[0], D[2] / gas_constant, D[1]) + + property is_sticking_coefficient: + """ + Get/Set a boolean indicating if the rate coefficient for this reaction + is expressed as a sticking coefficient rather than the forward rate + constant. + """ + def __get__(self): + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction + return r.is_sticking_coefficient + def __set__(self, stick): + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction + r.is_sticking_coefficient = stick + + property use_motz_wise_correction: + """ + Get/Set a boolean indicating whether to use the correction factor + developed by Motz & Wise for reactions with high (near-unity) sticking + coefficients when converting the sticking coefficient to a rate + coefficient. + """ + def __get__(self): + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction + return r.use_motz_wise_correction + def __set__(self, mw): + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction + r.use_motz_wise_correction = mw + + property sticking_species: + """ + The name of the sticking species. Needed only for reactions with + multiple non-surface reactant species, where the sticking species is + ambiguous. + """ + def __get__(self): + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction + return pystr(r.sticking_species) + def __set__(self, species): + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction + r.sticking_species = stringify(species) diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index b3c6b7ecf1b..8c37c2146b9 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -1082,6 +1082,81 @@ def test_chebyshev_bad_shape_yaml(self): - [0.3177, 0.26889, 0.094806, -7.6385e-03] - [-0.031285, -0.039412, 0.044375, 0.014458]''', gas) + def test_BlowersMasel(self): + r = ct.BlowersMaselReaction({'O':1, 'H2':1}, {'H':1, 'OH':1}) + r.rate = ct.BlowersMasel(3.87e1, 2.7, 6260*1000*4.184, 1e9*1000*4.184) + + gas1 = ct.Solution('BM_test.yaml') + + gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=gas1.species(), reactions=[r]) + + gas1.TP = self.gas.TP + gas2.TP = self.gas.TP + gas1.X = 'H2:0.1, H2O:0.2, O2:0.7, O:1e-4, OH:1e-5, H:2e-5' + gas2.X = 'H2:0.1, H2O:0.2, O2:0.7, O:1e-4, OH:1e-5, H:2e-5' + + self.assertNear(gas2.forward_rate_constants[0], + gas1.forward_rate_constants[0], rtol=1e-7) + self.assertNear(gas2.net_rates_of_progress[0], + gas1.net_rates_of_progress[0], rtol=1e-7) + + def test_Blowers_Masel_rate(self): + gas = ct.Solution('BM_test.yaml') + R = gas.reaction(0) + self.assertNear(R.rate(gas.T, gas.delta_enthalpy[0]), + gas.forward_rate_constants[0], rtol=1e-7) + + def test_negative_A_blowersmasel(self): + species = ct.Solution('BM_test.yaml').species() + r = ct.BlowersMaselReaction({'O':1, 'H2':1}, {'H':1, 'OH':1}) + r.rate = ct.BlowersMasel(-3.87e1, 2.7, 6260*1000*4.184, 1e9) + + self.assertFalse(r.allow_negative_pre_exponential_factor) + + with self.assertRaisesRegex(ct.CanteraError, 'negative pre-exponential'): + gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=species, reactions=[r]) + + r.allow_negative_pre_exponential_factor = True + gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=species, reactions=[r]) + + def test_Blowers_Masel_change_enthalpy(self): + gas = ct.Solution('BM_test.yaml') + r = gas.reaction(0) + E0 = r.rate.intrinsic_activation_energy + w = r.rate.bond_energy + A = r.rate.pre_exponential_factor + b = r.rate.temperature_exponent + vp = 2 * w * (w+E0) / (w - E0) + deltaH = gas.delta_enthalpy[0] + E = ((w + deltaH / 2) * (vp - 2 * w + deltaH) ** 2 / + (vp ** 2 - 4 * w ** 2 + deltaH ** 2)) + self.assertNear(E, gas.reaction(0).rate.activation_energy(deltaH)) + + deltaH_high = 10 * gas.reaction(0).rate.intrinsic_activation_energy + deltaH_low = -20 * gas.reaction(0).rate.intrinsic_activation_energy + index = gas.species_index('OH') + species = gas.species('OH') + perturbed_coeffs = species.thermo.coeffs.copy() + perturbed_coeffs[6] += deltaH_high / ct.gas_constant + perturbed_coeffs[13] += deltaH_high / ct.gas_constant + species.thermo = ct.NasaPoly2(species.thermo.min_temp, species.thermo.max_temp, + species.thermo.reference_pressure, perturbed_coeffs) + gas.modify_species(index, species) + self.assertEqual(deltaH_high, gas.reaction(0).rate.activation_energy(deltaH_high)) + self.assertNear(A*gas.T**b*np.exp(-deltaH_high/ct.gas_constant/gas.T), gas.forward_rate_constants[0]) + + perturbed_coeffs = species.thermo.coeffs.copy() + perturbed_coeffs[6] += deltaH_low / ct.gas_constant + perturbed_coeffs[13] += deltaH_low / ct.gas_constant + species.thermo = ct.NasaPoly2(species.thermo.min_temp, species.thermo.max_temp, + species.thermo.reference_pressure, perturbed_coeffs) + gas.modify_species(index, species) + self.assertEqual(0, gas.reaction(0).rate.activation_energy(deltaH_low)) + self.assertNear(A*gas.T**b*np.exp(0/ct.gas_constant/gas.T), gas.forward_rate_constants[0]) + def test_interface(self): surf_species = ct.Species.listFromFile('ptcombust.xml') gas = ct.Solution('ptcombust.xml', 'gas') @@ -1108,6 +1183,35 @@ def test_interface(self): self.assertNear(surf1.net_rates_of_progress[1], surf2.net_rates_of_progress[0]) + def test_BlowersMaselinterface(self): + gas = ct.Solution('gri30.yaml') + gas.TPX = 300, ct.one_atm, {"CH4": 0.095, "O2": 0.21, "AR": 0.79} + surf1 = ct.Interface('BM_test.yaml', 'Pt_surf', [gas]) + r1 = ct.BlowersMaselInterfaceReaction() + r1.reactants = 'H(S):2' + r1.products = 'H2:1, PT(S):2' + r1.rate = ct.BlowersMasel(3.7e20, 0, 67.4e6, 1e9) + r1.coverage_deps = {'H(S)': (0, 0, -6e6)} + + self.assertNear(r1.coverage_deps['H(S)'][2], -6e6) + + surf_species = [] + for species in surf1.species(): + surf_species.append(species) + surf2 = ct.Interface(thermo='Surface', species=surf_species, + kinetics='interface', reactions=[r1], adjacent=[gas]) + + surf2.site_density = surf1.site_density + surf1.coverages = surf2.coverages = 'PT(S):0.7, H(S):0.3' + gas.TP = surf2.TP = surf1.TP + + for T in [300, 500, 1500]: + gas.TP = surf1.TP = surf2.TP = T, 5*ct.one_atm + self.assertNear(surf1.forward_rate_constants[0], + surf2.forward_rate_constants[0]) + self.assertNear(surf1.net_rates_of_progress[0], + surf2.net_rates_of_progress[0]) + def test_modify_invalid(self): # different reaction type tbr = self.gas.reaction(0) @@ -1214,6 +1318,28 @@ def test_modify_chebyshev(self): kf = gas.forward_rate_constants self.assertNear(kf[4], kf[5]) + def test_modify_BlowersMasel(self): + gas = ct.Solution('BM_test.yaml') + gas.X = 'H2:0.1, H2O:0.2, O2:0.7, O:1e-4, OH:1e-5, H:2e-5' + gas.TP = self.gas.TP + R = gas.reaction(0) + A1 = R.rate.pre_exponential_factor + b1 = R.rate.temperature_exponent + Ta1 = R.rate.activation_energy(gas.delta_enthalpy[0]) / ct.gas_constant + T = gas.T + self.assertNear(A1* T ** b1 * np.exp(- Ta1 / T), + gas.forward_rate_constants[0]) + + # randomly modify the rate parameters of a Blowers-Masel reaction + A2 = 1.5 * A1 + b2 = b1 + 0.1 + Ta_intrinsic = R.rate.intrinsic_activation_energy * 1.2 + w = R.rate.bond_energy * 0.8 + R.rate = ct.BlowersMasel(A2, b2, Ta_intrinsic, w) + Ta2 = R.rate.activation_energy(gas.delta_enthalpy[0]) / ct.gas_constant + gas.modify_reaction(0, R) + self.assertNear(A2*T**b2*np.exp(-Ta2/T), gas.forward_rate_constants[0]) + def test_modify_interface(self): gas = ct.Solution('ptcombust.xml', 'gas') surf = ct.Interface('ptcombust.xml', 'Pt_surf', [gas]) @@ -1272,3 +1398,70 @@ def test_motz_wise(self): # M&W toggled on (locally) for reaction 9 self.assertNear(2.0 * k1[9], k2[9]) # sticking coefficient = 1.0 + + def test_modify_BMinterface(self): + gas = ct.Solution('gri30.yaml') + gas.TPX = 300, ct.one_atm, {"CH4": 0.095, "O2": 0.21, "AR": 0.79} + surf = ct.Interface('BM_test.yaml', 'Pt_surf', [gas]) + surf.coverages = 'O(S):0.1, PT(S):0.5, H(S):0.4' + gas.TP = surf.TP + + R = surf.reaction(0) + R.coverage_deps = {'O(S)': (0.0, 0.0, -3e6)} + surf.modify_reaction(0, R) + + # Rate constant should now be independent of H(S) coverage, but + # dependent on O(S) coverage + Ek = surf.reaction(0).coverage_deps['O(S)'][-1] + k1 = surf.forward_rate_constants[0] + O2_theta_k1 = surf.coverages[-1] + surf.coverages = 'O(S):0.2, PT(S):0.4, H(S):0.4' + k2 = surf.forward_rate_constants[0] + O2_theta_k2 = surf.coverages[-1] + O2_delta_theta_k = O2_theta_k1 - O2_theta_k2 + surf.coverages = 'O(S):0.2, PT(S):0.6, H(S):0.2' + k3 = surf.forward_rate_constants[0] + + self.assertNear(k1 / k2, np.exp(-O2_delta_theta_k * Ek / ct.gas_constant / surf.T)) + self.assertNear(k2, k3) + + def test_modify_BMsticking(self): + gas = ct.Solution('gri30.yaml') + gas.TPX = 300, ct.one_atm, {"CH4": 0.095, "O2": 0.21, "AR": 0.79} + surf = ct.Interface('BM_test.yaml', 'Pt_surf', [gas]) + surf.coverages = 'O(S):0.1, PT(S):0.5, H(S):0.4' + gas.TP = surf.TP + + R = surf.reaction(1) + R.rate = ct.BlowersMasel(0.25, 0, 0, 1000000) # original sticking coefficient = 1.0 + + k1 = surf.forward_rate_constants[1] + surf.modify_reaction(1, R) + k2 = surf.forward_rate_constants[1] + self.assertNear(k1, 4*k2) + + def test_BMmotz_wise(self): + # Motz & Wise off for all reactions + gas1 = ct.Solution('gri30.yaml') + gas1.TPX = 300, ct.one_atm, {"CH4": 0.095, "O2": 0.21, "AR": 0.79} + surf1 = ct.Interface('BM_test.yaml', 'Pt_surf', [gas1]) + surf1.coverages = 'O(S):0.1, PT(S):0.5, H(S):0.4' + gas1.TP = surf1.TP + + # Motz & Wise correction on for some reactions + gas2 = ct.Solution('BM-ptcombust-Motz-Wise.yaml', 'gas') + surf2 = ct.Interface('BM-ptcombust-Motz-Wise.yaml', 'Pt_surf', [gas2]) + surf2.TPY = surf1.TPY + + k1 = surf1.forward_rate_constants + k2 = surf2.forward_rate_constants + + # M&W toggled on (globally) for reactions 2 and 7 + self.assertNear(2.0 * k1[1], k2[1]) # sticking coefficient = 1.0 + self.assertNear(1.6 * k1[2], k2[2]) # sticking coefficient = 0.75 + + # M&W toggled off (locally) for reaction 4 + self.assertNear(k1[3], k2[3]) + + # M&W toggled on (locally) for reaction 5 + self.assertNear(2.0 * k1[4], k2[4]) # sticking coefficient = 1.0 diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index ef008ab112f..c664215e7ff 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -42,6 +42,11 @@ void GasKinetics::update_rates_T() } updateKc(); m_ROP_ok = false; + if (m_blowersmasel_rates.nReactions()) { + thermo().getPartialMolarEnthalpies(m_grt.data()); + getReactionDelta(m_grt.data(), m_dH.data()); + m_blowersmasel_rates.updateBlowersMasel(T, logT, m_rfn.data(), m_dH.data()); + } } if (T != m_temp || P != m_pres) { @@ -253,6 +258,8 @@ bool GasKinetics::addReaction(shared_ptr r) addPlogReaction(dynamic_cast(*r)); } else if (r->type() == "Chebyshev") { addChebyshevReaction(dynamic_cast(*r)); + } else if (r->type() == "Blowers-Masel") { + addBlowersMaselReaction(dynamic_cast(*r)); } else { throw CanteraError("GasKinetics::addReaction", "Unknown reaction type specified: '{}'", r->type()); @@ -324,6 +331,11 @@ void GasKinetics::addChebyshevReaction(ChebyshevReaction& r) m_cheb_rates.install(nReactions()-1, r.rate); } +void GasKinetics::addBlowersMaselReaction(BlowersMaselReaction& r) +{ + m_blowersmasel_rates.install(nReactions()-1, r.rate); +} + void GasKinetics::modifyReaction(size_t i, shared_ptr rNew) { // operations common to all bulk reaction types @@ -346,6 +358,8 @@ void GasKinetics::modifyReaction(size_t i, shared_ptr rNew) modifyPlogReaction(i, dynamic_cast(*rNew)); } else if (rNew->type() == "Chebyshev") { modifyChebyshevReaction(i, dynamic_cast(*rNew)); + } else if (rNew->type() == "Blowers-Masel") { + modifyBlowersMaselReaction(i, dynamic_cast(*rNew)); } else { throw CanteraError("GasKinetics::modifyReaction", "Unknown reaction type specified: '{}'", rNew->type()); @@ -380,6 +394,11 @@ void GasKinetics::modifyChebyshevReaction(size_t i, ChebyshevReaction& r) m_cheb_rates.replace(i, r.rate); } +void GasKinetics::modifyBlowersMaselReaction(size_t i, BlowersMaselReaction& r) +{ + m_blowersmasel_rates.replace(i, r.rate); +} + void GasKinetics::init() { BulkKinetics::init(); diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index a1834b7efdb..6946a1a8b84 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -54,6 +54,7 @@ void InterfaceKinetics::_update_rates_T() if (m_has_coverage_dependence) { m_surf->getCoverages(m_actConc.data()); m_rates.update_C(m_actConc.data()); + m_blowers_masel_rates.update_C(m_actConc.data()); m_redo_rates = true; } @@ -65,6 +66,13 @@ void InterfaceKinetics::_update_rates_T() // Calculate the forward rate constant by calling m_rates and store it in m_rfn[] m_rates.update(T, m_logtemp, m_rfn.data()); + for (size_t n = 0; n < nPhases(); n++) { + thermo(n).getPartialMolarEnthalpies(m_grt.data() + m_start[n]); + } + + // Use the stoichiometric manager to find deltaH for each reaction. + getReactionDelta(m_grt.data(), m_dH.data()); + m_blowers_masel_rates.updateBlowersMasel(T, m_logtemp, m_rfn.data(), m_dH.data()); applyStickingCorrection(T, m_rfn.data()); // If we need to do conversions between exchange current density @@ -501,49 +509,76 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base) if (!added) { return false; } + if (r_base->reaction_type == BMINTERFACE_RXN) { + BlowersMaselInterfaceReaction& r = dynamic_cast(*r_base); + BMSurfaceArrhenius rate = buildBMSurfaceArrhenius(i, r, false); + m_blowers_masel_rates.install(i, rate); - InterfaceReaction& r = dynamic_cast(*r_base); - SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, false); - m_rates.install(i, rate); + // Turn on the global flag indicating surface coverage dependence + if (!r.coverage_deps.empty()) { + m_has_coverage_dependence = true; + } + if (r.reversible) { + m_revindex.push_back(i); + } else { + m_irrev.push_back(i); + } - // Turn on the global flag indicating surface coverage dependence - if (!r.coverage_deps.empty()) { - m_has_coverage_dependence = true; - } + m_rxnPhaseIsReactant.emplace_back(nPhases(), false); + m_rxnPhaseIsProduct.emplace_back(nPhases(), false); - ElectrochemicalReaction* re = dynamic_cast(&r); - if (re) { - m_has_electrochem_rxns = true; - m_beta.push_back(re->beta); - m_ctrxn.push_back(i); - 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); + for (const auto& sp : r.reactants) { + size_t k = kineticsSpeciesIndex(sp.first); + size_t p = speciesPhaseIndex(k); + m_rxnPhaseIsReactant[i][p] = true; + } + for (const auto& sp : r.products) { + size_t k = kineticsSpeciesIndex(sp.first); + size_t p = speciesPhaseIndex(k); + m_rxnPhaseIsProduct[i][p] = true; } - } - if (r.reversible) { - m_revindex.push_back(i); } else { - m_irrev.push_back(i); - } + InterfaceReaction& r = dynamic_cast(*r_base); + SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, false); + m_rates.install(i, rate); + + // Turn on the global flag indicating surface coverage dependence + if (!r.coverage_deps.empty()) { + m_has_coverage_dependence = true; + } + ElectrochemicalReaction* re = dynamic_cast(&r); + if (re) { + m_has_electrochem_rxns = true; + m_beta.push_back(re->beta); + m_ctrxn.push_back(i); + 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); + } + } + if (r.reversible) { + m_revindex.push_back(i); + } else { + m_irrev.push_back(i); + } - m_rxnPhaseIsReactant.emplace_back(nPhases(), false); - m_rxnPhaseIsProduct.emplace_back(nPhases(), false); + m_rxnPhaseIsReactant.emplace_back(nPhases(), false); + m_rxnPhaseIsProduct.emplace_back(nPhases(), false); - for (const auto& sp : r.reactants) { - size_t k = kineticsSpeciesIndex(sp.first); - size_t p = speciesPhaseIndex(k); - m_rxnPhaseIsReactant[i][p] = true; - } - for (const auto& sp : r.products) { - size_t k = kineticsSpeciesIndex(sp.first); - size_t p = speciesPhaseIndex(k); - m_rxnPhaseIsProduct[i][p] = true; + for (const auto& sp : r.reactants) { + size_t k = kineticsSpeciesIndex(sp.first); + size_t p = speciesPhaseIndex(k); + m_rxnPhaseIsReactant[i][p] = true; + } + for (const auto& sp : r.products) { + size_t k = kineticsSpeciesIndex(sp.first); + size_t p = speciesPhaseIndex(k); + m_rxnPhaseIsProduct[i][p] = true; + } } - deltaElectricEnergy_.push_back(0.0); m_deltaG0.push_back(0.0); m_deltaG.push_back(0.0); @@ -555,10 +590,15 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base) void InterfaceKinetics::modifyReaction(size_t i, shared_ptr r_base) { Kinetics::modifyReaction(i, r_base); - InterfaceReaction& r = dynamic_cast(*r_base); - SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, true); - m_rates.replace(i, rate); - + if (r_base->reaction_type == BMINTERFACE_RXN) { + BlowersMaselInterfaceReaction& r = dynamic_cast(*r_base); + BMSurfaceArrhenius rate = buildBMSurfaceArrhenius(i, r, true); + m_blowers_masel_rates.replace(i, rate); + } else { + InterfaceReaction& r = dynamic_cast(*r_base); + SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, true); + m_rates.replace(i, rate); + } // Invalidate cached data m_redo_rates = true; m_temp += 0.1; @@ -656,6 +696,99 @@ SurfaceArrhenius InterfaceKinetics::buildSurfaceArrhenius( return rate; } +BMSurfaceArrhenius InterfaceKinetics::buildBMSurfaceArrhenius( + size_t i, BlowersMaselInterfaceReaction& r, bool replace) +{ + if (r.is_sticking_coefficient) { + // Identify the interface phase + size_t iInterface = npos; + size_t min_dim = 4; + for (size_t n = 0; n < nPhases(); n++) { + if (thermo(n).nDim() < min_dim) { + iInterface = n; + min_dim = thermo(n).nDim(); + } + } + + std::string sticking_species = r.sticking_species; + if (sticking_species == "") { + // Identify the sticking species if not explicitly given + bool foundStick = false; + for (const auto& sp : r.reactants) { + size_t iPhase = speciesPhaseIndex(kineticsSpeciesIndex(sp.first)); + if (iPhase != iInterface) { + // Non-interface species. There should be exactly one of these + if (foundStick) { + throw CanteraError("InterfaceKinetics::buildBMSurfaceArrhenius", + "Multiple non-interface species found" + "in sticking reaction: '" + r.equation() + "'"); + } + foundStick = true; + sticking_species = sp.first; + } + } + if (!foundStick) { + throw CanteraError("InterfaceKinetics::buildBMSurfaceArrhenius", + "No non-interface species found" + "in sticking reaction: '" + r.equation() + "'"); + } + } + + double surface_order = 0.0; + double multiplier = 1.0; + // Adjust the A-factor + for (const auto& sp : r.reactants) { + size_t iPhase = speciesPhaseIndex(kineticsSpeciesIndex(sp.first)); + const ThermoPhase& p = thermo(iPhase); + size_t k = p.speciesIndex(sp.first); + if (sp.first == sticking_species) { + multiplier *= sqrt(GasConstant/(2*Pi*p.molecularWeight(k))); + } else { + // Non-sticking species. Convert from coverages used in the + // sticking probability expression to the concentration units + // used in the mass action rate expression. For surface phases, + // the dependence on the site density is incorporated when the + // rate constant is evaluated, since we don't assume that the + // site density is known at this time. + double order = getValue(r.orders, sp.first, sp.second); + if (&p == m_surf) { + multiplier *= pow(m_surf->size(k), order); + surface_order += order; + } else { + multiplier *= pow(p.standardConcentration(k), -order); + } + } + } + + if (!replace) { + m_stickingData.emplace_back(StickData{i, surface_order, multiplier, + r.use_motz_wise_correction}); + } else { + // Modifying an existing sticking reaction. + for (auto& item : m_stickingData) { + if (item.index == i) { + item.order = surface_order; + item.multiplier = multiplier; + item.use_motz_wise = r.use_motz_wise_correction; + break; + } + } + } + } + + BMSurfaceArrhenius rate(r.rate.preExponentialFactor(), + r.rate.temperatureExponent(), + r.rate.activationEnergy_R0(), + r.rate.bondEnergy()); + + // Set up coverage dependencies + for (const auto& sp : r.coverage_deps) { + size_t k = thermo(reactionPhaseIndex()).speciesIndex(sp.first); + rate.addCoverageDependence(k, sp.second.a, sp.second.m, sp.second.E); + } + return rate; +} + void InterfaceKinetics::setIOFlag(int ioFlag) { m_ioFlag = ioFlag; diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 02874970c04..97b9486f52d 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -634,6 +634,7 @@ bool Kinetics::addReaction(shared_ptr r) m_ropr.push_back(0.0); m_ropnet.push_back(0.0); m_perturb.push_back(1.0); + m_dH.push_back(0.0); return true; } diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 10e699b5a6a..ba2899523cd 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -559,6 +559,14 @@ InterfaceReaction::InterfaceReaction(const Composition& reactants_, reaction_type = INTERFACE_RXN; } +void InterfaceReaction::calculateRateCoeffUnits(const Kinetics& kin) +{ + ElementaryReaction::calculateRateCoeffUnits(kin); + if (is_sticking_coefficient || input.hasKey("sticking-coefficient")) { + rate_units = Units(1.0); // sticking coefficients are dimensionless + } +} + void InterfaceReaction::getParameters(AnyMap& reactionNode) const { ElementaryReaction::getParameters(reactionNode); @@ -611,6 +619,96 @@ void ElectrochemicalReaction::getParameters(AnyMap& reactionNode) const } } +BlowersMaselReaction::BlowersMaselReaction() + : allow_negative_pre_exponential_factor(false) +{ + reaction_type = BLOWERSMASEL_RXN; +} + +void BlowersMaselReaction::validate() +{ + Reaction::validate(); + if (!allow_negative_pre_exponential_factor && + rate.preExponentialFactor() < 0) { + throw CanteraError("BlowersMaselReaction::validate", + "Undeclared negative pre-exponential factor found in reaction '" + + equation() + "'"); + } +} + +BlowersMaselReaction::BlowersMaselReaction(const Composition& reactants_, + const Composition& products_, + const BlowersMasel& rate_) + : Reaction(reactants_, products_) + , rate(rate_) + , allow_negative_pre_exponential_factor(false) +{ + reaction_type = BLOWERSMASEL_RXN; +} + +void BlowersMaselReaction::getParameters(AnyMap& reactionNode) const +{ + Reaction::getParameters(reactionNode); + if (allow_negative_pre_exponential_factor) { + reactionNode["negative-A"] = true; + } + AnyMap rateNode; + rate.getParameters(rateNode, rate_units); + reactionNode["rate-constant"] = std::move(rateNode); +} + +BlowersMaselInterfaceReaction::BlowersMaselInterfaceReaction() + : is_sticking_coefficient(false) + , use_motz_wise_correction(false) +{ + reaction_type = BMINTERFACE_RXN; +} + +BlowersMaselInterfaceReaction::BlowersMaselInterfaceReaction(const Composition& reactants_, + const Composition& products_, + const BlowersMasel& rate_, + bool isStick) + : BlowersMaselReaction(reactants_, products_, rate_) + , is_sticking_coefficient(isStick) + , use_motz_wise_correction(false) +{ + reaction_type = BMINTERFACE_RXN; +} + +void BlowersMaselInterfaceReaction::calculateRateCoeffUnits(const Kinetics& kin) +{ + BlowersMaselReaction::calculateRateCoeffUnits(kin); + if (is_sticking_coefficient || input.hasKey("sticking-coefficient")) { + rate_units = Units(1.0); // sticking coefficients are dimensionless + } +} + +void BlowersMaselInterfaceReaction::getParameters(AnyMap& reactionNode) const +{ + BlowersMaselReaction::getParameters(reactionNode); + if (is_sticking_coefficient) { + reactionNode["sticking-coefficient"] = std::move(reactionNode["rate-constant"]); + reactionNode.erase("rate-constant"); + } + if (use_motz_wise_correction) { + reactionNode["Motz-Wise"] = true; + } + if (!sticking_species.empty()) { + reactionNode["sticking-species"] = sticking_species; + } + if (!coverage_deps.empty()) { + AnyMap deps; + for (const auto& d : coverage_deps) { + AnyMap dep; + dep["a"] = d.second.a; + dep["m"] = d.second.m; + dep["E"].setQuantity(d.second.E, "K", true); + deps[d.first] = std::move(dep); + } + reactionNode["coverage-dependencies"] = std::move(deps); + } +} + Arrhenius readArrhenius(const XML_Node& arrhenius_node) { return Arrhenius(getFloat(arrhenius_node, "A", "toSI"), @@ -734,6 +832,32 @@ void readEfficiencies(ThirdBody& tbody, const AnyMap& node) } } +BlowersMasel readBlowersMasel(const Reaction& R, const AnyValue& rate, + const Kinetics& kin, const UnitSystem& units, + int pressure_dependence=0) +{ + double A, b, Ta0, w; + Units rc_units = R.rate_units; + if (pressure_dependence) { + Units rxn_phase_units = kin.thermo(kin.reactionPhaseIndex()).standardConcentrationUnits(); + rc_units *= rxn_phase_units.pow(-pressure_dependence); + } + if (rate.is()) { + auto& rate_map = rate.as(); + A = units.convert(rate_map["A"], rc_units); + b = rate_map["b"].asDouble(); + Ta0 = units.convertActivationEnergy(rate_map["Ea0"], "K"); + w = units.convertActivationEnergy(rate_map["w"], "K"); + } else { + auto& rate_vec = rate.asVector(4); + A = units.convert(rate_vec[0], rc_units); + b = rate_vec[1].asDouble(); + Ta0 = units.convertActivationEnergy(rate_vec[2], "K"); + w = units.convertActivationEnergy(rate_vec[3], "K"); + } + return BlowersMasel(A, b, Ta0, w); +} + void setupReaction(Reaction& R, const XML_Node& rxn_node) { // Reactant and product stoichiometries @@ -838,8 +962,8 @@ void setupReaction(Reaction& R, const AnyMap& node, const Kinetics& kin) R.allow_negative_orders = node.getBool("negative-orders", false); R.allow_nonreactant_orders = node.getBool("nonreactant-orders", false); - R.calculateRateCoeffUnits(kin); R.input = node; + R.calculateRateCoeffUnits(kin); } void setupElementaryReaction(ElementaryReaction& R, const XML_Node& rxn_node) @@ -1118,7 +1242,6 @@ void setupInterfaceReaction(InterfaceReaction& R, const AnyMap& node, R.rate = readArrhenius(R, node["rate-constant"], kin, node.units()); } else if (node.hasKey("sticking-coefficient")) { R.is_sticking_coefficient = true; - R.rate_units = Units(); // sticking coefficients are dimensionless R.rate = readArrhenius(R, node["sticking-coefficient"], kin, node.units()); R.use_motz_wise_correction = node.getBool("Motz-Wise", kin.thermo().input().getBool("Motz-Wise", false)); @@ -1187,6 +1310,53 @@ void setupElectrochemicalReaction(ElectrochemicalReaction& R, "exchange-current-density-formulation", false); } +void setupBlowersMaselReaction(BlowersMaselReaction& R, const AnyMap& node, + const Kinetics& kin) +{ + setupReaction(R, node, kin); + R.allow_negative_pre_exponential_factor = node.getBool("negative-A", false); + R.rate = readBlowersMasel(R, node["rate-constant"], kin, node.units()); +} + +void setupBlowersMaselInterfaceReaction(BlowersMaselInterfaceReaction& R, const AnyMap& node, + const Kinetics& kin) +{ + setupReaction(R, node, kin); + R.allow_negative_pre_exponential_factor = node.getBool("negative-A", false); + + if (node.hasKey("rate-constant")) { + R.rate = readBlowersMasel(R, node["rate-constant"], kin, node.units()); + } else if (node.hasKey("sticking-coefficient")) { + R.is_sticking_coefficient = true; + R.rate = readBlowersMasel(R, node["sticking-coefficient"], kin, node.units()); + R.use_motz_wise_correction = node.getBool("Motz-Wise", + kin.thermo().input().getBool("Motz-Wise", false)); + R.sticking_species = node.getString("sticking-species", ""); + } else { + throw InputFileError("setupBlowersMaselInterfaceReaction", node, + "Reaction must include either a 'rate-constant' or" + " 'sticking-coefficient' node."); + } + + if (node.hasKey("coverage-dependencies")) { + for (const auto& item : node["coverage-dependencies"].as()) { + double a, E, m; + if (item.second.is()) { + auto& cov_map = item.second.as(); + a = cov_map["a"].asDouble(); + m = cov_map["m"].asDouble(); + E = node.units().convertActivationEnergy(cov_map["E"], "K"); + } else { + auto& cov_vec = item.second.asVector(3); + a = cov_vec[0].asDouble(); + m = cov_vec[1].asDouble(); + E = node.units().convertActivationEnergy(cov_vec[2], "K"); + } + R.coverage_deps[item.first] = CoverageDependency(a, E, m); + } + } +} + std::vector > getReactions(const XML_Node& node) { std::vector > all_reactions; diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index fa741741f67..4ca20d400cf 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -143,6 +143,20 @@ ReactionFactory::ReactionFactory() [](Reaction* R, const AnyMap& node, const Kinetics& kin) { setupElectrochemicalReaction(*(ElectrochemicalReaction*)R, node, kin); }); + + // register Blowers Masel reactions + reg("Blowers-Masel", []() { return new BlowersMaselReaction(); }); + reg_AnyMap("Blowers-Masel", + [](Reaction* R, const AnyMap& node, const Kinetics& kin) { + setupBlowersMaselReaction(*(BlowersMaselReaction*)R, node, kin); + }); + + // register surface Blowers Masel reactions + reg("surface-Blowers-Masel", []() { return new BlowersMaselInterfaceReaction(); }); + reg_AnyMap("surface-Blowers-Masel", + [](Reaction* R, const AnyMap& node, const Kinetics& kin) { + setupBlowersMaselInterfaceReaction(*(BlowersMaselInterfaceReaction*)R, node, kin); + }); } bool isElectrochemicalReaction(Reaction& R, const Kinetics& kin) @@ -214,7 +228,7 @@ unique_ptr newReaction(const AnyMap& rxn_node, type = rxn_node["type"].asString(); } - if (kin.thermo().nDim() < 3) { + if (kin.thermo().nDim() < 3 && type == "elementary") { // See if this is an electrochemical reaction: type of // receiving reaction object is unimportant in this case ElementaryReaction testReaction; @@ -225,6 +239,10 @@ unique_ptr newReaction(const AnyMap& rxn_node, type = "interface"; } } + if (kin.thermo().nDim() < 3 && type == "Blowers-Masel") { + // Allow yaml file to specify "Blowers-Masel" for surface reactions + type = "surface-Blowers-Masel"; + } Reaction* R; try { diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index d1fe62efaf8..5571aca3042 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -68,6 +68,46 @@ void Arrhenius::getParameters(AnyMap& rateNode, const Units& rate_units) const rateNode.setFlowStyle(); } +BlowersMasel::BlowersMasel() + : m_logA(-1.0E300) + , m_b(0.0) + , m_A(0.0) + , m_w(0.0) + , m_E0(0.0) +{ +} + +BlowersMasel::BlowersMasel(double A, double b, double E0, double w) + : m_b(b) + , m_A(A) + , m_w(w) + , m_E0(E0) +{ + if (m_A <= 0.0) { + m_logA = -1.0E300; + } else { + m_logA = std::log(m_A); + } +} + +void BlowersMasel::getParameters(AnyMap& rateNode, const Units& rate_units) const +{ + if (rate_units.factor() != 0.0) { + rateNode["A"].setQuantity(preExponentialFactor(), rate_units); + } else { + rateNode["A"] = preExponentialFactor(); + // This can't be converted to a different unit system because the dimensions of + // the rate constant were not set. Can occur if the reaction was created outside + // the context of a Kinetics object and never added to a Kinetics object. + rateNode["__unconvertible__"] = true; + } + + rateNode["b"] = temperatureExponent(); + rateNode["Ea0"].setQuantity(activationEnergy_R0(), "K", true); + rateNode["w"].setQuantity(bondEnergy(), "K", true); + rateNode.setFlowStyle(); +} + SurfaceArrhenius::SurfaceArrhenius() : m_b(0.0) , m_E(0.0) @@ -198,4 +238,38 @@ ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, } } +BMSurfaceArrhenius::BMSurfaceArrhenius() + : m_b(0.0) + , m_A(0.0) + , m_E0(0.0) + , m_w(0.0) + , m_acov(0.0) + , m_ecov(0.0) + , m_mcov(0.0) +{ +} + +BMSurfaceArrhenius::BMSurfaceArrhenius(double A, double b, double Ta, double w) + : m_b(b) + , m_A(A) + , m_E0(Ta) + , m_w(w) + , m_acov(0.0) + , m_ecov(0.0) + , m_mcov(0.0) +{ +} + +void BMSurfaceArrhenius::addCoverageDependence(size_t k, double a, + double m, double e) +{ + m_sp.push_back(k); + m_ac.push_back(a); + m_ec.push_back(e); + if (m != 0.0) { + m_msp.push_back(k); + m_mc.push_back(m); + } +} + } diff --git a/test/data/BM-ptcombust-Motz-Wise.yaml b/test/data/BM-ptcombust-Motz-Wise.yaml new file mode 100644 index 00000000000..2749a1a2ede --- /dev/null +++ b/test/data/BM-ptcombust-Motz-Wise.yaml @@ -0,0 +1,115 @@ +description: |- + This is a test input file for Motz Wise correction for interface + Blowers Masel reaction tests. This file is tweaked from the model + made by Deutschmann, and it has no physical sense, which should only + be used for tests. + See https://www.detchem.com/mechanisms for original model. + Ref:- 1.) Deutschman et al., 26th Symp. (Intl.) on Combustion, 1996 + pp. 1747-1754 + ----------------------------------------------------------------------- +generator: cti2yaml +cantera-version: 2.5.0 +date: Wed, 11 Dec 2019 16:59:14 -0500 +input-files: [ptcombust.cti] + +units: {length: cm, quantity: mol, activation-energy: J/mol} + +phases: +- name: gas + thermo: ideal-gas + elements: [O, H, C, N, Ar] + species: + - gri30.yaml/species: [H2, H, O, O2, OH, H2O, HO2, H2O2, C, CH, CH2, CH2(S), + CH3, CH4, CO, CO2, HCO, CH2O, CH2OH, CH3O, CH3OH, C2H, C2H2, C2H3, + C2H4, C2H5, C2H6, HCCO, CH2CO, HCCOH, AR, N2] + skip-undeclared-elements: true + kinetics: gas + reactions: + - gri30.yaml/reactions: declared-species + transport: mixture-averaged + state: + T: 300.0 + P: 1.01325e+05 + X: {CH4: 0.095, O2: 0.21, AR: 0.79} +- name: Pt_surf + thermo: ideal-surface + elements: [Pt, H, O, C] + species: [PT(S), H(S), H2O(S), OH(S), O(S)] + Motz-Wise: true + kinetics: surface + reactions: all + state: + T: 900.0 + coverages: {O(S): 0.0, PT(S): 0.5, H(S): 0.5} + site-density: 2.7063e-09 + +species: +- name: PT(S) + composition: {Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +- name: H(S) + composition: {H: 1, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [-1.3029877, 5.4173199e-03, 3.1277972e-07, -3.2328533e-09, 1.136282e-12, + -4227.7075, 5.8743238] + - [1.0696996, 1.543223e-03, -1.5500922e-07, -1.6573165e-10, 3.8359347e-14, + -5054.6128, -7.1555238] +- name: H2O(S) + composition: {O: 1, H: 2, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [-2.7651553, 0.013315115, 1.0127695e-06, -7.1820083e-09, 2.2813776e-12, + -3.6398055e+04, 12.098145] + - [2.5803051, 4.9570827e-03, -4.6894056e-07, -5.2633137e-10, 1.1998322e-13, + -3.8302234e+04, -17.406322] +- name: OH(S) + composition: {O: 1, H: 1, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [-2.0340881, 9.3662683e-03, 6.6275214e-07, -5.2074887e-09, 1.7088735e-12, + -2.5319949e+04, 8.9863186] + - [1.8249973, 3.2501565e-03, -3.1197541e-07, -3.4603206e-10, 7.9171472e-14, + -2.6685492e+04, -12.280891] +- name: O(S) + composition: {O: 1, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [-0.94986904, 7.4042305e-03, -1.0451424e-06, -6.112042e-09, 3.3787992e-12, + -1.3209912e+04, 3.6137905] + - [1.945418, 9.1761647e-04, -1.1226719e-07, -9.9099624e-11, 2.4307699e-14, + -1.4005187e+04, -11.531663] + +reactions: +- equation: 2 H(S) => H2 + 2 PT(S) # Reaction 1 + type: Blowers-Masel + rate-constant: {A: 3.7e+21, b: 0, Ea0: 67400, w: 1000000} + coverage-dependencies: + H(S): {a: 0.0, m: 0.0, E: -6000.0} +- equation: H + PT(S) => H(S) # Reaction 2 + type: Blowers-Masel + sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 1000000} +- equation: H2O + PT(S) => H2O(S) # Reaction 3 + type: Blowers-Masel + sticking-coefficient: {A: 0.75, b: 0, Ea0: 0, w: 1000000} +- equation: O2 + 2 PT(S) => 2 O(S) # Reaction 4 + type: Blowers-Masel + sticking-coefficient: {A: 0.023, b: 0, Ea0: 0, w: 1000000} + Motz-Wise: false +- equation: OH + PT(S) => OH(S) # Reaction 5 + type: Blowers-Masel + sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 100000} + Motz-Wise: true diff --git a/test/data/BM_test.yaml b/test/data/BM_test.yaml new file mode 100644 index 00000000000..665bfb882c5 --- /dev/null +++ b/test/data/BM_test.yaml @@ -0,0 +1,116 @@ +description: |- + This is a test input file for interface Blowers Masel reaction tests. + This file is tweaked from the model made by Deutschmann, and it has + no physical sense, which should only be used for tests. + See https://www.detchem.com/mechanisms for original model. + Ref:- 1.) Deutschman et al., 26th Symp. (Intl.) on Combustion, 1996 + pp. 1747-1754 + ----------------------------------------------------------------------- +generator: cti2yaml +cantera-version: 2.5.0 +date: Wed, 11 Dec 2019 16:59:14 -0500 +input-files: [ptcombust.cti] + +units: {length: cm, quantity: mol, activation-energy: J/mol} + +phases: +- name: gas + thermo: ideal-gas + elements: [O, H, C, N, Ar] + species: + - gri30.yaml/species: [H2, H, O, O2, OH, H2O, CH4, AR, N2] + skip-undeclared-elements: true + kinetics: gas + reactions: [gas-reactions] + transport: mixture-averaged + state: + T: 300.0 + P: 1.01325e+05 +- name: Pt_surf + thermo: ideal-surface + elements: [Pt, H, O, C] + species: [{Pt-surf-species: all}] + kinetics: surface + reactions: [Pt-reactions] + state: + T: 900.0 + coverages: {O(S): 0.0, PT(S): 0.5, H(S): 0.5} + site-density: 2.7063e-09 + +Pt-surf-species: +- name: PT(S) + composition: {Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +- name: H(S) + composition: {H: 1, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [-1.3029877, 5.4173199e-03, 3.1277972e-07, -3.2328533e-09, 1.136282e-12, + -4227.7075, 5.8743238] + - [1.0696996, 1.543223e-03, -1.5500922e-07, -1.6573165e-10, 3.8359347e-14, + -5054.6128, -7.1555238] +- name: H2O(S) + composition: {O: 1, H: 2, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [-2.7651553, 0.013315115, 1.0127695e-06, -7.1820083e-09, 2.2813776e-12, + -3.6398055e+04, 12.098145] + - [2.5803051, 4.9570827e-03, -4.6894056e-07, -5.2633137e-10, 1.1998322e-13, + -3.8302234e+04, -17.406322] +- name: OH(S) + composition: {O: 1, H: 1, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [-2.0340881, 9.3662683e-03, 6.6275214e-07, -5.2074887e-09, 1.7088735e-12, + -2.5319949e+04, 8.9863186] + - [1.8249973, 3.2501565e-03, -3.1197541e-07, -3.4603206e-10, 7.9171472e-14, + -2.6685492e+04, -12.280891] +- name: O(S) + composition: {O: 1, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 3000.0] + data: + - [-0.94986904, 7.4042305e-03, -1.0451424e-06, -6.112042e-09, 3.3787992e-12, + -1.3209912e+04, 3.6137905] + - [1.945418, 9.1761647e-04, -1.1226719e-07, -9.9099624e-11, 2.4307699e-14, + -1.4005187e+04, -11.531663] + +gas-reactions: +- equation: O + H2 <=> H + OH # Reaction 1 + type: Blowers-Masel + rate-constant: {A: 38700, b: 2.7, Ea0: 2.619184e4, w: 4.184e9} + +Pt-reactions: +- equation: 2 H(S) => H2 + 2 PT(S) # Reaction 1 + type: Blowers-Masel + rate-constant: {A: 3.7e+21, b: 0, Ea0: 67400, w: 1000000} + coverage-dependencies: + H(S): {a: 0.0, m: 0.0, E: -6000.0} +- equation: H + PT(S) => H(S) # Reaction 2 + type: Blowers-Masel + sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 1000000} +- equation: H2O + PT(S) => H2O(S) # Reaction 3 + type: Blowers-Masel + sticking-coefficient: {A: 0.75, b: 0, Ea0: 0, w: 1000000} +- equation: O2 + 2 PT(S) => 2 O(S) # Reaction 4 + type: Blowers-Masel + sticking-coefficient: {A: 0.023, b: 0, Ea0: 0, w: 1000000} +- equation: OH + PT(S) => OH(S) # Reaction 5 + type: Blowers-Masel + sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 100000} +- equation: H2 + 2 PT(s) => 2 H(s) + type: Blowers-Masel + rate-constant: {A: 4.4579e10, b: 0.5, Ea0: 0.0, w: 1000000} + orders: {H2: 1.0, PT(s): 1.0} diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index c6a89bbf3a5..15584d87871 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -173,6 +173,39 @@ TEST(Reaction, ChebyshevFromYaml) EXPECT_NEAR(CR.rate.updateRC(std::log(T), 1.0/T), 130512.2773948636, 1e-9); } +TEST(Reaction, BlowersMaselFromYaml) +{ + auto sol = newSolution("gri30.yaml"); + AnyMap rxn = AnyMap::fromYamlString( + "{equation: O + H2 <=> H + OH," + " type: Blowers-Masel," + " rate-constant: [-3.87e+04 cm^3/mol/s, 2.7, 6260.0 cal/mol, 1e9 cal/mol]," + " negative-A: true}"); + + auto R = newReaction(rxn, *(sol->kinetics())); + EXPECT_EQ(R->reactants.at("H2"), 1); + EXPECT_EQ(R->products.at("OH"), 1); + EXPECT_EQ(R->reaction_type, BLOWERSMASEL_RXN); + + auto ER = dynamic_cast(*R); + doublereal E_intrinsic = 6260 / GasConst_cal_mol_K * GasConstant; // J/kmol + doublereal H_big = 5 * E_intrinsic; + doublereal H_small = -5 * E_intrinsic; + doublereal H_mid = 4 * E_intrinsic; + doublereal w = 1e9 / GasConst_cal_mol_K * GasConstant; // J/kmol + doublereal vp = 2 * w * ((w + E_intrinsic) / (w - E_intrinsic)); + doublereal Ea = (w + H_mid / 2) * (vp - 2 * w + H_mid) * (vp - 2 * w + H_mid) + / (vp * vp - 4 * w * w + H_mid * H_mid ); + EXPECT_DOUBLE_EQ(ER.rate.preExponentialFactor(), -38.7); + EXPECT_DOUBLE_EQ(ER.rate.activationEnergy_R0(), 6260 / GasConst_cal_mol_K); + EXPECT_DOUBLE_EQ(ER.rate.bondEnergy(), 1e9 / GasConst_cal_mol_K); + EXPECT_DOUBLE_EQ(ER.rate.activationEnergy_R(H_big), H_big / GasConstant); + EXPECT_DOUBLE_EQ(ER.rate.activationEnergy_R(H_small), 0); + EXPECT_NEAR(ER.rate.activationEnergy_R(H_mid), Ea / GasConstant, 1e-7); + EXPECT_TRUE(ER.allow_negative_pre_exponential_factor); + EXPECT_FALSE(ER.allow_negative_orders); +} + TEST(Kinetics, GasKineticsFromYaml1) { AnyMap infile = AnyMap::fromYamlFile("ideal-gas.yaml"); @@ -240,6 +273,31 @@ TEST(Kinetics, InterfaceKineticsFromYaml) EXPECT_TRUE(IR3->is_sticking_coefficient); } +TEST(Kinetics, BMInterfaceKineticsFromYaml) +{ + shared_ptr gas(newPhase("BM_test.yaml", "gas")); + shared_ptr surf_tp(newPhase("BM_test.yaml", "Pt_surf")); + shared_ptr surf = std::dynamic_pointer_cast(surf_tp); + std::vector phases{surf_tp.get(), gas.get()}; + auto kin = newKinetics(phases, "BM_test.yaml", "Pt_surf"); + EXPECT_EQ(kin->nReactions(), (size_t) 6); + EXPECT_EQ(kin->nTotalSpecies(), (size_t) 14); + auto R1 = kin->reaction(5); + auto IR1 = std::dynamic_pointer_cast(R1); + EXPECT_DOUBLE_EQ(R1->orders["PT(s)"], 1.0); + EXPECT_DOUBLE_EQ(IR1->rate.preExponentialFactor(), 4.4579e7); + + auto R2 = kin->reaction(0); + auto IR2 = std::dynamic_pointer_cast(R2); + EXPECT_DOUBLE_EQ(IR2->rate.preExponentialFactor(), 3.7e20); + EXPECT_DOUBLE_EQ(IR2->coverage_deps["H(S)"].E, -6e6 / GasConstant); + EXPECT_FALSE(IR2->is_sticking_coefficient); + + auto R3 = kin->reaction(1); + auto IR3 = std::dynamic_pointer_cast(R3); + EXPECT_TRUE(IR3->is_sticking_coefficient); +} + TEST(Kinetics, ElectrochemFromYaml) { shared_ptr graphite(newPhase("surface-phases.yaml", "graphite")); @@ -484,3 +542,25 @@ TEST_F(ReactionToYaml, unconvertible2) params.setUnits(U); EXPECT_THROW(params.applyUnits(), CanteraError); } + +TEST_F(ReactionToYaml, BlowersMasel) +{ + soln = newSolution("BM_test.yaml", "gas"); + soln->thermo()->setState_TPY(1100, 0.1 * OneAtm, "O:0.01, H2:0.8, O2:0.19"); + duplicateReaction(0); + EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); + compareReactions(); +} + +TEST_F(ReactionToYaml, BlowersMaselInterface) +{ + auto gas = newSolution("BM_test.yaml", "gas"); + soln = newSolution("BM_test.yaml", "Pt_surf", "None", {gas}); + gas->thermo()->setState_TPY(1100, 0.1 * OneAtm, "O:0.01, H2:0.8, O2:0.19"); + soln->thermo()->setState_TP(1100, 0.1 * OneAtm); + auto surf = std::dynamic_pointer_cast(soln->thermo()); + surf->setCoveragesByName("H(S):0.1, PT(S):0.8, H2O(S):0.1"); + duplicateReaction(0); + EXPECT_TRUE(std::dynamic_pointer_cast(duplicate)); + compareReactions(); +}