Modeling of Zirconium Atom Redistribution and Phase Transformation Coupling Behaviors in U-10Zr-Based Helical Cruciform Fuel Rods under Irradiation

: Uranium–zirconium metal-based Helical Cruciform Fuels (HCFs) have shown a promising prospect for their use in advanced nuclear reactors. However, during irradiation, dual-phase coexistence and the spatial heterogeneous distribution of zirconium atoms occur at higher powers, a ﬀ ecting the thermo-mechanical coupling behaviors and safety of fuel elements and assemblies. In this study, based on the phase-ﬁ eld approach, the coupled multi-ﬁ eld governing equations to describe the zirconium di ﬀ usion and phase evolution for U-Zr metallic fuels are improved. Furthermore, the corresponding numerical algorithms and procedures for multi-ﬁ eld coupling calculations are developed. The numerical predictions of zirconium atom fraction are in good agreement with the relevant experimental results, validating the developed models, algorithms and programs. The zirconium atom redistribution and phase transformation coupling behaviors in high-power U-10wt%Zr-based HCF rods are also obtained. Moreover, the complex evolution mechanisms of multi-ﬁ eld variables are analyzed. The results indicate the following: (1) the irradiation enhancement of the thermal mobility and chemical mobility plays a critical role in the redistribution of Zr atoms; (2) the multi-ﬁ eld results of HCF rods have helical symmetric characteristics; (3) the contribution competitions of the temperature gradient and chemical potential gradient within the α phase and γ phase signi ﬁ cantly in ﬂ uence the zirconium-atom redistribution, with the zirconium-rich zones formed in the elbow region and the zirconium-poor zones appearing inside. These research e ﬀ orts supply a foundation for the further involvement of mechanical ﬁ elds in multi-ﬁ eld coupling computation.


Introduction
To enhance the economic benefit of nuclear energy, it is imperative to achieve power updates or extend the refueling cycle in Pressurized Water Reactors [1].An advanced U-Zr metal-based helical cruciform fuel (HCF) has been developed recently by Lightbridge [2][3][4].In the HCF rod, the U-Zr metallic fuel core is metallurgically bonded with the zirconium alloy cladding.In comparison to oxide fuels like UO2, the U-Zr metal fuels generally exhibit superior thermal conductivity and facilitate easier fuel recovery [5,6].Furthermore, Zr-containing uranium alloys demonstrate excellent corrosion resistance and dimensional stability during thermal cycling [7].These inherent advantages position uranium-zirconium alloys as an exceptional choice for nuclear reactor fuels.In fact, the concept of HCF was first proposed by the Soviet Union [8,9].Compared to the traditional cylindrical fuels, the HCF rods possess a higher surface-to-volume ratio, enabling increased heat output and promoting enhanced cross-mixing characteristics among the liquid channels.Additionally, each HCF rod contacts the surrounding rods at multiple Citation: Chen, X.; Xie, Z.; Mao, X.; Ding, S. Modeling of Zirconium Atom Redistribution and Phase Transformation Coupling Behaviors in U-10Zr-Based Helical Cruciform Fuel Rods under Irradiation.Metals planes along the assembly length, eliminating the need for spacer grids due to their selfspacing traits [10].As a result, the U-Zr metal-based HCFs have a promising application prospect in future advanced nuclear reactors.
Before the actual application, the optimal design should be performed for HCFs, including the power density, the constituents of U-Zr metallic fuels and the rod geometrical structures.It should be mentioned that during irradiation the dual-phase coexistence and spatial distribution inhomogeneity of U-Zr metal fuels will occur at higher powers [11].Under the high-temperature conditions above 890 K, phase transformation will take place, and subsequently the zirconium atoms will be redistributed and the phase compositions will evolve with the increase in burnup [12].The diffusion coefficients of zirconium in different phases differ from each other [12], and the phase compositions depend on the temperature and zirconium concentration [13].Consequently, the coupling behaviors of zirconium atom diffusion, phase distribution and evolution are maintained during the inpile service.Thermo-mechanical coupling behaviors of fuel elements have received much attention [14][15][16][17].Due to the facts that the thermal conductivities, the creep and irradiation swelling performances are heavily affected by the constituents of U-Zr fuels [18][19][20][21][22], it can be understood that the thermo-mechanical coupling behaviors of HCF rod and assembly will be simultaneously impacted.The heterogeneities in the distribution and evolution of the zirconium atom concentration and phase constituent will pose a challenge to fuel safety.In order to realize the optimal design and efficiently control the operating conditions, it is critical to develop the theoretical models and numerical simulation methods for predicting the complex above-mentioned coupling behaviors.
So far, the modeling of the comprehensive coupling behaviors has not been reported, which pertains to the full coupling of zirconium atom diffusion, the phase distribution and evolution, the temperature field and the displacement/stress/strain fields.As for the zirconium redistribution kinetics, researchers have conducted some studies with a series of thermodynamic models on cylindrical fuel rods [12,23,24].These studies supply the basis for the modeling of the behaviors in U-Zr-based HCF rods because of the consistent thermodynamic properties of U-Zr alloys (e.g., phase composition) in both cylindrical fuel rods and HCF rods.Hofman et al. [12] used a diffusion model described directly with the gradients of Zr atom concentration and temperature, obtaining the radial concentration distribution within the U-10Zr fuel slug using a difference method.Cheol et al. [23] carried out calculations coupling the burnup-dependent temperature field and Zr atom redistribution field, with the effects of sodium intrusion into the porous fuel slug considered.Kim et al. [24] modeled the constituent redistribution in the U-Pu-Zr fuels, by simplifying the phase diagram of the U-Pu-Zr ternary system.It should be mentioned that in these previous studies, the used Zr atomic diffusion model was explicitly expressed with the gradients of Zr atom concentration and temperature; the model parameters were dynamically related to the phase distribution, which was assumed to be in the equilibrium state and depend on the current concentration and temperature.Due to the fact that the phase fractions could not be solved, the respective phase contributions could not be precisely reflected.The phase-field method has been well adapted to deal with interface problems and thermodynamic issues [25][26][27][28][29].In order to further improve the modeling ability, Hirschhorn et al. [30] developed the multi-field coupling models for U-10Zr alloys using the phase-field method, and the evolution results of phase transformation and Zr redistribution were acquired with the finite element method.Subsequently, the 1D simulation was extended to the 2D simulation for U-Pu-Zr fuels with BISON [31].Nevertheless, the Zr atom interdiffusion coefficients were generally scaled up by a certain number of times to match the experimental data.Wen et al. [32] perceived that the strong irradiation effects caused by the fission fragments should be responsible for this enhancement.So far, a quantitative relationship of the Zr atom diffusion coefficient with the fuel fission rate is absent.Moreover, a 3D numerical simulation of U-Zr redistribution behavior still needs to be developed.To assess the impacts of fuel fission rate on the comprehensive multi-field coupling behaviors in U-Zr-based HCFs, it becomes imperative to implement 3D simulations, with the fission-rate related multi-fields well modeled and fully coupled.
In this study, the irradiation enhancements are involved in the chemical mobility and thermal mobility models of the three U-Zr alloy phases.The mobility parameters are correlated directly with the fuel fission rate.The three-dimensional coupled theoretical models were presented for multi-fields, including the fields of the temperature, the zirconium atom concentration, the γ phase fraction, the zirconium atom fractions in different phases and the chemical potential.The equivalent integral weak forms for the corresponding differential governing equations are derived.Simultaneously, the fully coupled computational algorithms are developed.The iterative computation method of the nonlinear finite element equations for the multi-field variables are presented.The developed models, algorithms and procedures are validated with the experimental data.The coupling behaviors of zirconium atom diffusion and phase transformation in the U-10Zr-based HCF rods under irradiation are numerically obtained, and the effects of fuel fission rate enhancement are also investigated.This study lays a foundation for the further involvement of the irradiation-induced mechanical fields.

Multi-Field Coupling Models for Zirconium Atom Diffusion and Phase Transformations in U-10Zr Fuels
The geometry of the considered HCF rod in this study is similar to that in Ref. [33], as illustrated in Figure 1, with the core of U-10Zr perfectly bonded with the cladding of Zr-4 alloy.The dimensions are also shown in Figure 1.Compared to the traditional fuel rod, the external surface area of the HCF rod is amplified, facilitating heat exchange with the coolant.The outer surface of the fuel rod is denoted as F S , and the interface between the fuel core and the cladding is depicted as B S in Figure 1.
As mentioned in Section 1, higher temperatures of U-10Zr fuels will result in phase transformation and Zr redistribution.In this section, the theoretical models (differential governing equations) for the involved multi-fields are first integrated, in which the effects of the fuel fission rate are quantitatively described.

Temperature-Field Model
The differential governing equation for the temperature field is expressed as follows [34]: The boundary conditions are denoted as follows [35]: where ρ is the density; p c is the specific heat capacity; k is the thermal conductivity; q  is a heat generation rate; n  is the external normal of the external surface; w T demon- strates the coolant temperature; h denotes the heat transfer coefficient.It is noted that the temperature and Zr atom concentration dependence of the thermal conductivity for U-Zr fuels [18] is considered in this study.

Phase-Field Related Models
It has been mentioned in Section 1 that in some of the previous models for U-Zr component redistribution [12,23,24], the phase composition and evolution are determined according to the steady-state phase diagram.As shown in Figure 2, the phase composition can be obtained directly with the values of temperature and Zr atom concentration.However, the phase fractions in the binary regions are unable to be calculated.Moreover, the non-equilibrium phase evolution possibly occurs during in-service irradiation.Based on the thermodynamic theories, a phase-field model is used to describe precisely the dominant thermodynamic process, possessing the ability of considering the interactions of different parts.The U-Zr phase diagram in Figure 2 is plotted in this study without distinguishing the phases of 1 γ and 2 γ [12], regarded as a unified phase of γ.It can be understood that the phase compositions of U-Zr fuels depend on the temperature and Zr atom concentration, possibly occupied by the four single phase regions of γ, α, β, δ and the three binary phase regions.It is noted that the phase of α + δ will be modeled with the parameters of phase α [12].As simplified in Ref. [30], it is also assumed that only three phases exist in the whole system, implying 1 α β γ + + = , with γ, α and β representing the phase fractions of the three phases; the phase α and phase β are considered as one phase independent of phase γ.Then, the ternary phase system can be transformed into a binary phase system for calculation.Therefore, only the governing equations for the independent phase fraction γ should be developed.
The evolution equation of γ is expressed with the Allen-Cahn equation [30,36], namely with the boundary conditions denoted as where F δ δγ denotes the variational ration of the total free energy to the phase fraction of γ; L is the kinetic coefficient of phase boundary mobility to govern the evolution rate of the multi-phase region, described as where 0 L is a constant; T γ is the phase transition temperature of γ phase.It should be mentioned that Equation ( 5) is newly developed in this study, with which the phase transition could be precisely predicted during the process of temperature increase.
In fact, the phases of α and β do not co-exist, as displayed in Figure 2. A conversion function of h αβ is proposed to make a smooth transition between two phases, described as where , depicting the phase transition temperature from phase α to phase β.
Consequently, the phase fraction of phase α can be obtained as ( )( ) phase β being obtained as ( ) The total free energy of the system includes the contributions of both the bulk free energy and the gradient free energy, given as [30] ( ) where the bulk free energy density bulk f consists of the contributions from the γ phase, α phase and β phase, and the gradient free energy density grad f results from the concentration gradient and the phase fraction gradient.The bulk free energy density is derived from the CALPHAD free energy, expressed as [37] ( ) where is the reference molar volume [30]; ( ) g γ γ = − ; f αβ denotes the bulk free energy density function for the phases of α and β [32]; f γ is the bulk free energy density for the γ phase [32].Similar as that in Ref. [32], the gradient free energy density is defined as where γ κ is the energy coefficient of phase fraction gradient; c κ depicts the energy co- efficient of Zr atom concentration gradient c ∇ , needing to be identified.The parameters of γ κ and w characterize the size and shape of the diffusion inter- face, related to the interface width l and the interface energy σ [30]: In order to obtain the explicit expression of Equation ( 3), the variation in the total free energy in Equation ( 3) is firstly derived as ) Combined with the divergence theorem, integration by parts for the third item of the right-hand yields ( ) Simultaneously, Equation ( 12) can be written as Combining Equations ( 11)- (13) with Equation (4), we have It can be understood that the evolution equation of Equation (3) for the phase fraction of γ is a nonlinear equation, coupled with the Zr atom concentration of c .Furthermore, the free energy density of f αβ or f γ is expressed as a function of c αβ or c γ , with c αβ and c γ denoting the equilibrium atomic fractions of Zr in the respective phases for the region of α (or β) + γ.The new variables of c αβ and c γ should be further correlated to the Zr atom concentration of c and the phase fraction of γ .Generally, Kim-Kim-Suzuki (KKS) phase-field formulations are adopted as [37] ( ) with Equation (17) representing the required the equilibrium condition for c αβ and c γ at the phase boundaries.Here, Equations ( 16) and ( 17) should be satisfied in every point of the binary phase region.Actually, it can be understood that c αβ or c γ will be equal to c for the single-phase region.As a result, the fields of γ , c αβ , c γ and c are coupled with each other.

Zr Atom Concentration Field Related Models
The Zr atom concentration can be characterized with the atomic fraction.In this study, the independent variable c is chosen to denote the atomic fraction of Zr, meaning that the atomic fraction of U is obtained as 1 c − .The Zr atom concentration satisfies the mass conservation equation, with the diffusion essentially driven by the chemical potential gradient and temperature gradient, governed by the modified Cahn-Hilliard equation [38]: The interface conditions are expressed as where J denotes the Zr atom flux; c M is the chemical mobility; is the chemical potential of μ , depicting the variation ratio of the total free energy to the Zr atom concentration, similarly obtained as It is noted that the chemical potential of μ could not be explicitly described with the concentration c , due to the fact that the bulk free energy density is directly related to the variables of c αβ and c γ .The chemical potential also constitutes a field.If the chemical potential of μ is regarded as a field variable, the Zr atom concentration field is directly coupled with the fields of temperature T and μ .Simultaneously, the field of μ will directly depend on γ , c αβ , c γ and c .By bringing Equation ( 8) in and combining with the KKS phase-field model, the exact coupling relationship of the chemical potential can be derived as ( ) The chemical mobility in Equations ( 18) and ( 19) is considered as the homogenized contribution of the three phases, given as [30] ( )( ) ( )( ) In the phases of α and β, the chemical mobility is expressed as follows [30]: where i Q is the activation energy; R is the ideal gas constant; i D  is the diffusion coefficient.During irradiation, fission fragments enhance diffusion in the fuel medium.Similar to the enhancement of vacancy diffusion coefficient, the diffusion coefficient of i D  will be enhanced.This irradiation effect should be taken into account, especially at lower temperatures.In accordance with the treatment of the irradiation-enhanced effect of atomic migration in Ref. [39], the diffusion coefficient of i D  is correlated with the fuel fission rate, expressed as where f  is the fuel fission rate; In this study, the irradiation enhancement effect on the chemical mobility of the γ phase is also involved as where U β and Zr β are the thermal atomic mobilities of U and Zr, obtained as [40] 0 exp U,Zr where i H is the activation enthalpy and 0 β is the atomic mobility coefficient.
The thermal mobility of T M has a similar expression to the chemical mobility, de- picted as [30] ( )( ) ( )( ) where T M γ are the thermal mobilities of the corresponding phase.For the phases of α and β, the thermal mobility is obtained as [30] * exp , where * i Q  is the heat of transport of the corresponding phase.It should be mentioned that Zr atoms exhibit migration in both phase α and phase β against temperature gradients [12].When considering this in conjunction with Equation ( 18), it becomes imperative to incorporate a negative sign before Equation (28).
For the γ phase, based on the model in Ref. [30], the irradiation-enhanced thermal mobility model is developed as where * U Q  and * Zr Q  are the heat of transport of U and Zr atoms; v A γ is the irradiation enhancement coefficient.

Numerical Implementation of Multi-Field Governing Equations
The multi-field governing equations in Section 2 are numerically solved with the developed procedures in the platform of MOOSE (mooseframework.inl.gov), a multi-physics field coupling software based on the finite element method [41].Six field variables were considered in direct-coupling simulations: temperature T , c αβ , c γ , phase fraction γ , chemical potential μ and Zr atomic fraction c , implying that they are the nodal varia- bles of the finite elements.The application of the finite element method converts governing equations into a system of nonlinear equations: where the nonlinear system of equations of { } =0 u R depicts the equations that should be satisfied by the to-be-solved nodal variables, obtained by the weak forms of equivalent integral governing equations for the variable u.It is noted that { } u R represents an array of 1 N × , with N denoting the total node number; the expression for its component of ) will be derived in the following.
Owing to the nonlinearity and multi-field coupling characteristic of differential governing equations in this study, the arrays of { } u R depend on several correlated nodal variables.Consequently, the Newton iterative method is used to linearize this system of nonlinear equations to algebraic equations, with the iterative scheme expressed as where { } where ); u and v represent the corresponding two variables, and the derived ex- pressions will be obtained as follows.So, the array of { } R is defined as the residuals after the previous iterations, and the iterative convergence will be achieved when the array norm tends to be zero.It should be mentioned that the coercive boundary conditions should be introduced in order to eliminate the singularity of [ ] 6 6 m N N J × .

Residuals and Jacobian Matrix Components Related to the Temperature Field
The differential governing equations for the temperature field are shown in Equations (1) and (2).The equivalent integral form is derived as Combined with the divergence theorem and integration by parts, the equivalent integral weak form is obtained as With the finite element method, the domain is discretized as a sum of elements.Within a typical element, the test function can be described as where i T is the temperature value of the ith node, i T N is the corresponding shape function.It is noted that the dummy index of i is used here to indicate the summation of all items from 1 to N. The Crank-Nicolson method is used as the time integration scheme in this work.The Crank-Nicolson time integration scheme is a second-order implicit method and it is unconditionally stable [42].As a result, Equation ( 33) is rewritten as where t u is the value at time t , u without superscript represents the value at time N  depicts the shape functions for the corresponding surface element.So, we have The Jacobian matrix components are the partial derivatives of the residual terms for each coupled nodal variable.The governing equation for the temperature field is uncoupled with the other variables, so only the Jacobian matrix components related to the node temperatures are non-zero, obtained as

Residuals and Jacobian Matrix Components Related to the Phase Fraction Field
The differential governing equations for the phase fraction of γ are shown in Equa- tions (3) and (4).The corresponding equivalent integral form is expressed as where, the relations in Section 2.2 have been comprehensively considered and involved.
Combined with the divergence theorem and integration by parts, the equivalent integral weak form is obtained as Using the similar finite element interpolation method and Crank-Nicolson time integration scheme, we have ( ) ( ) As analyzed in Section 2, Equation ( 40) is coupled to the variables c αβ , c γ and T .
Therefore, non-zero off-diagonal Jacobian matrix components are provided as ( ) ( )

Residuals and Jacobian Matrix Components Related to the Equilibrium Atomic Fraction Field
The finite element method is also used to solve the KKS model, expressed as Equations ( 16) and (17), with the based equivalent integral form given as ( ) Using the similar finite element interpolation method, we have ( ) Because the KKS models are coupled with the variables c , γ and T , the compo- nents of the diagonal and off-diagonal Jacobian submatrices correlated to Equation (48) are derived as ) Similarly, the components of the diagonal and off-diagonal Jacobian submatrices correlated to Equation (49) are obtained as

Residuals and Jacobian Matrix Components Related to the Zr Atomic Fraction Field and Chemical Potential Field
The evolution of the Zr atomic fraction follows the Cahn-Hilliard equation, where the chemical potential μ is defined as Equation (21).In order to improve computational efficiency, the variables corresponding to these two governing equations were exchanged.Namely, Equation ( 18) is regarded as the equation that the Zr atomic fraction should satisfy.Then, combined with the boundary condition, the equivalent integral form is expressed as ( ) ( ) Integration by parts for the second item, the equivalent integral weak form is derived as ( ) Using the similar finite element interpolation method, we have It is noted that the nonlinear equation of Equation ( 59) is coupled with the variables c αβ , T , μ .Consequently, the components of the diagonal and off-diagonal Jacobian submatrices are supplied as As a result, the equivalent integral form for Equations ( 18) and ( 19) is similarly obtained as Integration by parts for the second and third items, the equivalent integral weak form is derived out as Using the similar finite element interpolation method, we have Due to the fact that Equation ( 65) is coupled with the variables of c and T , the components of the diagonal and off-diagonal Jacobian submatrices are described as It should be mentioned that the values of T in Equations ( 65) and ( 69) for the integration points are also calculated as

Validation of Multi-Field Coupling Models/Algorithms/Procedures
In order to validate the effectiveness of the multi-field models, algorithms and procedures, the coupling behaviors are simulated with the developed finite element procedures, for the multi-fields of T , c αβ , c γ , phase fraction γ , chemical potential μ and Zr atomic fraction c in the DP-81 fuels irradiated in an EBR-II reactor [12].The predictions will be compared with the relevant experimental data obtained by the researchers in Argonne National Laboratory [11].
Due to the axial-symmetric characteristics of this problem, the procedures are degenerated into a one-dimensional version to perform the simulation of the multi-field coupling behaviors.Essentially, the non-homogeneous temperature field is the dominant driving force in initiating the subsequent coupling behaviors.Consequently, the transient temperature field is simultaneously modeled in this study, which is different from the treatment in some references that use an unvaried radial temperature distribution [32].
It should be mentioned that the U-10Zr fuels were used for the DP-81 pin.The DP-81 pin geometry, power, and surface temperature data are taken from Hofman et al. [12], with an initial radius of 2.17 mm and a linear power of 24 kW/m.The possible outer surface temperatures range from 890 K to 916 K [43], and an outer surface temperature of 900 K is considered in this study, similar to those in some other references [30,32].Due to the fact of Zr atom redistribution during irradiation, the thermal conductivity model related to temperature and Zr atom concentration is used [18].Moreover, the effect of sodium permeation into the formed pores is considered, because the sodium infiltration into the pores will affect the local temperature distributions and the local temperature gradient [43].It is noted that the fuel porosity and the effects of sodium infiltration are involved in the local thermal conductivity [43].Therefore, different factors are multiplied over the calculated values by the model in Ref. [18] for each phase region, bringing the temperature distribution close to that in Ref. [43].In order to reasonably describe the irradiation effects on the redistribution behaviors of U-Zr alloys, the irradiation enhancement factor in Equations ( 24)-( 29) are identified in this study, as listed in Table 1, together with some other model parameters adopted in the simulation of the irradiation-induced coupling behaviors within the DP-81 pin.
The radial microprobe scans experimental data for the DP-81 pin were those at 3.6% FIMA (FIMA is the ratio of the fission atoms to the initial fissile atoms), providing an important reference for validation.As shown in Figure 3, distinct zirconium migration phenomena at 3.6% FIMA can be observed from the scattered black points, compared to the original Zr atom fraction of 10%.Simultaneously, Figure 3 shows the simulation results for the multiple fields.The predictions of Zr atom fractions plotted on the black line are found to be in a good agreement with the experimental results, validating the developed models, algorithms and procedures.It is noted that the dashed line in Figure 3a represents the predicted Zr atom fractions without considering the irradiation enhancement effects in Equations ( 24)-( 29), demonstrating that the irradiation enhancement plays a dominant role in the Zr atom redistribution behaviors.If the irradiation enhancement is not taken into account, the predictions will differ greatly from the experimental data.It can also be observed that the predicted phase distribution is close to those in Ref. [30], with a single phase of γ dominated in the high-temperature region within a distance of ~1.0 mm away from the pin center.In the other regions, the binary phases of β + γ or α + γ appear.In fact, the dynamic evolutions of multiple fields occur.In Figure 3b, the results of the Zr atom fraction at different burnup levels are exhibited.On the whole, the Zr atoms tend to diffuse towards the pin center from the external region.At 3.6% FIMA, a zirconium-poor region is formed in the middle of the radius.It should be mentioned that the temperature gradient induced by the sodium infiltration in this middle region near the interface between γ phase and β + γ phase plays a dominant role in the Zr atom distribution.In this study, the finite element model is firstly established for the whole helical cruciform fuel rod, with the mesh grid figure shown in Figure 4a.The rod is discretized as 24,000 elements.It should be mentioned that the end plugs are not established due to the fact that the multi-field coupling behaviors occur within the fuel meat, different from the geometry in Figure 1.The highest linear power of 80 kW/m given in Ref. [35] is chosen, corresponding to a heat generation rate of 4.56 × 10 9 W/m 3 and a fission rate of 1.41 × 10 20 fission/(m 3 •s).Uniform convective heat transfer boundary conditions are adopted for the side surfaces the fuel rod, with the heat transfer coefficient and the coolant temperature set as 65,900 W/(m 2 •K) [35] and 573 K, respectively.It is known that the helical cruciform fuel (HCF) developed recently by Lightbridge is referred to the U-50Zr metal-based ones.In this study, the multi-field coupling behaviors are investigated for U-10Zr metal-based HCFs.

Multi-Field Simulation Results for the HCF Rod
The simulated multi-field results for the HCF rod are shown in Figure 5.The crosssections of the HCF rod are known to be rotated around the rod axis, and it can be seen from Figure 5a that the results of the temperature and the Zr atom fraction at the side surface of the fuel rod have the helical symmetric distribution characteristics.Simultaneously, three typical helical paths are selected in the regions of rod elbow, blade and center, according to the contour plots.It can be found from Figure 5b that the corresponding results are almost unvaried along the three paths.Furthermore, the contour plots of the Zr atom fraction within several cross-sections at different heights are displayed in Figure 5c, demonstrating that the results within a certain helical cross-section can be regarded as the rotated ones of the other cross-sections.In fact, all the multi-field results in this study represent the helical symmetric distribution characteristics.Consequently, the finite element model can be further simplified to the 1/4 section of a thin slice, without allowing for the helical geometrical pattern.Thus, the adopted element number is heavily reduced, as illustrated in Figure 4b.The detailed results in the following sections are obtained with this simplified finite element model.It is noted that the mesh sensitivity analysis for the simplified finite element model is performed, demonstrating the convergences of multi-field variables.

The Multi-Field Results during the Reactor Startup
This section focuses on the multi-field result analysis for the HCF fuel rod with a linear power of 80 kW/m.It should be mentioned that the evolution process of temperature during the reactor startup is also simulated, corresponding to the solution of the transient temperature field.It is assumed that the coolant temperature increases from room temperature to the specified temperature within 160 min.As shown in Figure 6a, with the increasing of the temperatures on Path 1, the γ phase transition gradually occurs from the path origin towards the end, with the γ phase fraction achieving a value of ~0.4.It can be seen from the phase diagram in Figure 2 that the central region of the fuel meat is occupied by the α + γ phase.For the regions where the phase transition temperature is not achieved, the γ phase fraction is predicted to be zero and only the α phase appears.This heterogeneous phase transition process is successfully simulated in this study, depending on the phase mobility model proposed as Equation (5).It should be mentioned that the precious phase mobility model reflecting the actual phase transition rate should be developed according to the abundant experiment data.
Figure 6b depicts the contour plots of the temperature, the chemical potential, the γ phase fraction and the Zr atom fraction at the end of reactor startup.The highest temperature of ~909 K appears in the center of the fuel rod, and the lowest temperature of ~622 K locates in the fuel rod blade.It can be known that distinct temperature gradients are formed within the fuel meat.Figure 6a denotes that the magnitude of temperature gradient along Path 1 is enhanced near the path end, with the proceeding of temperature increase.Compared to those near the fuel rod elbow, the magnitude of temperature gradients surrounding the fuel rod blade are much lower.It can be also found that the higher values of the γ phase fraction only exist around the central region of the fuel meat.Correspondingly, the chemical potentials in this α + γ phase region are greatly reduced from the initial value of ~2.81 × 10 9 J/m 3 to the current value of ~0.78 × 10 9 J/m 3 .It is noted that near the boundary of the lowest chemical potential, the values of the local chemical potential are observed to increase sharply and then decrease, achieving a highest value of 3.5 × 10 9 J/m 3 .This region corresponds to that with a decrease in the γ phase fraction from ~0.4 to ~0 in Figure 6a.Actually, this is a result of the fact that the equilibrium atomic fractions c αβ decrease due to the presence of the γ phase in the central region.From Equation ( 21), it can be seen that the chemical potential is related to c αβ , and hence a sudden change is formed at the interface of α + γ phase and the pure α phase.Also, the chemical potential is temperature dependent, so in the pure α phase region the chemical potential decreases with decreasing temperature.Consequently, a high magnitude of chemical potential gradient is formed around this transition region.During the reactor startup, the Zr atom fraction is nearly the same as the initial value of ~0.22, due to the considerably short time period.

Redistribution of Zirconium Atom with Burnup and the Dominated Mechanisms
Figure 7 depicts the contour plots of the Zr atom fraction at different burnup levels.During the whole irradiation process, a significant redistribution of Zr atoms occurs.A zirconium-poor region with the Zr atom fraction of ~0 is formed at the α + γ phase boundary at 0.2% FIMA, and then the zirconium-poor region is expanding outward with burnup.At 1% FIMA, a zirconium-rich zone with the Zr atom fraction of ~1 is formed in the elbow region, which can be found to gradually extend along the boundary of the fuel to the blade region from the contour plots at 4% FIMA and 10% FIMA.At 10% FIMA, a zirconium-poor zone also appears around the center of the fuel.These evolution traits stem from the complex influences of multi-fields.Figure 8 displays the multi-field results at different burnup levels on Path 1, including the results of the Zr atom fraction, the temperature, the chemical potential and the γ phase fraction.Combined with Equation ( 18), it can be known that the migration of Zr atoms is mainly due to the combined effects of the chemical potential gradient as well as the temperature gradient.In order to better understand the above redistribution behaviors of Zr atoms, the contributions of the chemical potential gradient term and the temperature gradient term in Equation ( 18) are given in Figure 9.   8a.Simultaneously, it can be seen from Figure 9a that within a short time after the end of reactor startup, a part of the Zr atoms surrounding the α + γ phase boundary migrate towards the fuel center, creating a local Zr atom fraction increase to ~0.44 in Figure 8a.Different from the abnormal diffusion behaviors in the γ phase, in the α phase region Zr atoms will migrate against the temperature gradient, described with Equations ( 19) and (28).Therefore, the Zr atoms in the α phase region can be found from Figure 9a to diffuse away from the fuel center, induced by the temperature gradient there.Under the combined contributions of temperature gradient and chemical potential gradient (as shown in Figure 9a) to the Zr atom flux, the zirconium-poor region is firstly formed outside the α + γ phase boundary, as shown in Figures 7 and 8a,d.It can be found from Figure 8d that from the end of reactor startup to 0.2% FIMA the boundary of α + γ phase moves a little away from the fuel center, and the locally enhanced γ phase fraction can be observed near the α + γ phase boundary, similar to the distribution pattern of Zr atom fraction in Figure 8a.
As shown in Figure 8c,d, the chemical potential gradient around the α + γ phase boundary decreases with the increase in burnup, and the resulting diffusion flux of Zr atoms surrounding the phase boundary weakens.It can be seen that the region with higher magnitudes of chemical potential gradually moves towards the end of Path 1.It can be understood that the superposed contributions of the chemical potential gradient and temperature gradient mainly occur in the α phase region.As shown in Figure 9b, at 0.2% FIMA, the temperature gradient contribution to the diffusion flux of Zr atoms dominates over the whole region, especially near to the fuel elbow.As a result, the Zr atoms in the α phase region still migrate to the outside, leading to the expansion of the zirconium-poor zone, as depicted in Figures 7 and 8a.Due to the short diffusion distance and higher temperature gradient from the fuel center to the elbow region, at 1% FIMA the accumulation of Zr atoms can be seen at the elbow region, forming a zirconium-rich region in Figure 7. Afterwards, the width of the zirconium-poor region on Path 1 essentially stays ~0.5 mm, because the diffusion flux of Zr atoms along this path is severely lowered.As shown in Figure 9c, the contributions of the temperature gradient and chemical potential around Path 1 counteract significantly.Simultaneously, it can be found from Figure 9c that the temperature gradient contributions surrounding the previously formed zirconium-rich region dominate, resulting in the continuous expansion of this region in Figure 7.
Meanwhile, it can be seen from Figure 8a that the Zr atoms inside the α + γ phase region also change with burnup.The distribution characteristics of Zr atoms are more complicated.According to Equations ( 22)-( 29), it can be known that the Zr atom redistribution behaviors in the α + γ phase region depend on the individual contributions of the α phase and γ phase.It can be seen that after 0.2% FIMA the chemical gradient contribution to the Zr atom migration can be neglected in the α + γ phase region, replaced by the dominant contribution of temperature gradient.It should be mentioned that the Zr atoms within the γ phase tend to diffuse along the direction of temperature gradient, and on the contrary the Zr atoms within the α phase will migrate against the temperature gradient.As depicted in Figure 8d, the γ phase fraction on the left side of α + γ phase is locally enlarged around the phase boundary with the α phase on the right, with the maximum value of ~0.55 at 0.2% FIMA.Subsequently, the γ phase fraction varies greatly.Accordingly, the Zr atom fractions have a similar evolution pattern, resulting from the competition of the temperature gradient contributions of α phase and γ phase.As shown in Figure 8a, at 6% FIMA the Zr atom fraction in the α + γ phase region is already smaller than the original value of 0.22.With the decrease in the γ phase fraction in Figure 8d, the temperature gradient contribution from the α phase gradually dominates, resulting in an accelerated outward diffusion of Zr atoms in the α + γ phase region.Meanwhile, as shown in Equations ( 23), ( 26) and ( 28), the migration rate is positively correlated with temperature, so that at 10% FIMA a zirconium-poor region appears around the fuel center.It is noted that the results of γ phase fraction and Zr atom fraction in this zirconium-poor region coincide with the phase diagram at the corresponding temperature values in Figure 8b.It should be mentioned that the adopted thermal conductivity model in this study is dependent on the Zr atom fraction, causing the variations in the temperature field as a function of burnup.
As emphasized previously, the fission rate enhanced models for the thermal mobility and chemical mobility of various phases are developed in this study.The predictions in Figure 10a indicate that the irradiation enhancement plays an important role on the migration of Zr atoms.Here, the predictions of the Zr atom fraction within the HCF fuel are also obtained, without considering the fission rate enhanced terms in Equations ( 24), ( 25) and (29).
Figure 10a exhibits the contour plot of Zr atom fraction at 10% FIMA, and the predictions on Path 1 are given in Figure 10b for several burnup levels.It can be seen that when the accelerating effects of the fission rate on the mobility of Zr atoms are not taken into account, the redistribution of Zr atoms is not obvious even at 10% FIMA, and the redistribution region exists only nearby the α + γ phase boundary.As shown in Figure 10b, the Zr atom fractions are almost unchanged in the α phase region; in the other regions, the Zr atom fractions vary at a much slower rate compared to the results in Figure 8a, with the highest value of ~0.28 and the lowest value of ~0.056 at 10% FIMA.As shown in Figure 11, the diffusion flux magnitudes of Zr atoms at 10% FIMA are generally 200 times smaller than those with irradiation enhancement involved.In fact, the diffusion flux is majorly contributed by the large chemical potential gradient at the end of reactor startup.

Discussion
In this study, an improved model system and 3D multi-field coupling computational method for Zr atom redistribution behaviors are proposed for U-10Zr fuels, involving the irradiation-enhanced effects.The Zr atom redistribution behaviors within the U-10Zr fuels of HCF rods are investigated, under the linear power condition of 80 kW/m.It was discovered that the Zr atom redistribution occurs distinctly under these conditions, induced by the high temperatures, large temperature gradients and a higher fission rate.As shown in Figure 10, the irradiation enhancement can be found to be significant.Various fission rates can exist in real reactors [35], and the fission rate of a fuel rod may also vary with burnup.Therefore, the behaviors of Zr atom redistribution within U-Zr fuels deserve further investigation with the irradiation-enhanced models, for different cases of fission rate.
As mentioned in the part of Introduction, the mechanical fields have not been coupled with the other fields in most of the current studies.Actually, the irradiation-induced deformations will greatly affect the temperature field and the related Zr atom redistribution.It is known that the fission-induced irradiation swelling of U-Zr fuels will gradually change the macroscale configurations of fuel pellets or fuel meats.Especially, the differed behaviors of fission gas swelling and release will result in dissimilar internal microstructures [12,44], varying the local porosities and volumetric deformations in different phase regions.As a result, the macroscale average heat generation rate and the homogenized material parameters will change with the irradiation time and location, affecting the temperature field.
Currently, the multi-field simulation of Zr atom redistribution is generally based on the initial configuration of a fuel pellet, with the dynamically varying effects of configuration and local porosity considered specially.In Ref. [30], the thermal conductivity was multiplied with a time-dependent factor.In Ref. [32], some other material parameters were also adjusted to reflect the influences of irradiation swelling and microstructure variation.In this study, different adjustment factors are also employed to the thermal conductivities of integration points, as pointed out in Section 3. Due to the differences in microstructure and macroscale deformation in the three phase regions, the factors of 0.9, 0.45 and 0.75 are, respectively, assigned to the central region, the middle region and the external region across the radius of the considered fuel pellet.
Figure 3 indicates that the numerical predictions of Zr atom fraction in this study agree well with the experimental data.In fact, owing to the scattered experimental data around the zirconium-poor region, different predictions in Refs.[12,30,32] also seem to match the experimental results.It is found in this study that the diverse treatments to the thermal conductivity could lead to similar temperature values, but the local temperature gradients could differ distinctly.The distribution and evolution patters of Zr atom fraction would mainly depend on the local temperature gradients.Consequently, the predicted results of the Zr atom fraction in the references differ from each other, especially those around the zirconium-poor region.
In future work, the multi-field coupling modeling based on the dynamically updated configurations needs to be developed, with the irradiation-induced mechanical fields simultaneously involved.Thus, the effects of microstructure and irradiation-induced deformation [44,45] could be physically taken into account, instead of using adjustment factors.It should be mentioned that the U-50Zr fuels were proposed for HCF rods [16].In order to implement the corresponding multi-field coupling simulation for U-50Zr fuels under the conditions of a high fission rate or accidents, new models of free energy should be further developed.

Conclusions
In this study, the fission rate-related models for the thermal mobility and chemical mobility of U-10Zr fuel are developed.Based on the phase-field models of U-Zr redistribution, the comprehensive theoretical frame is integrated to describe the complicated multi-field coupling behaviors during the phase transformation and redistribution of Zr atoms.Correspondingly, the equivalent integral weak forms for every field are derived and presented.Based on the finite element method, the nonlinear iterative scheme for the direct-coupling computation of multi-fields is proposed.The developed models, algorithms and procedures are validated with the experimental data in the reference.Numerical simulations of the phase transformation and Zr atom redistribution behaviors within the HCF rod are carried out, with the predictions analyzed.The dominated mechanisms of zirconium-poor and zirconium-rich zones within the fuel core are revealed.The main conclusions are as follows: (1) The irradiation enhancement of the thermal mobility and chemical mobility plays a critical role in the redistribution of Zr atoms within both of the traditional fuel pellets and the HCFs; without considering this effect the experimental results could not be predicted out.(2) The obtained multi-field results have the helical symmetric characteristics; the multi-field results within a certain helical cross-section can be regarded as the rotated ones of the other cross-sections; the efficient quasi-three-dimensional finite element model can be adopted to perform the numerical analysis in this study.
(3) The U-10Zr alloy compositions are redistributed under the chemical potential gradient and temperature gradient.Under the irradiation conditions of this study, the highest temperature of the HCF rod can achieve 909 K, and the lowest temperature is ~622 K.The high temperature and the large temperature gradient lead to formation of a zirconium-rich zone in the elbow region and a zirconiumpoor zone inside the fuel core at 10% FIMA.(4) The dynamitic redistributions of Zr atoms are closely related to the heterogeneous high-temperature phase transformation after the end of reactor startup and the evolution of the phase fractions, depending on the competitive relations between the contributions of the chemical potential gradients and the temperature gradients within the α phase and γ phase.At the end of the start-up, the maximum γ phase fraction is 0.4; with the increase in burnup, the maximum γ phase fraction increases to 0.55 and then decreases to 0.38.

Figure 1 .
Figure 1.Geometry of a helical cruciform fuel rod.

M
γ are the chemical mobilities of the Zr atom in the corresponding phase.

u d stands for the array of 1 N
× for the node variable of u; [ ] matrix calculated with the known nodal values after the iterative step m, given as

Figure 3 .Table 1 . 4 . 1 .
Figure 3. Predictions of multi-field coupling U-Zr redistribution behaviors for DP-81 fuel.(a) The predicted results for DP-81 fuel at 3.6% FIMA together with the experimental data; (b) the simulation results of Zr atom fraction within the DP-81 fuel pin at different burnup levels.Table 1. Parameters required in calculation.Parament Value Refs.0 D α 

Figure 4 .
Figure 4. Mesh grids for (a) the whole helical cruciform fuel rod and (b) a 1/4 section of a slice.

Figure 5 .
Figure 5.The multi-field results along three typical helical paths of HCF rods at 5% FIMA.(a) Contour plots of the temperature and the Zr atom fraction in HCF rods; (b) the corresponding results along the three output paths across the rod elbow, blade and center; (c) distribution of Zr atomic fraction within the cross-sections at different heights.

Figure 6 .
Figure 6.(a) Predictions of temperature, phase fraction and Zr atom fraction on Path 1 during the reactor startup, with the gray lines representing the process results during startup; (b) the contour plots of multi-field results at the end of reactor startup.

Figure 9 .
Figure 9.The diffusion flux of Zr atoms and the corresponding contributions of temperature gradient and chemical potential gradient (a) after reactor startup (Arrow Magnification: 1), (b) at 0.2% FIMA (Arrow Magnification: 2) and (c) at 1% FIMA (Arrow Magnification: 4).Due to the high magnitude of chemical potential gradient along Path 1 shown in Figure 8c, Zr atoms around the α + γ phase boundary migrate rapidly within 0.2% FIMA, as depicted in Figure 8a.Simultaneously, it can be seen from Figure 9a that within a short time after the end of reactor startup, a part of the Zr atoms surrounding the α + γ phase boundary migrate towards the fuel center, creating a local Zr atom fraction increase to ~0.44 in Figure8a.Different from the abnormal diffusion behaviors in the γ phase, in the α phase region Zr atoms will migrate against the temperature gradient, described with Equations (19) and(28).Therefore, the Zr atoms in the α phase region can be found from Figure9ato diffuse away from the fuel center, induced by the temperature gradient there.Under the combined contributions of temperature gradient and chemical potential gradient (as shown in Figure9a) to the Zr atom flux, the zirconium-poor region is firstly formed outside the α + γ phase boundary, as shown in Figures7 and 8a,d.It can be found from Figure8dthat from the end of reactor startup to 0.2% FIMA the boundary of α + γ phase moves a little away from the fuel center, and the locally enhanced γ phase fraction can be

Figure 10 .
Figure 10.(a) Contour plots of Zr atom fraction at 10% FIMA without the irradiation enhancement effects and (b) the evolution results along Path1 from the end of reactor startup to 10% FIMA.