Recombination and dissociation flux BCs in a multi-species context

Assume you have 2 species A and B, associated with concentrations $c_A$ and $c_B$ in the bulk (in m-3).

At a surface in contact with a gas, you can have these three 2nd order reactions

A + A <-> A2
A + B <-> AB
B + B <-> B2

Species A2, AB, and B2 are associated with partial pressures $P_{A2}$, $P_{AB}$ and $P_{B2}$.

The reactions are associated with forward rate constants $K_{rAA}$, $K_{rAB}$ and $K_{rBB}$ in $\mathrm{m^{4} \ s^{-1}}$, and backwards rate constants $K_{dA2}$, $K_{dAB}$ and $K_{dB2}$ expressed in $\mathrm{m^{-2} \ s^{-1} \ Pa^{-1}}$.

The flux of species A and B at the surface is therefore

$$ \varphi_{A} = -2 K_{rAA} \ c_A \ c_A + 2 K_{dA2} P_{A2} - K_{rAB} \ c_A \ c_B + K_{dAB} P_{AB} $$

$$ \varphi_{B} =- 2 K_{rBB} \ c_B \ c_B + 2 K_{dB2} P_{B2} - K_{rAB} \ c_A \ c_B + K_{dAB} P_{AB} $$

Simple implementation

This is a simple design similar to Reaction. The product pressure here (eg. $P_{AB}$) is not dynamically calculated from the particle flux.

import festim as F

A = F.Species(name="A")
B = F.Species(name="B")

surf_reaction_1 = F.SurfaceReactionBC(
    reactant=[A, A], product_pressure=2.0, k_r=0.1, k_d=1, surface=...
)
surf_reaction_2 = F.SurfaceReactionBC(
    reactant=[A, B], product_pressure=0.0, k_r=0.1, k_d=1, surface=...
)
surf_reaction_3 = F.SurfaceReactionBC(
    reactant=[B, B], product_pressure=1.0, k_r=0.1, k_d=1, surface=...
)

Future proof implementation

But what if we want to dynamically calculate the pressure of the gas species. Then maybe we create a new GasSpecies object? this would allow us later on to implement enclosures (see this verification case).

However, a full implicit implementation of this would require Real finite elements which are not currently supported in dolfinx.

import festim as F

A = F.Species(name="A")
B = F.Species(name="B")

A2 = F.GasSpecies(name="A2", initial_pressure=2.0)
AB = F.GasSpecies(name="AB", initial_pressure=0)
B2 = F.GasSpecies(name="B2", initial_pressure=1.0)

surf_reaction_1 = F.SurfaceReactionBC(reactant=[A, A], product=[A2], k_r=0.1, k_d=1, surface=...)
surf_reaction_2 = F.SurfaceReactionBC(reactant=[A, B], product=[AB], k_r=0.1, k_d=1, surface=...)
surf_reaction_3 = F.SurfaceReactionBC(reactant=[B, B], product=[B2], k_r=0.1, k_d=1, surface=...)