- Email: [email protected]
Simulation of the fission-induced swelling and creep in the CERCER fuel pellets
Materials and Design 89 (2016) 183–195
Contents lists available at ScienceDirect
Materials and Design journal homepage: www.elsevier.com/locate/jmad
Simulation of the fission-induced swelling and creep in the CERCER fuel pellets Yunmei Zhao a, Xin Gong a, Yi Cui b, Shurong Ding a,c,⁎ a b c
Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China Department of Mechanical Engineering, University of Alberta, Edmonton, AB T6G 2G8, Canada Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institution of China, Chengdu, 610041, Sichuan, China
a r t i c l e
i n f o
Article history: Received 20 June 2015 Received in revised form 20 September 2015 Accepted 23 September 2015 Available online 26 September 2015 Keywords: CERCER Simulation Recrystallization Fission gas swelling Irradiation-induced creep Thermo-mechanical coupling
a b s t r a c t A ceramic-ceramic (CERCER) fuel is selected as a dedicated composite fuel for an Accelerator Driven System (ADS), with minor-actinide-enriched ceramic fuel particles dispersed in a ceramic MgO matrix. The multi-scale simulation method of the irradiation-induced thermo-mechanical behavior evolution in the CERCER fuel pellets is developed. The semi-analytical fission-induced gaseous swelling model is established with considering recrystallization process and gas atom resolution in fuel grain scale. It is involved in the macro-scale three-dimensional large-deformation incremental constitutive relation for fuel particles. The stress update algorithms are developed respectively for the fuel, and the MgO matrix in which different fission-induced creep effects are taken into account. The user subroutines to define their mechanical constitutive relations are correspondingly written and validated. The simulation results indicate that (1) recrystallization leads to accelerated intergranular bubble swelling, and the maximum bubble swelling locates in the central particle where the temperature is the highest; (2) the creep coefficient in the MgO fission-induced creep model has a remarkable effect on its mechanical behavior evolution; and it can be determined by the fracture pattern in the matrix; (3) the matrix can effectively shorten the radial deformation of the composite pellet. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction Minor actinides (MA, including Np, Am and Cm) have a 0.1% share of nuclear wastes in the spent fuels [1]. They have long-term radiotoxicity. Their continuous accumulation constitutes a huge threat to the environment and also affects the sustainable development of nuclear energy. ADS [2] is proposed as an effective means to transmute these actinides. Its subcritical reactor is equipped with a neutron source. As a kind of advanced nuclear fission energy system, it allows generating nuclear power and achieving significant decrease of the inventory of MA simultaneously. The nuclear fuel element is the basic component in the subcritical reactor of ADS, loading with MA and fissile nuclides. Its structure is similar to the nuclear fuel rod used in the current nuclear plants, with an original gap designed between the fuel pellets and the cladding. The aim is to avoid the intensive pellet-cladding mechanical interaction during its lifetime. The pellet in the ADS fuel rod is no longer a homogeneous one, but a composite one. In the EUROTRANS project, two kinds of candidates are highly favored [3–5]: one is CERCER with ceramic fuel particles dispersed in the ceramic MgO matrix and the other is CERMET with ⁎ Corresponding author at: Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China. E-mail address: [email protected] (S. Ding).
http://dx.doi.org/10.1016/j.matdes.2015.09.135 0264-1275/© 2015 Elsevier Ltd. All rights reserved.
ceramic fuel particles embedded in the metallic Mo92 matrix. For these innovative fuel pellets, it is an important issue to investigate their irradiation-induced thermo-mechanical behavior evolution. It can supply a basis for structural optimization of ADS fuel rods. For example, the radial deformation with increasing burnup should be understood in advance in order to design the original pellet-cladding gap. In the earlier research, UO2 fuel particles were used to replace the ones containing MA. In the framework of EFRTTRA [1,6,7] a large number of irradiation experiments were implemented for UO2-based CERCER pellets. These experiments focused on the irradiation-induced thermo-mechanical behavior, including the size stability, the effects of different matrix materials, the effect of particle sizes and the fracture characteristics. MgO was selected as a superior matrix for CERCER pellets because of its irradiation stability. The nuclear fuels in the CERCER pellets undergo sophisticated behavior in the extreme irradiation environment. In addition to generating fission heat, fission-induced irradiation swelling stems from solid and gaseous fission products. Distinguished from the homogeneous UO2 pellet, the produced fission gases can't be easily released since the fuel particles are strictly restrained by the MgO matrix. Irradiation experiments demonstrated that recrystallization [8,9] occurred in the fuel grains, by which more grain boundaries are formed to accommodate fission gas atoms. Recrystallization initiates from the boundary of the original grain and proceeds toward its center until the grain is totally
184
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
consumed. Large grains turn into small fine-grains so as to absorb gas atoms in the grain boundary effectively, thus large intergranular bubbles will produce and make a major contribution to the remarkable swelling. The considerable fission-induced swelling will cause enhanced mechanical interaction between the particles and the matrix, which results in large deformation. Simultaneously, the MgO matrix engenders fission-induced creep [7], which is perceived as an important mechanism to trigger the observed cracking [7]. In addition to the irradiation experimental researches [10], the corresponding theoretical models and numerical simulation method should also be developed to model the irradiation-induced thermomechanical behavior in the CERCER fuels, so as to give better comprehension of their in-pile behavior and benefit their optimal design. Above all, the fission gas swelling model should be established with the recrystallization process involved. J. Rest [9,11] established a recrystallization theory for UO2 and UMo fuel. Y. Cui [12] further considered the resolution effect [13] to develop a mechanistic model for the fission gas behavior in UMo fuels, and the semi-analytical expressions for the gaseous swelling were obtained. Simulation research on the threedimensional thermo-mechanical behavior is limited. The thermomechanical coupling behavior at the initial stage [10,14] were investigated, without considering the fuel swelling. S. Ding [15] simulated evolution of the in-pile micro-thermo-mechanical behavior, in which the irradiation creep was considered in the MgO matrix with the backward integration algorithm; however, the empirical swelling model suitable for the homogeneous UO2 fuel pellet was adopted. As mentioned above, the recrystallization-affected gaseous swelling should be considered in the fuel particles. In this study, the CERCER pellet [7] is chosen as the research object, in which UO2 is taken as the fuel particles. Based on Y. Cui's work [12], the theoretical fission gas swelling model for UO2 is built. The macroscale three-dimensional mechanical constitutive relations and stress update algorithms are established, with the grain-scale gas swelling model involved in the fuel particles and the irradiation creep considered in the matrix. Especially, the mid-value integration algorithm for the creep strain increments is adopted in the matrix. The formed userdefined subroutines [16,17] are validated and introduced to the finite element method (FEM) simulation of the multi-scale thermomechanical behavior within the CERCER pellet. The evolution rules are investigated and analyzed. The effect of the creep coefficient is examined.
2.1. Irradiation swelling in UO2 For the UO2 fuel particles in the CERCER pellet, the total swelling includes two contributions. One comes from the solid fission products, whose linear relationship with the burnup is given as [19] ΔV solid ¼ 2:5 10−29 F d V
ð1Þ
which is the volumetric variation relative to the original volume, that is the volumetric engineering strain; Fd depicts the fission density in fission/m3. The other swelling contribution results from the fission gas products, whose relation with burnup is of complexity. In this study, according to the theoretical model for the fission gas diffusion with the resolution effect and recrystallization effect involved [12], the theoretical model for the fission gas swelling is obtained and adopted.
2=15
Before recrystallization (F ≤ Fdx, F dx ¼ 4 1024 ðf Þ
, whereFdx is
the critical fission density; f is the fission rate in fissions/m3s), the fission gas swelling takes the form as the sum of the intragranular bubble swelling and the intergranular bubble swelling ΔV gas 4π 3 2πR3b C b ¼ : r b cb þ V r gr0 3
ð2Þ
3 In Eq. (2), 4π 3 r b cb depicts the intragranular bubble swelling, where r b
¼
bv nb 1=3 ð3hs4π Þ
¼
16π f n r g D0 c2g b0 nb
is the radius of the intragranular bubble [20], and cb
is the average concentration of intragranular bubbles. And
cgis the volume-averaged total gas content inside the grain, which is formulated as
cg ¼
−1 þ
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 64π f n r g D0 =b0 Y f t−3N=2rgr
32π f n r g D0 =b0
ð3Þ
where f is the fission rate, D0 =10−39m5is the proportionality constant
such that Dg ¼ D0 f , b0 = 2 × 10−23m3 is the proportionality constant
such that b ¼ b0 f , N is the concentration of the gas atoms in the grain boundary with the obtained solution as 2. Material properties For the considered CERCER pellet, the fission gas swelling model with consideration of recrystallization and gas atom resolution in grain scale is given for UO2 particles and the irradiated creep model for the MgO matrix is presented in this section. The other material properties used in the FEM simulation are the same as those in Ref. [15,18].
pffiffi pffiffi 8 2 3 2^ 2^ > es2 t erfc s2 ^t es1 t erfc s1 ^t ^ > r gr 5 ^ 2 2 2 β > > 4 ^ tbr gr =π t þ þ þ > > ^ ^ s2 ðs1 −s2 Þ s1 ðs2 −s1 Þ bλ < bλ ^ N t ¼ ð4Þ ! ^ 3 ^ 3 e−u21^t=r2gr > ^t 4βr > 2βr 1 gr gr > 2 2 ^ > tbr þ − =π > gr > 2 : 3 þ k1 r 2gr 5ð3 þ k1 Þ u21 u21 þ 3k1 þ k1
Fig. 1. (a) schematic figure of the fuel pellet with periodically-distributed particles (b) the selected RVE model (c) mesh grid.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
185
After recrystallization (F ≥Fdx), the original grain is divided into two regions including the recrystallized one and the unrecrystallized one. The total fission gas swelling is expressed as "
#
ΔV gas ΔV intra ΔV inter ΔV inter þ Vr ¼ ð1−V r Þ þ V V V V r gr r gr r grx
ð8Þ
where Vris the volume fraction of the recrystallized region; rgris the current grain radius excluding the recrystallized outside; rgrx = 0.1 μm is the average radius of the fine grains [29]. And the specified formulations for Vr and rgr can be found in Ref. [12]. In the recrystallized area, in consideration of the gas depletion [25, 26] taking place in the recrystallized area, the gas atom number on the grain boundary is obtained as [11] Fig. 2. Selected nodes in (a) the central particle and (b) the matrix.
^ are the expressions with the transforma^ and b where the parameters ^t, β Y f ^ ¼ b0 , the specific details can be found in ^ tion induced as ^t ¼ Dg t; β ¼ ; b
Nbx ¼
Y f t 3C bx =2r grx þ 1=8r 3grx
ð9Þ
Dg
Dg
Ref. [12]; b′ is the intergranular resolution rate with the employed value of 1×10−5s−1. Other parameters can be referred in Ref. [12]. Another item in Eq. (2)
2πR3b C b r gr0 ,
describes the intergranular bubble
swelling, where rgr0 is the original grain radius; Cbis the grainboundary bubble concentration, which is expressed as 0 qffiffiffiffi3 11=2 a dN 8z 12 ^ C b ¼ @ 2 dt A π ξδ
ð5Þ
3
a where 12 is the average atomic volume [28] with a is set as 5.47 × 10− 10m, z = 4 is the number of sites explored per gas-atom jump [11], ξ = 15is the grain-boundary diffusion enhancement factor [15], δ =2 × 10−9mis the width of the boundary [27], dN is the diffusion d^t
flux from grain interior to the grain boundary with the obtained form as
8 pffiffi pffiffi ^ 2β s2 s1 2^ 2^ > > ^tbr2 =π 2 > 1þ es2 t erfc s2 ^t þ es1 t erfc s1 ^t gr < ^ s1 −s2 s2 −s1 dN bλ ¼ 2^ 2 t=r −u ^ ^ > 2βr gr 4βrgr e 1 gr d^t > ^tNr2 =π 2 > − : gr 3 þ k1 u2 þ 3k1 þ k2 1
1
ð6Þ ^
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
In Eq. (4) and Eq. (6), s1 ¼ bλ 2 −
^2 λ2 b 4
^
^
þ rbλ , s2 ¼ bλ 2 þ gr
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^2 λ2 b 4
^
þ rbλ , u1 ¼ gr
^ gr . and k1 ¼ bλr 4:493− k14:883 þ3:612 Rb in Eq. (2) is the radius of the intergranular bubble, which satisfies the modified Van der Waals gas law [11], and neglecting the effect of hydrostatic pressure yields
2γ 4πR3b −hs bv Nb 3 Rb
where Cbx is the intergranular bubble density with the following expression as [11]
! ¼ Nb kT
ð7Þ
where γ = 1N/m is the surface tension [26], and Nb ¼ Nð^tÞ=C b is the gas atom number per intergranular bubble. It is known that Rbshould be obtained through the nonlinear iteration.
C bx
0 11=2 8zadN d^t r¼r grx C B ¼ @ 1=3 A : 12 π 2 ξδ
ð10Þ
Owing to the gas depletion, the fission gas atoms all migrate to the grain boundaries. As a consequence, the gaseous swelling is only contributed by the intergranular bubble swelling in the recrystallized area. According to Ref. [11], the volumetric gas swelling in the recrystallized region is obtained as !
ΔV intra 4πR3bx 3C bx 1 ¼ þ V 3 2r grx 8r 3grx r grx
ð11Þ
where Rbx is the radius of the intergranular bubble in the recrystallized area and it also satisfies the modified Van der Waals gas law as Eq. (7) with Nb replaced by Nbx and Rb replaced by Rbx. One can see that in the unrecrystallized part the gas swelling is still composed of the intragranular one and the intergranular one; while in the recrystallized part the gas swelling is only derived from the intergranular one because the gas depletion takes place in the recrystallized area [20] and almost all the generated gas atoms migrate to the grain boundaries directly. The detailed calculation formula is described in this section, and it needs to be pointed out that the intergranular bubble radius still should be obtained through the nonlinear iteration. 2.2. Irradiation creep in MgO Continuously attacked by fast neurons and other high energy fragments, irradiation creep occurs in the MgO matrix. The specific creep model for MgO hasn't been found in references. Here an irradiation creep model similar to CVD-SiC [21] is adopted as cr
ε ¼ c0 ψσ
ð12Þ
Table 1 Validation of the thermal constitutive relation for node M in the fuel particle. Time (day)
(K/mm) The temperature gradient ∂T ∂z
The heat flux qz obtained in ABAQUS (W/mm2)
The conductivity coefficient k used in ABAQUS (W/mm·K)
The conductivity from the model in Ref. [15] (W/mm·K)
The relative errors (%)
57.5 115 172.5
−68.14623 −82.50341 −96.03734
0.178414 0.1837 0.190923
0.002618 0.002227 0.001988
0.0026593 0.0022411 0.0020097
1.553 0.6291 1.080
186
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
Table 2 Validation of the thermal constitutive relation for node N in the matrix. Time (day)
The temperature gradient ∂T ∂z (K/mm)
The heat flux qz obtained in ABAQUS (W/mm2)
The conductivity coefficient k (W/mm·K)
The conductivity from the model in Ref. [15] (W/mm·K)
The relative errors (%)
57.5 115 172.5
−24.8468 −25.46282 −28.90954
0.313413 0.330782 0.354707
0.01261 0.012508 0.01227
0.012793 0.012699 0.012346
1.431 1.504 0.619
which is given based on the fact that a linear relationship exists between the irradiation creep and the Mises stress for ceramic materials [22]; cr
where ε is the equivalent irradiation creep strain rate(/s); c0is the creep coefficient in (MPa ⋅ n ⋅ cm−2)−1; σdenotes the equivalent Mises Stress in MPa and ψdepicts the fast neutron flux in n ⋅ cm−2s−1. In this study, the effects of different creep coefficients on the calculation results are investigated.
Then, the constitutive relation at time t + Δt can be obtained as
eðtþΔtÞ
eðtþΔt Þ
σ tþΔt ¼ 2GðT þ ΔT; t þ Δt Þε ij þ λðT þ ΔT; t þ Δt Þεkk δij ij h i h i eðt Þ eðt Þ ¼ 2GðT þ ΔT; t þ Δt Þ ε ij þ Δε eij þ λðT þ ΔT; t þ Δt Þ εkk þ Δεekk δij eðt Þ
¼ 2GðT þ ΔT; t þ Δt Þεij
þ 2GðT þ ΔT; t þ Δt ÞΔε eij
þ λðt þ Δt; T þ
þ λðT þ ΔT; t þ Δt ÞΔε ekk δij :
eðt Þ ΔT Þεkk δij
ð14Þ 3. Three-dimensional stress update algorithms To simulate the thermo-mechanical coupling behavior in ADS fuel pellets accurately, it is of importance to develop the three-dimensional irradiation-effect-involved mechanical constitutive relations. The whole calculation is divided into a number of time increments. The incremental mechanical constitutive relations and stress update algorithms for UO2 and MgO are established in this section. The userdefined subroutines to define the mechanical constitutive relations can be written in accordance with the algorithms, and ultimately be used in the finite element simulation. Being a large-deformation problem, the three-dimensional constitutive relations and stress-update algorithms are derived in a rotating coordinate system. Considering a time increment[t, t + Δt], the threedimensional incremental constitutive relation is first given. At each integration point, the relationship between the Cauchy
An incremental constitutive relation for an integration point yields Δσ ij ¼ σ ijtþΔt −σ tij ¼ 2GðT þ ΔT; t þ Δt ÞΔεeij þ λðT þ ΔT; t þ Δt ÞΔεekk δij eðt Þ
þ 2ΔGε ij
eðt Þ
ð15Þ
þ Δλεkk δij
where Δσij are the stress increments and Δεeij are the elastic logarithmic strain increments; ΔG and Δλ are the increments of the Lame coefficients. Based on the incremental constitutive relation, the threedimensional stress update algorithms for the fuel particles and the matrix are to be developed in the following.
eðtÞ
stresses σ tij and the elastic logarithmic strains εij at time t is given as eðt Þ
σ tij ¼ 2GðT; t Þεij
eðt Þ
þ λðT; t Þεkk δij
ð13Þ
where λ and Gare the Lame coefficients, which are related to both temperature and time.
3.1. The three-dimensional stress update algorithm for UO2 For the fuel particles, assuming that the total strain increments consist of the elastic ones, the thermal expansion ones, the irradiation swelling ones and the plastic ones, then the elastic incremental strains in a logarithmic form are expressed as p th sw Δε eij ¼ Δεtotal ij −Δε ij −Δε ij −Δε ij
Fig. 3. (a) Contour plot of the temperature in the pellet on the 115th day; (b) evolution of the temperature along Path 1.
ð16Þ
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
187
Fig. 4. (a) Contour plot of the swelling strains on the 230th day in different fuel particles and (b) the selected node P in Particle 1.
where Δεtotal denote the total strain increments; Δεth ij ij are the logarithmic thermal strain increments which can be derived as Δε th ½ ln ð1 þ α TþΔT ðT þ ΔT−T 0 ÞÞ− ln ð1 þ α T ðT−T 0 ÞÞδij ij ¼
1 þ α TþΔT ðT þ ΔT−T 0 Þ δij ¼ ln 1 þ α T ðT−T 0 Þ
ð17Þ
where αis known as the secant thermal expansion coefficient in K−1 and T0 is the reference temperature in K. The irradiation swelling strain incrementsΔεsw ij are expressed as
Δε sw ij ¼
ln ð1 þ SW ðt þ ΔtÞÞ swðt Þ δij −ε ij 3
ð18Þ
Fig. 5. (a) Validation of the volumetric swelling of node P; (b) evolution of the gas swelling in the unrecrystallized region; (b) evolution of the radius of the intergranular bubble.
188
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
Fig. 6. (a) The volumetric swelling distribution and (b) temperature distribution on the 115th day in the selected particles.
where SW(t +Δt) represents the volumetric swelling engineering strain ΔV
ΔV solid gas with SWðt þ ΔtÞ ¼ ΔV V ¼ V þ V , and the detailed expressions can be found in Section 2.1, which are obtained according to the grainscale theoretical analysis; εsw(t) denote the irradiation swelling logarithij mic strains at time t, which are known for the concerned time increment. As shown in Eq. (2), the intergranular swelling is dependent on the intergranular bubble radius. While the intergranular bubble radius Rb satisfies the nonlinear equation, as depicted in Eq. (7). Thus Newton– Raphson iteration method is applied to solve Eq. (7) written in a unified form as
!
4πR3b Δ 2γ −hs bv N b −N b kT ¼ 0 gðRb Þ ¼ 3 Rb
ð19Þ
where Rb denotes the intergranular bubble radius in the recrystallized or unrecrystallized region, and the definitions of the other variables can be found in Section 2.1. Eq. (19) is solved using a scheme of the Newton–Raphson (N–R) type, i.e. ∂g Rkb ¼ Rkb −g Rkb = ð20Þ Rkþ1 b ∂Rkb where ∂g Rkb ∂Rkb
¼ 16πγRkb =3 þ
2γhs bv Nb 2 : Rkb
ð21Þ
Once convergence is achieved, SW(t + Δt) can be worked out, then the irradiation swelling strain incrementsΔεsw ij can be obtained. For the integration points whose Mises stresses satisfy the yield condition, the contribution of the plastic strain increments should be considered. According to the plastic constitutive theory [23], we have dεpij ¼
3sij p dε 2σ
ð22Þ
where dεpij denotes differential plastic strain components; σ depicts the qffiffiffiffiffiffiffiffiffiffiffi Mises stress with σ ¼ 32 sij sij ; sij is the deviatoric stress component with sij ¼ σ ij − σ3kk δij ; dεp is the differential equivalent plastic strain. Backward Euler integration of Eq. (22) yields
Δεpij ¼
3sijtþΔt 2σ tþΔt
Δεp
ð23Þ
which implies that the plastic strain increments are related to the unknown stresses at time t+ Δt. Furthermore, the equivalent plastic strain increment Δεp needs to meet the consistent condition as described in Ref. [15]. For an integration point, firstly the plastic strain increments are supposed to be zero, the trial stresses at time t + Δtare obtained; the corresponding Mises stress is judged whether the yield condition is satisfied, and if not, the trial stresses are the real ones; if the yield condition is met, the plastic strain increments should be further included, and the
Fig. 7. Radius of intergranular bubbles in (a) unrecrystallized region and (b) recrystallized region on the 115th day.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
189
Fig. 8. Contour plot of the equivalent creep strains in the matrix on the (a) 11.5th day and (b) 115th day.
subsequent stress update algorithm is similar to the one presented in Ref. [17]. 3.2. The three-dimensional stress update algorithm for MgO As for the matrix, the corresponding total strain increments Δε-ij are constituted of the elastic strain ones Δεeij, the thermal expansion cr ones Δεth ij and the creep strain ones Δεij . Hence, the elastic strain increments can be derived as total
th cr Δε eij ¼ Δεtotal ij −Δε ij −Δε ij :
ð24Þ
Mid-value integration of Eq. (26) yields 3 stþΔt þ stij ij Δεcr Δεcr ij ¼ 2 σ tþΔt þ σ t
ð27Þ
whereΔε cr represents the equivalent creep strain increment. According to Eq. (12), the equivalent creep strain increment can be similarly derived as Δεcr ¼ c0 ψ
1 tþΔt σ þ σ t Δt: 2
ð28Þ
The thermal expansion logarithmic strain increments are demonstrated as
Substituting Eq. (28) into Eq. (27), the creep strain increments can be obtained as
Δε th ij ¼ ½ ln ð1 þ α TþΔT ðT þ ΔT−T 0 ÞÞ− ln ð1 þ α T ðT−T 0 ÞÞδij 1 þ α TþΔT ðT þ ΔT−T 0 Þ δij ¼ ln 1 þ α T ðT−T 0 Þ
Δεcr ij ¼ ð25Þ
3 c0 ψ sijtþΔt þ stij Δt: 4
ð29Þ
Combining Eq. (15) with Eq. (24), we have where α is the secant thermal expansion coefficient and T0 in K is the reference temperature. Based on the creep constitutive theory [23], the differential creep strains dεcr ij are expressed as
dεcr ij ¼
3sij cr dε 2σ cr
th cr σ ijtþΔt ¼ σ tij þ 2Gðt þ Δt; T þ ΔT Þ Δεtotal ij −Δε ij −Δε ij eðt Þ eðt Þ th þ λðt þ Δt; T þ ΔT Þ Δεtotal kk −Δε kk δij þ 2ΔGε ij þ Δλε kk δij :
ð30Þ
ð26Þ
where dε is the differential equivalent creep strain.
Firstly, the trial stresses σpr ij can be obtained assuming that no creep strain increments take place, namely, the components Δεcr ij in Eq. (24)
Fig. 9. (a) Validation of the creep algorithm in the matrix; and (b) evolution of the equivalent creep strains along Path 2.
190
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
Fig. 10. Contour plots of (a) the displacement magnitude and (b) y-directional displacement on the 230th day.
are set as zero. So, the corresponding Cauchy stresses σtij+ Δt can be expressed as 3 pr cr tþΔt σ ijtþΔt ¼ σ pr þ stij c0 ψΔt ij −2GΔε ij ¼ σ ij − 2 G sij
ð31Þ
where G represents the shear modulus at time t + Δt. The deviatoric part of Eq. (31) becomes 3 tþΔt þ stij c0 ψΔt: sijtþΔt ¼ spr ij − 2 G sij
ð32Þ
Then, the deviatoric stresses can be obtained as
sijtþΔt ¼
3 t spr ij − 2 sij Gc0 ψΔt 3 1 þ Gc0 ψΔt 2
ð33Þ
pr t+Δt tþΔt 1 pr where stþΔt ¼ σ ijtþΔt − 13 σ kk δij , spr =σpr kk. ij ¼ σ ij − 3 σ kk δij , σkk ij
can be updated as Therefore, the Cauchy stresses σt+Δt ij σ ijtþΔt ¼ sijtþΔt þ
σ pr kk δij : 3
ð34Þ
4. Finite element model A composite pellet containing UO2 particles and MgO matrix is taken as the research object. In this simulation, the geometric model is set up
with the assumption that the fuel particles are periodically-distributed along the axial direction and the circumferential one, as illustrated in Fig. 1(a). The pellet has a radius of 4.15 mm, the particles are supposed to be spherical with the radius of 200 μm and the volume fraction 6.23%. Considering the periodicity and symmetry, the Representative Volume Element (RVE) is built, as shown in Fig. 1(b). A pressure of 10 MPa and a constant temperature of 873 K are applied on the outside surface of the RVE. The symmetric boundary conditions are applied on the other surfaces of the RVE, with considering that the main deformation is along the radial direction of a pellet in the ADS fuel rod. Similar to the irradiation conditions in the experiment [7], the fission rate of fuel particles and the fast neutron flux are set as 2.5 × 1020 fission/m3s and 2.5 × 1015n/cm2s, respectively. Thus the corresponding heat generation rate of fuel particles is 8 W/mm3. As shown in Fig. 1(c), the RVE is discretized into 88,617 elements and 124,547 nodes with the element C3D10MT in ABAQUS, under which the convergent calculation results can be obtained. The special thermo-mechanical constitutive relations of fuel particles and the matrix are defined with the user-defined subroutines. 5. Results and discussions The finite element simulation of the multi-scale thermo-mechanical behavior evolution is successfully performed. The effectiveness of the user-defined subroutines should be confirmed. Lacking a series of definite experimental results, in this section some numerical results and the relative theoretical ones are compared in order to validate the correct definition of the material performances.
Fig. 11. (a) Evolution of the radial displacement along Path 3 with increasing burnup; (b) y-directional strain in the matrix on the 230th day.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
191
Fig. 12. Contour plot of the first principal stresses (a) at the initial stage; (b) on the 34.5th day;(c) on the 115th day and (d) the ones with no creep on the 115th day.
Besides, the evolution laws of the multi-scale thermo-mechanical behavior will be given and analyzed, and the effects of different creep coefficients in Eq. (12) will be investigated. 5.1. Temperature field The thermal flux and the temperature gradient comply with Fourier's law described as qz ¼ −k
∂T ∂z
ð35Þ
where K is the conductivity in W/m·K, which varies with temperature and time for the fuel particles, and for the matrix it is just
temperature-dependent; T depicts temperature in K; ∂T denotes the tem∂z perature gradient along the z-direction in Fig. 1(b); qzin W/mm2 represents the z-direction heat flux component. In order to confirm whether the thermal conductivity has been correctly defined, node M in the central particle and N in the matrix are selected to analyze their results, as shown in Fig. 2. In ABAQUS the heat flux of a node is supplied, while the temperature gradient hasn't been output. Then, the z-direction temperature gradient of node M is ob−T M2 , where TM1andTM2 are the temperature results of the tained asTzM1 M1 −zM2
neighboring nodes of node M in Fig. 2(a); zM1 and zM2 are the z coordinates. Simultaneously, the z-direction temperature gradient of node N can be obtained similarly. Thus, the fluxes and temperature gradients
Fig. 13. Evolution of the first principal stress along Path 2 with increasing burnup for (a) the case with creep and (b) the case without creep in the matrix.
192
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
Fig. 14. Equivalent creep stains along Path 2 on the172.5th day.
at different burnups are given in Table 1 and Table 2. Hence, the actually used thermal conductivity of nodes M and N can be calculated through Eq. (35), and their theoretical ones can be simultaneously obtained according to the model in Ref. [15]. The relative errors between them are slight, as depicted in Table 1 and Table 2, which demonstrates that the defined thermal constitutive relations of fuel particles and matrix are correct. The contour plot in Fig. 3(a) gives the temperature distribution on the 115th day. Thus, a representative path-Path 2 is selected to output the temperatures at different burnups, as depicted in Fig. 3(b). It can be found that (1) the temperature decreases along the radial direction due to the heat generation in the fuel particles; (2) the fuel particles have higher temperatures than the close-by matrix, which stems from the much lower thermal conductivity of the particles; (3) the peak temperature value locates at the center of the pellet and it goes up as the burnup increases, and the maximum value on the 230th day is about 118 K higher than that at the initial stage of burnup. Additionally, one can expect that larger volume fraction of fuel particles would bring much higher temperatures. In this case, a constant temperature condition is applied on the outer surface. If the surface temperature also increases with burnup, the maximum temperature will also be higher at a certain burnup. Furthermore, compared to our previous study [15], the maximum temperature is higher in this study, which is due to the higher heat generation rate adopted for the fuel particles. 5.2. Multi-scale results of volumetric swelling In this study, the used fission gas swelling model is obtained from our analytical solution of the governing equations for fission gas diffusion in grain scale, which involves the resolution effect and recrystallization effect. The details are demonstrated in Section 2.1 and its
Fig. 16. Evolution of the axial stresses along Path 2 considering creep in the matrix.
correctness has been validated [12]. As mentioned in Section 2.1 that the irradiation swelling calculation is complex, as a result the userdefined subroutine UMAT is designed for updating stresses. Firstly, the effectiveness of the introduced swelling is to be benchmarked. The contour plot of the swelling strains in different fuel particles is depicted in Fig. 4(a), which shows that the swelling strains are higher in the particles closer to the center of the pellet. Node P in Fig. 4 (b) is picked out to output the volumetric swelling strains at different burnup levels. The obtained volumetric swelling strains from UMAT are labeled as red points in Fig. 5 (a). At the same time, with the intragranular bubble radius and the intergranular bubble radius at different burnups, the theoretical results are also calculated out according to the swelling model in Section 2.1, as illustrated with the solid line in Fig. 5 (a). One can see that the theoretical results agree well with the ones obtained from UMAT, which identifies that the swelling strains are successfully introduced into the stress update. Meanwhile, the details of the swelling contributions are also given in Fig. 5 (a). It can be noticed that: (1) the solid swelling holds a proportional relationship with the irradiation time, and before the critical fission density Fdxits contribution is similar as the gas swelling. (2) After the critical fission density, the total swelling accelerates, with the gas swelling in the recrystallized region increasing promptly and the gas swelling contribution in the unrecrystallized part declining with burnup. At high burnup the gas swelling in the recrystallized region is far more than the solid swelling. This phenomenon attributes to the recrystallization process. At the critical fission density the recrystallization process begins, and the fine grains occupy the original large grain completely when the recrystallization process ends. For a very fine grain, the fission gas atoms can diffuse to the grain boundary quickly and in this study the atoms are assumed to reach the grain boundary directly without consuming time.
Fig. 15. Contour plot of the axial stresses (a) at the initial stage and (b) on the 115th day.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
193
Fig. 17. Evolution of the first principal stress along Path 2 with the creep coefficient in the matrix set as (a) 6.0×10−26(MPa n cm−2)−1; (b) 6.0×10−27(MPa n cm−2)−1; (c) 2.3×10−27 (MPa n cm−2)−1; (d) 2.0×10−27(MPa n cm−2)−1.
In addition, the evolution of the gas swelling in the unrecrystallized region is given in Fig. 5 (b). It reveals that the intragranular bubble swelling keeps a very small quantity all the time; and at the moment the recrystallization starts, the intergranular bubble swelling has a 43.7% share of the total swelling, while the intragranular swelling only account for ~0.93%. With the recrystallization proceeding, the amount of intergranular swelling of the unrecrystallization region is increasing as a whole, while the total gas swelling decreases gradually to zero due to the shrinked volume fraction. One can also find that at high burnup the gas swelling in the unrecrystallized region disappears (after the 218th day), which means that the original grain is completely replaced by the fine grains and the recrystallization process ends. For a better understanding, the intergranular bubble radius evolution is investigated for the unrecrystallized region and the recrystallized region, as shown in Fig. 5(c). It can be found that the intergranular bubble radius Rb in the unrecrystallized region first increases to 0.0287 μm and then decreases. Nevertheless, Rbx in the recrystallized region quickly jumps to a large magnitude after the recrystallization initiation and then increases with time all along. This phenomenon is consistent
The average irradiation swelling at 17.4% FIMA was obtained as 47 ± 17% from the post-irradiation-examination (PIE) [7]. Our simulation result falls in this range. This means that the established fission gas swelling model presents a good description of irradiation swelling of UO2 microspheres in the ADS fuel pellet. In Fig. 4(a), one can see that different swelling presents in the fuel particles, which majorly results from the temperature. The irradiation swelling results on the 172.5th day are depicted in Fig. 6(a). The fuel particle (Particle 6 in Fig. 4(a)) close to the outer surface of the fuel pellet has lower swelling than that in the pellet center (Particle 1), with their temperature having the similar distribution in Fig. 6(b). In accordance with the modified Van der Waals gas law as Eq. (7) in Section 2.1, high temperature will result in large intergranular bubble radius. To confirm this, the results of the intergranular bubble radii in the unrecrystallized and recrystallized regions are given in Fig. 7. It is depicted that the lowest intergranular bubble radius exists in Particle 6 with the lowest temperature, which leads to the lowest swelling. 5.3. Irradiation creep in the matrix
with the fact that the contribution of intergranular swelling V r ðΔVVinter Þjrgrx in Eq. (8) increases with time because both the volume fraction Vrof the recrystallized region and the relative bubble swelling ðΔVVinter ðRbx ÞÞ all rise. As demonstrated in Fig. 5(a), the peak swelling reaches approximately 43% on the 200th day (17.4% FIMA). The average irradiation swelling in this fuel particle is calculated out as approximately 41.3%.
The contour plots of equivalent creep strains in the matrix on the 11.5th day and the 115th day are shown in Fig. 8(a) and (b), respectively. Nodes on Path 2 are picked to output their equivalent creep strains in the time interval [0.003d, 0.006d], thus the equivalent creep strain increments can be calculated and displayed as a series of red points in Fig. 9(a). The theoretical results are also obtained through Eq. (28) with the Mises stresses at the starting time and terminal
194
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
time. It can be observed that the simulation results and the theoretical ones comply with each other, which verifies the correct definition of the creep model. Also, the evolutions of the equivalent creep strains along Path 2 are shown in Fig. 9(b). It can be found that the middle position exhibits larger creep strain and the value increases with time as a whole, and the peak value becomes 17.2% on the 172.5th day. Compared to the previous study [15], the magnitude of creep strains is larger, which is for the reason that the fission gas swelling considering the recrystallization effect is much larger in this study and it strengthens the mechanical interaction between the particles and the matrix. 5.4. Evolution of the radial deformation of the pellet In the conceptual design of ADS fuel elements, an original gap locates between the fuel pellet and the cladding in order to prevent them from an intense pellet-cladding mechanical interaction at higher burnup. In order to achieve the safety design of nuclear fuel elements, firstly, evolution of the radial deformation of ADS fuel pellets should be fully investigated. Contour plots of the displacement magnitude and the y-directional displacement on the 230th day are shown in Fig. 10. Path 3 in Fig. 10 (b) is selected to output its magnitude displacement in the Cartesian coordinate system, because in this path there is the maximum radial displacement which just equals to its y-directional component. As depicted in Fig. 11, (1) the displacements in the interfaces between the fuel particles and the matrix are larger than the ones in the adjacent positions, and the differences between the maximum one and the adjacent minimum one increase with time. This is due to the fact that the increased particle swelling enhances the mechanical interaction with the surrounding matrix, which results in appearance of larger compressive strains in the surrounding matrix; (2) relative large compressive strains occur in the cladding-pellet interface, as illustrated Fig. 11(b), which accords with the radial displacement distribution; (3) the radial displacement increases with irradiation time, and the pellet radius has an increase of ~ 78.83 μm (1.88%) on the 230th day, while the maximum displacement reaches about 89.16 μm in the vicinity of the pellet surface; the radial deformation is degraded by the constraint of the matrix; (3) compared to our previous work [15], the radial deformation of the pellet enlarges because in this study we adopt a modified swelling model considering the grain recrystallization. As a result, the original gap between the fuel pellet and the cladding ought to be appropriately designed, especially for the pellet with a large particle volume fraction. Numerical simulation of the radial deformation can supply a basis for optimal design of the ADS fuel element. 5.5. Evolution of the mechanical field in the matrix As for the matrix in the CERCER pellets, evolution of the first principal stresses is focused on, because the first principal stress is related to the failure of ceramic materials. The creep coefficient C0 in Eq. (12) was employed with a value of C0 = 2.1× 10−27(MPa n cm−2)−1. The obtained first principal stresses at different time points are, respectively, given in Fig. 12 together with the results without considering irradiation creep. Accordingly, to capture the evolution rules in detail, Path 2 in Fig. 5(a) is selected to output the first principal stresses at different burnup stages, as depicted in Fig. 13(a). It is denoted that (1) large magnitudes exist at the initial stage with the maximum at the lowest region of Path 2; (2) the first principal stresses all decrease heavily on the 11.5th day;(3) however, on the 57.5th day, the first principal stresses shows a different distribution character with the maximum value appearing at the two ends;(4) afterwards, the stresses increase with burnup as a whole, and on the 115th day the maximum first principal stresses in the vicinity of the two ends exceed the maximum value at the initial stage.
The evolution rule of the first principal stresses is induced by both effects of particle swelling and irradiation creep in the matrix. On one hand, the accumulated particle swelling enhances the interaction between the fuel particles and the matrix, which has a trend to increase the first principal stresses at the lowest region of Path 2. On the other hand, as presented in Fig. 14, the intensive mechanical interaction there results in high irradiation creep strains, and large creep strains can relax the stresses there. Both effects compete with each other, which will affect the evolution rules of the first principal stresses and the fracture pattern. Meanwhile, the first principal stresses without considering irradiation creep are also given in Fig. 13(b), which indicates that the location with the maximum value will always maintain and the enhanced irradiation swelling will lead to occurrence of extremely large stress (without considering cracking). Hence, it is revealed that the fission-induced creep effect in the MgO matrix plays an important role in its cracking pattern. For the ceramic materials, the maximum tensile stress would be the driving force to induce fracture. In order to explore the mechanism of the radial cracks discovered in the irradiation experiment [7], it is of significance to study the evolution rules of tensile stresses there. The contour plots of the axial stress component S33 are shown in Fig. 15. One can find that the axial stress components S33 at all points are negative at the initial stage, while on the 115th day S33 changes to be very large tensile ones nearby the axial symmetric plane across the fuel particles. The axial stress components S33 along Path 2 at different burnup stages are depicted in Fig. 16. It can be observed that (1) the axial stresses are compressive stresses at the initial stage; (2) the compressive stresses nearby the two ends are gradually evolved into tensile stresses and increase with burnup; nevertheless, the axial stresses at the lowest point remain compressive stresses and become larger at higher burnup. At the same time, one can find that the maximum tensile stress and the maximum first principle stress have the same value at the two ends (see Fig. 13(a) and Fig. 16), which means that the corresponding axial symmetric plane is the principal plane with the maximum tensile stresses. As a result, to determine the cracking by the maximum tensile stress, one can obtain that the radial cracks occur and propagate along this plane, which is consistent with the examined cracks in the irradiation experiment.
5.6. Effects of the creep coefficient In the previous sections, the results considering creep are calculated with the creep coefficient set as C0 = 2.1× 10−27(MPa n cm−2)−1, and it is an assumption of ours for lack of definite data. As Kim did [24], a proper parameter can be obtained through FE simulation if some experimental results can be used for reference. In this section, the effects of the creep coefficient on the mechanical behavior evolution are to be investigated. Hence, simulations are implemented using several C0 values, and the evolution of the first principal stresses along Path 2 are plotted in Fig. 17. One can find from Fig. 17 that the stress peaks locate at the lowest point at the initial stage; as time proceeds, the peaks gradually move toward the path ends. Especially, it can be observed from Fig. 17(a) that for the case withC0 = 6.0 × 10−26(MPa n cm−2)−1all the stresses decrease quickly to a very small value. Fig. 17(b) denotes that the stresses decrease with time, then increase and again decrease. When the creep coefficient is set as 2.3 × 10−27(MPa n cm−2)−1, as shown in Fig. 11 (c), after the initial stage the stresses will never exceed the original maximum. With C0 = 2.0 × 10−27(MPa n cm− 2)− 1 the maximum value at higher burnup exceeds the peaks at the initial stage, and it locates nearby the two ends. As a whole, a large creep coefficient strengthens the creep effect, thus the stresses will be heavily relaxed. If the radial cracking happens, the coefficient C0 = (2.0~ 2.1) × 10−27(MPa n cm−2)−1 will be suitable to explain the fracture mechanism with the obtained results.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
6. Conclusions For the CERCER pellet, the fission gas swelling model for UO2 is obtained with the resolution effect and recrystallization effect considered in the grain scale. The three-dimensional stress update algorithm is developed for the fuel particles, in which the irradiation swelling strains are obtained with N–R iterations. With the irradiation creep involved, using the mid-value integration to calculate the creep strain increments the stress update algorithm is also established for the matrix. Besides, the ABAQUS user subroutines are written and verified to define the thermo-mechanical properties. According to the irradiation conditions in the experiment, FE simulation of the multi-scale thermomechanical behavior is implemented. The conclusions can be drawn as follows. (1) Temperature decreases along the radial direction as a whole and the central particle has the highest temperature. The maximum temperature increases with burnup, and on the 230th day it reaches 1200 K. (2) Large volumetric swelling occurs in the particles at a higher burnup level. The recrystallization process speeds up the gas swelling. After the original grain is completely replaced by fine grains, the volumetric swelling takes on a linear relation with the burnup. (3) The recrystallization process can be captured with the intragranular swelling evolution. When the intragranular swelling disappears, the recrystallization process ends. The intergranular bubble radius in the uncrystallized zone will still grow after the recrystallization initiation, after reaching a maximum value it gradually reduces. While, the intergranular bubble radius in the crystallized zone will continuously increase with burnup, which leads to continuous increase of the total gaseous swelling. (4) The gaseous swelling in the fuel particles is different, which is mainly due to the different temperatures there. High temperature will induce large intergranular bubble radius, as a result high gaseous swelling appears in the fuel particles near the center of the pellet. (5) Fission-induced irradiation creep in the MgO matrix has an important effect on the mechanical behavior evolution. With a creep coefficient value of 2.1 × 10−27(MPa n cm−2)−1, at a certain level of burnup the maximum tensile stresses exist on the
195
axial symmetric plane, which will result in the radial cracks there. If the definite experimental results can be supplied, the suitable creep coefficient can be obtained by analyzing the FEM results.
Acknowledgments The authors thank the supports of the Natural Science Foundation of China (11172068, 91226101, 11272092), and the Foundation from Science and Technology on Reactor System Design Technology Laboratory. References [1] R.J.M. Konings, R. Conrad, G. Dassel, B.J. Pijlgroms, J. Somers, J. Nucl. Mater. 282 (2000) 159. [2] F. Delage, R. Belin, X.-N. Chen, et al., Energy Procedia 7 (2011) 303. [3] V. Sobolev, W. Uyttenhove, R. Therford, W. Maschek, J. Nucl. Mater. 414 (2011) 257. [4] D. Haas, A. Fernandez, D. Staicu, et al., Energy Convers. Manag. 49 (2008) 1928. [5] D. Haas, A. Fernandez, D. Staicu, et al., Energy Convers. Manag. 47 (2006) 2724. [6] E.A.C. Neeft, K. Bakker, R.P.C. Schram, et al., J. Nucl. Mater. 320 (2003) 106. [7] E.A.C. Neeft, K. Bakker, R.L. Belvroy, et al., J. Nucl. Mater. 317 (2003) 217. [8] J. Lamontagne, S. Bejaoui, K. Hanifi, et al., J. Nucl. Mater. 413 (2011) 137. [9] J. Rest, Compre Nucl Mater 3 (2012) 579. [10] S.E. Lemehov, V.P. Sobolev, M. Verwerft, J. Nucl. Mater. 416 (2011) 179. [11] J. Rest, J. Nucl. Mater. 346 (2005) 226. [12] Y. Cui, S.R. Ding, Z.T. Chen, Y.Z. Huo, J. Nucl. Mater. 457 (2015) 157. [13] Y. Cui, S.R. Ding, Y.Z. Huo, C.L. Wang, L. Yang, J. Nucl. Mater. 443 (2013) 570. [14] X.-N. Chen, A. Rineiski, W. Maschek, et al., Prog. Nucl. Energy 53 (2011) 855. [15] S.R. Ding, Y.M. Zhao, J.B. Wan, et al., J. Nucl. Mater. 442 (2013) 90. [16] X. Gong, Y. Zhao, S. Ding, Mech. Mater. 77 (2014) 14. [17] Y. Zhao, X. Gong, S. Ding, Y. Huo, Int. J. Mech. Sci. 81 (2014) 174. [18] MATPRO-09, A Handbook of Materials Properties for Use in the Analysis of Light Water Reactor Fuel Rod Behavior, USNRC TREENUREG-1005, 1976. [19] Y.S. Kim, G.L. Hofman, J. Nucl. Mater. 419 (2011) 291. [20] J. Spino, J. Rest, W. Goll, C.T. Walker, J. Nucl. Mater. 346 (2005) 131. [21] A. El-Azab, N.M. Ghoniem, Mech. Mater. 20 (1995) 291. [22] T.D. Gulden, J. Am. Ceram. Soc. 52 (1969) 585. [23] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons Ltd., UK, 2000. [24] Y.S. Kim, G.L. Hofman, J.S. Cheon, A.B. Robinson, D.M. Wachs, J. Nucl. Mater. 437 (2013) 37. [25] M.V. Speight, Nucl. Sci. Eng. 37 (1969) 180. [26] J. Rest, G.L. Hofman, J. Nucl. Mater. 277 (2000) 231. [27] J. Rest, G.L. Hofman, Y.S. Kim, J. Nucl. Mater. 385 (2009) 563. [28] J. Spino, D. Papaioannou, J. Nucl. Mater. 281 (2000) 146. [29] J. Rest, G.L. Hofman, J. Nucl. Mater. 210 (1994) 187.
Contents lists available at ScienceDirect
Materials and Design journal homepage: www.elsevier.com/locate/jmad
Simulation of the fission-induced swelling and creep in the CERCER fuel pellets Yunmei Zhao a, Xin Gong a, Yi Cui b, Shurong Ding a,c,⁎ a b c
Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China Department of Mechanical Engineering, University of Alberta, Edmonton, AB T6G 2G8, Canada Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institution of China, Chengdu, 610041, Sichuan, China
a r t i c l e
i n f o
Article history: Received 20 June 2015 Received in revised form 20 September 2015 Accepted 23 September 2015 Available online 26 September 2015 Keywords: CERCER Simulation Recrystallization Fission gas swelling Irradiation-induced creep Thermo-mechanical coupling
a b s t r a c t A ceramic-ceramic (CERCER) fuel is selected as a dedicated composite fuel for an Accelerator Driven System (ADS), with minor-actinide-enriched ceramic fuel particles dispersed in a ceramic MgO matrix. The multi-scale simulation method of the irradiation-induced thermo-mechanical behavior evolution in the CERCER fuel pellets is developed. The semi-analytical fission-induced gaseous swelling model is established with considering recrystallization process and gas atom resolution in fuel grain scale. It is involved in the macro-scale three-dimensional large-deformation incremental constitutive relation for fuel particles. The stress update algorithms are developed respectively for the fuel, and the MgO matrix in which different fission-induced creep effects are taken into account. The user subroutines to define their mechanical constitutive relations are correspondingly written and validated. The simulation results indicate that (1) recrystallization leads to accelerated intergranular bubble swelling, and the maximum bubble swelling locates in the central particle where the temperature is the highest; (2) the creep coefficient in the MgO fission-induced creep model has a remarkable effect on its mechanical behavior evolution; and it can be determined by the fracture pattern in the matrix; (3) the matrix can effectively shorten the radial deformation of the composite pellet. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction Minor actinides (MA, including Np, Am and Cm) have a 0.1% share of nuclear wastes in the spent fuels [1]. They have long-term radiotoxicity. Their continuous accumulation constitutes a huge threat to the environment and also affects the sustainable development of nuclear energy. ADS [2] is proposed as an effective means to transmute these actinides. Its subcritical reactor is equipped with a neutron source. As a kind of advanced nuclear fission energy system, it allows generating nuclear power and achieving significant decrease of the inventory of MA simultaneously. The nuclear fuel element is the basic component in the subcritical reactor of ADS, loading with MA and fissile nuclides. Its structure is similar to the nuclear fuel rod used in the current nuclear plants, with an original gap designed between the fuel pellets and the cladding. The aim is to avoid the intensive pellet-cladding mechanical interaction during its lifetime. The pellet in the ADS fuel rod is no longer a homogeneous one, but a composite one. In the EUROTRANS project, two kinds of candidates are highly favored [3–5]: one is CERCER with ceramic fuel particles dispersed in the ceramic MgO matrix and the other is CERMET with ⁎ Corresponding author at: Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China. E-mail address: [email protected] (S. Ding).
http://dx.doi.org/10.1016/j.matdes.2015.09.135 0264-1275/© 2015 Elsevier Ltd. All rights reserved.
ceramic fuel particles embedded in the metallic Mo92 matrix. For these innovative fuel pellets, it is an important issue to investigate their irradiation-induced thermo-mechanical behavior evolution. It can supply a basis for structural optimization of ADS fuel rods. For example, the radial deformation with increasing burnup should be understood in advance in order to design the original pellet-cladding gap. In the earlier research, UO2 fuel particles were used to replace the ones containing MA. In the framework of EFRTTRA [1,6,7] a large number of irradiation experiments were implemented for UO2-based CERCER pellets. These experiments focused on the irradiation-induced thermo-mechanical behavior, including the size stability, the effects of different matrix materials, the effect of particle sizes and the fracture characteristics. MgO was selected as a superior matrix for CERCER pellets because of its irradiation stability. The nuclear fuels in the CERCER pellets undergo sophisticated behavior in the extreme irradiation environment. In addition to generating fission heat, fission-induced irradiation swelling stems from solid and gaseous fission products. Distinguished from the homogeneous UO2 pellet, the produced fission gases can't be easily released since the fuel particles are strictly restrained by the MgO matrix. Irradiation experiments demonstrated that recrystallization [8,9] occurred in the fuel grains, by which more grain boundaries are formed to accommodate fission gas atoms. Recrystallization initiates from the boundary of the original grain and proceeds toward its center until the grain is totally
184
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
consumed. Large grains turn into small fine-grains so as to absorb gas atoms in the grain boundary effectively, thus large intergranular bubbles will produce and make a major contribution to the remarkable swelling. The considerable fission-induced swelling will cause enhanced mechanical interaction between the particles and the matrix, which results in large deformation. Simultaneously, the MgO matrix engenders fission-induced creep [7], which is perceived as an important mechanism to trigger the observed cracking [7]. In addition to the irradiation experimental researches [10], the corresponding theoretical models and numerical simulation method should also be developed to model the irradiation-induced thermomechanical behavior in the CERCER fuels, so as to give better comprehension of their in-pile behavior and benefit their optimal design. Above all, the fission gas swelling model should be established with the recrystallization process involved. J. Rest [9,11] established a recrystallization theory for UO2 and UMo fuel. Y. Cui [12] further considered the resolution effect [13] to develop a mechanistic model for the fission gas behavior in UMo fuels, and the semi-analytical expressions for the gaseous swelling were obtained. Simulation research on the threedimensional thermo-mechanical behavior is limited. The thermomechanical coupling behavior at the initial stage [10,14] were investigated, without considering the fuel swelling. S. Ding [15] simulated evolution of the in-pile micro-thermo-mechanical behavior, in which the irradiation creep was considered in the MgO matrix with the backward integration algorithm; however, the empirical swelling model suitable for the homogeneous UO2 fuel pellet was adopted. As mentioned above, the recrystallization-affected gaseous swelling should be considered in the fuel particles. In this study, the CERCER pellet [7] is chosen as the research object, in which UO2 is taken as the fuel particles. Based on Y. Cui's work [12], the theoretical fission gas swelling model for UO2 is built. The macroscale three-dimensional mechanical constitutive relations and stress update algorithms are established, with the grain-scale gas swelling model involved in the fuel particles and the irradiation creep considered in the matrix. Especially, the mid-value integration algorithm for the creep strain increments is adopted in the matrix. The formed userdefined subroutines [16,17] are validated and introduced to the finite element method (FEM) simulation of the multi-scale thermomechanical behavior within the CERCER pellet. The evolution rules are investigated and analyzed. The effect of the creep coefficient is examined.
2.1. Irradiation swelling in UO2 For the UO2 fuel particles in the CERCER pellet, the total swelling includes two contributions. One comes from the solid fission products, whose linear relationship with the burnup is given as [19] ΔV solid ¼ 2:5 10−29 F d V
ð1Þ
which is the volumetric variation relative to the original volume, that is the volumetric engineering strain; Fd depicts the fission density in fission/m3. The other swelling contribution results from the fission gas products, whose relation with burnup is of complexity. In this study, according to the theoretical model for the fission gas diffusion with the resolution effect and recrystallization effect involved [12], the theoretical model for the fission gas swelling is obtained and adopted.
2=15
Before recrystallization (F ≤ Fdx, F dx ¼ 4 1024 ðf Þ
, whereFdx is
the critical fission density; f is the fission rate in fissions/m3s), the fission gas swelling takes the form as the sum of the intragranular bubble swelling and the intergranular bubble swelling ΔV gas 4π 3 2πR3b C b ¼ : r b cb þ V r gr0 3
ð2Þ
3 In Eq. (2), 4π 3 r b cb depicts the intragranular bubble swelling, where r b
¼
bv nb 1=3 ð3hs4π Þ
¼
16π f n r g D0 c2g b0 nb
is the radius of the intragranular bubble [20], and cb
is the average concentration of intragranular bubbles. And
cgis the volume-averaged total gas content inside the grain, which is formulated as
cg ¼
−1 þ
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 64π f n r g D0 =b0 Y f t−3N=2rgr
32π f n r g D0 =b0
ð3Þ
where f is the fission rate, D0 =10−39m5is the proportionality constant
such that Dg ¼ D0 f , b0 = 2 × 10−23m3 is the proportionality constant
such that b ¼ b0 f , N is the concentration of the gas atoms in the grain boundary with the obtained solution as 2. Material properties For the considered CERCER pellet, the fission gas swelling model with consideration of recrystallization and gas atom resolution in grain scale is given for UO2 particles and the irradiated creep model for the MgO matrix is presented in this section. The other material properties used in the FEM simulation are the same as those in Ref. [15,18].
pffiffi pffiffi 8 2 3 2^ 2^ > es2 t erfc s2 ^t es1 t erfc s1 ^t ^ > r gr 5 ^ 2 2 2 β > > 4 ^ tbr gr =π t þ þ þ > > ^ ^ s2 ðs1 −s2 Þ s1 ðs2 −s1 Þ bλ < bλ ^ N t ¼ ð4Þ ! ^ 3 ^ 3 e−u21^t=r2gr > ^t 4βr > 2βr 1 gr gr > 2 2 ^ > tbr þ − =π > gr > 2 : 3 þ k1 r 2gr 5ð3 þ k1 Þ u21 u21 þ 3k1 þ k1
Fig. 1. (a) schematic figure of the fuel pellet with periodically-distributed particles (b) the selected RVE model (c) mesh grid.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
185
After recrystallization (F ≥Fdx), the original grain is divided into two regions including the recrystallized one and the unrecrystallized one. The total fission gas swelling is expressed as "
#
ΔV gas ΔV intra ΔV inter ΔV inter þ Vr ¼ ð1−V r Þ þ V V V V r gr r gr r grx
ð8Þ
where Vris the volume fraction of the recrystallized region; rgris the current grain radius excluding the recrystallized outside; rgrx = 0.1 μm is the average radius of the fine grains [29]. And the specified formulations for Vr and rgr can be found in Ref. [12]. In the recrystallized area, in consideration of the gas depletion [25, 26] taking place in the recrystallized area, the gas atom number on the grain boundary is obtained as [11] Fig. 2. Selected nodes in (a) the central particle and (b) the matrix.
^ are the expressions with the transforma^ and b where the parameters ^t, β Y f ^ ¼ b0 , the specific details can be found in ^ tion induced as ^t ¼ Dg t; β ¼ ; b
Nbx ¼
Y f t 3C bx =2r grx þ 1=8r 3grx
ð9Þ
Dg
Dg
Ref. [12]; b′ is the intergranular resolution rate with the employed value of 1×10−5s−1. Other parameters can be referred in Ref. [12]. Another item in Eq. (2)
2πR3b C b r gr0 ,
describes the intergranular bubble
swelling, where rgr0 is the original grain radius; Cbis the grainboundary bubble concentration, which is expressed as 0 qffiffiffiffi3 11=2 a dN 8z 12 ^ C b ¼ @ 2 dt A π ξδ
ð5Þ
3
a where 12 is the average atomic volume [28] with a is set as 5.47 × 10− 10m, z = 4 is the number of sites explored per gas-atom jump [11], ξ = 15is the grain-boundary diffusion enhancement factor [15], δ =2 × 10−9mis the width of the boundary [27], dN is the diffusion d^t
flux from grain interior to the grain boundary with the obtained form as
8 pffiffi pffiffi ^ 2β s2 s1 2^ 2^ > > ^tbr2 =π 2 > 1þ es2 t erfc s2 ^t þ es1 t erfc s1 ^t gr < ^ s1 −s2 s2 −s1 dN bλ ¼ 2^ 2 t=r −u ^ ^ > 2βr gr 4βrgr e 1 gr d^t > ^tNr2 =π 2 > − : gr 3 þ k1 u2 þ 3k1 þ k2 1
1
ð6Þ ^
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
In Eq. (4) and Eq. (6), s1 ¼ bλ 2 −
^2 λ2 b 4
^
^
þ rbλ , s2 ¼ bλ 2 þ gr
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^2 λ2 b 4
^
þ rbλ , u1 ¼ gr
^ gr . and k1 ¼ bλr 4:493− k14:883 þ3:612 Rb in Eq. (2) is the radius of the intergranular bubble, which satisfies the modified Van der Waals gas law [11], and neglecting the effect of hydrostatic pressure yields
2γ 4πR3b −hs bv Nb 3 Rb
where Cbx is the intergranular bubble density with the following expression as [11]
! ¼ Nb kT
ð7Þ
where γ = 1N/m is the surface tension [26], and Nb ¼ Nð^tÞ=C b is the gas atom number per intergranular bubble. It is known that Rbshould be obtained through the nonlinear iteration.
C bx
0 11=2 8zadN d^t r¼r grx C B ¼ @ 1=3 A : 12 π 2 ξδ
ð10Þ
Owing to the gas depletion, the fission gas atoms all migrate to the grain boundaries. As a consequence, the gaseous swelling is only contributed by the intergranular bubble swelling in the recrystallized area. According to Ref. [11], the volumetric gas swelling in the recrystallized region is obtained as !
ΔV intra 4πR3bx 3C bx 1 ¼ þ V 3 2r grx 8r 3grx r grx
ð11Þ
where Rbx is the radius of the intergranular bubble in the recrystallized area and it also satisfies the modified Van der Waals gas law as Eq. (7) with Nb replaced by Nbx and Rb replaced by Rbx. One can see that in the unrecrystallized part the gas swelling is still composed of the intragranular one and the intergranular one; while in the recrystallized part the gas swelling is only derived from the intergranular one because the gas depletion takes place in the recrystallized area [20] and almost all the generated gas atoms migrate to the grain boundaries directly. The detailed calculation formula is described in this section, and it needs to be pointed out that the intergranular bubble radius still should be obtained through the nonlinear iteration. 2.2. Irradiation creep in MgO Continuously attacked by fast neurons and other high energy fragments, irradiation creep occurs in the MgO matrix. The specific creep model for MgO hasn't been found in references. Here an irradiation creep model similar to CVD-SiC [21] is adopted as cr
ε ¼ c0 ψσ
ð12Þ
Table 1 Validation of the thermal constitutive relation for node M in the fuel particle. Time (day)
(K/mm) The temperature gradient ∂T ∂z
The heat flux qz obtained in ABAQUS (W/mm2)
The conductivity coefficient k used in ABAQUS (W/mm·K)
The conductivity from the model in Ref. [15] (W/mm·K)
The relative errors (%)
57.5 115 172.5
−68.14623 −82.50341 −96.03734
0.178414 0.1837 0.190923
0.002618 0.002227 0.001988
0.0026593 0.0022411 0.0020097
1.553 0.6291 1.080
186
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
Table 2 Validation of the thermal constitutive relation for node N in the matrix. Time (day)
The temperature gradient ∂T ∂z (K/mm)
The heat flux qz obtained in ABAQUS (W/mm2)
The conductivity coefficient k (W/mm·K)
The conductivity from the model in Ref. [15] (W/mm·K)
The relative errors (%)
57.5 115 172.5
−24.8468 −25.46282 −28.90954
0.313413 0.330782 0.354707
0.01261 0.012508 0.01227
0.012793 0.012699 0.012346
1.431 1.504 0.619
which is given based on the fact that a linear relationship exists between the irradiation creep and the Mises stress for ceramic materials [22]; cr
where ε is the equivalent irradiation creep strain rate(/s); c0is the creep coefficient in (MPa ⋅ n ⋅ cm−2)−1; σdenotes the equivalent Mises Stress in MPa and ψdepicts the fast neutron flux in n ⋅ cm−2s−1. In this study, the effects of different creep coefficients on the calculation results are investigated.
Then, the constitutive relation at time t + Δt can be obtained as
eðtþΔtÞ
eðtþΔt Þ
σ tþΔt ¼ 2GðT þ ΔT; t þ Δt Þε ij þ λðT þ ΔT; t þ Δt Þεkk δij ij h i h i eðt Þ eðt Þ ¼ 2GðT þ ΔT; t þ Δt Þ ε ij þ Δε eij þ λðT þ ΔT; t þ Δt Þ εkk þ Δεekk δij eðt Þ
¼ 2GðT þ ΔT; t þ Δt Þεij
þ 2GðT þ ΔT; t þ Δt ÞΔε eij
þ λðt þ Δt; T þ
þ λðT þ ΔT; t þ Δt ÞΔε ekk δij :
eðt Þ ΔT Þεkk δij
ð14Þ 3. Three-dimensional stress update algorithms To simulate the thermo-mechanical coupling behavior in ADS fuel pellets accurately, it is of importance to develop the three-dimensional irradiation-effect-involved mechanical constitutive relations. The whole calculation is divided into a number of time increments. The incremental mechanical constitutive relations and stress update algorithms for UO2 and MgO are established in this section. The userdefined subroutines to define the mechanical constitutive relations can be written in accordance with the algorithms, and ultimately be used in the finite element simulation. Being a large-deformation problem, the three-dimensional constitutive relations and stress-update algorithms are derived in a rotating coordinate system. Considering a time increment[t, t + Δt], the threedimensional incremental constitutive relation is first given. At each integration point, the relationship between the Cauchy
An incremental constitutive relation for an integration point yields Δσ ij ¼ σ ijtþΔt −σ tij ¼ 2GðT þ ΔT; t þ Δt ÞΔεeij þ λðT þ ΔT; t þ Δt ÞΔεekk δij eðt Þ
þ 2ΔGε ij
eðt Þ
ð15Þ
þ Δλεkk δij
where Δσij are the stress increments and Δεeij are the elastic logarithmic strain increments; ΔG and Δλ are the increments of the Lame coefficients. Based on the incremental constitutive relation, the threedimensional stress update algorithms for the fuel particles and the matrix are to be developed in the following.
eðtÞ
stresses σ tij and the elastic logarithmic strains εij at time t is given as eðt Þ
σ tij ¼ 2GðT; t Þεij
eðt Þ
þ λðT; t Þεkk δij
ð13Þ
where λ and Gare the Lame coefficients, which are related to both temperature and time.
3.1. The three-dimensional stress update algorithm for UO2 For the fuel particles, assuming that the total strain increments consist of the elastic ones, the thermal expansion ones, the irradiation swelling ones and the plastic ones, then the elastic incremental strains in a logarithmic form are expressed as p th sw Δε eij ¼ Δεtotal ij −Δε ij −Δε ij −Δε ij
Fig. 3. (a) Contour plot of the temperature in the pellet on the 115th day; (b) evolution of the temperature along Path 1.
ð16Þ
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
187
Fig. 4. (a) Contour plot of the swelling strains on the 230th day in different fuel particles and (b) the selected node P in Particle 1.
where Δεtotal denote the total strain increments; Δεth ij ij are the logarithmic thermal strain increments which can be derived as Δε th ½ ln ð1 þ α TþΔT ðT þ ΔT−T 0 ÞÞ− ln ð1 þ α T ðT−T 0 ÞÞδij ij ¼
1 þ α TþΔT ðT þ ΔT−T 0 Þ δij ¼ ln 1 þ α T ðT−T 0 Þ
ð17Þ
where αis known as the secant thermal expansion coefficient in K−1 and T0 is the reference temperature in K. The irradiation swelling strain incrementsΔεsw ij are expressed as
Δε sw ij ¼
ln ð1 þ SW ðt þ ΔtÞÞ swðt Þ δij −ε ij 3
ð18Þ
Fig. 5. (a) Validation of the volumetric swelling of node P; (b) evolution of the gas swelling in the unrecrystallized region; (b) evolution of the radius of the intergranular bubble.
188
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
Fig. 6. (a) The volumetric swelling distribution and (b) temperature distribution on the 115th day in the selected particles.
where SW(t +Δt) represents the volumetric swelling engineering strain ΔV
ΔV solid gas with SWðt þ ΔtÞ ¼ ΔV V ¼ V þ V , and the detailed expressions can be found in Section 2.1, which are obtained according to the grainscale theoretical analysis; εsw(t) denote the irradiation swelling logarithij mic strains at time t, which are known for the concerned time increment. As shown in Eq. (2), the intergranular swelling is dependent on the intergranular bubble radius. While the intergranular bubble radius Rb satisfies the nonlinear equation, as depicted in Eq. (7). Thus Newton– Raphson iteration method is applied to solve Eq. (7) written in a unified form as
!
4πR3b Δ 2γ −hs bv N b −N b kT ¼ 0 gðRb Þ ¼ 3 Rb
ð19Þ
where Rb denotes the intergranular bubble radius in the recrystallized or unrecrystallized region, and the definitions of the other variables can be found in Section 2.1. Eq. (19) is solved using a scheme of the Newton–Raphson (N–R) type, i.e. ∂g Rkb ¼ Rkb −g Rkb = ð20Þ Rkþ1 b ∂Rkb where ∂g Rkb ∂Rkb
¼ 16πγRkb =3 þ
2γhs bv Nb 2 : Rkb
ð21Þ
Once convergence is achieved, SW(t + Δt) can be worked out, then the irradiation swelling strain incrementsΔεsw ij can be obtained. For the integration points whose Mises stresses satisfy the yield condition, the contribution of the plastic strain increments should be considered. According to the plastic constitutive theory [23], we have dεpij ¼
3sij p dε 2σ
ð22Þ
where dεpij denotes differential plastic strain components; σ depicts the qffiffiffiffiffiffiffiffiffiffiffi Mises stress with σ ¼ 32 sij sij ; sij is the deviatoric stress component with sij ¼ σ ij − σ3kk δij ; dεp is the differential equivalent plastic strain. Backward Euler integration of Eq. (22) yields
Δεpij ¼
3sijtþΔt 2σ tþΔt
Δεp
ð23Þ
which implies that the plastic strain increments are related to the unknown stresses at time t+ Δt. Furthermore, the equivalent plastic strain increment Δεp needs to meet the consistent condition as described in Ref. [15]. For an integration point, firstly the plastic strain increments are supposed to be zero, the trial stresses at time t + Δtare obtained; the corresponding Mises stress is judged whether the yield condition is satisfied, and if not, the trial stresses are the real ones; if the yield condition is met, the plastic strain increments should be further included, and the
Fig. 7. Radius of intergranular bubbles in (a) unrecrystallized region and (b) recrystallized region on the 115th day.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
189
Fig. 8. Contour plot of the equivalent creep strains in the matrix on the (a) 11.5th day and (b) 115th day.
subsequent stress update algorithm is similar to the one presented in Ref. [17]. 3.2. The three-dimensional stress update algorithm for MgO As for the matrix, the corresponding total strain increments Δε-ij are constituted of the elastic strain ones Δεeij, the thermal expansion cr ones Δεth ij and the creep strain ones Δεij . Hence, the elastic strain increments can be derived as total
th cr Δε eij ¼ Δεtotal ij −Δε ij −Δε ij :
ð24Þ
Mid-value integration of Eq. (26) yields 3 stþΔt þ stij ij Δεcr Δεcr ij ¼ 2 σ tþΔt þ σ t
ð27Þ
whereΔε cr represents the equivalent creep strain increment. According to Eq. (12), the equivalent creep strain increment can be similarly derived as Δεcr ¼ c0 ψ
1 tþΔt σ þ σ t Δt: 2
ð28Þ
The thermal expansion logarithmic strain increments are demonstrated as
Substituting Eq. (28) into Eq. (27), the creep strain increments can be obtained as
Δε th ij ¼ ½ ln ð1 þ α TþΔT ðT þ ΔT−T 0 ÞÞ− ln ð1 þ α T ðT−T 0 ÞÞδij 1 þ α TþΔT ðT þ ΔT−T 0 Þ δij ¼ ln 1 þ α T ðT−T 0 Þ
Δεcr ij ¼ ð25Þ
3 c0 ψ sijtþΔt þ stij Δt: 4
ð29Þ
Combining Eq. (15) with Eq. (24), we have where α is the secant thermal expansion coefficient and T0 in K is the reference temperature. Based on the creep constitutive theory [23], the differential creep strains dεcr ij are expressed as
dεcr ij ¼
3sij cr dε 2σ cr
th cr σ ijtþΔt ¼ σ tij þ 2Gðt þ Δt; T þ ΔT Þ Δεtotal ij −Δε ij −Δε ij eðt Þ eðt Þ th þ λðt þ Δt; T þ ΔT Þ Δεtotal kk −Δε kk δij þ 2ΔGε ij þ Δλε kk δij :
ð30Þ
ð26Þ
where dε is the differential equivalent creep strain.
Firstly, the trial stresses σpr ij can be obtained assuming that no creep strain increments take place, namely, the components Δεcr ij in Eq. (24)
Fig. 9. (a) Validation of the creep algorithm in the matrix; and (b) evolution of the equivalent creep strains along Path 2.
190
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
Fig. 10. Contour plots of (a) the displacement magnitude and (b) y-directional displacement on the 230th day.
are set as zero. So, the corresponding Cauchy stresses σtij+ Δt can be expressed as 3 pr cr tþΔt σ ijtþΔt ¼ σ pr þ stij c0 ψΔt ij −2GΔε ij ¼ σ ij − 2 G sij
ð31Þ
where G represents the shear modulus at time t + Δt. The deviatoric part of Eq. (31) becomes 3 tþΔt þ stij c0 ψΔt: sijtþΔt ¼ spr ij − 2 G sij
ð32Þ
Then, the deviatoric stresses can be obtained as
sijtþΔt ¼
3 t spr ij − 2 sij Gc0 ψΔt 3 1 þ Gc0 ψΔt 2
ð33Þ
pr t+Δt tþΔt 1 pr where stþΔt ¼ σ ijtþΔt − 13 σ kk δij , spr =σpr kk. ij ¼ σ ij − 3 σ kk δij , σkk ij
can be updated as Therefore, the Cauchy stresses σt+Δt ij σ ijtþΔt ¼ sijtþΔt þ
σ pr kk δij : 3
ð34Þ
4. Finite element model A composite pellet containing UO2 particles and MgO matrix is taken as the research object. In this simulation, the geometric model is set up
with the assumption that the fuel particles are periodically-distributed along the axial direction and the circumferential one, as illustrated in Fig. 1(a). The pellet has a radius of 4.15 mm, the particles are supposed to be spherical with the radius of 200 μm and the volume fraction 6.23%. Considering the periodicity and symmetry, the Representative Volume Element (RVE) is built, as shown in Fig. 1(b). A pressure of 10 MPa and a constant temperature of 873 K are applied on the outside surface of the RVE. The symmetric boundary conditions are applied on the other surfaces of the RVE, with considering that the main deformation is along the radial direction of a pellet in the ADS fuel rod. Similar to the irradiation conditions in the experiment [7], the fission rate of fuel particles and the fast neutron flux are set as 2.5 × 1020 fission/m3s and 2.5 × 1015n/cm2s, respectively. Thus the corresponding heat generation rate of fuel particles is 8 W/mm3. As shown in Fig. 1(c), the RVE is discretized into 88,617 elements and 124,547 nodes with the element C3D10MT in ABAQUS, under which the convergent calculation results can be obtained. The special thermo-mechanical constitutive relations of fuel particles and the matrix are defined with the user-defined subroutines. 5. Results and discussions The finite element simulation of the multi-scale thermo-mechanical behavior evolution is successfully performed. The effectiveness of the user-defined subroutines should be confirmed. Lacking a series of definite experimental results, in this section some numerical results and the relative theoretical ones are compared in order to validate the correct definition of the material performances.
Fig. 11. (a) Evolution of the radial displacement along Path 3 with increasing burnup; (b) y-directional strain in the matrix on the 230th day.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
191
Fig. 12. Contour plot of the first principal stresses (a) at the initial stage; (b) on the 34.5th day;(c) on the 115th day and (d) the ones with no creep on the 115th day.
Besides, the evolution laws of the multi-scale thermo-mechanical behavior will be given and analyzed, and the effects of different creep coefficients in Eq. (12) will be investigated. 5.1. Temperature field The thermal flux and the temperature gradient comply with Fourier's law described as qz ¼ −k
∂T ∂z
ð35Þ
where K is the conductivity in W/m·K, which varies with temperature and time for the fuel particles, and for the matrix it is just
temperature-dependent; T depicts temperature in K; ∂T denotes the tem∂z perature gradient along the z-direction in Fig. 1(b); qzin W/mm2 represents the z-direction heat flux component. In order to confirm whether the thermal conductivity has been correctly defined, node M in the central particle and N in the matrix are selected to analyze their results, as shown in Fig. 2. In ABAQUS the heat flux of a node is supplied, while the temperature gradient hasn't been output. Then, the z-direction temperature gradient of node M is ob−T M2 , where TM1andTM2 are the temperature results of the tained asTzM1 M1 −zM2
neighboring nodes of node M in Fig. 2(a); zM1 and zM2 are the z coordinates. Simultaneously, the z-direction temperature gradient of node N can be obtained similarly. Thus, the fluxes and temperature gradients
Fig. 13. Evolution of the first principal stress along Path 2 with increasing burnup for (a) the case with creep and (b) the case without creep in the matrix.
192
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
Fig. 14. Equivalent creep stains along Path 2 on the172.5th day.
at different burnups are given in Table 1 and Table 2. Hence, the actually used thermal conductivity of nodes M and N can be calculated through Eq. (35), and their theoretical ones can be simultaneously obtained according to the model in Ref. [15]. The relative errors between them are slight, as depicted in Table 1 and Table 2, which demonstrates that the defined thermal constitutive relations of fuel particles and matrix are correct. The contour plot in Fig. 3(a) gives the temperature distribution on the 115th day. Thus, a representative path-Path 2 is selected to output the temperatures at different burnups, as depicted in Fig. 3(b). It can be found that (1) the temperature decreases along the radial direction due to the heat generation in the fuel particles; (2) the fuel particles have higher temperatures than the close-by matrix, which stems from the much lower thermal conductivity of the particles; (3) the peak temperature value locates at the center of the pellet and it goes up as the burnup increases, and the maximum value on the 230th day is about 118 K higher than that at the initial stage of burnup. Additionally, one can expect that larger volume fraction of fuel particles would bring much higher temperatures. In this case, a constant temperature condition is applied on the outer surface. If the surface temperature also increases with burnup, the maximum temperature will also be higher at a certain burnup. Furthermore, compared to our previous study [15], the maximum temperature is higher in this study, which is due to the higher heat generation rate adopted for the fuel particles. 5.2. Multi-scale results of volumetric swelling In this study, the used fission gas swelling model is obtained from our analytical solution of the governing equations for fission gas diffusion in grain scale, which involves the resolution effect and recrystallization effect. The details are demonstrated in Section 2.1 and its
Fig. 16. Evolution of the axial stresses along Path 2 considering creep in the matrix.
correctness has been validated [12]. As mentioned in Section 2.1 that the irradiation swelling calculation is complex, as a result the userdefined subroutine UMAT is designed for updating stresses. Firstly, the effectiveness of the introduced swelling is to be benchmarked. The contour plot of the swelling strains in different fuel particles is depicted in Fig. 4(a), which shows that the swelling strains are higher in the particles closer to the center of the pellet. Node P in Fig. 4 (b) is picked out to output the volumetric swelling strains at different burnup levels. The obtained volumetric swelling strains from UMAT are labeled as red points in Fig. 5 (a). At the same time, with the intragranular bubble radius and the intergranular bubble radius at different burnups, the theoretical results are also calculated out according to the swelling model in Section 2.1, as illustrated with the solid line in Fig. 5 (a). One can see that the theoretical results agree well with the ones obtained from UMAT, which identifies that the swelling strains are successfully introduced into the stress update. Meanwhile, the details of the swelling contributions are also given in Fig. 5 (a). It can be noticed that: (1) the solid swelling holds a proportional relationship with the irradiation time, and before the critical fission density Fdxits contribution is similar as the gas swelling. (2) After the critical fission density, the total swelling accelerates, with the gas swelling in the recrystallized region increasing promptly and the gas swelling contribution in the unrecrystallized part declining with burnup. At high burnup the gas swelling in the recrystallized region is far more than the solid swelling. This phenomenon attributes to the recrystallization process. At the critical fission density the recrystallization process begins, and the fine grains occupy the original large grain completely when the recrystallization process ends. For a very fine grain, the fission gas atoms can diffuse to the grain boundary quickly and in this study the atoms are assumed to reach the grain boundary directly without consuming time.
Fig. 15. Contour plot of the axial stresses (a) at the initial stage and (b) on the 115th day.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
193
Fig. 17. Evolution of the first principal stress along Path 2 with the creep coefficient in the matrix set as (a) 6.0×10−26(MPa n cm−2)−1; (b) 6.0×10−27(MPa n cm−2)−1; (c) 2.3×10−27 (MPa n cm−2)−1; (d) 2.0×10−27(MPa n cm−2)−1.
In addition, the evolution of the gas swelling in the unrecrystallized region is given in Fig. 5 (b). It reveals that the intragranular bubble swelling keeps a very small quantity all the time; and at the moment the recrystallization starts, the intergranular bubble swelling has a 43.7% share of the total swelling, while the intragranular swelling only account for ~0.93%. With the recrystallization proceeding, the amount of intergranular swelling of the unrecrystallization region is increasing as a whole, while the total gas swelling decreases gradually to zero due to the shrinked volume fraction. One can also find that at high burnup the gas swelling in the unrecrystallized region disappears (after the 218th day), which means that the original grain is completely replaced by the fine grains and the recrystallization process ends. For a better understanding, the intergranular bubble radius evolution is investigated for the unrecrystallized region and the recrystallized region, as shown in Fig. 5(c). It can be found that the intergranular bubble radius Rb in the unrecrystallized region first increases to 0.0287 μm and then decreases. Nevertheless, Rbx in the recrystallized region quickly jumps to a large magnitude after the recrystallization initiation and then increases with time all along. This phenomenon is consistent
The average irradiation swelling at 17.4% FIMA was obtained as 47 ± 17% from the post-irradiation-examination (PIE) [7]. Our simulation result falls in this range. This means that the established fission gas swelling model presents a good description of irradiation swelling of UO2 microspheres in the ADS fuel pellet. In Fig. 4(a), one can see that different swelling presents in the fuel particles, which majorly results from the temperature. The irradiation swelling results on the 172.5th day are depicted in Fig. 6(a). The fuel particle (Particle 6 in Fig. 4(a)) close to the outer surface of the fuel pellet has lower swelling than that in the pellet center (Particle 1), with their temperature having the similar distribution in Fig. 6(b). In accordance with the modified Van der Waals gas law as Eq. (7) in Section 2.1, high temperature will result in large intergranular bubble radius. To confirm this, the results of the intergranular bubble radii in the unrecrystallized and recrystallized regions are given in Fig. 7. It is depicted that the lowest intergranular bubble radius exists in Particle 6 with the lowest temperature, which leads to the lowest swelling. 5.3. Irradiation creep in the matrix
with the fact that the contribution of intergranular swelling V r ðΔVVinter Þjrgrx in Eq. (8) increases with time because both the volume fraction Vrof the recrystallized region and the relative bubble swelling ðΔVVinter ðRbx ÞÞ all rise. As demonstrated in Fig. 5(a), the peak swelling reaches approximately 43% on the 200th day (17.4% FIMA). The average irradiation swelling in this fuel particle is calculated out as approximately 41.3%.
The contour plots of equivalent creep strains in the matrix on the 11.5th day and the 115th day are shown in Fig. 8(a) and (b), respectively. Nodes on Path 2 are picked to output their equivalent creep strains in the time interval [0.003d, 0.006d], thus the equivalent creep strain increments can be calculated and displayed as a series of red points in Fig. 9(a). The theoretical results are also obtained through Eq. (28) with the Mises stresses at the starting time and terminal
194
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
time. It can be observed that the simulation results and the theoretical ones comply with each other, which verifies the correct definition of the creep model. Also, the evolutions of the equivalent creep strains along Path 2 are shown in Fig. 9(b). It can be found that the middle position exhibits larger creep strain and the value increases with time as a whole, and the peak value becomes 17.2% on the 172.5th day. Compared to the previous study [15], the magnitude of creep strains is larger, which is for the reason that the fission gas swelling considering the recrystallization effect is much larger in this study and it strengthens the mechanical interaction between the particles and the matrix. 5.4. Evolution of the radial deformation of the pellet In the conceptual design of ADS fuel elements, an original gap locates between the fuel pellet and the cladding in order to prevent them from an intense pellet-cladding mechanical interaction at higher burnup. In order to achieve the safety design of nuclear fuel elements, firstly, evolution of the radial deformation of ADS fuel pellets should be fully investigated. Contour plots of the displacement magnitude and the y-directional displacement on the 230th day are shown in Fig. 10. Path 3 in Fig. 10 (b) is selected to output its magnitude displacement in the Cartesian coordinate system, because in this path there is the maximum radial displacement which just equals to its y-directional component. As depicted in Fig. 11, (1) the displacements in the interfaces between the fuel particles and the matrix are larger than the ones in the adjacent positions, and the differences between the maximum one and the adjacent minimum one increase with time. This is due to the fact that the increased particle swelling enhances the mechanical interaction with the surrounding matrix, which results in appearance of larger compressive strains in the surrounding matrix; (2) relative large compressive strains occur in the cladding-pellet interface, as illustrated Fig. 11(b), which accords with the radial displacement distribution; (3) the radial displacement increases with irradiation time, and the pellet radius has an increase of ~ 78.83 μm (1.88%) on the 230th day, while the maximum displacement reaches about 89.16 μm in the vicinity of the pellet surface; the radial deformation is degraded by the constraint of the matrix; (3) compared to our previous work [15], the radial deformation of the pellet enlarges because in this study we adopt a modified swelling model considering the grain recrystallization. As a result, the original gap between the fuel pellet and the cladding ought to be appropriately designed, especially for the pellet with a large particle volume fraction. Numerical simulation of the radial deformation can supply a basis for optimal design of the ADS fuel element. 5.5. Evolution of the mechanical field in the matrix As for the matrix in the CERCER pellets, evolution of the first principal stresses is focused on, because the first principal stress is related to the failure of ceramic materials. The creep coefficient C0 in Eq. (12) was employed with a value of C0 = 2.1× 10−27(MPa n cm−2)−1. The obtained first principal stresses at different time points are, respectively, given in Fig. 12 together with the results without considering irradiation creep. Accordingly, to capture the evolution rules in detail, Path 2 in Fig. 5(a) is selected to output the first principal stresses at different burnup stages, as depicted in Fig. 13(a). It is denoted that (1) large magnitudes exist at the initial stage with the maximum at the lowest region of Path 2; (2) the first principal stresses all decrease heavily on the 11.5th day;(3) however, on the 57.5th day, the first principal stresses shows a different distribution character with the maximum value appearing at the two ends;(4) afterwards, the stresses increase with burnup as a whole, and on the 115th day the maximum first principal stresses in the vicinity of the two ends exceed the maximum value at the initial stage.
The evolution rule of the first principal stresses is induced by both effects of particle swelling and irradiation creep in the matrix. On one hand, the accumulated particle swelling enhances the interaction between the fuel particles and the matrix, which has a trend to increase the first principal stresses at the lowest region of Path 2. On the other hand, as presented in Fig. 14, the intensive mechanical interaction there results in high irradiation creep strains, and large creep strains can relax the stresses there. Both effects compete with each other, which will affect the evolution rules of the first principal stresses and the fracture pattern. Meanwhile, the first principal stresses without considering irradiation creep are also given in Fig. 13(b), which indicates that the location with the maximum value will always maintain and the enhanced irradiation swelling will lead to occurrence of extremely large stress (without considering cracking). Hence, it is revealed that the fission-induced creep effect in the MgO matrix plays an important role in its cracking pattern. For the ceramic materials, the maximum tensile stress would be the driving force to induce fracture. In order to explore the mechanism of the radial cracks discovered in the irradiation experiment [7], it is of significance to study the evolution rules of tensile stresses there. The contour plots of the axial stress component S33 are shown in Fig. 15. One can find that the axial stress components S33 at all points are negative at the initial stage, while on the 115th day S33 changes to be very large tensile ones nearby the axial symmetric plane across the fuel particles. The axial stress components S33 along Path 2 at different burnup stages are depicted in Fig. 16. It can be observed that (1) the axial stresses are compressive stresses at the initial stage; (2) the compressive stresses nearby the two ends are gradually evolved into tensile stresses and increase with burnup; nevertheless, the axial stresses at the lowest point remain compressive stresses and become larger at higher burnup. At the same time, one can find that the maximum tensile stress and the maximum first principle stress have the same value at the two ends (see Fig. 13(a) and Fig. 16), which means that the corresponding axial symmetric plane is the principal plane with the maximum tensile stresses. As a result, to determine the cracking by the maximum tensile stress, one can obtain that the radial cracks occur and propagate along this plane, which is consistent with the examined cracks in the irradiation experiment.
5.6. Effects of the creep coefficient In the previous sections, the results considering creep are calculated with the creep coefficient set as C0 = 2.1× 10−27(MPa n cm−2)−1, and it is an assumption of ours for lack of definite data. As Kim did [24], a proper parameter can be obtained through FE simulation if some experimental results can be used for reference. In this section, the effects of the creep coefficient on the mechanical behavior evolution are to be investigated. Hence, simulations are implemented using several C0 values, and the evolution of the first principal stresses along Path 2 are plotted in Fig. 17. One can find from Fig. 17 that the stress peaks locate at the lowest point at the initial stage; as time proceeds, the peaks gradually move toward the path ends. Especially, it can be observed from Fig. 17(a) that for the case withC0 = 6.0 × 10−26(MPa n cm−2)−1all the stresses decrease quickly to a very small value. Fig. 17(b) denotes that the stresses decrease with time, then increase and again decrease. When the creep coefficient is set as 2.3 × 10−27(MPa n cm−2)−1, as shown in Fig. 11 (c), after the initial stage the stresses will never exceed the original maximum. With C0 = 2.0 × 10−27(MPa n cm− 2)− 1 the maximum value at higher burnup exceeds the peaks at the initial stage, and it locates nearby the two ends. As a whole, a large creep coefficient strengthens the creep effect, thus the stresses will be heavily relaxed. If the radial cracking happens, the coefficient C0 = (2.0~ 2.1) × 10−27(MPa n cm−2)−1 will be suitable to explain the fracture mechanism with the obtained results.
Y. Zhao et al. / Materials and Design 89 (2016) 183–195
6. Conclusions For the CERCER pellet, the fission gas swelling model for UO2 is obtained with the resolution effect and recrystallization effect considered in the grain scale. The three-dimensional stress update algorithm is developed for the fuel particles, in which the irradiation swelling strains are obtained with N–R iterations. With the irradiation creep involved, using the mid-value integration to calculate the creep strain increments the stress update algorithm is also established for the matrix. Besides, the ABAQUS user subroutines are written and verified to define the thermo-mechanical properties. According to the irradiation conditions in the experiment, FE simulation of the multi-scale thermomechanical behavior is implemented. The conclusions can be drawn as follows. (1) Temperature decreases along the radial direction as a whole and the central particle has the highest temperature. The maximum temperature increases with burnup, and on the 230th day it reaches 1200 K. (2) Large volumetric swelling occurs in the particles at a higher burnup level. The recrystallization process speeds up the gas swelling. After the original grain is completely replaced by fine grains, the volumetric swelling takes on a linear relation with the burnup. (3) The recrystallization process can be captured with the intragranular swelling evolution. When the intragranular swelling disappears, the recrystallization process ends. The intergranular bubble radius in the uncrystallized zone will still grow after the recrystallization initiation, after reaching a maximum value it gradually reduces. While, the intergranular bubble radius in the crystallized zone will continuously increase with burnup, which leads to continuous increase of the total gaseous swelling. (4) The gaseous swelling in the fuel particles is different, which is mainly due to the different temperatures there. High temperature will induce large intergranular bubble radius, as a result high gaseous swelling appears in the fuel particles near the center of the pellet. (5) Fission-induced irradiation creep in the MgO matrix has an important effect on the mechanical behavior evolution. With a creep coefficient value of 2.1 × 10−27(MPa n cm−2)−1, at a certain level of burnup the maximum tensile stresses exist on the
195
axial symmetric plane, which will result in the radial cracks there. If the definite experimental results can be supplied, the suitable creep coefficient can be obtained by analyzing the FEM results.
Acknowledgments The authors thank the supports of the Natural Science Foundation of China (11172068, 91226101, 11272092), and the Foundation from Science and Technology on Reactor System Design Technology Laboratory. References [1] R.J.M. Konings, R. Conrad, G. Dassel, B.J. Pijlgroms, J. Somers, J. Nucl. Mater. 282 (2000) 159. [2] F. Delage, R. Belin, X.-N. Chen, et al., Energy Procedia 7 (2011) 303. [3] V. Sobolev, W. Uyttenhove, R. Therford, W. Maschek, J. Nucl. Mater. 414 (2011) 257. [4] D. Haas, A. Fernandez, D. Staicu, et al., Energy Convers. Manag. 49 (2008) 1928. [5] D. Haas, A. Fernandez, D. Staicu, et al., Energy Convers. Manag. 47 (2006) 2724. [6] E.A.C. Neeft, K. Bakker, R.P.C. Schram, et al., J. Nucl. Mater. 320 (2003) 106. [7] E.A.C. Neeft, K. Bakker, R.L. Belvroy, et al., J. Nucl. Mater. 317 (2003) 217. [8] J. Lamontagne, S. Bejaoui, K. Hanifi, et al., J. Nucl. Mater. 413 (2011) 137. [9] J. Rest, Compre Nucl Mater 3 (2012) 579. [10] S.E. Lemehov, V.P. Sobolev, M. Verwerft, J. Nucl. Mater. 416 (2011) 179. [11] J. Rest, J. Nucl. Mater. 346 (2005) 226. [12] Y. Cui, S.R. Ding, Z.T. Chen, Y.Z. Huo, J. Nucl. Mater. 457 (2015) 157. [13] Y. Cui, S.R. Ding, Y.Z. Huo, C.L. Wang, L. Yang, J. Nucl. Mater. 443 (2013) 570. [14] X.-N. Chen, A. Rineiski, W. Maschek, et al., Prog. Nucl. Energy 53 (2011) 855. [15] S.R. Ding, Y.M. Zhao, J.B. Wan, et al., J. Nucl. Mater. 442 (2013) 90. [16] X. Gong, Y. Zhao, S. Ding, Mech. Mater. 77 (2014) 14. [17] Y. Zhao, X. Gong, S. Ding, Y. Huo, Int. J. Mech. Sci. 81 (2014) 174. [18] MATPRO-09, A Handbook of Materials Properties for Use in the Analysis of Light Water Reactor Fuel Rod Behavior, USNRC TREENUREG-1005, 1976. [19] Y.S. Kim, G.L. Hofman, J. Nucl. Mater. 419 (2011) 291. [20] J. Spino, J. Rest, W. Goll, C.T. Walker, J. Nucl. Mater. 346 (2005) 131. [21] A. El-Azab, N.M. Ghoniem, Mech. Mater. 20 (1995) 291. [22] T.D. Gulden, J. Am. Ceram. Soc. 52 (1969) 585. [23] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons Ltd., UK, 2000. [24] Y.S. Kim, G.L. Hofman, J.S. Cheon, A.B. Robinson, D.M. Wachs, J. Nucl. Mater. 437 (2013) 37. [25] M.V. Speight, Nucl. Sci. Eng. 37 (1969) 180. [26] J. Rest, G.L. Hofman, J. Nucl. Mater. 277 (2000) 231. [27] J. Rest, G.L. Hofman, Y.S. Kim, J. Nucl. Mater. 385 (2009) 563. [28] J. Spino, D. Papaioannou, J. Nucl. Mater. 281 (2000) 146. [29] J. Rest, G.L. Hofman, J. Nucl. Mater. 210 (1994) 187.