UKAEA FESTIM Workshop 2026

 

Kaelyn Dunnell & Tez Orr

 June 29, 2026

 

The FESTIM team at MIT and beyond

James Dark
Postdoctoral associate
Huihua Yang
Postdoctoral associate
Rémi Delaporte-Mathurin
Research Scientist
Kaelyn Dunnell
Graduate Student
Chirag Khurana
Graduate Student
MIT Team
Joe Dean
Research Associate
University of Cambridge
Jørgen Dokken
Senior Research Scientist
Simula Research Laboratory
FEniCS collaborators
FESTIM Community
31+ contributors

DFT

Multi-scale hydrogen transport

Y. Ferro et al 2023 Nucl. Fusion 63 036017

Length scale

Time scale

MD

Length scale

Time scale

DFT

potentials

Grigorev et al 2018 JNM, 508 451-458

Multi-scale hydrogen transport

Component scale modelling

Length scale

Time scale

MD

DFT

D, S, other coeffs.

Multi-scale hydrogen transport

Length scale

Time scale

MD

DFT

Component scale modelling

Fuel cycle modelling

Residency times, fluxes, ...

Multi-scale hydrogen transport

Length scale

Time scale

MD

DFT

Component scale modelling

Fuel cycle modelling

Abstraction

Multi-scale hydrogen transport

Tritium Transport in Materials

 

Diffusion

\varphi = - D \ \nabla c

\( \varphi\): diffusion flux

\( D \): diffusion coefficient

\( c \): mobile hydrogen concentration

\( S\): source term

Fick's 1st law of diffusion

\frac{\partial c}{\partial t} = - \nabla \cdot \varphi + S

Fick's 2nd law of diffusion

Trapping

H

Potential energy

Distance

Diffusion barrier

Trap binding energy

Trapping energy

Common assumption:

\( E_k = E_D \)

Trapping

\mathrm{H} + [\ \ \ ]_\mathrm{trap} \ \substack{p \\[-1em] \longleftarrow\\[-1em] \longrightarrow \\[-1em] k} \ [\mathrm{H}]_\mathrm{trap}

0D

\frac{\partial c_\mathrm{t}}{\partial t} = -\frac{\partial c_\mathrm{m}}{\partial t} = k \ c_\mathrm{m} \ (n_\mathrm{trap,0} - c_\mathrm{t}) - p c_\mathrm{t}
c_\mathrm{m}
n_\mathrm{free \ trap}
c_\mathrm{t}

Total concentration of traps

McNabb & Foster Model

\frac{\partial c_\mathrm{m}}{\partial t} = \nabla \cdot (D \nabla c_\mathrm{m}) + S - \sum \frac{\partial c_\mathrm{t,i}}{\partial t}
\frac{\partial c_\mathrm{t, i}}{\partial t} = k_i \ c_\mathrm{m} \ (n_\mathrm{trap,i} - c_\mathrm{t}) - p_i c_\mathrm{t}

Particle implantation

Ziegler et al. 2010. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 268 (11): 1818–23. https://doi.org/10.1016/j.nimb.2010.02.091.

Implantation range

Implantation range & width and reflection coefficient can be computed with SRIM, SDTRIM...

Mutzke et al, SDTrimSP Version 6.00 2019

S = (1-r) \ \Gamma_\mathrm{incident} \ f(x)

\(\Gamma_\mathrm{incident} \): incident flux (particle/m2/s)

\( f(x) \): Gaussian distribution (/m)

\(r \): reflection coefficient

or

Surface effects

H2 molecules

Metal lattice

Surface effects

\varphi_\mathrm{dissociation} = K_\mathrm{d} \ P

Dissociation coefficient (H/m2/s/Pa)

Partial pressure of H (Pa)

Adsorbed H

Metal lattice

Surface effects

\varphi_\mathrm{dissociation} = K_\mathrm{d} \ P
\varphi_\mathrm{recombination} = K_\mathrm{r} \ c^2

Metal lattice

Recombination coefficient (m4/s)

Concentration (H/m3)

Surface effects

\varphi_\mathrm{dissociation} = K_\mathrm{d} \ P
\varphi_\mathrm{recombination} = K_\mathrm{r} \ c^2

Metal lattice

\varphi_\mathrm{net} = \varphi_\mathrm{dissociation} - \varphi_\mathrm{recombination}

Waelbroeck model

Surface effects

\varphi_\mathrm{dissociation} = K_\mathrm{d} \ P
\varphi_\mathrm{recombination} = K_\mathrm{r} \ c^2

Metal lattice

At equilibrium:

\varphi_\mathrm{net} = 0
c = \sqrt{\frac{K_\mathrm{d}}{K_\mathrm{r}}} \ \sqrt{P}
\varphi_\mathrm{net} = \varphi_\mathrm{dissociation} - \varphi_\mathrm{recombination}
c = K_S \ \sqrt{P}

Sievert's law of solubility

Surface effects

\varphi_\mathrm{dissociation} = K_\mathrm{d} \ P
\varphi_\mathrm{recombination} = K_\mathrm{r} \ c^1

Non-metallic liquid

At equilibrium:

\varphi_\mathrm{net} = 0
c = \frac{K_\mathrm{d}}{K_\mathrm{r}} \ P
\varphi_\mathrm{net} = \varphi_\mathrm{dissociation} - \varphi_\mathrm{recombination}
c = K_H \ P

Henry's law of solubility

Interfaces

Material 1

Material 2

Case 1:

Metal-Metal

- D_1 \nabla c_1 = - D_2 \nabla c_2
P = \left( \frac{c}{K_S} \right) ^2
\left( \frac{c_1}{K_{S,1}} \right)^2 =\left( \frac{c_2}{K_{S,2}} \right)^2

Sievert's law

Flux continuity

Conservation of chemical potential

Interfaces

Material 1

Material 2

Case 2:

Non metal-non metal

- D_1 \nabla c_1 = - D_2 \nabla c_2
P = \frac{c}{K_H}
\frac{c_1}{K_{H,1}} = \frac{c_2}{K_{H,2}}

Henry's law

Flux continuity

Conservation of chemical potential

Interfaces

Material 1

Material 2

Case 3:

Metal-Non metal

- D_1 \nabla c_1 = - D_2 \nabla c_2
P = \left( \frac{c}{K_S} \right) ^2
\left( \frac{c_1}{K_{S,1}} \right)^2 = \frac{c_2}{K_{H,2}}

Sievert's law

P = \frac{c}{K_H}

Henry's law

\frac{\partial c_\mathrm{m}}{\partial t} = \nabla \cdot \varphi_\mathrm{diffusion} + S - \sum \frac{\partial c_\mathrm{t,i}}{\partial t}
\frac{\partial c_\mathrm{t, i}}{\partial t} = k_i \ c_\mathrm{m} \ (n_\mathrm{trap,i} - c_\mathrm{t}) - p_i c_\mathrm{t}

McNabb & Foster model

Challenges

  • Number of degrees of freedom
  • Interface discontinuities

Governing equations

- D_1 \nabla c_1 = - D_2 \nabla c_2
\left( \frac{c_1}{K_{1}} \right)^2 =\left( \frac{c_2}{K_{2}} \right)^{\{1,2\} }

We can solve these numerically

Component modelling

Experimental analysis

3 main numerical methods

Finite Difference Method (FDM)

Finite Element Method (FEM)

Finite Volume Method (FVM)

Let's not bother

\frac{\partial c }{\partial x} \approx \frac{c(x=x_{i+1}) - c(x=x_{i})}{x_{i+1} - x_{i}}

FESTIM today

See more FESTIM stats:
festim-dev.github.io/pose

See all FESTIM papers:

festim.readthedocs.io/en/latest/publications.html

FESTIM is used worldwide

Docs users

Fusion machine developers

Research organisations

What changed in FESTIM2?

Upgrade from legacy-FEniCS to FEniCSx

  Before

  • Change of variable required - did not scale well
  • Single solver method available
  • No submesh coupling
  • Limited mesh element types

 

 

 Now

  • Mixed-domain: no substitution needed
  • Multiple solver methods (Penalty, Nitches)
  • Submesh coupling: 1D to 2D problems
  • Mixed topology meshes
  • Improved parallel I/O 

 

~11x faster

  Before

  • Monolithic approach
  • Change of variable required
    • Poor scaling
  • Little control over interface equations
  • Adding advection was tricky

  Now

  • Mixed-domain: no substitution needed
  • No change of variable needed
  • Improved performance
  • Full control over interface conditions
  • Very suitable for multi-physics 

 

Multi-material handling

Multi-species and Reactions

  Before

  • Single mobile species only
  • No multi-isotope support
  • No multi-occupancy trapping
  • Rigid trapping reactions

  Now

  • Arbitrary species
    • Multi-isotope
    • Multi-occupancy trapping
    • Defect diffusion
  • Flexible reactions
  • Bulk and surface reactions

 

The mixed formulation is great for multiphysics

velocity, temperature, turbulent viscosity...

tritium production

Heat source

Temperature

Neutronics

Thermo-hydraulics

Tritium transport

Examples

Using FESTIM for estimating tritium contamination in a heat exchanger

ARC Breeding Blanket

ARC blanket geometry

T source

Heat source

Velocity

Temperature

Turbulent viscosity

Tritium concentration

TDS analysis: codeposits

Patino et al. Nuclear Materials and Energy 41 (2024) 101808

  • Simulation of W codeposited layers
     
  • Influence of partial pressure
     
  • 10 different traps!

Thank you for your attention!

Any questions?