Cantera/Surface Coverage

In Cantera, the surface coverage is equivalent to a mole fraction of a surface species.

For information on surface reactions, see Cantera/Surface Reactions

The heterogeneous reaction depends on a number of "moles" of solid species. The number of "moles" of solid species A is the product of the surface coverage (i.e., the "mole fraction") and the site density (i.e., the "total moles"):

$ N_{A} = \theta_{A} \varrho $

where $ N_{A} $ is the total "moles" of solid species A, $ \theta_{A} $ is the coverage of species A, and $ \varrho $ is the site density of the solid.

This page covers the details of how this surface coverage is dealt with in the Cantera source code, and gives details about how the surface coverage is computed.

Overview of Class Hierarchy

Coverage is actually stored by a SurfPhase object, which is a manager for all things related to the surface. However, rather than forcing Reactor objects to interface directly with SurfPhase objects, Reactor objects deal with Walls instead. Walls, in turn, wrap SurfPhase objects.

The SurfPhase object manages the memory directly. The Wall object wraps the SurfPhase object to eliminate the need for Reactors to deal with memory access issues.

In the Wall class, the coverages are dealt with in three steps: get, set, and calculate.

• Get coverages - these methods are used when the coverage is needed to initialize the surface species mole fractions.
• Calculate coverages - these methods are used to calculate the updated coverage, based on the conditions the surface sees changing.
• Set coverages - once the new coverages are calculated, they need to be set.

This page covers the #Calculating Coverages, #Getting Coverages, and #Setting Coverages methods below.

Getting Coverages

Reactor::getInitialConditions()

Starting from the Reactor class, we see in the getInitialConditions method that we're looping over each Wall object associated with the Reactor, and we set the surface coverages for the wall.

    // set the first component to the total internal
    // energy
    y[0] = m_thermo->intEnergy_mass() * mass;

    // set the second component to the total volume
    y[1] = m_vol;

    // set the remaining components to the surface species
    // coverages on the walls
    size_t loc = m_nsp + 2;
    SurfPhase* surf;
    for (size_t m = 0; m < m_nwalls; m++) {
        surf = m_wall[m]->surface(m_lr[m]);
        if (surf) {
            m_wall[m]->getCoverages(m_lr[m], y + loc);
            loc += surf->nSpecies();
        }
    }

This piece of code is going to the Wall class and calling getCoverages. The arguments it passes are m_lr[m], which indicates whether the wall is left-sided or right-sided, and y+loc, which is a pointer to an array. The code y+loc means "the array y indexed at loc," which is set to the number of species plus two (the two being mass and volume).

The coverages, stored internally by the Wall object, are copied to the vector passed in the getCoverages function call, namely, y+loc.

Now, looking around more at the y array and the role that it plays in the updateState method, we can see that y is the solution vector - it holds all the solutions. Which means that when we call the getCoverages() method on the Wall object, we're passing it a pointer to our solution vector, and telling the Wall object to populate our solution vector with the coverages - that is, with the mole fractions of the surface species associated with the Wall object.

Wall::getCoverages()

The Wall::getCoverages method is not particularly interesting: it just copies data from one array into another.

void Wall::getCoverages(int leftright, doublereal* cov)
{
    if (leftright == 0) {
        copy(m_leftcov.begin(), m_leftcov.end(), cov);
    } else {
        copy(m_rightcov.begin(), m_rightcov.end(), cov);
    }
}

The m_leftcov and m_rightcov are C++ float vectors that hold the coverages of surface species on the left and right sides of the Wall.

SurfPhase

The Cantera Wall object basically wraps a Cantera Kinetics object and a Cantera SurfPhase object. It uses these to handle the underlying surface chemistry. The Wall class handles the interaction between the surface kinetics and the reactor. This happens via the Wall class's m_chem member (an array with two Kinetics objects, one for each side of the wall) and its m_surf member (an array with two SurfPhase objects, one for each side of the wall).

Setting Coverages

More interesting than getting the wall coverages is the setting the wall coverages.

Reactor::updateState()

The Reactor class's updateState method updates the Wall's surface coverages. In the file Reactor.cpp, the Reactor sets the coverages in the following code block:

    size_t loc = m_nsp + 2;
    SurfPhase* surf;
    for (size_t m = 0; m < m_nwalls; m++) {
        surf = m_wall[m]->surface(m_lr[m]);
        if (surf) {
            m_wall[m]->setCoverages(m_lr[m], y+loc);
            loc += surf->nSpecies();
        }
    }

Again, we're passing a flag to indicate whether the wall is left-handed or right-handed, and then we're passing our solution vector, indexed to the surface species (that's y+loc).

Wall::setCoverages()

In the Wall class, the setCoverages method is not doing anything too fancy. It's basically setting the contents of the Wall's coverage arrays (m_leftcov or m_rightcov) equal to the contents of another array (cov):

void Wall::setCoverages(int leftright, const doublereal* cov)
{
    if (leftright == 0) {
        copy(cov, cov + m_nsp[0], m_leftcov.begin());
    } else {
        copy(cov, cov + m_nsp[1], m_rightcov.begin());
    }
}

Calculating Coverages

More interesting than getting or setting coverages, is calculating coverages.

Reactor::evalEqs()

The money shot is in the Reactor::evalEqs method, which computes the wall terms for the Reactor:

// 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; double wallarea = m_wall[i]->area(); for (size_t k = 0; k < m_nsp; k++) { m_sdot[k] += m_work[k]*wallarea; } } }

Let's break down this block of code.

Check for Surface Reactions

The code first loops over each wall, checking if there is a surface phase associated with each wall, and if there is a reaction network on that surface (this happens via calls to the kinetics() and surface() methods of the Wall class):

    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 we do have a surface, and we do have kinetics, we do two things:

1. Compute surface reaction production rates
2. Compute the right-hand side of our species coverage equation, to advance the coverages

Compute Surface Production Rates

If there is a surface with kinetics, the first thing that's done is to store the inverse site density and number of species, for assembling the RHS of the species coverage equation (the second of the two steps here):

        if (surf && kin) {
            double rs0 = 1.0/surf->siteDensity();
            size_t nk = surf->nSpecies();
            double sum = 0.0;

Next, we set the surface temperature; we sync our coverages (more on this in a moment); and we get our net production rates from the surface kinetics:

            surf->setTemperature(m_state[0]);
            m_wall[i]->syncCoverages(m_lr[i]);
            kin->getNetProductionRates(DATA_PTR(m_work));

The Wall class's syncCoverages method sets the coverages in the SurfPhase object (which Wall wraps). As with the get and set methods, this is simply the Wall object wrapping the SurfPhase object such that memory is managed properly. From the Wall class:

void Wall::syncCoverages(int leftright)
{
    if (leftright == 0) {
        m_surf[0]->setCoverages(DATA_PTR(m_leftcov));
    } else {
        m_surf[1]->setCoverages(DATA_PTR(m_rightcov));
    }
}

The getNetProductionRates is a virtual method in the Kinetics class, but to speak abstractly about it, what it does is, calculate the net production rate of each species (the number of species here is the TOTAL number of species across ALL phases). So, if we have a gas phase with 100 species and a surface phase with 5 sruface species, the getNetProductionRates must return 105 production rates. These production rates are all in units of $ \frac{\text{kmol}}{\text{m}^2 \text{s}} $.

From the InterfaceKinetics.cpp class:

void InterfaceKinetics::getNetProductionRates(doublereal* net)
{
    updateROP();
    m_rxnstoich.getNetProductionRates(m_kk,
                                      &m_kdata->m_ropnet[0],
                                      net);
}

Compute RHS of Coverage Equation

Now that the surface coverages have been updated, and surface rates of production have been computed, we can compute the RHS of the surface coverage equation. The surface coverage equation is an ordinary differential equation that has the same form as the differential equations for gas phase species. The surface species are part of the solution vector, and are part of the set of ODEs that are being solved.

The RHS for the surface coverage equation is a net reaction rate (not a volumetric reaction rate); it is found by multiplying the surface production rate, the inverse site density, and the number of surface sites occupied by species k (that's what surf->size(k) computes; see documentation for the SurfPhase::setCoverages() method):

            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;

This allows us to compute RHS for surface coverage ODEs, which goes into the sum variable.