Gas Bubble Evolution in Polycrystalline UMo Fuels Under Elastic-Plastic Deformation: A Phase-Field Model With Crystal-Plasticity

In monolithic UMo fuels, the interaction between the Al cladding and large gas bubble volumetric swelling causes both elastic-plastic and creep deformation. In this work, a phase-field model of gas bubble evolution in polycrystalline UMo under elastic-plastic deformation was developed for studying the dynamic interaction between evolving gas bubble/voids and deformation. A crystal plasticity model, which assumes that the plastic strain rate is proportional to resolved shear stresses of dislocation slip systems on their slip planes, was used to describe plastic deformation in polycrystalline UMo. Xe diffusion and gas bubble evolution are driven by the minimization of chemical and deformation energies in the phase-field model, while evolving gas bubble structure was used to update the mechanical properties in the crystal plasticity model. With the developed model, we simulated the effect of gas bubble structures (different volume fractions and internal gas pressures) on stress-strain curves and the effect of local stresses on gas bubble evolution. The results show that 1) the effective Young’s modulus and yield stress decrease with the increase of gas bubble volume fraction; 2) the hardening coefficient increases with the increase of gas bubble volume fraction, especially for gas bubbles with higher internal pressure; and 3) the pressure dependence of Xe thermodynamic and kinetic properties in addition to the local stress state determine gas bubble growth or shrinkage. The simulated results can serve as a guide to improve material property models for macroscale fuel performance modeling.


INTRODUCTION
The measured thickness increase of a UMo fuel meat in monolithic UMo plate-type fuel as a function of the local fission density is shown in Figure 1. The profile of an irradiated UMo fuel plate and gas bubble structures at different fission densities are also shown in Figure 1. Irradiation results in a nonuniform thickness of the UMo fuel meart and large changes in the gas bubble structure with respect to fission density. The interaction between gas bubble volumetric swelling of the UMo fuel meat and the Al cladding in monolithic UMo fuels results in a non-uniform stress field in both the UMo fuel meat and Al cladding, which drives elastic-plastic and/or creep deformation, defect diffusion, and microstructure evolution (Kim et al., 2013;Meyer et al., 2014;Meyer et al., 2017;Ozaltun and Rabin, 2017;Medvedev et al., 2018). The fission product induced swelling generates high stresses within both the fuel and the cladding, leading to plasticity and creep. The interaction of gas bubble swelling with this induced mechanical deformation will subsequently affect the nature of gas bubble evolution, due to accommodation of stress. On the other hand, since the gas bubbles have different thermo-mechanical properties from UMo matrix, the evolution of gas bubbles impacts the thermomechanical properties of the entire fuel system. Therefore, to predict the fuel performance including mechanical integrity, geometric stability, and stable irradiation behavior, it is crucial to validate the correlation among deformation, gas bubble evolution, and material properties, and develop physics-based constitutive equations of swelling, elastic-plastic deformation, and creep.
The existing constitutive models for swelling, elastic-plastic deformation, and creep utilized in monolithic UMo fuel performance modeling (Miller et al., 2010;Ozaltun and Rabin, 2017;Medvedev et al., 2018;Yan et al., 2019) are empirical and developed primarily by fitting experimental data Burkes and Wachs, 2010;Kim et al., 2013) and secondarily to fit the results of finite element analyses. The mechanisms of deformations such as plastic and creep deformation and the coupling of gas bubble structure, stress, and property evolution are not sufficiently considered, which limits the predictive capability of the existing models.
In the absence of irradiation damage, a number of microstructure-dependent plastic and creep deformation models have been developed (Deutchman et al., 2012;Cottura et al., 2016;Wu and Sandfeld, 2017;Yang et al., 2018). A homogenized crystal plasticity FEM Model (Deutchman et al., 2012), which uses crystal plasticity parameters (such as activation energy, passing stress and activation volume) provided by a dislocation-density based crystal plasticity modeling, was developed to study the effect of various microstructures (precipitate shape and volume fraction, and channel width) on plastic deformation in Ni based superalloys. An empirical constitutive model-based phase-field model (Tsukada et al., 2011), which uses the von Mises yield criterion for plastic deformation and a simple creep evolution equation, was developed to study the evolution of microstructure and inelastic (plastic and creep) deformation during high temperature creep in nickel-based superalloys. Yang et al. (2018) developed a phase-field model combining a simple creep damage model (Vladimirov et al., 2009) and a dislocation dynamics model (Wu and Sandfeld, 2017) for studying the evolution of c ′ /c phases during creep in nickelbased single crystal superalloys. Phase-field models of diffusive transformations with dislocation density based plastic deformation (Cottura et al., 2016), grain boundary sliding with crystal plasticity deformation (Cheng et al., 2019), and continuum dislocation dynamic-based crystal plasticity (Koslowski et al., 2002;Koslowski and Ortiz, 2004;Zeng et al., 2016), in addition to continuum dislocation dynamics-based crystal plasticity (Hochrainer et al., 2014;Hochrainer et al., 2016), have also been developed to study dislocation structure evolution during deformation. These models strengthen the simulation capability in studying the interaction between phase separation, dislocation structure evolution, and plastic deformation. Previously, the authors have developed a number of microstructure evolution models in UMo fuels including gas bubble superlattice formation (Hu et al., 2016), radiation-induced recrystallization , and defect clustering with non-equilibrium gas bubble association growth . In this work, we leverage the existing computational capability of gas bubble evolution and crystal plasticity to develop a phase-field model of gas bubble evolution in polycrystalline UMo fuels under elastic-plastic deformation to study the effect of gas bubble structures, such as volume fraction and internal pressures, on mechanical properties as well as the effect of local stress on gas bubble evolution.

PHASE-FIELD MODEL OF GAS BUBBLE EVOLUTION IN POLYCRYSTALLINE UMO UNDER ELASTIC-PLASTIC DEFORMATION
We consider gas bubble evolution in a representative volume of polycrystalline UMo in monolithic UMo fuels as shown in FIGURE 1 | Experimental results reproduced from references (Kim et al., 2013;Meyer et al., 2017). (A) UMo fuel meat thickness change vs. fission density (Meyer et al., 2017), (B) profile of radiated UMo fuel meat, and (C) evolution of recrystallized grains and gas bubbles in UMo fuel meat (Kim et al., 2013). Figure 2. The axes xyz denote the global coordinate while x β y β z β is the local coordinate of grain β. The orientation of grain β in the global coordinate is described by the Euler angle (ϕ β , θ β , ψ β ). The volumetric swelling associated with fission gas bubbles in a UMo fuel meat is constrained by the Al cladding, which produce stresses, and both plastic and creep deformation. The interaction between gas bubble swelling and inelastic deformation affects gas bubble growth as well as stress and strain, and hence, fuel performance. To describe the gas bubble evolution in the polycrystalline UMo under irradiation and inelastic (plastic plus creep) deformation, two sets of field variables are used to describe the microstructure. One set is the order parameters η β (r, t) and χ(r, t). η β (r, t) describes the grain orientations and χ(r, t) describes the gas bubbles. η β (r, t) is equal to 1 inside grain β, 0 outside the grain β, and varies from 0 to 1 across the grain boundaries. The order parameter χ(r, t) is equal to 1 inside gas bubbles, 0 outside gas bubbles, and varies from 0 to 1 across the interface between gas bubble and matrix. The other set of field variables includes defect concentrations c li (n, r, t) and c Xe (r, t). c li (n, r, t), describes the concentrations of vacancy and interstitial clusters; where li denotes a cluster with defect i which could be vacancy or interstitial, while l means cluster or loop, and n is the number of defect i in the cluster li. c Xe (r, t) describes the concentration of fission gas atoms, here treated as Xe, and r, t are global coordinate and time, respectively.
A phase-field model of grain growth (Chen, 2002) is used to generate polycrystalline structures. The grain boundaries are defined by a shape function η(r, t) 2 β 0 β 1 (1 − η β ) 2 , which has the value of 0 inside the grains and continuously varies to 1 at the center of grain boundaries, and β 0 is the total number of grains in the simulation cell. Grain boundaries and gas bubbles are structural defects, which are sinks for vacancies and interstitials. The spatial distribution of sinks can be defined by θ(r, t) η(r, t) + χ 2 (r, t).
In UMo fuels, 235 U fission generates high-energy neutrons and fission fragments that cause radiation damage. A cluster dynamics model (Mansur, 1994;Kohnert and Wirth, 2015;Brimbal et al., 2016;Kohnert et al., 2018) is used to describe the evolution of vacancies, interstitials, and their clusters in polycrystalline structures with distributed gas bubbles. The generation of gas atoms, vacancies, and interstitials are calculated with the fission product yields and the kinetic energy distribution of the fission products. Grain boundaries, gas bubbles, and dislocations are treated as sink and emission sites of defects. The description of the cluster dynamics model of defect evolution is given in . In this work, we focused on the static or dynamic interaction between elastic-plastic deformation and gas bubbles. For simplicity, we assumed that interstitial concentration is low and fission gas Xe atoms occupy the vacancy or vacancy cluster sites. For a given Xe concentration, the vacancy concentration affects thermodynamic and kinetic properties such as lattice mismatch and diffusivity. So only Xe concentration c Xe (r, t) is taken into account in this work.
The KKS model (Kim et al., 1999) is used to describe the gas bubble evolution in polycrystalline UMo. The total free energy G of the system is formulated as a functional of the order parameter field χ(r, t) and concentration field c Xe (r, t) as where V is the material volume of the simulation cell, c m and c b are the Xe concentration in matrix and as bubble phases, respectively, f m (c m ) and f b (c b ) are the free energy density of matrix and gas bubble, respectively, h(χ) is a shape function having values between 0 and 1, g(χ) is the double-well potential, w is the double-well potential height, κ is the gradient energy coefficient, and f def is the deformation energy density. The shape function h(χ) and the double-well potential are selected as: For simplicity, the free-energy densities are approximated as quadratic functions of concentrations: where A im and A ib , (i 1, 2, 3) are free energy coefficients and f 0 is a parameter which scales the chemical driving force, which is related to the deformation energy driving force of Xe diffusion in the matrix.
The KKS model assumes that 1) at any material point the chemical potentials of the matrix and gas bubble phases are the same; and 2) Xe concentration is conserved. These two assumptions are described as: From Eqs 6, 7, we obtain a single solution of c p m (r, t) and c p b (r, t). The temporal evolution of the field variables are given by Allen-Cahn (Cahn and Allen, 1977) and Cahn-Hilliard (Cahn, 1961) equations: where D(χ, η m ) is the diffusion of Xe, the order parameters χ and η m are used to describe an inhomogeneous Xe diffusivity, L is the interfacial mobility coefficient, and is the total free energy density.
There are three unknown model parameters in the evolution Eqs 8, 9, i.e., w, κ, and L. These parameters can be determined by material properties including interface thickness 2λ, interface energy c, and interface mobility ]. Analyzing the equilibrium properties of the kinetic Eqs 8, 9 via thin interface limit analysis (Kim et al., 1999;Hu et al., 2007), the relationships can be obtained: where α is a coefficient that depends on the definition of the interface and is set to α 2.2, and the parameter L can be obtained from thin interface analysis . But for a diffusion-controlled process like the gas bubble evolution in this work, L can be chosen to be a large value to ensure a stable solution. The deformation energy density, f def , in Eq. 2 depends on the local elastic-plastic energy. The calculation methods incorporating elastic-plastic energy will be described in Elastic-Plastic Deformation Under Crystal Plasticity Framework and Material Properties of UMo.

ELASTIC-PLASTIC DEFORMATION UNDER CRYSTAL PLASTICITY FRAMEWORK
With the assumption of small deformation, the deformation energy density can be calculated by: where C ijkl (r) is the elastic constant tensor, ε e ij is the elastic strain, σ appl ij is the applied stress tensor, and ε kl is the average strain tensor.
The elastic strain is expressed as where ε ij is the total strain and ε p ij is the total eigenstrain associated with lattice mismatch between the matrix and the pressured gas bubble and plastic deformation. The eigenstrain is defined as where ε 0 (c Xe ) is the eigenstrain of gas bubbles which can be estimated by the equation of state of Xe gas phase inside the gas bubble, δ ij is the Kronecker delta function, and ε p ij (r, t) is the plastic strain which is calculated from crystal plasticity theory.
According to crystal plasticity theory, the plastic strain rate at the point r inside grain β can be generally calculated as (Ma and Roters, 2004;Ma et al., 2006): where _ c s (r), τ s 0 (r), and m s (r), are the shear strain rate, the critical resolved shear stress, and the Schmid tensor, respectively. The superscript s denotes the slip system s at material point r, and N is the total number of the slip systems of the crystal at material point r. _ c 0 is a normalization factor and n is the stress exponent (inverse of the rate-sensitivity exponent). σ ′ (r) is the deviatoric stress tensor. The Schmid tensor is a dyadic tensor and is calculated using: where b s and n s are the Burger's vector and the normal direction of slip system s at point r inside grain β. Then, the total plastic strain rate is calculated as We use A β (a β ij ) denoting the rotation matrix of the local coordinate of grain β. The coordinate transfer of the second order tensor ε, such as stress, strain and diffusivity tensors, from local coordinate to global can be written as where A T β is the transpose of A β . With the Euler angles of grain β, the components of the rotation matrix are given as Assuming that gas bubble phase has isotropic elastic constants C b ijkl and the single crystal UMo has anisotropic elastic constant C M ijkl , the elastic constant tensor in the global coordinate, can be described as To calculate deformation energy density in Eq. 12, we need to solve the shear strain rate _ c s (r) from Eq. 15. In this work, a formulation based on FFTs (Lebensohn et al., 2012) is employed to solve for the shear strain rate _ c s (r). Here we summarize the method as follows. At time t + Δt, the total strain includes elastic strains and plastic strains at a material point r: where ε(r) represents the total strain tensor, ε e (r) is the elastic strain tensor, ε p (r) is the viscoplastic strain tensor, and _ ε p (r) is the viscoplastic strain rate tensor. The viscoplastic strain rate _ ε p (r) is constitutively related to the local deviatoric stress, σ ′ (r) σ(r) − p(r)I, where p(r) (1/3)[σ 11 (r) + σ 22 (r) + σ 33 (r)] and I being the hydrostatic stresses and a unit matrix, respectively, via a sum over the N active slip systems described by Eq. 15.
The Euler implicit time discretization scheme is used to solve the solution of Eq. 21. The expression, in small strains, of the stress tensor at material point r at t + Δt is given by where σ(r) is the Cauchy stress tensor, and c(r) {c ijkl (r)} is the elastic stiffness tensor. The stresses must satisfy the stress equilibrium equation and associated boundary conditions.
For known plastic deformation strain ε p,t (r) at step t, the stress σ t+Δt (r), strain ε t+Δt (r), and plastic strain rate _ ε p,t+Δt (r) at time step t + Δt can be obtained by the following two steps: Step 1. Seek solutions of τ t+Δt (r) and e t+Δt (r) for the following equations by removing the superscript t + Δt, and keep the previous time step superscript t. The stress, τ(r), satisfies the equilibrium Eq. 23: The strains, e(x), are related to the displacements, u(r), as follows: Combining Eqs 24, 25, we have We use iteration and FFT to solve Eqs 25, 27 and let the obtained stresses and displacements satisfy the given boundary condition. The boundary condition is satisfied in the concept of average values. For example, (1) For a polycrystal under uniaxial tensile stress along the x 1axis with a strain rate of _ ε 11 , the boundary condition is given by e 11 ε t 11 + _ ε 11 Δt and τ 22 τ 33 τ 23 τ 13 τ 12 0 where ε t 11 is the average value of ε t 11 from the previous step t and is known for the current step.
(2) For a polycrystal under a constant pressure stress σ 0 along the x 1 -axis with a shear strain rate _ ε 12 and a fixed side-boundary condition to mimic billet material inside a die chamber, the boundary condition can be expressed as τ 11 σ 0 , e 22 e 33 e 23 e 13 0 and e 12 ε t 12 + _ ε 12 Δt, where ε t 12 is known for the current step, similar to ε t 11 .
Step 2. To get the final solutions of σ t+Δt (r), ε t+Δt (r), and _ ε p,t+Δt (r), or σ(r), ε(r), and _ ε p (r) without the superscript t+Δt for Eqs 21, 22 under given boundary conditions, a residual function R ij (r) is introduced (Lebensohn et al., 2012): where τ ij (r) and e ij (r), have been obtained from Step 1, and σ ij (r) and ε ij (r) will be solved through nullification of R ij (r) coupled with Eqs. 15, 22. The nullification of Eq. 15 is solved using a Newton-Raphson (N-R) scheme, i.e., where the superscript l denotes the l-th iteration step. The iteration is stopped when R ij is less than a predetermined value. Through this step, we can finally get σ ij (r), ε ij (r), and _ ε p ij (r) for the given boundary condition and time step. For materials with strength hardening, τ s 0 (r) varies with _ c s (r). For example, linear strength hardening can be expressed by Δτ s 0 H N s 1 _ c s (r)Δt with H being a constant. In such a case, τ s 0 (r) in Eq. 16 is replaced by Through Steps 1 and 2, σ ij (r), ε ij (r), and, _ ε p ij (r) are obtained. With a known strength hardening law such as Eq. 30, the shear strain rate _ c s (r) can be obtained from Eq. 15.

MATERIAL PROPERTIES OF UMO
The thermal and mechanical properties of UMo that depend on temperature and neutron fluence (Medvedev et al., 2018) have been assessed by experiments (Farrell and King, 1979;Polkinghorne and Lacy, 1991;Kaufman, 1999;Farrell;Meyer et al., 2016). The temperature dependence of Young's modulus E (GPa) is expressed as (Beghi, 1968), where T is the temperature (K). The Poisson's ratio was adapted from  and it is constant 0.324. The temperature dependence of yield strength σ y is expressed as (Klein), Since we do not have the yield stress of single crystal UMo, Eq. 32 is used to estimate the critical resolved stress in the crystal plasticity model. Formation energy and migration energy of Xe are calculated by atomistic simulations (Smirnova et al., 2012;Smirnova et al., 2013;Smirnova et al., 2015).
Both experiments and atomistic simulations show that gamma UMo has isotropic elastic properties. Thus, two of the three elastic moduli, i.e., Young's modulus E, Shear modulus G, and Poisson's ratio v, can describe the elastic properties of UMo. In this work, the temperature is set to be 500 K, which is approximately the operation temperature of UMo fuels in high performance research reactors (Meyer et al., 2016). From Eqs 30, 32, E is 70 GPa and the yield strength σ y is 0.718 GPa at T 500 K. The constant H in Eq. 30 is set to be 5.0 × 10 6 Pa. γ−UMo has a bodycentered cubic (bcc) structure, where 24 slip systems are often activated during deformation; 12 · {110}111 and 12 · {211}111 systems. The crystal plasticity model is general and can consider all the slip systems. For model validation, only 12 · {110}111 slip systems are considered in the simulations. In addition, the evolution of radiation defects such as vacancies, interstitials, and their clusters are not considered in the current phase-field model, whereas only Xe diffusion and Xe gas bubble evolution are taken into account. Table 1 lists the model parameters in the simulations. Figure 2 illustrates the simulation cell with dimensions of 128l 0 × 32l 0 × 128l 0 , cylindrical grains along y-direction, and distributed gas bubbles. The average grain size in the xz plane is about 340 nm, which is on the order of the typical grain size observed in recrystallized grains in UMo fuels. Periodic boundaries conditions are applied in the x-, y-, and zdirections and a strain along the z-direction is applied to perform tensile or compressive deformation. In this work, the developed model was applied to: 1) study the effect of gas bubble structures on mechanical response; and 2) study the effect of stresses and stress-dependent thermodynamic and kinetic properties on gas bubble evolution.

Effect of Gas Bubble Structures on Mechanical Properties
Three gas bubble structures with gas bubble volume fractions (V f 3.5%, 6.7%, and 9.7% ) are generated with a phase-field model of gas bubble evolution in polycrystalline structures. Gas bubbles, which have an average gas bubble size of 100 nm in diameter are randomly distributed in the simulation cell. It is assumed that bubbles are pressurized and that the pressure is associated with the lattice mismatch between the gas phase and the matrix UMo phase. The lattice mismatch is described by an eigenstrain tensor ε gbp ij ε b0 δ ij c p b (r, t)h(χ), where c p b (r, t) is the Xe concentration inside the gas bubble, h(χ) is the shape function defined by Eq. 3, δ ij is the Kronecker-Delta function, and ε b0 is the mismatch strain. For the first order approximation, if the bulk modulus, pressure, and Xe equilibrium concentration inside the gas bubble are B gb , p gb and c eq b , respectively, the mismatch strain can be  (Harrison, 1969;Beeler et al., 2020) can be used to estimate the bulk modulus B gb and the equilibrium Xe concentration for given a pressure, hence, the mismatch strain ε 0 . From the EOS , when the internal pressure is about 2 GPa the bulk modulus is about 30 GPa. Here, we assigned the elastic constants of the gas phase to be c b 11 90E 0 GPa, c b 12 30E 0 GPa and c b 44 30E 0 GPa, the bulk modulus B gb is 50E 0 GPa, and the Poisson's ratio is 0.125, where E 0 is a parameter which depends on the pressure p gb and the concentration c eq b inside the gas bubble. Pressure, equilibrium concentration, and lattice mismatch inside gas bubbles change with the local stresses and chemistry (local vacancy and Xe concentrations) in UMo fuels in service. For evolving gas bubbles, the molar volume is calculated by the Xe concentration inside the gas bubble. With the molar volume, the pressure and bulk modulus can be calculated with the EOS. To study the effect of steady state gas bubble structures on mechanical response, we can prescribe fixed values of c eq b and ε 0 , which are listed in Table 1, and vary E 0 to describe the pressure inside the gas bubbles.

Stress Field Around Pressured Gas Bubbles
In the elastic-plastic deformation model, the iteration approach (Hu and Chen, 2001) is used to solve the mechanical equilibrium equations and the stress field in an elastic inhomogeneous material with a distribution of stress-free strains as described in Eq. 15. The stress field around gas bubbles with average radius of 50 nm and different internal pressures ((P gb 0.07, 0.60, 1.2 and 2.1 GPa)) under elastic deformation is calculated. Stress fields on the middle plane around gas bubbles are presented in Figure 3.
The light black lines denote the grain boundary while the white circles show the interfaces between gas bubble and matrix. It is found that the pressure (P −(σ 11 + σ 22 + σ 33 )/3) inside the gas bubbles is uniform which is in agreement with Eshelby's solution (Eshelby, 1957), and the shear stress (σ 13 ) around the gas bubble is larger than the yield stress of UMo (0.718 GPa) when the internal pressure is larger than 1 GPa. The internal pressure inside nanosized gas bubbles may reach a few GPa according to MD simulations (Xiao and Long, 2014;Beeler et al., 2020), but with the increase of gas bubble size, the pressure decreases. In addition, a stress field associated with the cladding constraint in UMo monolithic fuels might increase the stresses in the matrix. Therefore, the internal pressure and the cladding constraint may result in plastic deformation in UMo under service.

Effect of Gas Bubble Structures on Stress-Strain Curves
Gas bubble structures with different volume fractions (V f 3.5%, 6.7%, and 9.7%) and different initial internal pressures (P gb 0.07, and 1.2 GPa) are used to study the effect of gas bubble structures on stress-strain curves under elastic-plastic deformation. In the simulations, a strain rate of dε 33 /dt 3 × 10 − 4 (1/s) (the other strain components are zero, ε ij 0 ) is applied in z-direction for tensile deformation while dε 33 /dt −3 × 10 − 4 (1/s) is applied for compressive deformation. Xe concentration in the matrix is set to be 5 × 10 − 5 and the stress-free strain associated with Xe induced lattice change in the matrix is set to be 0.1. Figure 4A,B presents the effect of gas bubble structures on stress-strain curves under tensile and compress deformation. The black curves are the stress-strain curves in polycrystalline structures with Xe concentration of 5 × 10 − 5 , but without gas bubbles. The results in Figure 4A are stress strain curves for gas bubbles with a low initial internal pressure of P gb 0.07 GPa, while the results in Figure 4B are for gas bubbles with a higher initial internal pressure of P gb 1.2 GPa. Because of the lattice mismatch associated with distributed Xe in the matrix and the internal pressure inside gas bubbles, a residual stress field is present. The residual stress, which is a compressive stress field due to a positive stress-free strain, shifts the total stress at applied strain of zero to a negative value. The negative stress value is marked by the small circle in the stress-strain curves. It can be seen that the effect of gas bubble volume fraction on the stress shift at an applied strain of ε 33 0 is small, especially for the case of gas bubbles with a low internal pressure P gb 0.07 GPa. For gas bubbles with high pressure, the stress shift increases with the increase of gas bubble volume fraction which can be seen by zooming in on the stress-strain curves at ε 33 0 in Figure 4B.
For the simulation conditions, the residual stress is mainly determined by the distributed Xe (its concentration and stressfree strain) in the matrix. If the matrix vacancy is rich, the distributed vacancies and Xe results in the reduction of UMo lattice constant, and the stress-free strain is negative. It is expected that the residual stress is a tensile stress field due to a negative stress-free strain, and the stress-strain curves shift to a positive value at ε 33 0. Comparing the results in Figure 4, we can conclude that 1) for all the cases the effective Young's modulus, which is the slope at the linear part of the stress-strain curves, decreases with the increase of gas bubble volume fraction. This is expected because the gas phase has a lower Young's modulus than that of the matrix UMo phase.
2) The Young's modulus depends on both gas bubble structure (gas bubble volume fraction and internal pressure) and applied stress (tensile or compress). 3) the yield stress decreases with the increase of gas bubble volume fraction. The yield stress has a similar dependence on gas bubble structure and applied stress as that of the Young's modulus; and 4) the hardening coefficient increases with the increase of gas bubble volume fraction, especially for gas bubbles with higher internal pressures, which is indicated by the slop of stress-strain curves in the plastic deformation stage.
The strain hardening is determined by the plastic strain rate. The distributions of plastic strain ε p 13 on the center plane S in polycrystalline structures with gas bubble volume fraction 9.7% at different applied strain ε 33 are shown in Figure 5. The results in Figures 5A,B are for gas bubbles with pressure of P gb 0.07 GPa and P gb 1.2 GPa under tensile deformation, respectively. Before the applied strain reaches ε 33 0.02, the deformation is elastic and the plastic strain is zero as shown in Figure 5. It is observed that plastic deformation first takes place near the gas bubble interface, particularly at the interface region of two nearby gas bubbles as shown ε 33 0.054, where the stress concentration is higher than that at the interface of an isolate gas bubble, as shown in Figure 3A.
With the increase of applied strain, plastic strain increases. The plastic strain in regions with yellow color has a positive sign while the plastic strain in regions with green color has a negative sign. The flaky pattern of plastic strain (ε p 13 ) distribution at ε 33 0.1 indicates the formation of shear bands where shear strain has a uniform and high value. Figure 6 plots the distributions of the shear stress component σ 13 on the center plane S at ε 33 0.1. The white lines show the <101> directions. From the results in Figure 6, we can see that 1) most bands of shear stress σ 13 align along the <101> directions while the effect of grain orientation on shear stress σ 13 is minor. The isotropic elastic properties of UMo, which has the Zener ratio (2C 44 /(C 11 − C 12 )) of 1, can explains the grain orientation independence of shear stresses, and 2) the alignment of gas bubbles along the <101> direction enhances the shear stress bands for both cases of gas bubbles (with low and high initial pressures). Compared with the shear stress, the bands of shear plastic strain (ε p 13 ) shown in Figure 5 does not well align along the <101> directions. This is because dislocation sliding depends not only on the resolved shear stress but also grain orientations. The red and blue of the color bar in Figures 5, 6 present the maximum and minimum values of shear strain (or stress) in the simulation cell during deformation for a given gas bubble structure with low (or high) pressure. Comparing the maximum values in the color bars in Figures 5, 6 we can see that both the maximum plastic strain and shear stress for gas bubbles with low pressure are larger but more localized near the gas bubbles than that for gas bubbles with high pressure. In other words, the shear stress and strain fields around gas bubbles with a low gas pressure is more inhomogeneous than those around gas bubbles with a high gas pressure. We also calculated the evolution of total shear plastic strain in the simulation cell during the FIGURE 4 | Effect of gas bubble volume fraction and internal pressure on stress strain curves. Results shown are for gas bubbles with a pressure of (A) P gb 0.07 GPa, (B) P gb 1.2 GPa. Both tensile and compressive stresses are applied.
deformation. The results show that the total plastic strain for a system with low pressure gas bubbles is higher than that for a system with high pressure gas bubbles. Therefore, we can conclude that the more inhomogeneous a stress field is, the less strain hardening is. And the gas bubble dependence of hardening behavior showed in Figure 4 is attributed to the inhomogeneous stress induced inhomogeneous plastic deformation.
For a given gas bubble structure, the average pressures (P (σ 11 + σ 22 + σ 33 )/3) in the matrix and inside gas bubbles are calculated during deformation. Figure 7 shows the average pressure in terms of applied strain. The circle on the curves presents the pressure associated with the residual stress at applied strain ε 33 0.0. Under tensile deformation the pressure P in the matrix changes from negative to positive. This means that the average volume of atoms in matrix increases. As a result, the formation energies of defects with positive lattice misfit such as large solutes and self-interstitials should decrease, and their migration energy should decrease. In contrast, under compressive deformation, the pressure P in the matrix simply becomes more negative. This causes the formation energies and migration energy of defects with positive lattice misfit (such as Xe substitutionals or interstitials) to increase. The average pressure inside the gas bubbles also depends on the applied strain. Based on the equation of state of the gas phase, which describes the correlation between gas atom concentration and pressure, the gas concentration inside gas bubbles should decrease with a pressure decrease. Therefore, the average pressure changes in the matrix and inside the gas bubble affect the thermodynamic and kinetic properties of diffusive defects in the matrix and inside gas bubbles, as well as the gas bubble evolution.

Effect of Elastic-Plastic Deformation on Gas Bubble Evolution
The gas bubble evolution model presented in Section Phase-Field Model of Gas Bubble Evolution in Polycrystalline UMo Under Elastic-Plastic Deformation requires thermodynamic and kinetic properties including diffusivities and the equilibrium concentration of Xe in the matrix and inside gas bubbles. We describe the diffusivity and equilibrium concentration as a function of average pressure in the matrix and inside gas bubbles. In principle, these thermodynamic and kinetic properties can be obtained from atomic simulations (Valikova and Nazarov, 2008;Sueoka et al., 2012). For instance, the equation of state of Xe gas phase can be used to determine the correlation between the equilibrium Xe concentration and pressure inside gas bubbles. The pressure dependent formation energy of Xe in the matrix can be used to develop the chemical free energy of Xe in the matrix phase, and the pressure dependent migration energy can be used to assess the correlation between Xe diffusivity and pressure. Since we don't have the data for UMo, a linear correlation of the thermodynamic and kinetic properties in terms of pressure near the equilibrium state are used to validate the effect of elastic-plastic deformation on gas bubbles.
In the simulations, an initial Xe concentration of 0.005 and distributed gas bubbles with a volume fraction 3.5% and initial pressure 2.1 GPa are introduced into polycrystalline UMo. The stress-free strain, ε 0 , of Xe in the matrix is assumed to be 0.1 or −0.1. The positive ε 0 describes a state which is Xe-rich in the matrix, causing a compressive lattice environment, while the negative ε 0 describes a vacancy-rich state in the matrix, causing a tensile lattice environment. For a given initial structure, the pressure of the residual stress is chosen as the equilibrium state. The equilibrium Xe concentrations in the matrix and gas phase are calculated by c eq i (P i ) (c eq i (P i0 ) + ω i (P i − P i0 )c eq i (P i0 ))/(P i,max − P i0 ). P i0 is the pressure of the residual stress at the equilibrium state. P i,max is the maximum pressure in different gas bubble structures for a give applied strain, i b or m, which denote the gas bubble and the matrix, respectively. ω i is a coefficient which is set to be 0.3, allowing a maximum 30% change of equilibrium concentration for the pressure changes. The pressure P i,max is estimated from the data in Figure 7. For simplicity, the pressure dependence of Xe diffusivity assumes that the diffusivity of Xe in a vacancy-rich environment, i.e., ε 0 −0.1, is one order magnitude larger than that in a Xe-rich environment, ε 0 0.1. The rest of the model parameters of the phase-field model are listed in Table 1. The simulations are carried out in three steps: 1) apply a strain rate, dε 33 /dt 3 × 10 − 4 (1/s) or dε 33 /dt −3 × 10 − 4 (1/s), up to a total applied strain of ε 33 0.05 for a tensile stress state, or ε 33 −0.05 for a compressive stress state; 2) update the diffusivity and equilibrium concentrations in the matrix and inside gas bubbles based on the average pressures; 3) simulate the gas bubble evolution with the phase-field model, update gas bubble structures, elastic and plastic properties; update the elastic-plastic solution every 1,000 diffusion time steps, and return to step 2) to simulate the dynamic interaction between elastic deformation and gas bubble evolution. Figure 8 shows the final gas bubble structures at time t 240 h for four simulation cases, and Figure 9 plots the evolution of gas bubble volume fraction. The circle at t 0 shows the volume fraction of initial gas bubbles. The solid lines present the gas bubble evolution under tensile stress states while the dashed lines denote the cases under compressive stress states. The results demonstrate that 1) a tensile stress state leads to gas bubble growth, while a compressive stress states causes a slight shrink of gas bubbles; 2) a vacancy rich environment promotes gas bubble growth under a tensile stress state while it prevents gas bubble from shrinking under a compressive stress state. Figure 10 shows the Xe concentration distribution on the center plane S at time t 240 h for the four simulation cases. From the color bar it can be seen that under a compressive stress state the Xe concentration inside the gas bubble, which is shown in Figures  10A, C, is about 0.75, which is an increase from 0.6 at the initial state due to the increase of internal pressure. The Xe concentration inside the gas bubble under a tensile stress state, which is shown in Figures 10B, D, changes from 0.6 at the initial state to around 0.4 due to the reduction of internal pressure inside gas bubbles. The white lines show the contour of the order parameter, χ, in the phase-field model which is equal to 0.5, representing the interface of gas bubbles and matrix. It is interesting to find that the gas bubbles no longer have a spherical shape. The gas bubble shape change can be explained by the fact that an inhomogeneous stress field around the gas bubble results in an inhomogeneous plastic deformation, pressure, and non-uniform equilibrium concentration, hence, different interface evolution kinetics. Experimental data also show the gas bubbles have an irregular shape (Miller et al., 2015;Gan et al., 2017). Beside the factor of elastic-plastic deformation, the solid fission phase precipitation on the interior of the  gas bubble surface can also affects the gas bubble interface evolution kinetics, hence, the gas bubble shape.

CONCLUSION
In this work we developed a phase-field model of gas bubble evolution in polycrystalline UMo under elastic-plastic deformation. Plastic strain rate-based crystal plasticity is employed to describe the elastic-plastic deformation. The effect of local pressure on the thermodynamic and kinetic properties of Xe is taken into account. With the developed model, we simulated the effect of gas bubble structures (different volume fraction and internal pressure) on stress-strain curves and the effect of local stress fields on gas bubble evolution. The results show that 1) the effective Young's modulus decreases with the increase of gas bubble volume fraction; 2) the yield stress decreases with the increase of gas bubble volume fraction; 3) the hardening coefficient increases with the increase of gas bubble volume fraction, especially for gas bubbles with higher internal pressure; 4) a tensile stress state leads to gas bubble growth while a compressive stress states causes a slight shrink of gas bubbles; 5) a vacancy rich environment promotes gas bubble growth under a tensile stress state while preventing gas bubble shrinkage under a compressive stress state; and 6) an inhomogeneous stress field around the gas bubble results in an inhomogeneous plastic deformation, pressure, and non-uniform equilibrium concentration, hence, an irregular gas bubble shape. The results demonstrate that the developed model is capable of studying the effect of gas bubble structures on mechanical properties and studying the dynamic interaction between elastic-plastic deformation and evolving gas bubbles. However, additional thermodynamic and kinetic properties, such as the effect of pressure on chemical free energy and mobility of Xe in matrix and gas phases, are required for more quantitative simulations. The model also needs to be extended to take radiation defect accumulation and potential creep deformation into account.

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

AUTHOR CONTRIBUTIONS
SH: model and code development, simulation, visualization, and writing-original draft and editing. BB: conceptualization, writingreview and editing.