Cantera/Reaction Rate Source Terms

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 {dc_{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)

See Cantera/Surface_Coverage

Code

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 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);
}

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 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);
}

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;
}