- Email: [email protected]
Fuel burnup analysis for the Moroccan TRIGA research reactor
Annals of Nuclear Energy 51 (2013) 112–119
Contents lists available at SciVerse ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Fuel burnup analysis for the Moroccan TRIGA research reactor B. El Bakkari a,b,⇑, T. El Bardouni b, B. Nacir a, C. El Younoussi a,b, Y. Boulaich a,b, H. Boukhal b, M. Zoubair b a b
Unité Conduite Réacteur, Centre d’Etudes Nucléaires de la Maâmora CNESTEN/CENM B.P.1382, R.P.10001, Rabat, Morocco ERSN-LMR, Department of Physics, Faculty of Sciences, Tetuan 93002, Morocco
a r t i c l e
i n f o
Article history: Received 24 April 2012 Received in revised form 9 July 2012 Accepted 21 July 2012 Available online 9 October 2012 Keywords: TRIGA Burnup BUCAL1 MCNP NJOY99
a b s t r a c t The fundamental advantage and main reason to use Monte Carlo methods for burnup calculations is the possibility to generate extremely accurate burnup dependent one group cross-sections and neutron fluxes for arbitrary core and fuel geometries. Yet, a set of values determined for a material at a given position and time remains accurate only in a local region, in which neutron spectrum and flux vary weakly — and only for a limited period of time, during which changes of the local isotopic composition are minor. This paper presents the approach of fuel burnup evaluation used at the Moroccan TRIGA MARK II research reactor. The approach is essentially based upon the utilization of BUCAL1, an in-house developed burnup code. BUCAL1 is a FORTRAN computer code designed to aid in analysis, prediction, and optimization of fuel burnup performance in nuclear reactors. The code was developed to incorporate the neutron absorption reaction tally information generated directly by MCNP5 code in the calculation of fissioned or neutron-transmuted isotopes for multi-fueled regions. The fuel cycle length and changes in several core parameters such as: core excess reactivity, control rods position, fluxes at the irradiation positions, axial and radial power factors and other parameters are estimated. Besides, this study gives valuable insight into the behavior of the reactor and will ensure better utilization and operation of the reactor during its life-time and it will allow the establishment of strategic planning for fuel management such as shuffling and/or reloading schemes and its safe implementation. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Adequate knowledge of burnup levels of fuel elements within a research reactor is of great importance for its optimum operation. Such knowledge is required for the monitoring of reactivity parameters, fluxes and power distributions throughout the reactor core, the estimation of the radioactive source term needed in accidental situations analysis, the evaluation of the amount of fissile materials present at any moment within the fuel for safeguards purposes and the estimation of cooling and shielding requirements for interim storage or transport of spent fuel elements. Accurate prediction of reactor main neutronics parameters requires accurate simulation of particle transport process by solving the Boltzmann transport equation. The Monte Carlo method (Spanier and Gelbard, 1969) is one of the most suitable analysis methods because of its accurate modeling of physical phenomena. In the last couple of decades, its application has been expanded to burnup calculations (Jafarikia et al., 2010) and its practicability has been shown. Most of these calculations are based on existing standard multipurpose continuous ⇑ Corresponding author at: Unité Conduite Réacteur, Centre d’Etudes Nucléaires de la Maâmora CNESTEN/CENM B.P.1382, R.P.10001, Rabat, Morocco. Tel.: +212 668782702; fax: +212 537803326. E-mail address: [email protected] (B. El Bakkari). 0306-4549/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2012.07.030
energy Monte Carlo codes, such as RACER, MVP and MCNP (Los Alamos Monte Carlo Group, 2003). The computational issues to perform accurate Monte Carlo calculations are however continuously reduced due to improvements made in the basic Monte Carlo algorithms, due to the development of variance reduction techniques and due to developments in computer architecture (more powerful processors, the so-called brute force approach through parallel processors, networked systems, etc.). This evolution of computer architecture is going to continue in the future. This paper presents a burnup dependent neutronic analysis of the 2 MW TRIGA Mark II Moroccan research reactor at Centre d’Etudes Nucléaire de la Maâmorra (CENM); burnup calculations were made using the in-house developed code BUCAL1. Such studies are needed all the way, from the preliminary analysis of new concepts to decommissioning, for a safe and economic operation throughout reactor life-time. For the modeling of the Moroccan TRIGA research reactor, the general purpose 3-D Monte Carlo N-particle transport code MCNP (version 5 release 1.40) was chosen because of its general modeling capability, correct representation of detailed geometry, transport effects and continuous energy cross sections. All fresh fuels, control rods, and other in-core (e.g., central thimble, graphite elements, etc.) and out-core elements (e.g., neutrons beams, thermal column,
113
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
etc.) were modeled using the maximum details allowed by the manufacturer General Atomics (GA). The repeated structure capability of MCNP was used to create a full core three dimensional model of TRIGA. The MCNP input was prepared in such a way that a very quick setup of any desired core configuration with an adequate position of all control rods is possible. In this study, Continuous energy cross section data from the more recent nuclear data evaluation ENDF/B-VII.0 (Chadwick et al., 2006) as well as S(a, b) thermal neutron scattering functions from the ENDF/B-VI.8 library were used. The cross section libraries were generated by using the updated version of the NJOY99 system (MacFarlane, 1999).
2. Modeling of TRIGA reactor The TRIGA reactor in CENM is a light water cooled, graphite-reflected one, designed for operation at a steady state thermal power level of 2000 kW. An outstanding feature of the TRIGA reactor is its proven safety, which stems from the large instantaneous negative temperature coefficient of reactivity of its U-ZrH fuel moderatormaterial. The TRIGA core consists of 101 fuel elements, 17 graphite elements, central thimble, a pneumatic transfer tube and a water filled position (G25) in the outer ring of the core. The cross-sectional view of the present core configuration of the reactor is
shown in Fig. 1 which was achieved on June 2007 during reactor start-up at full power operation. Elements are arranged in seven concentric rings in hexagonal geometry and the spaces between the rods are filled with water that acts as coolant and moderator. The reactor was modeled in full 3-D details to minimize geometry approximations. The repeated structure capability of MCNP was used to create a full core 3-D model of TRIGA. The TRIGA lattice can be represented as a hexagonal prism, solids with eight faces. The fuel elements were modeled explicitly specifying the detailed structure of the rod to eliminate any homogenization effects. The tapered end fixtures of stainless steel were also modeled with a very little approximation. The power level of the reactor is controlled with five control rods: a regulating rod and four shim safety rods. The control rods were explicitly modeled along the active length containing three vertical sections of boron carbide, fuel follower and void region. The central thimble was considered to be filled with water and the pneumatic tube was considered to be voided. The graphite dummy elements are of the same general dimensions and construction as the fuel-moderator elements, except these elements are made of aluminum alloy and filled entirely with graphite. The model was extended up radially containing the graphite reflector and lead shield. An annular well on the inside diameter in the top of the graphite reflector that provides for the rotary specimen rack was also modeled along with the radial and tangential beam ports that serve for experimentations around the TRIGA core. To complete the modeling of all reactor facilities the thermal column, which is a 1219 mm squared assembly located in the side of the reactor shield structure, is also modeled. It is located between beam ports NB1, NB4 and the reactor tank that consists of an aluminum vessel installed in the reactor shield structure, this facility serves for the irradiation of large experimental specimens. All the geometric and material data were taken from the fabrication shipment documentation provided by the reactor manufacturer General Atomics of USA. Physical proprieties of fuel and fuel follower elements are presented in Table 1. Figs. 2 and 3 represent the radial and axial view of MCNP5 model of the 2 MW TRIGA MARK II research reactor of CENM. An essential aspect of developing an accurate reactor physics model is validation. The accuracy of both the neutron transport physics as represented in MCNP code and the user-defined model must be assessed. However, MCNP has been proven to simulate the physical interactions correctly. Since the Monte Carlo method simulates individual particle tracks through a given system, it can provide a very accurate probabilistic transport solution. That does not mean that the model of TRIGA will provide accurate answers. Therefore, to build confidence, all the neutronics parameters including effective multiplication factor, radial and axial power peaking distributions as well as the reactivity worth and integral reactivity curves of the control rods were performed for the fresh core with MCNP modeling of the TRIGA reactor. The MCNP calculated neutronics parameters showed a very good agreement with
Table 1 Physical proprieties of fuel and fuel follower elements of TRIGA reactor.
Fig. 1. Cross-sectional view of the present core configuration of the Moroccan TRIGA reactor.
Outer diameter Fuel diameter Fuel height Cladding thickness Diameter of zirconium rod Hydrogen/zirconium ratio Amount of uranium in U-ZrH (wt.%) Enrichment (%)
Fuel element (cm)
Fuel follower (FF) (cm)
3.76 3.65 38.1 0.0508 0.6350 1.65 8.5 20
3.44 3.33 38.1 0.0508 0.6350
114
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
Fig. 2. The radial view of the MCNP model of TRIGA reactor.
Fig. 3. The axial view of the MCNP model of TRIGA reactor.
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
the measured ones and those available on the final safety analysis report (FSAR) (El Bakkari et al., 2010). 3. Computational methodology 3.1. Overview of BUCAL1 burnup code BUCAL1 (El Bakkari et al., 2009a,b, 2012) is a FORTRAN computer code designed to aid in analysis, prediction, and optimization of fuel burnup performance in nuclear reactors. The code was developed to incorporate the neutron absorption reaction tally information generated directly by MCNP5 code in the calculation of fissioned or neutron-transmuted isotopes for multi-fueled regions. This allows us to benefit of the full capabilities provided by MCNP and to incorporate them into burnup calculations in the aim to perform more accurate and robust treatment of the problem. Neutron transmutation, fission, and radioactive decay are included in the modeling of the production and removal terms for each isotope of interest. For a fueled region, neutron transmutation, fuel depletion, fission-product poisoning, actinide generation, burnable poison loading and depletion effects are included in the calculation. The code uses the fourth order Rung Kutta (Quarteroni et al., 2007) method with the predictor–corrector approach for the resolution of the depletion equation for more than 900 isotopes. Using BUCAL1 one can do: (1) standard burnup calculation, (2) burnup calculation followed by a space of time of cooling, (3) burnup calculation with shuffling of fueled regions and (4) burnup calculation with reloading new fresh fuels. Continuation card in BUCAL1 allows users to restart calculations after unwanted stops. In the current version of BUCAL1, the two groups of nuclides under consideration are: 1. actinides (ACTs) that contain heavy metal nuclides with atomic number Z P 90 and their decay daughters; and 2. fission products (FPs) produced by fissions and their decay/capture daughters. Specifically, the calculated reaction rates of fission products and actinides are shown in Table 2. Only the neutron capture cross section is considered for fission products since neutron absorption by fission products is primarily via (n, c) reaction (Xu, 2003). For actinides, four types of reactions are considered including capture, fission, (n, 2n) and (n, 3n) because they lead to the production of significant quantities of fission products and higher-mass actinides. Fig. 4 gives the simplified follow diagram of BUCAL1. 3.2. Continuous cross section library generation In this research work, the evaluated neutron data based on ENDF/B-VII was chosen as the data source. The process to construct the continuous cross section data libraries from the data source files is typically performed by the NJOY system. The procedure to process point-wise cross sections by each module of NJOY is shown in Fig. 5.
Table 2 MCNP-tallied reactions used in BUCAL1. Reaction type Actinides
Fission products
(n, (n, (n, (n,
c) f) 2n) 3n)
(n, c)
MCNP reaction identifier 102 6 16 17 102
115
Continuous cross section libraries were elaborated for 227 isotopes: 75 structural material isotopes and 147 fission products and actinides needed by BUCAL1 for burnup calculations. Cross sections for nuclear fuel isotopes were prepared for the average fuel temperature of 618°K, 393°K for the average cladding temperature and 313 °K for the bulk water temperature. The rest of structural materials (reflector, graphite elements, etc.) were considered at 300°K. The S(a, b) thermal scattering cross sections of bound nuclei (i.e. H in H2O, C in graphite, Zr in ZrH and H in HZr) were taken directly from the standard MCNP5 (release 1.40) cross section library respecting the operation conditions of the TRIGA reactor. These thermal scattering data are essential to accurately model the neutron interactions at energies below 4 eV.
4. Results and discussion Criticality calculation gave confidence to perform a set of burnup calculations for the core loaded with fresh LEU fuel. Burnup calculations were performed for realistic operating conditions at a power of 2 MW using our new elaborated burnup code BUCAL1 and the ENDF/B-VII nuclear data library. The initial burnup steps in the calculations were taken smaller to consider for the Xe and Sm build-up poisoning. This study was done using an HP proliant server IntelÒ Xeon™ CPU 3.40 GHz, 2 Gb RAM and 2 Mb Cache Memory under Linux Mandriva operating system. Calculations were performed with 3000 cycles of iterations on a nominal source size of 1000 particles per cycle in order to decrease statistical error estimates. Initial 50 cycles were skipped to insure homogeneous neutrons source distribution. The estimated statistical errors (1r) were reduced below 50 pcm for keff values.
4.1. TRIGA core life time estimation The calculation of the keff and its relationship with core burnup is of primary importance to determine the core life time. First the excess reactivity for the beginning of the core life was calculated by MCNP at 2 MW reactor power and was found to be 6.77$ (keff = 1.04975). Then burnup calculations were performed for the TRIGA core without changing the loading pattern to determine the core life of the primary fuel cycle and for the primary core configuration. The variation of reactivity as a function of total core burnup (MW h) is presented in Fig. 6. At the initial burnup time a sharp loss of reactivity of 3.52$ is observed which is particularly due to the build-up of 135Xe and 149 Sm. The concentration of these fission products poisons strongly influence the reactivity and eventually reaches equilibrium at about 150 MW h and 1500 MW h for 135Xe and 149Sm respectively, as can be seen from Figs. 7 and 8. Then, the excess reactivity researches zero at 3360 MW h of reactor operating history which correspond to 70 days of continuum operation at full power which is considered to be the maximum operational life time of the reactor at 2 MW, with the actual core configuration. According to the reactor manufacturer General Atomics (GA), the loss of reactivity due to 135Xe and 149Sm is estimated to be 3.57$ (FSAR), which is consistent with the results obtained using BUCAL1. The individual burnup (%235U depletion) of TRIGA fuel and fuelfollower elements at the end of life time (3360 MW h) are tabulated in Table 3. As expected, the maximum values around 7% are obtained for the hottest fuel elements of B ring (Fig. 2). Then the 235U depletion decreases when passing from the core centre to the core periphery until it reaches a value of 2.8% for fuel element G18 of G ring. The core average burnup is found to be 4.54%.
116
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
Fig. 4. Simplified flow diagram of BUCAL1.
4.2. Control rods position versus burnup In the following sections, steady states neutronics analysis were performed to determine the radial and axial specific power profiles as well as the neutron fluxes at the in-core irradiation facilities. To do this, we adopted the procedure used by Candalino (2006), after each burnup step, BUCAL1 creates an MCNP deck with new nuclide inventory for the desired materials to be burned (fuel and fuel follower elements). A new MCNP input was written using the new isotopic inventory. Then, the control rods were adjusted until a keff of 1.000 (±10 pcm) was achieved. This analysis was
performed at several burnup steps: 0, 960, 1920, 2400, 2880 and 3360 MW h. Fig. 9 shows the flow chart process to achieve the desired results and Fig. 10 gives the percent of the control rods withdrawing versus burnup supposing that all the CRs are in banked position. It shows that, the CRs withdrawing as a function of burnup is almost linear. Using a linear fit one found that in order to reach the full power of the reactor at 3744 MW h all the CRs must be fully withdrawn, which consistent with the reactor life time estimated on paragraph (4.1). Practically and due to safety reasons, this operating condition cannot be achieved with the same core configuration.
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
Fig. 5. Flow diagram of NJOY99 processing system for ACE format library construction.
Fig. 8.
149
Sm build-up as a function of core burnup.
Fig. 6. Excess of reactivity ($) for TRIGA reactor as a function of burnup (MW h).
4.3. Specific power values versus burnup The relative radial and axial power factors of fuel and fuel follower elements are presented in Figs. 11 and 12, respectively. The calculations are carried out at 0 MW h (BOC) and 3360 MW h (EOC).
Fig. 7.
135
Xe build-up as a function of core burnup.
Fig. 9. Flowchart for determining new CRs positions.
117
118
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
Fig. 10. CRs positions versus burnup (MW h).
Fig. 13. Thermal neutron flux at the CT.
Fig. 14. Thermal neutron flux at the G25 irradiation position. Fig. 11. Power factors distribution within the fuel and fuel – follower elements of the TRIGA reactor core.
Radially, it seems that for both BOC and EOC distributions, the maximum power factors are found to be for fuel elements located in B ring, due to the relatively higher thermalization of neutrons in the central region of the core.
Fig. 12. Axial power factors within the hottest fuel element.
At BOC, for each ring, the power factors have a cosines shape that peaks at the fuel elements farther than the control rods, for example for B ring the maximum is found to be 1.72 at B3 and for C ring the maximum is found to be 1.52 at C6. The amplitude of the shape is higher in rings B to D due to the influence of the control rods which are withdrawn with 55%. At EOC, one can see that the power factors values of rings B to D decrease considerably in average value, due to high burnup values at this region of the core and the cosines shape disappear partially due to the new control rods positions that became roughly fully withdrawn. EOC/BOC ratios show that, the power factors of fuel elements nearest to CRs become higher than in BOC due to the neutron flux evolution caused by the CRs withdrawing. Axially, the hottest fuel element was divided into 40 axial subzones to provide of the accurate heating profile of the fuel element. Fig. 12 shows that the hottest fuel region at BOC can be found at 5 cm below the central region with a value of 1.34, the location of this highest zone tends to shift upward as a function of burnup, because of the new CRs position to achieve criticality. As the fuel is depleted, positive reactivity insertion in the form of the CR withdrawal is needed to remain critical. It is also interesting to note that, when the core undergoes further burnup, the CRs need to be withdrawn and a flattening of the axial power profile can be observed. Thus the axial power peak decreases as shown in Fig. 12. The maximum peak value at EOC is found to be 1.28 at the fuel centreline.
119
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119 Table 3 Individual TRIGA fuel burnup. Fuel ID
%235U depletion
Fuel ID
%235U depletion
Fuel ID
%235U depletion
Fuel ID
%235U depletion
Fuel ID
%235U depletion
B01 B02 B03 B04 B05 B06 C01 C02 C03 C04 C05 C06 C07 C08 C09 C10 C11 C12 D01 D02 D03
7.06 7.10 7.07 7.01 6.99 7.01 6.21 6.37 6.24 6.37 6.12 6.30 6.13 6.24 6.02 6.20 6.10 6.27 5.54 5.68 5.67
D04 D05 D06 D07 D08 D09 D10 D11 D12 D13 D14 D15 D16 D17 D18 E01 E02 E03 E04 E05 E06
5.59 5.68 5.60 5.16 5.51 5.55 5.41 5.51 5.48 5.32 5.44 5.46 5.35 5.54 5.60 4.17 4.40 4.51 4.41 4.23 4.43
E07 E08 E09 E10 E11 E12 E13 E14 E15 E16 E17 E18 E19 E20 E21 E22 E23 E24 F01 F02 F03
4.55 4.30 4.03 4.23 4.36 4.25 4.01 4.23 4.32 4.15 3.89 4.13 4.24 4.15 3.98 4.26 4.40 4.33 3.25 3.48 3.77
F04 F05 F06 F07 F08 F09 F10 F11 F12 F13 F14 F15 F16 F17 F18 F19 F20 F21 F22 F23 F24
3.9 3.54 3.21 3.58 3.71 3.78 3.53 3.14 3.37 3.53 3.43 3.23 3.02 3.34 3.51 3.56 3.34 2.93 3.16 3.46 3.61
F25 F26 F27 F28 F29 F30 G02 G05 G07 G08 G13 G14 G17 G18 G22 G27 G29
3.42 3.14 3.30 3.67 3.72 3.47 3.09 2.86 3.00 3.08 2.53 2.39 2.75 2.86 2.63 2.94 3.02
4.4. Neutron fluxes at the irradiation Figs. 13 and 14 present the thermal fluxes versus burnup at the central thimble (CT) and G25 irradiation positions of the TRIGA reactor, respectively. For CT, it can be shown that the thermal flux is flattening as a function of burnup producing a decrease of the peak value of about 7% between BOC and EOC, due to the high burnup values at the core centre and the CRs positions. The peak value shift upward from 5 cm to 1 cm below the core centre line. For G25 position, it seems that thermal flux is also flattening with a small increase of the peak value at EOC. This can be explained as follows: more of the fissile material (U235) is depleted at the inner region of the core, more power must be supplied by the fuel elements located at the core periphery since the total power remain the same (2 MW), hence a small increase of the thermal flux at the outer part of the core.
5. Conclusion This work aimed to study the burnup dependent neutronic parameters for the 2 MW Moroccan TRIGA research reactor. Burnup calculations were made by means of the in-house developed and validated burnup code BUCAL1. The computational results show that, the reactor life time was found to be 3360 MW h (considering full power reactor operations conditions). This result was confirmed by the CRs positions versus burnup, in which we found that at 3744 MW h all the CRs must be completely withdrawn in order to reach 2 MW. It is also shown that the radial and axial power distributions are strongly affected by burnup due to the CRs positions to compensate the loss of reactivity. Analysis of the
neutron fluxes at CT and G25 irradiation positions versus burnup shows that the thermal flux for CT can be reduced by 7% due to burnup and it can be increased slightly in G25, due to the fact that more power must be supplied by fuel elements in the core periphery in order to reach the full power at the end of life. References Candalino, R.W., 2006. Engineering analysis of low enriched uranium fuel using improved zirconium hybrid cross sections. Thesis, Texas A&M University, August 2006. Chadwick, M.B., Oblozinsky, P., Herman, M., et al., 2006. ENDF/B-VII.0: next generation evaluated nuclear data library for nuclear science and technology. Nuclear data Sheets 107, 2931–3060. El Bakkari, B. et al., 2009a. Development of an MCNP-tally based burnup code and validation through PWR benchmark exercises. Annals of Nuclear Energy 36, 626–633. El Bakkari, B. et al., 2009b. Validation of a new continuous Monte Carlo burnup code using a MOX fuel assembly. Nuclear Engineering and Design 239, 1828–1838. El Bakkari, B. et al., 2010. Monte Carlo modeling of TRIGA research reactor. Radiation Physics and Chemistry 79, 1022–1030. El Bakkari, B. et al., 2012. Accuracy assessment of a new Monte Carlo based burnup computer code. Annals of Nuclear Energy 45, 29–36. Los Alamos Monte Carlo Group, 2003. MCNP-A general Monte Carlo N-particle transport code version 5. Los Alamos National Laboratory, vols. I–III. Available from Radiation Safety Information Computational Center at Oak Ridge National Laboratory as CCC-740. MacFarlane, R.E., 1999. NJOY-99 nuclear data processing system. (update web site). Quarteroni, A., et Sacco, R., Saleri, F., 2007. Méthodes Numériques: Algorithmes, Analyse et Applications. Springer-Verlag Italia Via Decembrio 28 20138 Milano, Italia. ISBN 13 978-88-470-0495-5. Jafarikia, S. et al., 2010. Validation of IRBURN calculation code system through burnup benchmark analysis. Annals of Nuclear Energy 37, 325–331. Spanier, J., Gelbard, E.M., 1969. Monte Carlo Principals and Neutron Transport Problems. Addison-Wesley Publishing Company, Inc., Reading, Mass. Xu, Zh., 2003. Design strategies for optimizing high burnup fuel in pressurized water reactors. Doctor thesis, Massachusetts Institute of Technology (January 2003).
Contents lists available at SciVerse ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Fuel burnup analysis for the Moroccan TRIGA research reactor B. El Bakkari a,b,⇑, T. El Bardouni b, B. Nacir a, C. El Younoussi a,b, Y. Boulaich a,b, H. Boukhal b, M. Zoubair b a b
Unité Conduite Réacteur, Centre d’Etudes Nucléaires de la Maâmora CNESTEN/CENM B.P.1382, R.P.10001, Rabat, Morocco ERSN-LMR, Department of Physics, Faculty of Sciences, Tetuan 93002, Morocco
a r t i c l e
i n f o
Article history: Received 24 April 2012 Received in revised form 9 July 2012 Accepted 21 July 2012 Available online 9 October 2012 Keywords: TRIGA Burnup BUCAL1 MCNP NJOY99
a b s t r a c t The fundamental advantage and main reason to use Monte Carlo methods for burnup calculations is the possibility to generate extremely accurate burnup dependent one group cross-sections and neutron fluxes for arbitrary core and fuel geometries. Yet, a set of values determined for a material at a given position and time remains accurate only in a local region, in which neutron spectrum and flux vary weakly — and only for a limited period of time, during which changes of the local isotopic composition are minor. This paper presents the approach of fuel burnup evaluation used at the Moroccan TRIGA MARK II research reactor. The approach is essentially based upon the utilization of BUCAL1, an in-house developed burnup code. BUCAL1 is a FORTRAN computer code designed to aid in analysis, prediction, and optimization of fuel burnup performance in nuclear reactors. The code was developed to incorporate the neutron absorption reaction tally information generated directly by MCNP5 code in the calculation of fissioned or neutron-transmuted isotopes for multi-fueled regions. The fuel cycle length and changes in several core parameters such as: core excess reactivity, control rods position, fluxes at the irradiation positions, axial and radial power factors and other parameters are estimated. Besides, this study gives valuable insight into the behavior of the reactor and will ensure better utilization and operation of the reactor during its life-time and it will allow the establishment of strategic planning for fuel management such as shuffling and/or reloading schemes and its safe implementation. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Adequate knowledge of burnup levels of fuel elements within a research reactor is of great importance for its optimum operation. Such knowledge is required for the monitoring of reactivity parameters, fluxes and power distributions throughout the reactor core, the estimation of the radioactive source term needed in accidental situations analysis, the evaluation of the amount of fissile materials present at any moment within the fuel for safeguards purposes and the estimation of cooling and shielding requirements for interim storage or transport of spent fuel elements. Accurate prediction of reactor main neutronics parameters requires accurate simulation of particle transport process by solving the Boltzmann transport equation. The Monte Carlo method (Spanier and Gelbard, 1969) is one of the most suitable analysis methods because of its accurate modeling of physical phenomena. In the last couple of decades, its application has been expanded to burnup calculations (Jafarikia et al., 2010) and its practicability has been shown. Most of these calculations are based on existing standard multipurpose continuous ⇑ Corresponding author at: Unité Conduite Réacteur, Centre d’Etudes Nucléaires de la Maâmora CNESTEN/CENM B.P.1382, R.P.10001, Rabat, Morocco. Tel.: +212 668782702; fax: +212 537803326. E-mail address: [email protected] (B. El Bakkari). 0306-4549/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2012.07.030
energy Monte Carlo codes, such as RACER, MVP and MCNP (Los Alamos Monte Carlo Group, 2003). The computational issues to perform accurate Monte Carlo calculations are however continuously reduced due to improvements made in the basic Monte Carlo algorithms, due to the development of variance reduction techniques and due to developments in computer architecture (more powerful processors, the so-called brute force approach through parallel processors, networked systems, etc.). This evolution of computer architecture is going to continue in the future. This paper presents a burnup dependent neutronic analysis of the 2 MW TRIGA Mark II Moroccan research reactor at Centre d’Etudes Nucléaire de la Maâmorra (CENM); burnup calculations were made using the in-house developed code BUCAL1. Such studies are needed all the way, from the preliminary analysis of new concepts to decommissioning, for a safe and economic operation throughout reactor life-time. For the modeling of the Moroccan TRIGA research reactor, the general purpose 3-D Monte Carlo N-particle transport code MCNP (version 5 release 1.40) was chosen because of its general modeling capability, correct representation of detailed geometry, transport effects and continuous energy cross sections. All fresh fuels, control rods, and other in-core (e.g., central thimble, graphite elements, etc.) and out-core elements (e.g., neutrons beams, thermal column,
113
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
etc.) were modeled using the maximum details allowed by the manufacturer General Atomics (GA). The repeated structure capability of MCNP was used to create a full core three dimensional model of TRIGA. The MCNP input was prepared in such a way that a very quick setup of any desired core configuration with an adequate position of all control rods is possible. In this study, Continuous energy cross section data from the more recent nuclear data evaluation ENDF/B-VII.0 (Chadwick et al., 2006) as well as S(a, b) thermal neutron scattering functions from the ENDF/B-VI.8 library were used. The cross section libraries were generated by using the updated version of the NJOY99 system (MacFarlane, 1999).
2. Modeling of TRIGA reactor The TRIGA reactor in CENM is a light water cooled, graphite-reflected one, designed for operation at a steady state thermal power level of 2000 kW. An outstanding feature of the TRIGA reactor is its proven safety, which stems from the large instantaneous negative temperature coefficient of reactivity of its U-ZrH fuel moderatormaterial. The TRIGA core consists of 101 fuel elements, 17 graphite elements, central thimble, a pneumatic transfer tube and a water filled position (G25) in the outer ring of the core. The cross-sectional view of the present core configuration of the reactor is
shown in Fig. 1 which was achieved on June 2007 during reactor start-up at full power operation. Elements are arranged in seven concentric rings in hexagonal geometry and the spaces between the rods are filled with water that acts as coolant and moderator. The reactor was modeled in full 3-D details to minimize geometry approximations. The repeated structure capability of MCNP was used to create a full core 3-D model of TRIGA. The TRIGA lattice can be represented as a hexagonal prism, solids with eight faces. The fuel elements were modeled explicitly specifying the detailed structure of the rod to eliminate any homogenization effects. The tapered end fixtures of stainless steel were also modeled with a very little approximation. The power level of the reactor is controlled with five control rods: a regulating rod and four shim safety rods. The control rods were explicitly modeled along the active length containing three vertical sections of boron carbide, fuel follower and void region. The central thimble was considered to be filled with water and the pneumatic tube was considered to be voided. The graphite dummy elements are of the same general dimensions and construction as the fuel-moderator elements, except these elements are made of aluminum alloy and filled entirely with graphite. The model was extended up radially containing the graphite reflector and lead shield. An annular well on the inside diameter in the top of the graphite reflector that provides for the rotary specimen rack was also modeled along with the radial and tangential beam ports that serve for experimentations around the TRIGA core. To complete the modeling of all reactor facilities the thermal column, which is a 1219 mm squared assembly located in the side of the reactor shield structure, is also modeled. It is located between beam ports NB1, NB4 and the reactor tank that consists of an aluminum vessel installed in the reactor shield structure, this facility serves for the irradiation of large experimental specimens. All the geometric and material data were taken from the fabrication shipment documentation provided by the reactor manufacturer General Atomics of USA. Physical proprieties of fuel and fuel follower elements are presented in Table 1. Figs. 2 and 3 represent the radial and axial view of MCNP5 model of the 2 MW TRIGA MARK II research reactor of CENM. An essential aspect of developing an accurate reactor physics model is validation. The accuracy of both the neutron transport physics as represented in MCNP code and the user-defined model must be assessed. However, MCNP has been proven to simulate the physical interactions correctly. Since the Monte Carlo method simulates individual particle tracks through a given system, it can provide a very accurate probabilistic transport solution. That does not mean that the model of TRIGA will provide accurate answers. Therefore, to build confidence, all the neutronics parameters including effective multiplication factor, radial and axial power peaking distributions as well as the reactivity worth and integral reactivity curves of the control rods were performed for the fresh core with MCNP modeling of the TRIGA reactor. The MCNP calculated neutronics parameters showed a very good agreement with
Table 1 Physical proprieties of fuel and fuel follower elements of TRIGA reactor.
Fig. 1. Cross-sectional view of the present core configuration of the Moroccan TRIGA reactor.
Outer diameter Fuel diameter Fuel height Cladding thickness Diameter of zirconium rod Hydrogen/zirconium ratio Amount of uranium in U-ZrH (wt.%) Enrichment (%)
Fuel element (cm)
Fuel follower (FF) (cm)
3.76 3.65 38.1 0.0508 0.6350 1.65 8.5 20
3.44 3.33 38.1 0.0508 0.6350
114
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
Fig. 2. The radial view of the MCNP model of TRIGA reactor.
Fig. 3. The axial view of the MCNP model of TRIGA reactor.
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
the measured ones and those available on the final safety analysis report (FSAR) (El Bakkari et al., 2010). 3. Computational methodology 3.1. Overview of BUCAL1 burnup code BUCAL1 (El Bakkari et al., 2009a,b, 2012) is a FORTRAN computer code designed to aid in analysis, prediction, and optimization of fuel burnup performance in nuclear reactors. The code was developed to incorporate the neutron absorption reaction tally information generated directly by MCNP5 code in the calculation of fissioned or neutron-transmuted isotopes for multi-fueled regions. This allows us to benefit of the full capabilities provided by MCNP and to incorporate them into burnup calculations in the aim to perform more accurate and robust treatment of the problem. Neutron transmutation, fission, and radioactive decay are included in the modeling of the production and removal terms for each isotope of interest. For a fueled region, neutron transmutation, fuel depletion, fission-product poisoning, actinide generation, burnable poison loading and depletion effects are included in the calculation. The code uses the fourth order Rung Kutta (Quarteroni et al., 2007) method with the predictor–corrector approach for the resolution of the depletion equation for more than 900 isotopes. Using BUCAL1 one can do: (1) standard burnup calculation, (2) burnup calculation followed by a space of time of cooling, (3) burnup calculation with shuffling of fueled regions and (4) burnup calculation with reloading new fresh fuels. Continuation card in BUCAL1 allows users to restart calculations after unwanted stops. In the current version of BUCAL1, the two groups of nuclides under consideration are: 1. actinides (ACTs) that contain heavy metal nuclides with atomic number Z P 90 and their decay daughters; and 2. fission products (FPs) produced by fissions and their decay/capture daughters. Specifically, the calculated reaction rates of fission products and actinides are shown in Table 2. Only the neutron capture cross section is considered for fission products since neutron absorption by fission products is primarily via (n, c) reaction (Xu, 2003). For actinides, four types of reactions are considered including capture, fission, (n, 2n) and (n, 3n) because they lead to the production of significant quantities of fission products and higher-mass actinides. Fig. 4 gives the simplified follow diagram of BUCAL1. 3.2. Continuous cross section library generation In this research work, the evaluated neutron data based on ENDF/B-VII was chosen as the data source. The process to construct the continuous cross section data libraries from the data source files is typically performed by the NJOY system. The procedure to process point-wise cross sections by each module of NJOY is shown in Fig. 5.
Table 2 MCNP-tallied reactions used in BUCAL1. Reaction type Actinides
Fission products
(n, (n, (n, (n,
c) f) 2n) 3n)
(n, c)
MCNP reaction identifier 102 6 16 17 102
115
Continuous cross section libraries were elaborated for 227 isotopes: 75 structural material isotopes and 147 fission products and actinides needed by BUCAL1 for burnup calculations. Cross sections for nuclear fuel isotopes were prepared for the average fuel temperature of 618°K, 393°K for the average cladding temperature and 313 °K for the bulk water temperature. The rest of structural materials (reflector, graphite elements, etc.) were considered at 300°K. The S(a, b) thermal scattering cross sections of bound nuclei (i.e. H in H2O, C in graphite, Zr in ZrH and H in HZr) were taken directly from the standard MCNP5 (release 1.40) cross section library respecting the operation conditions of the TRIGA reactor. These thermal scattering data are essential to accurately model the neutron interactions at energies below 4 eV.
4. Results and discussion Criticality calculation gave confidence to perform a set of burnup calculations for the core loaded with fresh LEU fuel. Burnup calculations were performed for realistic operating conditions at a power of 2 MW using our new elaborated burnup code BUCAL1 and the ENDF/B-VII nuclear data library. The initial burnup steps in the calculations were taken smaller to consider for the Xe and Sm build-up poisoning. This study was done using an HP proliant server IntelÒ Xeon™ CPU 3.40 GHz, 2 Gb RAM and 2 Mb Cache Memory under Linux Mandriva operating system. Calculations were performed with 3000 cycles of iterations on a nominal source size of 1000 particles per cycle in order to decrease statistical error estimates. Initial 50 cycles were skipped to insure homogeneous neutrons source distribution. The estimated statistical errors (1r) were reduced below 50 pcm for keff values.
4.1. TRIGA core life time estimation The calculation of the keff and its relationship with core burnup is of primary importance to determine the core life time. First the excess reactivity for the beginning of the core life was calculated by MCNP at 2 MW reactor power and was found to be 6.77$ (keff = 1.04975). Then burnup calculations were performed for the TRIGA core without changing the loading pattern to determine the core life of the primary fuel cycle and for the primary core configuration. The variation of reactivity as a function of total core burnup (MW h) is presented in Fig. 6. At the initial burnup time a sharp loss of reactivity of 3.52$ is observed which is particularly due to the build-up of 135Xe and 149 Sm. The concentration of these fission products poisons strongly influence the reactivity and eventually reaches equilibrium at about 150 MW h and 1500 MW h for 135Xe and 149Sm respectively, as can be seen from Figs. 7 and 8. Then, the excess reactivity researches zero at 3360 MW h of reactor operating history which correspond to 70 days of continuum operation at full power which is considered to be the maximum operational life time of the reactor at 2 MW, with the actual core configuration. According to the reactor manufacturer General Atomics (GA), the loss of reactivity due to 135Xe and 149Sm is estimated to be 3.57$ (FSAR), which is consistent with the results obtained using BUCAL1. The individual burnup (%235U depletion) of TRIGA fuel and fuelfollower elements at the end of life time (3360 MW h) are tabulated in Table 3. As expected, the maximum values around 7% are obtained for the hottest fuel elements of B ring (Fig. 2). Then the 235U depletion decreases when passing from the core centre to the core periphery until it reaches a value of 2.8% for fuel element G18 of G ring. The core average burnup is found to be 4.54%.
116
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
Fig. 4. Simplified flow diagram of BUCAL1.
4.2. Control rods position versus burnup In the following sections, steady states neutronics analysis were performed to determine the radial and axial specific power profiles as well as the neutron fluxes at the in-core irradiation facilities. To do this, we adopted the procedure used by Candalino (2006), after each burnup step, BUCAL1 creates an MCNP deck with new nuclide inventory for the desired materials to be burned (fuel and fuel follower elements). A new MCNP input was written using the new isotopic inventory. Then, the control rods were adjusted until a keff of 1.000 (±10 pcm) was achieved. This analysis was
performed at several burnup steps: 0, 960, 1920, 2400, 2880 and 3360 MW h. Fig. 9 shows the flow chart process to achieve the desired results and Fig. 10 gives the percent of the control rods withdrawing versus burnup supposing that all the CRs are in banked position. It shows that, the CRs withdrawing as a function of burnup is almost linear. Using a linear fit one found that in order to reach the full power of the reactor at 3744 MW h all the CRs must be fully withdrawn, which consistent with the reactor life time estimated on paragraph (4.1). Practically and due to safety reasons, this operating condition cannot be achieved with the same core configuration.
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
Fig. 5. Flow diagram of NJOY99 processing system for ACE format library construction.
Fig. 8.
149
Sm build-up as a function of core burnup.
Fig. 6. Excess of reactivity ($) for TRIGA reactor as a function of burnup (MW h).
4.3. Specific power values versus burnup The relative radial and axial power factors of fuel and fuel follower elements are presented in Figs. 11 and 12, respectively. The calculations are carried out at 0 MW h (BOC) and 3360 MW h (EOC).
Fig. 7.
135
Xe build-up as a function of core burnup.
Fig. 9. Flowchart for determining new CRs positions.
117
118
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119
Fig. 10. CRs positions versus burnup (MW h).
Fig. 13. Thermal neutron flux at the CT.
Fig. 14. Thermal neutron flux at the G25 irradiation position. Fig. 11. Power factors distribution within the fuel and fuel – follower elements of the TRIGA reactor core.
Radially, it seems that for both BOC and EOC distributions, the maximum power factors are found to be for fuel elements located in B ring, due to the relatively higher thermalization of neutrons in the central region of the core.
Fig. 12. Axial power factors within the hottest fuel element.
At BOC, for each ring, the power factors have a cosines shape that peaks at the fuel elements farther than the control rods, for example for B ring the maximum is found to be 1.72 at B3 and for C ring the maximum is found to be 1.52 at C6. The amplitude of the shape is higher in rings B to D due to the influence of the control rods which are withdrawn with 55%. At EOC, one can see that the power factors values of rings B to D decrease considerably in average value, due to high burnup values at this region of the core and the cosines shape disappear partially due to the new control rods positions that became roughly fully withdrawn. EOC/BOC ratios show that, the power factors of fuel elements nearest to CRs become higher than in BOC due to the neutron flux evolution caused by the CRs withdrawing. Axially, the hottest fuel element was divided into 40 axial subzones to provide of the accurate heating profile of the fuel element. Fig. 12 shows that the hottest fuel region at BOC can be found at 5 cm below the central region with a value of 1.34, the location of this highest zone tends to shift upward as a function of burnup, because of the new CRs position to achieve criticality. As the fuel is depleted, positive reactivity insertion in the form of the CR withdrawal is needed to remain critical. It is also interesting to note that, when the core undergoes further burnup, the CRs need to be withdrawn and a flattening of the axial power profile can be observed. Thus the axial power peak decreases as shown in Fig. 12. The maximum peak value at EOC is found to be 1.28 at the fuel centreline.
119
B. El Bakkari et al. / Annals of Nuclear Energy 51 (2013) 112–119 Table 3 Individual TRIGA fuel burnup. Fuel ID
%235U depletion
Fuel ID
%235U depletion
Fuel ID
%235U depletion
Fuel ID
%235U depletion
Fuel ID
%235U depletion
B01 B02 B03 B04 B05 B06 C01 C02 C03 C04 C05 C06 C07 C08 C09 C10 C11 C12 D01 D02 D03
7.06 7.10 7.07 7.01 6.99 7.01 6.21 6.37 6.24 6.37 6.12 6.30 6.13 6.24 6.02 6.20 6.10 6.27 5.54 5.68 5.67
D04 D05 D06 D07 D08 D09 D10 D11 D12 D13 D14 D15 D16 D17 D18 E01 E02 E03 E04 E05 E06
5.59 5.68 5.60 5.16 5.51 5.55 5.41 5.51 5.48 5.32 5.44 5.46 5.35 5.54 5.60 4.17 4.40 4.51 4.41 4.23 4.43
E07 E08 E09 E10 E11 E12 E13 E14 E15 E16 E17 E18 E19 E20 E21 E22 E23 E24 F01 F02 F03
4.55 4.30 4.03 4.23 4.36 4.25 4.01 4.23 4.32 4.15 3.89 4.13 4.24 4.15 3.98 4.26 4.40 4.33 3.25 3.48 3.77
F04 F05 F06 F07 F08 F09 F10 F11 F12 F13 F14 F15 F16 F17 F18 F19 F20 F21 F22 F23 F24
3.9 3.54 3.21 3.58 3.71 3.78 3.53 3.14 3.37 3.53 3.43 3.23 3.02 3.34 3.51 3.56 3.34 2.93 3.16 3.46 3.61
F25 F26 F27 F28 F29 F30 G02 G05 G07 G08 G13 G14 G17 G18 G22 G27 G29
3.42 3.14 3.30 3.67 3.72 3.47 3.09 2.86 3.00 3.08 2.53 2.39 2.75 2.86 2.63 2.94 3.02
4.4. Neutron fluxes at the irradiation Figs. 13 and 14 present the thermal fluxes versus burnup at the central thimble (CT) and G25 irradiation positions of the TRIGA reactor, respectively. For CT, it can be shown that the thermal flux is flattening as a function of burnup producing a decrease of the peak value of about 7% between BOC and EOC, due to the high burnup values at the core centre and the CRs positions. The peak value shift upward from 5 cm to 1 cm below the core centre line. For G25 position, it seems that thermal flux is also flattening with a small increase of the peak value at EOC. This can be explained as follows: more of the fissile material (U235) is depleted at the inner region of the core, more power must be supplied by the fuel elements located at the core periphery since the total power remain the same (2 MW), hence a small increase of the thermal flux at the outer part of the core.
5. Conclusion This work aimed to study the burnup dependent neutronic parameters for the 2 MW Moroccan TRIGA research reactor. Burnup calculations were made by means of the in-house developed and validated burnup code BUCAL1. The computational results show that, the reactor life time was found to be 3360 MW h (considering full power reactor operations conditions). This result was confirmed by the CRs positions versus burnup, in which we found that at 3744 MW h all the CRs must be completely withdrawn in order to reach 2 MW. It is also shown that the radial and axial power distributions are strongly affected by burnup due to the CRs positions to compensate the loss of reactivity. Analysis of the
neutron fluxes at CT and G25 irradiation positions versus burnup shows that the thermal flux for CT can be reduced by 7% due to burnup and it can be increased slightly in G25, due to the fact that more power must be supplied by fuel elements in the core periphery in order to reach the full power at the end of life. References Candalino, R.W., 2006. Engineering analysis of low enriched uranium fuel using improved zirconium hybrid cross sections. Thesis, Texas A&M University, August 2006. Chadwick, M.B., Oblozinsky, P., Herman, M., et al., 2006. ENDF/B-VII.0: next generation evaluated nuclear data library for nuclear science and technology. Nuclear data Sheets 107, 2931–3060. El Bakkari, B. et al., 2009a. Development of an MCNP-tally based burnup code and validation through PWR benchmark exercises. Annals of Nuclear Energy 36, 626–633. El Bakkari, B. et al., 2009b. Validation of a new continuous Monte Carlo burnup code using a MOX fuel assembly. Nuclear Engineering and Design 239, 1828–1838. El Bakkari, B. et al., 2010. Monte Carlo modeling of TRIGA research reactor. Radiation Physics and Chemistry 79, 1022–1030. El Bakkari, B. et al., 2012. Accuracy assessment of a new Monte Carlo based burnup computer code. Annals of Nuclear Energy 45, 29–36. Los Alamos Monte Carlo Group, 2003. MCNP-A general Monte Carlo N-particle transport code version 5. Los Alamos National Laboratory, vols. I–III. Available from Radiation Safety Information Computational Center at Oak Ridge National Laboratory as CCC-740. MacFarlane, R.E., 1999. NJOY-99 nuclear data processing system.