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

Implemented parsing of P-log and Chebyshev reactions from CTML

Browse source
This commit is contained in:
Ray Speth 2012-03-30 23:47:31 +00:00
parent 8386df50b9
commit b4c07978e6
5 changed files with 156 additions and 97 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
include/cantera/kinetics/GasKinetics.h
View file

@ -452,6 +452,8 @@ private:
void addElementaryReaction(const ReactionData& r);
void addThreeBodyReaction(const ReactionData& r);
void addFalloffReaction(const ReactionData& r);
void addPlogReaction(const ReactionData& r);
void addChebyshevReaction(const ReactionData& r);
void installReagents(const ReactionData& r);

28
include/cantera/kinetics/reaction_defs.h
View file

@ -37,6 +37,20 @@ const int THREE_BODY_RXN = 2;
*/
const int FALLOFF_RXN = 4;
/**
* A pressure-dependent rate expression consisting of several Arrhenius rate
* expressions evaluated at different pressures. The final rate is calculated
* by logarithmically interpolating between the two rates that bracket the
* current pressure.
*/
const int PLOG_RXN = 5;
/**
* A general pressure-dependent reaction where k(T,P) is defined in terms of
* a bivariate Chebyshev polynomial.
*/
const int CHEBYSHEV_RXN = 6;
/**
* A chemical activation reaction. For these reactions, the rate falls
* off as the pressure increases, due to collisional stabilization of
@ -75,12 +89,14 @@ const int GLOBAL_RXN = 30;
*/
//@{
const int ARRHENIUS_REACTION_RATECOEFF_TYPE = 1;
const int LANDAUTELLER_REACTION_RATECOEFF_TYPE = 2;
const int TSTRATE_REACTION_RATECOEFF_TYPE = 3;
const int SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE = 4;
const int ARRHENIUS_SUM_REACTION_RATECOEFF_TYPE = 5;
const int EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE = 6;
const int ARRHENIUS_REACTION_RATECOEFF_TYPE = 1;
const int LANDAUTELLER_REACTION_RATECOEFF_TYPE = 2;
const int TSTRATE_REACTION_RATECOEFF_TYPE = 3;
const int SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE = 4;
const int ARRHENIUS_SUM_REACTION_RATECOEFF_TYPE = 5;
const int EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE = 6;
const int PLOG_REACTION_RATECOEFF_TYPE = 7;
const int CHEBYSHEV_REACTION_RATECOEFF_TYPE = 8;
//@}

29
src/kinetics/GasKinetics.cpp
View file

@ -657,13 +657,24 @@ getRevRateConstants(doublereal* krev, bool doIrreversible)
void GasKinetics::
addReaction(const ReactionData& r)
{
if (r.reactionType == ELEMENTARY_RXN) {
switch (r.reactionType) {
case ELEMENTARY_RXN:
addElementaryReaction(r);
} else if (r.reactionType == THREE_BODY_RXN) {
break;
case THREE_BODY_RXN:
addThreeBodyReaction(r);
} else if (r.reactionType == FALLOFF_RXN) {
break;
case FALLOFF_RXN:
addFalloffReaction(r);
break;
case PLOG_RXN:
addPlogReaction(r);
break;
case CHEBYSHEV_RXN:
addChebyshevReaction(r);
break;
default:
throw CanteraError("GasKinetics::addReaction", "Invalid reaction type specified");
}
// operations common to all reaction types
@ -761,6 +772,16 @@ addThreeBodyReaction(const ReactionData& r)
}
//====================================================================================================================
void GasKinetics::addPlogReaction(const ReactionData& r)
{
// @todo: Not yet implemented
}
void GasKinetics::addChebyshevReaction(const ReactionData& r)
{
// @todo: Not yet implemented
}
void GasKinetics::installReagents(const ReactionData& r)
{

15
src/kinetics/ReactionData.h
View file

@ -9,6 +9,7 @@
#define CT_REACTION_DATA_H
#include "cantera/kinetics/reaction_defs.h"
#include "cantera/kinetics/RxnRates.h"
namespace Cantera
{
@ -71,6 +72,20 @@ public:
bool global;
bool isReversibleWithFrac;
doublereal beta; // for electrochemical reactions
//! Arrhenius parameters for P-log reactions.
//! The keys are the pressures corresponding to each Arrhenius expression.
std::map<double, Arrhenius> plogParameters;
double chebTmin; //!< Minimum temperature for Chebyshev fit
double chebTmax; //!< Maximum temperature for Chebyshev fit
double chebPmin; //!< Minimum pressure for Chebyshev fit
double chebPmax; //!< Maximum pressure for Chebyshev fit
size_t chebDegreeT; //!< Degree of Chebyshev fit in T
size_t chebDegreeP; //!< Degree of Chebyshev fit in P
//! Chebyshev coefficients. length chebDegreeT * chebDegreeP
vector_fp chebCoeffs;
};
}

179
src/kinetics/importKinetics.cpp
View file

@ -68,19 +68,6 @@ public:
bool validate_rxn) ;
};
/*
* First we define a couple of typedefs that will
* be used throught this file
*/
//! typedef for a pointer to an XML_Node
typedef const vector<XML_Node*> nodeset_t;
//! typedef for an XML_Node
typedef XML_Node node_t;
/*
* Check a reaction to see if the elements balance.
*/
@ -395,7 +382,7 @@ static void getStick(const XML_Node& node, Kinetics& kin,
E /= GasConstant;
}
static void getCoverageDependence(const node_t& node,
static void getCoverageDependence(const XML_Node& node,
thermo_t& surfphase, ReactionData& rdata)
{
vector<XML_Node*> cov;
@ -428,7 +415,7 @@ static void getCoverageDependence(const node_t& node,
<falloff type="Troe"> 0.5 73.2 5000. 9999. </falloff>
@endverbatim
*/
static void getFalloff(const node_t& f, ReactionData& rdata)
static void getFalloff(const XML_Node& f, ReactionData& rdata)
{
string type = f["type"];
vector<string> p;
@ -489,7 +476,7 @@ static void getFalloff(const node_t& f, ReactionData& rdata)
* reaction mechanism is homogeneous, so that all species belong
* to phase(0) of 'kin'.
*/
static void getEfficiencies(const node_t& eff, Kinetics& kin, ReactionData& rdata)
static void getEfficiencies(const XML_Node& eff, Kinetics& kin, ReactionData& rdata)
{
// set the default collision efficiency
@ -514,89 +501,106 @@ static void getEfficiencies(const node_t& eff, Kinetics& kin, ReactionData& rdat
*
* @param kf Reference to the XML Node named rateCoeff
*/
void getRateCoefficient(const node_t& kf, Kinetics& kin,
void getRateCoefficient(const XML_Node& kf, Kinetics& kin,
ReactionData& rdata, int negA)
{
string type = kf.attrib("type");
if (type == "") {
type = "Arrhenius";
rdata.rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE;
}
if (type == "ExchangeCurrentDensity") {
rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
} else if (type == "Arrhenius") {
if (rdata.reactionType == PLOG_RXN) {
rdata.rateCoeffType = PLOG_REACTION_RATECOEFF_TYPE;
for (size_t m = 0; m < kf.nChildren(); m++) {
const XML_Node& node = kf.child(m);
double A = getFloat(node, "A", "toSI");
double b = getFloat(node, "b");
double E = getFloat(node, "E", "actEnergy") / GasConstant;
double p = getFloat(node, "P", "toSI");
rdata.plogParameters[p] = Arrhenius(A, b, E);
}
} else if (rdata.reactionType == CHEBYSHEV_RXN) {
rdata.rateCoeffType = CHEBYSHEV_REACTION_RATECOEFF_TYPE;
rdata.chebTmin = getFloat(kf, "Tmin", "toSI");
rdata.chebTmax = getFloat(kf, "Tmax", "toSI");
rdata.chebPmin = getFloat(kf, "Pmin", "toSI");
rdata.chebPmax = getFloat(kf, "Pmax", "toSI");
const XML_Node& coeffs = kf.child("floatArray");
rdata.chebDegreeP = atoi(coeffs["degreeP"].c_str());
rdata.chebDegreeT = atoi(coeffs["degreeT"].c_str());
getFloatArray(kf, rdata.chebCoeffs, false);
} else {
throw CanteraError("getRateCoefficient",
"Unknown type: " + type);
}
nodeset_t& kf_children = kf.children();
vector_fp clow(3,0.0), chigh(3,0.0);
for (size_t m = 0; m < kf.nChildren(); m++) {
const node_t& c = *kf_children[m];
string nm = c.name();
int highlow=0;
string type = kf.attrib("type");
if (type == "") {
type = "Arrhenius";
rdata.rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE;
}
if (type == "ExchangeCurrentDensity") {
rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
} else if (type == "Arrhenius") {
if (nm == "Arrhenius") {
vector_fp coeff(3);
if (c["type"] == "stick") {
getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]);
chigh = coeff;
} else {
getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]);
if (highlow == 1 || rdata.reactionType == THREE_BODY_RXN
|| rdata.reactionType == ELEMENTARY_RXN) {
} else {
throw CanteraError("getRateCoefficient", "Unknown type: " + type);
}
vector_fp clow(3,0.0), chigh(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;
if (nm == "Arrhenius") {
vector_fp coeff(3);
if (c["type"] == "stick") {
getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]);
chigh = coeff;
} else {
clow = coeff;
getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]);
if (highlow == 1 || rdata.reactionType == THREE_BODY_RXN
|| rdata.reactionType == ELEMENTARY_RXN) {
chigh = coeff;
} else {
clow = coeff;
}
}
if (rdata.reactionType == SURFACE_RXN) {
getCoverageDependence(c,
kin.thermo(kin.surfacePhaseIndex()), rdata);
}
}
if (rdata.reactionType == SURFACE_RXN) {
getCoverageDependence(c,
kin.thermo(kin.surfacePhaseIndex()), rdata);
}
if (coeff[0] <= 0.0 && negA == 0) {
throw CanteraError("getRateCoefficient",
"negative or zero A coefficient for reaction "+int2str(rdata.number));
if (coeff[0] <= 0.0 && negA == 0) {
throw CanteraError("getRateCoefficient",
"negative or zero A coefficient for reaction "+int2str(rdata.number));
}
} else if (nm == "Arrhenius_ExchangeCurrentDensity") {
vector_fp coeff(3);
getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]);
chigh = coeff;
rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
} else if (nm == "falloff") {
getFalloff(c, rdata);
} else if (nm == "efficiencies") {
getEfficiencies(c, kin, rdata);
} else if (nm == "electrochem") {
rdata.beta = fpValue(c["beta"]);
}
} else if (nm == "Arrhenius_ExchangeCurrentDensity") {
vector_fp coeff(3);
getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]);
chigh = coeff;
rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
} else if (nm == "falloff") {
getFalloff(c, rdata);
} else if (nm == "efficiencies") {
getEfficiencies(c, kin, rdata);
} else if (nm == "electrochem") {
rdata.beta = fpValue(c["beta"]);
}
/*
* Store the coefficients in the ReactionData object for return
* from this function.
*/
if (rdata.reactionType == CHEMACT_RXN) {
rdata.rateCoeffParameters = clow;
} else {
rdata.rateCoeffParameters = chigh;
}
if (rdata.reactionType == FALLOFF_RXN) {
rdata.auxRateCoeffParameters = clow;
} else if (rdata.reactionType == CHEMACT_RXN) {
rdata.auxRateCoeffParameters = chigh;
}
}
/*
* Store the coefficients in the ReactionData object for return
* from this function.
*/
if (rdata.reactionType == CHEMACT_RXN) {
rdata.rateCoeffParameters = clow;
} else {
rdata.rateCoeffParameters = chigh;
}
if (rdata.reactionType == FALLOFF_RXN) {
rdata.auxRateCoeffParameters = clow;
} else if (rdata.reactionType == CHEMACT_RXN) {
rdata.auxRateCoeffParameters = chigh;
}
}
/*
* This function returns true if two reactions are duplicates of
* one another, and false otherwise. The input arguments are two
@ -842,6 +846,10 @@ bool rxninfo::installReaction(int i, const XML_Node& r, Kinetics* k,
rdata.falloffType = SIMPLE_FALLOFF;
} else if (typ == "threeBody") {
rdata.reactionType = THREE_BODY_RXN;
} else if (typ == "plog") {
rdata.reactionType = PLOG_RXN;
} else if (typ == "chebyshev") {
rdata.reactionType = CHEBYSHEV_RXN;
} else if (typ == "surface") {
rdata.reactionType = SURFACE_RXN;
} else if (typ == "edge") {
@ -1208,7 +1216,4 @@ bool buildSolutionFromXML(XML_Node& root, std::string id, std::string nm,
}
}
}