antioch-0.4.0
kinetics_reactive_scheme_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 // $Id$
28 //
29 //--------------------------------------------------------------------------
30 //--------------------------------------------------------------------------
31 
32 // C++
33 #include <limits>
34 #include <iomanip>
35 #include <string>
36 #include <vector>
37 
38 // Antioch
39 #include "antioch/vector_utils.h"
40 
44 #include "antioch/reaction_set.h"
46 #include "antioch/cea_thermo.h"
48 
49 
50 template <typename Scalar>
51 int check_test(const Scalar &exact,const Scalar &cal,const std::string &words, unsigned int rxn)
52 {
53  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 100;
54 
55  if(std::abs(exact - cal)/exact > tol)
56  {
57  std::cout << std::scientific << std::setprecision(20)
58  << "Erreur in tests of " << words << "\n"
59  << "of reaction #" << rxn << "\n"
60  << "Calculated value is " << cal << "\n"
61  << "Exact value is " << exact << "\n"
62  << "Relative error is " << std::abs(exact - cal)/exact << "\n"
63  << "Tolerance is " << tol << std::endl;
64  return 1;
65  }
66  return 0;
67 }
68 
69 template <typename Scalar>
70 int tester(const std::string& input_name)
71 {
72  using std::abs;
73 
74  std::vector<std::string> species_str_list;
75  const unsigned int n_species = 5;
76  species_str_list.reserve(n_species);
77  species_str_list.push_back( "N2" );
78  species_str_list.push_back( "O2" );
79  species_str_list.push_back( "N" );
80  species_str_list.push_back( "O" );
81  species_str_list.push_back( "NO" );
82 
83  Antioch::ChemicalMixture<Scalar> chem_mixture( species_str_list );
84  Antioch::ReactionSet<Scalar> reaction_set( chem_mixture );
85  Antioch::CEAThermodynamics<Scalar> thermo( chem_mixture );
86 
87  Antioch::read_reaction_set_data_xml<Scalar>( input_name, true, reaction_set );
88 
89  const Scalar T = 1500.0; // K
90  const Scalar P = 1.0e5; // Pa
91 
92  const Antioch::KineticsConditions<Scalar> conditions(T);
93 
94  // Mass fractions
95  std::vector<Scalar> Y(n_species,0.2);
96 
97  const Scalar R_mix = chem_mixture.R(Y); // get R_tot in J.kg-1.K-1
98 
99  const Scalar rho = P/(R_mix*T); // kg.m-3
100 
101  std::vector<Scalar> molar_densities(n_species,0.0);
102  chem_mixture.molar_densities(rho,Y,molar_densities);
103 
104  std::vector<Scalar> h_RT_minus_s_R(n_species);
105 
106  typedef typename Antioch::CEAThermodynamics<Scalar>::template Cache<Scalar> Cache;
107  thermo.h_RT_minus_s_R(Cache(T),h_RT_minus_s_R);
108 
109  std::vector<Scalar> net_rates,kfwd_const,kbkwd_const,kfwd,kbkwd,fwd_conc,bkwd_conc;
110 
111  reaction_set.get_reactive_scheme(conditions,molar_densities,h_RT_minus_s_R,net_rates,
112  kfwd_const,kbkwd_const,kfwd,kbkwd,fwd_conc,bkwd_conc);
113 
114  std::vector<Scalar> net_rates_exact,
115  kfwd_const_exact,
116  kbkwd_const_exact,
117  kfwd_exact,
118  kbkwd_exact,
119  fwd_conc_exact,
120  bkwd_conc_exact;
121  net_rates_exact.resize(5,0.L);
122  kfwd_const_exact.resize(5,0.L);
123  kbkwd_const_exact.resize(5,0.L);
124  kfwd_exact.resize(5,0.L);
125  kbkwd_exact.resize(5,0.L);
126  fwd_conc_exact.resize(5,1.L);
127  bkwd_conc_exact.resize(5,1.L);
128 
129  const Scalar Rcal = Antioch::Constants::R_universal<Scalar>() * Antioch::Constants::R_universal_unit<Scalar>().factor_to_some_unit("cal/mol/K");
130  const Scalar P0_RT(1.0e5/Antioch::Constants::R_universal<Scalar>()/T);
131 //N2 + M <=> 2 N + M
132  kfwd_const_exact[0] = 7e15L * std::pow(T,-1.6L) * std::exp(-224801.3L/(Rcal * T)) * //Kooij
133  (4.2857L * (molar_densities[2] + molar_densities[3]) // N and O
134  + molar_densities[0] + molar_densities[1] + molar_densities[4] ); //N2, O2, NO
135  fwd_conc_exact[0] = molar_densities[0]; //N2
136  bkwd_conc_exact[0] = std::pow(molar_densities[2],2); //2 N
137 
138 //O2 + M <=> 2 O + M
139  kfwd_const_exact[1] = 2e15L * std::pow(T,-1.5L) * std::exp(-117881.7L/(Rcal * T)) * //Kooij
140  (5.0L * (molar_densities[2] + molar_densities[3]) // N and O
141  + molar_densities[0] + molar_densities[1] + molar_densities[4] ); //N2, O2, NO
142  fwd_conc_exact[1] = molar_densities[1]; //O2
143  bkwd_conc_exact[1] = std::pow(molar_densities[3],2); //2 O
144 
145 //NO + M <=> N + O + M
146  kfwd_const_exact[2] = 5e9L * std::exp(-149943.0L/(Rcal * T)) * //Kooij
147  (22.0L * (molar_densities[2] + molar_densities[3] + molar_densities[4]) // N, O and NO
148  + molar_densities[0] + molar_densities[1]); //N2, O2
149  fwd_conc_exact[2] = molar_densities[4]; //NO
150  bkwd_conc_exact[2] = molar_densities[2] * molar_densities[3]; //N + O
151 
152 //N2 + O => NO + N
153  kfwd_const_exact[3] = 5.7e9L * std::pow(T,0.42) * std::exp(-85269.6L/(Rcal * T)); //Kooij
154  fwd_conc_exact[3] = molar_densities[0] * molar_densities[3]; //N2 + O
155 
156 //NO + O => O2 + N
157  kfwd_const_exact[4] = 8.4e9L * std::exp(-38526.0L/(Rcal * T)); //Kooij
158  fwd_conc_exact[4] = molar_densities[4] * molar_densities[3]; //NO + O
159 
160  for(unsigned int rxn = 0; rxn < 5; rxn++)
161  {
162  if(rxn < 3)
163  {
164  kbkwd_const_exact[rxn] = kfwd_const_exact[rxn]/reaction_set.reaction(rxn).equilibrium_constant(P0_RT,h_RT_minus_s_R);
165  kbkwd_exact[rxn] = kbkwd_const_exact[rxn] * bkwd_conc_exact[rxn];
166  }
167  kfwd_exact[rxn] = kfwd_const_exact[rxn] * fwd_conc_exact[rxn];
168  net_rates_exact[rxn] = kfwd_exact[rxn] - kbkwd_exact[rxn];
169  }
170 
171  int return_flag(0);
172  for(unsigned int rxn = 0; rxn < 5; rxn++)
173  {
174  return_flag = check_test(net_rates_exact[rxn],net_rates[rxn],"net rates",rxn) ||
175  return_flag;
176  return_flag = check_test(kfwd_const_exact[rxn],kfwd_const[rxn],"rate constant forward",rxn) ||
177  return_flag;
178  return_flag = check_test(kfwd_exact[rxn],kfwd[rxn],"rate forward",rxn) ||
179  return_flag;
180  return_flag = check_test(fwd_conc_exact[rxn],fwd_conc[rxn],"concentrations forward",rxn) ||
181  return_flag;
182 
183  if(rxn < 3)
184  {
185  return_flag = check_test(kbkwd_const_exact[rxn],kbkwd_const[rxn],"rate constant backward",rxn) ||
186  return_flag;
187  return_flag = check_test(kbkwd_exact[rxn],kbkwd[rxn],"rate backward",rxn) ||
188  return_flag;
189  return_flag = check_test(bkwd_conc_exact[rxn],bkwd_conc[rxn],"concentrations backward",rxn) ||
190  return_flag;
191  }
192  }
193  return return_flag;
194 }
195 
196 int main(int argc, char* argv[])
197 {
198  // Check command line count.
199  if( argc < 2 )
200  {
201  // TODO: Need more consistent error handling.
202  std::cerr << "Error: Must specify reaction set XML input file." << std::endl;
203  antioch_error();
204  }
205 
206  return (tester<float>(std::string(argv[1])) ||
207  tester<double>(std::string(argv[1])));/* ||
208  tester<long double>(std::string(argv[1]))
209  ); */
210 }
This class encapsulates all the reaction mechanisms considered in a chemical nonequilibrium simulatio...
CoeffType R(const unsigned int s) const
Gas constant for species s in [J/kg-K].
const Reaction< CoeffType > & reaction(const unsigned int r) const
Definition: reaction_set.h:232
#define antioch_error()
Antioch::enable_if_c< Antioch::is_valarray< T >::value, typename Antioch::state_type< T >::type >::type pow(const T &in, const T2 &n)
int tester(const std::string &input_name)
int check_test(const Scalar &exact, const Scalar &cal, const std::string &words, unsigned int rxn)
X(const unsigned int species, const StateType &M, const StateType &mass_fraction) const template< typename VectorStateType > void X(typename Antioch VectorStateType void molar_densities(const StateType &rho, const VectorStateType &mass_fractions, VectorStateType &molar_densities) const
Class storing chemical mixture properties.
StateType h_RT_minus_s_R(const Cache< StateType > &cache, unsigned int species) const
We currently need different specializations for scalar vs vector inputs here.
Definition: cea_thermo.h:420
int main(int argc, char *argv[])
This class contains the conditions of the chemistry.
void get_reactive_scheme(const KineticsConditions< StateType, VectorStateType > &conditions, const VectorStateType &molar_densities, const VectorStateType &h_RT_minus_s_R, VectorStateType &net_rates, VectorStateType &kfwd_const, VectorStateType &kbkwd_const, VectorStateType &kfwd, VectorStateType &kbkwd, VectorStateType &fwd_conc, VectorStateType &bkwd_conc) const
Definition: reaction_set.h:802
We parse the file here, with an exhaustive unit management.

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