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} \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.

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

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

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().

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

The surface production rate (second term in the equation above) is computed in the loop over 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);
            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;

            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. This is the source term for the surface species, and nk is the number of surface species associated with this particular wall. This is used in the surface coverage equation, which is covered below.

The second is happening inside a loop over m_nsp, which is a global variable (hence the m_ prefix) that stores the total number of species in the gas phase. This second term is the surface rate of production that is going into the species balance equation.

Cantera/Reactor Equations

From charlesreid1

Contents