Analysis methodology

Refueling reactivity coefficient is expressed as follows:(1)where ρ is the reactivity, and Q is the number of refueling pebbles.

The calculation of refueling reactivity coefficient requires critical calculation before and after refueling. First, calculate the group constant of the core. Then calculate critically based on neutron diffusion equation to obtain effective multiplication factor k_eff before and after refueling. Next, calculate the ratio of and to obtain the refueling reactivity coefficient. The result from this method is the average refueling reactivity coefficient within the calculated refueling amount range. Calculation accuracy relates to the refueling amount - generally, smaller amounts yield higher accuracy, though very small amounts can affect calculation precision. Perturbation theory is more suitable for calculating refueling reactivity coefficient [30-33]. Previous work [29, 34] derived the theoretical expression of refueling reactivity coefficient. Perturbation theory analyzes small perturbations’ influence when there is no obvious distortion. From the one-group neutron diffusion equation, the reactivity introduced by the perturbation can be derived as(2)where V is the volume, ϕ is the neutron flux, ν is the neutron yield per fission, is the macroscopic fission cross-section, is the macroscopic absorption cross-section, and D is the diffusion coefficient.

Equation (2) presents the specific form of the perturbation equation derived from the one-group neutron diffusion equation. It shows that the contributions of the fission and absorption cross-sections to reactivity are weighted by the square of neutron flux, whereas the contribution of the diffusion coefficient is weighted by the squared flux gradient. Because the cross-section perturbations are expressed as differences, Eq. (2) is referred to as the difference formula. For small perturbations, the neutron flux generally remains undistorted; thus, according to Eq. (2), the reactivity change can be evaluated by simply determining the perturbations in the neutron flux and macroscopic cross-sections.

The expression αQ is derived based on the refueling characteristics of the PB-HTR. Assuming that the flow velocity of fuel pebbles in the core is v, travel distance in the period Δt is . The change in at the space position r is(3)where Ni is the atomic density of the nuclide i and is the microscopic fission cross-section of nuclide i.

Using the first-order Taylor expansion, we obtain:(4)The PB-HTR uses graphite as the reflector. The cladding and matrix of the fuel pebbles were composed of graphite. The absorption cross-section of graphite is relatively small. The diffusion length of the neutrons in graphite is long, and the distribution of the neutron flux in the core is relatively flat. Simultaneously, the fuel pebbles were distributed randomly and evenly inside the core. The control rods are arranged in a reflector outside the core. Therefore, the neutron flux and other parameters do not exhibit relatively large local peaks. Therefore, the PB-HTR exhibited better macroscopic uniformity in the reactor core. This is the theoretical reason for using Taylor expansion in Eq. (4) to preserve only first order term.

Some researchers have conducted experimental studies on pebble flow in PB-HTR. This result proved that our simplification was reasonable. In the study of Khane et al. [35-37], the pebble flow dynamics were studied in a scaled-down test reactor. Using a non-invasive radioactive particle tracking (RPT) technique used a cobalt-60 based tracer to mimic the shape, size, and density of pebbles. The results show a plug-type flow in the cylindrical region of the reactor at all seeding positions. In a plug-type flow, the flow velocity is assumed to be constant at any cross-section perpendicular to the pipe axis. There was no velocity gradient in the plug, and the entire flow state moved forward at the same speed. This study shows that the flow velocity can be considered uniform. In other words, the transient perturbations within the flow were sufficiently weak. In the study of Gatt [38] and Li et al. [39], the results indicated that the pebble motion paths were streamlined, and there was little interference or crossing between pebble trajectories. This also supports the simplification of Eq. (4).

Put the Eq. (4) into Eq. (3), and we can get(5)Similarly,(6)where is the microscopic absorption cross-section of nuclide i.

The diffusion coefficient D is mainly attributed to the graphite in the PB-HTR, and remained unchanged before and after refueling. So,(7)Put the Eq. (5) – (7) into Eq. (2), and we can get(8)Refueling speed q is defined asfollows:(9)Thus, the expression for αQ based on perturbation theory can be obtained as(10)During the operation of PB-HTR, the refueling changes the cross-section and other parameters, and causes a change in reactivity. Equation (10) can be used to directly solve the reactivity introduced by refueling according to the current state parameters of reactor, without critical calculations after refueling, which improved the calculation efficiency and ensured the accuracy of the calculation. In this study, perturbation theory was used to calculate the refueling reactivity coefficient, and the effects of the key parameters were studied.

The main innovation of this method is the first use of perturbation theory to analyze the continuous refueling of the PB-HTR. The fuel management scheme of the PB-HTR involves continuous online refueling. The fuel pebbles were loaded from the top of the reactor core, flowed to the bottom, and then unloaded. Positive reactivity can be replenished in the reactor core through continuous online refueling. The refueling speed of the equilibrium core in the HTR-PM was approximately 6,000 fuel pebbles per day. The average distance travelled by the fuel pebbles was approximately 16 cm per day, which is very small compared to the height of the reactor core (11 m). Therefore, the continuous refueling of the PB-HTR can be regarded as the continuous introduction of small disturbances. Hence, it is reasonable to use the perturbation theory to analyze the continuous refueling of the PB-HTR. The application of perturbation theory is simpler than that of prior methods, which only needs to calculate the reactivity of the reactor core status before refueling without calculating the status after refueling. Previous methods must calculate the reactivities before and after refueling.

Calculation model

HTR-PM demonstration power plant is a 2250 MWth thermal reactor based on two modules, designed by INET at Tsinghua University. The HTR-PM reactor uses inert helium as the coolant, with a helium temperature of 250 ℃ at reactor inlet and 750 ℃ at reactor outlet. The primary helium pressure and flow rate were 7 MPa and 96 kg/s [1-4, 40]. The HTR-PM plant, which is the world’s first modular PB-HTR, was connected to a grid in December 2021.

The HTR-PM reactor used a spherical fuel element of diameter 6 cm consisting of a 1 cm-thick graphite cladding and a 5 cm-diameter graphite matrix, and the densities of the graphite cladding and matrix are 1.74 g/cm³. TRISO coated particles were dispersed in a graphite matrix with a diameter of 0.92 mm. The fuel kernel for TRISO-coated particles was uranium dioxide with a diameter of 0.50 mm and density 10.4 g/cm³, that was successively coated with a porous carbon buffer (BufferC), inner pyrolytic carbon (Inner PyC), silicon carbide (SiC) and outer PyC from the inside to the outside (their thicknesses are 95, 40, 35, and 40 μm, respectively, and their density is 1.05, 1.90, 3.18, 1.90 g/cm³ in order) [8-10, 41].

The active zone of the HTR-PM unit was 3 m in diameter with an average height of 11 m, and could accommodate approximately 420,000 fuel pebbles randomly arranged in the core. The graphite reflectors, including the side, top and bottom reflectors, were placed immediately outside the core. top and bottom reflectors, respectively. The graphite reflector contains control rods and helium channels and carbon brick outside the graphite reflector. A schematic diagram of the core structure of the HTR-PM is shown in Fig. 1.

Fig. 1Full-screen View

(Color online) Core structure schematic diagram of the HTR-PM

The fuel pebbles in the HTR-PM unit were cycled 15 times through the reactor before reaching the specified burnup and the average discharged pebble burnup was 90 GWd/tHM. The fuel pebbles initially installed in the reactor core are basically not burnup. If fuel pebbles with high enrichment (8.5%) are directly loaded, the excess reactivity was relatively large and difficult to control. Therefore, the initial core loading of HTR-PM adopts a low-enrichment (4.2%) fuel pebble and graphite pebble mixed charging scheme. When the HTR-PM unit is in equilibrium and operates at full power, each reactor module discharges 6,000 high-enrichment (8.5%) fuel pebbles per day where approximately 400 fuel pebbles were spent fuel. Approximately 5,600 fuel pebbles were reloaded into the reactor core, and approximately 400 fresh fuel pebbles were replenished.

The VSOP code [42] can be used for the reactor physics calculations of the PB-HTR. According to the characteristics of the PB-HTR, VSOP code adopts the coupling of "energy calculation" for a certain space region and the whole reactor diffusion calculation, that is before each reactor diffusion calculation, the spectrum was calculated according to the current core state. The nuclear database of VSOP code is stored in libraries, including the nuclear data required for resonance calculation. spectrum calculation and burnup calculation.

First, VSOP-ZUT code [43] was used to create the resonance integral table. The neutron slowing-down equation of the inhomogeneous system is established, and the probability of neutrons escaping into the moderator is calculated using the cell’s geometric model. The narrow and wide resonance approximations simplify the neutron slowing-down equation, yielding the neutron slowing-down energy spectrum. The resonance integral is then solved. The effective resonance absorption cross-sections of ²³²Th, ²³⁸U, ²⁴⁰Pu, ²⁴²Pu and other nuclides at different densities and temperatures are calculated based on fuel element design and composition, considering fuel pebbles and coated particles’ double non-uniformity. The resonance integral table is added to the nuclear database. The effective resonance absorption cross-section is obtained through density and temperature interpolation, and added as background absorption cross-section to fast spectrum calculation. VSOP-MS code performs physical calculations, including spectrum, diffusion, burnup, and thermal hydraulics calculations. The spectrum calculation merges nuclides’ multi-group cross-sections into wide group cross sections, with 2–33 wide groups and 1 thermal group. PB-HTR uses 4 wide groups: 1 thermal and 3 fast groups. The finite-mesh finite difference method with CITATION code calculates diffusion equation. The neutron balance calculation determines neutron leakage for spectrum calculation feedback. The diffusion calculation provides neutron flux and power distributions for burnup and thermal hydraulics calculations. The burnup equation’s solution yields nuclide density distribution, while thermal hydraulic calculation determines temperature distribution using power distribution. Both the nuclide density and temperature distributions must be fed back to the spectrum calculation to update the spectrum. Convergence results were obtained through repeated iterations.

The HTR-PM equilibrium core was modeled to study the factors influencing the refueling reactivity coefficient. The HTR-PM core was equivalent to a cylinder with a diameter of 3 m and height of 11 m. A two-dimensional model was used to simulate the fuel cycle of HTR-PM. The reactor core was divided into five channels, each of which was divided into 20 layers. There are 100 regions, each of which has the same volume. Based on this model, continuous refueling of HTR-PM can be simulated using the step method.

Reactor thermal power

Reactor thermal power P is the most basic parameter of nuclear reactor. When P increases, negative reactivity enters the core, called power loss. Positive reactivity must be introduced to compensate for this power loss to maintain stable reactor operation at the new power level. For the PB-HTR, positive reactivity can be introduced by raising the control rod and increasing refueling speed. At 250 MWth, the HTR-PM operation is simulated from transition to equilibrium core. We vary P from 125 to 250 MWth (125, 150, 175, 200, 225 and 250 MWth) to study the effect of P on refueling reactivity coefficient, with other parameters remaining as designed. To maintain core criticality, control rods are raised and thermal hydraulic parameters adjusted according to power level.

Based on the current core state, the refueling reactivity coefficient was calculated using the following perturbation method, and the results are shown in Fig. 2. Evidently, with an increase in P, the refueling reactivity coefficient gradually increases primarily due to the distribution of the neutron flux.

Fig. 2Full-screen View

Change of refueling coefficient reactivity with reactor thermal power

It can be expressed as(11)where A(r) and B(r) are functions of the nuclear cross-section.

Equation (11) shows that the square of neutron flux, as the cross-section disturbance weight, correlates positively with reactivity disturbance. After increasing reactor thermal power, ¹³⁵Xe atomic density increases, requiring raised control rods to supplement reactivity, so the neutron flux peak moves up to the reactor core top with higher macroscopic cross-section disturbance, as shown in Fig. 3, increasing the refueling reactivity coefficient.

Fig. 3Full-screen View

(Color online) Neutron flux in the innermost channel and control rod positions for different powers

To explain the neutron flux modulation effect, we calculate the axial distribution of refueling reactivity with flat neutron flux, and compare it to normal neutron flux distribution. A significant difference exists between reactivity distributions in normal versus flat neutron flux cases. With identical nuclide atomic density, the neutron flux modulation makes axial refueling reactivity distribution closer to neutron flux distribution. Reactor thermal power positively affects the refueling reactivity coefficient by influencing neutron flux. As shown in Fig. 2, the refueling reactivity coefficient changes minimally with P - when P decreases from 250 MWth to 125 MWth, the coefficient reduces by only 6.2%. In power regulation, neutron flux modulation dominates the relationship between refueling reactivity coefficient and P, though its effect is limited. As P increases, neutron flux distribution becomes steeper with an upward-shifting peak, enhancing modulation and upper core reactivity contribution, thereby increasing the refueling reactivity coefficient.

Enrichment of fresh fuel

The HTR-PM is a typical thermal reactor, using the uranium-plutonium fuel cycle with ²³⁵U as the primary fission nuclide. The enrichment of ²³⁵U in fresh fuel ϵ determines the number of fissile nuclei in a reactor. We vary ϵ from 6% to 11% (6%, 6.5%, 7%, 7.5%, 8%, 8.5%, 9%, 9.5%, 10%, 10.5% and 11%) to study the effect of ϵ on the refueling reactivity coefficient. The HTR-PM reactor reaches equilibrium core at 250 MWth, and refueling reactivity coefficients with enrichment of fresh fuel at same and different discharge burnup are calculated, as shown in Fig. 4. The change amplitude of refueling reactivity coefficient and ϵ are similar in both cases, as changes in ϵ affect refueling reactivity mainly through nuclide atomic density gradient distribution, primarily of ²³⁵U. The analysis uses Eq. (10). With same discharge burnup, refueling reactivity coefficient decreases as ϵ increases, shown in Fig. 4a. When discharge burnup remains constant, fuel pebbles pass through the core 15 times with same cyclic burnup; thus, ²³⁵U atomic density gradient changes minimally with increasing ϵ. The numerator in Eq. (10) remains similar, while the denominator increases with ϵ, causing the refueling reactivity coefficient to decrease.

Fig. 4Full-screen View

Change in the refueling coefficient reactivity with enrichment of fresh fuel. (a) Discharge burnup is the same; (b) discharge burnup is different

When fresh fuel pebbles with different ϵ are used, discharge burnup increases with ϵ, enabling full nuclear fuel utilization. The change of refueling reactivity coefficient with ϵ in different discharge burnup is emphasized. Taking 90 GWd/tHM discharge burnup at ϵ of 8.5% as benchmark, the fuel pebble’s average residence time is adjusted according to Eq. (12) for critical reactor state.(12)where BU_dis is the average discharge burnup, T is the average residence time of the fuel pebble in the reactor core, P is the reactor thermal power, and m_HM is the heavy metal mass of the reactor core.

With different discharge burnup, the refueling reactivity coefficient increases with ϵ, as shown in Fig. 4b. Figure 5a shows the ²³⁵U atomic density distribution in the innermost channel under different ϵ, where higher ϵ leads to greater axial density gradient. According to Eq. (12), ignoring the modulation effect of neutron flux and the influence of other nuclides, the relationship between the refueling reactivity coefficient and ²³⁵U atomic density can be approximated as(13)where αQ is the refueling reactivity coefficient, N₅ is the average atomic density of ²³⁵U along the reactor core axis, ΔN₅ is the difference for the atomic density of ²³⁵U at the top and bottom of the reactor core, and C₁ is constant.

Fig. 5Full-screen View

(Color online) Change of atomic density of ²³⁵U in the innermost channel and with enrichment of fresh fuel (Discharge burnup is different). (a) Atomic density of ²³⁵U in the innermost channel; (b)

Figure 5(b) shows the relationship between the refueling reactivity coefficient and in Eq. (13). As ϵ increases, ΔN₅ rises faster than N₅, increasing the refueling reactivity coefficient.

A comprehensive analysis of the impact of enrichment on the refueling reactivity coefficient underlines the gradient of nuclide atomic density within the PB-HTR as the primary factor influencing the refueling reactivity introduced to the reactor core.

Loading fraction of fresh fuel

Fuel pebbles discharged from the reactor core bottom return to the top multiple times to reach specified burnup. Fresh and cycled fuel pebbles that haven’t reached specified burnup return to the core, where the fresh fuel pebble proportion is called the loading fraction f, which is reciprocal to the number of fuel pebble passing the core n:(14)We set n from 5 to 18 (5, 6, 7, 8, 10, 12, 15, and 18), and accordingly, f from 5.56% to 20% (5.56%, 6.67%, 8.33%, 10%, 12.5%, 14.29%, 16.67%, and 20%). Then, we adjust n to match the corresponding f, while maintaining the same discharge burnup.

The calculation results are presented in Fig. 6. The refueling reactivity coefficient increases with increasing f, which is expected.

Fig. 6Full-screen View

Change in the refueling coefficient reactivity with loading fraction of fresh fuel

It can be expressed as(15)The change in the refueling reactivity coefficient parallels the change in the f, increasing by 112% when f rises from 10% to 20%. The refueling reactivity coefficient represents reactivity change from unit refueling amount, including fresh and recycled fuel pebbles. When f increases, the introduced reactivity change increases, though not linearly. From Fig. 6, as f increases, the growth rate of refueling reactivity coefficient gradually increases, due to neutron flux’s modulating effect on refueling reactivity distribution. Figure 7 shows the neutron flux distribution in the reactor core’s axial direction under different f. As f increases, the neutron flux peak moves up toward the reactor core top with higher atomic density gradient, causing faster increase in refueling reactivity coefficient.

Fig. 7Full-screen View

(Color online) Neutron flux in the innermost channel of the HTR-PM core for different fresh fuel loading fractions

The analyzed effect of f on the refueling reactivity coefficient shows the absence of a completely linear relationship between f and the refueling reactivity coefficient. This result is attributed to the modulating effect of the neutron flux on the distribution of refueling reactivity, according to the Eq. (10), affecting the refueling reactivity coefficient. Thus, curve fitting is used instead of linear fitting in Fig. 6.

Heavy metal loading per fuel pebble

The fuel pebble in the HTR-PM contains 12,000 TRISO coated particles. The fuel kernel is uranium dioxide, surrounded by four coating layers for reactor safety. The design value of heavy metal loading per fuel pebble m_HM is 7 g. The refueling reactivity coefficient is calculated based on Eq. (10). We vary m_HM from 6 g to 10 g to study its effect on the refueling reactivity coefficient. The reactor thermal power and discharge burnup remain constant for different heavy metal loadings, with refueling speed adjusted to maintain core criticality.

The results are shown in Fig. 8a, showing the refueling reactivity coefficient gradually decreases as the m_HM increases, with relatively large changes. When the m_HM increases from 6 g to 10 g, the refueling reactivity coefficient decreases by 35%. This occurs because m_HM affects the neutron spectrum and atomic density of ²³⁵U. As m_HM increases, both the product of neutron yield with microscopic fission cross-section of ²³⁵U () and microscopic absorption cross-section of ²³⁵U (σ_a5) decrease due to neutron spectrum hardening. Consequently, N₅ at the reactor core bottom increases significantly at the same discharge burnup. The increase of ΔN₅ is less than the increase of N₅ as m_HM increases. According to Eq. (13), the gradually decreases. The results of in the innermost channel for different heavy metal loadings are shown in Fig. 8b, decreasing with increasing m_HM, matching theoretical analysis. Heavy metal loading negatively affects the refueling reactivity coefficient through neutron spectrum and atomic density of ²³⁵U. The driving force of refueling reactivity is the gradient of nuclide atomic density in PB-HTR. Additionally, neutron flux modulation affects the refueling reactivity coefficient according to Eq. (10), making the relationship between m_HM and refueling reactivity coefficient non-linear and thus necessitating curve fitting in Fig. 8a.

Fig. 8Full-screen View

Change in the refueling coefficient reactivity and with heavy metal loading per fuel pebble