Mechanical–chemical coupling phase-field modeling for inhomogeneous oxidation of zirconium induced by stress–oxidation interaction

A phase-field model is proposed to study the inhomogeneous growth of zirconia induced by the stress–oxidation interaction, which captures the complex interplay among diffusion, oxidation kinetics, interfacial morphology evolution, and stress variation in an oxidation process. Through this numerical model, many experimentally observed but insufficiently understood phenomena can be well explained. Specifically, the numerical simulations reveal quantitatively the causes of interface roughening or smoothening during the inward oxide growth, the roughness-dependent oxide growth rate, and the nucleation sites of premature cracking. These numerical findings can be used as the theoretical references for the improving the durability of oxide scale and prolonging the service life of zirconium-based alloy cladding used in the nuclear power plant.


INTRODUCTION
As of 2018, more than 400 light-water reactors (LWRs) had been put into operation, and currently supply around 8% of the total global demand for electricity 1 . They are by far the most common reactors in the world. However, the performance of LWRs is severely limited by the high-temperature corrosion of metallic components. Of particular concern is the oxidation of the zirconium-based alloy used in the cladding of the fuel rods that operate at about 270-320°C in water. In this aqueous environment, oxidation occurs due to water dissociation 2 .
To protect the zirconium alloy cladding from further oxidation, a dense protective oxide layer (mainly zirconia) must be formed. This oxide layer inhibits the early stage of oxidation kinetics (a period of several hundred days) 3 by acting as a barrier to oxygen diffusion and hydrogen uptake (HU). However, in the early stage of oxidation, the change in specific volume from zirconium to zirconia can result in significant volumetric eigenstrain and therefore considerable stress increase (generally called growth stress) in the oxide layer, which easily leads to cracking of the oxide scale 4 .
The growth stresses in zirconia can be determined through measuring the curvature change of a specimen or from the peak shift in X-ray diffraction patterns or Raman spectra. Such growth stresses have been reported as ranging from a few hundred MPa to over 5 GPa 2 . The stress is also correlated with the interface morphology of the oxide layer. Some experimental observations 5,6 and numerical studies 6 have shown that if the metal-oxide interface is rough, the local stress in the oxide will be tensile in the direction perpendicular to the mean interface plane, and its magnitude and distribution are affect by the mechanical properties and vary with the morphological evolution of the oxide scale. This tensile stress is the reason why cracks initiate and propagate along the oxide-metal interface, as experimentally observed in refs 5,6 . In addition, Platt et al. 7 found that an initially rough interface became gradually smoother and an initially planar interface became rougher, and that both approached a similar level of roughness after a few hundred days of oxidation. It was also noted that the morphology of metal surface affected oxide growth rate [7][8][9][10] . An initially rough metal surface results in rapid oxide growth, whereas an initially planar meal surface leads to slow growth, as observed by Platt et al. 7 . Platt et al. 7 found that it was the growth stress, not the oxide geometry itself, that affected the oxidation rate, because when the samples were stress relieved, the oxide growth became insignificantly dependent on interface roughness 10 .
Previous experimental studies 7-10 have clearly indicated the interaction between stresses and oxidation kinetics. In the course of studying this interaction, the theory of the thermochemical equilibrium of a stressed solid was established by Larché and Cahn 11 . They assumed that species diffusion was affected by the stresses arising from compositional change due to diffusion, phase transformation, and chemical reaction, which gives rise to a hydrostatic-stress-dependent chemical potential. Following Larché and Cahn 11 , Krishnamurthy et al. 12 assumed that oxidation rate was a function of hydrostatic stress, and established a onedimensional (1D) theoretical model to study the interaction between stresses and oxidation kinetics. Krishnamurthy et al.'s model was further developed by Zhou et al. 13 and Wang et al. 14 . Zhou et al. 13 assumed that the chemical potential was dependent on both hydrostatic and deviatoric stresses, and dealt with large deformation in their calculation of growth stress. Wang et al. 14 considered inelastic deformation such as plasticity and creep. However, these theoretical works were all limited to 1D problems and could not be extended to higher dimensions, because they did not treat the evolution of interface morphology. To simulate a two-dimensional (2D) chemical-mechanical coupled problem, Loeffel et al. 15 used a field variable varying from 0 to 1 and representing oxide volume fraction to describe the transformation from metal to oxide. In this model, the variation rate of the field variable was related to the oxidation rate; and a generalized force conjugated to elastic energy density was derived and incorporated into the oxidation rate formula to reflect the contribution of elastic energy to oxidation kinetics. Lin et al. 16 15 to incorporate a hydrostatic-stress-dependent chemical potential 11 , which resulted in the stress-dependent absorption of oxygen molecules at the oxide surface. In contrast with Larché and Cahn's theory, Evans 17 assumed that the influence of stress on oxidation only takes place at the oxide-metal interface, that the oxygen diffusion through an oxide scale is due to the motion of vacancies nucleated at the interface by the reaction, and that the vacancy concentration is affected by stresses. Following Evans's theory, Saillard et al. 18 simulated the stress-induced nonuniform inward growth of oxide based on the sharp interface assumption; the interfacial moving problem was dealt with in their finite element (FE) simulation. However, their treatment of the interface required ad hoc assumptions about interfacial geometry, which significantly affected the accuracy of the simulation.
A phase-field (PF) model makes it possible to deal with complex morphology change without explicitly tracking interface positions. In a PF model, the interface has a finite thickness, which is associated with a smooth change of one or several order parameters. The change of these parameters can be coupled with various physical, chemical, and mechanical fields; therefore, the PF model can be utilized to predict many material processes. A few PF models of oxidation have been proposed. For example, Ammar et al. 19 proposed a PF model to describe diffusionmediated oxidation based on the assumption that the migration of the oxide-metal interface is driven by oxygen diffusion. Sherman and Voorhees et al. 20 proposed another PF model that used a Poisson equation to simulate the initial oxidation and take the electrostatic effect into account. Zaeem et al. 21 dealt with the stress-oxidation interaction, incorporating the effect of elastic energy on oxygen diffusion and establishing a Cahn-Hilliard type governing equation to describe the oxygen diffusion and oxide-metal phase transformation. Zhao et al. 22 further developed Zaeem's model to deal with inelastic deformation during hightemperature oxidation. In their models, however, stress only affected diffusivity. To cope with the influence of stress on oxidation kinetics, Lin et al. 23 recently proposed a PF model that included elastic strain energy in the formula for oxidation kinetics, and studied the inhomogeneous outward growth and surface roughening of the NiO scale.
To the best of our knowledge, there has not been a theoretical or numerical model to systematically explain the abovementioned phenomena of zirconium oxidation including the interface roughening or smoothening, the roughness-dependent oxide growth rate, and the premature spalling induced by growth stress. In order to unveil the fundamental cause of these phenomena, we present a PF model based on our previous model of mechanochemical corrosion 23 . The model tackles the effects of volumetric eigenstrain resulting from the change in specific volume when a metal transforms into oxide, the compositional mismatch strain induced by species diffusion, and the high-temperature creep of the stressed. The mechanical strain energy associated with these deformations is included in the Helmholtz free energy, along with the chemical and interfacial energy, to describe diffusion and reaction. The oxidation rate is expressed as a function of chemical potentials of reactants and products leading to a generalized Allen-Cahn equation. Coupled with the Allen-Cahn, diffusion--reaction, and mechanical equilibrium equations, the complete set of governing equations of PF model are established. The numerical results reveal the complex interplay among diffusion, oxidation kinetics, interfacial morphology evolution, and stress variation in an oxidation process, which can be furtherly regarded as the indications of what combination of properties (for example the initial surface roughness, creep behavior) should be optimized for improving the durability of oxide scale and prolonging the service life of zirconium-based alloy cladding used in the nuclear power plant. Note that the proposed PF model is not limited to zirconium but applicable to simulate inward oxide growth process in other metallic materials if the model parameters can be quantified.

Influence of eigenstrain
As mentioned earlier, growth stress affects the oxidation kinetics. To investigate this influence, the variations of oxide thickness against time for the different eigenstrains ε eg in = 0.005, 0.01, 0.02 are shown in Fig. 1a-c. It can be observed that oxide growth slows down as ε g in increases, because the compressive stresses lead to an increase in the chemical potential of the reaction products (as expressed in Eq. (41)), resulting in a decrease in the dimensionless reaction driving force Λ (see Eq. (47)) and therefore the reaction rate. In addition, the in-plane compressive stresses are nonuniform owing to the initial sinusoidal shape of the oxide scale. The stress concentrates near the peak of the interface, which induces a nonuniform growth of oxide, as shown in Fig. 1b, c. Such a roughening phenomenon has also been experimentally reported by Platt et al. 7 .
The time dependences of the average thickness of the oxide scale h ox and the roughness of interface Ra for the different eigenstrain, ε eg in = 0, 0.005, 0.01, 0.02, are shown in Fig. 2a, b, respectively. Note that if ε eg in is very small, e.g., ε eg in ≤ 0.005, the time dependence of the oxide thickness is approximately parabolic and obeys Wagner's theory. If ε eg in is large, i.e., ε eg in = 0.01-0.02, the oxide growth is no longer parabolic but subparabolic, which has also been reported by Zhou et al. 13 and Zhao et al. 22 . The simulation results are compared with the experimental data in Fig. 2a. A good agreement demonstrates that the proposed PF model describing the complex interplay among diffusion-reaction, interfacial evolution and stress variation can lead to a correct trend revealed in experiments. This model can therefore be used to study other parametric effects. The results clearly show that the roughness increases with the increase in eigenstrain; more notably, for small and large eigenstrains the evolutions of roughness are disparate, as shown in Fig. 2b. If ε eg in is very small, i.e., ε eg in ≤ 0.005, the roughness decreases slowly with time. In this case, diffusion rather than stress-induced nonuniform oxidation is the main cause of morphology evolution for the zirconia-zirconium interface, which results in interface flattening. If the eigenstrain is large, i.e., ε eg in = 0.01-0.02, stress-induced nonuniform oxidation dominates the morphology evolution. As shown in Fig. 2b, rapid roughening of zirconia-zirconium interface occurs in the early stage, followed by a slow increase of roughness. The decrease in the rate of roughening in the later stage is mainly owing to the mutual interference of neighboring oxide humps during growth, as illustrated in Fig. 2c. After a long period of oxidation, a dynamic balance between stress-induced interface roughening and interference-induced interface flattening can be achieved, which leads to a stable interface morphology, i.e., the roughness approaches an approximately constant magnitude.
Influence of chemical expansion and stress-dependent diffusivity The chemical expansion induced by composition mismatch is generally much smaller than the eigenstrain induced by metal-tooxide transformation. However, the former increases the influence of stresses on the chemical potential of the reactants, as derived in Eq. (39), which is a non-negligible effect. Figure 3a, b shows the growth of the oxide scale for two different chemical expansion coefficients, η O ox = 0.002 and 0.02. It can be seen that the oxide becomes slightly thicker and rougher with the increase in the chemical expansion coefficient. The larger chemical expansion elevates the chemical potential of the reactants (Eq. (39)) and drives the oxygen to flow from the compressive stress region to the tensile stress region (see Eq. (48)). The elevated chemical potential results in an increase in the dimensionless driving force (Eq. (47)) for oxidation and a higher oxidation rate. The change in the flow of oxygen makes their distribution less uniform, resulting in increased roughness.
If the diffusivity is also stress dependent, the oxide scale will become thinner and rougher, as shown in Fig. 3c. As Eqs. (49) and (50) show, compressive hydrostatic stress suppresses diffusion, while tensile stress promotes it. As the overall oxide scale is under compression, the diffusion of oxygen is suppressed, leading to reduced thickness. The stress-dependent diffusivity further leads to the slower diffusion of oxygen to the peak of the interface, where the most significant compressive stresses appear, leading to an increase in roughness. The effects of η O ox and ξ on the evolutions of the average thickness h ox and roughness Ra of oxide scale are further shown in Fig. 4. Note that the coefficients, η O ox = 0.002 and ξ = 0.05, are derived from refs 24,25 , respectively. Therefore, the simulation results shown in Fig. 4 indicate that the effect of compositional mismatch strain is insignificant because oxygen are small. This contrasts with the scenario of lithiation, in which chemical stress plays a very important role 26 . In contrast, the stress dependence of diffusivity has a more significant impact on thickness and roughness.

Influence of creep
As the nonuniform growth of oxide scale is mainly due to growth stress, creep brings about stress relaxation and should therefore mitigate this stress effect. The creep behavior of oxide depends on composition and microstructure. For zirconia, it can span several orders of magnitude due to variation in grain size or alloy composition 27 . To investigate the effect of creep, we vary the creep rate coefficient of zirconia in the range B ox -100B ox with the reference value B ox given in Supplementary Table 1 (see Supplementary information). Figure 5 shows the evolution of oxide scale over time for different creep rates, where ε eg in is set as 0.01. We observe that if the creep rate is small, i.e., in the cases B ox and 10B ox , the compressive stresses near the interface are not considerably reduced, although a significant stress relaxation occurs in the region far away from the interface (see Fig. 5a, b). In this case, the stresses induced by the eigenstrain still play a dominant role and cause interface roughening. If the creep rate coefficient is increased to 100B ox , the stress near the interface can be fully relaxed, resulting in an approximately stress-free state and a flat interface (see Fig. 5c). The simulation results of different creep rates are further summarized in Fig. 6. In terms of the thickness change of the oxide scale, larger creep rates always lead to greater thickness, as shown in Fig. 6a. We can also observe that the roughness first increases and then decreases, especially when the creep rate coefficient is 100B ox and eigenstrain is 0.02, as shown in Fig. 6b. This phenomenon can be regarded as competition between the roughening induced by stress generation and the flattening resulting from creep relaxation. At the early stage, rapid oxide growth leads to a significant increase in the stresses, which results in interface roughening, while after a long period of oxidation, the oxide growth rate decreases, which tends to slow down stress generation and causes creep relaxation to gradually dominate, resulting in interface flattening.
Influence of initial geometry Platt et al. 7 reported that the initial geometry of oxide scale significantly affected oxide growth and interface morphology. Their account of the effect of growth stresses can now be corroborated by the simulation model. Figure 7 shows the variations in the average thickness of the oxide scale and the roughness of the interface with the variations in the amplitude, A = 0.25-3 μm, and half-wavelength, W = 3-6 μm, of initial sinusoidal interface. Note that the thickness of the oxide scale increases as the initial roughness increases (see Fig. 7a), which agrees with the experimental results of Platt et al. 7 . The reason is that the increase in initial roughness reduces the region constrained by surrounding materials in the oxide scale, resulting in a larger low-stress region, which promotes oxide growth. The increase in initial half-wavelength, on contrary, increasing the constrained region, results in the suppression in the oxide growth (see Fig. 7a). The roughness evolution of the interface is also consistent with the experimental results reported by Platt et al. 7 (solid lines in Fig. 7b); i.e., the initially rough interface (A = 3 μm, W = 4 μm) becomes smoother, the initially flat interface (A = 0.25, 1.5 μm, W = 4 μm) becomes rougher, and after a long period of  oxidation, interface roughness approaches a constant magnitude whether it starts from a higher or lower value, as shown in Fig. 7b. This scenario can be attributed to competition between stressinduced interface roughening and interference-induced interface flattening, as illustrated in Fig. 2c. It is also observed that if halfwavelength is larger (W = 5 μm), the interference-induced interface flattening would be less significant, which results in a rougher interface (dish-line in Fig. 7b). For a more general case, the  evolution of roughness with the variation in amplitude, A = 0.25-3 μm, and half-wavelength W = 3-6 μm, are investigated, as shown in Fig. 7c. It is observed that evolution patten of interface roughness, in fact, is affected by the ratio, A/W, rather than the single parameter A or W. For the initial oxide scale with similar geometric properties, i.e., the ratios, A/W, are identical, the evolutions of roughness would have same patten. And the above-mention three evolution pattens of interface roughness can be identified by the initial ratio, A/W, with three regions: 0 < A/W ≤ 0.25, 0.25 < A/W ≤ 0.5, 0.5 < A/W. For the initial oxide scale with similar geometric properties, i.e., the ratios, A/W, are identical, the evolutions of roughness would have same patten. Impact of interface evolution on stress While growth stress induces interface roughening, the latter in turn causes variation in the stress field. One of the most important consequences is that a rough interface can result in localized tensile stresses in the direction perpendicular to the mean interface plane, i.e., σ 22 in our simulation model, which is the primary cause of spalling 5,6 . Figure 8a shows the time-series contour plots of σ 22 for two creep rate coefficients. Note that the maximum tensile stress first occurs near the shoulder of the interface, and the magnitude of this tensile stress increases with roughness. After the rapid roughening of the interface, a pronounced stress redistribution takes place when the oxide  scale becomes thicker, which results in two separate tensile stress concentration regions above the valley and peak of the interface. The largest tensile stress may occur in either of these two regions, depending on the creep rate of the oxide, and these two tensile stress concentration regions can both lead to cracking, as shown in the experimental result presented 6 .
The variation of the maximum tensile stress in the oxide scale is further exhibited in Fig. 8b, which also includes the effect of eigenstrain. Note that the spalling of the oxide scale probably occurs within the first 50 days in these cases. For example, with a large eigenstrain ε eg in = 0.02 and small creep rate B ox , the maximum tensile stress increases from 0.25 to 0.9 GPa in the first 20 days, which could easily result in premature cracking along the mean interface plane. The failure stress of a brittle material is not a constant, as it depends on the size, geometry, and quantity of flaws. To assess the likelihood of fracture, the statistical approach proposed by Beremin et al. 28 can be used.
In Beremin's model, the cumulative fracture probability of a component is estimated based on weakest link statistics and Weibull distribution: where m is the Weibull exponent describing the scattering of the strength of a brittle material; σ u is the characteristic tensile strength of the material associated with a reference volume V 0 ; σ w  Schematics of zirconium oxidation. Schematics of zirconium oxidation, which starts with surface adsorption (reaction formula (i) in Fig. 9). The adsorbed oxygen atoms are then ionized by capturing the free electrons (reaction formula (ii) in Fig. 9), meanwhile the easily oxidized metal atoms are converted to ions (reaction formula (iii). Consequently, counter diffusion leads to both the oxygen anion and metal cation diffuse into the oxide scale and then the oxidation takes place (reaction formula (iv) in Fig. 9).
is the Weibull stress determined from numerical simulation results, which is given by: In Eq. (2), V j is the volume associated with the jth integration point experiencing σ 22 > 0. Following Vermaak et al. 6 , we set m = 10, σ u = 100 MPa, and V 0 = 0.1 μm 3 for zirconia. Figure 8c shows the variations of the cumulative fracture probability P R for the cases studied in Fig. 8b. Note that if the eigenstrain is large (ε eg in ¼ 0:02) and the creep rate is small (B ox ), the fracture probability quickly reaches 100%. If the eigenstrain decreases to 0.01 and the creep rate coefficient is still B ox , the fracture probability first increases due to rapid roughening of the oxide-metal interface, then decreases due to stress redistribution, and finally rises continuously to a magnitude even higher than the maximum magnitude of the first peak, due to thickening of the oxide scale. This indicates that the initial roughening and later thickening can both lead to fracture of the oxide scale. When the creep rate increases to 10B ox , the impact of creep on the fracture probability is significant. In this case, rapid creep can lead to considerable stress relaxation and even result in near-elimination of fracture probability. This indicates that fracture probability is sensitive to the creep rate of oxide scale. Thus, an elevation in the creep rate of oxide scale is highly beneficial in terms of providing increased resistance to cracking.

DISCUSSION
A PF model is proposed to study the stress-oxidation interaction in the zirconia-zirconium system, which involves volumetric eigenstrain, compositional mismatch strain, and creep. The mechanical, chemical, and interface energy terms are incorporated into the Helmholtz free energy expression. To allow for a detailed balance of the reaction system, the oxidation rate is expressed as the difference of forward and backward reaction rates, which are the Arrhenius functions of the chemical potentials of reactants and products, respectively. This rate function is recast in the form of a generalized Allen-Cahn equation, which in conjunction with the reaction-diffusion equation and mechanical equilibrium equation constitutes the PF model that describes the interplay among diffusion, oxidation, and stress variation.
We use the PF model to study a zirconia-zirconium bilayer with the initial thickness of zirconia. This 2D numerical model enables a systematic investigation of how interface roughness affects the growth rate and leads to cracking, which has been observed in experimental studies. The main points derived from the simulation results are as follows.
(1) The growth stress originates from the eigenstrain associated with the volumetric change from zirconium to zirconia, the compositional mismatch strain induced by oxygen diffusion, and the stress-dependent characteristic of diffusivity. Greater stress at the interface leads to a lower oxidation rate, resulting in increased interfacial roughness. However, this impact is limited by the growth of neighboring humps. In addition, creep leads to stress relaxation and interface flattening.
(2) Oxide growth is also significantly influenced by the initial interface topology. Our simulation results agree well with Platt et al.'s finding 7 that a rougher interface leads to a higher growth rate. They also agree with the observed phenomenon that oxide scales with different initial interfacial roughness evolve toward a similar roughness level 7 . (3) The initial rapid roughening leads to two separate tensile stress concentration regions above the valley and peak of the interface, which is consistent with the cracked regions observed experimentally by Tejland et al. 5 and Vermaak et al. 6 . Due to this validation, we incorporate Beremin et al.'s statistical model 28 into the PF model to assess the fracture probability of oxide scale during its growth. In general, a higher creep rate reduces the probability of fracture, because creep leads to stress relaxation and reduces the tendency for interface roughening.
These numerical findings indicate that the durability of oxide scale can be significantly improved by the variation in the initial surface roughness and creep properties of zirconium alloy which can prolong the service life of zirconium-based alloy cladding used in the nuclear power plant.

METHODS
We consider the scenario in which zirconium coated with a thin protective oxide film (zirconia) is exposed to a high-temperature, oxygen-enriched environment over an extended period (hundreds of days). The process starts from the adsorption of oxygen atoms (reaction formula (i) in Fig. 9) on the oxide surface 29 . After capturing free electrons, the oxygen atoms are ionized at the oxide-gas interface (reaction formula (ii) in Fig. 9) and diffuse inward from the oxide surface to the oxide-metal interface; at the same time, the zirconium atoms are converted to cations at the oxide-metal interface and diffuse outward (reaction formula (iii) in Fig.  9) 30 . Their combination leads to the growth of zirconia (reaction formula (iv) in Fig. 9). Due to the difference in the diffusion rates of the reactants, the combination reaction takes place either at the oxide-gas interface (outward oxide growth) or at the oxide-metal interface (inward oxide growth). In the high-temperature oxidation of zirconium, the growth of the oxide scale is dominated by the inward diffusion of oxygen anions rather than the outward diffusion of zirconium cations 2 . Thus, inward oxide growth is considered in this paper, with the simplification that the zirconium cations are not diffusible. As the growth of zirconia is rather slow, the rapid adsorption and ionization processes are negligible 12 . We do not consider the details of the adsorption and sub-reactions, i.e., reaction formula (i-iv) in Fig. 9, but adopt the schematic equation of oxidation reaction by Krishnamurthy and Srolovitz 12 as: Zr þ 2O ! ZrO 2 ; (3)

Helmholtz free energy
The total Helmholtz free energy of the oxide-metal system can be expressed as: where F is the total Helmholtz free energy of the oxide-metal system in the domain Ω; f is the free energy density; and dω is the volume infinitesimal element in Ω. In our model, the free energy density, f, is the summation of the chemical potential energy density, f chem , the interface energy density, f inter , and the mechanical energy density, f mech . The electric potential energy density is ignored in our formula because the oxide is assumed to be charge neutral during a long-term oxidation with the oxide thickness larger than 1 μm following Krishnamurthy and Srolovitz 12 . If the oxide scale is extremely thin, typically <20 nm (called the Debye length 31 ), the charge neutrality assumption is invalid. In this scenario, the electric field would be large enough to dominate the ionic transport through oxide scale, which may be described using Wagner's theory 32 or Cabrera and Mott's theory 30 . The chemical and interface energy densities are functions of a set of elemental concentrations expressed as c ¼ c ZrO2 ; c O ; c Zr À Á , in which the superscripts "ZrO 2 ," "O", and "Zr" represent zirconia, oxygen, and zirconium, respectively. The mechanical energy density depends not only on concentrations, c, but also on the displacement field, d. The concentrations can be further expressed in dimensionless form, indicates the reference concentration of species *, which is generally the saturation concentration of species * in a media (see Supplementary Note 2 for the determination of c Ã ref ). In the simulation in this paper, the dimensionless concentration, c ZrO2 , is unity in the oxide and zero in the metal. Therefore, C. Lin et al.
it can be used as the order parameter, ϕ, to differentiate the oxide and metal phases, namely ϕ ¼ c ZrO2 .
Summing the contributions of diffusible species and the oxide sublattice in the system, the chemical potential energy density can be expressed as: where the first term on the right-hand side, W g (ϕ), with g ϕ ð Þ ¼ ϕ ð Þ 2 1 À ϕ ð Þ 2 , is a double-well function that ensures that both the oxide (ϕ = 1) and metal phases (ϕ = 0) are stable; W is the height of the energy barrier, and is related to the interface energy (per unit area), S, and the interface thickness, ζ, in the form of W ¼ 18S=ζ 19 ; and f O chem and f Zr chem are the chemical potential energy densities of the oxygen and zirconium, respectively. Following Loeffel and Anand 15 , the chemical potential energy densities can be expressed as: where R is the ideal gas constant; T is the thermodynamic temperature; and μ Ã 0 is the standard chemical potential. Based on the diffuse-interface approach, the interface energy density can be expressed as: where ∇ is the gradient operator and λ the scale factor of the interface energy density. If the interface is isotropic, with interface energy (per unit area) S and interface thickness ζ, λ ¼ Sζ 19 . The mechanical energy density is expressed as: where D e is the stiffness matrix for the oxide-metal system and ε e is the elastic strain tensor. The stiffness matrix is expressed as Þ D e met , where the subscripts "ox" and "met" represent the oxide and metal phases, respectively. p(ϕ) is an interpolation function that mollifies the material discontinuity across the metal-oxide interface, generally taking the form p ϕ ð Þ ¼ ϕ 3 10 À 15ϕ þ 6ϕ 2 À Á 33 . As mentioned earlier, mechanical deformation is caused by volumetric eigenstrain, compositional mismatch strain, and creep. Under the assumptions of isotropy and small deformation, the total strain in the zirconia-zirconium system can be expressed as: where ε, ε eg , ε mis , and ε cr are the total strain, volumetric eigenstrain, compositional mismatch strain, and creep strain, respectively. The assumption of small deformation allows the following geometric relation between the total strain, ε, and displacement, d: The volumetric eigenstrain, ε eg , is the result of the volume change from metal to oxide, which is expressed as: where ε eg 11 ; ε eg 22 ; ε eg 33 are the components of the eigenstrain tensor, which are related to the Pilling-Bedworth ratio (PBR). p(ϕ) here ensures that volumetric eigenstrain only occurs in the oxide phase. If the volumetric change is isotropic, the components of the eigenstrain tensor can be expressed as: In most cases of metallic oxidation, the PBRs are greater than unity, resulting in local expansion and compressive stresses. Note that PBR may also be less than unity, for example in the oxidation of K, Mg, or Na. In these cases, the eigenstrains are negative, which leads to a fragmented oxide scale 32 .
During diffusion, the increasing fraction of diffusive element dissolved interstitially in the host lattice is accompanied by a gradual increase in lattice parameter. The continuous change in lattice parameter with composition is referred to as stoichiometric expansion or chemical expansion 11 . In analogy to the thermal expansion coefficient, in numerical studies, the chemical expansion ε mis can be expressed as: where Δx Ã Ã ¼ OorZr ð Þis the deviation of the molar fraction of specie * from its stoichiometric composition; η is the coefficient of compositional expansion (CCE), also known as the chemical expansion coefficient, which can be experimentally determined using thermogravimetric analysis (TGA) under a controlled atmosphere 24 . In our oxide-metal PF model, η is expressed as The deviation of molar fraction Δx* can be expressed in the form of the concentration of specie * and the molar volume V for the oxide-metal system, as follows: where the molar volume can be expressed as Þ V met , in which V ox and V met are the molar volumes of zirconia and zirconium substrate, respectively; c Ã sto is the concentration of diffusible specie * when the oxide-metal system has a stoichiometric composition expressed as c Ã sto ¼ p ϕ ð Þc Ã sto;ox þ 1 À p ϕ ð Þ ð Þ c Ã sto;met . As the zirconium are assumed not to be diffusible, the compositional change induced by zirconium diffusion is omitted. In addition, it is apparent that the stoichiometric concentration of diffusible oxygen should be zero in both zirconia and zirconium substrate, i.e., c O sto ¼ 0. Thus, Eq. (14) is reduced to: The incremental creep strain dε cr can be expressed as: where ε cr eq is the equivalent creep strain; σ eq is the equivalent stress; and σ is the stress tensor, expressed as: Following the Norton power law, the equivalent creep strain rate can be expressed as where σ eq ref is the reference equivalent stress; Þ B met is the creep rate coefficient; and n ¼ p ϕ ð Þn ox þ 1 À p ϕ ð Þ ð Þ n met is the stress exponent.

Oxidation kinetics
A generalized reaction expressed as: describes a one-way reaction from the reactants (R i ) to the products (P j ), with the lowercase letters n i and m j representing the stoichiometric coefficients. According to Birks et al.'s definition 32 , the reaction rate, r, with the unit mol/m 3 /s, can be defined as: where c Ri and c Pj denote the molar concentrations of the reactant R i and product P j , respectively. Note that the rate of a reaction is always positive. A negative sign in Eq. (21) indicates that the concentrations of reactants are decreasing.
In a system containing both reactants and products, the forward and backward reactions take place simultaneously. If the forward reaction proceeds at the same rate as the backward reaction, there will be no net changes in the concentrations of reactants and products, which is called a chemical equilibrium state or detailed balance. However, if the forward or backward reaction is more favorable, there is a net increase in the concentrations of products or reactants. In this case, the net reaction rate is the difference between the forward and backward reaction rates, expressed as: where r 1→2 and r 2→1 represent the forward and backward reaction rates, respectively; k 0 is the rate coefficient, with the unit of mol/m 3 /s; μ 1 and μ 2 are the total chemical potentials of the reactants and the products, respectively; and μ TS is the chemical potential at the transition state.
The total chemical potentials for the reactants and the products can be expressed as: where μ Ã Ã ¼ R i or P j À Á is the chemical potentials for the specie *, which is defined as variational derivatives of the total free energy, F , with respect to the concentration, c*, expressed as: Following Bazant 34 , μ * can be further divided into two parts: where a * is the activity of specie * that is only concentration dependent, defined as: and μ ex Ã is the excess chemical potential originating from the standard chemical potential, mechanical energy, electric potential, and so on, expressed as: By substituting Eq. (26) into Eqs. (23) and (24), the chemical potentials for the reactants and the products can, respectively, be recast as: and where a 1 and a 2 are the sums of the activities of the reactants and products, respectively, and μ ex 1 and μ ex 2 are the total excess chemical potentials for the reactants and products, respectively.
The excess chemical potential in the transition state, μ ex TS , is defined as 34 : where a TS is the sum of the activities of the reactants and products; ρ is the asymmetry parameter, which is approximately a constant between zero and one for numerous reactions 35 ; and μ ex TS is the excess chemical potential at the transition state, which is defined as the linear combination of μ ex 1 and μ ex 2 with ρ. Substituting Eqs. (29), (31), and (33) into Eq. (22), the generalized reaction rate can be expressed as: r ¼ k 0 a TS a 1 exp ð1 À ρÞ μ ex 1 À μ ex 2 À Á RT À a 2 exp À ρ μ ex 1 À μ ex 2 À Á RT : For the zirconium oxidation given by Eq. (3), the total chemical potentials of the reactants and products can be derived as: and respectively, where the excess chemical potentials of the reactants μ ex 1 and products μ ex 2 are expressed as: where σ 11 , σ 22 , and σ 33 are the normal stresses; tr σ ð Þ ¼ σ 11 þ σ 22 þ σ 33 is the trace of the stress tensor; μ mis is the chemical potential resulting from the mismatch strain; μ eg is the chemical potential resulting from the eigenstrain; and μ stif is the chemical potential resulting from the change in elastic properties during oxidation.
Substituting Eqs. (36)-(41) into Eq. (35), the reaction rate for the zirconium oxidation can be expressed as: where μ int ¼ V ox W∂g ϕ ð Þ=∂ϕ À λ∇ 2 ϕ À Á is the chemical potential resulting from the interface energy. As we assume that zirconium atoms are not diffusible, the concentration of zirconium in the metallic phase is set at a constant c Zr ref . Thus, the dimensionless concentration of zirconium c Zr is set to be unity. In addition, following Chen et al. 36 , we assume that the chemical potential resulting from the interface energy is much smaller than RT. The approximation r % rj μ int ¼0 þ ∂r=∂μ int ð Þ μ int ¼0 h i μ int can be adopted to recast the reaction rate as follows: Governing equations for oxidation, diffusion, and mechanical deformation The oxidation rate r indicates the rate of oxide phase formation, i.e., r ¼ ∂c ZrO2 =∂t ¼ 1=V ox ð Þ∂ϕ=∂t. Therefore, the governing equation for the order parameter ϕ can be expressed as: where the function ∂p ϕ ð Þ=∂ϕ is used to ensure that the oxidation reaction only occurs at the oxide-metal interface 0 < ϕ < 1 ð Þ 23 ; L S and L κ are the coefficients to scale the contributions of the interface energy and the oxidation kinetics to phase migration, respectively; and Λ is the dimensionless reaction driving force. These variables are expressed as: C. Lin et al.