Extension of Pin-Based Point-Wise Energy Slowing-Down Method for VHTR Fuel with Double Heterogeneity

For the resonance treatment of a very high temperature reactors (VHTR) fuel with the double heterogeneity, an extension of the pin-based pointwise energy slowing-down method (PSM) was developed and implemented into DeCART. The proposed method, PSM-double heterogeneity (DH), has an improved spherical unit cell model with an explicit tri-structural isotropic (TRISO) model, a matrix layer, and a moderator for reflecting the moderation effect. The moderator volume was analytically derived using the relation of the Dancoff factor and the mean chord length. In the first step, the pointwise homogenized cross-sections for the compact was obtained after solving the slowing down equation for the spherical unit cell. Then, the shielded cross-section for the homogenized fuel compact was generated using the original PSM. The verification calculations were performed for the fuel pins with various packing fractions, compact sizes, TRISO sizes, and fuel temperatures. Additionally, two fuel block problems with very different sizes were examined and the depletion calculation was carried out for investigating the accuracy of the proposed method. They revealed that the PSM-DH has a good performance in the VHTR problems.


Introduction
The Korea atomic energy research institute (KAERI) has been developing very high temperature reactors (VHTR) with various sizes for the purposes of hydrogen production application or electricity supply for remote locations. The cores consist of graphite moderator and the fuel compact which contains a graphite matrix and tri-structural isotropic (TRISO) particle fuels randomly dispersed in the matrix. It introduces a unique neutronic characteristic, the double heterogeneity (DH), which is the combination of two level of heterogeneities: one for the fuel compact and the moderator and the other for the TRISO and the matrix within the compact.
The DH effect causes considerable complexity for the resonance treatment in a lattice calculation step. The SCALE code system employed a two-step resonance treatment method [1] for the DH fuel. The method generates pointwise homogenized cross-sections for the fuel compact through pointwise slowing down calculation for a spherical unit cell model with a fuel kernel and a matrix layer. By using the homogenized cross-sections for the fuel compact, the DH problem can be transformed into a single heterogeneity problem. However, this method not only ignores the effects of coating layers and moderator but also needs excessive computation time. The verification report [2] shows that the spherical unit cell without coating layers and a moderator causes error over 300 pcm. Afterward, a new resonance treatment method [3] for DH fuel was proposed based on the equivalence theory with the intermediate resonance (IR) approximation [4,5] and the Dancoff factor [5] for the particle fuel without coating layer. Kim [6] extended the method to take the coating layer of TRISO into account using the analytical Dancoff factor for TRISO particle fuel [7].
Another approach to enhance the accuracy of the equivalence theory for the DH fuel is to estimate accurate Dancoff factor using Monte Carlo method [8].
However, recent papers [9,10] report the limitation of the equivalence theory, overestimation of 238U absorption cross-section. The pin-based pointwise energy slowing-down method (PSM) [10] was developed to overcome the limitation of the equivalence theory. It can generate problem dependent self-shielded cross-sections in the lattice calculation step and has better accuracy than the pre-generated library within reasonably acceptable calculation time.
The lattice code developed at KAERI, DeCART (Deterministic Core Analysis based on Ray Tracing) [11], uses the pre-generated multi-group library. The library is prepared using the subgroup method [12] and the direct resonance integral table method [13] which are based on the equivalence theory with the IR approximation. The library generation procedure developed for conventional light water reactor (LWR) is used for the generation of the multi-group library for VHTR with a modification for DH fuel. During the library generation procedure for VHTR, the resonance integral table is corrected by using a Monte Carlo code. However, this correction leads to the dependency of the library on the specific reactor which the correction is performed based on. In addition, for improving the inaccuracy and system-dependency of the pre-generated library of DeCART, an extension of PSM for the resonance treatment of a DH fuel, PSM-DH, was developed and implemented into DeCART [14]. For extension of the PSM to a DH fuel, it is necessary to homogenize the fuel compact region. This approach is similar to the two-step resonance treatment method of SCALE but PSM-DH has an improved spherical unit cell model with the explicit TRISO fuel and a moderator to reflect the effect of the coating layers and moderator. Additionally, it applies the PSM in place of the equivalence theory for the homogenized fuel compact.
In this work, for completing our previous work [14], a new method for increasing the calculation speed of PSM-DH was proposed and the procedure for depletion calculation was established and the verification result was presented. Additionally, the DH effect of the proposed method was examined for typical VHTR problems and the verification calculations with various compact size and TRISO radius were performed.

Methods
The PSM-DH is essentially based on the PSM which can generate the self-shielded cross-section for a homogenized fuel region. The method needs a pointwise homogenized cross-section for a DH fuel of VHTR prior to the resonance treatment by PSM.
In order to provide a detailed description for the PSM-DH, a simple review of the PSM is firstly stated in Section 2.1 and then the extension of the PSM to a DH region is described in the next section.

Simple Review of Pin-Based Point-Wise Slowing down Method
The PSM solves the slowing down equation [10,15] for tens of thousands of energy points within the resonance energy range in lethargy space as follows: where i, j: sub-region indexes, u: lethargy, Σ t,i (u): total cross-section at sub-region i, φ i (u): the flux at sub-region i, P ji (u): probability that a neutron escaped from sub-region j has its first collision at sub-region i, V i : volume of a sub-region i, and Q s,j (u): slowing down scattering source of a sub-region j.
Note that the subscript M is the index of the moderator region. The collision probability (CP) tables for the given system should be known prior to solving Equation (1). The CP P ij can be evaluated by combining the CP table for an isolated fuel pin, P iso ij and the shadowing effect correction factor of the fuel region, η F . To evaluate P iso ij for every energy points, it is necessary to solve the fixed source neutron transport equation, Equation (2), for an isolated pin for every energy points. where and However, there are tens of thousands of energy points and it is time consuming to solve Equation (2) for every energy points. The computation time can be saved if the number of solution calculation of Equation (2) is reduced. Because the only input parameter of Equation (2) is the total cross-section (XS), P iso ij is a function of fuel region total XS, Σ t,F and it can be tabulated by solving Equation (2) for hundreds of Σ t,F values using method of characteristics (MOC). The CP table for an isolated fuel pin can be tabulated as a function of Σ t,F as follow: where φ j Σ k t,F is the flux at the sub-region j with the k-th level total cross-section of the fuel region, Σ k t,F , in the isolated fuel pin problem. The shadowing effect correction factor of the fuel region η F is defined as the ratio between the fuel escape probability in a lattice problem and the fuel escape probability in an isolated pin problem as shown in Equation (4).
where P e,F (u): fuel escape probability in a lattice problem and P iso e,F (u): fuel escape probability in an isolated pin problem. The fuel escape probabilities for the two systems can be calculated from the Carlvik's two-term rational approximation [16] with Dancoff factor. The Dancoff factor is obtained by solving the fixed source problem for the lattice problem using MOC, applying the enhanced neutron current method [17].
With the assumption that the ratio of P ij and P iso ij can be approximated as shown in Equation (5), P ij is expressed in terms of the η F and P iso ij as shown in Equation (6). where the following relationships are used: The CP table for an isolated fuel pin, P iso ij (u), is evaluated as follow using the tabulated function as prepared in Equation (3): Finally, the flux at an energy point of the slowing down equation, Equation (1), can be calculated with the CPs. Additionally, the effective multi-group cross-section for a resonance energy region can be simply obtained by the energy condensation as follow: where ∆u g is the lethargy width for an energy group g. Detailed derivation and further information on the PSM can be found in the reference [10].

Pointwise Energy Slowing down Method for a DH Region
For a lattice problem with double heterogeneous compacts, it is necessary to obtain the pointwise homogenized cross-sections for the compact prior to the resonance treatment by the PSM. Figure 1 shows the new spherical unit cell model proposed in this work. The spherical unit cell consists of an explicit TRISO, a matrix layer, and a moderator in the outermost layer. The moderator layer is to reflect the effect of the moderator in the pin cell. It is noted that one uniform XS region including the DH region in MOC calculation has one corresponding spherical unit cell. Thus, the homogenized compact can consider the effect of the radial particle distribution. The radii of the TRISO layers are given in the problem and the radius of the graphite matrix layer can be easily determined as follows: where F is the packing fraction of the TRISO particles in the compact and T R is the radius of the TRISO. For determining the radius of the moderator layer in the unit cell, the volume of the moderator for a specific pin cell of the lattice problem can be analytically derived from the Dancoff factor and the average chord length of the moderator region. Generally, they have the following relation [4]: The radii of the TRISO layers are given in the problem and the radius of the graphite matrix layer can be easily determined as follows: where F is the packing fraction of the TRISO particles in the compact and R T is the radius of the TRISO. For determining the radius of the moderator layer in the unit cell, the volume of the moderator for a specific pin cell of the lattice problem can be analytically derived from the Dancoff factor and the average chord length of the moderator region. Generally, they have the following relation [4]: where D is Dancoff factor and l is the mean chord length of the moderator region.l is the minimum distance between two compacts and Σ t,M is the total cross-section in the moderator. The Dancoff factor is calculated during the process of the PSM as described in Section 2.1. After replacing the mean chord length l in Equation (10) with 2π(R p +R c )H , the moderator volume for a pin cell of the lattice problem can be determined as follows: where R c , R p , and H are the radius of the compact, effective radius of the pin, and the height of the compact, respectively. Assuming that the moderator equally affects all TRISO particles in the compact, the outer radius of the moderator layer can be derived as follows: where N T is the number of the TRISO particles in a fuel compact. The pointwise energy slowing down equation for the spherical unit cell is as follow: where l, k: layer indices, L: total number of layers in a unit cell, and Q s,k : scattering source in layer k. The collision probability for the layer of the unit cell, P k,l , can be straightforwardly calculated by applying Kavenoky technique [18] which is widely used for CPs of a spherical geometry. Note that the CPs should be calculated with the white boundary condition for the unit cell. Equation (16) can be rewritten using the reciprocity relation as follow: After obtaining the pointwise energy flux for all layers from Equation (17), the pointwise homogenized cross-sections for the compact can be generated. The pointwise homogenized microscopic cross-section for a nuclide in the compact region can be obtained by the spatial homogenization as follows: where c represents the compact region which consists of the TRISO and matrix layer and N n k is the number density of nuclide n in layer k.
In the second step, the pointwise energy slowing down equation, Equation (1), for the homogenized fuel compact and moderator region can be readily solved by the original PSM.

Improvement of Calculation Speed for Collision Probability of Spherical Unit Cell
In Equation (16), P k,l should be obtained for all fuel compacts and all energy points at every burnup points in depletion calculation, which is a heavy computational burden. However, the computation time for the CP can be considerably cut down by applying the idea from the table of a CP used in the PSM.
In generating the CP table for the spherical unit cell by Kavenoky technique, it is assumed that the total cross-sections of all the layers except the fuel layer are constant in the resonance energy range and the CPs depend only on the total cross-section of the fuel kernel. Based on this assumption, the CPs are tabulated as a function of the total cross-section of the fuel kernel. This simplification can be justified by the comparison between the CPs of the original model and the simple model with the assumption. Figure 2 compares the CPs from the two models. The reference values are the CPs from the original model and were calculated with the total cross-sections for each their layer at the resonance energy region. On the contrary, the CPs from the simple model were calculated with the total cross-section of the fuel kernel at the energy point and the constant cross-sections for the other layers. The constant values for the layers were simply determined at the middle point of the resonance energy range. Except the very small oscillation in the low energy region, they are in a very good agreement. Actually, it was confirmed that the multiplication factors for pin problems obtained from two models are nearly the same. In this work, the CPs from a layer to a layer for the spherical unit cell were calculated using Kavenoky technique for wide range of total cross-sections with 500 levels from 10 −2 to 10 4 cm −1 and were tabulated as a function of the total cross-section of the fuel kernel. In a practical problem, the CP table for the specific spherical unit cell is prepared only once prior to the resonance treatment, regardless of a pin position inside a fuel block. Therefore, In this work, the CPs from a layer to a layer for the spherical unit cell were calculated using Kavenoky technique for wide range of total cross-sections with 500 levels from 10 −2 to 10 4 cm −1 and were tabulated as a function of the total cross-section of the fuel kernel. In a practical problem, the CP table for the specific spherical unit cell is prepared only once prior to the resonance treatment, regardless of a pin position inside a fuel block. Therefore, this technique can significantly reduce the computation time by determining the CPs for all spherical unit cells in the whole problem domain using the linear interpolation from the CP table as below in place of Kavenoky technique: where Σ t,F (u) is the total cross-section of the fuel layer at an energy point.

Reconstruction of Nuclide Number Densities for Depletion Calculation in a DH Region
In order to perform depletion analysis using the PSM-DH, a nuclide reconstruction process is needed for the pointwise unit cell slowing down calculation at the next burnup step. After a transport calculation using MOC, a depletion calculation for nuclide number densities at the next burnup step is performed for the homogenized compact region. The updated nuclide number densities for the homogenized compact can be simply converted to those for the fuel kernel using the reverse process of the homogenization. Afterward, the pointwise slowing down calculation for the double heterogeneous compact at the new burnup step can be performed using PSM-DH. Figure 3 shows the calculation procedure of the DeCART with PSM-DH. First, the collision probabilities for all layers of the spherical unit cell are calculated from the CP table which is priory generated using Kavenoky technique. After solving the pointwise slowing down equation for the spherical unit cell, the pointwise flux for all layers can be obtained and the pointwise homogenized cross-sections for the compact can be generated. The self-shielded multi-group cross-sections for the homogenized compact region can be then determined by applying the PSM. Finally, DeCART uses the cross-sections for a transport calculation in a DH region. In a depletion calculation, nuclide number densities of the fuel kernel of the spherical unit cell at the new burnup step should be reconstructed from the homogenized compact.

Reconstruction of Nuclide Number Densities for Depletion Calculation in a DH Region
In order to perform depletion analysis using the PSM-DH, a nuclide reconstruction process is needed for the pointwise unit cell slowing down calculation at the next burnup step. After a transport calculation using MOC, a depletion calculation for nuclide number densities at the next burnup step is performed for the homogenized compact region. The updated nuclide number densities for the homogenized compact can be simply converted to those for the fuel kernel using the reverse process of the homogenization. Afterward, the pointwise slowing down calculation for the double heterogeneous compact at the new burnup step can be performed using PSM-DH. Figure 3 shows the calculation procedure of the DeCART with PSM-DH. First, the collision probabilities for all layers of the spherical unit cell are calculated from the CP table which is priory generated using Kavenoky technique. After solving the pointwise slowing down equation for the spherical unit cell, the pointwise flux for all layers can be obtained and the pointwise homogenized cross-sections for the compact can be generated. The self-shielded multi-group cross-sections for the homogenized compact region can be then determined by applying the PSM. Finally, DeCART uses the cross-sections for a transport calculation in a DH region. In a depletion calculation, nuclide number densities of the fuel kernel of the spherical unit cell at the new burnup step should be reconstructed from the homogenized compact.

Numerical Results
For verifying the performance of the PSM-DH, the resonance treatment and transport calculations for various HTGR problems with DH regions based on MHTGR-350 [19] were performed and the results were compared with the reference values by McCARD [20] with ENDF/B-VII.1 library. DeCART used the multi-group cross-section based on

Numerical Results
For verifying the performance of the PSM-DH, the resonance treatment and transport calculations for various HTGR problems with DH regions based on MHTGR-350 [19] were performed and the results were compared with the reference values by McCARD [20] with ENDF/B-VII.1 library. DeCART used the multi-group cross-section based on ENDF/B-VII.1 for the fast and thermal energy range, which was processed by GROUPR of NJOY [21]. For the resonance energy groups, the PSM-DH module integrated in DeCART directly reads point-wise cross-sections processed by BROADR of NJOY and generates the self-shielded cross-sections. DeCART used the ray space of 0.02 cm, the number of polar angles of 4, and the number of azimuthal angles of 8 for the MOC calculation option. Figure 4 shows the configuration of the DH fuel pin used in the verification calculation and Table 1 lists nuclide number densities of the materials. The pin problem is taken from the OECD/NEA MHTGR-350 benchmark exercise III [19].  Figure 4 shows the configuration of the DH fuel pin used in the verification calculation and Table 1 lists nuclide number densities of the materials. The pin problem is taken from the OECD/NEA MHTGR-350 benchmark exercise III [19].   Figure 4. Configuration of the MHTGR fuel pin with DH region. Tables 2-5 compare the multiplication factors and the DH effect by McCARD and DeCART with PSM-DH for various packing fractions (PF) at 300 K, 600 K, 900 K, and 1100 K, respectively. The DH effect is defined with the multiplication factors of the homogeneous and double heterogeneous compact as follow:          The differences of the k inf in all cases are within 203 pcm and those of the double heterogeneous effect are within 146 pcm. McCARD uses the JT algorithm [22] for the random distribution of the TRISO particles and the random distribution of the TRISO particle causes additional statistical uncertainty in multiplication factor over 100 pcm. Additionally, it is known that the PSM has error of about 100 pcm [10]. Besides that, there is another error cause. The PSM-DH can consider the DH effect only for the resonance energy range. However, it was found that the error caused by the spatial self-shielding effect of DH in the thermal energy range and fast energy range is about 30 pcm in the problem.

MHTGR-350 Benchmark Problems
Hence, considering the statistical uncertainty caused by the random distribution in the McCARD and the inherent error of the PSM, it is clear that they are in good agreement. Figure 5 compares the multi-group absorption cross-sections in the resonance energy range for the homogenized compact with a packing fraction of 35% at 1100 K. There are differences of the maximum 4% between those by the two codes. Figure 6 shows the disadvantage factors of the fuel kernel and the matrix calculated by PSM-DH method.
dom distribution of the TRISO particles and the random distribution of the TRISO particle causes additional statistical uncertainty in multiplication factor over 100 pcm. Additionally, it is known that the PSM has error of about 100 pcm [10]. Besides that, there is another error cause. The PSM-DH can consider the DH effect only for the resonance energy range. However, it was found that the error caused by the spatial self-shielding effect of DH in the thermal energy range and fast energy range is about 30 pcm in the problem. Hence, considering the statistical uncertainty caused by the random distribution in the McCARD and the inherent error of the PSM, it is clear that they are in good agreement. Figure 5 compares the multi-group absorption cross-sections in the resonance energy range for the homogenized compact with a packing fraction of 35% at 1100K. There are differences of the maximum 4% between those by the two codes. Figure 6 shows the disadvantage factors of the fuel kernel and the matrix calculated by PSM-DH method.
Additionally, for verifying the accuracy of the proposed method for the moderation effect inside a fuel block, two fuel block problems with different size were examined. They are a normal size fuel block with 210 fuel pins and a small size fuel block with 12 fuel pins which have the same configuration with the pin problem as shown in Figure 4 and Table  1. Figure 7 shows the configuration of the two problems.   Additionally, for verifying the accuracy of the proposed method for the moderation effect inside a fuel block, two fuel block problems with different size were examined. They are a normal size fuel block with 210 fuel pins and a small size fuel block with 12 fuel pins which have the same configuration with the pin problem as shown in Figure 4 and Table 1. Figure 7 shows the configuration of the two problems. Tables 6 and 7 compare the results for the problems. They reveal that the differences of all cases are under 100 pcm. In the small size block, it is noted that there is a rapid change of the moderation effect between the outer fuel pin and the inner fuel pin due to the outmost moderator region. Nevertheless, the code produced the accurate results in the small size block.   Table 8 shows the burnup calculation results for the MHTGR pin problem. It was calculated under the condition of the hot full power for the DH pin with the packing fraction of 35% and the burnup at the 1500 effective full power day is about 98 GWD/MTU. The result shows that the differences between the two codes are within ±120 pcm. It can be seen that the PSM-DH/DeCART also has a good performance in the depletion calculation.  Table 9 shows the results for the pin problems with very large size compact. The half pitch size and the PF of the pin are fixed by 2.0 cm and 35%, respectively. The material compositions are the same with those of the previous problems. Considering the compact radius of the MHTGR-350 pin, 0.6225 cm, it is noted that they can cover large compact problem. The results reveal that the discrepancy between two codes are within 250 pcm and they are similar to the previous comparisons. Table 9. k inf of pin problems with large size compact and pitch.   Table 10 shows the results for the pin problems with various fuel kernel and TRISO radius. The compact radius was fixed by 0.6225 cm. In all cases, the results by two codes are similar within 150 pcm. It is clear that they are in a good agreement considering the statistical uncertainty by the random distribution of the particles. Table 10. k inf of pin problems with various fuel kernel size.   Table 11 compares the calculation time in DeCART with the pre-generated library and the PSM-DH for MHTGR single block problem as shown in Figure 5a. Actually, the existing DH treatment method in DeCART is the Sanchez-Pomraning method [23]. It needs much more calculation time than a homogenous case, because it performs the DH treatment for the explicit compact geometry with TRISO and matrix. However, DeCART with PSM-DH treats DH region as a homogeneous region and the calculation time was cut down to 72% compared with that by the existing DH treatment module for the DH problem.

Conclusions
In this paper, the extension of the PSM was proposed for performing the resonance treatment of a DH fuel of VHTR and the verification results for various VHTR fuel pins and fuel blocks were presented.
The pointwise homogenized cross-sections for the compact was obtained after solving the slowing down equation for a spherical unit cell composed of an explicit TRISO, a matrix layer, and a moderator layer. The moderator volume was analytically derived using the relation of the Dancoff factor and the mean chord length. Then, the shielded cross-section for the homogenized fuel compact was obtained using the original PSM.
The verification calculations were performed for the fuel pins with various packing fractions, compact sizes, TRISO sizes, and fuel temperatures. Additionally, two fuel block problems with very different sizes were examined and the depletion calculation was carried out for investigating the accuracy of the proposed method. They revealed that the PSM-DH has a good performance.
The PSM-DH has an advantage that problem dependent self-shielded cross-sections can be generated in the lattice calculation step for a VHTR problem. Thus, it is expected that the method could be applied for the development of the various type VHTRs. In the near future, it is planned that an improved model for considering multiple TRISO types in a compact will be examined.