Methodology for Discontinuity Factors Generation for Simplified P 3 Solver Based on Nodal Expansion Formulation

: The Simplified Spherical Harmonic (SP N ) approximation was first introduced as a three ‐ dimensional (3D) extension of the plane ‐ geometry Spherical Harmonic (P N ) equations. A third order SP N (SP 3 ) solver, recently implemented in the Nodal Expansion Method (NEM), has shown promis ‐ ing performance in the reactor core neutronics simulations. This work is focused on the develop ‐ ment and implementation of the transport ‐ corrected interface and boundary conditions in an NEM SP 3 solver, following recent published work on the rigorous SP N theory for piecewise homogeneous regions. A streamlined procedure has been developed to generate the flux zero and second or ‐ der/moment discontinuity factors (DFs) of the generalized equivalence theory to minimize the error introduced by pin ‐ wise homogenization. Moreover, several colorset models with varying sizes and configurations are later explored for their capability of generating DFs that can produce results equivalent to that using the whole ‐ core homogenization model for more practical implementations. The new developments are tested and demonstrated on the C5G7 benchmark. The results show that the transport ‐ corrected SP 3 solver shows general improvements to power distribution prediction compared to the basic SP 3 solver with no DFs or with only the zeroth moment DF. The complete equivalent calculations using the DFs can almost reproduce transport solutions with high accuracy. The use of equivalent parameters from larger size colorset models show a slightly reduced predic ‐ tion error than that using smaller colorset models in the whole ‐ core calculations.


Introduction
Obtaining solutions to the neutron transport equation with consistent angular discretization, such as discrete ordinates (S N ) or spherical harmonics (P N ), for the threedimensional (3D) transport problem can be challenging, even with the rapid increase in computing power.Thus, low-order approximations to the transport equation that can be solved at a significantly reduced computational cost are of special interest.The coarse mesh diffusion equations are traditionally the first choice to the program developer; however, such methods may be inadequate for advanced reactor designs, which usually demonstrate a high level of heterogeneity, or in core regions near material boundaries and strong absorbers.
The Simplified P N (SP N ) equations were first proposed to introduce additional transport effects into the standard P 1 equations without introducing the complexities and undesired increase in the computational cost that a full transport theory solution would entail.The P N equations are obtained from inserting the truncated spherical harmonics expansion (to some order n) of the angular flux and differential scattering cross section into the transport equation.
By replacing the derivatives of the one-dimensional (1D) P N equation in planar geometry with the 3D gradient and divergence operators, one arrives at a 3D generalization of Energies 2021, 14, 6478 2 of 17 the 1D P N equations [1].Such a substitution was performed in an ad hoc manner, leaving the method to lack a theoretical foundation.As asymptotic and variational derivations later assisted to establish the theoretical basis for the SP N method [2][3][4], the application of the method became more popular.The numerical results show that the SP N approximation can yield accurate solutions of the transport problems, which outperforms the diffusion method, but with considerably less computational expense than those for the S N or the full P N methods [5,6].For example, the computational time for SP 3 is only about twice that of the diffusion method.At locations in the core with a high flux gradient, such as material boundaries, the SP 3 has a higher precision due to the use of the higher Legendre moments P N scattering library and angular flux distribution.It is suggested that the SP 3 method should be utilized to a homogenized pin cell (other than a large assembly node) level calculation using a few groups (instead of two groups) to retain its superiority in accuracy over diffusion approximation [7].
Since spatial homogenization is applied in pin homogenized fine-mesh calculations, a mitigation method is required to eliminate or reduce the homogenization error.In this regard, discontinuity factors based on the generalized equivalence theory (GET) have been investigated, such as [8,9].Kozlowski et al. proposed the pin cell discontinuity factors (CDFs) to recover the error introduced by pin cell homogenization [6,10].The method showed the potential to improve the accuracy of the pin power prediction; however, its application for practical core problems was limited because of the considerable amount of data required for whole core calculations.Yamamoto developed a practical estimation method of the second angular moment φ 2 and the third angular moment J 3 from the results of method of characteristic (MOC), based on which two types of DFs are proposed [11].The independent DF preserves angular moments of J 1 and J 3 using two DFs at each surface, while the common DF preserves J 1 using one DF at each surface to reduce the memory storage for DF.There were also attempts to apply DFs to interfaces between assemblies in addition to using super homogenization (SPH) corrected cross sections in practical reactor core calculations, such as in [12].All studies have shown improved accuracy with varying degrees from the pin-by-pin SP 3 calculation without DFs.
The NEM (Nodal Expansion Method) code has been developed for 3D whole-core steady state and transient analysis with cartesian, hexagonal, and cylindrical geometries.A SP 3 solver was established in the NEM by taking advantage of the original diffusion-based solver to achieve a pin cell resolution of highly heterogeneous reactor cores [13].Further development and enhancement have since been conducted to improve its performance, such as incorporating the higher order scattering cross sections and discontinuity factors (DFs) [14,15].However, this implementation adopted the ad hoc interface and boundary conditions based on the assumption of 1D behavior near a surface, which prevents the angular flux from being represented by the SP N flux components.Therefore, an SP N solver cannot provide the exact scalar flux solution as the P N solver.
This work is focused on the development and implementation of the interface and boundary conditions in an NEM SP 3 solver, following recent published work on the rigorous SP N theory for piecewise homogeneous regions [16,17].The new theory proves that the SP N solutions can represent a particular set of the angular flux solution by the P N theory.The resulting interface and boundary conditions now involve terms of higher order gradients of flux as well as tangential gradients, although the SP N equations are identical to the conventional ones.The transported corrected interface and boundary conditions are formulated in the nodal expansion expression of the flux components and the solution method involving the response matrix utilized in the NEM.Then, a side-dependent DF generation approach is devised based on the GET for homogenization with the intention to eliminate the error introduced by pin-wise homogenization.Both zeroth-and secondmoment DFs can be generated using an in-house developed code using the solution from the transport solutions using Method of Characteristics (MOC).Several colorset models with varying sizes and configurations are later explored for their capability of generating DFs that can produce results equivalent to that using the whole-core model for more practical implementations.The new developments are tested and demonstrated on the C5G7 benchmark.
This paper is organized in the following way.In Section 2, the SP 3 equations based on the nodal expansion formulation are derived and the solution method of SP 3 equations using the response matrix is introduced.Section 3 is focused on the derivation of the interface and boundary conditions in the new SP N theory and their implementation in the NEM SP 3 solver.Section 4 introduces the equivalent calculation scheme developed for the SP 3 solver with the emphasis placed on the generation of equivalence parameters using different lattice models to facilitate practical applications.The transport-corrected SP3 solver is tested on a series of problems from the C5G7 benchmark, and the results are shown in Section 5. Finally, conclusions and future perspectives are given in Section 6.

SP 3 Method Based on Nodal Expansion Formulation
For the general case with odd n, the SP N equations can reduce to a (n + 1)-th order differential equation for the lowest order scalar flux, which is equivalent to (n + 1)/2 coupled second order differential equations.This indicates that being diffusion theory type equations, the SP N solvers have commonly been cast in such a way as to leverage existing diffusion machinery by iterating over the SP N moment equations.The steady-state SP 3 equations used in the NEM are the following [13]: where the unknowns are the synthesized zeroth moment flux approximation Φ 0 = φ 0 + 2φ 2 and the second moment flux φ 2 in all energy groups, Σ r is the transport-corrected removal macroscopic cross section, and the source terms on the right-hand side of the equations S 0,g is the sum of the scattering and fission source associated with the zeroth angular moment.
S 0,g = Σ G g =1,g =g Σ s,0,g →g φ 0,g + where Σ s and Σ f are the scattering and fission cross sections, respectively, and χ is the fission spectrum.Note that the scattering matrix Σ s,0,g →g includes only the isotropic scattering source, while the higher moments are approximated by the transport-corrected removal cross section Σ r .The two diffusion coefficients are defined as D 0 = 1/3Σ r,0,g and D 2 = 9/35Σ r,2,g .The Marshak boundary conditions in terms of surface fluxes and incoming (j + ) and outward (j − ) partial currents, as used in the NEM, are as follows: In the NEM, the 3D gradient operator is split into three directions, i.e., ∂y 2 + ∂ 2 ∂z 2 and the same equations are used in all directions.In an arbitrary node with constant neutronic properties and dimensions ∆x, ∆y, and ∆z in the Cartesian geometry, apply the transverse integration approximation, where Equation ( 1) is spatially integrated over the two dimensions transverse to the particular direction of interest.The resulting equations in the x direction take the following form: Energies 2021, 14, 6478 4 of 17 where the transverse integrated zeroth and second moment flux are as follows: Φ 0,g (x, y, z)dzdy and first moment transverse leakage term is as follows: and the transverse integrated source term is as follows: S 0,g (x, y, z)dzdy. (7) Replacing Φ 0 with φ 2 in Equation ( 5) yields the third moment transverse leakage terms L y,3,g (x) and L z,3,g (x).
In the context of the NEM, the intra-node flux moments Φ 0 (x) and φ 2 (x) are expanded in series within each node using fourth-order polynomial basis functions as follows: where Φ 0 and φ 2 are the node average flux moments and the polynomials are as follows: Here, we drop the energy group index g and focus on the mono-energetic expression for simplicity.
The expansion coefficients a and b can be determined by solving Equation (3) at specific locations, i.e., left node edge x = −∆x/2 and right node edge x = ∆x/2.Introducing the node edge flux moments, including Φ 0,L and φ 2,L on the left edge, and Φ 0,R and φ 2,R on the right edge yields the following: Energies 2021, 14, 6478 5 of 17 Next, multiply f 1 and f 2 to Equation (3), perform the transverse integration 1 ∆x ∆x/2 −∆x/2 dx, define the following flux-like terms: and the rest expansion coefficient can be derived as follows: Note that the expression of the flux-like terms can be derived from the 1D SP 3 equations shown in Equation ( 4) by multiplying f 1 and f 2 , and then performing the transverse integration along one direction as follows: where α = Σ r,2 + 4 5 Σ r,0 , and the J n,x,L and J n,x,R are the node edge net currents on the left (L) and right (R) edge, respectively.
In the above, L i,n,k = 1 ∆x ∆x/2 −∆x/2 L i,n f k dx is the transverse integrated transverse leakage term, where i refers to the direction (y or z), n is the flux moments (1 or 3), and k is the index integration function or polynomial (1 or 2).The new source term S 0,k = 1 ∆x ∆x/2 −∆x/2 S 0 f k dx is obtained in the same way.
The calculation mechanism implemented in NEM requires all the quantities of interest to be derived in terms of the partial currents and formulated in the response matrix, which expresses the outgoing partial currents as a function of incoming partial currents as well as intra-node sources and sinks.The response matrix is derived from Fick's Law for the partial currents on the node boundaries, using the x-axis as an example, where j + 1 and j − 1 are the outgoing and incoming partial currents, respectively.Substituting the polynomial expansion in Equation (3) into the equation above and performing differen-Energies 2021, 14, 6478 6 of 17 tiation yields four equations in each of the three directions.Take the right edge of the node x = ∆x/2, for example, the partial currents are as follows: while similar expressions can be obtained for j + 1,L and j + 3,L .Recall that the expansion coefficients a 1 -a 4 and b 1 -b 4 are functions of the nodeaverage flux, node edge fluxes, and flux-like terms, as shown in Equations ( 4)- (7), which can be expressed in terms of the partial currents.The node-average flux Φ 0 and φ 2 required in Equation ( 4) can be obtained by integrating the neutron balance equation (Equation ( 2)) over the node volume.Substituting for the node-average flux and fluxlike terms (Equation ( 7)) in Equation ( 9) yields the four partial current equations in the x-direction.This derivation is straight-forward but the algebra is quite involved, and thus not shown in detail here.More information can be found in [18].Note that all 12 partial current equations in the three directions are coupled on the right-hand side through the leakage transverse leakage and source terms.Rewriting them in the matrix form, we arrive at the following response matrix equation: where, in principle, J + , J − , S, and L are 12 × 1 vectors, and the coefficient matrices are the size of 12 × 12, although the actual implementation differs slightly from this representation [18].The eigenvalue problem is solved using the NEM using the traditional inner/outer iteration method.For each group, inner iterations, or multiple sweeps through the mesh with a known internal source are performed to invert the within-group SP 3 removal matrix.Outgoing partial currents are computed using the incoming partial currents and the node-dependent response functions.These outgoing partial currents become the incoming partial currents in the neighboring nodes.Outer fission source iterations are then performed around the inner iterations to calculate values for the problem multiplicative eigenvalue (k eff ) and the space-and energy-dependent fission neutron source distribution.

Derivations of Interface and Boundary Conditions
Homogenization techniques have been used in nuclear reactor analysis to reduce the spatial and angular domain complexity of a nuclear reactor by replacing pre-calculated heterogeneous subdomains by homogeneous ones and using a low order solver to solve the homogeneous problem.To guarantee this equivalence, especially for the fission rate and leakage rate in neutron balance equation, the following two main factors should be guaranteed: (1) homogeneous cross sections are calculated based on the spectrum of neutron fluxes to preserve the fission rate; (2) at node interfaces, first order flux (J 1 ) should be preserved to make a consistent leakage rate between heterogeneous and homogeneous calculations.
To achieve the second objective, DFs should be generated at the boundaries of the homogenized nodes as defined in the GET [19], which was derived from the equivalence theory [20].However, performing the complete equivalent calculation using the DF for SP 3 has long been problematic due to the ad hoc interface and boundary conditions adopted in the early developments of the SP N theory, which involves the diffusion theory type of first order derivatives in the surface normal direction.It prevents the angular flux from being represented by the SP 3 solution in the conventional SP 3 theory or being compared to the reference angular flux solution from the heterogeneous transport calculation.Simply put, it is not clear which higher order quantities should be continuous across the node surface, which makes it difficult to rigorously define and calculate the high order DF in SP 3 .
In a recent work, Chao rederived the SP N equations from the P N , proved that φ 0 solutions in 3D SP 3 are one set of solution to 3D P 3 theory, and provided an exact expression of the angular flux using the SP 3 flux solutions [16,17].This method has been adopted, modified, and implemented in the NEM SP 3 solver.For the convenience of the readers, we provide a concise review summary of Chao's work regarding the derivation of the transport-corrected interface and boundary conditions.
The derivation starts by decomposing the transport equation is into an even parity part and an odd parity part as follows: It was shown that the SP N equations can be derived via the variational method by introducing the following trial function into the even parity transport equation (Equation ( 18)): The operator ] and the coefficients a n,k are the same as those in Legendre polynomial P n .The functions F n are called auxiliary functions with the property as follows: Note that substituting Equation ( 21) into the SPn equation will yield the expression of the auxiliary function F n .By defining a dimensionless "unit vector" ∇ = ∇ ∇ 2 −1/2 , the operator L n (Ω, ∇) can be rewritten in terms of Legendre polynomials as follows: Therefore, the trial function Equation ( 19) can be rewritten as follows: Equation ( 23) is important not only because it can lead to the exact SPn equations (by plugging it into Equation ( 18)), but also because it gives the angular flux for the SPn theory because Equation ( 23) can be reconstructed in terms of the SP N solution functions, as shown below.
To derive the boundary and interface conditions, use the angular flux on a node surface.The angular partial currents going out or in through a surface with the normal vector n are defined as follows: (24) To calculate the k-th even order moment of the partial currents, Equation ( 22) is multiplied by the k-th even order Legendre polynomial of the cosine of the polar angle with respect to the normal vector n, and is then integrated over the angular space as follows: Replacing the even parity angular flux by Equation ( 22) and making use of the relationship in Equation ( 19) leads to the following: and It should be noted that Ψ k is the k-th moment of the generic even parity flux, not the n-th scalar moment φ n in the SPn equation.Now, one can derive the exact expression of Ψ k and J k using F n .
In the special case of SP3, we have the following terms on the y-z surface: By comparing with the conventional SP3 theory, the new theory contains additional terms involving the part of the divergence in the tangential directions of the surface.The continuity of Equations ( 30) and (32) imply that the following two quantities should be continuous on the node surface: (33) It is worth mentioning that the two terms φ 2 and J 3 are different from the usual ad hoc continuity conditions for φ 2 and J 3 , which will no longer be continuous on the surface in the new theory.

Implementation of Interface and Boundary Conditions
Converting the tangential divergence to the normal component of the divergence by applying the relationship 33) and (34), we arrive at the following interface and boundary conditions for the SP 3 theory: Energies 2021, 14, 6478 9 of 17 where (37) Again, the δ term above is a transport correction to the traditional ad hoc conditions and it represents the effect of the transverse terms.
Using Equation ( 26) and the definition of transport corrected terms in Equations ( 33) and (34), Equations ( 29) and (30) become the following: Note that J 3 = j + 3 − j − 3 and the corrected partial currents j + 3 and j − 3 are different from the previously defined j + 3 and j − 3 .Rearranging the above equations yields the interface and boundary conditions expressed by the partial currents required by the following response matrix: Since the nodal expansion in the base SP 3 equations is not impacted during the derivation in Section 3.1 and, thus, still valid for the new SP 3 theory, the corrected second moment flux term and third moment net partial current become the following: Take the right node surface in the x-direction for example, the corresponding node interface flux and net current terms as appeared in Equations ( 13) and ( 14) can be derived as the following: Follow the same procedure described in Section 2, a set of equations similar to the original response matrix based on the partial currents on the node surface can be derived.We can use the same inner/outer iteration scheme to solve the response matrix and the corrected terms φ 2 and J 3 will be obtained in addition to Φ 0 and J 1 .

Equivalent Calculation Scheme
As shown previously, the nodal SP 3 solver performs the finite volume integration, which requires material parameters to be constant in the node.Fuel pin homogenized cross sections are prepared using the flux-volume weighting procedure to preserve reaction rates between lattice calculations and SP 3 calculations.Here, we will only focus on the DF generation approach developed and implemented for the SP 3 solver aiming to preserve the leakage rate between nodes.

SP 3 Discontinuity Factors
The derivation above shows that the DFs needed to be used in the homogenized nodes can be generated without ambiguity by using the particular SP N angular flux solution.The work presented below is similar to a previous effort [11] in the sense that the DFs are generated for both zeroth and second order flux moments; however, angular moments in the SP 3 method are different and so is the application of DFs.It should be emphasized that the DF is defined for and applied to Ψ 0 and Ψ 2 in the following way: instead of Φ 0 and φ 2 , the unknowns in the SP 3 equations.In other words, given a surface limiting two adjacent homogenized regions, the DFs enforce the continuity for the heterogeneous reconstructed flux DF − 0 Ψ hom,− 0 The flux solutions from both heterogeneous and homogeneous models are required to compute DFs.The former can be obtained from transport calculations, while the latter needs to be calculated by imposing the conditions of conserving the total net current from the reference heterogeneous transport calculation.It can be seen in Equations ( 29)-(32) that the reference value of Ψ 0 , Ψ 2 , J 1 , and J 3 can all be readily computed in the transport calculation.The resulting values of J 1 and J 3 will then be imposed on each surface of a node (pin cell) as the boundary condition to uniquely determine the SP 3 solution (Ψ hom 0 and Ψ hom 2 ) inside the homogenized cell.In the context of a response matrix in the NEM, the updated incoming currents must be adjusted by the net current obtained in the transport solution, that is, by adding the net current to the outgoing current of the neighboring nodes.This leads to solving a fixed interface problem using the homogenized cross sections while also fixing the eigenvalue, which is similar to the calculation scheme described in [21].For this purpose, an independent program DF Generation Routine (DFGR) is developed to generate the homogeneous solutions.It can take the eigenvalue, cell interface net currents, and cell averaged scalar fluxes from the transport calculation and solve the fixed interface problem.Unlike the eigenvalue problem, no interactions between nodes are assumed here and each pin cell can be solved individually.
As depicted in Figure 1, the equivalent calculation process devised for the SP 3 solver can be carried out in the following steps: Energies 2021, 14, 6478 11 of 17 pendent program DF Generation Routine (DFGR) is developed to generate the homogeneous solutions.It can take the eigenvalue, cell interface net currents, and cell averaged scalar fluxes from the transport calculation and solve the fixed interface problem.Unlike the eigenvalue problem, no interactions between nodes are assumed here and each pin cell can be solved individually.As depicted in Figure 1, the equivalent calculation process devised for the SP3 solver can be carried out in the following steps:

Practical Approach to Generate SP3 Discontinuity Factors
The equivalent calculation procedure presented above relies on the whole-core transport solutions, which is not practical in routine calculations and, in a way, defeats the purpose of using low order solvers such as SP3.Therefore, efforts are also made to

1.
Perform the transport calculations and generate quantities including k eff , cell homogenized cross sections, cell averaged scalar flux, and side-dependent surface fluxes and currents.

2.
Solve the fixed interface problem for each cell sequentially by taking the reference values obtained in step one, generate the SP 3 surface fluxes and currents, then compute the DFs according to their definitions.

3.
Execute a normal SP3 calculation on a pin-by-pin level using the homogenized cross sections from step one and DFs from step two to obtain the SP 3 solution.

Practical Approach to Generate SP 3 Discontinuity Factors
The equivalent calculation procedure presented above relies on the whole-core transport solutions, which is not practical in routine calculations and, in a way, defeats the purpose of using low order solvers such as SP 3 .Therefore, efforts are also made to explore the feasibility of colorset lattice models to generate equivalent parameters for the core simulation.The three colorset models under consideration are as follows: 1.
Single-pin model: cross sections are homogenized over the pin cell and DFs are calculated for each of the four surfaces of a pin cell.

2.
Double-pin model: two pin cells of different type (i.e., material) are placed next to each other and DFs are calculated for each of the four surfaces in each pin cell.The cross sections are taken from the first model.

3.
Assembly model: both homogenized cross sections and DFs are location dependent, i.e., they are generated for each of the pin cells.
The reflective boundary condition is applied to all three models.Among the three models, the first one is the simplest one and does not capture the leakage between different fuel pins and, thus, is expected to result in a worse performance.The second one captures the effect from neighboring effects to a certain extend but the total number of models increases significantly with the types of fuel pin.The last model supports location-dependent pin-wise homogenizations and DF generation in an environment that is similar to that in the core.

Verification of Transported Corrected SP 3 Solver
The newly developed transported corrected SP 3 solver and the corresponding equivalent calculation scheme are verified using the C5G7 benchmark [22].It is a miniature light water reactor (LWR) with sixteen fuel assemblies (mini core): eight uranium oxide (UO 2 ) assemblies and eight mixed oxide (MOX) assemblies, surrounded by a water reflector.It features a quarter-core radial symmetry in the 2D configuration, as shown in Figure 2. On the fuel pin level, there are one UO 2 pin, three MOX pins, one guide tube, and one fission chamber.The three MOX fuels pins are 4.3, 7.0, and 8.7% plutonium weight enriched.
In this study, only the core geometry is adopted, and a new set of cross sections is generated in a seven-group energy structure because the original cross sections prepared by the benchmark do not satisfy the unique requirement of the SP 3 solver [15].
The transport code OpenMOC [23] is selected to perform the reference transport calculation because of the convenience of the code to extract node interfaces angular fluxes for DFs generations.However, because the current version only accepts isotropic scattering cross sections, transport corrections are applied to modify isotropic scattering cross sections.All the transport calculations shown below are conducted using the long characteristic method with 32 azimuthal angles, 3 polar angles, the ray trace width of 0.03 cm, and the eigenvalue convergence criterion of 1 × 10 −7 .
alent calculation scheme are verified using the C5G7 benchmark [22].It is a miniature light water reactor (LWR) with sixteen fuel assemblies (mini core): eight uranium oxide (UO2) assemblies and eight mixed oxide (MOX) assemblies, surrounded by a water reflector.It features a quarter-core radial symmetry in the 2D configuration, as shown in Figure 2. On the fuel pin level, there are one UO2 pin, three MOX pins, one guide tube, and one fission chamber.The three MOX fuels pins are 4.3, 7.0, and 8.7% plutonium weight enriched.In this study, only the core geometry is adopted, and a new set of cross sections is generated in a seven-group energy structure because the original cross sections prepared by the benchmark do not satisfy the unique requirement of the SP3 solver [15].
The transport code OpenMOC [23] is selected to perform the reference transport calculation because of the convenience of the code to extract node interfaces angular fluxes for DFs generations.However, because the current version only accepts isotropic scattering cross sections, transport corrections are applied to modify isotropic scattering cross sections.All the transport calculations shown below are conducted using the long characteristic method with 32 azimuthal angles, 3 polar angles, the ray trace width of 0.03 cm, and the eigenvalue convergence criterion of 1 10 .

Transport-Corrected SP3 Method
First, we compare the performance of the base-and transport-corrected SP3 solver without applying DFs to reveal the impact of the updated interface and boundary conditions.Three cases are selected from the C5G7 benchmark for this purpose, including the

Transport-Corrected SP 3 Method
First, we compare the performance of the base-and transport-corrected SP 3 solver without applying DFs to reveal the impact of the updated interface and boundary conditions.Three cases are selected from the C5G7 benchmark for this purpose, including the UO 2 assembly, MOX assembly, and C5G7 core.In the single assembly cases, pin-wise cross sections are generated using the assembly model with the infinite lattice approximation.For the last case, they are prepared in a whole-core transport calculation, including the cross sections of the water reflector.
The comparison of the eigenvalue yields negligible differences (less than 10 pcm in all cases), which indicates the correction term δ in Equation ( 37) is small and, thus, has a trivial impact on integral parameters such as k eff .The focus has been shifted to local quantities.The absolute relative error in the pin power at all locations is first computed for both solvers as d = |p SP3 /p MOC − 1|, where p denotes the normalized pin power.Then, the difference between the two solvers is calculated as ∆ = d t − d b , where the subscript t means "transport corrected" and b refers to "base".The value of ∆ is plotted for the C5G7 core problem according to the pin location, as shown in Figure 3.
The blue regions in the distribution indicate the improved prediction accuracy due to the correction to the interface and boundary conditions.It can be seen that the transportcorrected SP 3 solver yields slightly better pin power distributions, especially at the center of the UO 2 assembly (upper left assembly), regions close to the water reflector (right and bottom surfaces), as well as pin cells next to the fission chamber and guide tubes.As expected, the transport-corrected solver cannot eliminate the prediction error because the spatial homogenization errors still exist.
trivial impact on integral parameters such as  .The focus has been shifted to local quan tities.The absolute relative error in the pin power at all locations is first computed for both solvers as  = |  ⁄ − 1|, where p denotes the normalized pin power.Then, the dif ference between the two solvers is calculated as Δ =  −  , where the subscript t mean "transport corrected" and b refers to "base".The value of Δ is plotted for the C5G7 core problem according to the pin location, as shown in Figure 3.The blue regions in the distribution indicate the improved prediction accuracy due to the correction to the interface and boundary conditions.It can be seen that the transport corrected SP3 solver yields slightly better pin power distributions, especially at the cente of the UO2 assembly (upper left assembly), regions close to the water reflector (right and bottom surfaces), as well as pin cells next to the fission chamber and guide tubes.As ex pected, the transport-corrected solver cannot eliminate the prediction error because the spatial homogenization errors still exist.

Test and Comparison of DFs
Next, the performance of DFs is tested and analyzed in full equivalent calculations In addition to the three cases used in the previous sub-section, a C5 core problem is also introduced.It contains only fuel assemblies with reflective boundary conditions applied to the outer surface of the geometry.This case would help reveal the impact of DFs on th assembly interfaces by eliminating the effects of the reflector.Following the reference transport calculation, the procedure depicted in Figure 1 is carried out to generate equiv alent parameters, inducing the pin-wise cross sections and side-dependent DFs.The equivalent parameters are generated for and applied to each case independently.
It is worth mentioning that numerical instability emerges in regions experiencing a high flux gradient after applying DFs to the NEM response matrix, which forces the spatia discretization to be further refined to achieve convergence.As a result, each pin cell is sub divided into four nodes in equal size and the DFs are only imposed on the outer surface o the pin cell.The refinement allows the flux gradient to occur inside the pin cell; however

Test and Comparison of DFs
Next, the performance of DFs is tested and analyzed in full equivalent calculations.In addition to the three cases used in the previous sub-section, a C5 core problem is also introduced.It contains only fuel assemblies with reflective boundary conditions applied to the outer surface of the geometry.This case would help reveal the impact of DFs on the assembly interfaces by eliminating the effects of the reflector.Following the reference transport calculation, the procedure depicted in Figure 1 is carried out to generate equivalent parameters, inducing the pin-wise cross sections and side-dependent DFs.The equivalent parameters are generated for and applied to each case independently.
It is worth mentioning that numerical instability emerges in regions experiencing a high flux gradient after applying DFs to the NEM response matrix, which forces the spatial discretization to be further refined to achieve convergence.As a result, each pin cell is subdivided into four nodes in equal size and the DFs are only imposed on the outer surface of the pin cell.The refinement allows the flux gradient to occur inside the pin cell; however, the flux is assumed to be constant across surfaces within each pin cell when the DFs are generated.This inconsistency will inevitably introduce errors to the SP 3 calculation.
Table 1 compares the performance of the transport-corrected SP 3 solver in terms of its prediction of eigenvalue and pin power distribution.It can be seen that the DFs of the GET help reduce the prediction error significantly in all four test cases.The largest improvement is observed in the C5G7 core case, where the deviation from the reference value is reduced by ~200 pcm in eigenvalue and over halved in the root-mean-square (RMS) error of pin power distribution.The prediction accuracy is also drastically improved in the MOX assembly case where the local flux gradient is more profound than that in the UO 2 assembly.In principle, if the equivalent calculation is conducted properly, the low order operator should reproduce the transport solutions.The non-zero deviations from the reference solution, as shown in Table 1, indicate that the inconsistency between the DF generation and application introduced by the refined 2 × 2 meshing prevent the full equivalent calculation from being performed.Another source of deviation in the last test case stems from the fact that DFs for the nodes in the water reflector are all assumed to be one due to the limitation of the DFGR program.It would tilde the flux distribution in regions across the core/reflector boundary.
Next, various colorset models are tested for their capability to facilitate practical equivalent calculations, because generating equivalent parameters from the whole-core transport solution defeats the purpose of homogenization.The three colorset models described previously are as follows: 1.
Single-pin model: There are seven models each corresponding to one type of pins.
The non-fissionable node, such as guide tubes and water reflector cells, is placed in the center of a 3 × 3 configuration surrounded by UO 2 fuel pins.

2.
Double-pin model: Eight sets of 2 × 1 pin colorset models are developed for different combination of pin cells.These DFs will be used in the whole core on the interface between neighboring different cells.

3.
Assembly model: Single UO 2 and MOX assembly models.
A fourth model is also considered, which is a mix of options two and three.The DFs from the assembly colorset are used on node interfaces inside an assembly, while those from the double-pin colorset are used on the assembly interfaces.The comparison of the results of the C5G7 core calculation is listed in Table 2.We observe that the calculation accuracy in terms of the eigenvalue and pin power distribution is improved slightly and consistently as the size of the colorset model increases.The assembly model shows improved results compared to the single-pin and double-pin models and the addition of DFs on the assembly interfaces from the double-pin model (option four) further reduces the prediction error.It is expected that the use of DFs from the single-pin and double-pin models would not yield noticeably improved solutions, if any, compared with the results obtained without DFs, because DFs and associated cross sections are only applicable for the environment (spectrum and boundary conditions) within which they are generated.However, it is somewhat unexpected that none of the colorset models significantly outperform the approach without DFs in Table 1.That is to use only the location-dependent cross sections generated from the whole-core reference transport solution, although the associated computational cost in the latter case would be massively elevated.Figure 4 shows the relative error of the predicted pin power in the C5G7 whole-core calculation using the cross sections and DFs generated in option three.It can be seen that the main discrepancies occur in regions near the interface between the active core and reflector, where the SP 3 solution largely underestimates the fission power and neutron flux, while a uniform overestimation is observed in the central UO 2 assembly.An almost identical trend is also found in the comparison of the SP3 solution with option four against the transport results, which is expected because the DFs for nodes along the core/reflector interface are missing in both options three and four.
Assembly + double pin 272 0.047 We observe that the calculation accuracy in terms of the eigenvalue and pin power distribution is improved slightly and consistently as the size of the colorset model increases.The assembly model shows improved results compared to the single-pin and double-pin models and the addition of DFs on the assembly interfaces from the doublepin model (option four) further reduces the prediction error.
It is expected that the use of DFs from the single-pin and double-pin models would not yield noticeably improved solutions, if any, compared with the results obtained without DFs, because DFs and associated cross sections are only applicable for the environment (spectrum and boundary conditions) within which they are generated.However, it is somewhat unexpected that none of the colorset models significantly outperform the approach without DFs in Table 1.That is to use only the location-dependent cross sections generated from the whole-core reference transport solution, although the associated computational cost in the latter case would be massively elevated.Figure 4 shows the relative error of the predicted pin power in the C5G7 whole-core calculation using the cross sections and DFs generated in option three.It can be seen that the main discrepancies occur in regions near the interface between the active core and reflector, where the SP3 solution largely underestimates the fission power and neutron flux, while a uniform overestimation is observed in the central UO2 assembly.An almost identical trend is also found in the comparison of the SP3 solution with option four against the transport results, which is expected because the DFs for nodes along the core/reflector interface are missing in both options three and four.The results above imply that the impact of two factors in the equivalent calculation, the homogenized cross sections, and DFs, almost contribute equally to the prediction accuracy with regard to the eigenvalue and power distribution.The largest impact is found in regions in the vicinity of material boundaries, especially the core boundary, suggesting The results above imply that the impact of two factors in the equivalent calculation, the homogenized cross sections, and DFs, almost contribute equally to the prediction accuracy with regard to the eigenvalue and power distribution.The largest impact is found in regions in the vicinity of material boundaries, especially the core boundary, suggesting that colorset models including non-multiplying regions must be properly developed and tested.Effort should also be made in the future to determine the optimal approach to perform equivalent calculations for the SP 3 solver.

Conclusions and Outlook
In this work, we focus on the development of a transport-corrected nodal SP 3 solver for the equivalent reactor core calculation.The response matrix in the NEM SP 3 solver has been reformulated based on the recently published new SP 3 theory with rigorously derived interface and boundary conditions, while the computation scheme in the NEM is kept intact.A streamlined process of generating DFs of the GET for the correction of pinby-pin homogenization error has been implemented, which incorporates the high-fidelity transport code and an in-house developed DF generation program.We also propose a few models with varying sizes for the generation of equivalent parameters aiming to further reduce the computational burden and explore their feasibility to practical applications.The transport-corrected SP 3 solver and various colorset models are tested on the mini-core C5G7 benchmark problem.
It was found that the transport-corrected SP 3 solver can compute the inter-node leakage rate correctly due to the updated interface and boundary conditions and, thus, improve the prediction of the pin power distribution.In the whole-core calculation, the regions benefiting the largest improvement are located at the center of the core, close to the water reflector, and next to the fission chamber and guide tubes.
The performance of DFs is tested and analyzed in full equivalent calculations.We can see that the DFs of the GET can reduce the prediction error significantly in all the test cases from the single assembly, active core, and core with reflectors.In the last case, the deviation from the reference transport solutions is reduced by ~200 pcm in eigenvalue and over halved in the RMS error of pin power distribution.In order to increase the feasibility of the equivalent calculation scheme to practical applications, a number of colorset models are introduced and analyzed, including the single-pin, double-pin and assembly models.With DFs, the assembly colorset model slightly outperforms the most accurate homogenization approach without DFs while maintaining a superior computational efficiency.However, using DFs from various colorset models fails to achieve the prediction accuracy at the same level of the case when DFs are produced using the whole-core transport solution.
This study has also identified a few discrepancies in the implementation of the transport-corrected SP 3 methodology.For example, it has experienced numerical instabilities that cause an issue to achieve a converged solution when the DFs are imposed.Our preliminary investigation points to the use of the 1D forth order polynomial nodal expansion method implemented in the NEM solution method.We plan to implement and test a semi-analytical nodal expansion method in the future work to resolve this issue.
The inability to generate DFs for the water nodes in the reflectors is also found to be partially responsible for the deviations of the SP 3 results from the reference transport solution.Thus, the DF generation program will be revisited and updated to solve the fixed interface problem for the reflector region.At last, effort will continue to investigate colorset models that can maximize the benefit of the equivalent calculations while maintaining the computational cost to a reasonable level for practical whole-core simulations.

1 .
Perform the transport calculations and generate quantities including  , cell homogenized cross sections, cell averaged scalar flux, and side-dependent surface fluxes and currents.2. Solve the fixed interface problem for each cell sequentially by taking the reference values obtained in step one, generate the SP3 surface fluxes and currents, then compute the DFs according to their definitions.3. Execute a normal SP3 calculation on a pin-by-pin level using the homogenized cross sections from step one and DFs from step two to obtain the SP3 solution.

Figure 1 .
Figure 1.Procedure of equivalent core calculation using SP3 solver.

Figure 1 .
Figure 1.Procedure of equivalent core calculation using SP 3 solver.

Figure 2 .
Figure 2. OpenMOC model for C5G7 core.Materials are distinguished by different colors.

Figure 2 .
Figure 2. OpenMOC model for C5G7 core.Materials are distinguished by different colors.

Figure 3 .
Figure 3. Differences in relative pin power using the base and transport-corrected SP3 solvers.Blu regions indicate improved predictions by the transport-corrected SP3 solver, while red regions mean the opposite.

Figure 3 .
Figure 3. Differences in relative pin power using the base and transport-corrected SP 3 solvers.Blue regions indicate improved predictions by the transport-corrected SP 3 solver, while red regions mean the opposite.

Figure 4 .
Figure 4. Relative differences in pin power distribution predicted by a transport-corrected SP3 solver from the reference transport solution.The cross sections and DFs are generated using the assembly model (option 3).

Figure 4 .
Figure 4. Relative differences in pin power distribution predicted by a transport-corrected SP 3 solver from the reference transport solution.The cross sections and DFs are generated using the assembly model (option 3).

Table 1 .
Performance of the transport-corrected NEM SP 3 solver on various benchmark problems.The deviation from the reference results is given in pcm for the eigenvalue and RMS error for the pin power distribution.

Table 2 .
Performance of the transport-corrected SP 3 solver using different DF generation models.The deviation from the reference results is given in pcm for the eigenvalue and RMS error for the pin power distribution.