Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Theory and UG for the SurfaceKinetics BC #808

Merged
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 56 additions & 2 deletions docs/source/bibliography/references.bib
Original file line number Diff line number Diff line change
Expand Up @@ -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}},
Expand Down Expand Up @@ -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}},
Expand All @@ -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}
}
Binary file added docs/source/images/potential_diagram.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
94 changes: 82 additions & 12 deletions docs/source/theory.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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:

Expand All @@ -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
Expand All @@ -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:

Expand All @@ -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:

Expand All @@ -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
Expand All @@ -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.

Heat transfer
^^^^^^^^^^^^^^
Expand All @@ -122,7 +125,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
Expand Down Expand Up @@ -272,9 +276,75 @@ 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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Vladimir, do you not want to also mention your papers ? I think you implemented and used this kind of model in you laser desorption papers (but it is up to you).

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ehodille, well, I used, yes. However, I intended to mention here the pioneers, whereas my papers on LID were mainly based on your works.

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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

does J_vs contains the desorbing flux and all abstraction fluxes (ELy Rideal recombination etc) ?

Copy link
Collaborator Author

@KulaginVladimir KulaginVladimir Jul 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In general, yes, but it should be defined by a user. Here is an example from one validation case on how it can be defined:

def J_vs(T, surf_conc, t):
    phi_atom = SP * Gamma_atom * (1 - surf_conc / n_surf) # adsorption
    phi_exc = Gamma_atom * sigma_exc * surf_conc # abstraction
    phi_des = 2 * nu0 * (lambda_des * surf_conc) ** 2 * f.exp(-2 * E_des / F.k_B / T) # desorption
    return phi_atom - phi_exc - phi_des

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe this could be explicited then? something like

is the net flux of hydrogen atoms from the vacuum onto the surface (eg recombination, dissociation, Eley-Rideal recombination, etc.).

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this case, should we somehow clarify that the flux comes onto the surface: J_vs = +Dissoc. - Recomb. - Eley-Rideal recomb ...?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe just say J_vs = J_in - J_out where J_in is the sum of all fluxes coming from the vacuum to the surface and J_out from the surface to the vacuum

Or something like that and then give a few examples

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few examples on different processes that can be included?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes just like you had in your previous comment

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I updated this part. Let me know what you think.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me!

atoms from the vacuum onto the surface. 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
35 changes: 35 additions & 0 deletions docs/source/userguide/boundary_conditions.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
^^^^^^^^^^^^^^^^^^^^^^^^^^^
Expand Down
6 changes: 5 additions & 1 deletion docs/source/userguide/export_post_processing.rst
Original file line number Diff line number Diff line change
Expand Up @@ -167,7 +167,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`.

The data can be accessed in three different ways:
- directly from the :class:`festim.DerivedQuantities` (plural) object:
Expand Down
Loading