Merge pull request #753 from RemDelaporteMathurin/surface-reactions

Surface reactions

examples/surface_reaction_example.py

@@ -0,0 +1,169 @@
+import festim as F
+import numpy as np
+
+import dolfinx.fem as fem
+import ufl
+
+
+class FluxFromSurfaceReaction(F.SurfaceFlux):
+    def __init__(self, reaction: F.SurfaceReactionBC):
+        super().__init__(
+            F.Species(),  # just a dummy species here
+            reaction.subdomain,
+        )
+        self.reaction = reaction.flux_bcs[0]
+
+    def compute(self, ds):
+        self.value = fem.assemble_scalar(
+            fem.form(self.reaction.value_fenics * ds(self.surface.id))
+        )
+        self.data.append(self.value)
+
+
+my_model = F.HydrogenTransportProblem()
+my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000))
+my_mat = F.Material(name="mat", D_0=1, E_D=0)
+vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat)
+left = F.SurfaceSubdomain1D(id=1, x=0)
+right = F.SurfaceSubdomain1D(id=2, x=1)
+
+my_model.subdomains = [vol, left, right]
+
+H = F.Species("H")
+D = F.Species("D")
+my_model.species = [H, D]
+
+my_model.temperature = 500
+
+surface_reaction_hd = F.SurfaceReactionBC(
+    reactant=[H, D],
+    gas_pressure=0,
+    k_r0=0.01,
+    E_kr=0,
+    k_d0=0,
+    E_kd=0,
+    subdomain=right,
+)
+
+surface_reaction_hh = F.SurfaceReactionBC(
+    reactant=[H, H],
+    gas_pressure=lambda t: ufl.conditional(ufl.gt(t, 1), 2, 0),
+    k_r0=0.02,
+    E_kr=0,
+    k_d0=0.03,
+    E_kd=0,
+    subdomain=right,
+)
+
+surface_reaction_dd = F.SurfaceReactionBC(
+    reactant=[D, D],
+    gas_pressure=0,
+    k_r0=0.01,
+    E_kr=0,
+    k_d0=0,
+    E_kd=0,
+    subdomain=right,
+)
+
+my_model.boundary_conditions = [
+    F.DirichletBC(subdomain=left, value=2, species=H),
+    F.DirichletBC(subdomain=left, value=2, species=D),
+    surface_reaction_hd,
+    surface_reaction_hh,
+    surface_reaction_dd,
+]
+
+H_flux_right = F.SurfaceFlux(H, right)
+H_flux_left = F.SurfaceFlux(H, left)
+D_flux_right = F.SurfaceFlux(D, right)
+D_flux_left = F.SurfaceFlux(D, left)
+HD_flux = FluxFromSurfaceReaction(surface_reaction_hd)
+HH_flux = FluxFromSurfaceReaction(surface_reaction_hh)
+DD_flux = FluxFromSurfaceReaction(surface_reaction_dd)
+my_model.exports = [
+    F.XDMFExport("test.xdmf", H),
+    H_flux_left,
+    H_flux_right,
+    D_flux_left,
+    D_flux_right,
+    HD_flux,
+    HH_flux,
+    DD_flux,
+]
+
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True)
+
+my_model.settings.stepsize = 0.1
+
+my_model.initialise()
+my_model.run()
+
+import matplotlib.pyplot as plt
+
+plt.stackplot(
+    H_flux_left.t,
+    np.abs(H_flux_left.data),
+    np.abs(D_flux_left.data),
+    labels=["H_in", "D_in"],
+)
+plt.stackplot(
+    H_flux_right.t,
+    -np.abs(H_flux_right.data),
+    -np.abs(D_flux_right.data),
+    labels=["H_out", "D_out"],
+)
+plt.legend()
+plt.xlabel("Time (s)")
+plt.ylabel("Flux (atom/m^2/s)")
+plt.figure()
+plt.stackplot(
+    HD_flux.t,
+    np.abs(HH_flux.data),
+    np.abs(HD_flux.data),
+    np.abs(DD_flux.data),
+    labels=["HH", "HD", "DD"],
+)
+plt.legend(reverse=True)
+plt.xlabel("Time (s)")
+plt.ylabel("Flux (molecule/m^2/s)")
+
+
+plt.figure()
+plt.plot(H_flux_right.t, -np.array(H_flux_right.data), label="from gradient (H)")
+plt.plot(
+    H_flux_right.t,
+    2 * np.array(HH_flux.data) + np.array(HD_flux.data),
+    linestyle="--",
+    label="from reaction rates (H)",
+)
+
+plt.plot(D_flux_right.t, -np.array(D_flux_right.data), label="from gradient (D)")
+plt.plot(
+    D_flux_right.t,
+    2 * np.array(DD_flux.data) + np.array(HD_flux.data),
+    linestyle="--",
+    label="from reaction rates (D)",
+)
+plt.xlabel("Time (s)")
+plt.ylabel("Flux (atom/m^2/s)")
+plt.legend()
+plt.show()
+
+# check that H_flux_right == 2*HH_flux + HD_flux
+H_flux_from_gradient = -np.array(H_flux_right.data)
+H_flux_from_reac = 2 * np.array(HH_flux.data) + np.array(HD_flux.data)
+assert np.allclose(
+    H_flux_from_gradient,
+    H_flux_from_reac,
+    rtol=0.5e-2,
+    atol=0.005,
+)
+# check that D_flux_right == 2*DD_flux + HD_flux
+D_flux_from_gradient = -np.array(D_flux_right.data)
+D_flux_from_reac = 2 * np.array(DD_flux.data) + np.array(HD_flux.data)
+assert np.allclose(
+    D_flux_from_gradient,
+    D_flux_from_reac,
+    rtol=0.5e-2,
+    atol=0.005,
+)

src/festim/__init__.py

@@ -18,6 +18,9 @@
 )
 from .boundary_conditions.flux_bc import FluxBCBase, HeatFluxBC, ParticleFluxBC
 from .boundary_conditions.henrys_bc import HenrysBC
+from .boundary_conditions.flux_bc import FluxBCBase, ParticleFluxBC, HeatFluxBC
+from .boundary_conditions.surface_reaction import SurfaceReactionBC
+
 from .boundary_conditions.sieverts_bc import SievertsBC
 from .exports.average_surface import AverageSurface
 from .exports.average_volume import AverageVolume

src/festim/boundary_conditions/__init__.py

@@ -1,4 +1,11 @@
 from .dirichlet_bc import FixedConcentrationBC, FixedTemperatureBC
 from .flux_bc import HeatFluxBC, ParticleFluxBC
+from .surface_reaction import SurfaceReactionBC
 
-__all__ = ["FixedConcentrationBC", "FixedTemperatureBC", "ParticleFluxBC", "HeatFluxBC"]
+__all__ = [
+    "FixedConcentrationBC",
+    "FixedTemperatureBC",
+    "ParticleFluxBC",
+    "SurfaceReactionBC",
+    "HeatFluxBC",
+]

src/festim/boundary_conditions/surface_reaction.py

@@ -0,0 +1,127 @@
+from festim.boundary_conditions import ParticleFluxBC
+from festim import k_B
+from dolfinx import fem
+import ufl
+
+
+class SurfaceReactionBCpartial(ParticleFluxBC):
+    """Boundary condition representing a surface reaction
+    A + B <-> C
+    where A, B are the reactants and C is the product
+    the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T))
+    and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T))
+    The reaction rate is:
+    K = K_r * C_A * C_B - K_d * P_C
+    with C_A, C_B the concentration of species A and B,
+    P_C the partial pressure of species C at the surface.
+
+    This class is used to create the flux of a single species entering the surface
+    Example: The flux of species A entering the surface is K.
+
+
+    Args:
+        reactant (list): list of F.Species objects representing the reactants
+        gas_pressure (float, callable): the partial pressure of the product species
+        k_r0 (float): the pre-exponential factor of the forward reaction rate
+        E_kr (float): the activation energy of the forward reaction rate (eV)
+        k_d0 (float): the pre-exponential factor of the backward reaction rate
+        E_kd (float): the activation energy of the backward reaction rate (eV)
+        subdomain (F.SurfaceSubdomain): the surface subdomain where the reaction occurs
+        species (F.Species): the species to which the flux is applied
+    """
+
+    def __init__(
+        self,
+        reactant,
+        gas_pressure,
+        k_r0,
+        E_kr,
+        k_d0,
+        E_kd,
+        subdomain,
+        species,
+    ):
+        self.reactant = reactant
+        self.gas_pressure = gas_pressure
+        self.k_r0 = k_r0
+        self.E_kr = E_kr
+        self.k_d0 = k_d0
+        self.E_kd = E_kd
+        super().__init__(subdomain=subdomain, value=None, species=species)
+
+    def create_value_fenics(self, mesh, temperature, t: fem.Constant):
+        kr = self.k_r0 * ufl.exp(-self.E_kr / (k_B * temperature))
+        kd = self.k_d0 * ufl.exp(-self.E_kd / (k_B * temperature))
+        if callable(self.gas_pressure):
+            gas_pressure = self.gas_pressure(t=t)
+        else:
+            gas_pressure = self.gas_pressure
+
+        product_of_reactants = self.reactant[0].concentration
+        for reactant in self.reactant[1:]:
+            product_of_reactants *= reactant.concentration
+
+        self.value_fenics = kd * gas_pressure - kr * product_of_reactants
+
+
+class SurfaceReactionBC:
+    """Boundary condition representing a surface reaction
+    A + B <-> C
+    where A, B are the reactants and C is the product
+    the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T))
+    and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T))
+    The reaction rate is:
+    K = K_r * C_A * C_B - K_d * P_C
+    with C_A, C_B the concentration of species A and B,
+    P_C the partial pressure of species C at the surface.
+
+    The flux of species A entering the surface is K.
+    In the special case where A=B, then the flux of particle entering the surface is 2*K
+
+
+    Args:
+        reactant (list): list of F.Species objects representing the reactants
+        gas_pressure (float, callable): the partial pressure of the product species
+        k_r0 (float): the pre-exponential factor of the forward reaction rate
+        E_kr (float): the activation energy of the forward reaction rate (eV)
+        k_d0 (float): the pre-exponential factor of the backward reaction rate
+        E_kd (float): the activation energy of the backward reaction rate (eV)
+        subdomain (F.SurfaceSubdomain): the surface subdomain where the reaction occurs
+    """
+
+    def __init__(
+        self,
+        reactant,
+        gas_pressure,
+        k_r0,
+        E_kr,
+        k_d0,
+        E_kd,
+        subdomain,
+    ):
+        self.reactant = reactant
+        self.gas_pressure = gas_pressure
+        self.k_r0 = k_r0
+        self.E_kr = E_kr
+        self.k_d0 = k_d0
+        self.E_kd = E_kd
+        self.subdomain = subdomain
+
+        # create the flux boundary condition for each reactant
+        self.flux_bcs = [
+            SurfaceReactionBCpartial(
+                reactant=self.reactant,
+                gas_pressure=self.gas_pressure,
+                k_r0=self.k_r0,
+                E_kr=self.E_kr,
+                k_d0=self.k_d0,
+                E_kd=self.E_kd,
+                subdomain=self.subdomain,
+                species=species,
+            )
+            for species in self.reactant
+        ]
+
+    @property
+    def time_dependent(self):
+        return False  # no need to update if only using ufl.conditional objects

src/festim/hydrogen_transport_problem.py

@@ -508,10 +508,26 @@ def assign_functions_to_species(self):
             spe.test_function = sub_test_functions[idx]
 
     def define_boundary_conditions(self):
+        # @jhdark this all_bcs could be a property
+        # I just don't want to modify self.boundary_conditions
+
+        # create all_bcs which includes all flux bcs from SurfaceReactionBC
+        all_bcs = self.boundary_conditions.copy()
         for bc in self.boundary_conditions:
+            if isinstance(bc, boundary_conditions.SurfaceReactionBC):
+                all_bcs += bc.flux_bcs
+                all_bcs.remove(bc)
+
+        for bc in all_bcs:
             if isinstance(bc.species, str):
                 # if name of species is given then replace with species object
                 bc.species = _species.find_species_from_name(bc.species, self.species)
+            if isinstance(bc, boundary_conditions.ParticleFluxBC):
+                bc.create_value_fenics(
+                    mesh=self.mesh.mesh,
+                    temperature=self.temperature_fenics,
+                    t=self.t,
+                )
 
         super().define_boundary_conditions()
 
@@ -681,6 +697,13 @@ def create_formulation(self):
                     * bc.species.test_function
                     * self.ds(bc.subdomain.id)
                 )
+            if isinstance(bc, boundary_conditions.SurfaceReactionBC):
+                for flux_bc in bc.flux_bcs:
+                    self.formulation -= (
+                        flux_bc.value_fenics
+                        * flux_bc.species.test_function
+                        * self.ds(flux_bc.subdomain.id)
+                    )
 
         # check if each species is defined in all volumes
         if not self.settings.transient:
@@ -714,8 +737,15 @@ def update_time_dependent_values(self):
             self.temperature_fenics.interpolate(self.temperature_expr)
 
         for bc in self.boundary_conditions:
-            if bc.temperature_dependent:
-                bc.update(t=t)
+            if isinstance(
+                bc,
+                (
+                    boundary_conditions.FixedConcentrationBC,
+                    boundary_conditions.ParticleFluxBC,
+                ),
+            ):
+                if bc.temperature_dependent:
+                    bc.update(t=t)
 
         for source in self.sources:
             if source.temperature_dependent: