antioch-0.4.0
Functions
lindemann_falloff_unit.C File Reference
#include <limits>
#include "antioch/vector_utils.h"
#include "antioch/reaction.h"
#include "antioch/falloff_reaction.h"

Go to the source code of this file.

Functions

template<typename Scalar >
int tester (const std::string &type)
 
int main ()
 

Function Documentation

int main ( )

Definition at line 241 of file lindemann_falloff_unit.C.

242 {
243  return (tester<double>("double") ||
244  tester<long double>("long double") ||
245  tester<float>("float"));
246 }
template<typename Scalar >
int tester ( const std::string &  type)

Definition at line 36 of file lindemann_falloff_unit.C.

References Antioch::Reaction< CoeffType >::add_forward_rate(), Antioch::KineticsModel::ARRHENIUS, Antioch::KineticsModel::BERTHELOT, Antioch::KineticsModel::BHE, Antioch::FalloffReaction< CoeffType, FalloffType >::compute_forward_rate_coefficient(), Antioch::FalloffReaction< CoeffType, FalloffType >::compute_forward_rate_coefficient_and_derivatives(), Antioch::KineticsModel::HERCOURT_ESSEN, Antioch::KineticsModel::KOOIJ, Antioch::ReactionType::LINDEMANN_FALLOFF, std::pow(), and Antioch::KineticsModel::VANTHOFF.

37 {
38  using std::abs;
39  using std::exp;
40  using std::pow;
41 
42 
43 //values for 2 CH3 (+M) <=> C2H6 (+M) for the Kooij model, Ds are made up
44 
45  const Scalar Cf1 = 1.135e36L * 1e6L * 1e-12L; //(cm3/mol)^2/s -> kmol -> m3
46  const Scalar beta1 = 1.246L; //true value is -5.246
47  const Scalar Ea1 = 1704.8L / 1.9858775L; //cal/mol
48  const Scalar D1 = -4e-2L; // K^-1
49  const Scalar Cf2 = 6.22e16L * 1e3L * 1e-12L; //cm3/mol/s -> kmol -> m3
50  const Scalar beta2 = -1.174L;
51  const Scalar Ea2 = 635.8L / 1.9858775L; //cal/mol
52  const Scalar D2 = -5e-3L;
53 
54  const std::string equation("A + B -> C + D");
55  const unsigned int n_species(4);
56 
57  int return_flag = 0;
58  std::vector<Scalar> mol_densities;
59  mol_densities.push_back(1e-2);
60  mol_densities.push_back(1e-2);
61  mol_densities.push_back(1e-2);
62  mol_densities.push_back(1e-2);
63  Scalar M = mol_densities[0];
64  for(unsigned int i = 1; i < n_species; i++)
65  {
66  M += mol_densities[i];
67  }
68 
69  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 80;
70  std::cout << type << ", tolerance = " << tol;
71  Scalar max_diff(-1.L);
72 
73  for(Scalar T = 300.1; T <= 1500.1; T += 10.)
74  {
75 
76  const Antioch::KineticsConditions<Scalar> conditions(T);
77 
78  for(unsigned int ikinmod = 0; ikinmod < 6; ikinmod++)
79  {
80 
81  Antioch::KineticsType<Scalar> *rate_kinetics1(NULL);
82  Antioch::KineticsType<Scalar> *rate_kinetics2(NULL);
83  Scalar k0,kinf,dk0_dT,dkinf_dT;
84  Scalar rate_exact;
85  Scalar derive_exact;
86  std::vector<Scalar> derive_dX_exact;
87  derive_dX_exact.resize(n_species);
89 
90  switch(ikinmod)
91  {
92  case 0:
93  {
95  rate_kinetics1 = new Antioch::HercourtEssenRate<Scalar>(Cf1,beta1,1.);
96  rate_kinetics2 = new Antioch::HercourtEssenRate<Scalar>(Cf2,beta2,1.);
97  k0 = Cf1 * pow(T,beta1); kinf = Cf2 * pow(T,beta2);
98  dk0_dT = Cf1 * pow (T,beta1) * beta1/T; dkinf_dT = Cf2 * pow (T,beta2) * beta2/T;
99  break;
100  }
101  case 1:
102  {
104  rate_kinetics1 = new Antioch::BerthelotRate<Scalar>(Cf1,D1);
105  rate_kinetics2 = new Antioch::BerthelotRate<Scalar>(Cf2,D2);
106  k0 = Cf1 * exp(D1*T); kinf = Cf2 * exp(D2*T);
107  dk0_dT = Cf1 * exp(D1*T) * D1; dkinf_dT = Cf2 * exp(D2*T) * D2;
108 
109  break;
110  }
111  case 2:
112  {
114  rate_kinetics1 = new Antioch::ArrheniusRate<Scalar>(Cf1,Ea1,1.);
115  rate_kinetics2 = new Antioch::ArrheniusRate<Scalar>(Cf2,Ea2,1.);
116  k0 = Cf1 * exp(-Ea1/T); kinf = Cf2 * exp(-Ea2/T);
117  dk0_dT = Cf1 * exp(-Ea1/T) * Ea1/pow(T,2); dkinf_dT = Cf2 * exp(-Ea2/T) * Ea2/pow(T,2);
118  break;
119  }
120  case 3:
121  {
122  kin_mod = Antioch::KineticsModel::BHE;
123  rate_kinetics1 = new Antioch::BerthelotHercourtEssenRate<Scalar>(Cf1,beta1,D1,1.);
124  rate_kinetics2 = new Antioch::BerthelotHercourtEssenRate<Scalar>(Cf2,beta2,D2,1.);
125  k0 = Cf1 * pow(T,beta1) * exp(D1 * T); kinf = Cf2 * pow(T,beta2) * exp(D2 * T);
126  dk0_dT = Cf1 * pow(T,beta1) * exp(D1 * T) * (beta1/T + D1); dkinf_dT = Cf2 * pow(T,beta2) * exp(D2 * T) * (beta2/T + D2);
127  break;
128  }
129  case 4:
130  {
132  rate_kinetics1 = new Antioch::KooijRate<Scalar>(Cf1,beta1,Ea1,1.,1.);
133  rate_kinetics2 = new Antioch::KooijRate<Scalar>(Cf2,beta2,Ea2,1.,1.);
134  k0 = Cf1 * pow(T,beta1) * exp(-Ea1/T); kinf = Cf2 * pow(T,beta2) * exp(-Ea2/T);
135  dk0_dT = Cf1 * pow(T,beta1) * exp(-Ea1/T) * (beta1/T + Ea1/pow(T,2)); dkinf_dT = Cf2 * pow(T,beta2) * exp(-Ea2/T) * (beta2/T + Ea2/pow(T,2));
136  break;
137  }
138  case 5:
139  {
141  rate_kinetics1 = new Antioch::VantHoffRate<Scalar>(Cf1,beta1,Ea1,D1,1.,1.);
142  rate_kinetics2 = new Antioch::VantHoffRate<Scalar>(Cf2,beta2,Ea2,D2,1.,1.);
143  k0 = Cf1 * pow(T,beta1) * exp(-Ea1/T + D1 * T); kinf = Cf2 * pow(T,beta2) * exp(-Ea2/T + D2 * T);
144  dk0_dT = Cf1 * pow(T,beta1) * exp(-Ea1/T + D1 * T) * (D1 + beta1/T + Ea1/pow(T,2));
145  dkinf_dT = Cf2 * pow(T,beta2) * exp(-Ea2/T + D2 * T) * (D2 + beta2/T + Ea2/pow(T,2));
146  break;
147  }
148  }
149  rate_exact = k0 / (1.L/M + k0/kinf);
150  derive_exact = rate_exact * (dk0_dT/k0 - dk0_dT/(kinf/M + k0) + k0 * dkinf_dT/(pow(kinf,2)/M + k0*kinf));
151 
152  for(unsigned int i = 0; i < n_species; i++)
153  {
154  derive_dX_exact[i] = rate_exact/(M + pow(M,2)*k0/kinf);
155  }
156 
159 
160  fall_reaction->add_forward_rate(rate_kinetics1);
161  fall_reaction->add_forward_rate(rate_kinetics2);
162  Scalar rate1 = fall_reaction->compute_forward_rate_coefficient(mol_densities,conditions);
163  Scalar rate;
164  Scalar drate_dT;
165  std::vector<Scalar> drate_dx;
166  drate_dx.resize(n_species);
167  fall_reaction->compute_forward_rate_coefficient_and_derivatives(mol_densities,conditions,rate,drate_dT,drate_dx);
168 
169  for(unsigned int i = 0; i < n_species; i++)
170  {
171  Scalar diff = abs( (drate_dx[i] - derive_dX_exact[i])/derive_dX_exact[i] );
172  if(diff > max_diff)max_diff = diff;
173  if( diff > tol )
174  {
175  std::cout << std::scientific << std::setprecision(16)
176  << "\nError: Mismatch in rate values." << std::endl
177  << "Kinetics model (see enum) " << kin_mod << std::endl
178  << "species " << i << std::endl
179  << "drate_dc(T) = " << drate_dx[i] << std::endl
180  << "drate_dc_exact = " << derive_dX_exact[i] << std::endl
181  << "relative error = " << abs((drate_dx[i] - derive_dX_exact[i])/derive_dX_exact[i]) << std::endl
182  << "tolerance = " << tol << std::endl;
183  return_flag = 1;
184  }
185  }
186  Scalar diff = abs( (rate1 - rate_exact)/rate_exact );
187  if(diff > max_diff)max_diff = diff;
188  if( diff > tol )
189  {
190  std::cout << std::scientific << std::setprecision(16)
191  << "\nError: Mismatch in rate values." << std::endl
192  << "Kinetics model (see enum) " << kin_mod << std::endl
193  << "T = " << T << " K" << std::endl
194  << "rate(T) = " << rate1 << std::endl
195  << "rate_exact = " << rate_exact << std::endl
196  << "relative error = " << abs((rate1 - rate_exact)/rate_exact) << std::endl
197  << "tolerance = " << tol << std::endl;
198 
199  return_flag = 1;
200  }
201  diff = abs( (rate - rate_exact)/rate_exact );
202  if(diff > max_diff)max_diff = diff;
203  if( diff > tol )
204  {
205  std::cout << std::scientific << std::setprecision(16)
206  << "\nError: Mismatch in rate values." << std::endl
207  << "Kinetics model (see enum) " << kin_mod << std::endl
208  << "T = " << T << " K" << std::endl
209  << "rate(T) = " << rate << std::endl
210  << "rate_exact = " << rate_exact << std::endl
211  << "relative error = " << abs((rate - rate_exact)/rate_exact) << std::endl
212  << "tolerance = " << tol << std::endl;
213 
214  return_flag = 1;
215  }
216  diff = abs( (drate_dT - derive_exact)/derive_exact );
217  if(diff > max_diff)max_diff = diff;
218  if( diff > tol )
219  {
220  std::cout << std::scientific << std::setprecision(16)
221  << "\nError: Mismatch in rate derivative values." << std::endl
222  << "Kinetics model (see enum) " << kin_mod << std::endl
223  << "T = " << T << " K" << std::endl
224  << "drate_dT(T) = " << drate_dT << std::endl
225  << "derive_exact = " << derive_exact << std::endl
226  << "relative error = " << abs((drate_dT - derive_exact)/derive_exact) << std::endl
227  << "tolerance = " << tol << std::endl;
228 
229  return_flag = 1;
230  }
231 
232  delete fall_reaction;
233  }
234  }
235 
236  std::cout << " and maximum difference = " << max_diff << std::endl;
237 
238  return return_flag;
239 }
Berthelot rate equation.
base class for kinetics models
Definition: kinetics_type.h:82
Berthelot Hercourt-Essen rate equation.
void compute_forward_rate_coefficient_and_derivatives(const VectorStateType &molar_densities, const KineticsConditions< StateType, VectorStateType > &conditions, StateType &kfwd, StateType &dkfwd_dT, VectorStateType &dkfkwd_dX) const
Antioch::enable_if_c< Antioch::is_valarray< T >::value, typename Antioch::state_type< T >::type >::type pow(const T &in, const T2 &n)
Base class for falloff processes.
void add_forward_rate(KineticsType< CoeffType, VectorCoeffType > *rate)
Add a forward rate object.
Van't Hoff rate equation.
Definition: kinetics_type.h:62
Hercourt-Essen rate equation.
StateType compute_forward_rate_coefficient(const VectorStateType &molar_densities, const KineticsConditions< StateType, VectorStateType > &conditions) const
Kooij rate equation.
Definition: kinetics_type.h:59
This class contains the conditions of the chemistry.

Generated on Thu Jul 7 2016 11:09:47 for antioch-0.4.0 by  doxygen 1.8.8