Cantera  3.2.0a2
Loading...
Searching...
No Matches
InterfaceKinetics.h
Go to the documentation of this file.
1/**
2 * @file InterfaceKinetics.h
3 */
4
5// This file is part of Cantera. See License.txt in the top-level directory or
6// at https://cantera.org/license.txt for license and copyright information.
7
8#ifndef CT_IFACEKINETICS_H
9#define CT_IFACEKINETICS_H
10
11#include "Kinetics.h"
12
13namespace Cantera
14{
15
16// forward declarations
17class SurfPhase;
18class ImplicitSurfChem;
19
20//! A kinetics manager for heterogeneous reaction mechanisms. The reactions are
21//! assumed to occur at a 2D interface between two 3D phases.
22/*!
23 * There are some important additions to the behavior of the kinetics class due
24 * to the presence of multiple phases and a heterogeneous interface. If a
25 * reactant phase doesn't exists, that is, has a mole number of zero, a
26 * heterogeneous reaction can not proceed from reactants to products. Note it
27 * could perhaps proceed from products to reactants if all of the product phases
28 * exist.
29 *
30 * In order to make the determination of whether a phase exists or not actually
31 * involves the specification of additional information to the kinetics object.,
32 * which heretofore has only had access to intrinsic field information about the
33 * phases (for example, temperature, pressure, and mole fraction).
34 *
35 * The extrinsic specification of whether a phase exists or not must be
36 * specified on top of the intrinsic calculation of the reaction rate. This
37 * class carries a set of booleans indicating whether a phase in the
38 * heterogeneous mechanism exists or not.
39 *
40 * Additionally, the class carries a set of booleans around indicating whether a
41 * product phase is stable or not. If a phase is not thermodynamically stable,
42 * it may be the case that a particular reaction in a heterogeneous mechanism
43 * will create a product species in the unstable phase. However, other reactions
44 * in the mechanism will destruct that species. This may cause oscillations in
45 * the formation of the unstable phase from time step to time step within a ODE
46 * solver, in practice. In order to avoid this situation, a set of booleans is
47 * tracked which sets the stability of a phase. If a phase is deemed to be
48 * unstable, then species in that phase will not be allowed to be birthed by the
49 * kinetics operator. Nonexistent phases are deemed to be unstable by default,
50 * but this can be changed.
51 *
52 * @ingroup kineticsmgr
53 */
55{
56public:
57 //! Constructor
58 InterfaceKinetics() = default;
59
60 ~InterfaceKinetics() override;
61
62 void resizeReactions() override;
63
64 string kineticsType() const override {
65 return "surface";
66 }
67
68 //! Set the electric potential in the nth phase
69 /*!
70 * @param n phase Index in this kinetics object.
71 * @param V Electric potential (volts)
72 */
73 void setElectricPotential(int n, double V);
74
75 //! @name Reaction Rates Of Progress
76 //! @{
77
78 //! Equilibrium constant for all reactions including the voltage term
79 /*!
80 * Kc = exp(deltaG/RT)
81 *
82 * where deltaG is the electrochemical potential difference between
83 * products minus reactants.
84 */
85 void getEquilibriumConstants(double* kc) override;
86
87 void getDeltaGibbs(double* deltaG) override;
88
89 void getDeltaElectrochemPotentials(double* deltaM) override;
90 void getDeltaEnthalpy(double* deltaH) override;
91 void getDeltaEntropy(double* deltaS) override;
92
93 void getDeltaSSGibbs(double* deltaG) override;
94 void getDeltaSSEnthalpy(double* deltaH) override;
95 void getDeltaSSEntropy(double* deltaS) override;
96
97 //! @}
98 //! @name Reaction Mechanism Informational Query Routines
99 //! @{
100
101 void getActivityConcentrations(double* const conc) override;
102 void getFwdRateConstants(double* kfwd) override;
103 void getRevRateConstants(double* krev, bool doIrreversible=false) override;
104
105 //! @}
106 //! @name Reaction Mechanism Construction
107 //! @{
108
109 //! Add a thermo phase to the kinetics manager object.
110 /*!
111 * This must be done before the function init() is called or
112 * before any reactions are input. The lowest dimensional phase, where reactions
113 * occur, must be added first.
114 *
115 * This function calls Kinetics::addThermo). It also sets the following
116 * fields:
117 *
118 * m_phaseExists[]
119 *
120 * @param thermo Reference to the ThermoPhase to be added.
121 */
122 void addThermo(shared_ptr<ThermoPhase> thermo) override;
123
124 void init() override;
125 void resizeSpecies() override;
126 bool addReaction(shared_ptr<Reaction> r, bool resize=true) override;
127 void modifyReaction(size_t i, shared_ptr<Reaction> rNew) override;
128 void setMultiplier(size_t i, double f) override;
129 //! @}
130
131 //! Internal routine that updates the Rates of Progress of the reactions
132 /*!
133 * This is actually the guts of the functionality of the object
134 */
135 void updateROP() override;
136
137 //! Update properties that depend on temperature
138 /*!
139 * Current objects that this function updates:
140 * m_kdata->m_rfn
141 * updateKc();
142 */
143 void _update_rates_T();
144
145 //! Update properties that depend on the electric potential
146 void _update_rates_phi();
147
148 //! Update properties that depend on the species mole fractions and/or
149 //! concentration,
150 /*!
151 * This method fills out the array of generalized concentrations by calling
152 * method getActivityConcentrations for each phase, which classes
153 * representing phases should overload to return the appropriate quantities.
154 */
155 void _update_rates_C();
156
157 //! Advance the surface coverages in time
158 /*!
159 * This method carries out a time-accurate advancement of the
160 * surface coverages for a specified amount of time.
161 *
162 * @f[
163 * \dot {\theta}_k = \dot s_k (\sigma_k / s_0)
164 * @f]
165 *
166 * @param tstep Time value to advance the surface coverages
167 * @param rtol The relative tolerance for the integrator
168 * @param atol The absolute tolerance for the integrator
169 * @param maxStepSize The maximum step-size the integrator is allowed to take.
170 * If zero, this option is disabled.
171 * @param maxSteps The maximum number of time-steps the integrator can take.
172 * If not supplied, uses the default value in CVodeIntegrator (20000).
173 * @param maxErrTestFails the maximum permissible number of error test failures
174 * If not supplied, uses the default value in CVODES (7).
175 */
176 void advanceCoverages(double tstep, double rtol=1.e-7,
177 double atol=1.e-14, double maxStepSize=0,
178 size_t maxSteps=20000, size_t maxErrTestFails=7);
179
180 //! Solve for the pseudo steady-state of the surface problem
181 /*!
182 * This is the same thing as the advanceCoverages() function,
183 * but at infinite times.
184 *
185 * Note, a direct solve is carried out under the hood here,
186 * to reduce the computational time.
187 *
188 * @param ifuncOverride One of the values defined in @ref solvesp_methods.
189 * The default is -1, which means that the program will decide.
190 * @param timeScaleOverride When a pseudo transient is selected this value
191 * can be used to override the default time scale for
192 * integration which is one. When SFLUX_TRANSIENT is used, this
193 * is equal to the time over which the equations are integrated.
194 * When SFLUX_INITIALIZE is used, this is equal to the time used
195 * in the initial transient algorithm, before the equation
196 * system is solved directly.
197 */
198 void solvePseudoSteadyStateProblem(int ifuncOverride = -1,
199 double timeScaleOverride = 1.0);
200
201 void setIOFlag(int ioFlag);
202
203 //! Update the standard state chemical potentials and species equilibrium
204 //! constant entries
205 /*!
206 * Virtual because it is overridden when dealing with experimental open
207 * circuit voltage overrides
208 */
209 virtual void updateMu0();
210
211 //! Update the equilibrium constants and stored electrochemical potentials
212 //! in molar units for all reversible reactions and for all species.
213 /*!
214 * Irreversible reactions have their equilibrium constant set
215 * to zero. For reactions involving charged species the equilibrium
216 * constant is adjusted according to the electrostatic potential.
217 */
218 void updateKc();
219
220 //! Set the existence of a phase in the reaction object
221 /*!
222 * Tell the kinetics object whether a phase in the object exists. This is
223 * actually an extrinsic specification that must be carried out on top of
224 * the intrinsic calculation of the reaction rate. The routine will also
225 * flip the IsStable boolean within the kinetics object as well.
226 *
227 * @param iphase Index of the phase. This is the order within the
228 * internal thermo vector object
229 * @param exists Boolean indicating whether the phase exists or not
230 */
231 void setPhaseExistence(const size_t iphase, const int exists);
232
233 //! Set the stability of a phase in the reaction object
234 /*!
235 * Tell the kinetics object whether a phase in the object is stable.
236 * Species in an unstable phase will not be allowed to have a positive
237 * rate of formation from this kinetics object. This is actually an
238 * extrinsic specification that must be carried out on top of the
239 * intrinsic calculation of the reaction rate.
240 *
241 * While conceptually not needed since kinetics is consistent with thermo
242 * when taken as a whole, in practice it has found to be very useful to
243 * turn off the creation of phases which shouldn't be forming. Typically
244 * this can reduce the oscillations in phase formation and destruction
245 * which are observed.
246 *
247 * @param iphase Index of the phase. This is the order within the
248 * internal thermo vector object
249 * @param isStable Flag indicating whether the phase is stable or not
250 */
251 void setPhaseStability(const size_t iphase, const int isStable);
252
253 //! Gets the phase existence int for the ith phase
254 /*!
255 * @param iphase Phase Id
256 * @return The int specifying whether the kinetics object thinks the phase
257 * exists or not. If it exists, then species in that phase can be
258 * a reactant in reactions.
259 */
260 int phaseExistence(const size_t iphase) const;
261
262 //! Gets the phase stability int for the ith phase
263 /*!
264 * @param iphase Phase Id
265 * @return The int specifying whether the kinetics object thinks the phase
266 * is stable with nonzero mole numbers. If it stable, then the
267 * kinetics object will allow for rates of production of of
268 * species in that phase that are positive.
269 */
270 int phaseStability(const size_t iphase) const;
271
272 //! Gets the interface current for the ith phase
273 /*!
274 * @param iphase Phase Id
275 * @return The double specifying the interface current. The interface current
276 * is useful when charge transfer reactions occur at an interface. It
277 * is defined here as the net positive charge entering the phase
278 * specified by the Phase Id. (Units: A/m^2 for a surface reaction,
279 * A/m for an edge reaction).
280 */
281 double interfaceCurrent(const size_t iphase);
282
283 void setDerivativeSettings(const AnyMap& settings) override;
284 void getDerivativeSettings(AnyMap& settings) const override;
285 Eigen::SparseMatrix<double> fwdRatesOfProgress_ddCi() override;
286 Eigen::SparseMatrix<double> revRatesOfProgress_ddCi() override;
287 Eigen::SparseMatrix<double> netRatesOfProgress_ddCi() override;
288
289protected:
290 //! @name Internal service methods
291 //!
292 //! @note These methods are for internal use, and seek to avoid code duplication
293 //! while evaluating terms used for rate constants, rates of progress, and
294 //! their derivatives.
295 //! @{
296
297
298 //! Multiply rate with inverse equilibrium constant
299 void applyEquilibriumConstants(double* rop);
300
301 //! Process mole fraction derivative
302 //! @param stoich stoichiometry manager
303 //! @param in rate expression used for the derivative calculation
304 //! @return a sparse matrix of derivative contributions for each reaction of
305 //! dimensions nTotalReactions by nTotalSpecies
306 Eigen::SparseMatrix<double> calculateCompositionDerivatives(StoichManagerN& stoich,
307 const vector<double>& in);
308
309 //! Helper function ensuring that all rate derivatives can be calculated
310 //! @param name method name used for error output
311 //! @throw CanteraError if coverage dependence or electrochemical reactions are
312 //! included
313 void assertDerivativesValid(const string& name);
314
315 //! @}
316
317 //! Temporary work vector of length m_kk
318 vector<double> m_grt;
319
320 bool m_redo_rates = false;
321
322 //! Array of concentrations for each species in the kinetics mechanism
323 /*!
324 * An array of generalized concentrations @f$ C_k @f$ that are defined
325 * such that @f$ a_k = C_k / C^0_k, @f$ where @f$ C^0_k @f$ is a standard
326 * concentration/ These generalized concentrations are used by this
327 * kinetics manager class to compute the forward and reverse rates of
328 * elementary reactions. The "units" for the concentrations of each phase
329 * depend upon the implementation of kinetics within that phase. The order
330 * of the species within the vector is based on the order of listed
331 * ThermoPhase objects in the class, and the order of the species within
332 * each ThermoPhase class.
333 */
334 vector<double> m_conc;
335
336 //! Array of activity concentrations for each species in the kinetics object
337 /*!
338 * An array of activity concentrations @f$ Ca_k @f$ that are defined
339 * such that @f$ a_k = Ca_k / C^0_k, @f$ where @f$ C^0_k @f$ is a standard
340 * concentration. These activity concentrations are used by this
341 * kinetics manager class to compute the forward and reverse rates of
342 * elementary reactions. The "units" for the concentrations of each phase
343 * depend upon the implementation of kinetics within that phase. The order
344 * of the species within the vector is based on the order of listed
345 * ThermoPhase objects in the class, and the order of the species within
346 * each ThermoPhase class.
347 */
348 vector<double> m_actConc;
349
350 //! Vector of standard state chemical potentials for all species
351 /*!
352 * This vector contains a temporary vector of standard state chemical
353 * potentials for all of the species in the kinetics object
354 *
355 * Length = m_kk. Units = J/kmol.
356 */
357 vector<double> m_mu0;
358
359 //! Vector of chemical potentials for all species
360 /*!
361 * This vector contains a vector of chemical potentials for all of the
362 * species in the kinetics object
363 *
364 * Length = m_kk. Units = J/kmol.
365 */
366 vector<double> m_mu;
367
368 //! Vector of standard state electrochemical potentials modified by a
369 //! standard concentration term.
370 /*!
371 * This vector contains a temporary vector of standard state electrochemical
372 * potentials + RTln(Cs) for all of the species in the kinetics object
373 *
374 * In order to get the units correct for the concentration equilibrium
375 * constant, each species needs to have an RT ln(Cs) added to its
376 * contribution to the equilibrium constant Cs is the standard concentration
377 * for the species. Frequently, for solid species, Cs is equal to 1.
378 * However, for gases Cs is P/RT. Length = m_kk. Units = J/kmol.
379 */
380 vector<double> m_mu0_Kc;
381
382 //! Vector of phase electric potentials
383 /*!
384 * Temporary vector containing the potential of each phase in the kinetics
385 * object. length = number of phases. Units = Volts.
386 */
387 vector<double> m_phi;
388
389 //! Pointer to the single surface phase
390 SurfPhase* m_surf = nullptr;
391
392 //! Pointer to the Implicit surface chemistry object
393 /*!
394 * Note this object is owned by this InterfaceKinetics object. It may only
395 * be used to solve this single InterfaceKinetics object's surface problem
396 * uncoupled from other surface phases.
397 */
399
400 bool m_ROP_ok = false;
401
402 //! Current temperature of the data
403 double m_temp = 0.0;
404
405 //! Int flag to indicate that some phases in the kinetics mechanism are
406 //! non-existent.
407 /*!
408 * We change the ROP vectors to make sure that non-existent phases are
409 * treated correctly in the kinetics operator. The value of this is equal
410 * to the number of phases which don't exist.
411 */
413
414 //! Vector of booleans indicating whether phases exist or not
415 /*!
416 * Vector of booleans indicating whether a phase exists or not. We use this
417 * to set the ROP's so that unphysical things don't happen. For example, a
418 * reaction can't go in the forwards direction if a phase in which a
419 * reactant is present doesn't exist. Because InterfaceKinetics deals with
420 * intrinsic quantities only normally, nowhere else is this extrinsic
421 * concept introduced except here.
422 *
423 * length = number of phases in the object. By default all phases exist.
424 */
425 vector<bool> m_phaseExists;
426
427 //! Vector of int indicating whether phases are stable or not
428 /*!
429 * Vector of booleans indicating whether a phase is stable or not under
430 * the current conditions. We use this to set the ROP's so that
431 * unphysical things don't happen.
432 *
433 * length = number of phases in the object. By default all phases are stable.
434 */
435 vector<int> m_phaseIsStable;
436
437 //! Vector of vector of booleans indicating whether a phase participates in
438 //! a reaction as a reactant
439 /*!
440 * m_rxnPhaseIsReactant[j][p] indicates whether a species in phase p
441 * participates in reaction j as a reactant.
442 */
443 vector<vector<bool>> m_rxnPhaseIsReactant;
444
445 //! Vector of vector of booleans indicating whether a phase participates in a
446 //! reaction as a product
447 /*!
448 * m_rxnPhaseIsReactant[j][p] indicates whether a species in phase p
449 * participates in reaction j as a product.
450 */
451 vector<vector<bool>> m_rxnPhaseIsProduct;
452
453 int m_ioFlag = 0;
454
455 //! Number of dimensions of reacting phase (2 for InterfaceKinetics, 1 for
456 //! EdgeKinetics)
457 size_t m_nDim = 2;
458
459 //! Buffers for partial rop results with length nReactions()
460 vector<double> m_rbuf0;
461 vector<double> m_rbuf1;
462
463 //! A flag used to neglect rate coefficient coverage dependence in derivative
464 //! formation
466 //! A flag used to neglect electrochemical contributions in derivative formation
468 //! Relative tolerance used in developing numerical portions of specific derivatives
469 double m_jac_rtol_delta = 1e-8;
470 //! A flag stating if the object uses electrochemistry
472 //! A flag stating if the object has coverage dependent rates
474};
475
476}
477
478#endif
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:431
Advances the surface coverages of the associated set of SurfacePhase objects in time.
A kinetics manager for heterogeneous reaction mechanisms.
bool addReaction(shared_ptr< Reaction > r, bool resize=true) override
Add a single reaction to the mechanism.
int phaseExistence(const size_t iphase) const
Gets the phase existence int for the ith phase.
double interfaceCurrent(const size_t iphase)
Gets the interface current for the ith phase.
SurfPhase * m_surf
Pointer to the single surface phase.
vector< vector< bool > > m_rxnPhaseIsReactant
Vector of vector of booleans indicating whether a phase participates in a reaction as a reactant.
string kineticsType() const override
Identifies the Kinetics manager type.
vector< int > m_phaseIsStable
Vector of int indicating whether phases are stable or not.
Eigen::SparseMatrix< double > revRatesOfProgress_ddCi() override
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
size_t m_nDim
Number of dimensions of reacting phase (2 for InterfaceKinetics, 1 for EdgeKinetics)
void getFwdRateConstants(double *kfwd) override
Return the forward rate constants.
vector< bool > m_phaseExists
Vector of booleans indicating whether phases exist or not.
void getDeltaSSGibbs(double *deltaG) override
Return the vector of values for the reaction standard state Gibbs free energy change.
double m_temp
Current temperature of the data.
void setPhaseStability(const size_t iphase, const int isStable)
Set the stability of a phase in the reaction object.
void getDeltaGibbs(double *deltaG) override
Return the vector of values for the reaction Gibbs free energy change.
bool m_jac_skip_coverage_dependence
A flag used to neglect rate coefficient coverage dependence in derivative formation.
void _update_rates_phi()
Update properties that depend on the electric potential.
vector< double > m_phi
Vector of phase electric potentials.
void assertDerivativesValid(const string &name)
Helper function ensuring that all rate derivatives can be calculated.
double m_jac_rtol_delta
Relative tolerance used in developing numerical portions of specific derivatives.
void setMultiplier(size_t i, double f) override
Set the multiplier for reaction i to f.
void getDeltaSSEntropy(double *deltaS) override
Return the vector of values for the change in the standard state entropies for each reaction.
void advanceCoverages(double tstep, double rtol=1.e-7, double atol=1.e-14, double maxStepSize=0, size_t maxSteps=20000, size_t maxErrTestFails=7)
Advance the surface coverages in time.
InterfaceKinetics()=default
Constructor.
void getDeltaSSEnthalpy(double *deltaH) override
Return the vector of values for the change in the standard state enthalpies of reaction.
void resizeSpecies() override
Resize arrays with sizes that depend on the total number of species.
vector< double > m_grt
Temporary work vector of length m_kk.
void getActivityConcentrations(double *const conc) override
Get the vector of activity concentrations used in the kinetics object.
void getRevRateConstants(double *krev, bool doIrreversible=false) override
Return the reverse rate constants.
void getEquilibriumConstants(double *kc) override
Equilibrium constant for all reactions including the voltage term.
bool m_has_electrochemistry
A flag stating if the object uses electrochemistry.
void updateKc()
Update the equilibrium constants and stored electrochemical potentials in molar units for all reversi...
void setPhaseExistence(const size_t iphase, const int exists)
Set the existence of a phase in the reaction object.
void updateROP() override
Internal routine that updates the Rates of Progress of the reactions.
void applyEquilibriumConstants(double *rop)
Multiply rate with inverse equilibrium constant.
void getDerivativeSettings(AnyMap &settings) const override
Retrieve derivative settings.
void getDeltaEnthalpy(double *deltaH) override
Return the vector of values for the reactions change in enthalpy.
Eigen::SparseMatrix< double > fwdRatesOfProgress_ddCi() override
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
void modifyReaction(size_t i, shared_ptr< Reaction > rNew) override
Modify the rate expression associated with a reaction.
bool m_jac_skip_electrochemistry
A flag used to neglect electrochemical contributions in derivative formation.
void _update_rates_C()
Update properties that depend on the species mole fractions and/or concentration,.
vector< double > m_rbuf0
Buffers for partial rop results with length nReactions()
int phaseStability(const size_t iphase) const
Gets the phase stability int for the ith phase.
vector< double > m_mu0_Kc
Vector of standard state electrochemical potentials modified by a standard concentration term.
ImplicitSurfChem * m_integrator
Pointer to the Implicit surface chemistry object.
void _update_rates_T()
Update properties that depend on temperature.
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, double timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
bool m_has_coverage_dependence
A flag stating if the object has coverage dependent rates.
virtual void updateMu0()
Update the standard state chemical potentials and species equilibrium constant entries.
void addThermo(shared_ptr< ThermoPhase > thermo) override
Add a thermo phase to the kinetics manager object.
vector< double > m_mu
Vector of chemical potentials for all species.
vector< double > m_actConc
Array of activity concentrations for each species in the kinetics object.
vector< double > m_mu0
Vector of standard state chemical potentials for all species.
void getDeltaElectrochemPotentials(double *deltaM) override
Return the vector of values for the reaction electrochemical free energy change.
void setElectricPotential(int n, double V)
Set the electric potential in the nth phase.
void init() override
Prepare the class for the addition of reactions, after all phases have been added.
void resizeReactions() override
Finalize Kinetics object and associated objects.
int m_phaseExistsCheck
Int flag to indicate that some phases in the kinetics mechanism are non-existent.
Eigen::SparseMatrix< double > calculateCompositionDerivatives(StoichManagerN &stoich, const vector< double > &in)
Process mole fraction derivative.
Eigen::SparseMatrix< double > netRatesOfProgress_ddCi() override
Calculate derivatives for net rates-of-progress with respect to species concentration at constant tem...
void getDeltaEntropy(double *deltaS) override
Return the vector of values for the reactions change in entropy.
vector< vector< bool > > m_rxnPhaseIsProduct
Vector of vector of booleans indicating whether a phase participates in a reaction as a product.
void setDerivativeSettings(const AnyMap &settings) override
Set/modify derivative settings.
vector< double > m_conc
Array of concentrations for each species in the kinetics mechanism.
Public interface for kinetics managers.
Definition Kinetics.h:126
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition Kinetics.h:240
This class handles operations involving the stoichiometric coefficients on one side of a reaction (re...
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition SurfPhase.h:98
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595