# FalloffFactory.cpp

Go to the documentation of this file.
00001 /**
00002  *  @file FalloffFactory.cpp
00003  */
00004
00005 /*  $Author: hkmoffa$
00006  *  $Date: 2009-12-10 18:20:44 -0500 (Thu, 10 Dec 2009)$
00007  *  $Revision: 309$
00008  */
00009
00010 // Copyright 2001  California Institute of Technology
00011
00012 #ifdef WIN32
00013 #pragma warning(disable:4786)
00014 #pragma warning(disable:4503)
00015 #endif
00016
00017
00018 #include "FalloffFactory.h"
00019 #include "ctexceptions.h"
00020
00021 #include <cmath>
00022
00023 namespace Cantera {
00024
00025   FalloffFactory* FalloffFactory::s_factory = 0;
00027   boost::mutex FalloffFactory::falloff_mutex ;
00028 #endif
00029
00030
00031   //! The 3-parameter Troe falloff parameterization.
00032   /*!
00033    * The falloff function defines the value of \f$F \f$ in the following
00034    * rate expression
00035    *
00036    *  \f[
00037    *         k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
00038    *  \f]
00039    *  where
00040    *  \f[
00041    *         P_r = \frac{k_0 [M]}{k_{\infty}}
00042    *  \f]
00043    *
00044    * This parameterization is  defined by
00045    * \f[
00046    *        F = F_{cent}^{1/(1 + f_1^2)}
00047    * \f]
00048    * where
00049    * \f[
00050    *       F_{cent} = (1 - A)\exp(-T/T_3) + A \exp(-T/T_1)
00051    * \f]
00052    *
00053    * \f[
00054    *       f_1 = (\log_{10} P_r + C) / \left(N - 0.14
00055    *             (\log_{10} P_r + C)\right)
00056    * \f]
00057    *
00058    * \f[
00059    *          C = -0.4 - 0.67 \log_{10} F_{cent}
00060    * \f]
00061    *
00062    * \f[
00063    *          N = 0.75 - 1.27 \log_{10} F_{cent}
00064    * \f]
00065    *
00066    *  There are a few requirements for the parameters
00067    *
00068    *   T_3 is required to greater than or equal to zero. If it is zero,
00069    *   then the term is set to zero.
00070    *
00071    *   T_1 is required to greater than or equal to zero. If it is zero,
00072    *   then the term is set to zero.
00073    *
00074    * @ingroup falloffGroup
00075    */
00076   class Troe3 : public Falloff {
00077   public:
00078
00079     //! Default constructor.
00080     Troe3() : m_a (0.0), m_rt3 (0.0), m_rt1 (0.0) {}
00081
00082     //! Destructor. Does nothing.
00083     virtual ~Troe3() {}
00084
00085     /**
00086      * Initialize.
00087      * @param c Coefficient vector of length 3,
00088      * with entries \f$(A, T_3, T_1) \f$
00089      */
00090     virtual void init(const vector_fp& c) {
00091       m_a  = c[0];
00092
00093       if (c[1] <= 0.0) {
00094         if (c[1] == 0.0) {
00095           m_rt3 = 1000.;
00096         } else {
00097           throw CanteraError("Troe3::init()", "T3 parameter is less than zero");
00098         }
00099       } else {
00100         m_rt3 = 1.0/c[1];
00101       }
00102       if (c[2] <= 0.0) {
00103         if (c[2] == 0.0) {
00104           m_rt1 = 1000.;
00105         } else {
00106           throw CanteraError("Troe3::init()", "T1 parameter is less than zero");
00107         }
00108       } else {
00109         m_rt1 = 1.0/c[2];
00110       }
00111     }
00112
00113     //! Update the temperature parameters in the representation
00114     /*!
00115      *   The workspace has a length of one
00116      *
00117      *   @param T         Temperature (Kelvin)
00118      *   @param work      Vector of working space representing
00119      *                    the temperature dependent part of the
00120      *                    parameterization.
00121      */
00122     virtual void updateTemp(doublereal T, workPtr work) const {
00123       doublereal Fcent = (1.0 - m_a) * exp(- T * m_rt3 )
00124         + m_a * exp(- T * m_rt1 );
00125       *work = log10( fmaxx( Fcent, SmallNumber ) );
00126     }
00127
00128     //! Function that returns <I>F</I>
00129     /*!
00130      *  @param pr   Value of the reduced pressure for this reaction
00131      *  @param work Pointer to the previously saved work space
00132      */
00133     virtual doublereal F(doublereal pr, const_workPtr work) const {
00134       doublereal lpr,f1,lgf, cc, nn;
00135       lpr = log10( fmaxx(pr,SmallNumber) );
00136       cc = -0.4 - 0.67 * (*work);
00137       nn = 0.75 - 1.27 * (*work);
00138       f1 = ( lpr + cc )/ ( nn - 0.14 * ( lpr + cc ) );
00139       lgf = (*work) / ( 1.0 + f1 * f1 );
00140       return pow(10.0, lgf );
00141     }
00142
00143     //! Utility function that returns the size of the workspace
00144     virtual size_t workSize() { return 1; }
00145
00146   protected:
00147
00148     //! parameter a in the  4-parameter Troe falloff function
00149     /*!
00150      * This is unitless
00151      */
00152     doublereal m_a;
00153
00154     //! parameter 1/T_3 in the  4-parameter Troe falloff function
00155     /*!
00156      * This has units of Kelvin-1
00157      */
00158     doublereal m_rt3;
00159
00160     //! parameter 1/T_1 in the  4-parameter Troe falloff function
00161     /*!
00162      * This has units of Kelvin-1
00163      */
00164     doublereal m_rt1;
00165
00166   };
00167
00168   //! The 4-parameter Troe falloff parameterization.
00169   /*!
00170    * The falloff function defines the value of \f$F \f$ in the following
00171    * rate expression
00172    *
00173    *  \f[
00174    *         k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
00175    *  \f]
00176    *  where
00177    *  \f[
00178    *         P_r = \frac{k_0 [M]}{k_{\infty}}
00179    *  \f]
00180    *
00181    * This parameterization is defined by
00182    *
00183    * \f[
00184    *         F = F_{cent}^{1/(1 + f_1^2)}
00185    * \f]
00186    *    where
00187    * \f[
00188    *         F_{cent} = (1 - A)\exp(-T/T_3) + A \exp(-T/T_1) + \exp(-T_2/T)
00189    * \f]
00190    *
00191    * \f[
00192    *        f_1 = (\log_{10} P_r + C) /
00193    *              \left(N - 0.14 (\log_{10} P_r + C)\right)
00194    * \f]
00195    *
00196    * \f[
00197    *         C = -0.4 - 0.67 \log_{10} F_{cent}
00198    * \f]
00199    *
00200    * \f[
00201    *          N = 0.75 - 1.27 \log_{10} F_{cent}
00202    * \f]
00203    *
00204    *
00205    *  There are a few requirements for the parameters
00206    *
00207    *   T_3 is required to greater than or equal to zero. If it is zero,
00208    *   then the term is set to zero.
00209    *
00210    *   T_1 is required to greater than or equal to zero. If it is zero,
00211    *   then the term is set to zero.
00212    *
00213    *   T_2 is required to be greater than  zero.
00214    *
00215    * @ingroup falloffGroup
00216    */
00217   class Troe4 : public Falloff {
00218   public:
00219     //! Constructor
00220     Troe4() : m_a (0.0), m_rt3 (0.0), m_rt1 (0.0),
00221               m_t2 (0.0) {}
00222
00223     //! Destructor
00224     virtual ~Troe4() {}
00225
00226
00227     //! Initialization of the object
00228     /*!
00229      * @param c Vector of four doubles: The doubles are the parameters,
00230      *          a,, T_3, T_1, and T_2 of the SRI parameterization
00231      */
00232     virtual void init(const vector_fp& c) {
00233       m_a  = c[0];
00234       if (c[1] <= 0.0) {
00235         if (c[1] == 0.0) {
00236           m_rt3 = 1000.;
00237         } else {
00238           throw CanteraError("Troe4::init()", "T3 parameter is less than zero");
00239         }
00240       } else {
00241         m_rt3 = 1.0/c[1];
00242       }
00243       if (c[2] <= 0.0) {
00244         if (c[2] == 0.0) {
00245           m_rt1 = 1000.;
00246         } else {
00247           throw CanteraError("Troe4::init()", "T1 parameter is less than zero");
00248         }
00249       } else {
00250         m_rt1 = 1.0/c[2];
00251       }
00252       if (c[3] < 0.0) {
00253         throw CanteraError("Troe4::init()", "T2 parameter is less than zero");
00254       }
00255       m_t2 = c[3];
00256     }
00257
00258     //! Update the temperature parameters in the representation
00259     /*!
00260      *   The workspace has a length of one
00261      *
00262      *   @param T         Temperature (Kelvin)
00263      *   @param work      Vector of working space representing
00264      *                    the temperature dependent part of the
00265      *                    parameterization.
00266      */
00267     virtual void updateTemp(doublereal T, workPtr work) const {
00268       doublereal Fcent = (1.0 - m_a) * exp(- T * m_rt3 )
00269         + m_a * exp(- T * m_rt1 )
00270         + exp(- m_t2 / T );
00271       *work = log10( fmaxx( Fcent, SmallNumber ) );
00272     }
00273
00274     //! Function that returns <I>F</I>
00275     /*!
00276      *  @param pr   Value of the reduced pressure for this reaction
00277      *  @param work Pointer to the previously saved work space
00278      */
00279     virtual doublereal F(doublereal pr, const_workPtr work) const {
00280       doublereal lpr,f1,lgf, cc, nn;
00281       lpr = log10( fmaxx(pr,SmallNumber) );
00282       cc = -0.4 - 0.67 * (*work);
00283       nn = 0.75 - 1.27 * (*work);
00284       f1 = ( lpr + cc )/ ( nn - 0.14 * ( lpr + cc ) );
00285       lgf = (*work) / ( 1.0 + f1 * f1 );
00286       return pow(10.0, lgf );
00287     }
00288
00289     //! Utility function that returns the size of the workspace
00290     virtual size_t workSize() { return 1; }
00291
00292   protected:
00293
00294     //! parameter a in the  4-parameter Troe falloff function
00295     /*!
00296      * This is unitless
00297      */
00298     doublereal m_a;
00299
00300     //! parameter 1/T_3 in the  4-parameter Troe falloff function
00301     /*!
00302      * This has units of Kelvin-1
00303      */
00304     doublereal m_rt3;
00305
00306     //! parameter 1/T_1 in the  4-parameter Troe falloff function
00307     /*!
00308      * This has units of Kelvin-1
00309      */
00310     doublereal m_rt1;
00311
00312     //! parameter T_2 in the  4-parameter Troe falloff function
00313     /*!
00314      * This has units of Kelvin
00315      */
00316     doublereal m_t2;
00317   };
00318
00319
00320   //! The 3-parameter SRI falloff function for <I>F</I>
00321   /*!
00322    * The falloff function defines the value of \f$F \f$ in the following
00323    * rate expression
00324    *
00325    *  \f[
00326    *         k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
00327    *  \f]
00328    *  where
00329    *  \f[
00330    *         P_r = \frac{k_0 [M]}{k_{\infty}}
00331    *  \f]
00332    *
00333    *  \f[
00334    *         F = {\left( a \; exp(\frac{-b}{T}) + exp(\frac{-T}{c})\right)}^n
00335    *  \f]
00336    *      where
00337    *  \f[
00338    *         n = \frac{1.0}{1.0 + {\log_{10} P_r}^2}
00339    *  \f]
00340    *
00341    *   \f$c \f$ s required to greater than or equal to zero. If it is zero,
00342    *   then the corresponding term is set to zero.
00343    *
00344    * @ingroup falloffGroup
00345    */
00346   class SRI3 : public Falloff {
00347
00348   public:
00349
00350     //! Constructor
00351     SRI3() {}
00352
00353     //! Destructor
00354     virtual ~SRI3() {}
00355
00356     //! Initialization of the object
00357     /*!
00358      * @param c Vector of three doubles: The doubles are the parameters,
00359      *          a, b, and c of the SRI parameterization
00360      */
00361     virtual void init(const vector_fp& c) {
00362       m_a = c[0];
00363       m_b = c[1];
00364       m_c = c[2];
00365     }
00366
00367     //! Update the temperature parameters in the representation
00368     /*!
00369      *   The workspace has a length of one
00370      *
00371      *   @param T         Temperature (Kelvin)
00372      *   @param work      Vector of working space representing
00373      *                    the temperature dependent part of the
00374      *                    parameterization.
00375      */
00376     virtual void updateTemp(doublereal T, workPtr work) const {
00377       *work = m_a * exp( - m_b / T);
00378       if (m_c != 0.0) *work += exp( - T/m_c );
00379     }
00380
00381     //! Function that returns <I>F</I>
00382     /*!
00383      *  @param pr   Value of the reduced pressure for this reaction
00384      *  @param work Pointer to the previously saved work space
00385      */
00386     virtual doublereal F(doublereal pr, const_workPtr work) const {
00387       doublereal lpr = log10( fmaxx(pr,SmallNumber) );
00388       doublereal xx = 1.0/(1.0 + lpr*lpr);
00389       doublereal ff = pow( *work , xx);
00390       return ff;
00391     }
00392
00393     //! Utility function that returns the size of the workspace
00394     virtual size_t workSize() { return 1; }
00395
00396   protected:
00397
00398     //! parameter a in the  3-parameter SRI falloff function
00399     /*!
00400      * This is unitless
00401      */
00402     doublereal m_a;
00403
00404     //! parameter b in the  3-parameter SRI falloff function
00405     /*!
00406      * This has units of Kelvin
00407      */
00408     doublereal m_b;
00409
00410     //! parameter c in the  3-parameter SRI falloff function
00411     /*!
00412      * This has units of Kelvin
00413      */
00414     doublereal m_c;
00415   };
00416
00417
00418   //! The 5-parameter SRI falloff function.
00419   /*!
00420    * The falloff function defines the value of \f$F \f$ in the following
00421    * rate expression
00422    *
00423    *  \f[
00424    *         k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
00425    *  \f]
00426    *  where
00427    *  \f[
00428    *         P_r = \frac{k_0 [M]}{k_{\infty}}
00429    *  \f]
00430    *
00431    *  \f[
00432    *         F = {\left( a \; exp(\frac{-b}{T}) + exp(\frac{-T}{c})\right)}^n
00433    *              \;  d \; exp(\frac{-e}{T})
00434    *  \f]
00435    *      where
00436    *  \f[
00437    *         n = \frac{1.0}{1.0 + {\log_{10} P_r}^2}
00438    *  \f]
00439    *
00440    *   \f$c \f$ s required to greater than or equal to zero. If it is zero,
00441    *   then the corresponding term is set to zero.
00442    *
00443    *   m_c is required to greater than or equal to zero. If it is zero,
00444    *   then the corresponding term is set to zero.
00445    *
00446    *   m_d is required to be greater than zero.
00447    *
00448    * @ingroup falloffGroup
00449    */
00450   class SRI5 : public Falloff {
00451
00452   public:
00453
00454     //! Constructor
00455     SRI5() {}
00456
00457     //! Destructor
00458     virtual ~SRI5() {}
00459
00460     //! Initialization of the object
00461     /*!
00462      * @param c Vector of five doubles: The doubles are the parameters,
00463      *          a, b, c, d, and e of the SRI parameterization
00464      */
00465     virtual void init(const vector_fp& c) {
00466       m_a = c[0];
00467       m_b = c[1];
00468       m_c = c[2];
00469       m_d = c[3];
00470       m_e = c[4];
00471     }
00472
00473     //! Update the temperature parameters in the representation
00474     /*!
00475      *   The workspace has a length of two
00476      *
00477      *   @param T         Temperature (Kelvin)
00478      *   @param work      Vector of working space representing
00479      *                    the temperature dependent part of the
00480      *                    parameterization.
00481      */
00482     virtual void updateTemp(doublereal T, workPtr work) const {
00483       *work = m_a * exp( - m_b / T);
00484       if (m_c != 0.0) *work += exp( - T/m_c );
00485       work[1] = m_d * pow(T,m_e);
00486     }
00487
00488     //! Function that returns <I>F</I>
00489     /*!
00490      *  @param pr   Value of the reduced pressure for this reaction
00491      *  @param work Pointer to the previously saved work space
00492      */
00493     virtual doublereal F(doublereal pr, const_workPtr work) const {
00494       doublereal lpr = log10( fmaxx(pr,SmallNumber) );
00495       doublereal xx = 1.0/(1.0 + lpr*lpr);
00496       return pow( *work, xx) * work[1];
00497     }
00498
00499     //! Utility function that returns the size of the workspace
00500     virtual size_t workSize() { return 2; }
00501
00502   protected:
00503
00504     //! parameter a in the  5-parameter SRI falloff function
00505     /*!
00506      * This is unitless
00507      */
00508     doublereal m_a;
00509
00510     //! parameter b in the  5-parameter SRI falloff function
00511     /*!
00512      * This has units of Kelvin
00513      */
00514     doublereal m_b;
00515
00516     //! parameter c in the  5-parameter SRI falloff function
00517     /*!
00518      * This has units of Kelvin
00519      */
00520     doublereal m_c;
00521
00522     //! parameter d in the  5-parameter SRI falloff function
00523     /*!
00524      * This is unitless
00525      */
00526     doublereal m_d;
00527
00528     //! parameter d in the  5-parameter SRI falloff function
00529     /*!
00530      * This is unitless
00531      */
00532     doublereal m_e;
00533
00534   };
00535
00536
00537   //! Wang-Frenklach falloff function.
00538   /*!
00539    *
00540    * The falloff function defines the value of \f$F \f$ in the following
00541    * rate expression
00542    *
00543    *  \f[
00544    *         k = k_{\infty} \left( \frac{P_r}{1 + P_r} \right) F
00545    *  \f]
00546    *  where
00547    *  \f[
00548    *         P_r = \frac{k_0 [M]}{k_{\infty}}
00549    *  \f]
00550    *
00551    *  \f[
00552    *         F = 10.0^{Flog}
00553    *  \f]
00554    *  where
00555    *  \f[
00556    *         Flog = \frac{\log_{10} F_{cent}}{\exp{(\frac{\log_{10} P_r - \alpha}{\sigma})^2}}
00557    *  \f]
00558    *    where
00559    *
00560    *  \f[
00561    *       F_{cent} = (1 - A)\exp(-T/T_3) + A \exp(-T/T_1) + \exp(-T/T_2)
00562    *  \f]
00563    *
00564    *  \f[
00565    *      \alpha = \alpha_0 + \alpha_1 T + \alpha_2 T^2
00566    *  \f]
00567    *
00568    *  \f[
00569    *      \sigma = \sigma_0 + \sigma_1 T + \sigma_2 T^2
00570    *  \f]
00571    *
00572    *
00573    * Reference: Wang, H., and
00574    *            Frenklach, M., Chem. Phys. Lett. vol. 205, 271 (1993).
00575    *
00576    *
00577    * @ingroup falloffGroup
00578    */
00579   class WF93 : public Falloff {
00580
00581   public:
00582
00583     //! Default constructpr
00584     WF93() {}
00585
00586     //! Destructor
00587     virtual ~WF93() {}
00588
00589     //! Initialization routine
00590     /*!
00591      *  @param c  Vector of 10 doubles
00592      *            with the following ordering:
00593      *               a, T_1, T_2, T_3, alpha0, alpha1, alpha2
00594      *               sigma0, sigma1, sigma2
00595      */
00596     virtual void init(const vector_fp& c) {
00597       m_a = c[0];
00598       m_rt1 = 1.0/c[1];
00599       m_t2 = c[2];
00600       m_rt3  = 1.0/c[3];
00601       m_alpha0 = c[4];
00602       m_alpha1 = c[5];
00603       m_alpha2 = c[6];
00604       m_sigma0 = c[7];
00605       m_sigma1 = c[8];
00606       m_sigma2 = c[9];
00607     }
00608
00609     //! Update the temperature parameters in the representation
00610     /*!
00611      *   The workspace has a length of three
00612      *
00613      *   @param T         Temperature (Kelvin)
00614      *   @param work      Vector of working space representing
00615      *                    the temperature dependent part of the
00616      *                    parameterization.
00617      */
00618     virtual void updateTemp(doublereal T, workPtr work) const {
00619       work[0] = m_alpha0 + (m_alpha1 + m_alpha2*T)*T; // alpha
00620       work[1] = m_sigma0 + (m_sigma1 + m_sigma2*T)*T; // sigma
00621       doublereal Fcent = (1.0 - m_a) * exp(- T * m_rt3 )
00622         + m_a * exp(- T * m_rt1 ) + exp(-m_t2/T);
00623       work[2] = log10(Fcent);
00624     }
00625
00626     //! Function that returns <I>F</I>
00627     /*!
00628      *  @param pr   Value of the reduced pressure for this reaction
00629      *  @param work Pointer to the previously saved work space
00630      */
00631     virtual doublereal F(doublereal pr, const_workPtr work) const {
00632       doublereal lpr = log10( fmaxx(pr, SmallNumber) );
00633       doublereal x = (lpr - work[0])/work[1];
00634       doublereal flog = work[2]/exp(x*x);
00635       return pow( 10.0, flog);
00636     }
00637
00638     //! Utility function that returns the size of the workspace
00639     virtual size_t workSize() { return 3; }
00640
00641   protected:
00642
00643     //! Value of the \f$\alpha_0 \f$ coefficient
00644     /*!
00645      *  This is the fifth coefficient in the xml list
00646      */
00647     doublereal m_alpha0;
00648
00649     //! Value of the \f$\alpha_1 \f$ coefficient
00650     /*!
00651      *  This is the 6th coefficient in the xml list
00652      */
00653     doublereal m_alpha1;
00654
00655     //! Value of the \f$\alpha_2 \f$ coefficient
00656     /*!
00657      *  This is the 7th coefficient in the xml list
00658      */
00659     doublereal m_alpha2;
00660
00661     //! Value of the \f$\sigma_0 \f$ coefficient
00662     /*!
00663      *  This is the 8th coefficient in the xml list
00664      */
00665     doublereal m_sigma0;
00666
00667     //! Value of the \f$\sigma_1 \f$ coefficient
00668     /*!
00669      *  This is the 9th coefficient in the xml list
00670      */
00671     doublereal m_sigma1;
00672
00673     //! Value of the \f$\sigma_2 \f$ coefficient
00674     /*!
00675      *  This is the 10th coefficient in the xml list
00676      */
00677     doublereal m_sigma2;
00678
00679     //! Value of the \f$a \f$ coefficient
00680     /*!
00681      *  This is the first coefficient in the xml list
00682      */
00683     doublereal m_a;
00684
00685     //! Value of inverse of the \f$t1 \f$ coefficient
00686     /*!
00687      *  This is the second coefficient in the xml list
00688      */
00689     doublereal m_rt1;
00690
00691     //! Value of the \f$t2 \f$ coefficient
00692     /*!
00693      *  This is the third coefficient in the xml list
00694      */
00695     doublereal m_t2;
00696
00697     //! Value of the inverse of the  \f$t3 \f$ coefficient
00698     /*!
00699      *  This is the 4th coefficient in the xml list
00700      */
00701     doublereal m_rt3;
00702
00703   private:
00704
00705   };
00706
00707   // Factory routine that returns a new Falloff parameterization object
00708   /*
00709    *  @param type  Integer type of the falloff parameterization. These
00710    *               integers are listed in reaction_defs.h
00711    *
00712    *  @param c     Vector of input parameterizations for the Falloff
00713    *               object. The function is initialized with this vector.
00714    *
00715    *  @return  Returns a pointer to a newly malloced Falloff object
00716    */
00717   Falloff* FalloffFactory::newFalloff(int type, const vector_fp& c) {
00718     Falloff* f;
00719     switch(type) {
00720     case TROE3_FALLOFF:
00721       f = new Troe3(); break;
00722     case TROE4_FALLOFF:
00723       f = new Troe4(); break;
00724     case SRI3_FALLOFF:
00725       f = new SRI3(); break;
00726     case SRI5_FALLOFF:
00727       f = new SRI5(); break;
00728     case WF_FALLOFF:
00729       f = new WF93(); break;
00730     default: return 0;
00731     }
00732     f->init(c);
00733     return f;
00734   }
00735
00736 }
00737

Generated by  1.6.3