antioch-0.4.0
troe_falloff_unit.C
Go to the documentation of this file.
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 // C++
28 #include <limits>
29 // Antioch
30 #include "antioch/vector_utils.h"
31 
32 #include "antioch/reaction.h"
34 
35 template <typename Scalar>
36 int tester(const std::string & type)
37 {
38 
39  using std::abs;
40  using std::exp;
41  using std::pow;
42  using std::log;
43 
44 //values for 2 CH3 (+M) <=> C2H6 (+M) for the Kooij model, Ds are made up
45 
46  const Scalar Cf1 = 1.135e36L * 1e6L * 1e-12L; //(cm3/mol)^2/s -> kmol -> m3
47  const Scalar beta1 = 1.246L; //true value is -5.246
48  const Scalar Ea1 = 1704.8L / 1.9858775L; //cal/mol
49  const Scalar D1 = -4e-2L; // K^-1
50  const Scalar Cf2 = 6.22e16L * 1e3L * 1e-12L; //cm3/mol/s -> kmol -> m3
51  const Scalar beta2 = -1.174L;
52  const Scalar Ea2 = 635.8L / 1.9858775L; //cal/mol
53  const Scalar D2 = -5e-3L;
54  const Scalar alpha = 0.405L;
55  const Scalar T3 = 1120L; //K
56  const Scalar T1 = 69.6L; //K
57 // T2 too high, negligible
58 
59  const std::string equation("A + B -> AB");
60  const unsigned int n_species(3);
61 
62  int return_flag = 0;
63  std::vector<Scalar> mol_densities;
64  mol_densities.push_back(1e-2L);
65  mol_densities.push_back(1e-2L);
66  mol_densities.push_back(1e-2L);
67  Scalar M = mol_densities[0];
68  for(unsigned int i = 1; i < n_species; i++)
69  {
70  M += mol_densities[i];
71  }
72 
73  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 1500;
74 
75  std::cout << type << ", tolerance = " << tol;
76  Scalar max_diff(-1.L);
77  for(Scalar T = 300.1L; T <= 1500.1L; T += 10.L)
78  {
79 
80  const Antioch::KineticsConditions<Scalar> conditions(T);
81 
82  for(unsigned int ikinmod = 0; ikinmod < 6; ikinmod++)
83  {
84 
85  Antioch::KineticsType<Scalar> *rate_kinetics1(NULL);
86  Antioch::KineticsType<Scalar> *rate_kinetics2(NULL);
87  Scalar k0,kinf,dk0_dT,dkinf_dT;
88  Scalar rate_exact;
89  Scalar derive_exact;
90  std::vector<Scalar> derive_dX_exact;
91  derive_dX_exact.resize(n_species);
92 //Troe pre-calculations
93  Scalar Fcent = (1.L - alpha) * exp(-T/T3) + alpha * exp(-T/T1);
94  Scalar dFcent_dT = (alpha - 1.L)/T3 * exp(-T/T3) - alpha/T1 * exp(-T/T1);
95  Scalar dlog10Fcent_dT = Antioch::Constants::log10_to_log<Scalar>()/Fcent * dFcent_dT;
96  Scalar n = 0.75L - 1.27L * Antioch::Constants::log10_to_log<Scalar>()*log(Fcent);
97  Scalar c = - 0.40L - 0.67L * Antioch::Constants::log10_to_log<Scalar>()*log(Fcent);
98  Scalar d = 0.14L;
99  Scalar dn_dT = -1.27L * dlog10Fcent_dT;
100  Scalar dc_dT = -0.67L * dlog10Fcent_dT;
101 
103 
104  switch(ikinmod)
105  {
106  case 0:
107  {
109  rate_kinetics1 = new Antioch::HercourtEssenRate<Scalar>(Cf1,beta1,1.L);
110  rate_kinetics2 = new Antioch::HercourtEssenRate<Scalar>(Cf2,beta2,1.L);
111  k0 = Cf1 * pow(T,beta1); kinf = Cf2 * pow(T,beta2);
112  dk0_dT = Cf1 * pow (T,beta1) * beta1/T; dkinf_dT = Cf2 * pow (T,beta2) * beta2/T;
113  break;
114  }
115  case 1:
116  {
118  rate_kinetics1 = new Antioch::BerthelotRate<Scalar>(Cf1,D1);
119  rate_kinetics2 = new Antioch::BerthelotRate<Scalar>(Cf2,D2);
120  k0 = Cf1 * exp(D1*T); kinf = Cf2 * exp(D2*T);
121  dk0_dT = Cf1 * exp(D1*T) * D1; dkinf_dT = Cf2 * exp(D2*T) * D2;
122  break;
123  }
124  case 2:
125  {
127  rate_kinetics1 = new Antioch::ArrheniusRate<Scalar>(Cf1,Ea1,1.L);
128  rate_kinetics2 = new Antioch::ArrheniusRate<Scalar>(Cf2,Ea2,1.L);
129  k0 = Cf1 * exp(-Ea1/T); kinf = Cf2 * exp(-Ea2/T);
130  dk0_dT = Cf1 * exp(-Ea1/T) * Ea1/pow(T,2); dkinf_dT = Cf2 * exp(-Ea2/T) * Ea2/pow(T,2);
131  break;
132  }
133  case 3:
134  {
135  kin_mod = Antioch::KineticsModel::BHE;
136  rate_kinetics1 = new Antioch::BerthelotHercourtEssenRate<Scalar>(Cf1,beta1,D1,1.L);
137  rate_kinetics2 = new Antioch::BerthelotHercourtEssenRate<Scalar>(Cf2,beta2,D2,1.L);
138  k0 = Cf1 * pow(T,beta1) * exp(D1 * T); kinf = Cf2 * pow(T,beta2) * exp(D2 * T);
139  dk0_dT = Cf1 * pow(T,beta1) * exp(D1 * T) * (beta1/T + D1); dkinf_dT = Cf2 * pow(T,beta2) * exp(D2 * T) * (beta2/T + D2);
140  break;
141  }
142  case 4:
143  {
145  rate_kinetics1 = new Antioch::KooijRate<Scalar>(Cf1,beta1,Ea1,1.L,1.L);
146  rate_kinetics2 = new Antioch::KooijRate<Scalar>(Cf2,beta2,Ea2,1.L,1.L);
147  k0 = Cf1 * pow(T,beta1) * exp(-Ea1/T); kinf = Cf2 * pow(T,beta2) * exp(-Ea2/T);
148  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));
149  break;
150  }
151  case 5:
152  {
154  rate_kinetics1 = new Antioch::VantHoffRate<Scalar>(Cf1,beta1,Ea1,D1,1.L,1.L);
155  rate_kinetics2 = new Antioch::VantHoffRate<Scalar>(Cf2,beta2,Ea2,D2,1.L,1.L);
156  k0 = Cf1 * pow(T,beta1) * exp(-Ea1/T + D1 * T); kinf = Cf2 * pow(T,beta2) * exp(-Ea2/T + D2 * T);
157  dk0_dT = Cf1 * pow(T,beta1) * exp(-Ea1/T + D1 * T) * (D1 + beta1/T + Ea1/pow(T,2));
158  dkinf_dT = Cf2 * pow(T,beta2) * exp(-Ea2/T + D2 * T) * (D2 + beta2/T + Ea2/pow(T,2));
159  break;
160  }
161  }
162  // Troe calculations
163  Scalar Pr = M * k0/kinf;
164  Scalar dPr_dT = Pr * (dk0_dT/k0 - dkinf_dT/kinf);
165  Scalar log10Pr = Antioch::Constants::log10_to_log<Scalar>() * log(Pr);
166  Scalar dlog10Pr_dT = Antioch::Constants::log10_to_log<Scalar>() / Pr * dPr_dT;
167  std::vector<Scalar> dlog10Pr_dX(n_species,0);
168  for(unsigned int i = 0; i < n_species; i++)
169  {
170  dlog10Pr_dX[i] = Antioch::Constants::log10_to_log<Scalar>()/M;
171  }
172  Scalar logF = log(Fcent)/(1.L + pow(((log10Pr + c)/(n - d*(log10Pr + c) )),2));
173  Scalar dlogF_dT = logF * (dlog10Fcent_dT / Fcent
174  - 2.L *pow((log10Pr + c)/(n - d * (log10Pr + c)),2)
175  * ((dlog10Pr_dT + dc_dT)/(log10Pr + c) -
176  (dn_dT - d * (dlog10Pr_dT + dc_dT))/(n - d * (log10Pr + c))
177  )
178  / (1.L + pow((log10Pr + c)/(n - d * (log10Pr + c)),2))
179  );
180  Scalar F = exp(logF);
181  Scalar dF_dT = F * dlogF_dT;
182  std::vector<Scalar> dF_dX;
183  dF_dX.resize(n_species);
184  for(unsigned int i = 0; i < n_species; i++)
185  {
186  dF_dX[i] = - F * logF * logF/log(Fcent) * dlog10Pr_dX[i] * (1.L - 1.L/(n - d * (log10Pr + c))) * (log10Pr + c);
187  }
188 
189  rate_exact = k0 / (1.L/M + k0/kinf);
190 
191  derive_exact = rate_exact * (dk0_dT/k0 - dk0_dT/(kinf/M + k0) + k0 * dkinf_dT/(kinf*(kinf/M + k0))) * F
192  + dF_dT * rate_exact;
193 
194  for(unsigned int i = 0; i < n_species; i++)
195  {
196  derive_dX_exact[i] = rate_exact/(M + pow(M,2)*k0/kinf) * F + dF_dX[i] * rate_exact;
197  }
198 
199  rate_exact *= F;
200 
203  fall_reaction->add_forward_rate(rate_kinetics1);
204  fall_reaction->add_forward_rate(rate_kinetics2);
205  fall_reaction->F().set_alpha(alpha);
206  fall_reaction->F().set_T1(T1);
207  fall_reaction->F().set_T3(T3);
208  Scalar rate1 = fall_reaction->compute_forward_rate_coefficient(mol_densities,conditions);
209  Scalar rate;
210  Scalar drate_dT;
211  std::vector<Scalar> drate_dx;
212  drate_dx.resize(n_species);
213  fall_reaction->compute_forward_rate_coefficient_and_derivatives(mol_densities,conditions,rate,drate_dT,drate_dx);
214 
215  for(unsigned int i = 0; i < n_species; i++)
216  {
217  Scalar diff = abs( (drate_dx[i] - derive_dX_exact[i])/derive_dX_exact[i] );
218  if(diff > max_diff)max_diff = diff;
219 
220  if(diff > tol )
221  {
222  std::cerr << std::scientific << std::setprecision(16)
223  << "\nError: Mismatch in rate values." << std::endl
224  << "Kinetics model (see enum) " << kin_mod << std::endl
225  << "species " << i << std::endl
226  << "T = " << T << " K" << std::endl
227  << "drate_dc(T) = " << drate_dx[i] << std::endl
228  << "drate_dc_exact = " << derive_dX_exact[i] << std::endl
229  << "tolerance is " << tol << std::endl
230  << "criteria is " << abs( (drate_dx[i] - derive_dX_exact[i])/derive_dX_exact[i]) << std::endl
231  << "parameters are\nrate: " << rate_exact << "\n and " << rate_exact << std::endl
232  << "total density: " << M << " species " << mol_densities[i] << std::endl;
233  return_flag = 1;
234  }
235  }
236  Scalar diff = abs( (rate1 - rate_exact)/rate_exact );
237  if(diff > max_diff)max_diff = diff;
238  if(diff > tol )
239  {
240  std::cerr << std::scientific << std::setprecision(16)
241  << "\nError: Mismatch in rate values." << std::endl
242  << "Kinetics model (see enum) " << kin_mod << std::endl
243  << "T = " << T << " K" << std::endl
244  << "rate(T) = " << rate1 << std::endl
245  << "rate_exact = " << rate_exact << std::endl
246  << "tolerance is " << tol << std::endl
247  << "criteria is " << abs( (rate1 - rate_exact)/rate_exact ) << std::endl;
248 
249  return_flag = 1;
250  }
251  diff = abs( (rate - rate_exact)/rate_exact );
252  if(diff > max_diff)max_diff = diff;
253  if( diff > tol )
254  {
255  std::cerr << std::scientific << std::setprecision(16)
256  << "\nError: Mismatch in rate values." << std::endl
257  << "Kinetics model (see enum) " << kin_mod << std::endl
258  << "T = " << T << " K" << std::endl
259  << "rate(T) = " << rate << std::endl
260  << "rate_exact = " << rate_exact << std::endl
261  << "tolerance is " << tol << std::endl
262  << "criteria is " << abs( (rate - rate_exact)/rate_exact ) << std::endl;
263 
264  return_flag = 1;
265  }
266  diff = abs( (drate_dT - derive_exact)/derive_exact );
267  if(diff > max_diff)max_diff = diff;
268  if( diff > tol )
269  {
270  std::cerr << std::scientific << std::setprecision(16)
271  << "\nError: Mismatch in rate derivative values." << std::endl
272  << "Kinetics model (see enum) " << kin_mod << std::endl
273  << "T = " << T << " K" << std::endl
274  << "drate_dT(T) = " << drate_dT << std::endl
275  << "derive_exact = " << derive_exact << std::endl
276  << "tolerance is " << tol << std::endl
277  << "criteria is " << abs( (drate_dT - derive_exact)/derive_exact ) << std::endl;
278 
279  return_flag = 1;
280  }
281  delete fall_reaction;
282  }
283  }
284 
285  std::cout << " and maximum difference = " << max_diff << std::endl;
286 
287  return return_flag;
288 }
289 
290 int main()
291 {
292  return (tester<double>("double") ||
293  tester<long double>("long double") ||
294  tester<float>("float"));
295 }
Berthelot rate equation.
int tester(const std::string &type)
base class for kinetics models
Definition: kinetics_type.h:82
Berthelot Hercourt-Essen rate equation.
Scalar F(const Scalar &x)
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.
int main()
void add_forward_rate(KineticsType< CoeffType, VectorCoeffType > *rate)
Add a forward rate object.
const FalloffType & F() const
Return const reference to the falloff 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