Core design and source characteristics

The shielding design was performed based on a small helium xenon-cooled reactor, as shown in Figure 2. The internal geometric unit of the core has a repetitive hexagonal structure. The entire core is built with 1040 fuel rods with a diameter of 1.5 cm, surrounded by six coolant channels with a diameter of 0.8 cm. The overall radius of the core is 44 cm, the total height is 100 cm, and the design power is 20 MW.

Fig. 2Full-screen View

(Color online) Small helium xenon-cooled reactor. (a) Top view of the core. (b) Structural unit of a fuel rod. (c) Front view of the core.

The core is made of graphite as the base material, the fuel rods are filled with UC composed of uranium with an enrichment ratio of 19.75%, and the coolant cladding is made of Titnaium-Zirconium-Molybdenum (TZM) alloy. The core is wrapped with a 20 cm thick Be reflective layer. B₄C is a functional ceramic material with excellent properties, especially a high thermal neutron absorption capacity, the neutron shielding material selected in this study is B₄C. Stainless steel contains iron and is widely used for photon shielding. We selected SS-316 stainless steel as the photon-shielding material.

We defined the reactor source term as a volume source. For a reactor with a specified design power, when operating at full power conditions, the energy released by each fission is 180.912 MeV, where 1 MeV = 1.602×10^-13 J. Thus, the amount of nuclear fission required to produce 1 J of heat was calculated as follows: $\frac{1\text{ Mev}}{1.602\times {10}^{-13}\text{ J}}\times \frac{1\text{ fission}}{180.902\text{ Mev}}=3.45041\times {10}^{10}\text{ fission/J}\text{.}$ (1)

We assumed that the number n of neutrons produced by each fission is 2.45. If k_eff = 1, then every 2.45 source neutrons provided by the core will undergo a fission (the remaining fission neutrons are absorbed or escaped). The physical meaning of the above process is the number of fission neutrons required for every 1 J of energy released by the system per unit of time. If the power is P, then the source intensity A is calculated as follows: $A=P\times 3.4504\times {10}^{10}\times n=1.69\times {10}^{18}.$ (2)

Basic shielding schemes with empirical method

As shown in Figure 3(a), the total thickness of the shielding layer was set to 100 cm, and the thickness of each layer remained the same. We increased the number (limited to 3-16) of shielding layers to explore the variation law of radiation dose rate.

Fig. 3Full-screen View

(Color online) Determination of empirical shielding design scheme. (a) Initial shielding scheme. (b) Distribution law of radial dose along the axial direction. (c) Areas where dose extrema may locate. (d) Variation rule of dose with the increasing number of layers.

The schemes were further classified according to the sequence of shielding order for neutrons and photons. The statistical results for a large number of particles were the same as those for a large-scale repeated simulation of a single particle under the same conditions. However, with a continuous increase in the number of particles, it is not guaranteed that all particles will eventually be captured or absorbed by the shielding material, and some of the particles will pass through the shielding layer to reach the outer vacuum boundary. Therefore, by setting the detection area and counting the radiation dose within the area, the design scheme of the shielding layer can be evaluated. The dose was not uniformly distributed along the axis. To improve the statistical accuracy, a circle of air with a thickness of 5 cm was added at the periphery of the shielding layer to be the detection area, and 20 areas were uniformly divided along the axis, with a total height of 20×5 cm = 100 cm. As shown in Fig. 3(b), we calculated the dose in each area and determined the distribution law of the radial dose (mrem/h) in 20 areas along the axial direction of the six trial plans. Statistics show that the maximum dose value appears in the middle area, and the minimum value appears at the two ends, as shown in Fig. 3(c). Shielding calculations were performed using the Monte Carlo software RMC. The variation law of the maximum/minimum radiation dose in the obtained detection area with the number of layers and order of particle shielding is shown in Fig. 3(d).

As shown in Fig. 3(d), with the same total thickness, if the neutrons were shielded first (NSF), the dose rate and number of shielding layers were negatively correlated and ultimately converge to a certain range. Moreover, the radiation dose rate is lower when the total number of layers is even. If the photons were shielded first (PSF), with an increase in the number of shielding layers, the dose rate did not evidently decrease or converge. The dose rate was relatively low when the number of layers was odd. Compared to the scheme with the same number of layers and shielding neutrons first, the dose was reduced by 3–5 times. The dose difference between the odd and even numbers of shielding layers must be explained in conjunction with the particle transport process in the reactor. The scattering of neutrons in the reactor generates photons. If neutrons can be effectively shielded, it is equivalent to pre-shielding a portion of photons. For example, when neutrons are shielded first, and the number of layers is set to four, because the B₄C of the first and third layers have shielded most of the neutrons, the secondary photons generated by neutrons are small, and these photons can be absorbed by applying SS-316 to the outermost layer of the shielding body, and the dose is evidently lower than that of the three-layer shielding design with the same total thickness.

Therefore, two empirical shielding schemes were proposed: (1) Scheme I: shield photons first with an odd number of layers, and (2) Scheme II: shield neutrons first with an even number of layers.

Relative thickness optimization with GA

The empirical method determines the number of shielding layers and neutron/photon shielding order. These two parameters were not considered in the subsequent optimization and could be directly used as the modeling basis for the shielding body. However, the empirical method cannot optimize the relative thickness of the shielding layer because it only provides a scheme with uniform thickness. To further reduce the weight, the GA method was used to optimize the relative thickness of the shielding layers (Fig. 4).

Fig. 4Full-screen View

GA flowchart.

The genetic manipulation of radiation shielding scheme primarily includes population selection, crossover and mutation, which are discussed separately below.

(1) Selection

Suppose the group size is m and the fitness of individual i is F_i; then, the probability P of individual i being selected is $P=\frac{F_{\text{i}}}{{\displaystyle {\sum }_{i=1}^m{F}_i}}.$ (3)

(2) Crossover

We used the single-point crossover operator, randomly selected a position in the paired chromosomes, and then performed locus transformation on the paired chromosomes at the selected position.

(3) Mutation

To ensure the diversity of individuals, it is necessary to mutate the individual continuously in the iterative process. In this study, the single-point mutation method was adopted. That is, only a certain bit in the gene sequence needs to be mutated using binary coding as an example.

The mathematical function model for solving the optimal thickness of each layer under different layer-arrangement schemes is as follows: $\{\begin{array}{l}\min y\text{ = }F\text{(}x\text{) =}\left[f_3\text{(}x\text{), }{f}_4\text{(}x\text{)}\cdots \cdots f_{16}\text{(}x\text{)}\right]{}^{\text{T}}\\ \text{s}\text{.t}\text{. }{g}_i\text{(}x\text{) }\le \text{ 0}，i\text{ = }1,2,\cdots \cdots q\\ h_j\text{(}x\text{) = }r\text{, }j=1,2,\cdots \cdots p\end{array}$ (4) where x is the n-dimensional design variable, x = (x₁,x₂,…,x_n) $\in$ X, which includes the selection of shielding materials, the number of shielding layers and initial thickness of each layer. X is the design variable value space. A hybrid coding method was adopted. In other words, binary coding was used to express the order of each shielding layer in each scheme, and symbol coding was used to express the type of material assigned to a single shielding layer. Here, we chose only B₄C and SS-316 stainless steels as the shielding material input. y is the m-dimensional design objective. In this study, y includes the lowest level achieved by the maximum and minimum radiation doses in each detection area as well as the total weight of the shielding body. Simultaneously, the fitness function during the operation of the GA was determined according to the above design goals. f₃(x)–f₁₆(x) represent different schemes with 3–16 shielding layers, respectively. g_i(x) ≤ 0 defines q inequality constraints, h_j(x) = r defines p equality constraints.

When the iteration stop condition is satisfied, the optimal solution can be outputted according to the fitness ranking. The shielding design problem is essentially a multi-objective optimization problem; there will inevitably be conflicts between the various sub-objectives. To better reflect the importance of some sub-objectives, it is necessary to reasonably set the weight of each sub-objective in the fitness function. The Analytic Hierarchy Process (AHP) is adopted in this study to determine the weight of each index and select the optimal solution [25].

AHP is an efficient way to define the weights of each sub-objective in a multi-objective optimization problem. This method has the advantages of fewer information requirements and shorter decision-making time. The AHP steps are as follows.

(1) Create a positive reciprocal matrix $A={({\text{a}}_{ij})}_{\text{m}\times \text{m}}=\left[\begin{array}{l}1{\text{ a}}_{12}\cdots {\text{ a}}_{1\text{m}}\\ {\text{a}}_{21}1\cdots {\text{ a}}_{2\text{m}}\\ ⋮⋮\ddots ⋮\\ {\text{a}}_{\text{m}1}{\text{ a}}_{\text{m}2}\cdots 1\end{array}\right],$ (5) where ${\text{a}}_{ji}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\text{a}}_{ij}$}\right.$, and its value is determined using Table 1. For example, the total weight of the shielding body is compared to the lowest level that can be achieved by the maximum dose in the detection area. If the former is significantly more important than the latter, then ${\text{a}}_{12}=5$ and ${\text{a}}_{21}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$5$}\right.$. Each indicator was compared with other indicators, and the reciprocal matrix A was determined.

(2) Perform a consistency check.

Normalize each column vector of matrix A as follows: ${\tilde{w}}_{ij}={\text{a}}_{ij}/{\displaystyle {\sum }_{i=1}^m{\text{a}}_{ij}}.$ (6)

Sum ${\tilde{w}}_{ij}$ in rows to obtain ${\tilde{w}}_i={\displaystyle {\sum }_{j=1}^m{\tilde{w}}_{ij}}$.

Normalize ${\tilde{w}}_i$ to ${\tilde{w}}_i={\displaystyle {\sum }_{\text{i}=1}^{\text{m}}{\tilde{w}}_{ij}}$, where w = (w₁, w₂, …, w_m)^T is the eigenvector.

Calculate $\lambda =\frac{1}{m}{\displaystyle {\sum }_{i=1}^m\frac{{(A\text{w})}_i}{{\text{w}}_i}}$.

The consistency factor is given by $\text{CI}=\frac{\lambda -m}{m-1}$, where m is the order of the positive reciprocal matrix A. CI is compared with the same-order random consistency factor RI in Table 1; that is, the consistency ratio is $CR\text{=}\frac{CI}{RI}$. When CR < 0.1, it indicates that the positive reciprocal matrix A meets the consistency standard; otherwise, matrix A must be properly modified. When the positive reciprocal matrix A passes the consistency test, eigenvector $w={\left(w_1,w_2,\cdots ,w_m\right)}^T$ is the required weight vector.

To improve the probability of obtaining the optimal solution, the population size was set to 50, and the evolutionary generation was 100. Previous related research have shown that when the mutation probability is greater than 0.04, the result becomes unstable; we set the mutation probability to 0.02. Simultaneously, to obtain a better local search ability, the crossover probability was set to 0.5 [26]. Schemes with 7–13 shielding layers were selected as the calculation samples, and the statistical doses under different shielding layer schemes are summarized in Table 2 (unit: mrem/h) and Fig. 5.

Fig. 5Full-screen View

(Color online) Shielding and weight-loss effect analysis before and after optimization: (a) Radiation dose rate before and after optimization of two arrangements (NSF). (b) Radiation dose rate before and after optimization of two arrangements (PSF). (c) Analysis of weight-loss effect before and after optimization.

From Table 2 and Figs. 5(a) and 5(b), if the neutrons are shielded first, the dose rate output by the GA-optimized schemes is significantly lower than the initial series of schemes shown in Fig. 3(a). Furthermore, the maximum value can be reduced by 57%. If the photons are shielded first, as the number of layers increases, the difference between the dose rate obtained by optimization and that obtained by the initial scheme gradually decreases.

The relative thickness distributions of the optimized schemes are summarized in Table 2. Combined with the thickness distribution before and after optimization, the results in Fig. 5 can be divided into two types: the first type: neutrons are shielded first, and the total number of layers is odd; photons are shielded first, and the total number of layers is even; the second type: neutrons are shielded first, and the total number of layers is even; photons are shielded first, and the total number of layers is odd, that is, the two schemes presented at the end of Chapter 2. Considering the results in Table 2 and the particle transport process, most of the particles are absorbed because the shielding layer near the core is thicker. Simultaneously, in the first type of case, because the outermost layer is B₄C, which absorbs some of the remaining neutrons, the relative change rate of the first type of case before and after optimization will be more evident, as shown in Table 2. However, the overall dose level of the first type is higher than that of the second type. Therefore, we did not consider the two cases included in the first type. Also, when the number of layers in Fig. 5(b) increases, the dose before and after optimization does not change, which implies that the relative thickness of each layer decreases, and the absorption capacity of the corresponding particles worsens. This is similar to the scheme in which each layer is the same thickness.

The weight of each shielding scheme was then calculated, and the results are shown in Fig. 5(c).

Referring to Fig. 3(d), if scheme I is followed, the weight will also be reduced to a certain extent; if scheme II is followed, although the dose is reduced, the weight will be increased to a certain extent. Therefore, we recommend the use of Scheme I. Moreover, relative thickness optimization with GA helps reduce the weight of the shielding layers.

Determination of total thickness with quantitative correlation

The research results of relevant studies show that when the GA is applied to reactor shielding design, the general method for determining the genetic operators (population size, crossover probability, and mutation probability) is to preset different initial values [26]. Then, appropriate values of the genetic operators are determined by counting the average number of individuals required to achieve the optimal plan under different combinations. Finally, they are inputted into the algorithm framework. The above process involves a large amount of calculation and strong repeatability, which seriously affects the overall efficiency of the shielding design. To avoid repetition of the above process, this section proposes a method for constraining the algorithm directly and quantitatively by providing the total thickness of the shielding layer through the dose.

In this study, the GA-optimized scheme provides the relative thickness of each shielding layer determined according to the design target y in Equation (4). Because the setting of the total weight of the shielding body and dose in y are estimated values, that is, the upper limit of the corresponding value range, any solution within the range that satisfies this value may be regarded as the optimal solution. In addition, the algorithm exhibits the characteristics of a forwarding iteration. It is generally impossible to determine whether the total weight obtained by the sum of the relative thicknesses of each layer is the global optimal solution in the shielding design. Therefore, the total thickness can be determined by the shielding calculation as the quantitative constraint of the GA optimization process to avoid local optimal solutions [27].

The total thickness is related to the radiation dose requirement. In general, the best approach is to accurately perform Monte Carlo shielding calculations to determine the total thickness based on dose requirements. However, under the condition that the approximate value range cannot be obtained, the process of obtaining ideal results using the Monte Carlo method for the shielding calculation is usually random. To improve the efficiency of the design, we established a quantitative relationship between the total thickness and dose requirements through mass calculations in advance. The total thickness of the shielding layer can be directly determined according to the dose requirements to determine the complete shielding scheme.

As shown in Fig. 6(a), the thickness and dose exhibited an exponential change law. The thickness–dose design curves in the logarithmic coordinate system are shown in Fig. 6(b).

Fig. 6Full-screen View

(Color online) Establishment and validation of quantitative relationships: (a) Variation of dose with thickness for different shielding layers. (b) Thickness–dose design curve. (c) Comparison of dose values obtained by quantitative equation with those calculated by RMC under different thicknesses (mrem/h).

We first established a set of equations between the dose rate and total thickness under a constant power (20 MW):

The maximum equation is calculated as follows: $T=148.06-15.15D.$ (7)

And the minimum equation is calculated as follows: $T=140.34-14.93D.$ (8)

In the equation, T is the total thickness of the shielding layer (cm) and D = lga, where a = dose (mrem/h).

To verify the calculation accuracy of the above equation, three sets of thicknesses (T₁–T₃) were randomly selected, and the maximum and minimum dose results (D_1Max–D_3Max; D_1Min–D_3Min) were calculated using Eqs. (7) and (8). T₁–T₃ are substituted as quantitative constraints into the GA toolbox, that is h₁(x) =T₁, h₂(x) = T₂, h₃(x) = T₃. The total shield weight and dose range were set simultaneously. The design variables were entered into the x-dimensional vector space of Eq. (4), and the fitness function was set according to the AHP method. The obtained calculation results, that is, the relative thickness of each layer, together with other design variables, were used as input parameters of the RMC. The results of comparing the dose calculated by RMC with the dose calculated by Eqs. (7) and (8) are shown in Fig. 6(c).

As shown in Fig. 6(c), compared with the results calculated by RMC, the relative error of the extreme values of the dose in all detection areas calculated by Eqs. (7) and (8) is within 10%. The validity and accuracy of the quantitative solution given by the above equation were verified and can be applied to the rapid calculation of the shielding design of small reactors.

Generally, reactors with power levels in the range 5 to 200 MW are referred to as small reactors. Based on the above calculation results, to improve the universality of the equation, the power range was extended to 5–200 MW, and the thickness–dose relationship was specified for the different power levels. This process is illustrated in Fig. 7. Under actual design conditions, designers tend to focus on the maximum dose in the detection area; the relationship between the thickness and maximum dose is provided.

Fig. 7Full-screen View

Flow chart of quantitative shielding design method.

Combined with the conclusions obtained in Chapter 2 and Sect. 3.1, assuming that the particle shielding order is determined, it is evident from Table 2 that when the number of shielding layers is set to 7, the GA-optimized shielding scheme has the smallest dose value compared to those of the other schemes. Therefore, the number of shielding layers was set to 7. The results are shown in Fig. 8.

Fig. 8Full-screen View

Thickness–dose relationship curve at different power levels.

As observed, the thickness–dose relationship at different power levels continues to change exponentially. The quantitative relationship can be summarized as follows: $T=\text{a}\times {\text{D}}^{\text{b}}.$ (9)

As shown in Fig. 7, a represents the dose accumulation factor, and b represents the dose correction factor.

When the design power of the reactor is given, first, the specific values ​​of the accumulation factor a and correction factor b are determined through the relationship between them. Further, a and b are substituted into Eq. (9). Subsequently, the corresponding total thickness of the shielding layer is quantitatively determined according to the dose limit. Combined with the empirical design method and multi-objective optimization method proposed above, the final optimized design scheme can be obtained.

Evaluation of shielding and weight-loss effect

To verify the effect of the semi-empirical and semi-quantitative shielding design method on radiation shielding and weight reduction of the shielding body simultaneously, we set up three groups of samples for comparison. First, the power of the three groups of samples was randomly set to P = 37.5 MW, the dose D₁ = 100 mrem/h, and the total thickness of the shielding layer was obtained according to Eq. (9), that is, T₁ = 119.7 cm. Further, based on the empirical shielding design method, the shielding body model is established based on Scheme Ⅰ. The material, number of shielding layers, and initial thickness of each layer were entered as design variables in the GA model contained in Eq. (4). With D₁ and T₁ as quantitative constraints in the optimization process of the GA, corresponding to Eq. (4), h₁(x) = 100 and h₂(x) = 119.7. Finally, the optimized relative thickness from the inner to outer radius is shown in Figure 9(a) and Table 3.

Fig. 9Full-screen View

(Color online) Effect verification. (a) Validation model for shielding and weight-loss effects. (b) Dose comparison between D₂ and D₁.

Compared with the shielding scheme with equal total thickness, that is, T₂ = T₁ = 119.7 cm, and equal relative thickness of each layer, that is, t₂ = 17.1 cm, the neutron/photon coupling transport calculation is performed. The total weight of the shielding layer and dose in each detection area of the two schemes are listed in Table 3 and Fig. 9(b). It can be seen that under the condition of the same total thickness, the semi-empirical and semi-quantitative shielding design method can effectively help reduce the radiation dose by approximately 2.3 times. The maximum dose in the detection area decreased from D₂ = 235.23 mrem/h to D₁ = 101.04 mrem/h. Similarly, as listed in Table 3, when the dose requirements are equal, that is, D₃ = D₁ = 100 mrem/h, the semi-empirical and semi-quantitative shielding design method can effectively reduce the total weight of the shielding body from 60.71 t to 53.52 t. This is a reduction difference of 11.84% compared to the normal method (same dose requirement).

In addition, we provide a set of two-layer design schemes. To ensure the same order of particle shielding, SS-316 and B₄C were set up with one layer. The thicknesses of the two layers are given by referring to the scheme in the second column of Table 3; that is, the sum of the thicknesses of the two layers is 119.7 cm. The maximum dose calculated by this scheme is 4527.29 mrem/h, which is 42 times higher than the scheme optimized by the semi-empirical and semi-quantitative lightweight shielding design method.

To summarize, the power is considered an independent design variable, and the accumulation factor a and correction factor b are derived. Finally, the thickness value was obtained by substituting the two factors into Eq. (9). This method can restrict the total thickness before the specific optimization of the thickness of each layer by the GA. This can effectively avoid determining the local optimal solution owing to its large calculation scale.

Shield design optimization program SDIC1.0

We encapsulated the semi-empirical and semi-quantitative lightweight shielding design method into SDIC1.0.

The semi-empirical and semi-quantitative lightweight shielding design method is divided into three steps: (1) The number of shielding layers and neutron/photon shielding sequence are determined first by trial calculations, and an empirical scheme is established; (2) The thickness–dose relationship is established through the statistics of mass Monte Carlo calculation results, so as to determine the corresponding total thickness of the shielding layer through the dose requirements, and use this thickness as a constraint to prevent determining a local optimal solution when optimizing the relative thickness of each layer. (3) The relative thickness of the shield layer was optimized by the GA, and the weight of the shield layer was reduced. To encapsulate the semi-empirical and semi-quantitative lightweight shielding design method, we developed the shield optimization design program SDIC1.0. The program process is illustrated in Fig. 10.

Fig. 10Full-screen View

Process involved in SDIC1.0.

The accuracy of the procedure was verified by randomly setting different groups of power levels and dose limits based on a consistent core scheme. As shown in Table 4, the three groups of random power are 60 MW, 140 MW, and 110 MW, respectively. First, the accumulation factor a and correction factor b were obtained from the input power. Subsequently, three groups of dose limits at different power levels were correspondingly set, that is, D₁ = 10000 mrem/h; D₂ = 5000 mrem/h; D₃ = 1000 mrem/h. The corresponding total thickness of the shielding layer is given according to the dose limits. Furthermore, the content in Sect. 4.1 was repeated, the design variables and constraints were entered, the algorithm was run, and the output, including the relative thickness of each layer, was obtained. Finally, the thickness of each layer optimized by GA was used as the input of the RMC, and the calculated dose was compared with the initially set doses D₁–D₃. The results are presented in Table 4.

As shown in Table 4, the relative error between the results of SDIC1.0 and RMC can be effectively controlled at 9.63%, 10.40%, and 11.40%, respectively. By analyzing the reason for the error, it can be speculated that in determining the relationship between power and the two factors a and b, fewer sample points were used, resulting in a certain deviation in the fitting degree between the three. Therefore, in the future, the error caused by the calculation results of the program can be reduced by increasing the density of the initial sample points.

In addition, in terms of the calculation time, as shown in Table 4, T₁ is the time required to obtain the optimal solution using the method proposed in this paper, and T₂ is the time required to obtain a result similar to D_t by simply using the Monte Carlo code RMC. The calculation time can be shortened by approximately 6.3 times by using this method.