The perturbation theory could be adopted for the cases that there is no obvious change of neutron flux distribution for a specific neutronic system before and after a perturbation. Based on the neutron flux distribution before the perturbation, the neutron flux distribution after the perturbation could be obtained by solving the neutron adjoined equation (shown as Eq. 2). − 1 v ∂ ϕ ∗ ∂ t − Ω ⋅ ∇ ϕ ∗ + Σ ϕ ∗ − ∫ 0 ∞ d E ′ ∫ Ω ' Σ s ( r ; E , Ω → E ′ , Ω ′ ) ϕ ∗ ( r , E ′ , Ω ′ , t ) d Ω ′ , = v Σ f ( r , E ) 4 π ∫ 0 ∞ d E ′ ∫ Ω ' χ ( E ' ) ϕ ∗ ( r ′ , E ′ , Ω ′ , t ) d Ω ′ . (2)

Φ fnlowast(r, E, t) is the adjoined neutron flux which is the distribution of neutron value.

The perturbation calculation has been widely used in reactor physics calculation for fission reactors, including the core Doppler coefficient calculation, the differential value calculation of control rods, coolant cavitation value calculation, and so on. The type of perturbation includes density perturbation, cross-section perturbation, temperature perturbation, and so on.

The perturbation theory could also be put to use for the rapid calculation of neutron flux distribution of the solid breeder TBB of fusion reactors. As for a density perturbation, the TBR could be rapidly calculated through Eq. 1 based on the perturbation calculation. As for a geometric perturbation, the TBR could be rapidly calculated through the following flowchart (shown as Figure 3).

The density perturbation calculation for the change of the TBR on the TBB was verified. Comprehensively, the study for the geometry perturbation calculation for the change of the TBR on the TBB based on the virtual density theory could be performed (Hess et al., 1998).

According to the expression for isotropic deformation calculation based on the virtual density theory, the equivalent coefficient can be calculated as follows: ε = N 2 − N 1 N 1 = κ d κ l N 1 − N 1 N 1 = κ d κ l − 1. (3)

ε is the equivalent coefficient of a specific kind of material that indicates the rate of change of atom density (shown as the following formula); N ₁ is the atom density before the deformation, and N ₂ is the equivalent atom density; κ _d is the density change coefficient; κ _l is the linear scale change coefficient; δ _N and δ _ρ are the variation of the atom density and mass density of the specific material, separately. ε = δ N N = δ ρ ρ . (4)

According to the isotropic expansion case of a sphere, the density change coefficient and the linear scale change coefficient are expressed as follows: κ d = 1 f 3 . (5) κ l = f , (6) f is the change coefficient of radium. If f = 1.01, the radium will be increased by 1%. Thus, the equivalent coefficient of an isotropic expansion case can be calculated as follows: ε = κ d κ l − 1 = 1 f 2 − 1. (7)

As for the geometry adjustment (radial expansion or compression of each tritium breeder region or the neutron multiplier region for the HCCB TBB) of the TBB for the neutronic optimization, it is not an isotropic deformation case but an anisotropic deformation one.

In this article, the CFETR HCCB TBB typical module with a sandwich-like breeder zone was also selected for the geometry perturbation method study. According to Figure 1, both the tritium breeder regions and neutron multiplier regions are rectangular solids, which can be described with the radial thickness, toroidal width, and poloidal length. The toroidal width and poloidal length keep as constant during the neutronic optimization and the geometry adjustment for the neutronic optimization can be regarded as the radial expansion or compression case. In this case, there will be a change for the radial thickness (κ _d  = f) and the density change (κ _l  = 1/f) in a specific region, so the radial equivalent coefficient can be calculated as follows: ε r = f f − 1 = 0. (8)

Meanwhile, there is no change in the toroidal width and the poloidal length (d _o  = 1) and a reduction in the density (κ _l  = 1/f); consequently, the toroidal and poloidal equivalent coefficient can be calculated as follows: ε t = ε p = 1 f − 1. (9)

Figure 4 shows the local layout of the HCCB TBB before and after the geometric adjustment. The local layout model consists of the tritium breeder regions, CPs, and neutron multiplier regions. The total radial thickness of the local layout model is assumed to be constant during the geometry adjustment. According to the local layout, l ^CP is the radial thickness of the CP, which is also be assumed to be a constant during the neutronic optimization. l ₁ ^Li4SiO4 and l ₂ ^Li4SiO4 are the radial thickness of the tritium breeder region before and after the geometry adjustment, separately, and Δl is the increment of the tritium breeding region. Correspondingly, l ₁ ^Be and l ₂ ^Be are the radial thickness of the neutron multiplier region before and after the geometry adjustment, separately, and Δl is the decrement of the Be region. ρ ₁ ^Li4SiO4 and ρ ₁ ^Be are the densities of the tritium breeder region and the neutron multiplier region before the geometry adjustment separately, and ρ ₂ ^Li4SiO4 and ρ ₂ ^Be are the densities of the tritium breeder region and the neutron multiplier region after the geometry adjustment. The equivalent coefficients of the tritium breeder region and the neutron multiplier region can be calculated as follows: ε Li 4 SiO 4 = 1 f − 1 = − Δ l l 1 Li 4 SiO 4 + Δ l , (10) ε Be = 1 f − 1 = Δ l l 1 Be − Δ l . (11)

The equivalent density change of these two regions can be calculated as follows: δ ρ Li 4 SiO 4 = − Δ l ⋅ ρ 1 Li 4 SiO 4 l 1 Li 4 SiO 4 + Δ l = − Δ l ⋅ ρ 1 Li 4 SiO 4 l 2 Li 4 SiO 4 , (12) δ ρ Be = Δ l ⋅ ρ 1 Be l 1 Be − Δ l = Δ l ⋅ ρ 1 Be l 2 Be . (13)

ε ^Li4SiO4 and ε ^Be are the equivalent coefficients of the tritium breeder region and the neutron multiplier region, respectively. δρ ^Li4SiO4 and δρ ^Be are the equivalent density change of the tritium breeder region and the neutron multiplier, individually.

Then, the conclusions are made from the local layout model to the HCCB TBB typical model, which can be shown as Figure 5.

In the HCCB TBB typical module with the sandwich-like breeder zone, there are m tritium breeder regions, n neutron multiplier regions, and k CPs (k = m + n − 1). Before the geometry adjustment, the radial thickness and the density of each tritium breeder region are expressed by l _i ^Li4SiO4 and ρ _i ^Li4SiO4 (i = 1,2 … m), respectively. The radial thickness and the density of each neutron multiplier region are defined as l _j ^Be and ρ _j ^Be (j = 1,2 … n), separately. After the geometry adjustment, the change of the radial thickness and the change of the density of each tritium breeder region are Δl _i ^Li4SiO4 and Δρ _i ^Li4SiO4 (i = 1,2 … m), individually. The change of the radial thickness and the change of the density of each neutron multiplier region are Δl _j ^Be and Δρ _j ^Be (j = 1,2 … n). The total radial thickness of the breeder zone of the TBB remains unchanged (shown as the following equation). ∑ i = 1 m Δ l i Li 4 SiO 4 + ∑ j = 1 n Δ l j Be = 0. (14)

Based on the conclusions above, the equivalent coefficients and equivalent change of the density of each tritium breeder region and neutron multiplier region can be calculated as follows: ε i Li 4 SiO 4 = − Δ l i Li 4 SiO 4 l i Li 4 SiO 4 + Δ l i Li 4 SiO 4 ( i = 1,2 ⋯ m ) , (15) ε j Be = − Δ l j Be l j Be + Δ l j Be ( j = 1,2 ⋯ n ) . (16) δ ρ i Li 4 SiO 4 = − Δ l i Li 4 SiO 4 ⋅ ρ i Li 4 SiO 4 l i Li 4 SiO 4 + Δ l i Li 4 SiO 4 ( i = 1,2 ⋯ m ) , (17) δ ρ j Be = − Δ l j Be ⋅ ρ j Be l j Be + Δ l j Be ( j = 1,2 ⋯ n ) . (18)

The neutronic perturbation calculation for the change of the TBR and nuclear heat of each part on the TBB of three cases was performed, and the details of the three cases are listed in Table 2 (Mckinney and Iverson, 1996; Schnabel et al., 2021). Ten percent is selected for the interval of density perturbation (5% is selected for the TBR comparison under case 1).

Two kinds of the results were calculated for comparison: the transport results, and the 1st + 2nd order of the neutronic perturbation results. The TBR comparison of the HCCB TBB of cases 1, 2, and 3 are shown as Figure 6. The nuclear heat comparison (F6 tally results are compared) is shown as Figure 7. The curve of relative bias is also shown in each figure (in red). Regions of each relative bias band are marked with different colors. The variation of TBR and F6 calculated by transport code and perturbation method under each case is shown as Figure 8 and Figure 9 individually. The MC transport results were regarded as the reference. In this article, 1E7 particles are simulated, and some variance reduction techniques (such as weight windows, forced collision, energy splitting, and roulette) are used in the MCNP-4C calculation. In this way, the Monte-Carlo relative deviation could reduce to ∼0.01%. Therefore, the change of TBB neutronics performances were large enough (an order of magnitude larger) compared with the standard deviation of MCNP-4C code, and the perturbation problem could be verified using the MC method.

Analysis toward the above results indicates that the 1st + 2nd order of the neutronic perturbation calculation results (including TBR and nuclear heat) is consistent with the transport results under a density perturbation of -15% to +15% under each case (with a relative bias <0.2%). Meanwhile, the neutronic perturbation calculation is much faster than the transport calculation (the efficiency could be improved by more than 100 times in conservative estimating). Therefore, the perturbation calculation can be a substitute for the transport calculation, which will be a better choice for the rapid neutronic optimization for the TBB.

The CFETR HCCB TBB typical module was chosen for the verification of the geometry perturbation calculation. Geometry adjustments toward the CFETR HCCB TBB of five cases were performed, and the radial adjustment of each region of the HCCB TBB typical module of each case are listed in Table 3. Case 0 is the initial scheme of the CFETR HCCB TBB, which the radial lengths of each region can be found in Table 1. The maximum geometry adjustment of all regions is less than 15% (case 5 in Li-1 region). According to the virtual density theory, the equivalent mass densities of each tritium breeder region and neutron multiplier region can be summarized in Table 4.

Three methods (methods A, B, and C) were used for the neutronic calculation for the HCCB TBB. The model with the geometry adjustment (shown in Table 3) of each case was adopted, and 3D neutronic transport calculation was performed in method A; 3D neutronic transport calculation based on the model with the equivalent density adjustment (shown in Table 4) of each case was made in method B; 3D neutronic transport calculation of the initial scheme was performed, and perturbation calculation (the 1st + 2nd order) based on the virtual density theory of each case was made in model C. The neutron flux of each region calculated by methods A and B is listed and compared in Table 5 (the tritium breeding region) and Table 6 (the neutron multiplier region), separately. According to the design parameters of CFETR phase II, a fusion power of 1.5 GW was assumed, based on which a neutron wall load of 1.69 MW/m² was adopted in the calculations for a single TBB module (Cao et al., 2021). The TBRs calculated by using the three methods are listed in Table 7. The MCNP-4C code was adopted for the 3D neutronic and perturbation calculation based on FENDL-3.2. All the results calculated by method A are regarded as the references. The analysis and conclusions can be portrayed as follows:

(1) Relative deviations that the neutron flux of each region and the TBR are listed in Tables 5, 6, 7, separately for using methods A and B generally showing an increasing trend as the radial geometry variation of each region goes up.

Therefore, the geometry perturbation calculation can be adopted for the solid breeder TBB of the fusion reactor for the rapid neutronic optimization based on the virtual density theory.