Cantera  2.4.0
RxnRates.h
Go to the documentation of this file.
1 /**
2  * @file RxnRates.h
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at http://www.cantera.org/license.txt for license and copyright information.
7 
8 #ifndef CT_RXNRATES_H
9 #define CT_RXNRATES_H
10 
14 
15 #include <iostream>
16 
17 namespace Cantera
18 {
19 
20 class Array2D;
21 
22 //! Arrhenius reaction rate type depends only on temperature
23 /**
24  * A reaction rate coefficient of the following form.
25  *
26  * \f[
27  * k_f = A T^b \exp (-E/RT)
28  * \f]
29  */
30 class Arrhenius
31 {
32 public:
33  //! return the rate coefficient type.
34  static int type() {
35  return ARRHENIUS_REACTION_RATECOEFF_TYPE;
36  }
37 
38  //! Default constructor.
39  Arrhenius();
40 
41  /// Constructor.
42  /// @param A pre-exponential. The unit system is
43  /// (kmol, m, s). The actual units depend on the reaction
44  /// order and the dimensionality (surface or bulk).
45  /// @param b Temperature exponent. Non-dimensional.
46  /// @param E Activation energy in temperature units. Kelvin.
47  Arrhenius(doublereal A, doublereal b, doublereal E);
48 
49  //! Update concentration-dependent parts of the rate coefficient.
50  /*!
51  * For this class, there are no concentration-dependent parts, so this
52  * method does nothing.
53  */
54  void update_C(const doublereal* c) {
55  }
56 
57  /**
58  * Update the value of the natural logarithm of the rate constant.
59  */
60  doublereal updateLog(doublereal logT, doublereal recipT) const {
61  return m_logA + m_b*logT - m_E*recipT;
62  }
63 
64  /**
65  * Update the value the rate constant.
66  *
67  * This function returns the actual value of the rate constant. It can be
68  * safely called for negative values of the pre-exponential factor.
69  */
70  doublereal updateRC(doublereal logT, doublereal recipT) const {
71  return m_A * std::exp(m_b*logT - m_E*recipT);
72  }
73 
74  //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending
75  //! on the reaction order)
76  double preExponentialFactor() const {
77  return m_A;
78  }
79 
80  //! Return the temperature exponent *b*
81  double temperatureExponent() const {
82  return m_b;
83  }
84 
85  //! Return the activation energy divided by the gas constant (i.e. the
86  //! activation temperature) [K]
87  doublereal activationEnergy_R() const {
88  return m_E;
89  }
90 
91 protected:
92  doublereal m_logA, m_b, m_E, m_A;
93 };
94 
95 
96 /**
97  * An Arrhenius rate with coverage-dependent terms.
98  *
99  * The rate expression is given by [Kee, R. J., Coltrin, M. E., & Glarborg, P.
100  * (2005). Chemically reacting flow: theory and practice. John Wiley & Sons.
101  * Eq 11.113]:
102  * \f[
103  * k_f = A T^b \exp \left(
104  * \ln 10 \sum a_k \theta_k
105  * - \frac{1}{RT} \left( E_a + \sum E_k\theta_k \right)
106  * + \sum m_k \ln \theta_k
107  * \right)
108  * \f]
109  * or, equivalently, and as implemented in Cantera,
110  * \f[
111  * k_f = A T^b \exp \left( - \frac{E_a}{RT} \right)
112  * \prod_k 10^{a_k \theta_k} \theta_k^{m_k}
113  * \exp \left( \frac{- E_k \theta_k}{RT} \right)
114  * \f]
115  * where the parameters \f$ (a_k, E_k, m_k) \f$ describe the dependency on the
116  * surface coverage of species \f$k, \theta_k \f$.
117  */
119 {
120 
121 public:
122  static int type() {
123  return SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE;
124  }
125 
127  explicit SurfaceArrhenius(double A, double b, double Ta);
128 
129  //! Add a coverage dependency for species *k*, with pre-exponential
130  //! dependence *a*, rate constant exponential dependency *m*, and activation
131  //! energy dependence *e*, where *e* is in Kelvin, i.e. energy divided by
132  //! the molar gas constant.
133  void addCoverageDependence(size_t k, doublereal a,
134  doublereal m, doublereal e);
135 
136  void update_C(const doublereal* theta) {
137  m_acov = 0.0;
138  m_ecov = 0.0;
139  m_mcov = 0.0;
140  size_t k;
141  doublereal th;
142  for (size_t n = 0; n < m_ac.size(); n++) {
143  k = m_sp[n];
144  m_acov += m_ac[n] * theta[k];
145  m_ecov += m_ec[n] * theta[k];
146  }
147  for (size_t n = 0; n < m_mc.size(); n++) {
148  k = m_msp[n];
149  th = std::max(theta[k], Tiny);
150  m_mcov += m_mc[n]*std::log(th);
151  }
152  }
153 
154  /**
155  * Update the value the rate constant.
156  *
157  * This function returns the actual value of the rate constant. It can be
158  * safely called for negative values of the pre-exponential factor.
159  */
160  doublereal updateRC(doublereal logT, doublereal recipT) const {
161  return m_A * std::exp(std::log(10.0)*m_acov + m_b*logT -
162  (m_E + m_ecov)*recipT + m_mcov);
163  }
164 
165  //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending
166  //! on the reaction order) accounting coverage dependence.
167  /*!
168  * Returns reaction prexponent accounting for both *a* and *m*.
169  */
170  doublereal preExponentialFactor() const {
171  return m_A * std::exp(std::log(10.0)*m_acov + m_mcov);
172  }
173 
174  //! Return effective temperature exponent
175  doublereal temperatureExponent() const {
176  return m_b;
177  }
178 
179  //! Return the activation energy divided by the gas constant (i.e. the
180  //! activation temperature) [K], accounting coverage dependence.
181  doublereal activationEnergy_R() const {
182  return m_E + m_ecov;
183  }
184 
185 protected:
186  doublereal m_b, m_E, m_A;
187  doublereal m_acov, m_ecov, m_mcov;
188  std::vector<size_t> m_sp, m_msp;
189  vector_fp m_ac, m_ec, m_mc;
190 };
191 
192 
193 //! Pressure-dependent reaction rate expressed by logarithmically interpolating
194 //! between Arrhenius rate expressions at various pressures.
195 /*!
196  * Given two rate expressions at two specific pressures:
197  *
198  * * \f$ P_1: k_1(T) = A_1 T^{b_1} e^{E_1 / RT} \f$
199  * * \f$ P_2: k_2(T) = A_2 T^{b_2} e^{E_2 / RT} \f$
200  *
201  * The rate at an intermediate pressure \f$ P_1 < P < P_2 \f$ is computed as
202  * \f[
203  * \log k(T,P) = \log k_1(T) + \bigl(\log k_2(T) - \log k_1(T)\bigr)
204  * \frac{\log P - \log P_1}{\log P_2 - \log P_1}
205  * \f]
206  * Multiple rate expressions may be given at the same pressure, in which case
207  * the rate used in the interpolation formula is the sum of all the rates given
208  * at that pressure. For pressures outside the given range, the rate expression
209  * at the nearest pressure is used.
210  */
211 class Plog
212 {
213 public:
214  //! return the rate coefficient type.
215  static int type() {
216  return PLOG_REACTION_RATECOEFF_TYPE;
217  }
218 
219  //! Default constructor.
220  Plog() {}
221 
222  //! Constructor from Arrhenius rate expressions at a set of pressures
223  explicit Plog(const std::multimap<double, Arrhenius>& rates);
224 
225  //! Update concentration-dependent parts of the rate coefficient.
226  //! @param c natural log of the pressure in Pa
227  void update_C(const doublereal* c) {
228  logP_ = c[0];
229  if (logP_ > logP1_ && logP_ < logP2_) {
230  return;
231  }
232 
233  auto iter = pressures_.upper_bound(c[0]);
234  AssertThrowMsg(iter != pressures_.end(), "Plog::update_C",
235  "Pressure out of range: {}", logP_);
236  AssertThrowMsg(iter != pressures_.begin(), "Plog::update_C",
237  "Pressure out of range: {}", logP_);
238 
239  // upper interpolation pressure
240  logP2_ = iter->first;
241  ihigh1_ = iter->second.first;
242  ihigh2_ = iter->second.second;
243 
244  // lower interpolation pressure
245  logP1_ = (--iter)->first;
246  ilow1_ = iter->second.first;
247  ilow2_ = iter->second.second;
248 
249  rDeltaP_ = 1.0 / (logP2_ - logP1_);
250  }
251 
252  /**
253  * Update the value the rate constant.
254  *
255  * This function returns the actual value of the rate constant.
256  */
257  doublereal updateRC(doublereal logT, doublereal recipT) const {
258  double log_k1, log_k2;
259  if (ilow1_ == ilow2_) {
260  log_k1 = rates_[ilow1_].updateLog(logT, recipT);
261  } else {
262  double k = 1e-300; // non-zero to make log(k) finite
263  for (size_t i = ilow1_; i < ilow2_; i++) {
264  k += rates_[i].updateRC(logT, recipT);
265  }
266  log_k1 = std::log(k);
267  }
268 
269  if (ihigh1_ == ihigh2_) {
270  log_k2 = rates_[ihigh1_].updateLog(logT, recipT);
271  } else {
272  double k = 1e-300; // non-zero to make log(k) finite
273  for (size_t i = ihigh1_; i < ihigh2_; i++) {
274  k += rates_[i].updateRC(logT, recipT);
275  }
276  log_k2 = std::log(k);
277  }
278 
279  return std::exp(log_k1 + (log_k2-log_k1) * (logP_-logP1_) * rDeltaP_);
280  }
281 
282  //! Check to make sure that the rate expression is finite over a range of
283  //! temperatures at each interpolation pressure. This is potentially an
284  //! issue when one of the Arrhenius expressions at a particular pressure
285  //! has a negative pre-exponential factor.
286  void validate(const std::string& equation);
287 
288  //! Return the pressures and Arrhenius expressions which comprise this
289  //! reaction.
290  std::vector<std::pair<double, Arrhenius> > rates() const;
291 
292 protected:
293  //! log(p) to (index range) in the rates_ vector
294  std::map<double, std::pair<size_t, size_t> > pressures_;
295 
296  // Rate expressions which are referenced by the indices stored in pressures_
297  std::vector<Arrhenius> rates_;
298 
299  double logP_; //!< log(p) at the current state
300  double logP1_, logP2_; //!< log(p) at the lower / upper pressure reference
301 
302  //! Indices to the ranges within rates_ for the lower / upper pressure, such
303  //! that rates_[ilow1_] through rates_[ilow2_] (inclusive) are the rates
304  //! expressions which are combined to form the rate at the lower reference
305  //! pressure.
306  size_t ilow1_, ilow2_, ihigh1_, ihigh2_;
307 
308  double rDeltaP_; //!< reciprocal of (logP2 - logP1)
309 };
310 
311 //! Pressure-dependent rate expression where the rate coefficient is expressed
312 //! as a bivariate Chebyshev polynomial in temperature and pressure.
313 /*!
314  * The rate constant can be written as:
315  * \f[
316  * \log k(T,P) = \sum_{t=1}^{N_T} \sum_{p=1}^{N_P} \alpha_{tp}
317  * \phi_t(\tilde{T}) \phi_p(\tilde{P})
318  * \f]
319  * where \f$\alpha_{tp}\f$ are the constants defining the rate, \f$\phi_n(x)\f$
320  * is the Chebyshev polynomial of the first kind of degree *n* evaluated at
321  * *x*, and
322  * \f[
323  * \tilde{T} \equiv \frac{2T^{-1} - T_\mathrm{min}^{-1} - T_\mathrm{max}^{-1}}
324  * {T_\mathrm{max}^{-1} - T_\mathrm{min}^{-1}}
325  * \f]
326  * \f[
327  * \tilde{P} \equiv \frac{2 \log P - \log P_\mathrm{min} - \log P_\mathrm{max}}
328  * {\log P_\mathrm{max} - \log P_\mathrm{min}}
329  * \f]
330  * are reduced temperature and reduced pressures which map the ranges
331  * \f$ (T_\mathrm{min}, T_\mathrm{max}) \f$ and
332  * \f$ (P_\mathrm{min}, P_\mathrm{max}) \f$ to (-1, 1).
333  *
334  * A Chebyshev rate expression is specified in terms of the coefficient matrix
335  * \f$ \alpha \f$ and the temperature and pressure ranges. Note that the
336  * Chebyshev polynomials are not defined outside the interval (-1,1), and
337  * therefore extrapolation of rates outside the range of temperatures and
338  * pressures for which they are defined is strongly discouraged.
339  */
341 {
342 public:
343  //! return the rate coefficient type.
344  static int type() {
345  return CHEBYSHEV_REACTION_RATECOEFF_TYPE;
346  }
347 
348  //! Default constructor.
350 
351  //! Constructor directly from coefficient array
352  /*
353  * @param Tmin Minimum temperature [K]
354  * @param Tmax Maximum temperature [K]
355  * @param Pmin Minimum pressure [Pa]
356  * @param Pmax Maximum pressure [Pa]
357  * @param coeffs Coefficient array dimensioned `nT` by `nP` where `nT` and
358  * `nP` are the number of temperatures and pressures used in the fit,
359  * respectively.
360  */
361  ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax,
362  const Array2D& coeffs);
363 
364  //! Update concentration-dependent parts of the rate coefficient.
365  //! @param c base-10 logarithm of the pressure in Pa
366  void update_C(const doublereal* c) {
367  double Pr = (2 * c[0] + PrNum_) * PrDen_;
368  double Cnm1 = 1;
369  double Cn = Pr;
370  double Cnp1;
371  for (size_t j = 0; j < nT_; j++) {
372  dotProd_[j] = chebCoeffs_[nP_*j] + Pr * chebCoeffs_[nP_*j+1];
373  }
374  for (size_t i = 2; i < nP_; i++) {
375  Cnp1 = 2 * Pr * Cn - Cnm1;
376  for (size_t j = 0; j < nT_; j++) {
377  dotProd_[j] += Cnp1 * chebCoeffs_[nP_*j + i];
378  }
379  Cnm1 = Cn;
380  Cn = Cnp1;
381  }
382  }
383 
384  /**
385  * Update the value the rate constant.
386  *
387  * This function returns the actual value of the rate constant.
388  */
389  doublereal updateRC(doublereal logT, doublereal recipT) const {
390  double Tr = (2 * recipT + TrNum_) * TrDen_;
391  double Cnm1 = 1;
392  double Cn = Tr;
393  double Cnp1;
394  double logk = dotProd_[0] + Tr * dotProd_[1];
395  for (size_t i = 2; i < nT_; i++) {
396  Cnp1 = 2 * Tr * Cn - Cnm1;
397  logk += Cnp1 * dotProd_[i];
398  Cnm1 = Cn;
399  Cn = Cnp1;
400  }
401  return std::pow(10, logk);
402  }
403 
404  //! Minimum valid temperature [K]
405  double Tmin() const {
406  return Tmin_;
407  }
408 
409  //! Maximum valid temperature [K]
410  double Tmax() const {
411  return Tmax_;
412  }
413 
414  //! Minimum valid pressure [Pa]
415  double Pmin() const {
416  return Pmin_;
417  }
418 
419  //! Maximum valid pressure [Pa]
420  double Pmax() const {
421  return Pmax_;
422  }
423 
424  //! Number of points in the pressure direction
425  size_t nPressure() const {
426  return nP_;
427  }
428 
429  //! Number of points in the temperature direction
430  size_t nTemperature() const {
431  return nT_;
432  }
433 
434  //! Access the Chebyshev coefficients.
435  /*!
436  * \f$ \alpha_{t,p} = \mathrm{coeffs}[N_P*t + p] \f$ where
437  * \f$ 0 <= t < N_T \f$ and \f$ 0 <= p < N_P \f$.
438  */
439  const vector_fp& coeffs() const {
440  return chebCoeffs_;
441  }
442 
443 protected:
444  double Tmin_, Tmax_; //!< valid temperature range
445  double Pmin_, Pmax_; //!< valid pressure range
446  double TrNum_, TrDen_; //!< terms appearing in the reduced temperature
447  double PrNum_, PrDen_; //!< terms appearing in the reduced pressure
448 
449  size_t nP_; //!< number of points in the pressure direction
450  size_t nT_; //!< number of points in the temperature direction
451  vector_fp chebCoeffs_; //!< Chebyshev coefficients, length nP * nT
452  vector_fp dotProd_; //!< dot product of chebCoeffs with the reduced pressure polynomial
453 };
454 
455 }
456 
457 #endif
size_t ilow1_
Indices to the ranges within rates_ for the lower / upper pressure, such that rates_[ilow1_] through ...
Definition: RxnRates.h:306
void update_C(const doublereal *c)
Update concentration-dependent parts of the rate coefficient.
Definition: RxnRates.h:54
doublereal updateRC(doublereal logT, doublereal recipT) const
Update the value the rate constant.
Definition: RxnRates.h:257
Plog()
Default constructor.
Definition: RxnRates.h:220
doublereal preExponentialFactor() const
Return the pre-exponential factor A (in m, kmol, s to powers depending on the reaction order) account...
Definition: RxnRates.h:170
void addCoverageDependence(size_t k, doublereal a, doublereal m, doublereal e)
Add a coverage dependency for species k, with pre-exponential dependence a, rate constant exponential...
Definition: RxnRates.cpp:51
size_t nP_
number of points in the pressure direction
Definition: RxnRates.h:449
double Tmin() const
Minimum valid temperature [K].
Definition: RxnRates.h:405
static int type()
return the rate coefficient type.
Definition: RxnRates.h:215
Pressure-dependent rate expression where the rate coefficient is expressed as a bivariate Chebyshev p...
Definition: RxnRates.h:340
size_t nPressure() const
Number of points in the pressure direction.
Definition: RxnRates.h:425
double Pmin() const
Minimum valid pressure [Pa].
Definition: RxnRates.h:415
double logP2_
log(p) at the lower / upper pressure reference
Definition: RxnRates.h:300
std::vector< std::pair< double, Arrhenius > > rates() const
Return the pressures and Arrhenius expressions which comprise this reaction.
Definition: RxnRates.cpp:110
void update_C(const doublereal *c)
Update concentration-dependent parts of the rate coefficient.
Definition: RxnRates.h:366
Arrhenius()
Default constructor.
Definition: RxnRates.cpp:11
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:31
doublereal temperatureExponent() const
Return effective temperature exponent.
Definition: RxnRates.h:175
size_t nTemperature() const
Number of points in the temperature direction.
Definition: RxnRates.h:430
This file defines some constants used to specify reaction types.
double preExponentialFactor() const
Return the pre-exponential factor A (in m, kmol, s to powers depending on the reaction order) ...
Definition: RxnRates.h:76
doublereal updateRC(doublereal logT, doublereal recipT) const
Update the value the rate constant.
Definition: RxnRates.h:70
vector_fp dotProd_
dot product of chebCoeffs with the reduced pressure polynomial
Definition: RxnRates.h:452
double TrDen_
terms appearing in the reduced temperature
Definition: RxnRates.h:446
doublereal activationEnergy_R() const
Return the activation energy divided by the gas constant (i.e.
Definition: RxnRates.h:87
double temperatureExponent() const
Return the temperature exponent b
Definition: RxnRates.h:81
void validate(const std::string &equation)
Check to make sure that the rate expression is finite over a range of temperatures at each interpolat...
Definition: RxnRates.cpp:91
const vector_fp & coeffs() const
Access the Chebyshev coefficients.
Definition: RxnRates.h:439
Pressure-dependent reaction rate expressed by logarithmically interpolating between Arrhenius rate ex...
Definition: RxnRates.h:211
doublereal activationEnergy_R() const
Return the activation energy divided by the gas constant (i.e.
Definition: RxnRates.h:181
void update_C(const doublereal *c)
Update concentration-dependent parts of the rate coefficient.
Definition: RxnRates.h:227
doublereal updateLog(doublereal logT, doublereal recipT) const
Update the value of the natural logarithm of the rate constant.
Definition: RxnRates.h:60
Arrhenius reaction rate type depends only on temperature.
Definition: RxnRates.h:30
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:264
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
double Pmax_
valid pressure range
Definition: RxnRates.h:445
doublereal updateRC(doublereal logT, doublereal recipT) const
Update the value the rate constant.
Definition: RxnRates.h:389
size_t nT_
number of points in the temperature direction
Definition: RxnRates.h:450
std::map< double, std::pair< size_t, size_t > > pressures_
log(p) to (index range) in the rates_ vector
Definition: RxnRates.h:294
double Tmax() const
Maximum valid temperature [K].
Definition: RxnRates.h:410
const doublereal Tiny
Small number to compare differences of mole fractions against.
Definition: ct_defs.h:143
double logP_
log(p) at the current state
Definition: RxnRates.h:299
ChebyshevRate()
Default constructor.
Definition: RxnRates.h:349
Contains declarations for string manipulation functions within Cantera.
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
doublereal updateRC(doublereal logT, doublereal recipT) const
Update the value the rate constant.
Definition: RxnRates.h:160
static int type()
return the rate coefficient type.
Definition: RxnRates.h:344
double PrDen_
terms appearing in the reduced pressure
Definition: RxnRates.h:447
double rDeltaP_
reciprocal of (logP2 - logP1)
Definition: RxnRates.h:308
An Arrhenius rate with coverage-dependent terms.
Definition: RxnRates.h:118
double Tmax_
valid temperature range
Definition: RxnRates.h:444
static int type()
return the rate coefficient type.
Definition: RxnRates.h:34
vector_fp chebCoeffs_
Chebyshev coefficients, length nP * nT.
Definition: RxnRates.h:451
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
double Pmax() const
Maximum valid pressure [Pa].
Definition: RxnRates.h:420