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

Authors: Junxiao Zheng; Bin Zhang; Shengchun Shi; Yixue Chen
Journal: Science and Technology of Nuclear Installations
 (2016)

---

## 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 (E
>
1.0 MeV or E
>
0.1 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 X
Y plane) leads to approximately 23.7% calculation deviation compared to those of refined mesh (1 × 1 cm meshes in X
Y 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.


---

## Body

## 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 ofE
>
1.0 MeV or E
>
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 (S
N method) is completely feasible [2]. S
N 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. Section2 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 Table1.Table 1
AP1000 reactor characteristics.
Geometry
Material
Reactor core
Number of fuel assemblies
157
UO2 enriched to 2.7287%, density 10.42 × 103 kg/m3 Zircaloy-4, Inconel, stainless steel SS-304, and water
Fuel assembly type
17 × 17
Fuel assembly pitch
21.504 cm
Fuel assembly active core height
425.7 cm
Core shroud
DimensionsThickness
See Figure1
3 cm
Stainless steel SS-304
Core barrel
Inner radiusThickness
170.28 cm4.5 cm
Stainless steel SS-304
Neutron pad
Inner radiusThickness
175.78 cm8 cm
Stainless steel SS-304
Pressure vessel cladding
Inner radiusThickness
202.0 cm0.5 cm
Stainless steel SS-304
Pressure vessel
Inner radiusThickness
202.5 cm21.28 cm
Steel A533B, low carbon steelFigure1 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].Figure 1
Horizontal and vertical cross section of the AP1000 reactor.


### 2.2. Methodology

A multigroup cross section library based on ENDF/B-VII is used for the ARES transport calculation. Table2 lists the energy group boundaries in the library for E
>
0.1 MeV.Table 2
Cross section library energy group structure for energy rangeE > 0.1 MeV.
Group number
Upper energy (MeV)
Group number
Upper energy (MeV)
Group number
Upper energy (MeV)
Group number
Upper energy (MeV)
Group number
Upper energy (MeV)
1
19.640
22
7.788
43
2.466
64
0.962
85
0.297
2
17.333
23
7.408
44
2.385
65
0.907
86
0.295
3
16.905
24
7.047
45
2.365
66
0.863
87
0.287
4
16.487
25
6.703
46
2.346
67
0.821
88
0.273
5
15.683
26
6.592
47
2.307
68
0.781
89
0.247
6
14.918
27
6.376
48
2.231
69
0.743
90
0.235
7
14.550
28
6.065
49
2.123
70
0.707
91
0.224
8
14.191
29
5.770
50
2.019
71
0.672
92
0.213
9
13.840
30
5.488
51
1.921
72
0.639
93
0.202
10
13.499
31
5.221
52
1.827
73
0.608
94
0.193
11
12.840
32
4.966
53
1.738
74
0.578
95
0.183
12
12.523
33
4.724
54
1.653
75
0.550
96
0.174
13
12.214
34
4.493
55
1.572
76
0.523
97
0.166
14
11.618
35
4.066
56
1.496
77
0.498
98
0.158
15
11.052
36
3.679
57
1.423
78
0.450
99
0.150
16
10.513
37
3.329
58
1.353
79
0.408
100
0.143
17
10.000
38
3.166
59
1.287
80
0.388
101
0.136
18
9.512
39
3.012
60
1.225
81
0.369
102
0.129
19
9.048
40
2.865
61
1.165
82
0.334
103
0.123
20
8.607
41
2.725
62
1.108
83
0.302
104
0.117
21
8.187
42
2.592
63
1.003
84
0.299
105
0.111The 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 of235U 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.P
3
-
S
8 approximation is introduced. P
3 corresponds to the order of the expansion in Legendre polynomials of the scattering cross section matrix and S
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, P
5
-
S
16 Legendre/quadrature order may be needed. In this paper, from the engineering application perspective, P
3
-
S
8 can already obtain an accurate and converged result. According to the analysis of calculation results, the calculation precision of P
3
-
S
8 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 S
N 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 R
θ
Z 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.


## 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 Table1.Table 1
AP1000 reactor characteristics.
Geometry
Material
Reactor core
Number of fuel assemblies
157
UO2 enriched to 2.7287%, density 10.42 × 103 kg/m3 Zircaloy-4, Inconel, stainless steel SS-304, and water
Fuel assembly type
17 × 17
Fuel assembly pitch
21.504 cm
Fuel assembly active core height
425.7 cm
Core shroud
DimensionsThickness
See Figure1
3 cm
Stainless steel SS-304
Core barrel
Inner radiusThickness
170.28 cm4.5 cm
Stainless steel SS-304
Neutron pad
Inner radiusThickness
175.78 cm8 cm
Stainless steel SS-304
Pressure vessel cladding
Inner radiusThickness
202.0 cm0.5 cm
Stainless steel SS-304
Pressure vessel
Inner radiusThickness
202.5 cm21.28 cm
Steel A533B, low carbon steelFigure1 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].Figure 1
Horizontal and vertical cross section of the AP1000 reactor.


## 2.2. Methodology

A multigroup cross section library based on ENDF/B-VII is used for the ARES transport calculation. Table2 lists the energy group boundaries in the library for E
>
0.1 MeV.Table 2
Cross section library energy group structure for energy rangeE > 0.1 MeV.
Group number
Upper energy (MeV)
Group number
Upper energy (MeV)
Group number
Upper energy (MeV)
Group number
Upper energy (MeV)
Group number
Upper energy (MeV)
1
19.640
22
7.788
43
2.466
64
0.962
85
0.297
2
17.333
23
7.408
44
2.385
65
0.907
86
0.295
3
16.905
24
7.047
45
2.365
66
0.863
87
0.287
4
16.487
25
6.703
46
2.346
67
0.821
88
0.273
5
15.683
26
6.592
47
2.307
68
0.781
89
0.247
6
14.918
27
6.376
48
2.231
69
0.743
90
0.235
7
14.550
28
6.065
49
2.123
70
0.707
91
0.224
8
14.191
29
5.770
50
2.019
71
0.672
92
0.213
9
13.840
30
5.488
51
1.921
72
0.639
93
0.202
10
13.499
31
5.221
52
1.827
73
0.608
94
0.193
11
12.840
32
4.966
53
1.738
74
0.578
95
0.183
12
12.523
33
4.724
54
1.653
75
0.550
96
0.174
13
12.214
34
4.493
55
1.572
76
0.523
97
0.166
14
11.618
35
4.066
56
1.496
77
0.498
98
0.158
15
11.052
36
3.679
57
1.423
78
0.450
99
0.150
16
10.513
37
3.329
58
1.353
79
0.408
100
0.143
17
10.000
38
3.166
59
1.287
80
0.388
101
0.136
18
9.512
39
3.012
60
1.225
81
0.369
102
0.129
19
9.048
40
2.865
61
1.165
82
0.334
103
0.123
20
8.607
41
2.725
62
1.108
83
0.302
104
0.117
21
8.187
42
2.592
63
1.003
84
0.299
105
0.111The 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 of235U 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.P
3
-
S
8 approximation is introduced. P
3 corresponds to the order of the expansion in Legendre polynomials of the scattering cross section matrix and S
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, P
5
-
S
16 Legendre/quadrature order may be needed. In this paper, from the engineering application perspective, P
3
-
S
8 can already obtain an accurate and converged result. According to the analysis of calculation results, the calculation precision of P
3
-
S
8 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 S
N 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 R
θ
Z 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 ofX
Y 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.Figure 2
Three-dimensional distribution of fast neutron fluence (E
>
0.1 MeV).Figure 3
Radial profile of fast neutron fluence (E
>
1.0 MeV and E
>
0.1 MeV, 0 deg. azimuth at the core midplane).Here focus on the fast neutron fluence in pressure vessel. As illustrated in Figures4 and 5, the azimuthal flux has no strong attenuation on the order of magnitude within pressure vessel. In addition, the fast neutron fluence (E
>
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 E
>
0.1 MeV and E
>
1.0 MeV. Therefore, the following statements only analyze the results at 1/4T PV for E
>
0.1 MeV and this will no longer be repeated next.Figure 4
Azimuthal fast neutron fluence (E
>
0.1 MeV) distribution of PV clad, 1/2T, 1/4T, and 3/4T at core midplane.Figure 5
Azimuthal fast neutron fluence distribution at PV midplane 1/4T forE
>
0.1 MeV and E
>
1.0 MeV.


### 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 Figure6, 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.Figure 6
Azimuthal flux distribution of different differencing schemes: (a) DD, DZ, and TW and (b) ZW, DTW, and EDW.
(a)
(b)The azimuthal shape of flux ratios to the reference EDW solution is presented in Figure7. 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.Figure 7
Deviation between different differencing schemes on azimuthal flux distribution.As illustrated in Figure7, 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%.


### 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 θ
=
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.Figure 8
Effect ofθ parameter in TW scheme on azimuthal flux distribution: (a) θ
=
0.1,0.3,0.5 and (b) θ
=
0.7,0.9 and DZ.
(a)
(b)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.Figure 9
Deviation between different theta value adoptions in TW on azimuthal flux distribution.


### 3.4. Adaptation of Mesh Size

Apart from the differencing scheme, the effect of mesh size is investigated for AP1000 RPV as follows. Figures6, 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 X
Y 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.Figure 10
Calculation with 2 cm × 2 cmX
Y plane mesh spacing: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different differencing schemes.
(a)
(b)Figure 11
Calculation with 4 cm × 4 cmX
Y plane mesh spacing: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different differencing schemes.
(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. Figures12–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 X
Y 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%.Figure 12
Effect of mesh size for TW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.
(a)
(b)Figure 13
Effect of mesh size for ZW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.
(a)
(b)Figure 14
Effect of mesh size for DTW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.
(a)
(b)Figure 15
Effect of mesh size for EDW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.
(a)
(b)


## 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 ofX
Y 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.Figure 2
Three-dimensional distribution of fast neutron fluence (E
>
0.1 MeV).Figure 3
Radial profile of fast neutron fluence (E
>
1.0 MeV and E
>
0.1 MeV, 0 deg. azimuth at the core midplane).Here focus on the fast neutron fluence in pressure vessel. As illustrated in Figures4 and 5, the azimuthal flux has no strong attenuation on the order of magnitude within pressure vessel. In addition, the fast neutron fluence (E
>
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 E
>
0.1 MeV and E
>
1.0 MeV. Therefore, the following statements only analyze the results at 1/4T PV for E
>
0.1 MeV and this will no longer be repeated next.Figure 4
Azimuthal fast neutron fluence (E
>
0.1 MeV) distribution of PV clad, 1/2T, 1/4T, and 3/4T at core midplane.Figure 5
Azimuthal fast neutron fluence distribution at PV midplane 1/4T forE
>
0.1 MeV and E
>
1.0 MeV.


## 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 Figure6, 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.Figure 6
Azimuthal flux distribution of different differencing schemes: (a) DD, DZ, and TW and (b) ZW, DTW, and EDW.
(a)
(b)The azimuthal shape of flux ratios to the reference EDW solution is presented in Figure7. 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.Figure 7
Deviation between different differencing schemes on azimuthal flux distribution.As illustrated in Figure7, 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%.


## 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 θ
=
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.Figure 8
Effect ofθ parameter in TW scheme on azimuthal flux distribution: (a) θ
=
0.1,0.3,0.5 and (b) θ
=
0.7,0.9 and DZ.
(a)
(b)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.Figure 9
Deviation between different theta value adoptions in TW on azimuthal flux distribution.


## 3.4. Adaptation of Mesh Size

Apart from the differencing scheme, the effect of mesh size is investigated for AP1000 RPV as follows. Figures6, 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 X
Y 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.Figure 10
Calculation with 2 cm × 2 cmX
Y plane mesh spacing: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different differencing schemes.
(a)
(b)Figure 11
Calculation with 4 cm × 4 cmX
Y plane mesh spacing: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different differencing schemes.
(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. Figures12–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 X
Y 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%.Figure 12
Effect of mesh size for TW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.
(a)
(b)Figure 13
Effect of mesh size for ZW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.
(a)
(b)Figure 14
Effect of mesh size for DTW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.
(a)
(b)Figure 15
Effect of mesh size for EDW differencing scheme: (a) azimuthal flux distribution and (b) azimuthal profile of deviation between different mesh sizes.
(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.



---

*Source: 3461290-2016-09-15.xml*