Cantera  3.1.0a1
Loading...
Searching...
No Matches
InterfaceRate.cpp
Go to the documentation of this file.
1//! @file InterfaceRate.cpp
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
11#include "cantera/base/AnyMap.h"
13
14namespace Cantera
15{
16
18{
19 throw CanteraError("InterfaceData::update",
20 "Missing state information: 'InterfaceData' requires species coverages.");
21}
22
23void InterfaceData::update(double T, const vector<double>& values)
24{
25 warn_user("InterfaceData::update",
26 "This method does not update the site density.");
28 sqrtT = sqrt(T);
29 if (coverages.size() == 0) {
30 coverages = values;
31 logCoverages.resize(values.size());
32 } else if (values.size() == coverages.size()) {
33 std::copy(values.begin(), values.end(), coverages.begin());
34 } else {
35 throw CanteraError("InterfaceData::update",
36 "Incompatible lengths of coverage arrays: received {} elements while "
37 "{} are required.", values.size(), coverages.size());
38 }
39 for (size_t n = 0; n < coverages.size(); n++) {
40 logCoverages[n] = std::log(std::max(coverages[n], Tiny));
41 }
42}
43
44bool InterfaceData::update(const ThermoPhase& phase, const Kinetics& kin)
45{
46 int mf = 0;
47 for (size_t n = 0; n < kin.nPhases(); n++) {
48 mf += kin.thermo(n).stateMFNumber();
49 }
50
51 double T = phase.temperature();
52 bool changed = false;
53 const auto& surf = dynamic_cast<const SurfPhase&>(kin.thermo(0));
54 double site_density = surf.siteDensity();
55 if (density != site_density) {
56 density = surf.siteDensity();
57 changed = true;
58 }
59 if (T != temperature) {
61 sqrtT = sqrt(T);
62 changed = true;
63 }
64 if (changed || mf != m_state_mf_number) {
65 surf.getCoverages(coverages.data());
66 for (size_t n = 0; n < coverages.size(); n++) {
67 logCoverages[n] = std::log(std::max(coverages[n], Tiny));
68 }
69 for (size_t n = 0; n < kin.nPhases(); n++) {
70 size_t start = kin.kineticsSpeciesIndex(0, n);
71 const auto& ph = kin.thermo(n);
72 electricPotentials[n] = ph.electricPotential();
73 ph.getPartialMolarEnthalpies(partialMolarEnthalpies.data() + start);
74 ph.getStandardChemPotentials(standardChemPotentials.data() + start);
75 size_t nsp = ph.nSpecies();
76 for (size_t k = 0; k < nsp; k++) {
77 // only used for exchange current density formulation
78 standardConcentrations[k + start] = ph.standardConcentration(k);
79 }
80 }
82 changed = true;
83 }
84 return changed;
85}
86
87void InterfaceData::perturbTemperature(double deltaT)
88{
89 throw NotImplementedError("InterfaceData::perturbTemperature");
90}
91
92InterfaceRateBase::InterfaceRateBase()
93 : m_siteDensity(NAN)
94 , m_acov(0.)
95 , m_ecov(0.)
96 , m_mcov(0.)
97 , m_chargeTransfer(false)
98 , m_exchangeCurrentDensityFormulation(false)
99 , m_beta(0.5)
100 , m_deltaPotential_RT(NAN)
101 , m_deltaGibbs0_RT(NAN)
102 , m_prodStandardConcentrations(NAN)
103{
104}
105
107{
108 if (node.hasKey("coverage-dependencies")) {
110 node["coverage-dependencies"].as<AnyMap>(), node.units());
111 }
112 if (node.hasKey("beta")) {
113 m_beta = node["beta"].asDouble();
114 }
115 m_exchangeCurrentDensityFormulation = node.getBool(
116 "exchange-current-density-formulation", false);
117}
118
120{
121 if (!m_cov.empty()) {
122 AnyMap deps;
124 node["coverage-dependencies"] = std::move(deps);
125 }
126 if (m_chargeTransfer) {
127 if (m_beta != 0.5) {
128 node["beta"] = m_beta;
129 }
130 if (m_exchangeCurrentDensityFormulation) {
131 node["exchange-current-density-formulation"] = true;
132 }
133 }
134}
135
137 const UnitSystem& units)
138{
139 m_cov.clear();
140 m_ac.clear();
141 m_ec.clear();
142 m_mc.clear();
143 m_lindep.clear();
144 for (const auto& [species, coeffs] : dependencies) {
145 double a, m;
146 vector<double> E(5, 0.0);
147 if (coeffs.is<AnyMap>()) {
148 auto& cov_map = coeffs.as<AnyMap>();
149 a = cov_map["a"].asDouble();
150 m = cov_map["m"].asDouble();
151 if (cov_map["E"].isScalar()) {
152 m_lindep.push_back(true);
153 E[1] = units.convertActivationEnergy(cov_map["E"], "K");
154 } else {
155 m_lindep.push_back(false);
156 auto& E_temp = cov_map["E"].asVector<AnyValue>(1, 4);
157 for (size_t i = 0; i < E_temp.size(); i++) {
158 E[i+1] = units.convertActivationEnergy(E_temp[i], "K");
159 }
160 }
161 } else {
162 auto& cov_vec = coeffs.asVector<AnyValue>(3);
163 a = cov_vec[0].asDouble();
164 m = cov_vec[1].asDouble();
165 if (cov_vec[2].isScalar()) {
166 m_lindep.push_back(true);
167 E[1] = units.convertActivationEnergy(cov_vec[2], "K");
168 } else {
169 m_lindep.push_back(false);
170 auto& E_temp = cov_vec[2].asVector<AnyValue>(1, 4);
171 for (size_t i = 0; i < E_temp.size(); i++) {
172 E[i+1] = units.convertActivationEnergy(E_temp[i], "K");
173 }
174 }
175 }
176 addCoverageDependence(species, a, m, E);
177 }
178}
179
181{
182 for (size_t k = 0; k < m_cov.size(); k++) {
183 AnyMap dep;
184 dep["a"] = m_ac[k];
185 dep["m"] = m_mc[k];
186 if (m_lindep[k]) {
187 dep["E"].setQuantity(m_ec[k][1], "K", true);
188 } else {
189 vector<AnyValue> E_temp(4);
190 for (size_t i = 0; i < m_ec[k].size() - 1; i++) {
191 E_temp[i].setQuantity(m_ec[k][i+1], "K", true);
192 }
193 dep["E"] = E_temp;
194 }
195 dependencies[m_cov[k]] = std::move(dep);
196 }
197}
198
199void InterfaceRateBase::addCoverageDependence(const string& sp, double a, double m,
200 const vector<double>& e)
201{
202 if (std::find(m_cov.begin(), m_cov.end(), sp) == m_cov.end()) {
203 m_cov.push_back(sp);
204 m_ac.push_back(a);
205 m_ec.push_back(e);
206 m_mc.push_back(m);
207 m_indices.clear();
208 } else {
209 throw CanteraError("InterfaceRateBase::addCoverageDependence",
210 "Coverage for species '{}' is already specified.", sp);
211 }
212}
213
214void InterfaceRateBase::setSpecies(const vector<string>& species)
215{
216 m_indices.clear();
217 for (size_t k = 0; k < m_cov.size(); k++) {
218 auto it = find(species.begin(), species.end(), m_cov[k]);
219 if (it != species.end()) {
220 m_indices.emplace(k, it - species.begin());
221 } else {
222 throw CanteraError("InterfaceRateBase:setSpeciesIndices",
223 "Species list does not contain '{}'.", m_cov[k]);
224 }
225 }
226}
227
229 if (shared_data.ready) {
230 m_siteDensity = shared_data.density;
231 }
232
233 if (m_indices.size() != m_cov.size()) {
234 // object is not set up correctly (setSpecies needs to be run)
235 m_acov = NAN;
236 m_ecov = NAN;
237 m_mcov = NAN;
238 return;
239 }
240 m_acov = 0.0;
241 m_ecov = 0.0;
242 m_mcov = 0.0;
243 for (auto& [iCov, iKin] : m_indices) {
244 m_acov += m_ac[iCov] * shared_data.coverages[iKin];
245 if (m_lindep[iCov]) {
246 m_ecov += m_ec[iCov][1] * shared_data.coverages[iKin];
247 } else {
248 m_ecov += poly4(shared_data.coverages[iKin], m_ec[iCov].data());
249 }
250 m_mcov += m_mc[iCov] * shared_data.logCoverages[iKin];
251 }
252
253 // Update change in electrical potential energy
254 if (m_chargeTransfer) {
256 for (const auto& [iPhase, netCharge] : m_netCharges) {
258 shared_data.electricPotentials[iPhase] * netCharge;
259 }
261 }
262
263 // Update quantities used for exchange current density formulation
264 if (m_exchangeCurrentDensityFormulation) {
265 m_deltaGibbs0_RT = 0.;
267 for (const auto& [k, stoich] : m_stoichCoeffs) {
269 shared_data.standardChemPotentials[k] * stoich;
270 if (stoich > 0.) {
272 shared_data.standardConcentrations[k];
273 }
274 }
276 }
277}
278
280{
282
284 if (!m_chargeTransfer) {
285 return;
286 }
287
288 m_stoichCoeffs.clear();
289 for (const auto& [name, stoich] : rxn.reactants) {
290 m_stoichCoeffs.emplace_back(kin.kineticsSpeciesIndex(name), -stoich);
291 }
292 for (const auto& [name, stoich] : rxn.products) {
293 m_stoichCoeffs.emplace_back(kin.kineticsSpeciesIndex(name), stoich);
294 }
295
296 m_netCharges.clear();
297 for (const auto& [k, stoich] : m_stoichCoeffs) {
298 size_t n = kin.speciesPhaseIndex(k);
299 size_t start = kin.kineticsSpeciesIndex(0, n);
300 double charge = kin.thermo(n).charge(k - start);
301 m_netCharges.emplace_back(n, Faraday * charge * stoich);
302 }
303}
304
305StickingCoverage::StickingCoverage()
306 : m_motzWise(false)
307 , m_explicitMotzWise(false)
308 , m_stickingSpecies("")
309 , m_explicitSpecies(false)
310 , m_surfaceOrder(NAN)
311 , m_multiplier(NAN)
312 , m_factor(NAN)
313{
314}
315
317{
318 m_motzWise = node.getBool("Motz-Wise", false);
319 m_explicitMotzWise = node.hasKey("Motz-Wise");
320 m_stickingSpecies = node.getString("sticking-species", "");
321 m_explicitSpecies = node.hasKey("sticking-species");
322}
323
325{
326 if (m_explicitMotzWise) {
327 node["Motz-Wise"] = m_motzWise;
328 }
329 if (m_explicitSpecies) {
330 node["sticking-species"] = m_stickingSpecies;
331 }
332}
333
335{
336 // Ensure that site density is initialized
337 const ThermoPhase& phase = kin.thermo(0);
338 const auto& surf = dynamic_cast<const SurfPhase&>(phase);
339 m_siteDensity = surf.siteDensity();
340 if (!m_explicitMotzWise) {
341 m_motzWise = kin.thermo().input().getBool("Motz-Wise", false);
342 }
343
344 string sticking_species = m_stickingSpecies;
345 if (sticking_species == "") {
346 // Identify the sticking species if not explicitly given
347 vector<string> gasSpecies;
348 vector<string> anySpecies;
349 for (const auto& [name, stoich] : rxn.reactants) {
350 size_t iPhase = kin.speciesPhaseIndex(kin.kineticsSpeciesIndex(name));
351 if (iPhase != 0) {
352 // Non-interface species. There should be exactly one of these
353 // (either in gas phase or other phase)
354 if (kin.thermo(iPhase).phaseOfMatter() == "gas") {
355 gasSpecies.push_back(name);
356 }
357 anySpecies.push_back(name);
358 }
359 }
360 if (gasSpecies.size() == 1) {
361 // single sticking species in gas phase
362 sticking_species = gasSpecies[0];
363 } else if (anySpecies.size() == 1) {
364 // single sticking species in any phase
365 sticking_species = anySpecies[0];
366 } else if (anySpecies.size() == 0) {
367 throw InputFileError("StickingCoverage::setContext",
368 rxn.input, "No non-interface species found "
369 "in sticking reaction: '{}'", rxn.equation());
370 } else {
371 throw InputFileError("StickingCoverage::setContext",
372 rxn.input, "Multiple non-interface species ({})\nfound in sticking "
373 "reaction: '{}'.\nSticking species must be explicitly specified.",
374 fmt::format("'{}'", fmt::join(anySpecies, "', '")), rxn.equation());
375 }
376 }
377 m_stickingSpecies = sticking_species;
378
379 double surface_order = 0.0;
380 double multiplier = 1.0;
381 // Adjust the A-factor
382 for (const auto& [name, stoich] : rxn.reactants) {
383 size_t iPhase = kin.speciesPhaseIndex(kin.kineticsSpeciesIndex(name));
384 const ThermoPhase& p = kin.thermo(iPhase);
385 size_t k = p.speciesIndex(name);
386 if (name == sticking_species) {
387 multiplier *= sqrt(GasConstant / (2 * Pi * p.molecularWeight(k)));
388 } else {
389 // Non-sticking species. Convert from coverages used in the
390 // sticking probability expression to the concentration units
391 // used in the mass action rate expression. For surface phases,
392 // the dependence on the site density is incorporated when the
393 // rate constant is evaluated, since we don't assume that the
394 // site density is known at this time.
395 double order = getValue(rxn.orders, name, stoich);
396 if (&p == &surf) {
397 multiplier *= pow(surf.size(k), order);
398 surface_order += order;
399 } else {
400 multiplier *= pow(p.standardConcentration(k), -order);
401 }
402 }
403 }
404 m_surfaceOrder = surface_order;
405 m_multiplier = multiplier;
406}
407
408}
Header for reaction rates that occur at interfaces.
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase,...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
Definition AnyMap.h:630
bool getBool(const string &key, bool default_) const
If key exists, return it as a bool, otherwise return default_.
Definition AnyMap.cpp:1515
const string & getString(const string &key, const string &default_) const
If key exists, return it as a string, otherwise return default_.
Definition AnyMap.cpp:1530
A wrapper for a variable whose type is determined at runtime.
Definition AnyMap.h:86
double & asDouble()
Return the held value as a double, if it is a double or a long int.
Definition AnyMap.cpp:824
Base class for exceptions thrown by Cantera classes.
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition AnyMap.h:738
double m_mcov
Coverage term in reaction rate.
void setCoverageDependencies(const AnyMap &dependencies, const UnitSystem &units=UnitSystem())
Set coverage dependencies based on AnyMap node information.
double m_beta
Electrochemistry only.
double m_ecov
Coverage contribution to activation energy.
vector< pair< size_t, double > > m_stoichCoeffs
Pairs of species index and multipliers to calculate enthalpy change.
vector< double > m_ac
Vector holding coverage-specific exponential dependence.
void setParameters(const AnyMap &node)
Perform object setup based on AnyMap node information.
void setSpecies(const vector< string > &species)
Set association with an ordered list of all species associated with a given Kinetics object.
vector< string > m_cov
Vector holding names of coverage species.
double m_prodStandardConcentrations
Products of standard concentrations.
double m_deltaGibbs0_RT
Normalized standard state Gibbs free energy change.
virtual void addCoverageDependence(const string &sp, double a, double m, const vector< double > &e)
Add a coverage dependency for species sp, with exponential dependence a, power-law exponent m,...
void getCoverageDependencies(AnyMap &dependencies) const
Store parameters needed to reconstruct coverage dependencies.
bool m_chargeTransfer
Boolean indicating use of electrochemistry.
vector< double > m_mc
Vector holding coverage-specific power-law exponents.
map< size_t, size_t > m_indices
Map from coverage dependencies stored in this object to the index of the coverage species in the Kine...
vector< pair< size_t, double > > m_netCharges
Pairs of phase index and net electric charges (same order as m_stoichCoeffs)
double m_acov
Coverage contribution to pre-exponential factor.
vector< bool > m_lindep
Vector holding boolean for linear dependence.
void getParameters(AnyMap &node) const
Store parameters needed to reconstruct an identical object.
vector< vector< double > > m_ec
Vector holding coverage-specific activation energy dependence as a 5-membered array of polynomial coe...
double m_siteDensity
Site density [kmol/m^2].
void setContext(const Reaction &rxn, const Kinetics &kin)
Build rate-specific parameters based on Reaction and Kinetics context.
void updateFromStruct(const InterfaceData &shared_data)
Update reaction rate parameters.
double m_deltaPotential_RT
Normalized electric potential energy change.
Public interface for kinetics managers.
Definition Kinetics.h:125
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition Kinetics.h:242
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition Kinetics.h:184
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition Kinetics.h:276
size_t speciesPhaseIndex(size_t k) const
This function takes as an argument the kineticsSpecies index (that is, the list index in the list of ...
Definition Kinetics.cpp:281
An error indicating that an unimplemented function has been called.
double temperature() const
Temperature (K).
Definition Phase.h:562
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:129
int stateMFNumber() const
Return the State Mole Fraction Number.
Definition Phase.h:761
const vector< string > & speciesNames() const
Return a const reference to the vector of species names.
Definition Phase.cpp:148
double molecularWeight(size_t k) const
Molecular weight of species k.
Definition Phase.cpp:383
double 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:538
Abstract base class which stores data about a reaction and its rate parameterization so that it can b...
Definition Reaction.h:25
Composition orders
Forward reaction order with respect to specific species.
Definition Reaction.h:119
string equation() const
The chemical equation for this reaction.
Definition Reaction.cpp:345
bool usesElectrochemistry(const Kinetics &kin) const
Check whether reaction uses electrochemistry.
Definition Reaction.cpp:687
Composition products
Product species and stoichiometric coefficients.
Definition Reaction.h:114
Composition reactants
Reactant species and stoichiometric coefficients.
Definition Reaction.h:111
AnyMap input
Input data used for specific models.
Definition Reaction.h:139
void setStickingParameters(const AnyMap &node)
Perform object setup based on AnyMap node information.
bool m_explicitSpecies
Boolean flag.
string m_stickingSpecies
string identifying sticking species
bool m_explicitMotzWise
Correction cannot be overriden by default.
void getStickingParameters(AnyMap &node) const
Store parameters needed to reconstruct an identical object.
double m_multiplier
multiplicative factor in rate expression
bool m_motzWise
boolean indicating whether Motz & Wise correction is used
double m_surfaceOrder
exponent applied to site density term
void setContext(const Reaction &rxn, const Kinetics &kin)
Build rate-specific parameters based on Reaction and Kinetics context.
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition SurfPhase.h:98
double siteDensity() const
Returns the site density.
Definition SurfPhase.h:216
Base class for a phase with thermodynamic properties.
virtual double standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
virtual string phaseOfMatter() const
String indicating the mechanical phase of the matter in this Phase.
const AnyMap & input() const
Access input data associated with the phase description.
Unit conversion utility.
Definition Units.h:169
double convertActivationEnergy(double value, const string &src, const string &dest) const
Convert value from the units of src to the units of dest, allowing for the different dimensions that ...
Definition Units.cpp:673
R poly4(D x, R *c)
Evaluates a polynomial of order 4.
Definition utilities.h:153
const double Faraday
Faraday constant [C/kmol].
Definition ct_defs.h:131
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition ct_defs.h:120
const double Pi
Pi.
Definition ct_defs.h:68
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition global.h:267
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const double Tiny
Small number to compare differences of mole fractions against.
Definition ct_defs.h:173
const U & getValue(const map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a map.
Definition utilities.h:190
int m_state_mf_number
integer that is incremented when composition changes
vector< double > partialMolarEnthalpies
partial molar enthalpies
bool ready
boolean indicating whether vectors are accessible
double density
used to determine if updates are needed
Data container holding shared data for reaction rate specification with interfaces.
vector< double > logCoverages
logarithm of surface coverages
vector< double > electricPotentials
electric potentials of phases
bool update(const ThermoPhase &bulk, const Kinetics &kin) override
Update data container based on thermodynamic phase state.
vector< double > coverages
surface coverages
vector< double > standardChemPotentials
standard state chemical potentials
vector< double > standardConcentrations
standard state concentrations
double sqrtT
square root of temperature
virtual void update(double T)
Update data container based on temperature T
double temperature
temperature
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...