Analysis for the Effect of Spatial Discretization Method on AP1000 Reactor Pressure Vessel Fluence Calculation

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 (E > 1.0MeV or E > 0.1MeV) 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 cmmeshes inXY plane) leads to approximately 23.7% calculation deviation compared to those of refinedmesh (1 × 1 cmmeshes inXY plane). Comprehensive study on the deviation introduced by differencing scheme andmesh size has great significance for reasoned evaluation of RPV fast neutron fluence calculation results.


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  > 1.0 MeV or  > 0.1 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 2 Science and Technology of Nuclear Installations  code [4][5][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.

Pressure Vessel Fast Neutron
Fluence Calculation 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].

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  > 0.1 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 235 U and 239 Pu to the fission neutron source, the conversion factor of power to neutron source is calculated.This study takes the average of 235 U and 239 Pu fission spectra as the source energy spectrum. 3 - 8 approximation is introduced. 3 corresponds to the order of the expansion in Legendre polynomials of the scattering cross section matrix and  8 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,  5 - 16 Legendre/quadrature order may be needed.In this paper, from the engineering application perspective,  3 - 8 can already obtain an accurate and converged result.According to the analysis of calculation results, the calculation precision of  3 - 8 is indeed enough to obtain the fully converged reference to comparing with different cases.

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.

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 ( > 0.1 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  > 0.1 MeV and  > 1.0 MeV.Therefore, the following statements only analyze the results at 1/4T PV for  > 0.1 MeV and this will no longer be repeated next.

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.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 ( = 0.9) 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%.

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  = 0.9 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.Further,  = 0.9 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.

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.
From the overall trend presented above, EDW, which assumes an exponential distribution for angular flux within each cell, obtains bigger discrepancy compared to other Azi mu th (de g.) 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.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%.

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.

Figure 1 :
Figure 1: Horizontal and vertical cross section of the AP1000 reactor.

Figure 7 :
Figure 7: Deviation between different differencing schemes on azimuthal flux distribution.

Figure 9 :Figure 10 :
Figure 9: Deviation between different theta value adoptions in TW on azimuthal flux distribution.

Figure 11 :Figure 12 :
Figure 11: Calculation with 4 cm × 4 cm  plane mesh spacing: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different differencing schemes.

Figure 13 :Figure 14 :
Figure 13: Effect of mesh size for ZW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.

Figure 15 :
Figure 15: Effect of mesh size for EDW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.

Table 2 :
Cross section library energy group structure for energy range  > 0.1 MeV.