antioch-0.4.0
Functions
parsing_xml.C File Reference
#include "antioch/vector_utils_decl.h"
#include "antioch/physical_constants.h"
#include "antioch/reaction_set.h"
#include "antioch/chemical_species.h"
#include "antioch/chemical_mixture.h"
#include "antioch/read_reaction_set_data.h"
#include "antioch/units.h"
#include "antioch/vector_utils.h"
#include <cmath>
#include <limits>

Go to the source code of this file.

Functions

template<typename Scalar >
Scalar HE (const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Tf=1.)
 
template<typename Scalar >
Scalar Bert (const Scalar &T, const Scalar &Cf, const Scalar &D)
 
template<typename Scalar >
Scalar BHE (const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &D, const Scalar &Tf=1.)
 
template<typename Scalar >
Scalar Arrh (const Scalar &T, const Scalar &Cf, const Scalar &Ea, const Scalar &R=Antioch::Constants::R_universal< Scalar >())
 
template<typename Scalar >
Scalar Kooij (const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Ea, const Scalar &Tf=1., const Scalar &R=Antioch::Constants::R_universal< Scalar >())
 
template<typename Scalar >
Scalar VH (const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Ea, const Scalar &D, const Scalar &Tf=1., const Scalar &R=Antioch::Constants::R_universal< Scalar >())
 
template<typename Scalar >
Scalar FcentTroe (const Scalar &T, const Scalar &alpha, const Scalar &T3, const Scalar &T1, const Scalar &T2=-1.)
 
template<typename Scalar >
Scalar coeffTroe (const Scalar &coef1, const Scalar &coef2, const Scalar &Fcent)
 
template<typename Scalar >
Scalar cTroe (const Scalar &Fcent)
 
template<typename Scalar >
Scalar nTroe (const Scalar &Fcent)
 
template<typename Scalar >
Scalar FTroe (const Scalar &Fcent, const Scalar &Pr)
 
template<typename VectorScalar >
Antioch::value_type
< VectorScalar >::type 		k_photo (const VectorScalar &solar_lambda, const VectorScalar &solar_irr, const VectorScalar &sigma_lambda, const VectorScalar &sigma_sigma)
 
template<typename Scalar >
int tester (const std::string &kin_file, const std::string &solar_file, const std::string &CH4_cs_file)
 
int main (int argc, char *argv[])
 

Function Documentation

template<typename Scalar >
Scalar Arrh ( const Scalar &  T,
const Scalar &  Cf,
const Scalar &  Ea,
const Scalar &  R = Antioch::Constants::R_universal<Scalar>() 
)

Definition at line 61 of file parsing_xml.C.

Referenced by tester().

62 {
63  return Cf * std::exp(-Ea /(R * T));
64 }
template<typename Scalar >
Scalar Bert ( const Scalar &  T,
const Scalar &  Cf,
const Scalar &  D 
)

Definition at line 49 of file parsing_xml.C.

Referenced by tester().

50 {
51  return Cf * std::exp(D * T);
52 }
template<typename Scalar >
Scalar BHE ( const Scalar &  T,
const Scalar &  Cf,
const Scalar &  eta,
const Scalar &  D,
const Scalar &  Tf = 1. 
)

Definition at line 55 of file parsing_xml.C.

References std::pow().

Referenced by tester().

56 {
57  return Cf * std::pow(T/Tf,eta) * std::exp(D * T);
58 }
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
template<typename Scalar >
Scalar coeffTroe ( const Scalar &  coef1,
const Scalar &  coef2,
const Scalar &  Fcent 
)

Definition at line 88 of file parsing_xml.C.

Referenced by cTroe(), and nTroe().

89 {
90  return coef1 + coef2 * std::log10(Fcent);
91 }
template<typename Scalar >
Scalar cTroe ( const Scalar &  Fcent)

Definition at line 94 of file parsing_xml.C.

References coeffTroe().

Referenced by FTroe().

95 {
96  Scalar c1(-0.40);
97  Scalar c2(-0.67);
98  return coeffTroe(c1,c2,Fcent);
99 }
coeffTroe
Scalar coeffTroe(const Scalar &coef1, const Scalar &coef2, const Scalar &Fcent)
Definition: parsing_xml.C:88
template<typename Scalar >
Scalar FcentTroe ( const Scalar &  T,
const Scalar &  alpha,
const Scalar &  T3,
const Scalar &  T1,
const Scalar &  T2 = -1. 
)

Definition at line 81 of file parsing_xml.C.

Referenced by tester().

82 {
83  return (T2 < Scalar(0.))?(Scalar(1.) - alpha) * std::exp(-T/T3) + alpha * std::exp(-T/T1):
84  (Scalar(1.) - alpha) * std::exp(-T/T3) + alpha * std::exp(-T/T1) + std::exp(-T2/T);
85 }
template<typename Scalar >
Scalar FTroe ( const Scalar &  Fcent,
const Scalar &  Pr 
)

Definition at line 110 of file parsing_xml.C.

References cTroe(), nTroe(), and std::pow().

Referenced by tester().

111 {
112  Scalar d(0.14L);
113  return std::pow(10.L,std::log10(Fcent)/(Scalar(1.) +
114  std::pow( (std::log10(Pr) + cTroe(Fcent)) /
115  ( nTroe(Fcent) - d * (std::log10(Pr) + cTroe(Fcent)) ),2)));
116 }
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
nTroe
Scalar nTroe(const Scalar &Fcent)
Definition: parsing_xml.C:102
cTroe
Scalar cTroe(const Scalar &Fcent)
Definition: parsing_xml.C:94
template<typename Scalar >
Scalar HE ( const Scalar &  T,
const Scalar &  Cf,
const Scalar &  eta,
const Scalar &  Tf = 1. 
)

Definition at line 44 of file parsing_xml.C.

References std::pow().

Referenced by tester().

45 {
46  return Cf * std::pow(T/Tf,eta);
47 }
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
template<typename VectorScalar >
Antioch::value_type<VectorScalar>::type k_photo ( const VectorScalar &  solar_lambda,
const VectorScalar &  solar_irr,
const VectorScalar &  sigma_lambda,
const VectorScalar &  sigma_sigma 
)

Definition at line 119 of file parsing_xml.C.

References Antioch::SigmaBinConverter< VectorCoeffType >::y_on_custom_grid().

Referenced by tester().

121 {
122  Antioch::SigmaBinConverter<VectorScalar> bin;
123  VectorScalar sigma_rescaled(solar_lambda.size());
124  bin.y_on_custom_grid(sigma_lambda,sigma_sigma,solar_lambda,sigma_rescaled);
125 
126  typename Antioch::value_type<VectorScalar>::type _k(0.L);
127  for(unsigned int il = 0; il < solar_irr.size() - 1; il++)
128  {
129  _k += sigma_rescaled[il] * solar_irr[il] * (solar_lambda[il+1] - solar_lambda[il]);
130  }
131 
132  return _k;
133 }
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
template<typename Scalar >
Scalar Kooij ( const Scalar &  T,
const Scalar &  Cf,
const Scalar &  eta,
const Scalar &  Ea,
const Scalar &  Tf = 1.,
const Scalar &  R = Antioch::Constants::R_universal<Scalar>() 
)

Definition at line 67 of file parsing_xml.C.

References std::pow().

Referenced by tester().

69 {
70  return Cf * std::pow(T/Tf,eta) * std::exp(-Ea /(R * T));
71 }
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
int main ( int  argc,
char *  argv[] 
)

Definition at line 563 of file parsing_xml.C.

564 {
565  return (tester<float>(std::string(argv[1]),std::string(argv[2]),std::string(argv[3])) ||
566  tester<double>(std::string(argv[1]),std::string(argv[2]),std::string(argv[3])) ||
567  tester<long double>(std::string(argv[1]),std::string(argv[2]),std::string(argv[3])));
568 }
template<typename Scalar >
Scalar nTroe ( const Scalar &  Fcent)

Definition at line 102 of file parsing_xml.C.

References coeffTroe().

Referenced by FTroe().

103 {
104  Scalar c1(0.75);
105  Scalar c2(-1.27);
106  return coeffTroe(c1,c2,Fcent);
107 }
coeffTroe
Scalar coeffTroe(const Scalar &coef1, const Scalar &coef2, const Scalar &Fcent)
Definition: parsing_xml.C:88
template<typename Scalar >
int tester ( const std::string &  kin_file,
const std::string &  solar_file,
const std::string &  CH4_cs_file 
)

Elementary, + Kooij - Arrhenius conversion tested

Duplicate

Lindemann Falloff

Definition at line 136 of file parsing_xml.C.

References Antioch::KineticsModel::A, Antioch::KineticsConditions< StateType, VectorStateType >::add_particle_flux(), Arrh(), Bert(), BHE(), Antioch::Reaction< CoeffType, VectorCoeffType >::compute_forward_rate_coefficient(), Antioch::KineticsModel::D, Antioch::Units< T >::factor_to_some_unit(), FcentTroe(), FTroe(), Antioch::Units< T >::get_SI_factor(), HE(), k_photo(), Kooij(), Antioch::ReactionSet< CoeffType >::reaction(), and VH().

137 {
138 
139  std::vector<std::string> species_str_list;
140  species_str_list.push_back( "N2" );
141  species_str_list.push_back( "O2" );
142  species_str_list.push_back( "N" );
143  species_str_list.push_back( "O" );
144  species_str_list.push_back( "NO" );
145  species_str_list.push_back( "C" );
146  species_str_list.push_back( "C2" );
147  species_str_list.push_back( "CN" );
148  species_str_list.push_back( "CH4" );
149  species_str_list.push_back( "CH3" );
150  species_str_list.push_back( "H" );
151  unsigned int n_species = species_str_list.size();
152 
153  Antioch::ChemicalMixture<Scalar> chem_mixture( species_str_list );
154  Antioch::ReactionSet<Scalar> reaction_set( chem_mixture );
155  Antioch::read_reaction_set_data_xml<Scalar>( kin_file, true, reaction_set );
156 
157 //photochemistry set here
158  std::vector<Scalar> hv,lambda;
159  std::ifstream solar_flux(solar_file);
160  std::string line;
161 
162 
164 // all the steps, if ever someone needs a reference
165  getline(solar_flux,line);
166  Antioch::Units<Scalar> solar_wave("nm");
167  Antioch::Units<Scalar> solar_irra("W/m2/nm");
168  Antioch::Units<Scalar> i_unit = solar_irra - (Antioch::Constants::Planck_constant_unit<Scalar>() + Antioch::Constants::light_celerity_unit<Scalar>() - solar_wave); //photons.s-1 = irradiance/(h*c/lambda)
169  i_unit += Antioch::Units<Scalar>("nm"); //supress bin in unit calculations
170 
171  while(!solar_flux.eof())
172  {
173  Scalar l,i,di;
174  solar_flux >> l >> i >> di;
175 
176  hv.push_back(i /(Antioch::Constants::Planck_constant<Scalar>() * Antioch::Constants::light_celerity<Scalar>() / l) // irr/(h*c/lambda): power -> number of photons.s-1
177  * i_unit.get_SI_factor()); //SI for cs, keep nm for bin
178  lambda.push_back(l * solar_wave.factor_to_some_unit("nm")); //nm
179  if(lambda.size() == 796)break;
180  }
181  solar_flux.close();
182 
183  std::vector<Scalar> CH4_s,CH4_lambda;
184  std::ifstream CH4_file(CH4_cs_file);
185 
186  Scalar T = 2000.L;
187  Scalar Tr = 1.;
188  Antioch::Units<Scalar> unitA_m1("kmol/m3/s"),unitA_0("s-1"),unitA_1("m3/kmol/s"),unitA_2("m6/kmol2/s");
189 
190  Scalar Rcal = Antioch::Constants::R_universal<Scalar>() * Antioch::Constants::R_universal_unit<Scalar>().factor_to_some_unit("cal/mol/K");
191  getline(CH4_file,line);
192 
193  Antioch::Units<Scalar> cs_input("cm2");
194  Antioch::Units<Scalar> lambda_input("ang");
195  Scalar factor_cs = cs_input.get_SI_factor() / lambda_input.factor_to_some_unit("nm");
196  while(!CH4_file.eof())
197  {
198  Scalar l,s;
199  CH4_file >> l >> s;
200  CH4_s.push_back(s * factor_cs);
201  CH4_lambda.push_back(l * lambda_input.factor_to_some_unit("nm"));
202  if(CH4_s.size() == 137)break;
203  }
204  CH4_file.close();
205 
206  Antioch::ParticleFlux<std::vector<Scalar> > photons(lambda,hv);
207 
208  Antioch::KineticsConditions<Scalar,std::vector<Scalar> > conditions(T);
209 
210 //
211  // Molar densities
212  std::vector<Scalar> molar_densities(n_species,5e-4);
213  Scalar tot_dens((Scalar)n_species * 5e-4);
214 
216  std::vector<Scalar> k;
217  Scalar A,beta,Ea,D;
218 
219 // N2 -> 2 N
220  A = 7e18 * unitA_0.get_SI_factor();
221  beta = -1.6;
222  k.push_back(HE(T,A,beta));
223 
224 // O2 -> 2 O
225  A = 2e18 * unitA_0.get_SI_factor();
226  D = -5e-3;
227  k.push_back(Bert(T,A,D));
228 
229 //NO -> N + O
230  A = 5e12 * unitA_0.get_SI_factor();
231  Ea = 149943.0;
232  k.push_back(Arrh(T,A,Ea,Rcal));
233  beta = 0.42;
234  k.push_back(Kooij(T,A,beta,Ea,Tr,Rcal));
235 
236 //N2 + O -> NO + N
237  A = 5.7e9 * unitA_1.get_SI_factor();
238  beta = 0.42;
239  k.push_back(BHE(T,A,beta,D));
240 
241 //NO + O -> NO + N
242  A = 8.4e9 * unitA_1.get_SI_factor();
243  beta = 0.4;
244  Ea = 38526.0;
245  k.push_back(Kooij(T,A,beta,Ea,Tr,Rcal));
246  k.push_back(Arrh(T,A,Ea,Rcal));
247 
248 //C2 -> 2 C
249  A = 3.7e11 * unitA_0.get_SI_factor();
250  beta = -0.42;
251 // Ea = 138812.8; // cal/mol
252  Ea = 69900; // K
253  k.push_back(Kooij(T,A,beta,Ea,Tr,Scalar(1)));
254 
255 //CN -> C + N
256  A = 2.5e11 * unitA_0.get_SI_factor();
257  beta = 0.40;
258  Ea = 174240.9;
259  D = 0.05;
260  k.push_back(VH(T,A,beta,Ea,D,Tr,Rcal));
261 
263  Scalar A2,beta2,Ea2,D2,A3,beta3,Ea3,D3,A4,Ea4;
264 
265 // N2 -> 2 N
266  A = 7e18 * unitA_0.get_SI_factor();
267  beta = -1.6;
268  A2 = 5e17 * unitA_0.get_SI_factor();
269  beta2 = 0.5;
270  A3 = 3e18 * unitA_0.get_SI_factor();
271  beta3 = -0.6;
272  k.push_back(HE(T,A,beta) + HE(T,A2,beta2) + HE(T,A3,beta3));
273 
274 // O2 -> 2 O
275  A = 2e18 * unitA_0.get_SI_factor();
276  D = -5e-2;
277  A2 = 2e+16 * unitA_0.get_SI_factor();
278  D2 = 0.003;
279  k.push_back(Bert(T,A,D) + Bert(T,A2,D2));
280 
281 // NO -> N + O
282  A = 5e+12 * unitA_0.get_SI_factor();
283  Ea = 149943.0;
284  A2 = 3.5e+10 * unitA_0.get_SI_factor();
285  Ea2 = 1943.0;
286  A3 = 1.5e+8 * unitA_0.get_SI_factor();
287  Ea3 = 149.0;
288  A4 = 5.5e+8 * unitA_0.get_SI_factor();
289  Ea4 = 943.0;
290  k.push_back(Arrh(T,A,Ea,Rcal) + Arrh(T,A2,Ea2,Rcal) + Arrh(T,A3,Ea3,Rcal) + Arrh(T,A4,Ea4,Rcal));
291 
292 // N2 + O -> NO + N
293  A = 5.7e+9 * unitA_1.get_SI_factor();
294  beta = 0.42;
295  D = -5e-3;
296  A2 = 7e+7 * unitA_1.get_SI_factor();
297  beta2 = 0.5;
298  D2 = 2.5e-5;
299  k.push_back(BHE(T,A,beta,D) + BHE(T,A2,beta2,D2));
300 
301 //NO + O -> NO + N
302  A = 8.4e+09 * unitA_1.get_SI_factor();
303  beta = 0.40;
304  Ea = 38526.0;
305  A2 = 4e+07 * unitA_1.get_SI_factor();
306  beta2 = 0.50;
307  Ea2 = 40500.0;
308  A3 = 5e+10 * unitA_1.get_SI_factor();
309  beta3 = 0.10;
310  Ea3 = 15000.0;
311  k.push_back(Kooij(T,A,beta,Ea,Tr,Rcal) + Kooij(T,A2,beta2,Ea2,Tr,Rcal) + Kooij(T,A3,beta3,Ea3,Tr,Rcal));
312 
313 //C2 -> 2 C
314  A = 3.7e+11 * unitA_0.get_SI_factor();
315  beta = -0.42;
316  Ea = 138812.8;
317  A2 = 5.0e+10 * unitA_0.get_SI_factor();
318  beta2 = 1.32;
319  Ea2 = 150500.8;
320  k.push_back(Kooij(T,A,beta,Ea,Tr,Rcal) + Kooij(T,A2,beta2,Ea2,Tr,Rcal));
321 
322 //CN -> C + N
323  A = 2.5e+11 * unitA_0.get_SI_factor();
324  beta = 0.40;
325  D = -5e-3;
326  Ea = 174240.9;
327  A2 = 5e+10 * unitA_0.get_SI_factor();
328  beta2 = 0.50;
329  D2 = -1.5e-2;
330  Ea2 = 4240.9;
331  A3 = 3.2e+10 * unitA_0.get_SI_factor();
332  beta3 = 1.20;
333  D3 = -2.5e-5;
334  Ea3 = 174.9;
335  k.push_back(VH(T,A,beta,Ea,D,Tr,Rcal) + VH(T,A2,beta2,Ea2,D2,Tr,Rcal) + VH(T,A3,beta3,Ea3,D3,Tr,Rcal));
336 
337 //three body
338 // N2 -> 2 N
339  A = 7e18 * unitA_1.get_SI_factor();
340  beta = -1.6;
341  Ea = 149943.0;
342  k.push_back(HE(T,A,beta) * (Scalar(n_species) - 2. + 4.2857 + 4.2857) * 5e-4);
343 
344 // O2 -> 2 O
345  A = 2e18 * unitA_1.get_SI_factor();
346  D = -5e-3;
347  k.push_back(Bert(T,A,D) * (Scalar(n_species) - 2. + 5.0 + 5.0) * 5e-4);
348 
349 //NO -> N + O
350  A = 5e12 * unitA_1.get_SI_factor();
351  k.push_back(Arrh(T,A,Ea,Rcal) * (Scalar(n_species) - 3. + 22.0 + 22.0 + 22.0) * 5e-4);
352 
353 //N2 + O -> NO + N
354  A = 5.7e9 * unitA_2.get_SI_factor();
355  beta = 0.42;
356  D = -5e-3;
357  k.push_back(BHE(T,A,beta,D) * (Scalar(n_species) - 3. + 22.0 + 22.0 + 22.0) * 5e-4);
358 
359 //NO + O -> NO + N
360  A = 8.4e9 * unitA_2.get_SI_factor();
361  beta = 0.4;
362  Ea = 38526.0;
363  k.push_back(Kooij(T,A,beta,Ea,Tr,Rcal) * (Scalar(n_species) - 3. + 22.0 + 22.0 + 22.0) * 5e-4);
364 
365 //C2 -> 2 C
366  A = 3.7e11 * unitA_1.get_SI_factor();
367  beta = -0.42;
368  Ea = 138812.8;
369  k.push_back(Kooij(T,A,beta,Ea,Tr,Rcal) * Scalar(n_species) * 5e-4);
370 
371 //CN -> C + N
372  A = 2.5e11 * unitA_1.get_SI_factor();
373  beta = 0.40;
374  Ea = 729372.4;
375  D = 5e-3;
376  k.push_back(VH(T,A,beta,Ea,D,Tr,Rcal) * Scalar(n_species) * 5e-4);
378 // falloff is k(T,[M]) = k0*[M]/(1 + [M]*k0/kinf) * F = k0 * ([M]^-1 + k0 * kinf^-1)^-1 * F
379 // F = 1
380 
381 // N2 -> 2 N
382  A = 7e18 * unitA_1.get_SI_factor();
383  beta = -1.6;
384  A2 = 5e15 * unitA_0.get_SI_factor();
385  beta2 = 0.5;
386  k.push_back(HE(T,A,beta) / (1./tot_dens + HE(T,A,beta)/HE(T,A2,beta2)) );
387 
388 // O2 -> 2 O
389  A = 5e17 * unitA_1.get_SI_factor();
390  D = -2.5e-5;
391  A2 = 2e18 * unitA_0.get_SI_factor();
392  D2 = -5e-3;
393  k.push_back(Bert(T,A,D) / (1./tot_dens + Bert(T,A,D)/Bert(T,A2,D2)) );
394 
395 //NO -> N + O
396  A = 5.e+12 * unitA_1.get_SI_factor();
397  Ea = 149943.0;
398  A2 = 3e+15 * unitA_0.get_SI_factor();
399  Ea2 = 200000.0;
400  k.push_back(Arrh(T,A,Ea,Rcal) / (1./tot_dens + Arrh(T,A,Ea,Rcal)/Arrh(T,A2,Ea2,Rcal)) );
401 
402 //N2 + O -> NO + N
403  A = 5e+9 * unitA_2.get_SI_factor();
404  beta = 0.6;
405  D = -5e-4;
406  A2 = 5.7e+9 * unitA_1.get_SI_factor();
407  beta2 = -0.42;
408  D2 = -5e-3;
409  k.push_back(BHE(T,A,beta,D) / (1./tot_dens + BHE(T,A,beta,D)/BHE(T,A2,beta2,D2)) );
410 
411 //NO + O -> NO + N
412  A = 8.4e+09 * unitA_2.get_SI_factor();
413  beta = 0.40;
414  Ea = 38526.0;
415  A2 = 8.4e+05 * unitA_1.get_SI_factor();
416  beta2 = 0.02;
417  Ea2 = 3526.0;
418  k.push_back(Kooij(T,A,beta,Ea,Tr,Rcal) / (1./tot_dens + Kooij(T,A,beta,Ea,Tr,Rcal)/Kooij(T,A2,beta2,Ea2,Tr,Rcal)) );
419 
420 //C2 -> 2 C
421  A = 3.7e+11 * unitA_1.get_SI_factor();
422  beta = -0.42;
423  Ea = 138812.8;
424  A2 = 3.7e+12 * unitA_0.get_SI_factor();
425  beta2 = -0.52;
426  Ea2 = 135000.8;
427  k.push_back(Kooij(T,A,beta,Ea,Tr,Rcal) / (1./tot_dens + Kooij(T,A,beta,Ea,Tr,Rcal)/Kooij(T,A2,beta2,Ea2,Tr,Rcal)) );
428 
429 //CN -> C + N
430  A = 5e+10 * unitA_1.get_SI_factor();
431  beta = -0.10;
432  D = 1.5e-3;
433  Ea = 150240.9;
434  A2 = 2.5e+11 * unitA_0.get_SI_factor();
435  beta2 = 0.40;
436  D2 = -0.005;
437  Ea2 = 174240.9;
438  k.push_back(VH(T,A,beta,Ea,D,Tr,Rcal) / (1./tot_dens + VH(T,A,beta,Ea,D,Tr,Rcal)/VH(T,A2,beta2,Ea2,D2,Tr,Rcal)) );
439 //Troe falloff
440 //falloff is k(T,[M]) = k0*[M]/(1 + [M]*k0/kinf) * F = k0 * ([M]^-1 + k0 * kinf^-1)^-1 * F
441 // F is complicated...
442  Scalar Pr,k0,kinf;
443  Scalar Fc,alpha,T1,T2,T3;
444  alpha = 0.562;
445  T1 = 5836;
446  T2 = 8552;
447  T3 = 91;
448  Fc = FcentTroe(T,alpha,T3,T1,T2);
449 
450 // N2 -> 2 N
451  A = 7.e+18 * unitA_1.get_SI_factor();
452  beta = -1.6;
453  A2 = 5.e+15 * unitA_0.get_SI_factor();
454  beta2 = 0.5;
455  k0 = HE(T,A,beta);
456  kinf = HE(T,A2,beta2);
457  Pr = tot_dens * k0/kinf;
458  k.push_back(k0 / (1./tot_dens + k0/kinf) * FTroe(Fc,Pr));
459 
460 // O2 -> 2 O
461  A = 5e17 * unitA_1.get_SI_factor();
462  D = -2.5e-5;
463  A2 = 2e18 * unitA_0.get_SI_factor();
464  D2 = -5e-3;
465  k0 = Bert(T,A,D);
466  kinf = Bert(T,A2,D2);
467  Pr = tot_dens * k0/kinf;
468  k.push_back(k0 / (1./tot_dens + k0/kinf) * FTroe(Fc,Pr));
469 
470 //NO -> N + O
471  A = 5.e+12 * unitA_1.get_SI_factor();
472  Ea = 149943.0;
473  A2 = 3e+15 * unitA_0.get_SI_factor();
474  Ea2 = 200000.0;
475  k0 = Arrh(T,A,Ea,Rcal);
476  kinf = Arrh(T,A2,Ea2,Rcal);
477  Pr = tot_dens * k0/kinf;
478  k.push_back(k0 / (1./tot_dens + k0/kinf) * FTroe(Fc,Pr));
479 
480 //N2 + O -> NO + N
481  A = 5e+9 * unitA_2.get_SI_factor();
482  beta = 0.6;
483  D = -5e-4;
484  A2 = 5.7e+9 * unitA_1.get_SI_factor();
485  beta2 = -0.42;
486  D2 = -5e-3;
487  k0 = BHE(T,A,beta,D);
488  kinf = BHE(T,A2,beta2,D2);
489  Pr = tot_dens * k0/kinf;
490  k.push_back(k0 / (1./tot_dens + k0/kinf) * FTroe(Fc,Pr));
491 
492 //NO + O -> NO + N
493  A = 8.4e+09 * unitA_2.get_SI_factor();
494  beta = 0.40;
495  Ea = 38526.0;
496  A2 = 8.4e+05 * unitA_1.get_SI_factor();
497  beta2 = 0.02;
498  Ea2 = 3526.0;
499  k0 = Kooij(T,A,beta,Ea,Tr,Rcal);
500  kinf = Kooij(T,A2,beta2,Ea2,Tr,Rcal);
501  Pr = tot_dens * k0/kinf;
502  k.push_back(k0 / (1./tot_dens + k0/kinf) * FTroe(Fc,Pr));
503 
504 //C2 -> 2 C
505  A = 3.7e+11 * unitA_1.get_SI_factor();
506  beta = -0.42;
507  Ea = 138812.8;
508  A2 = 3.7e+12 * unitA_0.get_SI_factor();
509  beta2 = -0.52;
510  Ea2 = 135000.8;
511  k0 = Kooij(T,A,beta,Ea,Tr,Rcal);
512  kinf = Kooij(T,A2,beta2,Ea2,Tr,Rcal);
513  Pr = tot_dens * k0/kinf;
514  k.push_back(k0 / (1./tot_dens + k0/kinf) * FTroe(Fc,Pr));
515 
516 //CN -> C + N
517  A = 5e+10 * unitA_1.get_SI_factor();
518  beta = -0.10;
519  D = 1.5e-3;
520  Ea = 150240.9;
521  A2 = 2.5e+11 * unitA_0.get_SI_factor();
522  beta2 = 0.40;
523  D2 = -0.005;
524  Ea2 = 174240.9;
525  k0 = VH(T,A,beta,Ea,D,Tr,Rcal);
526  kinf = VH(T,A2,beta2,Ea2,D2,Tr,Rcal);
527  Pr = tot_dens * k0/kinf;
528  k.push_back(k0 / (1./tot_dens + k0/kinf) * FTroe(Fc,Pr));
529 //
530 //photochemistry
531  k.push_back(k_photo(lambda,hv,CH4_lambda,CH4_s));
532  conditions.add_particle_flux(photons,k.size()-1);
533 //Constant
534  k.push_back(2.5e11);
535 
536 
537  const Scalar tol = (std::numeric_limits<Scalar>::epsilon() < 1e-17L)?
538  std::numeric_limits<Scalar>::epsilon() * 5000:
539  std::numeric_limits<Scalar>::epsilon() * 100;
540  int return_flag(0);
541  for(unsigned int ir = 0; ir < k.size(); ir++)
542  {
543  const Antioch::Reaction<Scalar> * reac = &reaction_set.reaction(ir);
544  if(std::abs(k[ir] - reac->compute_forward_rate_coefficient(molar_densities,conditions))/k[ir] > tol)
545  {
546  std::cout << *reac << std::endl;
547  std::cout << std::scientific << std::setprecision(16)
548  << "Error in kinetics comparison\n"
549  << "reaction #" << ir << "\n"
550  << "temperature: " << T << " K" << "\n"
551  << "theory: " << k[ir] << "\n"
552  << "calculated: " << reac->compute_forward_rate_coefficient(molar_densities,conditions) << "\n"
553  << "relative error = " << std::abs(k[ir] - reac->compute_forward_rate_coefficient(molar_densities,conditions))/k[ir] << "\n"
554  << "tolerance = " << tol
555  << std::endl;
556  return_flag = 1;
557  }
558  }
559 
560  return return_flag;
561 }
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
Antioch::ReactionSet
This class encapsulates all the reaction mechanisms considered in a chemical nonequilibrium simulatio...
Definition: kinetics_evaluator.h:42
Arrh
Scalar Arrh(const Scalar &T, const Scalar &Cf, const Scalar &Ea, const Scalar &R=Antioch::Constants::R_universal< Scalar >())
Definition: parsing_xml.C:61
Antioch::ParticleFlux
Stores the incoming flux of particles.
Definition: particle_flux.h:38
Kooij
Scalar Kooij(const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Ea, const Scalar &Tf=1., const Scalar &R=Antioch::Constants::R_universal< Scalar >())
Definition: parsing_xml.C:67
VH
Scalar VH(const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Ea, const Scalar &D, const Scalar &Tf=1., const Scalar &R=Antioch::Constants::R_universal< Scalar >())
Definition: parsing_xml.C:74
HE
Scalar HE(const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &Tf=1.)
Definition: parsing_xml.C:44
BHE
Scalar BHE(const Scalar &T, const Scalar &Cf, const Scalar &eta, const Scalar &D, const Scalar &Tf=1.)
Definition: parsing_xml.C:55
Bert
Scalar Bert(const Scalar &T, const Scalar &Cf, const Scalar &D)
Definition: parsing_xml.C:49
FcentTroe
Scalar FcentTroe(const Scalar &T, const Scalar &alpha, const Scalar &T3, const Scalar &T1, const Scalar &T2=-1.)
Definition: parsing_xml.C:81
FTroe
Scalar FTroe(const Scalar &Fcent, const Scalar &Pr)
Definition: parsing_xml.C:110
Antioch::ChemicalMixture
Class storing chemical mixture properties.
Definition: chemical_mixture.h:75
Antioch::Units
An advanced unit class.
Definition: units.h:111
Antioch::KineticsConditions
This class contains the conditions of the chemistry.
Definition: kinetics_conditions.h:64
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_xml.C:119
template<typename Scalar >
Scalar VH ( const Scalar &  T,
const Scalar &  Cf,
const Scalar &  eta,
const Scalar &  Ea,
const Scalar &  D,
const Scalar &  Tf = 1.,
const Scalar &  R = Antioch::Constants::R_universal<Scalar>() 
)

Definition at line 74 of file parsing_xml.C.

References std::pow().

Referenced by tester().

76 {
77  return Cf * std::pow(T/Tf,eta) * std::exp(-Ea /(R * T) + D * T);
78 }
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
Generated on Thu Jul 7 2016 11:09:47 for antioch-0.4.0 by  doxygen 1.8.8