Implementation of 2D/1D Gamma Transport and Gamma Heating Capability in MPACT

This paper presents the most recent progress on the development of gamma transport capability for the CASL neutronics code MPACT: (1) 3D gamma transport and (2) explicit gamma heating capabilities. The 3D gamma calculation capability was implemented by leveraging the 2D/1D solver originally developed for neutron calculations. The results were verified by MCNP6 on a small assembly with 5 × 5 pins. Generally, errors were lower than 0.5% on each axial mesh as long as MPACT was running with enough axial meshes. The gamma heating calculation considered the energy deposition from photoelectric absorption, Compton scattering, and pair production. Verification with MCNP6 for both 2D and 3D benchmarks showed that the errors of energy depositions are comparable with those of gamma fluxes, proving the proper implementation of the energy deposition.


INTRODUCTION
The CASL (Consortium for Advanced Simulation of Light Water Reactors) developed the code MPACT, which is a three-dimensional (3D) full-core neutron transport code that solves the Boltzmann transport equation for the neutron flux distribution. MPACT commonly uses the 2D/1D method that treats the radial (x and y) variables differently from the axial (z) variable. The radial problem is treated by transport theory and is calculated by the 2D method of characteristics (MOC) on a fine spatial and angular mesh within each axial slice, and the axial problem is treated by NEM-P3 method [1]. Previous work successfully implemented the 2D gamma flux calculation capability [2], and this paper presents the recent progress on 3D gamma transport calculations using 2D/1D method in MPACT. In addition, explicit gamma energy deposition is developed and presented in this paper. The distribution of energy deposition in a reactor could be noticeably different from the fission rate distribution primarily due to the effect of global gamma transport [3]. The current energy deposition model in MPACT uses an empirical gamma smearing formula developed for PWR applications [4]. To improve the accuracy of heating calculation and extend MPACT capabilities for other reactor types, a general heating model based on gamma transport is desired.

Gamma Transport
Similar to neutron transport, gamma transport is focused on numerically solving gamma transport Eq. (1): where γ /N is gamma or neutron, and h/g is gamma or neutron groups.
ψ γ h , gamma angular flux in group h Σ tr,h , gamma transport cross section in group h Σ s,h →h , gamma transport-corrected scattering cross section from group h to h Σ N γ,g→h , neutron induced gamma production ((n,γ) source) cross section.
The equation has a form similar to the neutron transport equation, which makes it possible to reuse many of the solvers implemented for the neutron transport, including the 2D/1D solver with the radial 2D MOC sweeper and the axial NEM-P3 method, which will be discussed in the next subsection. Scattering and (n,γ) sources are assumed to be isotropic. A new procedure was developed for the neutron-induced gamma source calculation, and the gamma calculation was performed as a fixed-source problem.

2D/1D Method
The 2D/1D method for neutron calculations in MPACT is presented in detail by Collins et al. [1], and this subsection briefly summarizes some key points of the method. The 2D/1D method decomposes 3D problems into 2D planes in which the transport equation is solved. The planes are coupled together through a leakage source, and the axial variation in the flux is modeled using a lower-order transport approximation discretized by either a nodal or finite difference method.
The radial equation is obtained by averaging Eq. (1) axially in plane k: which includes an axial transverse leakage (TL) term in the source. For simplicity in calculation, the TL term is usually approximated to be isotropic and spatially flat within a pin cell in practice.
The axial equation is obtained by integrating Eq. (1) radially over both x and y: where the radial transverse leakage T L XY h (z, Ω) is usually assumed to be isotropic and can be constructed using the net currents at the right and left boundaries of each coarse radial cell.

Data library
As in previous work [2], the HELIOS library obtained from the University of Michigan was used for this development. The HELIOS library has 47 neutron groups and 18 gamma groups; the ngamma production data in the HELIOS library include all prompt, capture and delayed gamma rays, forming a big neutron-to-gamma group matrix(47 × 18).
Transport-corrected Compton scattering cross section integrates the contribution from pair production and forms the scattering kernel like the following: where c stands for Compton, pp stands for pair production, h is incoming gamma group, h is outgoing gamma group, and hPair represents the energy group containing the pair productioninduced photon (0.511 MeV). Gamma cross sections are temperature independent since all these interactions are between photons and the electrons surrounding the nucleus.
However, the HELIOS gamma scattering matrix does not include photon fluorescence reactions. The lack of this data can cause gamma flux to be underpredicted in the group where k-shell emission peak is located, e.g., about 0.1 MeV for uranium isotopes. For flux comparison, MCNP6 references were run with the photon fluorescence turned off to be consistent with the HELIOS library. This is discussed in the MCNP Setup subsection.
The transport cross section is defined as where H represents all gamma groups, PE is photoelectric, PP is pair production, and C is transport corrected Compton scattering.

Heating Calculation Scheme
The original scheme is presented by Luthi [5]. The total heat deposition H γ (r) is calculated as follows: where Φ(E, r) is gamma flux, and N i (r) is the number density of element i. K i,x (E, r) is the KERMA(Kinetic Energy Released to Material) factor for element i and reaction x at incident energy E and is defined as follows: where p.e. stands for photoelectric effect, c.s. for coherent scattering, i.s. for incoherent (Compton) scattering, and p.p. for pair production.
MPACT uses multigroup cross section libraries. HELIOS combined Compton scattering and pair production into a single matrix. Thus, Eqs. (6) and (7) are rewritten as follows for MPACT: Note that the H on the left side is heat, and the H on the right is total gamma group number H.

MCNP Setup
To be consistent with the available HELIOS library at University of Michigan, ENDF/B-VI data are used in most MCNP runs except those with isotopes that ENDF/B-VI in MCNP does not have gamma data like gadolinium. MCNP does not support direct delayed gamma tally. In this work it is scaled from prompt gamma as follows: where Q delayed and Q prompt are delayed and prompt energy per fission from the ENDF library respectively, and φ h,prompt (r) is prompt gamma flux in group h. Liu et al. [6] discusses the validity of this approximation.
As mentioned above, the HELIOS gamma scattering matrix does not include photon fluorescence reactions. MCNP does include photon fluorescence reactions, so MCNP references for flux distribution were run with the photon fluorescence effect turned off to match MPACT's spectra. However, MCNP references for energy depositions were run with this effect turned on. When this effect is turned off in MCNP, energy calculations are significantly impacted. For example, energy by prompt gamma per fission for 3.1% uranium fuel, which is supposed to be about 6.9 MeV according to ENDF.VI data and is actually 6.64 MeV in MCNP6 runs, becomes 6.17 MeV when the photon fluorescence effect is turned off.

2D Assembly
The first case is based on VERA Progression problem 2b [7]. Both the fuel and the moderator were set to 293K. Neutron and gamma flux distributions of this case are presented in Wang et al. [2], and this section focuses on the energy deposition distributions. As shown in Figure 1, depositions in guide tubes are about 15-20% of the values in the surrounding fuel pins. Gamma energy deposition distribution calculated by MPACT is consistent with that calculated by MCNP. Differences in all pins are under 1%. CASL problem 2h is chosen to test the capability to calculate control rod cases. Neutron and gamma flux distributions of this run are presented in Wang et al. [2]. As shown in Figure 2, gamma depositions in B4C control rods are about 25-30% of the values of fuel rod pins. This is higher than that in empty guide tubes from case 2b because B4C has stronger gamma sources compared to moderators. The errors are < 1% in most pins compared to MCNP. The maximum difference of 3% only occurs at the central guide tube pins, where the absolute deposition value is small.

3D Mini Assembly
To test the newly developed 3D gamma transport capability, a 5×5 mini assembly case was created. The geometry is shown in Figure 3. The assembly is 100 cm tall and has 20 fuel pins and 5 guide tubes. Everything is at 293K. It has vacuum boundaries on the top and bottom and reflective radial boundaries. The fuel pins are identical to those in CASL problem 2b and 2h, each with a 3.1% fuel rod with a radius of 0.4096 cm, cladding with an inner radius of 0.418 cm, and an outer radius of 0.475 cm with moderators outside. The guide tubes have an inner radius of 0.561 cm and an outer radius of 0.602 cm. A B4C control rod is inserted in the center guide tube. Results of the mini assembly case in which the control rod is half inserted are shown in Figure  4. Neutron and gamma fluxes are radially integrated and presented as axial distributions. Neutron flux calculated by MPACT is consistent with that calculated by MCNP, with errors less than 1%. Such consistency verifies the effectiveness of this modeling. Gamma flux results in MPACT are noticeably different compared with results from MCNP. However, the accuracy can be improved by running the case with finer axial meshes and integrating the result to 10 axial levels. As shown in Figure 4b, if the calculations are performed with 40 axial layers and integrated into 10, then the differences with MCNP results are all within 1%. The gamma energy deposition has a similar pattern as the gamma flux distribution, as shown in Figure 5. The errors are 2-3% on average when running with 10 planes, and using finer meshes can improve the results as in the flux modeling. When running with 40 meshes and integrating into 10 to compare with MCNP, errors are less than 0.5% in 9 out of 10 meshes. One exception is the bottom of control rod, which shows a 2% difference. Overall, the 2D/1D calculation gives accurate results for gamma flux and energy deposition.

CONCLUSIONS
In this work, the 2D/1D method for 3-D neutron transport in MPACT was leveraged to solve the gamma transport, and the gamma energy deposition was implemented in MPACT. Several 2D/1D functions written for neutron transport, such as nodal axial solvers and coarse/fine mesh homogenizing/projecting functions, were reused in the gamma code due to their similar transport behaviors. Energy deposition by photoelectric absorption, Compton scattering and pair production were calculated explicitly in this work, and benchmark results show that the algorithm and implementation work well on both flux and energy deposition calculations.

FUTURE WORK
Thermal feedback: To further improve thermal feedback calculations in MPACT, we will add gamma heat depositions to the thermal feedback scheme next.
ORNL Gamma Library: As indicated in the results of this work, some errors are thought to be the result of inconsistencies in the HELIOS library compared to MCNP's ENDF6. Therefore, future work will include runnning the cases again with the more current ENDF7-based ORNL gamma library once it is available.