The core idea of the XFEM is to deal with discontinuities in the computational domain with an extension function of the discontinuous property shape function base.

In ABAQUS software, the basic principle of the XFEM is based on the method of unit decomposition, which adds some special functions to the displacement function, reflecting the existence of discontinuity. The approximate displacement vector function based on the unit decomposition expansion can be expressed as Eq. 1 (Chen, 2013): u ¯ = ∑ I = 1 N N I x u I + H x a I + ∑ α = 1 4 F α x b I α , (1) where N I x is the displacement function of the ordinary node, u I is the continuous part of the displacement, a I and b I α are the vectors of the extended degrees of freedom of the node, H x is the interrupted jump function, and F α x is the asymptotic function.

In ABAQUS software, the XFEM crack propagation analysis is achieved by defining a damage model of the material. It consists of the initiation criterion and the evolution law. The former is used to determine whether the XFEM cracking mechanism is activated, and the latter is used to determine whether the crack is formally formed.

The damage initiation criterion is represented by an indicator value f , and the XFEM cracking mechanism is activated when Eq. 2 is satisfied: 1.0 ≤ f ≤ 1.0 + f tol , (2) where f tol is used to adjust the incremental step time to ensure that the damage criterion can meet the condition within a certain incremental step, generally 0.05. ABAQUS software provides six damage initiation criteria. In this paper, the maximum principal stress criterion is adopted as shown in Eq. 3 and 4: f = σ max σ max 0 , (3) σ max = σ max , σ max > 0 0 , σ max ≤ 0 , (4) where σ max is the maximum principal stress and σ max 0 is regarded as the tensile strength, indicating that there will be no cracks under compression.

When the criterion is satisfied, the material begins to enter the fracture stage. It introduces the damage variable D to represent the stiffness softening rate as shown in Eq. 5, 6 and 7: t n = 1 − D T n , T n ≥ 0 T n , T n < 0 , (5) t s = 1 − D T s , (6) t t = 1 − D T t , (7) where t n , t s , and t t are normal stress and two shear stresses after damage, and T n , T s , and T t are those before damage, respectively.

Energy-based linear damage evolution is considered in this paper. The damage variable D is shown in Eq. 8: D = G eff G effC , G eff < G effC 1 , G eff ≥ G effC , (8) where G eff is the equivalent energy release rate and G effC is the critical equivalent energy release rate. When G eff reaches G effC , the crack will be unstable and propagated.

According to Hooke’s law in general, in the linear elastic range of materials, the tensor form expresses the relationship between stress σ and elastic strain ε as Eq. 9, 10, 11 and 12 (Zhao et al., 2016): σ i j t = 2 G T , t ε i j e t + λ T , t ε k k e t δ i j , (9) G = E 2 1 + v , (10) λ = E v 1 + v 1 − 2 v , (11) δ i j = 0   i ≠ j 1   i = j , (12) where t is the time, G and λ are the Lame constants, E is the modulus of elasticity, v is the Poisson’s ratio, T is the temperature, and δ i j is the Kronecker function. Then, the stress increment will be expressed as Eq. 13: Δ σ i j = σ i j t + Δ t − σ i j t = 2 G T + Δ T , t + Δ t Δ ε i j e + λ T + Δ T , t + Δ t Δ ε k k e δ i j + 2 Δ G ε i j e t + Δ λ ε k k e t δ i j , (13) where elastic strain delta Δ ε e is expressed as Eq. 14: Δ ε e = Δ ε t − Δ ε th − Δ ε sw − Δ ε p , (14) where Δ ε t , Δ ε th , Δ ε sw , and Δ ε p represent the total, thermal, irradiation swelling, and plastic strain increment, respectively.

The thermal strain increment is expressed as Eq. 15: Δ ε th = ⁡ ln 1 + α T + Δ T T + Δ T − T 0 − ln 1 + α T T − T 0 δ i j . (15)

The irradiation swelling strain increment is expressed as Eq. 16: Δ ε sw = 1 3 ⁡ ln 1 + S W T + Δ T , t + Δ t − ln 1 + S W T , t δ i j . (16)

The plastic strain increment is expressed as Eq. 17: Δ ε i j p = 3 s i j t + Δ t 2 σ ¯ t + Δ t Δ ε ¯ p , (17) where s is the deviator stress, σ ¯ is the Mises stress, and Δ ε ¯ p is the increment of the equivalent plastic strain. The plastic strain increment is derived based on an iterative algorithm.

The above stress update algorithm and the variable properties of the material in Section 2.3 are implemented by programming the ABAQUS subroutine UMAT.

In this paper, the UMo-Al monolithic fuel element numbered L1P460 in the RERTR-12 experiment was selected as the research object, and its fuel was the U-10Mo alloy and the cladding was the Al6061-O alloy.

The results of the L1P460 annealing experiment after irradiation (Meyer, 2012) show nine blisters on the cladding surface, including six on the front and three on the back, and they are numbered from 1 to 9. The threshold temperature of the nine blisters is 673 K. The experimental results show that the blisters are mostly located at the junction. The cracking occurs when the UMo-Al monolithic fuel element fails. The failure process is that the fuel meat produces cracks and propagates to the cladding under the action of fission gas pressure during the annealing experiment. Therefore, it causes the blistering on the cladding surface.

Geometric modeling is based on the UMo-Al monolithic fuel element numbered L1P460 in the RERTR-12 experiment (Miller and Ozaltun, 2012). The fuel element geometric dimension is 101.473 mm × 25.4 mm × 1.397 mm, and the geometric dimensions of the UMo fuel meat is 82.55 mm × 19.05 mm × 0.254 mm, and its structure is shown in Figure 1. The XFEM attribute is given to the No. 4 blistering area, and the mini-crack of 1 mm² is set at the junction. Considering the amount of calculation and the convergence, the plasticity is only set in the area near the No. 4 blistering on the cladding surface.

In order to calculate the physical properties changing with irradiation of U10-Mo fuel, the burn-up of each node needs to be obtained. In this paper, based on the discrete burn-up of L1P460 after irradiation of the ATR reactor in Perez et al. (2012), the interpolation algorithm is programming by FORTRAN to allocate the burn-up to each node, and Figure 2 shows the U10-Mo burn-up distribution contour.

The boundary conditions for numerical simulation are that the xy plane, the xz plane, and the yz plane fix the displacement in the z direction, y direction, and x direction respectively, and the rest are free boundary conditions. A model is established with the element type of an eight-node linear hexahedral element and the non-coordinated mode of C3D8I. To determine the number of elements, the mesh independence verification is expanded, as shown in Figure 3, the final number of selected grid elements is 40,824, and the number of nodes is 46,870. The model is shown in Figure 4.

The steps of the actual annealing experiment are to heat, hold, cool to room temperature, and visually inspect. If no failure occurs, the experimental sample is raised to a higher temperature, and the above steps are repeated until failure occurs. Therefore, the temperature is changed periodically.

In this paper, the process of cooling and then heating to a higher temperature is simplified. This process is accompanied by the unloading and reloading of the thermal loads. In this paper, the cooling process is omitted and directly heated to a higher temperature, and the final effect is the same as that before simplification. At the same time, the simplification will greatly reduce the calculation amount of numerical simulation.

For the convergence, the first 1 min was applied to the irradiation swelling strain, and no thermal load was applied. Then, the annealing experiment was simulated, starting from 273 K. The strategy of 100 K heating for 15 min and holding for 45 min was selected to raise the temperature to 673 K, and then, the rate of decreasing 100 K for 15 min was selected to reduce the temperature to 273 K. The temperature curve is shown in Figure 5.

Since the XFEM Analysis Module in ABAQUS is located in the static analysis step, it cannot be used in the temperature–displacement coupling analysis. For the failure behavior, the stress is affected by the temperature, but the temperature is weakly affected by the stress. This paper chooses thermo–mechanical sequential coupling, which is less accurate but more efficient. In this paper, we first performed the heat transfer analysis through the ABAQUS subroutine UMATHT and then used the temperature field as the predefined field in the XFEM analysis and combined it with the ABAQUS subroutine UMAT to numerically simulate the failure behavior.

During the failure process, the crack propagates with the action of fission gas pressure. However, no load can be applied directly to the crack surface at the CAE interface of ABAQUS. In order to apply the gas pressure load on the propagating crack surface, this paper achieves three functions by modifying the input file of the ABAQUS model: (1) real-time creation of the XFEM crack surface, (2) application of the fission gas pressure load on the real-time XFEM crack surface and the amplitude of which is calculated by programming the ABAQUS subroutine UAMP, and (3) field output XFEM crack opening displacement, that is, the crack height.

For the calculation of the fission gas pressure load, the yield of fission gas (Burkes et al., 2016b) of UMo fuel is considered to be 26%, mainly the inert gas mixture of Xe and Kr, accounting for 93% and 7%, respectively. The gas equation of state adopts the ideal gas equation of state as shown in Eq. 27: P = n R T V , (27) where P is the fission gas pressure, n is the fission gas release, R is the gas constant, and V is the crack volume. Considering that when the crack passes through the fuel, all the fission gas in the thickness direction is released into the crack, then the molar amount of gas released by the crack through the fuel per unit area q can be expressed as Eq. 28: q = 0.26 × f d h f N 0 , (28) where N 0 is Avogadro’s constant, h f is the fuel meat thickness, and f d = 2.8 × 10 21  fissions / cm 3 for the No.4 blistering (Meyer, 2012). Then, the fission gas release will be expressed as Eq. 29: n = q ⋅ S , (29) where S is the crack area, which obtains Eq. 30: P = q R T ¯ h ¯ . (30) In other words, the fission gas pressure load P is related to the average crack opening height h ¯ and the average gas temperature T ¯ . The gas pressure load amplitude can be defined by programming the ABAQUS subroutine UAMP, and the fission gas pressure and crack propagation can be coupled and calculated.

For the average gas temperature, the historical output request is selected, the temperature is obtained in the subroutine UAMP through the GETSENSORVALUE function, and the arithmetic mean value is taken as the average gas temperature.

Since the GETSENSORVALUE function in the subroutine UAMP can only obtain the historical output and the crack height obtained by modifying the input file is the field output, it cannot be directly obtained in the subroutine UAMP. In order to expand the functions of ABAQUS software, this paper completes the calculation of fission gas pressure and realizes the coupling of fission gas pressure and crack propagation by programming Python to secondary develop ABAQUS with the FORTRAN subroutine and combining with the restart analysis of ABAQUS. Figure 6 shows the Python program flowchart.

Due to the good thermal physical properties, the temperature of each node of the UMo-Al monolithic fuel element calculated by the ABAQUS subroutine UMATHT heat transfer is basically the same when the external heat source is heated.

The maximum height (z-direction) displacement and temperature curves of the fuel meat and cladding surfaces obtained by numerical simulation are shown in Figure 7A. Before the temperature rises to 573 K, the maximum height displacements increase slowly with the increase in temperature during the heating period, while the maximum height displacements remain basically unchanged during the holding period, indicating that the mechanical load at this time is mainly thermal stress, the strain is mainly thermal expansion strain, and the height displacement is determined by the thermal expansion characteristics. Before 573 K, the fuel is always in the elastic stage. The maximum height displacement of the cladding is greater because the thermal expansion coefficient of the cladding is larger and the elastic modulus is smaller.

In the heating stage from 573 K to 673 K, it is found that the increase in the maximum height displacement is larger than before. The state of ① in Figure 7B is just heated up to 673 K. At this time, the maximum height displacement of the cladding shell is 0.01601 mm, and the state of crack propagation corresponds to Figure 8B. Contrasting with the initial state of the crack in Figure 8A, the increase in the magnitude becomes larger due to the beginning of crack propagation and the beginning of the action of fission gas pressure.

In the first 1,667 s of 673 K insulation, the crack basically did not propagate and the height displacement of the fuel was basically unchanged. After 1,667 s, there are two soaring of fuel height displacement. As shown in state ② in Figure 7B, the height displacement reached 0.04251 mm at 13426.99902 s, and the crack propagation state corresponds to Figure 8C. At this time, the crack further propagated, causing the height displacement to soar. As shown in the ③ state in Figure 7B, the height displacement reaches 0.04472 mm at 13501.95020 s, and the crack propagation state corresponds to Figure 8D, at which the crack shape is consistent, with the final shape of the crack at the end of annealing, as shown in Figure 8E, that is, at this time, the crack propagates sufficiently, and the height displacement soars and remains unchanged in the subsequent holding process.

The temperature at full crack propagation is designated as the blistering threshold temperature. The simulation results in this paper show that fuel failure occurs when the UMo-Al monolithic fuel is held at 673 K for approximately 1,742 s (∼29 min). The threshold temperature of 673 K obtained in this paper is consistent with the threshold temperature of 673 K for the L1P460 test in RERTR-12. The annealing test is difficult to visualize, and it is difficult to know which stage of blistering occurs during heating, holding, and cooling. The simulation results in this paper indicate that the failure occurs during the holding process and that the holding process cannot be omitted from the annealing test. It is also shown that the failure is related to the temperature duration.

In the cooling stage from 673 K to 273 K, the maximum height displacements decrease overall with the decrease in temperature but do not decrease to 0 since the fuel still has a crack with fission gas and has swelled. However, there is a height displacement rise, which is due to the negative feedback effect of a decrease in crack height, which, in turn, leads to an increase in fission gas pressure.

Figure 9A shows the crack opening displacement at the end of annealing, that is, the crack height. The maximum crack height at this point is 0.03162 mm. Figure 9B shows the height displacement of the crack at the end of annealing, where the crack height is the height displacement of the upper surface of the crack minus the height displacement of the lower surface. It can be seen that the crack propagates in all directions and basically only the center region bulges upward. Figure 9C shows a 2D sample of the fuel at the end of the annealing at the cracked surface, which is consistent in shape when compared to the experimental optical image (Meyer, 2012). Since the annealing experiment can only obtain the failure characteristics after cooling, the blistering height of the cladding surface obtained by the experiment is not the maximum height displacement. The height characteristics obtained by the experiment are not conservative enough. The maximum height displacement of the cladding surface can be obtained by numerical simulation, and the maximum height displacement and final height displacement of the cladding surface are 0.04472 mm and 0.02216 mm, as shown in Figures 10A, B, respectively, and the results are more conservation.

The blistering position is located at the junction, consistent with the crack position, which is mainly due to the fact that the fracture performance of Al6061-O is weaker than that of U-10Mo. Although the crack is generated in the fuel meat, the propagation direction is toward the cladding. In order to inhibit fuel failure, on the premise of not losing a large amount of heat transfer performance, the cladding material with better fracture performance can be considered instead of Al6061-O.

The maximum Mises stress, fission gas pressure, and temperature curves of the cladding surface are shown in Figure 11A. Before crack propagation, the Mises stress increases slowly with the increase in temperature and is basically unchanged during the holding period. The fission gas pressure does not act on the whole process. As shown in Figure 11B, with the crack just starting to propagate, the fission gas pressure starts to act, when the crack height is very small, leading to a rapid rise in the fission gas pressure. Before the crack propagates to a certain state, the Mises stress on the cladding surface rises slowly with the crack propagation. When the crack propagates to a certain state, the Mises stress on the surface of the clamshell rises rapidly. When the crack is fully propagated, the Mises stress on the cladding surface reaches the maximum value.

The left side (A, C, E) of Figure 12 shows the fuel meat stress distribution, and the right side (B, D, F) shows the cladding surface stress distribution. Figures 12A, B show the effects of irradiation swelling. Since only the fuel meat is swollen by irradiation, the stress in the fuel meat and cladding are concentrated in the areas of high burnup (long edges and corners, as shown in Figure 2. Figures 12C, D show the state before crack propagation, with the effect of thermal expansion superimposed on the irradiated swelling. Due to the different coefficients of thermal expansion of the fuel core and cladding, the stress concentration at the contact surface is caused. Figures 12E, F show the effect of crack expansion. For the fuel meat, the stress is concentrated at the tip of the crack propagation. For the cladding, there is more pressure load from the fission gas in the crack, so the stress is concentrated at the surface blistering.

Figures 13A, B show the Mises stress distribution at the blistering of the cladding surface at the pre-crack propagation and post-crack propagation stages, respectively. In the pre-crack propagation stage, the stress is concentrated at the highest point of the bulge. In the post-crack propagation stage, the stress concentration point is transferred to the contact position between the bottom surface of the bulge and the cladding. Therefore, during the actual operation of the reactor, the highest point of the bulge is most likely to rupture if the temperature rise transient, and during normal operation after blistering, the bottom surface of the bulge is most likely to rupture.

Based on the height displacement distribution at the blistering on the cladding surface in Figure 14, the blistering is approximated to be spherical crowns with a diameter of 2.4274 mm, and the blister area is calculated to be 4.63 mm². The experimental result for L1P460 in RERTR-12 is 4 mm² (Meyer, 2012). The simulation is slightly larger than the experiment, possibly due to the conservative fission gas release assumptions.