Study on Higher-Order Harmonic Calculation of Neutron Diffusion Equation and Its Application in Core Power Monitoring of the Gas-Cooled Microreactor

1  Introduction

The rapid development of gas-cooled microreactors (GMRs) has garnered significant attention for remote power supply and modular energy applications due to their inherent safety, high energy density, and structural versatility. As a cornerstone of next-generation nuclear energy, the integration of digital twins, autonomous control, and intelligent fault diagnosis is becoming a definitive trend to ensure reliable operation in unmanned or remote environments [1–5]. However, the unique design characteristics of GMRs—such as their extremely compact core, high-temperature helium cooling, and the absence of complex in-core instrumentation—pose stringent challenges for real-time core power monitoring. Accurate reconstruction of the three-dimensional (3D) power distribution is not merely a requirement for maintaining operational safety margins; it is the fundamental precursor for optimizing fuel burnup and predicting the thermo-mechanical integrity of the core under high-temperature gradients.

In recent years, the academic community has explored a diverse range of methodologies to achieve efficient power reconstruction. Traditional high-fidelity Monte Carlo (MC) simulations and fine-mesh deterministic codes offer superior accuracy but are computationally prohibitive for online applications. To bridge this gap, various reconstruction algorithms have been proposed. For instance, the Least-Squares method has been successfully applied to PWRs for power mapping [6], while Radial Basis Function (RBF) interpolation has demonstrated excellence in processing scattered measurement data [7]. Moreover, high-precision online monitoring systems, such as SOPHORA for Hualong One, have validated the use of self-powered neutron detectors (SPNDs) for real-time applications [8]. For microreactors specifically, recent research has gravitated toward dimensionality reduction techniques like Principal Component Analysis (PCA) to reconstruct power from limited ex-core detector signals [9].

Despite these advancements, existing methods face critical bottlenecks. Data-driven approaches, including PCA and Bayesian Neural Networks (BNN), typically rely on the construction of extensive “snapshot” libraries. These libraries must encompass a vast array of operational states, including various control rod configurations and temperature distributions. For GMRs, which may exhibit complex transient behaviors or “grey-box” uncertainties, the data dependency of these methods poses a significant risk: if the actual operational state deviates from the pre-computed database, the reconstruction error can escalate rapidly, leading to potential safety misjudgments. Furthermore, while the harmonic expansion method (HEM) provides a robust physics-based alternative [10], recent iterations utilizing fission matrix-based harmonics [11–13] still encounter substantial computational overhead during the matrix generation phase, particularly for specific 3D hexagonal geometries.

Currently, there is a lack of integrated frameworks that combine multi-group diffusion theory with higher-order harmonic expansion specifically optimized for the 3D hexagonal-z geometry of GMRs. This study aims to fill this research gap by developing a dedicated higher-order harmonic calculation module within a diffusion framework. Unlike purely data-driven methods, this approach is rooted in the physical governing equations, eliminating the need for exhaustive pre-calculations while ensuring higher adaptability to diverse and unforeseen operational transients. The contribution of this paper is threefold: (1) it establishes a theoretical formulation for higher-order harmonics in hexagonal geometry; (2) it demonstrates the implementation of an inverse solving algorithm using limited ex-core detectors; and (3) it provides a systematic sensitivity analysis to quantify the impact of measurement uncertainties on reconstruction precision. This work establishes a critical technical foundation for the engineering realization of autonomous monitoring systems in advanced microreactors.

2  Theoretical Model

2.1 Multi-Group Neutron Diffusion Equations in 3D Hexagonal-z Geometry

The core power monitoring of gas-cooled microreactors is fundamentally rooted in the accurate description of neutron transport and reaction rates. In this study, the multi-group neutron diffusion theory is adopted as the physical basis due to its computational efficiency and sufficient accuracy for global power distribution analysis. For a 3D reactor core, the steady-state multi-group diffusion equation is expressed as a set of coupled partial differential equations, expressed as follows (the spatial variable r of the neutron flux is omitted for brevity):

−Dg▽2ϕg+Σr,gϕg=χgkeff∑g′=1G(νΣf)g′ϕg′+∑g′=1g′≠gGΣs,g′→gϕg′(1)

In this equation, g and g′ denote the incident and outgoing energy groups, respectively; D is the diffusion coefficient; ▽2 is the Laplacian operator, representing the sum of second-order derivatives along all spatial coordinates; ϕg and ϕg′ are the neutron fluxes of the incident and outgoing energy groups; Σr is the removal cross section, which is calculated using Eq. (2).

Σr,g=Σa,g+∑g=1g≠g′ngΣs,g→g′(2)

χ is the fission spectrum; keff represents the fundamental eigenvalue, corresponding to the critical state of the reactor. However, for power reconstruction, we are not only interested in the fundamental mode but also in the higher-order spatial modes that capture the flux perturbations caused by control rod movements or localized xenon oscillations; G is the number of energy groups; ν is the average number of neutrons released per fission; Σf is the fission cross section; and Σs,g′→g is the scattering matrix.

The core assemblies of gas-cooled microreactors are hexagonal prisms. Thus, the hexagonal node method is adopted for geometric discretization. Using the transverse integration method [14,15] to separate spatial variables, for each node, the 3D neutron flux is decomposed into a 2D radial component and a 1D axial component. For the radial geometry, the governing equation can be written as:

−Dg(∂2∂x2+∂2∂y2)ϕgRm(x,y)+ΣrgϕgRm(x,y)−∑g′<gΣg′gϕg′Rm(x,y)=QgRm(x,y)−LgZm(x,y),m=1…6(3)

where R and Z denote the radial and axial directions, respectively; m is the index of the triangular subregion; QgRm(x,y) represents the radial source term, and LgZm(x,y) represents the axial leakage term.

For the axial geometry,

−Dg∂2∂z2ϕgZ(z)+∑rgϕgZ(z)−∑g′<gΣg′gϕg′Z(z)=QgZ(z)−LgR(z)(4)

where QgZ(z) is the axial source and LgR(z) is the radial leakage.

The source term is calculated as follows:

QgU(u)=χgkeff∑g′vΣfg′ϕg′R, old (u)+∑g′>gΣg′gϕg′R, old (u),u=x,y,z(5)

and the leakage is calculated as follows:

LgZm(x,y)=1hz(JgTZm(x,y)−JgBZm(x,y))LgR(z)=239hR∑u=13(Jg(u+3)x(z)−Jgux(z))(6)

where JgTZm(x,y) and JgBZm(x,y) represent the top-face and bottom-face fluxes of the subregion, respectively.

2.2 Higher-Order Harmonic Expansion Method

The core idea of the Harmonic Expansion Method [16–18] is to treat the actual fission rate distribution as a superposition of a series of orthogonal spatial eigenfunctions (harmonics) and transform the 3D monitoring problem into a coefficient estimation problem.

The harmonic modes {ψn} are defined by the generalized eigenvalue problem derived from the diffusion operator:

Lψn=1knMψn(7)

where L is the operator representing neutron loss (diffusion, absorption, and out-scattering) and M is the fission production operator. The eigenvalues kn are ordered such that k0>k1>k2>…, where k0 corresponds to the fundamental mode.

Any perturbed power distribution P can be reconstructed using a truncated series of these harmonics:

P≅∑n=0Nanψn(8)

In GMRs, where the core is compact and the flux shape is relatively “stiff,” a limited number of harmonics is sufficient to capture the majority of power distribution characteristics. This dimensionality reduction is the key to enabling real-time online monitoring.

To prevent the higher-order modes from collapsing back to the fundamental mode due to numerical round-off errors, a Gram-Schmidt Orthogonalization process is executed after each outer iteration. For any new mode ψn, orthogonality to all previously calculated modes {ψ0,…,ψn−1} with respect to the fission production operator M is ensured.

ψn(new)=ψn−∑i=0n−1⟨ψn,Mψi⟩⟨ψi,Mψi⟩ψi(9)

where the inner product ⟨⋅,⋅⟩ denotes integration over the entire core volume.

Specifically, the implementation is carried out based on the adjoint flux method. The adjoint neutron flux is used to characterize the “value” of neutrons, i.e., the contribution or importance of neutrons to physical quantities in the system, such as reaction rate, power, or reactivity. Its form is the same as that of Eq. (1), and the corresponding equation is as follows:

−Dg▽2ϕg∗+Σr,gϕg∗=(νΣf)gkeff∑g′=1Gχg′ϕg′∗+∑g=1g′≠gGΣs,g→g′ϕg′∗(10)

where ϕg∗ denotes the adjoint flux. The differences from Eq. (1) are: (1) transpose the scattering matrix energy group (i.e., change g→g′ to g′→g); (2) exchange the fission source and the fission spectrum.

To ensure the fundamental accuracy of the global power distribution, the fundamental mode (zeroth-order) forward and adjoint fluxes are calculated using a refined 25-group energy structure. However, extending this multi-group framework to higher-order modes would lead to a linear escalation in computational time with diminishing marginal returns on reconstruction precision. Therefore, a two-group approximation is employed for the higher-order forward and adjoint fluxes. This approach significantly accelerates the modal expansion process while preserving the essential spectral characteristics and orthogonality required for power reconstruction. Although the two-group approximation introduces minor energy-discretization biases compared to direct multi-group calculations, the resulting higher-order modes remain physically robust, and the reconstructed power distribution consistently satisfies the required accuracy standards.

Based on the zeroth-order flux and its adjoint flux, the first-order flux and adjoint flux can be obtained. Using recurrence relations, the calculation formulas for neutron flux and adjoint flux of any order can then be derived.

The calculation method for higher-order flux and higher-order adjoint flux are as follows:

ϕs,gn,(k)(r)=ϕs,gn,(k−1)(r)−∑l=0n−1∑r=1nz⋅nxyϕs,gl,∗(r)ψn,(k−1)(r)∑r=1nz⋅nxyϕs,gl,∗(r)ψl(r)ϕgl(r),g=1,2(11)

ϕs,gn,∗,(k)(r)=ϕs,gn,∗,(k−1)(r)−∑l=0n−1∑r=1nz⋅nxyϕs,1n,∗,(k−1)(r)ψl(r)∑r=1nz⋅nxyϕs,1l,∗(r)ψl(r)ϕs,gl,∗(r),g=1,2(12)

where l and n denote the flux order and the maximum order, respectively; (k) is the iteration index, and after each outer iteration the low-order components introduced by the outer iteration need to be removed; s represents the neutron spectrum that has not yet been normalized during the outer iteration process; r denotes the block location; and nz, nxy are the numbers of axial layers and radial blocks, respectively. The calculation method for the higher-order fission reaction rate ψn of the node is as follows:

ψn=V∑g=12ϕgn(νΣf)g(13)

where V is the volume of a node. It should be noted that the impact of the average number of neutrons generated per fission can be ignored and thus is not separated from the fission cross section.

2.3 Inverse Problem for Power Reconstruction

Due to the compact core volume and the harsh operating environment (high temperature and intense radiation) characteristic of gas-cooled micro-reactors, in-core instrumentation is significantly limited. Consequently, an ex-core monitoring strategy is employed. The ex-core axial detection system is configured with twelve detectors arranged in four layers, with three locations radially, as seen in Fig. 1, providing a spatially distributed signal field to capture the core’s leakage neutron information.

Figure 1: Radial (a) and axial (b) detectors layouts of GMR.

Assuming that the energy released per fission is consistent across different fuel regions [11–13], the three-dimensional power distribution P can be accurately represented as a superposition of the fundamental mode and higher-order harmonics of the fission reaction rate:

P≈∑n=0∞fnsn=Mx(14)

where sn and fn denote the harmonic vectors and their corresponding coefficients, respectively; M and x represent the harmonic matrix and the coefficient vector, respectively.

For the core power distribution P, after being transformed by the detector response matrix D, the ex-core detector signal R is obtained

DP=R(15)

The relationship between the internal power distribution P and the signal R measured by the ex-core detector is established through a response matrix D. Physically, each element of D represents the contribution of a unit power source at a specific core location to the detector response, which is typically pre-calculated using adjoint-based or fixed-source methods.

By expanding the power distribution using harmonic functions and combining it with Eq. (14), we obtain:

R=DP=D(Mx)=(DM)x(16)

Therefore, the problem of reconstructing the core power distribution from detector readings can be transformed into solving for the higher-order harmonic coefficients x. When the number of harmonics is smaller than the number of detectors (12), the optimal fitting coefficients x can be obtained using the least-squares method. Substituting these coefficients back into Eq. (14) then yields the core power distribution.

The integrated online monitoring workflow is illustrated in Fig. 2. The process is divided into two phases: (1) Offline Preparation: The high-fidelity Monte Carlo code RMC is utilized to model the three-dimensional core and accurately construct the response matrix, accounting for the complex transport effects between the core and ex-core detectors. (2) Online Reconstruction: During operation, the higher-order harmonic calculation module performs rapid three-dimensional diffusion calculations using real-time reactor state data (e.g., control rod positions, burnup, and temperature feedback). For the numerical implementation, a hexagonal-z nodal method is employed. The core is discretized into 90 blocks (30 radial fuel assemblies and 3 axial layers) to ensure a balance between computational speed and spatial resolution. By combining the real-time harmonics with live detector readings, the core power distribution is reconstructed instantaneously.

Figure 2: Core online monitoring workflow.

3  Numerical Results and Analysis

3.1 Online Core Monitoring and Verification of the Power Reconstruction Method

To validate the efficacy of the proposed power reconstruction methodology, a typical operating condition was defined by a specific set of parameters, including control rod position, burnup, and core temperature. The reference power distribution for this state was generated using the high-fidelity Monte Carlo code RMC. Subsequently, the ex-core detector signals were synthesized by applying the pre-calculated detector response matrix to this reference distribution.

In practical online monitoring scenarios, precise core state parameters are often unavailable or subject to measurement uncertainties. To simulate this, we assumed “approximate” state information: the control rod position, burnup, and core temperature were estimated to be near 80 cm, zero burnup, and 900 K, respectively. This approximate state serves as the input for the diffusion solver to generate the expansion basis (harmonics), which are then integrated with real-time detector readings to reconstruct the actual power distribution.

Fig. 3 illustrates the eigenvalues for the 0th to 12th orders of the diffusion equation under this specific state, while Fig. 4 displays the corresponding spatial distributions of the power harmonics. Due to the inherent radial symmetry of the gas-cooled microreactor core under nominal conditions, the harmonics exhibit distinct geometric characteristics. In Fig. 3, hollow markers denote eigenvalues associated with radially symmetric harmonic distributions. Conversely, solid markers indicate eigenvalues linked to degenerate harmonic modes, which correspond to the asymmetric harmonic shapes shown in Fig. 4. For the reconstruction process, the program automatically selected the 0th, 1st, and 6th order harmonics as the optimal expansion basis, determined by their ability to capture both the fundamental power level and potential dominate axial power variation.

Figure 3: 0th–12th order eigenvalues of the gas-cooled microreactor.

Figure 4: 0th–12th order harmonic distributions of the gas-cooled microreactor.

The reconstruction performance is evaluated through a systematic comparison as illustrated in Figs. 5–7. Fig. 5, from top to bottom, presents the reference power distribution from RMC, the predictive distribution from the Advanced Reactor Core Simulator (ARCS) diffusion solver, and the spatial distribution of their relative errors. Due to the initial parameter approximations, the standalone diffusion calculation exhibits a maximum deviation of 3.12% and a root-mean-square (RMS) deviation of 1.49%. To isolate the approximation capability of the generated harmonics, Fig. 6 displays the error obtained by directly fitting the RMC reference with the selected harmonics. The resulting maximum and RMS deviations are 1.50% and 0.75%, respectively. This outcome confirms that harmonics derived from diffusion theory possess sufficient spatial fidelity to serve as a robust basis for capturing the true power distribution under correlated operating conditions.

Figure 5: Power distribution of RMC and ARCS and their deviation.

Figure 6: Deviation of the harmonic fitted power distribution (without detector readings).

Figure 7: Deviation of the harmonic reconstructed power distribution (with detector readings).

The final power reconstruction, achieved by synthesizing the selected harmonics with real-time ex-core detector signals via the least-squares method, is shown in Fig. 7. This integrated approach suppresses the maximum deviation to 1.82% and the RMS deviation to 0.79%. A quantitative comparison reveals that, relative to the initial ARCS diffusion results, the harmonic-plus-detector reconstruction reduces the maximum and RMS deviations by 41.6% and 47.0%, respectively. Notably, the reconstruction precision achieved here significantly outperforms the approximately 2.0% RMS error threshold commonly reported in existing literature [11].

Furthermore, the robustness of this methodology was tested against uncertainties in reactor state observation, such as control rod positioning. As demonstrated in Fig. 8, while a ±3 cm mismatch in rod position can escalate the direct diffusion calculation error to 3.83%, the harmonic-based reconstruction effectively compensates for this discrepancy, maintaining the maximum deviation within 2.2%. This underscores the method’s ability to leverage ex-core measurements to rectify model-based predictions, ensuring highly reliable online monitoring even when core state parameters are imperfectly known.

Figure 8: MSE and RMSE of different control rod position.

3.2 Sensitivity Analysis: Impact of Detector Response Matrix and Measurement Errors on Reconstruction Accuracy

In real-world reactor environments, the parameters utilized in power reconstruction are subject to various sources of uncertainty. Specifically, the detector response matrix may undergo subtle variations due to localized changes in shielding or core state, while the precision of detector readings is inherently limited by the manufacturing and calibration accuracy of the instrumentation. To rigorously quantify the sensitivity of the HEM to these uncertainties, a series of stochastic numerical experiments were conducted.

For each sensitivity scenario, a simulation comprising 20,000 independent trials was performed. In each trial, random perturbations were introduced to every element of the response matrix and the detector reading vector. The resulting maximum relative deviations in the reconstructed power distribution were then statistically analyzed to establish confidence intervals and error propagation patterns.

Figs. 9–11 categorize these findings into three dimensions: the impact of response matrix perturbations on detector signals, the sensitivity of the reconstructed power distribution to detector reading fluctuations, and the combined effect of simultaneous perturbations in both components. The results demonstrate a clear linear relationship between the magnitude of individual perturbations and the resulting error in both detector signals and the final power distribution. Notably, while the combined effect of simultaneous perturbations is greater than that of any single source, it remains strictly lower than the arithmetic sum of the two individual effects. E.g., 2% perturbation of the response matrix may cause around 0.17% variation on average to the detector readings, which may further cause around 2.4% deviation of the reconstructed power; simultaneously and separately, 0.2% perturbation of the detector readings may cause around 2.8% deviation of the reconstructed power; all together around 5.2% deviation, higher than the actual mean value 3.3% with simultaneous perturbation. This “sub-linear” error accumulation suggests a degree of numerical robustness within the HEM framework, where certain randomized errors across the twelve detectors tend to partially offset one another during the least-squares fitting process.

Figure 9: The impact of perturbations in the detector response matrix on detector readings.

Figure 10: The impact of perturbations in detector readings on the power distribution.

Figure 11: The combined impact of simultaneous perturbations in both the response matrix and detector readings on the power distribution.

Given that the maximum variation in the detector response matrix across diverse operating conditions is approximately 1%–2%—though further research into the precise uncertainty bounds of these matrices is warranted—the simultaneous-perturbation analysis provides a critical engineering benchmark. The results indicate that to maintain the power reconstruction error within acceptable engineering design limits (typically <5%), the maximum allowable deviation in detector readings (i.e., detector accuracy) must be strictly controlled within 0.2%. This stringent requirement underscores the necessity for high-precision ex-core instrumentation and regular calibration to ensure the reliability of the online core monitoring system.

4  Conclusion

In this study, an advanced higher-order harmonic expansion methodology based on three-dimensional neutron diffusion theory was investigated and successfully integrated into the reactor core monitoring framework. By combining the detector response matrix with real-time measurement data, a robust system for the online core power distribution monitoring was established.

The numerical results demonstrate that for typical operating conditions, the power distribution reconstructed via higher-order harmonics yields a significant accuracy enhancement, exceeding 40% compared to conventional direct diffusion calculations. Furthermore, a comprehensive sensitivity analysis was conducted to quantify the impact of uncertainties in the detector response matrix and measurement readings on the reconstruction fidelity. The findings indicate that to satisfy stringent engineering design requirements and ensure operational safety, the precision of detector instrumentation must be maintained within a 0.2% tolerance. This research provides a high-performance and reliable computational tool for the real-time monitoring and digital twin development of advanced nuclear reactor systems.

Acknowledgement: Not applicable.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: Methodology and supervision, Peng Zhang; software and writing, Kui Hu; validation, Xiang Xiao; formal analysis, Yuan Yuan and Zhiyuan Feng; investigation, Yuan Xu and Yunhuang Zhang. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the Corresponding Author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Lan T, Long B, Peng X, Sheng F. A study on the seismic test of a planar model of the gas-cooled microreactor. Nucl Eng Technol. 2026;58(2):103984. doi:10.1016/j.net.2025.103984. [Google Scholar] [CrossRef]

2. Zhang Y, Lu D, Cao Q, Wang X, Xu C, Wang Y, et al. Experimental study on flow and heat transfer characteristics of the core rod bundle region in a space gas-cooled microreactor based on similarity analysis method. Prog Nucl Energy. 2026;190(1):105998. doi:10.1016/j.pnucene.2025.105998. [Google Scholar] [CrossRef]

3. Lan T, Feng T, Sheng F, Tan W. Internal stress analysis of irradiated graphite cores in a gas-cooled microreactor. Nucl Eng Des. 2024;429(7):113647. doi:10.1016/j.nucengdes.2024.113647. [Google Scholar] [CrossRef]

4. Doyle DL, Duchnowski EM, Trelewicz JR, Brown NR. Ceramic composite inert matrix fuel forms in high-temperature gas-cooled microreactors. Nucl Sci Eng. 2025;199(12):2234–53. doi:10.1080/00295639.2025.2475645. [Google Scholar] [CrossRef]

5. Guo J, Wang Y, Zhang H, Cui M, Li F. Challenges and progress of uncertainty analysis for the pebble-bed high-temperature gas-cooled reactor. Prog Nucl Energy. 2021;138(3):103827. doi:10.1016/j.pnucene.2021.103827. [Google Scholar] [CrossRef]

6. Lee K, Kim CH. The least-squares method for three-dimensional core power distribution monitoring in pressurized water reactors. Nucl Sci Eng. 2003;143(3):268–80. doi:10.13182/NSE03-A2335. [Google Scholar] [CrossRef]

7. Peng X, Wu Q, Cai Y, Lou L, Yu Y, Li Q. The application of radial basis function interpolation in reactor core power distribution on-line monitoring. Ann Nucl Energy. 2019;132(1):752–62. doi:10.1016/j.anucene.2019.06.059. [Google Scholar] [CrossRef]

8. Zhang X, Li W, Dang Z, Wang J, Li X, Li J, et al. Implementation and preliminary verification of 3D on-line core monitoring system: SOPHORA. High Power Laser Part Beams. 2017;29:016005. (In Chinese). doi:10.11884/HPLPB201729.160200. [Google Scholar] [CrossRef]

9. Mohsen MYM, Nagla TF, Elsamahy M, Abdel-Rahman MAE. Multivariate image analysis for core monitoring in PWRs. J Phys: Conf Ser. 2023;2616(1):012059. doi:10.1088/1742-6596/2616/1/012059. [Google Scholar] [CrossRef]

10. Meng ZY, Pan Y, Wang JP, Zhao X. Optical temporal interference model for investigation and manipulation of non-integer high-order harmonic generation. Chin Phys B. 2026;35(2):024208. doi:10.1088/1674-1056/adfdc5. [Google Scholar] [CrossRef]

11. Blaise P, Laureau A, Ros P, Leconte P, Routsonis K. Transient fission matrix approach for assessing complex kinetics behavior in the ZEPHYR ZPR coupled core configurations. Ann Nucl Energy. 2019;128:390–7. doi:10.1016/j.anucene.2019.01.031. [Google Scholar] [CrossRef]

12. Terlizzi S, Kotlyar D. Fission matrix decomposition method for criticality calculations: theory and proof of concept. Nucl Sci Eng. 2019;193(9):948–65. doi:10.1080/00295639.2019.1583948. [Google Scholar] [CrossRef]

13. He D, Li Y, Pan Q, Liu X. Development of a predictor-corrector quasi-static transient fission matrix combination method for the neutronics transient simulation. Nucl Eng Des. 2024;424(1):113301. doi:10.1016/j.nucengdes.2024.113301. [Google Scholar] [CrossRef]

14. Ikeda H, Takeda T. A new nodal SN transport method for three-dimensional hexagonal geometry. J Nucl Sci Technol. 1994;31(6):497–509. doi:10.1080/18811248.1994.9735184. [Google Scholar] [CrossRef]

15. Li Z. Study on the discrete ordinates axial-moment hexagonal nodal method and CMFD acceleration. At Energy Sci Technol. 2019;53(02):271–6. (In Chinese). [Google Scholar]

16. Urase Y, Yamamoto A, Endo T. Verification of in-core power distribution reconstruction based on POD for high-temperature gas-cooled reactors using ex-core detectors. J Nucl Sci Technol. 2026;63(1):66–81. doi:10.1080/00223131.2025.2533530. [Google Scholar] [CrossRef]

17. Wang C, Wu H, Cao L, Yang P. On-line reconstruction of in-core power distribution by harmonics expansion method. Nucl Eng Des. 2011;241(8):3042–50. doi:10.1016/j.nucengdes.2011.05.013. [Google Scholar] [CrossRef]

18. Liang J, Li J, Wang C, Chen J, Lu H, Yang Y, et al. Calculation of adjoint flux in the lattice code PINE and validation on neutron kinetics parameters. Ann Nucl Energy. 2025;213:111170. doi:10.1016/j.anucene.2024.111170. [Google Scholar] [CrossRef]