Cantera/Surface Coverage
From charlesreid1
In Cantera, the surface coverage is equivalent to a mole fraction of a surface species.
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, #Get Coverages, and #Set 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.
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;
}
}
}Now let's break down this block of code.
First, we check if we have a surface, and whether there are reactions on that surface - that's the kinetics() and surface() methods of 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++) {
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]);Now, if we do have a surface, and we do have kinetics, we do two things: 1. Advance surface coverages 2. Compute the right-hand side of our species coverage equation, to advance the coverages
First, we compute some stuff that will be used in finding the RHS of our species coverage equation:
if (surf && kin) {
double rs0 = 1.0/surf->siteDensity();
size_t nk = surf->nSpecies();
double sum = 0.0;Next, we set the surface temperature, sync our coverages, and 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));Finally, we compute the RHS for each of our surface species. This is a net reaction rate, and is found by multiplying production rate, 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.