Acta Polytechnica CTU Proceedings doi:10.14311/AP.2016.4.0062 Acta Polytechnica CTU Proceedings 4:62–67, 2016 © Czech Technical University in Prague, 2016 available online at http://ojs.cvut.cz/ojs/index.php/app NUMERICAL MULTIGROUP TRANSIENT ANALYSIS OF SLAB NUCLEAR REACTOR WITH THERMAL FEEDBACK Filip Osuský∗, Branislav Vrban, Peter Ballo, Štefan Čerba, Jakub Lüley, Gabriel Farkas Institute of Nuclear and Physical Engineering, Faculty of Electrical Engineering and Information Technology, Slovak University of Technology in Bratislava, Bratislava, Slovakia ∗ corresponding author: filip.osusky@stuba.sk Abstract. The paper describes a new numerical code for multigroup transient analyses with a thermal feedback. The code is developed at Institute of Nuclear and Physical Engineering. It is necessary to carefully investigate transient states of fast neutron reactors, due to recriticality issues after accident scenarios. The code solves numerical diffusion equation for 1D problem with possible neutron source incorporation. Crank-Nicholson numerical method is used for the transient states. The investigated cases are describing behavior of PWR fuel assembly inside the spent fuel pool and with the incorporated neutron source for better illustration of the thermal feedback. Keywords: diffusion equation, neutronics, transient state, Crank-Nicholson method, thermal feedback. 1. Introduction Nuclear reactors are most reliable source of energy at present time. The advantage of new fast neutron reac- tor concepts is the possibility of the transmutation of the 238U into the fissionable 239Pu, so the worldwide resources of fissionable materials could be enlarged when the fast neutron reactors will be commercialized. However, it is necessary to carefully investigate the transient states of these new concepts, because the re- criticality after accident scenario may be achieved [1]. For this purpose, a new numerical code has been developed within the frame of neutronic calculation group at the Institute of Nuclear and Physical En- gineering (INPE). The code has been used as the basic step in understanding of transient state and fast reactor R&D and we refer to it as TRANSOS. The TRANSOS code is able to solve a time depen- dent neutron flux distribution of slab systems with respect to the thermal feedback in multigroup ap- proximation. The power method with infinity norm (Manhattan norm) is used for a calculation of station- ary state [2]. The transient state is calculated by the Crank-Nicholson numerical method [3]. It is possible to define the geometry in the form of coarse mesh intervals with material parameters assigned to each of the mesh interval. Point and linear neutron sources can be incorporated to the design with different flu- ency for each investigated energy group of neutrons. The boundary conditions have to be specified in con- junction with design parameters such as temperature of the system and power generation rate. Initial state of neutron flux should be given for successful run. The developed code can be used for academic purposes at the faculty and for investigation of simple systems. This paper follows the previous work [4] and next sections describe calculation methodology of time de- pendent and in-dependent solution, accuracy of the methodology and analysis of two nuclear systems with different 235U enrichment. Higher impact of the ther- mal feedback can be observed on microscopic cross section properties for low energy neutron interaction with matter, so the thermal systems are investigated to show this effect. 2. Calculation Methodology 2.1. Stationary state (time in-dependent solution) The code solves common neutron diffusion equation with internal sources for the stationary state given by [5] Dg∆φg − Σa,gφg − ∑ h Σg→hφg+ + ∑ h Σh→gφh + χg ∑ h νΣf,hφh = 0. (1) Lower index g in all equations represents the neu- tron energy group and D stands for the associated diffusion coefficient. The first term represents the leakage of neutrons from the system, the second term represents absorptions (radiation capture), the third represents scattering from the group to another group h, the fourth represents contributions of scattering from another group and the last term describes the contribution of fission source to specified energy group (χg is the proportion of neutrons emitted by fission in group g, ν is an average yield of neutrons from the fission). The Box Scheme was used for the dis- cretization of examined area and therefore (1) can be written in a numerical form (2) for two energy groups 62 http://dx.doi.org/10.14311/AP.2016.4.0062 http://ojs.cvut.cz/ojs/index.php/app vol. 4/2016 Numerical Analysis of Slab Nuclear Reactor with Thermal Feedback of neutrons [2, 4], −D̃k−1g φ k−1 g + ( D̃k−1g +D̃ k g +Σ k rghk ) φk−D̃kgφ k+1 g = = λχkghk ( 2∑ g′=1 νΣkfg′φ k g′ ) + hk ( 2∑ g′=1 g′ 6=g Σks,g′gφ k g′ ) , (2) where Σr represents the removal macroscopic cross section (3), k represents spatial node, hk is the length of the spatial interval, λ stands for the fundamental eigenvalue, Σf is the fission macroscopic cross section, Σs is the scattering macroscopic cross section between energy groups of neutrons, and D̃ arbitrary coupling coefficient (3) [2]. Σkrg = Σ k ag + 2∑ g′=1 g′ 6=g Σks,gg′ D̃ k k = 2 Dkg hk Dk+1g hk+1 Dkg hk + D k+1 g hk+1 (3) 2.2. Transient state (time dependent solution) 2.2.1. Crank-Nicholson numerical method The Crank-Nicholson method is often applied to the diffusion problems and achieves good numerical sta- bility for the transient states. The finite difference method is used as for the stationary state. This method is able to solve heat equation and similar partial differential equations [3] M ∂2U(x,t) ∂x2 = ∂U(x,t) ∂t , (4) where M represents material parameters and U is so- lution of the differential equation. Applying the finite difference method, the central node of the common diffusion equation can be rewritten numerically for the transient state as − M̃k−1ht 2hk Uk−1,j − M̃kht 2hk Uk+1,j + ( 1 + M̃k−1ht 2hk + M̃kht 2hk ) Uk,j = M̃k−1ht 2hk Uk−1,j+1 + M̃kht 2hk Uk+1,j+1 + ( 1 − M̃k−1ht 2hk − M̃kht 2hk ) Uk,j+1 + f(hk), (5) where the index k represents spatial node, index j represents step time ( U·,j stands for initial condition of the function at present time step, U·,j+1 is initial condition of the function at next time step), the hk length of the spatial interval, ht length of the time in- terval, f(hk) is increment of the function U (particular solution). The variables and unknowns from (5) can be rewritten in the form (6) for the neutron diffusion equation (1) and (2). − D̃k−1g ht 2hk φk−1,jg − D̃kght 2hk φk+1,jg + ( 1 + D̃k−1g ht 2hk + D̃kght 2hk ) φk,jg = D̃k−1g ht 2hk φk−1,j+1g + D̃kght 2hk φk+1,j+1g + ( 1 − D̃k−1g ht 2hk − D̃kght 2hk ) φk,j+1g + χkght ( 2∑ g′=1 νΣkfg′φ k,j g′ ) − Σkrghtφ k,j g + ht ( 2∑ g′=1 g′ 6=q Σks,g′gφ k,j g′ ) + Sght, (6) where the last term S represents the external neutron source. The heat diffusion equation K ρc ∂2T(x,t) ∂x2 + q̇ = ∂T(x,t) ∂t (7) can be similarly rewritten to [6] − K̃k−1ht 2hk Tk−1,j − K̃kht 2hk Tk+1,j + ( 1 + K̃k−1ht 2hk + K̃kht 2hk ) Tk,j = K̃k−1ht 2hk Tk−1,j+ + K̃kht 2hk Tk+1,j+1 + ( 1 − K̃k−1ht 2hk − K̃kht 2hk ) Tk,j+1 + ( 2∑ g′=1 Σkfg′φ k,j g′ ) E ht ρkck e, (8) where K̃k = 2 Kk kkρkck Kk+1 kk+1ρk+1ck+1 Kk kkρkck + Kk+1 kk+1ρk+1ck+1 , and where K represents material conductivity, ρ is material density, c stands for material heat capacity, q̇ is power generation, T stands for temperature, E is an average released energy per fission in eV and parameter e is an electron charge. Accuracy of the above mentioned method is discussed in the next subsection. 2.2.2. Accuracy of the Crank-Nicholson method The accuracy of the method for the initial and bound- ary conditions (more in Section 3) is presented in Fig. 1. The analytical solution of the differential equa- tion for a homogenous slab system is sinus, so the problem can be approximated by function f(x) f(x) = a sin(bx) + c. 63 Filip Osuský, Branislav Vrban, Peter Ballo et al. Acta Polytechnica CTU Proceedings Figure 1. Accuracy of Crank-Nicholson method. In this case, the whole system is homogeneous sub- critical nuclear system. It is necessary to insert the external neutron source because the neutron flux will never stabilize on non-zero value without it. The external neutron source is linear, located between 217–227 mm, with fluency 1.0 × 106 s−1 in thermal energy spectrum. Deviations between the modeled problem and analytical solution are presented in Ta- ble 1. Parameter Mean value Dispersion Rel. error a 111 935 126.1 0.1127 % b 7.076 99 × 10−3 2.015 × 10−6 0.028 47 % c −785.941 97.17 12.36 % Table 1. Deviation between numerical and analytical solution. Amplitude (a) and shape parameter (b) are in accor- dance with the analytical solution (Table 1). Deviance in location parameter (c) is caused by the external neutron source and by zero flux boundary condition. 3. Geometry of the Nuclear Slab System The investigated geometry represents middle cut of one heterogenic VVER440 fuel assembly that is placed in the middle. One edge of the investigated geometry is represented by homogenous representation of the heterogeneous fuel assembly. Control rod assembly is placed on the opposite edge. This case was intro- duced due to simulation of the transient state during operational conditions. Geometry is also simplified into one dimension and some lengths are changed for faster convergence of numerical calculation (Fig. 2). Red material repre- sents UO2 (dimension d1 is 7 mm and is the same for every red and orange material). Orange material is the representation of UO2 with 3.35 % of 153Gd. The two cases of the enrichment are investigated (Table 2). Green is B4C (dimension d5 is 142 mm), Figure 2. Graphical representation of investigated nuclear reactor core. Enrichment Material Case 1 Case 2 UO2 5 % 22 % UO2 + 153Gd 4.4 % 21.4 % Table 2. Enrichment of UO2. blue stands for H2O (dimension d2 is 5 mm, d3 = 19 mm, d4 = 10 mm). Purple is a mixture of the whole area between x–y points (d5 is 142 mm). All cross sections are obtained from ENDF/B-VII.1 library for the temperatures of 296.3 K, 600 K and 1800 K, for the discrete energies of 0.0253 eV and 2 MeV [7]. Isotropic scattering is assumed and the diffusion coefficient can be calculated according to (9) [8], where Σt stands for the total macroscopic cross section, Σtr is the transport macroscopic cross section, and µ̄0 represents unsymmetrical parameter of neutron scattering (if 64 vol. 4/2016 Numerical Analysis of Slab Nuclear Reactor with Thermal Feedback isotropic scattering is applied µ̄0 = 0). In this case, a molecular cross section of H2O is not considered, even if the molecular cross section is higher than calculated cross sections of particular nuclides [8]. D ' λtr 1 3 = 1 3Σtr = 1 3 ( Σt − µ̄0Σs ) ≈ 1 3Σt . (9) More details about geometry input can be found in [4]. It is necessary to note that this is the only the- oretical approximation and that it is mainly because the values of the macroscopic scattering cross section were simplified with certain degree of uncertainty. 4. Initial and Boundary Conditions The zero incoming current boundary condition is used in all cases of the neutron diffusion α = J φ = 0.5, where J is the neutron current and φ is the neutron flux. This can be interpreted as that the environment behind the boundary of the nuclear system is vacuum. According to (2), the stationary state is represented by proportional neutron flux distribution that is normed to maximal value 1. It is necessary to set the initial condition in the form of total power generation P = ∫ V E ( 2∑ g′=1 Σfg′φg′ ) e dV. (10) The absolute value of the neutron flux distribution is calculated from the proportional neutron flux dis- tribution (2) and from the initial power generation condition (10). The power generation rate is calcu- lated for the case 1 from the typical power generation of VVER reactor (1375 MWth) and it is approximately 1.576 kW/m per the assembly P [W/m] = Pth ASSEMBLY NR · HEIGHT = = 1375 × 106 349 · 2.5 = 1576, where ASSMEBLY NR is the total number of fuel as- semblies in the reactor and HEIGHT represents height of the core. The power generation rate is the same for the case 2 with higher enrichment to lower impact of the thermal feedback to ensure the nuclear system’s subcritical state. Temperature is kept at constant 296.3 K at the boundaries that correspond to Dirich- let boundary condition. Temperature of the whole system is also set to 296.3 K before the transient for all cases. No thermo-mechanical physics is applied (the densities of the particular materials are constant). Also no fluid mechanics is applied (the heat diffusion equation solves just thermal conduction equation) and this simplification results in total temperature increase of the whole nuclear system. 5. Results 5.1. Case 1 keff = 0.282 601 is calculated from the stationary state of TRNASOS code for this scenario. This condition may be achieved in spent fuel pool. The worst case is simulated by the neutron flux corresponding to the operational state before the transient (the power generation is 1.576 kW/m). It is possible to observe (Fig. 3) that the thermal energy spectrum of neutrons decreases rapidly due to deep subcritical state of the nuclear system. The same evolution can be seen in case of fast energy spectrum of neutrons in Fig. 4. The temperature is set to 296.3 K at the beginning of the transient state. The power is generated by (10) according to the neutron flux and the temperature evolution is shown in Fig. 5. The temperature in- creases up to 460 K and according to the temperature change, the microscopic cross sections are modified. The whole system is below melting point temperature of the fuel. Figure 3. Neutron flux distribution for energy group 0.0253 eV in Case 1 (the neutron flux decreases with time evolution due to deep subcritical state). Figure 4. Neutron flux distribution for energy group 2 MeV in Case 1 (the neutron flux decreases with time evolution due to deep subcritical state). 65 Filip Osuský, Branislav Vrban, Peter Ballo et al. Acta Polytechnica CTU Proceedings Figure 5. Temperature evolution in Case 1 with keff = 0.282601 (the temperature decreases with time evolution – the nuclear system is brought into safe- conditions). 5.2. Case 2 The keff is equal to 0.982 693 in this case, due to higher enrichment of 235U (Table 2). The neutron source is incorporated within the core and its total emission of neutrons is 1.1 × 106 s−1 in the area be- tween 217–227 mm on x-axis. It corresponds to the concept of super-thermal liquid-helium (4He) source (UCN) [9] and portion of 1.0 × 106 s−1 is in the ther- mal energy spectrum. Subcritical multiplication can be derived from N∞ = S 1 1 −keff = 1.1 × 106 1 1 − 0.982 693 ≈ 108. The result of subcritical multiplication formula is consistent with the calculated results that are shown in Fig. 6 and Fig. 7. The whole system starts from zero power level. The transient state of the neutron flux is stabilized approximately after 50 000 s. However, the transient of the temperature evolution is represented by higher inertia and the transient state is stabilized in the region of 140 000 s (Fig. 8). The temperature does not exceed 380 K. Figure 6. Neutron flux distribution for energy group 0.0253 eV in Case 2 (nuclear system is stabilized after subcritical multiplication for keff = 0.982 693). Figure 7. Neutron flux distribution for energy group 2 MeV in Case 2 (nuclear system is stabilized after subcritical multiplication for keff = 0.982 693). Figure 8. Temperature evolution in Case 2 (temper- ature is stabilized below melting point of fuel material for keff = 0.982 693). 6. Conclusions The successful application of the TRANSOS code was demonstrated. It has been found that the maximum value of the neutron flux and its temperature is located in the middle of the system within the heterogeneous region. On the other hand, the diffusion theory is more appropriate for the homogenous problems and calculation error increases with greater heterogeneity. In all cases, the temperature does not exceed melting point of fuel material. The code that was written in C++ can be used for academic purposes at the institute and for solution of simple nuclear systems. It has been confirmed that the diffusion theory is suitable for fast calculations and it is used during development of theoretical reactor design. Therefore, development of simple codes increases the knowledge level of neutron physics and it can result in enhancing of nuclear safety. The contribution of 238U thermal feedback with thermo-hydraulic properties of the system has to be investigated in the future. Acknowledgements This work was financially supported by grant of Science and Technology Assistance Agency No. APVV-0123-12 and STU Grant scheme for Support of Young Researchers. 66 vol. 4/2016 Numerical Analysis of Slab Nuclear Reactor with Thermal Feedback References [1] H. Ninokata, T. Sawada, H. Tomozoe, et al. A study on recriticality characteristics of fast reactors in pursuit of recriticality-accident-free concepts. Progress in Nuclear Energy 29:387–393, 1995. doi:10.1016/0149-1970(95)00067-T. [2] J. Han-gyu. Solution of one-dimensional, one-group neutron diffusion equation, Lecture Note 1. Seoul National University, Reactor Physics Laboratory, 2008. [3] P. B. Patil, U. P. Verma. Numerical Computational Methods. Alpha Science International, Ltd, 2006. [4] F. Osuský, B. Vrban, Š. Čerba, et al. Two-group numerical analysis of one dimensional table nuclear reactor. In Proceedings of the 20th International Conference on Applied Physics of Condensed Matter (APCOM2014), pp. 320–323. 2014. [5] P. Reuss. Neutron physics. EDP Sciences, 2008. [6] M. Kalousek, B. Hučko. Prenos tepla. Vydavateľstvo STU, Bratislava, 1996. [7] Nuclear Data Services. POINT2015 Data. Temperature dependent ENDF/B-VII.1 cross section data, https://www-nds.iaea.org/point/. [8] J. R. Lamarsh. Introduction to Nuclear Reactor Theory. Amer Nuclear Society, 2002. [9] E. Lychagin, V. Mityukhlyaev, A. Muzychka, et al. UCN sources at external beams of thermal neutrons. an example of PIK reactor. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 823:47–55, 2016. doi:10.1016/j.nima.2016.04.008. 67 http://dx.doi.org/10.1016/0149-1970(95)00067-T https://www-nds.iaea.org/point/ http://dx.doi.org/10.1016/j.nima.2016.04.008 Acta Polytechnica CTU Proceedings 4:62–67, 2016 1 Introduction 2 Calculation Methodology 2.1 Stationary state (time in-dependent solution) 2.2 Transient state (time dependent solution) 2.2.1 Crank-Nicholson numerical method 2.2.2 Accuracy of the Crank-Nicholson method 3 Geometry of the Nuclear Slab System 4 Initial and Boundary Conditions 5 Results 5.1 Case 1 5.2 Case 2 6 Conclusions Acknowledgements References