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