Theory

Bulk physics

H transport

The model developed by McNabb & Foster [1] is used to model hydrogen transport in materials in FESTIM. The principle is to separate mobile hydrogen \(c_\mathrm{m}\,[\mathrm{m}^{-3}]\) and trapped hydrogen \(c_\mathrm{t}\,[\mathrm{m}^{-3}]\). The diffusion of mobile particles is governed by Fick’s law of diffusion where the hydrogen flux is

(1)\[J = -D \nabla c_\mathrm{m}\]

where \(D=D(T)\,[\mathrm{m}^{2}\,\mathrm{s}^{-1}]\) is the diffusivity. Each trap \(i\) is associated with a trapping and a detrapping rate \(k_i\,[\mathrm{m}^{3}\,\mathrm{s}^{-1}]\) and \(p_i\,[\mathrm{s}^{-1}]\), respectively, as well as a trap density \(n_i\,[\mathrm{m}^{-3}]\).

The temporal evolution of \(c_\mathrm{m}\) and \(c_{\mathrm{t}, i}\) are then given by:

(2)\[\frac {\partial c_\mathrm{m}} { \partial t} = \nabla \cdot (D\nabla c_\mathrm{m}) - \sum_i \frac{\partial c_{\mathrm{t},i}} { \partial t} + \sum_j S_j\]

(3)\[\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 \(S_j=S_j(x,y,z,t)\,[\mathrm{m}^{-3}\,\mathrm{s}^{-1}]\) is a source \(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.

Drift terms

Hydrogen can be carried by something other than its own concentration gradient. FESTIM writes all of these as a drift velocity added to the flux:

(4)\[J = -D \nabla c_\mathrm{m} + c_\mathrm{m} \mathbf{v}\]

The transport equation (2) becomes \(\partial c_\mathrm{m} / \partial t = -\nabla \cdot J + \dots\), so a drift term is assembled in divergence form, \(\nabla \cdot (c_\mathrm{m} \mathbf{v})\). This leaves behind the natural boundary condition (i.e. it’s automatically satisfied by the weak formulation of the problem). Correspondingly, the flux boundary condition constrains the total flux, drift included, and a boundary with no condition on it is genuinely no-flux (\(\nabla c_\mathrm{m} = 0\)).

What sets \(\mathbf{v}\) depends on the physics.

Advection (festim.AdvectionTerm) takes the velocity directly, for hydrogen carried by a moving fluid.

Note

Before FESTIM 2.2, festim.AdvectionTerm was assembled as \(\mathbf{v} \cdot \nabla c_\mathrm{m}\), which is implicitly divergence-free (\(\nabla \cdot \mathbf{v} = 0\)). However, it puts no flux through a boundary that carries no BC, so nothing was carried out through an outlet. Models with \(\nabla \cdot \mathbf{v} \neq 0\), or with an untagged boundary the flow crosses, will give different results.

Note

FESTIM does not stabilise the advection-diffusion form. When the cell Péclet number \(\mathrm{Pe} = |\mathbf{v}| h / D\) becomes large the solution can oscillate. Refine the mesh where the drift dominates to mitigate this issue.

Soret effect (festim.SoretTerm) [2, 3], also called thermophoresis, temperature-assisted diffusion, or thermodiffusion, drives hydrogen along a temperature gradient:

(5)\[J = -D \nabla c_\mathrm{m} - D\frac{Q^* c_\mathrm{m}}{k_B T^2} \nabla T\]

where \(Q^*\,[\mathrm{eV}]\) is the Soret coefficient (also called heat of transport) and \(k_B\) is the Boltzmann constant. For a positive \(Q^*\) hydrogen is driven along the temperature gradient (high to low). In a steady state closed container, the cocncentration profile follows \(\exp(Q^*/k_B T)\).

Electromigration (festim.ElectromigrationTerm) drives a charged species along an electric potential gradient \(\nabla\varphi\), the drift term of the Nernst-Planck equation:

(6)\[J = -D \nabla c_\mathrm{m} - \frac{z D c_\mathrm{m}}{k_B T} \nabla \varphi\]

where \(z\) is the charge number and \(\varphi\,[\mathrm{V}]\) is the electric potential. At equilibrium in a closed system, \(c_\mathrm{m} \propto \exp(-z\varphi/k_B T)\) (the Boltzmann distribution).

Conservation of chemical potential at interfaces

Continuity of local partial pressure \(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:

(7)\[P = \left(\frac{c_\mathrm{m}}{K_S}\right)^2\]

where \(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:

(8)\[P = \frac{c_\mathrm{m}}{K_H}\]

where \(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:

(9)\[\begin{split}\frac{c_\mathrm{m}^-}{K_S^-}&=&\frac{c_\mathrm{m}^+}{K_S^+} \\ &\mathrm{or}& \\ \frac{c_\mathrm{m}^-}{K_H^-}&=&\frac{c_\mathrm{m}^+}{K_H^+}\end{split}\]

where exponents \(+\) and \(-\) denote both sides of the interface.

At the interface between a Sievert and a Henry material:

(10)\[\left(\frac{c_\mathrm{m}^-}{K_S^-}\right)^2 = \frac{c_\mathrm{m}^+}{K_H^+}\]

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 [4]. The variable \(\theta\,[\mathrm{Pa}]\) is introduced and:

(11)\[\begin{split}\theta = \begin{cases} \left(\dfrac{c_\mathrm{m}}{K_S}\right)^2 & \text{in Sievert materials} \\ \\ \dfrac{c_\mathrm{m}}{K_H} & \text{in Henry materials} \end{cases}\end{split}\]

The variable \(\theta\) is continuous at interfaces.

Equations (2) and (3) are then rewritten and solved for \(\theta\). Note, the boundary conditions are also rewritten. Once solved, the discontinuous \(c_\mathrm{m}\) field is obtained from \(\theta\) and the solubilities by solving Equation (11) for \(c_\mathrm{m}\).

Arrhenius law

Many processes involved in hydrogen transport (e.g., diffusion, trapping/detrapping, desorption, etc.) are thermally activated. The coefficients characterising these processes (e.g., diffusivity, trapping/detrapping rates, recombination coefficients, etc.) are usually assumed to be temperature dependent and follow the Arrhenius law. According to the latter, the rate \(k(T)\) of a thermally activated process can be expressed as:

(12)\[k(T) = k_0 \exp \left[-\frac{E_k}{k_B T} \right]\]

where \(k_0\) is the pre-exponential factor, \(E_k\,[\mathrm{eV}]\) is the process activation energy, and \(T\,[\mathrm{K}]\) is the temperature.

Radioactive decay

Radioactive decay can be included in FESTIM. Neglecting other terms (diffusion, trapping…), the temporal evolution of the concentration of a species \(c\) is governed by:

(13)\[\frac{\partial c}{ \partial t} = -\lambda c\]

where \(\lambda\,[\mathrm{s}^{-1}]\) is the decay constant. In FESTIM it is possible to include radioactive decay in the evolution of the mobile hydrogen concentration and/or the trapped hydrogen concentrations.

Heat transfer

To properly account for the temperature-dependent parameters, an accurate representation of the temperature field is often required. FESTIM can solve a heat transfer problem governed by the heat equation:

(14)\[\rho C_p \frac{\partial T}{\partial t} = \nabla \cdot (\lambda \nabla T) + \sum_i Q_i\]

where \(T\) is the temperature, \(C_p\,[\mathrm{J}\,\mathrm{kg}^{-1}\,\mathrm{K}^{-1}]\) is the specific heat capacity, \(\rho\,[\mathrm{kg}\,\mathrm{m}^{-3}]\) is the material’s density, \(\lambda\,[\mathrm{W}\,\mathrm{m}^{-1}\,\mathrm{K}^{-1}]\) is the thermal conductivity and \(Q_i\,[\mathrm{W}\,\mathrm{m}^{-3}]\) is a volumetric heat source \(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

To fully pose the hydrogen transport problem and optionally the heat transfer problem, boundary conditions are required. Boundary conditions are separated in three categories: 1) enforcing the value of the solution at a boundary (Dirichlet’s condition) 2) enforcing the value of gradient of the solution (Neumann’s condition) 3) enforcing the value of the gradient as a function of the solution itself (Robin’s condition).

Dirichlet BC

In FESTIM, users can fix the mobile hydrogen concentration \(c_\mathrm{m}\) and the temperature \(T\) at boundaries \(\delta \Omega\) (Dirichlet):

(15)\[c_\mathrm{m} = f(x,y,z,t)~\text{on}~\delta\Omega\]

(16)\[T = f(x,y,z,t)~\text{on}~\delta\Omega\]

where \(f\) is an arbitrary function of coordinates \(x,y,z\) and time \(t\).

FESTIM has built-in Dirichlet’s boundary conditions for Sievert’s condition, Henry’s condition (see Equations (17) and (18), respectively).

(17)\[c_\mathrm{m} = K_S \sqrt{P}~\text{on}~\delta\Omega\]

(18)\[c_\mathrm{m} = K_H P~\text{on}~\delta\Omega\]

Plasma implantation approximation

Dirichlet’s boundary conditions can also be used to approximate plasma implantation in near surface regions to be more computationally efficient [5]. Let us consider a volumetric source term of hydrogen \(\Gamma=\varphi_{\mathrm{imp}}f(x)\), where \(f(x)\) is a narrow Gaussian distribution. The concentration profile of mobile species can be approximated by a triangular shape [6] with maximum at \(x=R_p\) (see the figure below).

Concentration profie with recombination flux and volumetric source term at \(x=R_p\). Dashed lines correspond to the time evolution

The expression of maximum concentration value \(c_{\mathrm{m}}\) can be obtained by expressing the flux balance at equilibrium:

(19)\[\varphi_{\mathrm{imp}} = \varphi_{\mathrm{recomb}} + \varphi_{\mathrm{bulk}}\]

where \(\varphi_{\mathrm{recomb}}\) is the recombination flux and \(\varphi_{\mathrm{bulk}}\) is the migration flux into the bulk. \(\varphi_{\mathrm{bulk}}\) can be expressed as:

(20)\[\varphi_{\mathrm{bulk}} = D \cdot \frac{c_{\mathrm{m}}}{R_d(t)-R_p}\]

with \(R_d\) the diffusion depth and \(R_p\) the implantation range. When \(R_d \gg R_p\), \(\varphi_{\mathrm{bulk}} \rightarrow 0\). Equation (19) can therefore be written as:

(21)\[\begin{split}\varphi_{\mathrm{recomb}} &=& D \cdot \frac{c_{\mathrm{m}} - c_0}{R_p} = \varphi_{\mathrm{imp}}\\ \Leftrightarrow c_{\mathrm{m}} &=& \frac{\varphi_{\mathrm{imp}} R_p}{D} + c_0\end{split}\]

Assuming second order recombination, \(\varphi_{\mathrm{recomb}}\) can also be expressed as a function of the recombination coefficient \(K_r\):

(22)\[\begin{split}\varphi_{\mathrm{recomb}} &=& K_r c_0^2 = \varphi_{\mathrm{imp}}\\ \Leftrightarrow c_0 &=& \sqrt{\frac{\varphi_{\mathrm{imp}}}{K_r}}\end{split}\]

By substituting Equation (22) into (21) one can obtain:

(23)\[c_\mathrm{m} = \frac{\varphi_{\mathrm{imp}} R_p}{D} + \sqrt{\frac{\varphi_{\mathrm{imp}}}{K_r}}\]

Similarly, dissociation can be accounted for:

(24)\[c_\mathrm{m} = \frac{\varphi_{\mathrm{imp}} R_p}{D} + \sqrt{\frac{\varphi_{\mathrm{imp}}+K_d P}{K_r}}\]

where \(K_d\) is the dissociation coefficient.

When recombination is fast (i.e. \(K_r\rightarrow\infty\)), Equation (24) can be reduced to:

(25)\[c_\mathrm{m} = \frac{\varphi_{\mathrm{imp}} R_p}{D}\]

Since the main driver of for the diffusion is the value \(c_{\mathrm{m}}\), when \(R_p\) is negligible compared to the dimension of the simulation domain, one can simply impose Equations (24) and (25) at boundaries \(\delta \Omega\).

Neumann BC

One can also impose hydrogen fluxes or heat fluxes at boundaries (Neumann). The condition is on the total flux \(J\) of Equation (4); writing it out for the case of no drift term, where \(J = -D\nabla c_\mathrm{m}\):

(26)\[J \cdot \mathrm{\textbf{n}} = -D\nabla c_\mathrm{m} \cdot \mathrm{\textbf{n}} =f(x,y,z,t)~\text{on}~\delta\Omega\]

With a drift term the same condition constrains \((-D\nabla c_\mathrm{m} + c_\mathrm{m}\mathbf{v}) \cdot \mathrm{\textbf{n}}\) instead, and festim.SurfaceFlux reports that same total. A boundary with no condition on it is therefore a wall. Where the drift should instead carry the species out of the domain, festim.OutflowBC cancels the drift part, leaving \(-D\nabla c_\mathrm{m} \cdot \mathrm{\textbf{n}} = 0\).

(27)\[-\lambda\nabla T \cdot \mathrm{\textbf{n}} = f(x,y,z,t)~\text{on}~\delta\Omega\]

where \(\mathrm{\textbf{n}}\) is the normal vector of the boundary.

Robin BC

Recombination and dissociation fluxes can also be applied:

(28)\[J \cdot \mathrm{\textbf{n}} = -D\nabla c_\mathrm{m} \cdot \mathrm{\textbf{n}} = K_d P - K_r c_\mathrm{m}^{\{1,2\}} ~\text{on}~\delta\Omega\]

In Equation (28), the exponent of \(c_\mathrm{m}\) is either 1 or 2 depending on the reaction order. These boundary conditions are Robin boundary conditions since the gradient is imposed as a function of the solution.

Finally, convective heat fluxes can be applied to boundaries:

(29)\[-\lambda\nabla T \cdot \mathrm{\textbf{n}} = h (T-T_{\mathrm{ext}})~\text{on}~\delta\Omega\]

where \(h\) is the heat transfer coefficient and \(T_{\mathrm{ext}}\) is the external temperature.

References

[1]

A. McNabb and P. K. Foster. A New Analysis of the Diffusion of Hydrogen in Iron and Ferritic Steels. Trans. Metall. Soc. AIME, 227:618, 1963.

[2]

J. H. Pendergrass. Temperature-dependent ordinary and thermal diffusion of hydrogen isotopes through thermonuclear reactor components. Technical Report, Los Alamos National Laboratory (LANL), 1976. URL: http://www.osti.gov/servlets/purl/7333557/.

[3]

G. R. Longhurst. The soret effect and its implications for fusion reactors. J. Nucl. Mater., 131(1):61–69, 1985. URL: https://linkinghub.elsevier.com/retrieve/pii/0022311585904258.

[4]

R. Delaporte-Mathurin, E. A. Hodille, J. Mougenot, Y. Charles, G. De Temmerman, F. Leblond, and C. Grisolia. Influence of interface conditions on hydrogen transport studies. Nucl. Fusion, 61(3):036038, 2021. URL: https://iopscience.iop.org/article/10.1088/1741-4326/abd95f.

[5]

R. Delaporte-Mathurin. Hydrogen transport in tokamaks : Estimation of the ITER divertor tritium inventory and influence of helium exposure. Theses, Université Paris-Nord - Paris XIII, 2022. URL: https://theses.hal.science/tel-04004369.

[6]

K. Schmid. Diffusion-trapping modelling of hydrogen recycling in tungsten under ELM-like heat loads. Phys. Scr., T167:014025, 2016. URL: https://iopscience.iop.org/article/10.1088/0031-8949/T167/1/014025.