Thermal-Hydraulic and Solid Mechanics Multiphysics Safety Analysis of a Heavy Water Reactor with Thorium-Based Fuel

1. Introduction

Thorium has been viewed as a prospective reactor fuel since the nuclear power sector emerged in the 1950s, yet—despite worldwide initiatives—its advancement has consistently lagged behind that of uranium-based fuels. Today, interest in thorium has been risen again due to the drive to develop safe, clean nuclear fuels capable of addressing modern power generation challenges, and its abundant reserves further underscore its potential to meet the rising global demand for energy [1,2]. Consequently, significant research has focused on investigating the replacement of uranium fuel with thorium as a means to reduce reliance on the scarce uranium fuel. Monte Carlo N-Particle eXtended (MCNPX) code version 2.7, with the ENDF.VII cross-section library, was used to study the neutronic performance of thorium-based fuel types [(_0.944Th, U)O₂, (_0.955Th, ²³³U)O₂, and (_0.934Th, rgPu)O₂] as alternatives to UO₂. The study demonstrated good viability by achieving acceptable safety parameter values and providing a favorable power distribution throughout the fuel assembly compared to UO₂ [3]. Incorporating 15% TRUO₂ into an 85% Th-²³³UO₂ fuel for Pebble Bed Modular Reactor (PBMR-400) optimizes neutronic performance by extending the fuel cycle, offer potential waste reduction and show good environmental benefits [4]. Using a 3D heterogeneous MCNPX model, the study evaluated thorium fuel for the Gas Cooled Fast Reactor (GCFR-400) found that the proposed fuel exhibits promising neutronic performance with satisfactory control and safety parameters throughout its burnup [5]. Thorium-based fuels—whether blended with High Enriched Uranium (HEU), Low Enriched Uranium ( LEU), or ²³³U or formulated as MOX/mononitride—deliver acceptable thermal neutron fractions for flat power density, lower fuel temperatures, improved reactivity coefficients, reduced plutonium and minor actinide production, and enhanced proliferation resistance [6,7,8,9,10]. The analysis of (Th_xU₁₋ₓ)N fuel in both infinite and compact reactor models indicates that thorium additions maintain criticality, extend reactor service life by reducing k_eff decline and enhancing ²³³U production and introduce a potential diversion risk due to proliferation vulnerabilities [11]. Based on state-of-the-art reactor designs, it is concluded that thorium-based molten salt small modular reactors represent the most promising technical pathway to advance nuclear energy into a new era—meeting market demands while ensuring maximum safety at an affordable capital cost [12].

Designing a reactor core that meets operational requirements remains one of the most significant challenges in reactor physics, where the quest for effective core optimization methods that ensure safety and improve economic efficiency has long been a central focus. Additionally, simulating nuclear fuel behavior is crucial for licensing and interpreting experiments with new fuels, as it involves evaluating temperature distributions, velocity profile and stresses. Finite Element Modeling (FEM) plays an improtant role in the design, safety analysis, and operation of nuclear reactors, because it incorporates multiple physics modules such as thermal, structural, and fluid flow. A fully coupled, time-dependent multi-physics model was developed using COMSOL Multiphysics to analyze the single-channel behavior of a Lead Fast Reactor (LFR) under both steady-state and transient conditions [13]. The model integrated neutronics, thermal-fluid dynamics, and thermal-elastic effects, enabled accurate prediction of reactivity feedbacks without relying on corrective factors, and was validated against Monte Carlo simulations using SERPENT. SERPENT is a 3D, continuous-energy Monte Carlo code for neutron and photon transport, developed at VTT Technical Research Centre of Finland since 2004. It supports stand-alone burnup simulations, lattice physics, group-constant generation, and multi-physics coupling for reactor analysis, shielding, fusion applications, and more. A steady-state thermal–mechanical safety analysis of a 0.5 MWth heat pipe-cooled, graphite-moderated micro reactor were provided using ANSYS Fluent and ANSYS Mechanical, focusing on temperature, heat pipe loads, and structural stresses during normal operation and heat pipe failure scenarios. The analysis showed that the design remains within safety limits under most conditions, and two design modifications were proposed—thermal paste and non-uniform heat sink conditions—to reduce thermal stresses in critical failure cases [14]. A two-dimensional thermomechanical-coupling model of a heat-pipe-cooled space-reactor core was built in COMSOL multiphysics to calculate key heat-transfer and structural-mechanics parameters under both steady-state and accident conditions and to evaluate the core’s thermal safety performance. The analysis shows that failure of a single heat pipe at any location poses no safety risk. However, thermal stress exceeds the allowable limit if two heat pipes fail at the core edge, and with three failed heat pipes it surpasses the limit even in the well-cooled central region, raising the local maximum temperature to 1959.8 K [15,16]. Coupled neutronics, thermal-hydraulics, and thermo-mechanics steady-state COMSOL analysis of a molten salt fast reactor enabled an accurate assessment of reactivity feedbacks and fuel behavior in a fully liquid-fueled fast reactor configuration[17,18]. Although extensive thermal-hydraulic and solid mechanics safety analyses have been conducted for CANDU-6 reactors—as shown in studies on pressure tube creep effects, critical channel power reduction, and high-temperature deformation behavior [19,20,21,22]—there remains a compelling need to provide new and updated safety studies. For instance, existing models often rely on historical creep correlations and coolant behavior under idealized or steady-state conditions. However, evolving reactor operation strategies and potential introduction of alternative fuels such as thorium demand a reassessment of heat transfer, mechanical stress, and safety margins. Studies like those evaluating thorium-based fuel channels (e.g., TACR-1300) demonstrate promise but remain limited to preliminary or simplified models, lacking full multi-physics coupling or transient condition analysis. Furthermore, modern computational tools such as COMSOL Multiphysics or advanced Computational Fluid Dynamics solvers can now incorporate temperature-dependent material properties and real-time coolant property variations more accurately.

In this study, COMSOL multiphysics was employed to carry out coupled thermal-hydraulic and structural analyses, evaluating the safe operating limits of a heavy-water reactor fueled with thorium. A transient model combined the Heat Transfer in Solids and Fluid Flow interfaces for heat-conduction and fluid-dynamics calculations, and the Solid Mechanics interface for the structural response. The simulations produced temperature fields, pressure and velocity distributions, as well as displacement, strain, and geometric deformation results.

2. Materials and Methods

Thermal-hydraulic and solid-mechanics calculations can be performed using the radial power distribution obtained from the neutronic analysis. The results of the neutronic calculations reported by [23] were used as input for our calculations. The reactor core is housed within a horizontal cylindrical vessel made of stainless steel. Inside this vessel are 380 fuel channels submerged in heavy water, which serves as a neutron moderator. The fuel channels are arranged in a square lattice with a spacing of 28.575 cm (Figure 1). Each channel includes a pressure tube made of zirconium–niobium alloy, encased within a zircaloy-2 calandria tube, and connected at both ends by stainless-steel fittings. Spacers are used to support the pressure tube and maintain its separation from the calandria tube. The space between these tubes is filled with carbon dioxide gas, serving as thermal insulation to reduce heat loss. Every fuel channel holds 12 identical fuel bundles, each 49.53 cm long, with each bundle containing 37 fuel rods arranged in a circular ring configuration. Each fuel rod consists of 30 fuel pellets enclosed within a cladding tube made of zirconium alloy. The design parameters and dimensions of reactor core can be found from Table 1.

Because all fuel bundles and channels in the core are symmetrical, a single fuel bundle was modeled to represent the neutronic behavior of the fuels under investigation. Figure 2 illustrates Pin power distribution through fuel bundle, which is obtained from neutronics calculations.

From the power distribution, we can identify the hottest rod, where the most limiting and safety-critical conditions in the reactor core are present. The aim is to demonstrate these conditions through simulations. For comparison, single fuel rods operating at the maximum, median, and minimum power levels were individually analyzed (Table 2). The initial step in the solution process involves breaking down the geometry into smaller elements, a procedure known as meshing. After that, the partial differential equations (PDEs) are mapped onto this mesh. Free triangular fine mesh was used for simulations and presented in Figure 3. The reason for using triangular elements is that they can conform more easily to complex geometries, such as curved boundaries, fuel rods, cladding layers, and coolant channels. Therefore, triangular meshes adapt better than quadrilateral meshes.

(Th, ²³³U)O₂ was selected as the fuel, and its relevant physical properties were taken from the literature. COMSOL Multiphysics can integrate heat-transfer, fluid-dynamics, and solid-mechanics physics into a single model. This multiphysics coupling offers the most faithful representation of real-world behavior, requiring far fewer simplifying assumptions than typical analytical approaches. Three physics modules are used as the following:

2.1. Heat Transfer in Solids

The heat-generation PDE is addressed with the FEM. Since a thermal power used for the heating cross surface of the core, the transient Fourier heat conduction equation was solved to find the spatial distribution of temperature. The PDE is then discretized over the mesh, and solving it provides the temperature distribution within the fuel, cladding, and coolant. Except the top, bottom and cross surfaces, all other boundaries are assumed to be thermally insulated.

For the thermal properties of fuel[24] and clad materials[25] the following expressions were used:

Specific heat capacity of (Th, ²³³U)O₂:

$${\mathrm{C}}_p=239.9+8.512\times {10}^{-5}T-5.754\times {10}^{-8}{T}^2+1.528 \times {10}^{-11}{T}^3 \left[{\displaystyle \frac{J}{mole K}}\right]$$

(1)

Thermal conductivity of (Th, ²³³U)O₂:

$$k=14.242-0.022T+1.764\times {10}^{-5}{T}^2-6.380\times {10}^{-9}{T}^3+8.038\times {10}^{-13}{T}^4 [W/m K]$$

(2)

Density of (Th, ²³³U)O₂:

$$\rho = 9326 – 0.0004\mathrm{T}+2.54\times {10}^{-7}{T}^2-7.897\times {10}^{-11}{T}^3 [\mathrm{k}\mathrm{g}/{\mathrm{m}}^3]$$

(3)

Specific heat capacity of zircalloy-4:

$$C_p=7.1\cdot {10}^{-2}+1.7\cdot {10}^{-5}\cdot T-0.89\cdot {10}^{-3}\cdot T^2 [{\displaystyle \frac{cal}{g}}\cdot K]$$

(4)

Thermal conductivity of zircalloy-4:

$$0.11804+0.00012\cdot T-4.33441\cdot {10}^{-8}\cdot T^2+5.90277\times {10}^{-11}\cdot T^3 [W/cm °C]$$

(5)

Density of zircalloy-4:

$$\rho = 10950 - 0.35\cdot \mathrm{T} - 1.2\cdot 10^(-4) \cdot T^2 [\mathrm{k}\mathrm{g}/{\mathrm{m}}^3]$$

(6)

2.2. Fluid Dynamics

The Navier–Stokes PDE for momentum and mass conservation are solved, yielding the pressure drop of the coolant along the axial (z) direction and the corresponding velocity field. When this turbulent-flow (k–ω) model is coupled with the heat-transfer physics, two key outputs emerge: 1)The convective heat-transfer coefficient h, which depends on the coolant velocity, the channel geometry (velocity-independent), and the temperature-dependent coolant properties, 2)The axial pressure drop obtained from the Navier–Stokes solution. Combining this pressure drop with the temperature profile and local heat flux from the heat-generation equation makes it possible to evaluate the coolant’s distance from nucleate-boiling conditions at any location along the z-axis. The water inlet and outlet boundaries are the top and bottom surfaces, respectively.

2.3. Solid Mechanics

Within the solid-mechanics interface, the model solves the governing elasticity equations (analogous to Navier–Stokes for solids) to determine the thermal-induced deformation, displacement magnitude, and strain in the fuel and cladding, fully coupled to the heat-transfer physics. Top and bottom boundaries are fixed to reflect the physical constraints imposed by the reactor’s structural support system.

For mechanical properties, the following thermal expansion expressions for fuel [26] and clad [24] were used:

$$\begin{gathered} a_{fuel}=5.416\times {10}^{-6}+3.47\times {10}^{-9}T+3.0107\times {10}^{-15}{T}^2-1.2212\times {10}^{-15}{T}^3+4.181\times {10}^{-19}{T}^4 [1/\mathrm{K}] \\ {a}_{clad}=\left(5.22+1.82\times {10}^{-3}T\right)\times {10}^{-6}[1/\mathrm{K}] \end{gathered}$$

(7)

Young’s Modulus expressions for fuel and clad materials [24]:

$$\begin{gathered} E_{fuel}=259.87\times {10}^9\times {(1-P}^2) \left[\mathrm{P}\mathrm{a}\right] \\ \mathrm{P} \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{p}\mathrm{o}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{y} \\ {E}_{clad}=9.21\times {10}^{10}-4.05\times {10}^7\times T \left[\mathrm{P}\mathrm{a}\right] \end{gathered}$$

(8)

The effects of burn-up, plutonium content, pellet cracking, swelling, oxygen-to-metal ratio, and gaseous fission product release on the thermal properties of the fuel are neglected. It is known that elementary particle irradiation affects the mechanical and structural properties of materials [27,28,29]. In this work, irradiation effects are also neglected due to the lack of available literature data.

3. Results and Discussion

To ensure safe reactor operation, the fuel’s centerline temperature must remain below its melting point, and its average temperature must stay under 1400 °C [30]. The spatial distribution of temperature after application of the maximum heat source for different times(from 0 to 3 seconds) is shown in Figure 1a. When the heat is transferred from the inner surface of the fuel material to its outer, the temperature changes parabolically. The temperature profile along the y–x direction is presented in Figure 1b. As evident from the figure, the reactor core exposed to a heat source for a longer period attains a higher temperature than when the exposure time is shorter. A similar temperature increase in the fuel and cladding materials as a function of time was also observed in the coupled neutronic, thermo-elastic, and fluid-dynamic time-dependent model for the single-channel transient analysis of a Lead Fast Reactor [13].

Figure 1. The temperature distribution through different materials of the fuel rod after application of the maximum heat source for different times(from 0 to 3 seconds) (a). The temperature profile along the y-x direction (b). 0 s means the beginning of the simulation time.

Figure 1. The temperature distribution through different materials of the fuel rod after application of the maximum heat source for different times(from 0 to 3 seconds) (a). The temperature profile along the y-x direction (b). 0 s means the beginning of the simulation time.

Figure 2 demonstrates the temperature values at 3 s after applying the maximum, median, and minimum heat sources. Overall, the temperature decreases radially from the center of the fuel rod to the coolant boundary. The centerline and average fuel temperatures are below the safety limit. A similar radial temperature distribution in the fuel rod, with a further decrease in the cladding, was observed in the thermo-mechanical analysis of a heat pipe–cooled space reactor [15].

Figure 2. Radial temperature distribution at 3 s after applying the maximum, median, and minimum heat sources.

Figure 2. Radial temperature distribution at 3 s after applying the maximum, median, and minimum heat sources.

Figure 3 shows the maximum volumetric strain for the fuel material in axial direction. After 3 s of heating with volumetric source of 245 MW m⁻³, the 2-D strain map shows a characteristic high-strain core surrounded by a low-strain rim (Figure 3a). Axial line profiles (Figure 3b,c) taken through the pellet mid-radius and at the fuel–cladding interface confirm that the strain is remarkably uniform along the central ~ 0.5 m rod, reflecting the nearly homogeneous internal power distribution. Peak volumetric strain at the center rises from 0.006 to 0.014, respectively, increasing almost linearly with heat generation. Only near the pellet end faces does the strain decrease slightly, a result of end-face mechanical constraints and local heat-flux gradients. An increased volumetric strain in the fuel compared to the cladding material was also observed in the literature [31].

Figure 3. Volumetric strain for the fuel after 3s of heating at three different heat sources. (a) 2D volumetric strain map, (b) strain at the center and (c) near the end.

Figure 3. Volumetric strain for the fuel after 3s of heating at three different heat sources. (a) 2D volumetric strain map, (b) strain at the center and (c) near the end.

The two-dimensional map of total displacement shows minimal expansion in the central core (Figure 4a) where the cladding constrains radial motion, and maximum displacement along the outer surface where thermal expansion is less restricted. Profiles extracted across the diameter confirm that surface displacement increases systematically with the applied heat source—rising from ~0.004 mm at 119 MW m⁻³ to ~0.012 mm at 245 MW m⁻³—while the center remains comparatively immobile. Axial centerline profiles demonstrate nearly uniform elongation along the heated length, but with magnitudes that also scale with heat source, from ~0.015 mm to ~0.04 mm. The displacement profile suggests non-uniform deformation along the radial length, which is influenced by local thermal and mechanical loading distributions [32]. The displacement magnitudes remain within mechanical safety limits, providing critical input for assessing fuel-cladding mechanical interaction and ensuring structural integrity under elevated reactor-power conditions.

At t = 3 s the fuel rod exhibits an annular von Mises stress distribution with maxima at the outer surface and near-uniform stress along the active height (Figure 5a). Increasing volumetric heat generation from 119 to 245 MW m⁻³ raises centerline stress from ~0.6 to ~1.5 GPa (Figure 5b) and outer surface stress from ~0.4 to ~1.0 GPa (Figure 5c). These results confirm that differential thermal expansion between the hot core and cooler cladding drives the dominant mechanical loading. The stresses remain within structural safety margins but provide critical input for evaluating fuel–cladding mechanical interaction and long-term integrity of the reactor core under high-power operation.

Figure 6 shows thermo-mechanical simulation results of clad material, where higher internal heat generation (119 → 163 → 245 MW m⁻³) causes the cladding to expand almost uniformly along its 0.5 m length while remaining safely elastic. Axial displacement grows steadily from about 0.012 mm to 0.03 mm with increasing heat source, with only small end-effects. However, average strain stays close to zero (≈10⁻⁴), rising slightly with power and showing minor local oscillations, confirming gentle, homogeneous volumetric growth. Mean stress rises from roughly 0.15 GPa to 0.35–0.4 GPa, with small peaks near the ends, remaining well below the yield strength of typical zirconium alloys.

Changes in the velocity profile are evident in Figure 7a, over 3 s of operation, the coolant maintains a stable upward flow through the cladding channel. Surface velocities fall from ~1 m s⁻¹ at the inlet to ~0.85 m s⁻¹ mid-height, while centerline velocities decrease smoothly to ~0.75 m s⁻¹ at the outlet. Increasing the internal heat source from 119 to 245 MW m⁻³ produces only slight reductions in surface velocity and negligible changes along the centerline. However, no change in pressure is observed over the 3 s period, even with the increase in heat sources. These results confirm that the hydrodynamic field remains stable and well within design limits, ensuring reliable convective heat removal under the investigated power densities.

4. Conclusions

Thermo-mechanical and hydrodynamic simulations of a thorium-based fuel rod subjected to transient volumetric heating (119 → 245 MW m⁻³ over 3 s) reveal stable structural and coolant responses under elevated power densities. The fuel experiences a characteristic high-strain core surrounded by a low-strain rim, with centerline volumetric strain increasing nearly linearly from 0.006 to 0.014, while remaining axially uniform along the heated 0.5 m segment. Total displacement grows from ~0.004 mm to ~0.012 mm at the surface and from ~0.015 mm to ~0.04 mm along the axis, confirming that radial constraint by the cladding limits central motion. Mechanical loading is dominated by differential thermal expansion, producing annular von Mises stress fields that rise from 0.6 → 1.5 GPa at the center and 0.4 → 1.0 GPa at the surface, yet remain comfortably below zirconium-alloy yield strengths. The cladding itself expands almost uniformly, with axial displacement increasing from ~0.012 mm to ~0.03 mm, average strain near 10⁻⁴, and mean stress limited to 0.35–0.4 GPa, all well within elastic and design margins. Coolant hydrodynamics remain robust throughout the transient. Upward flow is stable, with only small velocity reductions (surface: 1.0 → 0.85 m s⁻¹, centerline: 0.75 m s⁻¹ at outlet). Collectively, these results confirm that the combined fuel–cladding system maintains mechanical integrity and coolant stability under rapid power ramps to at least 245 MW m⁻³, providing critical design data for next-generation heavy-water reactors employing thorium-based fuels.