2.1.1 Thermal conductivity

The thermal conductivity is given as a function of temperature [17],

where k in W/(K⋅m) is the coefficient of thermal conductivity, T is the temperature in K.

2.1.2 Thermal expansion

Ozaltun et al. listed a series of data for the coefficient of thermal expansion [17], the correlation between the coefficient of thermal expansion α^T and temperature T was fitted as [15]

Logarithmic strains are adopted in this study considering that large strains may occur in the UMo meat at higher burnup. The logarithmic thermal expansion strain components can be expressed as

where T₀ is the reference temperature with a value of 300 K, δ_ij denotes the Kronecker delta.

2.1.3 Irradiation-induced creep

According to the creep constitutive theory [17], creep strain components can be calculated as

where, s_ij depicts the deviatoric stress components at the current time, σ¯=32sijsij is the Mises stress, εijcr denotes the creep strain components, ε̄^cr is the equivalent creep strain.

For Eq.(4), using the trapezoidal integration for a time increment, we have

which indicates that the creep strain increments depend on the stresses at time t and time t + Δt.

The irradiation-induced uniaxial creep strain rate for UMo is presented as [7]

where, B is the coefficient of creep rate equaling 300×10^–37 m³/Pa; σ depicts the uniaxial stress in Pa; ḟ denotes the fission rate in fission/(m³s), which means the nuclear fission number per second happens in a unit volume of fuel meat.

Similarly, the trapezoidal integration for Eq. (6) yields

Substituting Eq.(7) into Eq.(5) yields

Obviously,ΔΔεkkcr = 0.

2.1.4 Irradiation-induced swelling

Irradiation-induced swelling stems from fission gas swelling and fission solid swelling.

Here, an empirical equation to calculate the fission solid swelling was used [6]:

where f in fission/m³ describes the fission density, which describes the nuclear fission number occurred in a unit volume of fuel meat.

The fission gas model follows Yi Cui et al. [12] considering the recrystallization, hydrostatic pressure dependence and intergranular gas atom re-solution in grain scale. The swelling model in this study is based on the Booth model [18], in which the fuel grain is regarded as a spherical one. Fission gas swelling includes two parts, i.e., intragranular swelling and intergranular swelling.

Before recrystallization, the fission gas swelling can be calculated by the following equation:

where the intragranular swelling can be expressed as [19]:

where, r_b in m is the radius of the intragranular bubble; c_b in n/m³ is the fission gas bubble density, representing the number of fission gas bubbles in a unit volume of fuel meat. How to obtain the variables of r_b and c_b was presented in [12], which were calculated with the concentration of intergranular gas atoms.

The intergranular swelling takes the form as [20]:

where, C_b in n/m² is the intergranular bubble density, depicting the number of intergranular bubbles in a unit area of grain boundary; R_b in m is the radius of the intergranular bubble; r_gr in m is the grain radius.

Considering the influence of external hydrostatic pressure, the modified Van der Waals gas law is used [20]:

where, γ is the surface tension in N/m, p is the external hydrostatic pressure in Pa, N_b = N/C_b is the gas atom number per intergranular bubble and N depicts the concentration of intergranular fission gas atoms, h_s is the fitting parameter to make the Van der Waals equation equivalent to the hard-sphere equation of state [19], b_v is the Van der Waals constant for Xe. N is calculated with the analytical solution [12], which was obtained through solving the governing equations for diffusion of fission gas atoms in the equivalent spherical grain, where re-solution of the intergranular fission gas atoms is considered.

So,

with which the hydrostatic pressure can be solved through nonlinear iterations, as stated in the next section.

Considering Eq. (13) is a nonlinear equation, for a certain magnitude of p the Newton’s method to get R_b should be performed.

After recrystallization, the original grain is divided into two parts including the recrystallized region and the unrecrystallized region. In the unrecrystallized region, the calculation of the gas swelling is similar to the calculation before recrystallization; while in the recrystallized region, it depends only on the intergranular swelling, because almost all fission gas atoms will transfer to the intergranular bubbles. The intergranular swelling can be expressed as [20]:

where, r_grx is the recrystallized grain radius, R_bx is the radius of the intergranular bubble in the recrystallized region; C_bx is the intergranular bubble density [12].

Similar to that for R_b, according to the modified Van der Waals gas law, ∂Rbx∂p can be obtained as:

with which the hydrostatic pressure can be calculated through nonlinear iterations, as described in the next section.

Subsequently, the fission gas swelling can be calculated according to

where, V_r is the volume fraction of the recrystallized region [21].

So, the volumetric swelling strain can be expressed as:

Then, the logarithmic volumetric swelling strain can be obtained as:

Obviously, the logarithmic swelling strains are obtained as

which depicts that the irradiation swelling strain in all directions is the same.

Some of used parameters are listed in Table.1, while the other parameters can be found in Ref. [12].

2.1.5 Heat generation rate

The heat generation in the UMo meat is proportional to the fission rate:

where, q is the heat generation rate, c = 3.204 × 10^–11 J/fission is the generated heat energy for every fission event.

2.2.1 Thermal conductivity

The thermal conductivity is related to temperature [17], given as

where 273 K ≤ T ≤ 854 K.

2.2.2 Thermal expansion

Based on some experimental data [17], we have the coefficient of thermal expansion α^T by polynomial fitting:

where 373 K ≤ T ≤ 573 K.

2.2.3 The plasticity model

The strain-hardening curve is described as:

where, σ in MPa is the true stress, ε is the true strain, the strength coefficient K₀ = 569.6 MPa [7], the strain-hardening exponent n = 0.13 [7].

3.1.1 Stress update algorithm for the spherical part

Firstly, to update the spherical stresses, it’s necessary to set i = j in Eq. (29), and then multiply both sides by –1/3, we have

where, K = 2G3 is the bulk modulus.

Substitution of Eq. (27) and Eq. (28) into Eq. (30) yields

where, the known θtsw = θ^sw(t,T,p^t), the unknown logarithmic volumetric strain at time t + Δtθ^sw depends on the unkown hydrostatic pressure, the temperature and the irradiation time as θ^sw(p^t+Δt,T + ΔT,t + Δ t), and p=−σkkt+Δt3 depicts the hydrostatic pressure to be determined.

Eq. (31) is a nonlinear equation that the logarithmic volumetric swelling strain and the hydrostatic pressure satisfy. Using the hydrostatic-pressure-dependent equations for the logarithmic volumetric swelling strain, the Newton-Raphson iteration method is firstly adopted to get p:

where, p_i is the value of p after the i^th iteration. And it should be noted that θ^sw in (p_i) should be obtained with the iterative solution of Eq. (13) to get R_b for a certain value of p_i.

In Eq. (33), ∂θsw∂p can be obtained by implementing partial derivative operations with respect to p for the both sides of Eq. (19) together with Eq. (18), with the form described as

As given in Section 2.1.4, only the intergranular swelling is related to p. Before recrystallization, by combining Eq. (12) with Eq. (14) the partial derivative ∂∂pΔVgasV0 can be obtained as

Then using the partial derivative ∂θsw∂p,p before recrystallization will be worked out and θ^sw can be determined simultaneously.

While during the recrystallization, the partial derivative ∂∂pΔVgasV0 can be derived by using Eq. (14), Eq. (15) and Eq. (16) as

Also, p and θ^sw during the recrystallization can be obtained through nonlinear iterations. Considering p=−σkkt+Δt3, so the spherical stresses at time t + Δt can be obtained.

3.1.2 Stress update algorithm for the deviatoric part

As for the deviatoric stresses sijt+Δt, we introduce the trial stresses σijpr which refers to the stresses without involving the contribution of irradiation creep strain increments, then we have

As is known that the creep strain tensor is a deviatoric tensor, so

where, sijpr is the deviatoric part of σijpr.

Substituting Eq. (8) into Eq. (39), we have

So,

where, sijpr can be determined by the obtained stresses without considering the irradiation creep strain increments. As mentioned above, the increments of thermal expansion and irradiation swelling will not contribute to the deviatoric stresses. Thus, the deviatoric stresses in the fuel meat can be simply updated.

As a result, the Cauchy stress tensor can be updated with the contributions of the spherical part and the deviatoric part.

The flow chart of the stress-update algorithm is shown in Fig. 1. In every time increment [t,t + Δt], the swelling strain is firstly calculated with the hydrostatic pressure at time t, after which in every iteration step the hydrostatic pressure p at time t + Δt is obtained according to Eq. (34). It is noted that for a certain hydrostatic pressure, the radius of intergranular bubbles is obtained with Newton iteration method based on Eq. (13). As depicted in the flow chart, when the irradiation swelling strain and hydrostatic pressure p are convergent, the deviatoric stress components at time t + Δt can be obtained according to Eq. (41), with the creep effect taken into account. Finally, σijt+Δt=sijt+Δt−pδij, the stress update is finished.

Figure 1

Flow chart of stress-update algorithm

4.3.1 Thermal boundary conditions

As illustrated in Fig. 2(b), the convection heat transfer boundary refers to the lower surface of the finite element model, which meets the convection boundary condition: – k∂T∂n = h(T – T_f), where the temperature of the coolant T_f = 323K. The heat transfer coefficient h is related to parameters such as the heat transfer area and Reynolds number [23,24,25]. In this study, the heat transfer coefficient h = 0.035W/mm²K;

The other surfaces satisfy the heat insulation condition, that is, – k∂T∂n = 0.

4.3.2 Mechanical boundary conditions

1. Two symmetrical boundary conditions are respectively set on the surfaces, shown in Fig. 2(b);

2. The side surfaces are fixed to restrain the rigid body movement, as shown in Fig. 2(b);

3. The other surfaces are set as free boundaries because the pressure of the coolant is small [26].

5.1.1 Validation of the developed models

To validate the developed models and algorithms, the calculated thickness variations for the UMo fuel meat are compared with the measured data from RERTR-9A experiment [7].

From Fig. 4(b) it can be seen that the obtained numerical results for the relative variations of fuel meat thickness matches compariably well with the experimental results. And compared with the predicted results based on the swelling empirical formulas [7], the numerical results in this study agree better with the experimental ones near the heavily irradiated edge. This is for the reason that large hydrostatic pressures undergo at these areas, which will lead to smaller local fission gas swelling. As a result, the developed models and algorithms in this study are effective.

Figure 4

a) The output path; b) validation of the simulated fuel swelling

5.1.2 Verification of UMAT and UMATHT

To verify the accuracy and effectiveness of UMAT and UMATHT, the calculation results in the UMo monolithic plate are analyzed.

According to the logarithmic elastic strain components calculated by ABAQUS, the “theoretical” Cauchy stresses can be calculated by using the elastic constitutive relation.

Comparing the “theoretical” Cauchy stresses with the numerical results of the Cauchy stresses outputted by ABAQUS, the elastic constitutive relation can be verified. Fig. 5 gives the stress component σ₁₁ on the 98^th day. One can see that the numerical results are in good agreement with the “theoretical” ones, in fact, all of the other stress components also obey the elastic constitutive relations, which demonstrate that the elastic constitutive relations in the meat and cladding are correctly defined in UMAT.

Figure 5

Verification of elastic constitutive relation

Using Eq. (2), the thermal expansion coefficient can be calculated according to the temperature T + ΔT, with which substituted into Eq. (3) the “theoretical” logarithmic thermal expansion strain along one direction can be determined; comparing it with the logarithmic thermal expansion strain output by ABAQUS, the correct definition of thermal expansion strain can be verified.

Fig. 6 shows the comparison between the “theoretical” results and the numerical results on the 98^th day along the output path in Fig. 6, which indicates that they fit well with each other. So, it can be concluded that the thermal expansion strains in the meat and cladding are correctly defined in UMAT.

Figure 6

The verification of thermal expansion

The heat flux obeys:

which is introduced into the FEM simulation by the subroutine UMATHT.

According to the temperature T + ΔT on the 11^th day, the thermal conductivity can be calculated; using the outputted ∂T∂y defined as a state variable in the UMATHT, then the “theoretical” heat flux q_y is calculated with Eq. (53); comparing it with the numerical results of the heat flux output by ABAQUS, the defined thermal conductivity will be verified.

The comparison between the heat flux obtained above by Eq. (53) and the numerical one on the 98^th day can be found in Fig. 7, which shows that these results agree well with each other. Consequently, the thermal conductivities of the fuel meat and cladding are well defined in UMATHT and the corresponding programs are written correctly.

Figure 7

Verification of heat flux

The correctness is confirmed for introduction of the hydrostatic-pressure-dependent irradiation swelling into the mechanical constitutive relation, as shown in Fig. 8. The numerical results for the swelling strains on the 98^th day are obtained from the output files of ABAQUS. The “theoretical” true swelling strain including the fission solid swelling and fission gas swelling is obtained as

Figure 8

Verification of irradiation swelling strains

while ΔVV can be obtained with Eq. (9), (11), (12), (15) and (17) based on the fem-obtained bubble radius. One can find that the “theoretical” results are consistent with the numerical results, and the swelling values are lower near the heavily irradiated edge due to the large hydrostatic pressure there, which can be found in Section 5.2.3.

The “theoretical” Mises stress σ̄^t+Δt can be determined by [22]:

where, ε¯pt+Δt is output by ABAQUS. The theoretical results of σ̄^t+Δt are calculated through nonlinear iterative solution of Eq. (55)

Fig. 9 gives the “theoretical” results of the Mises stress calculated by the equivalent elasto-plastic strains ε¯pt+Δt+σ¯t+ΔtE. One can see that the points locate at the hardening curve, which depicts that the Mises stresses and the equivalent elasto-plastic strain obey the functional relation of the hardening curve. It is noted that the investigated point is experiencing the plastic loading process. As a result, the plasticity model is verified to be correctly involved in the mechanical constitutive relation for the cladding.

Figure 9

Verification of plasticity constitutive relation

The “theoretical” equivalent creep strain ε¯crt+Δt is calculated using Eq. (7) with substitution of σ̄^t+Δt, σ̄^t, ḟ and ε¯crt output by ABAQUS, where the Mises stress σ¯=32sijsij. As shown in Section 3.1.2, the deviatoric stresses s_ij depends on the equivalent creep strain. The “theoretical” results on the 1.1^th day are compared with the numerical ones from ABAQUS, as depicted in Fig. 10. One can observe that both results are consistent with each other, thus the irradiation induced creep in the meat is well involved in UMAT.

Figure 10

Verification of fission-induced creep effect

Simultaneously, it is confirmed that the calculation of the deviatoric stresses is correct.

5.2.1 Temperature field

Fig. 11 shows the distribution of temperature on the 98^th day. It can be seen that the temperature in the fuel meat is much higher than that in the cladding, and the highest temperature of 487.9 K appears at the location near the heavily irradiated fuel edge. According to the experimental result, the temperature at the middle of the fuel is about 425 K [7], and in this simulation the predicted temperature is 423K, which matches well with the experimental result. The distribution and evolution of temperatures along Path 1 (depicted in Fig. 11) are plotted in Fig. 12. As is known, the maximum fission rate takes place at the heavily irradiated fuel edge corresponding to the starting point of Path 1. One can see clearly that the highest temperature exists in the position a distance away from the starting point of the Path 1, which results from the fact that the cladding has a very high thermal conductivity and can effectively remove the generated heat in the fuel meat near the meat/cladding interface. Besides, it can be observed that the highest temperature increases remarkably with the irradiation time. The highest temperature on the 98^th day achieves a magnitude about 61.7 K higher than that on the 1.1^th day. It is the main reason that the fuel meat around the peak temperature zone is thickened with the irradiation time, which is reflected by the strain results in Section 5.2.2.

Figure 11

The contour plot of temperature

Figure 12

The evolution of temperature along Path 1

5.2.2 The strain field

In this study, the deformation contributions of thermal expansion, irradiation creep and irradiation swelling are considered in the fuel meat. Meanwhile, in the cladding the total deformation consists of the thermal expansion and elastic-plastic strains.

As depicted in Fig. 13, the distribution characters of thermal expansion strains in the fuel meat and cladding are similar to that of the temperature. It can also be found that the maximum thermal expansion strain is larger in the cladding than that in the fuel meat, which is due to the greater coefficient of thermal expansion of the cladding. The differences in coefficients of thermal expansion together with the non-homogeneous temperature field and mechanical constraints result in thermal stresses in the fuel plate.

Figure 13

The contour plots of thermal expansion strains on the 98^th day: (a) in the fuel meat; (b) in the cladding

Fig. 14 provides the distributions of equivalent creep strains and swelling strains on the 98^th day, and Fig. 15 presents their evolution results along Path 2 (illustrated above in Fig. 10) with variation of the irradiation time. It should be noted the fuel meat in Fig. 2 has been rotated for 180 degrees around the x-axis here.

Figure 14

The contour plots of (a) equivalent creep strains and (b) irradiation swelling strains on the 98^th day in the fuel meat

Figure 15

The evolution of (a) equivalent creep strains and (b) irradiation swelling strains along Path 2

One can see from Fig. 14(a) that the equivalent creep strains are generally higher at the regions at a distance away from the interface between the fuel meat and cladding. The equivalent creep strains are much higher around the heavily irradiated edge, and the maximum equivalent creep strains reach a very large magnitude of 236.4% at the meat corner. It is known that the equivalent creep strains depend on the Mises stress, the fission rate and the irradiation time, as expressed in Eq. (7). The highest fission rate in the fuel meat appears at locations corresponding to the starting point of Path 2 (depicted in Fig. 14(b)). Hence, it can be obtained that the distribution pattern stems from the Mises stresses in the fuel meat. Fig. 15(a) clearly depicts that (1) the equivalent creep strains increase with irradiation time; in the vicinity of the path origin, the equivalent creep strain decreases first, then increases to a peak value. The distribution character of the equivalent creep strains along Path 2 can be explained by the distribution of the Mises stresses, as given in Section 5.2.3.

As shown in Fig. 14(b), under the heterogeneous irradiation condition, the maximum irradiation swelling strain of 15.79% exists near the heavily irradiated edge and decreases along the direction of Path 2 as a whole. Owing to consideration of the hydrostatic-pressure dependence of the fission gas swelling, the irradiation swelling strains at the region close to the length ends of the fuel meat are observed to be lower. This can be interpreted since the hydrostatic pressures there must be higher because of intense mechanical interactions between the fuel meat and the cladding. From Fig. 15(b), it can be clearly seen that the irradiation swelling strains increase with irradiation time; and the swelling strains increase more from the 48^th day to 98^th day around the path origin than the other regions, which is also related to grain recrystallization occurring earlier there to accelerate the fission gas swelling despite the high irradiation there.

Superposition of the contributions of thermal expansion, irradiation creep and irradiation swelling in the fuel meat results in existence of higher thickness increment at the location with a distance away from the heavily irradiated edge. Thus, a higher temperature can be interpreted there.

Fig. 16 displays the distribution of equivalent plastic strains on the 98^th day in the cladding. According to the plastic constitutive theory, the plastic strains are also coupled with the stresses. Besides, the equivalent plastic strains are dependent on the Mises stresses. One can find from Fig. 16 that the equivalent plastic strains are larger in the vicinity of the heavily irradiated region than the interior part, especially at the corner of the meat/cladding interface. This is easy to understand because of the more intense mechanical interaction there. Besides, the thickness of the fuel meat along Path 2 is distributed similarly to equivalent creep strains in Fig. 15(a); thus the locally thickened fuel meat intensifies the extrusion between the fuel meat and the cladding. Path 3 in Fig. 16 goes through the regions having larger equivalent plastic strains. Fig. 17 shows evolution of equivalent plastic strains along Path 3. It is shown that equivalent plastic strains increase with irradiation time; the maximum magnitudes appear near the path end; the distribution and evolution patterns stem from irradiation swelling and creep effects occurring in the fuel meat, as explained above.

Figure 16

The contour plot of equivalent plastic strains on the 98^th day

Figure 17

The evolution of equivalent plastic strains along Path 3

5.2.3 The stress field

The Mises stresses are coupled with the creep and plastic strains, and they can be used to judge whether the material fails. What’s more, the hydrostatic pressure in the fuel meat is considered to be associated with fission gas swelling in this research. So, the Mises stresses and hydrostatic pressure will be discussed in the following.

Fig. 18 provides the distributions of the Mises stresses in the fuel meat and cladding on the 98^th day. Fig. 19 presents their evolution histories along path 2 in the fuel meat and Path 3 in the cladding. It should be noted that the contour plot in Fig. 18(a) is obtained by rotating the fuel meat in Fig. 1 for 180 degrees around the x-axis.

Figure 18

The contour plots of the Mises stresses (a) in the fuel meat and (b) the cladding

Figure 19

The evolution of Mises stresses: (a) Mises stresses in path 2; (b) Mises stresses in path 3; (c)Normalized Mises stress in path 3.

Comparing the contour plot in Fig. 18(a) with that in Fig. 15(a), we see that the distribution of equivalent creep strains is similar to that of Mises stresses in the fuel meat. Due to the creep relaxation effect, the Mises stress on the 98^th day in the fuel meat are quite small, much smaller than that in the cladding. The enlarged irradiation swelling has a trend of giving rise to the Mises stress, while irradiation creep tends to lessen the Mises stress. They both contribute to the evolution of the Mises stresses. As depicted in Fig. 19(a), the highest Mises stress increases quickly from the 1.1^th day to the 18^th day. One can also observe that there are high stress points at the locations at a distance away from the heavily irradiated edge, which is consistent with the distribution of irradiation creep.

For the Mises stresses in the cladding, one can find from Fig. 18(b) that comparably large Mises stresses mainly emerge in the interface between the fuel meat and the cladding, and the maximum magnitude on the 98^th day achieves 403.1 MPa. Fig. 19 shows that the Mises stresses increase with irradiation time, and the distribution is consistent with that of equivalent plastic strain in Fig. 17. Fig. 19(c) shows the Mises stresses normalized by 280 MPa, the initial yield stress of Al6061, which clearly depicts the strain hardening effect. On the 1.1^th day, the normalized Mises stresses are below 1.0, which means the material points have not met the yield condition.

In this study, the hydrostatic pressure differs significantly in the fuel plate, as shown in Fig. 20 which is obtained by rotating the fuel meat in Fig. 2 for 180 degrees around the x-axis. The edge effect is remarkable in Fig. 20(a), and the maximum value appears at the meat corner around the heavily irradiated region. This can explain the distribution pattern of swelling strains in Fig. 14(b). The edge effect of hydrostatic pressure is evidently reflected by the results in Fig. 20(b). One can find that the hydrostatic pressures are much higher near the two ends of Path 2 than those in the other locations, which is induced by enhanced mechanical interactions around the edges. With increasing irradiation time, the maximum hydrostatic pressure also rises, which tends to reduce fission gas swelling. The total irradiation swelling strain curve on the 98^th day in Fig. 15(b) has lower values in the vicinity of the path origin, which should result from enlarged hydrostatic pressure there.

Figure 20

(a) The contour plot of the hydrostatic pressure on the 98^th day in the fuel meat; (b) the evolution history of the hydrostatic pressures along Path 2