FESTIM: overview of hydrogen transport simulation capabilities

Rémi Delaporte-Mathurin^1,2, James Dark², Ian Wojtowicz³, Jørgen S. Dokken⁴

and the FESTIM community

¹NeutroGen Technologies ²MIT Plasma Science and Fusion Center ³University of Minnesota ⁴Simula Research Laboratory

4th RCM, CRP Hydrogen Permeation in Fusion-relevant Materials | IAEA, Vienna | 28 October 2026

Diffuse, trap, react, release

Hydrogen transport from plasma-facing wall through breeder to coolant and gas enclosure

One framework for solids, liquids, gases and the interfaces between them

Zoom in: the wall at atomic scale

Atomic-scale processes: dissociation, diffusion, trapping, recombination, implantation, isotope exchange, hydrides

\partial_t c_\mathrm{m} = \nabla \cdot \left( D \nabla c_\mathrm{m} - c_\mathrm{m}\, \mathbf{u} \right) + S - \sum_i \partial_t c_{\mathrm{t},i}

\partial_t c_{\mathrm{t},i} = k_i\, c_\mathrm{m} \left( n_i - c_{\mathrm{t},i} \right) - p_i\, c_{\mathrm{t},i}

mobile hydrogen: diffusion, drift, sources

one equation per species

FESTIM hello-world

1D, 2D, 3D, up to CAD-meshed components
Multi-material, multi-isotope, arbitrary species and reactions
Traps, surface kinetics, Sieverts / Henry, coupled heat transfer
Python API: physics vocabulary, no finite elements exposed
Built on FEniCSx, MPI-parallel
Open source (Apache 2.0), on conda-forge

import festim as F

model = F.HydrogenTransportProblem()
model.mesh = F.Mesh1D(np.linspace(0, 1e-3, 1001))

steel = F.Material(D_0=1e-7, E_D=0.2)
bulk = F.VolumeSubdomain1D(1, borders=[0, 1e-3], material=steel)
left = F.SurfaceSubdomain1D(1, x=0)
right = F.SurfaceSubdomain1D(2, x=1e-3)
model.subdomains = [bulk, left, right]

T = F.Species("T")
model.species = [T]
model.temperature = 600
model.boundary_conditions = [
    F.FixedConcentrationBC(left, value=1e20, species=T),
    F.FixedConcentrationBC(right, value=0, species=T),
]
model.exports = [F.SurfaceFlux(field=T, surface=right)]

model.initialise()
model.run()

FESTIM across this meeting

Atomistic modelling

Each trap energy gives its own TDS peak

Ion- and plasma-driven

Retention in an ITER divertor monoblock

Gas-driven permeation

Tritium permeator probe in a flowing liquid

Retention

Automated fitting of trap parameters to TDS

Retention: microstructure

Hydrogen in a polycrystalline microstructure

TDS of self-damaged tungsten, FESTIM versus experiment

Neutron-irradiated materials

Self-damaged W, 0 to 2.5 dpa: FESTIM vs TDS

Tritium inventory in ITER

HISP, an open-source workflow built on FESTIM, with the ITER Organization: how much tritium stays in the wall, and how well can we get it out?

ITER first wall and divertor split into bins

■ tungsten ■ boron co-deposit
HISP wall bins, from the open data

Plasma code outputs
fluxes, heat loads

Wall split in bins
first wall + divertor

1D FESTIM model
per bin, W or B

2-week DT scenarios
baking, GDC, DD

→

K. Dunnell, A. Lleal, E. A. Hodille, J. Dufour, R. Delaporte-Mathurin, T. Wauters, Hydrogen Inventory Simulations for PFCs (HISP), Nucl. Mater. Energy (2026)

Multiphysics, multi-code coupling

OpenMC: tritium source

OpenFOAM: velocity

FESTIM: tritium

Neutronics and CFD fields passed straight into FESTIM, no file conversion

ARC breeding blanket: ~243 mg tritium inventory | Dark, Sircar, Bae, Dokken & Delaporte-Mathurin, arXiv:2608.05398

And up to the fuel cycle: FESTIM components plug into PathSim system models

Gas enclosures and pumps

NEW

Permeation rigs, purged gaps, gloveboxes: a wall feeding a gas volume whose pressure feeds back on the wall.

\frac{\mathrm{d}P}{\mathrm{d}t} = \frac{R\,T}{V} \int_{\partial \Omega} \varphi \, \mathrm{d}s \; - \; \text{pumping}

Pressure is one extra unknown, solved together with the wall
Surfaces exchange with the gas via Sieverts or surface reactions
Pumps and purge streams are openings on the enclosure

Two membranes with a pumped gap. Pump on: gap pressure 25 Pa → 0, hydrogen reaching the second membrane down ~5 orders of magnitude.

Codimensional coupling

NEW

Subdomains can now be one dimension lower than the mesh: grain boundaries, thin films, barrier coatings.

Here: hydrogen reaches depth through the boundary network long before the bulk front arrives.

Microstructure effects on retention and permeation, without meshing every boundary in 3D

Much of this built by FESTIM Fellow Ian Wojtowicz (U. Minnesota)

Can be used to simulate trapping at interfaces!

Grain boundary tutorials

A fast path through a slow bulk

Concentration in the grain without and with a 1D grain boundary manifold

One grain boundary as a 1D manifold, no cells of its own: 4.5× the flux through the top of the square

Diffusion in a polycrystal

Hydrogen concentration in a 79-grain Voronoi polycrystal

79 Voronoi grains, D_GB = 1000 D_grain: hydrogen runs ahead along the boundaries

Both run in the browser, nothing to install: festim-workshop.readthedocs.io

Drift terms

NEW

Hydrogen does not only follow its own gradient: it drifts with temperature, electric potential, stress, and flow.

J = -D \nabla c + c\, \mathbf{u}, \qquad \mathbf{u} = \sum_j \mathbf{u}_j

Soret (thermodiffusion)

\mathbf{u} = -\frac{D\, Q^{*}}{k_B T^{2}} \nabla T

Electromigration

\mathbf{u} = -\frac{z\, D}{k_B T} \nabla \varphi

Stress-assisted diffusion

\mathbf{u} = \frac{D\, \bar{V}_\mathrm{H}}{k_B T} \nabla \sigma_h

Advection

\mathbf{u} = \mathbf{u}_\mathrm{fluid} \quad (\text{e.g. from OpenFOAM})

Long requested by users. Each drift is one more term in the weak form, so it works with multi-material, manifolds, enclosures and coupled heat transfer.

Multi-occupancy trapping

One defect, several hydrogen atoms, each level with its own binding energy.

\mathrm{H} + [n\mathrm{H}] \;\rightleftharpoons\; [(n{+}1)\mathrm{H}]

Each level is a species, linked by reactions. Works for mixed isotopes in the same trap too.

Binding energy per occupancy level: straight from atomistic modelling

Tutorial: festim-workshop.readthedocs.io, Reactions

Can you trust it?

Verification

650+ tests on every push and pull request, on two DOLFINx versions
~60 system tests against manufactured and exact solutions, 1D to 3D
The TMAP benchmark suite
Cross-checked against TMAP8, and independently against MHIMS and mHIT

Validation

8 experimental campaigns: permeation, retention in damaged W, Be co-deposition, H in Ti, EUROFER
BABY-1L breeding campaign (MIT, ClLiF): calibrated on 4 runs, predicts 2 withheld runs without recalibration

H. Yang et al., in preparation

Every case is a runnable script. The report rebuilds from them, in public.

festim-vv-report.readthedocs.io

LIBRA: breeding in BABY-1L

BACKUP

BABY-1L configuration: ClLiF salt in an Inconel 625 vessel under D-T neutrons, with inner-vessel and outer-vessel tritium collection

Predicted versus measured cumulative tritium recovery for withheld runs 4 and 6

ClLiF salt, 14 MeV D-T neutrons

TBR = 2.4×10^-3 per source neutron

2D FESTIM: calibrated on 4 runs, predicts 2 withheld runs

Salt / Inconel interface far from LTE: Da < 0.02

87 to 95% of outer-vessel tritium (runs 3, 4) was bred in earlier runs

R. Delaporte-Mathurin et al., Nucl. Fusion 66 (2026) 076046 | H. Yang et al., preprint (2026) | data: github.com/LIBRA-project

A community project

Map of FESTIM documentation users worldwide, from Google Analytics

30+ contributors

30+ institutions

8k forum views / month

Apache 2.0

NumFOCUS affiliated

Interactive tutorials, nothing to install: festim-workshop.readthedocs.io

conda install -c conda-forge festim | Monthly dev meetings, open to all | Discourse and Slack

Is local equilibrium enough?

FESTIM, TMAP8, TESSIM-X: interfaces modelled by local thermodynamic equilibrium (LTE)

\frac{c_A}{K_{S,A}} = \frac{c_B}{K_{S,B}}

c_B = K_{H,B} \left( \frac{c_A}{K_{S,A}} \right)^2

Sieverts / Sieverts (metal / metal)

Sieverts / Henry (metal / salt)

Three hidden assumptions:

A1 fast interfacial equilibration

A2 a single exchange pathway

A3 carrier species known on each side

Motivated by W/Cu (DFT + rate equations): under a net flux, the interface steady state is not equilibrium

Y. Silva-Solís, J. Denis, E. A. Hodille, Y. Ferro, Phys. Rev. Mater. 10 (2026) 045404

A2 and A3 fail at metal / molten salt interfaces

and they matter more than A1

From a law to reaction channels

Four interface models: LTE, first-order exchange, recombination, competing channels

w_r = k_r^+ \prod_i c_i^{\nu_{ir}} - k_r^- \prod_\alpha c_\alpha^{\mu_{\alpha r}}

one reversible, mass-action rate per channel

\frac{k_r^+}{k_r^-} = K_r(T)

detailed balance: same data as LTE

LTE = fast-kinetics limit of a single channel

Only new input: the magnitude of k⁺

Metal / metal: fast means LTE

Interfacial ratio and departure from LTE versus Damkohler number

\mathrm{Da} = \frac{k^+ L}{D}

bulk resistance over interface resistance

Departure from LTE scales as 1/Da
1.5% at Da = 10², 1.5×10^-6 at Da = 10⁶
Matches the analytical solution to 10^-12

Metal / metal with Da ≫ 1: LTE is a safe approximation

Metal / FLiBe: two exit channels

R: 2 H ⇄ H₂ quadratic

F: H + F ⇄ HF linear

fluorine activity set by the salt redox state

\mathcal{B} = \frac{w_\mathrm{F}}{2\, w_\mathrm{rec}}

n = \frac{2 + \mathcal{B}}{1 + \mathcal{B}}

branching ratio, kinetic

apparent exponent, 1 to 2

Regime map of LTE failure indicator over Damkohler number and branching ratio

Near B ≈ 1 no algebraic interface law exists, however fast the kinetics

Metal / salt: from √P to P

Metal/salt interface: hydrogen leaves the metal as H2 (recombination) or HF (fluorination)

Pressure exponent of the flux versus interfacial H2 pressure for three fluorine activities, FESTIM versus analytical

Same metal, same salt: the exponent sits anywhere between 0.5 and 1

A measured exponent reflects the salt redox state, not an intrinsic salt property

R. Delaporte-Mathurin & J. Dark, When chemical potential continuity fails: kinetic interface models for hydrogen isotope transport, arXiv:2608.26231