Improve ReactorNet handling of bulk phase interactions

[speth]

Abstract

The handling of surface reactions in the ReactorNet framework is mostly focused on catalytic surfaces that only interact with a single bulk (gas) phase, or involve deposition/ablation with an effectively infinite reservoir. Surfaces that define reactions involving net mass transfer between two bulk phases can be defined and mostly work, but they have some rough edges. The goals of this enhancement are:

1. To simplify setup of reactor networks involving multiple bulk phases connected by surfaces
2. To ensure that such systems are integrated correctly and without duplicating calculations

Motivation

Consider the following example, using the input file CaCO3-calcination-surface.yaml.txt (which needs to be renamed due to GitHub's file type limitations) adapted from this Users' Group post:

import numpy as np
import cantera as ct

p = ct.one_atm # pressure
T = 1080.0 # surface temperature

gas = ct.Solution("CaCO3-calcination-surface.yaml", "gas")
CaCO3_phase = ct.Solution("CaCO3-calcination-surface.yaml", "CaCO3_s")
CaO_phase = ct.Solution("CaCO3-calcination-surface.yaml", "CaO_s")
surf_phase = ct.Interface("CaCO3-calcination-surface.yaml", "CaO_surf",
                          adjacent=[gas, CaCO3_phase, CaO_phase])

surf_phase.TPX = T, p, 'CaO(S):1, CaCO3(S):1, CO2(S):0'
gas.TPX = T, p, 'H2:1, CO2:0'
CaCO3_phase.TP = T, p
CaO_phase.TP = T, p

r = ct.IdealGasConstPressureReactor(gas, energy='off')
rCaCO3 = ct.ConstPressureReactor(CaCO3_phase, energy='off')
rCaO = ct.ConstPressureReactor(CaO_phase, energy='off')

rCaCO3.volume = 0.0005
rCaO.volume = 0.0005

surf = ct.ReactorSurface(surf_phase, r=r, A=1)
surf = ct.ReactorSurface(surf_phase, r=rCaCO3, A=1)
surf = ct.ReactorSurface(surf_phase, r=rCaO, A=1)
net = ct.ReactorNet([r, rCaCO3, rCaO])

states = ct.SolutionArray(gas, extra=['t','surf_coverages'])

t_end = 0.001
t = 0

masses = []
while t < t_end:
    t = net.step()
    states.append(r.thermo.state, t=t, surf_coverages=surf.coverages)
    masses.append((r.mass, rCaCO3.mass, rCaO.mass))

The issues to be resolved here are:

(1) The need to define 3 ReactorSurface objects. Currently, a ReactorSurface always "belongs" to a single Reactor object, and is responsible for determining its own surface species' rates of production and modifying the rates of production for the Reactor that owns it. It cannot modify the rates of species production for other phases/reactors. This approach means that the system of equations includes the surface species three times each. Also, it requires the phases adjacent to each surface always having their state synchronized when the surface rates are evaluated, which will only happen if each reactor has an independent Solution object, which relates this to #201.

(2) The surface reaction rates are not limited by the mass of one of the adjacent bulk phases reaching zero.

Possible Solutions

From an API perspective, I think it would be preferable to simply define the ReactorSurface as being adjacent to several Reactors, for example:

surf = ct.ReactorSurface(surf_phase, r=[r, rCaCO3, rCaO], A=1)

Internally, I'm not quite sure what the best way to implement this change would be.

The least invasive approach would be to try to leave the current structure, where the ReactorSurface is managed by one of the bulk phases, but then needs to update components of the RHS/LHS arrays that affect other adjacent Reactors. However, this may not be possible because Reactor::eval is only provided with its own "local" segment of the overall RHS/LHS arrays.

Another approach would be to elevate ReactorSurface by a level in the evaluation of the equations, where it would be handled directly by ReactorNet. Here, ReactorSurface::eval would be responsible for computing the governing equations for each surface species, and for storing an array with the net production rates of the bulk phase species for each adjacent bulk phase. These would then be used by ReactorXyz::eval to compute the RHS/LHS terms for its bulk phase species.

I think either of these approaches require #201 to be resolved first.

I haven't yet thought about how to handle the mass of one of the phases going to zero.

References

This was motivated by the example posted on the Users' Group

[ischoegl]

Internally, I'm not quite sure what the best way to implement this change would be.

That's indeed a tricky question.

I haven't yet thought about how to handle the mass of one of the phases going to zero.

This makes me think of one of the (currently abandoned) things I had pondered a while ago #61 (especially what was outlined for diamond_cvd.py.) It's not necessarily directly related, but appears relevant.

With respect to #61, my instinct would be to decouple ReactorSurface from its associated bulk phase, which is more intrusive but imho a more future-proof approach?

[speth]

Thanks for the reminder of #61. I think resolving this feature using the second approach I described, which would decouple the rate evaluation for each surface from the evaluation of a specific Reactor, would also resolve several of the issues that were raised in #61. Specifically, I think it would make it possible to use a combination of a Reservoir and a ReactorSurface to efficiently replace InterfaceKinetics::advanceCoverages, and, with the eventual completion of Cantera/cantera#1021, to replace InterfaceKinetics::solvePseudoSteadyStateProblem.