1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // Antioch - A Gas Dynamics Thermochemistry Library
5 //
6 // Copyright (C) 2014-2016 Paul T. Bauman, Benjamin S. Kirk,
7 //                         Sylvain Plessis, Roy H. Stonger
8 //
9 // Copyright (C) 2013 The PECOS Development Team
10 //
11 // This library is free software; you can redistribute it and/or
12 // modify it under the terms of the Version 2.1 GNU Lesser General
13 // Public License as published by the Free Software Foundation.
14 //
15 // This library is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
23 // Boston, MA  02110-1301  USA
24 //
25 //-----------------------------------------------------------------------el-
26 
27 
28 #ifndef ANTIOCH_THREEBODY_REACTION_H
29 #define ANTIOCH_THREEBODY_REACTION_H
30 
31 // Antioch
32 #include "antioch/reaction.h"
33 #include "antioch/kinetics_conditions.h"
34 
35 //C++
36 #include <string>
37 #include <vector>
38 #include <iostream>
39 
40 namespace Antioch
41 {
42   //!A single reaction mechanism.
43   /*!\class ThreeBodyReaction
44  *
45     This class encapsulates a three-body reaction process. A three-body process
46     rate constant is defined by the equation
47     \f[
48         k(T,[M]) = \alpha(T)\times \sum_i\epsilon_ic_i
49     \f]
50     with \f$\alpha(T)\f$ being a kinetics model (see base classes Reaction and KineticsType), \f$[M]\f$
51     the mixture concentration (or pressure, it's equivalent, \f$[M] = \frac{P}{\mathrm{R}T}\f$
52     in an ideal gas model) and \f$c_i\f$ the concentration of species \f$i\f$.  All reactions are assumed to be reversible.
53     By default, the KooijRate kinetics model is used.
54 
55     We have:
56     \f[
57         \begin{split}
58            \frac{\partial k(T,[M])}{\partial T}   & = \frac{\partial \alpha(T)}{\partial T} \sum_i\epsilon_ic_i \\[10pt]
59            \frac{\partial k(T,[M])}{\partial c_i} & = \alpha(T) \epsilon_i
60         \end{split}
61     \f]
62     with \f$c_i\f$ the concentration of species \f$i\f$.
63   */
64   template <typename CoeffType=double>
65   class ThreeBodyReaction: public Reaction<CoeffType>
66   {
67   public:
68 
69     //! Construct a single reaction mechanism.
70     ThreeBodyReaction( const unsigned int n_species,
71                        const std::string &equation,
72                        const bool &reversible = true,
73                        const KineticsModel::KineticsModel kin = KineticsModel::KOOIJ);
74 
75     ~ThreeBodyReaction();
76 
77     //!
78     template <typename StateType, typename VectorStateType>
79     StateType compute_forward_rate_coefficient( const VectorStateType& molar_densities,
80                                                 const KineticsConditions<StateType,VectorStateType>& conditions ) const;
81 
82     //!
83     template <typename StateType, typename VectorStateType>
84     void compute_forward_rate_coefficient_and_derivatives( const VectorStateType& molar_densities,
85                                                            const KineticsConditions<StateType,VectorStateType>& conditions,
86                                                            StateType& kfwd,
87                                                            StateType& dkfwd_dT,
88                                                            VectorStateType& dkfwd_dX) const;
89 
90   };
91 
92   /* ------------------------- Inline Functions -------------------------*/
93   template <typename CoeffType>
94   inline
ThreeBodyReaction(const unsigned int n_species,const std::string & equation,const bool & reversible,const KineticsModel::KineticsModel kin)95   ThreeBodyReaction<CoeffType>::ThreeBodyReaction( const unsigned int n_species,
96                                                    const std::string &equation ,
97                                                    const bool &reversible,
98                                                    const KineticsModel::KineticsModel kin)
99     :Reaction<CoeffType>(n_species,equation,reversible,ReactionType::THREE_BODY,kin)
100   {
101     Reaction<CoeffType>::_efficiencies.resize(n_species);
102     std::fill (Reaction<CoeffType>::_efficiencies.begin(), Reaction<CoeffType>::_efficiencies.end(), 1.);
103     return;
104   }
105 
106   template <typename CoeffType>
107   inline
~ThreeBodyReaction()108   ThreeBodyReaction<CoeffType>::~ThreeBodyReaction()
109   {
110     return;
111   }
112 
113 
114 
115   template <typename CoeffType>
116   template<typename StateType, typename VectorStateType>
117   inline
compute_forward_rate_coefficient(const VectorStateType & molar_densities,const KineticsConditions<StateType,VectorStateType> & conditions)118   StateType ThreeBodyReaction<CoeffType>::compute_forward_rate_coefficient( const VectorStateType& molar_densities,
119                                                                             const KineticsConditions<StateType,VectorStateType>& conditions  ) const
120   {
121     //k(T,[M]) = (sum eff_i * C_i) * ...
122     StateType kfwd = (this->efficiency(0) * molar_densities[0] );
123 
124     for (unsigned int s=1; s<this->n_species(); s++)
125       {
126         kfwd += ( this->efficiency(s) * molar_densities[s] );
127       }
128 
129     //... alpha(T)
130     kfwd *= (*this->_forward_rate[0])(conditions);
131 
132     antioch_assert(!has_nan(kfwd));
133 
134     return kfwd;
135   }
136 
137 
138   template <typename CoeffType>
139   template<typename StateType, typename VectorStateType>
140   inline
compute_forward_rate_coefficient_and_derivatives(const VectorStateType & molar_densities,const KineticsConditions<StateType,VectorStateType> & conditions,StateType & kfwd,StateType & dkfwd_dT,VectorStateType & dkfwd_dX)141   void ThreeBodyReaction<CoeffType>::compute_forward_rate_coefficient_and_derivatives( const VectorStateType &molar_densities,
142                                                                                        const KineticsConditions<StateType,VectorStateType>& conditions,
143                                                                                        StateType& kfwd,
144                                                                                        StateType& dkfwd_dT,
145                                                                                        VectorStateType& dkfwd_dX) const
146   {
147     antioch_assert_equal_to(dkfwd_dX.size(),this->n_species());
148 
149     //dk_dT = dalpha_dT * [sum_s (eps_s * X_s)]
150     //dk_dCi = alpha(T) * eps_i
151     this->_forward_rate[0]->compute_rate_and_derivative(conditions,kfwd,dkfwd_dT);
152 
153     dkfwd_dX[0] = kfwd;
154     StateType coef = (this->efficiency(0) * molar_densities[0]);
155 
156     for (unsigned int s=1; s<this->n_species(); s++)
157       {
158         coef += ( this->efficiency(s) * molar_densities[s] );
159         dkfwd_dX[s] = kfwd;
160       }
161 
162     kfwd *= coef;
163     dkfwd_dT *= coef;
164 
165     for (unsigned int s=0; s<this->n_species(); s++)
166       {
167         dkfwd_dX[s] *= this->efficiency(s);
168       }
169 
170     antioch_assert(!has_nan(kfwd));
171 
172     return;
173   }
174 
175 } // namespace Antioch
176 
177 #endif // ANTIOCH_ELEMENTARY_REACTION_H
178