diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 089b60754..6c4f7ef40 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -17,6 +17,8 @@ namespace Cantera { +class Plog; + //! Arrhenius reaction rate type depends only on temperature /** * A reaction rate coefficient of the following form. @@ -104,6 +106,8 @@ public: protected: doublereal m_logA, m_b, m_E, m_A; + + friend class Plog; }; @@ -303,6 +307,9 @@ public: //! Constructor from ReactionData. explicit Plog(const ReactionData& rdata); + //! Constructor from Arrhenius rate expressions at a set of pressures + explicit Plog(const std::multimap& rates); + //! Update concentration-dependent parts of the rate coefficient. //! @param c natural log of the pressure in Pa void update_C(const doublereal* c) { diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index a47672958..89fd2e036 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -181,6 +181,62 @@ Plog::Plog(const ReactionData& rdata) } } +Plog::Plog(const std::multimap& rates) + : logP_(-1000) + , logP1_(1000) + , logP2_(-1000) + , m1_(npos) + , m2_(npos) + , rDeltaP_(-1.0) + , maxRates_(1) +{ + size_t j = 0; + size_t rateCount = 0; + // Insert intermediate pressures + for (std::multimap::const_iterator iter = rates.begin(); + iter != rates.end(); + iter++) { + double logp = std::log(iter->first); + if (pressures_.empty() || pressures_.rbegin()->first != logp) { + // starting a new group + pressures_[logp] = std::make_pair(j, j+1); + rateCount = 1; + } else { + // another rate expression at the same pressure + pressures_[logp].second = j+1; + rateCount++; + } + maxRates_ = std::max(rateCount, maxRates_); + + j++; + A_.push_back(iter->second.m_A); + n_.push_back(iter->second.m_b); + Ea_.push_back(iter->second.m_E); + } + + // For pressures with only one Arrhenius expression, it is more + // efficient to work with log(A) + for (pressureIter iter = pressures_.begin(); + iter != pressures_.end(); + iter++) { + if (iter->second.first == iter->second.second - 1) { + A_[iter->second.first] = std::log(A_[iter->second.first]); + } + } + + // Duplicate the first and last groups to handle P < P_0 and P > P_N + pressures_.insert(std::make_pair(-1000.0, pressures_.begin()->second)); + pressures_.insert(std::make_pair(1000.0, pressures_.rbegin()->second)); + + // Resize work arrays + A1_.resize(maxRates_); + A2_.resize(maxRates_); + n1_.resize(maxRates_); + n2_.resize(maxRates_); + Ea1_.resize(maxRates_); + Ea2_.resize(maxRates_); +} + void Plog::validate(const ReactionData& rdata) { double T[] = {200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0}; diff --git a/test/data/kineticsfromscratch.cti b/test/data/kineticsfromscratch.cti index f2e6b5511..38add76cc 100644 --- a/test/data/kineticsfromscratch.cti +++ b/test/data/kineticsfromscratch.cti @@ -18,3 +18,9 @@ falloff_reaction('2 OH (+ M) <=> H2O2 (+ M)', kf0=[2.300000e+12, -0.9, -1700.0], efficiencies='AR:0.7 H2:2.0 H2O:6.0', falloff=Troe(A=0.7346, T3=94.0, T1=1756.0, T2=5182.0)) + +pdep_arrhenius('H2 + O2 <=> 2 OH', + [(0.01, 'atm'), 1.212400e+16, -0.5779, 10872.7], + [(1.0, 'atm'), 4.910800e+31, -4.8507, 24772.8], + [(10.0, 'atm'), 1.286600e+47, -9.0246, 39796.5], + [(100.0, 'atm'), 5.963200e+56, -11.529, 52599.6]) diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index 08a46b939..7ba8c2335 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -101,3 +101,26 @@ TEST_F(KineticsFromScratch, add_falloff_reaction) kin.finalize(); check_rates(2); } + +TEST_F(KineticsFromScratch, add_plog_reaction) +{ + // reaction 3: + // pdep_arrhenius('H2 + O2 <=> 2 OH', + // [(0.01, 'atm'), 1.212400e+16, -0.5779, 10872.7], + // [(1.0, 'atm'), 4.910800e+31, -4.8507, 24772.8], + // [(10.0, 'atm'), 1.286600e+47, -9.0246, 39796.5], + // [(100.0, 'atm'), 5.963200e+56, -11.529, 52599.6]) + Composition reac = parseCompString("H2:1, O2:1"); + Composition prod = parseCompString("OH:2"); + std::multimap rates; + typedef std::multimap::value_type item; + rates.insert(item(0.01*101325, Arrhenius(1.212400e+16, -0.5779, 10872.7 / GasConst_cal_mol_K))); + rates.insert(item(1.0*101325, Arrhenius(4.910800e+31, -4.8507, 24772.8 / GasConst_cal_mol_K))); + rates.insert(item(10.0*101325, Arrhenius(1.286600e+47, -9.0246, 39796.5 / GasConst_cal_mol_K))); + rates.insert(item(100.0*101325, Arrhenius(5.963200e+56, -11.529, 52599.6 / GasConst_cal_mol_K))); + + shared_ptr R(new PlogReaction(reac, prod, Plog(rates))); + kin.addReaction(R); + kin.finalize(); + check_rates(3); +}