Research on the automatic generation code for nuclear fuel reloading patterns in pressurized water-cooled reactors

ARTICLE INFO A method for automated generation program for nuclear fuel reloading patterns in Pressurized Water Reactor (PWR) has been developed. This newlydeveloped method consists of six different steps to minimize the maximum FΔH value, and maximize the reactor cycle length. Step 1 is initial fuel placement that is expected to produce the longest cycle length possible with the selected Fuel Assemblies (FAs) for the current cycle. Step 2 is aiming to decrease the FΔH value of the FA with the maximum FΔH. Step 3 aims to increase the FΔH value of the old FA with the lowest FΔH. Step 4 is rotating FA with the lowest FΔH value to increase its FΔH value, and rotating several old FAs in the neighboring FA with the maximum FΔH value to decrease the maximum FΔH value. Step 5 is aiming to increase the FΔH value of FA with the lowest FΔH value. The last step or step 6, will try to move FAs that have high k∞ in the periphery zone, inward to increase the cycle length of the reactor. These steps are translated into code in the Python programming language to enable automatic execution in a computer. A 3D nuclear reactor core neutronic code, COCO, is used for the calculation of FΔH value and reactor cycle length. Every nuclear power plant designer company will have their FΔH peaking factor safety limit in accordance with their DNB experiments and calculations, and the FΔH value safety limit used in this research is 1.46. A PWR loading pattern model is used to test this method. During the test, all the steps in this method are successfully executed in a total of 25 iterations plus one initialization calculation and produced acceptable results. The results of this method are all of the loading patterns found in all steps which have the maximum FΔH value below the defined criterion values. In the mentioned PWR loading pattern model, four optimized loading patterns are found using this method, all of which can be selected in the PWR refueling loading pattern design. This is an open-access article under the CC–BY-SA license. S Article history Received: 7 April 2021 Revised: 28 June 2021 Accepted: 9 July 2021


Introduction
Besides the safety characteristic, another important aspect of a nuclear power plant is its economics. Both aspects have to be balanced to ensure the successful operation of a nuclear reactor. An effort to balance the safety and economics of a nuclear reactor is by employing incore fuel management. The objectives of in-core fuel management usually are maximizing cycle length or time between each reloading and flattening the power distribution in the core. Maximizing cycle length will improve the economics of the reactor because the reactor will produce more power before being subcritical. Flattening power distribution will improve the safety of the reactor since fuel element structural integrity will likely be compromised in the region with high power, and will have the potential to release fission products.
In the present day, the economics and safety of pressurized water reactor (PWR), is greatly improved by the employment of more sophisticated loading pattern optimization strategies. These new strategies are taking advantage of the computation technology and are mainly divided into two categories; deterministic method and stochastic method. Deterministic methods have the advantage of their effective computational resources with acceptable loading pattern results, but they might not produce the most optimum solution. Examples of deterministic approaches are using mixed-integer nonlinear optimization method [1], and Lagrange multiplier and direct adjoining approach optimization [2]. In contrast, the stochastic method will need an extensive computational resource, but it is a robust method to solve a complex optimization problem, and resulting in an optimum, or at least a near-optimum solution [3]. One disadvantage of stochastic methods is during its application, the solution of stochastic methods might be trapped in local optimum, so that multiple initial conditions have to be tried, thus adding the computational cost of these methods. Examples of stochastic approaches are using dynamic programming models [4], simulated annealing, and genetic algorithms [5]. Reviewing methods and approaches from various research is useful to design another loading pattern generation method. Lesson learned and optimization parameters from past research could be applied in the new optimization method.
The newly-developed automatic loading pattern generation described in this article is preliminary research and could be categorized as a deterministic approach, thus shares its advantages and disadvantages of the deterministic method mentioned above. However, compared to the deterministic methods mentioned in the paragraph above which taking into account the microscopic parameter such as neutron diffusion, and fuel mass distribution, the methods described in this article is a rule-based method that based on a macroscopic parameter such as effective multiplication factor (keff), and maximum enthalpy rise hot channel factor (FΔH). The rule-based methods have been developed for decades, and every researcher and expert on fuel reloading has their own ideas and rules in their methods. Examples of these methods are using 'IF-THEN' rules developed by a fuel manager [6], and an integrated expert system called Efficient and Safe Optimum In-core Fuel Management (ESOIFM Package) [7]. The rule-based methods or often also called heuristic methods is fast to calculate and could also be combined with stochastic methods to make a hybrid method to reduce calculation time and find more optimum result [2], [8].
Calculations of loading pattern optimization are determined at the end-of-cycle condition. Important parameters to compare and optimize, generally are keff, power distribution, soluble boron concentration, FΔH, uranium-235 inventory, and cycle length in effective full power days (efpd). Optimization effort needs an objective function to be achieved. The objective function will determine which parameters to be maximized or minimized. However, some constraints limiting certain parameters or variables, and they cannot be exceeded. For example, optimization by maximizing the discharged burnup will give considerable fuel cost-saving, but there is a burnup constraint that will limit how much burnup value is allowed to achieve before the fuel assembly structure could no longer withstand conditions in a full-power reactor core [9]. In actual conditions, maximizing discharged burnup does not always mean maximizing reactor cycle length, since these two parameters are not proportional to each other. Sometimes in a single cycle, discharge burnup shows a reverse trend against cycle length. Therefore, discharge burnup is considered as a long-term fuel management optimization index, while cycle length is considered as a short-term or single cycle optimization index [10].
The above methods are examples of how to balance nuclear reactor safety and economics in neutronic aspects. Balancing nuclear reactor safety and economics could also be achieved from its thermohydraulic aspects. The Economics of a nuclear reactor could be improved by operating it as near as possible to its safety limits. To do that, a departure from nucleate boiling (DNB) point calculation during normal and accident situations should be very accurate and precise. Besides, predicting DNB during accident conditions is very important to mitigate reactor damage, thus minimalize economic loss during minor or even severe accidents. Many experiments and computer simulations had been done to precisely and accurately predict DNB. One of the efforts is described in the reference [11], which uses a computer to simulate the coupling of neutron transport and thermal-hydraulics calculation to estimate DNB and also analyze its uncertainty. Another nuclear reactor safety concern when increasing its power is the thermal limits of the current fuels. These thermal limits are one factor that is limiting future power level increases since further local power increases in the core could result in fuel damage under anticipated operational occurrences. A new type of fuel, that had been envisioned to replace typical uranium oxide fuel and zirconium cladding, is called accident tolerant fuels (ATFs). ATFs are currently under research and development and designed to overcome current fuel thermal limits. The idealized goal for ATFs would be to provide substantially improved safety response to a design-basis, or beyond-design-basis accident while maintaining good operational characteristics under normal conditions [12].
During reactor operation, neutron flux and power distribution in the core are not uniform. This condition is causing many fuel assemblies in the core-periphery to operate far below the license limit, while fuel assemblies near the center of the core are operating near the license limit. The potential of localized power peaking needs to be considered during the fuel loading pattern design phase. One of the most important peaking factors is FΔH. FΔH is chosen as the safety criterion over the average fuel assembly power factor because limiting the average fuel assembly power factor does not ensure DNB will not occur. DNB depends not only on the local power density but also on the coolant flow rate and the local enthalpy. Maximum heat flux may occur in several core elevations in the same coolant channel or fuel assembly, so that heat flux will be limited at each location. However, the coolant enthalpy increase in that coolant channel or fuel assembly will be more than in any other in the reactor, thus DNB will be more likely to occur in that channel because heat flux that causing DNB is lower in the same location with the occurrence of higher enthalpy. To protect the reactor from such conditions, a peaking factor, FΔH, which define as maximum integrated rod power divided by average integrated rod power, is needed [13]. Every nuclear power plant designer company will have its FΔH peaking factor safety limit in accordance with their DNB experiments and calculation.

Reactor Core Model
A reactor core model is needed to design a new loading pattern optimization. In this research, a basic PWR core is chosen as a model. This reactor is consisting of 157 fuel assemblies. However, only one-eighth of the core will be used in the optimization process, and fuel assemblies, outside this one-eighth region, will be treated symmetrically and mirrored. Although the optimization process only modifies the position of fuel assemblies in the one-eighth region of the core, calculation by the neutronic code will be done in the quarter region of the core. Fig. 1 shows the illustration of the model core and the position of the mentioned one-eighth region in the core. In the early phase of refueling loading pattern optimization design, individual fuel assembly data for at least one cycle is needed. Old fuel or fuel assemblies that have already been used in the previous cycle can be differentiated from one another by looking at its individual k∞ and burnup value in each of its corners. Fresh fuel with the arrangement and enrichment of the same rods will have a similar k∞ value and zero burnups. However, different fuel rod arrangements, enrichment, and burnable poison composition will be reflected in its different k∞ values. Fuel assemblies k∞ and burn up in corners data, which are used in this research are shown in Fig.2.

COCO: Neutronic 3D Reactor Core Calculation Code
COCO is a 3D reactor core calculation code developed by China Nuclear Power Technology Research Institute (CNPRI), a subsidiary of China General Nuclear (CGN) Company. This code is simulating one-quarter of the reactor core and give reflective boundary condition for sides that have contact with another part of the reactor in the real core. COCO is able to calculate steadystate reactor condition parameters, including the effective multiplication factor (keff), the power distribution map, and the flux map. This code could also be used in in-core fuel management, since it is capable to calculate rod-by-rod power distribution and predict precisely safety parameters such as hot channel factor (FQ), and enthalpy rise hot channel factor (FΔH).
COCO Calculation algorithm starts with collecting input data from assembly calculation. After input data is collected, the cross-section reconstruction module will use interpolation processing to build the input cross-section tables. The nodal solver module will use the tables in the calculation of various reactor core parameters. The output from the nodal solver can be used by other modules such as the control rod insertion module, thermal-hydraulic feedback module, rod-by-rod power reconstruction module, and burnup calculation module [14].
The utilization example of COCO code is in the safety analysis of the Hualong Pressurized Reactor 1000 (HPR1000) project in the United Kingdom. In the pre-construction safety report of this project, the PCM code package, which includes COCO and PINE code is utilized for safety analysis. PINE can generate two-group parameter tables for macroscopic cross-sections and the assembly discontinuity factors, while COCO uses it to calculate these parameters. The PCM nuclear design code package has been intensively validated against Nuclear Power Plant (NPP) data, experimental data, and benchmarks [15].

General Rules
The first rule for this optimization method is that only old fuel assemblies (FAs), or FAs that have already been used in the previous cycle, are allowed to be swapped with another old FAs. The position of new FAs or fresh FAs will follow the new FAs position in the reference. Optimizing new FAs positions using this method is currently not possible, and needs further research. Besides this rule, new FAs will also not be placed in the core-periphery, to ensure long cycle length.
The second rule is that FA positions will be divided into two symmetries: 1/4 symmetry, and 1/8 symmetry. FAs that are included in one symmetry cannot be swapped with fuel assemblies from another symmetry. Symmetry definition is needed because only one eighth reactor region will be modified, and FAs position in other regions will follow the position move from this one eighth reactor region according to its position symmetry. Moreover, FA chosen for the reactor core center will not be moved in all optimization steps in this method. The illustration of the FAs position symmetry definition is shown in Fig.3.
The last general rule that will be applied in the optimization method is about FAs rotation. The rule is FAs in 1/8 can be rotated 90 o , 180 o , and 270 o , but FAs in 1/4 symmetry can only be rotated 180 o . Old FAs are already irradiated in the previous cycle, but due to their position in the reactor are not irradiated evenly. Every fuel rod in the assembly will have a different irradiation level depending on its position toward the core center. Fuel rods that have a position closer to the core center will interact with higher neutron flux than other fuel rods with a slightly farther position. The difference in irradiation level leads to a difference in burnup value of the fuel rods, and the distance between a single fuel rod with reactor core center will be proportional with neutron flux it received. However, rotating FAs directly during the refueling process is not possible. To achieve the same desired FA rotation, FAs in the modifiable one-eighth region can be filled with FAs from other regions with the same symmetry position as illustrated in Fig.4.   Fig. 3. Illustration of the FAs position symmetry. Red, and green squares, respectively, represent 1/4, and 1/8 symmetry FAs in the modifiable region. Light red and light green squares, respectively, represent 1/4, and 1/8 symmetry FAs in mirrored regions. The white square represents FA in the reactor core center.

Step 1: Fuel Placement / Initializing
Fuel placement begins with the selection of a FA for the core center. FA in the core center is not allowed to be a new fuel to avoid excessive power peaking in the center. FA for core center should have enough k∞ value, more than 1.01, and moderate burnup, less than 30,000 MWd/T. After FA in the core center is in place, new FAs will be placed next in the position according to the selected reference. Before old FAs placement, it needs to be sorted according to its k∞ value. Then, old FAs with the lowest k∞ in each symmetry will be placed in the reactor periphery zone, and old FAs with higher k∞ will be gradually placed in a position closer to the reactor core center. After all the FAs are in their place, old FAs will be rotated, so that its lowest burnup fuel rods will face the center of the reactor core. The placement of a high k∞ value in position near the core center and the old FAs rotation is expected to have resulted in the longest cycle length possible with the selected FAs for the current cycle. The final result is expected to not be far less than the longest cycle length obtained in this initialization step but with an acceptable peaking factor FΔH, which is less than 1.46. The flowchart of the initialization step or step 1 is shown in Fig.5.

Step 2: Swapping the Neighbor of the Maximum FΔH FA
After the initialization step, the loading pattern will be calculated using neutronic code to determine the value of FΔH and power for every FA. The location and value of FA with the maximum FΔH in the core will be identified and recorded. The location and value of FA with the highest power will also be identified and recorded. Then, the location of the swappable FA with the highest FΔH in the maximum FΔH FA's neighbor will be identified and recorded. The next FA identification and data record are for one swappable FA with the lowest FΔH, and with the same symmetry as FA with the highest FΔH the maximum FΔH FA's neighbor. Swappable FAs in this step are defined as old FAs, both 1/4 and 1/8 symmetry, that are not in the reactor core-periphery zone. The maximum FΔH FA's neighbor is defined as FAs that have direct contact with the FA that has maximum FΔH value, and are still in the one-eighth reactor core modifiable region. How step 2 modifies FA positions is illustrated in Fig.6.

Step 3: Swapping the Lowest FΔH Within Defined Zone
When the maximum FΔH reaches a value less than desired FΔH criterion value plus 0.02, step 2 sequences will still be continued. Then, if the maximum FΔH value increases again into more than desired FΔH criterion value plus 0.02, the loading pattern with the lowest maximum FΔH value will be taken from the memory as the initial pattern for step 3.
Step 3 is swapping FAs within a defined zone. The defined zone and illustration of step 3 example on how it modifies the loading pattern is shown in Fig.7.

Step 4: Rotating the Lowest and Highest FΔH FA Neighbors
Step 4 will start with the identification and data recording of the old FA with the lowest FΔH value in the swappable zone. The Swappable zone in step 4 is identical to the swappable zone in step 2. After being identified and recorded, the old FA with the lowest FΔH value is rotated, with the mechanism mentioned in the general rules, so that the lower burnup corner is facing the FA with the lowest FΔH value in its neighbor. A newly modified pattern then will be calculated with the neutronic code. If the old FA with the lowest FΔH value is changed, the sequence above will be repeated.
The next sequence in step 4 is the identification and data recording of the FA with the maximum FΔH value. Then, the old FA neighbors of the FA with the maximum FΔH value will be rotated, so that the higher burnup corner is facing the maximum FΔH value. This sequence will be repeated if the FA with the maximum FΔH value is changed, but if the maximum FΔH value becomes more than desired FΔH criterion + 0.01 again, the sequence will be undone to the pattern with the lowest maximum FΔH then the method will continue to step 5. The illustration of an example of how the FA position was modified in this step is shown in Fig.8.

Step 5: Swapping and Rotating Lowest FΔH FA
Step 5 starts with the identification and data recording of the old FA with the lowest FΔH value in the swappable zone. The important parameter for the step 5 sequence is the k∞ of the FA. The old FA, with slightly higher k∞ than the old FA with the lowest FΔH value in its symmetry and all zone (except the center core zone) will be identified. After that, two identified FAs above will be swapped, then rotated so that its highest burnup corner will face the FA neighbor with the highest FΔH value in its new position. The newly modified pattern then will be calculated using neutronic code. If the maximum FΔH value is less than desired FΔH criterion plus 0.01, the sequence will be repeated. If the maximum FΔH becomes more than the FΔH criterion plus 0.01, the sequence will be undone to the saved loading pattern with the lowest maximum FΔH. If the maximum FΔH value is less than desired FΔH criterion minus 0.01, the method will also continue to the next step. Fig. 9 shows the illustration of an example of how this step modifies FA positions. Step 6 will only be used if FAs in the core-periphery zone have been moved inward in the previous steps, thus it is an optional step. This step starts with the identification and data recording of the old FA with the highest k∞ value in the periphery zone. Then, the old FA with the highest k∞ value in the swappable zone, but has a slightly lower k∞ than the previously identified FA in the periphery zone will be identified and its data will be recorded. The two identified FAs will be swapped. The modified pattern then will be calculated using neutronic code. If the maximum FΔH value is still below the criterion, the sequence will be repeated. The illustration of an example of how step 6 modify FA position is shown in Fig.10.

Results and Discussion
All of the steps in the optimization method have been written as an automated computer program using Python language. The whole steps from initialization to the last step have also been tested using a reactor core model as described above It took in total 25 iterations plus one initialization calculation for this method to produce the result. There are four unique loading patterns with FΔH below the FΔH criterion, which is 1.46, that have been found in this method. Fig.  11 will summarize the optimization result by displaying the maximum FΔH value, and cycle length for every loading pattern calculated in this method. After all of the steps in this method are finished, the code will eliminate all loading patterns from the memory that have the maximum FΔH value above the defined criterion. Then, the code will sort the remaining loading pattern based on its maximum FΔH value. Table 1 summarizes the maximum FΔH value, and cycle length from the four loading patterns in the final result.

Conclusion
A new method for automated computation programs for nuclear fuel reloading pattern generation in PWR has been developed. This method consists of six steps, and tries to minimize the maximum peaking factor FΔH value, while at the same time maximize the reactor cycle length. This method is written in Python language and could operate automatically by modifying input and reading results from a 3D nuclear reactor core simulation: COCO. All of the steps were executed successfully with a reactor loading pattern model. The result of this method is all of the loading patterns found in all steps which have the maximum peaking factor FΔH value below the defined criterion value. The loading patterns from this method could be the reference for the fuel reloading design decision.
Power Technology Research Institute Co., Ltd., whom without their help, support, and tools, this research could not be completed. We would also like to thank the Department of Engineering Physics, Tsinghua University for their help to make this research better.