Pebble dynamics and thermal-fluid analysis of high-temperature gas-cooled pebble bed reactors using DEM and CFD simulations

This study presents a multiphysics computational simulation framework for analyzing pebble dynamics and thermal-fluid behavior in High-Temperature Gas-Cooled Pebble Bed Reactors (HTG-PBR). The pebble circulation and intermixing effects are predicted using Discrete Element Method (DEM) implemented in LIGGGHTS, while the thermal-fluid behavior is simulated with computational fluid dynamics (CFD) in OpenFOAM. The CFD model employs a porous-media formulation with a local thermal non-equilibrium model to capture the energy exchange between the helium coolant and pebbles. Integrating the DEM-based mixing effects into the porous CFD model enables a more physically representative and scalable approach for full-core reactor analysis. Both DEM and CFD solvers are validated using established pebble-bed benchmark problems to confirm the viability of the developed computational models. A HTG-PBR-like conical model reactor is employed as a test problem to evaluate the developed method. The simulation results confirm the predictive capability of the developed models for HTG-PBR performance analysis and provide insight for future multiphysics coupling strategies for reactor design optimization.

1 Introduction

High-Temperature Gas-Cooled Pebble Bed Reactors (HTG-PBR) are advanced nuclear reactors that utilize small spherical fuel elements (pebbles) packed within a reactor vessel. These reactors rely on helium gas as the coolant, which flows through the core and interacts thermally with the pebble bed to remove the heat generated in the reactor. Therefore, efficient heat transfer is vital for reactor operational safety and performance. However, accurately predicting thermal-hydraulic behavior, particularly the axial and radial temperature distributions, remains a critical challenge for HTG-PBR. This is because the complex interactions between coolant flow and randomly packed pebbles result in significant local temperature variations, which must be well understood to enhance reactor safety and efficiency (Zou et al., 2022). To address this challenge, high-fidelity computational simulations are essential for resolving the detailed temperature fields and heat transfer mechanisms inside the reactor core.

On the other hand, the coupled multiphysics nature of HTG-PBR—encompassing thermal-fluid dynamics, neutronics, and pebble motion—makes simulation and analysis highly complex and computationally expensive. To account for these effects, both resolved and unresolved thermal-hydraulics modeling approaches (Lee et al., 2007; Janse van Rensburg and Kleingeld, 2011; Ge et al., 2016; Wu et al., 2018; Wu et al., 2019; Sharma et al., 2020; Mulder and Boyes, 2020; Ahmed et al., 2021; Jaradat et al., 2023; Latifi et al., 2024) have been employed to study thermal-fluid behavior in HTG-PBR. While resolved approaches can capture detailed local phenomena such as flow channeling, contact resistance, and hot spots around individual pebbles, they become computationally infeasible for full-core analysis involving hundreds of thousands of pebbles. To surpass these limitations, Unresolved approaches like porous-media models are widely adopted owing to their significantly reduced computational demand. A comparative study by Sharma et al. (Sharma et al., 2020) showed that both resolved and unresolved approaches yield nearly identical predictions for global heat transfer characteristics. However, most porous-media studies (Wu et al., 2018; Jaradat et al., 2023; Latifi et al., 2024) assume spatially uniform porosity and impose a constant heat source across the domain—an oversimplification of actual reactor conditions. In practice, both porosity and heat generation vary spatially and temporally due to the stochastic nature of pebble packing, continuous fuel addition, spent fuel removal, and temperature-dependent nuclear reactions.

Additionally, although subregion-based neutronic modeling has been recognized as an efficient and practical method for full-scale HTG-PBR simulations, most existing neutronics studies neglect pebble intermixing between subregions and assume simplified or random mixing of pebbles. To address this gap, the present study enhances the fidelity of such models by incorporating a DEM-based simulation to explicitly capture pebble circulation and intermixing effects in the core. Moreover, the present work uses an unresolved porous-media model in OpenFOAM (Weller et al., 1998) for thermal-fluid analysis, offering a more realistic and computationally feasible approach for simulating HTG-PBR behavior under near-operational conditions.

To improve resolution, researchers have increasingly adopted the CFD-DEM coupling method for pebble bed simulations (Wu et al., 2018). By integrating CFD-DEM, the models enhance the understanding of heat transfer mechanisms, particularly in characterizing velocity profiles and temperature gradients, within the reactor (Janse van Rensburg and Kleingeld, 2011). This CFD-DEM coupling approach provides a more detailed representation of fluid-solid interactions, bridging the gap between purely continuum-based and particle-based models. During the preparation of this paper, we noticed Li et al. (2025) recently published a high-fidelity simulation work of HTG-PBR that integrates OpenMC (Romano et al., 2015) for neutronics, OpenFOAM (Weller et al., 1998) for thermal-fluid analysis, and LIGGGHTS (Kloss et al., 2012) for DEM-based pebble dynamics. Their findings demonstrated that pebble reshuffling and temperature-dependent neutronic behavior significantly influenced local power and burnup distributions.

To enable efficient neutronic analysis of HTG-PBR, our earlier work (Lu and Wu, 2022) proposed a computationally economical strategy to approximate the equilibrium state using a multi-fuel region model. In this approach, pebble motion and fresh fuel loading were not explicitly simulated. Instead, burnup values were adjusted based on a parabolic velocity profile and modified according to the effective multiplication factor of the reactor. Stewart et al. (2024) introduced a grouping-based approach in which pebbles were explicitly tracked and replaced upon reaching a burnup threshold or completing two circulations. Meanwhile, Robert et al. (2023) performed a hyper-fidelity simulation that computed burnup for each individual pebble. However, the associated computational cost makes such an approach impractical for full-core, long-term analysis. These efforts underscore a key modeling gap: the lack of accurate representation of pebble intermixing during pebble circulation. As a result, most previous studies (Stewart et al., 2024; Fratoni and Greenspan, 2010; Wu et al., 2013; Jiang et al., 2012; Reger et al., 2023; Tano and Ragusa, 2021) relied on assumed or randomized mixing models, which may not reflect actual pebble motions in HTG-PBR.

In this study, we aim to bridge this gap by conducting two separate analyses: (1) a pebble dynamics study using DEM simulations to predict pebble circulation and quantify intermixing across subregions, and (2) a thermal-fluid analysis that employs the unresolved porous-media approach in OpenFOAM to predict the temperature profiles of key components in the reactor. Through these analyses, we ensure both computational tractability and physical realism in capturing the dynamic and thermal behavior of the pebbles under representative HTG-PBR operating conditions.

In this work, the pebble dynamics are analyzed using LIGGGHTS, a DEM-based solver, to simulate realistic packing structure, circulation, and intermixing behavior of the pebbles under gravitational flow. To facilitate the analysis, the reactor core is divided into axial and radial subregions and the intermixing between these subregions is quantified based on pebble movement. This allows for a more physically representative assessment of pebble reshuffling, which is critical for evaluating fuel behavior over time. The thermal-fluid simulation is conducted using a CFD solver in OpenFOAM with an unresolved porous-media model. The governing equations for mass, momentum, and energy conservation are solved using a local thermal non-equilibrium (LTNE) model, where energy equations for the solid (pebbles) and fluid (helium) phases are treated separately. The model encompasses convection and radiation heat transfer, as well as conduction within the pebble phase, to ensure accurate thermal predictions across the reactor core.

Before applying the developed approaches to a model reactor problem, a couple of model verification and validation efforts are carried out on the DEM and CFD solvers using established pebble-bed type benchmark problems. The pebble dynamics simulation is validated against the experimental work of Jiang et al. (2012), which utilized a test vessel with similar geometry to HTR-10 design (Wu et al., 2013). Reger et al. (2023) extracted radial velocity profiles at various axial positions from Jiang’s experiments in a related study, which is also used here for benchmarking the DEM predictions. The thermal-fluid solver is validated using data from the SANA experimental facility, a well-known benchmark for natural convection and radiative heat transfer in packed beds (Stocker and Nieben, 1996; Zou et al., 2017).

While this study outlines a generic multiphysics coupled computational framework for HTG-PBR simulation, discussing potential coupling among DEM, CFD, and neutronics, the actual implementation in this work is limited to coupling DEM and CFD modules. Neutronics modeling, although discussed among the context, is not fully integrated into the models used in this study. Specifically, this study primarily focused on modeling efforts that accurately capture pebble dynamics and thermal-fluid behavior in the reactor by integrating DEM-based pebble mixing into a porous-media CFD model. In contrast to most traditional modeling approaches that rely on static core representations or simplified random pebble mixing assumptions, the present approach introduces a realistic pebble grouping and intermixing strategy based on DEM simulations. This approach replaces commonly assumed static or random mixing strategies with physics-based migration paths, capturing localized asymmetries in pebble circulation that affect power distribution and burnup modeling. In addition, the presented modeling approaches were validated part by part against experimental data, providing confidence in their predictive capabilities. The predictions of dynamic fuel movement and thermal fluid behavior on a model reactor problem using the developed approaches distinguishes this work from existing models and justifies them as a viable tool for design analysis, performance optimization, and safety evaluation of next-generation HTG-PBR systems.

The remainder of this paper is organized as follows. Section 2 presents theoretical details of the computational modeling approaches with the emphasis given to the mathematical models of DEM and CFD simulations. Section 3 summarizes the model verification and validation efforts, where both DEM and CFD solvers developed in this work are benchmarked against experimental and reference data to ensure accuracy and feasibility. Section 4 presents some results and discussion of pebble intermixing and thermal-fluid analysis on one model reactor problem that is akin to typical HTG-PBR. Section 5 offers some concluding remarks drawn from the study, summarizing the key findings and future research directions on high-fidelity multiphysics computational modeling of HTG-PBR.

2 Methodology

HTG-PBR employ a unique fuel recirculation system that enables continuous operation and uniform fuel utilization. Figure 1 depicts a typical pebble circulation path in HTG-PBR with the dynamic motion and forces acting on the pebbles indicated. As demonstrated by the labels in Figure 1, the reactor core is subdivided into multiple regions to enable a grouping strategy for fuel pebbles. This approach offers efficient modeling of pebble circulation, fuel mixing, and fresh fuel insertion, while maintaining mass conservation and preserving the essential physics governing reactor behavior.

Figure 1

Figure 1. Schematic demonstration of the circulation flow of pebbles in HTG-PBR.

In our earlier work (Mehta et al., 2024), a high-fidelity neutronics analysis of HTG-PBR was performed by coupling the Monte Carlo neutron transport method with DEM using OpenMC and LIGGGHTS, respectively. This approach enabled detailed tracking of neutron transport and fuel depletion in moving pebbles. In that work, the HTG-PBR core was discretized by grouping pebbles into both axial and radial regions (see Figure 1) to improve computational efficiency while preserving the fidelity of pebble circulation modeling. The study successfully established an efficient equilibrium core search process and rendered a preliminary neutronically equilibrium steady-state HTG-PBR core. The results of the study provided insights on the required fresh fuel loading rate and the number of pebble-passes necessary to sustain criticality for a typical HTG-PBR over long-term operation (Mehta et al., 2024). Following these efforts, the present study shifts the analysis focus from neutronics to the pebble movements and thermal-hydraulic behavior of HTG-PBR by employing a DEM-CFD coupling method. Given the strong interplay between pebble motion and heat transfer within the HTG-PBR, accurately modeling the fluid-solid interaction becomes essential for predicting local temperature gradients and flow dynamics within the reactor.

In the present study, a simplified reactor model problem is adopted to evaluate the developed method. This model is a scaled-down HTG-PBR model with a reactor-to-pebble diameter ratio of 7.33 and a total of 789 pebbles packed inside a conical vessel (Reger et al., 2024). The pebble dynamics and thermal-fluid analysis were carried out in a two-step procedure. First, the pebble dynamics simulation was performed using DEM in LIGGGHTS to investigate the circulation and intermixing behavior of pebbles across defined axial and radial subregions, This approach provides realistic insight into how pebbles redistribute within the reactor during operation. Second, the thermal-fluid analysis was performed using the CFD model in OpenFOAM. The mathematical models employed in this study, along with the DEM-CFD-neutronics coupling framework, are described in the following subsections.

Although this model problem contains only 789 pebbles, the developed DEM–CFD framework is designed to scale to full-core reactors containing more than 100,000 pebbles. This scalability is achieved by using a porous-media CFD formulation and a DEM grouping/recycling strategy, and is supported by the validation efforts involving ∼30,000 pebbles from Jiang and Reger benchmark (Jiang et al., 2012; Reger et al., 2023) and ∼9,500 pebbles from the SANA experiment (Stocker and Nieben, 1996; Zou et al., 2017). The DEM solver in LIGGGHTS has been verified in prior studies (Jiang et al., 2012; Reger et al., 2023) to efficiently simulate up to 30,000 pebbles with linear scaling of memory and runtime using parallel domain decomposition. Full-core HTG-PBR configurations involving over 100,000 pebbles are feasible on high-performance computing clusters with sufficient node distribution. The CFD solver in OpenFOAM also supports scalability through the use of unresolved porous-media modeling, which eliminates the need to resolve flow around each pebble. As such, the overall methodology is suitable for extension to steady-state full-core analyses with acceptable computational cost.

2.1 DEM–CFD-neutronics coupled approach

To enable consistent DEM-CFD–neutronic coupled analysis, the reactor domain is first divided into multiple subregions, with pebbles grouped accordingly as described earlier. The identical reactor geometry is independently constructed in LIGGGHTS for DEM simulation, OpenFOAM for the CFD analysis, and OpenMC for neutronics analysis, as illustrated in Figure 2. Furthermore, the OpenMC model shown in Figure 2c includes not only the fuel pebbles and core structure but also the surrounding graphite reflector and helium coolant channels indicated with sky blue. In addition to all side views of the reactor, a top-down view of the middle plane of the reactor is included for the OpenMC model.

Figure 2

Figure 2. Geometrically consistent reactor model created in (a) LIGGGHTS, (b) OpenFOAM, and (c) OpenMC.

A multiphysics coupling scheme between DEM, CFD, and neutronics simulation is illustrated in Figure 3. The power distribution is computed from OpenMC and exported as input to OpenFOAM, which serves as the CFD solver. The CFD model was constructed using porosity profiles obtained from DEM simulations, in which helium was modeled as the coolant flowing through the pebble bed. Parametric studies were conducted for different helium inlet velocities to assess their impact on thermal performance. Temperature distributions were computed under LTNE conditions, allowing evaluation of solid-fluid temperature differences and axial temperature gradients within the core.

Figure 3

Figure 3. Diagram of the coupling scheme for CFD-DEM-Neutronics analysis.

In this work, the coupling between DEM and CFD is established by first running the DEM simulation in LIGGGHTS to determine the spatial distribution of pebbles and their intermixing behavior across axial and radial subregions. The resulting local porosity profiles are then extracted and mapped into the CFD model in OpenFOAM. The porosity values used in the CFD model were obtained by averaging the local void fraction in each of the four subregions defined in the DEM simulation. For each region, the porosity was computed as the unoccupied volume divided by the total volume, based on the final packing configuration after steady circulation. These scalar porosity values were then assigned to their respective zones in the CFD mesh using the setFields utility in OpenFOAM. No interpolation or smoothing was applied, as each subregion was treated as a uniformly porous domain for this model problem.

These porosity values influence key thermal-fluid properties, including flow resistance and heat transfer characteristics, in each subregion of the porous-media model. While the DEM and CFD simulations are executed sequentially rather than in a tightly coupled loop, this approach allows the pebble dynamics—particularly mixing patterns and local packing structures—to inform and enhance the realism of the thermal-fluid model. While this sequential DEM–CFD approach provides a computationally efficient pathway for incorporating realistic pebble packing and mixing into thermal-fluid analysis, it does not account for two-way or transient feedback effects. Specifically, any influence of local temperature fields on pebble motion—such as thermally induced structural expansion, frictional property changes, or temperature-dependent packing rearrangement—is not modeled. This assumption is considered acceptable for steady-state performance analysis but may introduce limitations when modeling transient events or feedback-dominated phenomena.

This represents the core focus of the present work: using DEM-informed geometric and transport parameters to improve the fidelity of CFD-based heat transfer simulations in HTG-PBR systems. To incorporate the effect of pebble movement effect into neutronics analysis, the pebble positions and circulation dynamics are modeled using the DEM solver in LIGGGHTS, which tracks the gravitationally driven motion and mixing of pebbles in each subregion. Given the low rate of pebble movement during reactor operation (pebble moving speed ∼15.6 cm/day (Jaradat et al., 2023)), the dominant force considered in the DEM simulation is gravity. The calculated pebble positions and mixing information are then exported to the neutronics model developed with OpenMC using a custom Python script that generates the required pebble location information for OpenMC. At each position, OpenMC places fuel-containing pebbles to conduct neutronic analysis.

For the thermal-fluid simulation, the widely used CFD code OpenFOAM is employed using a homogenized porous-media model. A non-uniform porosity distribution is allowed for each subregion in the model. This local porosity is calculated from the DEM and imported into OpenFOAM. The local porosity is crucial for determining some key model parameters such as the flow resistance based on the fluid-solid interaction model and the thermal transport properties within each region of the reactor.

Regarding neutronics and thermal-hydraulics coupling, since neutronics analysis is not the main focus of this work, only the forward (one-way) coupling from OpenMC to OpenFOAM is implemented, which means the present work does not implement full multiphysics couplings. Only a one-way transfer of power distribution from a static neutronics solution is applied to the CFD model to enable spatially resolved thermal analysis. No feedback from temperature or flow fields to neutronics is included in this study. The full multiphysics coupling framework illustrated in Figure 3 represents a planned capability of future work building upon the present study.

Specifically, a steady-state neutronic simulation in OpenMC assuming an equilibrium core, and then the resulting spatially resolved power distribution is normalized and mapped as a volumetric heat source into OpenFOAM using the setFields utility. OpenFOAM solves the steady-state energy equations to obtain the temperature distribution in the reactor domain. However, Figure 3 above illustrates a full two-way coupling scheme that represents a planned extension of the present work, where temperature-dependent cross-sections will be used in OpenMC and iterative coupling with OpenFOAM will be implemented. As clarified earlier, although the overall simulation framework emphasizes a multiphysics coupling mechanism encompassing neutronics, CFD, and DEM calculations, the present work primarily focuses on the integration of DEM and CFD to evaluate thermal-fluid behavior under realistic pebble mixing effects. The one-way coupling between the neutronics and thermal-hydraulics models is implemented to provide a more physically meaningful and spatially resolved heat source for the thermal-hydraulic analysis. No iterative or time-dependent coupling with neutronics is performed in this study. This clarification is important to correctly position the scope of the present work.

2.2 Mathematical models

The detailed neutronics model employed in this study has been clearly outlined in the previous work (Mehta et al., 2024), and thus is not repeated here. The following parts only include computational models used for the DEM and CFD calculations, respectively.

2.2.1 DEM model

The DEM model employs a Lagrangian approach that involves the explicit tracking and solving of trajectories for all pebbles in the HTG-PBR, which is essentially governed by the linear momentum and angular momentum equations of pebbles (Latifi et al., 2024; Kloss et al., 2012) given by Equations 1, 2, respectively

$m_p\frac{d{\mathbf{U}}_p}{dt}={\mathbf{F}}^c+{\mathbf{F}}^g+{\mathbf{F}}^b+{\mathbf{F}}^d+{\mathbf{F}}^o(1)$

$I\frac{d{\mathbf{\omega }}_p}{dt}={\mathbf{M}}^c(2)$

Where ${\mathbf{F}}^c$ is the contact force acting on a pebble due to interaction with other pebbles or wall, respectively; ${\mathbf{F}}^g$ is the gravitational force; ${\mathbf{F}}^b$ and drag force ${\mathbf{F}}^d$ are the buoyancy and drag forces resulting from interactions with the surrounding fluid, respectively; ${\mathbf{F}}^o$ collectively represents all other additional forces including van der Waals, electrostatic, magnetic, cohesiveness, or any external forces. ${\mathbf{M}}^c$ and ${\mathbf{\omega }}_p$ are the torque and rotational velocity of the pebbles; The dominant forces considered in this study include pebble-pebble interactions, pebble-reactor interactions, gravitational, buoyancy, and drag forces.

2.2.2 CFD model

The flow of helium gas through a porous pebble bed is governed by the continuity and momentum equations of the fluid given in Equations 3, 4, respectively as follows (Jaradat et al., 2023; Nield and Bejan, 2006)

$\frac{\mathit{\partial }\epsilon {\rho }_f}{\mathit{\partial }t}+\nabla ·\left(\epsilon {\rho }_f\mathit{u}\right)=0(3)$

$\frac{\mathit{\partial }\epsilon {\rho }_f\mathbf{u}}{\mathit{\partial }t}+\nabla ·\left(\epsilon {\rho }_f\mathbf{\text{uu}}\right)-\nabla ·\left(\mu \nabla \mathbf{u}\right)=-\nabla \epsilon p-{\rho }_f{\mathbf{F}}_{pf}+\epsilon {\rho }_f\mathbf{g}(4)$

where ${\rho }_f$ and $\mu$ represent the fluid density and viscosity, respectively; $\mathbf{u}$ and p represent the fluid velocity and pressure fields, respectively; ${\mathbf{F}}_{pf}$ represents the momentum exchange between the solid structure (pebbles) and the fluid; $\epsilon$ stands for fluid medium porosity.

Note that the momentum exchange between the solid structure and the gaseous fluid, denoted as ${\mathbf{F}}_{pf}$ in Equation 4, is described by the Ergun equation established in 1952, in which the momentum exchange term ${\mathbf{F}}_{pf}$ is estimated with a drag coefficient $K_{pf}$ multiplied by the relative velocity between the gaseous fluid and the pebble $\left(\mathbf{u}-{\mathbf{U}}_{\mathit{p}}\right)$, namely by Equation 5,

${\mathbf{F}}_{pf}=K_{pf}\left(\mathbf{u}-{\mathbf{U}}_{\mathit{p}}\right)(5)$

where the drag coefficient is determined by the following empirical correlation (Ergun, 1952) given by Equation 6

$\begin{array}{l}{K}_{pf}=150\frac{{\left(1-\epsilon \right)}^2\vartheta }{\epsilon {\left(d_p\right)}^2}+1.75\frac{\left(1-\epsilon \right)\left|\left(\mathbf{u}-{\mathbf{U}}_{\mathit{p}}\right)\right|}{d_p}\end{array}(6)$

Here $\vartheta$ is the kinematic viscosity of the fluid; $d_p$ and ${\mathbf{U}}_{\mathit{p}}$ are the diameter and moving velocity of the pebble, respectively. However, the pebble movement is nearly negligible compared to the velocity of helium gas when the reactor is at normal operation. Therefore, the pebble movement is considered to be zero in all calculations in the present study. It is worth noting that while the Ergun model is widely used for predicting momentum exchange in porous beds, its accuracy can degrade near the reactor walls due to porosity gradients and flow channeling effects. Although the DEM-informed porosity field partially accounts for this spatial variation, the drag formulation itself remains uniform. This may introduce localized deviations in near-wall velocity and heat transfer predictions. However, the close agreement observed with the SANA benchmark—which includes radial heat dissipation to the sidewall—suggests that these limitations are not significant at the global scale for steady-state thermal analysis.

To computationally model the heat energy interactions within the HTG-PBR, the energy conservation equation in fluid and solid (pebbles) phases needs to be solved along with the momentum equations described earlier. The heat energy conservation equations for the fluid and solid phases are given by Equations 7, 8 respectively as follows (Jaradat et al., 2023; Nield and Bejan, 2006)

$\epsilon {\rho }_f{C}_{p,f}\frac{\mathit{\partial }{T}_f}{\mathit{\partial }t}+\epsilon {\rho }_f{C}_{p,f}\nabla \mathbf{u}·\left({\mathrm{T}}_f\right)-\nabla ·\left(\epsilon k_f\nabla {\mathrm{T}}_f\right)+\alpha \left(T_f-T_s\right)={\dot{q}}_f^{‴}(7)$

$\left(1-\epsilon \right){\rho }_s{C}_{p,s}\frac{\mathit{\partial }{T}_s}{\mathit{\partial }t}-\nabla ·\left(\epsilon k_s\nabla {\mathrm{T}}_s\right)+\alpha \left(T_s-T_f\right)={\dot{q}}_s^{‴}(8)$

where $C_{p,f}$ is the constant pressure specific heat, respectively. ${\rho }_s,C_{p,s}$ is the particle or pebble density and specific heat, respectively. $k_f$ and $k_s$ are the effective thermal conductivities of the fluid and solid, respectively. $T_f$ and $T_s$ are the temperature of the fluid and temperature of porous bed of pebbles, respectively.

Equations 7, 8 are coupled by the heat transfer coefficient ($\alpha$), and thus must be solved simultaneously. In solid (pebbles), the heat source (${\dot{q}}_f^{‴}$) counts the heat generated from the nuclear fuel. The coupling of the fluid energy equation (e.g., Equation 7) and the solid energy equation (e.g., Equation 8) through the heat transfer coefficient can be achieved with the Kerntechnischer Ausschuss (KTA) model (Jaradat et al., 2023; Nield and Bejan, 2006). This model provides a volumetric source coefficient through which, the solid (pebbles) releases the energy and fluid (helium) gains this energy due to the temperature difference between fluid and solid. The KTA correlation employs the Nusselt number to derive the heat transfer coefficient.

$\text{Nu}=1.27\frac{{\mathrm{Pr}}^{1/3}{\text{Re}}^{0.36}}{{\epsilon }^{1.18}}+0.033\frac{{\mathrm{Pr}}^{0.5}{\text{Re}}^{0.86}}{{\epsilon }^{1.07}}(9)$

Equation 9 is valid for a bed height greater than four times the pebble diameter and 0.36 < $\epsilon$ <0.42, 100 < Re < ${10}^5$, and Pr = 0.7.

Precise modeling of the energy and heat conduction equations also relies significantly on the accurate estimate of the thermal conductivities ($k_f$ and $k_s$) appearing in Equations 7, 8 because ensuring accuracy in these properties is crucial for assessing the temperature distribution within the reactor. In fluid flow through porous media, the effective thermal conductivity ($k_f^{\prime }$) is used as a function of the Reynolds number (Re) and Prandtl number (Pr) given by Equation 10. The correlation is developed based on the Péclet number (Pe), which is defined as the product of Reynolds and Prandtl numbers (Pe = Re × Pr), representing the ratio of advective to diffusive heat transport

$k_f^{\prime }=\epsilon k_f+C_0\text{Pe}{k}_f(10)$

where $C_0$ is a constant coefficient to be determined with additional conditions.

Likewise, obtaining an effective conductivity in the solid region ($k_s^{\prime }$) to use in Equation 11 is also crucial, yet it is more intricate to achieve. Various effective thermal conductivity models for two-phase flow have been extensively discussed in the literature (Van Antwerpen et al., 2010; Hsu et al., 1994). Among these well-studied models, the Zehner-Bauer-Schlünder (ZBS) model is considered in this work because it incorporates the heat transfer effect contributed from all three heat transfer modes

$k_s^{\prime }=k_{\text{radiation}}+k_{\text{solid}}+k_{\text{fluid}}(11)$

The radiation component of the solid effective conductivity is modeled using the Breitbach and Barthels model (Breitbach and Barthels, 1980) In the present model, these contributions are evaluated through the radiation, contact-conduction, and fluid-conduction components given in Equations 12–14.

$k_{\text{radiation}}=\left[\left(1-\sqrt{1-\epsilon }\right)\epsilon +\frac{\sqrt{1-\epsilon }}{2/e-1}\frac{B+1}{B}\frac{1}{1+\frac{1}{\left(2/e-1\right)\mathrm{\Lambda }}}\right]4\sigma d_p{\stackrel{—}{T}}^3(12)$

where $e$ represents the spectral emissivity of the pebbles, $B$ is the function of porosity, $\mathrm{\Lambda }$ is the geometric factor calculated with particle arrangement, and $\sigma$ is the Stefan–Boltzmann constant The pebble-to-pebble effective heat conduction effect is handled with the Chan and Tien model (Janse van Rensburg and Kleingeld, 2011; Latifi et al., 2024)

$k_{\text{solid}}=\frac{1}{2\times 0.53}\frac{N_A}{N_L}\left(\frac{d_c}{d_p}\right)d_p{k}_p(13)$

where $N_A/N_L$ is the ratio of actual contact points to the ideal contact points per pebbles, $d_c/d_p$ is the ratio of contact diameter to the pebble diameter, and $k_p$ is the thermal conductivity of the pebble. The heat conduction of pebbles through fluid component is handled with the ZBS model (Jaradat et al., 2023; Van Antwerpen et al., 2010; Hsu et al., 1994).

$k_{\text{fluid}}=(1-\left(\sqrt{1-\epsilon }\right)k_f+\sqrt{1-\epsilon }{K}_{SF}{k}_f(14)$

where $K_{SF}$ is a correction factor accounting for the interaction between the solid and fluid phases. The ZBS effective conductivity model accounts for numerous interaction modes between pebbles but surpasses the heat transfer modes observed in fluid flow across pebbles. The accounted modes in ZBS model involve radiation between pebbles, direct conduction at contact points, and pebble-conducted heat through the fluid present between pebbles. Therefore, the representation of solid effective conductivity encompasses these varied interaction modes. It is worth noting that the International Atomic Energy Agency recommends the use of the ZBS effective thermal conductivity model in pebble bed reactors (Romano et al., 2015).

The use of empirical correlations such as Ergun (for momentum exchange), KTA (for interphase heat transfer), and ZBS (for solid-phase conductivity) correlations inevitably introduces uncertainties, particularly when they are applied to spatially varying porosity fields from DEM. While these correlations are widely validated for packed-bed flows, their predictive accuracy may degrade at porosity extremes or near-wall gradients. Although a full uncertainty propagation analysis is beyond the scope of this study, care was taken to apply the correlations within their recommended ranges for Reynolds number and porosity. Additionally, turbulence effects are not modeled in the current CFD implementation, as the helium flow Reynolds numbers remain in the laminar-to-transitional range (typically Re < 1,000) for the operating conditions considered. This assumption simplifies the analysis while remaining consistent with prior reactor-scale porous media simulations. Moreover, the strong agreement with the SANA benchmark (presented in the next section) provides confidence in their applicability for reactor-scale thermal predictions under steady-state conditions.

3 Model validations

Before performing the pebble dynamics and thermal-fluid analysis of the HTG-PBR reactor, a couple of verification and validation efforts were carried out to benchmark the developed DEM and thermal models against established computational and experimental data.

3.1 Validation of the DEM model for pebble dynamics

In the first validation work, the DEM model for pebble dynamics was evaluated by comparing simulation results with experimental data measured by Jiang et al. (2012) and computational results provided by Reger et al. (2023). Jiang and his co-workers conducted a series of experiments to examine pebble flow behavior in cylindrical vessels of varying sizes, providing qualitative insights into flow patterns and packing configurations. Building on Jiang’s work, Reger and his co-workers utilized the same experimental dataset to extract radial velocity profiles by tracking colored pebbles, thereby enabling a quantitative comparison with their own DEM-based simulations.

In this study, the radial velocity profiles reported by Jiang et al. (2012) and Reger et al. (2023) were used to verify our DEM model. Comparisons are made against both the original experimental measurements and Reger’s DEM simulation results to confirm the consistency and accuracy of our simulation. Table 1 summarizes the key parameters of the specific cylindrical vessel and pebble used in the validation work.

Table 1

Table 1. Parameters of the experimental vessel and pebble.

Figure 4 compares the vertical velocity distribution of pebbles obtained from Reger’s experimental study and the present simulation. The velocity magnitudes are represented by colorful contours in the units of m/s, in which the blue color indicates downward motion and the red color indicates upward. The right plot shows axial locations 8d, 15d, 30d, and 70d from the pebble bed entrance, used to extract radial profiles as plotted in Figure 4.

Figure 4

Figure 4. Average velocity profile of pebbles from the present (left) and Reger’s study (right).

In the experimental setup, pebbles are discharged at a rate of 150 pebbles per minute, whereas in the simulation, the corresponding velocity profiles are generated for a discharge rate of 15,000 pebbles per minute to accelerate steady-state flow development. This 100× faster discharge rate was intentionally used as a numerical acceleration strategy to reach the steady pebble-flow pattern more quickly. The underlying assumption of this setup is that the fully developed velocity field is independent of the transient circulation speed, and the simulation results confirm this - the simulated velocity profiles (see Figure 5) match both Jiang’s experimental data and Reger’s DEM results despite the accelerated discharge rate. Importantly, this acceleration was a numerical strategy and does not imply a physically realistic circulation rate, but rather, it assumes that the quasi-steady flow field (the pebble velocity profile) is the same once it is fully developed, regardless of how fast we drive the circulation initially. Additional trial simulations at intermediate discharge rates (e.g., 50× and 75× of the experimental value) were conducted to confirm the robustness of this assumption. The results demonstrated that the resulting steady-state velocity profiles and mixing metrics remained effectively unchanged across this range, reinforcing the validity of using an accelerated circulation rate without compromising physical fidelity. The validity of this assumption is further supported by results discussed below.

Figure 5

Figure 5. Comparison of the velocity profile with Reger’s computational result and Jiang’s experimental data at different levels: (a) 8d, (b) 15d, (c) 30d, and (d) 70d.

Figure 5 displays the velocity contour plots from the simulation and compares them with the radial velocity profiles reported by Reger et al. (2023). In addition, Ref (Reger et al., 2023). presents velocity profiles at four axial positions along the reactor height, ranging from 8 times the pebble diameter (8d) above the cone to 70d, where 0d refers to the base of the conical section. The x-axis represents the horizontal position (in terms of pebble diameters), while the y-axis indicates the vertical velocity in pebble diameters per second. The profiles are plotted for computational results from the present study (solid line), numerical data from Reger et al. (2023) (blue circles), and experimental data from Jiang et al. (2012) (orange crosses). The present study shows excellent agreement with both experimental and previous computational results, capturing the characteristic “V-shaped” velocity distribution. These comparisons verify the accuracy and reliability of the present computational approach for modeling vertical velocity in the pebble bed reactor.

The simulated velocity profiles at 70 and 30d exhibit good agreement with both the experimental data and Reger’s DEM results. However, in the lower region near the cone (8d and below), the simulation slightly overpredicts the velocity magnitude compared to experimental observations, although it remains consistent with the velocity trends observed in Reger’s DEM simulation.

To further quantify the level of agreement, the mean absolute percent error (MAPE) were computed between the DEM-simulated and benchmark radial velocity profiles at each axial level. The errors remained below 5.0% across all locations, confirming the high accuracy of the DEM model. These results are summarized in Table 2.

Table 2

Table 2. Quantitative error metrics for DEM model validation.

3.2 Validation of the CFD model for thermal analysis

The second validation was performed using the SANA experimental benchmark (Stocker and Nieben, 1996; Zou et al., 2017) to ensure the accuracy of the developed thermal model. The SANA experiment provides well-characterized data on heat transfer in a pebble bed under steady-state conditions, allowing a direct comparison between simulated and measured temperature behavior. Using this experiment as a validation case helps verify that the model can predict real-world thermal responses in a pebble bed reactor, thereby building confidence in its use for reactor analysis and design.

The SANA experiment was designed to investigate heat transfer within a cylindrical pebble bed, with a focus on temperature profiles and heat dissipation mechanisms (conduction, convection, and radiation). The test facility consists of a stainless-steel cylinder (approximately 1.0 m tall and 1.5 m in diameter) filled with about 9,500 spherical pebbles. A centrally located electrical heater (0.13 m in diameter) provides a controlled heat source (around 10 kW) to simulate reactor-core conditions. Helium gas fills the void spaces between pebbles, serving as the coolant. Heat transfer in the experiment occurs via solid conduction through pebbles, gas-phase convection, and inter-surface thermal radiation within the bed. The boundary conditions of the setup are configured to mimic a high-temperature pebble-bed reactor: the top and bottom surfaces of the cylinder are insulated (adiabatic), while the cylindrical side wall is exposed, allowing natural convective cooling to the surrounding helium and radiative heat loss to the environment. To perform natural circulation analysis scenario, the gravitational buoyancy term in the CFD momentum equations is activated and imposed no forced inlet flow, allowing the helium coolant to circulate freely due to density differences. A schematic view of the SANA experiment is shown in Figure 6.

Figure 6

Figure 6. Schematic view of the SANA experiment.

The simulation results of the SANA experiment were evaluated using one-quarter of the reactor geometry by taking advantage of its geometric symmetry. Figure 7a gives a cut-away view of the three-dimensional steady-state temperature distribution of the pebble bed. The temperature is highest in regions close to the heat source and gradually decreases toward the outer cylindrical wall due to radial heat transfer. Figure 7b compares the radial temperature profile obtained from the present study to the recorded SANA experimental data and the predictions from the Pronghorn pebble-bed reactor code (Zou et al., 2017). The radial coordinate spans from the center (pebble bed core) to the outer wall. The model’s predicted temperatures (solid line) show excellent agreement with the SANA experiment data (discrete points) across the radius. Both in the core region and through the mid-radius of the bed, the simulation closely reproduces the experimental temperatures, indicating that the key heat transfer mechanisms are accurately captured. In contrast, the Pronghorn model (dashed line) underestimates the temperature in the mid-to-outer radial regions (Zou et al., 2017), deviating below the experimental curve. This discrepancy highlights the improved accuracy of the present model in representing radial heat dissipation in the pebble bed. The strong alignment between the current simulation and the SANA measurements validates the model’s ability to predict the steady-state thermal profile in a pebble bed reactor.

Figure 7

Figure 7. (a) Cut-away view of the temperature distribution in the pebble bed; (b) Comparison of temperature profiles obtained from the present study, the Pronghorn model, and the SANA experiment.

Overall, the validation results demonstrate that the computational model can reliably reproduce the thermal behavior observed in the SANA experiment. The close agreement between simulated and experimental temperature profiles confirms that the model accurately captures key heat transfer mechanisms in the pebble bed, including pebble-pebble conduction and multi-directional radiation. Deviations from the experimental data remain within acceptable error margins, supporting the model’s credibility.

With the SANA benchmark validation completed, the model is now applied to a representative high-temperature reactor scenario (e.g., an Xe-100 core). This includes evaluating temperature distributions under varying helium inlet velocities and power densities using DEM-derived porosity profiles and LTNE modeling, demonstrating the solver’s capability to simulate realistic high-temperature gas-cooled reactor conditions.

Although the axial temperature distribution may appear visually uniform in the contour plots, simulations reveal a 20–30 K temperature difference from bottom to top. This gradient is less apparent due to the color scale resolution. The thermal boundary conditions include adiabatic top and bottom walls and convective heat loss through the lateral surfaces, which collectively produce this moderate axial variation under steady-state conditions.

For the CFD model validation, the root mean square error (RMSE) and MAPE values were calculated between the simulated temperature profile and SANA experimental data shown in Figure 7b. The RMSE was 6.2 K and the MAPE was 2.7% across the radial domain, indicating strong agreement with measured thermal behavior.

4 Results and discussion

The modeling approaches outlined in Section 2 are employed to analyze pebble intermixing and thermal-fluid performance for a HTG-PBR-like model problem described in Reger et al. (2024). The DEM and CFD solvers are coupled by first performing a DEM simulation to predict the realistic pebble packing structure and intermixing across subregions of the model problem. The resulting spatial porosity distribution from the DEM model is then mapped into the CFD solver, where it is used to define flow resistance, heat source regions, and thermal transport properties within the porous media. This coupling scheme ensures that the local variations in pebble arrangement and void fraction—driven by realistic gravity-induced motion—are directly reflected in the fluid flow and temperature fields predicted by the CFD model. The DEM simulation of pebble recirculation in a conical-shaped core reveals the progressive intermixing process among the defined subregions, while the CFD simulation reaches a steady-state temperature field, integrating the impact of pebble dynamics on the thermal-fluid behavior of the system.

4.1 Model reactor problem

The model reactor problem, despite not representing an actual HTG-PBR, is widely used to investigate fundamental pebble mixing and flow behavior in pebble beds. The pebble mixing effect during pebble circulation is analyzed within a conical reactor geometry that maintains a reactor-to-pebble diameter ratio (D/d) of 7.33. The model reactor has a total height (H) of 0.6 m and is filled with spherical pebbles, each of which has a diameter (d) of 0.0381 m. The reactor is initially filled with 789 pebbles. Figure 8 shows a schematic illustration of the model problem configuration with an average porosity of 0.417 for the packed bed. The average porosity calculated from the DEM simulation (0.417) is consistent with reported values for random packing of spherical pebbles (typically in a range of 0.36–0.42). This agreement with established packing correlations reinforces the physical realism of the porosity fields used in the CFD model and supports the accuracy of the subsequent thermal-fluid predictions.

Figure 8

Figure 8. Schematic view of the model reactor problem used in this study.

The conical reactor vessel and the pebbles are explicitly modelled in LIGGGHTS. To facilitate the message coupling between DEM and CFD, the reactor domain is divided into four distinct subregions to monitor the redistribution of pebbles throughout reactor operation. In the course of simulation, the DEM and CFD solvers are executed individually. In the DEM simulation, pebbles are initially loaded into the reactor by defining solid wall boundaries and closing the outlet. Once the pebbles fill the center region, the outlet is opened, and a recycling boundary condition is applied—ensuring that each pebble exiting the bottom is reintroduced at the top—until all pebbles have completed about two or three circulation cycles. After this process, the degree of pebble mixing within each subregion is evaluated. In the CFD simulation, a no-slip boundary condition is applied at the reactor walls for fluid flow, and helium gas enters the domain through the inlet. This establishes a forced-circulation flow, unlike the buoyancy-driven natural-circulation regime used in the SANA validation case. For the thermal analysis, the inlet helium temperature is set to 330 K. The outlet temperature remains nearly unchanged (i.e., it has a negligible temperature gradient). At the reactor walls, convective heat transfer to the environment is modeled. For the fuel pebbles, both convective and radiative heat transfer mechanisms are considered to account for energy exchange with the surrounding medium.

4.2 Pebbles intermixing simulation

During the pebble mixing simulation, the reactor domain is subdivided into multiple regions by grouping pebbles spatially. Figure 9a illustrates the result of pebble groupings within the model reactor, where each group contains an equal number of pebbles. Figure 9b provides a cropped view of the grouping to show the inner subregions. After the reactor is fully filled with pebbles, pebble recirculation is initiated under the assumption of periodic boundary conditions in the axial direction. The number of pebbles transitioning between each subregion is recorded throughout the simulation, enabling a detailed assessment of mixing behavior inside the reactor.

Figure 9

Figure 9. The full view (a) and cropped view (b) of the pebble groupings in the reactor.

Figure 10 shows the percentage of total pebbles circulating inside the reactor along with the dimensionless time variable - t/T, where T is total time required for a pebble passing through the core and t is the current time the pebble is circulating on. The choice of time step affects the granularity of pebble displacement updates. A larger time step captures more displacement per interval, resulting in higher measured circulation percentages over the same normalized time. In this study, two different time steps (Δt = 0.02 and Δt = 0.04) are considered during the DEM simulation.

Figure 10

Figure 10. Pebble circulation percentage over normalized time for two DEM time steps (Δt = 0.02 and 0.04), showing consistent steady-state flow behavior despite differences in sampling resolutions.

As shown in Figure 10, each data point represents the cumulative fraction of pebbles that have circulated through the reactor domain at a given time step. For the case of Δt = 0.02, the average pebble circulation percentage is approximately 2.5%, while for Δt = 0.04, this number increases to nearly 5%. This shows that doubling the time step results in approximately double the number of pebbles circulating, indicating that the pebbles are circulating at a nearly constant velocity under gravitational forces within the simulation. Although the circulation percentage varies with time step due to the discrete nature of sampling intervals, the underlying velocity field and dominant intermixing pathways remain consistent. Additional simulations using intermediate time steps (e.g., Δt = 0.03) confirmed that the relative pebble exchange between subregions is stable and not sensitive to time-step variation. Thus, while the resolution of displacement tracking changes, the physical mixing behavior and directional migration trends are robust, supporting the reliability of the DEM-derived intermixing metrics.

Given that the reactor is divided into four subregions, each containing 25% of the total pebbles, the 2.5% circulation implies that roughly 10% of the pebbles in a given subregion have transitioned into its neighboring regions. Notably, approximately two complete subregions undergo a full exchange at around t/T = 0.4. This observation is consistent with the expected gravity-driven recirculation behavior of fuel pebbles in HTG-PBR.

Further insights can be gained by examining the pebble intermixing between subregions during circulation. Figure 11 shows the interchanges of pebbles between the four subregions (SR-P1 through SR-P4) at the level of approximately 2.5% total pebble circulation evaluated with Δt = 0.02. Each subplot (a–d) corresponds to one source subregion, with coloured lines showing the percentage of pebbles that have moved into the other three regions. Figure 11a indicates that pebbles from SR-P1 predominantly mix into SR-P3, due to the central upward position of SR-P1. There is also noticeable mixing with SR-P2, while mixing with SR-P4 is negligible. Figure 11b indicates that pebbles from SR-P2 primarily migrate into SR-P1 (∼15%), with minimal interaction with other regions. Figure 11c indicates that pebbles from SR-P3 mix significantly into SR-P4 (∼12%) and SR-P1 (∼7%). The reason for this redistribution is that as recirculated (or fresh) pebbles are inserted into SR-P3, the number of pebbles in this subregion increases. To maintain a constant number of pebbles in each subregion, some of the SR-P3 pebbles are redistributed into SR-P4. Without this balancing mechanism, the dominant mixing from SR-P3 would likely remain toward SR-P1, with only minor mixing toward SR-P2. Lastly, Figure 11d indicates intermixing of SR-P4 into SR-P3 and SR-P2, both around ∼7%, while migration to SR-P1 is minimal.

Figure 11

Figure 11. Pebble intermixing between subregions at Δt = 0.02. Subplots show pebble migration originating from (a) SR-P1, (b) SR-P2, (c) SR-P3, and (d) SR-P4. These plots highlight dominant mixing from SR-P3 to SR-P1 and SR-P4, revealing asymmetric and non-random migration patterns.

These pebble intermixing results generated by DEM simulations capture realistic pebble movement and mixing dynamics, which can be subsequently integrated into the neutronic simulation through subregion updates. For a full-core analysis, the HTG-PBR model would need to be divided into a greater number of axial and radial zones to account for the more frequent mixing that occurs near the conical outlet and the uppermost sections of the core. Integrating DEM-informed intermixing into neutronic modeling offers a more accurate and scalable strategy for the whole-core HTG-PBR analysis. Unlike typical random reshuffling assumptions, where each pebble is reassigned uniformly across regions after each pass, the DEM-informed intermixing captures asymmetric and localized transitions between subregions, such as the dominant flow from SR-P3 to SR-P1 and SR-P4 seen in Figure 11c. These physically driven mixing patterns result in a more accurate representation of fuel residence times and migration paths, which are essential for realistic power and burnup modeling in full-core simulations. Although the present thermal analysis uses a uniform power source for the model problem, in a fully coupled setting, the intermixing trends predicted by DEM would affect subregion-wise fuel composition and thus power distribution. For example, the asymmetric migration observed from SR-P3 to SR-P1 and SR-P4 (see Figure 11) would result in locally higher heat generation in those zones over time. This could alter local thermal gradients and impact peak temperature predictions. Future coupled simulations will incorporate this feedback to assess the sensitivity of thermal profiles to realistic mixing behavior.

Although no direct experimental data exist to validate the subregion-level pebble intermixing behavior, the results are physically consistent with gravity-driven granular flow in conical vessels. Given that the underlying velocity field has been benchmarked against experimental and numerical references, we consider the predicted intermixing behavior to be a reasonable and credible reflection of actual pebble motion under HTG-PBR-like conditions.

4.3 Thermal-fluid simulation

The thermal-fluid analysis of the model problem is performed by examining the temperature field of components within the reactor core. A constant volumetric heat source of 2 W/m³ is assumed in this analysis, along with a constant helium inlet velocity (2 m/s) at the reactor inlet. The heat source of 2 W/m³ used in the model problem is significantly lower than typical HTG-PBR volumetric power densities, which generally range from 1 to 3 MW/m³. This lower value was chosen to ensure numerical stability and to demonstrate the thermal-fluid solver’s behavior under simplified but representative flow and heat transfer conditions. Although absolute temperature magnitudes would increase under realistic power loads, the qualitative trends in coolant flow, thermal gradients, and solid-fluid temperature differences would remain consistent. The coupling framework and modeling strategy are not sensitive to this scaling, and future work will incorporate reactor-relevant heat loads alongside neutronic feedback.

This inlet boundary condition represents a pumped, forced-circulation flow, in contrast to the purely buoyancy-driven natural-circulation regime used in the SANA validation. Although the LTNE model was validated under natural circulation conditions (SANA), the underlying heat transfer parameters—such as the KTA correlation for interphase exchange—are based on local Reynolds and Prandtl numbers, allowing them to adapt to both buoyancy-driven and forced-convection regimes. For the reactor simulation, the flow remains within the acceptable Reynolds range for these correlations, ensuring their continued applicability in the forced-flow context. The total reactor height is 0.6 m, and the pebbles are packed up to a height of 0.4 m in the core.

Figure 12 shows the calculated temperature distribution of the helium coolant inside the core. As shown in the figure, the highest temperatures are observed in the bottom region of the reactor core, where the heat generation and accumulation is most significant. The coolant temperature gradually decreases both radially and axially as the heat dissipates through the fluid and reactor walls due to convection, resulting in cooler regions near the boundaries and the upper sections. The visualization highlights the thermal performance and heat transfer characteristics of the reactor, which are critical for the reactor operational safety assessment.

Figure 12

Figure 12. Steady-state helium temperature distribution in the model reactor.

To further examine the thermal-fluid characteristics of the components inside the reactor, Figure 13 displays the radial and axial temperature profiles for both the helium and pebble phases in the core region. The radial temperature is evaluated at the mid-plane of packed bed in the reactor core, while the axial temperature is measured along the central axis of the core.

Figure 13

Figure 13. Radial (a) and axial (b) temperature profiles of helium and pebble phases in the reactor.

The radial temperature profiles indicate that pebble temperatures are consistently higher than helium temperatures, due to the internal heat generation in pebbles and the imperfect heat exchange with the coolant. A near-parabolic temperature profile is observed, with the peak in the central region and a gradual decrease toward the walls, driven by radial conductive and convective transport. The axial temperature profiles indicate both helium and pebble temperatures rise with height due to the accumulation of heat generated in the lower sections. The helium temperature begins to level off in the upper region, while the pebble temperature shows a sharper decline near the outlet, where convective cooling becomes dominant. These profiles confirm that the LTNE assumption effectively captures the separate thermal responses of solid and fluid phases. The insights gained from the thermal-fluid analysis will contribute to optimizing reactor performance by providing a more detailed understanding of thermal transport mechanisms within the pebble bed of a full-scale HTG-PBR analysis.

5 Conclusions and future work

This study investigates the pebble dynamics and thermal-fluid behavior in HTG-PBR through a development of multiphysics-coupled computation simulation framework. A representative conical HTG-PBR core geometry was used in the DEM simulations, which accurately captured pebble circulation and intermixing phenomena across prescribed axial and radial subregions. Thermal-fluid simulations were carried out using a CFD approach in OpenFOAM with an unresolved porous-media formulation and an LTNE model. The accuracy of the DEM and CFD models was demonstrated through separate benchmark studies that validate their suitability for HTG-PBR analysis. Which are applied in this work to simulate a model reactor problem that mimics a realistic HTG-PBR environment.

The pebble intermixing results yielded from the model reactor problem reveal the dominant migration pathways that align with gravitational flow dynamics and the reactor’s geometrical configuration. In parallel, the thermal-fluid simulations of the model reactor indicate that the fuel pebbles maintain higher temperatures in the bottom region of the reactor and consistently remain hotter than the helium coolant. This reflects the distinct thermal responses of the solid and fluid phases, confirming that the model captures the key heat transfer characteristics of a generic HTG-PBR system.

One key novelty of the present work lies in its DEM-based pebble intermixing simulations, which introduce a more physically representative mechanism for modeling pebble migration during reactor operation. This allows for realistic fuel reshuffling without the computational burden of resolving individual pebble trajectories in full-core simulations. The approach suggests a possible way for future full-core thermal-fluid analyses that incorporate both spatially and temporally varying heat sources obtained from neutronic calculations. Moreover, the neutronic analysis can be enhanced to include the impacts of pebble intermixing and circulation on fuel burnup, enabling a more realistic representation of the reactor core in fuel cycle analysis. These improvements are planned for future work, building on the results presented here.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

KM: Data curation, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft. BG: Funding acquisition, Resources, Supervision, Writing – review and editing. ZW: Conceptualization, Funding acquisition, Resources, Supervision, Validation, Writing – review and editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This work is performed with the support of the U.S. Department of Energy’s Nuclear Energy University Program (NEUP) with the Award No. DE-NE0009304.

Acknowledgements

The authors acknowledge the College of Engineering at Virginia Commonwealth University for providing high-performance computing resources that contributed to the research results reported in this paper.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ahmed, F., Ara, N., Deshpande, V., Mollah, A. S., and Bhowmik, P. K. (2021). CFD validation with optimized mesh using benchmarking data of pebble-bed high-temperature reactor. Prog. Nucl. Energy 134, 103653. doi:10.1016/j.pnucene.2021.103653

CrossRef Full Text | Google Scholar

Breitbach, G., and Barthels, H. (1980). The radiant heat transfer in the high-temperature reactor core after failure of the afterheat removal systems. Nucl. Technol. 49, 392–399. doi:10.13182/nt80-a17687

CrossRef Full Text | Google Scholar

Ergun, S. (1952). Fluid flow through packed columns. Chem. Eng. Prog. 48, 89–94.

Google Scholar

Fratoni, M., and Greenspan, E. (2010). Equilibrium core composition search methodologies for pebble bed reactors. Nucl. Sci. Eng. 166, 1–16. doi:10.13182/nse09-66

CrossRef Full Text | Google Scholar

Ge, J., Wang, C., Xiao, Y., Tian, W., Qiu, S., Su, G. H., et al. (2016). Thermal-hydraulic analysis of a fluoride-salt-cooled pebble-bed reactor with CFD methodology. Prog. Nucl. Energy 91, 83–96. doi:10.1016/j.pnucene.2016.01.011

CrossRef Full Text | Google Scholar

Hsu, C. T., Cheng, P., and Wong, K. W. (1994). Modified zehner–schlünder models for stagnant thermal conductivity of porous media. Int. J. Heat Mass Transf. 37, 2751–2759. doi:10.1016/0017-9310(94)90392-1

CrossRef Full Text | Google Scholar

Janse van Rensburg, J. J., and Kleingeld, M. (2011). CFD applications in the pebble bed modular reactor project: a decade of progress. Nucl. Eng. Des. 241, 3683–3696. doi:10.1016/j.nucengdes.2011.05.018

CrossRef Full Text | Google Scholar

Jaradat, M., Schunert, S., and Ortensi, J. (2023). Gas-cooled high-temperature pebble-bed reactor reference plant model. Idaho Falls, ID: Idaho National Laboratory. INL/RPT-23-72192.

Google Scholar

Jiang, S., Yang, X., Tang, Z., Wang, W., Tu, J., Liu, Z., et al. (2012). Experimental and numerical validation of a two-region-designed pebble bed reactor with dynamic core. Nucl. Eng. Des. 246, 1–9. doi:10.1016/j.nucengdes.2012.02.005

CrossRef Full Text | Google Scholar

Kloss, C., Goniva, C., Hager, A., Amberger, S., and Pirker, S. (2012). Models, algorithms and validation for opensource DEM and CFD-DEM. Prog. Comput. Fluid Dyn. An Int. J. 12, 140–152. doi:10.1504/pcfd.2012.047457

CrossRef Full Text | Google Scholar

Latifi, M. S., du Toit, C. G., and Alonso, G. (2024). A comparative numerical investigation of the thermal-fluid phenomena between helium and nitrogen cooled pebble bed reactors. Powder Technol. 433, 119234. doi:10.1016/j.powtec.2023.119234

CrossRef Full Text | Google Scholar

Lee, J. J., Park, G. C., Kim, K. Y., and Lee, W. J. (2007). Numerical treatment of pebble contact in the flow and heat transfer analysis of a pebble bed reactor core. Nucl. Eng. Des. 237, 2183–2196. doi:10.1016/j.nucengdes.2007.03.046

CrossRef Full Text | Google Scholar

Li, R., Chen, J., Zhu, A., Liang, J., She, D., and Zhang, H. (2025). High-fidelity neutronics/thermal hydraulics/pebble flow coupling simulation of pebble bed reactor HTR-PM. Nucl. Sci. Eng. 201, 1–7. doi:10.1080/00295639.2025.2471712

CrossRef Full Text | Google Scholar

Lu, C., and Wu, Z. (2022). A new method to efficiently estimate the equilibrium state of pebble bed reactors. Nucl. Technol. 208, 1577–1590. doi:10.1080/00295450.2022.2049966

CrossRef Full Text | Google Scholar

Mehta, K. S., Goddard, B., and Wu, Z. (2024). Neutronics analysis on high-temperature gas-cooled pebble bed reactors by coupling monte carlo method and discrete element method. Energies 17 (20), 5188. doi:10.3390/en17205188

CrossRef Full Text | Google Scholar

Mulder, E. J., and Boyes, W. A. (2020). Neutronics characteristics of a 165 MW_th Xe-100 reactor. Nucl. Eng. Des. 357, 110415. doi:10.1016/j.nucengdes.2019.110415

CrossRef Full Text | Google Scholar

Nield, D. A., and Bejan, A. (2006). Convection in porous media. New York: Springer.

Google Scholar

Reger, D., Merzari, E., Balestra, P., Stewart, R., and Strydom, G. (2023). Discrete element simulation of pebble bed reactors on graphics processing units. Ann. Nucl. Energy 190, 109896. doi:10.1016/j.anucene.2023.109896

CrossRef Full Text | Google Scholar

Reger, D., Merzari, E., Balestra, P., Schunert, S., Hassan, Y., and King, S. (2024). Direct numerical simulation and large eddy simulation of a 67-pebble-bed experiment. Nucl. Technol. 210, 1258–1278. doi:10.1080/00295450.2023.2218245

CrossRef Full Text | Google Scholar

Robert, Y., Siaraferas, T., and Fratoni, M. (2023). Hyper-fidelity depletion with discrete motion for pebble bed reactors. Sci. Rep. 13, 12711. doi:10.1038/s41598-023-39186-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Romano, P. K., Horelik, N. E., Herman, B. R., Nelson, A. G., Forget, B., and Smith, K. (2015). OpenMC: a state-of-the-art monte carlo code for research and development. Ann. Nucl. Energy 82, 90–97. doi:10.1016/j.anucene.2014.07.048

CrossRef Full Text | Google Scholar

Sharma, A., Saha, S. K., Sharma, A., Sharma, D., and Chaudhuri, P. (2020). Thermal-hydraulic characteristics of purge gas in a rectangular packed pebble bed of a fusion reactor using DEM-CFD and porous medium analyses. Fusion Eng. Des. 160, 111848. doi:10.1016/j.fusengdes.2020.111848

CrossRef Full Text | Google Scholar

Stewart, R., Balestra, P., Reger, D., Merzari, E., and Strydom, G. (2024). High-fidelity simulations of the run-in process for a pebble-bed reactor. Ann. Nucl. Energy 195, 110193. doi:10.1016/j.anucene.2023.110193

CrossRef Full Text | Google Scholar

Stocker, B., and Nieben, H. Data sets of the SANA experiment 1994–1996. Technical Report, Forschungszentrum Jülich, Germany (1996).

Google Scholar

Tano, M. E., and Ragusa, J. C. (2021). Coupled computational fluid dynamics–discrete element method study of bypass flows in a pebble bed reactor. Nucl. Technol. 207, 1599–1614. doi:10.1080/00295450.2020.1820830

CrossRef Full Text | Google Scholar

Van Antwerpen, W., Du Toit, C. G., and Rousseau, P. G. (2010). A review of correlations to model the packing structure and effective thermal conductivity in packed beds of mono-sized spherical particles. Nucl. Eng. Des. 240, 1803–1818. doi:10.1016/j.nucengdes.2010.03.009

CrossRef Full Text | Google Scholar

Weller, H. G., Tabor, G., Jasak, H., and Fureby, C. (1998). A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput. Phys. 12, 620–631. doi:10.1063/1.168744

CrossRef Full Text | Google Scholar

Wu, S.-C., Sheu, R.-J., Peir, J.-J., and Liang, J.-H. (2013). “Burnup computation for HTR-10 using layer-to-layer movement,” in Proceedings of the international conference on nuclear engineering (ICONE), 55836. New York, NY, USA: American Society of Mechanical Engineers. doi:10.1115/icone21-16485

CrossRef Full Text | Google Scholar

Wu, Z., Wu, Y., Tang, S., Liu, D., Qiu, S., Su, G. H., et al. (2018). DEM-CFD simulation of helium flow characteristics in randomly packed bed for fusion reactors. Prog. Nucl. Energy 109, 29–37. doi:10.1016/j.pnucene.2018.07.010

CrossRef Full Text | Google Scholar

Wu, Z., Wu, Y., Wang, C., Tang, S., Liu, D., Qiu, S., et al. (2019). Experimental and numerical study on helium flow characteristics in randomly packed pebble bed. Ann. Nucl. Energy 128, 268–277. doi:10.1016/j.anucene.2019.01.016

CrossRef Full Text | Google Scholar

Zou, L., Novak, A. J., Martineau, R. C., and Gougar, H. D. (2017). Validation of pronghorn with the SANA experiments. Idaho Falls, ID, USA: Idaho National Laboratory Report. INL/EXT-17-44085.

Google Scholar

Zou, L., Hu, G., O'Grady, D., and Hu, R. (2022). Explicit modeling of pebble temperature in the porous-media model for pebble-bed reactors. Prog. Nucl. Energy 146, 104175. doi:10.1016/j.pnucene.2022.104175

CrossRef Full Text | Google Scholar