Application and validation of a higher-fidelity architecture for coupled neutronic and thermal-hydraulic analysis in the XE-prime nuclear thermal propulsion system

1 Introduction

The primary challenges of longer-term space missions are the stringent requirements imposed on the transport vehicle’s propulsion and life support systems. Solutions began to be investigated during the peak years of the space race when nuclear technology was rapidly evolving. It is from this time in history that nuclear thermal propulsion (NTP) first began to show its promise.

The rapid development of human space exploration capabilities by private companies in the United States and abroad in recent years has renewed interest in NTP systems, especially for crewed Mars missions. The power density of nuclear materials is sufficient to greatly reduce the amount of time required for interplanetary transit in such missions, thereby lowering the crew’s exposure to bone-embrittling microgravity and high-energy ionizing radiation. Conventional propulsion systems require the transport of both fuel and oxidizer; nuclear systems require only the transport of a low-mass working fluid such as hydrogen.

Testing of NTP designs is challenging due to environmental concerns and the consequent regulatory constraints. Many of these challenges can be addressed with computational tools. These tools are several orders of magnitude more powerful than those in existence during NTP’s infancy and are already widely used by today’s nuclear scientists and engineers. Parametric design studies do not require separate prototypes when designing a simulation; rather, only a few variables must be changed. However, these tools still have their limits. Even the smallest models require hours or days to solve the discretized momentum and energy transport equations using supercomputers and the underlying assumptions must be validated with experimental results. An enabling factor in the development of early test articles was relaxed regulations; much of the as-built experimental data comes from this time.

NTP development began in the United States in 1953 with ROVER, a program that developed thirteen total reactors, with each iteration improving the performance of the last. Several of the material selection and flow instability problems in these reactors were solved throughout the duration of the program (1953–1964) that informed the design of those systems (a total of seven) developed under the NERVA program that followed. NERVA intended to develop and demonstrate NTP systems that maintained their structural integrity, could be restarted, had predictable behavior, and could be reliable enough for crewed space operations (Walton, 1991).

The XE-Prime NTP engine was a joint development effort primarily between NASA and the Westinghouse Astronuclear Laboratory (WANL) under the NERVA program. Comprehensive reports from this time (Westinghouse Astronuclear Nuclear Laboratory WANL and Reese, 1967) can and should serve as the basis for validating modern codes and the analysis presented herein. This engine was chosen as it was the most flight-ready design of its time.

Analyzing the thermal-hydraulic behavior of the XE-Prime system can be approached using two general philosophies. High-fidelity models optimize the accuracy of the system’s physics at a substantial cost in terms of the computational power required and the time to complete single-state analyses. For all but the simplest thermal-hydraulic models, the system of interest must be reduced to the furthest possible extent and analyzed in a piecewise manner. Simulations are performed and tuned to accurately match experimental data, if available. This validation step is crucial. If experimental data is unavailable, code-to-code model verification can be used as proxy during initial design phases.

In contrast to high-fidelity models, the reduced-order approach seeks to optimize accuracy and simulation time. Neither accuracy nor simulation time shall be sacrificed to the extent that results are unusable, or the use of the reduced order model provides no benefit over higher fidelity models. The primary use case of reduced order models is to characterize the “bulk” behavior of large systems that cannot be feasibly represented using CFD because of extraordinary mesh and computational requirements. Parametric design can be applied at all scales of a system with reduced order codes due to their efficiency for this reason. Analogous to the high-fidelity models, the accuracy, or even the quantification of inaccuracy, must be understood for the reduced-order model prior to its usage in design applications.

Therefore, for both approaches, validation is a necessary condition when using results at a project management level. Verification is a sufficient condition if experimental data is unavailable or nonexistent.

The objectives of this paper can be broadly stated as follows: (1) Develop a framework that enables generic coupling between an MC code and a commercial CFD software package. (2) Validate the CFD model using representative XE-Prime cases. (3) Couple the CFD model with a neutronic model of an XE-Prime supercell and use the results for cross-verification of a reduced-order model’s thermal-hydraulic and neutronic behavior. Finally, (4) show how computational time can be reduced by using the reduced-order model in a predictor-corrector scheme, i.e., for greatly speeding up the high-fidelity solution.

In this analysis, the reduced-order thermal-hydraulic code is CoRE-ThM (Krecicki et al., 2022/08), which was developed in the Computational Reactor Engineering (CoRE) laboratory at the Georgia Institute of Technology and is the subject of ongoing verification efforts. The CFD analysis is automated with ANTHMM after the CAD model has been meshed. The Monte Carlo code, Serpent (Leppänen et al., 2025), is used for the coupled verification effort. The neutronic model is a reduced-geometry supercell of a full-core XE-Prime benchmark published recently (Eck and Kotlyar, 2026).

2 Description of the high-fidelity coupling framework

The ANTHMM package is a group of functions implemented in the Python programming language that is used to facilitate data transfer between a neutronic solver and Computational Fluid Dynamics software. Specifically, ANTHMM uses the Serpent MC code (Leppänen et al., 2025) to obtain power distributions and Ansys Fluent (ANSYS, 2025, Sec. 23.6) to solve the steady-state RANS equations to obtain thermal-hydraulic solutions.

The package uses these two solvers specifically, but much of the data transfer is standardized using the comma separated value (CSV) file format. The temperature/density interface files used to inform the neutronic simulations are generated by parsing the CSV files from the CFD simulations using cell-specific coordinates, volumes, temperatures, and densities.

3 Geometry of the XE-Prime supercell

The full geometry of the XE-Prime core and the general XE-Prime supercell considered in this analysis are shown in Figure 1. The supercell consists of a hexagonal outer domain enclosing a central moderating element and six primary fueled elements. An additional six partial fueled elements are included in the domain at each outer vertex. The coupled Serpent analysis uses symmetry boundary conditions in the radial directions and vacuum boundary conditions at the upper and lower axial boundaries.

The thermal-hydraulic models are only applied in the active region of the supercell. In the coupled analysis, the supercell model used in Serpent also includes hydrogen plenum regions, shielding, and structural plates which are assumed to remain at constant temperatures and densities.

The central moderating element (Figure 2) comprises of seven radial layers and five materials, extruded axially with a constant cross section for its entire length (Table 1). The Inconel tie rod is included as a structural element to counteract the high axial loads induced on each element that result from the large total pressure drop over the length of the active region of the core (Walton, 1991). Flow is introduced to the annular coolant channel (green color) in a single-pass design from the same plenum that supplies the fueled element’s coolant channels. The inlet pressure for the moderating element’s coolant channel is significantly lower than the fueled element’s channels due to more extreme orificing.

TABLE 1. Global dimensions of the XE-Prime Supercell.

Moving radially outward, there is a thin layer of stainless steel that separates the main flow region from a hydrogen-filled void layer. Pyrolytic graphite insulates the flow region and Inconel tie rod from the extreme heat transfer through the unfueled graphite matrix, which has direct-contact interface walls with the fueled elements. Dimensions for the moderating element are defined in [2, p. I-29] and summarized in Table 2.

TABLE 2. Data for each region of the XE-Prime moderating element.

The fueled elements are structurally simpler than the moderating element. A total of nineteen coolant channels direct hydrogen from the low-temperature plenum to the nozzle chamber at a rate of 32.4 k g / s. The fuel itself is a graphite-uranium matrix.

In contrast to the moderating element, the fueled element has two axial regions, each with the cross section shown in Figure 3. As indicated in Table 1, the final 8.255 cm of the fueled element is an unfueled tip, wherein the graphite-uranium matrix is replaced with the same unfueled graphite material that is used in the moderating element. Dimensions for the fueled element are stated in [2, pp. I21-I24, C-25]. All materials are explicitly included in the Serpent model and reported in Table 3.

TABLE 3. Dimensions for the fueled element.

4 Materials

Appropriate definition of the simulation’s material behavior is challenging due to the range of values reported in the literature and the actual properties of materials used in test articles. In the simulations presented herein, material properties are defined as functions of temperature, subject to the constraints of each solver’s capabilities. Ansys Fluent supports material definitions using polynomial correlations and field functions, whereas CoRE-ThM supports a variety of definitions including polynomial and exponential correlations, as well as properties provided as data points. The differences in the material models supplied to each TH solver are negligible; high-order polynomials are used to fit the available data within CoRE-ThM throughout the expected range of temperatures of each material. In the CFD model, other methods are used for the reasons provided in Section 4.3 to match the reference data with the same level of accuracy.

4.1 Solid material properties

The steady form of the general energy equation is defined in [7, p. 13] in equation 1:Where in v → is the velocity vector, ρ is the material density, h is the enthalpy, k is the material’s thermal conductivity, T is the temperature of the material, and S h represents thermal source terms. For the stationary solid regions, this reduces to the pure conduction as shown in equation 2:

∇ · v → ρ h = ∇ · k ∇ T + S h (1)

∇ · k ∇ T + S h = 0 (2)

Only thermal conductivities are required in solid zones for the steady-state simulations. Temperature-dependent material properties for all solid materials used in this analysis and their sources are given in Figure 4 and Table 4.

TABLE 4. Reference publications for material properties.

4.2 Fluid material properties

The nonzero velocity vector in Equation 1 for the flowing hydrogen coolant requires the inclusion of the specific heat, C P J k g · K, the thermal conductivity, k W m · K, and the density, ρ k g m 3 in the hydrogen’s material specification. The steady momentum equations must also be solved within the fluid domains of the XE-Prime geometry. This introduces the viscosity term, μ k g m · s, in addition the properties specified for the solid domains. Hydrogen material properties were evaluated from (National Aeronautics and Space Administration NASA, 2024). CoRE-ThM automatically calculates hydrogen densities as a function of temperature and pressure; this is not the case for the CFD model.

4.3 Extra considerations

In Ansys Fluent, a combination of field functions, polynomial fits, and piecewise-linear functions are used to define different materials. Among these material models, it is essential that each property’s temperature-dependent function is positive and well behaved within the examined temperature ranges.

An additional point must be made that is specific to Ansys Fluent. Polynomial fits of material properties with complex temperature-dependent behavior are truncated if the polynomials require more than eight coefficients (ANSYS, 2025, Sec. 23.4.2.1). This is the primary reason that piecewise linear functions are used for the fuel region. A field function is used in Ansys Fluent for the pyrolytic graphite to properly include the conductivities of the stainless steel and hydrogen inner gap without explicitly requiring separate mesh zones.

The temperature-dependent density of hydrogen is modeled in Ansys Fluent with the following general field function:

ρ k g / m 3 = A T + B + C T (3)

Note that Equation 3 expresses the hydrogen density as a function of temperature only. Although hydrogen density is generally dependent on both temperature and pressure, separate temperature-dependent correlations were developed (Table 5) for each analyzed operating condition (i.e., power levels) using the corresponding channel inlet pressure. Consequently, the effect of the absolute system pressure is implicitly incorporated into the density model. Within the analyzed cases, the pressure drop along the flow channel is relatively small compared with the absolute pressure level, and therefore its impact on the local hydrogen density is limited. As a result, the density variation is dominated by temperature effects and is adequately captured by Equation 3. The validity of this approximation was verified through comparisons with predictions obtained using the reduced-order CoRE-ThM model, which employs pressure-dependent thermal-hydraulic calculations and produced consistent density distributions throughout the channel.

TABLE 5. Coefficients for the pressure-specific hydrogen density model.

The coefficients in Equation 3 are determined for use in Ansys Fluent for each of the three power factors described in Section 7 and the moderating element by evaluating the densities reported in (National Aeronautics and Space Administration NASA, 2024) at their respective inlet pressures.

These coefficients were determined using the curve fitting function implemented in the SciPy.optimize Python package and documented by Vugrin et al. (2007). Note that these correlations give density values identical to those used by CoRE-ThM throughout the expected temperature range.

The pyrolytic graphite region is assumed to conduct heat isotropically using the temperature-dependent thermal conductivity relation for the basal-plane perpendicular direction given in [9, p. 547].

5 Description of CoRE-ThM

This work also used the thermal-hydraulic solver, CoRE-ThM (Krecicki et al., 2022/08), that provided an additional set of results and was also used to speed up the CFD solution. The primary objective of CoRE-ThM is to enable rapid, computationally inexpensive analysis of a wide range of reactor architectures, with a specific focus on NTP reactors. The most basic component within any CoRE-ThM model is the Element, which is a two-dimensional cross-section of a region within the core. This element contains dimensional and geometric information such as the element shape (e.g., circle, hexagon, etc.) and the number and diameter of coolant orifices. Each region of the element is assigned a material and a power fraction, which dictates the ratio of the heat source that is allocated to the material. Elements are then used to construct a Channel. A channel is composed of an arbitrary number of axial layers of user-defined length, with each layer being a three-dimensional extrusion of a given element.

Each channel is assigned an axial power profile as well as correlations for heat transfer and friction factor calculations. Channels are then dynamically linked into a Flow Scheme. A flow scheme contains all the boundary conditions shared between channels. For example, connected flows can share inlet/outlet pressure and temperature, adjacent flows can share a wall temperature, or flows that are adjacent but antiparallel can share inversed wall temperatures. It also contains other, externally driven boundary conditions, such as the pressure, temperature, and mass flow rates of channel inlets that do not have connected outlets. The flow scheme also determines the order in which the channels are solved as well as the convergence criteria for the model. Figure 5 the basic construction of a CoRE-ThM model.

The computational sequence implemented in CoRE-ThM relies on reducing a model into a 1.5D conduction-convection solution, where radiative heat transfer is not considered. Thus, each channel is approximated as an equivalent cylinder, with each material region in each axial layer being represented as an annulus. Each equivalent cylinder model conserves the element’s key hydraulic quantities, which includes heated and wetted perimeter, flow area, and hydraulic diameter. Previous publications have explored this analysis philosophy in more depth (Krecicki et al., 2022/08).

The CoRE-ThM supercell model is composed of three channels. The first of these is the fueled element, which contains a fueled length of 63 axial layers followed by an unfueled tip of four axial layers. The fueled length and unfueled length each have their own element definitions. The second channel represents the central element minus the tie rod, and the third channel is the tie rod on its own. These were separated because CoRE-ThM, for all channels containing flow regimes, requires flow to be either the innermost or outermost material region, as this simplifies the resistance matrix solution. The dimensions and material properties used for each element match those covered in the previous sections.

In this solution, the inlet temperature, pressure, and mass flow rate are fixed for both the fueled element and unfueled element flow. The fueled element is first solved without any heat flux term from the unfueled element, obtaining the shared wall temperature at the border between the fueled and unfueled element. The unfueled element is then solved with this wall temperature as a boundary condition, which outputs the heat flux into the unfueled element from all surrounding fueled elements as well as the bulk coolant temperature and coolant heat transfer coefficient. These latter two quantities allow for the tie rod solution to calculate the heat flux from the unfueled element into the rod, which is then passed back to the unfueled element as a power sink term. The unfueled element and tie rod solutions are then run until the sink term converges. Upon reaching convergence, the heat flux into the unfueled element is divided by six to calculate the heat flux out of the fueled element. This accounts for the fact that each unfueled element is bordered by six fueled elements. This scaled heat flux term is then passed back to the fueled element, allowing for the next iteration of the full loop. The loop can be rerun until convergence of this heat flux term. Figure 6 shows the solution sequence used for this model.

6 The CFD model

Ansys Fluent is a widely used commercial CFD software that is available at Georgia Tech. Fluent uses a finite-volume based approach to discretizing the transport equations.

6.1 Overview of thermal-hydraulic analysis using Ansys Fluent

For some volume within the mesh, the transport of some variable ϕ can be written in integral form:In Equation 4, ρ is the density of the material, ϕ is the generalized transport variable, v → is the velocity vector, Γ ϕ is the diffusion coefficient for the transported quantity ϕ, and S ϕ is the source term.

∫ V ∂ ρ ϕ ∂ t d V + ∮ ρ ϕ v → · d A → = ∮ Γ ϕ ∇ ϕ · d A → + ∫ V S ϕ d V (4)

Discretization of the above yields the following:By default, the gradient term ∇ ϕ is evaluated by Fluent using a least-squares cell-based method (ANSYS, Inc. 2025). As this analysis seeks steady-state solutions, the time derivative in Equations 4, 5 are replaced with a pseudo-time term to under-relax the equations in each cell.

∂ ρ ϕ ∂ t V + ∑ f N f a c e s ρ f v → f ϕ f · A → f = ∑ f N f a c e s Γ ϕ ∇ ϕ f · A → f + S ϕ V (5)

Diagonal dominance and numerical stability in cell equations is improved using the default global time-step method within Ansys Fluent to relax the governing equations. The pseudo-transient formulation of Equation 5 in a given cell is given in ANSYS, 2025 (Sec. 23.6)

ρ p Δ V ϕ P − ϕ P 0 Δ τ + a P ϕ P = ∑ n b a n b ϕ n b + b (6)

The linearization of Equation 6 employs the a P and a n b coefficients to describe the transport properties of the current cell and neighbor cells respectively. The coefficients themselves vary depending on the discretization scheme (Patankar, 1980).

Rearranging Equation 6 yields the following:

ϕ p = ∑ n b a n b ϕ n b + b ρ p Δ V Δ τ + a p + ϕ p 0 1 + Δ τ a p ρ p Δ V (7)

One sees that for very small values of Δ τ, the current value of ϕ p becomes mostly dependent on ϕ p 0, the previous-iteration value of ϕ. Similarly for very large values of Δ τ, Equation 7 becomes the steady-state solution for ϕ (Patankar, 1980). By default, the pseudo time step itself is calculated automatically within Fluent using the conservative time step calculation method; this behavior is changed in the present analysis to the aggressive time step calculation method to reduce computation time (ANSYS, 2025, Sec. 23.4.2.1).

The temperature-dependent density term is evaluated in each cell using arithmetic averaging of neighboring cells’ densities (ANSYS, 2025, Sec. 23.4.2.1).

Equation 7 is cast into the Reynold’s-Averaged Navier-Stokes (RANS) formulation, requiring additional turbulence modeling equations to capture the turbulent fluctuations of velocities and temperature fields in the fluid domain [(ANSYS, 2025, Sec. 23.3.3.4; Tu et al., 2024), pp. 98–99]. The turbulence equations applied in this analysis are the shear stress transport (SST) k − ω equations. Closure of the turbulence equations is obtained through the use of near-wall models to adjust the near-wall viscosity and determine the viscous stress terms. In the near-wall viscous layer of the fluid domains, the enhanced wall treatment (EWT) formulation for the momentum and energy equations is used to ensure y + insensitivity. Equation 8 stands for the momentum equations,Wherein u l a m + and u t u r b + , described in equation 9, are the laminar and turbulent nondimensional velocities given in the standard two-layer wall functions:

u + = e Γ u l a m + + e 1 / Γ u t u r b + (8)

u l a m + = y + (9)

u t u r b + = 1 κ ln E y +

The blending function, Γ, is defined in references (Tu et al., 2024; Kader, 1981) as Equation 10:

Γ = − 0.01 y + 4 1 + 5 y + (10)

The energy equation receives similar treatment (Tu et al., 2024; Kader, 1981; ANSYS, 2025, Sec. 23.6.2).

Two fundamental solver types are available within Ansys Fluent, namely, the pressure-based and density-based solvers. In the pressure-based solver, the pressure field is solved directly; the density-based solver uses the continuity equation to obtain the density field and solve for pressure using equations of state. The fully coupled pressure-based solver is the default solver type within Ansys Fluent. Of the two models, the pressure-based model is applicable to incompressible or low compressible flow, whereas the density model is generally used for compressible flow (ANSYS, 2025, Sec. 4.1.1). This analysis assumes all flow regions are incompressible, as the in-channel-shocks would damage the core itself. Solver settings are given in Table 6.

TABLE 6. Settings and numerical schemes for the CFD model.

6.2 The reduced CFD geometry

Symmetry planes reduce the total number of mesh cells required; these symmetry planes assume zero flux in the wall-normal directions. Figure 7 shows a cross-sectional view of the 3-Dimensional CFD mesh. The total CFD mesh accounts for one-sixth of the moderating element and one-sixth of a surrounding fuel element. Power is downscaled to conserve power densities. Mass fluxes are held constant among all fueled element channels by taking the inlet area fractions of each channel inlet to that of the full channel’s inlet area. In the moderating element, the mass flow rate is similarly scaled by a factor of six.

The y + insensitivity of the blending function used in the enhanced wall treatment is validated experimentally (Kader, 1981). In addition, a mesh sensitivity study was performed for the flow channels using multiple y + values and mesh densities during the development of the final CFD model. For the final geometry and material property functions adopted in this work, the predicted thermal-hydraulic quantities remained unchanged within engineering significance as the mesh was refined. This behavior is consistent with the expected y + insensitivity of the enhanced wall treatment and confirms its applicability to the present analysis. The good agreement observed in the temperature profile comparisons presented in the Results section further supports the adequacy of the selected fine mesh for the analyses performed in this work.

The temperature-dependent material properties of the pyrolytic graphite and the addition of source terms yield large thermal gradients and, consequently, a large range of thermal conductivities within the very small thicknesses of the pyrolytic graphite, inner hydrogen gap, and stainless-steel region as shown in Figure 8.

Thin material regions within the moderating element cause some difficulty, as cell refinement in these regions can greatly increase the total number of cells in the mesh. Due to the lack of available documentation characterizing heat transfer within the hydrogen gaps, the hydrogen in these gaps was assumed to be stagnant. Consequently, both the CFD and CoRE-ThM models considered heat transfer across the hydrogen gaps solely through conduction, while radiative heat transfer was neglected. Inclusion of radiative heat transfer would increase the effective thermal conductance across the gaps, thereby reducing the thermal resistance between adjacent components and potentially lowering local temperature gradients and peak temperatures. However, the impact is expected to be localized to the gap regions and is not anticipated to alter the overall conclusions of this work.

The CFD mesh assumes that the pyrolytic graphite zone extends into the spatial region where the inner gap regions are located, but material properties are adjusted at the cell level to maintain the appropriate spatial and temperature-dependent conductivities across the gap region. Mesh sensitivity studies were performed as part of the model development process to ensure that the numerical solution was independent of the spatial discretization. Since further mesh refinement produced only negligible changes in the quantities of interest, these results are not presented here for brevity.

Uniform inflation layers at the interface between the pyrolytic graphite and moderator coolant channel zones are used to capture these two materials, placing two cell centroids inside the regions occupied by the hydrogen gap and the stainless steel liner. This approach properly models the thermal gradient at this interface without requiring independent zones for these two thin regions. During postprocessing, the hydrogen gap and stainless steel layers are filtered out of the resulting data to preserve an accurate representation of the pyrolytic graphite’s volume-averaged temperatures. Material properties in the pyrolytic graphite zone are specified as a function of radial distance from the Inconel tie rod’s centerline as described in equation 15:

k p y r o l y t i c = i f R c e l l > R H 2 G a p , k p y r o l y t i c , i f R c e l l < R H 2 G a p , k s t a i n l e s s , k H 2 (15)

One-dimensional conduction zones are created on wall boundaries to approximate the location of the source terms in the coupled verification. The CFD mesh consists of 3,968,070 cells, twelve cell zones, and fifteen shell conduction regions.

6.3 Coupling scheme with CoRE-ThM

During system coupling with Serpent, the CFD model uses fully-converged coupled results from CoRE-ThM to initialize the required temperature and density fields. Figure 9 shows CoRE-ThM outside of the main iteration loop; however it uses a similar loop outside of that shown to produce first-iteration temperature and density predictions in each material for use in the coupled CFD analysis. This approach, using CoRE-ThM as an initializer, is hypothesized to decrease the total number of iterations and time to obtain the higher-fidelity coupled reference solution. However; fully-converged results can be obtained without the use of CoRE-ThM.

All steps shown in Figure 9 are automated by the defining Python simulation file. Detector tallies, power profiles, and CFD results are stored as ASCII Python dictionaries for postprocessing in addition to being used as input parameters in the simulation itself. Power profiles are determined for input into Ansys Fluent using the method described in the boundary conditions section of this dissertation. Once the CFD simulation has completed, cell data for each material zone is exported to.csv files. For purely axial temperature/density specification, the CFD data is volume-averaged and Serpent IFC files (Leppänen et al., 2025) are written for each material.

7 Boundary conditions

Experimental results for the fueled elements are given in Westinghouse Astronuclear Nuclear Laboratory WANL and Reese (1967). for three locations within the core. Each of these regions (see Figure 10) are characterized by their radial power factors, namely, 0.8, 1.0, and 1.2. The same axial power profile is applied to each of these cases, but the total power within the system is scaled from the average using the radial power factor. Thus, the 1.0 radial power factor case represents the steady-state, core-averaged behavior.

In the moderating element, full axial temperature distributions are documented only for the 1.0 radial power factor case; outlet temperature correlations are, however, provided in [2, p. I-31]. Boundary conditions for solid regions are given in the form of inlet face temperatures and shown in Table 7. In the moderating element, inlet face temperatures are assumed to be those of the 1.0 radial power factor case for all cases.

TABLE 7. Summary of boundary conditions for solid domains.

Mass flow rates are scaled directly (Table 8) at channel inlets in the fueled element by the value of the radial power factor. In the moderating element, inlet mass flow rates remain constant among all cases. Inlet pressures and temperatures (Table 8) are included for all channel inlets, both in the moderator and the fueled element.

TABLE 8. Summary of coolant inlet boundary conditions.

The axial power profiles of Figure 11 are directly provided in [2, pp. E−49, I-16]. Power source terms are only applied in the fueled element and moderating element’s unfueled graphite region in the uncoupled validation.

Coupled simulations collect neutron and photon heating tallies within Serpent and use these tallies to determine the power shares in each material. Power distributions are assumed to be radially invariant. The total power within the supercell is determined by the average full-core conditions given in (Westinghouse Astronuclear Nuclear Laboratory WANL and Reese, 1967). As power shares in the fueled and moderating elements are determined by Serpent, the total power in the supercell is conserved, but these power shares deviate slightly from the published data.

Tallies are collected for all materials in the supercell and averaged radially to find material-specific power shares. The CFD model requires power densities in units of watts per cubic meter. Axial power density distributions are obtained using the axial power profile, whether that profile is in the form of a normalized curve (in the uncoupled case) or from a normalized axial tally array. The peak power density value is found using the cross-sectional area of each material region by integrating over the length of the material. The peak power density known, the normalized data is converted to a profile file by the ANTHMM Python package for use in the Fluent simulation.

8 Results

Reference data is published explicitly in (Westinghouse Astronuclear Nuclear Laboratory WANL and Reese, 1967) for the fueled element’s regions at each of the radial power factors studied. Only the radial power factor of 1.0 is published for the moderating element, but the CFD model shows excellent agreement in the key regions for this radial power factor. Reference figures for materials in the moderating element show only the 1.0 radial power factor case.

The temperature profiles in the fueled element’s channels shown in Figure 12a vary only slightly between cases. Differences in outlet temperatures can be observed, though these differences are only a few Kelvin. Maximum deviations from the reference are in the near-inlet region, likely due to increased mixing and improved heat transfer that would be difficult to accurately directly measure or capture without a turbulence model. Some evidence of this comes from the CoRE-ThM curve - CoRE-ThM uses a heat transfer coefficient that accounts for turbulence in principle but is still reduced-order fidelity in nature.

The fueled region (Figure 12b) also shows insensitivity to the radial power factor, temperatures between the three cases are nearly identical. Deviations from the reference are most extreme in the section that extends from 20 cm to 35 cm. It is believed that the 10% uncertainty in the axial power distribution [2, p. I-6] could be the cause of this difference, or that the flow is not yet fully developed. During this analysis, the near-inlet bump seen in this region was shown to be sensitive to the power discontinuity at the inlet face; this suggests that the power profile may be reported with a downstream shift.

As the unfueled graphite has a relatively large thermal conductivity and is in direct contact with the fueled element, its axial temperature Figure 13a, profile is largely dependent on fuel temperatures. CoRE-ThM and the CFD model accurately capture the temperatures near the outlet and across the entire length of the graphite region with a similar temperature deviation from the reference that was shown previously in the fueled region.

The pyrolytic graphite region (Figure 13b) is intended to insulate the inner tie rod and is therefore characterized by extreme radial temperature gradients. At the outlet plane, a 0.25 mm distance corresponds to a temperature difference of 50 K. This requires the use of the composite, piecewise material property specification in the CFD mesh described previously, as neither the inner gap nor the stainless steel layer can be explicitly meshed as independent zones without incurring an enormous cell count penalty.

The moderating element’s coolant channel Figure 13c, is thermally downstream of all other regions except for the Inconel tie rod. Because of this, any errors in the temperature fields of the surrounding materials propagate inwards. Validating the temperature profile in this region is key to ensuring the heat transfer from outer regions is accurately modeled. The CFD model and CoRE-ThM are within 5 and 10 K (respectively) of the reference curve, indicating that the applied material models and heat transfer methods have a propagated thermal error that is well within the range of the uncertainty of the reference data.

Finally, the Inconel tie rod region is the region that is most susceptible to propagated error as shown in Figure 13d. The CFD model underpredicts the temperature at the tie rod’s outlet plane by 30 K, while CoRE-ThM matches the outlet temperature after a period of significant underprediction. This uncoupled validation assumes negligible power in the Inconel tie rod due to a vacancy in power information in (Westinghouse Astronuclear Nuclear Laboratory WANL and Reese, 1967) for this region; however, the reference data implicitly includes this missing source term. Comparison of the temperature profiles in the moderating element’s coolant channel and the tie rod in the CoRE-ThM model reveals that these temperatures are reported as identical. Because CoRE-ThM cannot capture the thermal boundary layer on the interfaces of the moderating element’s channel, it overpredicts the heat transfer through the coolant channel’s thickness.

The CFD model and CoRE-ThM are both within the margin of error specified in (Westinghouse Astronuclear Nuclear Laboratory WANL and Reese, 1967) and have been proven to predict nearly-identical temperature profiles for every material region of interest.

8.1 Coupled verification

Differences in power distributions between the coupled CFD and CoRE-ThM simulations are given in Table 10, Tables 11, 12 for each of the radial power factors studied. The total power of the entire system is preserved, with the powers in each fueled element and each moderating element summing to their expected values given their radial location within the XE-Prime core.

Reported power is for each hexagonal element. The first iteration of each coupled case uses CoRE-ThM’s output as a temperature and density initializer in Serpent. As the CoRE-ThM simulation is fully converged, the statistical error columns in Table 10, Tables 11, 12 provide an estimate of statistical uncertainty. Relative error percentages are calculated with respect to the CFD case–positive values represent more power in the CFD simulation and negative values represent more power in the CoRE-ThM simulation.

TABLE 9. Power shares in the fueled and moderating elements.

TABLE 10. Differences in power (kW) for each material; 0.8 RPF.

TABLE 11. Differences in power (kW) for each material; 1.0 RPF.

TABLE 12. Differences in power (kW) for each material; 1.2 RPF.

Power shares are well preserved across all cases. Differences occur in those materials wherein the statistical uncertainty from the MC simulations are high, such as the regions containing hydrogen and the regions that are located several mean free paths from the materials wherein most of the power is generated.

The bulk power shares resulting from the MC simulation are shifted slightly to the fueled regions when compared to the true full-core values given in Table 9. This effect has implications for the thermal-hydraulic analysis.

Similar trends are seen in the coupled simulation as in the uncoupled simulation when comparing axial temperature profiles as shown in Figure 14.

Excellent agreement is displayed between CoRE-ThM and the CFD model in the fueled element (Figure 14). Note that the near-inlet temperature gradient in the fuel is much greater than that shown in the nearly linear temperature rise of the uncoupled case. This occurs because of the upstream shift of the fuel region’s power profile. The characteristic bump in the fuel’s axial temperature profile of the reference curve is more accurately captured by the CFD model, especially in the 1.2 RPF case. CoRE-ThM represents this bump also, but with less magnitude.

Temperatures in the radially intermediate regions of the unfueled and pyrolytic graphite are equally well-matched as demonstrated in Figure 15.

In the innermost regions (shown in in Figure 16), CoRE-ThM reports the same outlet temperature overpredictions (shown in in Figure 16) as displayed in the uncoupled cases for the same reasons therein described. The presence of the source term in the coupled analysis amplifies this error, as there is no path for heat to propagate out of the Inconel. Recall that CoRE-ThM sets the wall temperature of the Inconel equal to the bulk temperature of the moderating element’s channel region; in the coupled case, CoRE-ThM stores more energy within the Inconel, causing a temperature rise in both regions.

This result can be confirmed by including the reference curve for these regions (Figure 17). The addition of source terms in the Inconel slightly increases the temperature in the moderating element’s flow channel for the CFD case, but the coupled CoRE-ThM curve shows a dramatically higher temperature than its uncoupled counterparts. The reference curves show the behavior of the full core and thus do not indicate the as-modeled effect of the power shift away from the moderating materials. The CoRE-ThM temperature overprediction can therefore be attributed to a breakdown in assumptions about heat transfer through the annular channel.

Across all domains, the maximum temperature deviation between the two thermal-hydraulic models is in the Inconel region, with a 55K overprediction of the outlet temperature in CoRE-ThM’s results for the 1.2 RPF case. As this is a very small portion of the total domain, it has been shown that there is no effect on the neutronic results.

The primary benefit of reduced-order solvers such as CoRE-THM is the avoidance of repeated CFD evaluations, thereby substantially accelerating the coupled neutronic–thermal-hydraulic calculation. The execution time of CoRE-THM is on the order of several seconds and is therefore negligible compared to the CFD execution time. Consequently, the computational cost of the coupled iteration is effectively determined by the neutronic and CFD simulations, with CoRE-THM contributing an insignificant computational overhead.

Similarity of the thermal-hydraulic feedback from CoRE-ThM allows its solutions to be used as a first-iteration predictor of densities and temperatures when coupling the CFD model. Two cases were compared using the 1.0 RPF case. The non-initialized case assumes constant temperatures of 400 K for every material and constant, nominal densities in the hydrogen regions. The initialized case uses the results of the 1.0 RPF CoRE-ThM simulation as an initial guess in an attempt to reduce the total number of iterations required.

Figure 18 indicates the magnitude of the error in the initialized case for the first iteration is several times less than that of the un-initialized case. Additional iterations would likely not be necessary for simulations using CoRE-ThM as an initializing tool. Results are compared using the axial variation of the FE coolant channel density. For each coupled iteration of the two cases, the densities are normalized by the fifth and final iteration’s density profile. As the number of iterations increases, this value approaches unity for the entire domain.

The convergence of the coupled calculations was evaluated by monitoring the axial hydrogen density distribution between successive Picard iterations. Since hydrogen density influences the neutronic solution, which in turn determines the power distribution and temperature field, convergence of the density profile implies convergence of the overall coupled neutronic–thermal-hydraulic solution. As shown in Figure 18, the maximum deviation from the converged density profile decreases from approximately 6% in the first iteration to below 0.3% after three iterations, indicating rapid convergence of the coupled framework.

The final, steady-state solution that is independent of CoRE-ThM influence. Review of Table 11 shows that power deviations in the fuel material at the 1.0 radial power factor are within about 6 W (out of ∼0.69 MW); this deviation corresponds to the small density errors shown in the fourth iteration of Figure 18.

9 Conclusions and future work

The objective of this work has been to develop a framework for modeling coupled thermal-hydraulic feedback on neutronic solutions obtained with an MC code. The intended use of the ANTHMM Python package is to provide reference solutions where data is scarce, to increase the analytical capabilities of its users, and to find potential weaknesses in assumptions. This publication provides an example of how this package can be used. Once a CFD model has been generated for a geometry of interest, it can be validated using explicit boundary condition definitions. Coupling the CFD model to a neutronic code allows the thermal feedback mechanisms to be accurately captured at runtime during the simulation. Results from coupled simulations can be used to identify weaknesses in the assumptions used for reduced-order codes. Finally, large simulation run times can be decreased if results from reduced-order codes are used to initialize the thermal fields for the first iteration of the neutronic simulation.

The work presented herein sits within the broader field of computational reactor engineering, where modeling challenges are often unique to each reactor concept and system design. Within this field, CFD remains an essential tool for resolving complex thermal-fluid phenomena; however, its high computational cost continues to be a significant limitation. The methodologies described in this study are intended to alleviate this burden through the use of reduced-order approaches while maintaining predictive capability. In parallel, recent studies (Guo et al., 2026; Lu et al., 2024) have demonstrated the successful application of porous-media representations for modeling complex reactor components that would otherwise require prohibitively expensive geometric resolution. Additional strategies, including reduced-order and surrogate models (e.g., simplified pump models), can further complement CFD analyses through the use of periodic boundary conditions, preprocessing techniques, and other approximations that reduce the need for extensive mesh refinement while preserving the key physical behavior of the system.

The ANTHMM simulation framework and that of CoRE-ThM are validated against the experimental reference data and then these codes are cross-verified. The expected error in CoRE-ThM has been quantified with the ANTHMM reference solution. Neutronic behavior in the coupled case is shown to be insensitive to differences between the models, enabling the use of CoRE-ThM as an initializing tool in a predictor-corrector scheme. Further expansion of the ANTHMM Python package will extend its capabilities and applications. Radiative heat transfer is not included in the current model, but further development of ANTHMM package may increase the fidelity of results. As interest in NTP grows both in industry and academia, an increased number of institutions will seek higher-fidelity solutions using larger CFD models and higher resolution of coupled effects. Working with common and capable software such as ANSYS Fluent enables implementation of thermomechanical models which require an extension of current methods to neutronic solvers that support deformable geometries.

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.

Footnotes

Edited by: Mark D. DeHart, Abilene Christian University, United States

Reviewed by: Heng Xie, Tsinghua University, China

Yeming Lu, Dalian University of Technology, China