antioch-0.4.0
Functions
kinetics_partial_order_unit.C File Reference
#include <limits>
#include <iomanip>
#include <string>
#include <vector>
#include <sstream>
#include "antioch/vector_utils.h"
#include "antioch/antioch_asserts.h"
#include "antioch/chemical_species.h"
#include "antioch/chemical_mixture.h"
#include "antioch/reaction_set.h"
#include "antioch/read_reaction_set_data.h"
#include "antioch/nasa_mixture_parsing.h"
#include "antioch/nasa_evaluator.h"
#include "antioch/kinetics_evaluator.h"

Go to the source code of this file.

Functions

template<typename Scalar >
std::vector< Scalar > temp_thermo ()
 
template<typename Scalar >
std::vector< Scalar > H_thermo (const Scalar &T)
 
template<typename Scalar >
std::vector< Scalar > H2_thermo (const Scalar &T)
 
template<typename Scalar >
std::vector< Scalar > H2O2_thermo (const Scalar &T)
 
template<typename Scalar >
std::vector< Scalar > HO2_thermo (const Scalar &T)
 
template<typename Scalar >
Scalar g (const std::vector< Scalar > &thermo, const Scalar &T)
 
template<typename Scalar >
Scalar G (const Scalar &T)
 
template<typename Scalar >
int check_test (const Scalar &exact, const Scalar &cal, const std::string &words)
 
template<typename Scalar >
int checker (const Scalar &net_rates_exact, const Scalar &kfwd_const_exact, const Scalar &kfwd_exact, const Scalar &fwd_conc_exact, const Scalar &kbkwd_const_exact, const Scalar &kbkwd_exact, const Scalar &bkwd_conc_exact, const Scalar &net_rates, const Scalar &kfwd_const, const Scalar &kfwd, const Scalar &fwd_conc, const Scalar &kbkwd_const, const Scalar &kbkwd, const Scalar &bkwd_conc, const Scalar &Temp)
 
template<typename Scalar >
int test_type (const std::string &input_name, const Antioch::ParsingType &inputType)
 
template<typename Scalar >
int tester (const std::string &input_name_xml, const std::string &input_name_ck)
 
int main (int argc, char *argv[])
 

Function Documentation

template<typename Scalar >
int check_test ( const Scalar &  exact,
const Scalar &  cal,
const std::string &  words 
)

Definition at line 224 of file kinetics_partial_order_unit.C.

Referenced by checker().

225 {
226  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 250;
227 
228  if(std::abs(exact - cal)/exact > tol)
229  {
230  std::cout << std::scientific << std::setprecision(20)
231  << "Erreur in tests of " << words << "\n"
232  << "Calculated value is " << cal << "\n"
233  << "Exact value is " << exact << "\n"
234  << "Relative error is " << std::abs(exact - cal)/exact << "\n"
235  << "Tolerance is " << tol << std::endl;
236  return 1;
237  }
238  return 0;
239 }
template<typename Scalar >
int checker ( const Scalar &  net_rates_exact,
const Scalar &  kfwd_const_exact,
const Scalar &  kfwd_exact,
const Scalar &  fwd_conc_exact,
const Scalar &  kbkwd_const_exact,
const Scalar &  kbkwd_exact,
const Scalar &  bkwd_conc_exact,
const Scalar &  net_rates,
const Scalar &  kfwd_const,
const Scalar &  kfwd,
const Scalar &  fwd_conc,
const Scalar &  kbkwd_const,
const Scalar &  kbkwd,
const Scalar &  bkwd_conc,
const Scalar &  Temp 
)

Definition at line 242 of file kinetics_partial_order_unit.C.

References check_test().

Referenced by test_type().

249 {
250 
251  int return_flag(0);
252 
253  std::stringstream os;
254  os << Temp << "K";
255 
256  return_flag = check_test(fwd_conc_exact,fwd_conc,"concentrations forward at " + os.str()) ||
257  return_flag;
258 
259  return_flag = check_test(kfwd_const_exact,kfwd_const,"rate constant forward at " + os.str()) ||
260  return_flag;
261 
262  return_flag = check_test(kfwd_exact,kfwd,"rate forward at " + os.str()) ||
263  return_flag;
264 
265  return_flag = check_test(bkwd_conc_exact,bkwd_conc,"concentrations backward at " + os.str()) ||
266  return_flag;
267 
268  return_flag = check_test(kbkwd_const_exact,kbkwd_const,"rate constant backward at " + os.str()) ||
269  return_flag;
270 
271  return_flag = check_test(kbkwd_exact,kbkwd,"rate backward at " + os.str()) ||
272  return_flag;
273 
274  return_flag = check_test(net_rates_exact,net_rates,"net rate at " + os.str()) ||
275  return_flag;
276 
277  return return_flag;
278 }
int check_test(const Scalar &exact, const Scalar &cal, const std::string &words)
template<typename Scalar >
Scalar g ( const std::vector< Scalar > &  thermo,
const Scalar &  T 
)

Definition at line 204 of file kinetics_partial_order_unit.C.

Referenced by G(), and test_g().

205 {
206  return - thermo[0] / ( 2 * T * T)
207  + thermo[1] * ( 1 + Antioch::ant_log(T) ) / T
208  + thermo[2] * ( 1 - Antioch::ant_log(T) )
209  - thermo[3] * T / 2
210  - thermo[4] * T * T / 6
211  - thermo[5] * T * T * T / 12
212  - thermo[6] * T * T * T * T / 20
213  + thermo[8] / T
214  - thermo[9];
215 }
template<typename Scalar >
Scalar G ( const Scalar &  T)

Definition at line 218 of file kinetics_partial_order_unit.C.

References g(), H2_thermo(), H2O2_thermo(), H_thermo(), and HO2_thermo().

Referenced by test_type().

219 {
220  return g(H2O2_thermo(T),T) + g(H_thermo(T),T) - g(HO2_thermo(T),T) - g(H2_thermo(T),T);
221 }
std::vector< Scalar > H_thermo(const Scalar &T)
std::vector< Scalar > HO2_thermo(const Scalar &T)
std::vector< Scalar > H2_thermo(const Scalar &T)
std::vector< Scalar > H2O2_thermo(const Scalar &T)
Scalar g(const std::vector< Scalar > &thermo, const Scalar &T)
template<typename Scalar >
std::vector<Scalar> H2_thermo ( const Scalar &  T)

Definition at line 103 of file kinetics_partial_order_unit.C.

Referenced by G().

104 {
105  std::vector<std::vector<Scalar> > out(3,std::vector<Scalar>(10,0));
106 
107  out[0][0] = 4.07832281e+04;
108  out[0][1] = -8.00918545e+02;
109  out[0][2] = 8.21470167e+00;
110  out[0][3] = -1.26971436e-02;
111  out[0][4] = 1.75360493e-05;
112  out[0][5] = -1.20286016e-08;
113  out[0][6] = 3.36809316e-12;
114  out[0][7] = 0.00000000e+00;
115  out[0][8] = 2.68248438e+03;
116  out[0][9] = -3.04378866e+01;
117 
118  out[1][0] = 5.60812338e+05;
119  out[1][1] = -8.37149134e+02;
120  out[1][2] = 2.97536304e+00;
121  out[1][3] = 1.25224993e-03;
122  out[1][4] = -3.74071842e-07;
123  out[1][5] = 5.93662825e-11;
124  out[1][6] = -3.60699573e-15;
125  out[1][7] = 0.00000000e+00;
126  out[1][8] = 5.33981585e+03;
127  out[1][9] = -2.20276405e+00;
128 
129  out[2][0] = 4.96671613e+08;
130  out[2][1] = -3.14744812e+05;
131  out[2][2] = 7.98388750e+01;
132  out[2][3] = -8.41450419e-03;
133  out[2][4] = 4.75306044e-07;
134  out[2][5] = -1.37180973e-11;
135  out[2][6] = 1.60537460e-16;
136  out[2][7] = 0.00000000e+00;
137  out[2][8] = 2.48835466e+06;
138  out[2][9] = -6.69552419e+02;
139 
140  return (T <= temp_thermo<Scalar>()[1])?out[0]:(T <= temp_thermo<Scalar>()[2])?out[1]:out[2];
141 }
template<typename Scalar >
std::vector<Scalar> H2O2_thermo ( const Scalar &  T)

Definition at line 144 of file kinetics_partial_order_unit.C.

Referenced by G().

145 {
146  std::vector<std::vector<Scalar> > out(2,std::vector<Scalar>(10,0));
147 
148  out[0][0] = -9.279533580E+04;
149  out[0][1] = 1.564748385E+03;
150  out[0][2] = -5.976460140E+00;
151  out[0][3] = 3.270744520E-02;
152  out[0][4] = -3.932193260E-05;
153  out[0][5] = 2.509255235E-08;
154  out[0][6] = -6.465045290E-12;
155  out[0][7] = 0.000000000E+00;
156  out[0][8] = -2.494004728E+04;
157  out[0][9] = 5.877174180E+01;
158 
159  out[1][0] = 1.489428027E+06;
160  out[1][1] = -5.170821780E+03;
161  out[1][2] = 1.128204970E+01;
162  out[1][3] = -8.042397790E-05;
163  out[1][4] = -1.818383769E-08;
164  out[1][5] = 6.947265590E-12;
165  out[1][6] = -4.827831900E-16;
166  out[1][7] = 0.000000000E+00;
167  out[1][8] = 1.418251038E+04;
168  out[1][9] = -4.650855660E+01;
169 
170  return (T <= temp_thermo<Scalar>()[1])?out[0]:out[1];
171 }
template<typename Scalar >
std::vector<Scalar> H_thermo ( const Scalar &  T)

Definition at line 63 of file kinetics_partial_order_unit.C.

Referenced by G().

64 {
65  std::vector<std::vector<Scalar> > out(3,std::vector<Scalar>(10,0));
66  out[0][0] = 0.00000000e+00;
67  out[0][1] = 0.00000000e+00;
68  out[0][2] = 2.50000000e+00;
69  out[0][3] = 0.00000000e+00;
70  out[0][4] = 0.00000000e+00;
71  out[0][5] = 0.00000000e+00;
72  out[0][6] = 0.00000000e+00;
73  out[0][7] = 0.00000000e+00;
74  out[0][8] = 2.54737080e+04;
75  out[0][9] = -4.46682853e-01;
76 
77  out[1][0] = 6.07877425e+01;
78  out[1][1] = -1.81935442e-01;
79  out[1][2] = 2.50021182e+00;
80  out[1][3] = -1.22651286e-07;
81  out[1][4] = 3.73287633e-11;
82  out[1][5] = -5.68774456e-15;
83  out[1][6] = 3.41021020e-19;
84  out[1][7] = 0.00000000e+00;
85  out[1][8] = 2.54748640e+04;
86  out[1][9] = -4.48191777e-01;
87 
88  out[2][0] = 2.17375769e+08;
89  out[2][1] = -1.31203540e+05;
90  out[2][2] = 3.39917420e+01;
91  out[2][3] = -3.81399968e-03;
92  out[2][4] = 2.43285484e-07;
93  out[2][5] = -7.69427554e-12;
94  out[2][6] = 9.64410563e-17;
95  out[2][7] = 0.00000000e+00;
96  out[2][8] = 1.06763809e+06;
97  out[2][9] = -2.74230105e+02;
98 
99  return (T <= temp_thermo<Scalar>()[1])?out[0]:(T <= temp_thermo<Scalar>()[2])?out[1]:out[2];
100 }
template<typename Scalar >
std::vector<Scalar> HO2_thermo ( const Scalar &  T)

Definition at line 174 of file kinetics_partial_order_unit.C.

Referenced by G().

175 {
176  std::vector<std::vector<Scalar> > out(2,std::vector<Scalar>(10,0));
177 
178  out[0][0] = -7.598882540E+04;
179  out[0][1] = 1.329383918E+03;
180  out[0][2] = -4.677388240E+00;
181  out[0][3] = 2.508308202E-02;
182  out[0][4] = -3.006551588E-05;
183  out[0][5] = 1.895600056E-08;
184  out[0][6] = -4.828567390E-12;
185  out[0][7] = 0.000000000E+00;
186  out[0][8] = -5.873350960E+03;
187  out[0][9] = 5.193602140E+01;
188 
189  out[1][0] = -1.810669724E+06;
190  out[1][1] = 4.963192030E+03;
191  out[1][2] = -1.039498992E+00;
192  out[1][3] = 4.560148530E-03;
193  out[1][4] = -1.061859447E-06;
194  out[1][5] = 1.144567878E-10;
195  out[1][6] = -4.763064160E-15;
196  out[1][7] = 0.000000000E+00;
197  out[1][8] = -3.200817190E+04;
198  out[1][9] = 4.066850920E+01;
199 
200  return (T <= temp_thermo<Scalar>()[1])?out[0]:out[1];
201 }
int main ( int  argc,
char *  argv[] 
)

Definition at line 400 of file kinetics_partial_order_unit.C.

References antioch_error.

401 {
402  // Check command line count.
403  if( argc < 3 )
404  {
405  // TODO: Need more consistent error handling.
406  std::cerr << "Error: Must specify reaction set XML input file and reaction set ChemKin input file." << std::endl;
407  antioch_error();
408  }
409 
410  return (tester<float>(std::string(argv[1]), std::string(argv[2]) ) ||
411  tester<double>(std::string(argv[1]), std::string(argv[2]) ) /*||
412  tester<long double>(std::string(argv[1]), std::string(argv[2]) )*/
413  );
414 }
#define antioch_error()
template<typename Scalar >
std::vector<Scalar> temp_thermo ( )

Definition at line 53 of file kinetics_partial_order_unit.C.

54 {
55  std::vector<Scalar> Ts(3,0);
56  Ts[0] = 200;
57  Ts[1] = 1000;
58  Ts[2] = 6000;
59  return Ts;
60 }
template<typename Scalar >
int test_type ( const std::string &  input_name,
const Antioch::ParsingType inputType 
)

Definition at line 281 of file kinetics_partial_order_unit.C.

References checker(), G(), Antioch::NASAEvaluator< CoeffType, NASAFit >::h_RT_minus_s_R(), Antioch::ChemicalMixture< CoeffType >::molar_densities, Antioch::ChemicalMixture< CoeffType >::R(), and Antioch::read_nasa_mixture_data().

282 {
283  using std::abs;
284 
285  std::vector<std::string> species_str_list;
286  const unsigned int n_species = 4;
287  species_str_list.reserve(n_species);
288  species_str_list.push_back( "H" );
289  species_str_list.push_back( "H2" );
290  species_str_list.push_back( "H2O2" );
291  species_str_list.push_back( "HO2" );
292 
293  // kinetics
294  Antioch::ChemicalMixture<Scalar> chem_mixture( species_str_list );
295  Antioch::ReactionSet<Scalar> reaction_set( chem_mixture );
296  Antioch::read_reaction_set_data<Scalar>( input_name, true, reaction_set, inputType );
297 
298  // thermo
300  Antioch::read_nasa_mixture_data( thermo_mixture ); // default
302 
303  const Scalar T = 1500.0; // K
304  const Scalar P = 1.0e5; // Pa
305 
306  // Mass fractions
307  std::vector<Scalar> Y(n_species,0.25);
308 
309  const Scalar R_mix = chem_mixture.R(Y); // get R_tot in J.kg-1.K-1
310 
311  std::vector<Scalar> molar_densities(n_species,0);
312 
313  std::vector<Scalar> h_RT_minus_s_R(n_species);
314 
315  Scalar net_rates_exact(0),
316  kfwd_const_exact(0),
317  kbkwd_const_exact(0),
318  kfwd_exact(0),
319  kbkwd_exact(0),
320  fwd_conc_exact(0),
321  bkwd_conc_exact(0);
322 
323  std::vector<Scalar> net_rates,kfwd_const,kbkwd_const,kfwd,kbkwd,fwd_conc,bkwd_conc;
324 
325  // golden values from bc
326  kfwd_const_exact = 3341803.012517167957974472239343341385324934527797794;
327  fwd_conc_exact = 0.13473900180907885567483372493084259843941997896096;
328  kfwd_exact = 450271.20214913586326479266723410049955538508784467327734;
329  bkwd_conc_exact = 0.06329787652743180276417606259788367794531198101580;
330  kbkwd_const_exact = 4327.11114439947889901291861278625576793065735934819023;
331  kbkwd_exact = 273.89694693867234146591036506821238311545088830214564;
332  net_rates_exact = 449997.30520219719092332675686903228717226963695635680885;
333 
334 
335  Scalar rho = P/(R_mix*T); // kg.m-3
336  chem_mixture.molar_densities(rho,Y,molar_densities);
337  const Antioch::KineticsConditions<Scalar> conditions(T);
338  const Antioch::TempCache<Scalar> Cache(T);
339  thermo.h_RT_minus_s_R(Cache,h_RT_minus_s_R);
340  reaction_set.get_reactive_scheme(conditions,molar_densities,h_RT_minus_s_R,net_rates,
341  kfwd_const,kbkwd_const,kfwd,kbkwd,fwd_conc,bkwd_conc);
342 
343  int return_flag = checker(net_rates_exact, kfwd_const_exact, kfwd_exact, fwd_conc_exact, kbkwd_const_exact, kbkwd_exact, bkwd_conc_exact,
344  net_rates[0], kfwd_const[0], kfwd[0], fwd_conc[0], kbkwd_const[0], kbkwd[0], bkwd_conc[0], T);
345 
346  const Scalar Rcal = Antioch::Constants::R_universal<Scalar>() * Antioch::Constants::R_universal_unit<Scalar>().factor_to_some_unit("cal/mol/K");
347  const Scalar fac(1e-6);
348 
349  for(Scalar Temp = 210; Temp < 5990; Temp += 10){
350 
351  rho = P/(R_mix*Temp); // kg.m-3
352 
353  molar_densities.resize(4,0);
354  chem_mixture.molar_densities(rho,Y,molar_densities);
355 
356  const Antioch::KineticsConditions<Scalar> conditionsTemp(Temp);
357  const Antioch::TempCache<Scalar> CacheTemp(Temp);
358  h_RT_minus_s_R.resize(4,0);
359  thermo.h_RT_minus_s_R(CacheTemp,h_RT_minus_s_R);
360 
361  net_rates.clear();
362  kfwd_const.clear();
363  kbkwd_const.clear();
364  kfwd.clear();
365  kbkwd.clear();
366  fwd_conc.clear();
367  bkwd_conc.clear();
368  reaction_set.get_reactive_scheme(conditionsTemp,molar_densities,h_RT_minus_s_R,net_rates,
369  kfwd_const,kbkwd_const,kfwd,kbkwd,fwd_conc,bkwd_conc);
370 
371 // H2O2 + H <=> HO2 + H2
372 // k = Arrhenius(4.82e13 (cm3/mol/s), 7.95e3 (cal/mol) )
373 // m_H2O2 = 1.5
374 // m_H = 0.5
375 // m_HO2 = 2
376 // m_H2 = 1
377  fwd_conc_exact = Antioch::ant_pow(molar_densities[2],1.5) * Antioch::ant_pow(molar_densities[0],0.5);
378  bkwd_conc_exact = Antioch::ant_pow(molar_densities[3],2) * molar_densities[1];
379  kfwd_const_exact = 4.82e13 * fac * std::exp(-7.95e3/(Rcal * Temp)); // Arrhenius
380  kfwd_exact = kfwd_const_exact * fwd_conc_exact;
381  kbkwd_const_exact = kfwd_const_exact / Antioch::ant_exp(G(Temp));
382  kbkwd_exact = kbkwd_const_exact * bkwd_conc_exact;
383  net_rates_exact = kfwd_exact - kbkwd_exact;
384 
385  return_flag = checker(net_rates_exact, kfwd_const_exact, kfwd_exact, fwd_conc_exact, kbkwd_const_exact, kbkwd_exact, bkwd_conc_exact,
386  net_rates[0], kfwd_const[0], kfwd[0], fwd_conc[0], kbkwd_const[0], kbkwd[0], bkwd_conc[0], Temp) ||
387  return_flag;
388  }
389 
390  return return_flag;
391 }
This class encapsulates all the reaction mechanisms considered in a chemical nonequilibrium simulatio...
int checker(const Scalar &net_rates_exact, const Scalar &kfwd_const_exact, const Scalar &kfwd_exact, const Scalar &fwd_conc_exact, const Scalar &kbkwd_const_exact, const Scalar &kbkwd_exact, const Scalar &bkwd_conc_exact, const Scalar &net_rates, const Scalar &kfwd_const, const Scalar &kfwd, const Scalar &fwd_conc, const Scalar &kbkwd_const, const Scalar &kbkwd, const Scalar &bkwd_conc, const Scalar &Temp)
void read_nasa_mixture_data(NASAThermoMixture< NumericType, CurveType > &thermo, const std::string &filename=DefaultSourceFilename::thermo_data(), ParsingType=ASCII, bool verbose=true)
Class storing chemical mixture properties.
Scalar G(const Scalar &T)
This class contains the conditions of the chemistry.
template<typename Scalar >
int tester ( const std::string &  input_name_xml,
const std::string &  input_name_ck 
)

Definition at line 394 of file kinetics_partial_order_unit.C.

References Antioch::CHEMKIN, and Antioch::XML.

395 {
396  return test_type<Scalar>(input_name_xml, Antioch::ParsingType::XML) ||
397  test_type<Scalar>(input_name_ck, Antioch::ParsingType::CHEMKIN);
398 }

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