Improvement to NEM SP3 Modelling and Simulation

Accurate reactor core steady state safety analysis requires coupling between thermalhydraulics and three dimensional multigroup pin by pin neutronics. Concerning the neutronics modeling, the Nodal Expansion Method (NEM) code is developed at North Carolina State University in the framework of high fidelity multiphysics coupling with CTF. NEM includes a simplified third-order Spherical Harmonic (SP3) solver. In this work, the solver has been improved by incorporating higher order scattering matrix library. The boundary conditions were corrected with one dimensional P3 theory and a consistent coupling coupling between zerothand second-order flux moments was established. two methods for generating second order discontinuity factors (DFs) has ben developed, one based on the Generalized Equivalence Theory (GET) and one based on Parial Current Equivalence Theory (PCET). DFs were generated with three lattice sizes: single pin, 2 pins and assembly level. These developments were tested using the C5G7 benchmark. The results of the SP3 solver improvement, by using P2 and P3 scattering cross sections, show a 50% decrease in the eigenvalue (keff) prediction error compared to the reference transport solution. The GET DFs are applied in the C5G7 core pin by pin calculation and are compared with PCET DFs. The results show that PCET have a better performance in global results (eigenvalue). Concerning the different lattice sizes studies, the results show that DFs generated in smalll colorsets can improve local solutions. However, in order to reveal strong global trends, DFs should be generated in a larger corloset representative of the whole core. For the core calculations, DFs generated with the three colorsets together with an additional mixed type DFs were tested. For the mixed type, DFs generated from assembly size lattice were used for the internal interfaces and DFs generated from 2 pins size lattice were used for the assemblies boundary interfaces. These mixed DFs outperformed all the other configurations indicating that they manage to accomplish a satisfying compromise between global and local trends.


INTRODUCTION
Accurate reactor core steady state safety analysis requires coupling between thermal-hydraulics and three dimensional multigroup pin by pin neutronics [1,2]. Monte Carlo Methods are considered for high fidelity coupling safety analysis [3,4] but compared to lower order determinsitic (Diffusion) methods they are computationally expensive. Simplified P3 theory is believed to be more compatible than diffusion theory in advanced reactors, mini cores, and other complicated whole core three dimensional transport calculations at pin level resolutions [5]. At locations in the core with high flux gradient, such as material boundaries, the SP3 has high precision due to the use of the higher order Pn scattering library. Diffusion, on the other hand, can only describe linear scattering, something that limits its applicability in modern reactor designs.
The NEM (Nodal Expansion Method) code has been developed for three-dimensional whole-core steady-state and transient analysis with cartesian, hexagonal, and cylindrical geometries. It is used in the high fidelity multiphysics CTF-NEM coupling framework. NEM has a diffusion based neutronic solver. This solver was recently extended to solve SP3 equations in order to achieve a pin cell resolution of highly heterogeneous reactor cores [6]. However, the SP3 solver does not outperform the diffusion solver by a large margin, in terms of the solution accuracy, without using discontinuity factors (DFs). The current method for generating DFs is based on the assembly lattice calculations. However, the lattice size could be reduced to few pins model and still preserve the accuracy of the solver. This paper summarizes some recent developments for the NEM SP3 solver enhancement. In Section 2 we briefly discuss the SP3 equations of the solver. The advantages of such higher order solver requires corresponding scattering cross sections. Hence, in Section 3 we present the the development of an extension in NEM to take into account higher moments of scattering cross sections and corrections regarding the boundary conditions. This creates an improved solver modeling (NEM Version 2) that we test and verify in the C5GC7 benchmark core. The eigenvalue (k eff ) of the core is calculated and compared with the previous SP3 version of the solver (NEM Version 1) and the Diffusion solver. The correct calculation of local quantities of a reactor core such as the flux and power in the different pin cells requires a consistent calculation of Discontinuities Factors (DFs) in order to account for the various heterogeneities. Reproducing whole core transport solutions with DFs generated from corresponding transport calculations is impractical. To this purpose a new DFs method was implemented in NEM and described in Section 4.1. The method generates DFs based on both Partial Current Equivalence Theory (PCET) and Generalized Equivalence Theory (GET) to eliminate the error introduced by pin-wise homogenization. In this approach, side-dependent zeroth and second order DFs are generated for each type of pin cell neighbored by other types of pin cells within different lattice size (colorset) calculations. There are three types of lattice size being studied. In the single pin model, lattice calculations are performed in 1×1 pin cell geometries with reflective boundary conditions. In 2 pins model, lattice calculations are performed in 1×2 pin cells geometries with reflective boundary conditions. DFs are generated at 2 pins' interface and applied between corresponding different pin cells in the whole core calculations. In the assembly model, lattice calculations are performed in single assembly geometries with reflective boundary conditions. In section 4.2, the generated DFs are verified with comparison to the reference transport solution in both assembly and core calculations. In Section 4.3, GET DFs and PCET DFs generated from different colorsets are compared in the whole C5G7 benchmark core. Finally, conclusions and future prospects are provided in Section 5.
In the NEM solver, the same equations are used in each direction and they are coupled through transverse leakage terms illustrated for the X direction by the following equations: with The Δi is the node cell length in the corresponding directions.

Modification To NEM SP3 Solver
Simplified P3 theory is believed to be more compatible than diffusion theory in advanced reactors, mini cores, and other highly heterogeneous whole core three dimensional transport calculations. Compared to practical diffusion application, SP3 is able to make most of higher order Pn scattering matrix library up to P 3 , which shows significant precision even at high flux gradient locations such as materials boundaries. The Version 2 NEM SP3 solver was developed in this work to include the high order cross sections using the outgoing approximation presented in Eq. (6).
The tightly coupling between Φ 0 and φ 2 is essential to demonstrate the improvements from higher order cross sections. All flux terms, for instance Φ 0 and φ 2 , in Eq.(1) are solved simultaneously through source iterations. The Marshak boundary condition from one dimension P3 theory replaces the original one used in NEM SP3 solver. Between neighboring nodes the correct boundary conditions of Eq.( 7) are implemented.

Verification of Improved NEM SP3 Solver
The steady-state two-dimension mini core C5G7 benchmark [7] was used to verify the Version 2 of SP3 solver. This benchmark is widely used for code verification purpose and is the one used to test the whole CTF-NEM multiphysics coupling framework. The core consist of four different fuel assemblies as illustrated in Fig.(1). The original C5G7 benchmark provides a scattering matrix only up to P1. For this work, where we want to test the higher order scattering up to P3 this is not enough. To this purpose, a modified version of the benchmark is considered, where Serpent2 [8] was used to generate the scattering cross sections. The pin-wise homogeneous cross sections feeding NEM were calculated by Dragon5 [9], using the method of characteristic solver with anisotropic source.  The steady state calculation of the C5G7 was performed with the different versions of the NEM solver (Diffusion, Version 1 SP3 and Version 2 SP3). The comparative results concerning the k eff are gathered in Table.1. Version 2 SP3 solver solutions demonstrate increasing accuracy at higher order cross sections. The same cross sections are used by the Version 1 SP3 solver, however no benefits is observed. The implementation of higher order cross sections in Version 1 and Version 2 solver are the same, but the effect of them are neglected because the higher moment fluxes are not fully coupled with low moment fluxes. Without at least P 1 cross sections, the SP3 method does not show any advantage compared to the diffusion method.

IMPROVEMENTS TO SP3 PRACTICAL SIMULATION WITH DFS
The SP3 solver does not outperform the diffusion solver by a large margin in terms of the solution accuracy without discontinuity factor technology. High order DFs were recently developed in the NEM SP3 solver to preserve node interface scalar flux and net currents. The high order DFs were applied to preserve cell interface flux φ 2 and showed expected improvements in whole core calculations [10]. The second order flux φ 2 can be very small and introduce numerical errors that will be amplified depending on the DFs definitions. To avoid this problem, a new DFs function was implemented in NEM. Different types of DFs were implemented in SP3 solver with different characteristics. Generalized Equivalence Theory (GET) DFs, defined as Eq. (8), preserve interface scalar flux and net currents.
Partial Current Equivalence Theory (PCET) DFs as Eq.(9) preserve interface partial currents and net currents.
The performance of the two types of DFs in different colorsets is studied as well and discussed in the following subsection. Three different types of colorsets are considered.

Generation of DFs in Colorsets
The previous DFs generation method is based on assembly lattice calculations, in which side dependent DF 0 and DF 2 are generated for each pin cell. For practical purposes, few pins colorsets could be used to generate DFs with reasonable precision. In this section single pin model, 2 pins model, and assembly model are proposed. The different studied DFs are tested using the C5G7 benchmark [10] as in section 3.2, in which there are: one type of UO 2 pin, three types of MOX pins, one type of guide tube, and one type of fission chamber. The three MOX fuels pins are 4.3%, 7.0%, and 8.7% Plutonium weight enriched. In this section, both single assembly problems and the whole problem are tested.
In this section, all the transport calculations are performed using the OpenMOC code [11]. The different materials cross sections are computed by Serpent as in Section 3.2. Because OpenMOC does not do higher order scattering calculation, only P 0 isotropic cross sections are used. Transport corrections to the cross sections are not performed, as it will result in negative in-group scattering cross sections.
For this work a dedicated program was developed to process the output of OpenMOC code. Different pin cell information are extracted: interfaces angular fluxes, cell averaged scalar flux and eigenvalue. Using these data the interface DFs are computed for each pin cell. Although Dr Chao [12] proposed two new methods describing the boundary conditions in a better way, to keep consistence of boundary conditions between reference solutions and NEM, one dimension P3 theory is used in this method. Verification of GET DFs is made by equivalence calculations based on assembly and whole core colorsets with the following process: 1. Run reference colorsets calculations with OpenMOC and generate k eff , cell-homogenized cross-sections and side-depndent heterogeneous surface fluxes/currents. 2. For each cell apply data processing program using the information including k eff from step 1 to generate homogeneous surface fluxes and currents. Compute the DFs using these fluxes and currents. 3. Run pin-by-pin homogeneous calculations using the homogenized cross-sections from 1 and DFs from 2. 4. Compare results from 3 to reference results from 1.
Performance of different DFs generated from different colorsets are tested in the C5G7 whole core problem. Reference solution is calculated by OpenMOC with options: 32 azimuthal angles, 3 polar angles, ray trace width of 0.03 cm and k eff criterion of 1E-7. Comparisons between GET and PCET DFs are performed together with comparisons between three different colorsets sizes. The three different colorset sizes are: 1. Single pin with all reflective boundary conditions. The unity DFs are calculated for each side of the different type of pins: UO 2 , MOX and Fission champers. Non-fissionable guide tubes and water cells are surrounded by UO 2 fuel pins in 3 × 3 colorsets. 2. Eight sets of 2 × 1 pin colorsets with reflective boundary conditions were defined to perturb the partial current and generate neighbor corrected interface DFs. These DFs will be used in the whole core on the interface between neighboring different cells. 3. Single UO 2 and MOX assembly transport calculations with reflective boundary conditions. DFs are generated for the four different interfaces.
The two different DFs types and the three different colorset creates a total of 6 different DFs configurations. Additionally, a mixed type of DFs is considered. For this type, DFs generated from the assembly colorset are used for the assemblies internal interfaces while DFs generated from 2 pins colorset are used for the boundary interfaces between assemblies. This creates a total of 8 different DFs configurations.

Equivalence Calculations
Before analyzing the performance of DFs, verification is required for the whole procedure including reference calculation, generating DFs and implementing DFs in NEM. Cell interface GET DFs are generated for each pin cell of the UO 2 assembly, MOX assembly and the quarter core of C5G7 benchmark. For all the calculations reflective boundary conditions are used. GET DFs reproduce the OpenMOC solutions including k eff , power distribution and flux distribution. The equivalence calculation results regarding the eigenvalue (k eff ) are presented in Table.2 and show a 0 pcm discrepancy for both assemblies and the core. The power distribution error in the whole core is presented in Fig.2. Power distribution and total flux distribution have been normalized to average value 1. The maximum power and flux error are 0.13% and 0.15% respectively and occur in the MOX assembly. The comparative results in both assemblies and the core verify that GET DFs can reproduce the OpenMOC reference solution with very high precision.

Comparison of DFs and colorsets in C5G7
C5G7 benchmark is used to test the GET and PCET DFs generated in single pin and two pins colorsets as described in Section 4.1. The DFs in single pin cell will be used in the whole core for each pin cell. Two pins based DFs will be applied at the interface where one of eight types of interface occurs. No DFs is applied between same pin cells. All NEM calculations with the 8 different DFs configuration described in the last section are performed with convergence criterion 1E-6 for k eff . In Table 3 we present the results using only the P0 cross sections while in Table 4 we present the results using the transport corrected P0 cross sections. The results of Table 3 show that PCET DFs offer a small improvement on the eigenvalues and power distribution compared to GET. Neighbor corrected DFs in larger colorsets have better precision at power distributions and peak power for both GET and PCET DFs.
By interpolating the DFs generated in whole core problem in Section 4.2, most of DFs increase leakage rate from center region to outside boundary. This is because the P 0 cross section have a smaller diffusion coefficients D 0 and D 2 , hence DFs play the key role of global transportation behaviors. PCET DFs could offer a better solution because they preserve partial currents that have a global effect. However, both of PCET and GET DFs are generated in colorsets with reflective boundary conditions and thus it is difficult to allow DFs to represent global trends.
The improvements to local solutions by DFs is demonstrated by applying transport corrected P 0 cross sections. They have relative larger diffusion coefficients which could match better with reference solutions in the sense of global leakage rate and therefore decrease the global trend information carried by the DFs. In the results of Table 4 for the corrected transport cross sections we see that the performance of GET and PCET are very similar. Single pin DFs without neighbor correction worsen the power distribution with both P 0 and transport corrected P 0 cross sections. In transport corrected cross sections calculations , DFs generated using larger colorsets increase the precision of power distribution. Finally, the assembly mixed with 2 pins model based DFs show improvements in the both eigenvalues and power distribution. This means that the mixed type DFs manage to find a good compromise between local and global trends.   show a 50% decrease in the eigenvalue error compared to the previous versions of the NEM code. This significant improvement is related to the higher than P1 cross sections and to some boundary conditions corrections.
A new high order DFs implementation method is developed in NEM for SP3 solver. The whole procedure of generating and utilizing high order DFs for SP3 solver has been verified in one assembly, and the whole core problem. Promising results for k eff , power distribution, and flux distribution are obtained. Reference transport solutions were reproduced by applying DFs.
Few pins colorsets are developed to correct pin cell interface discontinuity factors. GET and PCET theories are considered in this simplified method. PCET DFs show better performance in global results eigenvalue, with P0 cross sections. With transport corrected P0 cross sections GET and PCET DFs show similar behavior. DFs generated in smalll colorsets can offer improvements to local solutions. To reveal strong global trends, DFs should be generated in a larger corloset. The DFs were studied with three different colorsets: 1 pin , 2 pins and single assembly. An additional mixted DFs type was considered: the DFs from the assembly colorset were used in the internal interface while the DFs from the pins colorset were used for the assembly boundary interfaces. The mixed DFs type outpeform all other configurations in terms of both eigenvalue and power distribution. This result indicates that the mixed DFs manage to find a good compromise betweem local and global trends. Dr. Chao proposed rigorous boundary conditions of SP3 theory. In the future work, rigorous SP3 theory will be implemented in the NEM as one additional solver.
Corresponding rigorous DFs will be generated for this new solver.