diff --git a/include/cantera/kinetics/FalloffMgr.h b/include/cantera/kinetics/FalloffMgr.h index 831e8e84d..13c3f71bd 100644 --- a/include/cantera/kinetics/FalloffMgr.h +++ b/include/cantera/kinetics/FalloffMgr.h @@ -44,16 +44,18 @@ public: /* * @param rxn Index of the falloff reaction. This will be used to * determine which array entry is modified in method pr_to_falloff. - * @param type of falloff function to install. + * @param falloffType of falloff function to install. + * @param reactionType Either `FALLOFF_RXN` or `CHEMACT_RXN` * @param c vector of coefficients for the falloff function. */ - void install(size_t rxn, int type, + void install(size_t rxn, int falloffType, int reactionType, const vector_fp& c) { m_rxn.push_back(rxn); - Falloff* f = m_factory->newFalloff(type,c); + Falloff* f = m_factory->newFalloff(falloffType,c); m_offset.push_back(m_worksize); m_worksize += f->workSize(); m_falloff.push_back(f); + m_reactionType.push_back(reactionType); } //! Size of the work array required to store intermediate results. @@ -80,8 +82,15 @@ public: void pr_to_falloff(doublereal* values, const doublereal* work) { for (size_t i = 0; i < m_rxn.size(); i++) { double pr = values[m_rxn[i]]; - values[m_rxn[i]] *= - m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr); + if (m_reactionType[i] == FALLOFF_RXN) { + // Pr / (1 + Pr) * F + values[m_rxn[i]] *= + m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr); + } else { + // 1 / (1 + Pr) * F + values[m_rxn[i]] = + m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr); + } } } @@ -92,6 +101,9 @@ protected: vector_int m_loc; std::vector m_offset; size_t m_worksize; + + //! Distinguish between falloff and chemically activated reactions + vector_int m_reactionType; }; } diff --git a/include/cantera/kinetics/reaction_defs.h b/include/cantera/kinetics/reaction_defs.h index e8c696ac5..62c54d0cc 100644 --- a/include/cantera/kinetics/reaction_defs.h +++ b/include/cantera/kinetics/reaction_defs.h @@ -56,7 +56,6 @@ const int CHEBYSHEV_RXN = 6; * off as the pressure increases, due to collisional stabilization of * a reaction intermediate. Example: Si + SiH4 (+M) <-> Si2H2 + H2 * (+M), which competes with Si + SiH4 (+M) <-> Si2H4 (+M). - * @todo Implement chemical activation reactions. */ const int CHEMACT_RXN = 8; diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 2f88b48f5..43a59d9a3 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -301,3 +301,16 @@ class TestReactionPath(utilities.CanteraTest): # flux we're looking at for spec in species: self.assertTrue(gas.n_atoms(spec, element) > 0) + + +class TestChemicallyActivated(utilities.CanteraTest): + def test_rate_evaluation(self): + gas = ct.Solution('chemically-activated-reaction.xml') + P = [2026.5, 202650.0, 10132500.0] # pressure + + # forward rate of progress, computed using Chemkin + Rf = [2.851022e-04, 2.775924e+00, 2.481792e+03] + + for i in range(len(P)): + gas.TPX = 900.0, P[i], [0.01, 0.01, 0.04, 0.10, 0.84] + self.assertNear(gas.forward_rates_of_progress[0], Rf[i], 2e-5) diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index b87a62f4b..1f27e5288 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -359,7 +359,11 @@ void GasKinetics::processFalloffReactions() m_falloffn.pr_to_falloff(&pr[0], work); for (size_t i = 0; i < m_nfall; i++) { - pr[i] *= m_rfn_high[i]; + if (m_rxntype[m_fallindx[i]] == FALLOFF_RXN) { + pr[i] *= m_rfn_high[i]; + } else { // CHEMACT_RXN + pr[i] *= m_rfn_low[i]; + } } scatter_copy(pr.begin(), pr.begin() + m_nfall, @@ -487,6 +491,7 @@ addReaction(ReactionData& r) addThreeBodyReaction(r); break; case FALLOFF_RXN: + case CHEMACT_RXN: addFalloffReaction(r); break; case PLOG_RXN: @@ -533,7 +538,8 @@ addFalloffReaction(ReactionData& r) // install the falloff function calculator for // this reaction - m_falloffn.install(m_nfall, r.falloffType, r.falloffParameters); + m_falloffn.install(m_nfall, r.falloffType, r.reactionType, + r.falloffParameters); // forward rxn order equals number of reactants, since rate // coeff is defined in terms of the high-pressure limit @@ -541,7 +547,7 @@ addFalloffReaction(ReactionData& r) // increment the falloff reaction counter ++m_nfall; - registerReaction(reactionNumber(), FALLOFF_RXN, iloc); + registerReaction(reactionNumber(), r.reactionType, iloc); } void GasKinetics:: diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index f6c0b9cc0..c67eb0c4e 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -253,14 +253,16 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, * The Arrhenius expression is * \f[ k = A T^(b) exp (-E_a / RT). \f] */ -static void getArrhenius(const XML_Node& node, int& highlow, +static void getArrhenius(const XML_Node& node, int& labeled, doublereal& A, doublereal& b, doublereal& E) { if (node["name"] == "k0") { - highlow = 0; + labeled = -1; + } else if (node["name"] == "kHigh") { + labeled = 1; } else { - highlow = 1; + labeled = 0; } /* * We parse the children for the A, b, and E components. @@ -505,24 +507,24 @@ void getRateCoefficient(const XML_Node& kf, Kinetics& kin, throw CanteraError("getRateCoefficient", "Unknown type: " + type); } - vector_fp clow(3,0.0), chigh(3,0.0); + vector_fp c_alt(3,0.0), c_base(3,0.0); for (size_t m = 0; m < kf.nChildren(); m++) { const XML_Node& c = kf.child(m); string nm = c.name(); - int highlow=0; + int labeled=0; if (nm == "Arrhenius") { vector_fp coeff(3); if (c["type"] == "stick") { getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]); - chigh = coeff; + c_base = coeff; } else { - getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]); - if (highlow == 1 || rdata.reactionType == THREE_BODY_RXN + getArrhenius(c, labeled, coeff[0], coeff[1], coeff[2]); + if (labeled == 0 || rdata.reactionType == THREE_BODY_RXN || rdata.reactionType == ELEMENTARY_RXN) { - chigh = coeff; + c_base = coeff; } else { - clow = coeff; + c_alt = coeff; } } if (rdata.reactionType == SURFACE_RXN || rdata.reactionType == EDGE_RXN) { @@ -536,8 +538,8 @@ void getRateCoefficient(const XML_Node& kf, Kinetics& kin, } } else if (nm == "Arrhenius_ExchangeCurrentDensity") { vector_fp coeff(3); - getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]); - chigh = coeff; + getArrhenius(c, labeled, coeff[0], coeff[1], coeff[2]); + c_base = coeff; rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE; } else if (nm == "falloff") { getFalloff(c, rdata); @@ -551,17 +553,16 @@ void getRateCoefficient(const XML_Node& kf, Kinetics& kin, * Store the coefficients in the ReactionData object for return * from this function. */ - if (rdata.reactionType == CHEMACT_RXN) { - rdata.rateCoeffParameters = clow; + if (rdata.reactionType == FALLOFF_RXN) { + rdata.rateCoeffParameters = c_base; + rdata.auxRateCoeffParameters = c_alt; + } else if (rdata.reactionType == CHEMACT_RXN) { + rdata.rateCoeffParameters = c_alt; + rdata.auxRateCoeffParameters = c_base; } else { - rdata.rateCoeffParameters = chigh; + rdata.rateCoeffParameters = c_base; } - if (rdata.reactionType == FALLOFF_RXN) { - rdata.auxRateCoeffParameters = clow; - } else if (rdata.reactionType == CHEMACT_RXN) { - rdata.auxRateCoeffParameters = chigh; - } } } diff --git a/test/data/chemically-activated-reaction.xml b/test/data/chemically-activated-reaction.xml new file mode 100644 index 000000000..f6d3efa2e --- /dev/null +++ b/test/data/chemically-activated-reaction.xml @@ -0,0 +1,134 @@ + + + + + + + C H N O Ar + ch3 oh ch2o h2 n2 + + + 300.0 + 101325.0 + + + + + + + + + + + + H:3 C:1 + iu0702 + + + + 3.657179700E+00, 2.126597900E-03, 5.458388300E-06, -6.618100300E-09, + 2.465707400E-12, 1.642271600E+04, 1.673535400E+00 + + + + 2.978120600E+00, 5.797852000E-03, -1.975580000E-06, 3.072979000E-10, + -1.791741600E-14, 1.650951300E+04, 4.722479900E+00 + + + + + + + H:1 O:1 + iu3/03 + + + + 3.991984240E+00, -2.401066550E-03, 4.616640330E-06, -3.879163060E-09, + 1.363195020E-12, 3.368898360E+03, -1.039984770E-01 + + + + 2.838530330E+00, 1.107412890E-03, -2.940002090E-07, 4.206987290E-11, + -2.422898900E-15, 3.697808080E+03, 5.844946520E+00 + + + + + + + H:2 C:1 O:1 + g8/88 + + + + 4.793723120E+00, -9.908333220E-03, 3.732199900E-05, -3.792852370E-08, + 1.317726410E-11, -1.432278790E+04, 6.027980580E-01 + + + + 3.169526650E+00, 6.193205600E-03, -2.250563660E-06, 3.659756600E-10, + -2.201494580E-14, -1.449227560E+04, 6.042078980E+00 + + + + + + + H:2 + tpis78 + + + + 2.344331120E+00, 7.980520750E-03, -1.947815100E-05, 2.015720940E-08, + -7.376117610E-12, -9.179351730E+02, 6.830102380E-01 + + + + 2.932865750E+00, 8.266080260E-04, -1.464023640E-07, 1.541004140E-11, + -6.888048000E-16, -8.130655810E+02, -1.024328650E+00 + + + + + + + N:2 + g8/02 + + + + 3.531005280E+00, -1.236609880E-04, -5.029994330E-07, 2.435306120E-09, + -1.408812350E-12, -1.046976280E+03, 2.967470380E+00 + + + + 2.952576370E+00, 1.396900400E-03, -4.926316030E-07, 7.860101950E-11, + -4.607552040E-15, -9.239486880E+02, 5.871887620E+00 + + + + + + + + + ch3 + oh (+ M) [=] ch2o + h2 (+ M) + + + 2.823201E+02 + 1.46878 + -3270.564950 + + + 5.880000E-14 + 6.721 + -3022.227000 + + 1.671 434.782 2934.21 3919 + + ch3:1.0 oh:1 + ch2o:1.0 h2:1 + + +