Modification of Neutron Kinetic Code for Plate Type Fuel Nuclear Reactor

Abstract

The research is conducted on the modification of neutron kinetic code for the plate type fuel nuclear reactor. REMARK is a neutron kinetic code that works only for the cylindrical type fuel nuclear reactor. In this research, our main emphasis is on the modification of this code in order to be applicable for the plate type fuel nuclear reactor. For this purpose, detailed mathematical studies have been performed and are subjected to write the program in Fortran language. Since REMARK code is written in Fortran language, so we have developed the program in Fortran and then inserted it into the source library of the code. The main emphasis is on the modification of subroutine in the source library of the code for hexagonal fuel assemblies with plate type fuel elements in it. The number of steps involved in the modification of the code has been included in the paper. The verification studies were performed by considering the small modular reactor with hexagonal assemblies and plate type fuel in it to find out the power distribution of the reactor core. The purpose of the research is to make the code work for the hexagonal fuel assemblies with plate type fuel element.

1. Introduction

Real-time multigroup advanced reactor kinetics (REMARK) code is developed by the GSE power system used for the real-time simulation of nuclear reactor core. This code utilizes two-group, 3D time-dependent diffusion theory in the form of finite difference equations to simulate nuclear reactor core. This code has two neutron energy groups to accurately simulate the characteristics of fast and thermal neutron under normal and abnormal operating conditions of the nuclear reactor. The core geometry can be modeled as 3D mesh cell structure under the limited capacity of computer resources for real-time applications. Mesh cell sizes are chosen to ensure that the material properties within cell are homogenous or heterogeneous depending on the fuel, reflector, shielding, and so forth. In the REMARK code, there is improved quasistatic solution for obtaining the flux distribution and six delayed neutrons to be calculated in each mesh cell. Reactivity feedback is based on the core thermal hydraulic conditions like fuel temperature, coolant temperature, moderator density, and void fraction. This code can calculate the reactivity effect of control rod movement and soluble boron. It can also give the exact neutron flux readings at in-core and ex-core detector points and utilizes neutron cross-section for neutronic calculations [1, 2]. This is 3D neutron kinetic code, but it cannot examine the special kind of geometry such as hexagonal fuel assemblies. Therefore, the current research is focused on the modification of the code to make it applicable for the hexagonal fuel assemblies. The flow diagram of REMARK code is given in Figure 2.

The main concern of our research is to modify the nuclear reactor core for hexagonal fuel assemblies. The research on the new developed hexagonal fuel assembly design has already been done on the supercooled water reactor (SCWR) to show the moderation capability for the hexagonal fuel assemblies [3, 4].

The benchmark studies of different calculation schemes related to reactor design, collision probability, power reconstruction, and so forth were investigated by ten different institutes. The focus was to find out the difficulties caused by different neutron data reduction methods and for the homogenization of the fuel assemblies [5, 6].

2. Reactor Core Simulation Model

For the hexagonal fuel assemblies we consider the hexagon 𝐴𝐵𝐶𝐷𝐸𝐹 with sides “𝑎” and center “O” as shown in Figure 1. The distance between opposite sides is denoted by “Δ𝑥” and the height is denoted by “Δ𝑧” [3]. The area and volume of the hexagon can be written as 𝑆_{𝐴𝐵𝐶𝐷𝐸𝐹}=√32(Δ𝑥)²,𝑉_{𝐴𝐵𝐶𝐷𝐸𝐹}√=𝑆×Δ𝑧=32(Δ𝑥)²Δ𝑧. (1) The one-sided area of hexagon can be justified as 𝑆₁=𝑎Δ𝑧=Δ𝑥√3√×Δ𝑧=33Δ𝑥Δ𝑧. (2) We have considered eight sections of adjacent block of hexagon and calculated the neutron diffusion equations. Therefore the leakage of the neutron diffusion equation integral with volume of the node is given as _{(𝑖,𝑗,𝑘)}∇⋅𝐷∇𝜙𝑑𝑉=_𝑠𝐷𝜕𝜙𝜕𝑛𝑑𝑆=_𝑘𝑆𝑘𝐷𝜕𝜙=(),𝜕𝑛𝑑𝑆𝑇11+𝑇12+𝑇13+𝑇14+𝑇15+𝑇16+𝑇17+𝑇18 (3) where 𝑇11 can be calculated as 𝑇11=𝐷_{𝑖+1,𝑗,𝑘}𝜙_{𝑖+1,𝑗,𝑘}−𝜙_{𝑖+1/2,𝑗,𝑘}√Δ𝑥/23Δ𝑦Δ𝑧3, (4) where 𝜙_{𝑖+1/2,𝑗,𝑘} is neutrons flux density at the interface.

Figure 1 

Hexagonal plan.

Figure 2 

Flow diagram of REMARK code with modifications.

Since the net flux of neutrons at the adjacent interface is continuous, so 𝐷_{𝑖+1/2,𝑗,𝑘}𝜙_{𝑖+1/2,𝑗,𝑘}−𝜙_{𝑖,𝑗,𝑘}Δ𝑥/2=𝐷_{𝑖+1,𝑗,𝑘}𝜙_{𝑖+1,𝑗,𝑘}−𝜙_{𝑖+1/2,𝑗,𝑘},𝜙Δ𝑥/2_{𝑖+1/2,𝑗,𝑘}=𝐷_{𝑖,𝑗,𝑘}𝜙_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘}𝜙_{𝑖,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖+1,𝑗,𝑘}. (5) Substituting (5) in (4) we get 𝑇11=2𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖+1,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖+1,𝑗,𝑘}𝜙Δ𝑥^{𝑖+1,𝑗,𝑘}−𝜙_{𝑖,𝑗,𝑘}√3Δ𝑦Δ𝑧3,𝑇11𝑉=2𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖+1,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖+1,𝑗,𝑘}𝜙Δ𝑥^{𝑖+1,𝑗,𝑘}−𝜙_{𝑖,𝑗,𝑘}×√3Δ𝑦Δ𝑧/3√=3/2Δ𝑥Δ𝑦Δ𝑧4𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖+1,𝑗,𝑘}3𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖+1,𝑗,𝑘}Δ𝑥²𝜙_{𝑖+1,𝑗,𝑘}−𝜙_{𝑖,𝑗,𝑘}. (6) Similarly 𝑇12, 𝑇13, 𝑇14, 𝑇15, and 𝑇16 can be calculated as 𝑇11 and have the same value as 𝑇11. 𝑇17 is an integral on the top, whereas 𝑇18 is an integral on the bottom and can be written as𝑇17=𝐷_{𝑖,𝑗,𝑘}𝜙_{𝑖,𝑗,𝑘+1/2}−𝜙_{𝑖,𝑗,𝑘}√Δ𝑧/232Δ𝑥Δ𝑦, (7) where 𝜙_{𝑖,𝑗,𝑘+1/2} is neutrons flux density at the interface.

Here, the net neutrons flux at the adjacent interface is continuous; therefore 𝐷_{𝑖,𝑗,𝑘}𝜙_{𝑖,𝑗,𝑘+1/2}−𝜙_{𝑖,𝑗,𝑘}Δ𝑧/2=𝐷_{𝑖,𝑗,𝑘}𝜙_{𝑖,𝑗,𝑘+1}−𝜙_{𝑖,𝑗,𝑘+1/2},𝜙Δ𝑧/2_{𝑖,𝑗,𝑘+1/2}=𝐷_{𝑖,𝑗,𝑘}𝜙_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘+1}𝜙_{𝑖,𝑗,𝑘+1}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘+1}. (8) Substituting (8) in (7) we get 𝑇17=2𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘+1}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘+1}𝜙Δ𝑧^{𝑖,𝑗,𝑘+1}−𝜙_{𝑖,𝑗,𝑘}×√32Δ𝑥Δ𝑦,𝑇17𝑉=2𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘+1}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘+1}𝜙Δ𝑧^{𝑖,𝑗,𝑘+1}−𝜙_{𝑖,𝑗,𝑘}×√3/2Δ𝑥Δ𝑦√,3/2Δ𝑥Δ𝑦Δ𝑧𝑇17𝑉=2𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘+1}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘+1}Δ𝑧²𝜙_{𝑖,𝑗,𝑘+1}−𝜙_{𝑖,𝑗,𝑘}, (9) where 𝑇18 will have the same value as 𝑇17.

Now we will calculate the extrapolated length, that is, 𝑑=Δ𝑥_extd, to be considered as core boundary conditions: 𝐽=−𝐷_{𝑖,𝑗,𝑘}𝜙_{𝑖,𝑗,𝑘}=𝑑+Δ𝑥/2−2𝐷_{𝑖,𝑗,𝑘}𝜙_{𝑖,𝑗,𝑘},2𝑑+Δ𝑥_{𝑖,𝑗,𝑘}∇⋅𝐷∇𝜙𝑑𝑉=_𝑆𝐷𝜕𝜙𝜕𝑛𝑑𝑆=_𝑘𝑆𝑘𝐷𝜕𝜙𝜕𝑛𝑑𝑆,𝑇1=−𝐽∗𝑆𝑉=2𝐷_{𝑖,𝑗,𝑘}𝜙_{𝑖,𝑗,𝑘}×√2𝑑+Δ𝑥3Δ𝑦Δ𝑧/3√=3/2Δ𝑥Δ𝑦Δ𝑧4𝐷_{𝑖,𝑗,𝑘}3Δ𝑥²+2Δ𝑥Δ𝑥_extd𝜙_{𝑖,𝑗,𝑘}. (10) In three-dimensional hexagonal grid, 𝐷 is called diffusion coefficient and ∇⋅𝐷∇Ψ is the leakage term in four directions. The leakage term is calculated in four directions as 1∇⋅𝐷∇𝜙=𝑉_{𝑖,𝑗,𝑘}=()∇⋅𝐷∇𝜙𝑑𝑉𝑇11+𝑇12+𝑇13+𝑇14+𝑇15+𝑇16+𝑇17+𝑇18×𝑉⁻¹,∇⋅𝐷∇𝜙=𝑎⁻_𝑖𝜙_{𝑖−1,𝑗,𝑘}+𝑎⁺_𝑖𝜙_{𝑖+1,𝑗,𝑘}+𝑎⁻_𝑗𝜙_{𝑖,𝑗−1,𝑘}+𝑎⁺_𝑗𝜙_{𝑖,𝑗+1,𝑘}+𝑎⁻_𝑗2𝜙_{𝑖−1,𝑗−1,𝑘}+𝑎⁺_𝑗2𝜙_{𝑖+1,𝑗+1,𝑘}+𝑎⁻_𝑘𝜙_{𝑖−1,𝑗,𝑘}+𝑎⁺_𝑘𝜙_{𝑖+1,𝑗,𝑘}−𝑎_𝑖+𝑎_𝑗+𝑎_𝑗2+𝑎_𝑧𝜙_{𝑖,𝑗,𝑘}. (11) For the middle section block grid, 𝑎⁺_𝑖=4𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖+1,𝑗,𝑘}3𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖+1,𝑗,𝑘}Δ𝑥²,𝑎⁻_𝑖=4𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖−1,𝑗,𝑘}3𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖−1,𝑗,𝑘}Δ𝑥²,𝑎⁺_𝑘=2𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘+1}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘+1}Δ𝑧²,𝑎⁻_𝑘=2𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘−1}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘−1}Δ𝑧². (12) Similarly, 𝑎⁺_𝑗, 𝑎⁻_𝑗, 𝑎⁺_𝑗2, and 𝑎⁻_𝑗2 will have the same values.

For the grid block boundary section: 𝑋 direction left border 𝑎⁺_𝑖=4𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖+1,𝑗,𝑘}3𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖+1,𝑗,𝑘}Δ𝑥²,𝑎⁻_𝑖=4𝐷_{𝑖,𝑗,𝑘}3Δ𝑥²+2Δ𝑥Δ𝑥_extd; (13) 𝑋 direction right border 𝑎⁺_𝑖=4𝐷_{𝑖,𝑗,𝑘}3Δ𝑥²+2Δ𝑥Δ𝑥_extd,𝑎⁻_𝑖=4𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖−1,𝑗,𝑘}3𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖−1,𝑗,𝑘}Δ𝑥²; (14) 𝑍 direction border at the bottom 𝑎⁺_𝑘=2𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘+1}3𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘+1}Δ𝑧²,𝑎⁻_𝑘=2𝐷_{𝑖,𝑗,𝑘}Δ𝑧²+2Δ𝑧Δ𝑧_extd; (15) 𝑍 direction border at the top𝑎⁺_𝑘=2𝐷_{𝑖,𝑗,𝑘}Δ𝑧²+2Δ𝑧Δ𝑧_extd,𝑎⁻_𝑘=2𝐷_{𝑖,𝑗,𝑘}𝐷_{𝑖,𝑗,𝑘−1}𝐷_{𝑖,𝑗,𝑘}+𝐷_{𝑖,𝑗,𝑘−1}Δ𝑧². (16)

Similarly, 𝑎⁺_𝑗, 𝑎⁻_𝑗, 𝑎⁺_𝑗2,  and  𝑎⁻_𝑗2 will have the same values: 𝑎_𝑖=𝑎⁺_𝑖+𝑎⁻_𝑖,𝑎_𝑗=𝑎⁺_𝑗+𝑎⁻_𝑗,𝑎_𝑗/2=𝑎⁺_𝑗/2+𝑎⁻_𝑗/2,𝑎_𝑧=𝑎⁺_𝑧+𝑎⁻_𝑧. (17) In order to solve the diffusion and convection diffusion problems, we used alternating direction implicit method. This method was first discovered by Peaceman and Rechford [4] especially to solve parabolic and 2D elliptical problems.

The reactor core is a typical case of nonuniform problem and the parameters involved are not continuous. Therefore, ADI method was considered for calculating the neutron flux of each section in the block of reactor core. Neutron flux in whole reactor core is not continuous and there is a need to divide the coarse grid of core assembly. Hence, ADI method was employed to calculate the data and to give the accuracy of results, and, consistently, an algorithm should be adjusted hereafter [5].

3. Modification of REMARK Code

In real-time simulation, the computational time of each step is less than the simulation time, so it is difficult that the nodal expansion method, the nodal green function method, and other methods meet the requirements for the real-time simulation of the reactor core.

In REMARK code, neutron diffusion time-space dynamics equations with two-group and six-group delayed neutrons are utilized for the simulation [5]. These equations are given as below.

Fast Group with 𝐸≥1. We have 𝜕𝜙₁(𝑟,𝑡)𝑣₁𝜕𝑡=∇⋅𝐷₁(𝑟,𝑡)∇𝜙(𝑟,𝑡)−_𝑎1(𝑟,𝑡)𝜙₁(𝑟,𝑡)−₁₂(𝑟,𝑡)𝜙₁(𝑟,𝑡)+1−𝛽𝑣_𝑓1(𝑟,𝑡)𝜙₁()+𝑣𝑟,𝑡1−𝛽×_𝑓2(𝑟,𝑡)𝜙₂(𝑟,𝑡)+𝑆_𝑑(𝑟,𝑡)+𝑆(𝑟,𝑡). (18)

Thermal Group with E < 1. We have 𝜕𝜙₂(𝑟,𝑡)𝑣₂𝜕𝑡=∇⋅𝐷₂(𝑟,𝑡)∇𝜙₂(𝑟,𝑡)−₁₂(𝑟,𝑡)𝜙₁(𝑟,𝑡). (19) The six-group delayed neutron precursor concentration can be calculated as 𝜕𝐶_{𝑘(𝑟,𝑡)}𝜕𝑡=−𝜆_𝑘𝐶_𝑘(𝑟,𝑡)+𝛽_𝑘𝑣_𝑓1(𝑟,𝑡)𝜙₁(𝑟,𝑡)+𝛽_𝑘𝑣_𝑓2(𝑟,𝑡)𝜙₂(𝑟,𝑡). (20) In our research, we have experienced modelling of the real-time simulation technique for hexagonal fuel assemblies especially account to solve neutron diffusion equation. In hexagonal assembly, there are eight leakage directions which change the diffusion terms in diffusion equations. For the calculation of the reactor core, there are several steps that should be taken.(i)Grid homogenization (fuel/control grid in assemblies so as to make 1D model which consists of fuel, clad, and moderator) should be performed.(ii)Results from grid homogenization are used to calculate few-group or two-group diffusion calculations.(iii)The finally originating group calculations are used in diffusion equation to solve the diffusion calculation of the whole reactor core.(iv)Effective multiplication factor and power or neutron flux distribution will be calculated.(v)Neutron transport equations will be employed after a series of approximations to get the diffusion calculations.

3.1. Assumptions

There are some assumptions that need to be considered for modification in the REMARK code in terms of hexagonal assemblies.(i)Ignore transfer of neutrons from thermal to fast group; that is, transfer cross-section is zero.(ii)Delayed or prompt neutron should be considered as fast neutron.(iii)Poison concentrations like Xe and Sm from I-135 and Pm-149 should be considered in calculations.(iv)Taking into consideration other phenomena, fission power and decay power can be calculated.(v)Now the effective multiplication factor will be calculated by neutrons originating from fission divided by neutrons from absorption and leakage.In this paper, we have compiled all the data in separate program of hexagonal fuel assemblies into many subroutines like:(i)neutron flux distribution,(ii)group 1 and group 2 leakage,(iii)delayed neutrons,(iv)discontinuity factor,(v)method to solve diffusion equation.Since REMARK code was written in Fortran language, so we developed the program in Fortran language against each factor stated above (Figure 3). These programs are then added in the source code of REMARK code and they compile the whole subroutines. The code was then tested by considering simulation of hexagonal fuel assemblies and it gave good results, thus pointing to the accuracy of modification of REMARK code for hexagonal fuel assemblies.

Figure 3 

Modified REMARK code in Fortran language.

4. Verification Studies of REMARK Code

For verification of REMARK code, we have considered small modular nuclear reactor of power 220 MW designed by the College of Nuclear Science and Technology, Harbin Engineering University, China [6, 7]. The main design parameters of the reactor core are illustrated in Table 1.

Table 1 

Design parameter of reactor core.

For the verification studies, we have considered some of the parameters of the reactor core design and utilized in the lattice physic code, that is, HELIOS [8]. This code gives neutron cross-sectional data which is then used in the REMARK code for obtaining the power distribution of the reactor core. In the current research, the reactor core has hexagonal fuel assemblies with plate type fuel in it. The REMARK code in its original form cannot work for the plate type fuel nuclear reactor; therefore, the modification studies have been done in the REMARK code in order to be applicable for plate type fuel reactor. The cross-sectional data from HELIOS code exchanges with modified version of REMARK code which gives the power distribution of the reactor core and this verifies the modification of the code, although REMARK code needs some thermal hydraulic data which come from the plant design data and are used to get power distribution of the reactor core. The flow diagram of the process is shown in Figure 4.

Figure 4 

Flow diagram for verification studies.

In order to generate cross-sectional data, HELIOS code has been used which is two-dimensional (2D) lattice physics high-order transport code and is able to solve any kind of geometry. This code utilizes neutron and gamma groups for the burnup calculations, flux distribution, and micro/macro cross-sectional data. It works via two separate codes AURORA and ZENITH which are input and output processing codes. The data exchanged between the two codes are accessed by the subroutine package HERMIS [8].

The obtained power distribution of the reactor core from combination of HELIOS and modified version of the REMARK code is shown in Figure 5.

Figure 5 

Axial power distribution of the reactor core.

5. Conclusion

In this paper, the research is focused on the modification and verification of the neutron kinetics code, that is, REMARK code. Originally this code only works for the cylindrical type fuel elements, so we have modified this code to be applicable to plate type fuel nuclear reactor. The modification has been done by detailed mathematics and the modified files in the source library of the code have been computed. The modified code has been verified by coupling with HELIOS code to get the power distribution of the plate type fuel nuclear reactor. The result obtained from coupling analysis gives the accuracy of the modified version of the code.

Acknowledgment

The authors would like to extend their sincere appreciation to the Deanship of Scientific Research at King Saud University for its funding of this research through the Research Group Project no. RGP-VPP-255.