antioch-0.4.0
parsing_chemkin.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 //Antioch
28 #include "antioch/vector_utils_decl.h"
29 
30 #include "antioch/physical_constants.h"
31 #include "antioch/reaction_set.h"
32 #include "antioch/chemical_species.h"
33 #include "antioch/chemical_mixture.h"
34 #include "antioch/read_reaction_set_data.h"
35 #include "antioch/units.h"
36 
37 #include "antioch/vector_utils.h"
38 
39 //C++
40 #include <cmath>
41 #include <limits>
42 
43 template<typename Scalar>
44 Scalar HE(const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Tf = 1.L)
45 {
46  return Cf * std::pow(T/Tf,eta);
47 }
48 
49 template<typename Scalar>
50 Scalar Arrh(const Scalar &T, const Scalar &Cf, const Scalar &Ea)
51 {
52  return Cf * std::exp(-Ea /(Antioch::Constants::R_universal<Scalar>() * T));
53 }
54 
55 template<typename Scalar>
56 Scalar Kooij(const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Ea, const Scalar &Tf = 1.L)
57 {
58  return Cf * std::pow(T/Tf,eta) * std::exp(-Ea /(Antioch::Constants::R_universal<Scalar>() * T));
59 }
60 
61 template<typename Scalar>
62 Scalar FcentTroe(const Scalar &T, const Scalar &alpha, const Scalar &T3, const Scalar &T1, const Scalar &T2 = -1.L)
63 {
64  return (T2 < Scalar(0.))?(Scalar(1.) - alpha) * std::exp(-T/T3) + alpha * std::exp(-T/T1):
65  (Scalar(1.) - alpha) * std::exp(-T/T3) + alpha * std::exp(-T/T1) + std::exp(-T2/T);
66 }
67 
68 template<typename Scalar>
69 Scalar coeffTroe(const Scalar &coef1, const Scalar &coef2, const Scalar &Fcent)
70 {
71  return coef1 + coef2 * std::log10(Fcent);
72 }
73 
74 template<typename Scalar>
75 Scalar cTroe(const Scalar &Fcent)
76 {
77  Scalar c1(-0.40);
78  Scalar c2(-0.67);
79  return coeffTroe(c1,c2,Fcent);
80 }
81 
82 template<typename Scalar>
83 Scalar nTroe(const Scalar &Fcent)
84 {
85  Scalar c1(0.75);
86  Scalar c2(-1.27);
87  return coeffTroe(c1,c2,Fcent);
88 }
89 
90 template<typename Scalar>
91 Scalar FTroe(const Scalar &Fcent, const Scalar &Pr)
92 {
93  Scalar d(0.14L);
94  return std::pow(10.L,std::log10(Fcent)/(Scalar(1.) +
95  std::pow( (std::log10(Pr) + cTroe(Fcent)) /
96  ( nTroe(Fcent) - d * (std::log10(Pr) + cTroe(Fcent)) ),2)));
97 }
98 
99 template<typename VectorScalar>
100 typename Antioch::value_type<VectorScalar>::type k_photo(const VectorScalar &solar_lambda, const VectorScalar &solar_irr,
101  const VectorScalar &sigma_lambda, const VectorScalar &sigma_sigma)
102 {
103  Antioch::SigmaBinConverter<VectorScalar> bin;
104  VectorScalar sigma_rescaled;
105  bin.y_on_custom_grid(sigma_lambda,sigma_sigma,solar_lambda,sigma_rescaled);
106 
107  typename Antioch::value_type<VectorScalar>::type _k(0.L);
108  for(unsigned int il = 0; il < solar_irr.size() - 1; il++)
109  {
110  _k += sigma_rescaled[il] * solar_irr[il] * (solar_lambda[il+1] - solar_lambda[il]);
111  }
112 
113  return _k;
114 }
115 
116 template<typename Scalar>
117 int tester(const std::string &root_name)
118 {
119 
120  std::vector<std::string> species_str_list;
121  species_str_list.push_back( "O2" );
122  species_str_list.push_back( "OH" );
123  species_str_list.push_back( "H2" );
124  species_str_list.push_back( "H2O" );
125  species_str_list.push_back( "H2O2" );
126  species_str_list.push_back( "HO2" );
127  species_str_list.push_back( "O" );
128  species_str_list.push_back( "CH3" );
129  species_str_list.push_back( "CH4" );
130  species_str_list.push_back( "H" );
131  unsigned int n_species = species_str_list.size();
132 
133  Antioch::ChemicalMixture<Scalar> chem_mixture( species_str_list );
134  Antioch::ReactionSet<Scalar> reaction_set( chem_mixture );
135  Antioch::read_reaction_set_data_chemkin<Scalar>( root_name + "/test_parsing.chemkin", true, reaction_set );
136 
137  Scalar T = 2000.L;
138  Antioch::Units<Scalar> unitA_m1("(mol/cm3)/s"),unitA_0("s-1"),unitA_1("cm3/mol/s"),unitA_2("cm6/mol2/s"),
139  unitEa_cal("cal/mol");
140 
141 //
142  // Molar densities
143  std::vector<Scalar> molar_densities(n_species,5e-4);
144  Scalar tot_dens((Scalar)n_species * 5e-4);
145 
147  std::vector<Scalar> k;
148  Scalar A,b,Ea;
149 /*
150 ! Hessler, J. Phys. Chem. A, 102:4517 (1998)
151 H+O2=O+OH 3.547e+15 -0.406 1.6599E+4
152 */
153  A = 3.547e15L * unitA_1.get_SI_factor();
154  b = -0.406L;
155  Ea = 1.6599e4L * unitEa_cal.get_SI_factor();
156  k.push_back(Kooij(T,A,b,Ea));
157 
158 
159 /*
160 ! Sutherland et al., 21st Symposium, p. 929 (1986)
161 O+H2=H+OH 0.508E+05 2.67 0.629E+04
162 */
163  A = 0.508e5L * unitA_1.get_SI_factor();
164  b = 2.67L;
165  Ea = 0.629e4L * unitEa_cal.get_SI_factor();
166  k.push_back(Kooij(T,A,b,Ea));
167 
168 /*
169 ! Michael and Sutherland, J. Phys. Chem. 92:3853 (1988)
170 H2+OH=H2O+H 0.216E+09 1.51 0.343E+04
171 */
172  A = 0.216e9L * unitA_1.get_SI_factor();
173  b = 1.51L;
174  Ea = 0.343e4L * unitEa_cal.get_SI_factor();
175  k.push_back(Kooij(T,A,b,Ea));
176 
177 /*
178 ! Sutherland et al., 23rd Symposium, p. 51 (1990)
179 O+H2O=OH+OH 2.97e+06 2.02 1.34e+4
180 */
181  A = 2.97e6L * unitA_1.get_SI_factor();
182  b = 2.02L;
183  Ea = 1.34e4L * unitEa_cal.get_SI_factor();
184  k.push_back(Kooij(T,A,b,Ea));
185 
187 /*
188 ! Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986)
189 H2+M=H+H+M 4.577E+19 -1.40 1.0438E+05
190  H2/2.5/ H2O/12/
191 */
192  A = 4.577e19L * unitA_1.get_SI_factor();
193  b = -1.40L;
194  Ea = 1.0438e5L * unitEa_cal.get_SI_factor();
195  Scalar sum_eps = 5e-4L * (2.5L + 12.L + (Scalar)(species_str_list.size() - 2));
196  k.push_back(sum_eps * Kooij(T,A,b,Ea));
197 
198 /*
199 ! Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986)
200 O+O+M=O2+M 6.165E+15 -0.50 0.000E+00
201  H2/2.5/ H2O/12/
202 */
203  A = 6.165e15L * unitA_2.get_SI_factor();
204  b = -0.50L;
205  k.push_back(sum_eps * HE(T,A,b));
206 
207 /*
208 ! Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986)
209 O+H+M=OH+M 4.714E+18 -1.00 0.000E+00
210  H2/2.5/ H2O/12/
211 */
212  A = 4.714e18L * unitA_2.get_SI_factor();
213  b = -1.00L;
214  k.push_back(sum_eps * HE(T,A,b));
215 
216 /*
217 ! Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986)
218 !H+OH+M=H2O+M 2.212E+22 -2.00 0.000E+00
219 H+OH+M=H2O+M 3.800E+22 -2.00 0.000E+00
220  H2/2.5/ H2O/12/
221 */
222  A = 3.800e22L * unitA_2.get_SI_factor();
223  b = -2.00L;
224  k.push_back(sum_eps * HE(T,A,b));
225 
226 /*
227 !************** Formation and Consumption of HO2******************
228 
229 ! Cobos et al., J. Phys. Chem. 89:342 (1985) for kinf
230 ! Michael, et al., J. Phys. Chem. A, 106:5297 (2002) for k0
231 
232 !******************************************************************************
233 */
234 /*
235 ! MAIN BATH GAS IS N2 (comment this reaction otherwise)
236 !
237  H+O2(+M)=HO2(+M) 1.475E+12 0.60 0.00E+00
238  LOW/6.366E+20 -1.72 5.248E+02/
239  TROE/0.8 1E-30 1E+30/
240  H2/2.0/ H2O/11./ O2/0.78/
241 */
242  A = 6.366e20L * unitA_2.get_SI_factor();
243  b = -1.72L;
244  Ea = 5.248e2L * unitEa_cal.get_SI_factor();
245  Scalar k0 = Kooij(T,A,b,Ea);
246  A = 1.475e12L * unitA_1.get_SI_factor();
247  b = 0.60L;
248  Ea = 0.00L * unitEa_cal.get_SI_factor();
249  Scalar M = tot_dens + Scalar(5.39e-3L);
250  Scalar kinf = Kooij(T,A,b,Ea);
251  Scalar Pr = M * k0/kinf;
252  Scalar Fc = FcentTroe(T,(Scalar)0.8L,(Scalar)1e-30L,(Scalar)1e30L);
253  k.push_back(k0 / (1.L/M + k0/kinf) * FTroe(Fc,Pr));
254 
255 /*
256 ! Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) [modified]
257 HO2+H=H2+O2 1.66E+13 0.00 0.823E+03
258 */
259  A = 1.66e13L * unitA_1.get_SI_factor();
260  Ea = 0.823e3L * unitEa_cal.get_SI_factor();
261  k.push_back(Arrh(T,A,Ea));
262 
263 /*
264 ! Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) [modified]
265 HO2+H=OH+OH 7.079E+13 0.00 2.95E+02
266 */
267  A = 7.079e13L * unitA_1.get_SI_factor();
268  Ea = 2.95e2L * unitEa_cal.get_SI_factor();
269  k.push_back(Arrh(T,A,Ea));
270 
271 /*
272 ! Baulch et al., J. Phys. Chem. Ref Data, 21:411 (1992)
273 HO2+O=O2+OH 0.325E+14 0.00 0.00E+00
274 */
275  A = 0.325e14L * unitA_1.get_SI_factor();
276  k.push_back(A);
277 
278 /*
279 ! Keyser, J. Phys. Chem. 92:1193 (1988)
280 HO2+OH=H2O+O2 2.890E+13 0.00 -4.970E+02
281 */
282  A = 2.890e13L * unitA_1.get_SI_factor();
283  Ea = -4.97e2L * unitEa_cal.get_SI_factor();
284  k.push_back(Arrh(T,A,Ea));
285 
287 /*
288 ! Hippler et al., J. Chem. Phys. 93:1755 (1990)
289 HO2+HO2=H2O2+O2 4.200e+14 0.00 1.1982e+04
290  DUPLICATE
291 HO2+HO2=H2O2+O2 1.300e+11 0.00 -1.6293e+3
292  DUPLICATE
293 */
294  A = 4.200e14L * unitA_1.get_SI_factor();
295  Ea = 1.1982e4L * unitEa_cal.get_SI_factor();
296  Scalar A2 = 1.300e11L * unitA_1.get_SI_factor();
297  Scalar Ea2 = -1.6293e3L * unitEa_cal.get_SI_factor();
298  k.push_back(Arrh(T,A,Ea) + Arrh(T,A2,Ea2));
299 
300 /*
301 ! Brouwer et al., J. Chem. Phys. 86:6171 (1987) for kinf
302 ! Warnatz, J. in Combustion chemistry (1984) for k0
303 H2O2(+M)=OH+OH(+M) 2.951e+14 0.00 4.843E+04
304  LOW/1.202E+17 0.00 4.55E+04/
305  TROE/0.5 1E-30 1E+30/
306  H2/2.5/ H2O/12/
307 */
308  A = 1.202e17L * unitA_1.get_SI_factor();
309  Ea = 4.55e4L * unitEa_cal.get_SI_factor();
310  k0 = Arrh(T,A,Ea);
311  A = 2.951e14L * unitA_0.get_SI_factor();
312  Ea = 4.843e4L * unitEa_cal.get_SI_factor();
313  M = tot_dens + 6.25e-3;
314  kinf = Arrh(T,A,Ea);
315  Pr = M * k0/kinf;
316  Fc = FcentTroe(T,(Scalar)0.5L,(Scalar)1e-30L,(Scalar)1e30L);
317  k.push_back(k0 / (1.L/M + k0/kinf) * FTroe(Fc,Pr));
318 //
319 
320 /*
321 ! Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986)
322 H2O2+H=H2O+OH 0.241E+14 0.00 0.397E+04
323 */
324  A = 0.241e14L * unitA_1.get_SI_factor();
325  Ea = 0.397e4L * unitEa_cal.get_SI_factor();
326  k.push_back(Arrh(T,A,Ea));
327 
328 /*
329 ! Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986)
330 H2O2+H=HO2+H2 0.482E+14 0.00 0.795E+04
331 */
332  A = 0.482e14L * unitA_1.get_SI_factor();
333  Ea = 0.795e4L * unitEa_cal.get_SI_factor();
334  k.push_back(Arrh(T,A,Ea));
335 
336 /*
337 ! Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986)
338 H2O2+O=OH+HO2 9.550E+06 2.00 3.970E+03
339 */
340  A = 9.550e6L * unitA_1.get_SI_factor();
341  b = 2.00;
342  Ea = 3.970e3L * unitEa_cal.get_SI_factor();
343  k.push_back(Kooij(T,A,b,Ea));
344 
345 /*
346 ! Hippler and Troe, J. Chem. Phys. Lett. 192:333 (1992)
347 H2O2+OH=HO2+H2O 1.000E+12 0.00 0.000
348  DUPLICATE
349 H2O2+OH=HO2+H2O 5.800E+14 0.00 9.557E+03
350  DUPLICATE
351 */
352  A = 1.000e12L * unitA_1.get_SI_factor();
353  Ea = 0.000L * unitEa_cal.get_SI_factor();
354  A2 = 5.800e14L * unitA_1.get_SI_factor();
355  Ea2 = 9.557e3L * unitEa_cal.get_SI_factor();
356  k.push_back(Arrh(T,A,Ea) + Arrh(T,A2,Ea2));
357 
362  A = 1e12 * unitA_1.get_SI_factor();
363  b = 1.2;
364  Ea = 3.125e4 * unitEa_cal.get_SI_factor();
365  k.push_back(Kooij(T,A,b,Ea));
366  A = 1e8 * unitA_0.get_SI_factor();
367  b = 0.8;
368  Ea = 3.5e3 * unitEa_cal.get_SI_factor();
369  k.push_back(Kooij(T,A,b,Ea));
370 
371 
372  const Scalar tol = (std::numeric_limits<Scalar>::epsilon() < 1e-17L)?7e-16L:
373  std::numeric_limits<Scalar>::epsilon() * 100;
374  int return_flag(0);
375 
376  if(reaction_set.n_reactions() != k.size())
377  {
378  std::cerr << reaction_set << std::endl;
379  std::cerr << "Not the right number of reactions" << std::endl;
380  std::cerr << reaction_set.n_reactions() << " instead of " << k.size() << std::endl;
381  return_flag = 1;
382  }
383  {
384  for(unsigned int ir = 0; ir < k.size(); ir++)
385  {
386  const Antioch::Reaction<Scalar> * reac = &reaction_set.reaction(ir);
387 
388  if(std::abs(k[ir] - reac->compute_forward_rate_coefficient(molar_densities,T))/k[ir] > tol)
389  {
390  std::cout << *reac << std::endl;
391  std::cout << std::scientific << std::setprecision(16)
392  << "Error in kinetics comparison\n"
393  << "reaction #" << ir << "\n"
394  << "temperature: " << T << " K" << "\n"
395  << "theory: " << k[ir] << "\n"
396  << "calculated: " << reac->compute_forward_rate_coefficient(molar_densities,T) << "\n"
397  << "relative error = " << std::abs(k[ir] - reac->compute_forward_rate_coefficient(molar_densities,T))/k[ir] << "\n"
398  << "tolerance = " << tol
399  << std::endl;
400  return_flag = 1;
401  }
402  }
403  }
404 
405  return return_flag;
406 }
407 
408 int main(int argc, char* argv[])
409 {
410  return (tester<float>(std::string(argv[1])) ||
411  tester<double>(std::string(argv[1])) ||
412  tester<long double>(std::string(argv[1])));
413 }
FTroe
Scalar FTroe(const Scalar &Fcent, const Scalar &Pr)
Definition: parsing_chemkin.C:91
Antioch::Reaction::compute_forward_rate_coefficient
StateType compute_forward_rate_coefficient(const VectorStateType &molar_densities, const KineticsConditions< StateType, VectorStateType > &conditions) const
Definition: reaction.h:866
Antioch::Reaction
A single reaction mechanism.
Definition: reaction.h:108
tester
int tester(const std::string &root_name)
Definition: parsing_chemkin.C:117
Antioch::ReactionSet
This class encapsulates all the reaction mechanisms considered in a chemical nonequilibrium simulatio...
Definition: kinetics_evaluator.h:42
Antioch::ReactionSet::reaction
const Reaction< CoeffType > & reaction(const unsigned int r) const
Definition: reaction_set.h:232
Kooij
Scalar Kooij(const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Ea, const Scalar &Tf=1.L)
Definition: parsing_chemkin.C:56
main
int main(int argc, char *argv[])
Definition: parsing_chemkin.C:408
Arrh
Scalar Arrh(const Scalar &T, const Scalar &Cf, const Scalar &Ea)
Definition: parsing_chemkin.C:50
HE
Scalar HE(const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Tf=1.L)
Definition: parsing_chemkin.C:44
std::pow
Antioch::enable_if_c< Antioch::is_valarray< T >::value, typename Antioch::state_type< T >::type >::type pow(const T &in, const T2 &n)
Definition: valarray_utils.h:71
coeffTroe
Scalar coeffTroe(const Scalar &coef1, const Scalar &coef2, const Scalar &Fcent)
Definition: parsing_chemkin.C:69
Antioch::Units::get_SI_factor
T get_SI_factor() const
Multiplicative coefficient getter.
Definition: units.h:334
nTroe
Scalar nTroe(const Scalar &Fcent)
Definition: parsing_chemkin.C:83
Antioch::ReactionSet::n_reactions
unsigned int n_reactions() const
Definition: reaction_set.h:200
k_photo
Antioch::value_type< VectorScalar >::type k_photo(const VectorScalar &solar_lambda, const VectorScalar &solar_irr, const VectorScalar &sigma_lambda, const VectorScalar &sigma_sigma)
Definition: parsing_chemkin.C:100
cTroe
Scalar cTroe(const Scalar &Fcent)
Definition: parsing_chemkin.C:75
units.h
Advanced unit class.
FcentTroe
Scalar FcentTroe(const Scalar &T, const Scalar &alpha, const Scalar &T3, const Scalar &T1, const Scalar &T2=-1.L)
Definition: parsing_chemkin.C:62
Antioch::SigmaBinConverter::y_on_custom_grid
void y_on_custom_grid(const VectorCoeffType &x_old, const VectorCoeffType &y_old, const VectorStateType &x_new, VectorStateType &y_new) const
Definition: sigma_bin_converter.h:79
Antioch::ChemicalMixture
Class storing chemical mixture properties.
Definition: chemical_mixture.h:75
Antioch::Units
An advanced unit class.
Definition: units.h:111
read_reaction_set_data.h
We parse the file here, with an exhaustive unit management.
Generated on Thu Jul 7 2016 11:09:47 for antioch-0.4.0 by  doxygen 1.8.8