Merge pull request #941 from KulaginVladimir/surfacekinetics_refactoring

`SurfaceKinetics` minor update

docs/source/theory.rst

@@ -339,15 +339,14 @@ The fluxes for subsurface-to-surface and surface-to-subsurface transitions are d
 .. math::
     :label: eq_Jbs
 
-    J_\mathrm{bs} = k_\mathrm{bs} \lambda_\mathrm{abs} c_\mathrm{m} \left(1-\dfrac{c_\mathrm{s}}{n_\mathrm{surf}}\right)
+    J_\mathrm{bs} = k_\mathrm{bs} c_\mathrm{m} \left(1-\dfrac{c_\mathrm{s}}{n_\mathrm{surf}}\right)
 
 .. math::
     :label: eq_Jsb
 
     J_\mathrm{sb} = k_\mathrm{sb} c_\mathrm{s} \left(1-\dfrac{c_\mathrm{m}}{n_\mathrm{IS}}\right)
 
-where :math:`n_\mathrm{surf}\,[\mathrm{m}^{-2}]` is the surface concentration of adsorption sites, :math:`n_\mathrm{IS}\,[\mathrm{m}^{-3}]` is the bulk concentration of interstitial sites,
-:math:`\lambda_\mathrm{abs}=n_\mathrm{surf}/n_\mathrm{IS}\,[\mathrm{m}]` is the characteristic distance between surface and subsurface sites, :math:`k_\mathrm{bs}\,[\mathrm{s}^{-1}]` 
+where :math:`n_\mathrm{surf}\,[\mathrm{m}^{-2}]` is the surface concentration of adsorption sites, :math:`n_\mathrm{IS}\,[\mathrm{m}^{-3}]` is the bulk concentration of interstitial sites, :math:`k_\mathrm{bs}\,[\mathrm{m}\,\mathrm{s}^{-1}]` 
 and :math:`k_\mathrm{sb}\,[\mathrm{s}^{-1}]` are the rate constants for subsurface-to-surface and surface-to-subsurface transitions, respectively. 
 Usually, these rate constants are expressed in the Arrhenius form: :math:`k_i=k_{i,0}\exp(-E_{k,i} / kT)`. Both these processes are assumed to take place
 if there are available sites on the surface (in the subsurface). Possible surface/subsurface saturation is accounted for with terms in brackets.

docs/source/userguide/boundary_conditions.rst

@@ -116,13 +116,13 @@ The current class is supported for 1D simulations only. Refer to the :ref:`Kinet
     from festim import t
     import fenics as f
 
-    def k_bs(T, surf_conc, t):
+    def k_bs(T, surf_conc, mobile_conc, t):
         return 1e13*f.exp(-0.2/k_b/T)
 
-    def k_sb(T, surf_conc, t):
+    def k_sb(T, surf_conc, mobile_conc, t):
         return 1e13*f.exp(-1.0/k_b/T)
 
-    def J_vs(T, surf_conc, t):
+    def J_vs(T, surf_conc, mobile_conc, t):
 
         J_des = 2e5*surf_conc**2*f.exp(-1.2/k_b/T)
         J_ads = 1e17*(1-surf_conc/1e17)**2*f.conditional(t<10, 1, 0)

festim/boundary_conditions/fluxes/surface_kinetics.py

@@ -22,17 +22,14 @@ class SurfaceKinetics(FluxBC):
 
     .. math::
 
-         J_{\mathrm{bs}} = k_{\mathrm{bs}} c_{\mathrm{m}} \lambda_{\mathrm{abs}} \left(1 - \dfrac{c_\mathrm{s}}{n_{\mathrm{surf}}}\right),
+         J_{\mathrm{bs}} = k_{\mathrm{bs}} c_{\mathrm{m}} \left(1 - \dfrac{c_\mathrm{s}}{n_{\mathrm{surf}}}\right),
 
     the H flux from surface to subsurface :math:`J_{\mathrm{sb}}` (in :math:`\mathrm{m}^{-2} \ \mathrm{s}^{-1}`) is:
 
     .. math::
 
          J_{\mathrm{sb}} = k_{\mathrm{sb}} c_{\mathrm{s}} \left(1 - \dfrac{c_{\mathrm{m}}}{n_\mathrm{IS}}\right),
 
-    :math:`\lambda_{\mathrm{abs}}=n_{\mathrm{surf}}/n_{\mathrm{IS}}` is the characteristic distance between surface and
-    subsurface sites (:math:`\mathrm{m}`).
-
     For more details see:
         E.A. Hodille et al 2017 Nucl. Fusion 57 056002; Y. Hamamoto et al 2020 Nucl. Mater. Energy 23 100751
 
@@ -43,7 +40,7 @@ class SurfaceKinetics(FluxBC):
     Args:
         k_sb (float or callable): rate constant for the surface-to-subsurface transition (:math:`\mathrm{s}^{-1}`),
             can accept additional parameters (see example)
-        k_bs (float or callable): rate constant for the subsurface-to-surface transition (:math:`\mathrm{s}^{-1}`),
+        k_bs (float or callable): rate constant for the subsurface-to-surface transition (:math:`\mathrm{m} \ \mathrm{s}^{-1}`),
             can accept additional parameters (see example)
         lambda_IS (float): characteristic distance between two iterstitial sites (:math:`\mathrm{m}`)
         n_surf (float): surface concentration of adsorption sites (:math:`\mathrm{m}^{-2}`)
@@ -63,13 +60,13 @@ class SurfaceKinetics(FluxBC):
 
     Example::
 
-        def K_sb(T, surf_conc, prm1, prm2):
-            return 1e13 * f.exp(-2.0/F.k_B/T)
+        def K_sb(T, surf_conc, mobile_conc, prm1, prm2):
+            return 1e13 * f.exp(-2.0/F.k_B/T) + mobile_conc
 
-        def K_bs(T, surf_conc, prm1, prm2):
+        def K_bs(T, surf_conc, mobile_conc, prm1, prm2):
             return 1e13 * f.exp(-0.2/F.k_B/T)
 
-        def J_vs(T, surf_conc, prm1, prm2):
+        def J_vs(T, surf_conc, mobile_conc, prm1, prm2):
             return (1-surf_conc/5) ** 2 * fenics.exp(-2/F.k_B/T) + prm1 * prm2
 
         my_surf_model = SurfaceKinetics(
@@ -133,25 +130,22 @@ def create_form(self, solute, solute_prev, solute_test_function, T, ds, dt):
         lambda_IS = self.lambda_IS
         n_surf = self.n_surf
         n_IS = self.n_IS
-        lambda_abs = (
-            n_surf / n_IS
-        )  # characteristic distance between surface and subsurface sites
         self.form = 0
 
         for i, surf in enumerate(self.surfaces):
 
             J_vs = self.J_vs
             if callable(J_vs):
-                J_vs = J_vs(T.T, self.solutions[i], **self.prms)
+                J_vs = J_vs(T.T, self.solutions[i], solute, **self.prms)
             k_sb = self.k_sb
             if callable(k_sb):
-                k_sb = k_sb(T.T, self.solutions[i], **self.prms)
+                k_sb = k_sb(T.T, self.solutions[i], solute, **self.prms)
             k_bs = self.k_bs
             if callable(k_bs):
-                k_bs = k_bs(T.T, self.solutions[i], **self.prms)
+                k_bs = k_bs(T.T, self.solutions[i], solute, **self.prms)
 
             J_sb = k_sb * self.solutions[i] * (1 - solute / n_IS)
-            J_bs = k_bs * (solute * lambda_abs) * (1 - self.solutions[i] / n_surf)
+            J_bs = k_bs * solute * (1 - self.solutions[i] / n_surf)
 
             if dt is not None:
                 # Surface concentration form

test/system/test_system.py

@@ -1226,23 +1226,25 @@ def test_MMS_surface_kinetics():
     """
     MMS test for SurfaceKinetics BC
     """
-    exact_solution_cm = lambda x, t: 1 + 2 * x**2 + x + 2 * t
-    exact_solution_cs = (
-        lambda t: n_surf * (1 + 2 * t + 2 * lambda_IS - D) / (2 * n_IS - 1 - 2 * t)
-    )
 
     n_IS = 20
     n_surf = 5
     D = 7
     lambda_IS = 2
-    k_bs = n_IS / n_surf
+    k_bs = 3
     k_sb = 2 * n_IS / n_surf
 
+    exact_solution_cm = lambda x, t: 1 + 2 * x**2 + x + 2 * t
+    exact_solution_cs = (
+        lambda t: n_surf
+        * (3 * (1 + 2 * t) + 2 * lambda_IS - D)
+        / (2 * n_IS + 1 + 2 * t)
+    )
     solute_source = 2 * (1 - 2 * D)
 
-    def J_vs(T, surf_conc, t):
+    def J_vs(T, surf_conc, solute, t):
         return (
-            2 * n_surf * (2 * n_IS + 2 * lambda_IS - D) / (2 * n_IS - 1 - 2 * t) ** 2
+            2 * n_surf * (6 * n_IS - 2 * lambda_IS + D) / (2 * n_IS + 1 + 2 * t) ** 2
             + 2 * lambda_IS
             - D
         )

test/unit/test_boundary_conditions.py

@@ -620,14 +620,14 @@ def test_create_form_surf_kinetics():
     """
 
     # build
-    def k_sb(cs, T, prm1, prm2):
-        return 2 * T + cs**2 + prm1 - prm2
+    def k_sb(T, cs, cm, prm1, prm2):
+        return 2 * T + cs**2 / cm + prm1 - prm2
 
-    def k_bs(cs, T, prm1, prm2):
-        return 2 * T + 3 * cs + prm1 + prm2
+    def k_bs(T, cs, cm, prm1, prm2):
+        return 2 * T + 3 * cs + cm + prm1 + prm2
 
-    def J_vs(cs, T, prm1, prm2):
-        return 2 * T + 1
+    def J_vs(T, cs, cm, prm1, prm2):
+        return 2 * T + 5 * cm - 3 * cs
 
     lambda_IS = 1
     n_surf = 1
@@ -678,11 +678,11 @@ def J_vs(cs, T, prm1, prm2):
     # test
     p1 = my_bc.sub_expressions[0]
     p2 = my_bc.sub_expressions[1]
-    K_sb = k_sb(T.T, adsorbed, p1, p2)
-    K_bs = k_bs(T.T, adsorbed, p1, p2)
-    j_vs = J_vs(T.T, adsorbed, p1, p2)
+    K_sb = k_sb(T.T, adsorbed, solute, p1, p2)
+    K_bs = k_bs(T.T, adsorbed, solute, p1, p2)
+    j_vs = J_vs(T.T, adsorbed, solute, p1, p2)
     J_sb = K_sb * adsorbed * (1 - solute / n_IS)
-    J_bs = K_bs * (solute * n_surf / n_IS) * (1 - adsorbed / n_surf)
+    J_bs = K_bs * solute * (1 - adsorbed / n_surf)
 
     expected_form = (
         (adsorbed - adsorbed_prev) / dt.value * adsorbed_test_function * ds(1)