Cantera  2.3.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 InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right)
45 {
46  // Call the assignment operator
47  operator=(right);
48 }
49 
50 InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
51 {
52  // Check for self assignment.
53  if (this == &right) {
54  return *this;
55  }
56 
57  Kinetics::operator=(right);
58 
59  m_grt = right.m_grt;
60  m_revindex = right.m_revindex;
61  m_rates = right.m_rates;
62  m_redo_rates = right.m_redo_rates;
63  m_irrev = right.m_irrev;
64  m_conc = right.m_conc;
65  m_actConc = right.m_actConc;
66  m_mu0 = right.m_mu0;
67  m_mu = right.m_mu;
68  m_mu0_Kc = right.m_mu0_Kc;
69  m_phi = right.m_phi;
70  m_pot = right.m_pot;
71  deltaElectricEnergy_ = right.deltaElectricEnergy_;
72  m_surf = right.m_surf; //DANGER - shallow copy
73  m_integrator = right.m_integrator; //DANGER - shallow copy
74  m_beta = right.m_beta;
75  m_ctrxn = right.m_ctrxn;
76  m_ctrxn_BVform = right.m_ctrxn_BVform;
77  m_ctrxn_ecdf = right.m_ctrxn_ecdf;
78  m_StandardConc = right.m_StandardConc;
79  m_deltaG0 = right.m_deltaG0;
80  m_deltaG = right.m_deltaG;
81  m_ProdStanConcReac = right.m_ProdStanConcReac;
82  m_ROP_ok = right.m_ROP_ok;
83  m_temp = right.m_temp;
84  m_logtemp = right.m_logtemp;
85  m_has_coverage_dependence = right.m_has_coverage_dependence;
86  m_has_electrochem_rxns = right.m_has_electrochem_rxns;
87  m_has_exchange_current_density_formulation = right.m_has_exchange_current_density_formulation;
88  m_phaseExistsCheck = right.m_phaseExistsCheck;
89  m_phaseExists = right.m_phaseExists;
90  m_phaseIsStable = right.m_phaseIsStable;
91  m_rxnPhaseIsReactant = right.m_rxnPhaseIsReactant;
92  m_rxnPhaseIsProduct = right.m_rxnPhaseIsProduct;
93  m_ioFlag = right.m_ioFlag;
94 
95  return *this;
96 }
97 
99 {
100  warn_deprecated("InterfaceKinetics::type",
101  "To be removed after Cantera 2.3.");
102  return cInterfaceKinetics;
103 }
104 
105 Kinetics* InterfaceKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
106 {
107  InterfaceKinetics* iK = new InterfaceKinetics(*this);
108  iK->assignShallowPointers(tpVector);
109  return iK;
110 }
111 
113 {
115  m_redo_rates = true;
116 }
117 
119 {
120  // First task is update the electrical potentials from the Phases
123  m_surf->getCoverages(m_actConc.data());
124  m_rates.update_C(m_actConc.data());
125  m_redo_rates = true;
126  }
127 
128  // Go find the temperature from the surface
129  doublereal T = thermo(surfacePhaseIndex()).temperature();
130  m_redo_rates = true;
131  if (T != m_temp || m_redo_rates) {
132  m_logtemp = log(T);
133 
134  // Calculate the forward rate constant by calling m_rates and store it in m_rfn[]
135  m_rates.update(T, m_logtemp, m_rfn.data());
136  applyStickingCorrection(T, m_rfn.data());
137 
138  // If we need to do conversions between exchange current density
139  // formulation and regular formulation (either way) do it here.
142  }
145  }
146  m_temp = T;
147  updateKc();
148  m_ROP_ok = false;
149  m_redo_rates = false;
150  }
151 }
152 
154 {
155  // Store electric potentials for each phase in the array m_phi[].
156  for (size_t n = 0; n < nPhases(); n++) {
157  if (thermo(n).electricPotential() != m_phi[n]) {
158  m_phi[n] = thermo(n).electricPotential();
159  m_redo_rates = true;
160  }
161  }
162 }
163 
165 {
166  for (size_t n = 0; n < nPhases(); n++) {
167  const ThermoPhase* tp = m_thermo[n];
168  /*
169  * We call the getActivityConcentrations function of each ThermoPhase
170  * class that makes up this kinetics object to obtain the generalized
171  * concentrations for species within that class. This is collected in
172  * the vector m_conc. m_start[] are integer indices for that vector
173  * denoting the start of the species for each phase.
174  */
176 
177  // Get regular concentrations too
178  tp->getConcentrations(m_conc.data() + m_start[n]);
179  }
180  m_ROP_ok = false;
181 }
182 
184 {
185  _update_rates_C();
186  copy(m_actConc.begin(), m_actConc.end(), conc);
187 }
188 
190 {
191  fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
192 
193  if (m_revindex.size() > 0) {
194  /*
195  * Get the vector of standard state electrochemical potentials for
196  * species in the Interfacial kinetics object and store it in m_mu0[]
197  * and m_mu0_Kc[]
198  */
199  updateMu0();
200  doublereal rrt = 1.0 / thermo(0).RT();
201 
202  // compute Delta mu^0 for all reversible reactions
203  getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data());
204 
205  for (size_t i = 0; i < m_revindex.size(); i++) {
206  size_t irxn = m_revindex[i];
207  if (irxn == npos || irxn >= nReactions()) {
208  throw CanteraError("InterfaceKinetics", "illegal value: irxn = {}", irxn);
209  }
210  // WARNING this may overflow HKM
211  m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);
212  }
213  for (size_t i = 0; i != m_irrev.size(); ++i) {
214  m_rkcn[ m_irrev[i] ] = 0.0;
215  }
216  }
217 }
218 
220 {
221  // First task is update the electrical potentials from the Phases
223 
225  size_t ik = 0;
226  for (size_t n = 0; n < nPhases(); n++) {
228  for (size_t k = 0; k < thermo(n).nSpecies(); k++) {
229  m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k);
230  m_mu0_Kc[ik] -= thermo(0).RT() * thermo(n).logStandardConc(k);
231  ik++;
232  }
233  }
234 }
235 
237 {
238  warn_deprecated("InterfaceKinetics::checkPartialEquil",
239  "To be removed after Cantera 2.3.");
240  // First task is update the electrical potentials from the Phases
242 
243  vector_fp dmu(nTotalSpecies(), 0.0);
244  vector_fp rmu(std::max<size_t>(nReactions(), 1), 0.0);
245  if (m_revindex.size() > 0) {
246  cout << "T = " << thermo(0).temperature() << " " << thermo(0).RT() << endl;
247  size_t nsp, ik=0;
248  doublereal delta;
249  for (size_t n = 0; n < nPhases(); n++) {
250  thermo(n).getChemPotentials(dmu.data() + m_start[n]);
251  nsp = thermo(n).nSpecies();
252  for (size_t k = 0; k < nsp; k++) {
253  delta = Faraday * m_phi[n] * thermo(n).charge(k);
254  dmu[ik] += delta;
255  ik++;
256  }
257  }
258 
259  // compute Delta mu^ for all reversible reactions
260  getRevReactionDelta(dmu.data(), rmu.data());
261  updateROP();
262  for (size_t i = 0; i < m_revindex.size(); i++) {
263  size_t irxn = m_revindex[i];
264  writelog("Reaction {} {}\n",
265  reactionString(irxn), rmu[irxn]/thermo(0).RT());
266  writelogf("%12.6e %12.6e %12.6e %12.6e \n",
267  m_ropf[irxn], m_ropr[irxn], m_ropnet[irxn],
268  m_ropnet[irxn]/(m_ropf[irxn] + m_ropr[irxn]));
269  }
270  }
271 }
272 
274 {
275  updateMu0();
276  doublereal rrt = 1.0 / thermo(0).RT();
277  std::fill(kc, kc + nReactions(), 0.0);
278  getReactionDelta(m_mu0_Kc.data(), kc);
279  for (size_t i = 0; i < nReactions(); i++) {
280  kc[i] = exp(-kc[i]*rrt);
281  }
282 }
283 
285 {
286  // Calculate:
287  // - m_StandardConc[]
288  // - m_ProdStanConcReac[]
289  // - m_deltaG0[]
290  // - m_mu0[]
291 
292  // First collect vectors of the standard Gibbs free energies of the
293  // species and the standard concentrations
294  // - m_mu0
295  // - m_StandardConc
296  size_t ik = 0;
297 
298  for (size_t n = 0; n < nPhases(); n++) {
300  size_t nsp = thermo(n).nSpecies();
301  for (size_t k = 0; k < nsp; k++) {
303  ik++;
304  }
305  }
306 
307  getReactionDelta(m_mu0.data(), m_deltaG0.data());
308 
309  // Calculate the product of the standard concentrations of the reactants
310  for (size_t i = 0; i < nReactions(); i++) {
311  m_ProdStanConcReac[i] = 1.0;
312  }
313  m_reactantStoich.multiply(m_StandardConc.data(), m_ProdStanConcReac.data());
314 }
315 
317 {
318  // Compute the electrical potential energy of each species
319  size_t ik = 0;
320  for (size_t n = 0; n < nPhases(); n++) {
321  size_t nsp = thermo(n).nSpecies();
322  for (size_t k = 0; k < nsp; k++) {
323  m_pot[ik] = Faraday * thermo(n).charge(k) * m_phi[n];
324  ik++;
325  }
326  }
327 
328  // Compute the change in electrical potential energy for each reaction. This
329  // will only be non-zero if a potential difference is present.
331 
332  // Modify the reaction rates. Only modify those with a non-zero activation
333  // energy. Below we decrease the activation energy below zero but in some
334  // debug modes we print out a warning message about this.
335 
336  // NOTE, there is some discussion about this point. Should we decrease the
337  // activation energy below zero? I don't think this has been decided in any
338  // definitive way. The treatment below is numerically more stable, however.
339  for (size_t i = 0; i < m_beta.size(); i++) {
340  size_t irxn = m_ctrxn[i];
341 
342  // If we calculate the BV form directly, we don't add the voltage
343  // correction to the forward reaction rate constants.
344  if (m_ctrxn_BVform[i] == 0) {
345  double eamod = m_beta[i] * deltaElectricEnergy_[irxn];
346  if (eamod != 0.0) {
347  kf[irxn] *= exp(-eamod/thermo(0).RT());
348  }
349  }
350  }
351 }
352 
354 {
356  // Loop over all reactions which are defined to have a voltage transfer
357  // coefficient that affects the activity energy for the reaction
358  for (size_t i = 0; i < m_ctrxn.size(); i++) {
359  size_t irxn = m_ctrxn[i];
360 
361  // Determine whether the reaction rate constant is in an exchange
362  // current density formulation format.
363  int iECDFormulation = m_ctrxn_ecdf[i];
364  if (iECDFormulation) {
365  // If the BV form is to be converted into the normal form then we go
366  // through this process. If it isn't to be converted, then we don't
367  // go through this process.
368  //
369  // We need to have the straight chemical reaction rate constant to
370  // come out of this calculation.
371  if (m_ctrxn_BVform[i] == 0) {
372  // Calculate the term and modify the forward reaction
373  double tmp = exp(- m_beta[i] * m_deltaG0[irxn] / thermo(0).RT());
374  tmp *= 1.0 / m_ProdStanConcReac[irxn] / Faraday;
375  kfwd[irxn] *= tmp;
376  }
377  // If BVform is nonzero we don't need to do anything.
378  } else {
379  // kfwd[] is the chemical reaction rate constant
380  //
381  // If we are to calculate the BV form directly, then we will do the
382  // reverse. We will calculate the exchange current density
383  // formulation here and substitute it.
384  if (m_ctrxn_BVform[i] != 0) {
385  // Calculate the term and modify the forward reaction rate
386  // constant so that it's in the exchange current density
387  // formulation format
388  double tmp = exp(m_beta[i] * m_deltaG0[irxn] * thermo(0).RT());
389  tmp *= Faraday * m_ProdStanConcReac[irxn];
390  kfwd[irxn] *= tmp;
391  }
392  }
393  }
394 }
395 
397 {
398  updateROP();
399 
400  // copy rate coefficients into kfwd
401  copy(m_rfn.begin(), m_rfn.end(), kfwd);
402 
403  // multiply by perturbation factor
404  multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin());
405 }
406 
407 void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible)
408 {
409  getFwdRateConstants(krev);
410  if (doIrreversible) {
412  for (size_t i = 0; i < nReactions(); i++) {
413  krev[i] /= m_ropnet[i];
414  }
415  } else {
416  multiply_each(krev, krev + nReactions(), m_rkcn.begin());
417  }
418 }
419 
421 {
422  // evaluate rate constants and equilibrium constants at temperature and phi
423  // (electric potential)
424  _update_rates_T();
425  // get updated activities (rates updated below)
426  _update_rates_C();
427 
428  if (m_ROP_ok) {
429  return;
430  }
431 
432  // Copy the reaction rate coefficients, m_rfn, into m_ropf
433  m_ropf = m_rfn;
434 
435  // Multiply by the perturbation factor
436  multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin());
437 
438  // Copy the forward rate constants to the reverse rate constants
439  m_ropr = m_ropf;
440 
441  // For reverse rates computed from thermochemistry, multiply
442  // the forward rates copied into m_ropr by the reciprocals of
443  // the equilibrium constants
444  multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin());
445 
446  // multiply ropf by the activity concentration reaction orders to obtain
447  // the forward rates of progress.
448  m_reactantStoich.multiply(m_actConc.data(), m_ropf.data());
449 
450  // For reversible reactions, multiply ropr by the activity concentration
451  // products
452  m_revProductStoich.multiply(m_actConc.data(), m_ropr.data());
453 
454  for (size_t j = 0; j != nReactions(); ++j) {
455  m_ropnet[j] = m_ropf[j] - m_ropr[j];
456  }
457 
458  // For reactions involving multiple phases, we must check that the phase
459  // being consumed actually exists. This is particularly important for phases
460  // that are stoichiometric phases containing one species with a unity
461  // activity
462  if (m_phaseExistsCheck) {
463  for (size_t j = 0; j != nReactions(); ++j) {
464  if ((m_ropr[j] > m_ropf[j]) && (m_ropr[j] > 0.0)) {
465  for (size_t p = 0; p < nPhases(); p++) {
466  if (m_rxnPhaseIsProduct[j][p] && !m_phaseExists[p]) {
467  m_ropnet[j] = 0.0;
468  m_ropr[j] = m_ropf[j];
469  if (m_ropf[j] > 0.0) {
470  for (size_t rp = 0; rp < nPhases(); rp++) {
471  if (m_rxnPhaseIsReactant[j][rp] && !m_phaseExists[rp]) {
472  m_ropnet[j] = 0.0;
473  m_ropr[j] = m_ropf[j] = 0.0;
474  }
475  }
476  }
477  }
478  if (m_rxnPhaseIsReactant[j][p] && !m_phaseIsStable[p]) {
479  m_ropnet[j] = 0.0;
480  m_ropr[j] = m_ropf[j];
481  }
482  }
483  } else if ((m_ropf[j] > m_ropr[j]) && (m_ropf[j] > 0.0)) {
484  for (size_t p = 0; p < nPhases(); p++) {
485  if (m_rxnPhaseIsReactant[j][p] && !m_phaseExists[p]) {
486  m_ropnet[j] = 0.0;
487  m_ropf[j] = m_ropr[j];
488  if (m_ropf[j] > 0.0) {
489  for (size_t rp = 0; rp < nPhases(); rp++) {
490  if (m_rxnPhaseIsProduct[j][rp] && !m_phaseExists[rp]) {
491  m_ropnet[j] = 0.0;
492  m_ropf[j] = m_ropr[j] = 0.0;
493  }
494  }
495  }
496  }
497  if (m_rxnPhaseIsProduct[j][p] && !m_phaseIsStable[p]) {
498  m_ropnet[j] = 0.0;
499  m_ropf[j] = m_ropr[j];
500  }
501  }
502  }
503  }
504  }
505  m_ROP_ok = true;
506 }
507 
508 void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
509 {
510  // Get the chemical potentials of the species in the all of the phases used
511  // in the kinetics mechanism
512  for (size_t n = 0; n < nPhases(); n++) {
513  m_thermo[n]->getChemPotentials(m_mu.data() + m_start[n]);
514  }
515 
516  // Use the stoichiometric manager to find deltaG for each reaction.
517  getReactionDelta(m_mu.data(), m_deltaG.data());
518  if (deltaG != 0 && (m_deltaG.data() != deltaG)) {
519  for (size_t j = 0; j < nReactions(); ++j) {
520  deltaG[j] = m_deltaG[j];
521  }
522  }
523 }
524 
526 {
527  // Get the chemical potentials of the species
528  for (size_t n = 0; n < nPhases(); n++) {
530  }
531 
532  // Use the stoichiometric manager to find deltaG for each reaction.
533  getReactionDelta(m_grt.data(), deltaM);
534 }
535 
536 void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH)
537 {
538  // Get the partial molar enthalpy of all species
539  for (size_t n = 0; n < nPhases(); n++) {
541  }
542 
543  // Use the stoichiometric manager to find deltaH for each reaction.
544  getReactionDelta(m_grt.data(), deltaH);
545 }
546 
547 void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS)
548 {
549  // Get the partial molar entropy of all species in all of the phases
550  for (size_t n = 0; n < nPhases(); n++) {
552  }
553 
554  // Use the stoichiometric manager to find deltaS for each reaction.
555  getReactionDelta(m_grt.data(), deltaS);
556 }
557 
558 void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS)
559 {
560  // Get the standard state chemical potentials of the species. This is the
561  // array of chemical potentials at unit activity We define these here as the
562  // chemical potentials of the pure species at the temperature and pressure
563  // of the solution.
564  for (size_t n = 0; n < nPhases(); n++) {
566  }
567 
568  // Use the stoichiometric manager to find deltaG for each reaction.
569  getReactionDelta(m_mu0.data(), deltaGSS);
570 }
571 
573 {
574  // Get the standard state enthalpies of the species. This is the array of
575  // chemical potentials at unit activity We define these here as the
576  // enthalpies of the pure species at the temperature and pressure of the
577  // solution.
578  for (size_t n = 0; n < nPhases(); n++) {
579  thermo(n).getEnthalpy_RT(m_grt.data() + m_start[n]);
580  }
581  for (size_t k = 0; k < m_kk; k++) {
582  m_grt[k] *= thermo(0).RT();
583  }
584 
585  // Use the stoichiometric manager to find deltaH for each reaction.
586  getReactionDelta(m_grt.data(), deltaH);
587 }
588 
589 void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS)
590 {
591  // Get the standard state entropy of the species. We define these here as
592  // the entropies of the pure species at the temperature and pressure of the
593  // solution.
594  for (size_t n = 0; n < nPhases(); n++) {
595  thermo(n).getEntropy_R(m_grt.data() + m_start[n]);
596  }
597  for (size_t k = 0; k < m_kk; k++) {
598  m_grt[k] *= GasConstant;
599  }
600 
601  // Use the stoichiometric manager to find deltaS for each reaction.
602  getReactionDelta(m_grt.data(), deltaS);
603 }
604 
605 bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base)
606 {
607  size_t i = nReactions();
608  bool added = Kinetics::addReaction(r_base);
609  if (!added) {
610  return false;
611  }
612 
613  InterfaceReaction& r = dynamic_cast<InterfaceReaction&>(*r_base);
614  SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, false);
615  m_rates.install(i, rate);
616 
617  // Turn on the global flag indicating surface coverage dependence
618  if (!r.coverage_deps.empty()) {
620  }
621 
622  ElectrochemicalReaction* re = dynamic_cast<ElectrochemicalReaction*>(&r);
623  if (re) {
624  m_has_electrochem_rxns = true;
625  m_beta.push_back(re->beta);
626  m_ctrxn.push_back(i);
627  if (re->exchange_current_density_formulation) {
629  m_ctrxn_ecdf.push_back(1);
630  } else {
631  m_ctrxn_ecdf.push_back(0);
632  }
633 
637  r.reaction_type == GLOBAL_RXN) {
638  // Specify alternative forms of the electrochemical reaction
639  if (r.reaction_type == BUTLERVOLMER_RXN) {
640  m_ctrxn_BVform.push_back(1);
642  m_ctrxn_BVform.push_back(2);
643  } else {
644  // set the default to be the normal forward / reverse calculation method
645  m_ctrxn_BVform.push_back(0);
646  }
647  if (!r.orders.empty()) {
648  vector_fp orders(nTotalSpecies(), 0.0);
649  for (const auto& order : r.orders) {
650  orders[kineticsSpeciesIndex(order.first)] = order.second;
651  }
652  }
653  } else {
654  m_ctrxn_BVform.push_back(0);
655  if (re->film_resistivity > 0.0) {
656  throw CanteraError("InterfaceKinetics::addReaction()",
657  "film resistivity set for elementary reaction");
658  }
659  }
660  }
661 
662  if (r.reversible) {
663  m_revindex.push_back(i);
664  } else {
665  m_irrev.push_back(i);
666  }
667 
668  m_rxnPhaseIsReactant.emplace_back(nPhases(), false);
669  m_rxnPhaseIsProduct.emplace_back(nPhases(), false);
670 
671  for (const auto& sp : r.reactants) {
672  size_t k = kineticsSpeciesIndex(sp.first);
673  size_t p = speciesPhaseIndex(k);
674  m_rxnPhaseIsReactant[i][p] = true;
675  }
676  for (const auto& sp : r.products) {
677  size_t k = kineticsSpeciesIndex(sp.first);
678  size_t p = speciesPhaseIndex(k);
679  m_rxnPhaseIsProduct[i][p] = true;
680  }
681 
682  deltaElectricEnergy_.push_back(0.0);
683  m_deltaG0.push_back(0.0);
684  m_deltaG.push_back(0.0);
685  m_ProdStanConcReac.push_back(0.0);
686 
687  return true;
688 }
689 
690 void InterfaceKinetics::modifyReaction(size_t i, shared_ptr<Reaction> r_base)
691 {
692  Kinetics::modifyReaction(i, r_base);
693  InterfaceReaction& r = dynamic_cast<InterfaceReaction&>(*r_base);
694  SurfaceArrhenius rate = buildSurfaceArrhenius(i, r, true);
695  m_rates.replace(i, rate);
696 
697  // Invalidate cached data
698  m_redo_rates = true;
699  m_temp += 0.1;
700 }
701 
703  size_t i, InterfaceReaction& r, bool replace)
704 {
705  if (r.is_sticking_coefficient) {
706  // Identify the interface phase
707  size_t iInterface = npos;
708  size_t min_dim = 4;
709  for (size_t n = 0; n < nPhases(); n++) {
710  if (thermo(n).nDim() < min_dim) {
711  iInterface = n;
712  min_dim = thermo(n).nDim();
713  }
714  }
715 
716  std::string sticking_species = r.sticking_species;
717  if (sticking_species == "") {
718  // Identify the sticking species if not explicitly given
719  bool foundStick = false;
720  for (const auto& sp : r.reactants) {
721  size_t iPhase = speciesPhaseIndex(kineticsSpeciesIndex(sp.first));
722  if (iPhase != iInterface) {
723  // Non-interface species. There should be exactly one of these
724  if (foundStick) {
725  throw CanteraError("InterfaceKinetics::addReaction",
726  "Multiple non-interface species found"
727  "in sticking reaction: '" + r.equation() + "'");
728  }
729  foundStick = true;
730  sticking_species = sp.first;
731  }
732  }
733  if (!foundStick) {
734  throw CanteraError("InterfaceKinetics::addReaction",
735  "No non-interface species found"
736  "in sticking reaction: '" + r.equation() + "'");
737  }
738  }
739 
740  double surface_order = 0.0;
741  double multiplier = 1.0;
742  // Adjust the A-factor
743  for (const auto& sp : r.reactants) {
744  size_t iPhase = speciesPhaseIndex(kineticsSpeciesIndex(sp.first));
745  const ThermoPhase& p = thermo(iPhase);
746  const ThermoPhase& surf = thermo(surfacePhaseIndex());
747  size_t k = p.speciesIndex(sp.first);
748  if (sp.first == sticking_species) {
749  multiplier *= sqrt(GasConstant/(2*Pi*p.molecularWeight(k)));
750  } else {
751  // Non-sticking species. Convert from coverages used in the
752  // sticking probability expression to the concentration units
753  // used in the mass action rate expression. For surface phases,
754  // the dependence on the site density is incorporated when the
755  // rate constant is evaluated, since we don't assume that the
756  // site density is known at this time.
757  double order = getValue(r.orders, sp.first, sp.second);
758  if (&p == &surf) {
759  multiplier *= pow(p.size(k), order);
760  surface_order += order;
761  } else {
762  multiplier *= pow(p.standardConcentration(k), -order);
763  }
764  }
765  }
766 
767  if (!replace) {
768  m_stickingData.emplace_back(StickData{i, surface_order, multiplier,
770  } else {
771  // Modifying an existing sticking reaction.
772  for (auto& item : m_stickingData) {
773  if (item.index == i) {
774  item.order = surface_order;
775  item.multiplier = multiplier;
776  item.use_motz_wise = r.use_motz_wise_correction;
777  break;
778  }
779  }
780  }
781  }
782 
784  r.rate.temperatureExponent(),
785  r.rate.activationEnergy_R());
786 
787  // Set up coverage dependencies
788  for (const auto& sp : r.coverage_deps) {
789  size_t k = thermo(reactionPhaseIndex()).speciesIndex(sp.first);
790  rate.addCoverageDependence(k, sp.second.a, sp.second.m, sp.second.E);
791  }
792  return rate;
793 }
794 
795 void InterfaceKinetics::setIOFlag(int ioFlag)
796 {
797  m_ioFlag = ioFlag;
798  if (m_integrator) {
799  m_integrator->setIOFlag(ioFlag);
800  }
801 }
802 
804 {
806  m_phaseExists.push_back(true);
807  m_phaseIsStable.push_back(true);
808 }
809 
811 {
812  size_t ks = reactionPhaseIndex();
813  if (ks == npos) throw CanteraError("InterfaceKinetics::finalize",
814  "no surface phase is present.");
815 
816  // Check to see that the interface routine has a dimension of 2
817  m_surf = (SurfPhase*)&thermo(ks);
818  if (m_surf->nDim() != m_nDim) {
819  throw CanteraError("InterfaceKinetics::finalize",
820  "expected interface dimension = 2, but got dimension = {}",
821  m_surf->nDim());
822  }
823 }
824 
826 {
827  size_t kOld = m_kk;
829  if (m_kk != kOld && nReactions()) {
830  throw CanteraError("InterfaceKinetics::resizeSpecies", "Cannot add"
831  " species to InterfaceKinetics after reactions have been added.");
832  }
833  m_actConc.resize(m_kk);
834  m_conc.resize(m_kk);
835  m_StandardConc.resize(m_kk, 0.0);
836  m_mu0.resize(m_kk);
837  m_mu.resize(m_kk);
838  m_mu0_Kc.resize(m_kk);
839  m_grt.resize(m_kk);
840  m_pot.resize(m_kk, 0.0);
841  m_phi.resize(nPhases(), 0.0);
842 }
843 
844 doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const
845 {
846  for (size_t i = 0; i < m_ctrxn.size(); i++) {
847  if (m_ctrxn[i] == irxn) {
848  return m_beta[i];
849  }
850  }
851  return 0.0;
852 }
853 
855 {
856  if (m_integrator == 0) {
857  vector<InterfaceKinetics*> k{this};
860  }
861  m_integrator->integrate(0.0, tstep);
862  delete m_integrator;
863  m_integrator = 0;
864 }
865 
867  int ifuncOverride, doublereal timeScaleOverride)
868 {
869  // create our own solver object
870  if (m_integrator == 0) {
871  vector<InterfaceKinetics*> k{this};
874  }
875  m_integrator->setIOFlag(m_ioFlag);
876  // New direct method to go here
877  m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
878 }
879 
880 void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
881 {
882  if (iphase >= m_thermo.size()) {
883  throw CanteraError("InterfaceKinetics:setPhaseExistence", "out of bounds");
884  }
885  if (exists) {
886  if (!m_phaseExists[iphase]) {
888  m_phaseExistsCheck = std::max(m_phaseExistsCheck, 0);
889  m_phaseExists[iphase] = true;
890  }
891  m_phaseIsStable[iphase] = true;
892  } else {
893  if (m_phaseExists[iphase]) {
895  m_phaseExists[iphase] = false;
896  }
897  m_phaseIsStable[iphase] = false;
898  }
899 }
900 
901 int InterfaceKinetics::phaseExistence(const size_t iphase) const
902 {
903  if (iphase >= m_thermo.size()) {
904  throw CanteraError("InterfaceKinetics:phaseExistence()", "out of bounds");
905  }
906  return m_phaseExists[iphase];
907 }
908 
909 int InterfaceKinetics::phaseStability(const size_t iphase) const
910 {
911  if (iphase >= m_thermo.size()) {
912  throw CanteraError("InterfaceKinetics:phaseStability()", "out of bounds");
913  }
914  return m_phaseIsStable[iphase];
915 }
916 
917 void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable)
918 {
919  if (iphase >= m_thermo.size()) {
920  throw CanteraError("InterfaceKinetics:setPhaseStability", "out of bounds");
921  }
922  if (isStable) {
923  m_phaseIsStable[iphase] = true;
924  } else {
925  m_phaseIsStable[iphase] = false;
926  }
927 }
928 
929 void InterfaceKinetics::determineFwdOrdersBV(ElectrochemicalReaction& r, vector_fp& fwdFullOrders)
930 {
931  // Start out with the full ROP orders vector.
932  // This vector will have the BV exchange current density orders in it.
933  fwdFullOrders.assign(nTotalSpecies(), 0.0);
934  for (const auto& order : r.orders) {
935  fwdFullOrders[kineticsSpeciesIndex(order.first)] = order.second;
936  }
937 
938  // forward and reverse beta values
939  double betaf = r.beta;
940 
941  // Loop over the reactants doing away with the BV terms.
942  // This should leave the reactant terms only, even if they are non-mass action.
943  for (const auto& sp : r.reactants) {
944  size_t k = kineticsSpeciesIndex(sp.first);
945  fwdFullOrders[k] += betaf * sp.second;
946  // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet)
947  if (abs(fwdFullOrders[k]) < 0.00001) {
948  fwdFullOrders[k] = 0.0;
949  }
950  }
951 
952  // Loop over the products doing away with the BV terms.
953  // This should leave the reactant terms only, even if they are non-mass action.
954  for (const auto& sp : r.products) {
955  size_t k = kineticsSpeciesIndex(sp.first);
956  fwdFullOrders[k] -= betaf * sp.second;
957  // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet)
958  if (abs(fwdFullOrders[k]) < 0.00001) {
959  fwdFullOrders[k] = 0.0;
960  }
961  }
962 }
963 
964 void InterfaceKinetics::applyStickingCorrection(double T, double* kf)
965 {
966  if (m_stickingData.empty()) {
967  return;
968  }
969 
970  static const int cacheId = m_cache.getId();
971  CachedArray cached = m_cache.getArray(cacheId);
972  vector_fp& factors = cached.value;
973 
974  SurfPhase& surf = dynamic_cast<SurfPhase&>(thermo(reactionPhaseIndex()));
975  double n0 = surf.siteDensity();
976  if (!cached.validate(n0)) {
977  factors.resize(m_stickingData.size());
978  for (size_t n = 0; n < m_stickingData.size(); n++) {
979  factors[n] = pow(n0, -m_stickingData[n].order);
980  }
981  }
982 
983  for (size_t n = 0; n < m_stickingData.size(); n++) {
984  const StickData& item = m_stickingData[n];
985  if (item.use_motz_wise) {
986  kf[item.index] /= 1 - 0.5 * kf[item.index];
987  }
988  kf[item.index] *= factors[n] * sqrt(T) * item.multiplier;
989  }
990 }
991 
992 }
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Definition: Kinetics.cpp:622
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:913
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:982
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:943
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:251
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:276
void setElectricPotential(doublereal v)
Set the electric potential of this phase (V).
Definition: ThermoPhase.h:327
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:530
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
Definition: Kinetics.h:921
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:949
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.
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
Definition: ThermoPhase.h:494
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
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:179
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:504
virtual bool addReaction(shared_ptr< Reaction > r)
Add a single reaction to the mechanism.
Definition: Kinetics.cpp:522
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized concentrations.
Definition: ThermoPhase.h:403
virtual void getDeltaSSEntropy(doublereal *deltaS)
Return the vector of values for the change in the standard state entropies for each reaction...
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:54
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:262
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:438
vector_fp m_ropnet
Net rate-of-progress for each reaction.
Definition: Kinetics.h:985
const doublereal Pi
Pi.
Definition: ct_defs.h:51
size_t reactionPhaseIndex()
Phase where the reactions occur.
Definition: Kinetics.h:263
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:216
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:424
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:809
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:979
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:143
size_t surfacePhaseIndex()
This returns the integer index of the phase which has ThermoPhase type cSurf.
Definition: Kinetics.h:251
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:873
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:589
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.
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:290
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:288
A kinetics manager for heterogeneous reaction mechanisms.
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:821
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:973
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
Definition: Kinetics.h:910
size_t m_nDim
Number of dimensions of reacting phase (2 for InterfaceKinetics, 1 for EdgeKinetics) ...
virtual void assignShallowPointers(const std::vector< thermo_t *> &tpVector)
Reassign the pointers within the Kinetics object.
Definition: Kinetics.cpp:129
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:925
doublereal electricPotential() const
Returns the electric potential of this phase (V).
Definition: ThermoPhase.h:335
Public interface for kinetics managers.
Definition: Kinetics.h:111
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:310
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.
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:184
std::string reactionString(size_t i) const
Return a string representing the reaction.
Definition: Kinetics.h:671
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:520
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...
const U & getValue(const std::map< T, U > &m, const T &key)
Const accessor for a value in a std::map.
Definition: utilities.h:537
virtual doublereal logStandardConc(size_t k=0) const
Natural logarithm of the standard concentration of the kth species.
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...
Kinetics & operator=(const Kinetics &right)
Definition: Kinetics.cpp:41
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:976
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.
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:200
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.
InterfaceKinetics(thermo_t *thermo=0)
Constructor.
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:388
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
doublereal size(size_t k) const
This routine returns the size of species k.
Definition: Phase.h:414
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:579
void setPhaseExistence(const size_t iphase, const int exists)
Set the existence of a phase in the reaction object.
virtual int type() const
Identifies the kinetics manager type.
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:599
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:496
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: application.cpp:29
void getConcentrations(doublereal *const c) const
Get the species concentrations (kmol/m^3).
Definition: Phase.cpp:595
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:485
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:428
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
virtual Kinetics * duplMyselfAsKinetics(const std::vector< thermo_t *> &tpVector) const
Duplication routine for objects which inherit from Kinetics.
std::string equation() const
The chemical equation for this reaction.
Definition: Reaction.cpp:89