SP3-coupled global variance reduction method based on RMC code 增强出版

SP3-coupled global variance reduction method based on RMC code 1 corresp first-author panqingquan@sjtu.edu.cn panqingquan@sjtu.edu.cn panqingquan@sjtu.edu.cn https://orcid.org/0000-0002-9384-0070 1 https://orcid.org/0000-0002-6763-2534 1 1 2 https://orcid.org/0000-0003-0704-8069 School of Nuclear Science and Engineering, Shanghai JiaoTong University, Shanghai 200240, China School of Nuclear Science and Engineering, Shanghai JiaoTong University, Shanghai 200240, China Department of Engineering Physics, Tsinghua University, Beijing 100084, China Department of Engineering Physics, Tsinghua University, Beijing 100084, China A global variance reduction (GVR) method based on the SPN method is proposed. First, global multi-group cross-sections are obtained by Monte Carlo (MC) global homogenization. Then, the SP3 equation is solved to obtain the global flux distribution. Finally, the global weight windows are approximated by the global flux distribution, and the GVR simulation is performed. This GVR method is implemented as an automatic process in the RMC code. The SP3-coupled GVR method was tested on a modified version of the C5G7 benchmark with a thickened water shield. The results show that the SP3-coupled GVR method can improve the efficiency of the MC criticality calculation. RMC code Global homogenization Variance reduction SPN theory INTRODUCTION 1 The reactor criticality calculation is an essential step in the physical design of a reactor. The reactor burnup and neutronic/thermal (N-TH) coupling calculations [ 1 1 ] are based on the reactor criticality calculation. If the accuracy of the criticality calculation cannot be guaranteed, the accuracy of the reactor burnup and N-TH coupling calculations will also not be guaranteed. The Monte Carlo (MC) method is widely used in reactor criticality calculations because of its high fidelity. However, a non-uniform distribution of the relative errors has been found in MC criticality calculations. This non-uniform distribution of the relative errors has a negative influence on the efficiency and accuracy of some local tallies. In some cases, it even causes numerical instability in the MC multi-physics simulation [ 2 2 ] and MC kinetic calculation [ 3 3 ]. To improve the efficiency and accuracy of MC criticality calculations and ensure the numerical stability of MC multi-physics simulations, global variance reduction (GVR) methods should be studied. There are many variance reduction (VR) methods, such as Zheng’s method [ 4 4 , 5 5 ], the adaptive method [ 6 6 ], the FW-CADIS method [ 7 7 ], the automated weight windows method [ 8 8 ], and the importance function method [ 9 9 ]. These methods are well-known for their value in VR applications. However, these methods were initially developed for MC fixed-source calculations and not for MC criticality calculations. Some VR methods have also been d eveloped for MC criticality calculations. They include the UFS [ 10 10 ], UTD [ 11 11 ], and UVM [ 12 12 ] methods. These methods are easily implemented and have good universality. However, the UFS method assumes that the fission neutrons are uniformly distributed in the regions with fuel and have little influence on the non-fuel regions. Moreover, only the source bias technique and not the weight window technique is adopted in the UFS method. This limits its value for GVR. The UTD method can reduce the variance of various statistics, but because its basic principle is still the same as that of the UFS method, it faces the same limitations. The UVM method is newly proposed, and can improve the efficiency of the MC criticality calculation and flatten the global variance. However, the development of this method is still in the initial stages and needs to be further studied. A SP3-coupled GVR method is proposed here to improve the accuracy and efficiency of the MC criticality calculation. The SP3-coupled method can be used for GVR in the MC criticality calculation and allows the analysis of the reactor from its interior to its exterior through a single MC criticality calculation. Using the SP3-coupled GVR method, an accurate radiation source can be provided for reactor shielding design through the MC criticality. Moreover, the SP3-coupled GVR method can be coupled with the reactor burnup calculation to analyze the influence of reactor dynamics on reactor shielding. The SP3 method is a numerical method for solving the neutron transport equation with a simplified treatment of spatial anisotropy, and has been widely studied [ 13 13 , 14 14 ]. The traditional SP3 method starts from a one-dimensional model. The one-dimensional governing equation is first found and the one-dimensional operator is subsequently replaced with a three-dimensional operator to obtain the three-dimensional governing equation. Compared with traditional transport theory, the spatial anisotropy is simplified in the SP3 method to achieve higher calculation efficiency. Compared with traditional diffusion theory, the SP3 method provides a better treatment of spatial anisotropy to achieve higher calculation accuracy. Therefore, the SP3 method balances accuracy and efficiency, and has great potential for coupling with the MC method to achieve accelerated calculations and GVR. In this paper, a hybrid MC method based on SP3 theory is studied and a GVR method driven by the SP3 theory is developed for the RMC code [ 15 15 ]. Unlike the traditional two-step reactor physical calculation method, a one-step process of global MC homogenization and global reactor calculation is achieved in this method, and an automatic calculation framework for GVR is constructed. Unlike other hybrid MC methods, such as the fission matrix [ 16 16 ] and response matrix [ 17 17 ] methods, the neutron transport equation is solved instead of a matrix relationship. This provides advantages in terms of numerical stability, development potential, and application scope. The remainder of this paper is organized as follows: Section 2 introduces the theory, and provides the related formulas and details. Section 3 presents the numerical results with some basic analysis. A discussion is given in Section 4 before the conclusion. METHODS 1 The SP3-coupled GVR method is implemented as an automatic process in the RMC code. This implies that only the input card of the RMC code is required, and the GVR calculation will subsequently be automatically completed by the RMC code. The flowchart for the calculation is shown in Figure 1 Figure 1 . The calculation consists of three main steps: (1) The global cross-sections are determined by MC global homogenization. (2) The SP3 equation is solved to obtain the global flux distribution. (3) The global weight windows are determined by the global flux. Each of the three steps is described in the following subsections. MC Global Homogenization 2 The SP3-coupled GVR method is a hybrid MC method in which the MC method is coupled with the SP3 method. A deterministic reactor code [ 18 18 ] for solving the SP3 equation, NLSP3, is developed. In contrast to the discrete ordinate method (SN), global multi-group cross-sections (MGCS) are required to solve the SP3 equation. Therefore, the MGCS must first be determined. MC global homogenization was adopted to determine the MGCS for the following reasons: (1) Because a continuous-energy nuclear library is used, no resonance effect treatment [ 19 19 ] is required. (2) Because the homogenization process is not based on an assembly model, no leakage correction [ 20 20 ] is required. (3) Because the geometric model is not approximated, a higher accuracy can be achieved compared with traditional assembly homogenization. (4) Because the global homogenization is based on the entire reactor, the cross-sections of all the material regions can be determined through a single calculation. In particular, Point (4) implies that MC global homogenization is suitable for an automated process. Because the MGCS used for the MC global homogenization have a significant influence on the global flux distribution, a high accuracy of the MGCS must be ensured. Meanwhile, because the cross-sections in all the material regions are determined simultaneously, a non-uniform distribution of the relative errors will lead to different cross-section precisions in different material regions. To improve the accuracy of the cross-sections and ensure that the cross-sections in all material regions have the same precision, the following measures are taken [ 21 21 ]: (1)A fixed-source calculation with a flat neutron source is performed for MC global homogenization. The spatial distribution of neutrons is positively correlated with the neutron flux distribution if the VR technique is adopted in the MC simulation. In a large reactor, the neutron fluxes in different phase spaces may differ, especially when the reflective and shielding layers are considered. These differences result in a variation in the spatial distribution of neutrons. Therefore, although performing global homogenization is a reasonable approach in the MC criticality calculation, it is inevitable that the cross-sections in some regions will have higher precision compared to other regions. The spatial distribution of neutrons can be artificially set in a fixed-source calculation. Suppose a flat neutron source is provided in the fixed-source calculation. In this case, all the material regions have the same density of neutron histories, and the cross sections obtained will all have the same precision. Therefore, the MC global homogenization in the fixed-source calculation is performed with a flat neutron source. (2)The neutron fission reaction is restricted during neutron transport. Various reactions such as scattering, fission, and absorption occur during neutron transport. The reactor becomes critical when the loss of a neutron because of the absorption of the neutron or its transportation outside the system is compensated by the corresponding production of a neutron in a fission reaction. Typically, a fission reaction produces two or three neutrons. In a supercritical system, the loss of one neutron is accompanied by the production of more than one neutron. Therefore, if fission neutrons are produced during the fixed-source calculation, the number of simulated neutrons will increase, and the calculation will never be terminated. As stated earlier, global homogenization is performed during the MC fixed-source calculation. To ensure the feasibility of the fixed-source calculation, neutron fission is restricted during MC global homogenization: The neutrons only perform random walks in the system, and no fission reaction occurs. However, the various reaction probabilities are recorded during neutron transport. (3)The critical spectrum is approximated using the Watt fission spectrum. The MC homogenization involves a process of compressing and grouping the neutron spectra according to the energy. The neutron spectrum is essential in MC global homogenization. However, because the critical spectrum is unknown before the fixed-source calculation, it can only be approximated. As stated before, neutron fission is restricted during the neutron transport process. This restriction leads to errors. Therefore, the neutron spectrum can be artificially adjusted according to the fission spectrum when the flat neutron source is set during the fixed-source calculation. It is known that the fission spectrum of U-235 is described by the Watt spectrum [ 22 22 ]. Therefore, the external fission source is distributed according to the Watt spectrum <math display="block" id="M1"><mrow><mi>p</mi><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo><mo>=</mo><mi>C</mi><msup><mi>e</mi><mrow><mo>−</mo><mi>E</mi><mo>/</mo><mi>a</mi></mrow></msup><mi>sinh</mi><msqrt><mrow><mi>b</mi><mi>E</mi></mrow></msqrt><mo>,</mo></mrow></math> where C , a , and b are coefficients that are ordinarily equal to 0.453 MeV, 0.965 MeV, and 2.29 MeV -1 , respectively. It should be noted that although the Watt spectrum may not be the best one for the external source and is not universal for different fission nuclides, it is the best approximation that we can come up with. (4) Corrections are adopted to make the neutron spectrum asymptotic to the criticality spectrum. The main consideration in the Watt fission spectrum are the fission reactions that occur during the neutron transport process. The Watt fission spectrum is a good approximation in the reactor core region where many fission reactions occur. However, in regions without fissile materials, the neutron spectrum is mainly influenced by scattering reactions, and the error will be larger when the Watt spectrum is used to approximate the critical spectrum. Therefore, there should be as many scattering reactions as possible before the neutrons are lost. After many collisions, the neutron energy spectrum in the reflective and shielding layers become close to the critical spectrum. To make the neutron spectrum asymptotically approach the critical spectrum, the cutoff weight limit is reduced when more scattering reactions occur. In the SP3-coupled GVR method in this study, the weight cutoff threshold was set to 1.0E-8, which is the default value in the RMC code for the VR simulation. In the nonanalog MC method, the scattering reactions are effectively increased by reducing the cutoff weight windows. After multiple collisions in the reflective and shielding layers, the energy spectrum of the neutrons becomes close to the critical spectrum. However, at this stage, the weight of the neutrons is small. The neutrons therefore have little influence on the tallied results. To increase the influence of the neutrons after multiple collisions on the cross sections, the reaction rate is recorded based on the number of reactions rather than the weight of the neutrons. We calculated the MGCS using a fixed-source calculation, which is an approximation. Thus, the MGCS obtained inherently has errors, which implies that the efficiency is improved at the cost of compromised calculation accuracy. The SP3 solution only provides the parameters for VR, and the final result is given by the MC criticality calculation. Thus, even if the accuracy of the MGCS solution is compromised, the error in the MGCS solution does not affect the accuracy of the final result. This is one of the reasons why we adopt an approximate method to solve for the MGCS. The SP3 equation is solved to obtain the global flux distribution. The total grouped cross-sections <math id="M2"><mrow><msub><mi>Σ</mi><mrow><mtext>t,</mtext><mi>g</mi></mrow></msub></mrow></math> , grouped absorption cross-sections <math id="M3"><mrow><msub><mi>Σ</mi><mrow><mtext>a</mtext><mo>,</mo><mi>g</mi></mrow></msub></mrow></math> , grouped fission production cross-sections <math id="M4"><mrow><mi>ν</mi><msub><mi>Σ</mi><mrow><mtext>f</mtext><mo>,</mo><mi>g</mi></mrow></msub></mrow></math> , inter-group scattering matrix <math id="M5"><mrow><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow></msub></mrow></math> , and fission spectrum <math id="M6"><mrow><msub><mi>χ</mi><mi>g</mi></msub></mrow></math> need to be input to the NLSP3 code. Therefore, only these five grouped cross-sections are calculated in the MC global homogenization. It should be stated that it is not easy to obtain the anisotropic cross-sections during the MC global homogenization. To reduce the computational cost, only isotropic cross-sections are considered. The grouped flux distribution and grouped reaction rates are tallied to calculate the grouped cross-sections. The traditional volume-weighted method is employed in which the MGCS is defined as <math display="block" id="M7"><mrow><msubsup><mrow><mover accent="true"><mi>Σ</mi><mo stretchy="true">¯</mo></mover></mrow><mrow><mi>x</mi><mo>,</mo><mi>g</mi></mrow><mi>i</mi></msubsup><mo>≡</mo><mfrac><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><msub><mi>V</mi><mi>i</mi></msub></mrow></msub><mrow><msub><mi>Σ</mi><mrow><mi>x</mi><mo>,</mo><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><msub><mi>V</mi><mi>i</mi></msub></mrow></msub><mrow><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle></mrow></mfrac><mi>x</mi><mo>=</mo><mi>a</mi><mo>,</mo><mi>f</mi><mo>,</mo><mi>s</mi><mo>,</mo><mo>⋅</mo><mo>⋅</mo><mo>⋅</mo><mtext> </mtext><mi>g</mi><mo>=</mo><mn>1</mn><mo>,</mo><mo>⋅</mo><mo>⋅</mo><mo>⋅</mo><mo>,</mo><mi>G</mi><mo>,</mo></mrow></math> where r is a spatial variable, φ g ( r ) is the global flux distribution, and the subscript i in V i labels the different spatial regions. The track length is used to tally the flux distribution and reaction rates. The flux is calculated as <math display="block" id="M8"><mrow><msub><mi>ϕ</mi><mi>g</mi></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>E</mi><mi>g</mi></msub></mrow><mrow><msub><mi>E</mi><mrow><mi>g</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></msubsup><mrow><mtext>d</mtext><mi>E</mi></mrow></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mtext>d</mtext><mi>V</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><mi>W</mi><mo>⋅</mo><mi>T</mi><msubsup><mi>L</mi><mi>V</mi><mi>i</mi></msubsup><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></mrow></mstyle></mrow><mrow><mi>V</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msubsup><mi>W</mi><mn>0</mn><mi>i</mi></msubsup></mrow></mstyle></mrow></mfrac><mo>,</mo></mrow></math> and the reaction rate as <math display="block" id="M9"><mrow><msub><mi>R</mi><mi>g</mi></msub><mo>=</mo><msub><mi>Σ</mi><mi>g</mi></msub><msub><mi>ϕ</mi><mi>g</mi></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>E</mi><mi>g</mi></msub></mrow><mrow><msub><mi>E</mi><mrow><mi>g</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></msubsup><mrow><mtext>d</mtext><mi>E</mi></mrow></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mtext>d</mtext><mi>V</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><mi>W</mi><mo>⋅</mo><mi>T</mi><msubsup><mi>L</mi><mi>V</mi><mi>i</mi></msubsup><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo><mo>⋅</mo><mi>Σ</mi><mo stretchy="false">(</mo><mi>r</mi><mo>,</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></mrow></mstyle></mrow><mrow><mi>V</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msubsup><mi>W</mi><mn>0</mn><mi>i</mi></msubsup></mrow></mstyle></mrow></mfrac><mo>.</mo></mrow></math> The homogenized cross-section is therefore <math display="block" id="M10"><mrow><msub><mi>Σ</mi><mi>g</mi></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>E</mi><mi>g</mi></msub></mrow><mrow><msub><mi>E</mi><mrow><mi>g</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></msubsup><mrow><mtext>d</mtext><mi>E</mi></mrow></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mtext>d</mtext><mi>V</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><mi>W</mi><mo>⋅</mo><mi>T</mi><msubsup><mi>L</mi><mi>V</mi><mi>i</mi></msubsup><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo><mo>⋅</mo><mi>Σ</mi><mo stretchy="false">(</mo><mi>r</mi><mo>,</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>E</mi><mi>g</mi></msub></mrow><mrow><msub><mi>E</mi><mrow><mi>g</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></msubsup><mrow><mtext>d</mtext><mi>E</mi></mrow></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mtext>d</mtext><mi>V</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><mi>W</mi><mo>⋅</mo><mi>T</mi><msubsup><mi>L</mi><mi>V</mi><mi>i</mi></msubsup><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></mrow></mstyle></mrow></mfrac><mo>,</mo></mrow></math> where E g and E g -1 are the lower and upper bounds of the g th energy interval, respectively; <math id="M11"><mrow><mi>W</mi><mo>⋅</mo><mi>T</mi><msubsup><mi>L</mi><mtext>V</mtext><mi>i</mi></msubsup></mrow></math> is the product of the weight and track length of the i th neutron in region V ; and N is the total number of simulated neutrons. The calculation of the fission energy spectrum is relatively simple. The fission energy spectrum is calculated with reference to the number of fission neutrons produced in each energy group. The number of fission neutrons produced is then normalized to obtain the energy-dependent probability of neutron production in the energy interval E . The fission energy spectrum is equal to the energy-dependent probability of neutron production, that is, <math display="block" id="M12"><mrow><mi>χ</mi><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo><mo>=</mo><mi>p</mi><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where p ( E ) is the energy-dependent probability of neutron production. The process of calculating the fission matrix is the same as that for calculating the fission energy spectrum. The neutron energies before and after scattering are tallied to calculate the scattering probability: <math display="block" id="M13"><mrow><msub><mi>P</mi><mrow><mi>g</mi><mo>→</mo><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>E</mi><msup><mi>g</mi><mo>′</mo></msup></msub></mrow><mrow><msub><mi>E</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>−</mo><mn>1</mn></mrow></msub></mrow></msubsup><mrow><mtext>d</mtext><msup><mi>E</mi><mo>′</mo></msup></mrow></mrow></mstyle><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>E</mi><mi>g</mi></msub></mrow><mrow><msub><mi>E</mi><mrow><mi>g</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></msubsup><mrow><mtext>d</mtext><mi>E</mi></mrow></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mi>ϕ</mi><mo stretchy="false">(</mo><mi>r</mi><mo>,</mo><mi>E</mi><mo stretchy="false">)</mo><msub><mi>Σ</mi><mtext>s</mtext></msub><mo stretchy="false">(</mo><mi>r</mi><mo>,</mo><mi>E</mi><mo>→</mo><msup><mi>E</mi><mo>′</mo></msup><mo stretchy="false">)</mo><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>E</mi><mi>g</mi></msub></mrow><mrow><msub><mi>E</mi><mrow><mi>g</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></msubsup><mrow><mtext>d</mtext><mi>E</mi></mrow></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mi>ϕ</mi><mo stretchy="false">(</mo><mi>r</mi><mo>,</mo><mi>E</mi><mo stretchy="false">)</mo><msub><mi>Σ</mi><mtext>s</mtext></msub><mo stretchy="false">(</mo><mi>r</mi><mo>,</mo><mi>E</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle></mrow></mfrac><mo>.</mo></mrow></math> Using the scattering probability, the inter-group transfer cross-section is calculated as <math display="block" id="M14"><mrow><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>→</mo><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mo>=</mo><msub><mi>P</mi><mrow><mi>g</mi><mo>→</mo><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mo>⋅</mo><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi></mrow></msub><mo>.</mo></mrow></math> The MGCS is calculated using Equations (3)–(8) during the MC global homogenization. The MGCS of all material regions required by the NLSP3 code can thus be determined simultaneously in a single fixed-source calculation for subsequent use in the SP3 global calculation. Solution of the SP3 Equation 2 The SP3 equation is solved to obtain the global flux distribution. The SP3 equation is adopted because of the following reasons: (1) Compared with the diffusion and CMFD equations, the SP3 equation has higher accuracy [ 23 23 ]. (2) Solving the SP3 equation takes approximately the same amount of time as solving the diffusion equation [ 24 24 ]. (3) The SP3 equation is easier to apply to larger systems because it requires less memory and computing time compared to the SN equation. (4) To date, there has been no research on coupling the SP3 theory and the MC method for GVR. For academic interest, we would like to make an attempt to investigate this topic. The SPN theory [ 13 13 , 14 14 ] is a simplified spherical harmonics method for solving the neutron transport equation. It is widely studied because it has great potential for next-generation reactor physics calculations. There are two variants of the SPN theory, namely, the traditional and the rigorous SPN theories. In the traditional SPN theory, the one-dimensional governing equation is first found from a one-dimensional model, and the one-dimensional operator is subsequently replaced with a three-dimensional operator to obtain the three-dimensional governing equation. The traditional SPN theory lacks a rigorous mathematical basis. The rigorous SPN theory is derived based on the variational method and the concept of grouped angular flux. The angular flux is clearly defined in the rigorous SPN theory so that the physical fundamentals are completely considered in the treatment of the boundary conditions. These two SPN theories share the same SPN equation but use different boundary conditions. Although a complete mathematical derivation is provided in the rigorous SPN, the solution of the rigorous SPN equation suffers from numerical instability because of its incompatibility with the traditional transverse integral technique [ 25 25 ]. The traditional SPN equation is therefore solved in the SP3-coupled GVR method. The “ N ” in “SPN” represents the order of the polynomial expansion. When N is equal to 3, the SP3 equations are obtained as: <math display="block" id="M15"><mrow><mo>−</mo><msubsup><mi>D</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><msup><mo>∇</mo><mn>2</mn></msup><msubsup><mi>ϕ</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>+</mo><msubsup><mi>Σ</mi><mrow><mi>r</mi><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><msubsup><mi>ϕ</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><mn>2</mn><msubsup><mi>Σ</mi><mrow><mi>r</mi><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><msubsup><mi>ϕ</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>=</mo><msubsup><mi>S</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> <math display="block" id="M16"><mrow><mo>−</mo><msubsup><mi>D</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><msup><mo>∇</mo><mn>2</mn></msup><msubsup><mi>ϕ</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>+</mo><msubsup><mi>Σ</mi><mrow><mi>r</mi><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><msubsup><mi>ϕ</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><mfrac><mn>2</mn><mn>5</mn></mfrac><msubsup><mi>Σ</mi><mrow><mi>r</mi><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><msubsup><mi>ϕ</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>=</mo><mo>−</mo><mfrac><mn>2</mn><mn>5</mn></mfrac><msubsup><mi>S</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where the superscript k is the nodal index, the subscript g is the energy-interval index, <math id="M17"><mrow><msubsup><mi>ϕ</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup></mrow></math> is the zeroth-order scalar flux at the k th node in the g th energy interval, <math id="M18"><mrow><msubsup><mi>ϕ</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup></mrow></math> is the second-order scalar flux at the k th node in the g th energy interval, and <math display="block" id="M19"><mrow><msubsup><mi>D</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo>=</mo><mfrac><mn>1</mn><mrow><mn>3</mn><msubsup><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M20"><mrow><msubsup><mi>D</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo>=</mo><mfrac><mn>9</mn><mrow><mn>35</mn><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M21"><mrow><msubsup><mi>Σ</mi><mrow><mi>r</mi><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo>=</mo><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo>−</mo><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><msup><mi>g</mi><mo>′</mo></msup></mrow><mi>k</mi></msubsup><mo>,</mo></mrow></math> <math display="block" id="M22"><mrow><msubsup><mi>Σ</mi><mrow><mi>r</mi><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo>=</mo><mfrac><mn>9</mn><mn>5</mn></mfrac><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo>−</mo><mfrac><mn>4</mn><mn>5</mn></mfrac><msubsup><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><msup><mi>g</mi><mo>′</mo></msup></mrow><mi>k</mi></msubsup><mo>.</mo></mrow></math> The source is <math display="block" id="M23"><mrow><msubsup><mi>S</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>=</mo><munderover><mi>Σ</mi><mtable columnalign="left"><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>≠</mo><mi>g</mi></mrow></mtd></mtr></mtable><mi>G</mi></munderover><msubsup><mi>Σ</mi><mrow><mn>0</mn><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow><mi>k</mi></msubsup><mrow><mo>[</mo><mrow><msubsup><mi>ϕ</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><mn>2</mn><msubsup><mi>ϕ</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>K</mi><mrow><mtext>ef</mtext><mi>f</mi></mrow></msub></mrow></mfrac><msubsup><mi>χ</mi><mi>g</mi><mi>k</mi></msubsup><munderover><mi>Σ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mi>ν</mi><msubsup><mi>Σ</mi><mrow><mtext>f</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup></mrow><mi>k</mi></msubsup><mrow><mo>[</mo><mrow><msubsup><mi>ϕ</mi><mrow><mn>0</mn><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><mn>2</mn><msubsup><mi>ϕ</mi><mrow><mn>2</mn><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup></mrow><mi>k</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow><mo>.</mo></mrow></math> The boundary conditions [ 26 26 ] for the traditional SP3 equation are <math display="block" id="M24"><mrow><msubsup><mi>J</mi><mn>0</mn><mo>±</mo></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mn>4</mn></mfrac><msubsup><mi>ϕ</mi><mn>0</mn><mi>M</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>±</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mover accent="true"><mstyle mathvariant="bold" mathsize="normal"><mi>n</mi></mstyle><mo stretchy="true">⇀</mo></mover><mo>·</mo><msubsup><mi>J</mi><mn>0</mn><mi>M</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><mfrac><mn>3</mn><mrow><mn>16</mn></mrow></mfrac><msubsup><mi>ϕ</mi><mn>2</mn><mi>M</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> <math display="block" id="M25"><mrow><msubsup><mi>J</mi><mn>2</mn><mo>±</mo></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mn>21</mn></mrow><mrow><mn>80</mn></mrow></mfrac><msubsup><mi>ϕ</mi><mn>2</mn><mi>M</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>±</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mover accent="true"><mstyle mathvariant="bold" mathsize="normal"><mi>n</mi></mstyle><mo stretchy="true">⇀</mo></mover><mo>·</mo><msubsup><mi>J</mi><mn>2</mn><mi>M</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><mfrac><mn>3</mn><mrow><mn>80</mn></mrow></mfrac><msubsup><mi>ϕ</mi><mn>0</mn><mi>M</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where the superscript M represents the quantities at the node boundary, and <math id="M26"><mrow><msubsup><mi>J</mi><mn>0</mn><mo>±</mo></msubsup></mrow></math> and <math id="M27"><mrow><msubsup><mi>J</mi><mn>2</mn><mo>±</mo></msubsup></mrow></math> are the zeroth-order and second-order partial currents, respectively. The numerical instability problem was encountered when the SP3 equation was solved using the traditional nonlinear iterative method [ 27 27 ]. The zeroth-order flux is usually hundreds of times larger than the second-order flux. When the traditional coupling corrective relationship was used to calculate the corrective factors, there was also a hundred-fold difference between the zeroth-order and second-order corrective factors, which led to an ill-conditioned iteration matrix and convergence failure. The following new coupling corrective relationship is proposed to solve the numerical instability problem : <math display="block" id="M28"><mtable columnalign="left"><mtr><mtd><msubsup><mi>J</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><mtext>=</mtext><mo>−</mo><msubsup><mi>D</mi><mrow><mn>1</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mtext>FDM</mtext></mrow></msubsup><mo stretchy="false">(</mo><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">)</mo><mo>−</mo><msubsup><mi>D</mi><mrow><mn>1</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mi>N</mi><mi>O</mi><mi>D</mi></mrow></msubsup><mo stretchy="false">(</mo><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>           </mtext><mo>−</mo><msubsup><mi>D</mi><mrow><mn>3</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mtext>FDM</mtext></mrow></msubsup><mo stretchy="false">(</mo><msubsup><mi>f</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mi>f</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">)</mo><mo>−</mo><msubsup><mi>D</mi><mrow><mn>3</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mi>N</mi><mi>O</mi><mi>D</mi></mrow></msubsup><mo stretchy="false">(</mo><msubsup><mi>f</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msubsup><mi>f</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">)</mo><mo>,</mo></mtd></mtr></mtable></math> <math display="block" id="M29"><mtable columnalign="left"><mtr><mtd><msubsup><mi>J</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><mtext>=</mtext><mo>−</mo><msubsup><mi>D</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mtext>FDM</mtext></mrow></msubsup><mo stretchy="false">(</mo><msubsup><mi>f</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mi>f</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">)</mo><mo>−</mo><msubsup><mi>D</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mi>N</mi><mi>O</mi><mi>D</mi></mrow></msubsup><mo stretchy="false">(</mo><msubsup><mi>f</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msubsup><mi>f</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>2</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>           </mtext><mo>−</mo><msubsup><mi>D</mi><mrow><mn>4</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mtext>FDM</mtext></mrow></msubsup><mo stretchy="false">(</mo><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">)</mo><mo>−</mo><msubsup><mi>D</mi><mrow><mn>4</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mi>N</mi><mi>O</mi><mi>D</mi></mrow></msubsup><mo stretchy="false">(</mo><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mi>k</mi></msubsup><msubsup><mover accent="true"><mi>ϕ</mi><mo stretchy="true">¯</mo></mover><mrow><mn>0</mn><mo>,</mo><mi>g</mi></mrow><mi>k</mi></msubsup><mo stretchy="false">)</mo></mtd></mtr></mtable></math> where <math id="M30"><mrow><msubsup><mi>D</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mtext>FDM</mtext></mrow></msubsup></mrow></math> and <math id="M31"><mrow><msubsup><mi>D</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mi mathvariant="normal">FDM</mi></mrow></msubsup></mrow></math> are the zeroth-order and second-order pseudo-diffusion coefficients, respectively; <math id="M32"><mrow><msubsup><mi>D</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mtext>FDM</mtext></mrow></msubsup></mrow></math> and <math id="M33"><mrow><msubsup><mi>D</mi><mrow><mn>2</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>+</mo></mrow><mrow><mi>k</mi><mo>,</mo><mtext>FDM</mtext></mrow></msubsup></mrow></math> are the zeroth-order and second-order coupling corrective coefficients, respectively; and <math id="M34"><mrow><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup></mrow></math> and <math id="M35"><mrow><msubsup><mi>f</mi><mrow><mn>0</mn><mo>,</mo><mi>g</mi><mi>u</mi><mo>−</mo></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup></mrow></math> are the zeroth-order and second-order discontinuity factors, respectively. A new nonlinear iterative method [ 25 25 ] was proposed based on the above mentioned coupling corrective relationship, and the NLSP3 reactor core code was developed. Figure 2 Figure 2 illustrates the nonlinear iterative strategy. Using the NLSP3 code and global cross sections, the global flux distribution can be obtained easily. GVR for MC Criticality Calculation 2 The MC criticality calculation is composed of several inactive and active cycles. Each cycle is a fixed-source calculation with a cycle-estimated result. All the active cycles are averaged as follows to obtain the final result of the MC criticality calculation: <math display="block" id="M36"><mrow><mover accent="true"><mi>X</mi><mo stretchy="true">¯</mo></mover><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><mover accent="true"><mrow><msub><mi>Y</mi><mi>n</mi></msub></mrow><mo stretchy="true">¯</mo></mover></mrow></mstyle></mrow><mi>N</mi></mfrac><mo>,</mo></mrow></math> where the subscript n denotes the cycle index, <math id="M37"><mrow><mover accent="true"><mi>X</mi><mo stretchy="true">¯</mo></mover></mrow></math> the final result of the MC criticality calculation, N the total number of active cycles, and <math id="M38"><mrow><mover accent="true"><mrow><msub><mi>Y</mi><mi>n</mi></msub></mrow><mo stretchy="true">¯</mo></mover></mrow></math> the cycle-estimated result of the n th active cycle. <math id="M39"><mrow><mover accent="true"><mrow><msub><mi>Y</mi><mi>n</mi></msub></mrow><mo stretchy="true">¯</mo></mover></mrow></math> is given by <math display="block" id="M40"><mrow><mover accent="true"><mrow><msub><mi>Y</mi><mi>n</mi></msub></mrow><mo stretchy="true">¯</mo></mover><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></msubsup><mrow><msubsup><mi>Y</mi><mi>n</mi><mi>i</mi></msubsup></mrow></mstyle></mrow><mi>M</mi></mfrac><mo>,</mo></mrow></math> where <math id="M41"><mrow><msubsup><mi>Y</mi><mi>n</mi><mi>i</mi></msubsup></mrow></math> is the i th neutron history of the n th active cycle, and M is the number of fission sources in a single cycle. The standard deviation of the final result is calculated from the deviation between the individual active cycles: <math display="block" id="M42"><mrow><msub><mi>S</mi><mrow><mover accent="true"><mi>X</mi><mo stretchy="true">¯</mo></mover></mrow></msub><mo>=</mo><msqrt><mrow><mfrac><mrow><msup><mi>S</mi><mn>2</mn></msup></mrow><mi>N</mi></mfrac></mrow></msqrt><mo>=</mo><msqrt><mrow><mfrac><mrow><mover accent="true"><mrow><msup><mi>X</mi><mn>2</mn></msup></mrow><mo stretchy="true">¯</mo></mover><mo>−</mo><msup><mrow><mover accent="true"><mi>X</mi><mo stretchy="true">¯</mo></mover></mrow><mn>2</mn></msup></mrow><mrow><mi>N</mi><mo>−</mo><mn>1</mn></mrow></mfrac></mrow></msqrt><mo>,</mo></mrow></math> where <math display="block" id="M43"><mrow><mover accent="true"><mrow><msup><mi>X</mi><mn>2</mn></msup></mrow><mo stretchy="true">¯</mo></mover><mo>=</mo><mfrac><mrow><msup><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><mover accent="true"><mrow><msub><mi>Y</mi><mi>n</mi></msub></mrow><mo stretchy="true">¯</mo></mover></mrow></mstyle></mrow><mn>2</mn></msup></mrow><mi>N</mi></mfrac><mo>.</mo></mrow></math> Therefore, the standard deviation increases with the deviation between the individual active cycles. The deviation between the individual active cycles should therefore be reduced to reduce the standard deviation of the MC criticality calculation. As stated before, a fixed-source calculation is performed for each cycle, and the deviation between the individual cycles is reduced by various VR techniques. Suppose that all the active cycles are performed with the same weight windows. The individual cycle-estimated results will then be close to another. This leads to a reduction in the standard deviation. Therefore, the key to improving the MC criticality calculation is to obtain high-precision weight window parameters efficiently. In the SP3-coupled GVR method, only the forward SP3 calculation is performed. Because there is no known adjoint theory of the SP3 method, we cannot obtain the adjoint flux for the time being. We have therefore focused on GVR for the criticality calculations in the meantime. We believe that the forward flux distribution is sufficient for determining the weight window parameters in place of the fixed-source calculation. Therefore, we use the forward flux distribution to approximate the weight window parameters. The global flux distribution determined by the SP3 equation is denoted as <math id="M44"><mrow><msubsup><mi>ϕ</mi><mi>g</mi><mi>k</mi></msubsup></mrow></math> , where the superscript k denotes the spatial region and the subscript g denotes the energy interval. The SP3-coupled method can only bias the spatial position of the neutrons. Therefore, the independent energy variable should be eliminated by accumulating the energy-dependent flux in the spatial region: <math display="block" id="M45"><mrow><msub><mi>ϕ</mi><mi>k</mi></msub><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><msubsup><mi>ϕ</mi><mi>g</mi><mi>k</mi></msubsup></mrow></mstyle><mo>,</mo></mrow></math> where φ k is the space-dependent flux distribution, which is normalized before it is used to approximate the global weight windows, and <math display="block" id="M46"><mrow><msub><mover accent="true"><mi>ϕ</mi><mo>¯</mo></mover><mi>k</mi></msub><mo>=</mo><msub><mi>ϕ</mi><mi>k</mi></msub><mo>/</mo><mtext>max</mtext><mo stretchy="false">(</mo><msub><mi>ϕ</mi><mi>k</mi></msub><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where max( φ k ) is the largest space-dependent flux. Using the normalized global flux distribution, the global weight windows are determined as <math display="block" id="M47"><mrow><msub><mi>W</mi><mrow><mtext>s</mtext><mo>,</mo><mi>k</mi></mrow></msub><mo>=</mo><msub><mover accent="true"><mi>ϕ</mi><mo>¯</mo></mover><mi>k</mi></msub><mo>,</mo></mrow></math> where W s,k is the global space-dependent survival weight window in the k th space region. Using Eq. (20), the weight windows of all the spatial regions can be determined, following which the ceiling and floor weight windows can be determined using the following formulas: <math display="block" id="M48"><mrow><msub><mi>W</mi><mrow><mtext>c</mtext><mo>,</mo><mi>k</mi></mrow></msub><mo>=</mo><mfrac><mrow><mn>5</mn><msub><mi>W</mi><mrow><mtext>s</mtext><mo>,</mo><mi>k</mi></mrow></msub></mrow><mn>3</mn></mfrac><mo>,</mo></mrow></math> and <math display="block" id="M49"><mrow><msub><mi>W</mi><mrow><mtext>f</mtext><mo>,</mo><mi>k</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>W</mi><mrow><mtext>s</mtext><mo>,</mo><mi>k</mi></mrow></msub></mrow><mn>3</mn></mfrac><mo>.</mo></mrow></math> where W c ,k and W f ,k are the ceiling and floor weight windows in the k th space region, respectively. The standard weight window with the weight windows of all the spatial regions is used for GVR. NUMERICAL RESULTS 1 The SP3-coupled GVR method was tested on the C5G7 benchmark [ 28 28 ]. In the original C5G7 model, GVR can be easily achieved because the water shield layer is not sufficiently thick to pose a deep penetration problem. To better demonstrate the advantages of the SP3-coupled method for GVR, we artificially increased the thickness of the shield layer. The materials and geometry of the reactor core were derived from the C5G7 benchmark, and the core region was surrounded by a water shield layer, which was 64.362 cm thick and three times thicker than the side of the core assembly side as shown in Fig. 3 Fig. 3 . The superimposed mesh geometry was used for SP3-coupled GVR simulation. The C5G7 reactor core had a 4×4 assembly layout and the assemblies was surrounded by a water shield that was three times thicker as the assembly. The model was divided into a total of 10×10=100 material mesh in the modified C5G7 benchmark. The MGCS of the 100 mesh were determined through the MC global homogenization over the four energy intervals of ［0, 6.25×10 -7 ］, [6.25×10 -7 , 5.53×10 -3 ], [5.53×10 -3 , 8.21×10 -1 ], [8.21×10 -1 , 20] MeV [ 29 29 ]. A fixed-source calculation was performed for the MC global homogenization. The choice of the number of source neutrons per material mesh is a challenging issue. A greater number of source neutrons for the fixed-source calculation will increase the precision of the obtained cross sections but will also require more computation. To achieve a balance between efficiency and accuracy, the number of source neutrons per material mesh was chosen as 10,000 based on our experience. Therefore, a total of 1,000,000 source neutrons were simulated for the MC global homogenization. The SP3 solution only provided the weight windows, and the final result was given by the MC criticality calculation. Thus, even if there was an error in the MGCS results, the accuracy of the final result was not affected. This is why the number of source neutrons per material mesh was empirically determined for the solution of the MGCS. The RMC code called the SP3 solver for the global calculation with the cross-sections of the 10×10 material mesh. Each material mesh was divided into a total of 4×4=16 nodes with each node being 5 cm × 5 cm in size. Therefore, a total of 40×40=1600 nodes were used to solve the SP3 equation. After solving the SP3 equation, the flux distribution in the 1600 nodes was obtained. The three convergence criteria for the inner, outer, and k eff iterations in the NLSP3 code were set to 10 -3 , 10 -3 , and 10 -4 , respectively. The calculations were performed with 24 threads, and the calculation time taken for the MC global homogenization and SP3 solver was 1.14 minutes. Given the global flux distribution, the global weight windows can be determined using Eqs. (20) – (22). The MC criticality calculation was then performed using the global weight windows. There were 10,000 neutrons per cycle, 50 inactive cycles, and 450 active cycles in the SP3-coupled GVR criticality calculation. 15.67 minutes were required to complete the MC global homogenization, SP3 solver solution, and GVR criticality calculation. For comparison, the standard MC simulation method was also performed to calculate the flux distribution of the 1600 nodes. Because GVR was not performed in the standard simulation method, we call this the direct simulation. There were 1000,000 neutrons per cycle, 50 inactive cycles, and 450 active cycles in the simulation. The calculation time was 151.80 minutes. The distributions of the relative errors in the SP3-coupled GVR and the direct simulations are compared in Fig. 4 Fig. 4 . The distributions of the total flux in the SP3-coupled GVR and the direct simulations are compared in Figure 5 Figure 5 . As shown in Fig. 5 Fig. 5 , the distributions of the total flux obtained by the direct and the SP3-coupled GVR simulations are in good agreement. The deviation was calculated to be less than three times the standard deviation. Meanwhile, it can be seen from Figure 4 Figure 4 that when the modified C5G7 benchmark was simulated directly, the relative errors increased with the distance away from the reactor core. Compared with the direct simulation, the global variance distribution in the SP3-coupled GVR simulation was flattened in the entire system. High accuracy was achieved in both the reactor core and the shielding layer. Therefore, the SP3-coupled GVR method can improve the efficiency of MC criticality calculations. To better quantify its GVR performance, three quantities were compared: 1. The average figure of merit (AV.FOM), which is defined as <math display="block" id="M50"><mrow><mi>A</mi><mi>V</mi><mo>.</mo><mi>F</mi><mi>O</mi><mi>M</mi><mo>=</mo><mfrac><mi>N</mi><mrow><mi>T</mi><mo>⋅</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msubsup><mrow><mi>Re</mi></mrow><mi>i</mi><mn>2</mn></msubsup></mrow></mstyle></mrow></mfrac><mo>.</mo></mrow></math> 2. The standard deviation of the relative deviation, which is defined as <math display="block" id="M51"><mrow><msub><mi>σ</mi><mrow><mi>Re</mi></mrow></msub><mo>=</mo><msqrt><mrow><mfrac><mn>1</mn><mi>N</mi></mfrac><mo>⋅</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msubsup><mrow><mi>Re</mi></mrow><mi>i</mi><mn>2</mn></msubsup></mrow></mstyle><mo>−</mo><mfrac><mn>1</mn><mrow><msup><mi>N</mi><mn>2</mn></msup></mrow></mfrac><mo>⋅</mo><msup><mrow><mo stretchy="false">(</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msubsup><mrow><mi>Re</mi></mrow><mi>i</mi><mn>2</mn></msubsup></mrow></mstyle><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow></msqrt><mo>.</mo></mrow></math> 3. The ratio of the maximum to the minimum variances (σ max /σ min ), i.e., <math display="block" id="M52"><mrow><msub><mi>σ</mi><mrow><mi>max</mi></mrow></msub><mo>/</mo><msub><mi>σ</mi><mrow><mi>min</mi></mrow></msub><mo>=</mo><mfrac><mrow><mi>max</mi><mo stretchy="false">(</mo><msub><mrow><mi>Re</mi></mrow><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mi>min</mi><mo stretchy="false">(</mo><msub><mrow><mi>Re</mi></mrow><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mfrac><mo>.</mo></mrow></math> The quantitative results for the direct and SP3-coupled GVR simulations are compared in Table 1 Table 1 . 1 1 Quantities 1 1 Direct simulation 1 1 SP3-coupled GVR 1 1 AV.FOM 1 1 1.1364×10 1 1 1 3.8394×10 2 1 1 σ.Re 1 1 2.4070×10 -2 1 1 1.2893×10 -2 1 1 σ max /σ min 1 1 6.1998×10 2 1 1 1.2196×10 1 1 1 Times (min) 1 1 151.80 1 1 15.67 The results show that, compared with the direct simulation, the average figure of merit (AV. FOM) was improved 33 folds in the SP3-coupled GVR method. Meanwhile, the ratio of the maximum to minimum variances was reduced by a factor of five. The value of the SP3-coupled GVR method for GVR and its ability to improve the efficiency of the MC criticality calculation are hence clearly demonstrated. It has been argued that it is unnecessary to couple the SP3 method with the MC method for GVR, and that the diffusion method is sufficient. To demonstrate the advantage of coupling the SP3 method with the MC method for GVR over coupling the diffusion method, we developed a diffusion-coupled GVR method based on the RMC code. The flowcharts of the diffusion-coupled GVR and the SP3-coupled GVR methods are almost identical except for the solver used to solve the global flux distribution. In the SP3-coupled GVR method, the SP3 equation is solved, while in the diffusion-coupled method, the diffusion equation is solved. To solve the diffusion equation, the diffusion coefficient is required. Therefore, the total cross-section <math id="M53"><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow></msub></mrow></math> used in the SP3 equation was converted into the diffusion coefficient D g . Calculations were then performed on the same model using the diffusion-coupled GVR method with the same calculation parameters. The obtained AV.FOM from the results of the diffusion-coupled method was 136.81. As shown in Table 1 Table 1 , the AV.FOM of the result of the SP3-coupled method was 383.94, which is 2.81 times that of the diffusion-coupled method. Therefore, the SP3-coupled GVR method showed better GVR performance than the diffusion-coupled GVR method. We believe that the SP3 method has great potential to be coupled with the MC method for accelerated calculations and GVR. DISCUSSION AND CONCLUSION 1 An automatic process based on MC global homogenization, the SP3 global solver, and GVR criticality calculation was developed in the SP3-coupled GVR method. The numerical results show that this method can improve the efficiency of the MC criticality calculation. The SP3-coupled method proposed in this paper was developed for MC criticality calculations. It therefore superficially resembles the traditional two-step reactor physical calculation framework. However, the SP3-coupled method differs from the traditional framework in three aspects: (1) global homogenization is performed, instead of assembly homogenization; (2) the homogenization process is completed during the fixed-source calculation, and not during the criticality calculation; (3) global MC homogenization and the global reactor calculation are achieved in a one-step process instead of a two-step process. Moreover, only the weight window parameters for the MC criticality calculation are output in the SP3-coupled method. The final results are obtained from the MC calculation rather than the global calculation of the NLSP3 code. The SP3-coupled method has similarities and differences with other hybrid MC methods, such as the fission matrix and the response matrix methods. They share similar calculation processes in which some local information of the system is recorded in advance and used to build a global formula, and the global information is obtained by solving the global formula. In the fission matrix method, the fission coupling relationship between different regions is recorded first. In the response matrix method, the relationship between the currents and sources in different regions is recorded first. In the SP3-coupled method, the homogenized cross-sections in different regions are recorded first. The matrix relationship is solved in the fission matrix and response matrix methods to obtain global information. In contrast, in the SP3-coupled method, the neutron transport equation is solved. We therefore believe that the SP3-coupled method has advantages in terms of numerical stability, development potential, and application scope. In this study, the traditional SP3 equation is solved in the SP3-coupled method and the forward flux is used to approximate the weight windows. To further improve its performance for GVR, the rigorous SPN theory proposed by Prof. Yung-An Chao and the adjoint method for the SPN theory should be studied to obtain better weight window parameters. Because the SP3-coupled method is helpful for the GVR of the MC criticality calculation, we hope to apply it to reactor shielding calculations in the future. We are attempting to perform reactor shielding design using MC criticality and burnup calculations, and hope to have an accurate radiation source for reactor shielding calculations to analyze the influence of the reactor dynamics on the reactor shielding. REFERENCE Ma Y. , Liu S. , Luo Z. et al. , RMC/CTF multiphysics solutions to VERA core physics benchmark problem 9 . Ann. Nucl. Energy 133 , 837 - 852 ( 2019 ). doi: 10.1016/j.anucene.2019.07.033 http://doi.org/10.1016/j.anucene.2019.07.033 . RMC/CTF multiphysics solutions to VERA core physics benchmark problem 9 Kotlyar D. , Shwageraus E. , Numerical stable Monte Carlo-burnup-thermal-hydraulic coupling schemes . Ann. Nucl. Energy 63 , 371 - 381 ( 2014 ). doi: 10.1016/j.anucene.2013.08.016 http://doi.org/10.1016/j.anucene.2013.08.016 Numerical stable Monte Carlo-burnup-thermal-hydraulic coupling schemes Guo X.Y. , Shang X.T. , Song J. et al . Kinetic methods in Monte Carlo code RMC and its implementation to C5G7-TD benchmark . Ann. Nucl. Energy 151 , 107864 , ( 2021 ). doi: 10.1016/j.anucene.2020.107864 http://doi.org/10.1016/j.anucene.2020.107864 Kinetic methods in Monte Carlo code RMC and its implementation to C5G7-TD benchmark Zheng Z. , Wu H.C. , Cao L.Z. et al. , A deterministic and Monte Carlo coupling method for PWR cavity radiation streaming calculation . Ann. Nucl. Energy 63 , 162 - 171 ( 2014 ). doi: 10.1016/j.anucene.2013.07.040 http://doi.org/10.1016/j.anucene.2013.07.040 A deterministic and Monte Carlo coupling method for PWR cavity radiation streaming calculation Zheng Z. , Mei Q.L. , Deng L. , Study on variance reduction technique based on adjoint Discrete Ordinate method . Ann. Nucl. Energy 112 , 374 - 382 ( 2018 ). doi: 10.1016/j.anucene.2017.10.028 http://doi.org/10.1016/j.anucene.2017.10.028 Study on variance reduction technique based on adjoint Discrete Ordinate method Pan Q.Q. , Rao J.J. , Zhou Y.Y. et al. , Improved adaptive variance reduction algorithm based on RMC code for deep penetration problems . Ann. Nucl. Energy 137 , 107113 ( 2020 ). doi: 10.1016/j.anucene.2019.107113 http://doi.org/10.1016/j.anucene.2019.107113 Improved adaptive variance reduction algorithm based on RMC code for deep penetration problems Wagner J. , Peplow D. , Mosher S. , FW-CADIS method for global and regional variance reduction of Monte Carlo radiation transport calculations . Nucl. Sci. Eng. 176 , 37 - 57 ( 2014 ). doi: 10.13182/NSE12-33 http://doi.org/10.13182/NSE12-33 FW-CADIS method for global and regional variance reduction of Monte Carlo radiation transport calculations Cooper M. , Larsen E. , Automated weight windows for global Monte Carlo particle transport calculations . Nucl. Sci. Eng. 137 , 1 - 13 ( 2001 ). Automated weight windows for global Monte Carlo particle transport calculations Hendricks J. , A code-generated Monte Carlo importance function . T. Am. Nucl. Soc. 41 , 307 ( 1982 ). A code-generated Monte Carlo importance function Kelly D. , Sutton T. , Wilson S. . MC21 Analysis of the nuclear energy agency Monte Carlo performance benchmark problem . Proc. Advances in Reactor Physics—Linking Research, Industry, and Education, (PHYSOR 2012) , Knoxville, Tennessee , April 15-20 , American Nuclear Society , 2012 . MC21 Analysis of the nuclear energy agency Monte Carlo performance benchmark problem Shangguan D.H. , Li G. , Zhang B.Y. et al. , Uniform tally density based strategy for global tallying in Monte Carlo criticality calculation . Nucl. Sci. Eng. 182 ( 4 ), 555 - 562 ( 2016 ). doi: 10.13182/NSE15-32 http://doi.org/10.13182/NSE15-32 Uniform tally density based strategy for global tallying in Monte Carlo criticality calculation Pan Q.Q. , Wang K. , Uniform variance method for accelerated Monte Carlo criticality calculation . Prog. Nucl. Energy 139 , 103858 ( 2021 ). doi: 10.1016/j.pnucene.2021.103858 http://doi.org/10.1016/j.pnucene.2021.103858 Uniform variance method for accelerated Monte Carlo criticality calculation Chao Y. , A new and rigorous SPN theory – Part III: A succinct summary of the GSPN theory, the P3 equivalent GSP3 and implementation issues . Ann. Nucl. Energy 119 , 310 - 321 ( 2019 ). doi: 10.1016/j.anucene.2018.04.029 http://doi.org/10.1016/j.anucene.2018.04.029 A new and rigorous SPN theory – Part III: A succinct summary of the GSPN theory, the P3 equivalent GSP3 and implementation issues Chao Y. , Peng L. , Tang C. , A new and rigorous SPN theory – Part IV: Numerical qualification of GSP3 and the generalized transverse integration nodal method . Ann. Nucl. Energy 149 , 107768 ( 2020 ). doi: 10.1016/j.anucene.2020.107768 http://doi.org/10.1016/j.anucene.2020.107768 A new and rigorous SPN theory – Part IV: Numerical qualification of GSP3 and the generalized transverse integration nodal method Wang K. , Li Z.G. , She D. et al. , RMC – A Monte Carlo code for reactor core analysis . Ann. Nucl. Energy 82 , 121 - 129 ( 2015 ). doi: 10.1016/j.anucene.2014.08.048 http://doi.org/10.1016/j.anucene.2014.08.048 RMC – A Monte Carlo code for reactor core analysis Carney S. , Brown F. , Kiedrowski B. et al. , Theory and application of the fission matrix method for continuous-energy Monte Carlo . Ann. Nucl. Energy 73 , 423 - 431 ( 2014 ). doi: 10.1016/j.anucene.2014.07.020 http://doi.org/10.1016/j.anucene.2014.07.020 Theory and application of the fission matrix method for continuous-energy Monte Carlo Leppanen J. , Response matrix method–based importance solver and variance reduction scheme in the serpent 2 Monte Carlo code . Nucl. Technol. 205 ( 11 ), 1416 - 1432 ( 2019 ). doi: 10.1080/00295450.2019.1603710 http://doi.org/10.1080/00295450.2019.1603710 Response matrix method–based importance solver and variance reduction scheme in the serpent 2 Monte Carlo code Pan Q.Q. , Hao L. , Li D. et al. , A new nonlinear iterative method for SPN theory . Ann. Nucl. Energy 110 , 920 - 927 ( 2017 ). doi: 10.1016/j.anucene.2017.07.030 http://doi.org/10.1016/j.anucene.2017.07.030 A new nonlinear iterative method for SPN theory HidekiI T. , Yasushi M. , Accuracy of interpolation methods for resonance self-shielding factors . Nucl. Sci. Technol. 18 ( 2 ), 152 - 161 ( 1981 ). doi: 10.1080/18811248.1981.9733236I http://doi.org/10.1080/18811248.1981.9733236I . Accuracy of interpolation methods for resonance self-shielding factors Petrovic , Benoist P. . BN Theory: Advances and New Models for Neutron Leakage Calculation . New York : Plenum Press , 223 - 281 , 1996 . doi: 10.1093/molehr/1.1.36 http://doi.org/10.1093/molehr/1.1.36 BN Theory: Advances and New Models for Neutron Leakage Calculation Pan Q.Q. , Wang K. , Acceleration method of fission source convergence based on RMC code . Ann. Nucl. Energy 52 ( 7 ), 1347 - 1354 ( 2020 ). doi: 10.1016/j.net.2019.12.007 http://doi.org/10.1016/j.net.2019.12.007 Acceleration method of fission source convergence based on RMC code Duderstadt J. , Hamilton L. . Nuclear reactor theory . John Wiley & Sons , New York , 1976 . Nuclear reactor theory Beckert C. , Grundmann U. , Development and verification of a nodal approach for solving the multigroup SP3 equations . Ann. Nucl. Energy 35 , 75 - 86 ( 2008 ). doi: 10.1016/j.anucene.2007.05.014 http://doi.org/10.1016/j.anucene.2007.05.014 Development and verification of a nodal approach for solving the multigroup SP3 equations Smith K.S. , Multidimensional nodal transport using the simplified PL method . T. Am. Nucl. Soc. 52 , 427 ( 1986 ). Multidimensional nodal transport using the simplified PL method Pan Q.Q. , Lu H.L. , Li D.S. et al. , The rigorous SP3 theory and study on its numerical verification . Proceedings of 25th International Conference on Nuclear Engineering (ICONE 25) , 2017 . doi: 10.1115/ICONE25-67309 http://doi.org/10.1115/ICONE25-67309 The rigorous SP3 theory and study on its numerical verification Brantley P. , Larsen E. , The simplified P3 approximation . Nucl. Sci. Eng. 134 , 1 - 21 ( 2000 ). doi: 10.1016/S0168-583X(99)00603-5 http://doi.org/10.1016/S0168-583X(99)00603-5 The simplified P3 approximation Smith K. , Nodal method storage reduction by nonlinear iteration . T. Am. Nucl. Soc. 44 , 265 ( 1983 ). Nodal method storage reduction by nonlinear iteration Smith M. , Lewis E. , Na B. , Benchmark on deterministic transport calculations without spatial homogenization –a 2d/3d MOX fuel assembly benchmark, Tech. Rep. NEA/NSC/DOC 16, OECD/NEA , ( 2003 ). doi: 10.1016/j.pnueene.2004.09.003 http://doi.org/10.1016/j.pnueene.2004.09.003 Benchmark on deterministic transport calculations without spatial homogenization –a 2d/3d MOX fuel assembly benchmark, Tech. Rep. NEA/NSC/DOC 16, OECD/NEA Boyd W. , Reactor multi-group cross section generation for fine mesh deterministic neutron transport simulations. PhD Thesis , Massachusetts Institute of Technology , 2017 . Reactor multi-group cross section generation for fine mesh deterministic neutron transport simulations. PhD Thesis All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Qing-Quan Pan, Xiao-Jing Liu and Kan Wang. The first draft of the manuscript was written by Qing-Quan Pan and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-021-00973-0 This work was supported by the Shanghai Sailing Program, China (No. 21YF1421100) , and the Startup Fund for Youngman Research at SJTU. 2021-08-09 2021-09-26 2021-10-14 2021-11-16 2021-11-01 2021-12-06 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2021 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-021-00973-0 http://dx.doi.org/10.1007/s41365-021-00973-0 . 2021 Nuclear Science and Techniques 11 32