Abstract
Very high temperature reactor (VHTR) designs offer promising performance characteristics; they can provide sustainable energy, improved proliferation resistance, inherent safety, and high temperature heat supply. These designs also promise operation to high burnup and large margins to fuel failure with excellent fission product retention via the TRISO fuel design. The pebble bed reactor (PBR) is a design of gas cooled high temperature reactor, candidate for Generation IV of Nuclear Energy Systems. This paper describes the features of a detailed geometric computational model for PBR whole core analysis using the MCNPX code. The validation of the model was carried out using the HTR-10 benchmark. Results were compared with experimental data and calculations of other authors. In addition, sensitivity analysis of several parameters that could have influenced the results and the accuracy of model was made.
1. Introduction
Very high temperature reactor (VHTR) designs offer promising performance characteristics. If realized, these concepts can provide sustainable energy, offer improved proliferation resistance, and are more easily safeguarded than current light water reactors (LWR). These designs also promise operation to high burnup and large margins to fuel failure with excellent fission product retention via the TRISO (tristructural-isotropic) fuel design. The higher temperature of operation for these concepts can support industrial process applications that cannot be easily supported by LWR technology [1].
Two VHTR concepts have been studied: the prismatic reactor, with block-type fuel elements, and the pebble bed reactor (PBR), with spherical fuel elements, as shown in Figure 1.
In PBRs, the fuel is contained within graphite pebbles which form a randomly packed bed inside a graphite-walled cylindrical cavity. Due to the stochastic nature of this bed, the location of the individual pebbles is not well defined. The pebble bed in such a reactor for some type of calculations is often modeled as a homogeneous mixture of pebbles and coolants materials, with a uniform density throughout the core. Unfortunately, such a model does not include all the effects that the heterogeneity of the pebble bed entails, resulting in possible errors. Three of these effects are (1) the density fluctuations in the pebble bed near the wall, (2) the neutron streaming through the void space between the pebbles, and (3) the variations in the Dancoff factor near the edge of the pebble bed [2].
In the framework of a coordinated research program (CRP), the high temperature gas-cooled test reactor (HTR-10) benchmark was solved by several research organizations with the aim of validating analytical codes and performance models for actual operating conditions of HTRs. The calculations and experimental results were compared. The HTR-10, built in China, is an experimental facility whose construction was completed in 2000 and its initial criticality took place in December of that year. The HTR-10 should serve as a test facility to demonstrate the inherent safety features of the modular reactor design and to gain acceptance from regulatory institutions, utilities, and the public [3].
In the study by Taiwo et al. [4], an evaluation of the initial critical configuration of the HTR-10 pebble bed reactor was made, as part of an effort to identify experimental tests of the highest priority, recover the data for those cases, and then develop standard problems (benchmarks) that are of sufficient quality for use in the licensing of the VHTR analysis codes. The initial criticality measurement in HTR-10 was judged as an acceptable benchmark.
Monte Carlo techniques have been widely used to study complex systems. Several HTR-10 criticality analyses have been performed using codes of the MCNP family. In the study by Şeker and Çolak [5], three different geometrical models were employed to study the effect of geometric details on the criticality calculations. Results were compared with diffusion theory-based calculations as well as with experimental data.
Certain studies [6] have demonstrated that the high temperature reactor (HTR) performance is considerably influenced by the way its geometry is modeled with the Monte Carlo approach. The spatial and energy shielding of the neutron flux even in such small particles cannot be neglected for important isotopes with high resonance cross sections.
In previous studies [7], the design of a pebble bed transmuter (PBT) was done. It is a subcritical system cooled by helium and moderated by graphite, which uses as fuel small amounts of transuranic elements in the form of TRISO particles, confined in 3 cm radius graphite pebbles forming a pebble bed configuration. The PBT is a device designed for the transmutation of nuclear waste from the existing LWR. In the study by Abánades et al. [8], the conceptual design of a transmutation advanced device for sustainable energy applications (TADSEA) was presented. It could be used for simultaneous nuclear waste transmutation and hydrogen generation.
In the study by García et al. [9], a geometrical method for calculating the real number of pebbles that fit in a cylindrical ADS (accelerator driven system) core, according to its size and pebble configuration, is described. Based on its results, the packing fraction influence on the TADSEA’s main work parameters was studied, and the redesign of the previous configuration was done in order to maintain the thermal power that was established in the preliminary design. Results showed the capability of the system to reach coolant outlet temperatures high enough for its application to hydrogen production.
In order to improve the conceptual design of the TADSEA, a new computational model for whole core calculation of PBRs was developed. In this paper, the features of a detailed geometric computational model for PBR whole core analysis, using MCNPX code, are presented. The model considers an explicit representation of the graphite reflector and the double heterogeneity of the coated fuel particles and the graphite pebbles. The validation of model was carried out using the HTR-10 benchmark. Geometrical and material characteristics of the HTR-10 can be found in IAEA-TECDOC-1382 [3]. The characteristics of fuel and dummy pebbles are shown in Table 1. Results were compared with experimental data and calculations of other authors.
The HTR-10 benchmark was divided into two parts: the Original and the Deviated benchmarks. Additionally, a sensitivity analysis of the parameters that differ from the Deviated benchmark to the Original benchmark was performed to evaluate their influence on the differences between the calculations and the experimental results.
In PBRs, the helium gas flows among closely packed pebbles that are randomly loaded into the core region. It is difficult to model this distinguishing feature of arbitrary arrangement for pebbles. In general, most of the neutronic and thermohydraulic simulations assume an ordered arrangement of pebbles, including the simple cubic (SC), body-centered cubic (BCC), or the face-centered cubic (FCC) arrangement [10].
The effects of BCC and FCC arrangements on the critical height of the HTR-10 core are investigated using the proposed model based on MCNPX. In addition, the critical height parameter was calculated using two models for the simulation of the TRISO particles inside the pebbles: one model that considers the boundary effect, in which the particles cut by the surfaces are left out of the problem, and another model in which those particles are entirely considered.
The investigated geometric computational model is based on MCNPX code, version 2.6e [11]. MCNP was created in 1994 by Los Alamos National Laboratory. It is based on the Monte Carlo probabilistic method for computational modelling of time dependent transport of many kinds of particles, with a wide range of energies, in much different geometries. The MCNPX version 2.6e incorporates new capabilities with respect to previous versions and the available library in XSDIR, ENDF/B VI.2, was used.
2. Computational Model
In order to get more accurate results in the simulation of PBRs core, the double heterogeneity of the system must be considered. It consists of randomly located pebbles into the core and randomly located TRISO particles into the fuel pebbles. These features are often neglected due to the difficulty to model them with MCNP codes. The main reason is that the software allows defining a limited number of cells and surfaces (99999) [12]. For instance, inside the HTR-10 core, there are in average 27000 pebbles and 8335 TRISO particles per pebble to be modeled. This composition corresponds to a density of 10.4 g/cm3 of metal oxide and 5 g of pebble mass per coated fuel particles (CFPs) with a kernels radius of 250 μm.
In several studies, the geometric model is simplified homogenizing the material zones [7–9]. In other studies, the CFPs are distributed within the pebbles in a regular way, as well as the pebbles inside the core [5, 13]. In the study by Abedi et al. [14], the CFPs truncated by the boundary of fuel region inside fuel pebble were included in the study.
2.1. Models for the Simulation of the TRISO Particles inside the Pebbles
Distributing the CFPs in a regular way within the pebbles is a suitable way to simulate the fuel inside the core in complex geometries like in the HTR-10. In this work, a detailed modeling of TRISO fuel was carried out. The layers covering a TRISO particle were modeled using the exact densities and dimensions. In order to get a detailed geometric model equivalent to the homogeneous model, it was necessary to calculate the number of CFPs that guarantees the desired fuel mass per pebble.
In the model, the CFPs are uniformly located in a hexagonal lattice within the pebble. First of all, the number of CFPs that can fill the plane where the hexagonal lattice is located () is calculated (Figure 2(a)). Secondly, according to , a repetitive structure consisting of a CFP and its five layers is created. Then, the repetitive structure is spread all over the plane (Figure 2(b)). This model (in this paper called Model I) has two major approaches: (1) to consider that the CFPs are uniformly arranged inside the pebble, which ignores the real randomness of the CFPs location, and (2) not to include in the geometry the fraction of CFPs intersected with the pebble’s surface, causing the considered fuel mass in the pebble to be less than the real one.
(a) XY plane
(b) XZ plane
The last situation can be avoided by removing the CFPs located near the boundary. A model was built with a uniform distribution of the CFPs inside the pebble where those CFPs intercepted by the pebble’s surface are visually removed from the array and relocated in the successive plane, causing all the CFPs to be fully considered inside the pebble (Model II). However, this model has an extra difficulty: it is hard to construct.
2.2. Models for the Simulation of the Pebbles inside the Core
In the first load, the PBR core was filled with moderator pebbles and fuel pebbles. Fuel pebbles contain CFPs in their inner zone. The moderator pebbles are made only of graphite. The fuel-to-moderator pebble ratio (57%–43%) specified in the HTR-10 benchmark cannot be exactly modeled using repetitive structures of MCNPX [12]. To solve this problem, the size of the dummy pebbles was changed in order to obtain the designed fuel-to-moderator ratio in the HTR-10 core (Table 2). This approximation was considered previously in IAEA-TECDOC-1382 [3].
In order to simulate the core’s geometry, two types of lattice were used: one with a tetrahedral arrangement in two planes (FCC) and the other is the BCC mentioned above. To build the FCC arrangement, a tetrahedron composed of six spheres at its base forming an equilateral triangle and one on the top was assembled base to base with a similar tetrahedron rotated 60° with respect to the first one (Figure 3). The unit cell is the symmetric cube formed by spheres at the corners and spheres on the faces. Four pebbles () compose the cell. There are two dummy pebbles with the reduced diameter and the other two ones are fuel pebbles. In order to keep the volumetric packing fraction (PF) established in the reference benchmark, the pebbles were conveniently spaced into the unit cell, and the dimensions of the unit cell were calculated according to that (see Table 2).
To fill the core volume with unit cells, a square lattice for each plane was created and it was repeated inside the cylinder that delimits the reactor core. To perform the axial distribution of pebbles, as the number of unit cells in the axial direction should be integer, the top plane that delimits the height of the core does not cut the fuel pebbles, and the moderator-fuel ratio is preserved. In the radial direction, the surface of the cylinder cuts the unit cells so there is a boundary effect that can influence the accuracy of the results.
The height values, where the eigenvalue problems were solved, were obtained as integer multipliers of the unit cell size. The HTR-10 core heights for FCC model are displayed in Table 3.
A cubic unit cell composed of one pebble in the center and eight eights in the corners was considered for the BCC configuration. Therefore, the unit cell is composed of two pebbles (), a fuel pebble and a dummy pebble (Figure 4).
(a)
(b)
(c)
Like in FCC arrangement, the dummy pebble diameter is reduced in order to guarantee the fuel-to-moderator ratio (57%–43%). Again, the cubic unit cell size was calculated to keep the established PF in the core (0.61) (Table 2). With the calculated size of the unit cell, the square repetitive structure was applied to fill the core volume (Figure 4).
Similarly, the height values to solve the critical height problem with the BCC arrangement were taken as integers multiples of the cubic unit cell size, to avoid the boundary effect in the axial direction (Table 4). This way allows minimizing the influence of this approximation on the results, because the fuel pebbles are entirely taken into account in the simulation process. In both FCC and BCC arrangements, the pebbles in the discharge cone were simulated with the same repetitive structure, where all the pebbles are composed by graphite, to guarantee the material continuity.
3. Benchmark Solutions
In order to validate the proposed model, the neutronic benchmark for the HTR-10, included in IAEA-TECDOC-1382 [3], was calculated. The Original and Deviated problems for the initial criticality, control rod worth for full and initial core, and the temperature coefficient were solved.
3.1. Benchmarks Description
3.1.1. Original Benchmark Description
The critical core’s height is obtained calculating the effective multiplication factor () for several heights of the core and interpolating the results. To determine the temperature effects and the rod worth, is calculated for different reactor states varying the temperature in the core or modifying the positions of control rods, respectively.
The benchmark was divided into four different subproblems:(1)initial criticality (benchmark problem B1),(2)temperature coefficient (benchmark problem B2),(3)control rods worth for full core (benchmark problem B3),(4)control rods worth for initial core (benchmark problem B4).
3.1.2. Deviated Benchmark Description
After the Original HTR-10 physics benchmark problems had been defined and before the initial core loading, two conditions were changed from the benchmark, with respect to the initial core loading and the first critical experiment. First, dummy balls that were to be loaded in the initial core were different from those that had been defined in the Original benchmark. The graphite pebbles prepared for the first criticality had two main differences, which influence the neutronic calculation: density and impurities. The graphite pebbles prepared had a density of 1.84 g/cm3, larger than the value defined in the Original benchmark (1.73 g/cm3), and the equivalent boron impurities were 0.125 ppm instead of 1.3 ppm.
Second, the first criticality experiment was made under atmospheric air, instead of helium, which was indicated in the Original benchmark.
When the participants of the CRP knew the deviations of the Original benchmark, some calculations were repeated with new conditions and both results were analyzed in IAEA-TECDOC-1382 [3]. In order to compare our results with the published ones, these deviations were considered, causing changes in the dummy pebble’s composition and in the core’s atmosphere. Other conditions were kept with respect to the Original benchmark.
3.2. Simulation Results for the Original and Deviated Benchmarks
3.2.1. Initial Criticality
The benchmark problem B1 consists in estimating the core critical height of the HTR-10. The height values were measured over the upper surface of the discharge cone. In the Original and Deviated benchmark problems, it was considered that the temperature of the coolant (helium or atmospheric air) was 20°C, and the control rods and absorbent balls were completely removed.
The Original and Deviated benchmark problems B1 were simulated using the proposed MCNPX model. A detailed geometry of the HTR-10 pebbled bed, using a repetitive structure with FCC arrangement and the established PF of 0.61, was considered. A fuel-to-moderator ratio of 57%–43% was guaranteed reducing the size of dummy pebbles. The double heterogeneity, that is, fuel pebbles in the core and TRISO particles inside the pebbles (Model I), was included. Besides, an exact representation of all the reflector graphite blocks with the holes for the coolant channels, control rods, and absorbent balls, which surround the pebbled bed, was considered. In the case of the Deviated benchmark, the changes indicated in Section 3.1.2 were done.
The values were calculated for each loading height of Table 3 and the critical height was obtained using the least square method. In order to obtain the results, a million neutron histories were simulated, which ensures a statistical error () less than 0.1%. The calculations were done using the cross sections libraries ENDF/B-VI with a core’s temperature of 27°C (300 K).
Table 5 shows the experimental and calculated results of critical height for the Original and Deviated benchmarks. The participant countries in the CRP used different code systems (probabilistic and deterministic). In the Original benchmark, our results presented a better agreement with the experimental critical height than other Monte Carlo calculations. For the Deviated benchmark, only few calculations were made with Monte Carlo codes. However, a good agreement of our model’s results was also obtained.
3.2.2. Temperature Coefficients
The scope of this problem was to calculate the of the full core (5 m3), at core’s temperatures of 20°C (B21), 120°C (B22), and 250°C (B23) with no control rods inserted. The standard cross sections libraries of MCNPX do not have these values of core’s temperatures ready; however, they could be obtained with NJOY code.
In this work, the calculations were done using the available cross sections libraries for MCNPX, version 2.6e, obtained from ENDF/BVI library at core’s temperatures of 27°C (300 K), 326.3°C (600 K), and 526.3°C (800 K), with helium or atmospheric air as coolant and no control rods or absorbent balls inserted. The values for the B22 and B23 problems were obtained by polynomial interpolation of the calculated values for the available temperature libraries. In this way, a comparative analysis of the behavior with the core’s temperature was possible.
Table 6 shows our results of as a function of the temperature for the Original and Deviated benchmarks. In addition, results obtained with other code systems (probabilistic and deterministic) are presented.
In general, the obtained values present good agreement with the initial calculations for all the temperature values. However, the values present a better agreement with the Chinese results, especially for the Deviated benchmark.
3.2.3. Calculation of Control Rods Worth for Full Core
The reactivity worth of ten control rods fully inserted (B31) and one control rod fully inserted (B32, with the other remaining rods in withdrawal position), for full core and under helium or air atmosphere at 20°C, was calculated.
In our case, the calculation of was done with ENDF/B-VI library with a core’s temperature of 27°C. The reactivity worth was calculated by where is the reactivity worth and and represent the values in the reactor when the rods are completely inserted and when the rods are in withdrawal position, respectively.
Table 7 shows the results of the B3 problem for the Original and Deviated benchmarks, as well as the results obtained with other code systems (probabilistic and deterministic) used by the participant countries in the CRP. In general, a good agreement with initial calculations was obtained, especially with the Chinese results.
3.2.4. Calculation of Control Rod Worth for Initial Core
This problem includes calculating the reactivity worth of the ten fully inserted control rods (B41) under helium or air atmosphere, with a core’s temperature of 20°C and a loading core height of 126 cm. In the same way, the worth of one rod (B42) was calculated.
Table 8 shows the results of the B4 problem for both Original and Deviated benchmarks. A good agreement between our results and the Monte Carlo Chinese results can be observed.
4. Sensitivity Analysis
A sensitivity analysis of the parameters that differ from the Deviated benchmark to the Original benchmark was performed, to evaluate their possible influence on the differences between the calculations and the experimental results.
One difference between the Deviated benchmark and the Original benchmark was given by the variation of the dummy pebble’s composition.
4.1. Density Variations in the Dummy Pebbles
To evaluate the influence of the dummy pebble’s density variation on the results of the critical heights, the graphite density was changed keeping constant the boron concentration value. The Deviated benchmark was chosen for this study because in this benchmark the established conditions were closer to the experimental conditions than in the Original benchmark. The initial conditions for the study were the following:(i)dummy pebbles density = 1.84 g/cm3;(ii)boron impurities in dummy pebbles = 0.125 ppm;(iii)initial core’s atmosphere = air. Not all other parameters were varied.
The dummy pebble density was uniformly varied in the established range for both benchmarks (1.73 g/cm3 1.84 g/cm3). With the chosen values, the parameter was calculated for each height, and the corresponding critical height was obtained by linear interpolation methods. The results are plotted in Figure 5.
As observed in Figure 5, when the dummy pebble density increases, the critical height value is lower. Increasing the density of dummy pebbles, the moderator-to-fuel mass ratio also increases, improving the multiplicative properties of the fresh-fueled core. A difference of more than four cm was obtained for the critical height.
4.2. Boron Concentration Variations in the Dummy Pebbles
Obviously, increasing the boron concentration of dummy pebbles produces a decrease in . However, to analyze how this increase influences the critical height, the boron concentration was uniformly varied until it achieved the dummy pebble’s boron concentration of the Original benchmark (1.3 ppm), keeping constant the density value of the Deviated benchmark. The initial conditions for the analysis were the following:(i)boron impurities in dummy pebbles = 0.125 ppm;(ii)dummy pebble’s density = 1.84 g/cm3;(iii)initial core’s atmosphere = air. Not all other parameters were varied.
The critical height behavior is shown in Figure 6; it changed less than 3 cm. The critical height of core was more sensitive to the dummy pebble density variation than to the change in the boron concentration of dummy pebbles for the analyzed values.
4.3. Variations in the Composition of the Core’s Atmosphere
Another parameter that could affect the results of the calculated benchmark is the uncertainty of air humidity that replaced the helium coolant of the core. For this reason, a comparison between the critical height results with the air composition given by the Deviated benchmark and with dry air was made. The humid air composition was calculated according to the specifications given in IAEA-TECDOC-1382 [3] for the Deviated benchmark and the dry air composition was obtained from [15]. The compositions for both humid and dry air are given in Table 9. The rest of the parameters of the Deviated benchmark were kept unchanged. The values of and the critical height for both compositions of the air are given in Table 10. The results showed that the composition of the air in the studied range does not affect the critical height. In the study by Taiwo et al. [4], the change in from saturated air to totally dry air was found to be . This result is consistent with our study.
4.4. Comparison between Critical Height Values Using the FCC and BCC Arrangements
The aim of this study was to evaluate the possible influence of the arrangement used, FCC or BCC, on the critical height results. It is known that pebbles are randomly loaded in the reactor’s core, but in our computational model a fixed arrangement is assumed. The calculations were made for the Original benchmark.
The BCC and FCC unit cells have two and four pebbles, respectively, and they were simulated taking into account the parameters described in Section 2.
Table 11 shows the results of for each configuration. The critical height values were calculated using cubic spline polynomial interpolation. The obtained critical height for the BCC arrangement was 126.15 cm and it was 127.61 cm for the FCC.
The results showed that there are not significant differences between both models. The calculated critical heights have a relative error of 1.1%. This allows concluding that both models are equivalent for simulating the pebble bed whole core. It is important to notice that for the FCC arrangement the calculations were carried out for ten height values, whereas for the BCC arrangement eleven values were employed because of the smaller unit cell size. The mismatch of the calculated points makes it impossible to compare them one to one.
4.5. Comparison between Models I and II
Model II was made to avoid the boundary effect, which occurs when the surface of fuel region inside fuel pebble cuts the TRISO particles located near to boundary of fuel region inside fuel pebble (see Section 2). Model I does not consider this feature. With the aim of evaluating the influence of this approximation, the critical height results calculated by both models were compared.
The BCC arrangement was used to simulate the pebbles inside the core. The values were calculated for the core height values of Table 4, and the critical height was obtained using interpolation techniques. The behavior as function of the core’s height, for the two models, is presented in Figure 7. It can be seen that almost all values of the two models are very close one to one and the confidence intervals for 68% () overlap (except for 130.6687 cm of height). The values of were calculated ensuring a standard deviation value below 100 pcm. Accordingly, the critical height values calculated by both models are the same.
We can conclude that the pebble boundary effects for the whole core calculation of PBRs are negligible and are only important for cell calculations. This is an important result, taking into account the additional difficulties to generate this special feature in the model.
5. Conclusions
The features of a detailed geometric computational model for pebble bed whole core analysis using MCNPX code, 2.6e version, were discussed. The double heterogeneity of the pebble bed (fuel pebbles inside the core and the TRISO particles inside the pebbles) was taken into account. The validation of model was carried out using the HTR-10 benchmark.
The results obtained with the model, corresponding to the solution of the four benchmark problems, are consistent with the experimental results and initial calculations of the participant countries of the Coordinated Research Program. A very good agreement with the Chinese calculations was obtained.
A sensitivity analysis of the parameters that differ from the Deviated benchmark to the Original benchmark was performed to evaluate their possible influence on the differences between the calculations and the experimental results. The increase of the dummy pebble density and dummy pebble boron concentration produces opposite effects on the behavior of the critical height. This can make the Original and Deviated benchmarks results not significantly different. The critical height of core was more sensitive to the change of the dummy pebble density than to the change of the dummy pebble boron concentration in the analyzed value range.
Another parameter that could affect the results of the benchmark is the uncertainty of air humidity that replaced the helium coolant of the core. The analysis made showed that the variation of the composition of the air in the studied range does not affect the value of the critical height.
The FCC and BCC configurations for the distribution of the pebbles inside the core were compared. The corresponding critical heights were calculated and there was no significant difference between them, presenting a relative error of 1.1%.
Two models were evaluated for the CFPs distribution inside the pebbles. Model I ignores the boundary effect produced by the intersection of the pebble’s surface and the CFPs located near the edge. In Model II, the CFPs near the edge are manually relocated in other planes to mitigate the boundary effect. The critical height of the HTR-10 core was calculated for both models and the results were the same, showing that the boundary effect in the pebble is not an influential parameter in whole core calculations and both models are equivalent. This is an important result, taking into account the additional effort that requires the relocation of the CFPs.
The new detailed geometric computational model for whole core calculation of PBRs using MCNPX, 2.6e version, allows performing more accurate neutronic analysis to improve the accuracy of the design of the TADSEA.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was made in the framework of the MES-CAPES project: “Study of Advanced Nuclear Systems of Pebbled Bed Type, for the Transmutation of Long Life Waste and Hydrogen Production with High Temperature Technologies.” The authors would like to thank the Brazilian agencies CAPES (Higher Education Coordination Agency of Education Ministry) and CNPq (National Council for Scientific and Technological Development) for the financial support for this work.