From 5958578c401a51a952e439301e8b7d087f074b07 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 15 Nov 2014 00:47:25 +0000 Subject: [PATCH] [Kinetics] Check for negative and non-reactant reaction orders 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. --- doc/sphinx/cti/reactions.rst | 8 ++++ include/cantera/kinetics/Reaction.h | 18 +++++--- interfaces/cython/cantera/ctml_writer.py | 26 ++++++------ src/kinetics/Kinetics.cpp | 2 +- src/kinetics/Reaction.cpp | 47 ++++++++++++++++++--- test/data/frac.cti | 4 ++ test/data/frac.xml | 16 +++++++ test/kinetics/kineticsFromScratch.cpp | 53 ++++++++++++++++++++++++ test/kinetics/rates.cpp | 6 +-- 9 files changed, 154 insertions(+), 26 deletions(-) diff --git a/doc/sphinx/cti/reactions.rst b/doc/sphinx/cti/reactions.rst index f713f3fb9..6bd361f77 100644 --- a/doc/sphinx/cti/reactions.rst +++ b/doc/sphinx/cti/reactions.rst @@ -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 ======================================== diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 38ba218e3..b9510c88b 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -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; }; diff --git a/interfaces/cython/cantera/ctml_writer.py b/interfaces/cython/cantera/ctml_writer.py index f0bd8bf84..47ee358b0 100644 --- a/interfaces/cython/cantera/ctml_writer.py +++ b/interfaces/cython/cantera/ctml_writer.py @@ -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' diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 339853b8c..96a68b429 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -655,7 +655,7 @@ void Kinetics::addReaction(ReactionData& r) { void Kinetics::addReaction(shared_ptr 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 diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index fd1acc62f..e98470a42 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -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; diff --git a/test/data/frac.cti b/test/data/frac.cti index a338a42b3..e78717965 100644 --- a/test/data/frac.cti +++ b/test/data/frac.cti @@ -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') diff --git a/test/data/frac.xml b/test/data/frac.xml index 8b09e3c3c..53a0b40b5 100644 --- a/test/data/frac.xml +++ b/test/data/frac.xml @@ -153,5 +153,21 @@ H2:0.69999999999999996 O2:0.20000000000000001 OH:0.59999999999999998 H2O:1.0 + + + H2 + 0.5 O2 =] H2O + 1.0 + -0.25 + + + 3.981072E+00 + 0.0 + 0.000000 + + + H2:1.0 O2:0.5 + H2O:1.0 + + diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index ceacc2078..bcf1f789e 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -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 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 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 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 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 { diff --git a/test/kinetics/rates.cpp b/test/kinetics/rates.cpp index 7256d4967..384d8a839 100644 --- a/test/kinetics/rates.cpp +++ b/test/kinetics/rates.cpp @@ -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")]);