1. Introduction

Spent nuclear fuel (SNF) from nuclear reactors is stored for years to decades in storage pools for cooling prior to reprocessing or direct disposal [1]. As of the third quarter of 2024, the total amount of SNF stored in storage pools was 151,897 bundles according to data published by Korea Hydro & Nuclear Power (Fig. 1), which amounts to 89% of the total storage pool capacity of 170,207 bundles [2].

Fig. 1

Spent nuclear fuel storage status in wet storage in Q3 2024 [2].

SNF continues to be generated today, and the saturation of the storage pools at reactor sites is expected in near future. Thus, the preparations are underway for the deployment of interim storage facilities for SNF in Korea and the transportation of SNF from reactor sites to the interim storage facility seems an essential management option [3]. Maintaining the integrity of the cladding which encapsulates highly radioactive and toxic irradiated fuel is critical because the cladding acts as the first barrier to prevent release of radionuclides and maintains the geometric configuration of SNF [4-6]. The spent fuel rods can be exposed to various mechanical and thermal loads during the transportation. For example, the impact during a horizontal drop can significantly compromise the integrity of the spent fuel rods, creating a radiological risk. Thus, a comprehensive risk assessment is necessary, with the evaluation of the fuel damage ratio (FDR) under impact loads serving as a critical component of the process.

However, the structural evaluation of SNF assemblies and rods using CAE tools is challenging due to their complex characteristics. First, SNF assemblies and rods have highly complex geometries. A single transport cask can accommodate more than 30 fuel assemblies, each containing over 200 fuel rods, leading to a prohibitively large computational workload for detailed analysis. Second, the properties and characteristics of irradiated SNF are subject to significant uncertainties. The interfacial condition between cladding and pellets of a fuel rod is one source of the uncertainty. It is known that the complex interfacial interactions of the pellet and cladding have significant influence on the structural behavior of SNF rods. However, the mechanical properties governing these interactions remain poorly understood [7-12]. To address this challenges, simplified models of SNF rods and assemblies have been utilized in previous studies. In the work of Lee and Kim [13-17], the material properties of the equivalent beam model of SNF rods were found through optimization considering the effects of two pellet-cladding interaction (PCI) conditions. To evaluate the integrity of nuclear fuel rods in cask-level analyses, the Pacific Northwest National Laboratory (PNNL) and the Electric Power Research Institute (EPRI) have previously simplified nuclear fuel assembly into dummy shapes or single beam elements with equivalent mass and stiffness [18-20]. However, these simplified models do not account for the distinct characteristics of the structural parts within a fuel assembly, potentially reducing the accuracy of the analyses. Therefore, it is necessary to develop a simplified assembly model that considers the properties of each structural part during the simplification process.

Additionally, the process of quantifying the loads applied to the individual fuel rods in simplified fuel assembly models has not been extensively detailed in the literature. In this study, a detailed CE 16×16 fuel assembly finite element model was developed using the simplified fuel rod model developed by Kim and Lee [13-17]. The simplification for the assembly model was performed considering the stiffness of different parts within the fuel assembly as shown in Fig. 2. A procedure was developed to quantify the bending and pinch load during impact events using the proposed simplified fuel assembly model. Finally, the fuel damage ratio was evaluated by applying the failure criteria to the quantified bending and pinch loads.

Fig. 2

Flow chart of FDR calculation.

2. Development of Simplified SNF Assembly Model

2.1 Detailed Model of CE 16×16 Assembly

Fig. 3 shows a finite element model of a CE 16×16 fuel assembly containing 236 fuel rods and 5 guide tubes. Each fuel rod was modeled using a simplified equivalent beam model with the bonded interfacial condition between cladding and pellets developed in the study of Lee and Kim [13-15]. The fuel rod was simplified into a hollow beam with the same diameter with the actual fuel rod. The fuel rods and guide tubes are supported by 11 spacer grids modeled using shell elements. The springs and dimples in the unit cells of spacer grids are modeled as spring elements with different stiffness. At a spacer grid, each fuel rod in the assembly is connected to a total of four spring elements, consisting of two springs and two dimples. The top and bottom fixtures that hold the fuel rods and guide tubes in place were modeled using solid elements.

Fig. 3

CE 16×16 detailed model.

This CE 16×16 model was designed to quantify the bending load and pinch load exerted on the fuel rods due to impact acceleration. It will subsequently be simplified for cask-level analysis. This fuel assembly model in the form of CE 16×16 is referred to as detailed model in this study. Table 1 shows the mechanical properties for the main parts of the detailed model.

2.2 Equivalent Properties of Simplified Fuel Assembly Model

In this study, the detailed model introduced in the previous section was simplified into a composite structure composed of interconnected cuboids made of homogeneous elastic materials (Fig. 4). The process of determining the equivalent material properties for the cuboids corresponding to specific parts of the fuel assembly is also presented. Each cuboid was designed with external dimensions identical to the corresponding part of the fuel assembly to more accurately simulate the collision and contact behavior within the cask basket. The top and bottom fixtures were not simplified, while simplification was applied only to the regions containing spacer grids and the sections consisting solely of fuel rods.

Fig. 4

CE 16×16 detailed and simplified finite element models.

To obtain the equivalent mechanical properties of a particular section of the CE 16×16 detailed model, we performed structural analyses using ABAQUS/Explicit [21]. The spacer grid section includes fuel rods and guide tubes inserted into the spacer grid, while the fuel rod section consists solely of fuel rods and guide tubes. To simplify these sections into cuboids with homogeneous material properties, compression, shear, and torsion analyses were performed on two cross-sections of the detailed model shown in Figs. 5 and 6. The Young’s modulus of the fuel rods section in the simplified model was obtained by matching the flexural rigidity of the detailed model. The loading directions for the analyses mentioned above were determined to accurately reflect the resistance of the fuel assembly to the loads expected during the horizontal impact of the transport cask. The analysis results provided key physical properties, including elastic modulus, shear modulus, and Poisson’s ratio of the two sections.

Fig. 5

Compression and shear analyses of spacer grid assembly parts.

Fig. 6

Torsion analyses of fuel rod assembly parts.

Compression and shear analyses were conducted on the spacer grid section, as illustrated in Fig. 5. In both analyses, two rigid bodies are modeled on the top and bottom of the spacer grid and connected to the spacer grid using tie constraints. The bottom rigid body was then fully fixed and the top rigid body was subject to vertical and horizontal displacements. For the fuel rods section in Fig. 6, the shear analysis was replaced with a torsion analysis due to its geometry. The analysis similarly utilized rigid bodies with tie constraints applied to both ends of the fuel rod section. Torsional displacement conditions were then imposed on the rigid bodies.

Fig. 7 shows the nominal stress- strain curves obtained from the analyses of the spacer grid section. The nominal stress and nominal strain curves can be derived from the following equations:

Fig. 7

Stress-Strain curves for spacer grid section.

where σ is the nominal stress, F is the reaction force, and A₀ is the initial area of the cross section under load, ε is the nominal strain, L is the length of the section after deformation, and L₀ is the initial length before deformation. In both analyses, buckling is observed due to the structural features of the support grid section. Since the relationship between stress and strain is linear up to the critical stress point where buckling is observed, the modulus of elasticity was derived based on the point marked by the green circle in Fig. 7. Hook’s Law for deriving the modulus of elasticity is as follows:

where σ is the compressive stress, E is the elastic modulus, ε is the compressive strain, τ is the shear stress, G is the shear modulus of elasticity, and γ is the shear strain. Because the external dimensions of the simplified model were designed to match those of the detailed model, the calculated properties using Equation 3 and 4 can be regarded as the equivalent properties of the spacer grid section in the simplified model.

Fig. 8 shows the analysis results of the fuel rods section, showing the moment-angle curve obtained from the torsion analysis. The torsion analysis was performed to derive the equivalent shear properties of the fuel rod section. Through Hook’s law, the relation between the torsional angle and applied moment is given as in equation (5):

Fig. 8

Moment-Angle curve for fuel rods section.

where ϕ is the torsional angle and T is the applied torsional moment, L is the length of the axis, and I_P is the polar moments of inertia of the cross section of the fuel rods section in the simplified model. Since the cross-section of the fuel rod section has a square shape, the polar moment of inertia I_P can be calculated using the following equation:

where d is the width of the fuel rod section. From this relation, the effective shear modulus of the fuel rods section can be calculated. The elastic modulus of the fuel rod section can be derived from the condition that the flexural rigidity of this section should be identical in the detailed and simplified model. This relation can be formulated as follows:

where E₁I₁ is the flexural rigidity of the fuel rod section in detailed model and E₂I₂ is that of the simplified model. E₁I₁ can be determined by calculating the flexural rigidity of each component, including the fuel rods and guide tubes, that make up the cross-section of the detailed model. These individual contributions are summed using the parallel axis theorem with respect to the neutral axis of the entire crosssection. I₂ can be calculated using the following formula;

Then E₂ can be calculated using Eq. (7).

Finally, the Poisson’s ratio can be calculated using the obtained Young’s modulus and shear modulus, as expressed by the following formula:

where ν is the Poisson’s ratio. The equivalent material properties for the simplified models are summarized in Table 2.

2.3 Verification of Applicability of the Simplified Model

Before the simplified model developed in Section 2.2 can be used as a substitute for the CE 16×16 detailed model in drop impact analysis, its applicability must be verified. For this purpose, a finite element model of the transport cask, shown in Fig. 9, is utilized. The cask features an inner space capable of accommodating 21 nuclear fuel assemblies and includes a fuel basket and spacer disks to secure the assemblies. The cask body and the impact limiters attached to its top and bottom are also depicted in the figure. The material properties used for the cask model are summarized in Table 3. The detailed and simplified models of the fuel assembly were individually loaded into separate transport cask models, and horizontal drop impact analyses were performed from a height of 30 cm, as specified under normal transportation conditions in the regulations.

Fig. 9

Finite element model of transport cask.

To compare the dynamic behavior of the two models, structural responses such as the peak values of impact acceleration are compared because it is the dominating factor influencing the fuel rod failure. Acceleration data was extracted from the lower fixture of the fuel assembly models, a location less influenced by low-frequency vibrations due to its relatively stiff structural characteristics. The extracted acceleration is processed using the 4th order Butterworth filters to remove the high frequency signals generated by vibration and stress waves, leaving only the rigid body acceleration of the low frequency. The impact accelerations were extracted for three cases of loading position and the responses were filtered with three cutoff frequencies of 100, 200, and 300 Hz. Thus, the structural responses of the simplified assembly model were verified by comparing a total of nine results with those of the detailed model.

Fig. 10 shows the verification procedure, and Fig. 11 and Table 4 present the comparison of the impact accelerations measured at the bottom fixture of the detailed and simplified models, categorized by the cutoff frequency. The solid blue line shows the behavior of the CE16×16 detailed model and the dashed red line represents the behavior of the simplified model. In the results with the 100 Hz cutoff frequency, the impact accelerations of the two models show very close agreement. As the cutoff frequency increases, the discrepancies increase due to the influence of vibrations. In cases with relatively high error, the acceleration amplitude of the developed simplified model was observed to be higher than that of the detailed model. Based on the verification results, it can be concluded that the developed simplified model conservatively simulates the behavior of the CE 16×16 detailed model very closely and the developed simplified fuel assembly model can be utilized in place of the detailed model of the fuel assembly for extracting impact load for FDR calculations.

Fig. 10

Procedure for fuel assembly simplified model verification.

Fig. 11

Impact acceleration of the detailed and simplified models (blue line: detailed model, dashed red line: simplified model).

3. Procedure for Fuel Damage Ratio Calculation

3.1 9 m Drop Analysis Model

In this study, a total of 21 fuel assemblies were loaded in the cask to perform the drop impact analysis of the transport cask and to calculate the FDR. For efficient calculation, 20 simplified models and 1 CE 16×16 detailed model were used simultaneously to evaluate the FDR as shown in Fig. 12. The loading locations for the detailed model in the three cases are illustrated in Fig. 10 of Section 2.3. The drop impact analysis was performed based on the 9 m drop height of the hypothetical accident condition specified in the regulation. The 9 m height was converted to initial velocity and applied as initial condition for the dynamic simulation. To calculate the initial velocity from height, the following formula was used:

Fig. 12

Transport cask cross-section for 9-meter horizontal drop analysis.

where v is the impact velocity, g is the gravitational acceleration, and h is the drop height. The impacting surface was modeled as a rigid surface fixed in the space.

3.2 Quantification of Impact Loads

In order to evaluate the FDR, it is necessary to quantify the impact loads exerted on each fuel rod. In general, the loads experienced by a fuel rod in a drop impact can be categorized into two types. Bending loads due to the effects of inertia during the deceleration, as shown in Fig. 13, and pinch loads caused by the contact of the fuel rods with the other internal structure and fuel rods.

Fig. 13

Bending and pinch loads generated by impact loads.

The failure of fuel rods due to bending loads can be evaluated based on the membrane plus bending stress generated in beam model representing the fuel rod as demonstrated in the work of Lee and Kim [13-14]. Unlike the bending loads, which can be evaluated relatively simply through the stress response of beam, pinch loads can occur randomly at any location on the fuel rods, limiting their measurement by conventional methods. Pinch loads occur primarily in three cases:

For cases 1 and 2, the forces can be characterized as contact forces that occur at nodal points, primarily due to collisions between elements. In case 3, the forces are characterized as compressive or spring forces resulting from the compression of the spacer grid springs and dimples. To quantitatively extract these randomly occurring pinch loads, the Abaqus python scripting functionality was utilized.

The built-in python scripting in Abaqus enables the automatic extraction of the forces exerted by specific elements or nodes at every time increment. The python script we developed, shown in Fig. 14, allows us to extract the maximum contact forces and compressive spring forces acting on all specified fuel rods. Because the failure criterion for pinch loads on irradiated cladding is still under development, the extreme load values of unirradiated cladding, measured in the ring compression test, were used as an assumed failure criterion [22].

Fig. 14

Example of pinch load extraction using Python script code.

4. Results

4.1 Failure due to Bending Loads

The FDR evaluation was conducted for the three cases shown in Fig. 10 in Section 2.3. In the analysis, the detailed model was positioned at the red-marked locations, while the simplified model was placed at the unmarked locations.

The failure criterion due to bending loads is the membrane + bending stress of 749 MPa, calculated through the thickness of the cladding. Detailed explanation on this failure criterion can be found in [13-17]. The maximum von- Mises equivalent stress on the cross section of the beam is equivalent to the membrane plus bending component of von-Mises stress and it is multiplied by a stress correction factor of 1.78 to determine the failure of fuel rods [13-14]. The stress correction factor was introduced to reflect the stress concentration within the cladding due to pellet-cladding interaction.

Fig. 16 shows the results of the drop analysis for three different stacking positions, where the x-axis represents the identification numbers (Fig. 15) of the individual fuel rods inside the detailed fuel assembly and the y-axis represents the maximum von Mises stress experienced by the individual fuel rod multiplied by the stress correction factor. Different colors distinguish the three cases: Case 1 places the detailed fuel assembly at the bottom center of the cask, Case 2 at the center, and Case 3 at the top center. The results indicate that fuel rod failure due to bending loads was not observed in any of the three cases. Among the cases, the highest stress levels were experienced in Case 3, followed by Case 2, and then Case 1. The low stress values due to the bending load, despite the drop from a height of 9 meters, can be explained by the internal energy distribution shown in Fig. 17. This figure illustrates that the impact limiter absorbs more than 95% of the internal energy generated by impact, which is a typical characteristic of the design of type B transport cask.

Fig. 15

Load numbers for individual fuel rods in the detailed model.

Fig. 16

Maximum stress experienced during impact at all fuel rods.

Fig. 17

Internal energy distribution during 9 m drops impact.

4.2 Failure due to Pinch Loads

The failure criterion for pinch loads is still under development. Therefore, a provisional failure criterion of approximately 22,000 N was selected, based on the extreme load values observed in Ring Compression Tests (RCT) of unirradiated Zircaloy cladding [23]. Fig. 18 shows the quantified pinch loads applied to each fuel rod in the three cases. The pinch loads caused by the interaction with spring and dimple in spacer grid was not considered because their magnitude was much smaller than those of pinch load due to contacts with other fuel rods or basket. In Fig. 18, the x-axis of the graph represents the identification number of individual fuel rod (Fig. 15) within the detailed fuel assembly, and the y-axis represents the maximum contact force acting on the individual fuel rods. Significant contact forces were observed at the bottom, center and top of the cask, in ascending order. Correspondingly, the number of fuel rods exceeding the fracture criterion increased to 26, 37, and 45, respectively. Fig. 19 depicts various types of contact forces acting on the fuel rod derived from the drop analysis. In most cases, the largest contact forces resulted from collisions between fuel rods.

Fig. 18

Evaluating failure under contact forces.

Fig. 19

Types of contact forces exerting on fuel rods.

4.3 Calculation of Fuel Damage Ratio and Verification

The bending and pinch loads on each fuel rod were found to be the highest at the top-center position within the cask, referred to as Case 3. To confirm these results, a comparative analysis of key structural responses, including impact acceleration and internal energy, was performed for each case.

Fig. 20 shows the impact acceleration extracted from the bottom fixture filtered with a 300 Hz low-pass filter. The analysis revealed that the maximum acceleration occurred at the top center position, consistent with the trend observed during the quantification of pinch load acting on the fuel rods. Fig. 21 shows the internal energy absorbed by the detailed assembly model. It was observed that the fuel assembly at the top-center position absorbed the highest internal energy, which is consistent with the trends observed earlier.

Fig. 20

Impact acceleration by load position.

Fig. 21

Internal energy by load position.

Table 5 summarizes the fuel rod failure criteria and the corresponding number of failed fuel rods under impact load for each case. Out of a total of 236 fuel rods, the calculated damage ratios were 11.0%, 15.7%, and 19.1% for Case 1, 2 and 3, respectively.

One important point to note is that the accuracy of the pinch loads extracted using the Python script requires validation. The pinch loads analyzed in this study are contact forces occurring between beam elements or between beam and shell elements, which are known to be highly susceptible to numerical errors. Furthermore, it was observed that the magnitude of the calculated contact forces is significantly influenced by the size of the elements used. This issue necessitates further investigation and refinement in future studies. In this study, however, the focus is placed on presenting the complete evaluation procedures required for FDR calculation.

5. Discussion and Conclusions

This study aimed to evaluate the FDR in the event of a horizontal drop accident during transportation. To efficiently perform cask-level drop simulations, a simplified fuel rod model, previously developed by Kim and Lee [13-17], was employed. This model was used to create a detailed representation of a CE 16×16 fuel assembly, which was further simplified by considering its mass and stiffness properties. A procedure for calculating the FDR under bending and pinch loads generated by impact deceleration was proposed, with the following key findings:

The calculated FDR in this study may be very conservative, as the fuel rod model utilized beam elements which is prone to overestimating the contact force, and the pinch load failure criteria were based on un-irradiated cladding test results. Future research will address the refinement of pinch load failure criteria and the mesh resolution of beam elements to enable more realistic FDR assessments. Furthermore, to evaluate the FDR at all loading positions, we plan to directly apply position-specific acceleration data obtained from the simplified models in transportation cask drop analysis as acceleration loads to the detailed CE 16×16 model. This approach is expected to enhance the time efficiency of FDR calculations, contributing to more accurate and practical assessments.

Nomenclature