Solid heat conduction

Temperature is important for fuel rod performance analysis because of its influence on material properties, fission gas release, and irradiation behavior. Eq.1 describes the heat conduction in the fuel rod with heat source $S$, which is solved to obtain the temperature distribution of the fuel pellet and cladding. $\rho c_{\text{p}}\frac{\partial T}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(\lambda r\frac{\partial T}{\partial r}\right)+\frac{\partial }{\partial z}\left(\lambda \frac{\partial T}{\partial z}\right)+S$ (1)

Gap heat conduction

The change in gap heat conduction, which depends on the gap width, contact state, gas composition, and plenum pressure, causes a change in the temperature of the fuel rod because of its high thermal resistance. A modified Ross and Stoute model [15] was adopted to calculate the heat conduction ($h_{\text{gap}}$) in the gap under both normal and transient conditions, as shown in Eq.2. It is composed of three components: gas conduction; solid contact conduction, which is considered when there is contact between fuel pellet and cladding; and radiation between the fuel pellet and cladding surface. Gas conduction, $h_{\text{gas}}$, is related to surface roughness, and the model proposed by Ross and Stoute is adopted [15]. When pellet and cladding mechanical interaction (PCMI) occurs, solid contact conduction, $h_{\text{solid}}$, contributes to the gap heat transfer. This part is calculated using the modified Mikic-Todreas model considering the contact pressure and cladding Meyer hardness [16]. Radiation heat transfer, $h_{\text{rad}}$, is derived from the general radiation heat transfer equation based on the infinite-cylinder gray-body form [11]. $h_{\text{gap}}=h_{\text{gas}}+h_{\text{solid}}+h_{\text{rad}}$ (2)

Coolant heat transfer

A single fuel rod was used as the research object in this study; therefore, the coolant surrounding the fuel rod can be treated as a single-channel model. The temperature increased by convection between the cladding and coolant when it flowed axially. The coolant temperature was calculated based on the fluid energy equation in the axial direction. $\frac{\partial \left(\rho c_{\text{p}}T\right)}{\partial t}+\frac{\partial \left(\rho c_{\text{p}}Tw\right)}{\partial z}=q_v$ (3) where $q_v$ is the heat flux obtained from cladding surface; $w$ is the axial coolant velocity. The continuity equation for the coolant in the axial direction should be combined to complete the system. $\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho w\right)}{\partial z}=0$ (4)

The boundary condition was described by Newton's law in Eq.5. Considering that the research object is the fuel rod in LMFR, the heat transfer coefficient, $h_{\text{cool}}$, is calculated by the Mikityuk correlation [17], which is based on experiments of liquid metal in a bundle. $q_v=h_{\text{cool}}\left(T_{\text{co}}-T_{\text{c}}\right)$ (5) where $T_{\text{co}},T_{\text{c}}$ are the temperature of the cladding outer surface and coolant, respectively.

Parameters, such as thermal conductivity, heat capacity, and density are not constant; therefore, these equations must be solved by numerical iteration. The equations mentioned above were discretized using the finite volume or finite difference methods. The Gauss-Seidel iterative method was adopted to solve these equations over each mesh node to obtain the temperature distribution within the fuel rod.

Fission gas release analysis

The behavior of fission gas is essential for the fuel rod. The released fission gas may decrease the gap heat conduction, leading to a higher fuel pellet temperature. The plenum pressure may increase with burnup, which is an important boundary condition in the mechanical analysis of the fuel rod. Meanwhile, the fission gases in the matrix aggravate the swelling of the fuel pellet, thereby affecting the PCMI at a certain burnup.

To analyze the fission gas behaviors, much work has been done in the past decades. Booth et al. [18] first proposed an equivalent spherical model of the gas diffusion in a grain. Turnbull et al. [19] developed a diffusion coefficient model based on the experimental data. Speight et al. [20] introduced the effects of trapping and redissolving gas bubbles on gas atom diffusion. Based on previous work, Matthews and Wood [21],[22] proposed the operational gas release and swelling (OGRES) model, which physically describes the behavior of the fission gas produced in oxide nuclear fuel under operating conditions for fast reactors and has been incorporated into the KMC-Fueltra to simulate the transient behavior of fission gases. Based on a previous study [23], the URGAS algorithm, which was incorporated into the TRANSURANS code, was adopted to solve the production and release of fission gases.

Basic equations

Several assumptions were made to solve the complex mechanical problems. The first is that these deformations are small; second, the fuel rod is in quasi-static balance; and third, the planes perpendicular to the axial direction remain planar during deformation. The mechanical analysis of the fuel rod is controlled by mechanical principles as follows:

Geometric equation $\{\begin{array}{c}{\epsilon }_{\text{r}}=\frac{\text{d}u}{\text{d}r}\\ {\epsilon }_{\theta }=\frac{u}{r}\\ {\epsilon }_z=const\end{array}$ (6)

Equilibrium equation

The radial equilibrium equation is formulated as: $\frac{\text{d}{\sigma }_r}{\text{d}r}+\frac{{\sigma }_{\text{r}}-{\sigma }_{\theta }}{r}=0$ (7)

In addition to the radial equilibrium, the axial force balance for the fuel pellet is given by the integral equation: $\begin{array}{l}2\pi \underset{R_{\text{fi}}}{\overset{R_{\text{fo}}}{{\displaystyle \int }}}{\sigma }_z\left(r\right)\cdot r\text{d}r+\pi \left(R_{\text{fo}}^2-R_{\text{fi}}^2\right)\cdot P_{\text{plenum}}+{\displaystyle \sum _{k=i}^N\left(F_k+G_k\right)}\\ +Spr=0\end{array}$ (8)

Similarly, for the cladding: $\begin{array}{l}2\pi \underset{R_{\text{ci}}}{\overset{R_{\text{co}}}{{\displaystyle \int }}}{\sigma }_z\left(r\right)\cdot rdr+\pi R_{\text{co}}^2\left(P_{\text{c},\text{top}}+P_{\text{c},\text{bottom}}\right)+\pi R_{\text{ci}}^2\cdot P_{\text{plenum}}\\ +{\displaystyle \sum _{k=i}^N(F_k^{\text{'}}+G_k^{\text{'}})+Spr=0}.\end{array}$ (9)

Constitutive equation

The constitutive equation considering the in-pile behavior of the fuel rod is as follows: $\{\begin{array}{c}{\epsilon }_{\text{r}}=\frac{1}{E}\cdot \left[{\sigma }_{\text{r}}-v\left({\sigma }_{\theta }+{\sigma }_z\right)\right]+{\epsilon }_r^{\text{in}}\\ {\epsilon }_{\theta }=\frac{1}{E}\cdot \left[{\sigma }_{\theta }-v\left({\sigma }_r+{\sigma }_z\right)\right]+{\epsilon }_{\theta }^{\text{in}}\\ {\epsilon }_z=\frac{1}{E}\cdot \left[{\sigma }_z-v\left({\sigma }_r+{\sigma }_{\theta }\right)\right]+{\epsilon }_z^{\text{in}}\end{array}$ (10)

Eqs.11 and 12 provide the inelastic components in the constitutive equation. In addition to the elastic strain, thermal expansion, swelling, densification, relocation, plasticity, and thermal and irradiation creep strains for the fuel pellet and thermal expansion, swelling, plasticity, and thermal and irradiation creep strains for the cladding were considered. ${\epsilon }_{\text{p}}^{\text{in}}={\epsilon }_{\text{p}}^{\text{th}}+{\epsilon }_{\text{p}}^{\text{sw}}+{\epsilon }_{\text{p}}^{\text{den}}+{\epsilon }_{\text{p}}^{\text{relo}}+{\epsilon }_{\text{p}}^{\text{cr}}+{\epsilon }_{\text{p}}^{\text{plas}}$ (11) ${\epsilon }_{\text{c}}^{\text{in}}={\epsilon }_{\text{c}}^{\text{th}}+{\epsilon }_{\text{c}}^{\text{sw}}+{\epsilon }_{\text{c}}^{\text{cr}}+{\epsilon }_{\text{c}}^{\text{plas}}$ (12)

The fuel pellet plastic behavior was modeled by the ‘perfectly plastic' model, while the cladding was calculated using the Prandtl-Reuss flow rules: $\text{d}{\epsilon }_i^p=\frac{3S_i}{2{\sigma }_e}\text{d}{\overline{\epsilon }}^p,$ (13) where $E$ is Young's modulus; $v$ is Poisson's ratio; ${\epsilon }_r,{\epsilon }_{\theta },{\epsilon }_z$ are the total strain in radial, tangential and axial direction respectively; ${\epsilon }_r^{\text{in}},{\epsilon }_{\theta }^{\text{in}},{\epsilon }_z^{\text{in}}$ are the inelastic strain in radial, tangential and axial direction respectively; ${\sigma }_r,{\sigma }_{\theta },{\sigma }_z$ are the stress in radial, tangential and axial direction respectively; $u,w$ are the displacement in radial, and axial direction respectively; $R_{\text{fi}},R_{\text{fo}}$ are the inner and outer radius of the pellet; $R_{\text{ci}},R_{\text{co}}$ are the inner and outer radius of the cladding; $P_{\text{plenum}}$ is the plenum pressure; $F_k,F_k^{\text{"}}$ are the axial friction of the pellet and cladding respectively; $G_k,G_k^{\text{"}}$ are the gravity of the pellet and cladding respectively; $Spr$ is the spring force; $P_{\text{c},\text{top}},P_{\text{c},\text{bottom}}$ are the coolant pressure at the top and bottom respectively. ${\epsilon }_{\text{p}}^{\text{th}},{\epsilon }_{\text{p}}^{\text{sw}},{\epsilon }_{\text{p}}^{\text{den}},{\epsilon }_{\text{p}}^{\text{relo}},{\epsilon }_{\text{p}}^{\text{cr}},{\epsilon }_{\text{p}}^{\text{plas}}$ are the thermal, swelling, densification, relocation, creep, and plasticity strain of the fuel pellet respectively; ${\epsilon }_{\text{c}}^{\text{th}},{\epsilon }_{\text{c}}^{\text{sw}},{\epsilon }_{\text{c}}^{\text{cr}},{\epsilon }_{\text{c}}^{\text{plas}}$ are the thermal, swelling, creep and plasticity strain of the cladding respectively.

Cracking-healing model

Cracking is common in ceramic fuel pellets, which are inherently brittle even at low power. The low thermal conductivity induces high-temperature gradients, which leads to a high level of thermal stress and eventually cracks the fuel pellet into fragments. As long as cracking of the fuel pellet occurs, its mechanical properties can be modified as follows: $v_n={\left(\frac{1}{2}\right)}^n\cdot v$ (14) $E_n={\left(\frac{2}{3}\right)}^n\cdot E$ (15)

The criterion for cracking was that the fuel pellet cracks when the principal stress was larger than the fracture stress. The cracking number, $n$, in Eqs.14 and 15 increased by 1.

Cracking was considered to heal whenever the appropriate stress component was compressive at T >1400 °C and was applied continuously for at least 1 h, which is also used in LIFE [24]. If these cracks are healed, the mechanical properties of the fuel pellet will be restored to their original values.

Contact model

With increasing burnup, irradiation causes greater deformation of the fuel pellet and cladding, which results in a decrease in the gap width between them. PCMI occurs when the surfaces of the fuel pellet and cladding are in contact. This phenomenon is important for fuel rod performance because it can change the temperature distribution of the fuel rod by improving the gap heat conduction and changing the stress distribution as the boundary conditions. The plenum pressure was set as the boundary condition before the occurrence of PCMI, but it was later changed to the pellet-cladding interfacial pressure. The contact state between the pellet and the cladding is divided into three situations: no contact, soft contact, and hard contact. When there is no contact, PCMI does not occur and the contribution of friction is zero; when in soft contact, there is a relative slip between the pellet and the cladding and the friction between them is determined by $F_k=\mu \cdot F_N$; when in hard contact, PCMI occurs and there is no relative slip between the pellet and the cladding in the axial direction, and thus, $\cdot {\epsilon }_z^f=\cdot {\epsilon }_z^c$.

Based on the established equations, the stresses and deformations of the fuel rod are calculated with respect to each axial slice by imposing a continuous radial displacement and stress at the interface between mesh nodes and using the Gauss-Jordan elimination method to solve these equations by numerical iteration.

Thermal analysis

The thermal calculation of fuel rods is the basis of fission gas release and mechanical analyses. Analytical solutions and code results were chosen for comparison with the KMC-Fueltra. It is easy to find an analytical solution for steady-1D heat conduction by integration. For cladding, the inner surface temperature can be calculated using the following equation: $T_{\text{ci}}=T_{\text{co}}+\frac{q_l}{2\pi {\lambda }_{\text{c}}}ln\frac{R_{\text{co}}}{R_{\text{ci}}}$ (17) where $q_l$ is the linear power. A constant thermal conductivity was used in Eq.17 because it is not significantly dependent on the cladding temperature. However, the fuel pellets had the opposite effect. Therefore, the integrated thermal conductivity is used in Eq.18 to calculate the fuel pellet temperature at center, whereas the outer surface temperature can be obtained using Eq.19. Figure 4 shows the comparison results when the outer temperature of the cladding is set to 727.2 K and the linear power is set to 25.98 kW/m [30]. The maximum error from the analytical solution was approximately 10 K at the inner surface of the cladding, which was relatively small. $\underset{T_{\text{fo}}}{\overset{T_{\text{fi}}}{{\displaystyle \int }}}{\lambda }_{\text{f}}\text{d}T=\frac{q_l}{4\pi }$ (18) $T_{\text{fo}}=T_{\text{ci}}+\frac{q_l}{2\pi R_{\text{fo}}{h}_{\text{gap}}}$ (19)

Fig. 4Full-screen View

Comparison results of the fuel rod radial temperature

Simulation results from TRAC-AAA, developed by PSI, and EXCURS-M, developed by Japan Atomic Energy Research Institute (JAERI), which are widely used thermal hydraulics codes, were selected for comparison. An assumed beam interruption (BI) accident in an LBE-cooled and MOX-fueled accelerator-driven system has been simulated [30]. The axial power profile of the average fuel rod is shown in Fig. 5a. The power factor initially remained at one but dropped suddenly to zero at 1 ms. It later increased gradually at 6 s and then remained almost constant for 19 s. The history of the relative power is shown in Fig. 5b. The specifications of the fuel rod in the Preliminary Design Study of eXperimental Accelerator-Driven Systems (PDS-XADS) are listed in Table 1.

Fig. 5Full-screen View

Axial power profile (a) and transient relative power (b) during BI accident

Figure 6 shows the changes in the pellet center temperature and coolant outlet temperature. As can be seen from these figures, the temperature of the pellet and coolant changes with power during BI accidents, but the coolant temperature changes slightly later than the power because the process of heat transfer takes time. It can be concluded that the results of the simulations are consistent. These comparison results preliminarily illustrate the correctness of the numerical algorithm used in the thermal analysis, which provides a good foundation for other analysis modules.

Fig. 6Full-screen View

(a) Pellet center temperature and (b) coolant outlet temperature, during BI accident

Fission gas release analysis

Fission gas release behavior is essential for thermal and mechanical analysis by changing the gap conduction and acting as a boundary condition. The experimental data obtained from JOYO and the simulation results obtained by FEAST-OXIDE were selected to validate the correctness of this module. Experimental fast reactor of Japan, JOYO, was constructed to improve fast reactor safety and test advanced fuels and materials [31]. They have different core designs at different times. The MK-I core was critical in 1977, with an initial thermal power of 50 MW, which increased to 75 MW in 1979. The reactor core was then upgraded to MK-II in 1982, and its thermal power was increased to 100 MW. Table 2 lists the fuel rod specifications for cores MK-I and MK-II.

Figure 7 shows a comparison of the changes in the fission gas release fraction with burnup. It can be observed that the calculation results by KMC-Fueltra during the operation are in good agreement with most of the experimental data obtained from MK-Ⅰ and MK-Ⅱ and the calculation results by FEAST-OXIDE. However, there was still a slight difference on some occasions. It was concluded that these two factors could account for these differences. First, the joint oxide gain, which has not been considered in KMC-Fueltra, significantly affects the thermal conductivity of the fuel and results in a high release of fission gases. Second, the amount of fission gas released cannot be directly measured. It is calculated using other directly measured parameters, which means that the data itself has some uncertainties. The results from KMC-Fueltra agreed well with most of the experimental data within 6 at% burnup. In general, the accuracy of the fission gas release analysis module was valid within this range.

Fig. 7Full-screen View

Fission gas release fraction of JOYO MK-Ⅰ (a) and JOYO MK-Ⅱ (b)