Numerical modelling of modular high-temperature gas-cooled reactors with thorium fuel

The design of the nuclear reactor core requires the application of complex methods for reliable reconstruction of the geometry and material composition of the actual reactor core. The greater the complexity of the core geometry is, the longer the computation time of neutron transport will be, especially when using Monte Carlo methods. Additionally, the execution of many input files with varied key parameters is necessary for a comprehensive parametric study. The problem is particularly significant in the modelling of modular high-temperature gas-cooled reactors (MHTRs) with graphite fuel compacts containing spherical TRi-structural ISOtropic (TRISO) particles, called double heterogeneity [1]. In the paper, an alternative method called volumetric homogenization is presented for reducing complexity of the numerical model and the resulting computation time. The method assumes the elimination of double heterogeneity of TRISO particles embedded in the fuel compact. In this method, materials of TRISO particles are admixed to the graphite matrix of the compact and by this way, the complicated geometrical structure of multilayer TRISO particles is eliminated and replaced by one homogenized material. This, in turn, simplifies the numerical model and speeds up the simulations. The method was applied for the modelling of the MHTR core with thorium-uranium fuel [2, 3]. This kind of fuel mix is envisaged as an option for the currently designed nuclear reactors, including MHTR [4, 5]. Example results comprise the time evolutions of the effective neutron multiplication factor K_eff and of fissionable isotopes (²³³U, ²³⁵U, ²³⁹Pu, ²⁴¹Pu) for a few configurations of the initial reactor core. The parametric study includes a variation of geometrical and material parameters of TRISO particles, i.e. kernel radius, packing fraction, enrichment etc. The numerical model of the reactor was developed using the Monte Carlo continuous energy burnup code (MCB) [6, 7]. In the next section, the volumetric homogenization method is comprehensively described. Section “Results” focuses on the design of numerical model and example results. The study is summarized in the last section.

In this section, the volumetric homogenization method is described from the mathematical point of view. Table 1 shows the input parameters necessary to obtain the weight fractions of each isotope in the homogenized material mix and the density of the homogenized material for the Monte Carlo modelling of neutron transport. The presented mathematical set-up was transferred to the numerical script, which automatically calculates the aforementioned parameters using the input data from Table 1.

The calculations begin with the assessment of the volumes of each TRISO material in the fuel compact. The Appendix 1 contains list of symbols used in the applied equations. Initially, the volume of the spherical TRISO fuel kernel V^KR is calculated using Eq. (1), where r_KR is the kernel radius.

Next, the volumes V_i^TRISO of each TRISO material wrapping the kernel subsequently are calculated using Eq. (2). r_i, is the kernel radius r_KR increased by the material layer thickness t_i, i.e. buffer porous carbon (BPC), inner pyrolytic carbon (IPyC), silicon carbide (SiC) and outer pyrolytic carbon (OPyC), subsequently. For the calculating the volume of the BPC that is directly laid on the fuel kernel, it is known that r_i−1 equals r_KR.

The volume V^COM of the cylindrical fuel compact is calculated using Eq. (3), where r_COM is the compact radius and h_COM is the compact height.

In addition, Eq. (1) is used for the calculation of the TRISO particle volume V^TRISO, where r_KR is replaced by the radius of the whole TRISO particle r_TRISO.

The volumes calculated above may be used in Eq. (4) for calculations of the number of TRISO particles in the fuel compact N^TRISO for a given packing fraction P. The packing fraction is defined as the overall volume of all TRISO particles in the compact in relation to the compact volume (Eq. (5)).

Then, the volume of all TRISO particles in the compact V^TRISOinCOM and associated volume fractions of each material of TRISO in the compact V_i^TRISO,FR may be calculated by using Eqs. (6) and (7): (6) VTRISOinCOM=VTRISO⋅NTRISO {V^{{\rm{TRISOinCOM}}}} = {V^{{\rm{TRISO}}}} \cdot {N^{{\rm{TRISO}}}} (7) ViTRISO,FR=ViTRISOVCOM V_i^{{\rm{TRISO}},{\rm{FR}}} = {{V_i^{{\rm{TRISO}}}} \over {{V^{{\rm{COM}}}}}}

Then, with the given densities of each material in the compact ρ_i^COM, it is possible to calculate the mass of each material in the compact m_i^COM and its weight fraction w_i^COM,FR by means of Eqs. (8) and (9). m^COM is the total mass of the compact for a given number of TRISO particles; see Eq. (10). (8) miCOM=ViTRISO,FR⋅ρiCOM m_i^{{\rm{COM}}} = V_i^{{\rm{TRISO}},{\rm{FR}}} \cdot \rho _i^{{\rm{COM}}} (9) wiCOM,FR=miCOMmCOM w_i^{{\rm{COM}},{\rm{FR}}} = {{m_i^{{\rm{COM}}}} \over {{m^{{\rm{COM}}}}}} (10) mCOM=∑i=1nmiCOM {m^{{\rm{COM}}}} = \sum\limits_{i = 1}^n {m_i^{{\rm{COM}}}}

For the input of the Monte Carlo simulations, the density of the homogenized compact ρ^HOMO is also needed. It can be calculated using given densities ρ_i^COM and the calculated density of the (U,Th)O₂ (Eq. (17)), using Eq. (11). (11) aiU=wiU⋅M¯UMiU a_i^{\rm{U}} = w_i^{\rm{U}} \cdot {{{{\overline {\rm{M}}}^{\rm{U}}}} \over {{\rm{M}}_i^{\rm{U}}}}

In the next step, the material isotopic composition in the fuel compact is calculated. Initially, the isotopic composition of (U,Th)O₂ dioxide for given weight fractions of each dioxide w_i^(HM)O₂,FR in the dioxide mix is considered. The main input parameter is the isotopic composition of uranium with a defined enrichment in ²³⁵U. The isotopic composition of uranium is usually given as a weight fraction of each uranium isotope wiU w_i^{\rm{U}} . The isotopic compositions of light TRISO materials and graphite matrix are defined by natural abundances, which correspond to atom fractions. Therefore, the isotopic composition of uranium should be recalculated to atom fractions, which was performed using Eq. (12). In the case of thorium, the following recalculation is not necessary since it mainly contains one isotope: ²³²Th. aiU a_i^{\rm{U}} is the atom fraction of each uranium isotope, MiU M_i^{\rm{U}} is the atomic weight of each uranium isotope and M̄^U is the average atomic weight for given enrichment, which is calculated using Eq. (13). (12) aiU=wiU⋅M¯UMiU a_i^{\rm{U}} = w_i^{\rm{U}} \cdot {{{{\overline {\rm{M}}}^{\rm{U}}}} \over {{\rm{M}}_i^{\rm{U}}}} (13) M¯U=(∑i=1nwiUMiU) {\overline {\rm{M}} ^{\rm{U}}} = \left({\sum\limits_{i = 1}^n {{{w_i^U} \over {M_i^U}}}} \right)

Further, the atom fractions a_i^(HM)O₂ of each isotope in UO₂ and ThO₂ dioxides are calculated using Eq. (14), where the index HM corresponds to heavy metal (either U or Th). N_i ^(HM)O₂ is the average number of atoms of each isotope in the dioxide and N^{(HM) O₂} is the number of atoms in the dioxide molecule. (14) ai(HM)O2=Ni(HM)O2N(HM)O2 a_i^{({\rm{HM}}){{\rm{O}}_2}} = {{N_i^{({\rm{HM}}){{\rm{O}}_2}}} \over {{N^{({\rm{HM}}){{\rm{O}}_2}}}}}

Then, the atom fraction a_i^(HM)O₂,FR of each dioxide in the dioxide mix is calculated using Eq. (15). The weight fractions of each dioxide in the mix w_i^(HM)O₂,FR are given the input parameters. The atomic weights of each dioxide M^(HM)O₂ are calculated using Eq. (16), while the average atomic weight of the mixed dioxides M̄^(HM)O₂,MIX is calculated using Eq. (17), which is a modification of Eq. (13) for molecules. (15) ai(HM)O2,FR=wi(HM)O2,FRM¯(HM)O2,MIXMi(HM)O2 a_i^{({\rm{HM}}){{\rm{O}}_2},{\rm{FR}}} = w_i^{({\rm{HM}}){{\rm{O}}_2},{\rm{FR}}}{{{{\overline {\rm{M}}}^{({\rm{HM}}){{\rm{O}}_2},{\rm{MIX}}}}} \over {{\rm{M}}_i^{({\rm{HM}}){{\rm{O}}_2}}}} (16) M¯(HM)O2=MHM+2M0 {\overline {\rm{M}} ^{({\rm{HM}}){{\rm{O}}_2}}} = {{\rm{M}}^{{\rm{HM}}}} + 2{{\rm{M}}^0} (17) M¯(HM)O2,MIX=(∑i=1nwi(HM)O2,FRMi(HM)O2) {\overline {\rm{M}} ^{({\rm{HM}}){{\rm{O}}_2},{\rm{MIX}}}} = \left({\sum\limits_{i = 1}^n {{{w_i^{({\rm{HM}}){{\rm{O}}_2},{\rm{FR}}}} \over {{\rm{M}}_i^{({\rm{HM}}){{\rm{O}}_2}}}}}} \right)

Afterwards, the density of mixed dioxides ρ^(HM)O₂ for given densities of each dioxide ρ_i^(HM)O₂ is calculated using Eq. (18) and it is further used for the calculation of the overall density of the homogenized material mix ρ^(HM) (Eq. (11)). (18) ρ(HM)O2=(∑i=1nwi(HM)O2,FRρi(HM)O2) {\rho ^{({\rm{HM}}){{\rm{O}}_2}}} = \left({\sum\limits_{i = 1}^n {{{w_i^{({\rm{HM}}){{\rm{O}}_2},{\rm{FR}}}} \over {\rho _i^{({\rm{HM}}){{\rm{O}}_2}}}}}} \right)

Additionally, natural abundances of no-fuel compact and TRISO materials a_i^ISO should be recalculated in terms of weight fractions w_i^ISO, which was performed using Eq. (19). (19) wiISO=aiISOMiM w_i^{{\rm{ISO}}} = a_i^{{\rm{ISO}}}{{{{\rm{M}}_i}} \over {\rm{M}}}

The last stage of calculations considers the preparation of the weight fractions of the homogenized material mix w_i^HOMO,ISO and w_i^{HOMO,(HM)O₂} for neutron transport Monte Carlo simulations. Further, w_i^{HOMO,(HM)O₂} is calculated using Eq. (20), and w_i^HOMO,ISO using Eq. (22). (20) wiHOMO,(HM)O2=wiCOM,FR⋅wi(HM)O2,FR⋅wi(HM)O2 w_i^{{\rm{HOMO}},({\rm{HM}}){{\rm{O}}_2}} = w_i^{{\rm{COM}},{\rm{FR}}} \cdot w_i^{({\rm{HM}}){{\rm{O}}_2},{\rm{FR}}} \cdot w_i^{({\rm{HM}}){{\rm{O}}_2}} w_i^(HM)O₂ corresponds to the weight fractions of each isotope in the dioxide and it is calculated using Eq. (21). The w_i^COM,FR fractions were calculated using Eq. (9). (21) wi(HM)O2=ai(HM)O2⋅M(HM)O2M¯i(HM)O2 w_i^{({\rm{HM}}){{\rm{O}}_2}} = a_i^{({\rm{HM}}){{\rm{O}}_2}} \cdot {{{{\rm{M}}^{({\rm{HM}}){{\rm{O}}_2}}}} \over {\overline {\rm{M}} _i^{({\rm{HM}}){{\rm{O}}_2}}}} (22) wiHOMO,ISO=wiCOM,FR⋅wiISO w_i^{{\rm{HOMO}},{\rm{ISO}}} = w_i^{{\rm{COM}},{\rm{FR}}} \cdot w_i^{{\rm{ISO}}}

The weight fractions of carbon isotopes from different materials are finally summarized to simplify the input for numerical simulations and by this way, the final isotopic composition of the homogenized mix of fuel compact material was defined.

The example calculations were performed for 180 MW_th MHTR core, as depicted in Fig. 1. The core comprises of 310 hexagonal fuel blocks including 60 blocks with control rods (CR) holes. The numerical Monte Carlo simulations were performed for a two-year irradiation cycle. The MCB code was equipped with JEFF3.2 nuclear data libraries for neutron transport simulations. The reactor core was divided into four radial and 10 axial fuel zones, which gives a total of 40 fuel zones. The geometrical parameters of the core and the TRISO fuel are presented in Tables 2 and 3, respectively.

Fig. 1

Table 4 shows four versions of the homogenized thorium-uranium fuel composition. The main criterion is a similar mass of fissionable ²³⁵U in the fresh reactor core – about 205 kg. This may be achieved by a variation of the input parameters for the volumetric homogenization method, i.e. kernel radius, weight fraction of each dioxide in the dioxide mix, enrichment and packing fraction.

Figure 2 shows the evolution of K_eff for the investigated cases. The shapes of the curves are similar in each case. At the beginning, the decrease in K_eff due to the formation of absorbing fission products is observed. Then, K_eff increases to its peak value at about 200 days of irradiation due to the breeding of ²³³U, ²³⁹Pu and ²⁴¹Pu from residual ²³²Th and 238U. Afterwards, K_eff starts to decrease because of the ongoing fuel burnup. The highest values of K_eff were observed in the case C2 with the lowest mass of Th in the reactor core because of the lower fuel kernel radius. Although the mass of thorium for the subsequent cases C4 and C1 is similar, the relative values of K_eff are different; this difference originates from the variations in the initial ²³⁸U mass, which is lower in the case C4. The lowest values for K_eff were obtained in the case C3, characterized by the highest volume of ²³²Th and ²³⁸U.

Fig. 2

The evolution of ²³³U bred from ²³²Th is presented in Fig. 3. The mass of the produced ²³³U depends strictly on the initial mass of ²³²Th. The higher the initial mass of ²³²Th, the faster the production of ²³³U because of the higher absorption macroscopic cross section of ²³²Th. Therefore, the fastest production of ²³³U occurs in the case C3 and the lowest is observed in the case C2. The decrease of ²³⁵U due to the fuel burnup is shown in Fig. 4. The characteristics of all cases C1–C4 are similar due to the presence of almost the same amount of fissionable ²³⁵U in the reactor core. The initial mass of ²³⁵U decreases by about 60% during the irradiation period of two years. The production of ²³⁹Pu depends on the initial mass of ²³⁸U, as shown in Fig. 5. The case C4 with the lowest initial mass of ²³⁸U shows the lowest production of ²³⁹Pu and ²⁴¹Pu (Fig. 6). The production of ²³⁹Pu and ²⁴¹Pu is similar in other cases, where the initial mass of ²³⁸U is similar. However, during the ongoing fuel burnup, a difference in the shapes of the plutonium curves appears, which is related to the depletion of fissionable ²³⁵U and its replacement by ²³³U and plutonium isotopes as secondary fuel. This process is unambiguously observed in the case C2, where the mass of ²³⁹Pu starts to decrease after reaching its peak at about 500 days of irradiation. In this case, the production of ²³³U is the lowest; therefore, ²³⁵U is replaced by fissionable plutonium isotopes, mainly ²³⁹Pu.

Fig. 3

Fig. 4

Fig. 5

Fig. 6

The volumetric homogenization method for the modelling of the MHTR core was presented in this study. This method was applied for the numerical Monte Carlo modelling of neutron transport in the MHTR core with thorium-uranium fuel. The example results of the modelling were presented. The results prove the reliability of the method for the initial screening of the reactor core performance. The universal character of this method makes it suitable for the numerical modelling of any type of material fuel composition and geometry. Thus, the method can be used not only for the MHTR modelling but also for the modelling of any fissionable system with a complicated fuel geometry, especially using linear chain method [8]. Further study on this method will focus on benchmarking with a full 3D MHTR core model with double heterogeneity of the TRISO fuel. This will allow the verification of the simplified models developed and their use as a replacement of detailed models. In the benchmarking process, the condition under which a defined replacement may happen will be also determined. The benchmarking will define the areas of the reactor core physics that can be reliably modelled using the volumetric homogenization method, e.g. radiotoxicity of the spent nuclear fuel [9], which may facilitate the whole MHTR design methodology.