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

[Input] Parse YAML entries for Chebyshev reactions

Browse source
This commit is contained in:
Ray Speth 2018-12-13 17:21:26 -05:00
parent 7a97ff9917
commit aef101fee7
2 changed files with 68 additions and 4 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

43
src/kinetics/Reaction.cpp
View file

@ -284,9 +284,8 @@ Arrhenius readArrhenius(const XML_Node& arrhenius_node)
getFloat(arrhenius_node, "E", "actEnergy") / GasConstant);
}
Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate,
const Kinetics& kin, const UnitSystem& units,
int pressure_dependence=0)
Units rateCoeffUnits(const Reaction& R, const Kinetics& kin,
int pressure_dependence=0)
{
// Determine the units of the rate coefficient
double reaction_phase_ndim = static_cast<double>(
@ -313,9 +312,17 @@ Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate,
// pressure chemically-activated reaction limits
len_dim += pressure_dependence * reaction_phase_ndim;
quantity_dim -= pressure_dependence;
return Units(1.0, 0, len_dim, -1, 0, 0, quantity_dim);
}
Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate,
const Kinetics& kin, const UnitSystem& units,
int pressure_dependence=0)
{
auto& rate_vec = rate.asVector<AnyValue>(3);
double A = units.convert(rate_vec[0], Units(1.0, 0, len_dim, -1, 0, 0, quantity_dim));
Units rc_units = rateCoeffUnits(R, kin, pressure_dependence);
double A = units.convert(rate_vec[0], rc_units);
double b = rate_vec[1].asDouble();
double Ta = units.convertMolarEnergy(rate_vec[2], "K");
return Arrhenius(A, b, Ta);
@ -707,6 +714,30 @@ void setupChebyshevReaction(ChebyshevReaction& R, const XML_Node& rxn_node)
setupReaction(R, rxn_node);
}
void setupChebyshevReaction(ChebyshevReaction&R, const AnyMap& node,
const Kinetics& kin, const UnitSystem& units)
{
setupReaction(R, node);
R.reactants.erase("(+M)"); // remove optional third body notation
R.products.erase("(+M)");
const auto& T_range = node.at("temperature-range").asVector<AnyValue>(2);
const auto& P_range = node.at("pressure-range").asVector<AnyValue>(2);
auto& vcoeffs = node.at("data").asVector<vector_fp>();
Array2D coeffs(vcoeffs.size(), vcoeffs[0].size());
for (size_t i = 0; i < coeffs.nRows(); i++) {
for (size_t j = 0; j < coeffs.nColumns(); j++) {
coeffs(i, j) = vcoeffs[i][j];
}
}
Units rcUnits = rateCoeffUnits(R, kin);
coeffs(0, 0) += std::log10(units.convert(1.0, rcUnits));
R.rate = ChebyshevRate(units.convert(T_range[0], "K"),
units.convert(T_range[1], "K"),
units.convert(P_range[0], "Pa"),
units.convert(P_range[1], "Pa"),
coeffs);
}
void setupInterfaceReaction(InterfaceReaction& R, const XML_Node& rxn_node)
{
if (caseInsensitiveEquals(rxn_node["type"], "global")) {
@ -913,6 +944,10 @@ unique_ptr<Reaction> newReaction(const AnyMap& node, const Kinetics& kin,
unique_ptr<PlogReaction> R(new PlogReaction());
setupPlogReaction(*R, node, kin, units);
return unique_ptr<Reaction>(move(R));
} else if (type == "Chebyshev") {
unique_ptr<ChebyshevReaction> R(new ChebyshevReaction());
setupChebyshevReaction(*R, node, kin, units);
return unique_ptr<Reaction>(move(R));
} else {
throw CanteraError("newReaction", "Unknown reaction type '{}'", type);
}

29
test/kinetics/kineticsFromYaml.cpp
View file

@ -141,3 +141,32 @@ TEST(Reaction, PlogFromYaml)
EXPECT_DOUBLE_EQ(rates[0].second.preExponentialFactor(), 2.72e6);
EXPECT_DOUBLE_EQ(rates[3].second.preExponentialFactor(), 1.68e16);
}
TEST(Reaction, ChebyshevFromYaml)
{
IdealGasMix gas("gri30.xml");
AnyMap rxn = AnyMap::fromYamlString(
"equation: 'CH4 <=> CH3 + H'\n"
"type: Chebyshev\n"
"temperature-range: [290, 3000]\n"
"pressure-range: [0.0098692326671601278 atm, 98.692326671601279 atm]\n"
"data: [[-1.44280e+01, 2.59970e-01, -2.24320e-02, -2.78700e-03],\n"
" [ 2.20630e+01, 4.88090e-01, -3.96430e-02, -5.48110e-03],\n"
" [-2.32940e-01, 4.01900e-01, -2.60730e-02, -5.04860e-03],\n"
" [-2.93660e-01, 2.85680e-01, -9.33730e-03, -4.01020e-03],\n"
" [-2.26210e-01, 1.69190e-01, 4.85810e-03, -2.38030e-03],\n"
" [-1.43220e-01, 7.71110e-02, 1.27080e-02, -6.41540e-04]]\n");
UnitSystem U;
auto R = newReaction(rxn, gas, U);
EXPECT_EQ(R->reactants.size(), (size_t) 1);
auto CR = dynamic_cast<ChebyshevReaction&>(*R);
double logP = std::log10(2e6);
double T = 1800;
CR.rate.update_C(&logP);
EXPECT_EQ(CR.rate.nTemperature(), (size_t) 6);
EXPECT_EQ(CR.rate.nPressure(), (size_t) 4);
EXPECT_DOUBLE_EQ(CR.rate.Tmax(), 3000);
EXPECT_DOUBLE_EQ(CR.rate.Pmin(), 1000);
EXPECT_NEAR(CR.rate.updateRC(std::log(T), 1.0/T), 130512.2773948636, 1e-9);
}