Multiphysics tritium transport modelling of the ARC breeding blanket with FESTIM

OpenMC utilised the DAGMC library [47] for geometry modelling, in which volumes are represented by tessellated surfaces. The geometry can be created by meshing an existing CAD file and labelling the volumes bound by the surfaces. The labels are then read in OpenMC, where each volume is treated as an OpenMC cell with its material assignment. An initial uniform temperature of $800\text{\,}\mathrm{K}$ was assumed for the FLiBe. The DAGMC geometry is prepared using FERMI, which leverages the Python API of Cubit [6]. Additionally, in OpenMC, volumetric tetrahedral meshes can be provided for tallies. In the same workflow, a volumetric mesh is generated from the same CAD file and provided to OpenMC for spatially resolved tritium generation and heat deposition tallies.

The neutron source is modelled as a point cloud, with the space, energy, and angular distributions sampled using the tools from Fausser et al. [24]. The ENDF/B-VII.1 [11] nuclear data library was used. To ensure good statistics for the tetrahedral mesh tallies, a subroutine was set using the OpenMC Python API to monitor the maximum Monte Carlo relative error. The subroutine was set to continue running until the highest relative error across the tetrahedral mesh tallies fell below $1\text{\,}\%$. Nine hundred (900) million particles were simulated, with 300,000 particles per batch, using the OpenMC built-in Method of Automatic Generation of Importances by Calculation (MAGIC) [19] to update the weight windows of each batch for variance reduction. The material assignments for each layer, not considering any impurities, are shown in Table 1.

The thermal hydraulic behaviour of the ARC blanket coolant was modelled using the open-source CFD framework OpenFOAM [37], which provides a robust finite volume formulation for coupled momentum and heat transfer problems in complex geometries. OpenFOAM was selected for its established use in nuclear and fusion thermal hydraulics [20, 36, 10], its flexibility for incorporating volumetric source terms, and its suitability for integration within a broader multiphysics workflow [40].

The simulations were performed using a custom transient solver, rhoBlanketFoam, derived from the standard PIMPLE algorithm, and using a spatially resolved heat deposition field mapped from OpenMC neutronics. This custom solver was developed as part of the FERMI project. The solver advances the conservation equations for mass, momentum, and energy for a single-phase fluid. Although the equations are formulated in a mildly compressible form, the density was treated as constant, resulting in an effectively incompressible flow. The solved fields include velocity, pressure, temperature, and density. Turbulence was modelled using a Reynolds-averaged Navier–Stokes approach with the $k$-$\omega$ SST closure, including standard wall functions for near-wall treatment. Gravity was set to zero, and no buoyancy effects were included. Hybrid first- and second-order spatial discretisation schemes are used for flux calculations in the momentum and energy equations, and pure second-order schemes are used for derivatives, with non-orthogonal correctors to account for mesh non-orthogonality, given the arbitrary blanket profile.

Thermal transport was formulated in terms of sensible enthalpy, assuming constant thermophysical properties; see Table 2. The volumetric heat generation term, representing internal heating within the fluid domain, was obtained from OpenMC neutronics simulations and mapped using tools in the FERMI workflow. Turbulent heat transfer was accounted for through a turbulent thermal diffusivity with a constant turbulent Prandtl number.

The computational domain represents a reduced $5\text{\,}\mathrm{\SIUnitSymbolDegree}$ sector of the blanket coolant channel (see Figure 2), exploiting geometric symmetry where applicable. A symmetry condition was imposed on the sector’s side walls, consistent with the assumption of perfect geometric and flow symmetry among the 72 identical sectors that make up the full blanket. Velocity boundary conditions consisted of a fixed inlet velocity, no-slip walls, a zero-gradient outlet, and symmetry conditions at the toroidal faces; see Table 3. Pressure was prescribed at the outlet, with zero-gradient conditions elsewhere. The temperature boundary conditions included a fixed inlet temperature, prescribed temperatures at selected solid–fluid interfaces, adiabatic conditions on the remaining walls, and zero-gradient conditions at the outlet. Turbulent quantities were specified at the inlet to correspond to $20\text{\,}\%$ turbulence intensity and treated with wall functions at solid boundaries.

This modelling approach captures the dominant forced convection heat transfer behaviour in the blanket coolant but neglects density variations, buoyancy effects, and multiphase phenomena. The use of Reynolds-averaged turbulence modelling limits the resolution of transient and small-scale flow structures. In addition, volumetric source terms were prescribed based on the heat deposition from neutronics results. These assumptions are appropriate for the present study, in which fluid mechanics primarily provides temperature and flow fields for subsequent multiphysics analyses rather than being the primary focus. The mapping of neutronics heat deposition to the OpenFOAM mesh is handled within FERMI through .vtk files. A Python converter tool was developed, building on the
vtkUnstructuredToFoam capability and using mapFields within OpenFOAM to interpolate between OpenMC and OpenFOAM meshes [4]. The nearest-neighbour mapping method was found to conserve global energy between meshes to within $0.1\text{\,}\%$.

Tritium transport within the ARC breeding blanket is governed by the combined effects of diffusion, bulk advection due to molten salt flow, and volumetric generation from nuclear reactions. In the liquid immersion blanket concept, tritium produced in the FLiBe breeder is transported by the coolant before being extracted or interacting with surrounding structures. Capturing this behaviour, therefore, requires a transport model that accounts for spatially varying source terms, temperature-dependent material properties, and mixing induced by the coolant flow, while remaining compatible with the multiphysics workflow described above.

Under these assumptions, tritium transport in the breeder is modelled using a macroscopic advection–diffusion formulation, which is well suited to continuum-scale simulations of molten salts.

Equation (1) describes the rate equation used to model tritium transport, with a coupling advection term, $\nabla\cdot(\mathbf{u}c_{\mathrm{m}})$, and a diffusion term, $\nabla\cdot\left(D\nabla c_{\mathrm{m}}\right)$, to describe the diffusion of mobile hydrogen particle concentrations, $c_{\mathrm{m}}$, in units of m^-3. $D=D_{0}\cdot\exp(-E_{\text{D}}/(k_{B}\,T))$ is the Fickian diffusion coefficient of interstitial tritium in units of $\mathrm{m}^{2}\,\mathrm{s}^{-1}$, $k_{B}$ = $8.617\text{\times}{10}^{-5}\text{\,}\mathrm{e}\mathrm{V}\,\mathrm{K}^{-1}$ is the Boltzmann constant, and $S$ is the volumetric source term of mobile tritium in units of $\mathrm{m}^{-3}\,\mathrm{s}^{-1}$. The mean molecular diffusivity of tritium in FLiBe was evaluated from a hydrogen transport properties database (HTM) [21]; see Table 4. The Soret effect is not considered in this work. Hydrogen trapping is not considered in this work, as only the liquid breeder is present. The source term is obtained from the tritium production tally generated by OpenMC and ported to DOLFINx for use in the formulation via the openmc2dolfinx package. $\mathbf{u}$ represents the velocity field of the liquid FLiBe breeder material in units of $\mathrm{m}\,\mathrm{s}^{-1}$, and similarly to the source term, is ported to DOLFINx from OpenFOAM using the package foam2dolfinx.

Turbulence in the breeder flow enhances tritium mixing and transport well beyond what molecular diffusion alone can achieve. To account for this effect, the tritium transport model includes a turbulent diffusion term by defining an effective diffusivity:

where $D_{\text{turb}}$ is the turbulent diffusivity of hydrogen, $\nu_{t}$ is the turbulent kinematic viscosity (obtained from the fluid dynamics simulation), and $\mathrm{Sc}_{t}$ is the turbulent Schmidt number. The turbulent Schmidt number, $\mathrm{Sc}_{t}$, characterises the ratio between turbulent momentum diffusivity and mass diffusivity. For molten salt coolants such as FLiBe, experimentally validated values of the turbulent Schmidt number are scarce; in the absence of material-specific data, the adopted value falls within the commonly assumed range for turbulent mass transport in liquids and reflects enhanced scalar mixing relative to momentum diffusion. This formulation enables the hydrogen transport model to respond dynamically to local turbulence intensity, capturing enhanced tritium transport in regions of high shear or recirculating flow.

To ensure numerical stability when solving the advection-diffusion equation at high Péclet numbers, a stabilisation method was also employed. Solving advection-dominated transport problems with the finite element method (FEM) can result in oscillations unless additional stabilisation terms are introduced. A consistent stabilisation strategy was adopted, analogous to streamline upwind Petrov–Galerkin (SUPG) or artificial diffusion methods [13], which introduce additional diffusion aligned with the direction of flow to preserve physical accuracy without overly damping the solution. This improves the solver’s convergence, particularly in regions with strong gradients or recirculating flows where advection dominates. This is implemented by

where $D_{\mathrm{art}}$ is the artificial diffusion coefficient, $\delta$ is a dimensionless tuning parameter, $h$ is the mesh element size in $\mathrm{m}$, and $\|\mathbf{u}\|$ is the magnitude of the velocity. This value is used in tandem with the turbulent diffusivity term to evaluate the effective diffusivity term, $D_{\mathrm{eff}}$ given by equation (4), which is then used in place of the standard hydrogen diffusivity, $D$; see equation (1).

A homogeneous Dirichlet boundary condition is applied to the inlet, which assumes a perfect purification system.

The FESTIM code [15] was employed and extended to model tritium transport in the ARC breeding blanket, solving equation (1). It uses the FEM to simulate diffusion and trapping of hydrogen isotopes based on macroscopic rate equations [31].

The OpenMC results (figure 3) show a dramatic reduction with increasing distance from the neutron source (i.e., the plasma). The tritium generation rate and heat deposition values are converted to SI units and normalised by the volume of each mesh element to obtain volumetric source terms.

In the OpenFOAM simulations, the flow velocities in the cooling channel are much higher than in the blanket due to the channel’s much smaller cross-sectional area; flow velocities in the channel reach a maximum of $12\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$ compared with $2\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$ in the blanket. Although the channel lies closer to the plasma and therefore experiences much more intense nuclear heating (see figures 3(b) and 4(a)), this high velocity, together with the fixed inlet temperature of the incoming FLiBe, provides efficient heat removal, so the channel stays comparatively cool, with temperatures ranging from approximately $800\text{\,}\mathrm{K}$ at the inlet to $870\text{\,}\mathrm{K}$. In the blanket, by contrast, the much slower flow removes heat less effectively, and temperatures peak at around $1000\text{\,}\mathrm{K}$.

The flow field shows a region of very high mixing when the FLiBe flows through the channel and interconnects with the blanket. The behaviour is similar to that of a jet, since the fast-moving molten salt in the channel enters a region of slowly moving fluid in the blanket. This is potentially a region of high pressure drop and warrants further investigation in the ARC design. As the flow moves upwards, it carries heat and tritium produced by neutron interactions with the molten salt. Towards the top, given that the top divertor creates a flow restriction near the outlet, a recirculation zone is observed (see figures 4(b) and 6(d)). These flow features have important implications for tritium transport within the blanket, as discussed next.

Under steady-state conditions, the total tritium inventory within the liquid breeder is predicted to be approximately $243\text{\,}\mathrm{m}\mathrm{g}$. Peak local tritium concentrations reach values on the order of $2.5\text{\,}\mathrm{m}\mathrm{g}\,\mathrm{m}^{-3}$. At the blanket outlet, a steady-state tritium mass flux of approximately $1.0\text{\,}\mathrm{m}\mathrm{g}\,\mathrm{s}^{-1}$ is obtained, corresponding to an outlet-averaged tritium concentration of approximately $0.90\text{\,}\mathrm{m}\mathrm{g}\,\mathrm{m}^{-3}$. The blanket inventory and outlet flux approach steady state within approximately $30\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}$ (figure 8). This build-up time is of a similar order to the characteristic blanket residence times assumed in system-level analyses of the ARC fuel cycle [35]. The predicted inventory and build-up time differ from earlier ARC-relevant estimates, as expected given differences in geometry, dimensionality, and modelling assumptions across studies [25].

It should be emphasised that these values represent a lower-bound estimate. The homogeneous Dirichlet condition imposed at the inlet corresponds to an idealised, perfectly efficient extraction system; in practice, a finite inlet concentration would raise the equilibrium concentration throughout the breeder and hence increase the total inventory and outlet-averaged concentration. The reported inventory should therefore be interpreted as a conservative minimum for the assumed geometry and operating conditions.

The distribution of tritium concentration exhibits pronounced localisation effects, as shown in figure 6(e). Regions of elevated tritium concentration are strongly correlated with flow stagnation zones, indicating a strong sensitivity of tritium accumulation to the local flow structure within the breeder. In these regions, reduced convective transport leads to localised concentration hot spots. Conversely, areas characterised by high turbulence intensity exhibit enhanced effective diffusion and correspondingly lower tritium concentrations. This behaviour is consistent with the dominant contribution of turbulence-induced diffusion to tritium transport: comparison of the Fickian and effective diffusivity fields (figure 5) shows the latter exceeding the former by several orders of magnitude throughout the breeder. In such regions, increased turbulent mixing dissipates concentration gradients, thereby limiting local tritium accumulation. The interconnect region exhibits the converse effect, as shown in figure 7. FLiBe issuing from the coolant channel into the bulk blanket forms a jet of comparatively low tritium concentration, despite the volumetric tritium generation rate being highest in the channel owing to its proximity to the plasma (figure 3(a)). This apparent contradiction is explained by residence time: the flow velocity in the channel reaches $12\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$, an order of magnitude greater than in the blanket, so the salt traverses the high-generation region too rapidly for appreciable tritium to accumulate before entering the blanket. Combined with the homogeneous Dirichlet condition at the inlet, this produces a persistent low-concentration plume that penetrates well into the blanket before being dispersed by turbulent mixing in the shear layer at the jet boundary.

Concentration in the breeder is therefore governed by the interplay between local production and local residence time, rather than by the production rate alone. These results show that spatially resolved tritium transport simulations can be used to identify unfavourable flow features within the blanket. In particular, regions prone to stagnation or weak mixing may be targeted for geometric optimisation to enhance convective or turbulent transport.

The tritium mass flux at the blanket outlet provides a useful system-level metric for assessing overall blanket performance and tritium throughput. From the transient simulations, an outlet tritium mass flux of approximately $1.0\text{\,}\mathrm{m}\mathrm{g}\,\mathrm{s}^{-1}$ is obtained under steady-state conditions. However, the transient evolution of the outlet flux and blanket inventory, shown in figure 8, provides additional insight beyond the steady-state value alone. The results indicate that the outlet tritium flux approaches its steady-state value over approximately $30\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}$. This timescale is significant when compared with expected plasma pulse durations in fusion devices.

For reference, pulse lengths of order $30\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}$ are anticipated for ITER [41], while no definitive pulse length has yet been published for ARC-class devices. Although ARC-class reactors are generally envisaged to operate with long or quasi-steady-state plasma conditions [44], the predicted transient response time remains highly relevant. If the time required for the blanket to reach steady-state tritium throughput exceeds a single plasma pulse, multiple pulses may be needed for the blanket inventory and outlet flux to approach steady state. In such scenarios, the time for the blanket inventory and outlet flux to reach steady state would be longer than would be inferred solely from steady-state considerations. This transient behaviour has direct implications for tritium system design, as it informs both the time required to reach steady-state tritium throughput and the timescale over which the blanket inventory builds up. In particular, the predicted response time provides guidance on the doubling time of the tritium inventory during reactor start-up and early operation [35].

To verify the FESTIM predictions via code-to-code comparison, the same transport problem was solved independently in OpenFOAM using the passive scalar formulation described in Section 2.3. Since the two codes solve the same governing equations, they are expected to yield comparable results. The temporal evolution of the outlet flux and blanket inventory from both codes is shown in figure 8. The two solutions agree closely: The steady-state outlet flux is $1.0\text{\,}\mathrm{m}\mathrm{g}\,\mathrm{s}^{-1}$ in FESTIM and $1.01\text{\,}\mathrm{m}\mathrm{g}\,\mathrm{s}^{-1}$ in OpenFOAM, with corresponding inventories of $243\text{\,}\mathrm{m}\mathrm{g}$ and $240\text{\,}\mathrm{m}\mathrm{g}$. The outlet flux reaches $99.5\text{\,}\%$ of its steady-state value after $25.0\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}$ in FESTIM and $23.1\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}$ in OpenFOAM; the small difference is most likely attributable to the different time discretisation schemes used by the two codes. This close agreement between an independent finite volume solver and the finite element FESTIM model lends confidence to the predicted inventory and transient behaviour.

At steady state, the FESTIM outlet flux remains approximately $1.5\text{\,}\%$ below the neutronic production rate. By contrast, the finite volume OpenFOAM formulation conserves mass almost exactly, with a steady-state imbalance of only $-0.000\,75\text{\,}\%$, reflecting its inherently conservative construction. As global conservation requires that outlet flux and production balance once a steady state is reached, the FESTIM residual is attributed to the stabilisation scheme: the artificial diffusion term, although necessary to suppress oscillations in the continuous Galerkin formulation at high Péclet numbers, is not strictly conservative. A promising route to remove this discrepancy is to adopt a discontinuous Galerkin discretisation, which is locally conservative by construction and stabilises advection-dominated transport through cell-wise upwind numerical fluxes. Such a scheme would eliminate the need for the artificial diffusion term and its tuning parameter, $\delta$, altogether and is identified as a priority for the workflow’s future development.

From an extraction perspective, the tritium concentration in the outlet fluid is the key parameter, as the effectiveness of many extraction technologies depends directly on the local concentration, amongst other factors [27]. While the outlet flux characterises the overall tritium throughput, the outlet concentration governs the extraction system’s efficiency and operating conditions. The transient results indicate that the outlet concentration approaches its steady-state value on a timescale comparable to that of the outlet flux. This highlights the importance of considering both outlet flux and outlet concentration when evaluating tritium transport and recovery in breeding blanket systems, particularly during reactor start-up or changes in operating conditions.

These considerations become particularly important for reactor concepts operating in pulsed or partially steady-state regimes. Under such conditions, both the magnitude and temporal evolution of tritium concentrations and fluxes influence the performance of the extraction system. Spatially and temporally resolved transport simulations, therefore, provide an essential tool for assessing tritium recovery strategies under realistic operating scenarios, beyond what can be inferred from steady-state metrics alone.

To assess the robustness of the predicted tritium inventories and outlet behaviour, a limited parametric analysis was performed to examine the sensitivity of the results to key numerical and turbulence-related modelling parameters. In particular, variations were considered in both the turbulent Schmidt number, Sc_t, and the numerical tuning parameter, $\delta$. The tuning parameter, $\delta$, controls the magnitude of the artificial diffusion term introduced to ensure numerical stability when solving advection-dominated transport problems at high Péclet numbers. As described in Section 2.3, this tuning term introduces additional diffusion aligned with the local flow direction, analogous to SUPG or artificial diffusion methods, and is used in combination with the turbulent diffusivity to define the effective diffusion coefficient. While necessary for numerical convergence, the tuning term can introduce artificial diffusion and numerical bias. To assess its impact, $\delta$ was varied over a wide range, from 0.01 to 10, around a reference value of $\delta=0.1$. Particular attention was also given to the turbulent Schmidt number, Sc_t, for which very limited data are available for liquid FLiBe flows. Existing studies of turbulent scalar transport report a wide range of Sc_t values, typically between approximately 0.3 and 1.3, depending on flow conditions and fluid properties [38, 46, 48]. For comparison, values of Sc_t around 0.7 are commonly adopted for water flows [38]. In the absence of a well-established value for FLiBe, a range of Sc_t values from 0.3 to 1.3 was therefore investigated to bound the expected physical uncertainty.

The influence of the numerical tuning parameter, $\delta$, is considered first. Reducing $\delta$ reduces artificial diffusion and brings the solution closer to the physically correct behaviour; across the investigated range, it leads to increased tritium retention, together with higher outlet fluxes and concentrations (figure 9). Over this range, however, the resulting changes in inventory, outlet flux, and outlet concentration are modest, indicating that the predictions are relatively robust to this parameter. Reducing $\delta$ towards zero is nonetheless desirable for accuracy but was observed to introduce the first numerical instabilities near the interconnect and in flow stagnation zones, where concentration gradients are highest. This indicates that local mesh refinement in these regions is preferable to relying on artificial diffusion for stability.

The turbulent Schmidt number, $\mathrm{Sc}_{t}$, instead represents a physical uncertainty, as no consensus value exists for liquid FLiBe. Varying $\mathrm{Sc}_{t}$ across the range 0.3 to 1.3 again produces only modest changes in inventory, outlet flux, and outlet concentration (figure 9). The outlet flux in particular is only weakly sensitive to $\mathrm{Sc}_{t}$, indicating that transport at the outlet is dominated by advection, consistent with the high flow velocities in the breeder (figure 4(b)). Turbulent diffusion, therefore, mainly redistributes tritium within the breeder rather than governing its bulk removal.

Despite these insights, several limitations of the present model must be acknowledged. Most notably, the simulations neglect tritium diffusion into and trapping within solid structural materials, which is expected to contribute to the overall tritium inventory. Previous studies on the influence of trapping in structural materials (Eurofer) in a Water Cooled Lead Lithium blanket found a tritium inventory increase of approximately $20\text{\,}\%$ [14]. However, this work neglected the potential impact of neutron damage on tritium retention. This effect was found to increase model inventories in tungsten monoblocks by several orders of magnitude [32, 16]. Therefore, a dedicated study focusing on trapping in structural materials will be required in the future.

In addition, permeation through the inner wall is not explicitly modelled, potentially allowing tritium to migrate from the breeder into plasma-facing components. These effects are independent of the numerical and turbulence-related sensitivities discussed above and may contribute to tritium retention and redistribution. Addressing these mechanisms will be an important focus of future work.

The present results represent an initial step towards fully resolved tritium transport modelling in ARC-class breeding blankets. The analysis could be extended by incorporating a more detailed geometric description, including explicit treatment of solid walls and the associated tritium diffusion and trapping within structural materials. Introducing a nonzero inlet concentration would enable closer coupling with realistic tritium extraction systems and allow direct assessment of extraction performance under transient operating conditions. In addition, accounting for tritium implantation from the plasma into the first wall would provide a more complete representation of tritium sources and sinks within the blanket system. Beyond purely transport-driven effects, future extensions could include chemical reactions within the liquid breeder, allowing investigation of tritium speciation and the influence of redox chemistry on transport and retention.