Study on the Solution Method for Ultra-Fine Group Slowing-Down Equations Applicable to Stochastic Media

1. Introduction

As a low-carbon, clean, and high-energy-density baseload energy source, nuclear power is playing an irreplaceable strategic role in addressing global climate change and achieving carbon neutrality goals. According to the 2024 Energy Outlook Report [1] by the International Atomic Energy Agency (IAEA), nuclear power generation accounts for 10.2% of global electricity supply, with this proportion reaching 20-50% in advanced nuclear energy countries. As the mainstream reactor type in current nuclear power plants, the Pressurized Water Reactor (PWR) employs light water as both moderator and coolant. The accuracy of its core physics modeling directly impacts critical aspects such as core design optimization, operational strategy formulation, and safety margin evaluation [2].

In the field of reactor physics analysis, traditional deterministic methods, which are based on diffusion theory or transport theory for mathematical modeling, face significant challenges when dealing with strongly doubly heterogeneous media [3]. Particularly with the advancement of Generation IV nuclear energy system development, the engineering application of advanced fuel forms such as Fully Ceramic Micro-encapsulated fuel (FCM) and Mixed Oxide fuel (MOX) has introduced pronounced double heterogeneity (DH) effects that pose severe challenges to conventional analysis methods. This DH effect arises from the heterogeneous distribution characteristics at different spatial scales within the fuel elements: in FCM fuel, Tri-structural Isotropic (TRISO) fuel particles with diameters of approximately 0.8-1.0 mm are randomly dispersed in a graphite matrix; whereas in MOX fuel, plutonium-rich spots formed during the manufacturing process exhibit random spherical distributions at sub-millimeter scales [4].

From the perspective of neutron physics, this doubly heterogeneous structure affects core performance in three key aspects: First, in the resonance absorption region (0.1 eV–10 keV), the spatial self-shielding effect of fuel particles causes traditional homogenization methods to overestimate the effective resonance integral. Second, during neutron transport, the interface effects between randomly distributed fuel phases and matrix phases induce localized perturbations in neutron flux. Finally, at the core power distribution level, microscopic heterogeneity may lead to power peaking factors (PPF) exceeding safety limits due to neutron spectrum hardening effects [5].

The U.S. Nuclear Regulatory Commission (NRC) explicitly states in the NUREG-0800 Standard Review Plan that for DH fuel assemblies, validated high-fidelity modeling methods must be employed for safety analysis [6]. Among these challenges, effectively treating resonance cross-sections has become a key research difficulty. To overcome the limitations of traditional methods, the international academic community is currently advancing research along two main technical pathways:

(1) Developing multi-scale coupling algorithms based on the Monte Carlo method, which involve constructing high-resolution geometric models at the fuel particle scale while incorporating variance reduction techniques to improve computational efficiency [7].

(2) Enhancing deterministic methods by refining ultra-fine group slowing-down equation solutions and developing equivalent homogenization methods suitable for stochastic media [8].

For example, the APOLLO3 code developed by the French Alternative Energies and Atomic Energy Commission (CEA) employs a self-shielding correction model based on the subgroup method, converting the random distribution of TRISO particles into a probability density function for statistical treatment [9]. Meanwhile, the MC2-3 code developed by Argonne National Laboratory (ANL) in the U.S. achieves efficient coupling between Monte Carlo methods and multi-group transport equations [10].

Currently, deterministic neutron transport equation solvers remain the mainstream approach for reactor physics numerical simulations. Therefore, it is essential to develop neutron transport equation-solving methods suitable for randomly distributed media under doubly heterogeneous conditions.

To address this issue, Bende et al. [11] proposed the Dancoff method in 1999. This method first derives the neutron first-flight collision probability and transmission probability to calculate the probability of a neutron escaping from one fuel particle to another. Based on this escape probability, the average Dancoff factor is computed, which is then used to correct the cross-sections. However, the derivation of the Dancoff method was based on a non-coated fuel particle model for High-Temperature Gas-cooled Reactors (HTGRs), making it unsuitable for the multi-layered spherical structure of Fully Ceramic Micro-encapsulated (FCM) fuel. Ji et al. [12] improved the Dancoff method by proposing an average Dancoff factor calculation method based on a double-sphere model. By incorporating the chord length method to account for the multi-layered spherical structure, they extended the Dancoff method to TRISO particle calculations. Kim et al. [13] introduced a method using the average Dancoff factor to correct the background cross-section under double heterogeneity (DH) conditions and then inversely derived the effective resonance cross-section using resonance integral tables. However, due to the complexity of Dancoff factor calculations and the numerous approximations involved, this approach introduces additional biases. William [14,15] and She et al. [16] proposed the disadvantage factor method, initially applied to HTGRs. The core idea of this method is to simplify the doubly heterogeneous system into a conventional singly heterogeneous rod lattice problem, allowing the direct use of traditional resonance calculation methods to determine effective cross-sections. He [17] and Yin et al. [18] combined the disadvantage factor method with subgroup methods and ultra-fine group methods, applying it to FCM fuel resonance cross-section calculations. By homogenizing TRISO particles and matrix materials into a 1D model to compute the defect factor, this approach can be effectively integrated with traditional resonance calculation methods. However, it struggles to accurately resolve the fine spatial distribution of resonance cross-sections and neutron flux inside fuel particles.

Sanchez and Pomraning [19] proposed the Sanchez-Pomraning method to handle doubly heterogeneous geometries, which can be coupled with conventional transport methods such as the Method of Characteristics (MOC) [20]. This method iteratively computes the equivalent total cross-sections and source terms for each layer of the particle and the matrix material, then defines an update equation to correct the angular flux exiting the particle. After obtaining the matrix flux, the flux reconstruction method is used to recover the flux distribution inside each layer of the particle. The Sanchez-Pomraning method laid the foundation for solving subgroup fixed-source equations under double heterogeneity. Pogosbekyan et al. [21] applied this method to subgroup resonance calculations. However, resonance interference effects under double heterogeneity are more complex, and the background iteration method suffers from low accuracy. Since calculating interference factors under double heterogeneity is extremely challenging, traditional subgroup resonance interference methods fail to accurately account for these effects.

Currently, the most precise method for handling resonance interference is the ultra-fine group method [22]. Its fundamental approach involves extremely fine energy group division in the resonance energy range to avoid intra-group resonance peak treatment. Within each fine energy group, a simplified neutron transport equation (also called the neutron slowing-down equation) is established. The neutron flux distribution is then computed top-down for each energy group and flat source region, yielding a problem-dependent ultra-fine group neutron spectrum to achieve accurate resonance self-shielding treatment. The ultra-fine group method originated from Kier’s [23] approach of calculating two-region cell resonance integrals in fine energy groups. Building on this, the first practical ultra-fine group code, RABBLE [11], was developed. Since then, the ultra-fine group method has been widely applied in neutronics codes [24,25,26]. However, due to the need to solve the slowing-down equation in ultra-fine (or continuous-energy) groups, the computational cost is high. Codes like PEACO [13] significantly improve computational efficiency by precomputing collision probability interpolation tables, but their geometric modeling capability is limited. Alternatively, MOC-based fixed-source calculations can enhance geometric flexibility [27,28], though at the cost of reduced computational efficiency compared to PEACO-like methods

Therefore, it is necessary to investigate slowing-down equation solvers for stochastic media based on ultra-fine group methods. Theoretical innovations in neutron transport for randomly distributed media and the development of corresponding computational tools have become one of the most pressing challenges in reactor physics. The scientific value of this research lies not only in advancing fundamental neutron transport theory but also in its significant engineering implications for next-generation nuclear fuel development, reactor safety analysis, and improving the economic competitiveness of nuclear energy systems.

2. Materials and Methods

2.1. Ultra-fine group slowing-down equation

The energy range of the resolved resonance region typically spans 1 eV to 10 keV. Since fission neutrons have energies significantly higher than 10 keV, it can be assumed that no fission neutrons are generated within this energy range. Additionally, neutrons below 10 keV hardly undergo inelastic scattering, meaning the scattering source term can be simplified to account only for isotropic elastic scattering, while up-scattering can be neglected. Under these conditions, the neutron slowing-down equation can be expressed in the following form[29,30]:

$$\Omega \nabla \varphi \left(\overrightarrow{r},u,\overrightarrow{\Omega }\right)+{\Sigma }_t\left(\overrightarrow{r},u\right)\varphi \left(\overrightarrow{r},u\right)={\displaystyle {\int }_{u-\ln \left(1/\alpha \right)}^u\frac{{\Sigma }_s\left(\overrightarrow{r},u^{\prime }\right)\varphi \left(\overrightarrow{r},u^{\prime }\right)}{1-\alpha }}{e}^{u^{\prime }-u}\mathrm{d}{u}^{\prime }$$

(1)

Where α=[(A-1)/(A+1)]²，A represents the atomic mass, u is lethargy calculated as u = ln(E₀/E)，E₀=10MeV.

If there are more than one nuclides in region k, Eq. (1) could be re-written as:

$$\Omega \nabla \varphi \left(\overrightarrow{r},u,\overrightarrow{\Omega }\right)+{\displaystyle \sum _i{\Sigma }_{i,t}^k\left(\overrightarrow{r},u\right)\varphi \left(\overrightarrow{r},u\right)}=Q\left(\overrightarrow{r},u\right)$$

(2)

$$Q\left(\overrightarrow{r},u\right)={\displaystyle \sum _i{\displaystyle {\int }_{u-{\Delta }_i}^u{\Sigma }_{i,s}^k\left(\overrightarrow{r},u^{\prime }\right)\varphi \left(\overrightarrow{r},u^{\prime }\right)\frac{1}{1-\alpha }}{e}^{u^{\prime }-u}\mathrm{d}{u}^{\prime }}$$

(3)

Where ${\Sigma }_{i,x}^k$ represents the macroscopic cross-section for the i-th nuclide and reaction type x; Δ denotes the maximum logarithmic energy decrement that nuclide i can attain, and 1−α is the maximum energy fraction that nuclide i may lose in a single collision. By solving the slowing-down equations sequentially from the highest energy downward based on continuous-energy point cross-sections, the continuous neutron scalar flux can be obtained. Subsequently, using Equation (4), the within-group effective self-shielded cross-sections in multigroup format can be derived.

$${\sigma }_{x,G}={\displaystyle \frac{{\displaystyle {\int }_{\Delta u_G}{\sigma }_x\left(u\right)\varphi \left(u\right)du}}{{\displaystyle {\int }_{\Delta u_G}\varphi \left(u\right)du}}}$$

(4)

Where G is the group index.

When applying the above conclusions to practical resonance calculations, for a system containing multiple resonant and non-resonant materials, we first spatially discretize the system into multiple mesh regions. Within each mesh region, the source term is assumed to be uniform. Subsequently, the resonance energy range is subdivided into extremely fine energy groups, termed ultra-fine groups. To ensure that neutrons undergo only one collision within the energy range covered by a single ultra-fine group, the energy width of each ultra-fine group is set much smaller than the maximum logarithmic energy decrement (Δ) attainable by neutrons after a single collision with the heaviest nuclide in the system. Under these conditions, the scattering source in ultra-fine group g within a mesh region can be formulated as:

$$Q_g={\displaystyle \sum _{j=1}^J{\displaystyle \sum _{n=1}^{L_j}{P}_{nj}{\Sigma }_{s,g-n}^j{\varphi }_{g-n}}}\Delta u_f$$

(5)

Where P_nj represents the probability that a neutron, after undergoing a single collision with material jj, traverses nn energy groups and enters group g. L_j denotes the maximum number of energy groups that a neutron may traverse after a single collision with material j.

Since the number of reachable energy groups after neutron scattering is extremely large, directly solving Equation (5) would be computationally prohibitive. Therefore, a reformulation is necessary. First, consider P_nj the probability that a neutron’s energy decreases from u′to u after scattering can be expressed as:.

$$P\left(u^{\prime }-u\right)={\displaystyle \frac{1}{1-\alpha }}{e}^{-\left(u-u^{\prime }\right)}$$

(6)

Where u− ∆≤ u’ ≤ u.

The probability that neutrons pass through n ultra-fine groups after scattering is：

$$P_n\Delta u_f={\displaystyle \frac{1}{1-\alpha }}{\displaystyle {\int }_{u_0}^{u_0+\Delta u_f}du}{\displaystyle {\int }_{u_0-n\Delta u_f}^{u_0-\left(n-1\right)\Delta u_f}d{u}^{\prime }}{e}^{-\left(u-u^{\prime }\right)}$$

(7)

Substituting equation (6) into equation (7), we can obtain:

$$P_n\Delta u_f={\displaystyle \frac{1}{1-\alpha }}{\left(1-e^{-\Delta u_f}\right)}^2{e}^{-\left(n-1\right)\Delta u_f}$$

(8)

From equation (8), the following rule can be derived:

$$P_n\Delta u_f={\displaystyle \frac{1}{1-\alpha }}{\left(1-e^{-\Delta u_f}\right)}^2{e}^{-\left(n-1\right)\Delta u_f}$$

(9)

$$P_n=e^{-\Delta u_f}{P}_{n-1}$$

(10)

Substituting equations (9) and (10) into equation (5), we get:

$$\begin{array}{c}{Q}_g=P_1{\Sigma }_{g-1}{\varphi }_{g-1}+[P_2{\Sigma }_{g-2}{\varphi }_{g-2}+P_3{\Sigma }_{g-3}{\varphi }_{g-3}+\mathrm{...}+P_L{\Sigma }_{g-L+1}{\varphi }_{g-L+1}+\\ P_L{\Sigma }_{g-L}{\varphi }_{g-L}+P_{L+1}{\Sigma }_{g-L-1}{\varphi }_{g-L-1}]-P_{L+1}{\Sigma }_{g-L-1}{\varphi }_{g-L-1}\end{array}$$

(11)

According to equation (9), it can be seen that the term in parentheses in equation (11) is $Q_{g-1}{e}^{-\Delta u_f}$，where $Q_{g-1}$epresents the scattering source term of the (g-1)th ultra-fine group. Thus, the recursive relationship for the ultra-fine group scattering source term can be obtained:

$$Q_g=Q_{g-1}{e}^{-\Delta u_f}+P_1{\Sigma }_{g-1}{\varphi }_{g-1}-P_{L+1}{\Sigma }_{g-L-1}{\varphi }_{g-L-1}{e}^{-\Delta u_f}$$

(12)

For the recursive relationship in Equation (11), the first step is to calculate the scattering source of the leading ultra-fine group. The moderation neutrons in the leading ultra-fine group originate from the energy range above the resonance calculation region. In this higher energy range, the cross-sections vary relatively smoothly with energy. Therefore, the cross-sections of nuclides in this energy range can be approximated as constants, and the flux distribution within this region is assumed to follow the 1/E law. Based on Equation (5), the scattering source of the leading ultra-fine group can then be computed. Once the scattering source of the leading ultra-fine group is obtained, the neutron flux in this group can be determined through a single-group fixed-source calculation. The resulting neutron flux is then used to compute the scattering source for the next energy group. By iterating this process sequentially, a fine-resolution energy spectrum across the entire ultra-fine group range can be derived, ultimately enabling the computation of the resonance self-shielded cross-sections for a specific resonance group.

The accuracy of the fixed-source calculations directly affects the precision of the resonance self-shielded cross-sections. In this work, the method of characteristic is coupled to perform the intra-group fixed-source calculations for ultra-fine groups. The Method of Characteristics (MOC) transforms the neutron transport equation into equations along characteristic lines at specific angular directions using mathematical techniques. This approach simultaneously handles spatial and angular variables, enabling treatment of arbitrary geometries. The modular MOC approach primarily refers to the characteristic line generation methodology. This method generates characteristic lines on a modular basis, ensuring that lines in adjacent modules are seamlessly connected end-to-end.

2.2. Stochastic Medium Ultra-fine Group Equations Based on the Sanchez-Pomraning Method

In traditional pressurized water reactor (PWR) problems, only the heterogeneity arising from rod-and-lattice structures exists. However, for accident-tolerant fuels (ATFs) such as fully ceramic microencapsulated (FCM) fuel, the fuel particles are dispersedly distributed within a silicon carbide (SiC) matrix inside the fuel rod. As shown in Figure 1, the random distribution of fuel particles introduces a second layer of heterogeneity into core neutronics calculations. Additionally, for MOX fuel, due to manufacturing limitations, plutonium isotopes in the fuel rod form spherical agglomerations, where PuO₂ particles are dispersed within a UO₂matrix, also creating a dual heterogeneity phenomenon. Since the precise coordinates of the particles are difficult to determine, resonance calculations under double heterogeneity conditions cannot be performed directly. To address this, international approaches such as the volume-weighting method, Dancoff factor method, and disadvantage factor method have been proposed. These methods fundamentally aim to compute equivalent models by homogenizing fuel particles and the matrix, thereby converting the dual heterogeneity problem into a conventional single-heterogeneity problem, which is then solved using standard resonance and transport methods. A common limitation of these methods is their inability to resolve the fine spatial distributions of cross-sections and flux within individual particles, as they only provide homogenized resonance cross-sections for the mixed material. In recent years, the Sanchez-Pomraning method has been increasingly applied to MOC transport calculations for dual heterogeneity. This method iteratively updates effective cross-sections and computes the matrix flux using governing equations, then reconstructs the detailed flux distribution inside particles based on the matrix flux. The subgroup method exhibits intrinsic compatibility with transport solvers, making the Sanchez-Pomraning method equally applicable to subgroup resonance calculations. This section proposes a Sanchez-Pomraning-based ultrafine group equation solution framework. First, the traditional MOC method is enhanced using the Sanchez-Pomraning approach to enable transport calculations in dual-heterogeneity geometries. Subsequently, this improved method is applied to both ultrafine group equations and multigroup transport equations, achieving refined calculations for dual-heterogeneity problems.

Due to the difficulty in determining the coordinates of randomly distributed fuel particles, the conventional Method of Characteristics (MOC) cannot perform characteristic line scanning for fuel particles within the flat source region. Therefore, it is necessary to first homogenize the total cross-sections of each layer of the particles and the matrix. Unlike the conventional volume-weighted method, the Sanchez-Pomraning method iteratively calculates the equivalent total cross-section of the flat source region based on the probability of neutrons escaping from each layer of the fuel particles, and calculates the equivalent source term through the equivalent total cross-section and the double heterogeneity asymptotic flux. During the characteristic line scanning process, the renormalization factor is calculated through the update equation, and then the outgoing flux of the flat source region is corrected. The flux in the flat source region that converges iteratively through this process is the flux of the FCM fuel rod matrix, and the flux of each layer of the fuel particles can be calculated according to the reconstruction formula [31,32,33].

The trajectory of neutrons within the fuel particle is illustrated in Figure 2. Here, P_i→j represents the probability that a neutron located in the i-th spherical shell moves to the j-th shell and undergoes its first collision, while Ei denotes the probability that a neutron in the i-th spherical shell escapes the fuel particle without any collision. R is the radius of the spherical layer, and the length of each segment t can be directly determined based on geometric relationships. For an infinitely thin spherical layer dr at a distance r from the center of the sphere, P_dr→j can be calculated using Equation (13):

$$P_{dr\to j}={\displaystyle \frac{2\pi rdr}{{\Sigma }_i{V}_i}}\left(e^{-t_{i-1,j-1}^{+}}+e^{-t_{i,j-1}^{+}}\right){\left(1-e^{-t_i}\right)}^2$$

(13)

Where Σ_i and V_i are macro total cross section and the volume of region i respectively.

In Equation (13), integrating over dr from 0 to R_i yields the value of P_i→j. Since different spherical shells of the fuel particle have distinct material cross-sections, the integration must be performed segment-wise across the respective shells. The final expression for P_i→j is given by Equation (13), and the numerical integration can be solved using the Gauss-Seidel method.

$$P_{i\to j}={\displaystyle \frac{2\pi }{{\Sigma }_i{V}_i}}{\displaystyle \sum _{k=1}^i{\displaystyle {\int }_{R_{k-1}}^{R_k}rdr\left(e^{-t_{i-1,j-1}^{+}}+e^{-t_{i,j-1}^{+}}\right){\left(1-e^{-t_i}\right)}^2}}$$

(14)

Similar to P_i→j, the probability P_i→i that a neutron undergoes a collision within its current fuel layer can be calculated using Equation (15):

$$P_{i\to i}={\displaystyle \frac{2\pi }{{\Sigma }_i{V}_i}}{\displaystyle \sum _{k=1}^i{\displaystyle {\int }_{R_{k-1}}^{R_k}rdr\left[e^{-t_{i-1,i-1}^{+}}{\left(1-e^{-t_i}\right)}^2+2\left({\tau }_i+e^{-t_i}-1\right)\right]}}$$

(15)

Based on the reciprocity relation (14), the reverse collision probability P_j→i can be derived inversely using Equation (13).

$$V_i{\Sigma }_i{P}_{i\to j}=V_j{\Sigma }_j{P}_{j\to i}$$

(16)

The escape probability of neutrons from the i-th layer can be calculated using Equation (17) after obtaining the inter-layer collision probabilities.

$$E_i=1-{\displaystyle \sum _{j=1}^N{P}_{i\to j}}$$

(17)

Where N is the layer number of the fuel particle.

The transport calculation requires iterative solution of the equivalent total cross-section for the flat-source region, as given by Equation (18):

$${\Sigma }_t^{\left(n+1\right)}=\Sigma {}_{t,m}+{\displaystyle \frac{1}{p_m}}{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{i=1}^I{p}_{hi}}}\left[{\Sigma }_{t,hi}-{\Sigma }_t^{\left(n\right)}\right]{\widehat{E}}_{hi}$$

(18)

Where the subscript h represents particles; m represents the matrix; p_m and p_hi denote the volume fractions of the matrix and the i-th fuel layer of the h-th fuel particle, respectively; and $\widehat{E}$ is the reduced escape probability.

The calculation process of $\widehat{E}$ is similarly based on Equations (15)-(17), where the total cross-section value is taken as Σ_t,hi minus ${{\Sigma }_{\left(n\right)}}_t$. Starting from Equation (18), ${{\Sigma }_{\left(n\right)}}_t$ is first initialized to 0, followed by iterative computation of ${{\Sigma }_{(n+1)}}_t$ until convergence is achieved. Using the converged equivalent total cross-section Σ_t, the equivalent source term can then be calculated via Equation (19):

$$q={\Sigma }_t{\varphi }_{as}$$

(19)

Where ϕ_as denotes the asymptotic flux. Under doubly heterogeneous conditions, its calculation formula is given by Equation (20):

$${\varphi }_{as}=\left(q_m+{\displaystyle \frac{1}{p_m}}{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{i=1}^I{p}_{hi}}}{q}_{hi}{E}_{hi}\right)/{\tilde{\Sigma }}_t$$

(20)

Where ${\tilde{\Sigma }}_t$ is obtained by Equation (21):

$${\tilde{\Sigma }}_t=\Sigma {}_{t,m}+{\displaystyle \frac{1}{p_m}}{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{i=1}^I{p}_{hi}}}{\Sigma }_{t,hi}{E}_{hi}$$

(21)

Define q_m  and q_hi as the actual source terms for the matrix and each layer of the fuel particle, respectively. For energy group g, their calculation methods are given by Equation (22):

$$q_g=\left\{\begin{array}{cc}{\displaystyle \sum _{g\text{'}}^G{\Sigma }_{hi,s,g\text{'}\to g}{\varphi }_{hi,g\text{'}}}+{\displaystyle \frac{1}{k_{eff}}}{\chi }_g{\displaystyle \sum _{g\text{'}}^{G}v{\Sigma }_{hi,f,g\text{'}}{\varphi }_{hi,g\text{'}}}& Particle\\ {\displaystyle \sum _{g\text{'}}^G{\Sigma }_{m,s,g\text{'}\to g}{\varphi }_{m,g\text{'}}}+{\displaystyle \frac{1}{k_{eff}}}{\chi }_g{\displaystyle \sum _{g\text{'}}^{G}v{\Sigma }_{m,f,g\text{'}}{\varphi }_{m,g\text{'}}}& Matrix\end{array}\right.$$

(22)

where: ϕ_m and ϕ_hi denote the actual fluxes in the matrix and each layer of the fuel particle, respectively.

The update equation is calculated based on collision probabilities and escape probabilities as shown in Equation (23):

$${\varphi }_{out}\left(x\right)=\left\{\begin{array}{cc}{\displaystyle {\int }_0^{L}dy\left(S_h{Q}_h\left(y\right)+{\varphi }_{in,h}\left(x-y\right)f_h\left(y\right)\right)e^{-{\Sigma }_{h}y}}& Particle\\ {\varphi }_{out,m}\left(0\right)e^{-\left({\Sigma }_m+1/\overline{l_m}\right)x}+{\displaystyle {\int }_0^{x}dy\left(S_m+{\varphi }_{in,m}\left(x-y\right)/l_m\right)e^{-\left({\Sigma }_m+1/\overline{l_m}\right)y}}& Matrix\end{array}\right.$$

(23)

Where Q represents the cumulative chord length distribution function, f denotes the probability cumulative function, S is the source term, l stands for the chord length.

The Sanchez-Pomraning method simplifies Equation (23) through a boundary layer assumption, introducing a renormalized factor r_c to correct the outgoing neutron angular flux ϕ_out  in the Method of Characteristics (MOC) flat-source regions, which is shown in Equation (24):

$${\varphi }_{out}={\varphi }_{in}+r_c\left(1-e^{-{\Sigma }_{t}L}\right)\left({\varphi }_{as}-{\varphi }_{in}\right)$$

(24)

For flat-source regions containing fuel particles, r_c can be calculated using Equation (24). For regions without fuel particles (such as the moderator or cladding), r_c takes a value of 1, in which case Equation (24) becomes identical to the conventional MOC outgoing angular flux calculation formula. Under doubly heterogeneous conditions, the flux obtained through Equation (24) correction represents the actual matrix flux ϕ_m.

$$r_c=p_m+{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{i=1}^I{p}_{hi}{\widehat{E}}_{hi}}}$$

(25)

After obtaining the matrix flux, the flux for each layer of the fuel particles can then be calculated using Equation (26):

$${\varphi }_{hi}={\widehat{E}}_{hi}{\varphi }_m+\left(E_{hi}-{\widehat{E}}_{hi}\right){\varphi }_{as}+{\displaystyle \frac{1}{V_{hi}{\Sigma }_{hi,t}}}{\displaystyle \sum _{i=1}^I{V}_{hi}{q}_{hi}{P}_{h,i\to j}}$$

(26)

This study presents an enhanced ultra-fine group (UFG) approach incorporating modifications derived from the Sanchez-Pomraning (SP) technique, hereafter referred to as the UFGSP method. The removal cross-section formulation employed in the ultra-fine group slowing-down equations maintains consistency with Equation (18). By integrating Equations (12) and (22), we derive the comprehensive source terms for both particulate and matrix components within the ultra-fine group framework, as presented in Equation (27).

$$q_g=\left\{\begin{array}{cc}\begin{array}{l}{Q}_{hi,g-1}{e}^{-\Delta u_{g-1}}+{\displaystyle \sum _n^N{f}_{hi,1,n}{\Sigma }_{hi,n,s,g-1}}{\phi }_{hi,g-1}-\\ {\displaystyle \sum _n^N{f}_{hi,L,n}{\Sigma }_{hi,n,s,g-L-1}}{\phi }_{hi,g-L-1}{e}^{-\Delta u_{g-L-1}}\end{array}& -Particle\\ \begin{array}{l}{Q}_{m,g-1}{e}^{-\Delta u_{g-1}}+{\displaystyle \sum _n^N{f}_{m,1,n}{\Sigma }_{m,n,s,g-1}}{\phi }_{m,g-1}-\\ {\displaystyle \sum _n^N{f}_{m,L,n}{\Sigma }_{m,n,s,g-L-1}}{\phi }_{m,g-L-1}{e}^{-\Delta u_{g-L-1}}\end{array}& -Matrix\end{array}\right.$$

(27)

The computational procedure initiates with fast group flux determination using asymptotic flux values, followed by sequential resonance group flux calculations. The condensed effective resonance parameters for the predefined group structure is shown in Equation (28).

$${\sigma }_x={\displaystyle \frac{{\displaystyle \sum _{g_r=1}^{G_r}{\sigma }_{x,g_r}{\varphi }_{g_r}}}{{\displaystyle \sum _{g_r=1}^{G_r}{\varphi }_{g_r}}}}$$

(28)

Where x represents the reaction channel type, G_r denotes the number of ultra-fine groups corresponding to the current resonance energy interval in the multi-group structure.

The overall calculating process UFGSP is shown in Figure 3. Notably, the method requires just two applications of the Sanchez-Pomraning technique: first for doubly heterogeneous slowing-down equations and subsequently for multi-group transport solutions [34,35,36].

3. Results

In this section, various fully ceramic micro-encapsulated (FCM) fuel problems with different packing fractions and geometric configurations are selected for validation. Additionally, calculations and analyses are performed on issues involving poison particles and plutonium spots in traditional MOX fuels to verify the capability of the resonance module in the UFGSP code for handling doubly heterogeneous problems.

3.1. Typical FCM fuel problem

The basic structure of the FCM fuel is illustrated in Figure 1. The fuel kernel of the TRISO particles consists of UC fuel with an enrichment of 17.8 wt%, and the matrix material is SiC. The geometric configurations are listed in Table 1, and the material nuclide density information is provided in reference [37]. This section selects seven typical FCM fuel test cases. The first five cases have particle packing fractions of 1%, 10%, 20%, 30%, and 40%, respectively, with TRISO particles of Geometry A. The sixth and seventh cases correspond to Geometry B and Geometry C, both with a packing fraction of 30%. The computational parameters of MOC are set as follows: 16 azimuthal angles and 3 polar angles per octant, a ray-spacing of 0.01 cm, and a system temperature of 300 K. The reference solutions are provided by the Monte Carlo code OpenMC [38], where the calculations employ randomly packed particles and an explicit modeling approach.

Table 2 presents the basic information and effective multiplication factor k_eff results for different test cases. The computational accuracy is high across all packing fractions, with a maximum absolute deviation of 64 pcm. Taking Cases 2, 5, 6, and 7 as examples, Figure 4 shows the computational deviations of the absorption cross-sections for ²³⁵U and ²³⁸U in this series of problems. Test cases with lower packing fractions exhibit higher dilution levels and larger background cross-sections, resulting in relatively larger resonance cross-sections, whereas cases with higher packing fractions show smaller resonance cross-sections. UFGSP can accurately compute this phenomenon. For the 40% packing fraction case, the computational deviations for both ²³⁵U and ²³⁸U cross-sections remain within ±1%, demonstrating high accuracy. For the 10% packing fraction case, the deviation in the ²³⁸U cross-section increases, reaching relative deviations of 1.35% and 1.49% in the energy ranges of 13.71–29.02 eV and 78.89–130.07 eV, respectively. Since the resonance cross-sections in these energy ranges are relatively small, their impact on reaction rates is minor. Cases 6 and 7, which feature different geometric dimensions, also achieve high accuracy in resonance cross-section calculations, with most deviations staying within ±1% and the maximum deviation remaining around 1%. Table 3 summarizes the maximum (MAX), average (AVG), and root-mean-square (RMS )deviations in resonance cross-section calculations for different test cases, all of which exhibit high computational precision.

Taking Case 6 as an example, Figure 5 presents the radial distribution of the ²³⁸U absorption cross-section within the fuel kernel of TRISO particles. Similar to the spatial self-shielding effect observed in fuel rods, the fuel kernel also exhibits neutron shielding. Neutron reactions primarily concentrate near the surface of the fuel kernel, while the resonance cross-section in the inner region shows an increasing trend along the radial direction. The calculation of radial distribution is crucial for subsequent burnup analysis. However, traditional doubly heterogeneous resonance methods, such as the defect factor method or the Dancoff method, cannot resolve the internal resonance cross-section distribution within the fuel kernel. UFGSP based on the Sanchez-Pomraning method enables precise computation of the radial cross-section distribution in the fuel kernel. As shown in Figure 5, the computational deviations for the three resonant energy groups remain below ±1%, demonstrating high accuracy..

3.2. Burable poison and Pu spot problem

This section validates two special doubly heterogeneous problems: FCM fuel containing poison particles and the plutonium spot issue in conventional PWRs, as shown in Figure 6. For the poison particle problem, QUADRISO particles are conventional TRISO particles with an additional layer of burnable poison coated outside the fuel kernel, while BISO particles replace the TRISO fuel kernel entirely with a burnable poison kernel. The plutonium spot problem refers to the phenomenon where PuO₂ particles are randomly dispersed as spherical inclusions within UO₂ fuel rods in MOX fuel, where traditional methods typically homogenize the plutonium spots with the fuel matrix, whereas UFGSP can directly perform doubly heterogeneous calculations for plutonium spots. The MOC computational parameters for this series of problems are set with 16 azimuthal angles and 3 polar angles per octant, a ray spacing of 0.01 cm, and a system temperature of 300 K, with reference solutions provided by the Monte Carlo code OpenMC.。

3.2.1. Poison problems

The dimensions and material types of the QUADRISO and BISO particles selected in this section are shown in Table 4, with the poison materials including both B₄C and Gd₂O₃. The specific nuclide number densities are provided in references [37] and [39]. Four test cases were designed for this study: the first two cases involve fuel rods containing QUADRISO particles with a packing fraction of 30%, using B₄C and Gd₂O₃ as the poison materials respectively. The latter two cases feature BISO particles, where the fuel rods contain 33.8% TRISO particles and 10.8% BISO particles. The TRISO particle types are identical to those listed in Table 1, with the burnable poisons being B₄C and Gd₂O₃ respectively.

Figure 7(a) presents the computational deviations of ²³⁸U absorption cross sections in the fuel kernels for all four test cases, with maximum deviations remaining within ±1%, demonstrating high computational accuracy. Figure 7(b) shows the absorption cross-section deviations for ¹⁵⁵Gd and ¹⁵⁷Gd in the two Gd₂O₃ test cases. Compared with ²³⁸U, the maximum deviations for Gd cross-sections are slightly larger due to more pronounced flux depression within Gd-containing particles. Table 5 summarizes the cross-section calculation results for both B₄C and Gd₂O₃ poison problems, showing overall high computational accuracy. Table 6 presents the k_eff for all four poison test cases. For B₄C cases and the Gd₂O₃ BISO case, k_eff deviations remain within ±100 pcm. The deviation increases to 138 pcm for the Gd₂O₃ QUADRISO case, which features relatively small k_eff values and is rarely encountered in practical PWR applications.

3.2.2. Pu spots problem

In the plutonium spot problem, the matrix material is UO₂ while the fuel particles are uncoated PuO₂. Compared with FCM fuel problems, the main challenge of plutonium spot problems lies in the fact that the fuel rod matrix itself exhibits strong resonance effects, which conventional doubly heterogeneous resonance calculation methods cannot simultaneously handle for both the matrix and fuel particles. For the test cases in this section, the matrix consists of UO₂ with 0.16 wt% enrichment, while the particles are PuO₂ with 1.5 wt% enrichment. The plutonium spots have a radius of 100 μm with packing fractions of 0.5%, 1%, 1.5% and 2% respectively. Detailed material nuclide compositions and geometric information are provided in references [40,41].

Taking the plutonium spot problems with 1% or 2% packing fractions as examples, Figure 8(a) presents the calculation deviations of ²³⁸U absorption cross-sections in the UO₂ matrix for these two cases. Although the packing fraction of plutonium spots is extremely small, the strong resonance interference effects of Pu isotopes significantly influence the ²³⁸U resonance cross-sections. As the packing fraction increases, the dilution effect on ²³⁸U intensifies, leading to larger resonance cross-section values. The matrix resonance cross-section calculations for different cases demonstrate high accuracy, with deviations for most energy groups remaining within ±1%. Figure 8(b) shows the calculation deviation distributions of ²³⁹Pu and ²⁴¹Pu absorption cross-sections in the 1% packing fraction plutonium spots. Both Pu isotopes exhibit strong resonance peaks, with ²⁴¹Pu showing particularly intense resonance peaks below 10 eV. UFGSP accurately describes the resonance peak distributions of ²³⁹Pu and ²⁴¹Pu, with calculation deviations near resonance peaks all within ±1%. The maximum deviations across the entire energy range are -1.24% and -1.42% respectively. Table 7 summarizes the overall resonance cross-section calculation deviations for the plutonium spot problems, showing high accuracy for both UO₂ matrix and PuO₂ fuel particle resonance cross-section calculations. Table 8 presents the effective multiplication factor calculation deviations for the plutonium spot problems, with absolute deviations for all cases remaining within ±100 pcm.

4. Discussion

The UFGSP method advances deterministic neutron transport by addressing two critical challenges in DH systems: spatial self-shielding in stochastic media and resonance interference between matrix and particle phases. Traditional approaches, such as the Dancoff factor or defect factor methods, rely on homogenization techniques that fail to resolve intra-particle flux gradients or matrix-particle resonance coupling. In contrast, UFGSP leverages the Sanchez-Pomraning iterative framework to decouple particle and matrix transport, enabling precise calculation of equivalent cross-sections and angular flux corrections. However, computational costs remain a concern. While MOC-based angular discretization improves geometric flexibility, iterative cross-section homogenization and flux reconstruction increase runtime compared to conventional MOC. Future work could explore acceleration techniques, such as machine learning-aided collision probability precomputation or adaptive energy group condensation.

5. Conclusions

This study successfully develops and validates the UFGSP method, a high-precision solver for ultra-fine group slowing-down equations in doubly heterogeneous media. By integrating the Sanchez-Pomraning technique with MOC-based ultra-fine group slowing-down equation, the method overcomes the limitations of traditional homogenization approaches, enabling detailed resolution of resonance cross-sections in both stochastic particle distributions and matrix materials. Key achievements include:

(1) Accurate modeling of radial flux gradients within TRISO particles, critical for predicting power peaking and burnup effects.

(2) High-fidelity treatment of complex geometries, including BISO/QUADRISO poison particles and PuO₂ spots, with k_eff deviations consistently below 138 pcm.

(3) Demonstrated applicability to advanced fuel types, supporting the design and safety analysis of Generation IV reactors.

The method provides a practical tool for reactor physicists to address DH challenges in accident-tolerant fuels, MOX assemblies, and other advanced nuclear systems. Future extensions to multi-physics coupling and transient analysis could further solidify its role in next-generation reactor innovation.