Cantera  2.0
FalloffFactory.cpp
Go to the documentation of this file.
1 /**
2  * @file FalloffFactory.cpp
3  */
4 // Copyright 2001 California Institute of Technology
5 
6 #include "cantera/kinetics/FalloffFactory.h"
7 #include "cantera/base/ctexceptions.h"
8 #include "cantera/base/stringUtils.h"
9 
10 #include <cmath>
11 
12 namespace Cantera
13 {
14 
15 FalloffFactory* FalloffFactory::s_factory = 0;
16 mutex_t FalloffFactory::falloff_mutex;
17 
18 //! The 3-parameter Troe falloff parameterization.
19 /*!
20  * The falloff function defines the value of \f$ F \f$ in the following
21  * rate expression
22  *
23  * \f[
24  * k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
25  * \f]
26  * where
27  * \f[
28  * P_r = \frac{k_0 [M]}{k_{\infty}}
29  * \f]
30  *
31  * This parameterization is defined by
32  * \f[
33  * F = F_{cent}^{1/(1 + f_1^2)}
34  * \f]
35  * where
36  * \f[
37  * F_{cent} = (1 - A)\exp(-T/T_3) + A \exp(-T/T_1)
38  * \f]
39  *
40  * \f[
41  * f_1 = (\log_{10} P_r + C) / \left(N - 0.14
42  * (\log_{10} P_r + C)\right)
43  * \f]
44  *
45  * \f[
46  * C = -0.4 - 0.67 \log_{10} F_{cent}
47  * \f]
48  *
49  * \f[
50  * N = 0.75 - 1.27 \log_{10} F_{cent}
51  * \f]
52  *
53  * There are a few requirements for the parameters
54  *
55  * T_3 is required to greater than or equal to zero. If it is zero,
56  * then the term is set to zero.
57  *
58  * T_1 is required to greater than or equal to zero. If it is zero,
59  * then the term is set to zero.
60  *
61  * @ingroup falloffGroup
62  */
63 class Troe3 : public Falloff
64 {
65 public:
66 
67  //! Default constructor.
68  Troe3() : m_a(0.0), m_rt3(0.0), m_rt1(0.0) {}
69 
70  //! Destructor. Does nothing.
71  virtual ~Troe3() {}
72 
73  /**
74  * Initialize.
75  * @param c Coefficient vector of length 3,
76  * with entries \f$ (A, T_3, T_1) \f$
77  */
78  virtual void init(const vector_fp& c) {
79  m_a = c[0];
80 
81  if (c[1] <= 0.0) {
82  if (c[1] == 0.0) {
83  m_rt3 = 1000.;
84  } else {
85  throw CanteraError("Troe3::init()", "T3 parameter is less than zero");
86  }
87  } else {
88  m_rt3 = 1.0/c[1];
89  }
90  if (c[2] <= 0.0) {
91  if (c[2] == 0.0) {
92  m_rt1 = 1000.;
93  } else {
94  throw CanteraError("Troe3::init()", "T1 parameter is less than zero");
95  }
96  } else {
97  m_rt1 = 1.0/c[2];
98  }
99  }
100 
101  //! Update the temperature parameters in the representation
102  /*!
103  * The workspace has a length of one
104  *
105  * @param T Temperature (Kelvin)
106  * @param work Vector of working space representing
107  * the temperature dependent part of the
108  * parameterization.
109  */
110  virtual void updateTemp(doublereal T, doublereal* work) const {
111  doublereal Fcent = (1.0 - m_a) * exp(- T * m_rt3)
112  + m_a * exp(- T * m_rt1);
113  *work = log10(std::max(Fcent, SmallNumber));
114  }
115 
116  //! Function that returns <I>F</I>
117  /*!
118  * @param pr Value of the reduced pressure for this reaction
119  * @param work Pointer to the previously saved work space
120  */
121  virtual doublereal F(doublereal pr, const doublereal* work) const {
122  doublereal lpr,f1,lgf, cc, nn;
123  lpr = log10(std::max(pr,SmallNumber));
124  cc = -0.4 - 0.67 * (*work);
125  nn = 0.75 - 1.27 * (*work);
126  f1 = (lpr + cc)/ (nn - 0.14 * (lpr + cc));
127  lgf = (*work) / (1.0 + f1 * f1);
128  return pow(10.0, lgf);
129  }
130 
131  //! Utility function that returns the size of the workspace
132  virtual size_t workSize() {
133  return 1;
134  }
135 
136 protected:
137 
138  //! parameter a in the 4-parameter Troe falloff function
139  /*!
140  * This is unitless
141  */
142  doublereal m_a;
143 
144  //! parameter 1/T_3 in the 4-parameter Troe falloff function
145  /*!
146  * This has units of Kelvin-1
147  */
148  doublereal m_rt3;
149 
150  //! parameter 1/T_1 in the 4-parameter Troe falloff function
151  /*!
152  * This has units of Kelvin-1
153  */
154  doublereal m_rt1;
155 
156 };
157 
158 //! The 4-parameter Troe falloff parameterization.
159 /*!
160  * The falloff function defines the value of \f$ F \f$ in the following
161  * rate expression
162  *
163  * \f[
164  * k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
165  * \f]
166  * where
167  * \f[
168  * P_r = \frac{k_0 [M]}{k_{\infty}}
169  * \f]
170  *
171  * This parameterization is defined by
172  *
173  * \f[
174  * F = F_{cent}^{1/(1 + f_1^2)}
175  * \f]
176  * where
177  * \f[
178  * F_{cent} = (1 - A)\exp(-T/T_3) + A \exp(-T/T_1) + \exp(-T_2/T)
179  * \f]
180  *
181  * \f[
182  * f_1 = (\log_{10} P_r + C) /
183  * \left(N - 0.14 (\log_{10} P_r + C)\right)
184  * \f]
185  *
186  * \f[
187  * C = -0.4 - 0.67 \log_{10} F_{cent}
188  * \f]
189  *
190  * \f[
191  * N = 0.75 - 1.27 \log_{10} F_{cent}
192  * \f]
193  *
194  *
195  * There are a few requirements for the parameters
196  *
197  * T_3 is required to greater than or equal to zero. If it is zero,
198  * then the term is set to zero.
199  *
200  * T_1 is required to greater than or equal to zero. If it is zero,
201  * then the term is set to zero.
202  *
203  * @ingroup falloffGroup
204  */
205 class Troe4 : public Falloff
206 {
207 public:
208  //! Constructor
209  Troe4() : m_a(0.0), m_rt3(0.0), m_rt1(0.0),
210  m_t2(0.0) {}
211 
212  //! Destructor
213  virtual ~Troe4() {}
214 
215 
216  //! Initialization of the object
217  /*!
218  * @param c Vector of four doubles: The doubles are the parameters,
219  * a,, T_3, T_1, and T_2 of the SRI parameterization
220  */
221  virtual void init(const vector_fp& c) {
222  m_a = c[0];
223  if (c[1] <= 0.0) {
224  if (c[1] == 0.0) {
225  m_rt3 = 1000.;
226  } else {
227  throw CanteraError("Troe4::init()", "T3 parameter is less than zero");
228  }
229  } else {
230  m_rt3 = 1.0/c[1];
231  }
232  if (c[2] <= 0.0) {
233  if (c[2] == 0.0) {
234  m_rt1 = 1000.;
235  } else {
236  throw CanteraError("Troe4::init()", "T1 parameter is less than zero");
237  }
238  } else {
239  m_rt1 = 1.0/c[2];
240  }
241  m_t2 = c[3];
242  }
243 
244  //! Update the temperature parameters in the representation
245  /*!
246  * The workspace has a length of one
247  *
248  * @param T Temperature (Kelvin)
249  * @param work Vector of working space representing
250  * the temperature dependent part of the
251  * parameterization.
252  */
253  virtual void updateTemp(doublereal T, doublereal* work) const {
254  doublereal Fcent = (1.0 - m_a) * exp(- T * m_rt3)
255  + m_a * exp(- T * m_rt1)
256  + exp(- m_t2 / T);
257  *work = log10(std::max(Fcent, SmallNumber));
258  }
259 
260  //! Function that returns <I>F</I>
261  /*!
262  * @param pr Value of the reduced pressure for this reaction
263  * @param work Pointer to the previously saved work space
264  */
265  virtual doublereal F(doublereal pr, const doublereal* work) const {
266  doublereal lpr,f1,lgf, cc, nn;
267  lpr = log10(std::max(pr,SmallNumber));
268  cc = -0.4 - 0.67 * (*work);
269  nn = 0.75 - 1.27 * (*work);
270  f1 = (lpr + cc)/ (nn - 0.14 * (lpr + cc));
271  lgf = (*work) / (1.0 + f1 * f1);
272  return pow(10.0, lgf);
273  }
274 
275  //! Utility function that returns the size of the workspace
276  virtual size_t workSize() {
277  return 1;
278  }
279 
280 protected:
281 
282  //! parameter a in the 4-parameter Troe falloff function
283  /*!
284  * This is unitless
285  */
286  doublereal m_a;
287 
288  //! parameter 1/T_3 in the 4-parameter Troe falloff function
289  /*!
290  * This has units of Kelvin-1
291  */
292  doublereal m_rt3;
293 
294  //! parameter 1/T_1 in the 4-parameter Troe falloff function
295  /*!
296  * This has units of Kelvin-1
297  */
298  doublereal m_rt1;
299 
300  //! parameter T_2 in the 4-parameter Troe falloff function
301  /*!
302  * This has units of Kelvin
303  */
304  doublereal m_t2;
305 };
306 
307 
308 //! The 3-parameter SRI falloff function for <I>F</I>
309 /*!
310  * The falloff function defines the value of \f$ F \f$ in the following
311  * rate expression
312  *
313  * \f[
314  * k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
315  * \f]
316  * where
317  * \f[
318  * P_r = \frac{k_0 [M]}{k_{\infty}}
319  * \f]
320  *
321  * \f[
322  * F = {\left( a \; exp(\frac{-b}{T}) + exp(\frac{-T}{c})\right)}^n
323  * \f]
324  * where
325  * \f[
326  * n = \frac{1.0}{1.0 + {\log_{10} P_r}^2}
327  * \f]
328  *
329  * \f$ c \f$ s required to greater than or equal to zero. If it is zero,
330  * then the corresponding term is set to zero.
331  *
332  * @ingroup falloffGroup
333  */
334 class SRI3 : public Falloff
335 {
336 
337 public:
338 
339  //! Constructor
340  SRI3() {}
341 
342  //! Destructor
343  virtual ~SRI3() {}
344 
345  //! Initialization of the object
346  /*!
347  * @param c Vector of three doubles: The doubles are the parameters,
348  * a, b, and c of the SRI parameterization
349  */
350  virtual void init(const vector_fp& c) {
351  if (c[2] < 0.0) {
352  throw CanteraError("SRI3::init()",
353  "m_c parameter is less than zero: " + fp2str(c[2]));
354  }
355  m_a = c[0];
356  m_b = c[1];
357  m_c = c[2];
358  }
359 
360  //! Update the temperature parameters in the representation
361  /*!
362  * The workspace has a length of one
363  *
364  * @param T Temperature (Kelvin)
365  * @param work Vector of working space representing
366  * the temperature dependent part of the
367  * parameterization.
368  */
369  virtual void updateTemp(doublereal T, doublereal* work) const {
370  *work = m_a * exp(- m_b / T);
371  if (m_c != 0.0) {
372  *work += exp(- T/m_c);
373  }
374  }
375 
376  //! Function that returns <I>F</I>
377  /*!
378  * @param pr Value of the reduced pressure for this reaction
379  * @param work Pointer to the previously saved work space
380  */
381  virtual doublereal F(doublereal pr, const doublereal* work) const {
382  doublereal lpr = log10(std::max(pr,SmallNumber));
383  doublereal xx = 1.0/(1.0 + lpr*lpr);
384  doublereal ff = pow(*work , xx);
385  return ff;
386  }
387 
388  //! Utility function that returns the size of the workspace
389  virtual size_t workSize() {
390  return 1;
391  }
392 
393 protected:
394 
395  //! parameter a in the 3-parameter SRI falloff function
396  /*!
397  * This is unitless
398  */
399  doublereal m_a;
400 
401  //! parameter b in the 3-parameter SRI falloff function
402  /*!
403  * This has units of Kelvin
404  */
405  doublereal m_b;
406 
407  //! parameter c in the 3-parameter SRI falloff function
408  /*!
409  * This has units of Kelvin
410  */
411  doublereal m_c;
412 };
413 
414 
415 //! The 5-parameter SRI falloff function.
416 /*!
417  * The falloff function defines the value of \f$ F \f$ in the following
418  * rate expression
419  *
420  * \f[
421  * k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
422  * \f]
423  * where
424  * \f[
425  * P_r = \frac{k_0 [M]}{k_{\infty}}
426  * \f]
427  *
428  * \f[
429  * F = {\left( a \; exp(\frac{-b}{T}) + exp(\frac{-T}{c})\right)}^n
430  * \; d \; exp(\frac{-e}{T})
431  * \f]
432  * where
433  * \f[
434  * n = \frac{1.0}{1.0 + {\log_{10} P_r}^2}
435  * \f]
436  *
437  * \f$ c \f$ s required to greater than or equal to zero. If it is zero,
438  * then the corresponding term is set to zero.
439  *
440  * m_c is required to greater than or equal to zero. If it is zero,
441  * then the corresponding term is set to zero.
442  *
443  * m_d is required to be greater than zero.
444  *
445  * @ingroup falloffGroup
446  */
447 class SRI5 : public Falloff
448 {
449 
450 public:
451 
452  //! Constructor
453  SRI5() {}
454 
455  //! Destructor
456  virtual ~SRI5() {}
457 
458  //! Initialization of the object
459  /*!
460  * @param c Vector of five doubles: The doubles are the parameters,
461  * a, b, c, d, and e of the SRI parameterization
462  */
463  virtual void init(const vector_fp& c) {
464  if (c[2] < 0.0) {
465  throw CanteraError("SRI5::init()",
466  "m_c parameter is less than zero: " + fp2str(c[2]));
467  }
468  if (c[3] < 0.0) {
469  throw CanteraError("SRI5::init()",
470  "m_d parameter is less than zero: " + fp2str(c[3]));
471  }
472  m_a = c[0];
473  m_b = c[1];
474  m_c = c[2];
475  m_d = c[3];
476  m_e = c[4];
477  }
478 
479  //! Update the temperature parameters in the representation
480  /*!
481  * The workspace has a length of two
482  *
483  * @param T Temperature (Kelvin)
484  * @param work Vector of working space representing
485  * the temperature dependent part of the
486  * parameterization.
487  */
488  virtual void updateTemp(doublereal T, doublereal* work) const {
489  *work = m_a * exp(- m_b / T);
490  if (m_c != 0.0) {
491  *work += exp(- T/m_c);
492  }
493  work[1] = m_d * pow(T,m_e);
494  }
495 
496  //! Function that returns <I>F</I>
497  /*!
498  * @param pr Value of the reduced pressure for this reaction
499  * @param work Pointer to the previously saved work space
500  */
501  virtual doublereal F(doublereal pr, const doublereal* work) const {
502  doublereal lpr = log10(std::max(pr,SmallNumber));
503  doublereal xx = 1.0/(1.0 + lpr*lpr);
504  return pow(*work, xx) * work[1];
505  }
506 
507  //! Utility function that returns the size of the workspace
508  virtual size_t workSize() {
509  return 2;
510  }
511 
512 protected:
513 
514  //! parameter a in the 5-parameter SRI falloff function
515  /*!
516  * This is unitless
517  */
518  doublereal m_a;
519 
520  //! parameter b in the 5-parameter SRI falloff function
521  /*!
522  * This has units of Kelvin
523  */
524  doublereal m_b;
525 
526  //! parameter c in the 5-parameter SRI falloff function
527  /*!
528  * This has units of Kelvin
529  */
530  doublereal m_c;
531 
532  //! parameter d in the 5-parameter SRI falloff function
533  /*!
534  * This is unitless
535  */
536  doublereal m_d;
537 
538  //! parameter d in the 5-parameter SRI falloff function
539  /*!
540  * This is unitless
541  */
542  doublereal m_e;
543 
544 };
545 
546 
547 //! Wang-Frenklach falloff function.
548 /*!
549  *
550  * The falloff function defines the value of \f$ F \f$ in the following
551  * rate expression
552  *
553  * \f[
554  * k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
555  * \f]
556  * where
557  * \f[
558  * P_r = \frac{k_0 [M]}{k_{\infty}}
559  * \f]
560  *
561  * \f[
562  * F = 10.0^{Flog}
563  * \f]
564  * where
565  * \f[
566  * Flog = \frac{\log_{10} F_{cent}}{\exp{(\frac{\log_{10} P_r - \alpha}{\sigma})^2}}
567  * \f]
568  * where
569  *
570  * \f[
571  * F_{cent} = (1 - A)\exp(-T/T_3) + A \exp(-T/T_1) + \exp(-T/T_2)
572  * \f]
573  *
574  * \f[
575  * \alpha = \alpha_0 + \alpha_1 T + \alpha_2 T^2
576  * \f]
577  *
578  * \f[
579  * \sigma = \sigma_0 + \sigma_1 T + \sigma_2 T^2
580  * \f]
581  *
582  *
583  * Reference: Wang, H., and
584  * Frenklach, M., Chem. Phys. Lett. vol. 205, 271 (1993).
585  *
586  *
587  * @ingroup falloffGroup
588  */
589 class WF93 : public Falloff
590 {
591 
592 public:
593 
594  //! Default constructor
595  WF93() {}
596 
597  //! Destructor
598  virtual ~WF93() {}
599 
600  //! Initialization routine
601  /*!
602  * @param c Vector of 10 doubles
603  * with the following ordering:
604  * a, T_1, T_2, T_3, alpha0, alpha1, alpha2
605  * sigma0, sigma1, sigma2
606  */
607  virtual void init(const vector_fp& c) {
608  m_a = c[0];
609  m_rt1 = 1.0/c[1];
610  m_t2 = c[2];
611  m_rt3 = 1.0/c[3];
612  m_alpha0 = c[4];
613  m_alpha1 = c[5];
614  m_alpha2 = c[6];
615  m_sigma0 = c[7];
616  m_sigma1 = c[8];
617  m_sigma2 = c[9];
618  }
619 
620  //! Update the temperature parameters in the representation
621  /*!
622  * The workspace has a length of three
623  *
624  * @param T Temperature (Kelvin)
625  * @param work Vector of working space representing
626  * the temperature dependent part of the
627  * parameterization.
628  */
629  virtual void updateTemp(doublereal T, doublereal* work) const {
630  work[0] = m_alpha0 + (m_alpha1 + m_alpha2*T)*T; // alpha
631  work[1] = m_sigma0 + (m_sigma1 + m_sigma2*T)*T; // sigma
632  doublereal Fcent = (1.0 - m_a) * exp(- T * m_rt3)
633  + m_a * exp(- T * m_rt1) + exp(-m_t2/T);
634  work[2] = log10(Fcent);
635  }
636 
637  //! Function that returns <I>F</I>
638  /*!
639  * @param pr Value of the reduced pressure for this reaction
640  * @param work Pointer to the previously saved work space
641  */
642  virtual doublereal F(doublereal pr, const doublereal* work) const {
643  doublereal lpr = log10(std::max(pr, SmallNumber));
644  doublereal x = (lpr - work[0])/work[1];
645  doublereal flog = work[2]/exp(x*x);
646  return pow(10.0, flog);
647  }
648 
649  //! Utility function that returns the size of the workspace
650  virtual size_t workSize() {
651  return 3;
652  }
653 
654 protected:
655 
656  //! Value of the \f$ \alpha_0 \f$ coefficient
657  /*!
658  * This is the fifth coefficient in the xml list
659  */
660  doublereal m_alpha0;
661 
662  //! Value of the \f$ \alpha_1 \f$ coefficient
663  /*!
664  * This is the 6th coefficient in the xml list
665  */
666  doublereal m_alpha1;
667 
668  //! Value of the \f$ \alpha_2 \f$ coefficient
669  /*!
670  * This is the 7th coefficient in the xml list
671  */
672  doublereal m_alpha2;
673 
674  //! Value of the \f$ \sigma_0 \f$ coefficient
675  /*!
676  * This is the 8th coefficient in the xml list
677  */
678  doublereal m_sigma0;
679 
680  //! Value of the \f$ \sigma_1 \f$ coefficient
681  /*!
682  * This is the 9th coefficient in the xml list
683  */
684  doublereal m_sigma1;
685 
686  //! Value of the \f$ \sigma_2 \f$ coefficient
687  /*!
688  * This is the 10th coefficient in the xml list
689  */
690  doublereal m_sigma2;
691 
692  //! Value of the \f$ a \f$ coefficient
693  /*!
694  * This is the first coefficient in the xml list
695  */
696  doublereal m_a;
697 
698  //! Value of inverse of the \f$ t1 \f$ coefficient
699  /*!
700  * This is the second coefficient in the xml list
701  */
702  doublereal m_rt1;
703 
704  //! Value of the \f$ t2 \f$ coefficient
705  /*!
706  * This is the third coefficient in the xml list
707  */
708  doublereal m_t2;
709 
710  //! Value of the inverse of the \f$ t3 \f$ coefficient
711  /*!
712  * This is the 4th coefficient in the xml list
713  */
714  doublereal m_rt3;
715 
716 private:
717 
718 };
719 
720 // Factory routine that returns a new Falloff parameterization object
721 /*
722  * @param type Integer type of the falloff parameterization. These
723  * integers are listed in reaction_defs.h
724  *
725  * @param c Vector of input parameterizations for the Falloff
726  * object. The function is initialized with this vector.
727  *
728  * @return Returns a pointer to a newly malloced Falloff object
729  */
730 Falloff* FalloffFactory::newFalloff(int type, const vector_fp& c)
731 {
732  Falloff* f;
733  switch (type) {
734  case TROE3_FALLOFF:
735  f = new Troe3();
736  break;
737  case TROE4_FALLOFF:
738  f = new Troe4();
739  break;
740  case SRI3_FALLOFF:
741  f = new SRI3();
742  break;
743  case SRI5_FALLOFF:
744  f = new SRI5();
745  break;
746  case WF_FALLOFF:
747  f = new WF93();
748  break;
749  default:
750  return 0;
751  }
752  f->init(c);
753  return f;
754 }
755 
756 }
757 
Generated by   doxygen 1.8.2