Research Article  Open Access
The Verification of Coupled Neutronics ThermalHydraulics Code NODAL3 in the PWR Rod Ejection Benchmark
Abstract
A coupled neutronics thermalhydraulics code NODAL3 has been developed based on the fewgroup neutron diffusion equation in 3dimensional geometry for typical PWR static and transient analyses. The spatial variables are treated by using a polynomial nodal method while for the neutron dynamic solver the adiabatic and improved quasistatic methods are adopted. In this paper we report the benchmark calculation results of the code against the OECD/NEA CRP PWR rod ejection cases. The objective of this work is to determine the accuracy of NODAL3 code in analysing the reactivity initiated accident due to the control rod ejection. The NEACRP PWR rod ejection cases are chosen since many organizations participated in the NEA project using various methods as well as approximations, so that, in addition to the reference solutions, the calculation results of NODAL3 code can also be compared to other codes’ results. The transient parameters to be verified are time of power peak, power peak, final power, final average Doppler temperature, maximum fuel temperature, and final coolant temperature. The results of NODAL3 code agree well with the PHANTHER reference solutions in 1993 and 1997 (revised). Comparison with other validated codes, DYN3D/R and ANCK, shows also a satisfactory agreement.
1. Introduction
The National Nuclear Energy Agency of Indonesia (BATAN) has been operating three research reactors with the nominal thermal power of 100 kW, 2 MW, and 30 MW, respectively, for nuclear science, technology, and engineering research and development (R&D). The oldest 100 kW reactor has been operated since 1965. To the present date there is no nuclear power plant (NPP) due to strong dependency of the national primary energy on the fossil fuels [1]. However, recently, the evaluation of NPP reactor safety is becoming an important R&D activity at the agency since the nuclear option has been included into the policy of the national energy mix up to the fiscal year of 2025 [2]. One of the highly prioritized R&D activities in BATAN during the fiscal years of 2010 to 2014 is the capability to develop analytical tools for incore fuel management and transient analysis of typical NPPs.
The development of analytical tools in the agency was initiated 2 decades ago by the development of the 2dimensional (2D) and 3dimensional (3D) multigroup neutron diffusion codes for various types of research reactors, namely, the BATAN2DIFF and BATAN3DIFF codes, respectively [3–6]. There are several special features of those codes that do not exist in other similar codes, such as the capability of estimating the radial and axial power peaking factors based on the meshaveraged or meshedge approaches for each coregrid position. The satisfactory results of verification and validation of the codes have been achieved through several calculation benchmarks as well as against the experimental results using a critical assembly and a research reactor [7–11]. Furthermore, an incore fuel management code for research reactors has been developed based on those computer codes [12, 13]. The incore fuel management code, BATANFUEL, has been used for establishing the new equilibrium core of RSGGAS multipurpose reactor using silicide fuel with higher fuel loading. Currently the code is used for the routine incore fuel management of RSGGAS reactor [14, 15]. The BATANFUEL code has also several unique features such as the capability to search automatically an equilibrium core without performing lengthy and timeconsuming transient cores.
In the last several years, a 3D coupled neutronics and thermalhydraulic calculation code, MTRDYN, had been developed for safety analysis of a material testing research reactor (MTR) [16] such as RSG GAS. Several transient characteristics of the RSGGAS reactor, for example, reactivity insertion accident (RIA), reduced flow rate, and some combination of accident scenarios, had been analyzed by using the MTRDYN code [16, 17].
Based on the experiences in developing the static and transient calculation codes for research reactors, the development of the incore fuel management and 3D transient analysis codes is carried out for typical pressurized water reactors (PWR). PWRtype reactor is chosen based on the guidance of the national research programs in 2010–2014 [18]. These research programs are expected to assist the design evaluation of several PWR candidates which are expected to be offered and introduced by many international vendors.
As is well known, the PWR core dimension is considerably much larger compared to one of the research reactors, so that the neutron diffusion problem in PWRs is commonly solved by modern nodal methods [19]. Therefore, besides the finitedifference method adopted in the BATAN2DIFF, BATAN3DIFF, BATANFUEL, and MTRDYN codes, a nodal method has to be introduced for the spatial treatment. The nodal method adopted here is based on the polynomial nodal method proposed by Finnemann et al. [20]. The NODAL3 code, developed in this study, solves steady state as well as timedependent fewgroup neutron diffusion equations in 3D Cartesian geometry and is coupled with a simple thermalhydraulic model for PWRs.
The NODAL3 code has been verified with the steady state light water reactor (LWR) benchmarks, such as IAEA2D, KOERBERG, BIBLIS, and IAEA3D, and very satisfactory results were obtained [21]. However, the code has not been verified for transient LWR benchmark cases. In this paper, we will report the benchmark verification results of the code against the wellknown transient benchmarks of the NEACRP PWR rod ejection cases [22]. The objective of this work is to determine the accuracy of NODAL3 code in analysing the reactivity initiated accident (RIA) due to the control rod ejection which is one important aspect of PWR safety. The NEACRP PWR rod ejection cases are chosen since many organizations participated in the NEA project using various methods as well as approximations [22], so that, in addition to the reference solutions, the calculation results of NODAL3 code can also be compared to other codes’ results.
This paper is organized as follows. In Section 2, a brief numerical model in the NODAL3 code emphasized in the reactor dynamic solver is described. Section 3 presents the benchmarks and the selected cases for this work, followed by the results and discussion presented in Section 4. The last section, Section 5, presents the conclusions from this work and future activities that will be presented.
2. NODAL3 Code Brief Description
NODAL3 code consists of three modules; the first module deals with the nodal equation for the steady state problems; the second module deals with the thermalhydraulics model of a typical PWR fuel pin; and the third module is the timedependent solver for the reactor dynamics. In the first module, the fewgroup neutron diffusion equation in 3D Cartesian geometry is discretized spatially using the polynomial nodal method (PNM). A coarse mesh finite difference (CFMD) formulation is used to determine the nodeaveraged neutron fluxes and the eigenvalue, while the PNM method is used to estimate the accurate coupling between adjacent nodes in the core. Quadratic polynomial expansion for the transverseintegrated flux is adopted [23]. The detailed description of the PNM implementation in the NODAL3 code can be found in [21, 24].
In the second module, that is, the thermalhydraulic module, the heat conduction problem in the fuel rods is discretized in time and space using the conventional finitedifference method. Heat conduction is considered only in the radial direction. Fluid dynamic of the cooling water is modelled under a singlephase flow condition. The mass flow rate in each cooling channel is assumed to be known and specified by the code user. As a result, only the mass continuity and energy conservation equations are to be solved. These are discretized in space and time using finitedifference method and implicit scheme, respectively [24]. An appropriate steam table covering the operational temperature and pressure of PWR and BWR is included in the NODAL3 code [24]. The thermalhydraulics calculations in terms of fuel temperature, moderator (coolant) temperature, and density are fed back via appropriate macroscopic crosssections.
In the third module, two timedependent reactor dynamics models are available, that is, the adiabatic method (AM) and improved quasistatic methods (IQSM). These two methods are selected since they have a high accuracy [25]. This section describes briefly the application of these methods. The spatial and timedependent group neutron flux which is assumed can be factorized into [25] where and are the timedependent amplitude and shape functions, respectively. Under this assumption, the timedependent fewenergy group neutron diffusion equation is commonly written as follows: where prompt neutron energy group index; : total number of prompt neutron energy groups; : delayed neutron energy group index; : total number of delayed neutron energy groups; : prompt neutron index; : delayed neutron index; : neutron speed (cm s^{−1}); timedependent shape function (cm^{−2} s^{−1}); timedependent amplitude function; diffusion coefficient in time t (cm); macroscopic removal crosssection in time (cm^{−1}); macroscopic scattering cross section from group to in time ; number of prompt neutrons emitted per fission times macroscopic fission cross section in time ; fission spectrum of prompt neutron; fission spectrum of delayed neutron; decay constant of precursors (s^{−1}); : concentration of delayed neutron precursors in time (cm^{−3}).
In the AM, firstly, the difference between the neutron spectra of delayed neutrons and the ones of prompt neutrons is neglected. In other words, the delayed neutrons from their precursors are assumed to be born at the same time with the prompt neutrons. Secondly, all time derivatives of the amplitude and shape functions are neglected. Thus, (1) is simplified as where in (3) is the effective multiplication factor after the perturbation occurred. The usual eigenvalue criticality procedure can be readily applied to solve the shape function. The obtained again is used to calculate the new cross sections and other parameters.
In the IQSM, the time derivative of the shape function is approximated with the backward finite difference scheme: so that (1) is approximated as follows: Obviously IQSM is expected to give a better accuracy than AM. For both methods, the most time consuming part is the shape function calculations.
In addition to the three modules, what is not less important is the timestep adjustment. Different time steps for amplitude function (), shape function (), and thermalhydraulic calculation () are adopted in NODAL3 code in order to obtain accurate results with minimal computation times. In some cases, another time step for reactivity calculation () is required for simulating the movement of the control rod. There is no definitive rule to determine the optimal time steps and best relations among them. However, we usually adopt the following rule [24] to determine them manually. The should be taken small enough to maintain the accuracy and stability since the stability of the amplitude function is a necessary condition for the stability of other part of the calculations. In NODAL3 code, the amplitude function is solved using a fourth order explicit RungeKutta method with the typical time steps in the order of 1.0 × 10^{−5} seconds. On the other hand, and are determined through observation of the shape function transient rate. In the practical use of the AM, the temperature change gives a relatively smaller effect on the shape function compared to the composition change, such as a rod withdrawal. Thus, while composition changes occur, smaller and should be adopted. However, the existence of the solutions of the criticality problem in (3) is independent of . During the time intervals, where no composition change occurs, should not be taken less than or , and shape function calculation should wait until thermalhydraulics calculation is done. The is strongly dependent on the accumulation rate of energy (time integrated reactor power) or the transient rate of amplitude function. When the temperature increase/decrease produces no significant change in the reactivity, thermalhydraulics calculation can be delayed until needed. Of course, has to be checked to preserve the stability of the thermalhydraulics calculation itself. NODAL3 code is also equipped with an automatic algorithm for timestep adjustment based on the rule and users are recommended to use the option.
3. PWR Rod Ejection Benchmark
Control rod ejection event can occur as a consequence of the rupture of the control rod drive mechanism (CRDM) in a PWR. This event is followed by a significant localized perturbation of the neutronics and thermalhydraulics core parameters. Therefore, the PWR rod ejection benchmark prepared by the OECD/NEA can be used to evaluate the accuracy of a coupled neutronics thermalhydraulics code in analysing the transient characteristic of a PWR [22].
The transient events in the benchmark are initiated by a rapid ejection of control rod (CR) at HZP (hot zero power, 2775 W) and HFP (hot full power, 2775 MW) conditions. The core configuration and operational data, such as geometry and neutron cross sections, are derived from a real PWR. To allow the problem of a single rod ejection, a CR is added in the center of the core. As shown in [22], there are 6 cases in the benchmark, that is, cases A, B, and C, for both HZP and HFP conditions. In this work, only 4 cases are selected, namely, cases A and B (HZP and HFP), since the initial all control rod positions in the problem case C (HFP) are quite similar to case A and case B (HFP), although the position of the ejected CR in case A or B and case C is not the same.
There are 157 fuel assemblies, where 49 assemblies are with CR, with a radial dimension of 21.606 cm 21.606 cm per assembly. The axial zone of the reactor with the total height of 427.3 cm is divided into 18 layers of different thickness with the following configuration:(i)2 layers for the upper and lower axial reflector with thickness of 30 cm each;(ii)16 layers, from bottom to top, with heights of 7.7 cm, 11.0 cm, 15.0 cm, 30.0 cm (10 layers), 12.8 cm (2 layers), and 8.0 cm.The height of active core and the absorber length is 367.3 cm and 362.159 cm, respectively. The lower edge CR from the bottom of lower reflector is 37.7 cm (0 in unit of step) for fully inserted and 401.183 cm (228 in unit of step) for fully withdrawn, respectively [22]. Figure 1 describes the radial and axial core configuration.
(a)
(b)
The four selected cases are described in Table 1 and Figures 2, 3, 4, and 5. In this work, the benchmark cores are modelled in a symmetrical quarter core geometry with 2 × 2 nodes for a fuel assembly, radially, and 1 (one) node for each layer axially. Some mandatory assumptions in thermalhydraulic method are made in the benchmark specification, such as conductance which is constant and rod expansion and cross flow effects are not considered. Therefore, the thermal relations and properties in the NODAL3 code, such as heat conductivity and specific heat capacity, were identical with the thermal relations of the benchmark data.

4. Results and Discussions
Table 2 shows steady state and transient results of cases A1 and B1 at HZP condition and cases A2 and B2 at HFP condition. All calculation results of NODAL3 are compared to the reference calculation results of PANTHER which were published in 1993 (PANTHER1993) and revised in 1997 (PANTHER 1997) [24, 25]. In this paper, the NODAL3 results are compared to both references in terms of their relative deviation (%) as follows: The results of the steady state condition for the critical boron concentration using NODAL3 show the maximum deviations of 0.85% and 0.42% compared with the published and revised PANTHER solutions, respectively. It is noted that the deviation of 0.85% is equivalent to 10.1 ppm of boron concentration difference. For the critical boron concentration, the results of NODAL3 code are in a very good agreement with the revised reference results. Therefore, the feedback model handling the cross sections by their derivatives is correctly treated in the NODAL3 code [26].

The behaviour of reactor power and average Doppler temperature are shown in Figures 6 to 9. If the PANTHER1993 is used as the reference solution, the maximum deviation of 11.89% (using AM method) occurs in the calculated power peak of case B1, while, if the PANTHER1997 is used, the maximum deviation of 26.00% occurs in the calculated time of power peak of case B2. The deviation of 26.00% is equivalent to Δt = 0.026 s difference; that is, the power peak of NODAL3 occurs nearly 0.026 s later. The deviation of the calculated power peak of case B1 increased to 17.67% if it is compared to the PANTHER1997 result. On the other hand, compared with both references, the deviation of calculated final power (at 5 s) shows a good agreement by the maximum of 3.55%. Since there is no systematic difference, the numerical methods for treating the coupled neutronics thermalhydraulics in the NODAL3 code are considered correct. Probably these benchmark cases are sensitive, especially, concerning the control rod reactivity [27]. Therefore, sensitivity analysis needs to be done in the future to know clearly the cause of the deviations.
The maximum deviations of the fuel temperature parameters obtained by NODAL3 are lower compared to the final power parameters, since they are only 0.53% (°C) and 3.06% (°C) for the final average Doppler temperature and the maximum fuel temperature, respectively. For the final coolant outlet temperature, the maximum deviation of NODAL3 is 5.32% or equivalent to °C. Moreover, Table 2 and Figures 6, 7, 8, and 9 show that the AM and QSM methods are very similar for all cases and all transient parameters with the maximum deviations of 4.13%.
A comparison with the codes that have been validated for the same benchmark, DYN3D/R and ANCK codes, has been carried out as shown in Tables 3 and 4, respectively [25, 28]. It should be noted that the maximum deviation (%) for the DYN3D/R code was compared only to the PANTHER1993, while for the ANCK code was compared to the reference PANTHER1997. It is clearly shown in Table 3 that the maximum deviations of NODAL3 code are lower compared to DYN3D/R code, except for the final coolant outlet temperature parameter with 5.08% higher. It is not significant, since this 5.08% is equivalent to °C.


Furthermore, compared to the ANCK code, the maximum deviations of NODAL3 code are higher in the range of 1.02%–15.74%, except for the critical boron concentration parameter which is lower by 0.32%. The relatively large differences occur in two parameters, the time of power peak and the power peak, with the maximum deviations, that is, 15.74% and 12.50%, respectively. The difference of 15.74% is equivalent to Δt = 0.031 s, while the difference of 12.50% occurs in the reactor power of HZP condition (case B1). However, in the HFP condition of case B2, both codes have a difference of 0.85% for the reactor power parameter. We can conclude that the NODAL3 benchmark calculation results are very satisfactory since they have no significant difference from the validated codes of DYN3D/R and ANCK.
5. Conclusions and Future Works
A coupled neutronics thermalhydraulics code, NODAL3, based on the fewgroup neutron diffusion equation in 3dimensional geometry using the polynomial nodal method, has been verified with the OECD/NEACRP PWR rod ejection benchmark. The results of NODAL3 code show very good agreement with the PHANTHER reference solutions in 1993 and 1997 (revised).
As future works, the sensitivity analysis is needed to be carried out to analyse the cause of relatively higher deviation for the power peak parameters. In addition, the code will be verified to the PWR benchmark of uncontrolled withdrawal control rod to elaborate the accuracy in fast reactivity insertion.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work has been partially supported by the Finance Ministry of Indonesia for fiscal years of 2010–2012. The authors express a special gratitude to Dr. Kunihiko Nabeshima (JAEA) for providing an important reference for this paper.
References
 International Energy Agency, World Energy Outlook 2010, OECD/IEA, Paris, France, 2010.
 Republic of Indonesia, “Presidential Decree no. 5/2006 on National Energy Policy,” 2006. View at: Google Scholar
 P. H. Liem, “Development and verification of BATAN’s standard twodimensional multigroup neutron diffusion code (BATAN2DIFF),” Atom Indonesia, vol. 20, no. 2, pp. 1–19, 1994. View at: Google Scholar
 P. H. Liem, “Development of BATAN’s standard 3D multigroup diffusion code (BATAN3DIFF),” in Proceedings of the 5th Workshop of Computation in Nuclear Science and Technology, Jakarta, Indonesia, 1995. View at: Google Scholar
 P. H. Liem, “BATANADJOINT2D: BATAN's standard code for calculations of integral kinetic parameters of nuclear reactors,” Atom Indonesia, vol. 21, no. 2, pp. 1–18, 1995. View at: Google Scholar
 P. H. Liem, “Development and verification of BATAN's standard diffusion code modules for treating triangular mesh geometry,” Atom Indonesia, vol. 24, no. 2, pp. 75–93, 1998. View at: Google Scholar
 T. M. Sembiring and P. H. Liem, “Validation of BATAN's standard diffusion codes on IAEA benchmark static calculations,” Atom Indonesia, vol. 23, no. 2, pp. 73–91, 1997. View at: Google Scholar
 P. H. Liem and T. M. Sembiring, “Validation of BATAN's standard neutron diffusion codes for control rod worth analysis,” Atom Indonesia, vol. 23, no. 2, pp. 55–72, 1997. View at: Google Scholar
 Zuhair, T. M. Sembiring, and P. H. Liem, “BATAN2DIFF and3DIFF diffusion codes validation on Kyoto University Critical Assembly (KUCA),” Atom Indonesia, vol. 24, no. 1, pp. 15–32, 1998. View at: Google Scholar
 P. H. Liem, “Validation of BATAN standard 3D diffusion code, BATAN3DIFF, on first core of RSGGAS,” Atom Indonesia, vol. 25, no. 1, pp. 47–64, 1999. View at: Google Scholar
 T. M. Sembiring and P. H. Liem, “Validation of BATAN3DIFF code on 3D model of the IAEA 10 MWth benchmark core for partiallyinserted control rods,” Atom Indonesia, vol. 25, no. 2, pp. 91–100, 1999. View at: Google Scholar
 P. H. Liem, “BATANFUEL: a general incore fuel management code,” Atom Indonesia, vol. 22, no. 2, pp. 67–80, 1996. View at: Google Scholar
 P. H. Liem, “Development of an incore fuel management code for searching equilibrium core in 2D reactor geometry (BATANEQUIL2D),” Atom Indonesia, vol. 23, no. 1, pp. 1–19, 1997. View at: Google Scholar
 P. H. Liem, B. Arbie, T. M. Sembiring, P. Prayoto, and R. Nabbi, “Fuel management strategy for the new equilibrium silicide core design of RSG GAS (MPR30),” Nuclear Engineering and Design, vol. 180, no. 3, pp. 207–219, 1998. View at: Google Scholar
 T. M. Sembiring, Tukiran, S. Pinem, and Febrianto, “Neutronic design of mixed oxidesilicide cores for the core conversion of RSGGAS reactor,” Atom Indonesia, vol. 27, no. 2, pp. 85–101, 2001. View at: Google Scholar
 S. Pinem and T. M. Sembiring, “Application of neutronics modelling on the transient simulation of RSGGAS reactor,” in Proceedings of the International Conference on Mathematics and Natural Science, pp. 1056–1059, Institut Teknologi Bandung, Bandung, Indonesia, 2006. View at: Google Scholar
 S. Pinem and T. M. Sembiring, “Transient analysis of RSGGAS reactor core for coolant flow reduction using MTRDYN code,” Journal of Nuclear Reactor Technology, vol. 11, no. 3, pp. 153–1161, 2009. View at: Google Scholar
 Republic of Indonesia, “Presidential Regulation no. 5/2010 on National MediumTerm Development Plan 2010–2014”. View at: Google Scholar
 R. J. J. Stamm’ler and M. Abbate, Reactor Physics in Nuclear Design, Academic Press, London, UK, 1983.
 H. Finnemann, F. Bennewitz, and M. R. Wagner, “Interface current techniques for multidimensional reactor calculations,” Atomkernenergie, vol. 30, no. 2, pp. 123–128, 1977. View at: Google Scholar
 T. M. Sembiring and S. Pinem, “The validation of the NODAL3 code for static cases of the PWR benchmark core,” Journal of Nuclear Science and Technology, vol. 15, no. 2, pp. 82–92, 2012. View at: Google Scholar
 H. Finnemann and A. Galati, “NEACRP 3D LWR core transient benchmark—final specifications,” NEACRPL335 (Revision 1), 1992. View at: Google Scholar
 P. J. Turinsky, “A fewgroup neutron diffusion equation solver utilizing the nodal expansion method for eigenvalue, adjoint, fix source steadystate and transient problems,” Idaho National Engineering Laboratory EGGNRE11406, 1994. View at: Google Scholar
 “NODAL3: A Nodal Neutron Diffusion Code, version 2,” User’s Guide, 2010. View at: Google Scholar
 K. O. Ott and R. J. Neuhold, Introductory Nuclear Reactor Dynamic, American Nuclear Society, Cook County, Ill, USA, 1985.
 S. Aoki, T. Suemura, J. Ogawa, and T. Takeda, “The verification of 3 dimensional nodal kinetics code ANCK using transient benchmark problems,” Journal of Nuclear Science and Technology, vol. 44, no. 6, pp. 862–868, 2007. View at: Publisher Site  Google Scholar
 U. Grundmann and U. Rohde, “Verification of the code DYN3D/R with the help of international benchmarks,” Forschungszentrum Rossendorf FZR195, Research Center Rossendorf Inc., Dresden, Germany, 1997. View at: Google Scholar
 H. Finnemann, H. Bauer, A. Galati, and R. Martinelli, “Results of LWR core transient benchmark,” Tech. Rep., NEA/NSC/DOC(93)25, 1993. View at: Google Scholar
Copyright
Copyright © 2014 Surian Pinem et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.