Abstract
This paper reports the results of sensitivity analysis of the multidimension, multigroup neutron diffusion NODAL3 code for the NEACRP 3D LWR core transient benchmarks (PWR). The code input parameters covered in the sensitivity analysis are the radial and axial node sizes (the number of radial node per fuel assembly and the number of axial layers), heat conduction node size in the fuel pellet and cladding, and the maximum time step. The output parameters considered in this analysis followed the above-mentioned core transient benchmarks, that is, power peak, time of power peak, power, averaged Doppler temperature, maximum fuel centerline temperature, and coolant outlet temperature at the end of simulation (5 s). The sensitivity analysis results showed that the radial node size and maximum time step give a significant effect on the transient parameters, especially the time of power peak, for the HZP and HFP conditions. The number of ring divisions for fuel pellet and cladding gives negligible effect on the transient solutions. For productive work of the PWR transient analysis, based on the present sensitivity analysis results, we recommend NODAL3 users to use radial nodes per assembly, axial layers per assembly, the maximum time step of 10 ms, and 9 and 1 ring divisions for fuel pellet and cladding, respectively.
1. Introduction
The National Nuclear Energy Agency of Indonesia, as a Technical Support Organization (TSO), has long considered the capability of safety evaluation of a nuclear power plant (NPP), especially for pressurized water reactors (PWR), as one of the most important stages within the nation first NPP introduction program. To support the safety evaluation capability, an analytical code, NODAL3, has been developed [1]. The NODAL3 code solves steady-state as well as time-dependent few-group nodal neutron diffusion equations in 3D Cartesian geometry and coupled with a simple thermal-hydraulic model for typical PWRs. Three modules, that is, a steady-state problem, thermal-hydraulics model, and time-dependent solvers, are tailored in the code.
The verification and validation (V&V) works have been conducted and their results showed that the NODAL3 code produces very satisfactory results for the steady-state problems [2]. The code has also been verified for light water reactor (LWR) transient benchmarks. The transient benchmark verification works covered the well-known NEACRP 3D LWR core transient benchmark [3, 4]. The benchmark represents a severe reactivity initiated accident (RIA) due to single or multiple control rods ejection which is one important accident challenging the PWR safety. The benchmarks were selected since many organizations participated using various methods as well as approximations [3–6], so that in addition to the reference solutions (by the PANTHER code), the calculation results of the NODAL3 code can also be compared to other codes’ results. In the transient benchmarks, the initial conditions (critical boron concentrations), time of power peak, power peak, final power, final average Doppler temperature, final maximum fuel centerline temperature, final coolant outlet temperature, and so forth were evaluated. The NODAL3 transient benchmarks’ results also show an excellent agreement with the reference solutions (the PANTHER code’s solutions) for both the quarter- and full-core geometry models [7, 8].
During the above-mentioned static and transient benchmark calculations, we also found that the NODAL3 calculations’ accuracy is affected by some parameters related to the spatial and time variables treatment. This paper reports our sensitivity analysis of the NODAL3 code for the NEACRP 3D LWR core transient benchmarks. The parameters covered in the sensitivity analysis are the spatial (radial and axial) node sizes (number of nodes per fuel assembly and the number of axial layers), time step width on reactor power, and heat conduction node sizes in the fuel pellet and cladding.
Some other published researches on the sensitivity analysis for transient analysis of PWR reactor used a standalone neutronic and thermal-hydraulic codes which were coupled with other interface tools [9–11] while in the NODAL3 code, the neutronic and thermal-hydraulic modules are embedded in the code. Not found in those published works, in this study, the effects of axial zone width and time step width are also carried out. Through this sensitivity study, the parameters in the time-dependent few-group neutron diffusion and thermal-hydraulic equations which are sensitive can be recognized, and the ranges of the parameters, in which reliable and accurate results can be obtained, are quantized.
The paper is organized as follows. In Section 2, a brief description of the NEACRP 3D LWR core transient benchmark (PWR part) and the NODAL3 code solutions are given. In Section 3, sensitivity analysis results are presented and discussed. Section 4 provides the summary of the present work.
2. NEACRP 3D LWR Core Transient Benchmark and Reference Solutions
For the transient benchmark problem, the NEACRP 3D LWR core transient benchmark [3, 4] was selected. The detailed specifications of the benchmark problem are given in [3, 4], and readers should consult the reference. Control rod ejection events, stipulated in the benchmarks, may 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 fast perturbation of the neutronics and thermal-hydraulics core parameters which challenges the accuracy of any nodal code. The transient events in the benchmark are initiated by a rapid ejection of control rod (CR) at HZP (hot zero power, 0 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 [3, 4], there are 6 cases in the benchmark, that is, cases A, B, and C, for both HZP and HFP conditions. The calculation of core geometries and conditions of six cases is summarized as follows: A1: quarter-core geometry and ejection of the central CR at core condition HZP; A2: quarter-core geometry and ejection of the central CR at core condition HFP; B1: quarter-core geometry and ejection of the peripheral CR at core condition HZP; B2: quarter-core geometry and ejection of the peripheral CR at core condition HFP; C1: full-core geometry and ejection of one peripheral CR at core condition HZP; C2: full-core geometry and ejection of one peripheral CR at core condition HFP.
The reactor core configuration for cases A and B is shown in Figure 1, while for case C it is shown in Figure 2. There are 157 fuel assemblies in the core where 49 assemblies are with the control rod (CR) and 64 reflector elements (each has a width of 21.606 cm) surrounding the core. The axial zone of the reactor divided into 16 layers with the size 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 upper and lower axial reflectors (light water) have a thickness of 30 cm.
The reference solution of the above-mentioned benchmark is provided by the PANTHER code using fine temporal and spatial mesh. The reference solution contains several key steady states and transient parameters such as critical boron concentration, power peaking factor, time of power peak, power peak, averaged Doppler temperature at 5 s, maximum fuel centerline temperature at 5 s, and coolant outlet temperature at 5 s. The PANTHER reference solutions were calculated using the radial meshes per fuel assembly, the axial meshes, and the minimum time step of 2 ms [3]. Although there is no explicit endorsement that the PANTHER’s solutions are the most accurate, for the above-mentioned benchmark calculations, participants commonly used the PANTHER’s solutions as the reference solutions. For the present sensitivity analysis, we also adopted the PANTHER’s solutions as the reference solutions but for comparative purpose only. For example, as will be shown later, if we use radial meshes per fuel assembly instead of PHANTER’s radial meshes, our calculation result will deviate from PANTHER’s solutions, but they would probably be more accurate than the one of PANTHER’s. Our benchmark calculations reported in [7] were calculated using the following input data: radial meshes per fuel assembly, axial layers, the maximum time step of 5.5 ms, 9 fuel pellet radial divisions, and 2 fuel cladding radial divisions.
3. Sensitivity Analysis Results
As for the NODAL3 code sensitivity analysis, all calculations were carried out by using adiabatic method option. Among the three models adopted into NODAL3 code, point reactor model, adiabatic model (AM), and improved quasistatic model, the AM model is the optimal model from the point of view of speed and accuracy.
In the NODAL3 code modeling, a PWR fuel assembly can be divided into a node or several nodes in the radial direction, and many layers in the axial direction. For investigating the sensitivity of the number of node in the radial direction, we varied the number of radial nodes per fuel assembly from to (the original value used in [7] is ). As for the axial direction, the original number of the axial layers used in previous benchmark calculations is 18 [7]. We further subdivided the original axial layer from to . The NODAL3 code has a feature of automatic time step adjustment based on a maximum relative error of power, temperature and other parameters of interest. The user only needs to input the maximum time step. The effect of the maximum time step is investigated by varying it from 0.01 ms to 100 ms (the original value used in [7] is 5.5 ms). The number of radial ring divisions for the heat conduction calculation is varied from 5 to 36 for fuel pellet regions (the original value used in [7] is 9) and from 1 to 8 for cladding regions (the original value used in [7] is 2).
Principally, in the present sensitivity analysis, if one parameter is varied the others are fixed to the original values reported in [7] (cf. Section ). The relative differences shown hereafter are between the NODAL3 calculation result and the PANTHER’s reference solutions. However, as mentioned before, a smaller relative difference does not always represent better accuracy. It is also worthily noted that if a calculated parameter asymptotically converged then we could know that the input data is accurate enough and no further refinement is needed.
3.1. Radial Node Size
Sensitivity analysis results in the radial node size for selected cases A1, C1, and C2 are shown in Tables 1–3 (it is not possible to show tables for all cases). In all cases, time of power maximum is sensitive to the radial node size, as can be observed from the time difference in the range of 0 ms–268 ms; however, it converged to 0.530 s. The significant time differences occurred in cases A1 (later 268 ms) and C2 (later 100 ms). Furthermore, for cases with HZP initial condition (A1, B1, and C1), the calculation results with one node per assembly were not accurate enough regarding the time to power maximum and the power peak value.
First, we will discuss in particular cases A1 and C1 (HZP initial condition) since these cases involve a very fast power transient that may affect the calculation accuracy of reactor dynamics (kinetic) part. Figure 3 shows the reactor power peak versus time for case A1. For same radial node size with the PANTHER’s reference (), a good agreement with the reference solution is achieved where the relative difference is −0.68%. This is attributed to the same number of radial meshes per fuel assembly used in the PANTHER’s solutions. From the numerical method point of view, larger numbers of node per assembly (i.e., 4 × 4 and 8 × 8) are expected to provide more accurate results. From Figure 3, we can observe that larger numbers of node per assembly give higher, that is, more conservative, results. If the result of the size with unacceptable deviation is excluded, the power peak parameter for case A1 is sensitive to the radial node size since the relative difference is in the range of −53.05%–12.88%.
For case C1 (Figure 4), the result for the nodes per assembly could not be shown since the calculation did not converge. However, for the node per assembly, the relative difference is smaller than case A1 as can be observed from the relative difference for the power peak (−17.05%) and time of power peak (faster 10 ms). It should be noted that case C1 is more challenging since it involves control rod ejection located at the core periphery which is severer than case A1.
Next, the sensitivity on the fuel pellet and coolant temperatures will be discussed for cases with HFP initial condition, that is, cases A2 and C2. The two cases involve a slow transient, but a large amount of energy is accumulated in the fuel and coolant. The maximum fuel centerline temperature versus time for cases A2 and C2 is shown in Figures 5 and 6. The maximum fuel centerline temperature has a maximum difference of about 0.86%.
3.2. Axial Node Size
The sensitivity analysis results in the axial layer number are shown, for selected cases A2, B2, and C2, in Tables 4–6 and Figures 7–10. The sensitivity analysis results of the axial layer size show a slightly different trend compared with the one of the radial node size in terms of the power peak and the maximum fuel centerline temperature values. The transient solutions are not too sensitive to the range axial layer size investigated here. However, the time of power peak parameter is sensitive to the axial node size in the HFP condition with the maximum relative difference of 100 ms (case C2 with the layer numbers). The other relative differences worthily noted are 88 ms for case A2 (with layer numbers) and 74 ms for case B2 (with layer numbers).
Figures 9 and 10 show that the effect of axial layer size on the averaged Doppler temperature for cases A2 and C2 is not higher. However, the effect is relatively higher for the maximum fuel centerline temperature (−1.98%), as shown in Table 4. For the same case C2, the maximum fuel centerline temperature (0.80%) is slightly lower than the one of radial node number per assembly.
3.3. Maximum Time Step
Sensitivity analysis results in the maximum time step for selected cases A2, C1, and C2 are shown in Tables 7–9. The maximum time step gives impact when its value is greater than 1 ms for the time of power peak for cases with HFP initial condition. For cases with HFP initial condition, the maximum differences, up to 100 ms, occurred in the time step of 100 ms (longest). However, for maximum time step less than 1 ms, the relative differences are negligible. Tables 7–9 show that there is no impact of the maximum time step on the calculated maximum fuel centerline temperature for cases A2 and C2 since the relative differences are the same for all time steps, 0.44% and 0.70%, respectively. The impact of the maximum time step on the time of peak power of case A2 is shown in Figure 11. We found that for the maximum time step of 100 ms, case C1 did not converge since the calculated fuel temperatures exceeded the material properties range.
3.4. Heat Conduction Node in Fuel Pellet
Effect of heat conduction node size (number of the ring divisions, RD-FUEL) in the fuel for some selected cases A1 and C1 is presented in Tables 10 and 11. The pellet node number is varied for 5, 9, 18, and 36 nodes. This parameter is expected to give impact on the averaged Doppler temperature at 5 s and coolant outlet temperature at 5 s. However, as can be seen in Tables 10 and 11, for the maximum fuel centerline temperature, the maximum relative difference is only −2.21% (case A1). The sensitivity of heat conduction node size on the time to power peak is weak as can be observed in cases A1 and C1, where the maximum relative difference is less than 8 ms (cf. Figure 12). There is no specification or description of the number of nodes in the fuel pellet for the PANTHER’s solutions, so it is difficult to assess the differences. It is worthily noted that in the case C1 the ejected CR is located in the core periphery which should challenge the accuracy of the solutions of all codes (including PANTHER and NODAL3).
3.5. Heat Conduction Node Size in Cladding
The effect of heat conduction node size in the fuel cladding (number of the ring divisions, RD-CLAD) was investigated; however we found there was no significant impact of the node size for all transient parameters. Only the result of case A1 is shown in Table 12.
4. Concluding Remarks
The sensitivity analysis covering the spatial parameters (radial and axial node number per assembly, pellet and cladding ring division) and the temporal parameter (i.e., maximum time step) input on the important calculated transient parameters was conducted for multidimension, multigroup nodal neutron diffusion NODAL3 code, by taking the NEACRP LWR core transient problems as the representative benchmark. Qualitatively, the number of the radial nodes is the most sensitive input to accurately obtain the time to power peak for the HZP and HFP initial conditions. Furthermore, the effect of the maximum time step was found to be significant on the power peak and the time to power peak for HZP and HFP initial conditions. The number of ring divisions for fuel pellet and cladding gives negligible effect on the transient solutions. For productive work of PWR transient analysis, based on the present sensitivity analysis results, we recommend NODAL3 users to use radial nodes per assembly, axial layers per assembly, the maximum time step of 10 ms, and 9 and 1 ring divisions for fuel pellet and cladding.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work has been partially supported by the Finance Ministry of Indonesia for fiscal years of 2011–2013.