Equivalent low-order angular flux nonlinear finite difference equation of MOC transport calculation

Equivalent low-order angular flux nonlinear finite difference equation of MOC transport calculation 1 2 first-author https://orcid.org/0000-0002-6849-5156 1 corresp Corresponding author. E-mail address: haochen.heu@163.com Corresponding author. E-mail address: haochen.heu@163.com haochen.heu@163.com 3 https://orcid.org/0000-0002-4280-4117 Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin 150001, China. Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin 150001, China. Institute of Nuclear and New Energy Technology (INET), Tsinghua University, Beijing 100084, China Institute of Nuclear and New Energy Technology (INET), Tsinghua University, Beijing 100084, China Purdue University, West Lafayette, Indiana 47907, United States Purdue University, West Lafayette, Indiana 47907, United States The key issue in accelerating method of characteristics (MOC) transport calculations is in obtaining a completely equivalent low-order neutron transport or diffusion equation. Herein, an equivalent low-order angular flux nonlinear finite difference equation is proposed for MOC transport calculations. This method comprises three essential features: (1) the even parity discrete ordinates method is used to build a low-order angular flux nonlinear finite difference equation, and different boundary condition treatments are proposed; (2) two new defined factors, i.e., the even parity discontinuity factor and odd parity discontinuity factor, are strictly defined to achieve equivalence between the low-order angular flux nonlinear finite difference method and MOC transport calculation; (3) the energy group and angle are decoupled to construct a symmetric linear system that is much easier to solve. The equivalence of this low-order angular flux nonlinear finite difference equation is analyzed for two-dimensional (2D) pin, 2D assembly, and 2D C5G7 benchmark problems. Numerical results demonstrate that a low-order angular flux nonlinear finite difference equation that is completely equivalent to the pin-resolved transport equation is established. Key words: Angular flux Equivalence Even parity discrete ordinates method Nonlinear finite difference 1 Introduction 1 Significant advances in high-performance computing clusters (HPCs) have enabled three-dimensional (3D) whole-core high-fidelity heterogeneous transport simulations at the subpin level. This is becoming increasingly important for improving the reactor performance and safety as well as providing an efficient tool for more complex and novel designs. A numerical simulation method for solving the neutron transport equation, i.e., the method of characteristics (MOC), has garnered significant attention in whole-core high-fidelity simulations as the MOC can be applied to obtain detailed flux distributions without any restrictions on the geometry mesh. Hence, efficient performances can be achieved on various HPC systems. In fact, many famous high-fidelity neutron transport simulation codes were developed based on the MOC and have been used successfully in practical reactor applications, such as the CRX code [ 1 1 ] , DeCART code [ 2 2 , 3 3 ] , nTRACER [ 4 4 ] , MPACT [ 5 5 ] , NECP-X [ 6 6 ] , and STREAM [ 7 7 ] . However, the computational feasibility of whole-core MOC calculation depends on the acceleration technique used. The neutron transport equation contains seven independent variables associated with space, angle, energy, and time. The key to accelerating high-fidelity transport calculations is to establish a low-order equivalent equation relevant to the space, angle, energy, or time while minimizing the computational cost. However, previous studies focused more on space and energy accelerations. An equivalent low-order equation relevant to a coarse mesh or a coarse energy interval has been established to accelerate high-order transport calculations. The coarse mesh finite difference (CMFD) method is a popular method, and different types of CMFD methods [ 8 8 - 15 15 8-15 ] are widely utilized to accelerate MOC calculations. Furthermore, researchers have proposed a transient multilevel (TML) method to reduce the number of high-order transport calculations significantly [ 16 16 , 17 17 ] . Using the TML method, acceleration was achieved on the time scale; in this regard, large time steps were used for the MOC calculations, whereas small time steps were used for the low-order CMFD and point kinetics equations. In the high-fidelity core simulators established, the pin-homogenized multigroup CMFD (MG CMFD) method is widely used to accelerate the convergence of pin-resolved whole-core transport calculations. Excellent acceleration has been achieved. However, the MG CMFD method only updates the scalar flux and eigenvalue, and a bottleneck is encountered. Therefore, the MOC iterations will be reduced. Subsequently, the number of iterations cannot continue to decrease if only the MG CMFD is used for acceleration. In this case, the MOC calculations still incur a significant computational burden. Therefore, whole core pin-resolved MOC calculations require a more efficient acceleration technique. Fourier analysis [ 18 18 , 19 19 ] is a standard technique regularly used to investigate the stability of CMFD acceleration. Based on the Fourier analysis, the MG CMFD which based on flat source region or finer mesh can be used to further reduce the number of MOC iterations. Another potential acceleration technique is to directly update the incident angular flux and sources to further reduce the iterations of MOC calculations. Hence, an equivalent angular flux nonlinear finite difference (ANFD) equation must be established and applied to accelerate the MOC calculations. Furthermore, the domain decomposition technique is typically utilized for parallel computations, and the reactor core can be categorized into subdomains. Subsequently, the incident angular flux of each domain can be updated in time. However, the angular flux acceleration technique has not yet been investigated; hence, we will attempt to establish an equivalent ANFD equation in this study. The even parity discrete ordinates method can be used to derive a low-order angular flux nonlinear finite difference equation. However, ensuring the equivalence between this low-order equation and the high-order pin-resolved transport equation results in new challenges. The introduction of discontinuity factors used in popular CMFD methods is a good concept. Through the strict definition and calculation of even parity and odd parity discontinuity factors based on flux information obtained from MOC calculations, a completely equivalent angular flux nonlinear finite difference equation can be established. As relevant studies have not yet been performed, more details regarding the angular flux nonlinear finite difference equation based on the even parity discrete ordinates method, as well as the definition and calculation scheme of even and odd parity discontinuity factors including details for addressing different boundary conditions are discussed herein. Furthermore, an equivalent angular flux nonlinear finite difference equation with the pin-resolved transport equation was established in this study. This paper is organized as follows. In Sect. 2, we provide the necessary theory for establishing the equivalent low-order angular flux nonlinear finite difference equation. The numerical method to solve this equation is discussed in Sect. 3. A code based on this new method is developed, and two-dimensional (2D) pin, 2D assembly, and 2D C5G7 benchmark problems [ 20 20 ] are selected to verify the equivalence of the angular flux nonlinear finite difference method. The numerical results are discussed in Sect. 4, and important conclusions are provided. Sect. 5 presents the summary and conclusions. 2 Methodology of angular flux nonlinear finite difference equation 1 2.1 Even parity discrete ordinates method for neutron transport equation 2 To provide a detailed description of the angular flux nonlinear finite difference equation, we begin with the multigroup neutron transport equation in the Cartesian coordinate system, as follows: <math display="block" id="M1"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mo>⋅</mo><mo>∇</mo><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>S</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where vector r is the spatial position relative to the reference point, and vector Ω is the unit direction of flying with respect to the reference. Subscript g denotes the energy group index using multigroup approximation. This equation can represent the steady-state or transient neutron transport equation. When it represents a transient equation, the time variable is eliminated for brevity, as follows: <math display="block" id="M2"><mrow><msub><mi>S</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><mfrac><mn>1</mn><mi>v</mi></mfrac><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mn>1</mn><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><msub><mi>S</mi><mrow><mtext>f</mtext><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><mfrac><mn>1</mn><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><msub><mi>S</mi><mrow><mtext>d,</mtext><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><msub><mi>S</mi><mrow><mtext>e,</mtext><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where S s , S f , S d , and S e are the scattering, fission, delayed, and external sources, respectively. The external source is often disregarded, except when it must be considered. The fission and delayed neutron sources are isotropic. For the steady state, the time derivative term is zero, and the delayed neutron source is included in the fission source. If we discretize the time derivative term <math id="M3"><mrow><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow></mrow></math> with fluxes at the current and last time steps, the transient equation will become a transient fixed source equation, which can be solved in a manner similar to the steady-state equation. Therefore, we will only discuss the steady-state equation, and Eq. (2a) can be expressed as <math display="block" id="M4"><mrow><msub><mi>S</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mn>1</mn><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><msub><mi>S</mi><mrow><mtext>f</mtext><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>.</mo></mrow></math> The scattering and fission source can be expressed in terms of the angular flux as follows: <math display="block" id="M5"><mrow><msub><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mn>4</mn><mi>π</mi></mrow></msub><mrow><msub><mi>ψ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><msup><mi>Ω</mi><mo>′</mo></msup></mstyle></mrow><mo>)</mo></mrow><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><msub><mi>μ</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mtext>d</mtext><mstyle mathvariant="bold" mathsize="normal"><msup><mi>Ω</mi><mo>′</mo></msup></mstyle></mrow></mrow></mstyle></mrow></mstyle><mo>,</mo></mrow></math> <math display="block" id="M6"><mrow><msub><mi>S</mi><mrow><mtext>f</mtext><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msub><mi>χ</mi><mi>g</mi></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><mi>v</mi><msub><mi>Σ</mi><mtext>f</mtext></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mstyle><mo>,</mo></mrow></math> where ψ , φ , Σ , χ , μ 0 , and k eff are the angular flux, scalar flux, macro cross sections, fission spectrum, cosine of scattering angle in the center of mass system, and eigenvalue with a standard definition in nuclear reactor physics, respectively. Using the discrete ordinates method, we can represent the angular space with selected angles as follows: <math display="block" id="M7"><mrow><msub><mi>ϕ</mi><mi>g</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mn>4</mn><mi>π</mi></mrow></msub><mrow><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow><mo>)</mo></mrow><mtext>d</mtext><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow></mrow></mstyle><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>ω</mi><mi>i</mi></msub><msub><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mstyle><mo>,</mo></mrow></math> where <math id="M8"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub></mrow></math> and <math id="M9"><mrow><msub><mi>ω</mi><mi>i</mi></msub></mrow></math> are the selected angles and weights, respectively, and are known as the S N quadrature set; M is the total number of discrete directions. Here, we use subscript i to represent quantities at angle <math id="M10"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub></mrow></math> and subscript - i to represent quantities at angle <math id="M11"><mrow><mo>−</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub></mrow></math> . The equations for angular fluxes in opposite directions are as follows: <math display="block" id="M12"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msub><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msub><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> <math display="block" id="M13"><mrow><mo>−</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msub><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msub><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> The even and odd parity fluxes are expressed as follows: <math display="block" id="M14"><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><msub><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><msub><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>/</mo><mn>2</mn></mrow></math> <math display="block" id="M15"><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><msub><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>−</mo><msub><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>/</mo><mn>2</mn></mrow></math> The even and odd parity fluxes exhibit the following properties: <math display="block" id="M16"><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> <math display="block" id="M17"><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mo>−</mo><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> Using Eqs. (4) and (5), we can obtain the even and odd parity transport equations as follows: <math display="block" id="M18"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> <math display="block" id="M19"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M20"><mrow><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> and <math id="M21"><mrow><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> are the even and odd parity sources, respectively, expressed as <math display="block" id="M22"><mrow><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><msub><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><msub><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>/</mo><mn>2</mn></mrow></math> <math display="block" id="M23"><mrow><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><msub><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>−</mo><msub><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>/</mo><mn>2</mn></mrow></math> As the fission sources are isotropic, we have <math display="block" id="M24"><mrow><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msubsup><mi>S</mi><mrow><mtext>s,</mtext><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><mfrac><mn>1</mn><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><msub><mi>S</mi><mrow><mtext>f</mtext><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> <math display="block" id="M25"><mrow><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msubsup><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M26"><mrow><msubsup><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> and <math id="M27"><mrow><msubsup><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> are the even and odd parity scattering sources, expressed as <math display="block" id="M28"><mrow><msubsup><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><msub><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>+</mo><msub><mi>S</mi><mrow><mi mathvariant="normal">s</mi><mo>,</mo><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>/</mo><mn>2</mn><mo>,</mo></mrow></math> <math display="block" id="M29"><mrow><msubsup><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><msub><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>−</mo><msub><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>/</mo><mn>2.</mn></mrow></math> Here, <math id="M30"><mrow><msub><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> is the angle-dependent scattering source, expressed as <math display="block" id="M31"><mrow><msub><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mn>4</mn><mi>π</mi></mrow></msub><mrow><msub><mi>ψ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><msup><mi>Ω</mi><mo>′</mo></msup></mstyle></mrow><mo>)</mo></mrow><msub><mi>Σ</mi><mrow><mi>s</mi><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><msub><mi>μ</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mtext>d</mtext><mstyle mathvariant="bold" mathsize="normal"><msup><mi>Ω</mi><mo>′</mo></msup></mstyle></mrow></mrow></mstyle></mrow></mstyle><mo>.</mo></mrow></math> μ 0 is the cosine of the scattering angle in the center of mass system, where <math id="M32"><mrow><msub><mi>μ</mi><mn>0</mn></msub><mo>=</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><msup><mi>Ω</mi><mo>′</mo></msup></mstyle></mrow></math> . With a P n expansion of the scattering cross sections, we obtain <math display="block" id="M33"><mrow><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><msub><mi>μ</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn></mrow><mi>∞</mi></munderover><mrow><mfrac><mrow><mn>2</mn><mi>k</mi><mo>+</mo><mn>1</mn></mrow><mrow><mtext>4</mtext><mi>π</mi></mrow></mfrac><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mi>k</mi><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>P</mi><mi>k</mi></msub><mrow><mo>(</mo><mrow><msub><mi>μ</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></mstyle><mo>,</mo></mrow></math> where P k are the Legendre polynomials and <math id="M34"><mrow><msub><mi>Σ</mi><mrow><mi>s</mi><mo>,</mo><mi>k</mi></mrow></msub></mrow></math> are P k scattering moments. By expanding the even and odd parity scattering sources with a P 0 expansion, we have <math display="block" id="M35"><mrow><msubsup><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mtext>4</mtext><mi>π</mi></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>i</mi><mo>′</mo></msup></munder><mrow><msub><mi>ω</mi><msup><mi>i</mi><mo>′</mo></msup></msub><msubsup><mi>ψ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><msup><mi>i</mi><mo>′</mo></msup></mrow><mrow/></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mstyle></mrow></mstyle><mo>,</mo></mrow></math> <math display="block" id="M36"><mrow><msubsup><mi>S</mi><mrow><mtext>s</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mn>0.</mn></mrow></math> Combined with Eqs. (9) and (13), the even and odd parity sources can be expressed as follows: <math display="block" id="M37"><mrow><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><msub><mi>Σ</mi><mrow><mtext>s,</mtext><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>i</mi><mo>′</mo></msup></munder><mrow><msub><mi>ω</mi><msup><mi>i</mi><mo>′</mo></msup></msub><msubsup><mi>ψ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><msup><mi>i</mi><mo>′</mo></msup></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mstyle></mrow></mstyle><mo>+</mo><mfrac><mn>1</mn><mrow><mtext>4</mtext><mi>π</mi></mrow></mfrac><mfrac><mrow><msub><mi>χ</mi><mi>g</mi></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><mi>v</mi><msub><mi>Σ</mi><mrow><mtext>f</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>i</mi><mo>′</mo></msup></munder><mrow><msub><mi>ω</mi><msup><mi>i</mi><mo>′</mo></msup></msub><msubsup><mi>ψ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><msup><mi>i</mi><mo>′</mo></msup></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mstyle></mrow></mstyle><mo>,</mo></mrow></math> <math display="block" id="M38"><mrow><msubsup><mi>S</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mn>0.</mn></mrow></math> Integrate the even parity transport equation in Eq. (7a) over a region with constant a cross section as follows: <math display="block" id="M39"><mrow><mstyle displaystyle="true"><mrow><msub><mo>∬</mo><mi>Γ</mi></msub><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle><mtext>d</mtext><mi>Γ</mi></mrow></mrow></mstyle><mo>+</mo><msub><mi>V</mi><mi>r</mi></msub><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>=</mo><msub><mi>V</mi><mi>r</mi></msub><msubsup><mover accent="true"><mover accent="true"><mi>S</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>,</mo></mrow></math> where subscript r is the region index; <math id="M40"><mrow><msub><mi>V</mi><mi>r</mi></msub></mrow></math> is the volume of region r ; <math id="M41"><mi>Γ</mi></math> is the surface of region r ; <math id="M42"><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle></math> is the unit normal direction vector of the surface; <math id="M43"><mrow><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></math> and <math id="M44"><mrow><msubsup><mover accent="true"><mover accent="true"><mi>S</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></math> are the node average even parity flux and node average even parity source weighted by the region volume, respectively. <math id="M45"><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow></math> is the angle-dependent total macro cross section, defined as follows: <math display="block" id="M46"><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><mrow><msub><mo>∭</mo><mi>V</mi></msub><mrow><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><mrow><msub><mo>∭</mo><mi>V</mi></msub><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle></mrow></mfrac><mo>.</mo></mrow></math> It is noteworthy that the divergence theorem was used for deriving Eq. (15). <math display="block" id="M47"><mrow><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>⋅</mo><mtext>+</mtext><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mo>∇</mo><mo>⋅</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>=</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M48"><mrow><mstyle displaystyle="true"><mrow><msub><mo>∬</mo><mi>Γ</mi></msub><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle><mtext>d</mtext><mi>Γ</mi></mrow></mrow></mstyle><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∭</mo><mi>V</mi></msub><mrow><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub></mrow><mo>)</mo></mrow><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∭</mo><mi>V</mi></msub><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle><mo>.</mo></mrow></math> The surface of region r is divided into n r segments, where n r is the number of neighbors. Each segment is the interface between region r and one of its neighbors. Define the average odd parity flux on each segment as follows: <math display="block" id="M49"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><mfrac><mtext>1</mtext><mrow><msub><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><msub><mo>∬</mo><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></msub><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>−</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle><mtext>d</mtext><mi>Γ</mi></mrow></mrow></mstyle><mo>,</mo></mrow></math> <math display="block" id="M50"><mrow><msub><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∬</mo><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></msub><mrow><mrow><mo>|</mo><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle></mrow><mo>|</mo></mrow><mtext>d</mtext><mi>Γ</mi></mrow></mrow></mstyle><mo>,</mo></mrow></math> where n is the neighbor node of node r , and <math id="M51"><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></math> is the interface between r and n . The average odd parity flux on the surface can represent the net number of neutrons passing through the interface in a certain direction, and it exhibits the following property: <math display="block" id="M52"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><mo>−</mo><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></math> By rearranging Eq. (15), we can obtain the even parity flux equation as follows: <math display="block" id="M53"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>∈</mo><mtext>neighbors</mtext><mtext> </mtext></mrow><mrow/></munderover><mrow><msub><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mstyle><mo>+</mo><msub><mi>V</mi><mi>r</mi></msub><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>=</mo><msub><mi>V</mi><mi>r</mi></msub><msubsup><mover accent="true"><mover accent="true"><mi>S</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>.</mo></mrow></math> The even parity flux equation exactly preserves the neutron balance in each direction and in each coarse mesh. If the relationship between the average odd parity flux <math id="M54"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></math> and the node average even parity flux <math id="M55"><mrow><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></math> is obtained, then Eq. (20) can be solved further. This relationship can be obtained from the odd parity transport equation. Integrating the odd parity transport equation (Eq. (7b)) over each segment of the surface, we obtain <math display="block" id="M56"><mrow><mstyle displaystyle="true"><mrow><msub><mo>∬</mo><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></msub><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle><mtext>d</mtext><mi>Γ</mi></mrow></mrow></mstyle><mo>+</mo><msub><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><mn>0.</mn></mrow></math> The left hand of Eq. (21) is the average derivative of the even parity flux on interface <math id="M57"><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></math> in the <math id="M58"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub></mrow></math> direction and can be defined as follows: <math display="block" id="M59"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msup><mrow/><mo>′</mo></msup><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><msub><mo>∬</mo><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></msub><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle><mtext>d</mtext><mi>Γ</mi></mrow></mrow></mstyle><mo>.</mo></mrow></math> Substituting Eq. (22) into Eq. (21), we have <math display="block" id="M60"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msup><mrow/><mo>′</mo></msup><mo>=</mo><mo>−</mo><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>.</mo></mrow></math> Through a finite difference approximation for <math id="M61"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msup><mrow/><mo>′</mo></msup></mrow></math> , we obtain the following odd parity flux equation: <math display="block" id="M62"><mrow><mfrac><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow><mrow><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub></mrow></mfrac><mo>≈</mo><mo>−</mo><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>,</mo></mrow></math> where <math id="M63"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow></math> is the average even parity flux on the surface between r and n , and <math id="M64"><mrow><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub></mrow></math> is the thickness of the node center to the surface center in the direction of angle <math id="M65"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub></mrow></math> . <math display="block" id="M66"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><msub><mo>∬</mo><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></msub><mrow><msubsup><mi>ψ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>)</mo></mrow><mrow><mo>|</mo><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle></mrow><mo>|</mo></mrow><mtext>d</mtext><mi>Γ</mi></mrow></mrow></mstyle><mo>,</mo></mrow></math> <math display="block" id="M67"><mrow><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>=</mo><mrow><mo>|</mo><mrow><mfrac><mrow><mstyle displaystyle="true"><mrow><msub><mo>∬</mo><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></msub><mrow><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>⋅</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub></mrow><mo>)</mo></mrow><mrow><mo>|</mo><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle></mrow><mo>|</mo></mrow><mtext>d</mtext><mi>Γ</mi></mrow></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><mrow><msub><mo>∬</mo><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></msub><mrow><mrow><mo>|</mo><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle></mrow><mo>|</mo></mrow><mtext>d</mtext><mi>Γ</mi></mrow></mrow></mstyle></mrow></mfrac><mo>−</mo><mfrac><mrow><mstyle displaystyle="true"><mrow><msub><mo>∭</mo><mi>V</mi></msub><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>⋅</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><mrow><msub><mo>∭</mo><mi>V</mi></msub><mrow><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle></mrow></mfrac></mrow><mo>|</mo></mrow><mo>.</mo></mrow></math> The average even parity flux on the surface exhibits following priority: <math display="block" id="M68"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>=</mo><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mo>−</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>.</mo></mrow></math> 2.2 Definition and application of discontinuity factor 2 To render all <math id="M69"><mrow><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></math> , <math id="M70"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow></math> and <math id="M71"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></math> in the odd parity flux equation (Eq. (24)) satisfy the MOC transport calculation, we introduced a factor known as the even parity discontinuity factor, <math id="M72"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow></math> , as follows: <math display="block" id="M73"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>=</mo><mfrac><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow><mo>)</mo></mrow></mrow></mfrac></mrow></math> <math id="M74"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow></math> will be greater than zero when the denominator of Eq. (27) is positive. However, the denominator can less than or equal to zero. In this case, <math id="M75"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow></math> will be infinite or negative. To avoid this, we introduced another discontinuity factor to Eq. (27), i.e., the odd parity discontinuity factor, <math id="M76"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></math> . <math display="block" id="M77"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>=</mo><mfrac><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow><mo>)</mo></mrow></mrow></mfrac></mrow></math> An algorithm was designed to determine the even and odd parity discontinuity factors. This algorithm begins by evaluating the following conditional: <math display="block" id="M78"><mrow><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>></mo><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>.</mo></mrow></math> If the condition in Eq. (29) is true, then the even and odd parity discontinuity factors can be calculated from the MOC solutions as follows: <math display="block" id="M79"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>=</mo><mfrac><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mrow><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M80"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><mn>1.</mn></mrow></math> If the condition is false, then different equations are designed to calculate the even and odd parity discontinuity factors as follows: <math display="block" id="M81"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>=</mo><mfrac><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>+</mo><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow><mrow><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M82"><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>.</mo></mrow></math> Using this approach, the even and odd parity discontinuity factors will be positive at all times. Rearranging Eq. (28) yields <math display="block" id="M83"><mrow><mfrac><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow></mfrac><mo>+</mo><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>.</mo></mrow></math> If neighbor node n is within the computational domain, then a similar equation is obtained for the other side of the interface <math id="M84"><mrow><msub><mi>Γ</mi><mrow><mi>r</mi><mi>n</mi></mrow></msub></mrow></math> , as follows: <math display="block" id="M85"><mrow><mfrac><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></mfrac><mo>+</mo><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t,</mtext><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow><mo>+</mo></msubsup><mo>.</mo></mrow></math> Combining Eqs. (32) and (33) and the continuity condition <math id="M86"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>=</mo><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></math> , <math id="M87"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><mo>−</mo><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>−</mo></msubsup></mrow></math> , we obtain <math display="block" id="M88"><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>−</mo></msubsup></mrow></mfrac></mrow><mo>)</mo></mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow><mo>+</mo></msubsup><mo>.</mo></mrow></math> As for the vacuum boundary surface, the incident angular flux is zero. Therefore, we have the equation <math id="M89"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow></math> , which can be derived from Eqs. (18a) and (25a). Combining this equation with Eq. (32), we obtain <math display="block" id="M90"><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>.</mo></mrow></math> For the reflective boundary surface, the angular fluxes in the reflect directions are equal. Therefore, we have the continuity condition <math id="M91"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>=</mo><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow></math> and <math id="M92"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><mo>−</mo><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></math> , where subscript i* and i represent a pair of angles in the incident and reflective directions, respectively. Fig. 1 Fig. 1 shows an illustration of the reflective boundary condition. Combined with Eq. (32) for both the angle i and i* , we obtain <math display="block" id="M93"><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac></mrow><mo>)</mo></mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></math> In general, the relationships between the surface average odd parity flux and node average even parity flux for both the interior surface and boundary surface are established through Eqs. (34), (35), and (36). Moreover, these equations satisfy the high-order MOC solution through the strict definitions of the even and odd parity discontinuity factors. 2.3 Construction of angular flux nonlinear finite difference linear system 2 Based on Eqs. (34), (35), and (36), the odd parity fluxes on the interior and boundary surfaces can be expressed as a linear function of the neighbor average even parity fluxes using the coupling coefficient C . For the interior surface, we have <math display="block" id="M94"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow><mo>+</mo></msubsup><mo>,</mo></mrow></math> where <math display="block" id="M95"><mrow><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>=</mo><mfrac><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>−</mo></msubsup></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M96"><mrow><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub><mo>=</mo><mfrac><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>−</mo></msubsup></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mfrac><mo>.</mo></mrow></math> For the vacuum boundary surface, we have <math display="block" id="M97"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>,</mo></mrow></math> where <math display="block" id="M98"><mrow><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>=</mo><mfrac><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></mfrac><mo>.</mo></mrow></math> For the reflective boundary surface, we have <math display="block" id="M99"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>,</mo></mrow></math> where <math display="block" id="M100"><mrow><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>=</mo><mfrac><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M101"><mrow><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>=</mo><mfrac><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mfrac><mo>.</mo></mrow></math> Substituting Eqs. (37)–(39) into Eq. (20), we obtain the following angular flux nonlinear difference linear system: <math display="block" id="M102"><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>n</mi><mo>∈</mo><mtext>neighbors</mtext><mtext> </mtext></mrow></munder><mrow><msub><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mo>)</mo></mrow></mrow></mstyle><mo>+</mo><msub><mi>V</mi><mi>r</mi></msub><msub><mi>Σ</mi><mrow><mi>t</mi><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>=</mo><msubsup><mi>b</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mrow/></msubsup><mo>.</mo></mrow></math> If the homogenized cross sections and coupling coefficients are obtained from the results of the MOC transport calculation, then the ANFD linear system can be solved. As the even parity fluxes in opposite directions are equal, only the ANFD linear system need to be solved on a half angular space. 3 Solving angular flux nonlinear difference equation 1 To solve the angular flux nonlinear difference equation, different solving strategies for decoupling the energy group and angle are proposed herein. If the scattering sources are determined simultaneously with the even parity fluxes, then all groups and all angles are coupled together. The terms on the right side of the linear system can be expressed as follows: <math display="block" id="M103"><mrow><msubsup><mi>b</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mrow/></msubsup><mo>=</mo><msub><mi>V</mi><mi>r</mi></msub><mfrac><mn>1</mn><mrow><msub><mi>k</mi><mrow><mi>e</mi><mi>f</mi><mi>f</mi></mrow></msub></mrow></mfrac><mfrac><mrow><msub><mi>χ</mi><mi>g</mi></msub></mrow><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><mi>v</mi><msub><mover accent="true"><mi>Σ</mi><mo>¯</mo></mover><mrow><mtext>f</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><mi>r</mi></mrow></msub><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>i</mi><mo>′</mo></msup></munder><mrow><msub><mi>ω</mi><msup><mi>i</mi><mo>′</mo></msup></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><msup><mi>i</mi><mo>′</mo></msup><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></mstyle></mrow></mstyle><mo>.</mo></mrow></math> If all scattering sources other than the in-group scattering term are provided, then all angles are coupled together, and the linear system can still be solved by group separately, as follows: <math display="block" id="M104"><mrow><msubsup><mi>b</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mrow/></msubsup><mo>=</mo><msub><mi>V</mi><mi>r</mi></msub><mfrac><mn>1</mn><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mfrac><mrow><msub><mi>χ</mi><mi>g</mi></msub></mrow><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><mi>v</mi><msub><mover accent="true"><mi>Σ</mi><mo>¯</mo></mover><mrow><mtext>f</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><mi>r</mi></mrow></msub><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>i</mi><mo>′</mo></msup></munder><mrow><msub><mi>ω</mi><msup><mi>i</mi><mo>′</mo></msup></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><msup><mi>i</mi><mo>′</mo></msup><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></mstyle></mrow></mstyle><mo>+</mo><msub><mi>V</mi><mi>r</mi></msub><mfrac><mn>1</mn><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>≠</mo><mi>g</mi></mrow></munder><mrow><msub><mover accent="true"><mi>Σ</mi><mo>¯</mo></mover><mrow><mtext>s</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi><mo>,</mo><mi>r</mi></mrow></msub><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>i</mi><mo>′</mo></msup></munder><mrow><msub><mi>ω</mi><msup><mi>i</mi><mo>′</mo></msup></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><msup><mi>i</mi><mo>′</mo></msup><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></mstyle></mrow></mstyle><mo>.</mo></mrow></math> If the scattering sources are all treated as specified, then the linear system can be solved by group and angle separately, as follows: <math display="block" id="M105"><mrow><msubsup><mi>b</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mrow/></msubsup><mo>=</mo><msub><mi>V</mi><mi>r</mi></msub><mfrac><mn>1</mn><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mfrac><mrow><msub><mi>χ</mi><mi>g</mi></msub></mrow><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><mi>v</mi><msub><mover accent="true"><mi>Σ</mi><mo>¯</mo></mover><mrow><mtext>f</mtext><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><mi>r</mi></mrow></msub><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>i</mi><mo>′</mo></msup></munder><mrow><msub><mi>ω</mi><msup><mi>i</mi><mo>′</mo></msup></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><msup><mi>i</mi><mo>′</mo></msup><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></mstyle></mrow></mstyle><mo>+</mo><msub><mi>V</mi><mi>r</mi></msub><mfrac><mn>1</mn><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><msub><mover accent="true"><mi>Σ</mi><mo>¯</mo></mover><mrow><mi>s</mi><mo>,</mo><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi><mo>,</mo><mi>r</mi></mrow></msub><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>i</mi><mo>′</mo></msup></munder><mrow><msub><mi>ω</mi><msup><mi>i</mi><mo>′</mo></msup></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>,</mo><msup><mi>i</mi><mo>′</mo></msup><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow></mstyle></mrow></mstyle><mo>+</mo><msub><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>*</mo><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>.</mo></mrow></math> It is noteworthy that the last term in Eq. (41c) appears only when a reflect boundary condition exists on the node surface. 3.1 Decoupling of energy group and angle 2 To solve the linear system by group and angle separately, we would like to keep <math id="M106"><mrow><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>=</mo><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub></mrow></math> in the interior surfaces such that the linear system will be symmetric. A symmetric linear system is much easier to solve than an asymmetric system. To create a symmetric ANFD linear system, the coupling coefficients for the interior surface are redefined as follows: <math display="block" id="M107"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>=</mo><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mrow><mo>(</mo><mrow><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mo>)</mo></mrow><mo>−</mo><msub><mrow><menclose notation="downdiagonalstrike"><mi>ψ</mi></menclose></mrow><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>,</mo></mrow></math> <math display="block" id="M108"><mrow><msubsup><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mrow/></msubsup><mo>=</mo><msubsup><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mrow/></msubsup><mo>=</mo><mfrac><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><mo>+</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t,</mtext><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>−</mo></msubsup></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M109"><mrow><msub><mrow><menclose notation="downdiagonalstrike"><mi>ψ</mi></menclose></mrow><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mtext>=</mtext><mfrac><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup><mo>−</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow><mo>+</mo></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup></mrow></mfrac><mo>+</mo><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>+</mo></msubsup><msub><mi>h</mi><mrow><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub><mfrac><mrow><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi></mrow></msub></mrow><mrow><msubsup><mi>f</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow><mo>−</mo></msubsup></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M110"><mrow><msub><mrow><menclose notation="downdiagonalstrike"><mi>ψ</mi></menclose></mrow><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub></mrow></math> is the asymmetric term of Eq. (42a), which can be determined using even parity fluxes. Substituting Eq. (42) into Eq. (40) and placing the asymmetric term to the right side of the equation, the ANFD linear system is rearranged as follows: (43) where N , S , W , E , T , and B represent the north, south, west, east, top, and bottom surfaces of the node , respectively. The superscript k denotes the power iteration number. As the symmetric and diagonal dominant matrix is a symmetric positive definite matrix, Eq. (43) will become an expression for a symmetric positive definite linear system. Furthermore, the preconditioned conjugate gradient (CG) method [ 21 21 ] can be used to solve this linear system. It is noteworthy that decoupling the angle and energy group will increase the iteration number of the ANFD method. However, the ANFD computing costs are low compared with those of MOC calculations, and they can be reduced significantly using the multilevel acceleration method, which has been implemented in several established core simulators [ 22 22 - 24 24 22-24 ] . The main idea of the multilevel acceleration method is to use a few-group CMFD linear system, such as a one- or two-group CMFD, to reduce the computing burden of the MG CMFD linear system. The few-group CMFD linear system is used to provide an effective estimation of the eigenvalue at the beginning of each power iteration, thereby significantly reducing the outer iteration number of the MG CMFD linear system. Therefore, a new approach similar to the multilevel acceleration method can be used to reduce the iteration number of ANFD required for convergence. This will be investigated further in future studies. Meanwhile, if we wish to solve the linear system with all angles coupled or all energy groups coupled, it would be meaningless to maintain <math id="M111"><mrow><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>=</mo><msub><mi>C</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>n</mi><mi>r</mi></mrow></msub></mrow></math> . In this case, we will define the interior surface coupling coefficients directly from Eq. (37), and the preconditioned CG method will not be used. 3.2 Preconditioned Krylov subspace method 2 The direct solution of this linear system is unpractical as the matrix can become highly ill-conditioned, causing problems in the numerical solution. The Krylov subspace method has been proven to be extremely effective for solving large, nonsymmetric, sparse linear systems, particularly when a good preconditioner is applied. In this study, the right preconditioned generalized minimal residual (GMRES) method was used to solve the ANFD linear system. The right preconditioned GMRES algorithm can be summarized as follows: Meanwhile, if we decoupling the energy group and angle to construct a symmetric positive definite linear system, then the preconditioned conjugate gradient method can be used to solve the ANFD linear system. The left preconditioned conjugate gradient algorithm can be summarized as follows: As for precondition technology, a self-developed RSILU preconditioner has been used to improve parallel computing capability and minimize computational burden. The RSILU preconditioner does not require a multicolor ordering strategy for parallel computing and can provide convergence rates comparable to the standard incomplete lower–upper preconditioner. More details are available in the literature [ 25 25 ] . 3.3 Iteration scheme of ANFD method 2 Once the node and surface angular flux obtained from the MOC calculation become available, the ANFD homogenized cross sections as well as the even and odd parity discontinuity factors can be easily obtained using Eq. (16) and (30), respectively. The strategy of decoupling the energy group and angle was used to construct the ANFD linear system, and the coupling coefficients were calculated using Eqs. (37)–(39). Hence, an equivalent ANFD linear system can be established. It is noteworthy that the asymmetric term arising from the construction of a symmetric positive definite matrix should be updated with the even party angular flux when solving the ANFD linear system. The right preconditioned GMRES solver was used to solve the ANFD linear system in this study. The exit from the GMRES solver to the ANFD cycle occurs when the required reduction in the relative L2-norm of the residual vector is achieved. After solving the ANFD linear system, new even parity flux distributions and eigenvalues were obtained. To update the 2D MOC parameters, the average scalar flux on each node and the average angular flux on each surface were reconstructed as follows, indicated as <math id="M112"><mrow><msub><mover accent="true"><mi>ϕ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow></math> and <math id="M113"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub></mrow></math> , respectively: <math display="block" id="M114"><mrow><msub><mover accent="true"><mi>ϕ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>r</mi></mrow></msub><mo>=</mo><mn>2</mn><mi>π</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><mrow><mo>(</mo><mrow><msub><mi>w</mi><mi>i</mi></msub><msubsup><mover accent="true"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi></mrow><mo>+</mo></msubsup></mrow><mo>)</mo></mrow></mrow></mstyle></mrow></math> <math display="block" id="M115"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub><mo>=</mo><mfrac><mrow><mrow><mo>|</mo><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle></mrow><mo>|</mo></mrow></mrow><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mi>i</mi></msub><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>n</mi><mo>^</mo></mover></mstyle></mrow></mfrac><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>−</mo></msubsup><mo>+</mo><msubsup><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow><mo>+</mo></msubsup></mrow></math> The former was used to update the scalar flux in the flat source regions, whereas the latter was used to update the incident angular flux on the boundary surface. The M in Eq. (44a) denotes the total number of discrete directions, which is consistent with the number of MOC azimuthal angles. After updating the average node scalar flux <math id="M116"><mrow><msub><mover accent="true"><mi>ϕ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>r</mi></mrow></msub></mrow></math> and average surface angular flux <math id="M117"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>g</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>r</mi><mi>n</mi></mrow></msub></mrow></math> , 2D MOC calculations were performed, and the cell-averaged cross sections as well as the even and odd parity discontinuity factors were updated. Subsequently, ANFD was performed using the updated cross sections as well as the even and odd parity discontinuity factors. This process was repeated until all convergence criteria were satisfied. The iteration scheme of the ANFD method is depicted in Fig. 2 Fig. 2 . 4 Numerical results 1 To investigate the equivalence of the ANFD method proposed herein, a code written in the C language was developed and three different test problems were selected: 2D pin, 2D assembly, and 2D C5G7 benchmark problems. In addition, a self-defined problem was tested to compare the acceleration capabilities of ANFD and CMFD in the void region. For all these problems, the flux convergence criterion was 10 -5 by comparing the infinite norm of the scalar flux in the flat source region, and the eigenvalue convergence criterion was 10 -7 . For the ray tracing parameter, a 0.01 cm ray spacing with 72 azimuthal angles and a Tabuchi–Yamamoto polar quadrature [ 26 26 ] with three polar angles per half-space were used for all the problems. The computer used to execute the code comprised an Inter(R) Core (TM) i5-7200U CPU @ 2.50 GHz with an 8.0G RAM. The basic computing unit for the ANFD was the pin cell. For the 2D pin cell and 2D assembly problems, the eigenvalue and flux distribution error in the angle, energy group, and space were compared with those of the MOC calculations. In addition, three parameters were selected to assess the overall relative error of the even parity flux, i.e., the average relative error (AVG), root mean square of relative error (RMS), and maximum relative error (MAX), defined as follows: <math display="block" id="M118"><mrow><mtext>AVG</mtext><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>N</mi></munder><mrow><mrow><mo>|</mo><mrow><msub><mi>e</mi><mi>n</mi></msub></mrow><mo>|</mo></mrow></mrow></mstyle></mrow><mi>N</mi></mfrac></mrow></math> <math display="block" id="M119"><mrow><mtext>MAX</mtext><mo>=</mo><mi>max</mi><mrow><mo>{</mo><mrow><msub><mi>e</mi><mi>n</mi></msub></mrow><mo>}</mo></mrow></mrow></math> <math display="block" id="M120"><mrow><mtext>RMS</mtext><mo>=</mo><mfrac><mrow><msqrt><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>N</mi></munder><mrow><msub><mi>e</mi><mi>n</mi></msub><msup><mrow/><mn>2</mn></msup></mrow></mstyle></mrow></msqrt></mrow><mi>N</mi></mfrac><mo>,</mo></mrow></math> where e n is the relative error of the even parity flux compared with that of the MOC solution, and N is the total number of even parity fluxes. 4.1 2D pin cell problem 2 As indicated in Fig. 3(a) Fig. 3(a) , the side length of the cell was 1.26 cm, and the radius of the fuel pins was 0.54 cm for the 2D pin cell problem. The cross sections data for the UO2 fuel and moderator were obtained from the 2D C5G7 benchmark report. A vacuum boundary was applied to the east and south boundary surfaces, and the reflected boundary was applied to the north and west boundary surfaces. Three moderator rings and five fuel rings with 16 azimuthal divisions were selected for 128 flat-source regions. For the 2D pin cell problem, Fig. 5(a) Fig. 5(a) shows the angular distribution of the even parity flux in the first energy group, and Fig. 5(b) Fig. 5(b) shows the energy group distribution of the even parity flux in the first azimuth angle. It is clear that the even parity flux distribution of the ANFD was consistent with the MOC solution in both the angle and energy groups. Table 1 Table 1 shows the eigenvalue difference and overall even parity flux relative error. The difference in the eigenvalue was 0 pcm, and the flux maximum relative error was 3.683E-06. 3 1 k eff 3 1 Flux relative error 1 1 MOC 1 1 ANFD 1 1 Error (pcm) 1 1 AVG 1 1 RMS 1 1 MAX 1 1 0.01561 1 1 0.01561 1 1 0 1 1 3.495E-07 1 1 9.145E-07 1 1 3.683E-06 4.2 2D assembly problem 2 As indicated in Fig. 3(b) Fig. 3(b) , the 2D assembly problem comprised a 17 × 17 lattice of square pin cells containing a UO2 fuel pin cell, guide tube, and fission chamber. The cross sections data for these materials were obtained from the 2D C5G7 benchmark problem. The division of the flat-source region was the same as that of the 2D pin cell problem. As for the boundary conditions, vacuum boundary conditions were applied to the east and south boundary surfaces, whereas the reflected boundary conditions were applied to the north and west boundary surfaces. Fig. 6 Fig. 6 shows the space relative error of the even parity flux in the first energy group and first azimuth angle for the 2D assembly problem. As shown, the even parity flux distribution of the ANFD indicated a relatively good agreement in space compared with the MOC solution. Table 2 Table 2 shows the eigenvalue difference and overall even parity flux relative error for the 2D assembly problem. The difference in the eigenvalue was 0 pcm, whereas the flux maximum relative error was 1.162E-05. 3 1 k eff 3 1 Flux relative error 1 1 MOC 1 1 ANFD 1 1 Error (pcm) 1 1 AVG 1 1 RMS 1 1 MAX 1 1 0.96848 1 1 0.96848 1 1 0 1 1 1.030E-06 1 1 1.742E-06 1 1 1.162E-05 4.3 2D C5G7 problem 2 The 2D C5G7 benchmark problem was originally a pin-cell level heterogeneous benchmark for testing the accuracy of modern deterministic transport methods and codes without spatial homogenization. The 2D C5G7 problem configuration is shown in Fig. 4 Fig. 4 . The overall geometry of this problem was 64.26 cm × 64.26 cm, and each assembly measured 21.42 cm × 21.42 cm. Each fuel assembly comprised a 17 × 17 lattice of square pin cells, including a moderator, guide tube, fission chamber, MOX fuel pin cell with different enrichment, and UO2 fuel pin cell. For the reflector pin cases, the six additional pin cells of each reflector assembly used 1 × 1 sub mesh, whereas the 11 pin cells near the fuel used a 6 × 6 Cartesian sub mesh. For the other pin cells of the fuel and control rods, 128 flat-source regions comprising five fuel rings and three moderator rings with 16 azimuthal divisions were used. An accurate Monte-Carlo reference solution was provided in the 2D C5G7 benchmark report, which contains a precise eigenvalue solution and normalized pin power (NPP) that can be used as a reference solution. The NPP ( <math id="M121"><mrow><msub><mi>p</mi><mi>r</mi></msub></mrow></math> ) in each fuel pin is defined as follows: <math display="block" id="M122"><mrow><msub><mi>p</mi><mi>r</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>R</mi><mrow><mtext>f,</mtext><mi>r</mi></mrow></msub></mrow><mrow><msub><mover accent="true"><mi>R</mi><mo>¯</mo></mover><mtext>f</mtext></msub></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M123"><mrow><msub><mi>R</mi><mrow><mtext>f,</mtext><mi>r</mi></mrow></msub><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><msub><mi>Σ</mi><mrow><mtext>f</mtext><mo>,</mo><mi>g</mi><mo>,</mo><mi>r</mi></mrow></msub><msub><mi>ϕ</mi><mrow><mi>g</mi><mo>,</mo><mi>r</mi></mrow></msub><msub><mi>V</mi><mi>r</mi></msub></mrow></mstyle></mrow></math> <math display="block" id="M124"><mrow><msub><mover accent="true"><mi>R</mi><mo>¯</mo></mover><mtext>f</mtext></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>r</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msub><mi>R</mi><mrow><mtext>f,</mtext><mi>r</mi></mrow></msub></mrow></mstyle></mrow><mi>N</mi></mfrac><mo>,</mo></mrow></math> where <math id="M125"><mrow><msub><mi>R</mi><mrow><mtext>f</mtext><mo>,</mo><mi>r</mi></mrow></msub></mrow></math> is the fission rate in each fuel pin, <math id="M126"><mrow><msub><mover accent="true"><mi>R</mi><mo>¯</mo></mover><mtext>f</mtext></msub></mrow></math> the average fission rate of all fuel pins, and <math id="M127"><mi>N</mi></math> the number of all fuel pins. Figure 7 Figure 7 indicates the percent error of the NPP for the 2D C5G7 benchmark problem; as shown, the main error appeared in the fuel pin cell near the reflector. Table 3 Table 3 shows the eigenvalue difference and NPP error from our results and the MCNP reference solutions. Compared with the reference results, the error of the eigenvalue was 3 pcm, which was less than the stochastic uncertainty of the MCNP reference results. In general, the eigenvalue, NPP, etc. agreed well with the reference results. 1 1 Parameters 1 1 Reference MCNP (Statistic error) 1 1 ANFD 1 1 Error 1 1 Eigenvalue 1 1 k eff (error/pcm) 1 1 1.18655 (±9.5 pcm) 1 1 1.18652 1 1 -3 pcm 1 1 Specific NPP data 1 1 Maximum NPP (error) 1 1 2.498 (±0.16%) 1 1 2.494 1 1 -0.14% 1 1 Minimum NPP (error) 1 1 0.232 (±0.58%) 1 1 0.234 1 1 +1.04% 1 1 Normalized assembly power data 1 1 Inner UO2 power (error) 1 1 492.8 (±0.10%) 1 1 492.3 1 1 -0.10% 1 1 MOX power (error) 1 1 211.7 (±0.18%) 1 1 211.9 1 1 +0.10% 1 1 Outer UO2 power (error) 1 1 139.8 (±0.20%) 1 1 139.9 1 1 +0.10% 1 1 NPP error results 1 1 MAX error 1 1 -- 1 1 1.25% 1 1 -- 1 1 AVG error 1 1 0.32% 1 1 0.21% 1 1 -- 1 1 RMS error 1 1 0.34% 1 1 0.30% 1 1 -- 1 1 MRE error 1 1 0.27% 1 1 0.16% 1 1 -- 4.4 Void region problem 2 The ANFD method was developed from transport theory and does not require the calculation of the diffusion coefficient. Therefore, the ANFD method may exhibit a better acceleration than the CMFD method in certain strongly anisotropic cases, such as those involving void regions and anisotropic scattering. In this study, a lattice problem involving void regions was constructed and tested; 5 × 5 square pin cells were set to the void regions, whereas UO 2 fuel cells were used in the other regions. The geometry and cross sections of the UO 2 fuel cell were the same as those provided in the 2D C5G7 benchmark report. The reflected boundary conditions were applied to all boundary surfaces. The layout of this void region problem is shown in Fig. 8 Fig. 8 . The convergence history of the MOC, ANFD, and CMFD are shown in Table 4 Table 4 and Fig. 9 Fig. 9 . As expected, the CMFD failed to converge in this problem, whereas the ANFD converged and remained effective. 1 1 1 1 k eff 1 1 Error 1 1 MOC iteration number 1 1 Speedup 1 1 No acceleration 1 1 1.32592 1 1 -- 1 1 91 1 1 -- 1 1 CMFD 1 1 Divergence 1 1 -- 1 1 -- 1 1 -- 1 1 ANFD 1 1 1.32602 1 1 10 pcm 1 1 11 1 1 8.3 Void regions are typical in real reactors, such as TREAT [ 27 27 ] and HTR-10 [ 28 28 ] . In a previous study [ 29 29 ], the quasi-diffusion method was used to address the void region problem in TREAT core analysis. In this study, the quasi-diffusion method was implemented using the nodal expansion method. Meanwhile, the nodal wise average cross sections, Eddington factors, and discontinuity factors were generated from 3D whole-core Serpent Monte-Carlo simulations. By comparison, the ANFD is more concise and only requires even and odd parity discontinuity factors generated from MOC calculations. 5 Conclusion 1 In this study, an innovative equivalent low-order ANFD method for MOC transport calculations was developed and implemented. The even parity discrete ordinates method was used to build a low-order ANFD equation, and different boundary condition treatments were proposed. Even and odd parity discontinuity factors based on the flux information obtained from MOC calculations were strictly defined to ensure the equivalence between this low-order equation and the high-order pin-resolved transport equation. The energy group and angle were decoupled to construct an easier-to-solve symmetric linear system. The self-developed RSILU preconditioned GMRES solver was applied to solve this linear system. The equivalence of the ANFD method was verified for 2D pin, 2D assembly, and 2D C5G7 benchmark problems. The numerical results demonstrated that a low-order ANFD equation that was completely equivalent to the pin-resolved transport equation was established. This study provides an important foundation for further developing a novel acceleration technique based on the ANFD method. Only an equivalent equation was established, and a converged consistent solution was obtained to achieve acceleration. In the future, the acceleration performance of the ANFD will be tested. Meanwhile, we will compare the pin-homogenized ANFD method with the pin-homogenized MG CMFD method, which is widely used to accelerate MOC calculations. Our preliminary analysis indicated that the ANFD method achieved an excellent and stable acceleration performance compared with the CMFD method when void regions existed. Other anisotropic cases, such as anisotropic scattering, will be tested in our future studies. Subsequently, Fourier analysis will be used for convergence analysis. Additionally, we will investigate the use of the multilevel acceleration method to reduce the computing burden of ANFD and then attempt to extend this method to 3D cases. References N. Z. Cho , G. S. Lee , C. J. Park , Fusion of method of characteristics and nodal method for 3-d whole-core transport calculation . Trans. Am. Nucl. Soc . 87 , 322 ( 2002 ). Fusion of method of characteristics and nodal method for 3-d whole-core transport calculation H. G. Joo , J. Y. Cho , K. S. Kim et al , Methods and performance of a three dimensional whole-core transport code DeCART . PHYSOR 2004-The Physics of Fuel Cycles and Advanced Nuclear Systems , Chicago, Illinois, US. ( 2004 ) Methods and performance of a three dimensional whole-core transport code DeCART J. Y. Cho , H. G. Joo , Solution of the C5G7MOX benchmark three-dimensional extension problems by the DeCART direct whole core calculation code . Prog. Nucl. Energy 48 , 456 ( 2006 ). https://doi.org/10.1016/j.pnucene.2006.01.006 https://doi.org/10.1016/j.pnucene.2006.01.006 Solution of the C5G7MOX benchmark three-dimensional extension problems by the DeCART direct whole core calculation code Y. S. Jung , C. B. Shim , C. H. Lim et al ., Practical numerical reactor employing direct whole core neutron transport and subchannel thermal/hydraulic solvers . Ann. Nucl. Energy 62 , 357 ( 2013 ). https://doi.org/10.1016/j.anucene.2013.06.031 https://doi.org/10.1016/j.anucene.2013.06.031 Practical numerical reactor employing direct whole core neutron transport and subchannel thermal/hydraulic solvers B. Collins , S. Stimpson , B. W. Kelley et al ., Stability and accuracy of 3D neutron transport simulations using the 2D/1D method in MPACT . J. Comput. Phys . 326 , 612 ( 2016 ). https://doi.org/10.1016/j.jcp.2016.08.022 https://doi.org/10.1016/j.jcp.2016.08.022 Stability and accuracy of 3D neutron transport simulations using the 2D/1D method in MPACT Z. Liu , L. Liang , J. Chen et al ., Accuracy and analysis of the iteratively coupling 2D/3D method for the three-dimensional transport calculations . Ann. Nucl. Energy 113 , 130 ( 2018 ). https://doi.org/10.1016/j.anucene.2017.10.020 https://doi.org/10.1016/j.anucene.2017.10.020 Accuracy and analysis of the iteratively coupling 2D/3D method for the three-dimensional transport calculations J. Choe , S. Choi , P. Zhang , et al ., Verification and validation of STREAM/RASTK for PWR analysis . Nucl. Eng. Tech . 51 ( 2 ), 356 ( 2019 ). https://doi.org/10.1016/j.net.2018.10.004 https://doi.org/10.1016/j.net.2018.10.004 Verification and validation of STREAM/RASTK for PWR analysis K. S. Smith . Spatial Homogenization Methods for Light Water Reactors. Doctoral thesis , Massachusetts Institute of Technology , 1980 . Spatial Homogenization Methods for Light Water Reactors. Doctoral thesis A. Zhu , M. Jarrett , Y. Xu , et al ., An optimally diffusive coarse mesh finite difference method to accelerate neutron transport calculations . Ann. Nucl. Energy 95 , 116 ( 2016 ). https://doi.org/10.1016/j.anucene.2016.05.004 https://doi.org/10.1016/j.anucene.2016.05.004 An optimally diffusive coarse mesh finite difference method to accelerate neutron transport calculations S. Yuk , N. Z. Cho , Two-level convergence speedup schemes for p-CMFD acceleration in neutron transport calculation . Nucl. Sci. Eng . 188 , 1 ( 2017 ). https://doi.org/10.1080/00295639.2017.1332891 https://doi.org/10.1080/00295639.2017.1332891 Two-level convergence speedup schemes for p-CMFD acceleration in neutron transport calculation Y. Xu , T. J. Downar , Convergence Analysis of a CMFD Method Based on Generalized Equivalence Theory . PHYSOR 2012 , Knoxville, Tennessee , April 15-20. (2012) Convergence Analysis of a CMFD Method Based on Generalized Equivalence Theory D. Wang , S. Xiao . A linear prolongation approach to stabilizing CMFD . Nucl. Sci. Eng . 190 , 45 ( 2018 ). https://doi.org/10.1080/00295639.2017.1417347 https://doi.org/10.1080/00295639.2017.1417347 A linear prolongation approach to stabilizing CMFD X. Chai , D. Yao , K. Wang et al ., Generalized coarse-mesh finite difference acceleration for method of characteristics in unstructured meshes . Chinese Journal of Computational Physics , 27 ( 4 ), 541 ( 2010 ). Generalized coarse-mesh finite difference acceleration for method of characteristics in unstructured meshes C. Hao , L. Kang , Y. Xu , et al ., 3D whole-core neutron transport simulation using 2D/1D method via multi-level generalized equivalence theory based CMFD acceleration . Ann. Nucl. Energy 122 , 79 ( 2018 ). https://doi.org/10.1016/j.anucene.2018.08.014 https://doi.org/10.1016/j.anucene.2018.08.014 3D whole-core neutron transport simulation using 2D/1D method via multi-level generalized equivalence theory based CMFD acceleration Z. Liu , X. Zhou , L. Cao , et al ., A new three-level CMFD method based on the loosely coupled parallel strategy . Ann. Nucl. Energy 145 , 107529 ( 2020 ). https://doi.org/10.1016/j.anucene.2020.107529 https://doi.org/10.1016/j.anucene.2020.107529 A new three-level CMFD method based on the loosely coupled parallel strategy L. Kang , C. Hao , Q. Zhao , et al ., Two-level time-dependent GET based CMFD acceleration for 3D whole core transient transport simulation using 2D/1D method . Ann. Nucl. Energy 142 , 107405 ( 2020 ). https://doi.org/10.1016/j.anucene.2020.107405 https://doi.org/10.1016/j.anucene.2020.107405 Two-level time-dependent GET based CMFD acceleration for 3D whole core transient transport simulation using 2D/1D method A. Zhu . Transient methods for pin-resolved whole-core neutron transport. Doctoral thesis , The University of Michigan , 2016 . Transient methods for pin-resolved whole-core neutron transport. Doctoral thesis M. Jarrett , B. Kochunas , A. Zhu , et al ., Analysis of stabilization techniques for CMFD acceleration of neutron transport problems . Nucl. Sci. Eng . 184 , 208 ( 2016 ). https://doi.org/10.13182/NSE16-51 https://doi.org/10.13182/NSE16-51 Analysis of stabilization techniques for CMFD acceleration of neutron transport problems Q. Shen , Y. Xu , T. Downar , Stability analysis of the CMFD scheme with linear prolongation . Ann. Nucl. Energy 129 , 298 ( 2019 ). https://doi.org/10.1016/j.anucene.2019.02.011 https://doi.org/10.1016/j.anucene.2019.02.011 Stability analysis of the CMFD scheme with linear prolongation M. A. Smith , E. E. Lewis , C. B. Na , Benchmark on deterministic 2-D MOX fuel assembly transport calculations without spatial homogenization . Prog. Nucl. Energy 45 , 107 ( 2004 ). https://doi.org/10.1016/j.pnueene.2004.09.003 https://doi.org/10.1016/j.pnueene.2004.09.003 Benchmark on deterministic 2-D MOX fuel assembly transport calculations without spatial homogenization Y. Saad , Iteration Methods for Sparse Linear systems . ( PWS Publishing Company , Boston , 1996 ), pp. 245 - 247 . Iteration Methods for Sparse Linear systems C. Hao , Y. Xu , T. Downar , Multi-level Coarse Mesh Finite Difference Acceleration with Local Two-Node Nodal Expansion Method . Ann. Nucl. Energy 116 , 105 ( 2018 ). https://doi.org/10.1016/j.anucene.2018.02.002 https://doi.org/10.1016/j.anucene.2018.02.002 Multi-level Coarse Mesh Finite Difference Acceleration with Local Two-Node Nodal Expansion Method J. I. Yoon , H. G. Joo , Two-level coarse mesh finite difference formulation with multigroup source expansion nodal kernels . J. Nucl. Sci. Technol . 45 , 668 ( 2008 ). Two-level coarse mesh finite difference formulation with multigroup source expansion nodal kernels Z. Zhong , T. Downar , Y. Xu , et al ., Implementation of two-level coarse-mesh finite difference acceleration in an arbitrary geometry, two-dimensional discrete ordinates transport method . Nucl. Sci. Eng . 158 ( 3 ), 289 ( 2008 ). https://doi.org/10.13182/NSE06-24TN https://doi.org/10.13182/NSE06-24TN Implementation of two-level coarse-mesh finite difference acceleration in an arbitrary geometry, two-dimensional discrete ordinates transport method Y Xu , C. Hao , A novel and efficient hybrid RSILU preconditioner for the parallel GMRES solution of the coarse mesh finite difference equations for practical reactor simulations . Nucl. Sci. Eng . 194 ( 2 ), 104 ( 2019 ). https://doi.org/10.1080/00295639.2019.1657322 https://doi.org/10.1080/00295639.2019.1657322 A novel and efficient hybrid RSILU preconditioner for the parallel GMRES solution of the coarse mesh finite difference equations for practical reactor simulations A. Yamamoto , M. Tabuchi , N. Sugimura , et al ., Derivation of optimum polar angle quadrature set for the method of characteristics based on approximation error for the Bickley function . J. Nucl. Sci. Technol . 44 ( 2 ), 129 ( 2007 ). https://doi.org/10.1080/18811248.2007.9711266 https://doi.org/10.1080/18811248.2007.9711266 Derivation of optimum polar angle quadrature set for the method of characteristics based on approximation error for the Bickley function G.A. Freund , P. Elias , D. R. MacFarlane , J. D. Geier , J.F. Boland , “Design Summary Report on the Transient Reactor Test Facility (TREAT)”, ANL-6034 , Argonne National Laboratory , ( 1960 ). “Design Summary Report on the Transient Reactor Test Facility (TREAT)”, ANL-6034 Z. Wu , D. Lin , D. Zhong , The design features of the HTR-10 , Nucl. Eng. Des . 218 ( 3 ), 25 ( 2002 ). https://doi.org/10.1016/S0029-5493(02)00182-6 https://doi.org/10.1016/S0029-5493(02)00182-6 The design features of the HTR-10 Y. Xu , V. Seker , T. Downar , Quasi-diffusion method with 3-D cross sections for TREAT core analysis , Nucl. Technol . 206 ( 6 ), 825 , ( 2020 ). https://doi.org/10.1080/00295450.2019.1672451 https://doi.org/10.1080/00295450.2019.1672451 Quasi-diffusion method with 3-D cross sections for TREAT core analysis 10.1007/s41365-020-00834-2 This work was supported by the National Key R&D Program of China (No. 2018YFE0180900). 2020-07-24 2020-10-13 2020-10-16 2020-12-07 2020-12-01 2021-04-10 © 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. 2020 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-020-00834-2 http://dx.doi.org/10.1007/s41365-020-00834-2 . 2020 Nuclear Science and Techniques 12 31