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