Abstract
Maintaining the structural integrity of the reactor pressure vessel (RPV) is a critical concern related to the safe operation of nuclear power plants. To estimate the structural integrity over the designed lifetime and to support analyses for a potential plant life extension, an accurate calculation of the fast neutron fluence ( MeV or MeV) at the RPV is significant. The discrete ordinates method is one of the main methods to solve such problems. During the calculation process, many factors will affect the results. In this paper, the deviations introduced by different differencing schemes and mesh sizes on the AP1000 RPV fast neutron fluence have been studied, which are based on new discrete ordinates code ARES. The analysis shows that the differencing scheme (diamond difference with or without linear zero fix-up, theta weighted, directional theta weighted, and exponential directional weighted) introduces a deviation within 4%. The coarse mesh (4 × 4 cm meshes in plane) leads to approximately 23.7% calculation deviation compared to those of refined mesh (1 × 1 cm meshes in plane). Comprehensive study on the deviation introduced by differencing scheme and mesh size has great significance for reasoned evaluation of RPV fast neutron fluence calculation results.
1. Introduction
Reactor pressure vessel (RPV) is the key equipment which is positioned to withstand the enormous operating pressure in nuclear power plant. From the viewpoint of the safety over the designed lifetime, the integrity of RPV must be guaranteed. In the process of reactor operation, the material properties of RPV will be deteriorated gradually due to being continuously irradiated by high-energy neutrons. Usually, energy ranges of MeV or MeV are considered [1]. Hence, an accurate estimate of RPV fast neutron fluence to evaluate the structural material radiation damage becomes important.
The RPV fast neutron fluence calculation model is complex for its anisotropy scattering and deep penetration. To meet the demands mentioned above and obtain detailed fast neutron flux distribution in RPV, three-dimensional calculation with the discrete ordinates method ( method) is completely feasible [2]. method simulation of the reactor pressure vessel fast neutron fluence obtains a high degree of confidence in the results. However, in the calculation process, deviations introduced by alternative algorithms cannot be neglected. A relatively sufficient analysis and conclusion have been done both at home and abroad. In 1995, Haghighat, who has made great contributions to this area, has evaluated the related effects of various numerical techniques specifically for pressure vessel fluence calculations using DORT [1]. In 2001, the Regulatory Guide 1.190, published by NRC [3], puts forward the demands for a reasonable assessment of the biases caused by different calculation methods.
Recently, the fast development of AP1000 reactor has attracted extensive attention. This paper further studies the AP1000 pressure vessel fast neutron fluence calculation correspondingly. In this paper, fast neutron flux distribution along AP1000 pressure vessel is calculated with multidimensional parallel discrete ordinates transport code ARES. ARES, developed by School of Nuclear Science and Engineering in North China Electric Power University, is primarily applied for the nuclear devices’ shielding calculation and analysis. Sufficient verification and validation for ARES transport code system have been done by experimental benchmark and reference code [4–6]. This research is focused on the deviation caused by different algorithms and mainly discusses the effects of alternative spatial differencing methods and mesh sizes. Errors on account of these two aspects are illustrated by comparing calculation results from different test cases.
This study has modeled the AP1000 reactor and performed a series of fluence calculation cases. Section 2 describes the AP1000 reactor model and explains the methodology used for these calculations. The results are indicated and the effect of different parameters is discussed in Section 3. Section 4 summarizes the findings and conclusions.
2. Pressure Vessel Fast Neutron Fluence Calculation
2.1. Description of AP1000 Model
The AP1000 model for these studies is a 3400 MW (thermal power) pressurized water reactor designed by Westinghouse. The core consists of 157 fuel elements and is surrounded by the core shroud, core barrel, thermal shield, pressure vessel, and biological shield. The overall dimensions of the AP1000 model are 503 × 503 × 485.876 cm, while each assembly is 21.504 × 21.504 cm and each fuel assembly is composed of 17 × 17 array of fuel pins. However, in concrete shield region, neutron flux is difficult to converge for particular calculation case and the absence of whom does not affect the results in pressure vessel. Hence the concrete shield region is removed from the three-dimensional calculation model. The details of the geometry dimensions and material compositions are available in Table 1.
Figure 1 illustrates the geometry configuration modeled with GGTM. GGTM is a program which is included in the BOT3P software package and aimed at describing the geometry, the material, and the fixed neutron source for three-dimensional transport applications [7].
2.2. Methodology
A multigroup cross section library based on ENDF/B-VII is used for the ARES transport calculation. Table 2 lists the energy group boundaries in the library for MeV.
The studied problem is provided for a standard core loading. To describe the model accurately, the three-dimensional calculation of neutron source is significant. For the spatial source distribution, the cycle average assembly radial and axial powers are employed. Based on the contributions of 235U and 239Pu to the fission neutron source, the conversion factor of power to neutron source is calculated. This study takes the average of 235U and 239Pu fission spectra as the source energy spectrum.
approximation is introduced. corresponds to the order of the expansion in Legendre polynomials of the scattering cross section matrix and represents the order of the flux angular discretization. Fully symmetrical quadrature sets are applied. However, for the sake of obtaining precise reactor pressure vessel fluence calculation results under certain situations, Legendre/quadrature order may be needed. In this paper, from the engineering application perspective, can already obtain an accurate and converged result. According to the analysis of calculation results, the calculation precision of is indeed enough to obtain the fully converged reference to comparing with different cases.
2.3. Differencing Scheme and Mesh Dividing
As mentioned above, there are several factors that can lead to deviation in the RPV fast neutron fluence calculation. The discrepancy caused by alternative spatial differencing methods and mesh sizes here mainly examined is part of the overall biases. The detailed explication of differencing scheme employed in ARES is presented as follows.
The spatial differencing scheme applied in ARES includes the diamond difference with or without linear zero fix-up (DZ/DD), theta weighted (TW), directional theta weighted (DTW), and exponential directional weighted (EDW). The traditional diamond differencing scheme (DD), which assumes a linear relationship between the directional flux at the cell center and cell boundaries [8], requires small grid spacing to ensure the accuracy. This contributes to large amount of calculation and is time consuming. When the mesh interval is relatively large, it is easy to yield negative fluxes and these negative fluxes may cause oscillations in the iterative process, even iteration divergent [8]. Usually, “zero correction” approach (DZ) is adopted to eliminate the negative flux but this will lead to nonlinear differential equations, reduce the original second-order computational accuracy of DD, and cause load imbalance in parallel solving process [9]. Hence, to overcome this, the theta weighted (TW) is developed to ensure the nonnegativity of extrapolated fluxes [10]. The value of , depending on the practical issues and experiences, may also have a great impact on the results [1]. In multidimensional calculation, since the discrete direction does not match the spatial mesh, the phenomenon of nonphysical oscillation may appear even with the mesh refinement. Therefore, on the basis of the theta weighted differencing scheme, the directional theta weighted (DTW) and exponential directional weighted (EDW) differencing schemes which have adaptive features and can mitigate flux oscillation to get more accurate scalar flux are introduced [11, 12]. DTW is suitable for the area in which flux changes gently. It can be found that EDW, which assumes an exponential distribution for angular flux within each cell, can obtain higher accuracy in larger grid spacing and is more stable than DTW.
The guideline provided by the Regulatory Guide 1.190 [3] describes the principle of mesh dividing in geometry exactly. It is obvious that mesh spacing will introduce great influence on the calculation results. Hence, this paper studies the effect of mesh size based on Cartesian coordinates in detail and presents corresponding conclusions.
3. Results and Discussion
3.1. Results of AP1000 Model
In this section, the deviation introduced by differencing scheme and mesh size is examined. The objective is to obtain the relative differences between the solutions and further better estimate the fast neutron fluence calculation method of AP1000 pressure vessel. Differencing schemes include DD, DZ, TW, DTW, and EDW. EDW is chosen to be the base case. The uniform mesh interval of plane is 1 cm (original mesh), 2 cm, and 4 cm. Considering the model calculation scale, axial direction is divided into 42 grids referring to the axial power distributions. Figures 2 and 3 show the fast neutron fluence distribution. Figure 2 depicts the three-dimensional distribution of flux and Figure 3 presents the radial profile of flux at 0 deg. azimuth of the core midplane. From the core periphery to the pressure vessel along the radial, the fast neutron fluence attenuation is about four orders of magnitude. Further, the descent is steepest throughout the pressure vessel.
Here focus on the fast neutron fluence in pressure vessel. As illustrated in Figures 4 and 5, the azimuthal flux has no strong attenuation on the order of magnitude within pressure vessel. In addition, the fast neutron fluence ( MeV) distributions on four different positions are roughly similar and neutron energy greater than 1.0 MeV has the same behavior. For the distribution at 1/4T, the flux change trend is substantially alike for MeV and MeV. Therefore, the following statements only analyze the results at 1/4T PV for MeV and this will no longer be repeated next.
3.2. Effect of Differencing Scheme
This paper studies the DD, DZ, TW, ZW, DTW, and EDW differencing schemes individually and compares between them. All of the cases are calculated based on 1 cm × 1 cm radial mesh spacing. As it can be seen in Figure 6, different schemes exhibit similar flux distribution trend. The schemes in Figure 6(a) show somewhat severe flux oscillations, while ZW, DTW, and EDW differencing scheme mitigate the nonphysical oscillation. However, the oscillations are not discussed in-depth here and only the overall behavior is considered.
(a)
(b)
The azimuthal shape of flux ratios to the reference EDW solution is presented in Figure 7. As a result of the oscillation introduced by DD, DZ, and TW, compared to EDW which has smoother results, the corresponding deviations also have oscillation along the azimuth.
As illustrated in Figure 7, DZ and TW () have similar results. Although the flux distribution calculated by DD is close to these two schemes, the oscillation of which is stronger and causes bigger biases. The ZW differencing scheme, which has a higher accuracy in the source region but poor precision in the nonsource region [13], presents roughly the largest results and deviations. In some cases, such as when the neutron flux is relatively flat (central of the core, the source region, etc.), due to the small flux gradient, the truncation error of DTW will also decrease; hence the accuracy of the calculation is improved. For pressure vessel fast neutron fluence calculation, the flux attenuation gradient is not gentle. EDW is more suitable for solving such situation than DTW differencing scheme. Therefore, DTW provides a deviation within ±2% versus reference EDW method. In summary, different differencing schemes will affect the pressure vessel fluence to some extent. Based on the comparisons between EDW and other schemes for AP1000 model, the maximum bias is 4%.
3.3. Different Value Adoptions in TW Scheme
In TW differencing scheme, different adoptions of parameter, which could be any value between 0 and 1, will lead to different calculation results. Currently, DORT and TORT apply as a default value; thus the calculated results are always close to DZ and converge faster. For pressure vessel region, as Figure 8 shows, the increase of the value does not affect the flux distribution too much when is greater than 0.5. However, large value of will bring in obviously flux oscillation compared to small one.
(a)
(b)
Further, is taken as a reference and the relative differences between other values are provided. The comparison indicated in Figure 9 shows that the effect of different theta values on the calculation results can be up to 5.6%. The smaller the value is, the larger the neutron fluence calculation result will be. According to the analysis mentioned above, the reasonable choice of value in TW scheme should be on the basis of actual situation and experience, so as to avoid a great deviation and simultaneously decrease the oscillation as possible.
3.4. Adaptation of Mesh Size
Apart from the differencing scheme, the effect of mesh size is investigated for AP1000 RPV as follows. Figures 6, 7, 10, and 11 depict the azimuthal flux distribution of alternative spatial differencing methods with 1 cm × 1 cm, 2 cm × 2 cm, and 4 cm × 4 cm plane mesh spacing, respectively. For coarse mesh, negative flux appears when applying DD scheme. Moreover, DZ and TW always display similar behavior. Therefore, here we only show the results and biases upon TW, ZW, DTW, and EDW method for various mesh sizes. EDW differencing scheme is taken as the reference solution.
(a)
(b)
(a)
(b)
From the overall trend presented above, EDW, which assumes an exponential distribution for angular flux within each cell, obtains bigger discrepancy compared to other schemes. Furthermore, with the change of mesh size, coarse mesh makes larger deviation between different spatial differencing solutions than refined mesh. The maximum bias can achieve 17.1% difference relative to EDW methods.
This paper also has examined the sensibility of alternative differencing scheme to mesh size. From the perspective of CPU time, coarse mesh spends much less running time than refined mesh. Although CPU time is significantly reducing with the increase of mesh size, the influence of mesh size on the flux calculation results should not be ignored. Figures 12–15 exhibit the effect of mesh size for different differencing schemes. These figures show the azimuthal flux distribution and azimuthal profile of percent differences for 2 cm × 2 cm and 4 cm × 4 cm versus reference 1 cm × 1 cm plane mesh spacing. As shown when spatial grid width varies, TW, ZW, and DTW act out similar changing trend. These solutions are sensible to the size of mesh. Comparing to the original 1 cm mesh, different mesh spacing divisions cause the deviation up to 23.7%. The maximum bias mainly appears in the range of 5 degrees away from the azimuth of zero degree. It may be because the gradient of flux is large in this area; thus different grid step will greatly affect the calculation results. As it has been mentioned before, EDW is more suitable to calculate the problem with big mesh size and can obtain higher accuracy. As illustrated in Figure 15, EDW method gives more stable results than other differencing schemes, and the bias in the case of different mesh sizes is smaller. The maximum deviation is reduced to 8.9%.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
4. Conclusions
Comprehensive analysis of the influence factors is of great significance to evaluate the RPV fast neutron fluence calculation reasonably and ensure the reactor operated safely. Based on the AP1000 PV model and a new discrete ordinates code ARES, this paper has analyzed the effects of differencing scheme, mesh size, and their combinations. From the discussion mentioned above, it is found out that the differencing scheme (DZ, DD, TW, DTW, and EDW) introduces a deviation within 4%, and the discrepancy caused by TW scheme with an alterable value is less than 5.6%. Furthermore, it is observed that the bigger the mesh is, the greater the bias between different differencing schemes will be. As for the impact of variable mesh sizes employed in this calculation, the coarse mesh gives up to 23.7% flux deviation compared to those of the refined mesh.
These conclusions are specific to the particular AP1000 problem studied here. The results demonstrate the magnitude and distribution of the deviation between different differencing schemes and mesh sizes. The consequences and discussions presented here are part of an ongoing effort to access the overall deviation of AP1000 pressure vessel fast neutron calculation. When taking the effect of differencing scheme and mesh size on AP1000 pressure vessel fast neutron fluence calculation into account, the findings above can establish a relatively sufficient guideline.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (11505059 and 11575061) and the Fundamental Research Funds for the Central Universities (2016MS59).