diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 2cc362045..2ce74d357 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -377,10 +377,11 @@ public: protected: //! Build a SurfaceArrhenius object from a Reaction, taking into account //! the possible sticking coefficient form and coverage dependencies - //! @param i Reaction number. Set to npos if this reaction is replacing - //! an existing rate constant. + //! @param i Reaction number //! @param r Reaction object containing rate coefficient parameters - SurfaceArrhenius buildSurfaceArrhenius(size_t i, InterfaceReaction& r); + //! @param replace True if replacing an existing reaction + SurfaceArrhenius buildSurfaceArrhenius(size_t i, InterfaceReaction& r, + bool replace); //! Temporary work vector of length m_kk vector_fp m_grt; @@ -669,12 +670,17 @@ protected: */ std::vector > m_rxnPhaseIsProduct; - //! Pairs of (reaction index, total order) for sticking reactions, which are - //! needed to compute the dependency of the rate constant on the site - //! density. - std::vector > m_sticking_orders; + //! Values used for converting sticking coefficients into rate constants + struct StickData { + size_t index; //!< index of the sticking reaction in the full reaction list + double order; //!< exponent applied to site density term + double multiplier; //!< multiplicative factor in rate expression + }; - void applyStickingCorrection(double* kf); + //! Data for sticking reactions + std::vector m_stickingData; + + void applyStickingCorrection(double T, double* kf); int m_ioFlag; diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index f163a98dd..13d4a378d 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -139,7 +139,7 @@ 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()); - applyStickingCorrection(m_rfn.data()); + applyStickingCorrection(T, m_rfn.data()); // If we need to do conversions between exchange current density // formulation and regular formulation (either way) do it here. @@ -633,7 +633,7 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base) } InterfaceReaction& r = dynamic_cast(*r_base); - SurfaceArrhenius rate = buildSurfaceArrhenius(i, r); + SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, false); m_rates.install(i, rate); // Turn on the global flag indicating surface coverage dependence @@ -717,7 +717,7 @@ void InterfaceKinetics::modifyReaction(size_t i, shared_ptr r_base) { Kinetics::modifyReaction(i, r_base); InterfaceReaction& r = dynamic_cast(*r_base); - SurfaceArrhenius rate = buildSurfaceArrhenius(npos, r); + SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, true); m_rates.replace(i, rate); // Invalidate cached data @@ -726,11 +726,8 @@ void InterfaceKinetics::modifyReaction(size_t i, shared_ptr r_base) } SurfaceArrhenius InterfaceKinetics::buildSurfaceArrhenius( - size_t i, InterfaceReaction& r) + size_t i, InterfaceReaction& r, bool replace) { - double A_rate = r.rate.preExponentialFactor(); - double b_rate = r.rate.temperatureExponent(); - if (r.is_sticking_coefficient) { // Identify the interface phase size_t iInterface = npos; @@ -742,7 +739,6 @@ SurfaceArrhenius InterfaceKinetics::buildSurfaceArrhenius( } } - b_rate += 0.5; std::string sticking_species = r.sticking_species; if (sticking_species == "") { // Identify the sticking species if not explicitly given @@ -768,6 +764,7 @@ SurfaceArrhenius InterfaceKinetics::buildSurfaceArrhenius( } 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)); @@ -775,7 +772,7 @@ SurfaceArrhenius InterfaceKinetics::buildSurfaceArrhenius( const ThermoPhase& surf = thermo(surfacePhaseIndex()); size_t k = p.speciesIndex(sp.first); if (sp.first == sticking_species) { - A_rate *= sqrt(GasConstant/(2*Pi*p.molecularWeight(k))); + 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 @@ -785,19 +782,32 @@ SurfaceArrhenius InterfaceKinetics::buildSurfaceArrhenius( // site density is known at this time. double order = getValue(r.orders, sp.first, sp.second); if (&p == &surf) { - A_rate *= pow(p.size(k), order); + multiplier *= pow(p.size(k), order); surface_order += order; } else { - A_rate *= pow(p.standardConcentration(k), -order); + multiplier *= pow(p.standardConcentration(k), -order); } } } - if (i != npos) { - m_sticking_orders.emplace_back(i, surface_order); + + if (!replace) { + m_stickingData.emplace_back( + StickData{i, surface_order, multiplier}); + } else { + // Modifying an existing sticking reaction. + for (auto& item : m_stickingData) { + if (item.index == i) { + item.order = surface_order; + item.multiplier = multiplier; + break; + } + } } } - SurfaceArrhenius rate(A_rate, b_rate, r.rate.activationEnergy_R()); + SurfaceArrhenius rate(r.rate.preExponentialFactor(), + r.rate.temperatureExponent(), + r.rate.activationEnergy_R()); // Set up coverage dependencies for (const auto& sp : r.coverage_deps) { @@ -976,9 +986,9 @@ void InterfaceKinetics::determineFwdOrdersBV(ElectrochemicalReaction& r, vector_ } } -void InterfaceKinetics::applyStickingCorrection(double* kf) +void InterfaceKinetics::applyStickingCorrection(double T, double* kf) { - if (m_sticking_orders.empty()) { + if (m_stickingData.empty()) { return; } @@ -989,14 +999,15 @@ void InterfaceKinetics::applyStickingCorrection(double* kf) SurfPhase& surf = dynamic_cast(thermo(reactionPhaseIndex())); double n0 = surf.siteDensity(); if (!cached.validate(n0)) { - factors.resize(m_sticking_orders.size()); - for (size_t n = 0; n < m_sticking_orders.size(); n++) { - factors[n] = pow(n0, -m_sticking_orders[n].second); + factors.resize(m_stickingData.size()); + for (size_t n = 0; n < m_stickingData.size(); n++) { + factors[n] = pow(n0, -m_stickingData[n].order); } } - for (size_t n = 0; n < m_sticking_orders.size(); n++) { - kf[m_sticking_orders[n].first] *= factors[n]; + for (size_t n = 0; n < m_stickingData.size(); n++) { + const StickData& item = m_stickingData[n]; + kf[item.index] *= factors[n] * sqrt(T) * item.multiplier; } }