Coupled Modeling and Simulation of Phase Transformation in Zircaloy-4 Fuel Cladding Under Loss-of-Coolant Accident Conditions

Under loss-of-coolant conditions, the temperature on fuel cladding will increase rapidly (up to 1000–1500 K), which will not only cause a dramatic oxidation reaction of Zircaloy-4 and an increase in hydrogen concentration but also cause an allotropic phase transformation of Zircaloy-4 from hexagonal (α-pahse) to cubic (β-phase) crystal structure. As we all know, thermophysical properties have a close relationship with the microstructure of the material. Moreover, because of an important influence of the phase transformation on the creep resistance and the ductility of the fuel rod, studying the crystallographic phase transformation kinetics is pivotal for evaluating properties for fuel rod completeness. We coupled the phase transformation model together with the existing physical models for reactor fuel, gap, cladding, and coolant, based on the finite element analysis and simulation software COMSOL Multiphysics. The critical parameter for this transformation is the evolution of the volume fraction of the favored phase described by a function of time and temperature. Hence, we choose two different volume fractions (0 and 10%) of BeO for UO2-BeO enhanced thermal conductivity nuclear fuel and zircaloy cladding as objects of this study. In order to simulate loss-of-coolant accident conditions, five relevant parameters are studied, including the gap size between fuel and cladding, the temperature at the extremities of the fuel element, the coefficient of heat transfer, the linear power rate, and the coolant temperature, to see their influence on the behavior of phase transformation under non-isothermal conditions. The results show that the addition of 10vol%BeO in the UO2 fuel decreased the phase transformation effect a lot, and no significant phase transformation was observed in Zircaloy-4 cladding with UO2-BeO enhanced thermal conductivity nuclear fuel during existing loss-of-coolant accident conditions.


INTRODUCTION
With the increasing energy demand, human beings have never stopped exploring new energy sources. Among them, nuclear energy is one of the most popular and promising future energy sources. Although nuclear energy is considered a clean and efficient energy source, it deals with a serious hazard -nuclear radiation. Many studies have shown that complex lesions can easily occur in cellular DNA (Sutherland et al., 2000;Sutherland et al., 2001;Yang et al., 2004), which means that nuclear radiation can cause massive damage to the human body, causing people to suffer from cancer and even death.
Nearly 70°years have passed, nuclear energy has been developed from the 1950s to the present. In these years, nuclear energy's development process is not smooth sailing, which means that there have been accidents in nuclear energy development history. Some major nuclear power cases have attracted the whole world's attention on nuclear safety, such as the Chernobyl nuclear leak accident in 1986, rated as INES 7 and considered the worst nuclear power accident in history. More recently, the Fukushima Daiichi Nuclear Power Plant accident in 2011, which is also rated as INES 7, brought widespread attention. These accidents caused substantial economic losses and caused many people to be exposed to varying degrees of radiation pollution, affecting health and even losing their lives. What is more, huge adverse effects have been caused on nuclear energy development worldwide. Thus, safety issues occupy a vital position in the development of nuclear energy.
With respect that an entire reactor of a nuclear power station is incredibly complex, there are many ways to ensure nuclear power plants' safety. One of the critical aspects is how to avert the leakage of radioactive energy as much as possible in a loss-ofcoolant accident (LOCA). In LOCA, the pivotal step to controlling nuclear leakage is to maintain the nuclear fuel rod's integrity and the fuel cladding when they are inserted into cold water under the emergency cooling system. If the fuel cladding is not firm and tenacious enough, the nuclear fuel rod may crack easily, causing the radioactive material to break through the cladding and leak, resulting in severe consequences. Subsequently, accident tolerant fuel and cladding materials under LOCA conditions are studied by many researchers (Isobe and Suda, 1999;Forgeron et al., 2000a;Sawarn et al., 2014;Park et al., 2016;Suman et al., 2016;Gamble et al., 2017;Tang et al., 2017;Jailin et al., 2020), especially after the event at Fukushima Daiichi Nuclear Power Plant. Accident Tolerant Fuel, or ATF for short, is a new generation of the fuel system to improve the fuel element's ability to fight against severe accidents. Compared to the previous fuel system, this updated fuel system can resist grave accident conditions for a longer time and, in the meanwhile, maintain the same or have even better performance under normal operating conditions. The unique material characteristics of accident tolerant fuel can slow down the velocity of deterioration in severe cases, which would help people reclaim more valuable time to take emergency measures so that the radiation leakage risk of fuel is significantly reduced. To sum up, the critical point of accident tolerant fuel is the endurance capacity in a loss-of-coolant accident. For the moment, there are two effective methods to enhance this capacity, one is to find new material (such as SiC, FeCrAl) to replace the fuel cladding that we are using now, and the other is to enhance the properties of the existing fuel system (such as coating cladding, modifying fuel). For the former, the database of new cladding material on the physical properties and phenomena is scarce. There are still many properties and phenomena of existing cladding that are not well understood at present. Further research on the existing cladding material is beneficial to analyzing accident tolerant fuel and cladding as the object for comparison and validation (to keep its strong points and overcome the shortcomings). Our work has used UO 2 -BeO enhanced thermal conductivity nuclear fuel and zircaloy cladding as an accident tolerant fuel system, which means we have chosen the latter method to research the physical properties of cladding.
In today's nuclear energy industry, zirconium alloys (especially Zircaloy-4) are still widely used as structural materials for reactors due to their superior properties. The small thermal neutron capture cross-section of zirconium allows it to ensure sufficient thermal neutrons to sustain the reactor's normal operation. Also, zirconium alloy has the advantages of strong corrosion resistance (Isobe and Suda, 1999) and excellent mechanical properties, making the research of zirconium alloys internationally occupy an increasingly important position.
Under extreme conditions, such as in a loss-of-coolant accident, the fuel cladding will undergo a rapid temperature increase (up to 1000-1500 K) (Hales et al., 2016), which will not only cause a dramatic oxidation reaction of Zircaloy-4 and an increase in hydrogen concentration but also cause an allotropic phase transformation of Zircaloy-4 from hexagonal (α-pahse) to cubic (β-phase) crystal structure (Northwood and Lim, 1979). Thermophysical properties are closely related to the microstructure of the materials themselves. In other words, under LOCA conditions, the behavior of fuel rod cladding depends mainly on the metallurgical evolution at high temperatures. Researchers have also pointed out an important influence of the phase transformation on the creep resistance and ductility, two essential characteristics for fuel rod integrity (Forgeron et al., 2000a). Therefore, studying the crystallographic phase transformation kinetics is pivotal for evaluating the mechanical properties essential for fuel rod completeness (deformation and burst) to improve its performance during LOCA conditions. The essential parameter for the transformation kinetics is the evolution of the new phase's volume fraction as a function of time and temperature. This paper has selected the method for calculating the volume fraction of the advantageous phase in Zircaloy-4 as a function of time and temperature during phase transformation under non-isothermal conditions (Hales et al., 2016).
This study has implemented the physical model of phase transformation coupled with the existing physical models for reactor fuel, gap, cladding, and coolant, based on the finite element analysis and simulation software COMSOL Multiphysics. COMSOL Multiphysics originated from the Toolbox of MATLAB, is a finite element analysis and simulation software that is good at coupling multiple physical fields described by the PDEs. Besides, some physical models of the UO 2 -BeO-Zircaloy fuel-cladding system have already been implemented (Liu et al., 2015) into the COMSOL Multiphysics finite-element platform.

Model Geometry
The model used in this work adopted a 2D axisymmetric plane with UO 2 -BeO fuel rod and Zircaloy-4 cladding (see Figure 1A). For the reason of the periodic boundary condition in the axial direction, a single pellet is chosen to represent all the pellets with a mapped mesh (see Figure 1B) (Liu et al., 2015).

Transition Model
We choose a variable y as the volume fraction of the new transformed phase (β-phase) as a function of time t and temperature T. The value of y is in the range of 0-1. Following the research of Leblond and Devaux (Leblond and Devaux, 1984), we considered that the value of y and the steadystate or equilibrium value y s (T) has not much difference at a given temperature T. In this case, we write dy dt where τ c is the characteristic time of phase transition. We can see that Eq.
(1) has two external functions y s (T) and τ c (T), depending on temperature, with y s (T) being the volume fraction of phase newly formed at the temperature T after an infinitely long time and y s (T) ∈ [0, 1]. Both of these functions are temperature-dependent quantities and can be derived from the experimental data. Under non-isothermal conditions, the temperature T may change as a function of time, and therefore Eq. (1) needs to be solved numerically. By putting the rate parameter k(T) ≡ τ −1 c (T), we rewrite Eq. (1) in a reduced form To calculate volume fraction y as a function of time and temperature, we need to specify the two function y s (T) and k(T), which appear in Eq. (2) Firstly, let us consider the former function y s (T). The steady-state volume fraction experimental data under phase transformation show that y s (T) has an S-shaped or sigmoid form (Massih, 2009). Thus, the equilibrium volume fraction of β-phase is represented by a sigmoid function of temperature.
where T cent and T span are material-specific coefficients connected to the center and span of the mixed-phase region, respectively. These two parameters are determined by where T α and T β are measured phase boundary temperatures corresponding to 99% α-phase and 99% β-phase fractions, respectively. For Zircaloy-4, we used T cent 1159 − 0.096w(K) and T span 44 + 0.026w(K) (Massih, 2009), where w is the hydrogen concentration in the range of 0-1,000 weight parts per million hydrogen (wppm). We have chosen 500wppm for the value of w because of the imperfection of the hydride formation model. In addition to the steady-state volume fraction y s (T), we also need to express the rate parameter k(T), which is the inverse of τ c (T), in detail.
From Eq. 5 , we can see that the rate parameter k(T) is highly dependent on temperature. Here, k 0 and k b are used to represent respectively the kinetic prefactor and the Boltzmann constant. E is used to show the total effective activation energy, and k m is just a constant changing with the direction of the phase transition. In this model, E is the effective activation energy, which associates the activation energy of nucleus formation and development. The justifiability of this associated effect has been discussed by Mittemeijer and his colleagues in (Mittemeijer, 1992;Kempen et al., 2002). For Zircaloy-4, we use k 0 60457 + 18129|Q|(S − 1 ) and E/k b 16650(K) (Massih, 2009;, where Q dT/dt(Ks −1 ) is used to represent the heat rate with |Q| ∈ [0.1, 100] Ks −1 (Massih and Jernkvist, 2009). The phase transformation α → β is entirely diffusion-controlled, while the transformation of the opposite direction is partly martensitic. This difference is revealed by the constant k m at the end of Eq. (5), which is given in the form   Threshold Temperature Models For the material-dependent temperatures for the beginning of phase transformation, the experimental data on Zircaloy-4 show that this quantity depends on the heating or cooling rate Q (Forgeron et al., 2000b). In this model, we have used the relation below to calculate the starting temperature in Kelvin for α → α + β transition based on the experimental data in (Forgeron et al., 2000b;Brachet et al., 2002).
where Q is the heating rate in Kelvin per second and the hydrogen concentration w is in the range 0 ≤ w ≤ 1000 wppm. Also, we have another relationship for the opposite direction of phase transformation basing on the experimental data of Zircaloy-4 reported in (Forgeron et al., 2000b;Brachet et al., 2002), i.e., from β → α + β.
A temperature (time) lag of the start temperature of phase transition from β to α + β is observed in the cooling process. Eq. (8) is not symmetric with the equation on heating, i.e., Eq. (7).
All these material-dependent quantities above allow us to calculate the β-phase volume fraction as a function of time by numerical integration of Eq. (2). Using the above-mentioned phase transformation models, the Zircaloy-4 cladding phase formation and redistribution can be investigated.

MODELING RESULTS
The behavior of fuel and cladding are presented for a 2D axisymmetric LWR fuel rodlet in COMSOL Multiphysics. The 2D axisymmetric model simulates a simplified fuel pellet with a typical finite element mesh shown in Figure 1B. The model contains an individual fuel pellet and Zircaloy-4 cladding. It is important to note a minimal but nonnegligible gap between pellet and cladding considering the actual situation. A width of 80°μm is considered to be the nominal gap size in this model. This section has used two different fuel systems, i.e., UO 2 -10% BeO and UO 2 , following the parameters setting in (Liu et al., 2015). The UO2-BeO fuel properties are shown in Supplementary Appendix S1. The typical RWR operating conditions used are shown in Table 1 above.
After implementing the phase transformation model for Zircaloy-4 cladding, we first calculated the volume fraction of β-phase under the normal operating condition to test this model, and we obtain Figure 2 as a result for the UO 2 fuel and Zircaloy-4 cladding. Figure 2 shows that the value of volume fraction of β-phase attains a very small magnitude. The green line with circles represents the volume fraction of β-phase of cladding outersurface, while the blue line with asterisks represents the cladding inner-surface. Although the value of the inner surface of cladding is almost twice that of the outer surface, it is barely more than 1.1 × 10 − 7 when it reaches its highest point, which means the effect of phase transformation in the Zircaloy-4 cladding is so small that it can be ignored under normal operating condition. The result calculated by this model is consistent with the figure shown in (Hales et al., 2016). Then we estimated the effect of phase transformation under postulated loss-of-coolant accident conditions.
During a LOCA condition, on account of the decrease in heat-transfer capacity under loss-of-coolant accident conditions, the heat or energy generated in the fuel cannot be passed outside, which causes the changes in the boundary conditions between fuel elements. Hence it has a high possibility that the temperature at the extremities of the fuel element will rise higher than before. The boundary conditions will change between fuel elements and change between fuel, cladding, and coolant because heat transfer capability decreases. For example, the heat transfer coefficient between cladding and coolant will decrease in the transition boiling interval. If we consider the swelling effect of fission gas, the fuel diameter will grow a little, which will make the size of the gap between fuel and cladding smaller. The linear power rate will also rise from 20000 W/m (set for normal condition) to 25,000 W/m, even 30,000 W/m when it is in serious condition. According to NRC (U.S. Nuclear Regulatory Commission) in (The U.S. Nuclear Regulatory Commission, 2011), a peak linear power density of around 60,000 W/m is considered the safe limit for core in operation. Hence, we can also try to raise the linear power rate to 60,000 W/m to observe its effect on phase transformation. The water saturation pressure at 530 K is between 3.3469 and 4.6923 MPa, according to Engineering Toolbox's data in (Engineering ToolBox (2004). However, from Figure 2, we can see the value of coolant pressure is up to 15.5 MPa while the coolant temperature is only 530 K. Per the data for water saturation pressure, the corresponding temperature of 15.5 MPa is between 613 and 633 K. In the following simulation, we have taken 613 K as the coolant temperature into account during the calculations.
To simulate a relatively severe loss-of-coolant accident, we had reduced the physical size of the gap between fuel and cladding, raising the linear power generation rate, the coolant temperature, and the temperature at the extremities of the fuel element. For the sake of decreasing the heat transfer coefficient of cladding to coolant, we have multiplied the coefficient with a number less than one. All these conditions are taken into account simultaneously, and the results obtained are shown below.
To make the results more clear, they are divided into four groups. The first group is the volume fraction of β-phase varying with the linear power rate, with 1500 K -the temperature at the extremities of the fuel element, 8e−7 m -the gap size, and 0,6 -the coefficient multiplied with the heat transfer coefficient between cladding and coolant. Figures 3A,C showed that the volume fraction on cladding inner-surface was increased almost by an order of magnitude from 4.3 × 10 − 5 to 2.6 × 10 − 4 when the linear power rate varied from 20,000 W/m to 25,000 W/m. However, when the power rate is increased by an additional 5,000 W/m to 30,000 W/m in Figure 3E, the volume fraction on the cladding inner-surface has reached 0.05, which is much bigger than before. Compared with UO 2 , the results for UO 2 -10vol%BeO in the left, i.e., Figures  3B,D,F, are all smaller. The increment of linear power rate up to 30,000 W/m for UO 2 -10vol%BeO did not have the same influence as that for UO 2 on the volume fraction of β-phase, which just attained around 4.0 × 10 − 4 on the cladding innersurface, the same magnitude for UO 2 with the linear power rate, which equals to 25,000 W/m.
The second group is about the effect of the temperature at the extremities of fuel element on the volume fraction of β-phase, with 30,000 W/m -the linear power rate, 8e−7 m-the gap size, and 0,6-the coefficient multiplied with the heat transfer coefficient between cladding and coolant. Figure 4 showed that there was not much difference in these results when the temperature at the extremities of the fuel element changed. From 700 K to 1100 K, the volume fractions on cladding inner-surface and outer-surface had all increased a little for both UO 2 and UO 2 -10vol%BeO, where the increments are nearly 8 × 10 − 3 for UO 2 and 2 × 10 − 5 for UO 2 -10vol%BeO on cladding inner surface. From 1100 K to 1500 K, the increment is noticeable for UO 2 -10vol%BeO compared with that from 700 to 1100 K in Figures 4B,D. We could see the volume fraction has increased around 6 × 10 − 5 on the inner-surface for Figures 4D,F. The third group describes the variations of volume fraction of β-phase when the gap size changes. The other values are 0,6-the coefficient multiplied with the heat transfer coefficient between cladding and 30,000 W/m -the linear power rate, 1500 K-the temperature at the extremities of the fuel element.
In this group of calculations, we have selected three different magnitudes of gap size to observe the variations in the volume fraction of β-phase on cladding inner-surface and outer surface. We note that 8e−5 m in Figures 5A,B is the initial gap size. From Figures 5A,C, E, the gap size has increased tenfold and a 100 times, the volume fraction on cladding inner-surface has increased from 0.02 to 0.036 and from 0.02 to 0.052, respectively, while we nearly see no difference between Figures 5B,D. In Figures 5D,F, when the gap size had increased a hundred times, the volume fraction of β-phase on cladding inner-surface for UO 2 -10vol%BeO had increased a lot from 2.8 × 10 − 14 to 3.9 × 10 − 4 , which is much more apparent than the variation in Figures 5B,D. Moreover, for UO 2 -10vol%BeO, the gap of volume fractions between two cladding surfaces has increased significantly from Figures 5D-F. The last group is to observe the influence of the coefficient multiplied with the heat transfer coefficient between cladding and coolant on the volume fraction of β-phase, with 30,000 W/m -the linear power rate, 1500 K -the temperature at the extremities of the fuel element, and 8e−7 m-the gap size.
In Figures 6A,C,E, we can see the volume fractions of β-phase on both the cladding inner-surface and outer-surface for UO 2 increased rapidly when the heat transfer coefficient decreased. More specifically, when the coefficient decreased from 1 to 0.8, the volume fraction of β-phase for UO 2 on cladding inner-surface increased from 3.5 × 10 − 4 to 1.5 × 10 − 3 . When the coefficient decreased 0.2 again, i.e., from 0.8 to 0.6, the volume fraction of β-phase for UO 2 on cladding inner-surface increased from 1.5 × 10 − 3 to 5.2 × 10 − 2 . For Figures 6B, D, F, the volume fractions  Figures 7D-F, the volume fraction of β-phase on cladding inner-surface had increased from 1.8 × 10 − 7 to 2.1 × 10 − 6 when the coefficient decreased from 1.0 to 0.8. In Figures 6B,D, the volume fraction of β-phase for UO 2 -10vol%BeO on cladding innersurface had increased from 2.1 × 10 − 6 to 3.9 × 10 − 4 when the coefficient decreased from 0.8 to 0.6. Compared with UO 2 , the volume fractions for UO 2 -10vol%BeO are much lower, and they are off by three orders of magnitude.
From all the figures above (i.e., Figures 3-6), we found that the highest volume fractions of β-phase are only 0.052 for UO 2 and 3.9 × 10 − 4 for UO 2 -10vol%BeO, which are still not high enough to reach the significant phase transformation interval. To simulate a severe condition, we reset the four parameters mentioned above: the linear power rate, the temperature at the extremities of the fuel element, the gap size, and the heat transfer coefficient to simulate again.
This time we set the linear power rate to 60,000 W/m according to Westinghouse Technology Systems Manual by NRC (The U.S. Nuclear Regulatory Commission, 2011). We decreased the gap size from 8e-7 to 1e-7 m because we have seen a relatively significant influence of gap size on the volume fraction of the new phase (β-phase). We can also decrease the coefficient of heat transfer by multiplying a number smaller (this time, we choose 0.4) for the same reason. Nevertheless, we could keep the temperature at the extremities of fuel element 1500 K because it seems that this temperature has little influence on the volume fraction of the β-phase. The results are shown below.
As we can see, the volume fraction of the β-phase on cladding inner-surface for UO 2 in Figure 7A has reached nearly 0.4, which means the inner surface of cladding had reached the significant phase transformation interval.
The result for UO 2 -10vol%BeO in Figure 7B, comparing with the results before (i.e., Figures 3-6), showed a significant increase in the volume fraction of β-phase, reaching nearly 0.01. However, the volume fraction is still minimal, which showed that the phase transformation did not significantly occur on both the inner and outer surfaces of cladding for UO 2 -10vol%BeO.
To see more clearly the variation of the β-phase volume fraction with axial direction on cladding inner-surface, in Figures 8A,B, we found that the β-phase volume fraction decreases with the cladding radius increase.

ANALYSIS AND DISCUSSION
As seen in the results above, we found that all parameters have an evident influence on the volume fraction of the favored phase except the temperature at the extremities of the fuel element, which has little impact. In the results of Figures 3-6, the highest volume fractions of β-phase are 0.052 for UO 2 and 3.9 × 10 − 4 for UO 2 -10vol%BeO, respectively. These two values have a significant difference of two orders of magnitude, which showed the BeO addition's significant influence on the fuelcladding system. For the last case in Figure 7, the volume fractions of β-phase on cladding inner-surface are around 0.39 for UO 2 and 0.0094 for UO 2 -10vol%BeO, where there is still a significant difference between them.
According to the results in (Liu et al., 2015), the addition of BeO decreased fuel temperature by increasing the fuel thermal conductivity. However, we found that the temperature distribution in the cladding (not the fuel) for both UO 2 and UO 2 -10vol%BeO is very close to our results. Although the temperature is close, the volume fraction has a significant difference. According to Massih in (Massih, 2009), the heating rate Q affects the position of the transformed volume fraction. More specifically, when the absolute value of Q increases, the graph of volume fraction of β-phase will move to a higher temperature position. Thus from our results presented in Figures 3-7, we can deduce that the heating rate Q for UO 2 -10vol%BeO is higher than that for UO 2 because the volume fraction for UO 2 -10vol%BeO is lower than that for UO 2 under the similar temperature. We may also deduce that the addition of BeO affects the temperature variation, which is the heating rate Q.

CONCLUSION AND PROSPECTS
In summary, this model simulates the phase transformation by presenting the volume fraction of the favored phase (β − phase) on the cladding inner and outer surfaces. In the case of Figure 7, which has the most severe LOCA conditions in all the groups of calculations, the addition of 10vol%BeO in the UO 2 fuel decreased the phase transformation effect a lot, which shows that the addition of BeO may increase the rate parameter for temperature under heating. All the results for UO 2 -10vol%BeO showed that no significant phase transformation was observed in the Zircaloy-4 cladding with UO 2 -BeO enhanced thermal conductivity nuclear fuel during the existing loss-of-coolant accident conditions.
For further study on phase transformation behavior in Zircaloy-4 cladding with the UO 2 -BeO enhanced thermal conductivity nuclear fuel, adjusting parameters to simulate a more realistic loss-of-coolant accident would be the right choice. In this model, we have used a constant to represent the variation of hydrogen concentration. To analog a more realistic phase transformation, implanting the model of hydrogen pickup and ydride formation would be beneficial. The influence of excess oxygen due to oxidation is another important factor affecting phase transformation behavior in Zircaloy-4 cladding and needs to be modeled in the future work.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
The work is mainly done by WL. LQ contributed to the accomplishment of this work, and WZ is the advisor of WL.