Cantera  2.4.0
InterfaceKinetics.cpp
Go to the documentation of this file.
1 /**
2  * @file InterfaceKinetics.cpp
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at http://www.cantera.org/license.txt for license and copyright information.
7 
12 
13 #include <cstdio>
14 
15 using namespace std;
16 
17 namespace Cantera
18 {
19 
20 InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
21  m_redo_rates(false),
22  m_surf(0),
23  m_integrator(0),
24  m_ROP_ok(false),
25  m_temp(0.0),
26  m_logtemp(0.0),
27  m_has_coverage_dependence(false),
28  m_has_electrochem_rxns(false),
29  m_has_exchange_current_density_formulation(false),
30  m_phaseExistsCheck(false),
31  m_ioFlag(0),
32  m_nDim(2)
33 {
34  if (thermo != 0) {
35  addPhase(*thermo);
36  }
37 }
38 
39 InterfaceKinetics::~InterfaceKinetics()
40 {
41  delete m_integrator;
42 }
43 
44 void InterfaceKinetics::setElectricPotential(int n, doublereal V)
45 {
47  m_redo_rates = true;
48 }
49 
51 {
52  // First task is update the electrical potentials from the Phases
55  m_surf->getCoverages(m_actConc.data());
56  m_rates.update_C(m_actConc.data());
57  m_redo_rates = true;
58  }
59 
60  // Go find the temperature from the surface
61  doublereal T = thermo(surfacePhaseIndex()).temperature();
62  m_redo_rates = true;
63  if (T != m_temp || m_redo_rates) {
64  m_logtemp = log(T);
65 
66  // Calculate the forward rate constant by calling m_rates and store it in m_rfn[]
67  m_rates.update(T, m_logtemp, m_rfn.data());
68  applyStickingCorrection(T, m_rfn.data());
69 
70  // If we need to do conversions between exchange current density
71  // formulation and regular formulation (either way) do it here.
74  }
77  }
78  m_temp = T;
79  updateKc();
80  m_ROP_ok = false;
81  m_redo_rates = false;
82  }
83 }
84 
86 {
87  // Store electric potentials for each phase in the array m_phi[].
88  for (size_t n = 0; n < nPhases(); n++) {
89  if (thermo(n).electricPotential() != m_phi[n]) {
90  m_phi[n] = thermo(n).electricPotential();
91  m_redo_rates = true;
92  }
93  }
94 }
95 
97 {
98  for (size_t n = 0; n < nPhases(); n++) {
99  const ThermoPhase* tp = m_thermo[n];
100  /*
101  * We call the getActivityConcentrations function of each ThermoPhase
102  * class that makes up this kinetics object to obtain the generalized
103  * concentrations for species within that class. This is collected in
104  * the vector m_conc. m_start[] are integer indices for that vector
105  * denoting the start of the species for each phase.
106  */
108 
109  // Get regular concentrations too
110  tp->getConcentrations(m_conc.data() + m_start[n]);
111  }
112  m_ROP_ok = false;
113 }
114 
116 {
117  _update_rates_C();
118  copy(m_actConc.begin(), m_actConc.end(), conc);
119 }
120 
122 {
123  fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
124 
125  if (m_revindex.size() > 0) {
126  /*
127  * Get the vector of standard state electrochemical potentials for
128  * species in the Interfacial kinetics object and store it in m_mu0[]
129  * and m_mu0_Kc[]
130  */
131  updateMu0();
132  doublereal rrt = 1.0 / thermo(0).RT();
133 
134  // compute Delta mu^0 for all reversible reactions
135  getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data());
136 
137  for (size_t i = 0; i < m_revindex.size(); i++) {
138  size_t irxn = m_revindex[i];
139  if (irxn == npos || irxn >= nReactions()) {
140  throw CanteraError("InterfaceKinetics", "illegal value: irxn = {}", irxn);
141  }
142  // WARNING this may overflow HKM
143  m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);
144  }
145  for (size_t i = 0; i != m_irrev.size(); ++i) {
146  m_rkcn[ m_irrev[i] ] = 0.0;
147  }
148  }
149 }
150 
152 {
153  // First task is update the electrical potentials from the Phases
155 
157  size_t ik = 0;
158  for (size_t n = 0; n < nPhases(); n++) {
160  for (size_t k = 0; k < thermo(n).nSpecies(); k++) {
161  m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k);
162  m_mu0_Kc[ik] -= thermo(0).RT() * thermo(n).logStandardConc(k);
163  ik++;
164  }
165  }
166 }
167 
169 {
170  updateMu0();
171  doublereal rrt = 1.0 / thermo(0).RT();
172  std::fill(kc, kc + nReactions(), 0.0);
173  getReactionDelta(m_mu0_Kc.data(), kc);
174  for (size_t i = 0; i < nReactions(); i++) {
175  kc[i] = exp(-kc[i]*rrt);
176  }
177 }
178 
180 {
181  // Calculate:
182  // - m_StandardConc[]
183  // - m_ProdStanConcReac[]
184  // - m_deltaG0[]
185  // - m_mu0[]
186 
187  // First collect vectors of the standard Gibbs free energies of the
188  // species and the standard concentrations
189  // - m_mu0
190  // - m_StandardConc
191  size_t ik = 0;
192 
193  for (size_t n = 0; n < nPhases(); n++) {
195  size_t nsp = thermo(n).nSpecies();
196  for (size_t k = 0; k < nsp; k++) {
198  ik++;
199  }
200  }
201 
202  getReactionDelta(m_mu0.data(), m_deltaG0.data());
203 
204  // Calculate the product of the standard concentrations of the reactants
205  for (size_t i = 0; i < nReactions(); i++) {
206  m_ProdStanConcReac[i] = 1.0;
207  }
208  m_reactantStoich.multiply(m_StandardConc.data(), m_ProdStanConcReac.data());
209 }
210 
212 {
213  // Compute the electrical potential energy of each species
214  size_t ik = 0;
215  for (size_t n = 0; n < nPhases(); n++) {
216  size_t nsp = thermo(n).nSpecies();
217  for (size_t k = 0; k < nsp; k++) {
218  m_pot[ik] = Faraday * thermo(n).charge(k) * m_phi[n];
219  ik++;
220  }
221  }
222 
223  // Compute the change in electrical potential energy for each reaction. This
224  // will only be non-zero if a potential difference is present.
226 
227  // Modify the reaction rates. Only modify those with a non-zero activation
228  // energy. Below we decrease the activation energy below zero but in some
229  // debug modes we print out a warning message about this.
230 
231  // NOTE, there is some discussion about this point. Should we decrease the
232  // activation energy below zero? I don't think this has been decided in any
233  // definitive way. The treatment below is numerically more stable, however.
234  for (size_t i = 0; i < m_beta.size(); i++) {
235  size_t irxn = m_ctrxn[i];
236 
237  // If we calculate the BV form directly, we don't add the voltage
238  // correction to the forward reaction rate constants.
239  if (m_ctrxn_BVform[i] == 0) {
240  double eamod = m_beta[i] * deltaElectricEnergy_[irxn];
241  if (eamod != 0.0) {
242  kf[irxn] *= exp(-eamod/thermo(0).RT());
243  }
244  }
245  }
246 }
247 
249 {
251  // Loop over all reactions which are defined to have a voltage transfer
252  // coefficient that affects the activity energy for the reaction
253  for (size_t i = 0; i < m_ctrxn.size(); i++) {
254  size_t irxn = m_ctrxn[i];
255 
256  // Determine whether the reaction rate constant is in an exchange
257  // current density formulation format.
258  int iECDFormulation = m_ctrxn_ecdf[i];
259  if (iECDFormulation) {
260  // If the BV form is to be converted into the normal form then we go
261  // through this process. If it isn't to be converted, then we don't
262  // go through this process.
263  //
264  // We need to have the straight chemical reaction rate constant to
265  // come out of this calculation.
266  if (m_ctrxn_BVform[i] == 0) {
267  // Calculate the term and modify the forward reaction
268  double tmp = exp(- m_beta[i] * m_deltaG0[irxn] / thermo(0).RT());
269  tmp *= 1.0 / m_ProdStanConcReac[irxn] / Faraday;
270  kfwd[irxn] *= tmp;
271  }
272  // If BVform is nonzero we don't need to do anything.
273  } else {
274  // kfwd[] is the chemical reaction rate constant
275  //
276  // If we are to calculate the BV form directly, then we will do the
277  // reverse. We will calculate the exchange current density
278  // formulation here and substitute it.
279  if (m_ctrxn_BVform[i] != 0) {
280  // Calculate the term and modify the forward reaction rate
281  // constant so that it's in the exchange current density
282  // formulation format
283  double tmp = exp(m_beta[i] * m_deltaG0[irxn] * thermo(0).RT());
284  tmp *= Faraday * m_ProdStanConcReac[irxn];
285  kfwd[irxn] *= tmp;
286  }
287  }
288  }
289 }
290 
292 {
293  updateROP();
294 
295  // copy rate coefficients into kfwd
296  copy(m_rfn.begin(), m_rfn.end(), kfwd);
297 
298  // multiply by perturbation factor
299  multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin());
300 }
301 
302 void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible)
303 {
304  getFwdRateConstants(krev);
305  if (doIrreversible) {
307  for (size_t i = 0; i < nReactions(); i++) {
308  krev[i] /= m_ropnet[i];
309  }
310  } else {
311  multiply_each(krev, krev + nReactions(), m_rkcn.begin());
312  }
313 }
314 
316 {
317  // evaluate rate constants and equilibrium constants at temperature and phi
318  // (electric potential)
319  _update_rates_T();
320  // get updated activities (rates updated below)
321  _update_rates_C();
322 
323  if (m_ROP_ok) {
324  return;
325  }
326 
327  // Copy the reaction rate coefficients, m_rfn, into m_ropf
328  m_ropf = m_rfn;
329 
330  // Multiply by the perturbation factor
331  multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin());
332 
333  // Copy the forward rate constants to the reverse rate constants
334  m_ropr = m_ropf;
335 
336  // For reverse rates computed from thermochemistry, multiply
337  // the forward rates copied into m_ropr by the reciprocals of
338  // the equilibrium constants
339  multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin());
340 
341  // multiply ropf by the activity concentration reaction orders to obtain
342  // the forward rates of progress.
343  m_reactantStoich.multiply(m_actConc.data(), m_ropf.data());
344 
345  // For reversible reactions, multiply ropr by the activity concentration
346  // products
347  m_revProductStoich.multiply(m_actConc.data(), m_ropr.data());
348 
349  for (size_t j = 0; j != nReactions(); ++j) {
350  m_ropnet[j] = m_ropf[j] - m_ropr[j];
351  }
352 
353  // For reactions involving multiple phases, we must check that the phase
354  // being consumed actually exists. This is particularly important for phases
355  // that are stoichiometric phases containing one species with a unity
356  // activity
357  if (m_phaseExistsCheck) {
358  for (size_t j = 0; j != nReactions(); ++j) {
359  if ((m_ropr[j] > m_ropf[j]) && (m_ropr[j] > 0.0)) {
360  for (size_t p = 0; p < nPhases(); p++) {
361  if (m_rxnPhaseIsProduct[j][p] && !m_phaseExists[p]) {
362  m_ropnet[j] = 0.0;
363  m_ropr[j] = m_ropf[j];
364  if (m_ropf[j] > 0.0) {
365  for (size_t rp = 0; rp < nPhases(); rp++) {
366  if (m_rxnPhaseIsReactant[j][rp] && !m_phaseExists[rp]) {
367  m_ropnet[j] = 0.0;
368  m_ropr[j] = m_ropf[j] = 0.0;
369  }
370  }
371  }
372  }
373  if (m_rxnPhaseIsReactant[j][p] && !m_phaseIsStable[p]) {
374  m_ropnet[j] = 0.0;
375  m_ropr[j] = m_ropf[j];
376  }
377  }
378  } else if ((m_ropf[j] > m_ropr[j]) && (m_ropf[j] > 0.0)) {
379  for (size_t p = 0; p < nPhases(); p++) {
380  if (m_rxnPhaseIsReactant[j][p] && !m_phaseExists[p]) {
381  m_ropnet[j] = 0.0;
382  m_ropf[j] = m_ropr[j];
383  if (m_ropf[j] > 0.0) {
384  for (size_t rp = 0; rp < nPhases(); rp++) {
385  if (m_rxnPhaseIsProduct[j][rp] && !m_phaseExists[rp]) {
386  m_ropnet[j] = 0.0;
387  m_ropf[j] = m_ropr[j] = 0.0;
388  }
389  }
390  }
391  }
392  if (m_rxnPhaseIsProduct[j][p] && !m_phaseIsStable[p]) {
393  m_ropnet[j] = 0.0;
394  m_ropf[j] = m_ropr[j];
395  }
396  }
397  }
398  }
399  }
400  m_ROP_ok = true;
401 }
402 
403 void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
404 {
405  // Get the chemical potentials of the species in the all of the phases used
406  // in the kinetics mechanism
407  for (size_t n = 0; n < nPhases(); n++) {
408  m_thermo[n]->getChemPotentials(m_mu.data() + m_start[n]);
409  }
410 
411  // Use the stoichiometric manager to find deltaG for each reaction.
412  getReactionDelta(m_mu.data(), m_deltaG.data());
413  if (deltaG != 0 && (m_deltaG.data() != deltaG)) {
414  for (size_t j = 0; j < nReactions(); ++j) {
415  deltaG[j] = m_deltaG[j];
416  }
417  }
418 }
419 
421 {
422  // Get the chemical potentials of the species
423  for (size_t n = 0; n < nPhases(); n++) {
425  }
426 
427  // Use the stoichiometric manager to find deltaG for each reaction.
428  getReactionDelta(m_grt.data(), deltaM);
429 }
430 
431 void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH)
432 {
433  // Get the partial molar enthalpy of all species
434  for (size_t n = 0; n < nPhases(); n++) {
436  }
437 
438  // Use the stoichiometric manager to find deltaH for each reaction.
439  getReactionDelta(m_grt.data(), deltaH);
440 }
441 
442 void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS)
443 {
444  // Get the partial molar entropy of all species in all of the phases
445  for (size_t n = 0; n < nPhases(); n++) {
447  }
448 
449  // Use the stoichiometric manager to find deltaS for each reaction.
450  getReactionDelta(m_grt.data(), deltaS);
451 }
452 
453 void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS)
454 {
455  // Get the standard state chemical potentials of the species. This is the
456  // array of chemical potentials at unit activity We define these here as the
457  // chemical potentials of the pure species at the temperature and pressure
458  // of the solution.
459  for (size_t n = 0; n < nPhases(); n++) {
461  }
462 
463  // Use the stoichiometric manager to find deltaG for each reaction.
464  getReactionDelta(m_mu0.data(), deltaGSS);
465 }
466 
468 {
469  // Get the standard state enthalpies of the species. This is the array of
470  // chemical potentials at unit activity We define these here as the
471  // enthalpies of the pure species at the temperature and pressure of the
472  // solution.
473  for (size_t n = 0; n < nPhases(); n++) {
474  thermo(n).getEnthalpy_RT(m_grt.data() + m_start[n]);
475  }
476  for (size_t k = 0; k < m_kk; k++) {
477  m_grt[k] *= thermo(0).RT();
478  }
479 
480  // Use the stoichiometric manager to find deltaH for each reaction.
481  getReactionDelta(m_grt.data(), deltaH);
482 }
483 
484 void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS)
485 {
486  // Get the standard state entropy of the species. We define these here as
487  // the entropies of the pure species at the temperature and pressure of the
488  // solution.
489  for (size_t n = 0; n < nPhases(); n++) {
490  thermo(n).getEntropy_R(m_grt.data() + m_start[n]);
491  }
492  for (size_t k = 0; k < m_kk; k++) {
493  m_grt[k] *= GasConstant;
494  }
495 
496  // Use the stoichiometric manager to find deltaS for each reaction.
497  getReactionDelta(m_grt.data(), deltaS);
498 }
499 
500 bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base)
501 {
502  if (!m_surf) {
503  init();
504  }
505 
506  // Check that the number of surface sites is balanced
507  double reac_sites = 0.0;
508  double prod_sites = 0.0;
509  for (const auto& reactant : r_base->reactants) {
510  size_t k = m_surf->speciesIndex(reactant.first);
511  if (k != npos) {
512  reac_sites += reactant.second * m_surf->size(k);
513  }
514  }
515  for (const auto& product : r_base->products) {
516  size_t k = m_surf->speciesIndex(product.first);
517  if (k != npos) {
518  prod_sites += product.second * m_surf->size(k);
519  }
520  }
521  if (fabs(reac_sites - prod_sites) > 1e-5 * (reac_sites + prod_sites)) {
522  throw CanteraError("InterfaceKinetics::addReaction", "Number of surface"
523  " sites not balanced in reaction {}.\nReactant sites: {}\n"
524  "Product sites: {}", r_base->equation(), reac_sites, prod_sites);
525  }
526 
527  size_t i = nReactions();
528  bool added = Kinetics::addReaction(r_base);
529  if (!added) {
530  return false;
531  }
532 
533  InterfaceReaction& r = dynamic_cast<InterfaceReaction&>(*r_base);
534  SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, false);
535  m_rates.install(i, rate);
536 
537  // Turn on the global flag indicating surface coverage dependence
538  if (!r.coverage_deps.empty()) {
540  }
541 
542  ElectrochemicalReaction* re = dynamic_cast<ElectrochemicalReaction*>(&r);
543  if (re) {
544  m_has_electrochem_rxns = true;
545  m_beta.push_back(re->beta);
546  m_ctrxn.push_back(i);
547  if (re->exchange_current_density_formulation) {
549  m_ctrxn_ecdf.push_back(1);
550  } else {
551  m_ctrxn_ecdf.push_back(0);
552  }
553 
557  r.reaction_type == GLOBAL_RXN) {
558  // Specify alternative forms of the electrochemical reaction
559  if (r.reaction_type == BUTLERVOLMER_RXN) {
560  m_ctrxn_BVform.push_back(1);
562  m_ctrxn_BVform.push_back(2);
563  } else {
564  // set the default to be the normal forward / reverse calculation method
565  m_ctrxn_BVform.push_back(0);
566  }
567  if (!r.orders.empty()) {
568  vector_fp orders(nTotalSpecies(), 0.0);
569  for (const auto& order : r.orders) {
570  orders[kineticsSpeciesIndex(order.first)] = order.second;
571  }
572  }
573  } else {
574  m_ctrxn_BVform.push_back(0);
575  if (re->film_resistivity > 0.0) {
576  throw CanteraError("InterfaceKinetics::addReaction()",
577  "film resistivity set for elementary reaction");
578  }
579  }
580  }
581 
582  if (r.reversible) {
583  m_revindex.push_back(i);
584  } else {
585  m_irrev.push_back(i);
586  }
587 
588  m_rxnPhaseIsReactant.emplace_back(nPhases(), false);
589  m_rxnPhaseIsProduct.emplace_back(nPhases(), false);
590 
591  for (const auto& sp : r.reactants) {
592  size_t k = kineticsSpeciesIndex(sp.first);
593  size_t p = speciesPhaseIndex(k);
594  m_rxnPhaseIsReactant[i][p] = true;
595  }
596  for (const auto& sp : r.products) {
597  size_t k = kineticsSpeciesIndex(sp.first);
598  size_t p = speciesPhaseIndex(k);
599  m_rxnPhaseIsProduct[i][p] = true;
600  }
601 
602  deltaElectricEnergy_.push_back(0.0);
603  m_deltaG0.push_back(0.0);
604  m_deltaG.push_back(0.0);
605  m_ProdStanConcReac.push_back(0.0);
606 
607  return true;
608 }
609 
610 void InterfaceKinetics::modifyReaction(size_t i, shared_ptr<Reaction> r_base)
611 {
612  Kinetics::modifyReaction(i, r_base);
613  InterfaceReaction& r = dynamic_cast<InterfaceReaction&>(*r_base);
614  SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, true);
615  m_rates.replace(i, rate);
616 
617  // Invalidate cached data
618  m_redo_rates = true;
619  m_temp += 0.1;
620 }
621 
623  size_t i, InterfaceReaction& r, bool replace)
624 {
625  if (r.is_sticking_coefficient) {
626  // Identify the interface phase
627  size_t iInterface = npos;
628  size_t min_dim = 4;
629  for (size_t n = 0; n < nPhases(); n++) {
630  if (thermo(n).nDim() < min_dim) {
631  iInterface = n;
632  min_dim = thermo(n).nDim();
633  }
634  }
635 
636  std::string sticking_species = r.sticking_species;
637  if (sticking_species == "") {
638  // Identify the sticking species if not explicitly given
639  bool foundStick = false;
640  for (const auto& sp : r.reactants) {
641  size_t iPhase = speciesPhaseIndex(kineticsSpeciesIndex(sp.first));
642  if (iPhase != iInterface) {
643  // Non-interface species. There should be exactly one of these
644  if (foundStick) {
645  throw CanteraError("InterfaceKinetics::addReaction",
646  "Multiple non-interface species found"
647  "in sticking reaction: '" + r.equation() + "'");
648  }
649  foundStick = true;
650  sticking_species = sp.first;
651  }
652  }
653  if (!foundStick) {
654  throw CanteraError("InterfaceKinetics::addReaction",
655  "No non-interface species found"
656  "in sticking reaction: '" + r.equation() + "'");
657  }
658  }
659 
660  double surface_order = 0.0;
661  double multiplier = 1.0;
662  // Adjust the A-factor
663  for (const auto& sp : r.reactants) {
664  size_t iPhase = speciesPhaseIndex(kineticsSpeciesIndex(sp.first));
665  const ThermoPhase& p = thermo(iPhase);
666  size_t k = p.speciesIndex(sp.first);
667  if (sp.first == sticking_species) {
668  multiplier *= sqrt(GasConstant/(2*Pi*p.molecularWeight(k)));
669  } else {
670  // Non-sticking species. Convert from coverages used in the
671  // sticking probability expression to the concentration units
672  // used in the mass action rate expression. For surface phases,
673  // the dependence on the site density is incorporated when the
674  // rate constant is evaluated, since we don't assume that the
675  // site density is known at this time.
676  double order = getValue(r.orders, sp.first, sp.second);
677  if (&p == m_surf) {
678  multiplier *= pow(m_surf->size(k), order);
679  surface_order += order;
680  } else {
681  multiplier *= pow(p.standardConcentration(k), -order);
682  }
683  }
684  }
685 
686  if (!replace) {
687  m_stickingData.emplace_back(StickData{i, surface_order, multiplier,
689  } else {
690  // Modifying an existing sticking reaction.
691  for (auto& item : m_stickingData) {
692  if (item.index == i) {
693  item.order = surface_order;
694  item.multiplier = multiplier;
695  item.use_motz_wise = r.use_motz_wise_correction;
696  break;
697  }
698  }
699  }
700  }
701 
703  r.rate.temperatureExponent(),
704  r.rate.activationEnergy_R());
705 
706  // Set up coverage dependencies
707  for (const auto& sp : r.coverage_deps) {
708  size_t k = thermo(reactionPhaseIndex()).speciesIndex(sp.first);
709  rate.addCoverageDependence(k, sp.second.a, sp.second.m, sp.second.E);
710  }
711  return rate;
712 }
713 
714 void InterfaceKinetics::setIOFlag(int ioFlag)
715 {
716  m_ioFlag = ioFlag;
717  if (m_integrator) {
718  m_integrator->setIOFlag(ioFlag);
719  }
720 }
721 
723 {
725  m_phaseExists.push_back(true);
726  m_phaseIsStable.push_back(true);
727 }
728 
730 {
731  size_t ks = reactionPhaseIndex();
732  if (ks == npos) throw CanteraError("InterfaceKinetics::finalize",
733  "no surface phase is present.");
734 
735  // Check to see that the interface routine has a dimension of 2
736  m_surf = (SurfPhase*)&thermo(ks);
737  if (m_surf->nDim() != m_nDim) {
738  throw CanteraError("InterfaceKinetics::finalize",
739  "expected interface dimension = 2, but got dimension = {}",
740  m_surf->nDim());
741  }
742 }
743 
745 {
746  size_t kOld = m_kk;
748  if (m_kk != kOld && nReactions()) {
749  throw CanteraError("InterfaceKinetics::resizeSpecies", "Cannot add"
750  " species to InterfaceKinetics after reactions have been added.");
751  }
752  m_actConc.resize(m_kk);
753  m_conc.resize(m_kk);
754  m_StandardConc.resize(m_kk, 0.0);
755  m_mu0.resize(m_kk);
756  m_mu.resize(m_kk);
757  m_mu0_Kc.resize(m_kk);
758  m_grt.resize(m_kk);
759  m_pot.resize(m_kk, 0.0);
760  m_phi.resize(nPhases(), 0.0);
761 }
762 
763 doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const
764 {
765  for (size_t i = 0; i < m_ctrxn.size(); i++) {
766  if (m_ctrxn[i] == irxn) {
767  return m_beta[i];
768  }
769  }
770  return 0.0;
771 }
772 
774 {
775  if (m_integrator == 0) {
776  vector<InterfaceKinetics*> k{this};
779  }
780  m_integrator->integrate(0.0, tstep);
781  delete m_integrator;
782  m_integrator = 0;
783 }
784 
786  int ifuncOverride, doublereal timeScaleOverride)
787 {
788  // create our own solver object
789  if (m_integrator == 0) {
790  vector<InterfaceKinetics*> k{this};
793  }
794  m_integrator->setIOFlag(m_ioFlag);
795  // New direct method to go here
796  m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
797 }
798 
799 void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
800 {
801  if (iphase >= m_thermo.size()) {
802  throw CanteraError("InterfaceKinetics:setPhaseExistence", "out of bounds");
803  }
804  if (exists) {
805  if (!m_phaseExists[iphase]) {
807  m_phaseExistsCheck = std::max(m_phaseExistsCheck, 0);
808  m_phaseExists[iphase] = true;
809  }
810  m_phaseIsStable[iphase] = true;
811  } else {
812  if (m_phaseExists[iphase]) {
814  m_phaseExists[iphase] = false;
815  }
816  m_phaseIsStable[iphase] = false;
817  }
818 }
819 
820 int InterfaceKinetics::phaseExistence(const size_t iphase) const
821 {
822  if (iphase >= m_thermo.size()) {
823  throw CanteraError("InterfaceKinetics:phaseExistence()", "out of bounds");
824  }
825  return m_phaseExists[iphase];
826 }
827 
828 int InterfaceKinetics::phaseStability(const size_t iphase) const
829 {
830  if (iphase >= m_thermo.size()) {
831  throw CanteraError("InterfaceKinetics:phaseStability()", "out of bounds");
832  }
833  return m_phaseIsStable[iphase];
834 }
835 
836 void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable)
837 {
838  if (iphase >= m_thermo.size()) {
839  throw CanteraError("InterfaceKinetics:setPhaseStability", "out of bounds");
840  }
841  if (isStable) {
842  m_phaseIsStable[iphase] = true;
843  } else {
844  m_phaseIsStable[iphase] = false;
845  }
846 }
847 
848 void InterfaceKinetics::determineFwdOrdersBV(ElectrochemicalReaction& r, vector_fp& fwdFullOrders)
849 {
850  // Start out with the full ROP orders vector.
851  // This vector will have the BV exchange current density orders in it.
852  fwdFullOrders.assign(nTotalSpecies(), 0.0);
853  for (const auto& order : r.orders) {
854  fwdFullOrders[kineticsSpeciesIndex(order.first)] = order.second;
855  }
856 
857  // forward and reverse beta values
858  double betaf = r.beta;
859 
860  // Loop over the reactants doing away with the BV terms.
861  // This should leave the reactant terms only, even if they are non-mass action.
862  for (const auto& sp : r.reactants) {
863  size_t k = kineticsSpeciesIndex(sp.first);
864  fwdFullOrders[k] += betaf * sp.second;
865  // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet)
866  if (abs(fwdFullOrders[k]) < 0.00001) {
867  fwdFullOrders[k] = 0.0;
868  }
869  }
870 
871  // Loop over the products doing away with the BV terms.
872  // This should leave the reactant terms only, even if they are non-mass action.
873  for (const auto& sp : r.products) {
874  size_t k = kineticsSpeciesIndex(sp.first);
875  fwdFullOrders[k] -= betaf * sp.second;
876  // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet)
877  if (abs(fwdFullOrders[k]) < 0.00001) {
878  fwdFullOrders[k] = 0.0;
879  }
880  }
881 }
882 
883 void InterfaceKinetics::applyStickingCorrection(double T, double* kf)
884 {
885  if (m_stickingData.empty()) {
886  return;
887  }
888 
889  static const int cacheId = m_cache.getId();
890  CachedArray cached = m_cache.getArray(cacheId);
891  vector_fp& factors = cached.value;
892 
893  double n0 = m_surf->siteDensity();
894  if (!cached.validate(n0)) {
895  factors.resize(m_stickingData.size());
896  for (size_t n = 0; n < m_stickingData.size(); n++) {
897  factors[n] = pow(n0, -m_stickingData[n].order);
898  }
899  }
900 
901  for (size_t n = 0; n < m_stickingData.size(); n++) {
902  const StickData& item = m_stickingData[n];
903  if (item.use_motz_wise) {
904  kf[item.index] /= 1 - 0.5 * kf[item.index];
905  }
906  kf[item.index] *= factors[n] * sqrt(T) * item.multiplier;
907  }
908 }
909 
910 }
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Definition: Kinetics.cpp:572
virtual double size(size_t k) const
Returns the number of sites occupied by one molecule of species k.
Definition: SurfPhase.h:317
std::string sticking_species
For reactions with multiple non-surface species, the sticking species needs to be explicitly identifi...
Definition: Reaction.h:231
void applyVoltageKfwdCorrection(doublereal *const kfwd)
Apply modifications for the forward reaction rate for interfacial charge transfer reactions...
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
Definition: Kinetics.h:852
void setElectricPotential(int n, doublereal V)
Set the electric potential in the nth phase.
vector_fp m_deltaG0
Vector of delta G^0, the standard state Gibbs free energies for each reaction.
vector_fp m_ropr
Reverse rate-of-progress for each reaction.
Definition: Kinetics.h:921
std::vector< thermo_t * > m_thermo
m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator...
Definition: Kinetics.h:882
bool is_sticking_coefficient
Set to true if rate is a parameterization of the sticking coefficient rather than the forward rate co...
Definition: Reaction.h:221
int getId()
Get a unique id for a cached value.
Definition: ValueCache.cpp:21
int reaction_type
Type of the reaction.
Definition: Reaction.h:45
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase, assuming an ideal solution model (see Thermodynamic Properties and class SurfPhase).
virtual bool addReaction(shared_ptr< Reaction > r)
Add a single reaction to the mechanism.
doublereal electrochem_beta(size_t irxn) const
Return the charge transfer rxn Beta parameter for the ith reaction.
doublereal m_logtemp
Current log of the temperature.
virtual void getDeltaEnthalpy(doublereal *deltaH)
Return the vector of values for the reactions change in enthalpy.
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
size_t speciesIndex(const std::string &name) const
Returns the index of a species named &#39;name&#39; within the Phase object.
Definition: Phase.cpp:175
thermo_t & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism...
Definition: Kinetics.h:227
void setElectricPotential(doublereal v)
Set the electric potential of this phase (V).
Definition: ThermoPhase.h:294
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, doublereal timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
Definition: ThermoPhase.h:497
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
Definition: Kinetics.h:860
SurfPhase * m_surf
Pointer to the single surface phase.
std::vector< size_t > m_start
m_start is a vector of integers specifying the beginning position for the species vector for the n&#39;th...
Definition: Kinetics.h:888
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
vector_fp m_phi
Vector of phase electric potentials.
void advanceCoverages(doublereal tstep)
Advance the surface coverages in time.
vector_fp m_deltaG
Vector of deltaG[] of reaction, the delta Gibbs free energies for each reaction.
const int SURFACEAFFINITY_RXN
This is a surface reaction that is formulated using the affinity representation, common in the geoche...
Definition: reaction_defs.h:93
Composition orders
Forward reaction order with respect to specific species.
Definition: Reaction.h:56
virtual void getDeltaEntropy(doublereal *deltaS)
Return the vector of values for the reactions change in entropy.
vector_fp m_beta
Electrochemical transfer coefficient for the forward direction.
virtual void updateMu0()
Update the standard state chemical potentials and species equilibrium constant entries.
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
Definition: Kinetics.cpp:449
virtual bool addReaction(shared_ptr< Reaction > r)
Add a single reaction to the mechanism.
Definition: Kinetics.cpp:461
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized concentrations.
Definition: ThermoPhase.h:370
virtual void getDeltaSSEntropy(doublereal *deltaS)
Return the vector of values for the change in the standard state entropies for each reaction...
virtual void getEquilibriumConstants(doublereal *kc)
Equilibrium constant for all reactions including the voltage term.
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:266
STL namespace.
std::vector< std::vector< bool > > m_rxnPhaseIsProduct
Vector of vector of booleans indicating whether a phase participates in a reaction as a product...
vector_fp m_StandardConc
Vector of standard concentrations.
virtual void addPhase(thermo_t &thermo)
Add a phase to the kinetics manager object.
virtual void getRevReactionDelta(const doublereal *g, doublereal *dg)
Given an array of species properties &#39;g&#39;, return in array &#39;dg&#39; the change in this quantity in the rev...
Definition: Kinetics.cpp:383
vector_fp m_ropnet
Net rate-of-progress for each reaction.
Definition: Kinetics.h:924
const doublereal Pi
Pi.
Definition: ct_defs.h:51
size_t reactionPhaseIndex()
Phase where the reactions occur.
Definition: Kinetics.h:214
void convertExchangeCurrentDensityFormulation(doublereal *const kfwd)
When an electrode reaction rate is optionally specified in terms of its exchange current density...
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition: Kinetics.h:167
bool m_has_electrochem_rxns
Boolean flag indicating whether any reaction in the mechanism has a beta electrochemical parameter...
virtual doublereal standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
Definition: ThermoPhase.h:391
void multiply_each(OutputIter x_begin, OutputIter x_end, InputIter y_begin)
Multiply each entry in x by the corresponding entry in y.
Definition: utilities.h:162
int m_phaseExistsCheck
Int flag to indicate that some phases in the kinetics mechanism are non-existent. ...
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:748
vector_fp m_mu0_Kc
Vector of standard state electrochemical potentials modified by a standard concentration term...
vector_fp deltaElectricEnergy_
Storage for the net electric energy change due to reaction.
SurfaceArrhenius buildSurfaceArrhenius(size_t i, InterfaceReaction &r, bool replace)
Build a SurfaceArrhenius object from a Reaction, taking into account the possible sticking coefficien...
vector_fp m_ropf
Forward rate-of-progress for each reaction.
Definition: Kinetics.h:918
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
Values used for converting sticking coefficients into rate constants.
void updateKc()
Update the equilibrium constants and stored electrochemical potentials in molar units for all reversi...
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
doublereal m_temp
Current temperature of the data.
std::map< std::string, CoverageDependency > coverage_deps
Adjustments to the Arrhenius rate expression dependent on surface species coverages.
Definition: Reaction.h:217
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition: SurfPhase.h:142
size_t surfacePhaseIndex()
This returns the integer index of the phase which has ThermoPhase type cSurf.
Definition: Kinetics.h:202
bool m_has_coverage_dependence
Boolean flag indicating whether any reaction in the mechanism has a coverage dependent forward reacti...
ValueCache m_cache
Cache for saved calculations within each Kinetics object.
Definition: Kinetics.h:812
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
Definition: ThermoPhase.h:556
vector_int m_ctrxn_ecdf
Vector of booleans indicating whether the charge transfer reaction rate constant is described by an e...
virtual void getDeltaGibbs(doublereal *deltaG)
Return the vector of values for the reaction Gibbs free energy change.
void getElectrochemPotentials(doublereal *mu) const
Get the species electrochemical potentials.
Definition: ThermoPhase.cpp:86
bool reversible
True if the current reaction is reversible. False otherwise.
Definition: Reaction.h:63
std::vector< size_t > m_ctrxn_BVform
Vector of Reactions which follow the Butler-Volmer methodology for specifying the exchange current de...
void getCoverages(doublereal *theta) const
Return a vector of surface coverages.
Definition: SurfPhase.cpp:261
virtual void getFwdRateConstants(doublereal *kfwd)
Return the forward rate constants.
vector_int m_phaseIsStable
Vector of int indicating whether phases are stable or not.
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition: Kinetics.h:239
vector_fp m_conc
Array of concentrations for each species in the kinetics mechanism.
doublereal multiplier(size_t i) const
The current value of the multiplier for reaction i.
Definition: Kinetics.h:768
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, doublereal timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
std::vector< bool > m_phaseExists
Vector of booleans indicating whether phases exist or not.
vector_fp m_rfn
Forward rate constant for each reaction.
Definition: Kinetics.h:912
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
Definition: Kinetics.h:849
size_t m_nDim
Number of dimensions of reacting phase (2 for InterfaceKinetics, 1 for EdgeKinetics) ...
void _update_rates_T()
Update properties that depend on temperature.
A reaction occurring on an interface (i.e. a SurfPhase or an EdgePhase)
Definition: Reaction.h:206
vector_fp m_actConc
Array of activity concentrations for each species in the kinetics object.
vector_fp m_perturb
Vector of perturbation factors for each reaction&#39;s rate of progress vector.
Definition: Kinetics.h:864
doublereal electricPotential() const
Returns the electric potential of this phase (V).
Definition: ThermoPhase.h:302
double preExponentialFactor() const
Return the pre-exponential factor A (in m, kmol, s to powers depending on the reaction order) ...
Definition: RxnRates.h:76
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:261
virtual void getActivityConcentrations(doublereal *const conc)
Get the vector of activity concentrations used in the kinetics object.
void _update_rates_C()
Update properties that depend on the species mole fractions and/or concentration,.
vector_fp m_pot
Vector of potential energies due to Voltages.
const U & getValue(const std::map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a std::map.
Definition: utilities.h:504
int phaseExistence(const size_t iphase) const
Gets the phase existence int for the ith phase.
void setPhaseStability(const size_t iphase, const int isStable)
Set the stability of a phase in the reaction object.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
doublereal activationEnergy_R() const
Return the activation energy divided by the gas constant (i.e.
Definition: RxnRates.h:87
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition: Kinetics.h:135
Composition reactants
Reactant species and stoichiometric coefficients.
Definition: Reaction.h:48
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
Definition: ThermoPhase.h:487
double temperatureExponent() const
Return the temperature exponent b
Definition: RxnRates.h:81
bool m_has_exchange_current_density_formulation
Boolean flag indicating whether any reaction in the mechanism is described by an exchange current den...
virtual doublereal logStandardConc(size_t k=0) const
Natural logarithm of the standard concentration of the kth species.
Definition: ThermoPhase.cpp:65
Declarations for the implicit integration of surface site density equations (see Kinetics Managers an...
std::vector< std::vector< bool > > m_rxnPhaseIsReactant
Vector of vector of booleans indicating whether a phase participates in a reaction as a reactant...
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
std::vector< size_t > m_revindex
List of reactions numbers which are reversible reactions.
vector_fp m_rkcn
Reciprocal of the equilibrium constant in concentration units.
Definition: Kinetics.h:915
int phaseStability(const size_t iphase) const
Gets the phase stability int for the ith phase.
virtual void getRevRateConstants(doublereal *krev, bool doIrreversible=false)
Return the reverse rate constants.
std::vector< StickData > m_stickingData
Data for sticking reactions.
const int GLOBAL_RXN
A global reaction.
void _update_rates_phi()
Update properties that depend on the electric potential.
doublereal film_resistivity
Film Resistivity value.
Definition: Reaction.h:249
virtual void init()
Prepare the class for the addition of reactions, after all phases have been added.
std::vector< size_t > m_irrev
Vector of irreversible reaction numbers.
Advances the surface coverages of the associated set of SurfacePhase objects in time.
virtual void initialize(doublereal t0=0.0)
vector_fp m_mu
Vector of chemical potentials for all species.
void updateExchangeCurrentQuantities()
values needed to convert from exchange current density to surface reaction rate.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
virtual void updateROP()
Internal routine that updates the Rates of Progress of the reactions.
size_t speciesPhaseIndex(size_t k)
This function takes as an argument the kineticsSpecies index (i.e., the list index in the list of spe...
Definition: Kinetics.cpp:333
bool use_motz_wise_correction
Set to true if rate is a sticking coefficient which should be translated into a rate coefficient usin...
Definition: Reaction.h:227
doublereal beta
Forward value of the apparent Electrochemical transfer coefficient.
Definition: Reaction.h:252
Rate1< SurfaceArrhenius > m_rates
Templated class containing the vector of reactions for this interface.
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
std::vector< size_t > m_ctrxn
Vector of reaction indexes specifying the id of the charge transfer reactions in the mechanism...
vector_fp m_ProdStanConcReac
Vector of the products of the standard concentrations of the reactants.
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
Definition: Phase.h:585
const int BUTLERVOLMER_NOACTIVITYCOEFFS_RXN
This is a surface reaction that is formulated using the Butler-Volmer formulation and using concentra...
Definition: reaction_defs.h:81
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
Definition: ThermoPhase.h:546
void setPhaseExistence(const size_t iphase, const int exists)
Set the existence of a phase in the reaction object.
CachedArray getArray(int id)
Get a reference to a CachedValue object representing an array (vector_fp) with the given id...
Definition: ValueCache.h:171
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
Definition: ThermoPhase.h:566
const int BUTLERVOLMER_RXN
This is a surface reaction that is formulated using the Butler-Volmer formulation.
Definition: reaction_defs.h:87
Composition products
Product species and stoichiometric coefficients.
Definition: Reaction.h:51
doublereal molecularWeight(size_t k) const
Molecular weight of species k.
Definition: Phase.cpp:420
virtual void getDeltaSSGibbs(doublereal *deltaG)
Return the vector of values for the reaction standard state Gibbs free energy change.
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
void getConcentrations(doublereal *const c) const
Get the species concentrations (kmol/m^3).
Definition: Phase.cpp:519
virtual void getDeltaElectrochemPotentials(doublereal *deltaM)
Return the vector of values for the reaction electrochemical free energy change.
vector_fp m_mu0
Vector of standard state chemical potentials for all species.
vector_fp m_grt
Temporary work vector of length m_kk.
virtual void addPhase(thermo_t &thermo)
Add a phase to the kinetics manager object.
Definition: Kinetics.cpp:430
virtual void getDeltaSSEnthalpy(doublereal *deltaH)
Return the vector of values for the change in the standard state enthalpies of reaction.
ImplicitSurfChem * m_integrator
Pointer to the Implicit surface chemistry object.
An Arrhenius rate with coverage-dependent terms.
Definition: RxnRates.h:118
virtual void getReactionDelta(const doublereal *property, doublereal *deltaProperty)
Change in species properties.
Definition: Kinetics.cpp:373
An interface reaction which involves charged species.
Definition: Reaction.h:235
void integrate(doublereal t0, doublereal t1)
Integrate from t0 to t1. The integrator is reinitialized first.
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:577
doublereal siteDensity()
Returns the site density.
Definition: SurfPhase.h:312
std::string equation() const
The chemical equation for this reaction.
Definition: Reaction.cpp:89