Cantera/Reaction Rate Source Terms
From charlesreid1
Notes
- Reaction rates are a functon of TPX, thermodynamic state of gas
- phase object wraps an equation of state and wraps kinetic object
- use phase to obtain reaction rates
- Definition of reaction rate
- for a component: sum over each reaction rate
- Cantera phase object - calls
- production rate for component i
- production rate for reaction j
- production rate for ij
- species production rate for component i
- species production rate for reaction j
- species production rate for ij
Definition of Rates
Species Production Rate
In a differential, perfectly uniform control volume in which a chemical reaction is the only phenomena occurring, the species rate of production (species i) can be defined as the change in moles of i with time:
$ \frac{d c_i}{dt} = R_i \qquad i = 1 \dots N_{sp} $
where:
- $ c_i $ is molar volume (units of $ \frac{\text{mol}}{\text{m}^3} $)
- $ R_i $ is the volume-specific species rate of production (units of $ \frac{\text{mol}}{\text{m}^3 \text{s}} $)
Reaction Production Rate
A species production rate can be further broken up into its contributions from each reaction. These rates of production are stoichiometric and known as the reaction rates of production:
$ R_i = \sum_{j=1}^{N_{rxns}} \nu_{ij} \mathfrak{R}_j $
where:
- $ \nu_{ij} $ - stoichiometric coeff of species i in rxn j
- $ \mathfrak{R}_j $ - volume-specific stoichiometric reaction rate for reaction j
Dealing with Phases
Gas Only: Volumetric Reaction Rates
Volumetric reaction rate is: (volume) x (rate of production per unit volume)
Rate of production per unit volume is: (stoichiometric coefficient) x (rate of production of reaction)
Gas and Surface: Surface Reaction Rates
Surface reaction rates are: (area) x (rate of production per unit area)
Rate of production per area is: (mass cat) x (surface area per mass cat) x (surface density of active sites) x (fractional coverage of surface species) x (rate of production of surface reaction)
Code
Gas Reaction Rates
The surface reaction rates of production are computed in the GasKinetics class, in GasKinetics.cpp:
void GasKinetics::getNetProductionRates(doublereal* net)
{
updateROP();
m_rxnstoich.getNetProductionRates(m_kk, &m_ropnet[0], net);
}
Surface Reaction Rates
The surface reaction rates of production are computed in the InterfaceKinetics class, in InterfaceKinetics.cpp:
void InterfaceKinetics::getNetProductionRates(doublereal* net)
{
updateROP();
m_rxnstoich.getNetProductionRates(m_kk,
&m_kdata->m_ropnet[0],
net);
}
Code for Rate of Production
For both surfaces and gases, calculating the rates of production is a two-step process.
First, the updateROP() method is called.
Second, the reaction stoichiometry manager object, ReactionStoichMgr, is called.
Gas updateROP()
void GasKinetics::updateROP()
{
update_rates_C();
update_rates_T();
if (m_ROP_ok) {
return;
}
// copy rate coefficients into ropf
copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin());
// multiply ropf by enhanced 3b conc for all 3b rxns
if (!concm_3b_values.empty()) {
m_3b_concm.multiply(&m_ropf[0], &concm_3b_values[0]);
}
if (m_nfall) {
processFalloffReactions();
}
// multiply by perturbation factor
multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin());
// copy the forward rates to the reverse rates
copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin());
// for reverse rates computed from thermochemistry, multiply
// the forward rates copied into m_ropr by the reciprocals of
// the equilibrium constants
multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin());
// multiply ropf by concentration products
m_rxnstoich.multiplyReactants(&m_conc[0], &m_ropf[0]);
//m_reactantStoich.multiply(m_conc.begin(), ropf.begin());
// for reversible reactions, multiply ropr by concentration
// products
m_rxnstoich.multiplyRevProducts(&m_conc[0], &m_ropr[0]);
//m_revProductStoich.multiply(m_conc.begin(), ropr.begin());
for (size_t j = 0; j != m_ii; ++j) {
m_ropnet[j] = m_ropf[j] - m_ropr[j];
}
m_ROP_ok = true;
}
Surface updateROP
The updateROP method for a surface is found in the InterfaceKinetics class:
void InterfaceKinetics::updateROP()
{
_update_rates_T();
_update_rates_C();
if (m_kdata->m_ROP_ok) {
return;
}
const vector_fp& rf = m_kdata->m_rfn;
const vector_fp& m_rkc = m_kdata->m_rkcn;
vector_fp& ropf = m_kdata->m_ropf;
vector_fp& ropr = m_kdata->m_ropr;
vector_fp& ropnet = m_kdata->m_ropnet;
// copy rate coefficients into ropf
copy(rf.begin(), rf.end(), ropf.begin());
// multiply by perturbation factor
multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
// copy the forward rates to the reverse rates
copy(ropf.begin(), ropf.end(), ropr.begin());
// for reverse rates computed from thermochemistry, multiply
// the forward rates copied into m_ropr by the reciprocals of
// the equilibrium constants
multiply_each(ropr.begin(), ropr.end(), m_rkc.begin());
// multiply ropf by concentration products
m_rxnstoich.multiplyReactants(DATA_PTR(m_conc), DATA_PTR(ropf));
//m_reactantStoich.multiply(m_conc.begin(), ropf.begin());
// for reversible reactions, multiply ropr by concentration
// products
m_rxnstoich.multiplyRevProducts(DATA_PTR(m_conc),
DATA_PTR(ropr));
//m_revProductStoich.multiply(m_conc.begin(), ropr.begin());
// do global reactions
//m_globalReactantStoich.power(m_conc.begin(), ropf.begin());
for (size_t j = 0; j != m_ii; ++j) {
ropnet[j] = ropf[j] - ropr[j];
}
/*
* For reactions involving multiple phases, we must check that the phase
* being consumed actually exists. This is particularly important for
* phases that are stoichiometric phases containing one species with a unity activity
*/
if (m_phaseExistsCheck) {
for (size_t j = 0; j != m_ii; ++j) {
if ((ropr[j] > ropf[j]) && (ropr[j] > 0.0)) {
for (size_t p = 0; p < nPhases(); p++) {
if (m_rxnPhaseIsProduct[j][p]) {
if (! m_phaseExists[p]) {
ropnet[j] = 0.0;
ropr[j] = ropf[j];
if (ropf[j] > 0.0) {
for (size_t rp = 0; rp < nPhases(); rp++) {
if (m_rxnPhaseIsReactant[j][rp]) {
if (! m_phaseExists[rp]) {
ropnet[j] = 0.0;
ropr[j] = ropf[j] = 0.0;
}
}
}
}
}
}
if (m_rxnPhaseIsReactant[j][p]) {
if (! m_phaseIsStable[p]) {
ropnet[j] = 0.0;
ropr[j] = ropf[j];
}
}
}
} else if ((ropf[j] > ropr[j]) && (ropf[j] > 0.0)) {
for (size_t p = 0; p < nPhases(); p++) {
if (m_rxnPhaseIsReactant[j][p]) {
if (! m_phaseExists[p]) {
ropnet[j] = 0.0;
ropf[j] = ropr[j];
if (ropf[j] > 0.0) {
for (size_t rp = 0; rp < nPhases(); rp++) {
if (m_rxnPhaseIsProduct[j][rp]) {
if (! m_phaseExists[rp]) {
ropnet[j] = 0.0;
ropf[j] = ropr[j] = 0.0;
}
}
}
}
}
}
if (m_rxnPhaseIsProduct[j][p]) {
if (! m_phaseIsStable[p]) {
ropnet[j] = 0.0;
ropf[j] = ropr[j];
}
}
}
}
}
}
m_kdata->m_ROP_ok = true;
}