diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index b484ca8ad..11806f8b2 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -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( @@ -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(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(2); + const auto& P_range = node.at("pressure-range").asVector(2); + auto& vcoeffs = node.at("data").asVector(); + 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 newReaction(const AnyMap& node, const Kinetics& kin, unique_ptr R(new PlogReaction()); setupPlogReaction(*R, node, kin, units); return unique_ptr(move(R)); + } else if (type == "Chebyshev") { + unique_ptr R(new ChebyshevReaction()); + setupChebyshevReaction(*R, node, kin, units); + return unique_ptr(move(R)); } else { throw CanteraError("newReaction", "Unknown reaction type '{}'", type); } diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 7c69d0a9e..0fc265660 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -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(*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); +}