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

Volume Equation

The equation for reactor volume is:

$ \frac{dV}{dt} = \sum_{walls} K_w A_w \Delta P + F(t) $

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(t) $ is an arbitrary, user-specified function of time.

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 + 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, $ 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.

Energy Equation Code

As with the volume equation, the wall 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;
}