GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

[Kinetics] Implement evaluation of "chemically activated" reaction rates

Browse Source
This commit is contained in:
Ray Speth
2013-07-06 21:44:02 +00:00
parent 1cec7cc051
commit 77ccab92bd
6 changed files with 194 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
include/cantera/kinetics/FalloffMgr.h
View File

@@ -44,16 +44,18 @@ public:
/*
* @param rxn Index of the falloff reaction. This will be used to
* determine which array entry is modified in method pr_to_falloff.
* @param type of falloff function to install.
* @param falloffType of falloff function to install.
* @param reactionType Either `FALLOFF_RXN` or `CHEMACT_RXN`
* @param c vector of coefficients for the falloff function.
*/
void install(size_t rxn, int type,
void install(size_t rxn, int falloffType, int reactionType,
const vector_fp& c) {
m_rxn.push_back(rxn);
Falloff* f = m_factory->newFalloff(type,c);
Falloff* f = m_factory->newFalloff(falloffType,c);
m_offset.push_back(m_worksize);
m_worksize += f->workSize();
m_falloff.push_back(f);
m_reactionType.push_back(reactionType);
}
//! Size of the work array required to store intermediate results.
@@ -80,8 +82,15 @@ public:
void pr_to_falloff(doublereal* values, const doublereal* work) {
for (size_t i = 0; i < m_rxn.size(); i++) {
double pr = values[m_rxn[i]];
values[m_rxn[i]] *=
m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr);
if (m_reactionType[i] == FALLOFF_RXN) {
// Pr / (1 + Pr) * F
values[m_rxn[i]] *=
m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr);
} else {
// 1 / (1 + Pr) * F
values[m_rxn[i]] =
m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr);
}
}
}
@@ -92,6 +101,9 @@ protected:
vector_int m_loc;
std::vector<vector_fp::difference_type> m_offset;
size_t m_worksize;
//! Distinguish between falloff and chemically activated reactions
vector_int m_reactionType;
};
}

1
include/cantera/kinetics/reaction_defs.h
View File

@@ -56,7 +56,6 @@ const int CHEBYSHEV_RXN = 6;
* off as the pressure increases, due to collisional stabilization of
* a reaction intermediate. Example: Si + SiH4 (+M) <-> Si2H2 + H2
* (+M), which competes with Si + SiH4 (+M) <-> Si2H4 (+M).
* @todo Implement chemical activation reactions.
*/
const int CHEMACT_RXN = 8;

13
interfaces/cython/cantera/test/test_kinetics.py
View File

@@ -301,3 +301,16 @@ class TestReactionPath(utilities.CanteraTest):
# flux we're looking at
for spec in species:
self.assertTrue(gas.n_atoms(spec, element) > 0)
class TestChemicallyActivated(utilities.CanteraTest):
def test_rate_evaluation(self):
gas = ct.Solution('chemically-activated-reaction.xml')
P = [2026.5, 202650.0, 10132500.0] # pressure
# forward rate of progress, computed using Chemkin
Rf = [2.851022e-04, 2.775924e+00, 2.481792e+03]
for i in range(len(P)):
gas.TPX = 900.0, P[i], [0.01, 0.01, 0.04, 0.10, 0.84]
self.assertNear(gas.forward_rates_of_progress[0], Rf[i], 2e-5)

12
src/kinetics/GasKinetics.cpp
View File

@@ -359,7 +359,11 @@ void GasKinetics::processFalloffReactions()
m_falloffn.pr_to_falloff(&pr[0], work);
for (size_t i = 0; i < m_nfall; i++) {
pr[i] *= m_rfn_high[i];
if (m_rxntype[m_fallindx[i]] == FALLOFF_RXN) {
pr[i] *= m_rfn_high[i];
} else { // CHEMACT_RXN
pr[i] *= m_rfn_low[i];
}
}
scatter_copy(pr.begin(), pr.begin() + m_nfall,
@@ -487,6 +491,7 @@ addReaction(ReactionData& r)
addThreeBodyReaction(r);
break;
case FALLOFF_RXN:
case CHEMACT_RXN:
addFalloffReaction(r);
break;
case PLOG_RXN:
@@ -533,7 +538,8 @@ addFalloffReaction(ReactionData& r)
// install the falloff function calculator for
// this reaction
m_falloffn.install(m_nfall, r.falloffType, r.falloffParameters);
m_falloffn.install(m_nfall, r.falloffType, r.reactionType,
r.falloffParameters);
// forward rxn order equals number of reactants, since rate
// coeff is defined in terms of the high-pressure limit
@@ -541,7 +547,7 @@ addFalloffReaction(ReactionData& r)
// increment the falloff reaction counter
++m_nfall;
registerReaction(reactionNumber(), FALLOFF_RXN, iloc);
registerReaction(reactionNumber(), r.reactionType, iloc);
}
void GasKinetics::

41
src/kinetics/importKinetics.cpp
View File

@@ -253,14 +253,16 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp,
* The Arrhenius expression is
* \f[ k = A T^(b) exp (-E_a / RT). \f]
*/
static void getArrhenius(const XML_Node& node, int& highlow,
static void getArrhenius(const XML_Node& node, int& labeled,
doublereal& A, doublereal& b, doublereal& E)
{
if (node["name"] == "k0") {
highlow = 0;
labeled = -1;
} else if (node["name"] == "kHigh") {
labeled = 1;
} else {
highlow = 1;
labeled = 0;
}
/*
* We parse the children for the A, b, and E components.
@@ -505,24 +507,24 @@ void getRateCoefficient(const XML_Node& kf, Kinetics& kin,
throw CanteraError("getRateCoefficient", "Unknown type: " + type);
}
vector_fp clow(3,0.0), chigh(3,0.0);
vector_fp c_alt(3,0.0), c_base(3,0.0);
for (size_t m = 0; m < kf.nChildren(); m++) {
const XML_Node& c = kf.child(m);
string nm = c.name();
int highlow=0;
int labeled=0;
if (nm == "Arrhenius") {
vector_fp coeff(3);
if (c["type"] == "stick") {
getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]);
chigh = coeff;
c_base = coeff;
} else {
getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]);
if (highlow == 1 || rdata.reactionType == THREE_BODY_RXN
getArrhenius(c, labeled, coeff[0], coeff[1], coeff[2]);
if (labeled == 0 || rdata.reactionType == THREE_BODY_RXN
|| rdata.reactionType == ELEMENTARY_RXN) {
chigh = coeff;
c_base = coeff;
} else {
clow = coeff;
c_alt = coeff;
}
}
if (rdata.reactionType == SURFACE_RXN || rdata.reactionType == EDGE_RXN) {
@@ -536,8 +538,8 @@ void getRateCoefficient(const XML_Node& kf, Kinetics& kin,
}
} else if (nm == "Arrhenius_ExchangeCurrentDensity") {
vector_fp coeff(3);
getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]);
chigh = coeff;
getArrhenius(c, labeled, coeff[0], coeff[1], coeff[2]);
c_base = coeff;
rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
} else if (nm == "falloff") {
getFalloff(c, rdata);
@@ -551,17 +553,16 @@ void getRateCoefficient(const XML_Node& kf, Kinetics& kin,
* Store the coefficients in the ReactionData object for return
* from this function.
*/
if (rdata.reactionType == CHEMACT_RXN) {
rdata.rateCoeffParameters = clow;
if (rdata.reactionType == FALLOFF_RXN) {
rdata.rateCoeffParameters = c_base;
rdata.auxRateCoeffParameters = c_alt;
} else if (rdata.reactionType == CHEMACT_RXN) {
rdata.rateCoeffParameters = c_alt;
rdata.auxRateCoeffParameters = c_base;
} else {
rdata.rateCoeffParameters = chigh;
rdata.rateCoeffParameters = c_base;
}
if (rdata.reactionType == FALLOFF_RXN) {
rdata.auxRateCoeffParameters = clow;
} else if (rdata.reactionType == CHEMACT_RXN) {
rdata.auxRateCoeffParameters = chigh;
}
}
}

134
test/data/chemically-activated-reaction.xml Normal file
View File

@@ -0,0 +1,134 @@
<?xml version="1.0"?>
<ctml>
<validate reactions="yes" species="yes"/>
<!-- phase gas -->
<phase dim="3" id="gas">
<elementArray datasrc="elements.xml">C H N O Ar</elementArray>
<speciesArray datasrc="#species_data">ch3 oh ch2o h2 n2</speciesArray>
<reactionArray datasrc="#reaction_data"/>
<state>
<temperature units="K">300.0</temperature>
<pressure units="Pa">101325.0</pressure>
</state>
<thermo model="IdealGas"/>
<kinetics model="GasKinetics"/>
<transport model="None"/>
</phase>
<!-- species definitions -->
<speciesData id="species_data">
<!-- species ch3 -->
<species name="ch3">
<atomArray>H:3 C:1 </atomArray>
<note>iu0702</note>
<thermo>
<NASA Tmax="1000.0" Tmin="200.0" P0="100000.0">
<floatArray name="coeffs" size="7">
3.657179700E+00, 2.126597900E-03, 5.458388300E-06, -6.618100300E-09,
2.465707400E-12, 1.642271600E+04, 1.673535400E+00</floatArray>
</NASA>
<NASA Tmax="6000.0" Tmin="1000.0" P0="100000.0">
<floatArray name="coeffs" size="7">
2.978120600E+00, 5.797852000E-03, -1.975580000E-06, 3.072979000E-10,
-1.791741600E-14, 1.650951300E+04, 4.722479900E+00</floatArray>
</NASA>
</thermo>
</species>
<!-- species oh -->
<species name="oh">
<atomArray>H:1 O:1 </atomArray>
<note>iu3/03</note>
<thermo>
<NASA Tmax="1000.0" Tmin="200.0" P0="100000.0">
<floatArray name="coeffs" size="7">
3.991984240E+00, -2.401066550E-03, 4.616640330E-06, -3.879163060E-09,
1.363195020E-12, 3.368898360E+03, -1.039984770E-01</floatArray>
</NASA>
<NASA Tmax="6000.0" Tmin="1000.0" P0="100000.0">
<floatArray name="coeffs" size="7">
2.838530330E+00, 1.107412890E-03, -2.940002090E-07, 4.206987290E-11,
-2.422898900E-15, 3.697808080E+03, 5.844946520E+00</floatArray>
</NASA>
</thermo>
</species>
<!-- species ch2o -->
<species name="ch2o">
<atomArray>H:2 C:1 O:1 </atomArray>
<note>g8/88</note>
<thermo>
<NASA Tmax="1000.0" Tmin="200.0" P0="100000.0">
<floatArray name="coeffs" size="7">
4.793723120E+00, -9.908333220E-03, 3.732199900E-05, -3.792852370E-08,
1.317726410E-11, -1.432278790E+04, 6.027980580E-01</floatArray>
</NASA>
<NASA Tmax="6000.0" Tmin="1000.0" P0="100000.0">
<floatArray name="coeffs" size="7">
3.169526650E+00, 6.193205600E-03, -2.250563660E-06, 3.659756600E-10,
-2.201494580E-14, -1.449227560E+04, 6.042078980E+00</floatArray>
</NASA>
</thermo>
</species>
<!-- species h2 -->
<species name="h2">
<atomArray>H:2 </atomArray>
<note>tpis78</note>
<thermo>
<NASA Tmax="1000.0" Tmin="200.0" P0="100000.0">
<floatArray name="coeffs" size="7">
2.344331120E+00, 7.980520750E-03, -1.947815100E-05, 2.015720940E-08,
-7.376117610E-12, -9.179351730E+02, 6.830102380E-01</floatArray>
</NASA>
<NASA Tmax="6000.0" Tmin="1000.0" P0="100000.0">
<floatArray name="coeffs" size="7">
2.932865750E+00, 8.266080260E-04, -1.464023640E-07, 1.541004140E-11,
-6.888048000E-16, -8.130655810E+02, -1.024328650E+00</floatArray>
</NASA>
</thermo>
</species>
<!-- species n2 -->
<species name="n2">
<atomArray>N:2 </atomArray>
<note>g8/02</note>
<thermo>
<NASA Tmax="1000.0" Tmin="200.0" P0="100000.0">
<floatArray name="coeffs" size="7">
3.531005280E+00, -1.236609880E-04, -5.029994330E-07, 2.435306120E-09,
-1.408812350E-12, -1.046976280E+03, 2.967470380E+00</floatArray>
</NASA>
<NASA Tmax="6000.0" Tmin="1000.0" P0="100000.0">
<floatArray name="coeffs" size="7">
2.952576370E+00, 1.396900400E-03, -4.926316030E-07, 7.860101950E-11,
-4.607552040E-15, -9.239486880E+02, 5.871887620E+00</floatArray>
</NASA>
</thermo>
</species>
</speciesData>
<reactionData id="reaction_data">
<!-- reaction 0001 -->
<reaction reversible="yes" type="chemAct" id="0001">
<equation>ch3 + oh (+ M) [=] ch2o + h2 (+ M)</equation>
<rateCoeff>
<Arrhenius>
<A>2.823201E+02</A>
<b>1.46878</b>
<E units="cal/mol">-3270.564950</E>
</Arrhenius>
<Arrhenius name="kHigh">
<A>5.880000E-14</A>
<b>6.721</b>
<E units="cal/mol">-3022.227000</E>
</Arrhenius>
<falloff type="Troe">1.671 434.782 2934.21 3919 </falloff>
</rateCoeff>
<reactants>ch3:1.0 oh:1</reactants>
<products>ch2o:1.0 h2:1</products>
</reaction>
</reactionData>
</ctml>