FlowReactor does not conserve sum of mass fractions when last gas phase species is reacting

[speth]

Problem description

The FlowReactor can produce output that does not satisfy the constraint that $\sum Y_k = 1$ when the last species in the mechanism has a non-zero net rate of production from surface reactions. It's a little surprising that the last gas phase species is special here. While there is special handling of the last surface phase species for each attached surface related to this constraint, since the coverages are algebraic variables, there is no corresponding treatment for the gas, where the variables are differential.

Steps to reproduce

Run the following simplification / modification of the surf_pfr.py example, which swaps the last two gas phase species (CO2 and AR).

import cantera as ct

cm = 0.01
minute = 60.0
T0 = 1073.0
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

gas_def = """
phases:
    - name: gas
      thermo: ideal-gas
      elements: [O, H, C, N, Ar]
      species: [{methane_pox_on_pt.yaml/species: [H2, O2, H2O, CH4, CO, AR, CO2]}]
      skip-undeclared-elements: true
      state:
        T: 300.0
        P: 1.01325e+05
        X: {CH4: 0.095, O2: 0.21}
"""
gas = ct.Solution(yaml=gas_def)
surf = ct.Interface('methane_pox_on_pt.yaml', 'Pt_surf', adjacent=[gas])
surf.TP = T0, ct.one_atm
gas.TPX = T0, ct.one_atm, 'CH4:1, O2:1.5, H2:0.1, CO2:0.5'

# 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 = velocity * gas.density * area * porosity
r.energy_enabled = False

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

n = 0
print('    distance       Y_CH4        Y_H2        Y_CO       Y_CO2      sum(Y)')
print('  {:10f}  {:10f}  {:10f}  {:10f}  {:10f}  {:10f}'.format(
      0, *r.thermo['CH4', 'H2', 'CO', 'CO2'].Y, sum(r.thermo.Y)))

while sim.distance < length:
    dist = sim.distance * 1e3  # convert to mm
    sim.step()

    if n % 500 == 0 or (dist > 1 and n % 20 == 0):
        print('  {:10f}  {:10f}  {:10f}  {:10f}  {:10f}  {:10f}'.format(
              dist, *r.thermo['CH4', 'H2', 'CO', 'CO2'].Y, sum(r.thermo.Y)))
    n += 1

Behavior

The above code outputs the following:

    distance       Y_CH4        Y_H2        Y_CO       Y_CO2      sum(Y)
    0.000000    0.186014    0.002337    0.000000    0.255136    1.000000
    0.000000    0.186014    0.002337    0.000000    0.255136    1.000000
    0.000030    0.091973    0.000394    0.008276    0.375955    0.837798
    0.000073    0.052324    0.000170    0.005348    0.431867    0.762734
    0.000074    0.051099    0.000179    0.005947    0.433060    0.761132
    0.000084    0.044205    0.000105    0.008037    0.440768    0.750785
    0.000142    0.036697    0.000034    0.007743    0.451153    0.736842
    0.000334    0.036270    0.000138    0.007733    0.451738    0.736057
    0.015505    0.029418    0.002992    0.013701    0.456384    0.729820
    1.233854    0.006130    0.011906    0.051386    0.458689    0.726725
    1.464645    0.005405    0.012175    0.052745    0.458616    0.726823
    1.710634    0.004743    0.012421    0.053994    0.458544    0.726919
    2.178760    0.003800    0.012770    0.055784    0.458433    0.727068
    2.673283    0.003115    0.013023    0.057093    0.458346    0.727185
    2.988742    0.002776    0.013148    0.057742    0.458301    0.727246

Modifying the species definition to be:

species: [{methane_pox_on_pt.yaml/species: [H2, O2, H2O, CH4, CO, CO2, AR]}]

instead produces output that satisfies the $\sum Y_k = 1$ constraint:

    distance       Y_CH4        Y_H2        Y_CO       Y_CO2      sum(Y)
    0.000000    0.186014    0.002337    0.000000    0.255136    1.000000
    0.000000    0.186014    0.002337    0.000000    0.255136    1.000000
    0.000016    0.130539    0.000705    0.014398    0.384694    1.000000
    0.000073    0.068400    0.000223    0.007079    0.566651    1.000000
    0.000074    0.066891    0.000236    0.007960    0.569408    1.000000
    0.000084    0.058474    0.000135    0.010751    0.588111    1.000000
    0.000136    0.049908    0.000044    0.010521    0.611970    1.000000
    0.000254    0.049396    0.000129    0.010491    0.613422    1.000000
    0.002482    0.047102    0.001208    0.011507    0.618118    1.000000
    0.151922    0.023696    0.010654    0.043737    0.631686    1.000000
    1.136859    0.008976    0.016184    0.069706    0.631265    1.000000
    1.277229    0.008267    0.016445    0.071022    0.631142    1.000000
    1.551396    0.007056    0.016891    0.073281    0.630913    1.000000
    1.949381    0.005808    0.017350    0.075626    0.630653    1.000000
    2.257913    0.005086    0.017615    0.076991    0.630491    1.000000
    2.470240    0.004658    0.017771    0.077800    0.630391    1.000000
    2.733889    0.004187    0.017944    0.078695    0.630278    1.000000

System information

• Cantera version: 3.0.0 (conda-forge package)
• OS: macOS 14.2.1
• Python 3.11.5 (conda-forge)

Additional context

Based on an issue posted in the Cantera Users' Group.