Neutronic and safety evaluation of a hybrid UO 2 - U 3 Si 2 accident tolerant fuel assembly for VVER-1000 reactors

1 Introduction

In the continuous evolution of Generation III + Pressurized Water Reactors (PWRs) and VVERs, the development of Accident Tolerant Fuels (ATF) has emerged as a primary trajectory to simultaneously improve the economics of normal operation and enhance coping times during severe design basis accidents (Zinkle et al., 2014; IAEA, 2004). Following the paradigm shift in nuclear safety protocols, extensive research has been directed toward transitioning away from standard Uranium Dioxide (UO 2) and zirconium alloy systems toward robust alternatives that mitigate runaway oxidation and hydrogen generation during Loss of Coolant Accidents (LOCAs) (Carmack et al., 2013). Among these alternatives, high-density fuel matrices, particularly Uranium Silicide (U 3 Si 2), have gained significant attention. U 3 Si 2 offers a heavy metal uranium density of 11.3 gU/c m 3—approximately 12 % higher than standard UO 2—providing a direct mechanism for power uprates and cycle length extension (White et al., 2015). Furthermore, U 3 Si 2 transitions away from low-conductivity ceramic heat transfer mechanisms. It offers metallic thermal conductivities exceeding 15.0 W/m⋅K, fundamentally altering radial heat transfer dynamics and drastically reducing the sensible stored energy within the core (Xiong et al., 2024; Fink, 2000).

Despite these profound thermal-hydraulic advantages, the integration of U 3 Si 2 presents significant reactor physics challenges. Recent computational and experimental studies have extensively evaluated the neutronic penalties associated with full-core ATF implementation or homogeneously dispersed advanced fuels (Xiong et al., 2024; Mohsen et al., 2021). Because high-density matrices act as aggressive thermal neutron absorbers, the complete replacement of a standard assembly lattice with U 3 Si 2 drastically alters the neutronic spectrum and power distribution (Zinkle et al., 2014). This excessive absorption leads to severe spectrum hardening, which subsequently degrades the integral worth of safety control rods and complicates maintaining a deeply negative moderator temperature coefficient (Lewis, 2008; Stammler and Abbate, 1983).

While previous literature has thoroughly investigated full-core U 3 Si 2 implementations and their associated global penalties, the commercial nuclear industry has recently pivoted toward “hybrid” or partial ATF loading schemes as a pragmatic transitional strategy (Terrani, 2018). Recent operational milestones, such as the 2026 pilot loading of heterogeneous hybrid ATF assemblies in the Balakovo VVER-1000 (Rosatom, 2026), indicate that strategically mixing standard UO 2 with advanced fuel matrices within a single assembly can optimize the neutronic economy while mitigating the severe regulatory and spectral penalties of a sudden full-core transition. However, despite this growing interest in partial loading, no existing study has evaluated the feasibility of a highly localized, high-density silicide hybrid radial configuration specifically tailored for the neutronic environment of a 3D VVER-1000 geometry. This represents the core novelty of this work.

To mitigate global core spectral effects while retaining the economic and thermal-hydraulic benefits of U 3 Si 2, this paper proposes a novel hybrid 3D VVER-1000 assembly. By strategically restricting the U 3 Si 2 pins to the outer periphery, the design capitalizes on the highly thermalized inter-assembly water gap neutron flux while maintaining a standard UO 2 driving inner core. This paper systematically evaluates the multi-physics viability of this hybrid 3D geometry. Continuous-energy Monte Carlo simulations using OpenMC are utilized to quantify isotopic breeding, spatial self-shielding, safety feedback parameters, and kinetic penalties (Romano et al., 2015). Concurrently, OpenFOAM Computational Fluid Dynamics (CFD) is employed to evaluate conjugate heat transfer, proving that the ATF’s superior thermal diffusivity safely absorbs the extreme neutronic power peaking and drastically enhances Loss of Coolant Accident (LOCA) survivability.

2 Methodology

Neutronic assessments were undertaken using the OpenMC continuous-energy Monte Carlo transport code (Romano et al., 2015) with the ENDF/B-VIII.0 nuclear data library (Brown et al., 2018). The geometric model was developed as a finite 3D assembly to adequately compute axial neutron leakage, true 3D spatial distributions of power, and reflector savings. To ensure statistical reliability, the Monte Carlo transport solver was executed with 30,000 particles per batch for a total of 120 batches (with 20 inactive batches discarded for fission source convergence).

2.1 Geometric model and hybrid assembly layout

The proposed hybrid ATF arrangement is based on the structural platform of a standard VVER-1000 assembly (Figure 1). An x-oriented hexagonal lattice defines the spatial domain. The model represents the critical reactor physics parameters intrinsic to VVER designs: exactly 312 fuel pins, 18 guide tubes, 1 central instrumentation channel, a nominal pin pitch of 1.275 cm, and the standard Hydrogen-to-Heavy-Metal (H/HM) volume ratio (IAEA, 2004). To model an infinite radial core array, the exterior of the assembly is defined by a reflective hexagonal prism with an edge length of 13.62 cm. To transition to a finite 3D space, the active fuel length was extruded explicitly to 353.0 cm (covering the axial domain between Z = − 176.5 cm and Z = + 176.5 cm). To simulate axial neutron leakage and correctly capture flux flattening at the core extremities, 20.0 cm light water reflectors were attached directly above and below the active fuel region. Vacuum boundary conditions were then applied to the outer surfaces of these top and bottom reflector regions to mathematically represent the physical boundaries of the reactor system.

The hybrid assembly consists of a total of 331 lattice positions, distributed into 11 concentric hexagonal rings. The pins are modeled geometrically with the following radial dimensions:

• Fuel Pins: Fuel pellet radius of 0.386 cm with an inner helium gap (outer radius of 0.400 cm), surrounded by a zirconium alloy cladding with an outer radius of 0.455 cm.
• Guide/Instrumentation Tubes: Moderator-filled inner region with a radius of 0.550 cm, enclosed by zirconium cladding with an outer radius of 0.630 cm.

The hybrid loading pattern strictly restricts the ATF to the outermost periphery:

• Standard UO 2 Pins (252 total): Standard UO 2 pins ( ρ = 10.5 g/cm 3 ) at 3.3 % enrichment occupy the inner regions of the core (Rings 1 through 9).
• High-Density U 3 Si 2 Pins (60 total): The high-density U 3 Si 2 ATF matrix ( ρ = 12.2 g/cm 3 ) at 3.3 % enrichment populates the 60 positions of the outermost ring (Ring 10) to make use of the highly thermalized inter-assembly water gap neutron flux.
• Control & Instrumentation Tubes (19 total): 18 control rod guide tubes are distributed symmetrically across Rings 3, 6, and 9, with 1 central instrumentation tube placed at the origin (Ring 0).

To explicitly represent the axial power distribution and extract the Axial Peaking Factor ( K z ) for subsequent thermal-hydraulic analyses, a high-resolution 1D regular mesh tally consisting of 50 equidistant axial bins was overlaid exactly across the active fuel length. This mesh filter was configured to score recoverable fission energy (kappa-fission), enabling the precise resolution of the chopped cosine flux profile.

2.2 Computational Fluid Dynamics and meshing methodology

In order to determine the thermal-hydraulic safety margins that were dictated by the 3D neutronic power distribution, a matching steady-state CFD model was created. The VVER-1000 fuel assembly thermal-hydraulic analysis was carried out with the use of OpenFOAM, which is an open-source Computational Fluid Dynamics (CFD) framework (Weller et al., 1998; OpenCFD Ltd, 2023). The steady-state conjugate heat transfer solver chtMultiRegionSimpleFoam was used to describe the salient hydraulic nature of heat transfer between the solid fuel, the Zirconium cladding, and the fluid coolant. This solver couples the Navier-Stokes equations with the energy equations in the fluid region, to the heat conduction equations in the solid regions, and solves until the temperatures at the physical interfaces are at thermal equilibrium.

2.2.1 Geometry and mesh generation

A hexagonal lattice 3D computational domain model of a VVER-1000 fuel assembly has been used (Figure 2). Both the fluid and solid fuel domains were explicitly modeled to the full active fuel length of Z = 3.53 m. To accurately map the OpenMC spatial tallies to the OpenFOAM mesh, the coordinate system was centered at the core midplane, defining the active inlet at Z = − 1.765 m and the active outlet at Z = + 1.765 m. Simulating the full, un-truncated 3.53 m geometry ensures that the axial heat conduction and the absolute peak cladding temperatures are captured with complete mathematical fidelity. The pin geometry was specified with an outer radius of the fuel pellet of r fuel = 3.78 mm and an outer radius of the cladding of r clad = 4.55 mm.

The multi-region mesh has been created with a strictly Cartesian, castellated method through the OpenFOAM native blockMesh and topoSet programs with the primary goal of achieving computational robustness and orthogonality over the large pin array. To begin with, a single structured hexahedrally-defined background mesh was created using blockMesh. The X and Y dimensions were bound to − 0.165 m and 0.165 m respectively, which contained the full assembly, to produce a detailed computational grid of about 9,000,000 hexahedral cells.

Instead of using body-fitted snapping algorithms, topoSet was used to cut up regional cell zones directly out of the background mesh. The cylinderToCell topological source was used to test the spatial coordinates of the centers of cells. The cells that were within the physical coordinates of the standard pins were assigned to the UO2 cellSet, with certain peripheral positions that would interface with the inter-assembly water gap assigned to the U3Si2 cellSet. The same geometric inclusion logic was used to define the cladding sets by the 4.55 mm outer radius in order to get a stair-stepped approximation of the cylindrical geometries. Lastly, the program splitMeshRegions -cellZones -overwrite was run to physically divide the single background mesh into isolated, interacting mathematical entities (e.g., fuel, cladding, and coolant). The regionProperties dictionary was then adjusted such that the fluid and solid regions were connected in order that the solver could apply isolated thermophysical properties and governing equations to each discrete zone.

2.2.2 Boundary conditions and source terms

At the bottom inlet boundary ( Z = − 1.765 m ), coolant water was injected at a realistic VVER-1000 baseline operational temperature of T i n = 563 K (IAEA, 2004). The external boundaries of the hexagonal assembly were subjected to symmetry/periodic boundary conditions to model an infinite lattice environment. Any axial boundary conditions that the solid fuel pins (inlet and outlet patches) had with the rest of the system were rigorously defined to be zeroGradient (adiabatic) in order to correctly isolate the fuel volume. This type of boundary treatment is essential for the proper modeling of intrinsic axial heat conduction as a primary mechanism of heat removal in the highly conducting U 3 Si 2 pins.

To simplify the computational demands of meshing a microscopic volume of fluid, the helium gap was modeled explicitly as a continuous layer of 1D thermal resistance utilizing the compressible::turbulentTemperatureRadCoupledMixed boundary condition at the fuel-to-cladding physical interface. This allowed the precise fabrication of a 0.14 mm thick layer having a constant thermal conductivity of k = 0.25 W/m - K, successfully simulating the heat transfer penalty of the helium backfill gas (Todreas and Kazimi, 2012).

OpenFOAM was used to prescribe the generation of heat in the fuel pins through the fvOptions framework. A chopped-cosine axial power shape was used instead of a uniform profile. As this study prioritized evaluating the bounding thermal safety limits of a static 3D power distribution, a strict one-way multi-physics coupling scheme was employed, utilizing no Picard iterations for dynamic thermal feedback (Ivanov and Avramova, 2007). The power-to-mesh mapping was executed via a custom Python script that parsed the high-resolution OpenMC kappa-fission Cartesian mesh tally. The script geometrically mapped the neutronics grid data directly to the OpenFOAM body-fitted cellSet centroids, generating discrete spatial scalarSemiImplicitSource terms. This one-way extraction of the exact volumetric heat source distributions—including the radial peaking factor ( F x y ) and the axial peaking factor ( K z )—aligns with recent multi-physics coupling benchmarks for advanced fuel assessments (Baidoo et al., 2024). Thermophysical properties of the standard UO 2 fuel were modeled using temperature-dependent thermal conductivity ( k ) to adequately model the degradation of heat transfer in the centerline at high temperatures, whereas the U 3 Si 2 Accident Tolerant Fuel (ATF) was given properties that reflect its much smoother, enhanced metallic thermal conductivity profile (White et al., 2015).

2.3 Reactivity feedback coefficients

To ensure the hybrid assembly retains passive safety capabilities, continuous-energy stochastic sweeps were performed to extract the Fuel Temperature Coefficient (FTC), Moderator Temperature Coefficient (MTC), and Control Rod Worth (CRW) (Lewis, 2008; Glasstone and Sesonske, 1994). The generalized reactivity coefficient α_T is calculated using Equation 1:

α T = k 2 − k 1 k 1 k 2 T 2 − T 1 × 1 0 5 p c m / K (1)

For the MTC evaluations requiring cross-sections at intermediate temperatures not explicitly provided in the evaluated library (e.g., 550 K, 580 K), OpenMC was configured to utilize stochastic mathematical interpolation between bounding datasets. Because the updated geometric model incorporates finite axial boundaries, Control Rod Worth was evaluated by comparing the effective system eigenvalue ( k eff ) between the All Rods Out (ARO) and All Rods In (ARI) states, inherently accounting for axial neutron leakage. To accurately simulate the mechanics of the VVER-1000 Control Rod Cluster (CRC), the central instrumentation channel was locked as a continuous liquid moderator path. The 18 surrounding guide tubes were toggled between moderator (ARO) and enriched B 4C ( ρ = 2.52 g/cm 3 ) absorber pins (ARI). The integral worth is calculated using Equation 2:

C R W = k ARO − k ARI k ARO ⋅ k ARI × 1 0 5 p c m (2)

3 Results and discussion

3.1 Reactivity evolution and fission product poisoning

The overall goal of the hybrid assembly is to sustain a long 540-day full-power cycle of criticality. As the finite 3D geometry is considered, axial neutron leakage intrinsically reduces the absolute reactivity of the system in comparison to an infinite lattice. Figure 3 shows the evolution of the effective multiplication factor ( k eff ) as a function of depletion time and burnup. At the Beginning of Cycle (BOC), the fresh hybrid assembly is in a highly reactive state ( k eff ≈ 1.155 ). However, during the initial 15 days of operation, there is a sharp, non-linear reactivity penalty, which decreases rapidly during that period to an effective value of approximately k eff ≈ 1.115.

This initial suppressive reactivity is wholly controlled by the rapid accumulation of large absorption cross-section fission products (Glasstone and Sesonske, 1994). The isotope ¹³⁵Xe accumulates very quickly, driven by both its primary direct fission yield and the beta decay of ¹³⁵I, reaching its steady-state equilibrium concentration in approximately 72 hours. Concurrently, the accumulation of ¹⁴⁹Sm (through the beta decay chain of ¹⁴⁹Nd and 149Pm) adds a second reactivity penalty that becomes nearly saturated during the initial 300 days of the operating cycle (Lewis, 2008). This isotopic behavior is explicitly captured and visualized in Figure 4, which plots the total mass accumulation of these specific neutron poisons within the inner UO 2 core. The steep asymptotic rise of the Xenon and Samarium curves directly correlates with, and explains, the sharp initial drop in the k eff trajectory observed during the first month of operation.

3.2 Breeding plutonium and minor actinides

The fundamental mechanism enabling the cycle extension is the strategic placement and heavy metal loading of the U 3 Si 2 matrix. Because U 3 Si 2 possesses a density of 11.3 gU/cm 3—approximately 12 % higher than standard UO 2—it provides a significantly larger macroscopic fertile base of ²³⁸U atoms per unit volume (White et al., 2015). Confining this high-density matrix exclusively to a single, outermost peripheral row (Ring 10) is a deliberate neutronic design choice. The periphery is exposed to the highly thermalized neutron flux scattering back from the bulk inter-assembly water gap. The dense silicide pins act as an aggressive thermal neutron sink, maximizing localized radiative captures and rapidly accelerating the breeding of fissile ²³⁹Pu. By transitioning to the 3D finite geometry, the isotopic inventories physically integrate the full 3.53 m active core height. As shown in Figure 5, this peripheral configuration breeds ²³⁹Pu per full-length pin at a significantly higher rate than the inner UO 2 core. The sustained in situ burning of this newly bred Plutonium supplies the critical additional neutron source needed to offset the depletion of the original ²³⁵U inventory, thereby achieving the 540-day extended working life. Furthermore, this single-row arrangement serves a vital safety function: the “black” peripheral absorber ring successfully shields the inner core from severe spectral hardening, ensuring the central standard UO 2 pins retain optimal moderation for effective control rod insertion and stable temperature feedback margins.

Although this accelerated rate of Plutonium breeding provides a strict neutronic benefit in extending the cycle, it also inherently alters the minor actinide inventory at End of Cycle (EOC). The total 3D generation of the highly radiotoxic minor actinide of interest, ²⁴¹Am, is mapped in Figure 6. This inventory dictates spent nuclear fuel management and decay heat shielding requirements (Salvatores, 2005). The generation of ²⁴¹Am is a direct consequence of the Plutonium breeding chain; sequential radiative captures ultimately result in its formation via the reaction sequence shown in Equation 3:Because the peripheral ring breeds more parent Plutonium, it predictably exhibits a small, yet significant, increase in the total mass of ²⁴¹Am per pin by EOC (Glasstone and Sesonske, 1994). This implies that while the hybrid assembly is highly efficient at generating power, the long-term radiotoxicity profile of the discharged assembly will be dominated by the metallic periphery.

P u 239 → n , γ 240 P u → n , γ 241 P u → β − 241 A m (3)

To mitigate the radiotoxic burden of this elevated ²⁴¹Am inventory, the deployment of this hybrid assembly must be coupled with advanced backend fuel cycle strategies. The most effective mitigation strategy is Partitioning and Transmutation (P&T) (Salvatores, 2005). By reprocessing the discharged U 3 Si 2 pins, the separated minor actinides can be fabricated into specialized transmutation targets and consumed within the fast neutron spectrum of Generation IV sodium- or lead-cooled fast reactors. Alternatively, if a once-through cycle is maintained, the thermal limits of deep geological repository canisters must be re-evaluated to account for the increased decay heat, potentially requiring extended interim wet-pool cooling periods prior to final dry-cask loading.

3.3 Fissile inventory shift and macroscopic poison burden

The extended reactivity lifespan is directly related to the exceptional neutronic economy of the outer U 3 Si 2 pins. The primary mechanism sustaining the cycle is the fundamental shift in the macroscopic fissile inventory of the assembly. As time progresses throughout the operational cycle, the mass of the primary fissile isotope, ²³⁵U, decreases continuously. This depletion is a faithful representation of the volume-integrated axial depletion profile of the 3D core, as shown in Figure 7. Simultaneously, the metallic ATF matrix drives the continuous ²³⁸U(n, γ) radiative capture reaction, breeding large amounts of fissile ²³⁹Pu. This dynamic inventory crossover ensures that as the ²³⁵U reaction rate drops, the in situ burning of the bred ²³⁹Pu supplies an additional source of neutrons, leveling the depletion curve.

In addition to the short-lived equilibrium poisons (¹³⁵Xe and ¹⁴⁹Sm), the assembly develops a steady accumulation of long-lived macroscopic fission products. The cumulative mass generation of primary long-lived neutron absorbers (e.g., ¹⁴³Nd, ¹⁴⁷Pm, ¹³³Cs, and ¹⁰⁵Rh) is illustrated in Figure 8. The proportionately greater absolute poison burden observed in the high-density periphery is directly associated with the higher local fission density in that region. The strict linearity of this accumulation validates that, within the 3D spatial domain, fuel burnup occurs in a continuous, stable manner without erratic local spectral shifts.

3.4 Spatial self-shielding and isotopic breeding

The assumption of uniform burnup across the peripheral ATF matrix is physically invalid due to severe spatial resonance self-shielding (Stammler and Abbate, 1983). In the 3D finite geometry, a static reaction rate tally was performed to map the total volume-integrated ²³⁸U(n,γ) captures across the full 353 cm active length of the 60 peripheral pins at Beginning of Cycle (BOC). The findings reveal a profound azimuthal breeding gradient, as illustrated in Figure 9. The maximum ²³⁹Pu breeding rates occur at the flat face centers of the hexagonal geometry, where the pins possess a continuous solid angle facing the highly thermalized bulk inter-assembly water gap. At these face centers, the maximum capture rate peaks at 5.00037 × 1 0 − 3 reactions per source neutron. Conversely, the minimum capture rates are observed at the hexagonal vertices (corner pins), where the geometry pinches and the local moderator-to-fuel ratio drops, causing a localized spectrum hardening. The capture rate at these corner pins drops to 4.51260 × 1 0 − 3 reactions per source neutron. This represents a significant spatial gradient variance of 10.81 % across a single assembly ring.

3.5 Radial power peaking and the 3D Hot Channel Factor

The primary physical consequence of the peripheral thermalization and accelerated burnup is localized azimuthal power disparity. Because the metallic ATF ring acts as an aggressive thermal neutron sink, it drives a higher volumetric fission rate than the inner UO 2 core. Figure 10 visualizes the volume-integrated radial pin power distribution across the hybrid assembly. The extraction of the OpenMC statepoint data yields a total assembly power of 18.40 MWth. The localized peak volumetric power generation was recorded as 399.18 W/cm 3 against an assembly average of 372.18 W/cm 3, strictly yielding a Radial Peaking Factor ( F x y ) of 1.073. For a highly heterogeneous assembly configuration, a radial peaking factor of 1.073 indicates a remarkably well-balanced macroscopic flux distribution. To properly bound the thermal-hydraulic safety analysis, this radial component was multiplied by the chopped-cosine Axial Peaking Factor ( K z = 1.8858 ) extracted from the 3D mesh tally, mathematically resulting in a total 3D Hot Channel Factor ( F q ) of 2.023. The total multi-physics peaking factor is calculated using Equation 4:

F q = F x y × K z = 1.073 × 1.8858 = 2.023 (4)

This F q value of 2.023 dictates the absolute maximum localized volumetric heat generation ( q ‴ ) that must be dissipated by the primary coolant in order to maintain operational safety margins.

3.6 Neutron flux lethargy and spectrum hardening

To definitively quantify the spectral impact of the high-density fuel loading within the 3D finite geometry, the volume-integrated neutron flux per unit lethargy ( ϕ u ) was tallied independently for the inner UO 2 core and the peripheral U 3 Si 2 matrix across the full 3.53 m active fuel length. The explicit physics of spectrum hardening are revealed in Figure 11. The silicide matrix exhibits a visibly depressed thermal flux peak (in the 1 0 − 2 to 1 0 − 1 eV domain) compared to the standard UO 2 core. This confirms that the dense packing of Uranium atoms in the metallic fuel acts as a “black” absorber to thermalized neutrons entering from the inter-assembly water gap, successfully shielding the inner core but shifting the localized peripheral spectrum toward the epithermal and fast energy regimes.

3.7 Safety validation: reactivity feedback coefficients

To ensure the hybrid assembly retains passive safety capabilities, continuous-energy stochastic sweeps were performed to extract the fundamental temperature coefficients.

Doppler Defect (FTC): Sweeping the fuel temperature from 600 K to 1200 K while holding the moderator fixed yielded an FTC of − 1.895 pcm/K (Figure 12). It must be noted that this apparently diminished magnitude is a mathematical artifact of evaluating the defect in a highly reactive, unpoisoned single-assembly lattice at BOC ( k ∞ ≈ 1.33 ). The formulation of α T scales inversely with the k 1 k 2 denominator ( ≈ 1.75 ). If normalized to a critical whole-core state ( k ≈ 1.0 ), the hybrid assembly FTC scales to approximately − 3.31 pcm/K, perfectly aligning with standard VVER-1000 Doppler feedback expectations (Lewis, 2008). This confirms that the resonance broadening of ²³⁸U in both the oxide and silicide matrices provides prompt, negative reactivity suppression.

Moderator Expansion (MTC): Sweeping the moderator temperature from 550 K to 610 K, with corresponding density drops from 0.760 g/cm 3 to 0.630 g/cm 3 (incorporating an intermediate evaluation at 580 K and 0.710 g/cm 3) at an operating pressure of 15.5 MPa, yielded an MTC of − 36.533 pcm/K (Figure 13). The highly negative magnitude guarantees robust neutronic choking in response to coolant-driven transients, as the loss of moderation efficiently hardens the neutron spectrum into the ²³⁸U capture resonances (Glasstone and Sesonske, 1994).

3.8 Control Rod Worth and spectrum efficacy

The spectrum hardening visualized in the prior lethargy evaluation raises a critical safety concern: the potential degradation of 1 / v thermal absorbers such as boron carbide (B 4C). To ensure that the hardened spectrum does not overwhelm the available shutdown margin, the total Control Rod Worth (CRW) was evaluated. To accurately simulate the mechanical reality of a VVER-1000 Control Rod Cluster (CRC), the central instrumentation channel was locked as a continuous liquid moderator path. The 18 surrounding guide tubes were toggled between moderator (All Rods Out, ARO) and enriched B 4C ( ρ = 2.52 g/cm 3 ) absorber pins (All Rods In, ARI). The resulting eigenvalue changes are summarized in Table 1. The magnitude of the reactivity drop ( 65,323.56 pcm ) does not represent a standard operational shutdown margin. Rather, because the OpenMC geometry utilizes reflective radial boundary conditions, this value represents a theoretical maximum absolute bounding state: an infinite core array where every single adjacent assembly drops its control cluster simultaneously. This confirms that despite the localized thermal flux depression in the periphery, the 18-rod clusters maintain massive, overwhelming shutdown efficacy.

TABLE 1. Control Rod Worth evaluation for an 18-pin B 4C cluster in a 3D finite lattice at BOC.

3.9 Time-evolution of the axial peaking factor

The axial power profile within the VVER-1000 assembly is not a static parameter. As the operational cycle progresses, the spatial distribution of the neutron flux is continuously reshaped by isotopic depletion and the localized accumulation of macroscopic neutron absorbers. The time-evolution of the Axial Peaking Factor ( K z ) across the full 540-day depletion cycle is mapped in Figure 14.

The assembly exhibits a continuous, stable axial flux flattening phenomenon. At the Beginning of Cycle (BOC), the uniform fresh fuel loading dictates a strict chopped cosine flux distribution, resulting in a maximum K z of 1.8858. However, because the localized fission rate is highest at the core midplane, this region experiences the most aggressive differential burnup. During the first 30 days of full-power operation, the curve exhibits a steep, non-linear decline. This initial flattening is primarily governed by the rapid, localized buildup of equilibrium fission product poisons, specifically ¹³⁵Xe and ¹⁴⁹Sm. Because these isotopes possess massive thermal neutron absorption cross-sections, their concentrated accumulation at the midplane heavily suppresses the central thermal flux peak. Following the saturation of these short-lived poisons, the K z trajectory transitions into a steady, linear decline through to End of Cycle (EOC). This extended phase is driven by macroscopic differential burnup. The primary fissile inventory (²³⁵U) is consumed more rapidly at the core center than at the axial extremities. Although the metallic periphery actively breeds ²³⁹Pu, the net accumulation of long-lived parasitic fission products continuously forces the neutron flux to migrate axially outward toward the less-depleted regions of the fuel rod. By day 540, the Axial Peaking Factor reduces to 1.8539, representing a total absolute peak reduction of Δ K z = 0.0319. As the Zirconium cladding accumulates radiation damage over the prolonged cycle, the maximum localized linear heat generation rate it must endure is systematically relaxed, preserving structural integrity.

4 Thermal-hydraulic safety analysis

Applying the bounding values 3D Hot Channel Factor ( F q = 2.023 ) obtained from the OpenMC tallies (Section 3.5), the steady-state OpenFOAM model of conjugate heat transfer was run in order to evaluate the safety limits in terms of thermal-hydraulic conditions, The chief justification for the deployment of Accident Tolerant Fuels (ATF) is based on the extinction of the maximum fuel centerline temperature ( T C L ). The axial centerline temperature profiles for both the standard UO 2 inner pins and the peripheral U 3 Si 2 ATF pins, evaluated from the coolant inlet ( Z = − 1.765 m ) up to the active fuel outlet ( Z = + 1.765 m ).

4.1 Axial centerline temperature and thermal safety margins

Applying the bounding 3D Hot Channel Factor ( F q = 2.023 ) derived from the OpenMC tallies (Section 3.5), the steady-state OpenFOAM conjugate heat transfer model was executed to evaluate the thermal-hydraulic safety margins. The primary justification for the deployment of Accident Tolerant Fuels (ATF) relies on the suppression of the maximum fuel centerline temperature ( T C L ). Figure 15 illustrates the axial centerline temperature profiles for both the standard UO 2 inner pins and the peripheral U 3 Si 2 ATF pins, evaluated from the coolant inlet ( Z = − 1.765 m ) up to the active fuel outlet ( Z = + 1.765 m ).

The CFD analysis shows a large decrease in peak fuel temperatures. The standard UO 2 pin reaches a maximum centerline temperature ( T C L ) of approximately 740 K at the midplane elevation. In comparison, U 3 Si 2 pin has a maximum T C L of about 690 K, absolute thermal margin expansion of 50 K, at the hottest, localised spot in the fuel assembly. This thermal suppression is directly determined by the most basic differences in thermophysical properties between the two materials. Standard UO 2 is a ceramic oxide of operatively low thermal conductivity, which further decreases with increasing temperatures due to phonon (Umklapp) scattering (Fink, 2000; Todreas and Kazimi, 2012), This thermal resistance making sensible heat to be bottleneck in radial center of the pellet, some steep parabolic temperature gradient. Conversely, U 3 Si 2 shows high thermal diffusivity in order to quickly and efficiently transfer heat from the pellets outward across the helium gap and into the Zirconium cladding (White et al., 2015); Furthermore, the crossing of the temperature profiles in the lower half of the core ( Z < 0.0 m ) highlights a secondary, but critical, mechanism of heat transfer: macroscopic axial conduction. At the core inlet ( Z = − 1.765 m ), the standard UO 2 centerline temperature is tightly coupled to the 563 K coolant, starting near 570 K. Because UO 2 is a thermal insulator, it is incapable of efficiently conducting heat down the length of the rod. Consequently, the heat generated at the core midplane remains trapped there, driving the peak to 740 K. The U 3 Si 2 profile, however, initializes at a much higher temperature of approximately 590 K at the inlet. This is not an anomaly, but rather a direct manifestation of its superior metallic diffusivity. The silicide fuel essentially acts as an axial heat pipe, rapidly conducting thermal energy away from the high-power midplane and forcing it down toward the cooler extremities of the rod. By efficiently redistributing these severe thermal loads both radially and axially, the U 3 Si 2 pin eliminates the concentrated midplane hot-spot that structurally stresses standard oxide fuels. This 50 K reduction in the peak T C L translates directly to a massive reduction in localized sensible stored energy. During a postulated Loss of Coolant Accident (LOCA), this reduced initial energy state strictly limits the severity of the subsequent temperature transient, fundamentally increasing the reactor’s margin to melting and providing invaluable extended operator coping time (Zinkle et al., 2014).

4.2 Radial temperature gradient and margin to melting

While the absolute maximum temperatures in the core are controlled by the axial temperature profiles, internal heat transfer kinetics and localized margin to fuel melting are dictated by the radial distribution of temperature in the fuel pellet. An example of this is Figure 16, which shows the radial temperature gradient from the center line of pellet ( r = 0.0 mm ) to outer surface of pellet ( r = 3.78 mm ), (extracted at the bounding axial midplane elevation of Z = 0.0 m). This steady-state CFD And without delving into the chemical details of the system, the actual analysis should demonstrate steep contrast in the localized heat transfer mechanics in the two fuel types. The standard UO 2 pin has a steep parabolic temperature profile characteristic of low conductivity ceramic oxides, which creates a severe temperature gradient from the centerline to the edge (Δ T UO 2) of 69.4 K. In direct contact, the metallic ATF ring is an effective way to remove sensible heat rapidly toward the pellet periphery. This gives extremely flat radial temperature profile with a minimum gradient (Δ T U 3 Si 2) of only 12.7 K. This 81.7 % reduced the radial thermal gradient is an operational need for silicide fuels to survive. A known thermodynamic tradeoff of U 3 Si 2 is its lowest absolute melting point (around 1938 K) with respect to UO 2 standard (around 3120 K) (White et al., 2015). However, the results from the CFD calculation show that the highly conductive nature of the metallic lattice effectively neutralises this drawback”. By keeping the concentration of heat at the centerline, the U 3 Si 2 pellet has a near-isothermal radial profile. Consequentially, even under bounding 3D hot channel conditions Flatt radial gradient ATF centerline is more than 1200 K, melting is below melting infrastructure. Additionally, this low variation in temperatures suppresses temperature-dependent phenomena, e.g., localized fission gas diffusion, which ensures enhanced In literature reported thermal stability of the fuel matrix during the steady state operation (Carmack et al., 2013).

4.3 Azimuthal temperature gradient and thermal bowing

In addition to analysis of absolute maximum core temperatures, the multi-physics coupling analysis needs to account for the spatial location of the ATF. In this type of hybrid assembly design, the external fuel lattice is infiltrated by a highly asymmetric neutronic and thermal-hydraulic environment. The inter-assembly water gap is a localized thermal neutron trap surrounding the hexagonal assembly. This casts a temperature-localized moderator resulting in a sharp thermal current of neutrons peak, resulting in an asymmetric production of heat very strongly asymmetrical towards the gap-facing outer surface of the pins of the periphery. The endothermic response of the UO 2 in the fuel cells that is specifically designed for precipitation is set by the high degree of incoitance of the UO 2 cubes within the fuel element fabric. fuel cannot well redistribute the asymmetric load of this heat, leading to a severe azimuthal temperature gradient in the fuel rod. This graded thermal expansion causes thermal expansion over the diameter of the cladding, resulting in a process of thermomechanical degradation (thermal bowing which causes the rod to warp and bend (Olander, 1976). Causing extreme thermal bowing can severely colouring of the adjacent coolant subchanels, the local flow rates are slowed down and induce premature Departure from Nucleate Boiling [DNB] Point (Todreas and Kazimi, 2012). As shown in Figure 17, the diametrical temperature profile of peripheral ATF pin at the bounding axial midplane elevation ( Z = 0.0 m ), between the entire 9.1 mm diameter From the assembly-facing inside surface to the gap-facing outside surface Asymmetric penetrance The coupled CFD results show that the asymmetric neutronic penalty is completely counteracted by this very high conductivity of the silicide matrix. Instead of localized-pressing the wildfire, trapping the wildfire. power spike at the outside edge, the U 3 Si 2 matrix vigorously passes the thermal energy through the pellet diameter. Consequently, the (azimuthal) temperature tilt ( Δ T tilt ) across the opposite cladding surfaces is restricted to a minimum of microscopic 4.08 K (602.70 K on the gap side compared to 598.62 K, on the assembly side). That near isothermal azimuthal profile guarantees that the cladding uniformly expands radially. By essentially cancelling the forces for divergent thermal expansion the ATF memorizes the latitudinal structure of the fuel rod, to prevent subchannel restriction and maintain critical thermal-hydraulic safety margins throughout the core.

4.4 Coolant temperature and thermodynamic conservation

Although it is the solid state thermal analysis which proves the increase in safety margins within the fuel and the cladding it is essential to prove the thermodynamic response of the surrounding fluid. Because of the effective heat transfer by the highly conductive U 3 Si 2 ATF to the cladding, far superior to the heat transfer by stan-dard UO 2, the local heat flux boundary condition at the fluid-solid interface is changed. Figure 18 shows the axial increase of the bulk coolant temperature in the subchannels close to the standard UO 2 and ATF pins. Two different aspects of thermodynamics can be distinguished from the CFD results: the kinetics of the heat flow process as well as the conservation of energy at macroscopic level. First, from the path of the curves, we can see that the coolant surrounding the U 3 Si 2 pin (blue) gets much heat up in the bottom segment of the core (Z = − 1.765 to − 0.965 m). Since the silicide metallic matrix does not suffer from the thermal resistance at this stage that is so harsh with ceramic oxides, it rejects the sensible heat to the fluid almost instantaneously upon creation. On the other hand, the less thermally-resistant UO 2 pin delays the radial cooling, in addition to carrying the bulk of the coolant temperature increase further up the axial length of the rod. Although these two subchannels have a different heat transfer kinetic process, the coolant exits the active fuel region at the same peak temperature of 568.0 K. This physical confluence is as direct manifestation of the First Law of Thermodynamics of the conjugate heat transfer model. Because the mass flow rate (18.40 MWth, as calculated in Section 3.5), fluid specific heat capacity ( C p ) and cumulative integrated neutronic power generation (18.40 MWth, as calculated in Section 3.5) are equally conserved within the channels, the total enthalpy increase (18.40 MWth, as calculated in Section 3.5) of the bulk fluid must be equal at last. Most importantly this analysis provides a vessel supporting the thermal-hydraulic sustainability of the ATF assembly under nominal operating conditions. Standard VVER-1000 reactors are designed for the primary system pressure of about 15.7 MPa with corresponding coolant saturation (boiling) temperature of about 618 K (IAEA, 2004). The maximum localized coolant temperature, measured adjacent to the ATF pins at the outlet, reached 568.0 K. Evaluated against the system saturation (boiling) temperature of 618 K at 15.7 MPa, this confirms the methodology successfully maintained a 50 K subcooled margin. This categorically validates the assumption of single-phase flow throughout the domain, ensuring the solid-state heat transfer profiles were isolated without interference from Departure from Nucleate Boiling (DNB) mechanics.

5 Conclusion

The proposed hybrid UO 2-U 3 Si 2 VVER-1000 assembly successfully demonstrates an optimized balance between neutronic economy and thermal-hydraulic safety. Through high-fidelity, one-way multi-physics coupling using OpenMC and OpenFOAM, the following key outcomes were established:

• Cycle Length Extension: The peripheral deployment of high-density U 3 Si 2 leverages inter-assembly thermalization to rapidly breed ²³⁹Pu. This dynamic fissile inventory shift sustains criticality ( k eff ≈ 1.069 ) across an extended 540-day operating cycle.
• Passive Safety Margins: Despite severe physical spectrum hardening in the high-density periphery, the reactor maintains strong passive safety feedback. The assembly exhibits deeply negative temperature coefficients (FTC: − 1.895 pcm/K; MTC: − 36.533 pcm/K) and retains a massive shutdown margin ( 65,323.56 pcm ) utilizing standard B 4C control clusters.
• Thermal-Hydraulic Superiority: The multi-physics analysis confirms that the severe localized radial power peaking ( 399.18 W/cm 3 ) induced by the ATF is completely offset by its high metallic thermal conductivity. The U 3 Si 2 pins act as axial heat pipes, conducting thermal energy away from the high-power midplane.
• LOCA Survivability: By operating at a significantly lower absolute peak centerline temperature ( 690 K ) compared to the driving UO 2 core ( 740 K ), the localized sensible stored energy is drastically reduced. This 50 K reduction limits the severity of temperature transients during a postulated Loss of Coolant Accident (LOCA), minimizing thermal shock on the cladding and maximizing operator coping time.

To further advance this hybrid concept toward commercial implementation, future work will focus on fully coupled transient 3D subchannel CFD simulations. Assessing the localized Departure from Nucleate Boiling Ratio (DNBR) during dynamic flow excursions and characterizing the LOCA blowdown phase will serve as the ultimate validation of the enhanced safety margins provided by this design.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Generative AI (Large Language Models) was used exclusively for English language refinement, grammar correction, and LaTeX formatting assistance during the preparation of this manuscript. Generative AI was not used in the conceptualization, execution of the OpenMC/OpenFOAM simulations, data analysis, or generation of scientific conclusions. The author has meticulously reviewed all text and takes full, sole responsibility for the manuscript’s scientific integrity and content.

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: Lingfeng He, North Carolina State University, United States

Reviewed by: Shuvendu Shivam, Indian Institute of Technology Jammu, India

Dennis Solomon Balami, University of Maiduguri, Nigeria