Cantera/Reactor Equations

The thermochemical state of a gas inside a gas reactor is set by specifying the gas's energy state, pressure, and composition. Cantera does this by solving equations for a the gas internal energy U, the reactor volume V, the mass of each species, and the coverage of surface species on walls of the reactor.

Equations from Cantera Source

Volume Equation

The equation for reactor volume is:

${\frac {dV}{dt}}=\sum _{walls}\left(K_{w}A_{w}\Delta P+F_{w}(t)\right)$

where $K_{w}$ is a wall factor that controls the dependence of wall velocity on pressure gradient across the wall, $A_{w}$ is the wall area, $\Delta P$ is the pressure difference across the wall (so that the wall velocity is pressure-dependent), and $F_{w}(t)$ is an arbitrary, user-specified, time-dependent volumetric source term.

In words: the volume change is due to expansion of all walls connected to our reactor, plus a user-defined volume change function.

Volume Equation Code

The RHS of the volume change equation, above, is assembled by summing over each wall's ${\dot {V}}$. This term is computed in Wall::vdot():

/**
 * The volume rate of change is given by
 * \f[ \dot V = K A (P_{left} - P_{right}) + F(t) \f]
 * where \f$ F(t) \f$ is a specified function of time.
 *
 * This method is used by class Reactor to compute the
 * rate of volume change of the reactor.
 */
doublereal Wall::vdot(doublereal t)
{
    double rate1 = m_k * m_area *
                   (m_left->pressure() - m_right->pressure());
    if (m_vf) {
        rate1 += m_area * m_vf->eval(t);
    }
    return rate1;
}

The actual computation of these terms, and their assembly into the RHS of the volume change equation above, happens in the Reactor::evalEqs() function:

    m_vdot = 0.0;

    // compute wall terms
    size_t loc = m_nsp+2;
    fill(m_sdot.begin(), m_sdot.end(), 0.0);
    for (size_t i = 0; i < m_nwalls; i++) {
        int lr = 1 - 2*m_lr[i];
        double vdot = lr*m_wall[i]->vdot(time);
        m_vdot += vdot;

    [...]

    // volume equation
    ydot[1] = m_vdot;

Energy Equation

To track the energy state of the gas in the batch reactor, an internal energy differential equation is solved. This differential equation is of the form:

${\frac {dU}{dt}}=-P{\frac {dV}{dt}}-\sum _{walls}\left(h_{w}A_{w}\Delta T+A_{w}\varepsilon _{w}\sigma \Delta (T^{4})+G_{w}(t)\right)+\sum _{inlets}{\dot {m}}_{i}h_{i,in}+\left(\sum _{outlets}{\dot {m}}_{i}\right)h$

where $h_{w}$ is the heat transfer coefficient for the wall, $A_{w}$ is the wall area, $\Delta T$ is the temperature difference across the wall, $\varepsilon _{w}$ is the emissivity of the wall, $\sigma$ is the Stefan constant, $G_{w}(t)$ is an arbitrary, user-specified, time-dependent energy source term, ${\dot {m}}_{i}$ is the mass flowrate of all streams entering into the reactor, $h_{i,in}$ is the enthalpy of the streams entering into the reactor, and $h$ is the enthalpy of any stream leaving the reactor.

In words: the energy change is due to compression work, heat transfer to the system through walls in the reactor, heat transfer via wall radiation/emissivity, and net enthalpy being convected into/out of the reactor.

Note that this is a more specific form of the general energy balance derived from Reynolds Transport Theorem:

${\frac {\partial \rho U}{\partial t}}+\nabla \cdot \left(\rho uU\right)=-{\displaystyle \int }PdV-\nabla \cdot \mathbf {q}$

Where is Heat of Reaction?

In case you're wondering, "Where is the heat of reaction contained in the internal energy differential equation?"

The answer is, Cantera accounts for heat of reaction directly - not through an actual heat of reaction term.

The definition of the heat of reaction is, after all, defined simply as the difference in energy state between two mixtures,

$\Delta H_{rxn}=H_{prod}-H_{reac}$

So rather than compute the heat of reaction and substitute it explicitly into the internal energy balance, the heat of reaction is taken into account through the composition of the system, and through the conversion of internal energy U into temperature T.

For example, consider a reaction,

$A+B\rightarrow C+D$

that is highly exothermic - a large amount of heat is released when turning two molecules A and B into two molecules C and D. Let's say $H_{prod}=2000$ and $H_{reac}=3000$, so that $\Delta H=-1000$ (units are arbitrary). Then the internal energy change of the system undergoing this reaction is...

ZERO!

The internal energy change is zero, because the total amount of energy in the system is still the same - it's just changed form, from potential energy $A+B$ to thermal energy $C+D+{\text{thermal energy}}$. The internal energy balance is:

$A+B=C+D+{\text{thermal energy}}$

$3000=2000+1000$

This shakes out, when computing the system temperature, because although the energy level in the system is the same, the composition is different - the system is now composed of lower-energy materials, C and D. Therefore, the same amount of energy will give a higher temperature.

If the energy level of the system were increased or decreased by compression work, heat transfer through walls, or flow into or out of the system, then the internal energy level of the system would change, and this would affect the conversion of internal energy into a temperature.

Energy Equation Code

For the reactor energy equation, the RHS of the internal energy ODE is constructed in Reactor::evalEqs(), where all the ODE RHS are constructed:

void Reactor::evalEqs(doublereal time, doublereal* y,
                      doublereal* ydot, doublereal* params)
{
    [...]

It is important, in understanding the importance and layout of this method, to understand that the structure of the solution vector is as follows:

y[0] = energy ODE RHS
y[1] = volume ODE RHS
y[2] = Species 1 ODE RHS
...
y[2+n] = Species n ODE RHS

The RHS is constructed in three parts. The first part is to assemble the heat flux terms, excluding flow and compression energy terms. In the evalEqs method:

    // compute wall terms
    size_t loc = m_nsp+2;
    fill(m_sdot.begin(), m_sdot.end(), 0.0);

    for (size_t i = 0; i < m_nwalls; i++) {
        int lr = 1 - 2*m_lr[i];
        double vdot = lr*m_wall[i]->vdot(time);
        m_vdot += vdot;
        m_Q += lr*m_wall[i]->Q(time);

        [...]

Second, the compression energy term is added (later in the same method):

    /*
     *  Energy equation.
     *  \f[
     *  \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in}
     * - \dot m_{out} h.
     * \f]
     */
    if (m_energy) {
        ydot[0] = - m_thermo->pressure() * m_vdot - m_Q;
    } else {
        ydot[0] = 0.0;
    }

And finally, the in and out terms for flow work are added:

    // add terms for open system
    if (m_open) {
        const doublereal* mf = m_thermo->massFractions();
        doublereal enthalpy = m_thermo->enthalpy_mass();

        // outlets
        for (size_t i = 0; i < m_nOutlets; i++) {
            double mdot_out = m_outlet[i]->massFlowRate(time);
            [...]
            if (m_energy) {
                ydot[0] -= mdot_out * enthalpy;
            }
        }

        // inlets
        for (size_t i = 0; i < m_nInlets; i++) {
            double mdot_in = m_inlet[i]->massFlowRate(time);
            [...]
            if (m_energy) {
                ydot[0] += mdot_in * m_inlet[i]->enthalpy_mass();
            }
        }
    }

For computing the heat flux from each Wall object, the heat flux terms are computed in the Wall class, specifically Wall:Q():

/**
 * The heat flux is given by
 * \f[ Q = h A (T_{left} - T_{right}) + A G(t) \f]
 * where h is the heat transfer coefficient, and
 * \f$ G(t) \f$ is a specified function of time.
 */
doublereal Wall::Q(doublereal t)
{
    double q1 = (m_area * m_rrth) *
                (m_left->temperature() - m_right->temperature());
    if (m_emiss > 0.0) {
        double tl = m_left->temperature();
        double tr = m_right->temperature();
        q1 += m_emiss * m_area * StefanBoltz * (tl*tl*tl*tl - tr*tr*tr*tr);
    }
    if (m_qf) {
        q1 += m_area * m_qf->eval(t);
    }
    return q1;
}

Species Equations (Gas)

The mass of each species $M_{k}$ in the batch reactor is tracked with a differential equation of the form:

${\frac {dM_{k}}{dt}}=MW_{k}\left[V{\dot {\omega }}_{k}+\sum _{walls}A_{w}{\dot {s}}_{k,w}\right]+\sum _{inlets}{\dot {m}}_{i}Y_{k,i,in}-\left(\sum _{outlets}{\dot {m}}_{o}\right)Y_{k}$

where:

• $MW_{k}$ is molecular weight of species k
• ${\dot {\omega }}_{k}$ is the net production rate of species k
• ${\dot {s}}_{k,w}$ is the production rate of species k by the wall (surface)
• ${\dot {m}}_{i}$ is mass flowrate of stream i entering reactor
• $Y_{k,i,in}$ is mass fraction of species k in stream i entering reactor
• ${\dot {m}}_{o}$ is mass flowrate out of reactor
• $Y_{k}$ is mass fraction of species k in batch reactor (in stream exiting reactor)

The four terms going into the RHS of the species mass balance are:

• Volumetric rate of production
• Surface rate of production
• Mass flow into reactor
• Mass flow out of reactor

In words: the species are changed due to gas-phase (volumetric) reactions, surface-phase (surface) reactions, and net convection into/out of the reactor.

The net production rate of species k can be further broken down as:

${\dot {\omega }}_{k}=\sum _{j=1}^{N_{rxns}}\nu _{jk}r_{j}$

where $\nu _{jk}$ is the stoichiometric coefficient of species k in reaction j, and $r_{j}$ is the molar production rate for the reaction j.

Species Equation Code

The species equation RHS is assembled, as with the other equations, in Reactor::evalEqs().

void Reactor::evalEqs(doublereal time, doublereal* y,
                      doublereal* ydot, doublereal* params)
{
    [...]

The volumetric production rate (first term in the equation above) and the inflow and outflow terms are computed after the loop over each wall:

    /* species equations
     *  Equation is:
     *  \dot M_k = \hat W_k \dot\omega_k + \dot m_{in} Y_{k,in}
     *             - \dot m_{out} Y_{k} + A \dot s_k.
     */
    const vector_fp& mw = m_thermo->molecularWeights();
    if (m_chem) {
        m_kin->getNetProductionRates(ydot+2);   // "omega dot"
    } else {
        fill(ydot + 2, ydot + 2 + m_nsp, 0.0);
    }
    for (size_t n = 0; n < m_nsp; n++) {
        ydot[n+2] *= m_vol;     //           moles/s/m^3 -> moles/s
        ydot[n+2] += m_sdot[n];
        ydot[n+2] *= mw[n];
    }

Here, the line

        m_kin->getNetProductionRates(ydot+2);   // "omega dot"

gets the rates of production corresponding to each species, following the equation:

${\dot {\omega }}_{k}=\sum _{j=1}^{N_{rxns}}\nu _{jk}r_{j}$

The surface production rate (second term in the equation above) is reactor-specific, because although each surface has a given surface rate associated with it, the surface rate of production that a reactor sees depends on the composition in that reactor, the thermodynamic conditions in the reactor, etc.

Therefore, the Reactor class loops over each wall that it touches, and computes a surface rate of production for each surface of each wall:

    // compute wall terms
    size_t loc = m_nsp+2;
    fill(m_sdot.begin(), m_sdot.end(), 0.0);
    for (size_t i = 0; i < m_nwalls; i++) {
        [....]
        if (surf && kin) {
            [...]
            kin->getNetProductionRates(DATA_PTR(m_work));
            size_t ns = kin->surfacePhaseIndex();
            size_t surfloc = kin->kineticsSpeciesIndex(0,ns);

            // First RHS term:
            for (size_t k = 1; k < nk; k++) {
                ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
                sum -= ydot[loc + k];
            }
            ydot[loc] = sum;
            loc += nk;

            // Second RHS term:
            double wallarea = m_wall[i]->area();
            for (size_t k = 0; k < m_nsp; k++) {
                m_sdot[k] += m_work[k]*wallarea;
            }
        }
    }

It's important to notice here that there are two RHS terms being computed:

The first is happening inside a loop over nk. nk is the number of surface species associated with this particular wall. Thus, the first RHS term is the source term for the surface species - the change in surface coverage by a particular species due to it being produced or consumed by surface reactions. This is used in the surface coverage equation (i.e., the solid phase species material balance), which is covered below.

The second RHS term is happening inside a loop over m_nsp, which is a class-specific variable (hence the m_ prefix) that stores the total number of species in the gas phase. This second term is the rate of surface production of each gas species. This goes into the gas phase species material balance.

Surface Coverage Equation

See also: Cantera/Surface Coverage

The surface coverage equation for species m is given by:

${\frac {d\theta _{m}}{dt}}={\frac {{\dot {s}}_{m}\sigma _{m}}{\varrho }}$

where:

• ${\dot {s}}_{m}$ is the surface production rate of species m (the sum of surface reaction rates involving species m)
• $\sigma _{m}$ is the total number of surface sites occupied by species m
• $\varrho$ is the site density of the surface

In words: the surface coverage (that is, the concentration of surface species on the surface) only changes due to surface reactions.

In order to really understand this equation, you should understand how Cantera treats surface phase reactions.

In order to model surface reactions, Cantera treats them similar to the way it treats gas phase reactions. In a gas phase reaction, the reaction rate of a simple reaction like:

$A+B\rightarrow ^{k_{ABC}}C$

depends on the concentrations of A and B. Assuming a first order reaction, and calling this the ABC reaction, the reaction rate looks like:

$r_{ABC}=k_{ABC}n_{A}n_{B}$

where $n_{i}$ is the molar concentration of i, with units of moles/volume.

This approach can be extended, by analogy, to surface reactions. If a surface reaction proceeds as (asterisk denotes surface species):

$A+B*\rightarrow ^{k_{surf}}C*$

then the reaction rate looks like:

$r_{surf}=k_{surf}n_{A}n_{B*}$

The problem, though, is that $[B*]$ can't be defined in the same way as $[A]$.

The solution is to define $[B*]$ as a number of "moles" of a surface species, based on two quantities:

• Site density - the number of sites on the surface available to be "occupied" by any surface species
• Surface coverage - the fraction of all sites on the surface that are occupied by a particular species

This makes the total number of "moles" of surface species B:

$n_{B*}=\theta _{B}\varrho$

where $\theta _{B}$ is the surface coverage of species B, and $\varrho$ is the site density of the surface.

The surface coverage, then, is a "mole fraction" of a species on the surface. The surface coverages also evolve with time, due to changing reaction conditions and production rates, so the surface coverages must be evolved in time using a differential equation, just like the species.

Surface Coverage Equation Code

The surface coverage differential equation solution technique is to make the surface coverages part of the solution vector that is being advanced in the Reactor class. The Reactor class already solves a coupled set of ODEs for the volume, energy, and species, so the surface coverages come along for the ride.

This can be seen in the Reactor::evalEqs method:

void Reactor::evalEqs(doublereal time, doublereal* y,
                      doublereal* ydot, doublereal* params)
{
    [....]

There is a master loop where all wall terms are computed, and this is where the bulk of the surface coverage equation code goes:

    // compute wall terms
    size_t loc = m_nsp+2;
    fill(m_sdot.begin(), m_sdot.end(), 0.0);
    for (size_t i = 0; i < m_nwalls; i++) {
        int lr = 1 - 2*m_lr[i];
        double vdot = lr*m_wall[i]->vdot(time);
        m_vdot += vdot;
        m_Q += lr*m_wall[i]->Q(time);
        Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
        SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
        if (surf && kin) {
            double rs0 = 1.0/surf->siteDensity();
            size_t nk = surf->nSpecies();
            double sum = 0.0;
            surf->setTemperature(m_state[0]);
            m_wall[i]->syncCoverages(m_lr[i]);
            kin->getNetProductionRates(DATA_PTR(m_work));
            size_t ns = kin->surfacePhaseIndex();
            size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
            for (size_t k = 1; k < nk; k++) {
                ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
                sum -= ydot[loc + k];
            }
            ydot[loc] = sum;
            loc += nk;

            [...]
        }
    }

The for loop over nk, where the terms are being put into ydot, is where the actual RHS for the surface coverage differential equation is constructed. The m_work array is populated with the surface production rates (just a few lines above), and surf->size(k) gets the total number of surface sites occupied by species k, so the expression:

                ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);

is constructing the RHS of the equation:

${\frac {d\theta _{m}}{dt}}={\frac {{\dot {s}}_{m}\sigma _{m}}{\varrho }}$

Units

A word on units of important quantities.

Volume units are ${\text{m}}^{3}$, area units are ${\text{m}}^{2}$

Site density is ${\frac {\text{sites}}{{\text{m}}^{2}}}$

Equations from Aris

The following are (mostly) equivalent reactor equations from Aris's book The Mathematical Theory of Diffusion and Reaction in Permeable Catalysts (1975).

1D Model

Aris gives non-dimensional forms of the mass and energy balances in a catalyst pellet as follows:

1D Material Balance

For the material balance,

${\dfrac {\partial u_{i}}{\partial \tau }}=\sum _{k=1}^{Nsp}\nabla _{p}\cdot \left(\Delta _{ik}\nabla _{p}u_{k}\right)+\sum _{j=1}^{Nrxn}\phi _{j}^{2}\alpha _{ji}R_{j}\left(u_{i},v\right)$

where p is the dimensionless spatial direction, $\Delta _{ik}$ is a non-dimensionalized diffusivity,

$\Delta _{ik}={\dfrac {D_{ik,eff}}{D}}$

the reaction rate $R_{j}$ is non-dimensionalized using the reaction rate at reference conditions,

$R_{j}(u_{i},v)={\dfrac {{\hat {r}}_{j}(c_{if}u_{i},T_{f}v)}{{\hat {r}}(c_{if},T_{f})}}$

$\alpha _{ji}$ is the stoichiometric coefficient of species i for reaction j, the time is non-dimensionalized as

$\tau ={\dfrac {Dt}{a^{2}\varepsilon }}$

and a is the scaling constant for distance, $\phi _{j}^{2}$ is the Thiele modulus.

1D Energy Balance

For the energy balance,

$\mathbb {L} {\dfrac {\partial v}{\partial \tau }}=\nabla _{p}\cdot \left(\kappa \nabla _{p}v\right)+\sum _{j=1}^{Nrxn}\phi _{j}\beta _{j}R_{j}\left(u_{i},v\right)$

where $\mathbb {L}$ is the inverse Lewis number,

$\mathbb {L} ={\frac {k\varepsilon }{DC_{p}}}$

beta relates to release of energy via reaction,

$\beta _{j}={\dfrac {(-\Delta H_{j})Dc_{f}}{kT_{f}}}$

and $\kappa ={\frac {k_{eff}}{k}}$ is a non-dimensional thermal diffusivity.

Equations from Froment and Bischoff

The following are (mostly) equivalent reactor equations from Froment and Bischoff's book Chemical Reactor Analysis and Design (1990).

1D Model

From p. 403. Conservation eqn for steady state, single rxn carried out in cylindrical tube:

1D Material Balance

$-{\dfrac {d}{dz}}\left(u_{s}C_{A}\right)=r_{A}\rho _{B}$

1D Energy Balance

$u_{s}\rho _{g}c_{p}{\dfrac {dT}{dz}}=\left(-\Delta H\right)r_{A}\rho _{B}-r{\frac {U}{d_{t}}}(T-T_{r})$

1D Momentum Balance

$-{\frac {dp_{t}}{dz}}=f{\frac {\rho _{g}u_{s}^{2}}{d_{p}}}$

Friction factor dependent on shape factors, &c.

Reynolds number, a, b, hollow rings, packed bed of spheres, etc.

2D Equations

The 2 dimensional equations, accounting for radial mixing, are treated with an effective diffusivity method. The equations become coupled two-point boundary value problems. The governing equations become:

Material Balance

$\varepsilon D_{ea}{\dfrac {d^{2}C_{A}}{dz^{2}}}-u_{s}{\frac {dC_{A}}{dz}}-r_{A}\rho _{B}=0$

Energy Equation

$\lambda _{ea}{\frac {d^{2}T}{dz^{2}}}-\rho _{g}u_{s}c_{p}{\frac {dT}{dz}}+\left(-\Delta H\right)r_{A}\rho _{B}-{\frac {4U}{d_{t}}}\left(T-T_{r}\right)=0$

Boundary Conditions

$u_{s}(C_{A0}-C_{A})=-\varepsilon D_{ea}{\frac {dC_{A}}{dz}}\qquad z=0$

$\rho _{g}u_{s}c_{p}(T_{0}-T)=-\lambda _{ea}{\frac {dT}{dz}}$

${\frac {dC_{A}}{dz}}={\frac {dT}{dz}}=0\qquad z=L$