Calculation of effective point reactor kinetic parameters combined with the Galerkin finite element method

Bao-Xin Yuan· Wan-Kui Yang · Song-Bao Zhang · Bin Zhong ·Jun-Xia Wei · Yang-Jun Ying

Abstract Precise point reactor kinetic parameters are essential in the study of reactor dynamics. Point reactor kinetic parameters include the effective delayed neutron fraction and prompt neutron lifetime. In this work, effective point reactor kinetic parameters, which can be applied to unstructured grids, were calculated based on the Galerkin finite element method. First, two-dimensional and three-dimensional benchmarks were used to verify the calculation of steady-state neutronic parameters. Then, the Tehran research reactor core was divided into hexahedral meshes, and the forward flux and adjoint flux were calculated.Finally,the effective point reactor kinetic parameters of the Tehran research reactor were obtained by comprehensively processing the steady-state neutronic parameters.

Keywords Finite element method · Reactor · Kinetic parameters

1 Introduction

The finite element method can handle unstructured meshes and can thus be adapted to complex geometries.Since the 1970s, this method has been gradually applied and promoted in the field of neutronics computing [1–4].

Zhang et al.[5]proposed a strategy to accelerate solving the neutron diffusion equation via the finite element method, and the acceleration ratio reached 103. Hosseini et al. [6–9] proposed a finite element method for solving steady-state neutronics and neutron noise problems.Avvakumov et al. [10] also solved a typical neutron dynamics problem via the finite element method.The finite element platform FEniCS was used to develop the neutron dynamics code,and the extensible tool package SLEPc was used to solve the eigenvalue of the coefficient matrix.Gonza´lez-Pintor et al. [11] proposed a high-order finite element method to solve the lambda eigenvalue problem of hexagonal geometry.

For analyzing reactor dynamics, effective point reactor kinetic parameters need to be obtained in advance, and finite element methods can be used.

Saadi et al. [12] reviewed the calculation methods of effective point reactor kinetic parameters and compared deterministic and stochastic methods by using the same benchmark.Arkani et al.[13]calculated the effective point reactor kinetic parameters of the Tehran research reactor(TRR) by using the six-group weighting method. Lashkari et al. [14] calculated the effective point reactor kinetic parameters of the TRR mixed core through the core analysis tool MTR_PC. The effective delayed neutron fraction was found to decrease with burnup, while the prompt neutron lifetime increased.With increasing highly enriched uranium fuel assemblies, the prompt neutron lifetime increases,while the effective delayed neutron fraction does not change significantly.Hosseini et al.[15]calculated and measured the effective point reactor kinetic parameters of the TRR. MTR_PC was used to calculate the neutronic parameters of the TRR in cold and hot states. By using WIMS and CITATION, Zaker et al. [16] calculated the effective point reactor kinetic parameters of two TRR cores with both highly enriched uranium and low enriched uranium.

To improve the steady-state calculation in the current deterministic method of effective point reactor kinetic parameters, the calculation of these parameters and the Galerkin finite element method are combined in this work based on the existing framework of FEMN[17].This paper is organized as follows: the basic concepts of the Galerkin finite element method are outlined, and the 2D and 3D steady-state benchmarks based on the Galerkin finite element method are verified in Sect. 2. In Sect. 3, the TRR benchmarks used worldwide are described, and the effective point reactor kinetic parameters of the TRR are obtained by the improved strategy proposed above. In the final section, the calculation results are discussed, and a brief summary is provided.

2 Galerkin finite element method

2.1 Galerkin finite element formulation of the neutron diffusion problem

Table 1 Macroscopic cross sections of the TWIGL 2D benchmark

Table 2 k*eff of TWIGL

After integration by parts,the finite element form of the 3D steady-state neutronics equation can be obtained:

Table 5 Macroscopic cross sections of the IAEA 3D benchmark

Table 6 keff and k*eff of IAEA-3D PWR

The 3D steady-state adjoint equation of neutronics is as follows:

Table 7 Power distribution of the 3D benchmark obtained by FEMN

Similarly,the finite element form of the 3D steady-state neutronics adjoint equation can be obtained as follows:

2.2 Benchmark test

1. 2D steady-state neutronics benchmark.

The 2D benchmark TWIGL is used for the code test[18]. The geometric layout of the 1/4 core is shown in Fig. 1,and the macroscopic cross sections used are listed in Table 1. Both triangular and quadrilateral elements can be used for 2D geometry meshing, and Fig. 2 shows the 2D mesh adopted by the present study. The steady-state neutronic parameters are calculated by FEMN.k*effand adjoint flux are shown in Tables 2 and 3 and 4, respectively.

Compared with the reference value[18],the deviation of k*effis 0.000002, while the deviation of adjoint flux in the outermost assemblies is less than 4%. The deviation is mainly affected by the mesh size, and it can be further reduced by mesh refinement.

1. 3D steady-state neutronics benchmark.

The 3D benchmark IAEA-3D PWR is used for the code test [19]. A schematic of the 3D benchmark is shown inFig. 3,and the macroscopic cross sections used are listed in Table 5. Both triangular prism and hexahedral elements can be used for 3D geometric meshing.Figure 4 shows the 3D mesh adopted by the present study. The steady-state neutronic parameters are calculated by FEMN. keff, k*eff,and power distribution are shown in Tables 6 and 7.

Compared with the reference value[19],the deviation of keffis 0.00004. The deviation of heat deposition of the inner fuel assembly is less than 1%.Although the deviation of the outer fuel assembly close to the reflector is high,the overall deviation is less than 2.2%.The deviation is mainly affected by the mesh size,and it can be further reduced by mesh refinement.

The forward flux and adjoint flux can also be adequately determined by FEMN.The distribution of forward flux and adjoint flux for the 3D benchmark is shown in Figs. 5 and 6, respectively.

3 Results and discussion

3.1 Description of the TRR

The TRR is a 5-MW pool-type light-water moderated,heterogeneous solid fuel reactor,in which the water is also used for cooling and shielding.The TRR core is immersed in either section of a two-section, concrete pool filled with water. Utilization of the reactor is essential for research,training, and production of radioisotopes [12]. Substantial research has been performed on the point reactor kinetic parameters of the TRR;thus,it was chosen as a benchmark to facilitate the calculation comparison. The first core configuration of the TRR is shown in Fig. 7 [13].

The standard fuel element (SFE) consists of 19 fuel plates with 20% enrichment of uranium 235. Each fuel plate consists of two unit areas of fuel cladding and one fuel meat;the meat is made of U3O8 powder dispersed in a pure aluminum matrix, and the cladding and other structural materials are Al6061. The coolant channel is located between the fuel plates. The control fuel element (CFE)consists of 14 fuel plates with 20%enrichment of uranium 235, and each fuel plate consists of two unit areas of fuel cladding and one fuel meat.

The shim safety control rod is composed of two absorbing plates and a knuckle subassembly. The absorber plate is an alloy of silver, indium, and cadmium (80, 15,and 5%wt.,respectively),while the regulating rod is made of stainless steel. All absorbing rods are fork type. The main geometrical data of the CFE and SFE can be found in Ref. [12], and the configuration is shown in Fig. 8 [13].

The calculation of TRR’s effective point reactor kinematic parameters involves three steps: (1) Preparation of the six-group cross section and six-group neutron velocity by WIMS and 2DSN; (2) Calculation of the six-group forward flux and adjoint flux with the Galerkin finite element method;(3)Calculation of the effective point reactor kinetic parameters. Because the TRR core is asymmetric,the use of triangular prism elements for geometric meshing and six-group calculation will cause substantial memory consumption, which is impractical for personal computers.Therefore, this work uses a memory-saving hexahedron meshing method, and Fig. 9 shows the meshes used in the calculation.

To distinguish the energy spectrum between prompt and delayed neutrons, the total number of energy groups must be greater than two groups,because in a typical two-group structure, neutrons can only be produced in the fast group,so the effect of delayed neutrons is not reflected. The neutron group structure,transient neutron energy spectrum,delayed neutron energy spectrum, and delayed neutron parameters of235U and238U can be found in Ref. [12].

3.2 Results and discussion

The steady-state neutronic parameters of the TRR were calculated using the FEMN, and the excess reactivity in this state is 6344 pcm.The forward flux and adjoint flux in the axial plane of 51.5 cm are shown in Figs. 10 and 11,respectively.The forward flux and adjoint flux distribution of the TRR in this state is not given in relevant literature,so only the calculation results are given here for peer reference. For FEMN, the relative error of the steady-state calculation has been discussed, and interested readers can refer to Sect. 2. From the perspective of excess reactivity,the deviation between the calculated value(6344 pcm)and the reference value (6481 pcm) of Ref. [12] is 2%, which proves that the steady-state calculation is reliable from the side.

Table 8 Effective point reactor kinetic parameters of the TRR

Equations (1) and (2) proposed above are used to process the effective point reactor kinetic parameters of the TRR, and the calculated results are shown in Table 8.Compared with the recommended value of Ref. [12], the relative error of the effective delayed neutron fraction is 2%, while the relative error of the prompt neutron lifetime is 10%. The main reason for the relative error is that the cross sections and neutron energy spectra used by different researchers are inconsistent, which eventually leads to some slight differences in forward flux, adjoint flux, and group neutron velocity.

4 Conclusion

1. Based on the internally developed code FEMN, a method for determining effective point reactor kinetic parameters was developed and validated.

2. Compared with the reference values,the relative errors of the effective delayed neutron fraction and prompt neutron lifetime are 2% and 10%, respectively. The main factor that led to the errors was the inconsistency between the cross sections and the neutron energy spectra.