From 2fdf504c9972d71af243cf31d0dc399a3e667132 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 4 Mar 2021 16:00:43 -0500 Subject: [PATCH 01/20] [Kinetics/Cython] Add Blowers Masel reaction class in GasKinetics This commit inclues Blowers Masel rate class and Blowers Masel reaction class, which are aimed to be used in GasKinetics. The changes are implemented in both C++ and Python/Cython interfaces. --- include/cantera/kinetics/GasKinetics.h | 3 + include/cantera/kinetics/Kinetics.h | 3 + include/cantera/kinetics/RateCoeffMgr.h | 7 ++ include/cantera/kinetics/Reaction.h | 12 +++ include/cantera/kinetics/RxnRates.h | 103 +++++++++++++++++++++++ include/cantera/kinetics/reaction_defs.h | 5 ++ interfaces/cython/cantera/_cantera.pxd | 20 +++++ interfaces/cython/cantera/reaction.pyx | 85 +++++++++++++++++++ src/kinetics/GasKinetics.cpp | 20 +++++ src/kinetics/Kinetics.cpp | 1 + src/kinetics/Reaction.cpp | 60 +++++++++++++ src/kinetics/RxnRates.cpp | 24 ++++++ 12 files changed, 343 insertions(+) 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/Kinetics.h b/include/cantera/kinetics/Kinetics.h index a355f8a687c..78ecac9f144 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -944,6 +944,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..b4a163fafc9 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -416,6 +416,18 @@ 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 validate(); + BlowersMasel rate; + bool allow_negative_pre_exponential_factor; +}; + //! Create Reaction objects for all `` nodes in an XML document. //! //! The `` nodes are assumed to be children of the `` diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index e2a90951e94..8aae365b145 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -95,6 +95,109 @@ class Arrhenius doublereal m_logA, m_b, m_E, m_A; }; +//! Blowers Masel reaction rate type depends on enthalpy +/** + * The Blowers Masel approximation is written by [Paul Blowrs, + * Rich Masel(2004). Engineering Approximations For Activation Energies in Hydrogen + * Transfer Reactions. Eq (10)-(11)] to adjust the activation energy based on + * enthalpy change of a reaction: + * + * \f[ + * Ea = 0 if DeltaH < -4E_0 + * Ea = DeltaH if DeltaH > 4E_0 + * Ea = \frac{(w_0 + DeltaH / 2)(V_P - 2w_0 + DeltaH)^2}{(V_P^^2 - 4w_0^2 + DeltaH^2)} + * \f] + * where + * \f[ + * V_P = \frac{2w_0 (w_0 + E_0)}{w_0 - E_0} + * \f] + * and \f$ w_0 \f$ is theaverage of the bond energy of the + * bond breaking and that being formed, which can be approximated as + * arbitrary high value like 1000kJ/mol as long as \f$ w_0 >= 2 E_0 \f$ + * The rate constant then can be calculated by Arrhenius function. + * \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 bond energy of the bond being formed or broken, in temperature unts. Kelvin. + + BlowersMasel(doublereal A, doublereal b, doublereal E0, doublereal w); + + //! 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) { + 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 energy *w* divided by the gas constant[K] + doublereal bondEnergy() const { + return m_w; + } + +protected: + doublereal m_logA, m_b, m_E, m_A, m_w, m_E0; +}; /** * An Arrhenius rate with coverage-dependent terms. diff --git a/include/cantera/kinetics/reaction_defs.h b/include/cantera/kinetics/reaction_defs.h index 7a1aedee97a..45210c5f427 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 diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index cb9dce619fd..a772ec25e8f 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -450,6 +450,22 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxTestReaction "Cantera::TestReaction" (CxxReaction2): 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) @@ -1101,6 +1117,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/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index eb0135f1c68..e36b2ae5bc7 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -910,6 +910,91 @@ cdef class ChebyshevReaction(Reaction): r.rate.update_C(&logP) return r.rate.updateRC(logT, recipT) +cdef class BlowersMasel: + 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]. + """ + return self.rate.activationEnergy_R(dH) * gas_constant + + property bond_energy: + 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): """ diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index ef008ab112f..270c7636b22 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,12 @@ 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 +359,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() == "BlowersMasel") { + modifyBlowersMaselReaction(i, dynamic_cast(*rNew)); } else { throw CanteraError("GasKinetics::modifyReaction", "Unknown reaction type specified: '{}'", rNew->type()); @@ -380,6 +395,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/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..6b84f7ec363 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -611,6 +611,32 @@ void ElectrochemicalReaction::getParameters(AnyMap& reactionNode) const } } +BlowersMaselReaction::BlowersMaselReaction() + : Reaction(BLOWERSMASEL_RXN) + , allow_negative_pre_exponential_factor(false) +{ +} + +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(BLOWERSMASEL_RXN, reactants_, products_) + , rate(rate_) + , allow_negative_pre_exponential_factor(false) +{ +} + Arrhenius readArrhenius(const XML_Node& arrhenius_node) { return Arrhenius(getFloat(arrhenius_node, "A", "toSI"), @@ -734,6 +760,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 @@ -1187,6 +1239,14 @@ 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()); +} + std::vector > getReactions(const XML_Node& node) { std::vector > all_reactions; diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index d1fe62efaf8..c4cbcbfce53 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -68,6 +68,30 @@ void Arrhenius::getParameters(AnyMap& rateNode, const Units& rate_units) const rateNode.setFlowStyle(); } +BlowersMasel::BlowersMasel() + : m_logA(-1.0E300) + , m_b(0.0) + , m_E(0.0) + , m_A(0.0) + , m_w(0.0) + , m_E0(0.0) +{ +} + +BlowersMasel::BlowersMasel(doublereal A, doublereal b, doublereal E0, doublereal w) + : m_b(b) + , m_E(E0) + , 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); + } +} + SurfaceArrhenius::SurfaceArrhenius() : m_b(0.0) , m_E(0.0) From 2d6db7690259fe526505fa16a11de08e5bf0724c Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 4 Mar 2021 16:06:47 -0500 Subject: [PATCH 02/20] [Test/Cython] Add unit tests for Blowers Masel reaction and rate classes The C++ test is only for Blowers Masel rate test, and the tests in Python interface are for both the reaction and rate classes. Two minimal yaml files are added as the input file for the tests --- .../cython/cantera/test/test_kinetics.py | 97 +++++++++++++++++++ test/data/BM_example.yaml | 86 ++++++++++++++++ test/data/surface-phases-BM.yaml | 84 ++++++++++++++++ test/kinetics/kineticsFromYaml.cpp | 33 +++++++ 4 files changed, 300 insertions(+) create mode 100644 test/data/BM_example.yaml create mode 100644 test/data/surface-phases-BM.yaml diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index b3c6b7ecf1b..16097d635e5 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_example.yaml') + + gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=self.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_example.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.Species.listFromFile('BM_example.yaml') + 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_example.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') @@ -1214,6 +1289,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_example.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]) diff --git a/test/data/BM_example.yaml b/test/data/BM_example.yaml new file mode 100644 index 00000000000..8c25789004a --- /dev/null +++ b/test/data/BM_example.yaml @@ -0,0 +1,86 @@ +description: |- + This is the minimal example tweaked from Hydrogen-Oxygen + submechanism extracted from GRI-Mech 3.0. It should only + be used for testing the Blowers Masel reaction class in + gas phase. + +generator: ck2yaml +input-files: [h2o2.inp, gri30_tran.dat] +cantera-version: 2.5.0 +date: Wed, 11 Dec 2019 16:59:04 -0500 + +units: {length: cm, time: s, quantity: mol, activation-energy: cal/mol} + +phases: +- name: ohmech + thermo: ideal-gas + elements: [O, H, Ar, N] + species: [H2, H, O, OH, H2O, O2] + kinetics: gas + state: {T: 300.0, P: 1 atm} +species: +- name: H2 + composition: {H: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [2.34433112, 7.98052075e-03, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12, + -917.935173, 0.683010238] + - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14, + -950.158922, -3.20502331] +- name: H + composition: {H: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [2.5, 7.05332819e-13, -1.99591964e-15, 2.30081632e-18, -9.27732332e-22, + 2.54736599e+04, -0.446682853] + - [2.50000001, -2.30842973e-11, 1.61561948e-14, -4.73515235e-18, 4.98197357e-22, + 2.54736599e+04, -0.446682914] +- name: O + composition: {O: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [3.1682671, -3.27931884e-03, 6.64306396e-06, -6.12806624e-09, 2.11265971e-12, + 2.91222592e+04, 2.05193346] + - [2.56942078, -8.59741137e-05, 4.19484589e-08, -1.00177799e-11, 1.22833691e-15, + 2.92175791e+04, 4.78433864] +- name: OH + composition: {O: 1, H: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [3.99201543, -2.40131752e-03, 4.61793841e-06, -3.88113333e-09, 1.3641147e-12, + 3615.08056, -0.103925458] + - [3.09288767, 5.48429716e-04, 1.26505228e-07, -8.79461556e-11, 1.17412376e-14, + 3858.657, 4.4766961] +- name: H2O + composition: {H: 2, O: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [4.19864056, -2.0364341e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12, + -3.02937267e+04, -0.849032208] + - [3.03399249, 2.17691804e-03, -1.64072518e-07, -9.7041987e-11, 1.68200992e-14, + -3.00042971e+04, 4.9667701] +- name: O2 + composition: {O: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [3.78245636, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, + -1063.94356, 3.65767573] + - [3.28253784, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, + -1088.45772, 5.45323129] + +reactions: +- equation: O + H2 <=> H + OH # Reaction 1 + type: Blowers-Masel + rate-constant: {A: 3.87e+04, b: 2.7, Ea0: 6260.0, w: 1e9} diff --git a/test/data/surface-phases-BM.yaml b/test/data/surface-phases-BM.yaml new file mode 100644 index 00000000000..2000629e2b8 --- /dev/null +++ b/test/data/surface-phases-BM.yaml @@ -0,0 +1,84 @@ +phases: +- name: gas + thermo: ideal-gas + species: [{gas-species: all}] + state: {T: 900 K, P: 1 atm, X: {H2: 0.4, Ar: 0.6}} + +- name: Pt-surf + thermo: ideal-surface + kinetics: surface + reactions: [Pt-reactions] + species: [{Pt-surf-species: all}] + site-density: 2.7063e-9 mol/cm^2 + state: {T: 900 K, P: 1 atm, coverages: {Pt(s): 0.5, H(s): 0.4, O(s): 0.1}} + + +gas-species: +- name: H2 + composition: {H: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [2.34433112, 0.00798052075, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12, + -917.935173, 0.683010238] + - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14, + -950.158922, -3.20502331] +- name: H + composition: {H: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [2.5, 7.05332819e-13, -1.99591964e-15, 2.30081632e-18, -9.27732332e-22, + 25473.6599, -0.446682853] + - [2.50000001, -2.30842973e-11, 1.61561948e-14, -4.73515235e-18, 4.98197357e-22, + 25473.6599, -0.446682914] +- name: Ar + composition: {Ar: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 5000.0] + data: + - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.366] + + +Pt-surf-species: +- name: Pt(s) + composition: {Pt: 1} + thermo: + model: constant-cp +- name: H(s) + composition: {H: 1, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300, 1000, 3000] + data: + - [-1.3029877E+00, 5.4173199E-03, 3.1277972E-07, -3.2328533E-09, + 1.1362820E-12, -4.2277075E+03, 5.8743238E+00] + - [1.0696996E+00, 1.5432230E-03, -1.5500922E-07, -1.6573165E-10, + 3.8359347E-14, -5.0546128E+03, -7.1555238E+00] +- name: O(s) + composition: {O: 1, Pt: 1} + thermo: + model: NASA7 + temperature-ranges: [300, 1000, 3000] + data: + - [-9.4986904E-01, 7.4042305E-03, -1.0451424E-06, -6.1120420E-09, + 3.3787992E-12, -1.3209912E+04, 3.6137905E+00] + - [1.9454180E+00, 9.1761647E-04, -1.1226719E-07, -9.9099624E-11, + 2.4307699E-14, -1.4005187E+04, -1.1531663E+01] + +Pt-reactions: +- equation: H2 + 2 Pt(s) => 2 H(s) + type: surface-Blowers-Masel + rate-constant: {A: 4.4579e10 cm^3/mol/s, b: 0.5, Ea0: 0.0 J/mol, w: 1000000 J/mol} + orders: {H2: 1.0, Pt(s): 1.0} +- equation: 2 H(s) => H2 + 2 Pt(s) + type: surface-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}} +- equation: H + Pt(s) => H(s) + type: surface-Blowers-Masel + sticking-coefficient: [1.0, 0.0, 0.0, 1000000] + diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index c6a89bbf3a5..f47bd461e28 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"); From 814db54b8e14fd210c0d88b28e23096a2a6aeea6 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 4 Mar 2021 16:08:08 -0500 Subject: [PATCH 03/20] [Kinetics/Cython] Add Blowers Masel interface reaction and reaction rate classes in InterfaceKinetics The changes are made in both C++ and Cython interfaces. --- include/cantera/kinetics/InterfaceKinetics.h | 20 ++ include/cantera/kinetics/Reaction.h | 46 +++- include/cantera/kinetics/RxnRates.h | 135 ++++++++++++ include/cantera/kinetics/reaction_defs.h | 5 + interfaces/cython/cantera/_cantera.pxd | 6 + interfaces/cython/cantera/reaction.pyx | 69 ++++++ src/kinetics/GasKinetics.cpp | 2 +- src/kinetics/InterfaceKinetics.cpp | 212 +++++++++++++++---- src/kinetics/Reaction.cpp | 65 +++++- src/kinetics/ReactionFactory.cpp | 14 ++ src/kinetics/RxnRates.cpp | 36 ++++ 11 files changed, 565 insertions(+), 45 deletions(-) diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index fec3d801b59..a20c6876dc6 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 BMInterfaceReaction; //! 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, BMInterfaceReaction& r, + bool replace); + //! Temporary work vector of length m_kk vector_fp m_grt; @@ -401,6 +413,14 @@ 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/Reaction.h b/include/cantera/kinetics/Reaction.h index b4a163fafc9..71fe18fa897 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -424,10 +424,49 @@ class BlowersMaselReaction: public Reaction BlowersMaselReaction(const Composition& reactants, const Composition& products, const BlowersMasel& rate); 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 BMInterfaceReaction : public BlowersMaselReaction +{ +public: + BMInterfaceReaction(); + BMInterfaceReaction(const Composition& reactants, const Composition& products, + const BlowersMasel& rate, bool isStick=false); + + 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 `` @@ -499,7 +538,12 @@ void setupElectrochemicalReaction(ElectrochemicalReaction&, //! @internal May be changed without notice in future versions 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 setupBMInterfaceReaction(BMInterfaceReaction&, + const AnyMap&, const Kinetics&); } #endif diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 8aae365b145..3a6906e42d8 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -543,6 +543,141 @@ 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[ + * Ea = 0 if DeltaH < -4E_0 + * Ea = DeltaH if DeltaH > 4E_0 + * Ea = \frac{(w_0 + DeltaH / 2)(V_P - 2w_0 + DeltaH)^2}{(V_P^2 - 4w_0^2 + DeltaH^2)} + * \f] + * where + * \f[ + * V_P = 2w_0 (w_0 + E_0) / (w_0 - E_0) + * \f] + * and \f$ w_0 \f$ is the average of the bond energy of the + * bond breaking and that being formed, which can be approximated as an + * arbitrary high value like 1000kJ/mol as long as \f$ w_0 >= 2 E_0 \f$ + */ +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, doublereal a, + doublereal m, doublereal 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) + if (deltaH_R < -4 * m_E0) { + m_E = 0; + } else if (deltaH_R > 4 * m_E0) { + m_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) + m_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 m_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_E + m_ecov; + } + + //! Return the bond energy *w* divided by the gas constant[K] + doublereal bondEnergy() const { + return m_w; + } + +protected: + doublereal m_b, m_E, 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 45210c5f427..7e3050b6427 100644 --- a/include/cantera/kinetics/reaction_defs.h +++ b/include/cantera/kinetics/reaction_defs.h @@ -88,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 a772ec25e8f..9ebb1b23cbc 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -478,6 +478,12 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cbool is_sticking_coefficient cbool use_motz_wise_correction string sticking_species + + cdef cppclass CxxBMInterfaceReaction "Cantera::BMInterfaceReaction" (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 diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index e36b2ae5bc7..053fd6ff60c 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -1161,3 +1161,72 @@ cdef class InterfaceReaction(ElementaryReaction): def __set__(self, species): cdef CxxInterfaceReaction* r = self.reaction r.sticking_species = stringify(species) + +cdef class BMInterfaceReaction(BlowersMaselReaction): + """ A reaction occurring on an `Interface` (i.e. a surface or an edge) """ + 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 CxxBMInterfaceReaction* 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 CxxBMInterfaceReaction* 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 CxxBMInterfaceReaction* r = self.reaction + return r.is_sticking_coefficient + def __set__(self, stick): + cdef CxxBMInterfaceReaction* 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 CxxBMInterfaceReaction* r = self.reaction + return r.use_motz_wise_correction + def __set__(self, mw): + cdef CxxBMInterfaceReaction* 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 CxxBMInterfaceReaction* r = self.reaction + return pystr(r.sticking_species) + def __set__(self, species): + cdef CxxBMInterfaceReaction* r = self.reaction + r.sticking_species = stringify(species) diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 270c7636b22..6573f5bd7aa 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -359,7 +359,7 @@ 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() == "BlowersMasel") { + } else if (rNew->type() == "Blowers-Masel") { modifyBlowersMaselReaction(i, dynamic_cast(*rNew)); } else { throw CanteraError("GasKinetics::modifyReaction", diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index a1834b7efdb..eeea03152b8 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,64 +509,95 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base) if (!added) { return false; } + if (r_base->reaction_type == BMINTERFACE_RXN) { + BMInterfaceReaction& 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); - m_ProdStanConcReac.push_back(0.0); - + m_ProdStanConcReac.push_back(0.0); return true; } 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) { + BMInterfaceReaction& 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 +695,99 @@ SurfaceArrhenius InterfaceKinetics::buildSurfaceArrhenius( return rate; } +BMSurfaceArrhenius InterfaceKinetics::buildBMSurfaceArrhenius( + size_t i, BMInterfaceReaction& 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/Reaction.cpp b/src/kinetics/Reaction.cpp index 6b84f7ec363..bae362990c6 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -612,9 +612,9 @@ void ElectrochemicalReaction::getParameters(AnyMap& reactionNode) const } BlowersMaselReaction::BlowersMaselReaction() - : Reaction(BLOWERSMASEL_RXN) - , allow_negative_pre_exponential_factor(false) + : allow_negative_pre_exponential_factor(false) { + reaction_type = BLOWERSMASEL_RXN; } void BlowersMaselReaction::validate() @@ -631,10 +631,29 @@ void BlowersMaselReaction::validate() BlowersMaselReaction::BlowersMaselReaction(const Composition& reactants_, const Composition& products_, const BlowersMasel& rate_) - : Reaction(BLOWERSMASEL_RXN, reactants_, products_) + : Reaction(reactants_, products_) , rate(rate_) , allow_negative_pre_exponential_factor(false) { + reaction_type = BLOWERSMASEL_RXN; +} + +BMInterfaceReaction::BMInterfaceReaction() + : is_sticking_coefficient(false) + , use_motz_wise_correction(false) +{ + reaction_type = BMINTERFACE_RXN; +} + +BMInterfaceReaction::BMInterfaceReaction(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; } Arrhenius readArrhenius(const XML_Node& arrhenius_node) @@ -1247,6 +1266,46 @@ void setupBlowersMaselReaction(BlowersMaselReaction& R, const AnyMap& node, R.rate = readBlowersMasel(R, node["rate-constant"], kin, node.units()); } +void setupBMInterfaceReaction(BMInterfaceReaction& 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_units = Units(); // sticking coefficients are dimensionless + 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("setupBMInterfaceReaction", 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..b920f13fd38 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 BMInterfaceReaction(); }); + reg_AnyMap("surface-Blowers-Masel", + [](Reaction* R, const AnyMap& node, const Kinetics& kin) { + setupBMInterfaceReaction(*(BMInterfaceReaction*)R, node, kin); + }); } bool isElectrochemicalReaction(Reaction& R, const Kinetics& kin) diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index c4cbcbfce53..320d8fca239 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -222,4 +222,40 @@ ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, } } +BMSurfaceArrhenius::BMSurfaceArrhenius() + : m_b(0.0) + , m_E(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_E(Ta) + , 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, doublereal a, + doublereal m, doublereal 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); + } +} + } From 79329ec06d7f90e6831f103222ea109655c6b18d Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 4 Mar 2021 16:11:24 -0500 Subject: [PATCH 04/20] [Test/Cython] The unit tests for Blowers Masel interface reaction and reaction rates classes The tests are added in both C++ and Python interfaces. Two minimal yaml files are included for the tests. --- data/BM-ptcombust-Motz-Wise.yaml | 115 ++++++++++++++++++ data/BM_ptcombust.yaml | 111 +++++++++++++++++ .../cython/cantera/test/test_kinetics.py | 87 +++++++++++++ test/kinetics/kineticsFromYaml.cpp | 25 ++++ 4 files changed, 338 insertions(+) create mode 100644 data/BM-ptcombust-Motz-Wise.yaml create mode 100644 data/BM_ptcombust.yaml diff --git a/data/BM-ptcombust-Motz-Wise.yaml b/data/BM-ptcombust-Motz-Wise.yaml new file mode 100644 index 00000000000..ba4fa5acc8f --- /dev/null +++ b/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: surface-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: surface-Blowers-Masel + sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 1000000} +- equation: H2O + PT(S) => H2O(S) # Reaction 3 + type: surface-Blowers-Masel + sticking-coefficient: {A: 0.75, b: 0, Ea0: 0, w: 1000000} +- equation: O2 + 2 PT(S) => 2 O(S) # Reaction 4 + type: surface-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: surface-Blowers-Masel + sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 100000} + Motz-Wise: true \ No newline at end of file diff --git a/data/BM_ptcombust.yaml b/data/BM_ptcombust.yaml new file mode 100644 index 00000000000..4dfb09695cb --- /dev/null +++ b/data/BM_ptcombust.yaml @@ -0,0 +1,111 @@ +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, 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)] + 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: surface-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: surface-Blowers-Masel + sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 1000000} +- equation: H2O + PT(S) => H2O(S) # Reaction 3 + type: surface-Blowers-Masel + sticking-coefficient: {A: 0.75, b: 0, Ea0: 0, w: 1000000} +- equation: O2 + 2 PT(S) => 2 O(S) # Reaction 4 + type: surface-Blowers-Masel + sticking-coefficient: {A: 0.023, b: 0, Ea0: 0, w: 1000000} +- equation: OH + PT(S) => OH(S) # Reaction 5 + type: surface-Blowers-Masel + sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 100000} \ No newline at end of file diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 16097d635e5..0e3d0055423 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -1183,6 +1183,34 @@ def test_interface(self): self.assertNear(surf1.net_rates_of_progress[1], surf2.net_rates_of_progress[0]) + def test_BMinterface(self): + gas = ct.Solution('BM_ptcombust.yaml') + surf1 = ct.Interface('BM_ptcombust.yaml', 'Pt_surf', [gas]) + r1 = ct.BMInterfaceReaction() + 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) @@ -1369,3 +1397,62 @@ 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('BM_ptcombust.yaml', 'gas') + surf = ct.Interface('BM_ptcombust.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 + k1 = surf.forward_rate_constants[0] + surf.coverages = 'O(S):0.2, PT(S):0.4, H(S):0.4' + k2 = surf.forward_rate_constants[0] + surf.coverages = 'O(S):0.2, PT(S):0.6, H(S):0.2' + k3 = surf.forward_rate_constants[0] + self.assertNotAlmostEqual(k1, k2) + self.assertNear(k2, k3) + + def test_modify_BMsticking(self): + gas = ct.Solution('BM_ptcombust.yaml', 'gas') + surf = ct.Interface('BM_ptcombust.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('BM_ptcombust.yaml', 'gas') + surf1 = ct.Interface('BM_ptcombust.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/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index f47bd461e28..35807933f4b 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -273,6 +273,31 @@ TEST(Kinetics, InterfaceKineticsFromYaml) EXPECT_TRUE(IR3->is_sticking_coefficient); } +TEST(Kinetics, BMInterfaceKineticsFromYaml) +{ + shared_ptr gas(newPhase("surface-phases-BM.yaml", "gas")); + shared_ptr surf_tp(newPhase("surface-phases-BM.yaml", "Pt-surf")); + shared_ptr surf = std::dynamic_pointer_cast(surf_tp); + std::vector phases{surf_tp.get(), gas.get()}; + auto kin = newKinetics(phases, "surface-phases-BM.yaml", "Pt-surf"); + EXPECT_EQ(kin->nReactions(), (size_t) 3); + EXPECT_EQ(kin->nTotalSpecies(), (size_t) 6); + auto R1 = kin->reaction(0); + 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(1); + 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(2); + auto IR3 = std::dynamic_pointer_cast(R3); + EXPECT_TRUE(IR3->is_sticking_coefficient); +} + TEST(Kinetics, ElectrochemFromYaml) { shared_ptr graphite(newPhase("surface-phases.yaml", "graphite")); From 2e9a7c8217f34c152ec24f683a85ecf81ba3b9df Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 25 Mar 2021 15:56:07 -0400 Subject: [PATCH 05/20] reduce the test input files and modify the tests to keep them working --- .../cython/cantera/test/test_kinetics.py | 36 ++++---- .../data}/BM-ptcombust-Motz-Wise.yaml | 13 +-- test/data/BM_example.yaml | 86 ------------------- .../data/BM_test.yaml | 37 ++++---- test/data/surface-phases-BM.yaml | 84 ------------------ test/kinetics/kineticsFromYaml.cpp | 26 +++--- 6 files changed, 61 insertions(+), 221 deletions(-) rename {data => test/data}/BM-ptcombust-Motz-Wise.yaml (95%) delete mode 100644 test/data/BM_example.yaml rename data/BM_ptcombust.yaml => test/data/BM_test.yaml (82%) delete mode 100644 test/data/surface-phases-BM.yaml diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 0e3d0055423..4202b0d5a8b 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -1086,10 +1086,10 @@ 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_example.yaml') + gas1 = ct.Solution('BM_test.yaml') gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', - species=self.species, reactions=[r]) + species=gas1.species(), reactions=[r]) gas1.TP = self.gas.TP gas2.TP = self.gas.TP @@ -1102,13 +1102,13 @@ def test_BlowersMasel(self): gas1.net_rates_of_progress[0], rtol=1e-7) def test_Blowers_Masel_rate(self): - gas = ct.Solution('BM_example.yaml') + 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.Species.listFromFile('BM_example.yaml') + 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) @@ -1123,7 +1123,7 @@ def test_negative_A_blowersmasel(self): species=species, reactions=[r]) def test_Blowers_Masel_change_enthalpy(self): - gas = ct.Solution('BM_example.yaml') + gas = ct.Solution('BM_test.yaml') r = gas.reaction(0) E0 = r.rate.intrinsic_activation_energy w = r.rate.bond_energy @@ -1183,10 +1183,11 @@ def test_interface(self): self.assertNear(surf1.net_rates_of_progress[1], surf2.net_rates_of_progress[0]) - def test_BMinterface(self): - gas = ct.Solution('BM_ptcombust.yaml') - surf1 = ct.Interface('BM_ptcombust.yaml', 'Pt_surf', [gas]) - r1 = ct.BMInterfaceReaction() + 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) @@ -1318,7 +1319,7 @@ def test_modify_chebyshev(self): self.assertNear(kf[4], kf[5]) def test_modify_BlowersMasel(self): - gas = ct.Solution('BM_example.yaml') + 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) @@ -1399,8 +1400,9 @@ def test_motz_wise(self): self.assertNear(2.0 * k1[9], k2[9]) # sticking coefficient = 1.0 def test_modify_BMinterface(self): - gas = ct.Solution('BM_ptcombust.yaml', 'gas') - surf = ct.Interface('BM_ptcombust.yaml', 'Pt_surf', [gas]) + 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 @@ -1419,8 +1421,9 @@ def test_modify_BMinterface(self): self.assertNear(k2, k3) def test_modify_BMsticking(self): - gas = ct.Solution('BM_ptcombust.yaml', 'gas') - surf = ct.Interface('BM_ptcombust.yaml', 'Pt_surf', [gas]) + 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 @@ -1434,8 +1437,9 @@ def test_modify_BMsticking(self): def test_BMmotz_wise(self): # Motz & Wise off for all reactions - gas1 = ct.Solution('BM_ptcombust.yaml', 'gas') - surf1 = ct.Interface('BM_ptcombust.yaml', 'Pt_surf', [gas1]) + 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 diff --git a/data/BM-ptcombust-Motz-Wise.yaml b/test/data/BM-ptcombust-Motz-Wise.yaml similarity index 95% rename from data/BM-ptcombust-Motz-Wise.yaml rename to test/data/BM-ptcombust-Motz-Wise.yaml index ba4fa5acc8f..97aa1c8471d 100644 --- a/data/BM-ptcombust-Motz-Wise.yaml +++ b/test/data/BM-ptcombust-Motz-Wise.yaml @@ -95,21 +95,22 @@ species: reactions: - equation: 2 H(S) => H2 + 2 PT(S) # Reaction 1 - type: surface-Blowers-Masel + 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: surface-Blowers-Masel + type: Blowers-Masel sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 1000000} - equation: H2O + PT(S) => H2O(S) # Reaction 3 - type: surface-Blowers-Masel + 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: surface-Blowers-Masel + 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: surface-Blowers-Masel + type: Blowers-Masel sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 100000} - Motz-Wise: true \ No newline at end of file + Motz-Wise: true + diff --git a/test/data/BM_example.yaml b/test/data/BM_example.yaml deleted file mode 100644 index 8c25789004a..00000000000 --- a/test/data/BM_example.yaml +++ /dev/null @@ -1,86 +0,0 @@ -description: |- - This is the minimal example tweaked from Hydrogen-Oxygen - submechanism extracted from GRI-Mech 3.0. It should only - be used for testing the Blowers Masel reaction class in - gas phase. - -generator: ck2yaml -input-files: [h2o2.inp, gri30_tran.dat] -cantera-version: 2.5.0 -date: Wed, 11 Dec 2019 16:59:04 -0500 - -units: {length: cm, time: s, quantity: mol, activation-energy: cal/mol} - -phases: -- name: ohmech - thermo: ideal-gas - elements: [O, H, Ar, N] - species: [H2, H, O, OH, H2O, O2] - kinetics: gas - state: {T: 300.0, P: 1 atm} -species: -- name: H2 - composition: {H: 2} - thermo: - model: NASA7 - temperature-ranges: [200.0, 1000.0, 3500.0] - data: - - [2.34433112, 7.98052075e-03, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12, - -917.935173, 0.683010238] - - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14, - -950.158922, -3.20502331] -- name: H - composition: {H: 1} - thermo: - model: NASA7 - temperature-ranges: [200.0, 1000.0, 3500.0] - data: - - [2.5, 7.05332819e-13, -1.99591964e-15, 2.30081632e-18, -9.27732332e-22, - 2.54736599e+04, -0.446682853] - - [2.50000001, -2.30842973e-11, 1.61561948e-14, -4.73515235e-18, 4.98197357e-22, - 2.54736599e+04, -0.446682914] -- name: O - composition: {O: 1} - thermo: - model: NASA7 - temperature-ranges: [200.0, 1000.0, 3500.0] - data: - - [3.1682671, -3.27931884e-03, 6.64306396e-06, -6.12806624e-09, 2.11265971e-12, - 2.91222592e+04, 2.05193346] - - [2.56942078, -8.59741137e-05, 4.19484589e-08, -1.00177799e-11, 1.22833691e-15, - 2.92175791e+04, 4.78433864] -- name: OH - composition: {O: 1, H: 1} - thermo: - model: NASA7 - temperature-ranges: [200.0, 1000.0, 3500.0] - data: - - [3.99201543, -2.40131752e-03, 4.61793841e-06, -3.88113333e-09, 1.3641147e-12, - 3615.08056, -0.103925458] - - [3.09288767, 5.48429716e-04, 1.26505228e-07, -8.79461556e-11, 1.17412376e-14, - 3858.657, 4.4766961] -- name: H2O - composition: {H: 2, O: 1} - thermo: - model: NASA7 - temperature-ranges: [200.0, 1000.0, 3500.0] - data: - - [4.19864056, -2.0364341e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12, - -3.02937267e+04, -0.849032208] - - [3.03399249, 2.17691804e-03, -1.64072518e-07, -9.7041987e-11, 1.68200992e-14, - -3.00042971e+04, 4.9667701] -- name: O2 - composition: {O: 2} - thermo: - model: NASA7 - temperature-ranges: [200.0, 1000.0, 3500.0] - data: - - [3.78245636, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, - -1063.94356, 3.65767573] - - [3.28253784, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, - -1088.45772, 5.45323129] - -reactions: -- equation: O + H2 <=> H + OH # Reaction 1 - type: Blowers-Masel - rate-constant: {A: 3.87e+04, b: 2.7, Ea0: 6260.0, w: 1e9} diff --git a/data/BM_ptcombust.yaml b/test/data/BM_test.yaml similarity index 82% rename from data/BM_ptcombust.yaml rename to test/data/BM_test.yaml index 4dfb09695cb..47598ef9378 100644 --- a/data/BM_ptcombust.yaml +++ b/test/data/BM_test.yaml @@ -18,30 +18,26 @@ phases: 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] + - gri30.yaml/species: [H2, H, O, O2, OH, H2O, CH4, AR, N2] skip-undeclared-elements: true kinetics: gas - reactions: - - gri30.yaml/reactions: declared-species + reactions: [gas-reactions] 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)] + species: [{Pt-surf-species: all}] kinetics: surface - reactions: all + 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 -species: +Pt-surf-species: - name: PT(S) composition: {Pt: 1} thermo: @@ -91,21 +87,30 @@ species: - [1.945418, 9.1761647e-04, -1.1226719e-07, -9.9099624e-11, 2.4307699e-14, -1.4005187e+04, -11.531663] -reactions: +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: surface-Blowers-Masel + 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: surface-Blowers-Masel + type: Blowers-Masel sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 1000000} - equation: H2O + PT(S) => H2O(S) # Reaction 3 - type: surface-Blowers-Masel + 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: surface-Blowers-Masel + type: Blowers-Masel sticking-coefficient: {A: 0.023, b: 0, Ea0: 0, w: 1000000} - equation: OH + PT(S) => OH(S) # Reaction 5 - type: surface-Blowers-Masel - sticking-coefficient: {A: 1.0, b: 0, Ea0: 0, w: 100000} \ No newline at end of file + 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/data/surface-phases-BM.yaml b/test/data/surface-phases-BM.yaml deleted file mode 100644 index 2000629e2b8..00000000000 --- a/test/data/surface-phases-BM.yaml +++ /dev/null @@ -1,84 +0,0 @@ -phases: -- name: gas - thermo: ideal-gas - species: [{gas-species: all}] - state: {T: 900 K, P: 1 atm, X: {H2: 0.4, Ar: 0.6}} - -- name: Pt-surf - thermo: ideal-surface - kinetics: surface - reactions: [Pt-reactions] - species: [{Pt-surf-species: all}] - site-density: 2.7063e-9 mol/cm^2 - state: {T: 900 K, P: 1 atm, coverages: {Pt(s): 0.5, H(s): 0.4, O(s): 0.1}} - - -gas-species: -- name: H2 - composition: {H: 2} - thermo: - model: NASA7 - temperature-ranges: [200.0, 1000.0, 3500.0] - data: - - [2.34433112, 0.00798052075, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12, - -917.935173, 0.683010238] - - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14, - -950.158922, -3.20502331] -- name: H - composition: {H: 1} - thermo: - model: NASA7 - temperature-ranges: [200.0, 1000.0, 3500.0] - data: - - [2.5, 7.05332819e-13, -1.99591964e-15, 2.30081632e-18, -9.27732332e-22, - 25473.6599, -0.446682853] - - [2.50000001, -2.30842973e-11, 1.61561948e-14, -4.73515235e-18, 4.98197357e-22, - 25473.6599, -0.446682914] -- name: Ar - composition: {Ar: 1} - thermo: - model: NASA7 - temperature-ranges: [300.0, 5000.0] - data: - - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.366] - - -Pt-surf-species: -- name: Pt(s) - composition: {Pt: 1} - thermo: - model: constant-cp -- name: H(s) - composition: {H: 1, Pt: 1} - thermo: - model: NASA7 - temperature-ranges: [300, 1000, 3000] - data: - - [-1.3029877E+00, 5.4173199E-03, 3.1277972E-07, -3.2328533E-09, - 1.1362820E-12, -4.2277075E+03, 5.8743238E+00] - - [1.0696996E+00, 1.5432230E-03, -1.5500922E-07, -1.6573165E-10, - 3.8359347E-14, -5.0546128E+03, -7.1555238E+00] -- name: O(s) - composition: {O: 1, Pt: 1} - thermo: - model: NASA7 - temperature-ranges: [300, 1000, 3000] - data: - - [-9.4986904E-01, 7.4042305E-03, -1.0451424E-06, -6.1120420E-09, - 3.3787992E-12, -1.3209912E+04, 3.6137905E+00] - - [1.9454180E+00, 9.1761647E-04, -1.1226719E-07, -9.9099624E-11, - 2.4307699E-14, -1.4005187E+04, -1.1531663E+01] - -Pt-reactions: -- equation: H2 + 2 Pt(s) => 2 H(s) - type: surface-Blowers-Masel - rate-constant: {A: 4.4579e10 cm^3/mol/s, b: 0.5, Ea0: 0.0 J/mol, w: 1000000 J/mol} - orders: {H2: 1.0, Pt(s): 1.0} -- equation: 2 H(s) => H2 + 2 Pt(s) - type: surface-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}} -- equation: H + Pt(s) => H(s) - type: surface-Blowers-Masel - sticking-coefficient: [1.0, 0.0, 0.0, 1000000] - diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 35807933f4b..48c01a6432b 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -275,26 +275,26 @@ TEST(Kinetics, InterfaceKineticsFromYaml) TEST(Kinetics, BMInterfaceKineticsFromYaml) { - shared_ptr gas(newPhase("surface-phases-BM.yaml", "gas")); - shared_ptr surf_tp(newPhase("surface-phases-BM.yaml", "Pt-surf")); + 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, "surface-phases-BM.yaml", "Pt-surf"); - EXPECT_EQ(kin->nReactions(), (size_t) 3); - EXPECT_EQ(kin->nTotalSpecies(), (size_t) 6); - auto R1 = kin->reaction(0); - auto IR1 = std::dynamic_pointer_cast(R1); - EXPECT_DOUBLE_EQ(R1->orders["Pt(s)"], 1.0); + 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(1); - auto IR2 = std::dynamic_pointer_cast(R2); + 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_DOUBLE_EQ(IR2->coverage_deps["H(S)"].E, -6e6 / GasConstant); EXPECT_FALSE(IR2->is_sticking_coefficient); - auto R3 = kin->reaction(2); - auto IR3 = std::dynamic_pointer_cast(R3); + auto R3 = kin->reaction(1); + auto IR3 = std::dynamic_pointer_cast(R3); EXPECT_TRUE(IR3->is_sticking_coefficient); } From 6de06d0922c4bcd9aa55cacc937760f46bc56a9f Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 25 Mar 2021 16:07:39 -0400 Subject: [PATCH 06/20] change 'BMInterfaceReaction' to 'BlowersMaselInterfaceReaction' --- include/cantera/kinetics/InterfaceKinetics.h | 6 +++--- include/cantera/kinetics/Reaction.h | 8 ++++---- interfaces/cython/cantera/_cantera.pxd | 2 +- interfaces/cython/cantera/reaction.pyx | 18 +++++++++--------- src/kinetics/InterfaceKinetics.cpp | 8 ++++---- src/kinetics/Reaction.cpp | 8 ++++---- src/kinetics/ReactionFactory.cpp | 4 ++-- 7 files changed, 27 insertions(+), 27 deletions(-) diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index a20c6876dc6..581bb698315 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -19,7 +19,7 @@ namespace Cantera class SurfPhase; class ImplicitSurfChem; class InterfaceReaction; -class BMInterfaceReaction; +class BlowersMaselInterfaceReaction; //! A kinetics manager for heterogeneous reaction mechanisms. The reactions are //! assumed to occur at a 2D interface between two 3D phases. @@ -393,8 +393,8 @@ class InterfaceKinetics : public Kinetics //! @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, BMInterfaceReaction& r, - bool replace); + BMSurfaceArrhenius buildBMSurfaceArrhenius(size_t i, BlowersMaselInterfaceReaction& r, + bool replace); //! Temporary work vector of length m_kk vector_fp m_grt; diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 71fe18fa897..841e85d77de 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -436,11 +436,11 @@ class BlowersMaselReaction: public Reaction //! A reaction occurring on an interface (i.e. a SurfPhase or an EdgePhase) //! with the rate calculated with Blowers-Masel approximation. -class BMInterfaceReaction : public BlowersMaselReaction +class BlowersMaselInterfaceReaction : public BlowersMaselReaction { public: - BMInterfaceReaction(); - BMInterfaceReaction(const Composition& reactants, const Composition& products, + BlowersMaselInterfaceReaction(); + BlowersMaselInterfaceReaction(const Composition& reactants, const Composition& products, const BlowersMasel& rate, bool isStick=false); virtual std::string type() const { @@ -542,7 +542,7 @@ void setupElectrochemicalReaction(ElectrochemicalReaction&, void setupBlowersMaselReaction(BlowersMaselReaction&, const AnyMap&, const Kinetics&); //! @internal May be changed without notice in future versions -void setupBMInterfaceReaction(BMInterfaceReaction&, +void setupBlowersMaselInterfaceReaction(BlowersMaselInterfaceReaction&, const AnyMap&, const Kinetics&); } diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 9ebb1b23cbc..b42facfa691 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -479,7 +479,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cbool use_motz_wise_correction string sticking_species - cdef cppclass CxxBMInterfaceReaction "Cantera::BMInterfaceReaction" (CxxBlowersMaselReaction): + cdef cppclass CxxBlowersMaselInterfaceReaction "Cantera::BlowersMaselInterfaceReaction" (CxxBlowersMaselReaction): stdmap[string, CxxCoverageDependency] coverage_deps cbool is_sticking_coefficient cbool use_motz_wise_correction diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index 053fd6ff60c..36528dfa9a6 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -1162,7 +1162,7 @@ cdef class InterfaceReaction(ElementaryReaction): cdef CxxInterfaceReaction* r = self.reaction r.sticking_species = stringify(species) -cdef class BMInterfaceReaction(BlowersMaselReaction): +cdef class BlowersMaselInterfaceReaction(BlowersMaselReaction): """ A reaction occurring on an `Interface` (i.e. a surface or an edge) """ reaction_type = "surface-Blowers-Masel" @@ -1176,7 +1176,7 @@ cdef class BMInterfaceReaction(BlowersMaselReaction): and the activation energy [J/kmol], respectively. """ def __get__(self): - cdef CxxBMInterfaceReaction* r = self.reaction + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction deps = {} cdef pair[string,CxxCoverageDependency] item for item in r.coverage_deps: @@ -1184,7 +1184,7 @@ cdef class BMInterfaceReaction(BlowersMaselReaction): item.second.E * gas_constant) return deps def __set__(self, deps): - cdef CxxBMInterfaceReaction* r = self.reaction + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction r.coverage_deps.clear() cdef str species for species, D in deps.items(): @@ -1198,10 +1198,10 @@ cdef class BMInterfaceReaction(BlowersMaselReaction): constant. """ def __get__(self): - cdef CxxBMInterfaceReaction* r = self.reaction + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction return r.is_sticking_coefficient def __set__(self, stick): - cdef CxxBMInterfaceReaction* r = self.reaction + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction r.is_sticking_coefficient = stick property use_motz_wise_correction: @@ -1212,10 +1212,10 @@ cdef class BMInterfaceReaction(BlowersMaselReaction): coefficient. """ def __get__(self): - cdef CxxBMInterfaceReaction* r = self.reaction + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction return r.use_motz_wise_correction def __set__(self, mw): - cdef CxxBMInterfaceReaction* r = self.reaction + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction r.use_motz_wise_correction = mw property sticking_species: @@ -1225,8 +1225,8 @@ cdef class BMInterfaceReaction(BlowersMaselReaction): ambiguous. """ def __get__(self): - cdef CxxBMInterfaceReaction* r = self.reaction + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction return pystr(r.sticking_species) def __set__(self, species): - cdef CxxBMInterfaceReaction* r = self.reaction + cdef CxxBlowersMaselInterfaceReaction* r = self.reaction r.sticking_species = stringify(species) diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index eeea03152b8..c95e6309dfb 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -510,7 +510,7 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base) return false; } if (r_base->reaction_type == BMINTERFACE_RXN) { - BMInterfaceReaction& r = dynamic_cast(*r_base); + BlowersMaselInterfaceReaction& r = dynamic_cast(*r_base); BMSurfaceArrhenius rate = buildBMSurfaceArrhenius(i, r, false); m_blowers_masel_rates.install(i, rate); @@ -590,7 +590,7 @@ void InterfaceKinetics::modifyReaction(size_t i, shared_ptr r_base) { Kinetics::modifyReaction(i, r_base); if (r_base->reaction_type == BMINTERFACE_RXN) { - BMInterfaceReaction& r = dynamic_cast(*r_base); + BlowersMaselInterfaceReaction& r = dynamic_cast(*r_base); BMSurfaceArrhenius rate = buildBMSurfaceArrhenius(i, r, true); m_blowers_masel_rates.replace(i, rate); } else { @@ -696,9 +696,9 @@ SurfaceArrhenius InterfaceKinetics::buildSurfaceArrhenius( } BMSurfaceArrhenius InterfaceKinetics::buildBMSurfaceArrhenius( - size_t i, BMInterfaceReaction& r, bool replace) + size_t i, BlowersMaselInterfaceReaction& r, bool replace) { - if (r. is_sticking_coefficient) { + if (r.is_sticking_coefficient) { // Identify the interface phase size_t iInterface = npos; size_t min_dim = 4; diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index bae362990c6..8399ce2d3c5 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -638,14 +638,14 @@ BlowersMaselReaction::BlowersMaselReaction(const Composition& reactants_, reaction_type = BLOWERSMASEL_RXN; } -BMInterfaceReaction::BMInterfaceReaction() +BlowersMaselInterfaceReaction::BlowersMaselInterfaceReaction() : is_sticking_coefficient(false) , use_motz_wise_correction(false) { reaction_type = BMINTERFACE_RXN; } -BMInterfaceReaction::BMInterfaceReaction(const Composition& reactants_, +BlowersMaselInterfaceReaction::BlowersMaselInterfaceReaction(const Composition& reactants_, const Composition& products_, const BlowersMasel& rate_, bool isStick) @@ -1266,7 +1266,7 @@ void setupBlowersMaselReaction(BlowersMaselReaction& R, const AnyMap& node, R.rate = readBlowersMasel(R, node["rate-constant"], kin, node.units()); } -void setupBMInterfaceReaction(BMInterfaceReaction& R, const AnyMap& node, +void setupBlowersMaselInterfaceReaction(BlowersMaselInterfaceReaction& R, const AnyMap& node, const Kinetics& kin) { setupReaction(R, node, kin); @@ -1282,7 +1282,7 @@ void setupBMInterfaceReaction(BMInterfaceReaction& R, const AnyMap& node, kin.thermo().input().getBool("Motz-Wise", false)); R.sticking_species = node.getString("sticking-species", ""); } else { - throw InputFileError("setupBMInterfaceReaction", node, + throw InputFileError("setupBlowersMaselInterfaceReaction", node, "Reaction must include either a 'rate-constant' or" " 'sticking-coefficient' node."); } diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index b920f13fd38..7376c3a30c9 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -152,10 +152,10 @@ ReactionFactory::ReactionFactory() }); // register surface Blowers Masel reactions - reg("surface-Blowers-Masel", []() { return new BMInterfaceReaction(); }); + reg("surface-Blowers-Masel", []() { return new BlowersMaselInterfaceReaction(); }); reg_AnyMap("surface-Blowers-Masel", [](Reaction* R, const AnyMap& node, const Kinetics& kin) { - setupBMInterfaceReaction(*(BMInterfaceReaction*)R, node, kin); + setupBlowersMaselInterfaceReaction(*(BlowersMaselInterfaceReaction*)R, node, kin); }); } From 977ff3fde6d3b3e75d40d6911e8e371101ac8af2 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 25 Mar 2021 16:10:46 -0400 Subject: [PATCH 07/20] make the type of surface Blowers-Masel reaction to be specified without surface in the surface phase in yaml file --- src/kinetics/ReactionFactory.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index 7376c3a30c9..29b3a523a4b 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -228,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; @@ -239,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 { From d12dee78d1a7a3b67566c0e99b4f78c1c07007ad Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 25 Mar 2021 16:11:35 -0400 Subject: [PATCH 08/20] delete memeber variable 'm_E' --- include/cantera/kinetics/RxnRates.h | 19 ++++++++++--------- src/kinetics/RxnRates.cpp | 4 ---- 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 3a6906e42d8..11da0666ef0 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -150,7 +150,7 @@ class BlowersMasel * safely called for negative values of the pre-exponential factor. */ doublereal updateRC(doublereal logT, doublereal recipT, doublereal deltaH) { - m_E = activationEnergy_R(deltaH); + double m_E = activationEnergy_R(deltaH); return m_A * std::exp(m_b * logT - m_E * recipT); } @@ -196,7 +196,7 @@ class BlowersMasel } protected: - doublereal m_logA, m_b, m_E, m_A, m_w, m_E0; + doublereal m_logA, m_b, m_A, m_w, m_E0; }; /** @@ -644,25 +644,26 @@ class BMSurfaceArrhenius //! 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) { - m_E = 0; + E = 0; } else if (deltaH_R > 4 * m_E0) { - m_E = deltaH_R; + 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) - m_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 + 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 m_E + m_ecov; + 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_E + m_ecov; + return m_E0 + m_ecov; } //! Return the bond energy *w* divided by the gas constant[K] @@ -671,7 +672,7 @@ class BMSurfaceArrhenius } protected: - doublereal m_b, m_E, m_A, m_E0, m_w; + 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; diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 320d8fca239..31003414647 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -71,7 +71,6 @@ void Arrhenius::getParameters(AnyMap& rateNode, const Units& rate_units) const BlowersMasel::BlowersMasel() : m_logA(-1.0E300) , m_b(0.0) - , m_E(0.0) , m_A(0.0) , m_w(0.0) , m_E0(0.0) @@ -80,7 +79,6 @@ BlowersMasel::BlowersMasel() BlowersMasel::BlowersMasel(doublereal A, doublereal b, doublereal E0, doublereal w) : m_b(b) - , m_E(E0) , m_A(A) , m_w(w) , m_E0(E0) @@ -224,7 +222,6 @@ ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax, BMSurfaceArrhenius::BMSurfaceArrhenius() : m_b(0.0) - , m_E(0.0) , m_A(0.0) , m_E0(0.0) , m_w(0.0) @@ -236,7 +233,6 @@ BMSurfaceArrhenius::BMSurfaceArrhenius() BMSurfaceArrhenius::BMSurfaceArrhenius(double A, double b, double Ta, double w) : m_b(b) - , m_E(Ta) , m_A(A) , m_E0(Ta) , m_w(w) From f0dc1be37e2a25230ef6cd8ea1d28fd1afd0c9db Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 26 Mar 2021 15:22:37 -0400 Subject: [PATCH 09/20] fix the the syntax error in doc strings to make it rendered correctly by doxygen --- include/cantera/kinetics/RxnRates.h | 35 ++++++++++++++++------------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 11da0666ef0..08b17846194 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -97,24 +97,27 @@ class Arrhenius //! Blowers Masel reaction rate type depends on enthalpy /** - * The Blowers Masel approximation is written by [Paul Blowrs, - * Rich Masel(2004). Engineering Approximations For Activation Energies in Hydrogen - * Transfer Reactions. Eq (10)-(11)] to adjust the activation energy based on - * enthalpy change of a reaction: - * - * \f[ - * Ea = 0 if DeltaH < -4E_0 - * Ea = DeltaH if DeltaH > 4E_0 - * Ea = \frac{(w_0 + DeltaH / 2)(V_P - 2w_0 + DeltaH)^2}{(V_P^^2 - 4w_0^2 + DeltaH^2)} - * \f] + * 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_0 (w_0 + E_0)}{w_0 - E_0} + * V_P = \frac{2w (w_0 + E_0)}{w - E_0}, * \f] - * and \f$ w_0 \f$ is theaverage of the bond energy of the - * bond breaking and that being formed, which can be approximated as - * arbitrary high value like 1000kJ/mol as long as \f$ w_0 >= 2 E_0 \f$ - * The rate constant then can be calculated by Arrhenius function. + * \f$ w \f$ is the average of the bond energy of the + * bond breaking and that being formed. Since the expression is + * very insensitive to \f$ w \f$ for \f$ w >= 2 E_0 \f$, it + * 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] @@ -131,7 +134,7 @@ class BlowersMasel /// order and the dimensionality (surface or bulk). /// @param b Temperature exponent. Non-dimensional. /// @param E0 Intrinsic activation energy in temperature units. Kelvin. - /// @param w bond energy of the bond being formed or broken, in temperature unts. Kelvin. + /// @param w bond energy of the bond being formed or broken, in temperature units. Kelvin. BlowersMasel(doublereal A, doublereal b, doublereal E0, doublereal w); From e42204eb204cd63186d54be7ca64560121519c94 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 26 Mar 2021 15:50:03 -0400 Subject: [PATCH 10/20] add documentation about Blowers Masel in gas and interface phases --- AUTHORS | 1 + doc/sphinx/cython/kinetics.rst | 14 ++++++++ doc/sphinx/yaml/reactions.rst | 60 ++++++++++++++++++++++++++++++---- 3 files changed, 69 insertions(+), 6 deletions(-) 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/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..5d961037262 100644 --- a/doc/sphinx/yaml/reactions.rst +++ b/doc/sphinx/yaml/reactions.rst @@ -26,12 +26,15 @@ 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 +258,31 @@ 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 energy of the bond breaking and that being formed :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 +291,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 +349,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}} From def5e5dddb15b6315ed054a120f794965b11895b Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 29 Mar 2021 00:27:56 -0400 Subject: [PATCH 11/20] Adjust the definition of the average bond energy w and the input data type of Blowers-Masel reaction rates --- doc/sphinx/yaml/reactions.rst | 3 ++- include/cantera/kinetics/RxnRates.h | 40 +++++++++++++++-------------- src/kinetics/RxnRates.cpp | 6 ++--- 3 files changed, 26 insertions(+), 23 deletions(-) diff --git a/doc/sphinx/yaml/reactions.rst b/doc/sphinx/yaml/reactions.rst index 5d961037262..ca007ca0baf 100644 --- a/doc/sphinx/yaml/reactions.rst +++ b/doc/sphinx/yaml/reactions.rst @@ -271,7 +271,8 @@ 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 energy of the bond breaking and that being formed :math:`w`. + 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 diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 08b17846194..d7c6c89ba62 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -109,12 +109,12 @@ class Arrhenius * \f} * where * \f[ - * V_P = \frac{2w (w_0 + E_0)}{w - E_0}, + * V_P = \frac{2w (w + E_0)}{w - E_0}, * \f] - * \f$ w \f$ is the average of the bond energy of the - * bond breaking and that being formed. Since the expression is - * very insensitive to \f$ w \f$ for \f$ w >= 2 E_0 \f$, it - * can be approximated to an arbitrary high value like 1000 kJ/mol. + * \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. @@ -134,9 +134,9 @@ class BlowersMasel /// order and the dimensionality (surface or bulk). /// @param b Temperature exponent. Non-dimensional. /// @param E0 Intrinsic activation energy in temperature units. Kelvin. - /// @param w bond energy of the bond being formed or broken, 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(doublereal A, doublereal b, doublereal E0, doublereal w); + BlowersMasel(double A, double b, double E0, double w); //! Update concentration-dependent parts of the rate coefficient. /*! @@ -193,7 +193,7 @@ class BlowersMasel doublereal activationEnergy_R0() const { return m_E0; } - //! Return the bond energy *w* divided by the gas constant[K] + //! Return the bond dissociation energy *w* divided by the gas constant[K] doublereal bondEnergy() const { return m_w; } @@ -572,18 +572,20 @@ class ChebyshevRate * Transfer Reactions. Eq (10)-(11)], which can be used to change the activation * based on enthalpy of the reaction: * - * \f[ - * Ea = 0 if DeltaH < -4E_0 - * Ea = DeltaH if DeltaH > 4E_0 - * Ea = \frac{(w_0 + DeltaH / 2)(V_P - 2w_0 + DeltaH)^2}{(V_P^2 - 4w_0^2 + DeltaH^2)} - * \f] + * \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_0 (w_0 + E_0) / (w_0 - E_0) + * V_P = 2w (w + E_0) / (w - E_0), * \f] - * and \f$ w_0 \f$ is the average of the bond energy of the - * bond breaking and that being formed, which can be approximated as an - * arbitrary high value like 1000kJ/mol as long as \f$ w_0 >= 2 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 { @@ -595,8 +597,8 @@ class BMSurfaceArrhenius //! 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, doublereal a, - doublereal m, doublereal e); + void addCoverageDependence(size_t k, double a, + double m, double e); void update_C(const doublereal* theta) { m_acov = 0.0; diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 31003414647..4efac3e9a59 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -77,7 +77,7 @@ BlowersMasel::BlowersMasel() { } -BlowersMasel::BlowersMasel(doublereal A, doublereal b, doublereal E0, doublereal w) +BlowersMasel::BlowersMasel(double A, double b, double E0, double w) : m_b(b) , m_A(A) , m_w(w) @@ -242,8 +242,8 @@ BMSurfaceArrhenius::BMSurfaceArrhenius(double A, double b, double Ta, double w) { } -void BMSurfaceArrhenius::addCoverageDependence(size_t k, doublereal a, - doublereal m, doublereal e) +void BMSurfaceArrhenius::addCoverageDependence(size_t k, double a, + double m, double e) { m_sp.push_back(k); m_ac.push_back(a); From a640137a4994e469fce86cf5df9b6ecfc7372552 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 29 Mar 2021 00:28:46 -0400 Subject: [PATCH 12/20] add more docstrings about the Blowers-Masel classes --- doc/sphinx/yaml/reactions.rst | 2 +- include/cantera/kinetics/RxnRates.h | 17 +++++++++-------- interfaces/cython/cantera/reaction.pyx | 19 +++++++++++++++++-- 3 files changed, 27 insertions(+), 11 deletions(-) diff --git a/doc/sphinx/yaml/reactions.rst b/doc/sphinx/yaml/reactions.rst index ca007ca0baf..431c47eb0fe 100644 --- a/doc/sphinx/yaml/reactions.rst +++ b/doc/sphinx/yaml/reactions.rst @@ -271,7 +271,7 @@ 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 + and the average of the bond dissociation energy of the bond breaking and that being formed in the reaction :math:`w`. ``negative-A`` diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index d7c6c89ba62..13ab483d280 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -111,11 +111,11 @@ class Arrhenius * \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 + * \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[ @@ -134,7 +134,8 @@ class BlowersMasel /// 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. + /// @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); @@ -573,17 +574,17 @@ class ChebyshevRate * based on enthalpy of the reaction: * * \f{eqnarray*}{ - * E_a &=& 0\; \text{if }\Delta H < -4E_0 \\ + * 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 + + * 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 + * \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. */ diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index 36528dfa9a6..e5dfbf5fc8d 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -911,6 +911,12 @@ cdef class ChebyshevReaction(Reaction): 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) @@ -937,11 +943,17 @@ cdef class BlowersMasel: def activation_energy(self, float dH): """ - The activation energy *E* [J/kmol]. + 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 @@ -1163,7 +1175,10 @@ cdef class InterfaceReaction(ElementaryReaction): r.sticking_species = stringify(species) cdef class BlowersMaselInterfaceReaction(BlowersMaselReaction): - """ A reaction occurring on an `Interface` (i.e. a surface or an edge) """ + """ + 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: From 757bbc005f4349dbbf72ddcadaa1d1b47d07ab46 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 30 Mar 2021 13:11:10 -0400 Subject: [PATCH 13/20] Fix some formatting/whitespace issues --- doc/sphinx/yaml/reactions.rst | 16 ++++++++------ include/cantera/kinetics/InterfaceKinetics.h | 3 +-- include/cantera/kinetics/Kinetics.h | 2 +- include/cantera/kinetics/Reaction.h | 5 +++-- include/cantera/kinetics/RxnRates.h | 22 +++++++++---------- interfaces/cython/cantera/_cantera.pxd | 5 ++--- interfaces/cython/cantera/reaction.pyx | 6 ++--- .../cython/cantera/test/test_kinetics.py | 18 +++++++-------- src/kinetics/GasKinetics.cpp | 1 - src/kinetics/InterfaceKinetics.cpp | 7 +++--- src/kinetics/ReactionFactory.cpp | 2 +- test/data/BM-ptcombust-Motz-Wise.yaml | 9 ++++---- test/data/BM_test.yaml | 2 +- test/kinetics/kineticsFromYaml.cpp | 6 ++--- 14 files changed, 52 insertions(+), 52 deletions(-) diff --git a/doc/sphinx/yaml/reactions.rst b/doc/sphinx/yaml/reactions.rst index 431c47eb0fe..c3b5ee82d68 100644 --- a/doc/sphinx/yaml/reactions.rst +++ b/doc/sphinx/yaml/reactions.rst @@ -29,11 +29,13 @@ The fields common to all ``reaction`` entries are: - :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 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 `. ``duplicate`` @@ -269,7 +271,7 @@ approximation as `described here `__. -Includes the fields of a :ref:`sec-yaml-Blowers-Masel` reaction and +Includes the fields of a :ref:`sec-yaml-Blowers-Masel` reaction and the fields of an :ref:`sec-yaml-interface-reaction` reaction. Example:: diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 581bb698315..fe74ed39480 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -395,7 +395,7 @@ class InterfaceKinetics : public Kinetics //! future. BMSurfaceArrhenius buildBMSurfaceArrhenius(size_t i, BlowersMaselInterfaceReaction& r, bool replace); - + //! Temporary work vector of length m_kk vector_fp m_grt; @@ -420,7 +420,6 @@ class InterfaceKinetics : public Kinetics */ 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 78ecac9f144..620cb89c792 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -944,7 +944,7 @@ 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; diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 841e85d77de..aad41288360 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -538,12 +538,13 @@ void setupElectrochemicalReaction(ElectrochemicalReaction&, //! @internal May be changed without notice in future versions void setupElectrochemicalReaction(ElectrochemicalReaction&, const AnyMap&, const Kinetics&); + //! @internal May be changed without notice in future versions void setupBlowersMaselReaction(BlowersMaselReaction&, - const AnyMap&, const Kinetics&); + const AnyMap&, const Kinetics&); //! @internal May be changed without notice in future versions void setupBlowersMaselInterfaceReaction(BlowersMaselInterfaceReaction&, - const AnyMap&, const Kinetics&); + const AnyMap&, const Kinetics&); } #endif diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 13ab483d280..5099a3f0980 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -95,16 +95,16 @@ class Arrhenius doublereal m_logA, m_b, m_E, m_A; }; -//! Blowers Masel reaction rate type depends on enthalpy +//! 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 + * 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 &=& 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 + + * 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 @@ -183,7 +183,7 @@ class BlowersMasel // 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) / + 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; @@ -197,7 +197,7 @@ class BlowersMasel //! 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; @@ -646,7 +646,7 @@ class BMSurfaceArrhenius } //! Return the activation energy divided by the gas constant (i.e. the - //! activation temperature) [K] based on the input enthalpy + //! 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) @@ -671,12 +671,12 @@ class BMSurfaceArrhenius 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; diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index b42facfa691..c2ddaaf4dd3 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -450,7 +450,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cdef cppclass CxxTestReaction "Cantera::TestReaction" (CxxReaction2): CxxTestReaction() cbool allow_negative_pre_exponential_factor - + cdef cppclass CxxBlowersMasel "Cantera::BlowersMasel": CxxBlowersMasel() CxxBlowersMasel(double, double, double, double) @@ -461,7 +461,6 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": double activationEnergy_R0() double bondEnergy() - cdef cppclass CxxBlowersMaselReaction "Cantera::BlowersMaselReaction"(CxxReaction): CxxBlowersMaselReaction() CxxBlowersMasel rate @@ -478,7 +477,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": cbool is_sticking_coefficient cbool use_motz_wise_correction string sticking_species - + cdef cppclass CxxBlowersMaselInterfaceReaction "Cantera::BlowersMaselInterfaceReaction" (CxxBlowersMaselReaction): stdmap[string, CxxCoverageDependency] coverage_deps cbool is_sticking_coefficient diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index e5dfbf5fc8d..3db13faf693 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -914,8 +914,8 @@ 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`. + 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: @@ -960,7 +960,7 @@ cdef class BlowersMasel: property intrinsic_activation_energy: """ The intrinsic activation energy *E0* [J/kmol]. - """ + """ def __get__(self): return self.rate.activationEnergy_R0() * gas_constant diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 4202b0d5a8b..61b09f41b19 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -1104,7 +1104,7 @@ def test_BlowersMasel(self): 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]), + self.assertNear(R.rate(gas.T, gas.delta_enthalpy[0]), gas.forward_rate_constants[0], rtol=1e-7) def test_negative_A_blowersmasel(self): @@ -1121,7 +1121,7 @@ def test_negative_A_blowersmasel(self): 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) @@ -1131,10 +1131,10 @@ def test_Blowers_Masel_change_enthalpy(self): 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 / + 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') @@ -1142,7 +1142,7 @@ def test_Blowers_Masel_change_enthalpy(self): 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 = 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)) @@ -1151,7 +1151,7 @@ def test_Blowers_Masel_change_enthalpy(self): 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 = 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)) @@ -1188,13 +1188,13 @@ def test_BlowersMaselinterface(self): 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.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) @@ -1327,7 +1327,7 @@ def test_modify_BlowersMasel(self): 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), + self.assertNear(A1* T ** b1 * np.exp(- Ta1 / T), gas.forward_rate_constants[0]) # randomly modify the rate parameters of a Blowers-Masel reaction diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 6573f5bd7aa..c664215e7ff 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -336,7 +336,6 @@ 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 diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index c95e6309dfb..6946a1a8b84 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -542,11 +542,11 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base) 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; @@ -582,7 +582,8 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base) deltaElectricEnergy_.push_back(0.0); m_deltaG0.push_back(0.0); m_deltaG.push_back(0.0); - m_ProdStanConcReac.push_back(0.0); + m_ProdStanConcReac.push_back(0.0); + return true; } diff --git a/src/kinetics/ReactionFactory.cpp b/src/kinetics/ReactionFactory.cpp index 29b3a523a4b..4ca20d400cf 100644 --- a/src/kinetics/ReactionFactory.cpp +++ b/src/kinetics/ReactionFactory.cpp @@ -143,7 +143,7 @@ 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", diff --git a/test/data/BM-ptcombust-Motz-Wise.yaml b/test/data/BM-ptcombust-Motz-Wise.yaml index 97aa1c8471d..2749a1a2ede 100644 --- a/test/data/BM-ptcombust-Motz-Wise.yaml +++ b/test/data/BM-ptcombust-Motz-Wise.yaml @@ -1,9 +1,9 @@ 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 + 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. + See https://www.detchem.com/mechanisms for original model. Ref:- 1.) Deutschman et al., 26th Symp. (Intl.) on Combustion, 1996 pp. 1747-1754 ----------------------------------------------------------------------- @@ -113,4 +113,3 @@ reactions: 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 index 47598ef9378..665bfb882c5 100644 --- a/test/data/BM_test.yaml +++ b/test/data/BM_test.yaml @@ -2,7 +2,7 @@ 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. + See https://www.detchem.com/mechanisms for original model. Ref:- 1.) Deutschman et al., 26th Symp. (Intl.) on Combustion, 1996 pp. 1747-1754 ----------------------------------------------------------------------- diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 48c01a6432b..42ae28442c3 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -190,7 +190,7 @@ TEST(Reaction, BlowersMaselFromYaml) 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_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)); @@ -200,7 +200,7 @@ TEST(Reaction, BlowersMaselFromYaml) 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_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); @@ -292,7 +292,7 @@ TEST(Kinetics, BMInterfaceKineticsFromYaml) 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); From 770f1e75a67b8ba8d470f19c43f845d11c6f324e Mon Sep 17 00:00:00 2001 From: 12Chao Date: Tue, 6 Apr 2021 15:38:00 -0400 Subject: [PATCH 14/20] Add an example about use case of Blowers-Masel code --- .../examples/kinetics/blowers_masel.py | 67 +++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 interfaces/cython/cantera/examples/kinetics/blowers_masel.py 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..87937ee52ba --- /dev/null +++ b/interfaces/cython/cantera/examples/kinetics/blowers_masel.py @@ -0,0 +1,67 @@ +""" +A simplest example to demonstrate the difference between Blowers-Masel +reaction and elementary reaction. First we set up 2 same reactions with +Arrhenius and Blowers-Masel rate parameters. The Blowers-Masel parameters +are same as Arrhenius parameters with an additional value -- bond energy. +First we show that the forward rate constants of 2 different reactions +are different even if the rate parameters are same. Then we change the +enthalpy of a species involved in all the reactions, and it shows that +the Arrhenius forward rate constant does not change, and the Blowers-Masel +forward rate constants change because the enthalpy change of the reactions +are changed. +""" + +import cantera as ct + +#Create an elementary reaction +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 gas-phase Blowers-Masel reaction +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 +r3 = ct.BlowersMaselReaction({'H':1, 'CH4':1}, {'CH3':1, 'H2':1}) +r3.rate = ct.BlowersMasel(3.87e1, 2.7, 6260*1000*4.184, 1e9) + +gas1 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=ct.Solution('gri30.yaml').species(), reactions=[r1, r2, r3]) + +gas1.TP = 300, ct.one_atm + +print(gas1.forward_rate_constants[0]) + +# Comparing the second and the third rate constants shows that +# even the rate parameters are same, the rate constants are +# different due to the enthalpy change of the reactions are different +print(gas1.forward_rate_constants[1]) +print(gas1.forward_rate_constants[2]) + +def change_species_enthalpy(gas, species_name, dH): + """ + Find the species by name and change it's enthlapy 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) + +#change the enthalpy of species 'H' so that the reaction enthalpy changes +change_species_enthalpy(gas1, 'H', +1e8) + +# the rate constant of elementary should not change with enthalpy of the reaction +print(gas1.forward_rate_constants[0]) + +# the rate constant of Blowers-Masel reaction changes with enthalpy of the reaction +print(gas1.forward_rate_constants[1]) +print(gas1.forward_rate_constants[2]) From 597f6eb546b6570492db2d47f1899d583ca2a6a9 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Sat, 10 Apr 2021 18:18:14 -0400 Subject: [PATCH 15/20] Add 2 plots in the Blowers-Masel example --- .../examples/kinetics/blowers_masel.py | 94 +++++++++++++------ 1 file changed, 67 insertions(+), 27 deletions(-) diff --git a/interfaces/cython/cantera/examples/kinetics/blowers_masel.py b/interfaces/cython/cantera/examples/kinetics/blowers_masel.py index 87937ee52ba..3c50ffa43a3 100644 --- a/interfaces/cython/cantera/examples/kinetics/blowers_masel.py +++ b/interfaces/cython/cantera/examples/kinetics/blowers_masel.py @@ -1,43 +1,75 @@ """ A simplest example to demonstrate the difference between Blowers-Masel -reaction and elementary reaction. First we set up 2 same reactions with -Arrhenius and Blowers-Masel rate parameters. The Blowers-Masel parameters -are same as Arrhenius parameters with an additional value -- bond energy. -First we show that the forward rate constants of 2 different reactions -are different even if the rate parameters are same. Then we change the -enthalpy of a species involved in all the reactions, and it shows that -the Arrhenius forward rate constant does not change, and the Blowers-Masel -forward rate constants change because the enthalpy change of the reactions -are changed. +reaction and elementary reaction. The first two reactions have same reaction +equations with Arrhenius and Blowers-Masel rate parameters, respectively. +The Blowers-Masel parameters are same as 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 2 reactions have same Blowers-Masel parameters can have different +forward rate constants. The first plot generated shows the rate constant changes +with respect to temperature for elementary and Blower-Masel reactions are different. +The second plot shows the activation energy change of a Blowers-Masel reaction +with respect to the delta enthalpy of the reaction. """ import cantera as ct +import numpy as np +import matplotlib.pyplot as plt -#Create an elementary reaction +#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 gas-phase Blowers-Masel reaction +#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) -gas1 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', +gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=ct.Solution('gri30.yaml').species(), reactions=[r1, r2, r3]) -gas1.TP = 300, ct.one_atm +gas.TP = 300, ct.one_atm -print(gas1.forward_rate_constants[0]) +r1_rc = gas.forward_rate_constants[0] +r2_rc = gas.forward_rate_constants[1] +r3_rc = gas.forward_rate_constants[2] -# Comparing the second and the third rate constants shows that -# even the rate parameters are same, the rate constants are -# different due to the enthalpy change of the reactions are different -print(gas1.forward_rate_constants[1]) -print(gas1.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 cosntant 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.scatter(T_range, r1_kf, 1*1, label='Reaction 1 (Elementary)') +plt.scatter(T_range, r2_kf, 1*1, label='Reaction 2 (Blowers-Masel)') +plt.xlabel("Temperature(K)") +plt.ylabel("Forward Rate Constant $(kmol/(m^3.s))$") +plt.title("Comparison Of kf vs. Temperature For Reaction 1 and 2",y=1.1) +plt.legend() +plt.savefig("kf_to_T.png") +plt.clf() +# 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 enthlapy by dH (in J/kmol) @@ -50,18 +82,26 @@ def change_species_enthalpy(gas, species_name, dH): 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) -#change the enthalpy of species 'H' so that the reaction enthalpy changes -change_species_enthalpy(gas1, 'H', +1e8) +# 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 -# the rate constant of elementary should not change with enthalpy of the reaction -print(gas1.forward_rate_constants[0]) +Ea_list = [] +deltaH_list = np.arange(lower_limit_enthalpy, upper_limit_enthalpy, (upper_limit_enthalpy-lower_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])) -# the rate constant of Blowers-Masel reaction changes with enthalpy of the reaction -print(gas1.forward_rate_constants[1]) -print(gas1.forward_rate_constants[2]) +plt.scatter(deltaH_list, Ea_list, 1) +plt.xlabel("Enthalpy Change (J/kmol)") +plt.ylabel("Activation Energy Change (J/kmol)$") +plt.title(r"Ea vs. $\Delta H$ For O+H2<=>H+OH", y=1.1) +plt.savefig("Ea_to_H.png") From 11abb3d86e995de92770e5ca3fd887860a9751d3 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 13 Apr 2021 20:29:39 -0400 Subject: [PATCH 16/20] [Examples/Python] Reformat docstring for Blowers-Masel example Create a shorter leading paragraph which will be shown in the index of examples on the website. Document the dependencies. Modify plot formatting slightly and display them before the script exits. --- .../examples/kinetics/blowers_masel.py | 57 +++++++++++-------- 1 file changed, 32 insertions(+), 25 deletions(-) diff --git a/interfaces/cython/cantera/examples/kinetics/blowers_masel.py b/interfaces/cython/cantera/examples/kinetics/blowers_masel.py index 3c50ffa43a3..90db4174c69 100644 --- a/interfaces/cython/cantera/examples/kinetics/blowers_masel.py +++ b/interfaces/cython/cantera/examples/kinetics/blowers_masel.py @@ -1,16 +1,21 @@ """ -A simplest example to demonstrate the difference between Blowers-Masel -reaction and elementary reaction. The first two reactions have same reaction -equations with Arrhenius and Blowers-Masel rate parameters, respectively. -The Blowers-Masel parameters are same as 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 2 reactions have same Blowers-Masel parameters can have different -forward rate constants. The first plot generated shows the rate constant changes -with respect to temperature for elementary and Blower-Masel reactions are different. -The second plot shows the activation energy change of a Blowers-Masel reaction -with respect to the delta enthalpy of the reaction. +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 @@ -45,7 +50,7 @@ " 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 cosntant of second reaction is {0:.3f} kmol/(m^3.s)," + " 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)) @@ -59,20 +64,20 @@ gas.TP = temp, ct.one_atm r1_kf.append(gas.forward_rate_constants[0]) r2_kf.append(gas.forward_rate_constants[1]) -plt.scatter(T_range, r1_kf, 1*1, label='Reaction 1 (Elementary)') -plt.scatter(T_range, r2_kf, 1*1, label='Reaction 2 (Blowers-Masel)') +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("Forward Rate Constant $(kmol/(m^3.s))$") -plt.title("Comparison Of kf vs. Temperature For Reaction 1 and 2",y=1.1) +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") -plt.clf() + # 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 enthlapy by dH (in J/kmol) + Find the species by name and change it's enthalpy by dH (in J/kmol) """ index = gas.species_index(species_name) @@ -83,25 +88,27 @@ def change_species_enthalpy(gas, species_name, dH): perturbed_coeffs[6] += dx perturbed_coeffs[13] += dx - species.thermo = ct.NasaPoly2(species.thermo.min_temp, species.thermo.max_temp, + 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 +# 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.arange(lower_limit_enthalpy, upper_limit_enthalpy, (upper_limit_enthalpy-lower_limit_enthalpy)/100) +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.scatter(deltaH_list, Ea_list, 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"Ea vs. $\Delta H$ For O+H2<=>H+OH", y=1.1) +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() From da22bc0ec2705c1418709be9e04708043f3e4cf2 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 14 Apr 2021 11:53:38 -0400 Subject: [PATCH 17/20] [Kinetics] Use calculateRateCoeffUnits for sticking coefficients --- include/cantera/kinetics/Reaction.h | 2 ++ src/kinetics/Reaction.cpp | 20 +++++++++++++++++--- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index aad41288360..25781a2ad75 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 { @@ -442,6 +443,7 @@ class BlowersMaselInterfaceReaction : public BlowersMaselReaction BlowersMaselInterfaceReaction(); BlowersMaselInterfaceReaction(const Composition& reactants, const Composition& products, const BlowersMasel& rate, bool isStick=false); + virtual void calculateRateCoeffUnits(const Kinetics& kin); virtual std::string type() const { return "surface-Blowers-Masel"; diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 8399ce2d3c5..7e8c104ab48 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); @@ -656,6 +664,14 @@ BlowersMaselInterfaceReaction::BlowersMaselInterfaceReaction(const Composition& 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 + } +} + Arrhenius readArrhenius(const XML_Node& arrhenius_node) { return Arrhenius(getFloat(arrhenius_node, "A", "toSI"), @@ -909,8 +925,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) @@ -1189,7 +1205,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)); @@ -1276,7 +1291,6 @@ void setupBlowersMaselInterfaceReaction(BlowersMaselInterfaceReaction& R, const R.rate = readBlowersMasel(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 = readBlowersMasel(R, node["sticking-coefficient"], kin, node.units()); R.use_motz_wise_correction = node.getBool("Motz-Wise", kin.thermo().input().getBool("Motz-Wise", false)); From c8474f171cb31d881f3e03a44f4f4120906760e9 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 14 Apr 2021 12:21:32 -0400 Subject: [PATCH 18/20] [Input] Serialize BlowersMasel reactions --- include/cantera/kinetics/Reaction.h | 2 ++ include/cantera/kinetics/RxnRates.h | 2 ++ src/kinetics/Reaction.cpp | 37 +++++++++++++++++++++++++++++ src/kinetics/RxnRates.cpp | 18 ++++++++++++++ test/kinetics/kineticsFromYaml.cpp | 22 +++++++++++++++++ 5 files changed, 81 insertions(+) diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 25781a2ad75..c709b5b2b3e 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -424,6 +424,7 @@ class BlowersMaselReaction: public Reaction 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 { @@ -443,6 +444,7 @@ class BlowersMaselInterfaceReaction : public BlowersMaselReaction 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 { diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 5099a3f0980..5f526ebf220 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -139,6 +139,8 @@ class BlowersMasel 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 diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 7e8c104ab48..ba2899523cd 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -646,6 +646,17 @@ BlowersMaselReaction::BlowersMaselReaction(const Composition& reactants_, 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) @@ -672,6 +683,32 @@ void BlowersMaselInterfaceReaction::calculateRateCoeffUnits(const Kinetics& kin) } } +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"), diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 4efac3e9a59..5571aca3042 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -90,6 +90,24 @@ BlowersMasel::BlowersMasel(double A, double b, double E0, double w) } } +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) diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 42ae28442c3..15584d87871 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -542,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(); +} From 2200f6fdaab9f5d98acef686903e2984b4ab64bd Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 14 Apr 2021 12:27:22 -0400 Subject: [PATCH 19/20] Bump version to 2.6.0a2 --- README.rst | 2 +- SConstruct | 2 +- doc/doxygen/Doxyfile | 2 +- interfaces/cython/cantera/ck2yaml.py | 2 +- interfaces/cython/cantera/cti2yaml.py | 2 +- interfaces/cython/cantera/ctml2yaml.py | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) 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/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), } ) From 0f4bc3b17d07a8e54cfaa26f3eb3bbb971443fa7 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 19 Apr 2021 00:04:39 -0400 Subject: [PATCH 20/20] use assertNear to test the exact ratio of rate constants before and after change the coverage dependencies --- interfaces/cython/cantera/test/test_kinetics.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 61b09f41b19..8c37c2146b9 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -1412,12 +1412,17 @@ def test_modify_BMinterface(self): # 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.assertNotAlmostEqual(k1, k2) + + 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):