research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Kinetics] Refactor handling of sticking reactions

Browse source
This commit is contained in:
Ray Speth 2016-07-28 11:52:46 -04:00
parent cc2725ca37
commit 593ab8a573
2 changed files with 46 additions and 29 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
include/cantera/kinetics/InterfaceKinetics.h
View file

@ -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<std::vector<bool> > 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<std::pair<size_t, double> > 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<StickData> m_stickingData;
void applyStickingCorrection(double T, double* kf);
int m_ioFlag;

53
src/kinetics/InterfaceKinetics.cpp
View file

@ -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<Reaction> r_base)
}
InterfaceReaction& r = dynamic_cast<InterfaceReaction&>(*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<Reaction> r_base)
{
Kinetics::modifyReaction(i, r_base);
InterfaceReaction& r = dynamic_cast<InterfaceReaction&>(*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<Reaction> 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<SurfPhase&>(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;
}
}