The Verification of Coupled Neutronics Thermal-Hydraulics Code NODAL 3 in the PWR Rod Ejection Benchmark

A coupled neutronics thermal-hydraulics code NODAL3 has been developed based on the few-group neutron diffusion equation in 3-dimensional 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.


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 2-dimensional (2D) and 3-dimensional (3D) multigroup neutron diffusion codes for various types of research reactors, namely, the BATAN-2DIFF and BATAN-3DIFF codes, respectively [3][4][5][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 mesh-averaged or mesh-edge approaches for each core-grid 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][8][9][10][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, BATAN-FUEL, has been used for establishing the new equilibrium core of RSG-GAS multipurpose reactor using silicide fuel with higher fuel loading.Currently the code is used for the routine incore fuel management of RSG-GAS reactor [14,15].The BATAN-FUEL code has also several unique features such as the capability to search automatically an equilibrium core without performing lengthy and time-consuming transient cores.
In the last several years, a 3D coupled neutronics and thermal-hydraulic calculation code, MTR-DYN, had been developed for safety analysis of a material testing research reactor (MTR) [16] such as RSG GAS.Several transient characteristics of the RSG-GAS reactor, for example, reactivity insertion accident (RIA), reduced flow rate, and some combination of accident scenarios, had been analyzed by using the MTR-DYN 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).PWR-type 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 finite-difference method adopted in the BATAN-2DIFF, BATAN-3DIFF, BATAN-FUEL, and MTR-DYN 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 time-dependent fewgroup neutron diffusion equations in 3D Cartesian geometry and is coupled with a simple thermal-hydraulic model for PWRs.
The NODAL3 code has been verified with the steady state light water reactor (LWR) benchmarks, such as IAEA-2D, KOERBERG, BIBLIS, and IAEA-3D, 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 well-known 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.

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 thermal-hydraulics model of a typical PWR fuel pin; and the third module is the time-dependent solver for the reactor dynamics.In the first module, the few-group 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 node-averaged 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 transverse-integrated 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 thermal-hydraulic 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 single-phase 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 finite-difference 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 thermal-hydraulics calculations in terms of fuel temperature, moderator (coolant) temperature, and density are fed back via appropriate macroscopic cross-sections.
In the third module, two time-dependent 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 time-dependent group neutron flux which is assumed can be factorized into [25] where () and Ψ  (, ) are the time-dependent amplitude and shape functions, respectively.Under this assumption, the time-dependent few-energy 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; V  : neutron speed (cm s −1 ); Ψ  (, ): timedependent shape function (cm −2 s −1 ); (): time-dependent amplitude function;   (, ): diffusion coefficient in time t (cm); Σ  (, ): macroscopic removal cross-section in time  (cm −1 ); Σ   →  (, ): macroscopic scattering cross section from group   to  in time (cm −1 ); ]   Σ   (, ): number of prompt neutrons emitted per fission times macroscopic fission cross section in time (cm −1 );   : 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  eff 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 time-step adjustment.Different time steps for amplitude function (Δ  ), shape function (Δ  ), and thermal-hydraulic 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 Runge-Kutta 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 thermal-hydraulics 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, thermal-hydraulics calculation can be delayed until needed.Of course, Δ  has to be checked to preserve the stability of the thermal-hydraulics calculation itself.NODAL3 code is also equipped with an automatic algorithm for time-step adjustment based on the rule and users are recommended to use the option.Figure 1: The radial (a) and axial (b) PWR benchmark core configuration [22]. in a PWR.This event is followed by a significant localized perturbation of the neutronics and thermal-hydraulics core parameters.Therefore, the PWR rod ejection benchmark prepared by the OECD/NEA can be used to evaluate the accuracy of a coupled neutronics thermal-hydraulics 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.

PWR Rod Ejection Benchmark
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.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.condition.All calculation results of NODAL3 are compared to the reference calculation results of PANTHER which were published in 1993 (PANTHER-1993) 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:

Results and Discussions
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 PANTHER-1993 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 PANTHER-1997 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 PANTHER-1997 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 thermal-hydraulics 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% (Δ = 1.13 ∘ C)    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  PANTHER-1993, while for the ANCK code was compared to the reference PANTHER-1997.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 Δ = 14.90 ∘ 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.

Conclusions and Future Works
A coupled neutronics thermal-hydraulics code, NODAL3, based on the few-group 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.

Figure 2 :Figure 3 :
Figure 2: The initial control rod positions for case A1 (circled position is ejected).

Figure 4 :
Figure 4: The initial control rod positions for case B1 (circled position is ejected).

Figure 5 :
Figure 5: The initial control rod positions for case B2 (circled position is ejected).

Figure 6 :
Figure 6: Comparison between NODAL3 and references results for case A1.

Figure 7 :
Figure 7: Comparison between NODAL3 and references results for case A2.

and 3 .
06% (Δ = 20.8∘ 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 Δ = 15.6 ∘ C.Moreover,

Figure 8 :
Figure 8: Comparison between NODAL3 and references results for case B1.

Figure 9 :
Figure 9: Comparison between NODAL3 and references results for case B2.
Control rod ejection event can occur as a consequence of the rupture of the control rod drive mechanism (CRDM)

Table 2
shows steady state and transient results of cases A1 and B1 at HZP condition and cases A2 and B2 at HFP

Table 2 :
The calculation results of NODAL3 code for steady sate and transient problems.