© W.R. Kendrick and B. Forget, Published by EDP Sciences, 2025

1. Introduction

Microreactors have become a significant source of interest in the nuclear power community, with concepts being proposed and developed for commercial, space, and defense purposes [1, 2]. The designs for these reactors include a variety of fuel types, heat removal mechanisms, and both fast and thermal spectrum specializations. One of the key challenges for such small reactors seeking a thermal spectrum is the need for a very powerful moderator, specifically a moderator with a high macroscopic slowing down power (MDSP).

Table 1 quantifies some of the moderating properties of commonly selected materials. For a microreactor to remain “micro”, keeping the diffusion length short is important to focus thermalization and reduce neutron leakage. This constraint typically eliminates most moderators that do not include hydrogen in their molecule. While water is cheap and well understood, it requires high pressure to achieve adequate thermal conversion efficiency, which necessitates costly pressure vessels. Zirconium hydride, on the other hand, presents an incredibly appealing alternative. Zirconium is inexpensive, easily converted to hydride form, and is neutronically transparent from an absorption standpoint. These same properties are why zirconium has been used in the nuclear power industry for many years as cladding material for fuel [4].

The ability of zirconium to form hydrides is a reflection of the high diffusivity of hydrogen atoms within the zirconium lattice, which also leads to a well documented concern of keeping the hydrogen contained during high temperature operation [5, 6]. Hydrogen leakage presents a neutronic challenge, as less moderation means less thermalization and a potentially lower fission rate, as well as posing a safety concern, since hydrogen on its own is highly combustible. Even if well contained, this high diffusivity means that hydrogen will migrate within the zirconium [7–9]. While likely less concerning than outright moderator loss, strong enough redistribution of the moderator can change the reactor flux profile and impact reactor behavior.

Solving the movement of hydrogen in a reactor core requires solving a nonlinear problem. Neutron flux affects the power profile, which affects the temperature profile, which will affect the hydrogen profile, which will in turn affect neutron flux. In order to simulate the entirety of this feedback loop, this work couples OpenMC (neutron transport) and MOOSE (heat transport and hydrogen migration) as a neutronic-thermal simulation structure to calculate hypothetical hydrogen migration and its resulting effects.

The geometry used to explore this behavior is based on the epithermal reactor modeled by Choi et al. [10]. One defining feature of this reactor geometry is the use of zirconium hydride in a contiguous, monolithic form throughout which heat pipes and fuel pins are embedded. It is more common to see zirconium hydride moderator in pin form, where neutronic feedback is mainly driven by axial migration [11]. In the monolithic configuration, hydrogen can theoretically migrate from one end of the reactor to the other, exhibiting radial and azimuthal mobility that would not occur in pin-form moderators. Since heat transfer in the reactor primarily acts in the radial and azimuthal direction, temperature gradients are likely to develop within the moderator. This phenomenon provides the impetus for this work: to simulate the hydrogen migration in the moderator block and analyze the neutronic-thermal impacts of this redistribution.

This work’s novel results are (1) analyzing a case where hydrogen migration results in a positive reactivity effect for the core, and (2) using functional expansions and continually varying material tracking to perform Monte Carlo particle transport on continually varying hydrogen distributions within the moderator region.

2. Background

The following subsection discusses the physics behind the diffusion of hydrogen in zirconium, followed by a subsection detailing additional particle transport physics that are introduced in order to represent this redistribution.

2.1. Diffusion of hydrogen in zirconium hydride

Diffusion of hydrogen in the zirconium lattice is largely governed by two systems, Fick’s Law and the Soret effect [7]. Essentially, this is diffusion driven by concentration gradient and thermal gradient. Interest in the alpha phase of zirconium hydride for cladding purposes has been a significant source of characterizing hydrogen migration in zirconium. A good example of this is Courtey’s work [7] modeling hydrogen behavior in Zircaloy-4, which is the source for the following equation that describes hydrogen migration in zirconium, Equation (1). Here, J_Fick is hydrogen flux governed by Fick’s law, where C is the hydrogen concentration and D is the diffusion coefficient for hydrogen in the zirconium lattice. J_Soret is the hydrogen flux governed by the Soret effect, where Q^* is the heat of transport, R is the gas constant, and T is the temperature.

$$\begin{array}{c}\hfill J=J_{\mathrm{Fick}}+J_{\mathrm{Soret}}=-D\mathrm{\nabla }C-\frac{DC{Q}^{\ast }}{R{T}^2}\mathrm{\nabla }T\end{array}$$(1)

D is determined via the Arrhenius equation, Equation (2). In this equation, A_D is a pre-exponential factor, R is the ideal gas constant, and Q_D is the activation energy.

$$\begin{array}{c}\hfill D=A_D\times exp(-\frac{Q_D}{\mathit{RT}})\end{array}$$(2)

The heat of transport is assumed to be a value of 25.1 kJ/mol, in agreement with [7]. The specific pre-exponential factor and activation energy chosen were A_D = 1.53 × 10⁻³ cm²/s and Q_D = −14 000 J/mol, based on a similar use in [8, 12].

The partial differential equation solved in order to determine the time evolution of hydrogen concentration is

$$\begin{array}{c}\hfill \frac{\partial C}{\partial t}=\mathrm{\nabla }·J=\mathrm{\nabla }·(-D\mathrm{\nabla }C-\frac{DC{Q}^{\ast }}{R{T}^2}\mathrm{\nabla }T).\end{array}$$(3)

For the majority of the work in this study, the problem is assumed to be steady state, setting the time derivative to zero, and an initial condition of all cells having the same initial hydrogen concentration is assumed. All moderator regions are set to have a Neumann boundary condition of zero, referring to the gradient of hydrogen concentration at that boundary.

It is important to note that zirconium hydride is not a single-phase material, but has multiple phases depending on temperature and hydrogen concentration [13]. In order to simplify the complexity of the problem, no phase transitions will included in the coupled simulation; a single phase will be assumed when setting material properties. Delta phase zirconium hydride is a common phase assumed for moderating purposes [8, 12, 14], and will be what this paper assumes for material properties like density and thermal conductivity.

2.2. Functional expansion and CVMT

One distinct challenge in coupling Monte Carlo neutronics to a finite element-based solver is the difference in physical domains being solved. Monte Carlo neutronics codes typically employ Constructive Solid Geometry (CSG), where surfaces are explicitly defined and volumes are created through Boolean operations with those surfaces. For Monte Carlo codes, the vast majority of tallied statistics are volume-integrated values, where the volume in question will often be a specific cell or material region. In contrast, finite element solvers calculate node-to-node or quadrature point-to-quadrature point, with volume averages being a secondary operation performed by integrating this nodal data.

This difference in the output data format implies that traditional coupling approaches would require passing volume-averaged data between simulations to inform the behavior of each method. Refinement of the quality of the information being transferred requires refining the geometry itself, discretizing the CSG model, making already computationally expensive neutronic solve even costlier [15].

Instead, the use of functional expansions to represent spatial heterogeneity has been suggested as a solution to this challenge [16, 17]. By calculating expansion coefficients that fit orthonormal basis sets to an input, a series of easy-to-calculate, continuous functions can replace the often complicated point-wise data. The choice of function series used for this process strongly impacts the accuracy and efficiency with which the true spatial information can be reproduced. A common choice, and the choice that this work uses, is using Legendre polynomial moments to rebuild the axial variation of hydrogen in the moderator and heat deposition in the fuel.

Tallying the functional expansion moments across the cylindrical geometry is straightforward for both OpenMC and MOOSE and exists in an easily implementable form for users. The term “tallying” here means calculating expansion coefficients that, when applied to their corresponding polynomial basis set, rebuild a fit of the original data that was sampled against. In a Monte Carlo simulation, these expansion coefficients are estimated by

$$\begin{array}{c}\hfill {\widehat{a}}_n=\frac{1}{N}\sum _{i=1}^N{P}_n(z_i)\rho (z_i)\end{array}$$(4)

where ${\widehat{a}}_n$ is the estimated n^th expansion coefficient, N is the number of particles tallied, P_n is the n^th term in the basis set, and ρ is the weighting function. This ρ would be the flux multiplied by a scoring value for the case of an OpenMC simulation, and P_n is the Legendre polynomial series used in this work.

Using the expansion coefficients as input for MOOSE is easily achievable, but in the case of OpenMC, one is required to introduce another complication. MOOSE will return a set of expansion coefficients that rebuild a shape for the number density of hydrogen within a moderator cell. This presents a problem, since in the publicly available version of OpenMC, materials are defined to beconstant in their properties within a cell [18]. This is necessity because Monte Carlo particle transport codes sample particle distance-to-collision by evaluating the cumulative density function (CDF) of the probability of colliding during flight:

$$\begin{array}{c}\hfill F(s)=P_{\text{NC}}(s_{}\text{b})H(s-s_{}\text{b})+{\int }_0^s\frac{d\tau (s)}{\mathit{ds}}{e}^{-\tau (s)ds}\end{array}$$(5)

where s_b is the distance to the boundary, s is the distance along the particle path, and H is the Heaviside step function. P_NC(s_b) is the probability of not colliding until the boundary, calculated via P_NC(s)=e^−τ(s) and τ(s) is the optical depth, defined as

$$\begin{array}{c}\hfill \tau (s)={\int }_0^s{\sigma }_t(s^{\prime })N(s^{\prime })d{s}^{\prime }\end{array}$$(6)

where σ_t is the microscopic total cross-section and N is the number density.

Sampling distance-to-collision requires inverting the CDF in Equation (5), which includes performing an analytical integral along a characteristic line through the material. While Legendre polynomials are simple to evaluate at any point in the axial phase space, they are not easily analytically invertable along this line. Instead, this work utilizes Continually Varying Material Tracking (CVMT), developed by Brown and Martin [19], expounded upon by Griescheimer [17], and first utilized in the context of OpenMC by Ellis [15]. This method is in essence a way of estimating distance-to-collision by performing a series of numerical integrals along the particle flight path.

This represents a significant increase in computational steps during particle transport. What was once a simple evaluation by integrating a constant cross-section across a distance is now potentially a large number of numerical integrals. Beyond that, this also means that track-length estimation is no longer viable for tallying, as it would require a similar numerical integration process to accurately represent scoring across the track-length. Instead, collision-based tally estimators must be used for all tallies, including when tallying functional expansion moments described earlier in this section. While this restriction does not change the tally results, it slightly decreases the statistical convergence rate.

The benefit to implement this method is that it allows for particle transport in materials with spatially varying properties without requiring highly discretized geometry. In the context of this work, CVMT allows for performing particle transport in the moderator material with spatially varying hydrogen concentration.

CVMT is not the only method for handling spatially varying material properties in neutronics solvers. Serpent-2 features an implementation of a method similar to delta tracking, a method that samples collision distances based on a majorant cross-section, thereby avoiding the complexities described earlier in this section [20, 21]. This method was not pursued for this work as calculating the majorant cross-section can be difficult when using continuous distributions of number densities, OpenMC does not currently have a released branch with delta tracking enabled, and because CVMT has been previously implemented in the OpenMC framework [15].

3. Methods

This section discusses some of the more technical aspects of this work, beginning with the two primary codes used during coupled simulation, MOOSE and OpenMC. The coupling scheme that controls the flow of information between both simulations is also discussed.

3.1. OpenMC & MOOSE

OpenMC is utilized extensively for all neutronics simulations described in this work. OpenMC is an open source Monte Carlo particle transport code primarily developed by the Massachusetts Institute of Technology and Argonne National Laboratory [18]. Designed with a focus on parallel performance and featuring an extremely user-friendly Python API, the OpenMC code is very well suited for the work described in this article. Thermal scattering data for zirconium hydride is used for both zirconium and hydrogen scattering, sourced from ENDF-8 cross-sectional data. As mentioned at the end of Section 2.2, the CVMT methodology has been integrated into the transport code, triggering when a neutron enters a cell with a material defined using polynomial coefficients.

Functional expansions are available in OpenMC for tallying as a filter. In order to convert heating rates, in units of eV per source particle, to volumetric heating, each cell score was divided by the cell volume and total heating score of the geometry, and then multiplied by the total input power. Functional expansion tally results from OpenMC must also be orthonormalized before use in MOOSE. This orthonormalization process corrects the expansion coefficients in order to rebuild the original distribution [22].

Additionally, because the azimuthal and radial subdivision of the geometry studied is so dense (see Fig. 3, Direct Accelerated Geometry Monte Carlo (DAGMC) is utilized with OpenMC [23, 24]). This allows for simulation of particle transport directly on a skinned version of the mesh generated for MOOSE. This ensures there are no geometric differences between the OpenMC and MOOSE solves, and allows the radial and azimuthal subdivisions of the moderator to be handled during mesh generation, which is more straightforward than via CSG processes.

For solving heat transfer and hydrogen diffusion, MOOSE is utilized. MOOSE is a finite element physics framework developed by Idaho National Laboratories [25]. It features extremely parallelizable code for high efficiency on clustered computing nodes, as well as well-documented and user-friendly methods of implementing one’s own physics kernels. MOOSE is used for all thermal and hydrogen simulation in this work, relying on modules that are readily available as an open-source download.

The generation of the mesh model of the reactor studied in this work was processed by the Reactor module in MOOSE [26]. This module allows for regularized mesh generation for assemblies and has functionality to cut and stich meshes together, which was heavily leveraged to generate the geometry seen in Figure 3. Heat transfer in the model is generic heat conductance, except at gaps between the heat pipes/fuel pins and the moderator monolith. Heat transfer across the gaps was calculated via a simple radiative heat transfer model, assuming that the gaps are voided.

Hydrogen movement was solved as a generic diffusion problem, with diffusivity calculated using the relation described in Section 2.1. All MOOSE solves for this work are performed with a steady state assumption, except for those mentioned at the end of Section 4.2. As noted by Terlizzi [27], however, performing this hydrogen diffusion problem without a time derivative opens up the solution space to cases where the total hydrogen content is not preserved. This is obviously a large concern, as this would be an unphysical moderator loss and would make isolating the impact of hydrogen migration impossible. To avoid this issue, all hydrogen functional expansion moments are normalized after each MOOSE run to preserve the total number of hydrogen atoms in the system. This normalization does not impact the shape of the hydrogen distribution, only the magnitude, so the normalized steady state results match the transient solution results within a small margin. This is confirmed in Figure 30 where the root mean squared error is reported between all 0^th moments from 3 different steady state and transient MOOSE runs.

Finally, it should be noted again that leakage of hydrogen is not considered in the hydrogen diffusion solve. Hydrogen is an element well known for its ability to penetrate and escape its containment, so zirconium hydride moderator is expected to have a cladding of some kind [28]. While both OpenMC and MOOSE models do not feature any cladding on the moderator, it is assumed that no leakage occurs in order to isolate the effects of hydrogen migration.

3.2. Coupling scheme

Coupling for this process is distinctly separate, with each simulation process occurring stand-alone followed by a step for processing results and converting information into a format readable for the next simulation. The exchange of information and initialization of data containers is handled by a Python script. The general layout for this coupled process can be seen in Figure 1.

Fig. 1. Operational flowchart for coupled neutronic-thermal simulation with hydrogen diffusion.

The cycle begins with an initial parse of the pre-generated mesh. This is done to collect cell-id information for all cells in the mesh, identify what material corresponds to each cell, and tabulate all cell volumes. Following that, the “working” mesh is generated, clearing all block definitions and creating a unique block for every unique volume in the mesh. These processes of loading, parsing, and altering mesh information are handled by Coreform Cubit’s Python API [29].

Before the OpenMC simulation can begin, the mesh file is exported to a facet geometry that OpenMC can read via DAGMC implementation. This facet file replaces the classic Geometry block in OpenMC, so all materials and temperatures are assigned to every cell at this point. After that, the OpenMC solve occurs, tallying cell-wise energy deposition, Legendre polynomials for all fuel pins, and collecting various information like global spectrum and absorption rates.

After the OpenMC step concludes, the Python script loads the information from the output statepoint file, and energy deposition values are converted to volumetric heating rates by performing power normalization. The MOOSE input file is then generated, making sure to assign all material properties to the correct volume ids, and assigning correct heating rates to all volumes, except for fuel pins that will receive Legendre coefficients. The MOOSE run outputs a CSV file with all cell temperatures, a text file with all output hydrogen Legendre coefficients, as well as an output mesh file that contains temperature and hydrogen nodal data. Note that both hydrogen and axial power moments are tallied up to the 20^th order coefficient.

Note that as OpenMC is a Monte Carlo code, all tallied values have associated standard deviations. This includes the functional expansion moments that describe the axial heating rates for all fuel pins. To the best of the authors knowledge at this time, these function expansion standard deviations are unable to be used in MOOSE for propagation of error. This means that there should be some propagated uncertainties for the temperature and hydrogen results from MOOSE that are unable to be quantified for this work.

Convergence of this coupled process is determined by comparing convergence for both OpenMC and MOOSE. Convergence for OpenMC is determined by calculating the normalized root mean square (RMS) of all relative changes in cell-averaged energy deposition compared to the previous run, seen in Equation (7), where E_dep,i is the tallied energy deposition per source particle for cell i. Convergence for MOOSE is similarly based on the RMS of all cell-averaged relative temperature changes compared to the previous run, seen in Equation (8), where T is the cell temperature. After the final MOOSE simulation, if both norms are sufficiently small, the simulation ends; otherwise, a new OpenMC run is prepared with new temperatures and hydrogen coefficients.

$$\begin{array}{cc}\hfill {\epsilon }_{}\text{OpenMC}& =\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(\frac{\mathrm{\Delta }{\text{E}}_{\text{dep},i}}{E_{\text{dep},i}}\right)}^2}\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill {\epsilon }_{}\text{MOOSE}& =\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(\frac{\mathrm{\Delta }{T}_i}{T_i}\right)}^2}\hfill \end{array}$$(8)

4. Results

As first mentioned in Section 1, the geometry studied for this work is based on a heat pipe reactor design produced by the Korea Atomic Energy Research Institute [10, 30, 31], specifically the 5 kW_th epithermal design with cylindrical fuel pins of 1.5 cm diameter, referenced in Table 7 of [10].

This reactor geometry, XY and YZ slices of which are visible in Figure 2, features 19.5% enriched uranium oxide fuel, stainless steel Type-316 (SS316) sheathed potassium heat pipes, a zirconium hydride (ZrH_1.5) moderator block, and is surrounded by a beryllium oxide reflector, both radially and axially. The reactor is modeled as being surrounded by vacuum. The fuel pins and heat pipes are spaced with a pitch of 2.4 cm and both feature 0.005 cm thick gaps between their outer radius and the moderator. The moderator and fuel pins are 35 cm in axial length, and the reflector is 14.385 cm thick in both the radial and axial directions.

Fig. 2. Geometry of the 5 kWth reactor studied, colored by material. In the actual model simulated, heat pipes have a SS316 sheath and potassium inner fluid, they are simplified to one material for this visual. On the left is an XY plane cut, on the right is a YZ plane cut.

Fig. 3. Reactor geometry colored by cell, highlighting the radial and azimuthal discretizations of the zirconium hydride moderator. On the left is an XY plane cut focused around the center of the core, on the right is a YZ plane cut.

The inner wall of the SS316 sheath for the heat pipes are set as a constant temperature boundary of 1000 K, meant to be a similar constraint to that used in [30]. As discussed in Section 3.1, heat transfer between outer SS316 walls and the hydride monolith as well as between fuel pins and the monolith are simulated to be solely via radiation. The emissivities used for the UO₂, ZrH, and SS316 are 0.8, 0.158, and 0.6, respectively [32–34].

While axial variation of hydrogen concentration is approximated with calculated Legendre polynomial moments, radial and azimuthal movement is distinguished via rigorously splitting the zirconium hydride mesh block into many subdomains. These many unique cells, of both ZrH and all other materials, can be seen in Figure 3.

This work analyzes the reactor geometry mentioned above in two different case studies. The initial study is focused on replicating the parameters described in [10], with a reactor thermal power of 5 kilowatts and heat pipes as a boundary condition of 1000 K. The second study takes the same geometry and performs analysis with a range of eleven different thermal powers from 500 W to 50 kW, as well using two different heat pipe temperatures, 700 K and 1000 K. The first study acts as a “realistic” data point for temperature and powers suggested by the work of [10, 30]. The second study expands the input power and temperature parameters in order to explore the phase space for interesting results and to attempt to characterize broader trends.

4.1. 1000 K, 5 kW case study

As described above, this first case study analyzes the impacts of hydrogen migration on steady-state neutronic-thermal behavior when reactor power and heat pipe temperatures are as described in [10], namely 5 kW thermal power and heat pipes fluid held at 1000 K. The initial OpenMC simulation of the reactor geometry provides cell-average heating rates, as well as Legendre moments for the axial heating rates of the fuel pins. Note that the fuel pins are modeled as two concentric cylinders with outer radii of 0.4 cm and 0.75 cm, in order to capture some of the radial shape of heating in the fuel. An XY slice of the cell-average volumetric heating rates, plotted with a logged scale, can be seen in Figure 4. An axial profile of three fuel pins is shown in Figure 5, including both radial layers of each fuel pin.

Fig. 4. XY slice of reactor geometry showing cell-averaged volumetric heating rates in logged scale, produced by the initial OpenMC run before hydrogen migration.

Fig. 5. Reconstruction of axial heating rates of three fuel pins using OpenMC-generated heating moments. The dashed lines are the outer concentric cylinder cell and the solid lines are the inner cell. Starting at the origin in Figure 4 and moving in the positive Y direction, the “inner” pin is the first pin encountered, followed by the “middle” pin, and finally the “outer” pin at the reactor periphery.

Figure 6 is a visual of the normalized fuel pin powers in the reactor, showing a maximum peaking factor of 1.08 and the highest power pins at the outer edges of the reactor. This pattern of normalized powers visibly matches well with Figure 10 from [10], and the magnitudes of heating from Figures 4 and 5 visibly match fairly well with those in Figure 10 of [10]. It should be noted that the eigenvalue calculated by this work differs from that which is reported in Table 7 of [10]. Input material masses and geometry were confirmed to be identical and the power profile aligns well, so the cause of the difference is unknown at this time. Regardless, the original reactor design is meant to be an inspiration for a hypothetical reactor.

Fig. 6. Radial, normalized pin powers for reactor fuel pins, before hydrogen diffusion.

These heating rates are then passed to MOOSE to solve the temperature profile of the reactor. An XY slice showing temperature in the central portion of the core is visible in Figure 7, and Figure 8 shows the temperature profile in the YZ plane. Temperature peaks in the fuel pins at roughly 1025 K, with the coldest region of the core being around the heat pipes at 1000 K. This roughly 25 degree variation differs from the temperature range seen in Figure 4 of [30] which displays a variation of roughly 150 degrees. This difference may come from a different spatial representation of power in the core or a different heat transfer model between fuel pins and the rest of the geometry. Without additional information on the full methodology in [30], one can potentially treat this case study as a more conservative estimation of the hydrogen diffusion, given the dependence of J_soret on ΔT as seen in Equation (1).

Fig. 7. XY slice of the reactor geometry, focused on the inner portion of the core, showing temperature of the MOOSE output mesh.

Fig. 8. YZ slice of the reactor geometry showing temperature of the MOOSE output mesh.

The temperature profiles shown in Figures 7 and 8 inform the diffusion of hydrogen within the zirconium hydride mesh during the MOOSE solve. The cell-averaged hydrogen concentration (x in ZrH_x) can be seen in Figure 9, showing noticeable diffusion to the colder, heat pipe neighboring cells with a peak x of roughly 1.53. In contrast, the cells in the outer areas of the core saw a drop in average concentration from 1.5 to about 1.48. Axially, the hydrogen concentration exhibits a negative cosine shape, with the highest concentrations seen at the axial upper and lower regions. This variation can be seen in Figure 10, where the MOOSE-calculated Legendre moments of three ZrH cells are used to rebuild the estimated axial profile. The cosine shape of heating in the fuel pins (see Fig. 5) create a similar cosine temperature profile, and the hydrogen migrates from high to low temperatures, resulting in the shape seen in Figure 10.

Fig. 9. XY slice of zirconium hydride moderator showing cell-average hydrogen concentration after migration is complete. Areas of greatest concentration are centered around heat pipes.

Fig. 10. Axial profile of hydrogen concentration rebuilt using MOOSE-produced moments. Three regions are selected, the first being adjacent to a fuel pin, the second being adjacent to a heat pipe, and the last being on the outer edge of the monolith.

The hydrogen migration seen in Figures 9 and 10 have a measurable impact on neutron transport in the reactor, most notably increasing the multiplication factor by roughly 250 pcm, as recorded in Table 2. Other neutronic impacts recorded in Table 2 include a decrease in particle leakage by approximately 0.1% and an increase in the thermal fission factor (ThFF) by approximately 0.2%. The increased eigenvalue suggests a higher rate of fission due to hydrogen migration, and the decreased leakage and increased ThFF suggest that this increase in fission rate is due to an increase in thermalization of neutrons. The material-averaged absorption rates in Table 2 support this conclusion, with all materials seeing an increase in absorption rate except for the reflector material, as fewer neutrons leave the center of the core due to increased thermalization. Analysis of a fine-group, global neutron flux tally shows a decrease in the proportion of fast neutrons (greater than 0.1 MeV) of 0.54 ± 0.04%, further supporting this conclusion.

As explained in Section 3.2, providing MOOSE with volume-averaged and Legendre moment-based volumetric heating rates begins with tallying energy deposition per particle in OpenMC. This allows for summarizing material-wise energy deposition rates, which can be seen in Table 3. These energy deposition results show that every material in the reactor experiences more energy deposition per neutron after hydrogen migration, with the fuel and moderator most strongly exhibiting a 0.26% and 0.42% increase, respectively. These are normalized values, however, as is obvious from their units of “per source particle”.

Because reactor power is one of the constraints of the simulation (5 kilowatts thermal) the magnitude of these deposition increases matter less than the relative change in magnitude in comparison to the other materials, as this would indicate some movement of power in the reactor. These values do suggest, however, that per neutron there will be 0.255 ± 0.03 percent more energy deposited in the reactor. This may contribute to a transient increase in reactor power during and after hydrogen migration.

After power normalization, where the previously mentioned energy deposition results are used to spatially distribute assigned reactor power, volumetric heating values can be estimated for each cell. Like Table 3, material average heating rates and changes in heating rate due to hydrogen migration can be calculated, and are seen in Table 4. Looking at these tabulated statistics, what is most evident is that outside of a small increase in heating rate in the moderator and a decrease in the reflector, power remains concentrated in its original locations.

While there is no dramatic change in material-average heating rates, the fuel pins do see a slight change in spatial distribution of power. This change in power is small and can be difficult to distinguish from stochastic noise, but averaging the change in fuel heating rate by “layer” (where each layer represents a unique pin position), a shift can be seen. Figure 11 shows this behavior. The outer-most fuel pins decrease in heating rate by an average of 0.37 ± 0.08 percent, while the inner-most pins increase by an average of 0.14 ± 0.08 percent. This movement in power is opposite to the power shape in Figure 6, helping reduce power peaking in the core.

Fig. 11. Percent changes in fuel pin heating rate, averaged by layer.

This shift in power is interesting to note and will be a continued trend seen in Section 4.2, but the absolute values of this power shift are still very low, and therefore do not have a significant impact on thermal results. An XY slice of the inner reactor showing change in temperature due hydrogen migration can be seen in Figure 12. The lack of a clear pattern in temperature variation indicates that the impact of hydrogen migration is so low for this case that stochastic noise becomes an influencing factor. Additionally, the magnitudes of 0.3 degree difference at most are well below noticeable levels during operation.

Fig. 12. Slice in XY plane showing change in predicted temperature after hydrogen migration for the moderator and fuel. The lack of clear pattern suggests that these deltas are likely influenced by stochastic or numerical noise.

4.2. Multi-temperature, multi-power case study

As mentioned in Section 4.1 and recognizable from Equation (1), the migration of hydrogen in this system is driven both by the temperature of the medium and by the gradient of temperature. Because of this, the case study described in this section modulates both the thermal power simulated as well as the heat pipe temperature boundary condition. This provides a series of data points to characterize this reactor geometry’s neutronic and thermal response to hydrogen migration driven under varying temperatures and temperature gradients.

The eleven powers simulated range from 500 watts to 50 kilowatts, and the two heat pipe temperatures used are 700 and 1000 Kelvin. The difference between minimum and maxiumum cell-averaged temperature for the moderator block at all cases can be seen in Figure 13. Both heat pipe cases experience a nearly identical delta in temperature, and the distribution of temperatures looks similar to Figure 7, just with new maximums and minimums.

Fig. 13. Difference between maximum and minimum cell-average temperature reported for zirconium hydride. Includes different reactor powers and two different heat pipe temperatures as boundary conditions. The vertical dashed line indicates the power prescribed by [10].

The resulting maximum and minimum x in ZrH_x can be seen in Figure 14. The overall pattern of distribution remains similar to the one seen in Figure 9, just with the range of values scaled to new maximums and minimums.

Fig. 14. Minimum and maximum cell-average stoichiometric ratio of hydrogen in zirconium for different reactor powers and two different heat pipe temperatures. The vertical dashed line indicates the power prescribed by [10].

Note how the 700 K case sees more dramatic migration compared to the 1000 K case. This may seem counterintuitive, given the equation for the diffusion coefficient (see Sect. 2.1) is such that lower temperatures will have lower coefficient values. Interestingly, migration is greater in the 700 K case due to the T² value in the denominator of Equation (1). For the results presented in this section, it seems that the 1/T² contribution outweighs the reduced diffusion coefficient, resulting in stronger diffusion for the lower temperature cases. Additionally, note that the 700 K/50 kW case has a x value that becomes quite large (exceeding 2 in some regions) which cannot be properly modeled by our hydrogen model as these values would lead to significant phase change in the material and additional constraints on the migration.

Real-life zirconium hydride held under these temperatures would experience phase changes that make the earlier assumption of delta-phase behavior incorrect. Because this work continues to use delta-phase properties, the reader should consider these results as estimations of potential impact of hydrogen migration, and expect that higher accuracy phase-change models may show different hydrogen responses, particularly at the transition from delta to epsilon phase [35].

Figure 15 shows the reactivity effect of hydrogen migration for all set powers and heat pipe temperatures. Like the results from Section 4.1, hydrogen migration has a positive effect on core multiplication factor that increases as power increases and temperature gradients increase. Notably, the 700 K case reports a reactivity roughly 1500 pcm higher than the 1000 K case when simulated at 50 kW of power. The reactivity effect per power is shown in Figure 16. After some initial variation at low temperatures, the relationship between power and reactivity seems to become nearly constant, with the 700 K case seeing between 0.1 and 0.075 pcm effect per watt.

Fig. 15. Reactivity effect of hydrogen migration in the zirconium hydride at different reactor powers and two different heat pipe temperatures. Reactivity is calculated via (k2 − k1)/k2k1 where k2 is the eigenvalue in the case with hydrogen migration, and k1 is the eigenvalue in the case where hydrogen concentration is unchanged. The vertical dashed line indicates the power prescribed by [10].

Fig. 16. Reactivity per power plotted for different reactor powers and two different heat pipe temperatures. Values for reactivity come from Figure 15. The vertical dashed line indicates the power prescribed by [10].

Once again, these reactivity effects are driven by increased thermalization in the core centered around the heat pipes. Figures 17 and 18 show the neutronic spectral shifts due to hydrogen migration. With the strong decrease in fast flux comes a reduction in particle leakage for the core, seen in Figure 19, peaking at a reduction of 1.6% for the 700 K case and 0.8% for the 1000 K case.

Fig. 17. Percent difference in global neutron flux due to hydrogen migration at three different reactor powers and one heat pipe temperature. Shows clear decrease in fast flux as higher power drives stronger hydrogen migration which increases thermalization.

Fig. 18. Percent change in global fast flux (neutrons with energy greater than 0.1 MeV) due to hydrogen migration, shown at different reactor powers and two different heat pipe temperatures. This shows clear increasing thermalization as reactor power increases and hydrogen migration becomes more pronounced. The vertical dashed line indicates the power prescribed by [10].

Fig. 19. Percent change in OpenMC-reported leakage due to hydrogen migration at different reactor powers and two different heat pipe temperatures. The vertical dashed line indicates the power prescribed by [10].

As hydrogen migrates inwards in the moderator region and results in a more thermal flux, the fission rate of the core increases can be seen in Figure 20. The 700 K case experiences a roughly 4% higher likelihood of fissioning per neutron after hydrogen migration. Similarly to this behavior, with neutrons in the core becoming slower on average, the amount of total energy deposited per neutron in the core increases. By taking the total change in energy deposited per neutron Figure 21 can be generated, showing how much the fission rate would theoretically need to be reduced to account for this increase in energy.

Fig. 20. Percent change in fission rate due to hydrogen migration at different reactor powers and two different heat pipe temperatures. The vertical dotted line indicates the power prescribed by [10].

Fig. 21. The required change in fission rate to maintain reactor power after the effects of hydrogen migration, based on the change in energy deposition per neutron/photon. Values are for different reactor powers and two different heat pipe temperatures. The vertical dotted line indicates the power prescribed by [10].

After converting energy deposition tallies to heating rates during power normalization, analysis on how power has spatially shifted in the core can be done. Figure 22 is a plot of percent changes in material-wise heating rates. For all cases, the amount of power located in the fuel material remains essentially constant, while the moderator experiences an increase in heating rate and the reflector experiences a decrease in heating rate.

Fig. 22. Percent change in material-average heating rates due to hydrogen migration at different powers. The dashed lines represent the case where heat pipes are held at 700 K, the solid lines represent the case where heat pipes are held at 1000 K. The vertical dotted line indicates the power prescribed by [10].

By isolating a single power and heat pipe temperature, namely the “worst case” 700 K/50 kW case, a closer look at spatial shift in power is possible. Figure 23 shows a cell-averaged plot of percent change in heating rate, where a noticeable shift inwards in power is visible. This movement follows the movement of hydrogen, and is another result of the increased thermalization in the center of the core. An axial shift in power for three fuel pins can be seen in Figure 24, showing that the hydrogen profiles seen in Figure 5 can be impactful on fuel axial heating profiles when migration is exacerbated under higher temperature gradients. Note that the percent differences are based on the axial region, so a hypothetical 3% increase at one of the ends of the fuel rod would not be equal to a 3% decrease in one of the central axial regions.

Fig. 23. XY slice focused on the inner core region showing percent change in cell-average heating rates due to hydrogen migration. These results are from the case simulated with 700 K heat pipes and 50 kW reactor power, the cae with the stronger hydrogen migration modeled. Values with relative standard deviations above 50% were ignored.

Fig. 24. Axially averaged percent changes in fuel pin heat rates for the same three pins selected in Figure 5. These results are from the case simulated with 700 K heat pipes and 50 kW reactor power, the case that has the strongest hydrogen migration modeled.

This inward shift in power can be quantified as a change in pin powers, seen in Figure 25 where the 700 K/50 kW case results are shown. The strongest difference seen is a decrease in pin power at the outer-most pins. These strongest changes in pin powers are plotted for all powers and heat pipe temperatures in Figure 26. This shift helps to flatten the power profile in the core since the highest power fuel pins are located at the outer edges, seen in Figure 6. A plot of percent change to the pin power Peaking Factor (the ratio of the highest power fuel pin to the average fuel pin power) can be seen in Figure 27.

Fig. 25. Change in normalized pin power due to hydrogen diffusion for the 700 K/50 kW case. Increased thermalization in the center of the core causes power to move inwards. Note that pin powers are normalized such that the average fuel pin power is 1.

Fig. 26. Maximum values of the percent change in pin power due to hydrogen migration at different reactor powers and two different heat pipe temperatures. Note that these are absolute percent changes, and most often the highest percent change is a decrease in pin power at the outer fuel pins, as seen in Figure 25. The vertical dotted line indicates the power prescribed by [10].

Fig. 27. Percent difference in pin power peaking factor due to hydrogen migration at different reactor powers and two different heat pipe temperatures. The vertical dotted line indicates the power prescribed by [10].

The stronger temperature gradients modeled for this case study cause stronger hydrogen migration which in turn has a stronger resulting effect on reactor temperature profile. The strongest case of migration, the 700 K/50 kW case, causes noticeable change in temperature in the core, seen in Figures 28 and 29. The outer fuel pins experience a roughly 8 degree decrease in temperature, while the inner fuel pins experience a roughly 6 degree increase in temperature, specifically at the upper and lower axial regions, following the trends seen in Figure 24.

Fig. 28. XY slice focused on the inner core region showing the change in predicted temperature due to hydrogen migration for the 700 K/50 kW case.

Fig. 29. YZ slice of the core region showing the change in predicted temperature due to hydrogen migration for the 700 K/50 kW case.

Fig. 30. Root mean square error of hydrogen 0th moments over time at three different powers with heat pipes at 1000 K. RMSE is calculated via ${\text{RMSE}}_t=\sqrt{\left(\sum _N{({\mu }_{n,t}^0-{\mu }_{}^0\text{n,ss})}^2\right)/N}$ where μn, t0 is the 0th axial hydrogen moment for cell n at timestep t, μn,ss0 is the same value but from the steady state solution, and N is the total number of ZrH cells.

Finally, because all results presented so far in this work assume steady-state behavior, it is important to estimate the speed at which this hydrogen migration occurs. In order to do so, the MOOSE calculation was changed to transient mode and a time derivative was added to the hydrogen system of equations. Tested on three different powers with 1000 K heat pipe temperatures, the root mean squared error of the cell-average hydrogen concentrations is plotted over time in Figure 30. The 500 W, 5 kW, and 50 kW cases took respectively 5.8 days, 8.6 days, and 29.2 days to reach steady-state behavior. These values are in a similar range to the 50–200 hours reported for full diffusion in [8].

Taking all these results into account, the hydrogen migration for this core tends to have more of a global thermalization impact. The concentration of hydrogen in the center of the core around the heat pipes results in the neutron spectrum softening, leading to more energy deposited in the core and in general more reactions per neutron, including fissions. This can have a significant effect on eigenvalue as temperature gradients in the moderator increase. There is not a significant change to the spatial distribution of these now more thermal neutrons, however. The power remains concentrated in the fuel material and the maximum shift in pin power seen is around 3.5%, an inward movement that has a minor effect in flattening the power shape.

5. Conclusions

To summarize the results presented in Section 4, particularly Section 4.1, temperature gradients in the moderator material drive hydrogen migration towards the colder regions in the core, specifically the heat pipes. This concentration of hydrogen around the heat pipes results in increased thermalization of the neutron flux. This slower neutron population escapes the core less often, deposits more of its energy in the core, and is more likely to interact with core material including via fission. In the most extreme case (the 700 K/50 kW case), hydrogen migration has a reactivity effect of over 3500 pcm and increases fission rate by roughly 4%. While reactor eigenvalue seems to be relatively responsive to hydrogen migration, spatial effects on heating rates are mostly minor after power normalization has taken place. In the most extreme case identified, pin powers change by a maximum of roughly 3% and the maximum temperature difference caused by the hydrogen migration was 8 degrees Kelvin. In the case described in Section 4.1 where power and heat pipe temperature are 5 kW and 1000 K, the impacts of hydrogen migration are mild. Eigenvalue increases by around 250 pcm, but no meaningful spatial shifts in power or impacts to temperature are recorded.

One of the critical takeaways from this work is that the impact of hydrogen migration on neutronic-thermal behavior is highly dependent on core geometry, particularly heat pipe placement. In the geometry studied, heat pipes located at the center of the core cause migrating hydrogen to concentrate centrally. This central concentration is what drives a positive eigenvalue impact and the flattening behavior seen in Figures 25 and 27. If heat pipes were positioned on the outer periphery surrounding fuel pins, hydrogen migration would likely resemble moderator loss from the core center, resulting in a negative eigenvalue impact and an increased peaking factor.

Contemporary studies of hydrogen migration with neutronic-thermal simulation indeed report negative reactivity effects [11, 12, 27]. These negative reactivities stem from separated moderator regions rather than contiguous moderator blocks, which isolate hydrogen migration primarily to the axial direction (producing negative reactivity), and from heat removal locations that create temperature gradients moving hydrogen away from the core center. This study is unique in reporting a positive reactivity effect when using a contiguous moderator block core geometry.

When using zirconium hydride as a moderator for a nuclear reactor, the primary concern has been and should remain on characterizing the loss of hydrogen from the moderator. The authors believe that the results of this work show that when using zirconium hydride as a contiguous block there should also be attention paid to a potential redistribution of hydrogen. Hydrogen diffusion is largely beholden to temperature gradient and heat moves radially in most reactor designs (going from fuel pin to heat pipe, for example). By using a contiguous block moderator in the design, there is potential that the radial heat movement will manifest as a strong temperature gradient in the moderator, and in turn drive significant hydrogen migration. This work shows that hydrogen migration can change core thermalization behavior and significantly affect multiplication factor, as well as shift power profile under more extreme hydrogen migration. When studying the leakage of hydrogen from the moderator, it may be important to take this intra-lattice migration of hydrogen into account. The leakage rate from the surface of uniform ZrH_1.5 is likely to differ from the leakage rate at a surface that has a diffusion-driven local value of ZrH_1.9, for example.

There are many avenues for improvement that future work on this topic could take. One of the most obvious tasks is to improve the accuracy of the hydrogen diffusion coefficients and account for additional constraints that may occur at the boundary of phase changes. For this work, a single pre-exponential factor and activation energy was used for the Arrhenius equation in Equation (2). These values, also used in [8, 12], come from [9] who fit these constants for ZrH_1.58 at a temperature range of 600–970 K. More accurate diffusion behavior would be possible to model by adding concentration dependence to the Arrhenius equation variables.

The next avenue for further work could be to improve the thermodynamic accuracy of the model. Because thermodiffusion of hydrogen in the zirconium relies on temperature gradients, confidence in the accuracy of hydrogen migration is directly tied to confidence in the accuracy of the temperature gradients. Some areas of improvement for this work would be in adding temperature-dependent thermal conductivity values for the reactor materials, and using a heat pipe model to replace the constant temperature boundary conditions.

One area for future work that involves a large change in methodology would be using Cardinal [36] to couple neutronic response directly during MOOSE simulation. Cardinal includes OpenMC transport simulation in directly in the MOOSE framework, and features adaptive meshing capability in order to capture material property gradients (e.g. temperature). Having neutronic feedback included in the MOOSE framework would simplify many of the complicated data-processing and transfers described in Section 3.2. Some downsides may be that strong concentration gradients will require significant mesh discretization, which would be computationally costly during particle transport.

The final area for future work could be focused on simulating this behavior in a fully coupled transient process. Figure 30 indicates that hydrogen migration would likely occur over the course of multiple days, not a dramatically fast process. There still may be value in coupling this migration with a systems model, as perhaps this behavior will impact reactor operational regimes. This would also provide a good opportunity to include a model for hydrogen leakage, and measuring the impact of hydrogen migration on leakage rate.

Acknowledgments

Thank you to Guillaume Giudicelli for his guidance with MOOSE. This research made use of the resources of the High Performance Computing Center at Idaho National Laboratory, which is supported by the Office of Nuclear Energy of the U.S. Department of Energy and the Nuclear Science User Facilities under Contract No. DE-AC07-05ID14517.

Funding

This research was funded by the Department of Nuclear Science and Engineering at the Massachusetts Institute of Technology.

Conflicts of interest

The authors have nothing to disclose, nor any conflicts of interest.

Data availability statement

The data included in this article can be accessed by contacting the primary author.

Author contribution statement

Conceptualization, W.R. Kendrick and B. Forget; Methodology and Investigation, W.R. Kendrick; Supervision, B. Forget.

References

All Tables

All Figures

	Fig. 1. Operational flowchart for coupled neutronic-thermal simulation with hydrogen diffusion.
In the text

	Fig. 2. Geometry of the 5 kWth reactor studied, colored by material. In the actual model simulated, heat pipes have a SS316 sheath and potassium inner fluid, they are simplified to one material for this visual. On the left is an XY plane cut, on the right is a YZ plane cut.
In the text

	Fig. 3. Reactor geometry colored by cell, highlighting the radial and azimuthal discretizations of the zirconium hydride moderator. On the left is an XY plane cut focused around the center of the core, on the right is a YZ plane cut.
In the text

	Fig. 4. XY slice of reactor geometry showing cell-averaged volumetric heating rates in logged scale, produced by the initial OpenMC run before hydrogen migration.
In the text

Fig. 5. Reconstruction of axial heating rates of three fuel pins using OpenMC-generated heating moments. The dashed lines are the outer concentric cylinder cell and the solid lines are the inner cell. Starting at the origin in Figure 4 and moving in the positive Y direction, the “inner” pin is the first pin encountered, followed by the “middle” pin, and finally the “outer” pin at the reactor periphery.

In the text

	Fig. 6. Radial, normalized pin powers for reactor fuel pins, before hydrogen diffusion.
In the text

	Fig. 7. XY slice of the reactor geometry, focused on the inner portion of the core, showing temperature of the MOOSE output mesh.
In the text

	Fig. 8. YZ slice of the reactor geometry showing temperature of the MOOSE output mesh.
In the text

	Fig. 9. XY slice of zirconium hydride moderator showing cell-average hydrogen concentration after migration is complete. Areas of greatest concentration are centered around heat pipes.
In the text

	Fig. 10. Axial profile of hydrogen concentration rebuilt using MOOSE-produced moments. Three regions are selected, the first being adjacent to a fuel pin, the second being adjacent to a heat pipe, and the last being on the outer edge of the monolith.
In the text

	Fig. 11. Percent changes in fuel pin heating rate, averaged by layer.
In the text

	Fig. 12. Slice in XY plane showing change in predicted temperature after hydrogen migration for the moderator and fuel. The lack of clear pattern suggests that these deltas are likely influenced by stochastic or numerical noise.
In the text

	Fig. 13. Difference between maximum and minimum cell-average temperature reported for zirconium hydride. Includes different reactor powers and two different heat pipe temperatures as boundary conditions. The vertical dashed line indicates the power prescribed by [10].
In the text

	Fig. 14. Minimum and maximum cell-average stoichiometric ratio of hydrogen in zirconium for different reactor powers and two different heat pipe temperatures. The vertical dashed line indicates the power prescribed by [10].
In the text

	Fig. 15. Reactivity effect of hydrogen migration in the zirconium hydride at different reactor powers and two different heat pipe temperatures. Reactivity is calculated via (k2 − k1)/k2k1 where k2 is the eigenvalue in the case with hydrogen migration, and k1 is the eigenvalue in the case where hydrogen concentration is unchanged. The vertical dashed line indicates the power prescribed by [10].
In the text

	Fig. 16. Reactivity per power plotted for different reactor powers and two different heat pipe temperatures. Values for reactivity come from Figure 15. The vertical dashed line indicates the power prescribed by [10].
In the text

	Fig. 17. Percent difference in global neutron flux due to hydrogen migration at three different reactor powers and one heat pipe temperature. Shows clear decrease in fast flux as higher power drives stronger hydrogen migration which increases thermalization.
In the text

	Fig. 18. Percent change in global fast flux (neutrons with energy greater than 0.1 MeV) due to hydrogen migration, shown at different reactor powers and two different heat pipe temperatures. This shows clear increasing thermalization as reactor power increases and hydrogen migration becomes more pronounced. The vertical dashed line indicates the power prescribed by [10].
In the text

	Fig. 19. Percent change in OpenMC-reported leakage due to hydrogen migration at different reactor powers and two different heat pipe temperatures. The vertical dashed line indicates the power prescribed by [10].
In the text

	Fig. 20. Percent change in fission rate due to hydrogen migration at different reactor powers and two different heat pipe temperatures. The vertical dotted line indicates the power prescribed by [10].
In the text

	Fig. 21. The required change in fission rate to maintain reactor power after the effects of hydrogen migration, based on the change in energy deposition per neutron/photon. Values are for different reactor powers and two different heat pipe temperatures. The vertical dotted line indicates the power prescribed by [10].
In the text

	Fig. 22. Percent change in material-average heating rates due to hydrogen migration at different powers. The dashed lines represent the case where heat pipes are held at 700 K, the solid lines represent the case where heat pipes are held at 1000 K. The vertical dotted line indicates the power prescribed by [10].
In the text

	Fig. 23. XY slice focused on the inner core region showing percent change in cell-average heating rates due to hydrogen migration. These results are from the case simulated with 700 K heat pipes and 50 kW reactor power, the cae with the stronger hydrogen migration modeled. Values with relative standard deviations above 50% were ignored.
In the text

	Fig. 24. Axially averaged percent changes in fuel pin heat rates for the same three pins selected in Figure 5. These results are from the case simulated with 700 K heat pipes and 50 kW reactor power, the case that has the strongest hydrogen migration modeled.
In the text

	Fig. 25. Change in normalized pin power due to hydrogen diffusion for the 700 K/50 kW case. Increased thermalization in the center of the core causes power to move inwards. Note that pin powers are normalized such that the average fuel pin power is 1.
In the text

	Fig. 26. Maximum values of the percent change in pin power due to hydrogen migration at different reactor powers and two different heat pipe temperatures. Note that these are absolute percent changes, and most often the highest percent change is a decrease in pin power at the outer fuel pins, as seen in Figure 25. The vertical dotted line indicates the power prescribed by [10].
In the text

	Fig. 27. Percent difference in pin power peaking factor due to hydrogen migration at different reactor powers and two different heat pipe temperatures. The vertical dotted line indicates the power prescribed by [10].
In the text

	Fig. 28. XY slice focused on the inner core region showing the change in predicted temperature due to hydrogen migration for the 700 K/50 kW case.
In the text

	Fig. 29. YZ slice of the core region showing the change in predicted temperature due to hydrogen migration for the 700 K/50 kW case.
In the text

Fig. 30. Root mean square error of hydrogen 0th moments over time at three different powers with heat pipes at 1000 K. RMSE is calculated via ${\text{RMSE}}_t=\sqrt{\left(\sum _N{({\mu }_{n,t}^0-{\mu }_{}^0\text{n,ss})}^2\right)/N}$ where μn, t0 is the 0th axial hydrogen moment for cell n at timestep t, μn,ss0 is the same value but from the steady state solution, and N is the total number of ZrH cells.

In the text