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

[Kinetics] Check for negative and non-reactant reaction orders

Browse source 
Allow non-reactant orders for electrochemical reactions

Allow negative orders specifically requested, e.g. by setting the
'negative_orders' option in the CTI definition of the reaction.
This commit is contained in:
Ray Speth 2014-11-15 00:47:25 +00:00
parent e09db92756
commit 5958578c40
9 changed files with 154 additions and 26 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
doc/sphinx/cti/reactions.rst
View file

@ -92,6 +92,8 @@ Note that the ID string is only used when selectively importing reactions. If
all reactions in the local file or in an external one are imported into a phase
or interface, then the reaction ``ID`` field is not used.
.. _sec-reaction-options:
Options
-------
@ -136,6 +138,12 @@ should be handled.
positive, then negative *A* parameters are acceptable, as long as the
``'negative_A'`` option is specified.
``negative_orders``
Reaction orders are normally required to be non-negative, since negative
orders are non-physical and undefined at zero concentration. Cantera allows
negative orders for a global reaction only if the ``negative_orders``
override option is specified for the reaction.
Reactions with Pressure-Independent Rate
========================================

18
include/cantera/kinetics/Reaction.h
View file

@ -28,8 +28,9 @@ public:
virtual std::string productString() const;
std::string equation() const;
//! Ensure that the rate constant for this reaction is valid.
virtual void validateRateConstant() {}
//! Ensure that the rate constant and other parameters for this reaction are
//valid.
virtual void validate();
//! Type of the reaction. The valid types are listed in the file,
//! reaction_defs.h, with constants ending in `RXN`.
@ -55,6 +56,13 @@ public:
//! True if the current reaction is marked as duplicate
bool duplicate;
//! True if reaction orders can be specified for non-reactant species.
//Default is `false`.
bool allow_nonreactant_orders;
//! True if negative reaction orders are allowed. Default is `false`.
bool allow_negative_orders;
};
@ -66,7 +74,7 @@ public:
ElementaryReaction();
ElementaryReaction(const Composition& reactants, const Composition products,
const Arrhenius& rate);
virtual void validateRateConstant();
virtual void validate();
Arrhenius rate;
bool allow_negative_pre_exponential_factor;
@ -115,7 +123,7 @@ public:
const vector_fp& falloff_params);
virtual std::string reactantString() const;
virtual std::string productString() const;
virtual void validateRateConstant();
virtual void validate();
Arrhenius low_rate;
Arrhenius high_rate;
@ -148,7 +156,7 @@ public:
PlogReaction();
PlogReaction(const Composition& reactants, const Composition& products,
const Plog& rate);
virtual void validateRateConstant();
virtual void validate();
Plog rate;
};

26
interfaces/cython/cantera/ctml_writer.py
View file

@ -1128,8 +1128,9 @@ class reaction(object):
An optional identification string. If omitted, it defaults to a
four-digit numeric string beginning with 0001 for the first
reaction in the file.
:param options:
Processing options, as described in :ref:`sec-phase-options`.
:param options: Processing options, as described in
:ref:`sec-reaction-options`. May be one or more (as a list) of the
following: 'skip', 'duplicate', 'negative_A', 'negative_orders'.
"""
self._id = id
self._e = equation
@ -1218,11 +1219,12 @@ class reaction(object):
else:
r['reversible'] = 'no'
for s in self._options:
if s == 'duplicate':
r['duplicate'] = 'yes'
elif s == 'negative_A':
r['negative_A'] = 'yes'
if 'duplicate' in self._options:
r['duplicate'] = 'yes'
if 'negative_A' in self._options:
r['negative_A'] = 'yes'
if 'negative_orders' in self._options:
r['negative_orders'] = 'yes'
ee = self._e.replace('<','[').replace('>',']')
r.addChild('equation',ee)
@ -1337,8 +1339,8 @@ class three_body_reaction(reaction):
An optional identification string. If omitted, it defaults to a
four-digit numeric string beginning with 0001 for the first
reaction in the file.
:param options:
Processing options, as described in :ref:`sec-phase-options`.
:param options: Processing options, as described in
:ref:`sec-reaction-options`.
"""
reaction.__init__(self, equation, kf, id, '', options)
self._type = 'threeBody'
@ -1426,7 +1428,7 @@ class falloff_reaction(pdep_reaction):
four-digit numeric string beginning with 0001 for the first
reaction in the file.
:param options:
Processing options, as described in :ref:`sec-phase-options`.
Processing options, as described in :ref:`sec-reaction-options`.
"""
kf2 = (kf, kf0)
reaction.__init__(self, equation, kf2, id, '', options)
@ -1469,7 +1471,7 @@ class chemically_activated_reaction(pdep_reaction):
four-digit numeric string beginning with 0001 for the first
reaction in the file.
:param options:
Processing options, as described in :ref:`sec-phase-options`.
Processing options, as described in :ref:`sec-reaction-options`.
"""
reaction.__init__(self, equation, (kLow, kHigh), id, '', options)
self._type = 'chemAct'
@ -1599,7 +1601,7 @@ class surface_reaction(reaction):
four-digit numeric string beginning with 0001 for the first
reaction in the file.
:param options:
Processing options, as described in :ref:`sec-phase-options`.
Processing options, as described in :ref:`sec-reaction-options`.
"""
reaction.__init__(self, equation, kf, id, order, options)
self._type = 'surface'

2
src/kinetics/Kinetics.cpp
View file

@ -655,7 +655,7 @@ void Kinetics::addReaction(ReactionData& r) {
void Kinetics::addReaction(shared_ptr<Reaction> r)
{
r->validateRateConstant();
r->validate();
// If reaction orders are specified, then this reaction does not follow
// mass-action kinetics, and is not an elementary reaction. So check that it

47
src/kinetics/Reaction.cpp
View file

@ -16,6 +16,8 @@ Reaction::Reaction(int type)
: reaction_type(type)
, reversible(true)
, duplicate(false)
, allow_nonreactant_orders(false)
, allow_negative_orders(false)
{
}
@ -26,9 +28,36 @@ Reaction::Reaction(int type, const Composition& reactants_,
, products(products_)
, reversible(true)
, duplicate(false)
, allow_nonreactant_orders(false)
, allow_negative_orders(false)
{
}
void Reaction::validate()
{
if (!allow_nonreactant_orders) {
for (Composition::iterator iter = orders.begin();
iter != orders.end();
++iter) {
if (reactants.find(iter->first) == reactants.end()) {
throw CanteraError("Reaction::validate", "Reaction order "
"specified for non-reactant species '" + iter->first + "'");
}
}
}
if (!allow_negative_orders) {
for (Composition::iterator iter = orders.begin();
iter != orders.end();
++iter) {
if (iter->second < 0.0) {
throw CanteraError("Reaction::validate", "Negative reaction "
"order specified for species '" + iter->first + "'");
}
}
}
}
std::string Reaction::reactantString() const
{
std::ostringstream result;
@ -87,11 +116,12 @@ ElementaryReaction::ElementaryReaction()
{
}
void ElementaryReaction::validateRateConstant()
void ElementaryReaction::validate()
{
Reaction::validate();
if (!allow_negative_pre_exponential_factor &&
rate.preExponentialFactor() < 0) {
throw CanteraError("ElementaryReaction::validateRateConstant",
throw CanteraError("ElementaryReaction::validate",
"Undeclared negative pre-exponential factor found in reaction '"
+ equation() + "'");
}
@ -165,10 +195,11 @@ std::string FalloffReaction::productString() const {
}
}
void FalloffReaction::validateRateConstant() {
void FalloffReaction::validate() {
Reaction::validate();
if (low_rate.preExponentialFactor() < 0 ||
high_rate.preExponentialFactor() < 0) {
throw CanteraError("FalloffReaction::validateRateConstant", "Negative "
throw CanteraError("FalloffReaction::validate", "Negative "
"pre-exponential factor found for reaction '" + equation() + "'");
}
}
@ -342,6 +373,9 @@ void setupElementaryReaction(ElementaryReaction& R, const XML_Node& rxn_node)
if (rxn_node["negative_A"] == "yes") {
R.allow_negative_pre_exponential_factor = true;
}
if (rxn_node["negative_orders"] == "yes") {
R.allow_negative_orders = true;
}
setupReaction(R, rxn_node);
}
@ -421,8 +455,9 @@ void setupPlogReaction(PlogReaction& R, const XML_Node& rxn_node)
setupReaction(R, rxn_node);
}
void PlogReaction::validateRateConstant()
void PlogReaction::validate()
{
Reaction::validate();
rate.validate(equation());
}
@ -508,6 +543,7 @@ void setupElectrochemicalReaction(ElectrochemicalReaction& R,
R.orders.clear();
// Reaction orders based on species stoichiometric coefficients
R.allow_nonreactant_orders = true;
for (Composition::const_iterator iter = R.reactants.begin();
iter != R.reactants.end();
++iter) {
@ -524,6 +560,7 @@ void setupElectrochemicalReaction(ElectrochemicalReaction& R,
if (rxn_node.hasChild("reactionOrderFormulation")) {
Composition initial_orders = R.orders;
R.orders.clear();
R.allow_nonreactant_orders = true;
const XML_Node& rof_node = rxn_node.child("reactionOrderFormulation");
if (lowercase(rof_node["model"]) == "reactantorders") {
R.orders = initial_orders;

4
test/data/frac.cti
View file

@ -97,3 +97,7 @@ reaction( "H2O => 1.4 H + 0.6 OH + 0.2 O2", [1.0e13, 0.0, 0.0])
# coefficients.
reaction( "0.7 H2 + 0.6 OH + 0.2 O2 => H2O", [1.0e13, 0.0, 0.0],
order = "H2:0.8 OH:2 O2:1")
# A reaction with negative reaction orders
reaction( "H2 + 0.5 O2 => H2O", [1.0e9, 0.0, 0.0],
order = "H2:1.0 O2:-0.25", options='negative_orders')

16
test/data/frac.xml
View file

@ -153,5 +153,21 @@
<reactants>H2:0.69999999999999996 O2:0.20000000000000001 OH:0.59999999999999998</reactants>
<products>H2O:1.0</products>
</reaction>
<reaction reversible="no" id="0002" negative_orders="yes">
<equation>H2 + 0.5 O2 =] H2O</equation>
<order species="H2">1.0</order>
<order species="O2">-0.25</order>
<rateCoeff>
<Arrhenius>
<A>3.981072E+00</A>
<b>0.0</b>
<E units="cal/mol">0.000000</E>
</Arrhenius>
</rateCoeff>
<reactants>H2:1.0 O2:0.5</reactants>
<products>H2O:1.0</products>
</reaction>
</reactionData>
</ctml>

53
test/kinetics/kineticsFromScratch.cpp
View file

@ -261,6 +261,59 @@ TEST_F(KineticsFromScratch, invalid_reversible_with_orders)
ASSERT_EQ(0, kin.nReactions());
}
TEST_F(KineticsFromScratch, negative_order_override)
{
Composition reac = parseCompString("O:1 H2:1");
Composition prod = parseCompString("H:1 OH:1");
Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
shared_ptr<ElementaryReaction> R(new ElementaryReaction(reac, prod, rate));
R->reversible = false;
R->allow_negative_orders = true;
R->orders["H2"] = - 0.5;
kin.addReaction(R);
ASSERT_EQ((size_t) 1, kin.nReactions());
}
TEST_F(KineticsFromScratch, invalid_negative_orders)
{
Composition reac = parseCompString("O:1 H2:1");
Composition prod = parseCompString("H:1 OH:1");
Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
shared_ptr<ElementaryReaction> R(new ElementaryReaction(reac, prod, rate));
R->reversible = false;
R->orders["H2"] = - 0.5;
ASSERT_THROW(kin.addReaction(R), CanteraError);
ASSERT_EQ(0, kin.nReactions());
}
TEST_F(KineticsFromScratch, nonreactant_order_override)
{
Composition reac = parseCompString("O:1 H2:1");
Composition prod = parseCompString("H:1 OH:1");
Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
shared_ptr<ElementaryReaction> R(new ElementaryReaction(reac, prod, rate));
R->reversible = false;
R->allow_nonreactant_orders = true;
R->orders["OH"] = 0.5;
kin.addReaction(R);
ASSERT_EQ((size_t) 1, kin.nReactions());
}
TEST_F(KineticsFromScratch, invalid_nonreactant_order)
{
Composition reac = parseCompString("O:1 H2:1");
Composition prod = parseCompString("H:1 OH:1");
Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K);
shared_ptr<ElementaryReaction> R(new ElementaryReaction(reac, prod, rate));
R->reversible = false;
R->orders["OH"] = 0.5;
ASSERT_THROW(kin.addReaction(R), CanteraError);
ASSERT_EQ(0, kin.nReactions());
}
class InterfaceKineticsFromScratch : public testing::Test
{

6
test/kinetics/rates.cpp
View file

@ -116,10 +116,10 @@ TEST_F(FracCoeffTest, CreationDestructionRates)
EXPECT_DOUBLE_EQ(0.6*ropf[0], cdot[kOH]);
EXPECT_DOUBLE_EQ(0.2*ropf[0], cdot[kO2]);
EXPECT_DOUBLE_EQ(0.7*ropf[1], ddot[kH2]);
EXPECT_DOUBLE_EQ(0.7*ropf[1]+ropf[2], ddot[kH2]);
EXPECT_DOUBLE_EQ(0.6*ropf[1], ddot[kOH]);
EXPECT_DOUBLE_EQ(0.2*ropf[1], ddot[kO2]);
EXPECT_DOUBLE_EQ(ropf[1], cdot[kH2O]);
EXPECT_DOUBLE_EQ(0.2*ropf[1]+0.5*ropf[2], ddot[kO2]);
EXPECT_DOUBLE_EQ(ropf[1]+ropf[2], cdot[kH2O]);
EXPECT_DOUBLE_EQ(0.0, cdot[therm.speciesIndex("O")]);
EXPECT_DOUBLE_EQ(0.0, ddot[therm.speciesIndex("O")]);