Inconsistent Size of `net_production_rates` in `ReactorSurface()` Object

[dylanrubini]

Problem description

My assumption was that the net_production_rates vector stored in ReactorSurface() under kinetics, would only store the surface species. However, this is not the case, instead both gas and surface species are stored in concatonated form. Is this behaviour intentional? For the partial_molar_enthalpies only the surface species are stored.

Steps to reproduce

import csv

import cantera as ct

# unit conversion factors to SI
cm = 0.01
minute = 60.0

#######################################################################
# Input Parameters
#######################################################################

tc = 800.0  # Temperature in Celsius
length = 0.3 * cm  # Catalyst bed length
area = 1.0 * cm**2  # Catalyst bed area
cat_area_per_vol = 1000.0 / cm  # Catalyst particle surface area per unit volume
velocity = 40.0 * cm / minute  # gas velocity
porosity = 0.3  # Catalyst bed porosity

# input file containing the surface reaction mechanism
yaml_file = 'methane_pox_on_pt.yaml'
output_filename = 'surf_pfr2_output.csv'

#####################################################################

t = tc + 273.15  # convert to Kelvin

# import the model and set the initial conditions
surf = ct.Interface(yaml_file, 'Pt_surf')
surf.TP = t, ct.one_atm
gas = surf.adjacent['gas']
gas.TPX = t, ct.one_atm, 'CH4:1, O2:1.5, AR:0.1'

mass_flow_rate = velocity * gas.density * area * porosity

# create a new reactor
r = ct.FlowReactor(gas)
r.area = area
r.surface_area_to_volume_ratio = cat_area_per_vol * porosity
r.mass_flow_rate = mass_flow_rate
r.energy_enabled = False

# Add the reacting surface to the reactor
rsurf = ct.ReactorSurface(surf, r)

sim = ct.ReactorNet([r])

print(surf.n_species)
print(rsurf.kinetics.net_production_rates.shape)
print(rsurf.kinetics.partial_molar_enthalpies.shape)
print(gas.n_species)

The ouput to terminal is:

11

(18,)

(11,)

7

System information

• Cantera version: 3.2.0a1 (dev)
• OS: Apple Silicon M4 (OS 15.3.2)
• Python/MATLAB/other software versions: Python

[ischoegl]

@dylanrubini thank you for reporting ... this is mostly a quirk of how ReactorSurface objects are handled (they are attached to an owning Reactor). There are some open feature requests (e.g., Cantera/enhancements#202) that may directly address this behavior. I agree that the behavior is not intuitive.

[dylanrubini]

@ischoegl ah okay thanks! So then currently, the best way to get the net_production_rates for only the surfaces species at each timestep in the context of a zero dimensional reactor problem would be to slice like this ---> rsurf.kinetics.net_production_rates[n_species_gas:]? Or is there a more correct way?

[speth]

I don't think this is a "quirk" but a consequence of surface kinetics being inherently somewhat complicated.

The kinetics member of a ReactorSurface is an Interface object, which represents the thermodynamics of the surface species and the reactions occurring on that surface. Those reactions produce and consume not only surface species but also gas/bulk phase species, and so must be part of the full creation/destruction/net production rates arrays. This structure complemented by other methods/properties of the Kinetics and InterfaceKinetics classes that provide information on these structures:

>>> surf = ct.Interface('methane_pox_on_pt.yaml', 'Pt_surf')
>>> print(surf.n_species, surf.n_total_species)
11 18
>>> print(surf.kinetics_species_names)
['PT(S)', 'H(S)', 'H2O(S)', 'OH(S)', 'CO(S)', 'CO2(S)', 'CH3(S)', 'CH2(S)',
 'CH(S)', 'C(S)', 'O(S)', 'H2', 'O2', 'H2O', 'CH4', 'CO', 'CO2', 'AR']
>>> print(surf.kinetics_species_index('CO2'))
16
>>> print(surf.phase_index('Pt_surf'))
0

In the case where you are only interested in rates in one phase, methods like InterfaceKinetics.get_creation_rates accept a phase argument, which can be either the phase name or index.

Note that in the case of a steady-state plug flow reactor, as in your example, the net production rates for the surface species should be zero, within the tolerances of the numerical integrator.

[dylanrubini]

@speth, thanks for clarifying.

With this in mind, I want to store the heat_release_rate in a SolutionArray() object during each step of the FlowReactor() time integration. How can I calculate the heat release/absorbtion rate when both surface and gas phase species are involved?

Currently, my understanding is that only the gas-phase heat_release_rate is automatically calculated and stored in the SolutionArray()

For surface species, I guess the heat release rate would be a per unit area quantity(?)

[ischoegl]

@dylanrubini and @speth ... I think I'll stick with my assessment that the ReactorSurface implementation has quirks that, from my perspective, are implementation details that are exposed at the user-facing interface.

To be honest, querying the FlowReactor from the top post should expose gas species, whereas the ReactorSurface should only expose surface species. I still consider the current implementation of the API somewhat problematic, with

>>> r.n_vars  # should not include surface species
22
>>> len(rsurf.kinetics.net_production_rates)  # should not include gas phase reactions
18
>>> rsurf.n_vars  # throws error

representing various inconsistencies (I made sure to check this with stable Cantera 3.1).

From that perspective, @dylanrubini's comment about potentially conflating values with inconsistent units is just another example (which, however, currently throws errors).

[dylanrubini]

@ischoegl @speth just to clarify - i'm trying to calculate the last term on the RHS of the following code taken from the FlowReactor() .cpp source file

    if (m_energy) {
        // energy conservation, Kee 16.58, adiabatic
        // -sum(wk' * Wk * hk) - (perim / Ac) * sum(sk' * Wk * hk)
        // Note: the partial molar enthalpies are in molar units, while Kee uses mass
        //       units, where:
        //              h_mass = h_mole / Wk
        //       hence:
        //              h_mass * Wk = h_mol
        m_thermo->getPartialMolarEnthalpies(m_hk.data());
        for (size_t i = 0; i < m_nsp; ++i) {
            ydot[3] -= m_hk[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
        }
    } else {
        ydot[3] = 0;
    }

In other words, I want to calculate m_hk[i] * hydraulic * m_sdot[i] - which is a heat release/absorbtion rate per unit volume related to surface chemistry.

While I appreciate, the net production rates for the surface species are zero in steady state, there is still a net production rate of gas-phase species that results from surface reactions (and thus a heat release/absorbtion rate from surface-phase chemistry.

I am still a little unclear on how best to correctly calculate this heat absorbtion rate at each time step of a FlowReactor() simulation and store it in the SolutionArray(). Currently, the precalculated heat_release_rate variable is just for gas-phase reactions.

I apologise in advance for my confusion!

[dylanrubini]

Sorry, ignore me, I figured this out now

sdot = rsurf.kinetics.net_production_rates
hk = r.thermo.partial_molar_enthalpies
heat_release_rate_surf = sdot[surf.n_species:] @ hk * cat_area_per_vol

[speth]

@dylanrubini and @speth ... I think I'll stick with my assessment that the ReactorSurface implementation has quirks that, from my perspective, are implementation details that are exposed at the user-facing interface.

To be honest, querying the FlowReactor from the top post should expose gas species, whereas the ReactorSurface should only expose surface species. I still consider the current implementation of the API somewhat problematic, with

>>> r.n_vars  # should not include surface species
22
>>> rsurf.n_vars  # throws error

This must be the case as long as the evaluation of governing equations is handled entirely by Reactor objects. Refactoring things so that a ReactorSurface implements its own version of eval(t, LHS, RHS) and "owns" its own state variables is something I expect to be done as part of implementing Cantera/enhancements#202.

>>> len(rsurf.kinetics.net_production_rates)  # should not include gas phase reactions
18

Why not? These are all production rates on the surface, with units of kmol/m²/s. Where else would these rates be reported?

There is perhaps a case to be made for adding a property of the Reactor and ReactorSurface objects that would combine the production rates integrated over the bulk and all surfaces, giving rates in kmol/s. This would then give values only for the surface species on the ReactorSurface and bulk species in the Reactor.

[ischoegl]

@speth ... I believe the disconnect between expectations and Cantera's traditional implementation is that mechanism definitions do not cleanly map to zero-D Reactor topology.

Why not? These are all production rates on the surface, with units of kmol/m²/s. Where else would these rates be reported?

It is true that these production rates occur at the surface, but imho it is inconsistent to have rate expressions with a different length than the number of species defined for that entity (which is exactly what was reported on top).

Along the lines of Cantera/enhancements#202, I have long eyed some kind of ConnectorNode object (per #1848) to handle interactions between surface phase and bulk phase: this BiLayer (or whatever we may want to call it) would handle the transfer of gas species to the ReactorSurface, which would appear as source terms in the governing equations, with suitable conversions between surface production rates and volumetric equivalents.

This approach would resolve the one missing item of the 'bipartite graph' approach I was envisioning in #1848/#1864, and would also be in line with future extensions that address phase change at an interface (which is currently entirely missing).

[speth]

I think this is more fundamental than the reactor network topology. If you want all arrays for an object that involve species to have the same length, then this essentially becomes an argument against the Python Interface class. The InterfaceKinetics object must be defined over a set of species including all phases participating in reactions.

While you could truncate the arrays to only cover the surface species, that leaves a problem of how to make the information available regarding bulk phase species that are being created/destroyed on the surface. I'll again note that we do already provide convenience methods to this effect like InterfaceKinetics.get_creation_rates(phase).

What is clearly worthwhile is to better document the InterfaceKinetics class to note that all species-sized arrays are dimensioned as n_total_species and how to extract information for particular phases / species. This applies both to the Python class as well as the underlying C++ class.