2.1 Galerkin finite element formulation of the neutron diffusion problem

The effective delayed neutron fraction is defined as follows:

For group g, ${\varnothing }_g{}^{\text{*}}$ is adjoint flux, ${\chi }_{\text{d}i}{}^q$ is the delayed neutron spectrum of fissionable isotope $q$ in group i, ${\beta }_i{}^q$ is the delayed neutron yield for fissionable isotope $q$ in group i, $\text{ν}$ is the number of neutrons emitted per fission, ${\left({\sum }_{\text{f}}\right)}_{g^{\prime }}$ is the fission cross section, ${\varnothing }_g$ is forward flux, and ${\chi }_p{}^q$ is the prompt neutron spectrum of fissionable isotope $q$.

Similarly, the definition of prompt neutron lifetime can be given as follows:

where $\frac{1}{v_g}$ is the inverse velocity of group $g$.

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

For group g, $D_{\text{g}}$ is the diffusion coefficient, ${\sum }_{r,g}$ is the removal cross section, ${\sum }_{g^{\prime }\to g}$ is the scattering cross section from group $g^{\prime }$ to $g$, ${\text{χ}}_g$ is the fission spectrum, and $k_{\text{eff}}$ is the effective multiplication factor.

The boundary conditions are described as follows:

where ${\text{Г}}_1$ is the boundary, n is the normal direction of ${\text{Г}}_1$, and $\text{α}$ and $\text{β}$ are constants.

By following the Galerkin finite element method, the weak form of the diffusion equation can be written as follows:

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

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

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 neutronics parameters are calculated by FEMN. $k_{\text{eff}}^{*}$ and adjoint flux are shown in Table 2 and 3 and 4, respectively.

Fig. 1.Full-screen View

(Color online) Schematic of the 2D benchmark

Fig. 2.Full-screen View

(Color online) Meshes used by the 2D benchmark

Compared with the reference value [18], the deviation of $k_{\text{eff}}^{\text{*}}$ is 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.

2) 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 in Fig. 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 neutronics parameters are calculated by FEMN. k_eff, $k_{\text{eff}}^{*}$, and power distribution are shown in Table 6 and 7.

Fig. 3.Full-screen View

(Color online) Schematic of the 3D benchmark. (a) Top view; (b) Side view.

Fig. 4.Full-screen View

(Color online) Meshes used by the IAEA 3D benchmark

Table 5.Full-screen View

Macroscopic cross sections of the IAEA 3D benchmark

Material types		Energy group	$D_g \left({\text{cm}}^{-1}\right)$	${\sum }_{t,g} \left({\text{cm}}^{-1}\right)$	${\left(\nu {\sum }_f\right)}_g\left({\text{cm}}^{-1}\right)$	${\sum }_{g\to g^{\text{'}}}\left({\text{cm}}^{-1}\right)$
	Fuel1	1	1.500	0.010	0.000	0.020
		2	0.400	0.085	0.135	—
	Fuel1+rod	1	1.500	0.010	0.000	0.020
		2	0.400	0.130	0.135	—
	Fuel2	1	1.500	0.010	0.000	0.020
		2	0.400	0.080	0.135	—
	Reflector	1	2.000	0.000	0.000	0.040
		2	0.300	0.010	0.000	—
	Reflector+Rod	1	2.000	0.000	0.000	0.040
		2	0.300	0.055	0.000	—

Show more

Compared with the reference value [19], the deviation of $k_{\text{eff}}$ is 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 are shown in Fig. 5 and Fig. 6, respectively.

Fig. 5.Full-screen View

(Color online) Distribution of forward flux for the 3D benchmark obtained by FEMN (Z=19 cm, 99 cm, 189 cm, 289 cm, 363 cm). (a) Fast group forward flux; (b) Thermal group forward flux

Fig. 6.Full-screen View

(Color online) Distribution of adjoint flux for the 3D benchmark obtained by FEMN (Z=19 cm, 99 cm, 189 cm, 289 cm, 363 cm). (a) Fast group adjoint flux; (b) Thermal group adjoint flux.

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].

Fig. 7.Full-screen View

(Color online) First core configuration of the TRR

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 reference [12], and the configuration is shown in Fig. 8 [13].

Fig. 8.Full-screen View

Fork-type control rod and cross-sectional view of the CFE and SFE

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.

Fig. 9.Full-screen View

(Color online) Meshes used for the TRR benchmark

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 of ²³⁵U and ²³⁸U can be found in reference [12].

3.2 Results and discussion

The steady-state neutronics 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 reference [12] is 2%, which proves that the steady-state calculation is reliable from the side.

Fig. 10.Full-screen View

(Color online) Distribution of forward flux for the TRR obtained by FEMN (Z=51.5 cm). (a) Forward flux of group 1; (b) Forward flux of group 2; (c) Forward flux of group 3; (d) Forward flux of group 4; (e) Forward flux of group 5; (f) Forward flux of group 6.

Fig. 11.Full-screen View

(Color online) Distribution of adjoint flux for the TRR obtained by FEMN (Z=51.5 cm). (a) Adjoint flux of group 1; (b) Adjoint flux of group 2; (c) Adjoint flux of group 3; (d) Adjoint flux of group 4; (e) Adjoint flux of group 5; (f) Adjoint flux of group 6.

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 reference [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.