Phase-Field Model of Hydride Blister Growth Kinetics on Zirconium Surface

Introduction

Zirconium alloy has a low neutron absorption rate and good high-temperature corrosion resistance and is widely used in structural materials in fuel cladding of water-cooled reactors (Krishnan and Asundi, 1981; Bair et al., 2015; Motta et al., 2015). However, severe waterside corrosion reactions occur during use and storage, and its performance will also be limited (Hong et al., 2002; Han et al., 2019). During the corrosion process on the waterside, the outer layer of the zirconium alloy cladding is oxidized, releasing hydrogen atoms, and the hydrogen atoms migrate into the zirconium matrix. It moves along the hydrostatic stress gradient and moves against the thermal gradient (Kammenzind and Berquist, 2000; Reheman and Ståhle, 2018). In other words, hydrogen is attracted to high-stress areas and cold areas, or areas where the two interact together. When the hydrogen concentration exceeds the solid solution limit in the zirconium matrix, zirconium and solid solute hydrogen transform into the zirconium hydride phase (Ghosal et al., 2002; Bair et al., 2015; Reheman and Ståhle, 2018; Han et al., 2019). The formation of brittle hydride will not only cause the material to expand and increase its volume greatly, but also it is easy to initiate crack growth and penetrate the matrix (Kim et al., 2007; Puls, 2009; McRae et al., 2010). The hydride formed on the metal surface appears as a small bump, usually called a blister (Cheong et al., 2001; Singh et al., 2002). The growth of zirconium hydride blister is an important form of pitting corrosion. It is the main initiation point of fuel cladding failure in nuclear power plants, which will pose a serious threat to nuclear power safety. Therefore, the basic understanding of hydride blister growth behavior is of great significance for predicting material performance degradation and failure.

The size of the hydride blister is much larger than that of the single hydride, reaching the micrometer level (Cheong et al., 2001; Singh et al., 2002; Long et al., 2017; Kim et al., 2018). The formation of the hydride blister requires the existence of a cold spot in the surrounding environment, and the effect of cold spot aggregation on the growth of hydride blisters is also different (Singh et al., 2002). In the experiment of (Domizzi et al., 1996), the growth of blisters on the surface of the Zr-2.5% Nb alloy tube was studied. It divides the observed phenomenon into three areas (I) the main part of the zirconium hydride blister, whose growth is similar to that of a single hydride, and the anisotropic growth is ellipsoidal, (II) the radial hydride accumulation area, and (III) circumferential hydride accumulation area. The phase-field method has been widely used to study the microstructure changes such as hydride precipitation in Zr cladding materials (Bair et al., 2015; Bair et al., 2016; Bair and Zaeem, 2017; Bair et al., 2017; Han et al., 2019; Toghraee et al., 2021). Based on region (I) (Reheman and Ståhle, 2018), studies the growth and cracking of hydride blisters under elastic-plastic action by using the phase-field method, but its model is established based on isotropy, which can not represent the anisotropic growth of zirconium hydride blisters. In this work, we study the morphological evolution of zirconium hydride blisters by using the phase-field method coupled with anisotropic elastic strain energy. In addition, some studies (Cheong et al., 2001; Motta and Chen, 2012) show that during reactor operation, the fission gas and the helium filling gas in the cladding tube create a combined pressure that exposes the nuclear fuel cladding to constant stress. Even the swelling of nuclear fuel also produces a certain pressure on the zirconium alloy cladding tube (Rozhnov et al., 2011). Therefore, it is important to understand the evolution mechanism of zirconium hydride blisters under specific pressure. The evolution of hydride blisters under radial stress is considered in this work.

Phase-Field Model

Combining the zirconium hydrogen phase diagram (Yu et al., 2020) and experimental research (Santisteban and Domizzi, 2009; Vicente Alvarez et al., 2011; Lin et al., 2016; Motta et al., 2019), it can be found that the hydrides in the zirconium hydride blister are all δ-hydrides, and the ratio varies approximately in the range of 60–85%, and the rest is the Zr matrix and impurities. Since the chemical energy parameters of the zirconium hydride blister are not clear yet, we simulate the zirconium hydride blister by using the chemical energy (Bair et al., 2016; Bair and Zaeem, 2017; Bair et al., 2017; Toghraee et al., 2021) and elastic parameters (Bair et al., 2017; Han et al., 2019) of the δ-hydride. The size of the cladding tube in the model is similar to some studies (Hong et al., 2002; Hellouin de Menibus et al., 2014), and 11.25-degree section is simulated. Since the longitudinal plane of the blister is usually observed in the experiment, we convert the data in the z-y direction of the coordinate axis to x-y. The whole calculation process is carried out in two dimensions. Figure 1 shows the setting of the hydride blister growth model and boundary conditions. For the displacement field $u$, the top boundary is set as free boundary conditions, the left and right are periodic boundary conditions, and the bottom is fixed boundary conditions. The boundary settings of phase-field η and concentration field c are the same. All boundaries are set as zero flux boundary conditions, and hydride and the matrix phase are distinguished by the diffuse interface.

FIGURE 1. Schematic diagram of zirconium-zirconium hydride blister model and boundary condition settings.

Governing Equations and Bulk Free Energy

Interfacial reaction, hydrogen transport driven by diffusion potential, and internal stress caused by precipitation and growth of zirconium hydride blister. Multi-physical field coupling based on dynamic equations can simulate these effects and reflect the dynamic characteristics (Shi et al., 2021; Liu et al., 2022; Shi et al., 2022). The phase field, hydrogen concentration field, and elastic strain energy field were coupled to establish the zirconium hydride blister growth model in this study. These fields are necessary to track the movement of the diffusion interface, calculate the diffusion/redistribution of hydrogen (Yang et al., 2021), and reveal the effect of radial stress on hydride growth. The advantage of the phase-field model in this study is that the shape of the grown precipitates can be calculated (Gao et al., 2021; Gao et al., 2022; Xu et al., 2022). The phase transition from metal to hydride is given by a phase variable η that continuously covers a range of values, where the ends of the range represent pure metal or pure hydride, respectively. The analysis is based on totally free energy, including strain energy, chemical potential energy, and gradient energy describing the diffusion interface. In the phase-field model, both the diffusion of hydrogen atoms and the precipitation of hydrides are considered. The time evolution of the concentration c is determined by the Cahn-Hilliard diffusion equation (Kim et al., 2008), and the evolution of the microstructural parameter η is governed by the time-dependent Allen-Cahn equation (Allen and Cahn, 1979). In order to solve the evolution equation, the finite element method is adopted, and the mesh adaptive method is used to speed up the calculation.

The driving force of the zirconium-hydride blister phase-field model comes from the minimization of the total free energy ($F_{\text{total}}$) of the system, which includes the uniform bulk chemical free energy ($F_{\text{bulk}}$) the interface gradient energy ($F_{\text{int}}$) and the elastic energy ($F_{\text{el}}$). For simplicity, it is assumed that the diffusion coefficients of hydrogen in the matrix and precipitate are the same.

Here, all the parameters used in this simulation are listed in Table. 1 The bulk chemical free energy $F_{\text{bulk}}$ is a function of the concentration $c$ and the order parameter $\eta$ at T = 600 k, while $F_{\text{int}}$ is the gradient energy due to the inhomogeneity of the interface and $F_{\text{el}}$ is the elastic strain energy during the Zr matrix phase to the hydride phase transition.

TABLE 1. Model parameters and conditions used in the simulations.

The Kim-Kim-Suzuki (KKS) binary alloy model (Kim et al., 1999) is used to determine the expression form of bulk chemical free energy density. The KKS model assumes that each material point is considered to be a mixture of two or more phases with different concentrations, and the diffusion potentials of the same concentration between different phases are equal. The assumption can be expressed as follows:

where $h\left(\eta \right)=3{\eta }^2-2{\eta }^3$ is the interpolation function, the double potential well function $g\left(\eta \right)=\frac{1}{2}\left({\eta }^2\times {\left(1-\eta \right)}^2\right)$ , and $w$ is the potential barrier. The free energy density in this study is constructed according to the quadratic fitting method adopted by (Bair et al., 2016; Bair and Zaeem, 2017; Bair et al., 2017; Toghraee et al., 2021). Taking advantage of the equal slope at the two-phase equilibrium point, the two-phase free energy curve can be constructed by using the artificial parameters $A_{\alpha }$ and $c_{\alpha }$.

Two-phase molar free energy $G_{\alpha }$ and $G_{\delta }$ at T = 600 K can be expressed as:

Here, $c_{\alpha }$ and $c_{\delta }$ are the molar fractions of atomic H in the matrix and blister. $c_{\alpha }^{\prime }$ is the concentration of maximum hydrogen solubility in the matrix before precipitation occurs (Une et al., 2009); $c_{\delta }^{\prime }$ is the hydrogen concentration in the hydride blister; $\Delta G_{\delta }^0$ is the molar Gibbs free energy of hydride blister formation. The constant $A_{\alpha }$ and $A_{\delta }$ are the contribution of the free energy of the two phases to the interface, which are determined by constructing a common tangent between the molar Gibbs free energy curves of the two phases (Yang et al., 2020; Yang et al., 2021; Sheng et al., 2022), from which an approximate expression of the free energy can be obtained. $f_{\alpha }$ and $f_{\delta }$ are the free energy densities of the zirconium matrix and hydride blister, respectively, as shown in Figure 2 The blue dashed line is the common tangent when the two phases are in equilibrium.

FIGURE 2. Two-phase chemical free energy density curves of zirconium matrix and zirconium hydride blister.

The interfacial energy density can generally be written as a function of the gradient of the phase-field variable,

where ${\kappa }_{\eta }$ is the gradient energy coefficient related to the phase-field variable $\eta$.

Elastic Strain Energy and Applied Stress

When hydride blister precipitate and grow in the matrix, their elastic strain energy comes from the structural difference between hydride blister and matrix. The elastic strain energy $F_{e\text{l}}$ is expressed as follows:

where $f_{\text{el}}$ is the mismatch elastic strain energy density, ${\sigma }_{ij}$ and ${\mathit{ϵ}}_{ij}^{el}$ are the stress and elastic strain, respectively. Considering that the cladding tube will be subjected to the radial outward stress from the swelling of nuclear fuel and so on, we add a new cylindrical coordinate system and add the radial outward stress ${\sigma }_{app}$ to the model to study the evolution difference of the blister under different stresses. The radial stress ${\sigma }_{app}\mathit{ }$ is added through the stress coordinate conversion in elasticity.

${\sigma }_r$ and ${\sigma }_{\phi }$ are the radial stress and circumferential stress applied in the cylindrical coordinate system respectively. $\phi$ is the angle value in the cylindrical coordinate system. ${\sigma }_x$, ${\sigma }_y$, ${\tau }_{xy}$ are the stress value in the rectangular coordinate system.

The elastic strain ${\epsilon }_{ij}^{el}$ can be obtained from Equation 19, where ${\epsilon }_{ij}\left(r\right)$ is the total strain represented by the components $u_i$ and $u_j$ of the displacement field variable $u$ obtained by solving the mechanical equilibrium Equation 23 and ${\epsilon }_{ij}^{00}\left(r\right)$ is related to the lattice inherent strain ${\epsilon }_{ij}^{00}$. Tummala (Tummala et al., 2018) and Singh (Singh et al., 2007) give the eigenstrain of hydride at different temperatures. The eigenstrain of the hydride blister at T = 600 k can be obtained.

where $C_{ijkl}\left(\eta \right)$ is the elastic constant, a function of the order parameter $\eta$, $C_{ijkl}^{\alpha }$ and $C_{ijkl}^{\delta }$ are the elastic constants of the matrix and precipitate, respectively, and $\Delta C_{ijkl}$ is the difference between the elastic constants of $C_{ijkl}^{\alpha }$ and $C_{ijkl}^{\delta }$.

Results and Discussion

The effect of anisotropic elastic energy on the growth morphology of hydride blisters is explored, as shown in Figures 3A,C, which are the evolution results of hydride blisters after initial state t = 0 s and final state t = 1500 s, respectively. Due to the effect of elastic energy, the corrosion width of hydride blister is greater than the corrosion depth, the growth characteristics of anisotropy appear, and the morphology is semi ellipsoid. In addition, there is a slight bulge at the top of the blister. The growth of the hydride blister is closely related to the intrinsic strain, and to minimize the elastic properties, different growth rates will appear in different directions. Figures 3B,D show that the concentration in the hydride blister is relatively high, and there is a diffusion area at the interface between the matrix and the hydride blister. Figures 3E,F are the stress distribution results in the radial and circumferential directions of the coupled anisotropic elastic energy, respectively. The hydride blisters in both directions are in a compressed state, and there are stress concentration area near the phase interface and the tips on both sides. The stress concentration will redirect the δ-hydrides around the hydride blisters, form more harmful radial hydride clusters, and even cause the initiation and propagation of microcracks.

FIGURE 3. The evolution results of the coupled anisotropic elastic energy zirconium hydride blister, where (A) and (C) are the morphology evolution results (order parameter field) of the hydride blister growth at t = 0 s and t = 1,500 s, (B), and (D) is the distribution of hydrogen concentration field after t = 0 s and t = 1,500 s, (E) and (F) are the stress distribution of ${\sigma }_{rd}$ (Radial stress) and ${\sigma }_{\phi }$ (Circumferential stress).

Zirconium alloy cladding tube will be subjected to a certain radial outward stress during operation, which will have a great impact on the formation and growth of hydride. We applied stress of 1, 2, and 3 GPa along the radial direction in this model. Because this model is established under the ideal condition of no defects, the applied stress is much greater than that of the actual cladding tube. We can qualitatively study the effect of radial stress on the growth morphology of hydride blisters. The evolution results after t = 1,500 s are shown in Figure 4. Figures 4A–C shows that the width and depth of hydride blister increase with the increase of stress under the application of different radial stress. Figure 4D shows the contour comparison diagram of hydride blister under t = 1,500 s applied stress and no applied stress. The results in Figure 4 confirm that applying radial stress will change the morphology and dynamic evolution of blisters, promote the growth of blisters and enhance the process of pitting corrosion.

FIGURE 4. At t = 1,500 s, (A–C) are the growth morphologies of hydride blisters under 1, 2, 3 GPa radial stress, respectively (The direction of stress application is indicated by blue arrows). (D) is the contours of hydride blisters.

Figure 5A shows the change of corrosion depth of the hydride blister with time. It can be seen intuitively that the growth of hydrogenated blister is affected by different radial stresses, and the corrosion depth of hydride blister will increase with the increase of time and radial outward stress. Figure 5B shows the changes of the area of hydride blister and the average corrosion rate along the pitting depth under different radial stresses. It can be found that the greater the applied radial stress is, the more the corrosion area will be increased, which also shows that the radial stress will increase the influence of pitting corrosion. The corrosion rate changed from fast to slow, which may be due to the continuous consumption of hydrogen in the surrounding environment during the evolution process. The reduction of hydrogen content has become one of the key factors restricting its growth. After the stress is applied, the initial growth speed becomes larger, and the larger the value is, the more obvious the speed increase is, and then it will tend to be consistent. The matrix has a compressive effect on the expansion and growth of hydride blisters. The expansion state of the hydride blister in the simulation area is realized by the bending of the grid (Sheng et al., 2022). After the external stress is applied, the external stress will affect the overall stress state of the blister, and then affect its growth morphology. The growth of hydride bubbles is controlled by these three effects. Figure 5C,D show the distribution of displacement $u_x$, $u_y$ and η along the red arrow line after applying different radial stresses at t = 1,500 s. On both sides of the hydride blister, $u_x$ has the same size and changes in the opposite direction. The application of stress will change the initial distribution of displacement, but will not affect its change trend. Combined with the change of order parameters, it can be found that the displacement changes rapidly at the matrix-blister interface and reaches the maximum at the tip of the blister. The application of stress will promote the movement of the interface, thus changing the change of displacement at the same grid point, showing a larger expansion trend on both sides. In the depth direction of the hydride blister, when no stress is applied, as shown in the black line in the figure, the displacement $u_y$ first starts from zero (Due to the displacement of the lower surface $u_y$ is constrained) and then the value starts to increase near the phase interface, which means that this part of the matrix point moves along the z-axis and reaches a maximum at the upper free surface. The protrusion of hydride blister can also be realized by this part. This part can realize the bulge of hydride blister, which is the reason for the bulge in Figure 3. When the radial stress is applied, the stress distribution of the hydride blister and the matrix changes, and the displacement will change from the fixed bottom along the direction of the red arrow. The positive displacement changes after the stress are applied. A general trend is a form of “check”, and the inflection points appear at the interface.

FIGURE 5. (A) The corrosion pit depth, (B) the average corrosion rate and the area of hydride blister as a function of time under different radial stresses. Variation of (C) displacement component $u_x$ (solid line) and order parameter η (dotted line), (D) displacement component $u_y$ (solid line) and order parameter η (dotted line) along the red arrow in the diagram under different radial stresses at t = 1,500 s (The displacement of the green dotted line and the order parameter of the light blue dotted line are zero).

As mentioned above, applying different of radial outward stress to the whole simulation area will have a great impact on the growth of the hydride blister. Figure 6A and Figure 6B are t = 1,500 s respectively, under the application of 2 GPa radial stress ${\sigma }_{rd}$ and ${\sigma }_{\phi }$ stress distribution. Figure 6C and Figure 6D are t = 1,500 s respectively under different stress conditions ${\sigma }_{rd}$, ${\sigma }_{\phi }$ and η with the change of the black arrow in the diagram. Combined with the evolution process of order parameters, it can be found that when hydride blisters grow normally, the matrix bears compressive stress. The stress value changes rapidly at the interface and reaches the maximum, and there is a stress concentration area around the blister. The radial stress ${\sigma }_{rd}$ reverses at the interface, and part of the matrix bears tensile stress. The tensile effect caused by tensile stress is very weak, which can not prevent the expansion and compression of hydride. The circumferential stress ${\sigma }_{\phi }$ has a large stress concentration at the tip of the blister, which may lead to the cracking of the blister. It can be found from the change of the displacement $u_y$ in Figure 5D. When we apply stress and change the initial stress distribution, the overall stress of blister and matrix will change from compression to tension. No matter how the applied force changes, there will be stress concentration areas at the tip and interface of the hydride blister. The side of the upper surface of the blister close to the oxide layer will bear outward and upward stresses, accelerating the cracking of the upper surface. The growth process of zirconium hydride blisters is the result of the combined effect of chemical energy driving force and elastic energy driving force (Simon et al., 2021). The applied stress affects the growth of hydride blisters by affecting the elastic driving force, and the radial stress reduces the effect of the elastic driving force hindering the growth of hydrides. This leads to the difference in the morphology of hydride blisters under different stress states.

FIGURE 6. The distribution of (A) ${\sigma }_{rd}$ (Radial stress) and (B) ${\sigma }_{\phi }$ (Circumferential stress) under 2 GPa radial stress at t = 1,500 s. Variation of (C) ${\sigma }_{rd}$ and η, (D) ${\sigma }_{\phi }$ and η along the black arrow in the figure under different radial stresses at t = 1,500 s.

Conclusion

In this work, the anisotropic elastic properties are introduced into Zr-Hydride blister system, and the growth and evolution dynamics of zirconium hydride blister are studied by using phase-field model. The simulation results show that the natural growth morphology of hydride blister tends to be semi-ellipsoidal and in a compressed state, while the matrix near the interface is in a tensile state. There is a stress concentration area on the outer contour of the blister, and the maximum value appears at the end points on both sides of the hydride. The external stress affects the growth morphology of zirconium hydride blister by affecting the elastic stress distribution. Applying radial outward stress will promote the growth of blister, and its corrosion width and depth will increase. The greater the applied stress, the more obvious the promotion effect and the faster the growth rate. This model also has certain application value in other fields of pitting corrosion.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

HH and YL conceived and guided the work. SW conducted simulation and data integration. The results were analyzed and discussed by all authors. JS gave guidance in the theoretical part. SW, JS, HH and YL jointly wrote the manuscript, with contributions from all the authors.

Funding

This work was supported by the Science Challenge Project (Grant No. TZZT2019-D1-03) and the Foundation of LCP.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Authors would like to thank Changqing Guo, Yuechao Wang and Guomin Han for their help on the model and basic knowledge in this work.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References