Solid state phase transformation kinetics in Zr-base alloys

We present a kinetic model for solid state phase transformation (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \rightleftharpoons \beta$$\end{document}α⇌β) of common zirconium alloys used as fuel cladding material in light water reactors. The model computes the relative amounts of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α phase fraction as a function of time or temperature in the alloys. The model accounts for the influence of excess oxygen (due to oxidation) and hydrogen concentration (due to hydrogen pickup) on phase transformation kinetics. Two variants of the model denoted by A and B are presented. Model A is suitable for simulation of laboratory experiments in which the heating/cooling rate is constant and is prescribed. Model B is more generic. We compare the results of our model computations, for both A and B variants, with accessible experimental data reported in the literature covering heating/cooling rates of up to 100 K/s. The results of our comparison are satisfactory, especially for model A. Our model B is intended for implementation in fuel rod behavior computer programs, applicable to a reactor accident situation, in which the Zr-based fuel cladding may go through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \rightleftharpoons \beta$$\end{document}α⇌β phase transformation.


Experimental data
In this section, we provide a brief survey of published experimental data on solid-to-solid ( α ⇋ β ) phase transformation of zirconium alloys used widely in the nuclear industry. In particular, we go over the data that we have utilized in our modeling and analysis. The alloys discussed here are listed in Table 1 with their nominal chemical compositions. We first recount the phase equilibria properties in regard to the excess oxygen concentration in the alloy and the hydrogen content of the alloy. Then we peruse data on phase transformation of the alloys during time varying temperature, namely continuous heating and cooling. The phase boundary temperature versus excess oxygen concentration for Zircaloy-4 was determined by Chung and Kassner 14 . These authors studied Zircaloy/oxygen phase diagram by resistivity measurements and by metallographic examination of equilibrated specimens. Chung and Kassner determined the location of phase boundary of oxidized Zircaloy-4 specimens annealed to form uniform oxygen distribution in the α phase region. Similar specimens were brought out to equilibrium temperatures in β phase and in ( α + β ) mixed phase and then quenched to obtain the corresponding phase boundaries. We have extracted phase boundary data, i.e. α/(α + β) and β/(α + β) transus, determined by metallographic examination from table 2 of 14 . These data are displayed here in Fig. 1 by circles.
The phase boundary temperature versus excess oxygen concentration for certain ZrNbO alloys for β/(α + β) transus was determined by Hunt and Niessen 15 . The authors investigated this transus for ZrNbO alloys containing up to 5 wt% niobium by metallographic examination. The oxygenated specimens were heated into the β phase region then quenched in water. Thereafter, they examined the microstructure of the samples for retained α phase in the transformed β matrix. Their data on Zr0.5wt%NbO were used to our intent for Zr1NbO, which are shown in Fig. 1 by blue diamonds. Because oxygen is an alloying element in Zircaloy-4 and Zr1NbO (Table 1), to obtain the excess oxygen, we have subtracted 0.12 w% from the absolute values of the aforementioned data.
A few data points for low excess oxygen concentrations are available for Zr1NbO close to the α/(α + β) boundary 16 . These data are shown in Fig. 1 as green diamonds. Toffolon et al. 16 determined these data by calorimetry starting from a Zr1wt%Nb binary alloy. A more detailed assessment of the phase equilibria of the Zr-Nb-O ternary system can be found in 8 .
Kinetics of solid state phase transformation in Zr-based alloys has been investigated experimentally over the years. Early studies include references [17][18][19] , which used electrical resistivity versus temperature measurements to determine the α/β phase boundaries in Zircaloy-2. Corchia and Righini studied phase transformation 17 in the temperature range of 1000-1400 K with heating rates in the range 200-1400 K/s and cooling rates 10-60 K/s. They concluded that, in nonequilibrium conditions, the kinetics of the phase transformation is strongly influenced by the microstructure of the specimen and by inhomogeneous distribution of the fast diffusing elements (Cr, Fe, Ni). Arias and coworkers employed much slower heating/cooling ( ≈ 0.05 K/s) to establish the  Figure 1. Phase boundary temperatures versus metal excess oxygen concentration for Zircaloy-4≡Zry4 (circles) 14 , Zr-0.5%NbO (blue diamonds) 15 , and Zr-1%Nb (green diamonds) 16 . The lines are empirical correlations for Zry4 and Zr1NbO alloy.
α/β phase boundaries 19 . They also evaluated 18 the influence of tin, iron and oxygen on the ( α + β ) field. The aforementioned early studies provide us with some experimental facts and are not in a form that could be used for direct comparison with our modeling.
A systematic study of α/β phase transformation of stress-relieved annealed (SRA) Zircaloy-4 and recrystallized (RX) Zr1NbO alloy was carried out by Forgeron et al. 11 . The study includes phase transformation under quasi-equilibrium conditions, i.e. very low heating/cooling rates, ranging from 0.002 to 0.33 K/s. The researchers used DSC, i.e. Differential Scanning Calorimetry, to determine the fraction of β (or α ) phase as a function of temperature in the alloys. They also made direct observation of microstructure by optical microscopy on samples annealed for a few hours at different temperatures in the (α + β) two-phase domain. The experimental uncertainty of the phase fraction measurement was estimated to be around ±0.05 or less.
In order to study the kinetics of phase transformation, Forgeron et al. 11 utilized dilatometric technique on small samples (12 mm in length) machined from as-fabricated cladding tubes. The dilatometric measurements for Zircaloy-4 samples included data with heating/cooling rates of 10 and 100 K/s but only with heating rates 10 and 100 K/s for Zr1NbO alloy. As in their DSC measurements, the uncertainty of the phase fraction measurement is reported to be ≤ 0.05 and ±10 K for measured temperatures.
In a subsequent study, Brachet et al. 20 investigated the influence of hydrogen content in the alloy on α/β phase transformation temperatures of SRA Zircaloy-4 and to some extent that of RX Zr1NbO. Most of the data presented in 20 are for Zircaloy-4L (Table 1)  No continuous set of data on β phase volume fraction versus temperature at different concentrations of hydrogen in Zircaloy-4 is reported in 20 . However, computations with some verifications for β phase volume fraction versus temperature at hydrogen concentrations of [H] ≈ 0 , 500 and 1000 wppm are reported, which we have used to compare the results of our model calculations in section "Results: computations versus measurements". Regarding kinetic experiments, the α/β phase transformation temperatures were measured by dilatometric technique for Zircaloy-4L during heating at rates of 10 K/s and 100 K/s as a function hydrogen content (in the range 10-970 wppm) at β volume fractions 0.1, 0.5 and 0.90 20 . The measurement uncertainty for the volume fraction and temperature is around ±0.05 and ±10 K, respectively.
Frechinet 21 investigated the α/β phase transformation of Zircaloy-4 by electrical resistivity, dilatometry and DSC. His resistivity measurements comprised heating/cooling ratings of ±5 K/s and ±0.8 K/s for which the relative amount of β phase fraction as a function of temperature was determined. Frechinet's 21 dilatometric measurements were carried out at a heating rate of 0.8 K/s as a kind of verification of his resistivity measurements. Some discrepancy between the two methods can be observed in regard to the evolution of β phase. His DSC measurements covered heating rates of 0.02, 0.08 and 0.33 K/s, which can be considered as equilibrium or quasi-equilibrium conditions for phase transformation.
A comparative study of the phase transformation behavior of SRA Zircaloy-4 and a RX Zr1NbO alloy at low heating rates is reported in 9 . Kaddour et al. 9 used both resistivity and DSC to determine the evolution of β phase during heating as a function of temperature for Zircaloy-4, in which the heating rates utilized were ≈ 0.8 K/s (resistivity) and ≈ 0.3 K/s (DSC). For Zr1NbO alloy, in addition to resistivity and DSC, they used image analysis, i.e. optical microscopy to determine the volume fraction of β phase as a function of temperature. The results obtained by these different techniques are in fair agreement. The onset of the α → β transformation temperature determined by resistivity was about 1093 K for Zircaloy-4 and 1043 K for Zr1NbO alloy.
More recent investigations of α → β phase transformation kinetics include Nguyen et al's 22,23 in-situ synchrotron X-ray diffraction (SXRD) measurements of cold-worked (CW) Zircaloy-4L subjected to heating rates in the range of 10 to 100 K/s, and electrical resistivity measurements of CW Zircaloy-4L and RX Zircaloy-4L at 100 K/s 22,24 , supporting their SXRD measurements. In a new study, Jailin et al. 25 using dilatometric technique measured the α → β phase transformation kinetics of SRA Zircaloy-4L subjected to very high heating rates, i.e. from 50 to 2000 K/s. They found that from thermal equilibrium to 500 K/s, the temperature at the start of transformation shifts toward higher values with increasing heating rates. However, above 500 K/s, the effect of www.nature.com/scientificreports/ heating rate saturates. Moreover, they found that the temperature at which the transformation terminates, i.e. when the alloys is in complete β phase, is not much affected by heating rate. Some phase transformation kinetic measurements on related Zr-Nb alloys have also been reported in the literature 23,26 . Nguyen 23 has reported in-situ SXRD measurements of α to β phase transformation of RX ZIRLO (Table 1) at a heating rate of 10 K/s, which indicates that the starting α → β phase transformation temperature is around 1040 K compared to ≈ 1130 K for CW Zircaloy-4, measured at the same heating rate. However, at about 1170 K the two alloys' β fractions become roughly equal, i.e., ≈ 0.4 . Běnes et al. 26 examined the α/β transformation of a type of Zr1NbO alloy, namely the E110 alloy, with an oxygen concentration of about 0.05 wt% (cf. Table 1). They measured the temperature of the transformation at low heating rates in the range of 0.05 to 0.33 K/s by a DSC technique. Their measurements indicate a weak shift (increase) in the onset of α to ( α + β ) transition temperature with heating rate. The equilibrium onset transition temperature for this alloy was put as T α/(α+β) = 1020 K.

Models
In this section, we describe the models that are used to evaluate the experimental data discussed in the preceding section. The models comprise both thermal statics (equilibrium) and kinetics (time varying) for α ⇋ β phase transformation in Zr-alloys (Table 1) for the intend of computing the β (or α ) phase volume fractions as a function of time and temperature. Portions of the models were presented in our earlier publications 12, 13, 27 . Thermal equilibrium. We assume that the volume fraction β-phase material under thermal equilibrium, denoted as y eq , depends on the Zr-alloy under consideration, temperature T, and excess oxygen and hydrogen concentrations x = (x O , x H ) in the metal 12, 13 , namely where T mid and T span are material specific functions of the cladding metal layer excess oxygen concentration (>1200 wppm), and hydrogen concentration (> 10 wppm). They are related to the middle and the span of the mixed-phase temperature region at equilibrium. They are calculated from the mixed phase lower and upper boundary temperatures, T α and T β , through For the cladding materials considered (Table 1), the phase boundary temperatures T α and T β can be written as where x O and x H are the weight fraction of excess oxygen and hydrogen in the cladding metal respectively, and A 1 , A 2 and m are constants given in Table 2. Equation 4 treats x O and x H as independent variables with no interaction between them. Furthermore, the relation is appropriate for x O < 0.01 and x H < 0.001 . The phase boundary temperatures versus excess oxygen concentration, calculated through Eq. (4), are compared with experimental data for Zircaloy-4, Zr-0.5%NbO and Zr-1%Nb in Fig. 1.

Kinetics.
Our approach to kinetics of the α/β phase transformation in zirconium alloys is similar to the method used by Leblond and Devaux 28 for diffusion controlled transformations in steels. We assume that there are only two phases, α and β , and we denote the volume fraction of phase β by the variable y. We also posit unique phase transformations α ⇋ (α + β) ⇋ β occurring at specified temperatures. Furthermore, we assume that the kinetics equation furnishes the equilibrium proportion y eq discussed in the foregoing section.
The general form of the kinetics or the time derivative of y, ẏ = dy/dt , depends only on y and T, namely 3 .

Alloy type
Phase boundary temperature T 0 (K) A 1 (K) where f is a function of thermal history T(t) and y is the volume fraction of β phase, and in thermal equilibrium f (T, y = y eq ) = 0 . In modeling the kinetics of the α/β transformation, two variants of the model are presented here, which we call model A and model B. In model A, a key parameter entering the kinetic equation is temperature rate of change ( q = dT/dt ), i.e. in principle ẏ = f (q, T, y) , whereas in model B, q does not appear explicitly in this equation. In the forthcoming two subsections, we will delineate the equations used for these two models.
Model A. We consider a simplest equation of type (5), by assuming y = y eq (T) is a steady-state solution at constant temperature, i.e. f (T, y = y eq ) = 0 for each T. We obtain f by assuming y is not too far from thermal equilibrium; thereby we Taylor expand f up to the first power Since the first term on the right-hand side is null, we can write Eq. (5) as where τ is a relaxation time, which is temperature dependent and formally is defined as τ (T) is a positive quantity and characterizes the duration of phase transformation. Equation (7) assumes initial values y(0) = 0 for heating and y(0) = 1 for cooling. For model A during heat-up, τ is calculated through an empirical relation by 13 where τ is in s, T is in kelvin, and q = dT/dt is the heating rate in K/s. Furthermore, y eq (T) is given by Eq. (1). During cooling ( β → α ), however, as the temperature is lowered ( q < 0 ), τ rapidly increases using this formula. And at sufficiently low temperatures, it blows up thus rendering artificially ẏ = 0 in Eq. (7). Surely the formula (9) requires a constraint during cooling to get meaningful results. Hence, in our computer program, we have modified Eq. (9) according to that is, we have introduced a lower-bound temperature constraint T b to avoid unphysically slow phase transformation during rapid cooling to low temperature. Choosing T b = 1000 K, gives reasonable results for experiments performed with cooling rates from −5 to −100 K/s; other choices are also possible. In model A, the heating/ cooling rate q also enters the expressions used for the phase boundary temperatures T α and T β 13 . Model A may be used for an experiment or a fabrication process when the heating/cooling rate q is prescribed in advance. As will be shown in section "Kinetics" below, the model describes measured data satisfactorily, i.e. it is engineered to do so.
Model B. As noted in the foregoing subsection, the relaxation time in model A embraces the temperature rate q, thereby renders it ineffective for reactor fuel rod simulations, where temperature history is not known a priori and may vary nonuniformly and rapidly during transients. In order to remedy this shortcoming, in model B we have extended equation (7) to a higher order and modified the relaxation time τ in Model A to fit the measured data suitably. Moreover, we have introduced a cutoff on τ.
The basic equations for the volume fraction of the β phase in model B are where the +sign is for heat-up ( y < y eq ) while the −sign is for cool-down ( y > y eq ), and y eq (T) is given by Eq.
(1). The relaxation time τ (in s) is given as where B, b and E are constants; see Table 3. During cooling, we have introduced an upper cutoff limit on time constant τ ≤ τ * to avoid excessively long relaxation times, namely we take where τ * = 6 s is used in all our computations (Table 3). It will be shown in section "Kinetics" below that model B reproduces measured data fairly satisfactorily.

Results: computations versus measurements
In this section, we use the models described in the foregoing section to evaluate some of the data alluded in section "Experimental data". We compare model computations of the volume fraction of β phase with experimental data through a number of figures.
Thermal equilibrium. In Fig. 3, the temperature dependence of the equilibrium volume fraction of the β phase Zircaloy-4 and Zr1NbO alloy using Eq. (1) is shown against the experimental data (markers) extracted from ref. 11 . Solid lines are calculated with the model parameters listed in Table 2 whereas the dashed lines are those from ref. 13 ; here x O = x H = 0 . As can be seen, the agreement is quite good. It is noted that the transition temperatures for Zr1NbO are lower than for Zircaloy-4 and the volume fraction of β phase shifts to a lower temperature because niobium acts as a β phase stabilizer. The influence of hydrogen content on phase transformation of Zircaloy-4 is displayed in Fig. 4. It is seen that the results of our computations for hydrogenated Zircaloy-4 are in fair agreement with experimental data and analysis in 20 . However, the deviation for as-fabricated Zircaloy-4 (H≈ 0 wppm) is salient. This is because, our   Table 2 whereas the dashed lines are those from ref. 13 .  Table 2 whereas the dashed lines are those extracted from figure 12 of ref. 20 , which are based on analysis of experimental data. The disks show measured data for H ≈ 0 wppm from ref. 11 11 ; discs in the figure. Note that as hydrogen concentration is increased, the α → (α + β) transition temperature shifts to a lower value and the volume fraction of the favored phase is moved to lower temperatures. We should note that, to our knowledge, no similar measurements on the effect of excess oxygen on phase transformation of zirconium alloys have been reported in the literature.
Kinetics. Model A. The results of model A computations were presented earlier 13 , however, with some clarification included here that was missing in that reference; e.g. usage of Eq. (10), and additional data. We perform our computations using this model and parameter values as in 13 for the sake of a benchmark with model B, which is the main contribution of our paper. We evaluate some of the measurements discussed in section "Experimental data", namely data in 11 with heating/cooling rates q = dT/dt = ±10, ±100 K/s and in 21 with q = dT/dt = ±5 K/s. We solve Eq. (7) with initial conditions y(0) = 0 on heating and y(0) = 1 on cooling together with relations (9)-(10) for the relaxation times (with T b = 25|q| ) by using an explicit Runge-Kutta method. The results for the fraction of the transformed phase for Zircaloy-4 under heating/cooling are shown in Fig. 5. The markers denote measured data while the lines are calculations according to model A. The corresponding calculations and data for Zr1NbO 11 subject to heating rates q = 10 and 100 K/s are shown in Fig. 6. It is seen that as q is increased, the transformed volume fraction shifts to higher/lower temperatures under heating/cooling, respectively.

Model B.
Here, we solve Eq. (11) with initial conditions y(0) = 0 on heating and y(0) = 1 on cooling together with relations (12)-(13) for the relaxation time by using an explicit Runge-Kutta method. The results of model B computations, the volume fraction of β phase versus temperature, against measured data for Zircaloy-4 (markers) which were performed at thermal rates of ±0.8(≈ ±1) ± 5, ±10, ±100 K/s and at thermal equilibrium (cf. section "Experimental data") are shown in Fig. 7. As can be seen from this figure, the agreement between calculations and measurements are satisfactory, except for the case of heating at 100 K/s, at which computations underestimate the measurements at temperatures ≈ 1230 and ≈ 1280 K. However, the trend in α → β phase transition is qualitatively featured. The corresponding calculations and data for Zr1NbO 11 subject to heating rates of +10 , +100 K/s and at thermal equilibrium are plotted in Fig. 8; cf. Fig. 6. As can be seen, the concordance between computations and measurements is fair for this alloy.
To illustrate the phase transformation kinetics expected for reactor fuel cladding under a postulated loss-ofcoolant accident (LOCA) in a pressurized water reactor (PWR) 29 , we apply model B to a cladding temperature history measured in a test within the QUENCH-LOCA experimental program at the Karlsruhe Institute of Technology, Germany. This program comprised seven out-of-reactor LOCA-simulation tests on bundles with 21 electrically heated and extensively instrumented fuel rod simulators with various Zr-base cladding materials 30 . The rods in the considered test, QUENCH-L1, were clad with as-received Zircaloy-4 31 . More specifically, we consider the cladding temperature history measured for rod 4 at an elevation 450 mm from the bottom of the heated section (thermocouple TFS 4/8 31 ), which is shown in the left panel of Fig. 9. The temperature history is typical for a postulated large-break LOCA in a PWR: An initial heat-up phase, representing loss of primary coolant water at full reactor power, is followed by cooling in steam with a rate around 3-4 Ks −1 . At t=207 s, the axial steam flow through the fuel bundle is blocked by liquid water that reaches the lower (cold) part of the bundle during re-flood of the reactor pressure vessel by the emergency core cooling system. The fuel cladding temperature then rises to about 870 K until the water front reaches the thermocouple position and the cladding is rapidly quenched. For illustration, this temperature history is simulated with model B. The calculated results are shown in the right panel of Fig. 9 for two different cases: For the 'no oxidation' case, the effects of high-temperature oxidation are neglected, keeping x O =x H =0 throughout the simulated test. For the ' oxidation' case, the evolution of x O during the test is calculated through the correlation by Leistikow and Schanz 32 , while hydrogen pick-up during the test is neglected ( x H =0). The calculated results in Fig. 9 suggest significant kinetic effects on the cladding phase composition during a typical PWR large-break LOCA. The calculated effects of cladding high-temperature oxidation are, however, moderate. The calculated excess oxygen concentration in the cladding metal is merely 280 Equilibrium q =10 K/s q =100 K/s Figure 6. Volume fraction of β phase as a function of temperature calculated (lines) by using model A versus measured data (markers) 11 for Zr1NbO at heating rates q = dT/dt = 10 and 100 K/s and at thermal equilibrium.  Next, we consider the case of hydrogenated Zircaloy-4 and rely on experimental data reported in 20 . As mentioned in section "Experimental data", no continuous set of data on β phase volume fraction versus temperature at different concentrations of hydrogen in Zircaloy-4 have been reported in 20 . Brachet et al. 20 , besides equilibrium data, have presented data on heating at 10 K/s and 100 K/s for β phase volume fractions of 0.1, 0.5 and 0.9 for Zircaloy-4 containing hydrogen concentrations in the range ≈ 10 to 970 wppm. We have extracted these data from their Fig. 8 by grouping the data to four levels of hydrogen concentrations, namely, ≈ 12 wppm, 188 wppm, 500 wppm, and 920 wppm. The results are depicted as a contour plot in Fig. 10a for the heating rate of 100 K. In this figure, the β-fraction values are interpolated between the measured data points of 0.1, 0.5 and 0.9 for the sake of contiguity. The corresponding computations with model B are presented in Fig. 10b. It is seen that measurements and computations are in fair agreement, albeit the complexity of the phenomenon.
In order to recount the aforementioned experimental data on Zircaloy4+H, we have extended the relaxation time relation of phase transformation, i.e. Eq. (12), to account for hydrogen concentration as Figure 8. Volume fraction of β phase as a function of temperature calculated (lines) using model B versus measured data (markers) 11 for Zr1NbO at heating rates q = dT/dt = 10 and 100 K/s and at thermal equilibrium.   www.nature.com/scientificreports/ β phase as a function of temperature for several values of hydrogen concentration in Zircaloy-4 at heating rates of 10 K/s and 100 K/s. It is seen that as hydrogen content in the material is increased, the transition from α to β phase initiates at a lower temperature.
As has been noted the deviation between model B computations and measurements are most prominent for the case of Zircaloy-4 subjected to a heating rate of 100 K/s (cf. Fig. 7a) and even, to some extend, at a cooling rate of −100 K/s, cf. Fig. 7b. Thereupon, an uncertainty evaluation of measured data and their deviation from the model calculation are worthy of estimation. As alluded in section "Experimental data", the uncertainty in measured volume fraction of β phase as a function of temperature is estimated to be about σ = ±0.05 . Moreover, there is an uncertainty on temperature measurements of around ±10 K. In Fig. 12a,b, we have plotted the measured data for Zircaloy-4 for q = dT/dt = ±10 K/s and (b) q = dT/dt = ±100 K/s including uncertainties (bars) of σ = ±0.05 on β phase volume fraction and ±10 K on temperature. We have also included in these figures the corresponding model B computations for comparison. As can be seen the model B computations for the case of q = ±10 K/s are reasonably satisfactory. However, for the case of q = ±100 K/s, the model B deviates from measured β phase fractions up to around 0.15 in certain temperature ranges, Fig. 12b. For Zr1NbO at heating rates q = +10 and +100 K/s, Fig. 8, the data are reproduced fairly satisfactory.

Discussion
The phase boundary temperatures versus excess oxygen concentration, calculated through Eq. (4), were compared with experimental data for Zircaloy-4 and Zr-0.5%NbO cladding in Fig. 1. The data were obtained from metallographic analyses of cladding tube samples, quenched from different equilibrium temperatures. We should note that T β , calculated through Eq. (4), may exceed the Zircaloy or Zr1NbO melting temperatures at high oxygen concentrations. The melting (solidus) temperature of Zircaloy-4 increases from about 2025 K at zero excess oxygen concentration to about 2320 K for x O ≥ 0.04 36 .
Turning now to the thermal equilibrium relation (1), we stated this formula without any theoretical justification or explanation, albeit that it captures experimental data adequately. Let's write the hyperbolic tangent part of this equation (the first term 1/2 is just a constant shift on the ordinate) in a simplified form where T 0 ≡ T mid and ξ ≡ T span / √ 2 . Indeed, Eq. (15) is a stationary solution of the Cahn-Hilliard equation of phase field theory in the temperature domain, namely where f (y) = y 4 − y 2 /2 is a double-well function (potential) of y; see e.g. 37,38 . This can be verified either by direct substitution or by integrating Eq. (16) with appropriate boundary conditions, e.g. y 0 (T 0 ) = 0 and y ′ 0 (T 0 ) = √ 2/(4ξ) . Differentiation of Eq. (15) with respect to T gives www.nature.com/scientificreports/ It is seen that the parameter ξ appears both in the argument of the sech-function and as its coefficient in inverse.
In phase field theory parlance, ξ is the correlation length, here the temperature span of the coexisting phase. It plays an important role in the analysis of phase transitions. Figure 13 shows schematic plots of y 0 and ζ as a function of temperature. It is seen that ζ is highly localized in temperature with a maximum at the center of the bump, T = T 0 , and falling rapidly to zero for |T − T 0 | > √ 2ξ. Experimental data indicate that the level of impurities in the material affects the width of ξ . For example, as hydrogen concentration in Zircaloy-4 is increased, so does ξ , namely at x H ≈ 0 wppm, ξ ≈ 27 K, whereas at x H = 500 wppm, ξ ≈ 37 K; using relation (3) and data in Table 2 for Zircaloy-4 with ξ ≡ T span / √ 2 . In regard to the influence of oxygen on the Zr-O phase equilibrium, oxygen is an " α stabilizer". It expands the α phase region of the phase diagram by formation of an interstitial solid solution in zirconium 4 . That is, excess oxygen widens ξ . For example, at an excess oxygen concentration of 1wt%, ξ ≈ 80 K; see relation (3) and data in Table 2 for Zircaloy-4. We may also compare ξ for Zr1NbO versus Zircaloy-4. Again using data in Table 2, we see that for x H ≈ 0 wppm, ξ Zr1NbO ≈ 34 K, i.e. somewhat wider than that for Zircaloy-4. This is expected, since Nb is a β stabilizer, meaning that it widens the (α + β ) domain in Zr-based alloys 4 .
Regarding the issue of hydrogen in Zircaloy, we should say that in the studied range of temperature and hydrogen concentration, cf. Figs. 4, 10, 11, hydrogen is primarily in solid solution, i.e. not in zirconium hydride phases. This can be checked from data on the temperature dependence of hydrogen solubility limit in zirconium alloys 39,40 and the ensuing phase diagram 39 . It is, however, worth recalling that the temperature dependence of the hydrogen solubility limits for hydride dissolution (during heat-up) and precipitation (during cool-down) are different in Zircaloys 40 ; see Fig. 14.
Turning now into the issues in kinetics, the α ⇌ β phase transformation in zirconium base alloys comprise the process of nucleation and growth. For instance, on heating from α phase, the transformation involves β-nucleation, β-growth and β-completion and vice versa; i.e. α-nucleation, α-growth and α-completion. A commonly used macroscopic model for this process is the Kolmogorov-Johnson-Mehl-Avrami [41][42][43][44] or KJMA description 45,46 . The KJMA model assumes that tiny spherical grains nucleate at a constant rate per unit volume in the metastable phase and grow isotropically at constant velocity once formed. As first formulated by Kolmogorov 41 , the volume fraction of the favored phase at time t can be expressed as where R m ≡ vt is the grain radius obeying an algebraic growth law. The nucleation process in Kolmogorov's theory is considered as homogeneous, meaning that, it occurs uniformly throughout the metastable phase. In most cases, however, nucleation is heterogeneous. That is, nuclei primarily form at impurities and defects in the crystal, which are present before the transformation started. A simple model in the spirit of KJMA has been discussed in 47 , which places nuclei randomly with density ρ throughout the solid. When the phase transformation is initiated, spherical grains form at each nuclei and grow isotropically with velocity v. Both the KJMA model and its simplified heterogeneous nucleation variant 47 are a mean-field type theory. Meaning that, there is no spatial variation in φ through the solid, i.e. ∇φ = 0. In 12 , we made a detailed comparison between Eq. (18) and our model A. To sum up that account, we note that under isothermal condition, by putting φ = y/y eq , we can write Eq. (7)  where K is a function of temperature designating the nucleation and growth rates, and n is a numerical exponent (Avrami exponent) whose value may vary from ∼ 1 to 4. According to 48 if there is no change in the nucleation mechanism, n is independent of temperature. So, as can be inferred from Eq. (20), our model A corresponds to n = 1 and K = 1/τ.
The behavior of the model B, as described by Eq. (21), is slightly more complicated. Solution to Eq. (21) subject to φ(0) = 0 under isothermal condition reads  Fig. 15 for several values of b. As can be seen from these plots, the behavior for b = 0 and b = 0.3 is rather close. So, we expect even for b = 0.3 the Avrami exponent should be around one. It is worth mentioning here that the space dimension d (or n ≡ d + 1 ) appearing in the foregoing equations depends on genre of nucleation and/or growth. Avrami 44 attributes n = 1 to 1-d growth (needles), n = 2 to 2-d growth (plates), n = 3 to 3-d growth (spheres or polyhedra), and n = 4 to nucleation and 3-d growth 49 . Whereas Cahn 50 , in a subsequent theoretical study of the Avrami model, relations of type (23), concludes that after site saturation, the stage in a reaction at which all potential nucleation sites are consumed, the reaction proceeds by growth alone. In Cahn's analysis n = 1 describes nucleation on grain surfaces, n = 2 nucleation on grain edges and n = 3 nucleation on grain corners 49 .
We should note that our model B is tailored for fuel rod behavior computer program application, e.g. the FRAPTRAN program 51, 52 , where its parameters are independent of thermal rate dT/dt, contrary to model A. Its performance against experimental data is fairly good, considering the uncertainty in the measured measured phase fraction and temperature, which is about ±0.05 and ±10 K, respectively. The model, despite some theoretical justification, is basically empirically based with adjustable parameter constants. Its results perhaps can be improved further by more tune-up of those constants. The model should be applicable to thermal rating in the range dT/dt = −100 to +100 K/s and hydrogen contents in the alloys up to around 1000 wppm; however, with more uncertainty toward higher thermal rates.

Summary and conclusions
In this paper we have introduced a new model (model B) for the kinetics of α ⇋ β transformation in zirconium alloys used in the nuclear industry as fuel cladding material. We first surveyed pertinent experimental data for the alloys reported in the literature, section "Experimental data". They constitute an empirical basis for our model. We then presented our model B together with our previously published model A, the latter for the sake of comparison and completeness, section "Models". The main difference between models A and B is that model B is formulated such that it can be applied to arbitrary heating/cooling histories, while model A is restricted to cases with a constant heating/cooling rate. Moreover, the kinetic equation for model B extends to a higher order term in Taylor expansion about the equilibrium. Both models are extended to account for the effect of excess oxygen and hydrogen content in the alloy. Surplus in oxygen may occur under an accident situation due to the high temperature metal-water reaction as the cladding temperature raises say from roughly 650 K to 1470 K. The increase in the hydrogen content of cladding from its as-fabricated value of around 10 wppm occurs during normal reactor operation and is enhanced during an accident, again due to the metal-water reaction. All the experiments which we have evaluated with our models were performed on unirradiated materials. The effect of irradiation on α ⇋ β is believed not to be significant, at least for the rather low heating rates expected in light water reactor LOCA; cf. Fig. 9. The thermodynamic stability of crystal structures is the dominant factor for phase transformation. During α → β phase transformation as the cladding temperature raises, the irradiation-induced defects would mostly be removed as the recrystallizing front is swept through the lattice 53 .
We compared the results of our model computations with experimental data as regards the relative amounts of β (or α ) phase fractions as a function of temperature in the alloys, section "Results: computations versus measurements". Our computations are mainly on Zircaloy-4, including hydrogenated material, and to a limited degree on Zr1NbO material because of dearth of available measurements for this alloy. The agreement between calculations and measurements is fairly satisfactory considering the complexity in the details and the measurement uncertainties. Finally, in section "Discussion" we made a theoretical analysis of the thermal equilibrium expression used and qualified for the transformation. In that section, we also briefly discussed the kinetics in regard to nucleation and growth, and delineated the range of applicability of the putative model.
Our model B is basically empirically based, with adjustable parameter constants. It is intended for implementation in fuel rod behavior computer programs and reproduces available experimental data fairly satisfactory.