The effect of control rods on the reactivity and flux distribution of BWR 4 bundle using MCNPX Code

This paper has three main objectives related to the neutronic and burnup analysis of the BWR (Boiling Water Reactor) Four-Lattice. The first objective is to provide partial validation of the MCNPX code for this lattice by comparing its results with Scale-5.1 results. Validation of the MCNPX to calculate effective multiplication factor and reactivity rod worth for the F-Lattice is provided. This is carried out in case of instantly removing the control blade and replacing it with a graphite moderator. Moreover, spatial neutron flux distributions using F-mesh card over the bundle and the control blade are investigated at inserting and withdrawing the B4C. The second objective is to perform parametric design studies of the F-Lattice. Areas of particular interest are the effect of increased or decreased blade width on the neutron flux throughout the bundle. It is found that the presence of carbon in the control blade at withdrawing the B4C makes the reactor supercritical, (K-eff = 1.22206). On the other hand, the use of B4C blade presents (K-eff = 0.93521). Consequently, the reactivity of 10% B4C thinner case is higher that of 10% B4C thicker. The simulation also showed that the B4C blade had an effective role in decreasing the thermal flux at the periphery of the bundle. This is contrast to the effect of carbon that moderates fast to thermal neutrons. The third part of this work aims at studying the burnup calculations using MCNPX code for 30 days burn with 1 day time step then for 20 months burn with 2 week time steps for the lattice. At the end of the work, it is very important to determine the most proper bundle model that achieves a prolonged fuel burn and flatting thermal flux distribution. For reaching this goal, three cases (B4C, 10% thinner of B4C and 10% thicker) are simulated by MCNPX code till 70 GWd/ton. It is found that the B4C and 10% thicker are the appropriate models that can satisfy the safety considerations of the Compact Modular Boiling Water Reactor.

www.nature.com/scientificreports/ When F-Lattice is utilized in the ESBWR (Economically Simplified Boiling Water Reactor) design, only half the number of control rods is required in the core. Fewer control rods means a reduction in material and manufacturing costs, fewer drives and hydraulic control units. Furthermore, there would have been approximately half the number of control rod related components requiring maintenance or posing the risk of failure. This is accompanied by a reduction in cost and an increase in safety; the two key goals of all design modifications in a nuclear reactor. While the F-Lattice is no longer a part of the ESBWR, it still remains a key component for the CM-BWR design. Therefore, study and analysis of this new lattice design must be performed before the CM-BWR design can be completed 1,3 .
Generally, the production of electricity by boiling water reactors is one of the basic considerations of power plants because of reactor power raising and higher fuel assembly burnups. Highly reliable thermal margins and stability assessment are required for BWRs operating at higher powers. Higher burnup required higher average U-235 enrichment in the fresh assemblies and led to more heterogeneous cores, which stressed the thermal and shutdown margins of the reactor. The required better thermal and shutdown margins were gained by significantly improving BWR fuel assembly designs 4 .
The fuel characteristics of BWRs are studied less than PWRs. However, BWRs are widely deployed internationally. This is attributed to the complexity of the modeling associated with their highly heterogeneous fuel assemblies [5][6][7] . This has been recently confirmed for an advanced BWR with a published core design which had been modeled by experts in reactor physics 8 .
The infinite multiplication factor for the F-lattice is evaluated when the thickness of B 4 C blade is increased and decreased by 10% of its nominal thickness. It is expected that the decreased thickness of B 4 C leads to an increase in the effective multiplication factor due to the low boron content in this case.
The reactivity worth of control rod is another important safety parameter that is affected by the variation of isotopes during fuel burnup and the operational time of the rod and its position in the core 9 . This parameter can be calculated using the equation ρ = K−K 0 K , where K 0 is the multiplication factor with control rods inserted in, and K is the multiplication factor with the control rods pulled out and replaced with pure graphite. This means that when B 4 C is replaced by graphite, neutron absorption dramatically decreases, causing an increase in neutron population, increase in flux and a larger effective multiplication factor. All the obtained results are accompanied by calculating the standard deviation and the percent error relative to Scale-5.1 code.
Because of the strong dependence of reactor safety on the flux distribution, F-Mesh card is used for plotting the thermal, epithermal, fast and total flux distributions in xy-plane of the F-bundle. To accomplish this, the three dimensions of the assembly geometry are divided into 10 meshes. Spatial flux distribution is not only studied for the fuel F-bundle, but also plotted over the control blade at inserting and withdrawing the B 4 C. This is conducted for distinguishing between the role of both B 4 C as a strong absorber and graphite as a strong moderator in determining the flux profile.
The last section of this work, depletion calculations of BWR bundle, focused mainly on searching for which model among the three models under investigation (B 4 C(nominal case), 10% B 4 C thinner and 10% B 4 C thicker) can achieve a prolonged life cycle and flattened thermal flux distribution. These two features are one of the desired safety considerations of the BWR cores for enhancing the stabilization of reactor operation. As we will see, the B 4 C and 10% B 4 C thicker can satisfy these goals.

Methodology
The MCNPX code was developed by the Los Alamos National Laboratory. It is a general purpose Monte Carlo code, which facilitates independent or coupled neutron, photon and electron transport calculations 10 . All the neutronic and burnup calculations in this research are conducted by this code. It is capable of modeling complex 3D geometries and utilizes extensive point-wise cross-section data library on a continuous energy spectrum. The cross section data library used in the calculations is the updated library, endf71X. This data library provides the necessary probability distributions for simulating particle interactions through use of random number sampling 11 . Through a large number of neutrons and active cycles, MCNPX is capable of determining k-effective, neutron flux, current and many other neutronic parameters. Increasing the number of histories and the number of neutrons tracked in each history provides more detailed results, though at the expense of higher computational costs.
Flux calculations require the use of a mesh, which divides the model into different cells. Each individual cell describes a grouping of fuel, gap, clad and surrounding water as one geometric unit. Flux results are then averaged within each geometric unit. A more refined mesh leads to smaller cells, which reduces the effect of averaging. Smaller cells lead to more accurate flux profiles. The flux, or power, distribution can be determined with the average flux results for each cell. The superimposed mesh that covers only desired specific locations in the model can be specified in MCNP6 which is involved in MCNPX code. The FMESH tally is a general track-length tally which does not record tracks in pre-assigned geometry cells, but uses a mesh superimposed over the problem geometry 12,13 . The mesh can be defined in Cartesian or cylindrical coordinates (Supplementary Material).

The BWR bundle under investigation
The F-lattice was modeled by MCNPX code as a collection of unit cells. A unit fuel cell is composed of the fuel, gap, clad and surrounding moderator. This unit fuel cell was then repeated using the MCNPX "lattice fill" option to create an 8 × 8 unit bundle. A SS-304 canister was then modeled around the unit bundle. This unit bundle was then repeated four times with moderator spacing between them. This was done using the "like m but" command in MCNPX. This command allows a unit to be repeated with slight changes in material or physical location. The control rods were then modeled around the four bundles. This completed the 4-bundle cell. The parameters and dimensions of F-bundle are depicted in Table 1. A schematic diagram of F-lattice is given in Fig. 1 www.nature.com/scientificreports/ used was the "ksrc", which allows the user to define specific locations for the initial sources to be located. Sources were placed in the center of several different fuel rods throughout the system. MCNP6 automatically expands its source locations with each generation history, using information from each previous history to better approximate the next. Finally, the F4 tally was utilized for the flux analysis. Neutron energy levels can also be specified with the F4 tally (e0 card). The flux results are normalized to 1000 MWth.

Results and discussions
Studying the reactivity of different models of the F-lattice. Firstly, the effective multiplication factor results for the F-Lattice with B 4 C and C are obtained by MCNPX and then they are compared with SCALE-5.1 K-eff results as shown in the previous Table 2. It is expected that the K-eff in case of inserting B 4 C blade  www.nature.com/scientificreports/ (0.93521) is lower than that is produced at withdrawing the B 4 C blade (1.22206). This is due to the high absorption cross section of B-10 which in turn decreases the thermal absorption in the fuel, consequently decreases the number of fission events. Therefore, its K-eff is lower. In case of using carbon, the carbon is considered a good moderator that converts the fast neutrons into thermal neutrons which participate in increasing the fission processes in the fuel. The percent error is calculated relative to SCALE-5.1result. Relative percent differences are calculated through the equation, It is necessary to interfere that the 4-bundle MCNPX case used 100,000 particles for 250 inactive cycles and 550 active cycles. More particles and histories were used in this case for evaluation of the flux distribution, which requires more information than the k-effective calculations. The SCALE-5.1 case was run with 10,000 particles for 200 inactive cycles and 400 active cycles. The relative percent difference between MCNPX and Scale-5.1 is included to show the level of agreement between the results of the two codes (MCNPX and KENO). This relative difference is shown to be around one percent (0.87% for the B 4 C case and 1.19% for the graphite case).
Secondly, it is also noted that the rod worth of the B 4 C control blades was determined by MCNPX code to be 0.2347. On the other hand, the SCALE-5.1 result leads to a calculated rod worth of 0.2188. Overall, a rod worth of about 0.2347 indicates that the B 4 C control blades (fully inserted) have a very strong effect on the core reactivity. This is beneficial to overall core life, indicating that a gradual removal of the control blades from the core will allow for a prolonged fuel burn. The calculated relative percent error for this case is 7.26%.
Thirdly, the variation in the thickness of the B 4 C control blade has an impact of the bundle reactivity. Studying this parameter (variation of thickness) confirms that MCNPX code is capable of calculating the effective multiplication factor with less than one percent error for the F-Lattice when the thickness of the control blade is increased or decreased by 10%. In these two models, the bundle pitch remains constant but the thickness of the control blade is varied by ± 10%, thus affecting the size of the water gap between the fuel and the control blade. The variation of the thickness of the control blades also results in a significant change in the amount of control material present in the bundle.
According to the results, it is observed the 10% thinner model presents higher K-eff (0.93883) than that of 10% thinner (0.93152). The 10% thinner case leads to an increase in the effective multiplication factor by 0.362% from nominal case (B 4 C case). The 10% thicker case produces a decrease by nearly the same factor (0. 369%). The (±) 10% thickness variation leads to overall reactivity changes of about (±) 0.004. Relative error for MCNPX compared to SCALE-5.1 remained below one percent for this variation (0.76 and 0.82) for − 10% thinner and + 10% thicker respectively.
Fourthly, the rod worth calculations were carried out because of the significant increase or decrease in control material volume. From the tabulated data, it is obvious that rod worth is directly proportional to the ratio of control blade surface area to fuel volume. And the fuel volume is constant. Therefore, the rod worth is directly proportional to the control blade area. The rod worth is higher for the 10% thicker. The relative percent error is around 7% (6.95% for thinner case and 6.48% for thicker case).
Finally, Reactivity changes are expected to be small due to the fact that the control blades' absorption crosssection, the control blade surface area and the fuel volume are not significantly affected in this parametric study. The three cases of nominal control blade thickness of 8.3 mm, 10% thicker (9.13 mm) and 10% thinner (7.47 mm) are studied with MCNPX code with 50,000 particles, 250 inactive cycles and 350 active cycles.
Comparison between the thermal, epithermal, fast and total flux distribution of the F-bundle in case inserting and withdrawing the B 4 C control blade. In this section, the flux distribution is studied over the fuel bundle in case of using B 4 C (nominal case) and C in the control blade. This is carried out by  where P = Thermal reactor power in MW (1000 MWt), Q = Energy released per fission, MeV (200 meV), υ = Average number of neutrons released per fission, determined by MCNPX code, K-eff = Effective neutron multiplication number, determined by MCNPX. The neutron flux profile is an important characteristic in determining the behavior of the reactor core and understanding the isotopic composition and the contributions from Plutonium. For these two purposes, the thermal, epithermal and flux profiles for the BWR bundle are plotted and analyzed 15 .
From Fig. 2, there is a decrease in the thermal flux values at the bundle periphery due to the presence of B 4 C that absorbs thermal neutrons. The peak of the thermal flux is found in the middle of the bundle owing to the presence of water that moderate fast neutrons to thermal neutrons. On the other hand, in Fig. 3, the thermal flux peaks found in that figure are located near the graphite control blade that does not absorb thermal neutrons. A central thermal flux peak is also found in the middle of the bundle due to the reason mentioned before. The very sharp decrease in both figures towards the Y-axis in blue and violet colors is due to locating the mesh cells of this part outside the bundle. This part summarizes that replacing B 4 C with graphite would actually result in a larger power level, and also larger flux values. 1000 MWt was maintained in both cases to allow for a direct comparison of the effect the material change has on the flux distribution. The constant power normalization is the reason that the peak flux values in the thermal flux (Fig. 2) are larger for the nominal case, compared to the graphite case (Fig. 3). Figures 4 and 5 show that there is a significant difference in the epithermal flux distribution in the two cases under study. It is clear that the graphite case contributes to flat and lower epithermal flux profile through the bundle than the B 4 C case. This may be attributed to the interaction of fast neutrons with carbon and this resulting in lower epithermal flux over the F-lattice. For the B4C model, the absorption of thermal neurons causes plenty of epithermal and fast neurons over the bundle fuel. Figure 6 gives the fast flux distribution in case of B 4 C case. Four fast flux peaks are found in the figure. Each peak corresponds to a fuel bundle. This means that decreasing the fission reactions in the four bundles contributes in increasing the number of fast neutrons over the four bundles. A very sharp decrease in the fast flux profile is also found because of the void mesh cells located outside the F-lattice. www.nature.com/scientificreports/ The behavior of fast flux profile at withdrawing the B 4 C is nearly similar to that is obtained in case of inserting B 4 C but the values are lower in the carbon case. The scattering of fast neutrons by carbon around the bundle makes fast neutrons to be lower in the F-lattice (Fig. 7). Figures 8 and 9 give the total flux profile which is the combination between thermal, epithermal and fast flux. The total flux is essentially flat for Carbon model.   Thermal neutron flux shown in Fig. 10 provides the most relevant example of the model's rod worth. The strong absorber, B 4 C results in near zero thermal neutron flux through the control blade. The stronger moderator, graphite, is shown to maintain a high thermal neutron flux throughout the blade, which is directly related to the significantly larger k-effective associated with the graphite system. Figure 11 is a plot of the epithermal neutron flux. The graphite system presents a larger flux profile through the blade than the B 4 C system. This is expected as more of the fast neutrons are slowed down via interactions with the graphite than are by B 4 C. Figure 12 provides the flux profile for the fast neutrons. As expected, the graphite system again presents larger flux values, though it is not as significant as it is in the epithermal or thermal cases. Slightly more fast neutrons are present in the graphite case due to the increase in overall fission reactions that are taking place compared to the B 4 C control rod system.   16 . The production of poisonous isotopes, such as Xe-135, quickly accumulates and reaches to an equilibrium state, which in turn helps in regulating the reactivity of the core. At the same time, fission products, such as Cs-137, are produced. All of these reactions have an effect on core life, power distribution, flux distribution and reactivity. How the core changes with time is an important aspect of any reactor  Two tests were carried out on the F-Lattice, assuming a fresh core with an average burn rate of 9.80 MWt (control blades remain fully inserted over the entire duration of burn). The first burn-up calculation was performed for 1 month (31 days). MCNPX code performed 31 burn steps (1 day intervals), which refers to the time step in which neutron fluxes are actualized with time.
The second test was performed for 20 months with 40 steps (2 week time steps) and 40 internal steps. K-effective calculations are only presented in this test. (For a preliminary study, and due to time constraints, MCNPX simulations were run with reduced accuracy, using only 10,000 particles for 150 inactive cycles and 150 active cycles. The obtained results by MCNPX code are compared with those obtained by (MCNP5-Monteburns) combination. The idea of coupling MCNP5 with Monteburns in reference 1 depends on calculating the    The results of the 1 month burn-up calculation with 1 day time step. The large drop in k-effective on the first day of burn is a result of the production of Xe-135 present in the reactor (Fig. 14). This well known Xenon poisoning which in turn has a fast and dramatic effect on the core reactivity before it reaches an equilibrium state at which the mass of Xenon is nearly 0.08 gm as in Fig. 15. Fission product accumulation of Sm-151 and Pu-239 are linear as in Figs. 16 and 17. The mass of U-235 decreases linearly (Fig. 18). Results of K-eff for 20 month burn with 2 week time steps are depicted in Fig. 19.
The final k-effective of Reference results, at the end of 30 days, is 0.9072, a 3.3% decrease from the initial value of 0.9380. But for MCNPX results, the final value (0.9455) decreases by 3.1% from the initial value (0.93556). As observed in Fig. 20 and Table 3, the highest values of K-eff are obtained at withdrawing the B 4 C blade. The lowest values are for 10% thicker of B 4 C. This shows that as the content of absorbing material increases, the reactivity decreases. Figure 21 and Table 4 depict the variation of U-235 mass with time. In all cases, the mass decreases with time. From the tabulated results, the mass of U-235 at 70 GWd/ton for B 4 C and 10% thicker cases is higher than that of the other cases (C and 10% thinner). Therefore, those two models achieves the longest life cycle among the 4 models.
Flux distribution for C and B 4 C models at 70 GWd/ton (4.679E+03 days). At the end of reactor operation, it is very important to study how the thermal flux is distributed at inserting and withdrawing the B 4 C blade. Thermal flux profiles are plotted by the way illustrated before. Figure 22 show the spatial thermal neutron flux for the B 4 C model at the end of cycle. Based on these results, we can see the trend where maximum flux is in the centre of the bundle and gradually reduced towards the edges due to the presence of B 4 C that absorbs thermal neutrons. The significant peak at the bundle centre is also because it has larger water volume compare to other locations in the geometry. Hence more moderation process occur which increased thermal neutron flux at the centre of F-lattice. Therefore, it is obvious from the figure that the thermal neutron flux has no fluctuations in the bundle that is resulted from uniform power distribution over each bundle.     www.nature.com/scientificreports/ At withdrawing the B 4 C (Fig. 23), the thermal flux fluctuations are observable due to the presence of carbon that participates in decreasing the thermal flux over the four bundles and increasing it at the edges and the middle. So, the radial thermal flux profile was found to be non-cosine shape. Finally, it can be said that the B 4 C model achieves more safety and long operation for the Boiling Water Reactor. Increasing the ratio between the average flux and maximum flux is evidence that a flattening behavior of thermal neutrons is obtained for the B 4 C

Conclusion
MCNPX is found be capable of modeling and simulating the BWR F-lattice. Four cases are investigated. These cases include inserting, withdrawing, using 10% thinner and using 10% thicker of B 4 C. A comparison is held between the four cases from the point of view of multiplication factor, reactivity and rod worth. The results showed the role of B 4 C control blade in decreasing the reactivity of the bundle. Spatial thermal, epithermal and fast fluxes are studied for the B 4 C and C cases. This is carried out by MCNP6 found in MCNPX. The flux results summarized that thermal flux over the B 4 C material is lower at the bundle periphery and is peaked at the centre but for the C material, thermal flux peaks are found at the periphery due to the moderation effect of carbon. For investigating which model that provides long life cycle and flatting thermal flux at the end of operation, depletion calculations are conducted till 70 GWd/ton. It was clear that the B 4 C and 10% thicker models provide prolonged life cycles due to the high amounts of U-235. The results also concluded that the flattening of thermal flux is observed for B 4 C model and this was predicted via ɸ average /ɸ maximum . This gives high stability and safety for the BWR operation.