BSS and spectral unfolding methods

The BSS, manufactured by Centronic Ltd., UK, consists of an SP9 ³He thermal neutron detector and a set of high-density polyethylene spheres (0.938 ± 0.009 g/m³) with the following 12 different diameters: 3, 3.5, 4, 4.5, 5, 6, 7, 8, 9, 9.5, 10, and 12 inches (1 inch=2.54 cm). The internal diameter of SP9 was 3.2 cm, the stainless-steel shell was 1 mm thick, and the nominal pressures were 230 kPa ³He and 120 kPa Kr at 20 ℃. The incident neutrons moderated by polyethylene spheres were detected by SP9 at the center of each sphere and a coupled series of suitable electronic equipment [12,13].

The neutron spectrum unfolding method of the BSS is based on the theory of the development of few channel data, which is an underdetermined system of equations [2,14] that has infinite solutions. When the response function R_i(E) and each detector count rate C_i are known, the neutron fluence rate $\phi \left(E\right)$ can be determined using Equation (1): $C_i=\int R_i\left(E\right)\phi \left(E\right)\text{d}E,i=1,2,\dots ,m,$ (1) where i denotes the ith sphere. Equation (1) is discretized in the unfolding procedure, and can be transformed into Equation (2) as follows: $C_i={\displaystyle \sum _{j=1}^n{R}_{ij}{\phi }_j}, i=1,2,3,\mathrm{...},m,$ (2) where $R_{ij}$ is the response of the ith sphere in the jth energy interval, which can be calculated by a simulation and validated by a sufficient monoenergetic neutron field, and ${\phi }_j$ is the fluence rate in the jth energy interval.

The maximum entropy, iteration, and GA methods were used to unfold the neutron energy spectra and provide mutual verification. The maximum entropy and iteration methods achieved by FORTRAN integrated into UMG, called MAXED and GRAVEL, respectively, were employed to obtain the neutron energy spectrum. A priori default spectrum for the preliminary determination of the structure of the neutron energy spectrum was required when using the maximum entropy and iteration methods. The two methods aim to change the default spectrum into one that fits the data but remains “as close as possible” to the default spectrum. Therefore, the a priori default spectrum significantly affects the final results of the spectral unfolding process. Detailed descriptions and analyses of these two methods are available in previous studies [15,16]. For MAXED, the chi-squared per degree of freedom was 1.8, and the highest energy, temperature, and temp reductions were set to 20, 1, and 0.85, respectively. For GRAVEL, the chi-square per degree of freedom was 1.8, and the highest energy and maximum number of iterations were set to 20 and 2000, respectively. The selection of the default spectrum used in this study is described in Sect. 2.3.

The GA, which imitates the Darwinian evolution paradigm, also known as the “survival of the fittest” strategy, was applied to solve the problem of unfolding the neutron energy spectrum [17]. It is a new type of unfolding method that has recently emerged and can eliminate the limitation of requiring a default energy spectrum during calculations. In this study, the GA code was achieved using C++ within the ROOT framework. The main parameters used in the code were as follows: the upper and lower boundaries of the spectral unfolding zone were 20 MeV and 0, respectively, genetic algebra was 500, population size was 30000, crossover probability was 0.3, mutation probability was 0.001, and genetic probability was 0.2.

Model of response functions calculation

The response functions, which are the key parameters of a BSS, were calculated using Geant4, a Monte Carlo code library used for simulating particle transport, coupled with ENDF/B-VII [18]. Figure 2 shows the geometry of a 4.5 inches polyethylene sphere with SP9. The incident neutrons were assumed to be a monoenergetic parallel beam emitted from a circular plane with a diameter equal to that of the simulated sphere, and the circular plane was placed at a distance of 20 cm from the SP9 center. The structure, dimensions, material, and material density of the SP9 and polyethylene spheres were the same as the actual parameters described above. S(α, β) was applied to the neutron transport calculation in polyethylene to solve the problem of neutron thermalization [19]. Excluding the polyethylene sphere and SP9, all the spaces were set to a vacuum. A neutron was considered to escape, and the trace was stopped once it occurred in the vacuum zone during the simulation. The calculated response energy ranged from 10^-9 eV to 20 MeV, considering 20 energy blocks per order of magnitude from 10^-9 eV to 10 MeV, and five energy blocks from 10 MeV to 20 MeV. In addition, the responses of six experimental energies, 144 keV, 250 keV, 565 keV, 1.2 MeV, 2.5 MeV, and 14.8 MeV, which were determined from an accelerator-driven neutron source were also calculated.

Fig. 2Full-screen View

(Color online) Sketch of a Bonner sphere with incident neutrons

To validate the simulated BSS energy response values, the six aforementioned energy responses were calibrated in a monoenergetic neutron field using the shadow cone method, and the thermal neutron response of SP9 was calibrated in a thermal neutron reference irradiation field formed by a neutron source and graphite moderator [20] at PTB. We should also consider that the air gap between the SP9 holder and sphere is a key factor for BSS [21].

Neutron spectral analysis of reactor building

Neutron energy spectra have different structures for incident neutrons with various energies because the absorption cross-section of the moderator varies with the neutron energy. Neutrons above 20 MeV are produced by a proton-induced reaction on a lithium target in the laboratory; its characteristics can be described by a quantum molecular dynamics model [22]. However, in fission reactors or other common neutron fields, the neutron energy is lower than 20 MeV and the energy spectra can be expressed by the following three equations:

In the thermal region (<0.5 eV), the energy spectrum follows the Maxwellian distribution as follows (Eq. (3)) [23]: ${\phi }_1=\frac{2\pi }{{\left(\pi kT\right)}^{3/2}}\sqrt{\frac{2}{m}}E{e}^{-\frac{E}{kT}},$ (3) where ${\phi }_1$ is the thermal neutron fluence rate in cm^-2·s^-1, E is the neutron energy in eV, k is Boltzmann's constant (8.617 × 10 ^{− 5} eV/K), T is the temperature, and m is the neutron mass.

In the intermediate energy region (0.5 eV-0.1 MeV), the classical 1/E distribution was abandoned because the function diverges for E→0 and decreases very slowly for an energy of ~1 keV; the distribution was replaced by Eq. (4) as follows [24]: ${\phi }_2=A\left(1-e^{-\frac{E^2}{E_{\text{d}}^2}}\right)E^{{\alpha }_1-1}{e}^{-\frac{E}{{\beta }_1}},$ (4) where ${\phi }_2$ is the intermediate neutron fluence rate, E_d = 0.0707 eV limits the function at a lower energy, and α₁ (-0.5 < α₁ < 0.5) and β₁ (0 < β₁ < 1) are the parameters in the intermediate energy region.

Moreover, for the fission energy spectrum of U-235 (fast neutrons, 1-20 MeV), the distribution is approximated by the Watt distribution (Eq. (5)) [3,25], which can be explained by the nuclear evaporation model as follows [26]: ${\phi }_3=B{E}^{{\alpha }_2}{e}^{-E/{\beta }_2},$ (5) where ${\phi }_3$ is the fission neutron fluence rate and α₂ (0 < α₂ <1) and β₂ (1 < β₂ < 2) are the parameters of the fast component of the spectrum.

Normally, fission neutrons emitted from U-235 are moderated by both the coolant and moderator, such as light water and heavy water, and then become thermal neutrons and a small number of intermediate energy neutrons, epithermal neutrons, and fission neutrons [27].

According to the physical process, a hypothesis indicating that the spectrum of a reactor building can be explained as a superposition of these energy spectra with different parameters and weight factors can be proposed as follows (Eq.(6)): $\phi =C_1{\phi }_1+C_2{\phi }_3+C_3{\phi }_3,$ (6) where C₁, C₂, and C₃ are different weight factors that are positive or null.

Neutrons in the thermal energy region account for a high proportion of the final spectra, and the most probable energy of a fast neutron may move to a lower energy (<1 MeV) owing to the surrounding moderators. In addition, the neutron-absorbing material also affects the final superposition energy spectrum around the measurement locations. In addition, moderators with a high resonance absorption cross-section distort the energy spectrum in the intermediate energy region, although this effect can be neglected owing to the low energy resolution of the BSS in the intermediate energy region.

According to the analysis of the neutron energy spectrum of reactor buildings, a previously published [28] neutron spectrum at the workplace of a nuclear power plant (Fig. 3), with a covering thermal, intermediate energy, and fast neutrons, was chosen as the default spectrum for the MAXED and GRAVEL methods.

Fig. 3Full-screen View

Default neutron energy spectrum of reactor buildings

On-site measurement

The neutron energy spectra were measured at the following six locations: F1-A (305), F1-B (access route), F2-C (1R326), F2-D (1R511), F3-E (position between pumps 1 and 3 in 501), and F3-F (501 ladder). F1-A and F1-B are located in factory 1, F2-C and F2-D in factory 2, and F3-E and F3-F in factory 3. Table 1 summarizes the measurement devices used. Extra polyethylene holders should fill the holes of the Bonner spheres as much as possible (Fig. 4(a)) because the air gap between the holder and sphere has a significant influence on the response function. The direction of SP9 was parallel to the floor during the measurement (Fig. 4(b)). The height of each sphere support varied to make the geometric center of SP9 with different moderators 35 cm above the ground. According to the uncertainty formula of the linear model (Eq. (7)), the time and total count of a single measurement should be greater than 200 s and 13000, respectively, to ensure that the relative uncertainty of the count rate is less than 1%. $u_{\text{c}}^2={\left(\frac{1}{\sqrt{N}}\right)}^2+u_{\text{t}}^2$ (7) Here, $u_{\text{c}}$ and $u_{\text{t}}$ are the relative uncertainties of the count rate and measurement time, respectively, and N is statistical fluctuation.

Fig. 4Full-screen View

(Color online) (a) Photograph of an 8 inches polyethylene sphere and one holder of SP9. (b) Photograph of the on-site measurement setup

In addition, a neutron dose equivalent ratemeter, LB6411, a ³He detector located at the center of a 25 cm diameter polyethylene moderator, was employed to directly measure the neutron dose equivalent rate $\dot{H^{*}}$(10). The LB6411 is one of the best-known neutron dose equivalent ratemeters, and its measurement results were representative [29]. The calibration factor was 0.998, which was calibrated in a neutron reference radiation field based on a neutron reference with a ²⁴¹Am-Be source. The final results are the average of 10 measurement values.

The on-site process of measurement relied on the following three assumptions: (1) The fluence rate of the neutron radiation field is homogeneous in the local area of measurement, and the minor positional changes when replacing different spheres are negligible. (2) The fluence rate of the neutron radiation field was stable and unchanged during the measurement. (3) The disturbance induced by BSS can be ignored in the measurement area of the neutron radiation field.

Response functions and uncertainty

The response functions that were simulated using Geant4 are shown in Fig. 5. The thermal neutron response of SP9 was 3.2 cm² based on the calibration results in PTB, and the real gas pressure of ³He in SP9 was 248.7 kPa, which is 7.5% higher than the nominal value. Therefore, in subsequent simulations, the actual value of the gas pressure was used. Figure 6 shows the deviation of the simulation compared to the experiment, which indicates the D-values between the ratio of the simulation to the experimental value and 1; excluding the 3 inches and 3.5 inches spheres at a neutron energy of 14.8 MeV, all the results were within 0.9~1.1 at each energy point, indicating that they sufficiently coincided with one another. The largest relative deviation was observed at a high neutron energy (14.8 MeV) for the smallest sphere (3 inches), and the maximum was 14%. Table 2 lists the combined uncertainty (k=1) and components, including the ³He pressure, density of the polyethylene sphere, air gap, isotropy of the sphere, statistic of uncertainty by simulation, maximum deviation, and experimental uncertainty. The uncertainty of the ³He pressure was fixed at 1% for all the spheres and was mainly due to the absorption of the thermal energy by stainless steel. The uncertainty of the density of polyethylene caused by the polyethylene nonuniformity and measurement error ranged from 5% to 13%, which increased as the diameter of the sphere increased. A constant air gap between the SP9 holder and the sphere of 0.4 mm was used in the response function simulations, and the estimations of the deduced uncertainties caused by this factor were based on a 0.4 mm variation. The largest uncertainty (approximately 12 %) occurred for the 3 inches sphere responses. The isotropy of the sphere uncertainty and statistical uncertainty of the simulation changed monotonically, decreasing from 3% to 1% and increasing from 1% to 4%, respectively, when the sphere diameter varied from 3 inches to 12 inches The largest experimental uncertainty of each sphere was relatively stable, which was approximately 4-5% during the calibration based on a monoenergetic neutron field. For each sphere, the maximum deviation between the experimental and simulation results was also selected among the six energy points as an uncertainty component of the response function, which are representative because they cover the thermal neutrons, intermediate energy neutrons, and fast neutrons. Finally, the maximum combined uncertainty, namely 20% (k = 1), was considered as the final response function uncertainty owing to the inability to obtain the experimental values of all the energy points.

Fig. 5Full-screen View

(Color online) Response functions of the Bonner sphere spectrophotometer simulated by Geant4

Fig. 6Full-screen View

(Color online) Plot of the deviation of the simulation from the experimental results versus sphere diameter at different neutron energies (the two red dashed lines represent the region with a deviation of less than 10%)

Normalization

Figure 7 displays the normalized measurement results of the BSS relative to the average count rate with polyethylene spheres and a bare detector at the six locations. The normalization constants at the six locations were 3616, 729, 979, 343, 1354, and 1871 cps, respectively, which were the average count rates of each location. The peak package of F1-A was approximately 3 inches, which was lower than the other locations, indicating that F1-A would have the lowest neutron energy among these positions according to the response function shown in Fig. 5. The peak locations of F1-B and F3-E corresponded to a larger sphere, approximately 4-5 in, thus the neutron energy could be higher than that of the other locations. The peak locations were essentially the same and were located at 3.5 inches or 4 inches for F2-C, F2-D, and F3-F. Therefore, a preliminary conclusion indicates that the neutron energy at each location would be in the following order: F1-B and F3-E > F2-C, F2-D, and F3-F > F1-A.

Fig. 7Full-screen View

(Color online) Normalization of count rate versus sphere diameter at different measurement locations

Unfolding of energy spectra

Figure 8 shows the neutron energy spectra at six positions obtained by MAXED and GRAVEL; the two methods had nearly the same distribution. The thermal and fast neutron peak positions occurred at approximately 0.01 eV and 100 keV, respectively. Although the peak height of the GA unfolded spectrum was different from that of the other methods at certain locations, their corresponding peak positions were essentially in unity in the entire energy region, and the most probable energies of the fast neutrons were also near 10-100 keV. We also observed that the spectrum given by GA at F1-B and F3-E did not have a small thermal neutron peak compared to MAXED and GRAVEL at the same position. This may be caused by the constraint of the default spectrum required when using MAXED and GRAVEL, which presents an apparent peak for thermal neutrons. In addition, Fig. 9 shows the proportion of neutrons with different energy regions over the total neutrons at the six locations. Compared to MAXED and GRAVEL, there were larger proportions of intermediate-energy neutrons and smaller proportions of thermal neutrons and fast neutrons given by GA. This may be due to the influence of the default spectrum when using MAXED and GRAVEL.

Fig. 8Full-screen View

(Color online) Comparison of spectra given by the MAXED, GRAVEL, and GA unfolding methods at (a) F1-A, (b) F1-B, (c) F2-C, (d) F2-D, (e) F3-E, and (f) F3-F

Fig. 9Full-screen View

(Color online) Bar plot of the proportion of fast, intermediate, and thermal neutrons over the total neutrons at six measurement positions (MX, GV, and GA represent MAXED, GRAVEL, and the genetic algorithm, respectively, and the red arrows indicate the borderlines of the measurement positions)

According to the shapes and structures of the unfolded energy spectra (Figs. 8 and 9), three groups were identified; this grouping was consistent with the pre-analysis by the count rate. For group 1, including the spectra obtained at F2-C, F2-D, and F3-F, which was in good agreement with the aforementioned analysis in Section 2.3, most of the neutrons were concentrated in the thermal region, and there was an apparent spike in the fast neutron energy region. For F2-D, the energy of the thermal neutron peak given by GA (1 eV) was higher than that of MAXED and GRAVEL (0.1 eV). Thermal, intermediate energy, and fast neutrons accounted for approximately 40%, 50%, and 10% of the total, respectively. For group 2, including the spectra obtained at F1-A, nearly all the neutrons were located in the thermal (approximately 50%) and intermediate energy (approximately 50%) regions, and the proportion of fast neutrons was near 0%, as shown in Fig. 9. The mean energy given by the three methods was 0.3 eV, indicating that a thicker material with a higher neutron moderation efficiency may be present in the environment. For group 3, including the neutron spectra obtained at F1-B and F3-E, the intermediate energy and fast neutrons were the vast majority (> 80 %), whereas the thermal neutron proportions were only 10-20%. Their corresponding average energies were 49.8 keV and 91.6 keV, respectively, which was significantly higher than those of the other locations in the same buildings, 0.3 keV of F1-A (Fig. 8(a)) and 42.9 keV of F3-F (Fig. 8(f)). However, considering the high mean fluence rates at F1-A (3167 cm^-2·s^-1) and F3-F (1474 cm^-2·s^-1), the mean fluence rates at F1-B and F3-E were relatively low at 513 and 967 cm^-2·s^-1, respectively, as shown in Fig. 10. A higher neutron energy and lower fluence rate indicates that there should be an extra absorption material of thermal neutrons nearby; intermediate and fast neutrons that were not totally moderated were detected by the BSS.

Fig. 10Full-screen View

(Color online) Bar plot of the (a) total fluence rate, and (b) the average energy comparison given by MAXED, GRAVEL, and the genetic algorithm at different locations

Neutron ambient dose equivalent rate

It is easy to derive the neutron ambient dose equivalent rate $\dot{H^{*}}$(10) based on the neutron energy spectra and conversion coefficients that relate the neutron fluence rate to the dose equivalent rate. Considering that the results obtained by the three unfolding spectrum methods provide only an approximate value, the average value per location was considered as the reference value for a better comparison (Fig. 11). The maximum appeared at F3-E, which was approximately 209.3 μSv/h, owing to a higher average energy and fluence rate. Moreover, F1-A only had 1.7 times the $\dot{H^{*}}$(10) of the access route, although the total fluence rate of F1-A was far higher than that of F1-B in Factory 1. Neutrons above the keV range are prone to the reaction of recoil protons. The human body contains a large amount of water, thus the conversion coefficients rapidly rise when the neutron energy is above ~1 keV. Therefore, neutrons with a higher energy have a greater influence on $\dot{H^{*}}$(10) and the human body.

Fig. 11Full-screen View

(Color online) Bar plot of $\dot{H^{*}}$(10) given by the energy spectrum methods and LB6411 at different measurement locations

The $\dot{H^{*}}$(10) ratios of the direct measurement results (LB6411) to the energy spectrum methods exhibit large variations at different positions (Fig. 11, red point-and-line curve). The $\dot{H^{*}}$(10) measured by LB6411 was significantly greater than that derived from the spectrum methods at F1-A (more than 60%), whereas at F3-F, the $\dot{H^{*}}$(10) measured by LB6411 was lower than that derived from the spectrum methods. At F1-B and F3-E, the values obtained by the direct measurement and spectrum methods were not significantly different. These disparities resulted from the differences in the survey instrument responses in the various energy regions; the neutron spectra structure of the calibration laboratory was not the same as in practice. The $\dot{H^{*}}$(10) of the direct measurement was significantly higher than that of the spectrum methods, which may be due to the over-response of LB6411 for the thermal neutrons at F1-A. For F1-B and F3-E, the survey instrument responded properly with more intermediate-energy neutrons. However, for the neutron spectrum structures at F2-C, F2-D, and F3-F, the direct measurement results demonstrated a significant fluctuation, overresponse in certain locations, and under-response in others. In addition, according to the relative deviation, the calibration factors were 0.62, 0.94, 0.87, 0.78, 1.07, and 1.12, respectively, at the six locations and were different from 0.998 based on the ²⁴¹Am-Be neutron source. Considering a practical protection, these factors are more promising. Therefore, energy spectrum methods can effectively avoid the problems of over- or under-responses and help obtain a more accurate $\dot{H^{*}}$(10) in a complex field environment.

Evaluation of neutron spectra

There is no determined and unique spectrum; owing to the few-channel data of the BSS, certain criteria and conditions are required to verify the results of the unfolding spectra and ensure that they approximate the true values. However, all criteria are necessary and insufficient for the true spectrum. Three criteria, including the count rates of the inverse solutions, neutron spectrum structures, and key parameter comparison among the three methods, were used to evaluate the quality of the spectrum unfolding results.

3.5.1

Count rates of the inverse solutions and experiment

The count rate of the inverse solution of the neutron spectrum is one of the most important constraints. Although the spectra may not indicate the correct results despite the inverse solution values sufficiently matching the experimental results, poor results of the unfolding spectra can be excluded if they are inconsistent with one another. Figure 12 shows the absolute experimental count rate, inverse solution count rate, and their ratios for the three unfolding methods applied in this study. The inverse count rate of the GA approximated the measurement results more closely than MAXED and GRAVEL at most locations, and the maximum relative deviation was 18% for the 9 inches sphere at F1-A. It is difficult to quantify the uncertainty of the neutron spectrum involving a response function, measurement, experiment (temperature, scattering, etc.), unfolding spectrum, and other unknown factors, especially in an on-site complex environment. Therefore, a ±20% relative deviation between the measurement and inverse solution values may be considered satisfactory owing to the 20% response function uncertainty and 1% measurement uncertainty.

Fig. 12Full-screen View

(Color online) Plot of measurement and inverse solution neutron count rate against the sphere diameter at (a) F1-A, (b) F1-B, (c) F2-C, (d) F2-D, (e) F3-E, and (f) F3-F. The count rate ratio of the experimental results to the results derived by each of the three methods is plotted on the second y axis. The two red lines represent the region of deviation of the ratio that is less than 20%