GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

Initial commit for polynomial Ea coverage dependency

Browse Source 
Add documentation
This commit is contained in:
Jongyoon Bae 2023-06-02 15:10:22 -04:00 committed by Ray Speth
parent 952e465ec9
commit deb9bf87c4
5 changed files with 97 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
doc/sphinx/yaml/reactions.rst
View File

@ -388,7 +388,9 @@ Includes the fields of an :ref:`sec-yaml-elementary` reaction plus:
``E``
Activation energy dependence on coverage, which uses the same sign convention
as the leading-order activation energy term
as the leading-order activation energy term. This can be a scalar value for
the linear dependency or a list of four values for the polynomial dependency
given in the order of 1st, 2nd, 3rd, and 4th-order coefficients
or a list containing the three elements above, in the given order.
@ -403,12 +405,21 @@ Examples::
rate-constant: {A: 3.7e21 cm^2/mol/s, b: 0, Ea: 67400 J/mol}
coverage-dependencies: {H(s): {a: 0, m: 0, E: -6000 J/mol}}
- equation: 2 O(S) => O2 + 2 Pt(S)
rate-constant: {A: 3.7e+21, b: 0, Ea: 213200 J/mol}
coverage-dependencies: {O(S): {a: 0.0, m: 0.0,
E: [1.0e3 J/mol, 3.0e3 J/mol , -7.0e4 J/mol , 5.0e3 J/mol]}
- equation: CH4 + PT(S) + O(S) => CH3(S) + OH(S)
rate-constant: {A: 5.0e+18, b: 0.7, Ea: 4.2e+04}
coverage-dependencies:
O(S): [0, 0, 8000]
PT(S): [0, -1.0, 0]
- equation: 2 O(S) => O2 + 2 Pt(S)
rate-constant: {A: 3.7e+21, b: 0, Ea: 213200 J/mol}
coverage-dependencies:
O(S): [0, 0, [1.0e6, 3.0e6, -7.0e7, 5.0e6]]
.. _sec-yaml-interface-Blowers-Masel:

27
include/cantera/kinetics/InterfaceRate.h
View File

@ -84,9 +84,20 @@ struct InterfaceData : public BlowersMaselData
* It is evident that this expression combines a regular modified Arrhenius rate
* expression \f$ A T^b \exp \left( - \frac{E_a}{RT} \right) \f$ with coverage-related
* terms, where the parameters \f$ (a_k, E_k, m_k) \f$ describe the dependency on the
* surface coverage of species \f$ k, \theta_k \f$. The InterfaceRateBase class implements
* terms related to coverage only, which allows for combinations with arbitrary rate
* parameterizations (for example Arrhenius and BlowersMaselRate).
* surface coverage of species \f$ k, \theta_k \f$. In addition to the linear coverage
* dependence on the activation energy modifier \f$ E_k \f$, polynomial coverage
* dependence is also available. When the dependence parameter \f$ E_k \f$ is given as
* a scalar value, the linear dependency is applied whereas if a list of four values
* are given as \f$ [E^{(1)}_k-E^{(4)}_k] \f$, a polynomial dependency is applied as
* \f[
* k_f = A T^b \exp \left( - \frac{E_a}{RT} \right)
* \prod_k 10^{a_k \theta_k} \theta_k^{m_k}
* \exp \left( \frac{- E^{(1)}_k \theta_k - E^{(2)}_k \theta_k^2
* - E^{(3)}_k \theta_k^3 - E^{(4)}_k \theta_k^4}{RT} \right)
* \f]
* The InterfaceRateBase class implements terms related to coverage only, which allows
* for combinations with arbitrary rate parameterizations (for example Arrhenius and
* BlowersMaselRate).
*/
class InterfaceRateBase
{
@ -116,7 +127,7 @@ public:
//! Add a coverage dependency for species *sp*, with exponential dependence
//! *a*, power-law exponent *m*, and activation energy dependence *e*,
//! where *e* is in Kelvin, that is, energy divided by the molar gas constant.
virtual void addCoverageDependence(const string& sp, double a, double m, double e);
virtual void addCoverageDependence(const string& sp, double a, double m, vector_fp e);
//! Boolean indicating whether rate uses exchange current density formulation
bool exchangeCurrentDensityFormulation() {
@ -243,7 +254,11 @@ protected:
map<size_t, size_t> m_indices;
vector<string> m_cov; //!< Vector holding names of coverage species
vector_fp m_ac; //!< Vector holding coverage-specific exponential dependence
vector_fp m_ec; //!< Vector holding coverage-specific activation energy dependence
//! Vector holding coverage-specific activation energy dependence as a
//! 5-membered array of polynomial coeffcients starting from 0th-order to
//! 4th-order coefficients
std::vector<vector_fp> m_ec;
std::vector<bool> m_lindep; //!< Vector holding boolean for linear dependence
vector_fp m_mc; //!< Vector holding coverage-specific power-law exponents
private:
@ -434,7 +449,7 @@ public:
return RateType::activationEnergy() + m_ecov * GasConstant;
}
void addCoverageDependence(const string& sp, double a, double m, double e) override {
void addCoverageDependence(const string& sp, double a, double m, vector_fp e) override {
InterfaceRateBase::addCoverageDependence(sp, a, m, e);
RateType::setCompositionDependence(true);
}

52
src/kinetics/InterfaceRate.cpp
View File

@ -9,6 +9,7 @@
#include "cantera/thermo/ThermoPhase.h"
#include "cantera/thermo/SurfPhase.h"
#include "cantera/base/AnyMap.h"
#include "cantera/base/utilities.h"
namespace Cantera
{
@ -140,18 +141,38 @@ void InterfaceRateBase::setCoverageDependencies(const AnyMap& dependencies,
m_ac.clear();
m_ec.clear();
m_mc.clear();
m_lindep.clear();
for (const auto& [species, coeffs] : dependencies) {
double a, E, m;
double a, m;
vector_fp E(5, 0.0);
if (coeffs.is<AnyMap>()) {
auto& cov_map = coeffs.as<AnyMap>();
a = cov_map["a"].asDouble();
m = cov_map["m"].asDouble();
E = units.convertActivationEnergy(cov_map["E"], "K");
if (cov_map["E"].isScalar()) {
m_lindep.push_back(true);
E[1] = units.convertActivationEnergy(cov_map["E"], "K");
} else {
m_lindep.push_back(false);
auto& E_temp = cov_map["E"].asVector<AnyValue>(4);
for (size_t i = 0; i < E_temp.size(); i++) {
E[i+1] = units.convertActivationEnergy(E_temp[i], "K");
}
}
} else {
auto& cov_vec = coeffs.asVector<AnyValue>(3);
a = cov_vec[0].asDouble();
m = cov_vec[1].asDouble();
E = units.convertActivationEnergy(cov_vec[2], "K");
if (cov_vec[2].isScalar()) {
m_lindep.push_back(true);
E[1] = units.convertActivationEnergy(cov_vec[2], "K");
} else {
m_lindep.push_back(false);
auto& E_temp = cov_vec[2].asVector<AnyValue>(4);
for (size_t i = 0; i < E_temp.size(); i++) {
E[i+1] = units.convertActivationEnergy(E_temp[i], "K");
}
}
}
addCoverageDependence(species, a, m, E);
}
@ -166,22 +187,39 @@ void InterfaceRateBase::getCoverageDependencies(AnyMap& dependencies,
warn_deprecated("InterfaceRateBase::getCoverageDependencies",
"To be changed after Cantera 3.0: second argument will be removed.");
vector_fp dep(3);
if (m_lindep[k]) {
dep[2] = m_ec[k][1];
} else {
std::vector<AnyValue> dep(3);
vector_fp E_temp(4);
for (size_t i = 0; i < m_ec[k].size() - 1; i++) {
E_temp[i] = m_ec[k][i+1];
}
dep[2] = E_temp;
}
dep[0] = m_ac[k];
dep[1] = m_mc[k];
dep[2] = m_ec[k];
dependencies[m_cov[k]] = std::move(dep);
} else {
AnyMap dep;
dep["a"] = m_ac[k];
dep["m"] = m_mc[k];
dep["E"].setQuantity(m_ec[k], "K", true);
if (m_lindep[k]) {
dep["E"].setQuantity(m_ec[k][1], "K", true);
} else {
std::vector<AnyValue> E_temp(4);
for (size_t i = 0; i < m_ec[k].size() - 1; i++) {
E_temp[i].setQuantity(m_ec[k][i+1], "K", true);
}
dep["E"] = E_temp;
}
dependencies[m_cov[k]] = std::move(dep);
}
}
}
void InterfaceRateBase::addCoverageDependence(const std::string& sp,
double a, double m, double e)
double a, double m, vector_fp e)
{
if (std::find(m_cov.begin(), m_cov.end(), sp) == m_cov.end()) {
m_cov.push_back(sp);
@ -226,7 +264,7 @@ void InterfaceRateBase::updateFromStruct(const InterfaceData& shared_data) {
m_mcov = 0.0;
for (auto& [iCov, iKin] : m_indices) {
m_acov += m_ac[iCov] * shared_data.coverages[iKin];
m_ecov += m_ec[iCov] * shared_data.coverages[iKin];
m_ecov += poly4(shared_data.coverages[iKin], m_ec[iCov].data());
m_mcov += m_mc[iCov] * shared_data.logCoverages[iKin];
}

4
test/data/surface-phases.yaml
View File

@ -148,6 +148,10 @@ Pt-reactions:
Pt-multisite-reactions:
- equation: 2 O(s) <=> O2(s)
rate-constant: {A: 3.1e9 cm^2/mol/s, b: 0.0, Ea: 0.0}
# dummy reaction for polynomial Ea coverage dependence
coverage-dependencies:
{O(s): {a: 0, m: 0, E: [1.0e3 J/mol, 3.0e3 J/mol,
-7.0e4 J/mol, 5.0e3 J/mol]}}
graphite-anode-species:
- name: EC(e)

15
test/kinetics/kineticsFromYaml.cpp
View File

@ -605,8 +605,23 @@ TEST(Kinetics, InterfaceKineticsFromYaml)
auto& cov_map = coverage_deps["H(s)"].as<Cantera::AnyMap>();
EXPECT_DOUBLE_EQ(cov_map["E"].asDouble(), -6e6);
auto R3 = kin->reaction(2);
EXPECT_TRUE(std::dynamic_pointer_cast<StickingArrheniusRate>(R3->rate()));
auto soln2 = newInterface("surface-phases.yaml", "Pt-multi-sites");
auto kin2 = soln2->kinetics();
auto R4 = kin2->reaction(3);
vector_fp coeffs{1.0e6, 3.0e6, -7.0e7, 5.0e6};
const auto rate4 = std::dynamic_pointer_cast<InterfaceArrheniusRate>(R4->rate());
Cantera::AnyMap coverage_deps2;
rate4->getCoverageDependencies(coverage_deps2);
coverage_deps2.applyUnits();
auto& cov_map2 = coverage_deps2["O(s)"].as<Cantera::AnyMap>();
auto& poly_dep = cov_map2["E"].asVector<AnyValue>();
for (size_t i = 0; i < poly_dep.size(); i++) {
EXPECT_DOUBLE_EQ(poly_dep[i].asDouble(), coeffs[i]);
}
}
TEST(Kinetics, BMInterfaceKineticsFromYaml)