Abstract
The in-house coupled neutronic and thermal-hydraulic (N/T-H) code of BATAN (National Nuclear Energy Agency of Indonesia), NODAL3, based on the few-group neutron diffusion equation in 3-dimensional geometry using the polynomial nodal method, has been verified with static and transient PWR benchmark cases. This paper reports the verification of NODAL3 code in the NEA-NSC PWR uncontrolled control rods withdrawal at zero power benchmark. The objective of this paper is to determine the accuracy of NODAL3 code in solving the continuously slow and fast reactivity insertions due to single and group of control rod bank withdrawn while the power and temperature increment are limited by the Doppler coefficient. The benchmark is chosen since many organizations participated using various methods and approximations, so the calculation results of NODAL3 can be compared to other codes’ results. The calculated parameters are performed for the steady-state, transient core averaged, and transient hot pellet results. The influence of radial and axial nodes number was investigated for all cases. The results of NODAL3 code are in very good agreement with the reference solutions if the radial and axial nodes number is 2 × 2 and 2 × 18 (total axial layers), respectively.
1. Introduction
An in-house analytical tool for Pressurized Water Reactor (PWR) safety, NODAL3 code, has been developed by National Nuclear Energy Agency of Indonesia (BATAN) [1]. The capability of developing analytical tool for neutronic and thermal-hydraulic parameters of PWR is important to support the self-review of the proposed nuclear power plant in Indonesia. NODAL3 code is a coupled neutronic and thermal-hydraulic code based on the few-group nodal neutron diffusion equations in 3-dimensional (3-D) Cartesian geometry for steady and time-dependent problems. The code has been verified with the PWR static and transient benchmark cases [1–6]. The calculated results of NODAL3 code show very good agreement with the reference solutions [1, 5]. The PWR transient benchmark cases are reactivity initiated accident (RIA) [7, 8]. In those cases, the reactor is scrammed after the power reaches the limit condition; however, in the real operating condition, there is a delay time for the control rods (CRs) to shut down after the neutron flux or power reaches the limit condition.
The important reactivity insertion analysis in PWR beside the RIA is the uncontrolled withdrawal of single or group CRs at zero power. Fraikin and Finnemann proposed the PWR benchmark cases for the uncontrolled safety withdrawal of CRs at zero power accident [9]. In contrast to the RIA benchmark, the reactor will be scrammed after a certain delay time; it means the maximum power occurs before the scram. The uncontrolled withdrawal of CRs accident affects the reactivity and power distribution anomaly that should be postulated in the safety analysis of DBA (design basis accident) [10]. Therefore, the proposed benchmark is very useful to know the accuracy of an analytical tool in solving the PWR postulated initiating event of uncontrolled withdrawal of CRs at zero power. In addition, there were many participated institutions using various methods and approximation [11], including the reference solutions, so the calculation results of NODAL3 code can be compared directly.
In the benchmark of uncontrolled control rod withdrawal at zero power, the accidents are continuously slow and fast reactivity insertion, due to the single and the group control rod bank withdrawal, respectively. The power and temperature increments are limited by the feedback reactivity of Doppler. The objective of this work is to determine the accuracy of NODAL3 code in analyzing the transient parameters of continuous slow and fast reactivity accident due to the uncontrolled control rod withdrawal at zero power which is one important aspect in the DBA analysis of PWR.
This paper is organized as follows. In Section 2, a description of the NEA-NSC uncontrolled control rod withdrawal at zero power benchmark cases is given. Section 3 describes the core calculation method of NODAL3 code. Section 4 presents some results and discussion, followed by the conclusion of this work in Section 5.
2. NEA-NSC PWR 3D/1D Uncontrolled Control Rods Withdrawal at Zero Power Benchmark Cases
The NEA-NSC PWR uncontrolled CR withdrawal benchmark core is almost the same as the rod ejection benchmark except the CR layout [12] and dynamic evolution [9]. There are 157 fuel assemblies (FA) in the benchmark core, where 48 assemblies are with CR. As seen in Figure 1, there are 3 types of enrichment for FAs, 2.1 w/o, 2.6 w/o, and 3.1 w/o, with 4 CRs banks A, B, C, D, and one shutdown bank S. The banks of B, C, D and S (type 1 absorbers) are located in 2.1 w/o FAs and the bank of A (type 2 absorbers) is located in 3.1 w/o FAs [9]. The radial dimension is 21.606 cm × 21.606 cm while the active core height is 367.3 cm with the upper and lower axial reflector thickness of 30.0 cm. Some FAs are with burnable absorber rods (BA), such as 12 BAs, 16 BAs, and 20 BAs in the 2.6 w/o FAs, as well as 12 BAs and 16 BAs in the 3.1 w/o FAs.
Reference [7] shows some specifications of CR; total CR length is 362.159 cm. The full insertion of CR is 37.7 cm from the bottom of the lower reflector, so the height of 401.183 cm is the full withdrawal CR position. The CRs position for fully insertion and withdrawal correspond to 0 and 228 in unit of steps, respectively. The detailed thermal-hydraulic data and the macroscopic cross sections and derivatives are prescribed in [9].
In the uncontrolled CR withdrawal benchmark, there are 4 cases, Cases A, B, C and D; however, only three Cases, which are Case A, Case B, and Case D, are considered in this work. Case C is not considered since it is different from Case B only by the heat transfer coefficient between cladding and water which is set constant of 30,000 W/m2/K [9]. In addition, the reference results showed only small difference between Case B and Case C [12]. In each case, the transient parameters are analyzed until 10 s after reactor trip signal occurred. The calculated benchmark cases are as follows [9]:(i)Case A represents a single CRs bank withdrawal (bank D) and other CRs banks (C, B, A, and S) remain fully withdrawn until scram.(ii)Case B represents a group CRs bank withdrawal (banks B and C), CRs bank A and D remain fully inserted, and other CRs banks S are fully withdrawn until scram.(iii)Case D represents group CRs bank withdrawal (banks A and B), CRs banks C and D remain fully inserted, and other CRs banks S are fully withdrawn until scram.
3. Benchmark Core Calculations
All calculation results of NODAL3 code will be compared with the reference solution, PANTHER code, as described in [12]. Same as the previous works [1, 5], all transient cases are calculated by adiabatic (AM) and improved quasistatic (IQM) methods. The detail numerical of those methods are described in [1].
NODAL3 code used a 3D symmetrical quarter core configuration as described in Figure 1 for the core calculation. The core is divided into 18 layers axially, from bottom to top, 1 layer for lower reflector (30.0 cm thickness), 16 layers for active core (7.7 cm, 11.0 cm, 15.0 cm, 30.0 cm (10 layers), 12.8 cm (2 layers), 8.0 cm), and 1 layer for upper reflector (30.0 cm). The active core configuration of 2nd layer from bottom (7.7 cm) is described in Figure 2, while the remaining layers are configured as shown in Figure 1.
The calculation results of NODAL3 are compared with the reference solutions for [12]:Steady-state results(i)Critical boron concentration, ppm (result B1)(ii)Axial power peak factor, (result B2)(iii)Radial power peak factor, (result B3)(iv)Global (3D) power peak factor, (result B6) Note: values at 6th (result B4) and 13th (result B5) axial layer are not available to be presented since the NODAL3 code cannot give the output for the those active layers.Transient core averaged results(i)Fission power relative to nominal, %, and time of maximum power, s (result C1)(ii)Average coolant outlet temperature (maximum value), , °C (result C3)(iii)Average core Doppler temperature (maximum value), , °C (result C4) Note: NODAL3 code cannot give the output for result C2, coolant heating.Transient hot pellet results(i)Maximum fuel centreline temperature, °C (result D6) Note: NODAL3 code cannot give the output for result D1–D5 and D7. For information, NODAL3 calculates the maximum inner cladding temperature, however the benchmark result D7 is maximum outer cladding temperature.
It is assumed that the CRs are withdrawn with the velocity of 72 steps/minute, while the CRs drop with constant velocity of 228 steps in 2.2 s. The CRs begin to drop 0.6 s after the fission power reaches 35% of nominal power (2,775 MW). It is noted that the initial power is 10−13 of nominal power. The moderator inlet conditions (flow, pressure, temperature, and boron concentration) are constant during the transient [9].
The NODAL3 calculations were carried out by using the following input conditions: 2 × 2 radial nodes per FA, 1 × 18 axial nodes layers and the maximum time step of 5.5 ms. In this work, the influence of the radial and axial node number is investigated. We found in the previous research work that the number of nodes gives a significant influence in the calculation results [5]. In order to be consistent with the previous researches [1, 5], the NODAL3 results are compared to the references in terms of their absolute value of relative deviation (%) [1, 12].
4. Results and Discussion
4.1. Steady-State Results
Table 1 shows the steady-state parameter results (B1, B2, B3, and B6) for Cases A, B, and C. The calculation results of NODAL3 have the maximum deviations of 0.64% in the calculated and parameters for Case A, while the deviation of other calculated parameters is less than 0.37% for all cases. The NODAL3 steady-state benchmarks’ results show an excellent agreement with the reference’s results [12]. It is noted that the radial and axial node number used in NODAL3 code is different with the reference, PANTHER code. Therefore, the effect of the radial and axial node number is presented in Tables 2 and 3.
Tables 2 and 3 show the influence of the radial and axial node number per FA, respectively. It is clear that the radial node number has significant influence compared with the axial one, since the maximum deviation changes in the range of 0.64%–1.27% and 0.32%–0.64% for radial and axial nodes, respectively. Among those parameters, the critical boron concentration and the global (3D) power peak factor, , are sensitive to the radial node number. Those results show that the core configuration of Case A is sensitive to the nodes number, probably because of the local effect of the single CRs bank withdrawal. The calculation results of NODAL3 codes become closer to the reference’s results if it is used with the radial nodes number per FA of 2 × 2 and the axial node numbers of 2 × 18. By using the recommended node numbers, the deviation can be lowered to 50% or the maximum deviation is to be 0.32%. However, a further analysis needs to be carried out in the future to clarify which the node number is optimum to obtain a satisfy calculation results.
4.2. Transient Core Averaged Results
Table 4 and Figures 3–5 show the calculated transient core parameters by using NODAL3 code for Cases A, B, and C. More than 10% deviations occur in the fission power relative to nominal () of Cases B and D, with the range of 12.78%–17.90%, while the deviation range for Case A is < 0.48%. As noted in [12], the power increased rapidly because there is continuous fast reactivity insertion, so those deviations are observed. However, the deviation of NODAL3 is less than the standard deviation of the benchmark participants’ results in previous work [12], because their deviations are in the range of 21%–24.5%.
Table 4 shows that the deviation of time of maximum fission power is about 1.22%–1.44% for all cases. The deviation is equivalent to difference; that is, of NODAL3 occurs about 0.42 s–1.18 s later. Those are clearly shown in Figures 3–5; however, the deviations of NODAL3 are still smaller than the benchmark participants’ results [12].
Tables 5 and 6 show effect of radial and axial nodes number on the (fission power relative to nominal) and time of maximum fission power parameters. It is clear that the calculated results are sensitive with the radial and axial node number, since the deviation is in the range of 0.56%–17.90% and 0.48%–17.90% for radial and axial nodes number, respectively. On the other hand, the calculated time of maximum fission power results is not sensitive since the deviation is in the range of 0.0%–2.59% and 0.0%–1.34% for radial and axial node number, respectively. It is noted that the reference, PANTHER code, used the 3 × 3 nodes and 48 nodes for radial and axial, respectively [12], so the increment of axial nodes in NODAL3 calculations, 18 nodes (1 × 18) to 36 nodes (2 × 18), give closer result to the reference.
4.3. Transient Hot Pellet Results
Table 7 shows that the deviation of >1% occurs in Cases B and D, while the deviation in Case A is < 1%, except IQM method. The maximum deviation of fuel centerline temperature is 3.62% (Case B with AM method) or equivalent to (lower). However, the maximum deviation of NODAL3 is less than the standard deviation of the benchmark participants’ results in previous work [12], because their deviations are in the range of 19°C–30.8°C.
Tables 8 and 9 indicate that the maximum fuel temperature is more sensitive to the axial node number than the radial node number, since the results vary with the case. It is clear in Case B, since the deviation is in the range of 0.73%–3.62% and 3.53%–3.85% for axial and radial node number, respectively. By using some axial node numbers, the deviation can be lowered to 79.8% (2 × 18 nodes), 66.6% (4 × 18 nodes), and 69.3% (8 × 18 nodes).
5. Conclusions and Future Works
The in-house coupled neutronic and thermal-hydraulic code of BATAN, NODAL3, based on the few-group neutron diffusion equation in 3-dimensional geometry using the polynomial nodal method, has been verified with the NEA-NSC PWR uncontrolled control rods withdrawal at zero power benchmark. The results of NODAL3 code show very good agreement with the reference solutions, especially with the 2 × 2 of radial nodes and 2 nodes per axial layer (total 18 layers).
As future works, the NODAL3 code will be used to analyze the real experimental data of PWR such as described in the IAEA-TECDOC-815 [13], as well as transient analysis of an advanced PWR type reactor, such as AP1000 reactor.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work has been partially supported by the IAEA Research Contract no. 18731/R0 and the Finance Ministry of Indonesia (DIPA) for fiscal years of 2015-2016.