Analysis of Fission Products’ Release in Pebble-Bed High-Temperature Gas-Cooled Reactor Fuel Elements Using a Modified FRESCO II Numerical Model

Abstract

The radiation safety design and emergency analysis of an advanced nuclear system highly depends on the source term analysis results. In modular high-temperature gas-cooled reactors (HTGRs), the release rates of fission products (FPs) from fuel elements are the key issue of source term analysis. The FRESCO-II code has been established as a useful tool to simulate the accumulation and transport behaviors of FPs for many years. However, it has been found that the mathematical method of this code is not comprehensive, resulting in large errors for short-lived nuclides and large time step during calculations. In this study, we used the original model of TRISO particles and spherical fuel elements and provided a new method to amend the FRESCO-II code. The results show that, for long-lived radionuclides (Cs-137), the two methods are perfectly consistent with each other, while in the case of short-lived radionuclides (Cs-138), the difference can be more than 1%. Furthermore, the matrix method is used to solve the final release rates of FPs from fuel elements. The improved analysis code can also be applied to the source term analysis of other HTGRs.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Modular high-temperature gas-cooled reactors (HTGRs) are generally considered to have the technical characteristics of a Generation IV (Gen-IV) nuclear energy system [1]. When analyzing the emergency action and safety principles of an HTGR, it is crucial to know the source term under normal and accident conditions. Source term analysis could provide the generation, quantity, release, and radiation hazard of radionuclides in a nuclear power plant. Owing to the use of high performance fuel elements and TRISO particles, the fission products (FPs) can be effectively retained in the reactor core [2]. However, after decades of operation of the reactor, many FPs may escape from the fuel elements (with a lifetime of approximately 3 years) through permeation and diffusion effects, which has contributed to the main source of radioactivity in the primary circuit and other parts of the nuclear power plant [3, 4].

In order to calculate the amount of FPs in the primary circuit, the numerical simulation code FRESCO-II was developed at Forschungszentrum Jülich in Germany in the early 1980s by Krohn and Finken [5]. This code has been used for decades in HTGRs at Germany, China, and other countries; the results have been compared in many benchmark calculations and have been validated through experiments [6]. Based on the FRESCO-II code and other codes (FRESCO-I, PANAMA, SPATRA, and so on), the uniform and integrated code package STACY (Source Term Analysis Code System) was developed by Allelein et al. in 2010 [7]. Recently, the new simulation package FRAT for prediction and evaluation of full core fission product release was developed by Li et al. at Tsinghua University (China) [8].

All of these simulation codes are based on Fick’s law of diffusion [9, 10], and the results are obtained based on different discrete methods and numerical simulation processes. After using the original model and method of the FRESCO-II code, it has been found that there is a small error in the expression of the process coefficient. Consequently, the model is nonphysical when considering nuclide decay terms, resulting in large errors for short-lived nuclides and large time steps. In this study, all coefficients were strictly derived on the basis of the modified model. Additionally, the simulation results of the new and old models are compared herein.

2. Features of HTGR Fuel

The fuel element used in a HTGR is coated spherical fuel particles with a graphite substrate [11]. The research and design of coated fuel particles have gone through two stages. In the initial design, the coated fuel particle structure is of Bi-isotropic (BISO) type, i.e., a loose and a dense pyrolytic carbon layer on the outer coating of the active core. However, after many irradiation experiments and safety tests, the BISO fuel particle gradually exposed its structural shortcomings; these were mainly reflected in the high level of heavy metal contamination during the coating deposition and the low retention ability for certain metal FPs (such as strontium). Based on the BISO fuel particle, the tri-isotropic (TRISO) fuel particle is designed. The main advantage of TRISO fuel particles is their excellent retention ability for various radioactive FPs. Experiments have shown that the integrity of the TRISO structure can be maintained at temperatures below 1600°C.

Considering a pebble-bed modular HTGR as an example, a spherical fuel element with a diameter of 60 mm was used (shown in Figure 1). Each fuel element contains approximately 12,000 TRISO fuel particles consisting of five parts: (1) an enriched UO2 kernel surrounded by subsequent layers, (2) a loose pyrolytic carbon (buffer) layer, (3) an inner pyrolytic carbon (IPyC) layer, (4) silicon carbide (SiC), and finally, (5) an outer pyrocarbon (OPyC) layer [11].

[figure omitted; refer to PDF]

If the fuel element, TRISO particles, and graphite grains are regarded as uniform spheres, the diffusion of FPs in these media can be described by the same equation. According to Fick’s diffusion law, the release fractions of FPs from the spherical fuel elements can be calculated using the following diffusion equation [12]:$$\begin{array}{}\text{(1)}& \frac{\partial c}{\partial t}=DT\frac{{\partial }^{2}c}{\partial t^2}+\frac{2}{r}\cdot \frac{\partial c}{\partial t}+Q-\lambda C,\end{array}$$where C is the FP concentration in the diffusion medium (in atoms/cm³) and can be obtained from the inventory, t is the time (in s), r is the spatial coordinate (in cm), Q is the production rate of FPs (in atom/cm³ • s), λ is the decay constant of nuclide (in s⁻¹), and D (T) is the diffusion coefficient (in cm²/s) depending on temperature T (in K) and activation energy E (in J/mol) given by the following equation:$$\begin{array}{}\text{(2)}& DT=D_0\cdot \exp -\frac{E}{R\cdot T},\end{array}$$where R is the ideal gas constant (R = 8.3143 J/(mol·K)) and D₀ is the diffusion frequency factor (in cm²/s). The corresponding boundary condition can be expressed as$$\begin{array}{}\text{(3)}& {-D\frac{\partial c}{\partial r}|}_{r=r_p}=\beta C_{\text{rp}}-C_{\text{gr}},\end{array}$$where β is the interface diffusion between different materials, $C_{\text{rp}}$ is the concentration of FPs on the surface of the fuel element, and $C_{\text{gr}}$ is the average concentration of FPs in the helium outside the fuel element.

4. Numerical Simulation

Theoretically, the concentration at any time and position can be obtained by solving equations (1)–(3). However, the explicit expression of C(r,t) is usually very complex and difficult to solve [13].

To address this challenge, the numerical simulation method is used in the FRESCO-II code and other related codes [14]. In this method, the whole sphere is divided into N thin layers of spherical shell with a volume of V_i (i = 1/2, ... N) (presented in Figure 3(a)).

[figures omitted; refer to PDF]

Considering that the concentration at the i^th spherical shell is C_i and assuming that the concentration is a linear function of radius r, the concentration between two adjacent spherical shells can be expressed as follows (shown in Figure 3(b)):$$\begin{array}{}\text{(4)}& c=c_i+\frac{c_{i+1}-c_i}{r_{i+1}-r_i}r-r_i.\end{array}$$

In particular, assuming that the concentration distribution is a quadratic curve for the central sphere because the concentration gradient is zero at the central point, the concentration can be derived as follows:$$\begin{array}{}\text{(5)}& c=c_N-2\frac{c_{N+1}-c_N}{r_N}r-r_N-\frac{c_{N+1}-c_N}{r_N^2}{r-r_N}^2,\end{array}$$where c_N+1 represents the concentration at the central point. For the i^th spherical shell, the mass conservation equation is as follows:$$\begin{array}{}\text{(6)}& \frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{V_i}c\cdot \text{d}V}=Q\cdot V_i-\lambda {\displaystyle {\int }_{V_i}c\cdot \text{d}V}-J_i\cdot F_i+J_{i+1}\cdot F_{i+1},\end{array}$$where the left side of this equation represents the concentration change during dt time. The right side represents the various sources of concentration change, i.e., the generation of fission with rate Q, the disappearance of decay, and the diffusion from the inner layer into the target layer and from the target layer into the outer layer, respectively. According to Fick’s first diffusion law, the diffusion rate per unit area J can be expressed as$$\begin{array}{c}\text{(7)}& J_i=-D_i\frac{c_{i+1}-c_i}{r_{i+1}-r_i},\hspace{1em}i\ne N,J_N=-2D_N\frac{c_{N+1}-c_N}{r_{N+1}-r_N},\hspace{1em}i=N.\end{array}$$

Then, we can discretize the mass conservation equation (6) using the average change rates instead of instantaneous change rates. For the left side of equation (6), it can be approximated by replacing dt with Δt as follows:$$\begin{array}{}\text{(8)}& \frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{V_i}c\cdot \text{d}V}=\frac{1}{\Delta t}{\displaystyle {\int }_{V_i}{c}^{k+1}\cdot \text{d}V}-{\displaystyle {\int }_{V_i}{c}^k\cdot \text{d}V},\end{array}$$where the superscript index k represents the time step; a larger value corresponds to a later time point.

It should be emphasized that, until this step, all the symbols and equations are the same as the ones in FRESCO-II code. However, the next derivation proposes amendments to the FRESCO-II code for the right side of equation (6). In the FRESCO-II code, the concentration appearing on the right side is set to be a future time value, i.e., the top index is k + 1. Obviously, this assumption is nonphysical because the concentrations in each spherical shell are all changing dynamically. Therefore, a new method is proposed through using the time average concentrations instead of future time concentrations; this is also consistent with the treatment of the left side of equation (6). Then, the expression is derived as follows:$$\begin{array}{c}\text{(9)}& Q\cdot V_i-\lambda {\displaystyle {\int }_{V_i}c\cdot \text{d}V}-J_i\cdot F_i+J_{i+1}\cdot F_{i+1}=Q\cdot V_i-\frac{\lambda }{2}{\displaystyle {\int }_{V_i}{c}^{k+1}\cdot \text{d}V}+{\displaystyle {\int }_{V_i}{c}^k\cdot \text{d}V}-\frac{1}{2}{J}_i^{k+1}+J_i^k\cdot F_i+\frac{1}{2}{J}_{i+1}^{k+1}+J_{i+1}^k\cdot F_{i+1}.\end{array}$$

From the combination of equations (6)–(9), the new equation can be reorganized as$$\begin{array}{}\text{(10)}& a_{i+1}{c}_i^{k+1}+b_{i+1}{c}_{i+1}^{k+1}+d_{i+1}{c}_{i+2}^{k+1}=e_{i+1},\hspace{1em}i=\mathrm{1,2,3},\mathrm{...},N-1.\end{array}$$

The coefficients are given as follows:$$\begin{array}{}\text{(11)}& a_{i+1}={\gamma }_{i}1+\frac{1}{2}\lambda \Delta t+\frac{1}{2}{D}_i\frac{F_i}{V_i}\frac{\Delta t}{r_{i+1}-r_i},\text{(12)}& b_{i+1}=1-{\gamma }_i\times 1+\frac{1}{2}\lambda \Delta t-\frac{1}{2}{D}_i\frac{F_i}{V_i}\frac{\Delta t}{r_{i+1}-r_i}-\frac{1}{2}{D}_{i+1}\frac{F_{i+1}}{V_i}\frac{\Delta t}{r_{i+2}-r_{i+1}},\text{(13)}& d_{i+1}=\frac{1}{2}{D}_{i+1}\frac{F_{i+1}}{V_i}\frac{\Delta t}{r_{i+2}-r_{i+1}},\text{(14)}& e_{i+1}={\gamma }_{i}1-\frac{1}{2}\lambda \Delta t-\frac{1}{2}{D}_i\frac{F_i}{V_i}\frac{\Delta t}{r_{i+1}-r_i}{c}_i^k+1-{\gamma }_{i}1-\frac{1}{2}\lambda \Delta t+\frac{1}{2}{D}_i\frac{F_i}{V_i}\frac{\Delta t}{r_{i+1}-r_i}+\frac{1}{2}{D}_{i+1}\frac{F_{i+1}}{V_i}\frac{\Delta t}{r_{i+2}-r_{i+1}}1-{\gamma }_i{c}_{i+1}^k-\frac{1}{2}{D}_{i+1}\frac{F_{i+1}}{V_i}\frac{\Delta t}{r_{i+2}-r_{i+1}}{c}_{i+2}^k+Q_i\Delta t,\text{(15)}& {\gamma }_i=1-\frac{3}{r_{i+1}^3-r_i^3}\cdot \frac{1}{r_{i+1}-r_i}\frac{r_{i+1}^4-r_i^4}{4}-\frac{r_{i+1}^3{r}_i-r_i^4}{3}.\end{array}$$

Notably, equation (15) must be identical to that in the FRESCO-II code; nevertheless, the incorrect subscript index i-1 of the radius cubic in the rightmost fraction appeared in the FRESCO-II code.

The two boundary conditions for the outermost shell contact with the environment and the central shell can be expressed as follows:$$\begin{array}{c}\text{(16)}& -D_1\frac{c_2^{k+1}-c_1^{k+1}}{r_2-r_1}=\beta c_1^{k+1}-c_{\text{env}}⟹\beta -\frac{D_1}{r_2-r_1}{c}_1^{k+1}+D_1\cdot \frac{1}{r_2-r_1}\cdot c_2^{k+1}=\beta c_{\text{env}}\\ ⟹b_1{c}_1^{k+1}+d_1{c}_2^{k+1}=e_1,\\ a_{N+1}={\gamma }_{N}1+\frac{1}{2}\lambda \Delta t+\frac{1}{2}{D}_N\frac{F_N}{v_N}\cdot \frac{2\Delta t}{-r_N},\\ b_{N+1}=1-{\gamma }_{N}1+\frac{1}{2}\lambda \Delta t-\frac{1}{2}{D}_N\frac{F_N}{v_N},\\ d_{N+1}=0,\\ e_{N+1}={\gamma }_{N}1-\frac{1}{2}\lambda \Delta t-\frac{1}{2}{D}_N\frac{F_N}{v_N}\frac{2\Delta t}{-r_N}{c}_N^k+1-{\gamma }_{i}1-\frac{1}{2}\lambda \Delta t+\frac{1}{2}{D}_N\frac{F_N}{v_N}\frac{2\Delta t}{-r_N}{c}_{N+1}^k+Q_N\Delta t,\\ {\gamma }_N=\mathrm{0.6.}\end{array}$$

Equation (10) is a recurrence relation for both time and space dimensions and can be further represented in the matrix form as follows:$$\begin{array}{}\text{(17)}& \mathbb{M}\cdot {ℂ}^{k+1}={\mathbb{E}}^k,\text{(18)}& \mathbb{M}=\begin{array}{ccc}\begin{array}{ccc}{b}_1& d_1& 0\\ a_2& b_2& d_2\\ 0& a_3& b_3\end{array}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & \begin{array}{ccc}{b}_{N-1}& d_{N-1}& 0\\ a_N& b_N& d_N\\ 0& a_{N+1}& b_{N+1}\end{array}\end{array},\text{(19)}& {ℂ}^{k+1}={c_1,c_2,\dots ,c_{N+1}}^{\text{T}},\text{(20)}& {\mathbb{E}}^k={e_1,e_2,\dots ,e_{N+1}}^{\text{T}}.\end{array}$$

Equation (17) has a unique solution because the matrix $\mathbb{M}$ is full rank. Furthermore, the elements of the matrix $\mathbb{M}$ are all time independent, i.e., they are only related to the space index and can be determined by geometric parameters. The vector ${\mathbb{E}}^{\text{k}}$ is a function of the space index and the concentrations in past time points. The unknown vector ${\mathrm{ℂ}}^{\text{k}+1}$ is our target concentration in the future time point, which can be obtained by solving the matrix equation (17) as follows:$$\begin{array}{}\text{(21)}& {ℂ}^{k+1}={\mathbb{M}}^{-1}\cdot {\mathbb{E}}^k.\end{array}$$

This equation indicates that any concentrations in any future time point can be recursively solved when an initial concentration state is given.

5. Results and Discussion

The numerical simulation method introduced in Section 4 has been implemented in Python and is named Fuel Element Release Rate Analysis (FERRA) Code. Considering the HTR-PM (high-temperature gas-cooled reactor-pebble bed modules) fuel element as an example, a series of comparative calculations are performed in this study. The main parameters of the HTR-PM are listed in Table 1 [2].

Table 1

Parameters of fuel particles in HTR-PM.

Parameter	Normal data
Diameter of kernel	500 μm
Thickness of buffer layer	95 μm
Thickness of OPyC layer	40 μm
Thickness of IPyC layer	40 μm
Thickness of SiC layer	35 μm
Diameter of graphite grain	15 μm
Failure fraction of TRISO at objective burn-up	2.60E − 04
Natural uranium contamination of fuel elements	7.00E − 07

The kernel and the four-layer coating materials are separated into 39 thin spherical shells in both the FRESCO-II and FERRA codes. Two cesium isotopes, Cs-137 and Cs-138 (with a decay constant of approximately 7.3E − 10 s⁻¹ and 3.6E − 04 s⁻¹, respectively), were used as an example to compare the differences between the results of the two methods. The diffusion coefficients of Cs in each material are listed in Table 2 [11].

Table 2

Diffusion coefficients.

	Kernel	Buffer layer	IPyC	SiC	OPyC
D0 (cm2/s)	5.6E − 04	1.0E − 08	6.3E − 04	1.0E − 14	6.3E − 04
E (J/mol)	2.1E + 05	0	2.2E + 05	1.0E + 05	2.2E + 05

Based on these parameters, after a 25,400 h burn-up (3 years average life time of fuel element), the fractional release rates of Cs-137 and Cs-138 from a pebble over the temperature range from 600°C to 1600°C are calculated and presented in Tables 3 and 4.

Table 3

Comparison of two methods for Cs-137 fractional release rates from a pebble.

Temperature (°C)	600	700	800	900	1000	1100	1200	1300	1400	1500	1600
FRESCO-II (s−1)	3.9192E − 16	1.7100E − 15	6.5096E − 14	6.3373E − 13	2.4107E − 12	2.3764E − 10	1.5134E − 09	3.9838E − 09	7.1964E − 09	9.7246E − 09	1.0717E − 08
FERRA (s−1)	3.9198E − 16	1.7096E − 15	6.5130E − 14	6.3414E − 13	2.3990E − 12	2.3774E − 10	1.5145E − 09	3.9864E − 09	7.2005E − 09	9.7284E − 09	1.0718E − 08
Ratio	0.9998	1.0002	0.9995	0.9994	1.0049	0.9996	0.9993	0.9993	0.9994	0.9996	0.9999

Table 4

Comparison of two methods for Cs-138 fractional release rates from a pebble.

Temperature (°C)	600	700	800	900	1000	1100	1200	1300	1400	1500	1600
FRESCO-II (s−1)	1.6609E − 15	2.1412E − 14	1.6480E − 13	7.8264E − 13	2.3831E − 12	5.3521E − 12	1.0020E − 11	1.6930E − 11	3.1376E − 11	8.9169E − 11	3.0714E − 10
FERRA (s−1)	1.6609E − 15	2.1412E − 14	1.6480E − 13	7.8264E − 13	2.3831E − 12	5.3521E − 12	1.0027E − 11	1.7014E − 11	3.1707E − 11	8.9918E − 11	3.0952E − 10
Ratio	1.0000	1.0000	1.0000	1.0000	1.0000	1.0000	0.9993	0.9951	0.9896	0.9917	0.9923

The results show that, for the long-lived radionuclide Cs-137, the two methods are perfectly consistent with each other at all temperatures with an error less than 1%. However, the case for the short-lived radionuclide Cs-138 is slightly different. A maximum of more than 1% error occurred when the temperature was higher than 1200°C. Furthermore, the release rates are overall underestimated in the FRESCO-II code; this can be attributed to the underestimation of the interlayer transport rates.

6. Conclusions

Thus far, the FRESCO-II code and its successor provided powerful simulation tools for handling the FP release rates in HTGR. Based on reliable release rates, the source term analysis could be credible and safety design and emergency evaluation could be performed.

In order to improve the reliability of the source term analysis codes, we used the original model and derivation process of the FRESCO-II code. An index error has been identified, and a new treatment for transport rates is proposed. After the code implementation, the new simulation code FERRA is developed, and a series of tests are performed. The results show that the FRESCO-II and FERRA codes are perfectly consistent with each other for both long-lived and short-lived radionuclides. Another valuable finding is that the release rates are overall underestimated in the FRESCO-II code, which suggests that the FERRA results are more conservative. As a conclusion, continuous improvements to computing tools can make a great contribution to the design of pebble-bed HTGRs.

Disclosure

The contents of this manuscript are not now under consideration for publication elsewhere.

Authors’ Contributions

All authors of this manuscript have directly participated in planning, execution, and/or analysis of this study.

Acknowledgments

This work was supported by the National S & T Major Project (Grant no. 2018ZX06902013).

References

[1] C. Fang, R. Morris, F. Li, "Safety features of high temperature gas cooled reactor," Science and Technology of Nuclear Installations, vol. 2017,DOI: 10.1155/2017/9160971, 2017.

[2] Z. Zhang, Y. Dong, F. Li, "The Shandong Shidao Bay 200 MW e high-temperature gas-cooled reactor pebble-bed module (HTR-PM) demonstration power plant: an engineering and technological innovation," Engineering, vol. 2 no. 1, pp. 112-118, DOI: 10.1016/j.eng.2016.01.020, 2016.

[3] A. F. Michael, E. D’Agata, X. Raepsaet, "Is tritium an issue for high temperature reactors?," Nuclear Engineering and Design, vol. 306, pp. 160-169, 2016.

[4] Y. Liu, C. Jianzhu, "Fission product release and its environment impact for normal reactor operations and forrelevant accidents," Nuclear Engineering and Design, vol. 218 no. 1, pp. 81-90, DOI: 10.1016/S0029-5493(02)00200-5, 2002.

[5] H. Krohn, R. Finken, "FRESCO-II: a computer program for analysis of fission product release from spherical HTR fuel element in irradiation and annealing experiments," 1983. Internal Report JUEL-SPEZ--212

[6] K. Verfondern, J. Cao, T. Liu, H.-J. Allelein, "Conclusions from V and V studies on the German codes PANAMA and FRESCO for HTGR fuel performance and fission product release," Nuclear Engineering and Design, vol. 271, pp. 84-91, DOI: 10.1016/j.nucengdes.2013.11.015, 2014.

[7] H.-J. Allelein, S. Kasselmann, A. Xhonneux, S.-C. Herber, "Progress on the development of a fully integrated HTR code package," Nuclear Engineering and Design, vol. 251, pp. 400-406, DOI: 10.1016/j.nucengdes.2011.09.053, 2012.

[8] J. Li, D. She, L. Shi, "Status of development of an integrated source term analysis code package for HTGR," Proceedings of the HTR 2018, .

[9] J. Crank, The Mathematics of Diffusion, 1956.

[10] W. Jost, “Diffusion in Solids, Liquids, Gases” (3rd Printing with Addendum), 1960.

[11] IAEA, "Fuel performance and fission product behavior in gas cooled reactors," 1997. Report IAEA-TECDOC-978

[12] C. Yang, C. Fang, J. Zhang, "Study on cumulative fractional release of radionuclides in HTGR fuel particles," Nuclear Science and Techniques, vol. 63 no. 3,DOI: 10.7498/aps.63.032802, 2014.

[13] C. Tian, L. Sun, C. Fang, "Exact solutions to the diffusion process of fission products in the TRISO fuel particles of pebble bed HTRs: case study of137Cs and134Cs," Nuclear Science and Engineering, vol. 175 no. 2, pp. 204-211, DOI: 10.13182/nse12-51, 2013.

[14] H. Krohn, R. Finken, "FRESCO-II: Ein Rechenprogramm zur Berechnung der Spaltprodukt freisetzung aus Kugelförmigen HTR-Brennelementen in Bestrahlungs- und Ausheizexperimenten," , 1983.