Thermomechanical Microstructural Predictions of Fracture Nucleation of Zircaloy-4 Alloys With δ and ɛ Hydride Distributions

Material systems comprised of zirconium alloys are used as cladding protection against the release of fission gases and other products during the operation of pressurized water reactors (PWR). This is mainly due to its favorable thermomechanical behavior, ductility, and the insignificant neutron absorption of these zirconium alloys [1]. The most commonly used alloys, for these harsh thermal environments, are Zircaloy-2, Zircaloy-4, and Zirlo (zirconium low oxidation).

However, due to thermal gradients across the cladding and the corrosion of zirconium alloys due to reaction with the coolant water, hydrogen can be transported to the zirconium alloys, and this can result in the precipitation of hydrides at the edge of zirconium-oxide layers. These hydrides generally have crystalline structures that are different from the zirconium alloy h.c.p. matrix and have an atomic structure of zirconium and hydrogen. These hydrides can initially form with circumferential orientations, but due to cycles of stress and thermal loading conditions, these hydrides can reorient to mostly radial distributions.

The orientation and distribution of hydride populations can significantly affect the ductile to brittle transition temperature (DBTT) of the parent zirconium alloy. The DBTT of hydrided zirconium alloys has been reported to range from 50 °C to 260 °C [2–7]. Billone et al. [8] conducted ring compression tests (RCT) with high burn-up Zircaloy-4 and Zirlo cladding specimens for temperatures ranging between 25 °C and 200 °C that indicated the transition temperatures significantly affected the hydride morphology, orientation, and crystallography. Hydride reorientation in the Zircaloy-4 can also take place during long-term storage, in the absence of a corrosive environment, and this can also significantly reduce the ductility of the cladding system [9–11]. Kim et al. [12–15] experimentally investigated how hydrides affected the cladding mechanical behavior for different conditions, such as those thermal conditions associated with the loss of coolant accident (LOCA). Their results show that hydrides in the circumferential direction at room temperature are generally brittle with some plasticity, but can be ductile above 150 °C. In contrast, for radial hydrides, at 150 °C, micro-cracks were nucleated and as they can act as paths for crack propagation between the different crystalline structures of the matrix and the hydrides.

Chan et al. [16] conducted micro-fracture experiments to determine the local mechanical properties of multi-crystalline systems with a zirconium/zirconium hydride microstructure. Their experiments clearly indicate that different interphase crystalline structures between the δ, ZrH_1.66 (f.c.c.) hydrides, the ɛ, ZrH₂, and the h.c.p. matrix affect crack nucleation and propagation. How these different interfaces and interphases collectively affect fracture for the hydride population is not clearly understood due to the complex crystalline interactions between the matrix and hydride populations. For example, there are interfacial dislocations and misfit strains, due to semicoherent interfaces between the matrix and the hydrides. [17]. Due to the different lattice parameters and atomic arrangements between the parent and hydride phases, interfacial planes can form and result in heterogeneous crystalline behavior. This is because interfaces, between the closed packed planes of the matrix and the hydrides, have misorientation dislocation and resulting misfit strains. Hence, orientations relations (ORs) have to be physically represented to account for both coherent and semicoherent interfaces [18,19].

Hence, our main goal is to postulate a microstructurally based validated approach that can accurately predict and understand the thermomechanical failure behavior of the hydrided microstructure. This approach is predicated on an accurate physical representation of distributions of circumferential and radial precipitated hydride morphologies, which have distinct crystalline systems than the matrix, such that a fundamental understanding of how multiple fracture failure modes and surfaces nucleate and propagate because of hydride orientations, morphologies, and distributions and matrix interactions with a multiphase system that is comprised of h.c.p. matrix and f.c.c. and b.c.c. hydrides. The multiple-slip crystal plasticity dislocation-density-based formulation and the ORs are briefly outlined in Sec. 2, the microstructural fracture criterion is discussed in Sec. 3, the numerical implementation of the fracture approach is given in Sec. 4, the overall computational approach is given in Sec. 5, the results are presented and discussed in Sec. 6, and a summary of the results and conclusions are given in Sec. 7.

The ORs are a critical microstructural characteristic that ensures different crystalline interphases are accurately represented [20]. The ORs for the f.c.c. and b.c.c. hydrides are based on the formulation developed by Ref. [21], and in this investigation, we used an OR of (0001)_h.c.p.//(111)_f.c.c. for the f.c.c. δ hydrides and an OR of (0001)_h.c.p.//(110)_f.c.c. for the b.c.c. ɛ hydrides. There are numerous experimentally observed ORs, but as discussed in Ref. [21], these are the most commonly observed experimental ORs for hydrided zirconium alloys.

By splitting an element into two halves along the crack line as first proposed by Hansbo and Hansbo [25] and detailed by Wu and Zikry [26], cracks in the mesh can be obtained. The advantage of the overlap method is that it overcomes the limitations of methods, such as cohesive fracture and element enrichment. As noted by Wu and Zikry [26], the cohesive fracture can introduce compliant surfaces that can affect numerical stability. Enrichment methods are limited by elastic fracture solutions that may not be applicable for inelastic failure modes.

An approach that has been developed specifically for microstructural fracture and dislocation-density crystalline plasticity will be utilized in this investigation. Only a brief outline will be given here, the methodologies are detailed in Refs. [13,28]. An finite element (FE) method with Broyden–Fletcher–Goldfarb–Shanno (BFGS) iteration and one-point integration with hour-glass control is used to obtain the total deformation rate tensor. The dislocation-density evolution equations and the plastic deformation rate tensor is updated by an implicit–explicit fifth-order Runge–Kutta method that accounts for numerical stiffness that arises due to the field quantities associated with the different crystalline structures evolving at different rates. A Crank-Nicholson approach was used to update the local temperatures with both thermal conductivity and internal heat generation due to dissipated plastic.

The dislocation-density crystalline representation, the fracture nucleation framework, and the Voronoi scheme for hydride formation with OR interfaces were used to investigate the behavior of hydrided zirconium systems. The multiphase aggregate was assumed to have random low-angle misorientations, where the maximum misorientations were assumed to be less than 5 deg. Voronoi tessellations were used with different Euler angles to represent these misorientations, and plane-strain boundary conditions were used to represent cladding behavior.

Seven thousand plane-strain four-noded quadrilateral elements (Fig. 2) were used for an aggregate size of 200 μm by 200 μm, and the mesh was optimized based on a convergence analysis for different grain misorientations. For tensile plane-strain behavior, a quasi-static nominal strain rate of 200 × 10⁻¹ s⁻¹ was applied. The material properties that are used here were obtained for experimental measurements of Zircaloy-4, and δ and ɛ hydrides [6,7], and are listed in Table 1. The main differences between the hydrided aggregates are the b.c.c. and f.c.c. hydrides, which result in distinct fracture behavior due to the evolution of the different dislocation densities and slip-systems.

To fully understand and predict how δ and ɛ hydrides affect fracture nucleation and propagation, we investigated aggregates, for each crystalline hydride, with radial, circumferential, and mixed distributions with different HVF. We refer to distributions as mixed as those with high HVFs that result in aggregates that are then neither fully circumferential or radial. Furthermore, to understand the effects of temperature on deformation and local behavior, we investigated aggregate behavior with radial hydrides for three different temperatures of 20 °C, 250 °C, and 500 °C.

The microstructural fracture behavior of hydride zircaloy-4 for radial and circumferential hydride distributions, with different hydride crystallographic structures and different hydride volume fractions (HVFs), were investigated. The failure criterion is based on the normal components of the stress traction exceeding the experimentally obtained fracture stresses given in Table 1. The fracture formulation also accounts for the crystallographic structures of ORs that are uniquely inherent to f.c.c. and b.c.c. hydrides and the parent h.c.p. matrix microstructures [21].

We investigated three different fracture stresses to further understand the response of the Zircaloy-4 system without hydrides (Fig. 3), and these cases correspond to different fracture stresses of 720, 840, and 960 MPa, and it confirms the sensitivity of the matrix to different fracture stresses. Namely that as the fracture stress increases, the fracture is inhibited as a function of the applied strain, which is physically realistic. Based on this, we investigated the following cases to obtain a more detailed representative understanding of fracture:

• Different morphologies with different hydride volume fractions.

• Radial and circumferential distributions.

• Hydrides with different crystalline structures: δ f.c.c. and ɛ b.c.c. hydrides.

To better understand the mechanistic effects of HVF on the fracture behavior, we analyzed three different models with different HVFs for a fracture stress of 840 MPa and at a room temperature of 20 °C for a radial distribution of δ f.c.c. hydrides. As shown in the stress–strain response (Fig. 4), as the HVF value increases, the fracture is inhibited. For the higher HVFs, there are significant increases in plasticity, as indicated by global behavior beyond the apparent yield point, before the global unloading of the stresses. This can be further understood from the contours for the different values of the HVF, which are given for the nominal strains that correspond to the final point of unloading. The contours are for the normalized (normalized by the static yield stress) fracture stress, and the circled regions correspond to the different crack nucleation sites. All the contours correspond to nominal strains of 1.8%. As these results show, the aggregate with the higher HVF ruptures at a much higher nominal strain than the two other cases (Fig. 5) and much higher fracture stress is needed in comparison with the two other cases. This can be attributed to the significant plasticity associated with the higher nominal strains of the higher HVFs. Plasticity ahead of the crack fronts inhibits crack propagation, which requires higher fracture stresses to rupture the aggregate. This is physically consistent with fracture observations of δ radial hydrides within zirconium alloys [16].

We wanted to further understand the effects of hydride orientation on this local ductile fracture behavior, so we investigated the behavior of δ-hydrides in the radial and circumferential directions with the highest fracture stress of 960 MPa. Figure 6 shows the global nominal stress–strain response for the two models at an HVF of 7%. For the radial hydrides, there is rapid unloading, in comparison with the circumferential case. For the circumferential hydride case, there is more ductility, which would delay crack rupture as the predicted results indicate. Billone et al. [8] reported that radial hydrides promote unstable crack propagation, whereas circumferential hydrides tend to promote plasticity and/or stable crack growth, so the predicted results are consistent with these experimental results. The results also showed that the hydrides act as initiation sites for crack nucleation (Fig. 7) for both the radial and circumferential hydrides. There are multiple crack nucleation sites, but not all of them coalesce to form bigger cracks. All the contours correspond to nominal strains of 2.2%. As a single crack develops, others could be blocked by potential pile-ups and accumulations of dislocation densities (Fig. 8).

We investigated two radial hydride crystallographies: δ-hydrides, which are f.c.c., and ɛ hydrides, which are b.c.c., as outlined in Table 1. We examined the effects of these two types of hydrides, which span different fracture stresses (Table 1) for an HVF of 17%. At this high HVF, the contours indicate that the hydride distribution can be considered to be mixed: neither circumferential nor radial due to the high volume fraction. Figure 9 shows the nominal stress–strain response for the two hydrides types. As these results indicate, the aggregate with the δ f.c.c. hydrides unloads gradually in comparison with the ɛ b.c.c. aggregate. This gradual unloading is due to the higher plasticity associated with the δ f.c.c. aggregate in comparison with the b.c.c. aggregate, and this can be due to the different ORs associated with the h.c.p./b.c.c. and the h.c.p/f.c.c. systems as discussed by Ref. [29]. Figure 10 shows the normalized fracture stresses for the two hydride populations. As these predictions indicate, and based on the global nominal stress strain, the gradual unloading of the aggregates with the delta hydrides results in rupture at a higher nominal strain and a higher normalized fracture stress (Fig. 10). The contours correspond to a nominal strain of 2.8%. The higher plasticity associated with the δ hydrides results in the higher fracture stresses. This is consistent with experimental observations of how the different microstructures can affect fracture [30].

We next investigated the effects of the thermomechanical loading on the behavior of Zircaloy-4 deformation and cracking. The coupling between the heat conduction and the mechanical loading was done at each time-step. Here, we will focus on radial hydrides and present the effects of thermomechanical loading on ductility, crack nucleation, and propagation.

To understand the effects of temperature on deformation and local behavior, we subjected the aggregate with radial hydrides to three different temperatures of 20, 250, and 500 °C for an HVF of 10%, and we inhibited fracture to focus on ductility and how deformation modes evolve with increasing temperature. Figure 11 shows the global nominal stress–strain response for the three temperatures. As the temperature increases, the response is significantly more ductile. Furthermore, as the applied temperatures increase, the normalized (normalized by the static yield stress) decreases as shown in Fig. 12. This has also been observed by Kim et al. [15,31]. This ductility is due to the increased normalized immobile dislocation density as the temperature increases (Fig. 13). These dislocation-density systems are the most active normalized dislocation density for the pyramidal system, as shown.

To further understand the effects of thermal softening due to high temperature on the dislocation density, a section plot is shown in Fig. 14 for a temperature of 500 °C. This line was taken along a hydride and shows the evolution of the immobile dislocation density along the x-direction. The plot shows that the dislocation-density values for both the hydrides and the parent matrix at a nominal strain of 12%. The dislocation densities for both f.c.c. hydrides and the h.c.p. matrix are almost of the same magnitude at this high temperature, which further underscores the thermal softening associated with this high temperature.

To investigate the effect of temperature on fracture, we investigated two aggregates with two different temperatures of 20 °C and 500 °C. Figure 15 shows the global nominal stress–strain response for the two aggregates. As the temperature increases, the aggregates become more ductile. As noted earlier, this was also observed by Kim et al. [15,31] for δ hydrides. At room temperature, the aggregate has little ductility in comparison with the aggregate subjected to higher temperatures. Also, at room temperature, the specimens were fractured at the early stages of deformation and breached even before reaching the yield point. Similar behavior is also obtained for the circumferential hydrides (Fig. 16) at an HVF of 10%.

The effects of temperature on shear strain accumulation can be clearly seen in Figs. 17 and 18 for both radial and circumferential distributions. As the contours show, at room temperature, for both radial and circumferential distributions, rupture is significantly inhibited as the temperature increases. As the predictions indicate, at a nominal strain of 2.5% for both thermal cases, the shear strain for both radial and circumferential distributions, the shear strain at 500 °C is almost four times larger than what it is at 20 °C. This increased dislocation-density activity and shear strain accumulation are the mechanistic effects that inhibit fracture at high temperatures for these hydride populations and orientations. These high temperatures lead to thermal softening, shear localization, and ductile failure [2,32].

A fracture methodology specialized for simultaneous crack surface nucleation was implemented within a nonlinear finite-element framework to accurately represent a microstructural multiphase system with f.c.c. δ and b.c.c. ɛ hydrides within an h.c.p. Zircaloy-4 matrix. This validated predictive framework provides an understanding of fundamental microstructural mechanisms and characteristics that span the physical scales of crystalline structures of the hydrides, the parent matrix, and the ORs between the h.c.p. matrix and the f.c.c. and b.c.c. hydrides.

Radial, circumferential, and mixed hydride orientations were investigated to understand how inelastic mechanistic behavior of these multiphase systems evolves for radial, circumferential, and mixed orientations for different HVFs. Circumferential hydride distributions had higher strength and ductility than the radial distributions. A key aspect for both distributions is that the hydrides do contribute to ductile behavior due to dislocation-density accumulations at hydride-matrix interphases. As the HVF was increased beyond a volume fraction of 7%, the hydride distributions were mixed, in the sense that radial and circumferential hydrides could no longer be distinguished. As the HVF increased, the stresses and dislocation densities increased, which further indicates that hydrides both can strengthen the hydrided system and render it more ductile.

The hydrided ɛ b.c.c. system generally had more strength and ductility than the δ f.c.c. system, and this was mainly because of the availability of more slip-systems associated with the b.c.c. ɛ hydrides. This has resulted in the matrix with the b.c.c. ɛ hydrides having fracture modes that occurred at higher nominal strains than the matrix with the f.c.c. δ hydrides. A key aspect of the predictions was the incorporation of temperature changes that accounted for both thermal conductivity and the dissipation of plastic work. For hydrided aggregates subjected to high temperatures, thermal softening resulted in higher ductility due to increased dislocation-density activity, which led to higher shear strain accumulation and inhibited crack nucleation and propagation. The present analysis indicates how underlying thermomechanical mechanisms collectively related to slip-systems, dislocation densities, local stresses, volume fractions, phase boundaries, and interfaces between crystalline systems affect and dominate local microstructural behavior. Therefore, it is essential to account for and identify these dominant microstructural mechanisms to understand how multiple fracture surfaces nucleate and propagate in these multiphase crystalline systems.

Supports from the Consortium for Advanced Simulation of Light (CASL) Water Reactors an Energy Innovation Hub for Modeling and Simulation of Nuclear Reactors, U.S. Department of Energy, Contract No. DE-AC05-00OR227 and the DOE NEUP Integrated Research Project IRP-DOE 17-13708, Development of a Mechanistic Hydride Behavior Model for Spent Fuel Cladding Storage and Transportation are gratefully acknowledged.

There are no conflicts of interest.