[Kinetics] Add option for Motz & Wise correction to sticking reactions

This commit is contained in:
Ray Speth
2016-07-28 17:40:59 -04:00
parent 593ab8a573
commit 1e5eb8c871
9 changed files with 210 additions and 2 deletions

View File

@@ -675,6 +675,7 @@ protected:
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
bool use_motz_wise; //!< 'true' if Motz & Wise correction is being used
};
//! Data for sticking reactions

View File

@@ -217,6 +217,12 @@ public:
//! 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;

View File

@@ -330,6 +330,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera":
cdef cppclass CxxInterfaceReaction "Cantera::InterfaceReaction" (CxxElementaryReaction):
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":

View File

@@ -253,6 +253,17 @@ _valrxn = ''
_valexport = ''
_valfmt = ''
# default for Motz & Wise correction
_motz_wise = None
def enable_motz_wise():
global _motz_wise
_motz_wise = True
def disable_motz_wise():
global _motz_wise
_motz_wise = False
def export_species(filename, fmt = 'CSV'):
global _valexport
global _valfmt
@@ -354,6 +365,8 @@ def write(outName=None):
r = x.addChild('reactionData')
r['id'] = 'reaction_data'
if _motz_wise is not None:
r['motz_wise'] = str(_motz_wise).lower()
for rx in _reactions:
rx.build(r)
@@ -1067,6 +1080,15 @@ class Arrhenius(rate_expression):
addFloat(c, 'e', cov[3], fmt = '%f', defunits = _ue)
class stick(Arrhenius):
def __init__(self, *args, **kwargs):
"""
:param motz_wise: 'True' if the Motz & Wise correction should be used,
'False' if not. If unspecified, use the mechanism default (set using
the functions `enable_motz_wise` or `disable_motz_wise`).
"""
self.motz_wise = kwargs.pop('motz_wise', None)
Arrhenius.__init__(self, *args, **kwargs)
def build(self, p, name=''):
a = p.addChild('Arrhenius')
a['type'] = 'stick'
@@ -1077,6 +1099,8 @@ class stick(Arrhenius):
+ str(ngas) + ': ' + str(self.gas_species))
a['species'] = self.gas_species[0]
if self.motz_wise is not None:
a['motz_wise'] = str(self.motz_wise).lower()
self.unit_factor = 1.0
Arrhenius.build(self, p, name, a)

View File

@@ -704,6 +704,20 @@ cdef class InterfaceReaction(ElementaryReaction):
cdef CxxInterfaceReaction* r = <CxxInterfaceReaction*>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 CxxInterfaceReaction* r = <CxxInterfaceReaction*>self.reaction
return r.use_motz_wise_correction
def __set__(self, mw):
cdef CxxInterfaceReaction* r = <CxxInterfaceReaction*>self.reaction
r.use_motz_wise_correction = mw
property sticking_species:
"""
The name of the sticking species. Needed only for reactions with

View File

@@ -967,3 +967,28 @@ class TestReaction(utilities.CanteraTest):
surf.modify_reaction(2, R)
k2 = surf.forward_rate_constants[2]
self.assertNear(k1, 4*k2)
def test_motz_wise(self):
# Motz & Wise off for all reactions
gas1 = ct.Solution('ptcombust.xml', 'gas')
surf1 = ct.Interface('ptcombust.xml', '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('../data/ptcombust-motzwise.cti', 'gas')
surf2 = ct.Interface('../data/ptcombust-motzwise.cti', '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[2], k2[2]) # sticking coefficient = 1.0
self.assertNear(1.6 * k1[7], k2[7]) # sticking coefficient = 0.75
# M&W toggled off (locally) for reaction 4
self.assertNear(k1[4], k2[4])
# M&W toggled on (locally) for reaction 9
self.assertNear(2.0 * k1[9], k2[9]) # sticking coefficient = 1.0

View File

@@ -791,14 +791,15 @@ SurfaceArrhenius InterfaceKinetics::buildSurfaceArrhenius(
}
if (!replace) {
m_stickingData.emplace_back(
StickData{i, surface_order, multiplier});
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;
}
}
@@ -1007,6 +1008,9 @@ void InterfaceKinetics::applyStickingCorrection(double T, double* kf)
for (size_t n = 0; n < m_stickingData.size(); n++) {
const StickData& item = m_stickingData[n];
if (item.use_motz_wise) {
kf[item.index] /= 1 - 0.5 * kf[item.index];
}
kf[item.index] *= factors[n] * sqrt(T) * item.multiplier;
}
}

View File

@@ -234,6 +234,7 @@ ChebyshevReaction::ChebyshevReaction(const Composition& reactants_,
InterfaceReaction::InterfaceReaction()
: is_sticking_coefficient(false)
, use_motz_wise_correction(false)
{
reaction_type = INTERFACE_RXN;
}
@@ -481,6 +482,19 @@ void setupInterfaceReaction(InterfaceReaction& R, const XML_Node& rxn_node)
if (lowercase(arr["type"]) == "stick") {
R.is_sticking_coefficient = true;
R.sticking_species = arr["species"];
if (lowercase(arr["motz_wise"]) == "true") {
R.use_motz_wise_correction = true;
} else if (lowercase(arr["motz_wise"]) == "false") {
R.use_motz_wise_correction = false;
} else {
// Default value for all reactions
XML_Node* parent = rxn_node.parent();
if (parent && parent->name() == "reactionData"
&& lowercase((*parent)["motz_wise"]) == "true") {
R.use_motz_wise_correction = true;
}
}
}
std::vector<XML_Node*> cov = arr.getChildren("coverage");
for (const auto& node : cov) {

View File

@@ -0,0 +1,119 @@
units(length = "cm", time = "s", quantity = "mol", act_energy = "J/mol")
enable_motz_wise()
ideal_gas(name = "gas",
elements = "O H C N Ar",
species = """gri30: 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""",
transport = 'Mix',
reactions = 'gri30: all',
options = ['skip_undeclared_elements',
'skip_undeclared_species'],
initial_state = state(temperature = 300.0, pressure = OneAtm,
mole_fractions = 'CH4:0.095, O2:0.21, AR:0.79')
)
ideal_interface(name = "Pt_surf",
elements = " Pt H O C ",
species = """ ptcombust: PT(S) H(S) H2O(S) OH(S) CO(S)
CO2(S) CH3(S) CH2(S)s CH(S) C(S) O(S) """,
phases = "gas",
site_density = 2.7063e-9,
reactions = "all",
initial_state = state(temperature = 900.0,
coverages = 'O(S):0.0, PT(S):0.5, H(S):0.5')
)
# Reaction 1
surface_reaction("H2 + 2 PT(S) => 2 H(S)", [4.45790E+10, 0.5, 0],
order = "PT(S):1")
# Reaction 2
surface_reaction( "2 H(S) => H2 + 2 PT(S)",
Arrhenius(3.70000E+21, 0, 67400,
coverage = ['H(S)', 0.0, 0.0, -6000.0]))
# Reaction 3
surface_reaction( "H + PT(S) => H(S)", stick(1.00000E+00, 0, 0))
# Reaction 4
surface_reaction( "O2 + 2 PT(S) => 2 O(S)", Arrhenius(1.80000E+21, -0.5, 0),
options = 'duplicate')
# Reaction 5
surface_reaction( "O2 + 2 PT(S) => 2 O(S)", stick(2.30000E-02, 0, 0, motz_wise=False),
options = 'duplicate')
# Reaction 6
surface_reaction( "2 O(S) => O2 + 2 PT(S)",
Arrhenius(3.70000E+21, 0, 213200,
coverage = ['O(S)', 0.0, 0.0, -60000.0]) )
# Reaction 7
surface_reaction( "O + PT(S) => O(S)", stick(1.00000E+00, 0, 0))
# Reaction 8
surface_reaction( "H2O + PT(S) => H2O(S)", stick(7.50000E-01, 0, 0))
# Reaction 9
surface_reaction( "H2O(S) => H2O + PT(S)", [1.00000E+13, 0, 40300])
# Reaction 10
surface_reaction( "OH + PT(S) => OH(S)", stick(1.00000E+00, 0, 0, motz_wise=True))
# Reaction 11
surface_reaction( "OH(S) => OH + PT(S)", [1.00000E+13, 0, 192800])
# Reaction 12
surface_reaction( "H(S) + O(S) <=> OH(S) + PT(S)", [3.70000E+21, 0, 11500])
# Reaction 13
surface_reaction( "H(S) + OH(S) <=> H2O(S) + PT(S)", [3.70000E+21, 0, 17400])
# Reaction 14
surface_reaction( "OH(S) + OH(S) <=> H2O(S) + O(S)", [3.70000E+21, 0, 48200])
# Reaction 15
surface_reaction( "CO + PT(S) => CO(S)", [1.61800E+20, 0.5, 0],
order = "PT(S):2")
# Reaction 16
surface_reaction( "CO(S) => CO + PT(S)", [1.00000E+13, 0, 125500])
# Reaction 17
surface_reaction( "CO2(S) => CO2 + PT(S)", [1.00000E+13, 0, 20500])
# Reaction 18
surface_reaction( "CO(S) + O(S) => CO2(S) + PT(S)", [3.70000E+21, 0, 105000])
# Reaction 19
surface_reaction( "CH4 + 2 PT(S) => CH3(S) + H(S)", [4.63340E+20, 0.5, 0],
order = "PT(S):2.3")
# Reaction 20
surface_reaction( "CH3(S) + PT(S) => CH2(S)s + H(S)",
[3.70000E+21, 0, 20000])
# Reaction 21
surface_reaction( "CH2(S)s + PT(S) => CH(S) + H(S)", [3.70000E+21, 0, 20000])
# Reaction 22
surface_reaction( "CH(S) + PT(S) => C(S) + H(S)", [3.70000E+21, 0, 20000])
# Reaction 23
surface_reaction( "C(S) + O(S) => CO(S) + PT(S)", [3.70000E+21, 0, 62800])
# Reaction 24
surface_reaction( "CO(S) + PT(S) => C(S) + O(S)", [1.00000E+18, 0, 184000])
# Reaction 25 (12/28/2009 HKM added: This is a fictious rxn that is added for numerical stability.
# The issue is that if multiple surface species have a negative concentration, the
# Jacobian for the surface problem will go singular due to the way negative concentrations
# are truncated within Cantera. Adding in unimolecular desorption rxns with neglibigle real
# effects alleviates the problem.)
surface_reaction( "C(S) => C + PT(S)", [3.7E7, 0, 62800])