Merge pull request #808 from KulaginVladimir/SurfaceKinetics_docs

Theory and UG for the `SurfaceKinetics` BC

docs/source/bibliography/references.bib

@@ -4,7 +4,7 @@ @article{McNabb1963
 	year         = 1963,
 	journal      = {Trans. Metall. Soc. AIME},
 	volume       = 227,
-	pages        = 618,
+	pages        = 618
 }
 @article{Longhurst1985,
 	title        = {{The soret effect and its implications for fusion reactors}},
@@ -33,7 +33,7 @@ @article{Delaporte-Mathurin2021
 	number       = 3,
 	pages        = {036038},
 	issn         = {0029-5515},
-	url          = {https://iopscience.iop.org/article/10.1088/1741-4326/abd95f},
+	url          = {https://iopscience.iop.org/article/10.1088/1741-4326/abd95f}
 }
 @phdthesis{Delaporte-Mathurin2022,
 	title        = {{Hydrogen transport in tokamaks : Estimation of the ITER divertor tritium inventory and influence of helium exposure}},
@@ -54,4 +54,58 @@ @article{Schmid2016
 	pages        = {014025},
 	issn         = {0031-8949},
 	url          = {https://iopscience.iop.org/article/10.1088/0031-8949/T167/1/014025}
+}
+@article{Guterl2019,
+	title        = {Effects of surface processes on hydrogen outgassing from metal in desorption experiments},
+	author       = {Guterl, Jerome and Smirnov, RD and Snyder, P},
+	year         = 2019,
+	journal      = {Nuclear Fusion},
+	publisher    = {IOP Publishing},
+	volume       = 59,
+	number       = 9,
+	pages        = {096042},
+	url          = {https://iopscience.iop.org/article/10.1088/1741-4326/ab280a/meta}
+}
+@article{Pick1985,
+	title        = {A model for atomic hydrogen-metal interactions—application to recycling, recombination and permeation},
+	author       = {Pick, MA and Sonnenberg, K},
+	year         = 1985,
+	journal      = {Journal of Nuclear Materials},
+	publisher    = {Elsevier},
+	volume       = 131,
+	number       = {2-3},
+	pages        = {208--220},
+	url          = {https://www.sciencedirect.com/science/article/abs/pii/0022311585904593}
+}
+@article{Hodille2017,
+	title        = {Simulations of atomic deuterium exposure in self-damaged tungsten},
+	author       = {Hodille, EA and Zalo{\v{z}}nik, A and Markelj, S and Schwarz-Selinger, T and Becquart, CS and Bisson, R{\'e}gis and Grisolia, Christian},
+	year         = 2017,
+	journal      = {Nuclear Fusion},
+	publisher    = {IOP Publishing},
+	volume       = 57,
+	number       = 5,
+	pages        = {056002},
+	url          = {https://iopscience.iop.org/article/10.1088/1741-4326/aa5aa5/meta}
+}
+@article{Schmid2021,
+	title        = {On the use of recombination rate coefficients in hydrogen transport calculations},
+	author       = {Schmid, K and Zibrov, M},
+	year         = 2021,
+	journal      = {Nuclear Fusion},
+	publisher    = {IOP Publishing},
+	volume       = 61,
+	number       = 8,
+	pages        = {086008},
+	url          = {https://iopscience.iop.org/article/10.1088/1741-4326/ac07b2/meta}
+}
+@article{Hamamoto2020,
+	title        = {Comprehensive modeling of hydrogen transport and accumulation in titanium and zirconium},
+	author       = {Hamamoto, Yoshiki and Uchikoshi, Takeru and Tanabe, Katsuaki},
+	year         = 2020,
+	journal      = {Nuclear Materials and Energy},
+	publisher    = {Elsevier},
+	volume       = 23,
+	pages        = 100751,
+	url          = {https://www.sciencedirect.com/science/article/pii/S2352179120300272}
 }

docs/source/theory.rst

@@ -8,14 +8,17 @@ Bulk physics
 
 H transport
 ^^^^^^^^^^^
-The model developed by McNabb & Foster :cite:`McNabb1963` is used to model hydrogen transport in materials in FESTIM. The principle is to separate mobile hydrogen :math:`c_\mathrm{m}` and trapped hydrogen :math:`c_\mathrm{t}`. The diffusion of mobile particles is governed by Fick’s law of diffusion where the hydrogen flux is
+The model developed by McNabb & Foster :cite:`McNabb1963` is used to model hydrogen transport in materials in FESTIM. 
+The principle is to separate mobile hydrogen :math:`c_\mathrm{m}\,[\mathrm{m}^{-3}]` and trapped hydrogen :math:`c_\mathrm{t}\,[\mathrm{m}^{-3}]`. 
+The diffusion of mobile particles is governed by Fick’s law of diffusion where the hydrogen flux is
 
 .. math::
     :label: eq_difflux
     
     J = -D \nabla c_\mathrm{m}
 
-where :math:`D=D(T)` is the diffusivity. Each trap :math:`i` is associated with a trapping and a detrapping rate :math:`k_i` and :math:`p_i`, respectively, as well as a trap density :math:`n_i`.
+where :math:`D=D(T)\,[\mathrm{m}^{2}\,\mathrm{s}^{-1}]` is the diffusivity. 
+Each trap :math:`i` is associated with a trapping and a detrapping rate :math:`k_i\,[\mathrm{m}^{3}\,\mathrm{s}^{-1}]` and :math:`p_i\,[\mathrm{s}^{-1}]`, respectively, as well as a trap density :math:`n_i\,[\mathrm{m}^{-3}]`.
 
 The temporal evolution of :math:`c_\mathrm{m}` and :math:`c_{\mathrm{t}, i}` are then given by:
 
@@ -29,7 +32,7 @@ The temporal evolution of :math:`c_\mathrm{m}` and :math:`c_{\mathrm{t}, i}` are
 
     \frac {\partial c_{\mathrm{t}, i}} { \partial t} = k_i c_\mathrm{m} (n_i - c_{\mathrm{t},i}) - p_i c_{\mathrm{t},i}
 
-where :math:`S_j=S_j(x,y,z,t)` is a source :math:`j` of mobile hydrogen. In FESTIM, source terms can be space and time dependent. These are used to simulate plasma implantation in materials, tritium generation from neutron interactions, etc. 
+where :math:`S_j=S_j(x,y,z,t)\,[\mathrm{m}^{-3}\,\mathrm{s}^{-1}]` is a source :math:`j` of mobile hydrogen. In FESTIM, source terms can be space and time dependent. These are used to simulate plasma implantation in materials, tritium generation from neutron interactions, etc. 
 These equations can be solved in cartesian coordinates but also in cylindrical and spherical coordinates. This is useful, for instance, when simulating hydrogen transport in a pipe or in a pebble. FESTIM can solve steady-state hydrogen transport problems.
 
 Soret effect
@@ -41,18 +44,18 @@ FESTIM can include the Soret effect :cite:`Pendergrass1976,Longhurst1985` (also
 
     J = -D \nabla c_\mathrm{m} - D\frac{Q^* c_\mathrm{m}}{k_B T^2} \nabla T
 
-where :math:`Q^*` is the Soret coefficient (also called heat of transport) and :math:`k_B` is the Boltzmann constant.
+where :math:`Q^*\,[\mathrm{eV}]` is the Soret coefficient (also called heat of transport) and :math:`k_B` is the Boltzmann constant.
 
 Conservation of chemical potential at interfaces
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-Continuity of local partial pressure :math:`P` at interfaces between materials has to be ensured. In the case of a material behaving according to Sievert’s law of solubility (metals), the partial pressure is expressed as:
+Continuity of local partial pressure :math:`P\,[\mathrm{Pa}]` at interfaces between materials has to be ensured. In the case of a material behaving according to Sievert’s law of solubility (metals), the partial pressure is expressed as:
 
 .. math::
     :label: eq_Sievert   
 
     P = \left(\frac{c_\mathrm{m}}{K_S}\right)^2
 
-where :math:`K_S` is the material solubility (or Sivert's constant).
+where :math:`K_S\,[\mathrm{m}^{-3}\,\mathrm{Pa}^{-0.5}]` is the material solubility (or Sivert's constant).
 
 In the case of a material behaving according to Henry's law of solubility, the partial pressure is expressed as:
 
@@ -61,7 +64,7 @@ In the case of a material behaving according to Henry's law of solubility, the p
 
     P = \frac{c_\mathrm{m}}{K_H}
 
-where :math:`K_H` is the material solubility (or Henry's constant).
+where :math:`K_H\,[\mathrm{m}^{-3}\,\mathrm{Pa}^{-1}]` is the material solubility (or Henry's constant).
 
 Two different interface cases can then occur. At the interface between two Sievert or two Henry materials, the continuity of partial pressure yields:
 
@@ -85,7 +88,7 @@ At the interface between a Sievert and a Henry material:
 
 It appears from these equilibrium equations that a difference in solubilities introduces a concentration jump at interfaces.
 
-In FESTIM, the conservation of chemical potential is obtained by a change of variables :cite:`Delaporte-Mathurin2021`. The variable :math:`\theta` is introduced and:
+In FESTIM, the conservation of chemical potential is obtained by a change of variables :cite:`Delaporte-Mathurin2021`. The variable :math:`\theta\,[\mathrm{Pa}]` is introduced and:
 
 .. math::
     :label: eq_theta
@@ -111,7 +114,7 @@ According to the latter, the rate :math:`k(T)` of a thermally activated process
 
     k(T) = k_0 \exp \left[-\frac{E_k}{k_B T} \right]
 
-where :math:`k_0` is the pre-exponential factor, :math:`E_k` is the process activation energy, and :math:`T` is the temperature. 
+where :math:`k_0` is the pre-exponential factor, :math:`E_k\,[\mathrm{eV}]` is the process activation energy, and :math:`T\,[\mathrm{K}]` is the temperature. 
 
 Radioactive decay
 ^^^^^^^^^^^^^^^^^
@@ -135,7 +138,8 @@ To properly account for the temperature-dependent parameters, an accurate repres
 
     \rho C_p \frac{\partial T}{\partial t} = \nabla \cdot (\lambda \nabla T) + \sum_i Q_i
 
-where :math:`T` is the temperature, :math:`C_p` is the specific heat capacity, :math:`\rho` is the material's density, :math:`\lambda` is the thermal conductivity and :math:`Q_i` is a volumetric heat source :math:`i`. As for the hydrogen transport problem, the heat equation can be solved in steady state. In FESTIM, the thermal properties of materials can be arbitrary functions of temperature.
+where :math:`T` is the temperature, :math:`C_p\,[\mathrm{J}\,\mathrm{kg}^{-1}\,\mathrm{K}^{-1}]` is the specific heat capacity, :math:`\rho\,[\mathrm{kg}\,\mathrm{m}^{-3}]` is the material's density, 
+:math:`\lambda\,[\mathrm{W}\,\mathrm{m}^{-1}\,\mathrm{K}^{-1}]` is the thermal conductivity and :math:`Q_i\,[\mathrm{W}\,\mathrm{m}^{-3}]` is a volumetric heat source :math:`i`. As for the hydrogen transport problem, the heat equation can be solved in steady state. In FESTIM, the thermal properties of materials can be arbitrary functions of temperature.
 
 ---------------
 Surface physics 
@@ -285,9 +289,79 @@ Finally, convective heat fluxes can be applied to boundaries:
 
 where :math:`h` is the heat transfer coefficient and :math:`T_{\mathrm{ext}}` is the external temperature.
 
----------------
+Kinetic surface model
+^^^^^^^^^^^^^^^^^^^^^
+
+Modelling hydrogen retention or outgassing might require considering the kinetics of surface processes. 
+A representative example is the hydrogen uptake from a gas phase, when the energy of incident atoms/molecules is not high enough to 
+overcome the surface barrier for implantation. The general approach to account for surface kinetics :cite:`Pick1985, Hodille2017, Guterl2019, Schmid2021` consists in 
+introducing hydrogen surface species :math:`c_\mathrm{s}\,[\mathrm{m}^{-2}]`. 
+
+Evolution of hydrogen surface concentration is determined by the atomic flux balance at the surface, as sketched in the simplified energy diagram below.
+
+.. thumbnail:: images/potential_diagram.png
+    :align: center
+    :width: 800
+
+    Idealised potential energy diagram for hydrogen near a surface of an endothermic metal. Energy levels are measured from the :math:`\mathrm{H}_2` state
+
+The governing equation for surface species is:
+
+.. math::
+    :label: eq_surf_conc
+    
+    \dfrac{d c_\mathrm{s}}{d t} = J_\mathrm{bs} - J_\mathrm{sb} + J_\mathrm{vs}~\text{on}~\delta\Omega
+
+where :math:`J_\mathrm{bs}\,[\mathrm{m}^{-2}\,\mathrm{s}^{-1}]` is the flux of hydrogen atoms from the subsurface (bulk region just beneath the surface) onto the surface, 
+:math:`J_\mathrm{sb}\,[\mathrm{m}^{-2}\,\mathrm{s}^{-1}]` is the flux of hydrogen atoms from the surface into the subsurface, and :math:`J_\mathrm{vs}\,[\mathrm{m}^{-2}\,\mathrm{s}^{-1}]` 
+is the net flux of hydrogen atoms from the vaccuum onto the surface. The latter is defined as :math:`J_\mathrm{vs}=J_\mathrm{in}-J_\mathrm{out}`, where :math:`J_\mathrm{in}\,[\mathrm{m}^{-2}\,\mathrm{s}^{-1}]` 
+is the sum of all fluxes coming from the vacuum onto the surface and :math:`J_\mathrm{out}\,[\mathrm{m}^{-2}\,\mathrm{s}^{-1}]` is the sum of all fluxes coming from the surface to the vacuum.
+:math:`J_\mathrm{in}` can be used to set up adsorption fluxes from different processes such as molecular dissociation, adsorption of low-energy atoms, etc. Similarly, 
+:math:`J_\mathrm{out}` can be used to define desorption fluxes from various processes such as Langmuir-Hinshelwood recombination, Eley-Rideal recombination, sputtering, etc.
+It worth noticing that the current model does not account for possible surface diffusion and is available only for 1D hydrogen transport simulations.
+
+The connection condition between surface and bulk domains represents the Robin boundary condition for the hydrogen transport problem.
+
+.. math::
+    :label: eq_subsurf_conc
+    
+    -D \nabla c_\mathrm{m} \cdot \mathbf{n} = \lambda_{\mathrm{IS}} \dfrac{\partial c_{\mathrm{m}}}{\partial t} + J_{\mathrm{bs}} - J_{\mathrm{sb}}~\text{on}~\delta\Omega
+
+where :math:`\lambda_\mathrm{IS}\,[\mathrm{m}]` is the distance between two interstitial sites in the bulk. 
+
+.. note::
+
+    At the left boundary, the normal vector :math:`\textbf{n}` is :math:`-\vec{x}`. The steady-state approximation of eq. :eq:`eq_subsurf_conc` at the left boundary 
+    is, therefore, :math:`D\frac{\partial c_\mathrm{m}}{\partial x}=J_\mathrm{bs}-J_\mathrm{sb}` representing eq. (12) in the original work of M.A. Pick & K. Sonnenberg :cite:`Pick1985`.
+
+The fluxes for subsurface-to-surface and surface-to-subsurface transitions are defined as follows:
+
+.. 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)
+
+.. 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}]` 
+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.
+
+.. note::
+
+    In eq. :eq:`eq_Jsb`, the last term in brackets is usually omitted :cite:`Guterl2019, Pick1985, Hodille2017, Schmid2021`, 
+    since :math:`c_\mathrm{m} \ll n_\mathrm{IS}` is assumed. However, this term is included in some works (e.g. :cite:`Hamamoto2020`) 
+    to better reproduce the experimental results.   
+
+
+------------
 References
----------------
+------------
 
 .. bibliography:: bibliography/references.bib
     :style: unsrt

docs/source/userguide/boundary_conditions.rst

@@ -105,6 +105,41 @@ Where :math:`Kd` is the dissociation coefficient, :math:`P` is the partial press
 
     my_bc = DissociationFlux(surfaces=2, Kd_0=2, E_Kd=0.1, P=1e05)
 
+Kinetic surface model (1D)
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Kinetic surface model can be included to account for the evolution of adsorbed hydrogen on a surface with the :class:`festim.SurfaceKinetics` class.
+The current class is supported for 1D simulations only. Refer to the :ref:`Kinetic surface model` theory section for more details.
+
+.. testcode:: BCs
+
+    from festim import t
+    import fenics as f
+
+    def k_bs(T, surf_conc, t):
+        return 1e13*f.exp(-0.2/k_b/T)
+
+    def k_sb(T, surf_conc, t):
+        return 1e13*f.exp(-1.0/k_b/T)
+
+    def J_vs(T, surf_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)
+
+        return J_ads - J_des
+
+    my_bc = SurfaceKinetics(
+        k_bs=k_bs,
+        k_sb=k_sb,
+        lambda_IS=1.1e-10,
+        n_surf=1e17,
+        n_IS=6.3e28,
+        J_vs=J_vs,
+        surfaces=3,
+        initial_condition=0,
+        t=t
+        )
 
 Sievert's law of solubility
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^

docs/source/userguide/export_post_processing.rst

@@ -168,7 +168,11 @@ Finally, you can add the :class:`festim.DerivedQuantities` object to the simulat
     my_model.exports = [my_derived_quantities]
 
 
-The complete list of derived quantities can be found at: :ref:`Exports`.
+The complete list of derived quantities can be found at: :ref:`Exports`. 
+
+.. note::
+
+    There is a specific derived quantity :class:`festim.AdsorbedHydrogen` which can be used only with :class:`festim.SurfaceKinetics`.
 
 ^^^^^^^^^^^^^^^^^^
 Accessing the data