Surface reactions

festim/__init__.py

@@ -27,6 +27,7 @@
 from .boundary_conditions.sieverts_bc import SievertsBC
 from .boundary_conditions.henrys_bc import HenrysBC
 from .boundary_conditions.flux_bc import FluxBCBase, ParticleFluxBC, HeatFluxBC
+from .boundary_conditions.surface_reaction import SurfaceReactionBC
 
 from .material import Material
 

festim/boundary_conditions/surface_reaction.py

@@ -0,0 +1,127 @@
+import festim as F
+from dolfinx import fem
+import ufl
+
+
+class SurfaceReactionBCpartial(F.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 / (F.k_B * temperature))
+        kd = self.k_d0 * ufl.exp(-self.E_kd / (F.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

festim/hydrogen_transport_problem.py

@@ -248,7 +248,6 @@
         self.define_temperature()
         self.define_boundary_conditions()
         self.create_source_values_fenics()
-        self.create_flux_values_fenics()
         self.create_initial_conditions()
         self.create_formulation()
         self.create_solver()
@@ -512,14 +511,30 @@
         )
 
     def define_boundary_conditions(self):
-        """Defines the dirichlet boundary conditions of the model"""
+        """Create forms for DirichletBC and value_fenics for ParticleFluxBC"""
+        # @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, F.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 = F.find_species_from_name(bc.species, self.species)
             if isinstance(bc, F.DirichletBC):
                 form = self.create_dirichletbc_form(bc)
                 self.bc_forms.append(form)
+            if isinstance(bc, F.ParticleFluxBC):
+                bc.create_value_fenics(
+                    mesh=self.mesh.mesh,
+                    temperature=self.temperature_fenics,
+                    t=self.t,
+                )
 
     def create_dirichletbc_form(self, bc):
         """Creates a dirichlet boundary condition form
@@ -590,18 +605,6 @@
                     t=self.t,
                 )
 
-    def create_flux_values_fenics(self):
-        """For each particle flux create the value_fenics"""
-        for bc in self.boundary_conditions:
-            # create value_fenics for all F.ParticleFluxBC objects
-            if isinstance(bc, F.ParticleFluxBC):
-
-                bc.create_value_fenics(
-                    mesh=self.mesh.mesh,
-                    temperature=self.temperature_fenics,
-                    t=self.t,
-                )
-
     def create_initial_conditions(self):
         """For each initial condition, create the value_fenics and assign it to
         the previous solution of the condition's species"""
@@ -701,6 +704,13 @@
                     * bc.species.test_function
                     * self.ds(bc.subdomain.id)
                 )
+            if isinstance(bc, F.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:

surface_reaction_example.py

@@ -0,0 +1,52 @@
+import festim as F
+import numpy as np
+import ufl
+
+
+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 = 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_manual = F.ParticleFluxBC(
+    subdomain=right,
+    value=lambda c1, T: -2 * 0.01 * ufl.exp(-0 / F.k_B / T) * c1**2,
+    species=H,
+    species_dependent_value={"c1": H},
+)
+
+my_model.boundary_conditions = [
+    F.DirichletBC(subdomain=left, value=5, species=H),
+    F.DirichletBC(subdomain=left, value=5, species=D),
+    surface_reaction,
+    # surface_reaction_manual,
+]
+
+my_model.exports = [F.XDMFExport("test.xdmf", H)]
+
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=10, transient=True)
+
+my_model.settings.stepsize = 0.1
+
+my_model.initialise()
+my_model.run()

test/test_h_transport_problem.py

@@ -964,7 +964,7 @@ def test_update_time_dependent_values_flux(bc_value, expected_values):
     my_model.define_meshtags_and_measures()
     my_model.assign_functions_to_species()
     my_model.define_temperature()
-    my_model.create_flux_values_fenics()
+    my_model.define_boundary_conditions()
 
     for i in range(3):
         # RUN
@@ -1018,7 +1018,7 @@ def test_update_fluxes_with_time_dependent_temperature(
     my_model.define_function_spaces()
     my_model.assign_functions_to_species()
     my_model.define_meshtags_and_measures()
-    my_model.create_flux_values_fenics()
+    my_model.define_boundary_conditions()
 
     for i in range(3):
         # RUN
@@ -1031,8 +1031,8 @@ def test_update_fluxes_with_time_dependent_temperature(
             assert np.isclose(computed_value, expected_values[i])
 
 
-def test_create_source_values_fenics_multispecies():
-    """Test that the create_flux_values_fenics method correctly sets the value_fenics
+def test_define_boundary_conditions_multispecies():
+    """Test that the define_boundary_conditions method correctly sets the value_fenics
     attribute in a multispecies case"""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
@@ -1056,7 +1056,7 @@ def test_create_source_values_fenics_multispecies():
     my_model.define_temperature()
 
     # RUN
-    my_model.create_flux_values_fenics()
+    my_model.define_boundary_conditions()
 
     # TEST
     assert np.isclose(my_model.boundary_conditions[0].value_fenics.value, 5)