The steady-state neutron transport equation can be written as follows: Ω ⋅ ∇ ψ r , Ω , E + Σ t r , E ψ r , Ω , E = ∫ 0 ∞ ∫ 4 π Σ s r , Ω ′ → Ω , E ′ → E ψ r , Ω ′ , E ′ d Ω ′ d E ′ + χ r , E 4 π k e f f ∫ 0 ∞ ν Σ f r , E ′ ∫ 4 π ψ r , Ω ′ , E ′ d Ω ′ d E ′ , (1) where r , Ω , E are the space, angle, and energy; ψ is the angular neutron flux (cm⁻²s⁻¹); Σ t is the total macroscopic cross section (cm⁻¹); Σ s is the scattering macroscopic cross section (cm⁻¹); χ is the fission spectrum; ν Σ f is the production macroscopic cross section (cm⁻¹); k e f f is the effective multiplication factor.

Multi-group approximation is adopted for energy, and the discrete ordinate method is adopted for angular discretization. After that, a single group neutron transport equation in each direction can be obtained. Ω m ⋅ ∇ ψ g m r + Σ t , g r ψ g m r = Q s , g m r + Q f , g r , (2) where m represents the direction index, g represents the energy group index, and Q s , Q f represent scattering source and fission source, respectively. This equation can be solved by energy group in each direction, and the different energy groups and directions can be coupled through the right hand side source term.

Ignoring m and g indexes, the transport equation for the hexagonal prism in the calculated coordinate system (Figure 1) can be written as ∑ s = x , y , z μ s h s ∂ ∂ s ψ x , y , z + Σ t ψ x , y , z = Q x , y , z , (3) where x , z ∈ − 1 / 2,1 / 2 and y ∈ − y s x , y s x , y s x = 1 − x / 3 . h x = h y = h r is the actual opposite side distance of the hexagon, and h z is actual height of the hexagonal prism.

Three-dimensional transport equations can be decoupled into three radial and one axial equation by transverse integral technique. Integral Eq. 3 within y ∈ − y s x , y s x and z ∈ − 1 / 2,1 / 2 , and the transverse integral equation along x direction can be obtained: d d x y s x ψ x x + Σ t h r μ x y s x ψ x x = q x x . (4)

Let ε x ≡ e − Σ t h r / μ x solve the ordinary differential Eq. 4, and the nodal radial integral equation can be obtained. y s x ψ x x = 3 6 ψ x − ε x x + 1 / 2 + ∫ − 1 / 2 x q x ′ ε x x − x ′ d x ′ . (5)

Then, the outgoing integral flux can be obtained as ψ x + = ψ x − ε x + 2 3 ∫ − 1 / 2 1 / 2 q x ′ ε x 1 / 2 − x ′ d x ′ . (6)

Similar equations can be obtained for the direction of u, v, and z.

The transverse integral flux, source, and axial leakage term can be expanded by orthogonal polynomials with weights as ψ x x = ∑ i = 0 3 ψ x i h i x , (7) Q x x = ∑ i = 0 3 Q x i h i x . (8)

The first three order polynomials are defined as follows: h i x = 1 ,   sgn μ x x , x 2 − 5 / 72 , sgn μ x x 3 − 7 / 50 x . (9)

We insert Eq. 7 and Eqs 8–5, then multiply both sides of the equation by the basis function, and integrate over the nodal, so the nodal response relationship can be obtained. We insert Eq. 7 and Eqs 8–6, and the relationship of the outgoing surface flux with the incoming surface flux can be obtained.

In hexagonal prism geometry, for any given direction of the quadrature with 60° symmetry, there are three radial planes, one axial plane as the incoming plane, and the remaining planes as the outgoing planes. The response relationship of the outgoing surface flux with the nodal average flux and incoming surface flux can be expressed as ψ s + = A s ψ ¯ + 1 − B s ψ s − + C s , (10) where the coefficients and derivation process can be referred to Wang et al. (2020. Integrating Eq. 3 within a nodal, the neutron balance equation can be obtained. ∑ s = v , x , u , z α s μ s ψ s + − ψ s − + Σ t ψ ¯ = Q , ¯ (11) where α v = α x = α u = 2 / 3 h r , α z = 1 / h z . Combining Eqs 10, 11, the nodal average flux can be calculated from the incoming surface flux: ψ ¯ = ∑ s = v , x , u , z α s μ s B s ψ s − − C s + Q ¯ Σ t + ∑ s = v , x , u , z α s μ s A s . (12)

According to the sweeping order determined in advance, the average flux of the nodal is calculated by Eq. 12 with incoming boundary condition, and then, the outgoing surface flux is calculated by Eq. 10 and taken as the incoming boundary condition of the next nodal. The flux distribution and keff are obtained by iterative of the nodal, direction, and energy group, and then, the power distribution can be obtained as q r = ∑ g = 1 G Σ f − k a p a , g r ϕ g r , (13) where q is the power distribution, Σ f − k a p a is the energy release cross section, and ϕ is the scalar flux distribution. Once q is obtained, the TH calculation can be performed.

In the steady state and ignoring the axial heat transfer, the heat transfer equation and boundary condition in the cylindrical coordinate system are as follows: 1 ρ C p 1 r d d r r k d T d r + 1 ρ C p Q = 0 , (25) d T d r = 0 , r = 0 ; T = T c l a d _ s , r = R p i n , (26) where ρ is the material density (kg·m⁻³); C p is the specific heat capacity (J·kg⁻¹K.⁻¹); k is the heat conductivity (W·m⁻²K.⁻¹); Q is the heat source (W·m⁻³).

Integrating the equation over the space, the following equation can be obtained: r k d T d r + Q r 2 2 = 0 . (27)

The finite difference method (FDM) is used to solve Eq. (27).The fuel is divided into more than five segments, the gap is one segment, and the cladding is at least two segments (as shown in Figure 3). The radial temperature distribution of the fuel pin can be obtained by solving the linear equations generated by FDM.

The temperature-dependent cross-section library is generated by the external program, which uses 1D cylindrical geometry to equivalent the hexagonal assembly, superfine group method to do the resonance calculation, and super homogenization (SPH) method to do cross-section homogenization. Once the cross-section library of the fuel, cladding, and coolant at tabulated temperature is generated, the cross section at any temperature can be obtained by interpolation. The linear interpolation formula is as follows: σ = 1 − w σ 1 + w σ 2 , (28) where σ is the cross section at the temperature point T , σ 1 and σ 2 are the cross sections at two adjacent temperature points T 1 and T 2 T 1 < T 2 , and w is a weighting coefficient. For fuel materials, the temperature dependence of cross sections is dominated by the Doppler effect, which is proportional to the square-root of temperature; therefore, the cross sections are interpolated using the square root of temperature interpolation. Thus, w = T − T 1 T 2 − T 1 . (29)

For cladding and coolant materials, the cross sections are linearly interpolated by temperature. Thus, w = T − T 1 T 2 − T 1 . (30)

Based on the temperature distribution obtained from the TH calculation, the cross section data required by the NE calculation can be obtained.

The macroscopic cross section is related to the microscopic cross section and density of nuclides, both of which vary with temperature. The microscopic cross section can be obtained by interpolation described previously, while the variation in the nuclide density with temperature is mainly due to the variation in material density. Assuming that the material density at temperature T 0 is ρ 0 and the nuclide density is N 0 , the nuclide density at temperature T is N T = N 0 ρ T ρ 0 , (31) where ρ T is the material density at temperature T . The macroscopic cross section at temperature T can be calculated as Σ T = N T σ T . (32)

The NE/TH coupling calculation process can be divided into four parts. First, the homogeneous microscopic cross-section library of the hexagonal assembly at different state points is obtained by the assembly calculation. Then, the cross section is interpolated according to the temperature distribution (taking the design temperature for the first calculation) to obtain the cross section at the actual state point. Next, neutronics calculation is carried out to obtain the power distribution. Finally, the power distribution is set as the heat source for TH calculation and the temperature distribution of the fuel, cladding, and coolant in the reactor core is obtained. The calculation and data transfer process are shown in Figure 4. The Picard iteration method is adopted to do the NE/TH coupling calculation, which keeps iterating until the NE fields and TH fields have converged.

The thermal power of the reactor is 50 MW, the material of the fuel is 19.95% enrichment UO₂, the control-rod absorber material is B₄C, the reflector material is BeO, and the coolant material is the lead-bismuth eutectic (LBE) (44.5% Pb and 55.5% Bi). The reactor core consists of 144 fuel assemblies, 18 control-rod assemblies, and 48 reflector assemblies, each assembly consists of 61 rods, which are arranged in a triangular grid, and the configuration is shown in Figure 5. The detailed core parameters can be found in the work of Zhao et al. (2022, p. 50).

The physical properties of the materials in SLBR refer to Sobolev, (2007) and Roelofs, (2019), and the thermal conductivity of fuel material varies with temperature as k = 1.0 / 0.042 + 2.71 × 10 − 4 T + 6.9 × 10 − 11 T 3 . (33)

The thermal conductivity of cladding materials varies with temperature as k = 462.0 + 0.134 T . (34)

The melting point of the LBE coolant is 397 K and the boiling point is 1943 K. When using liquid metals as coolants and the fuel rods are arranged in triangular grids with 1.1 ≤ P / D ≤ 1.5 (where P is the pitch of fuel rods and D is the diameter of fuel rods) and 80 ≤ P e ≤ 4000 , the correlation (Mikityuk, 2009) for the Nusselt number is obtained as Nu = 0.72 R e 0.45 Pr 0.45 , (35) where Pe is the Peclet number, Re is the Reynolds number, and Pr is the Prandtl number. The heat transfer coefficient of LBE can be obtained by the Nusselt number. To avoid corrosion of the cladding by the liquid metal coolant, the velocity of LBE should not exceed 2 m/s.

When making the cross section library, the selected temperature points should cover all working states, thus, the microscopic cross sections at 600, 700, 800, 1000, and 1200 K are made for fuel, cladding, and coolant materials. For the NE calculation, the order of discrete ordinates is set at S6, the radial mesh is used to create hexagonal assembly, and the axial mesh is 5 cm per layer. For the TH calculation, only fuel assemblies are considered, the axial mesh is computed using the NE calculation, the radial mesh is divided into five segments for fuel, one segment for gap, and two segments for cladding. The power distribution result is shown in Figure 6.

The maximum power location is at the central fuel assembly, and the fuel–rod temperature distribution at the central assembly is shown in Figure 7. The maximum temperature of the fuel, cladding, and coolant is all presented in this assembly, which is shown in Table 1.

Figures 8, 9, 10, 11, 12 show the average fuel temperature, coolant temperature, coolant density, coolant flow rate, and coolant velocity distribution, which are outputted with the hexagonal nodal mesh, respectively. It can be seen that the average temperature distribution of the fuel is similar to the power distribution, and the maximum average temperature of the fuel is less than 1200 K. Since the single-channel model is applied to each parallel channels, the coolant temperature presents a one-dimensional distribution, the density of the coolant changes very little, and its distribution is the opposite of the temperature distribution. The flow rate and velocity distribution are similar to the two-dimensional power distribution, the maximum velocity is equal to 1.8 m/s, which is less than the limit 2 m/s.

The NE/TH coupling calculation converges after two iterations, and the results of neutronics parameters with and without the TH feedback are shown in Table 2. As we can see, whether TH feedback is considered has little influence on the results of the NE calculation. After considering TH feedback, the effective multiplication factor is reduced by 207 pcm, and the axial and radial power peak factors have small changes. Figure 13 compares the axial power distribution between with TH feedback and without TH feedback and the power distribution almost stays the same. It can be considered that, in steady-state calculation of the fast reactor, the TH calculation has little influence on the NE calculation, and TH feedback calculation is generally not required.