A new method to simulate the micro-thermo-mechanical behaviors evolution in dispersion nuclear fuel elements

Mechanics of Materials 77 (2014) 14–27

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

A new method to simulate the micro-thermo-mechanical behaviors evolution in dispersion nuclear fuel elements Xin Gong a, Yunmei Zhao a, Shurong Ding a,b,⇑ a b

Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China 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 27 January 2014 Received in revised form 17 May 2014 Available online 24 June 2014 Keywords: Dispersion nuclear fuel element Incremental constitutive relations Stress update algorithm Irradiation effect UMAT UMATHT

a b s t r a c t The large-deformation constitutive relations and stress update algorithms in a corotational framework are formulated respectively for the fuel particles, the matrix and cladding in dispersion nuclear fuel elements undergoing irradiation, with the main irradiation-induced effects within them considered. Their specific consistent tangent stiffness moduli are also developed. Correspondingly, the user subroutines UMAT have been programmed for definition of their mechanical constitutive relations. Besides, the user subroutines UMATHT have been written to define their thermal constitutive relations, in which degradation of the thermal conductivity of fuel particles are involved. An efficient method is established for modeling the irradiation-induced micro-thermo-mechanical behaviors evolution in dispersion nuclear fuel elements. The developed methodology is validated with the simulation results of the thermo-mechanical behaviors in fuel elements under an assumed irradiation condition. This study lays a foundation for optimal design of dispersion fuel elements. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction A dispersion nuclear fuel element comprises of a dispersion nuclear fuel meat encapsulated in the alloy cladding, and in the composite meat the fissile fuel particles are dispersively distributed in the non-fissile matrix. Compared to the current nuclear fuel elements in the nuclear plants, the dispersion nuclear fuel elements (Holden, 1967; Neeft et al., 2003) can be subject to much higher burnup. They are widely used in the research and test reactors and have a promising prospect (Xu, 2003; Duyn, 2003; Lombardi et al., 2008) to be used in the advanced nuclear reactors and disposal of nuclear wastes. UMo/Al dispersion nuclear ⇑ Corresponding author at: Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China. Tel./fax: +86 021 65642587. E-mail address: [email protected] (S. Ding). http://dx.doi.org/10.1016/j.mechmat.2014.06.004 0167-6636/Ó 2014 Elsevier Ltd. All rights reserved.

fuel plates (Kim and Hofman, 2011) have the potential to displace the highly-enriched Uranium fuel elements in the research and test reactors, because of their high density and stable irradiation performances. PuO2/Zircaloy dispersion nuclear fuel rods (Duyn, 2003) were studied to be used in the commercial power water reactors in order to burn the weapons-grade or reactor-grade plutonium. Plate-type PuO2/Zircaloy dispersion nuclear fuel elements should also be an alternative and are waiting to be studied. Dispersion fuels containing the minor Actinides (Ding et al., 2013) are the dedicated fuels in the future ADS systems. The lifetime of dispersion nuclear fuel elements is directly related to the micro-structure of the fuel meat (Neeft et al., 2003; Duyn, 2003; Vatulin et al., 1999; Schram and Klaassen, 2007). It is of significance to investigate the effects of micro-structures of the fuel meat on the in-reactor performances of fuel elements. Owing to the fact that the irradiation experiments are of long time

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

with high cost and hard to be in situ examined, it is very necessary to develop the appropriate theoretical models and numerical simulation methods to assess the irradiation-induced behaviors and lifetime issues. Simulation research on the micro-scale irradiation behaviors in dispersion nuclear fuel elements was paid great attention (Meyer, 2001; Verwerft, 2007). In the extreme irradiation environment of nuclear reactors, a dispersion nuclear fuel element experiences complex irradiation behaviors (Duyn, 2003; Suzuki and Saitou, 2005): (1) in the fuel particles fission heat is generated; the produced solid and gaseous fission products result in irradiation swelling, which leads to the strengthened mechanical interaction between the fuel particles, the matrix and cladding so that large deformation occurs within them; simultaneously, the thermal conductivity of fuel particles is degraded with burnup; (2) especially, the metal matrix and cladding undergo the irradiation damage effects (Rodriguez et al., 1984; Rowcliffe et al., 2009; Rodriguez and Fellow, 2005) induced by the high-energy fission fragments and fast neutrons, including the irradiation-hardening, irradiation creep and irradiation growth effects for Zircaloy (MATPRO-09, 1976). It is most critical to build the special thermal and mechanical constitutive relations considering the synergic irradiation effects within the fuel particles, the metal matrix and cladding. For the in-reactor thermo-mechanical behaviors evolution in dispersion nuclear fuel elements, some new theoretical models and simulation methods have been found in the references Duyn (2003), Rest and Hofman (1999), Taboada et al. (2002), Saliba et al. (2003) and Marelle et al. (2004). In these studies, the constitutive relations for the nuclear fuels and the metal structural materials were simplified and should be improved with the irradiation effects involved. In our previous works, some simulation methods were developed. The irradiation hardening effect in the matrix and cladding and the irradiation growth effect in the cladding were gradually realized with the virtual temperature increase method (Ding et al., 2008, 2009; Wang et al., 2011; Jiang et al., 2011), where the irradiation-time dependant mechanical constitutive behaviors were regarded as the virtual temperature dependant behaviors. Because the real temperature variations after the initial stage were ignored, the simulation method were further improved (Gong et al., 2013) with several simple user subroutines in ABAQUS programmed for simulation of the irradiation swelling and degradation of the thermal conductivity in the fuel particles, irradiation hardening and creep effects in the matrix and cladding. For the heterogeneous irradiation condition (Kim and Hofman, 2011), the above simulation method is still inconvenient since the time–temperature-dependent Young’s modulus of the metal matrix and cladding has to be discretized with the defined field variables. As in the references Roh and Bae (2010), Nanthikesan and Shyam Sunder (1995) and Mashayekhi et al. (2007), the special constitutive relations should be better to be defined in the user subroutines UMAT. In this study, for the PuO2/Zircaloy dispersion nuclear fuel elements, the stress update algorithms and the consistent tangent stiffness moduli are derived out based

15

on the incremental constitutive relations for large deformation problems in the co-rotational reference frame (Andrade-Campos et al., 2006; Belytschko et al., 2000; Voyiadjis et al., 2006; Simo and Hughes, 1998; Ponthot, 2002), with the irradiation effects involved. Accordingly, the user subroutines UMAT are programmed to define the mechanical constitutive behaviors of the fuel particles, the matrix and cladding. And the user subroutines UMATTH are written to define their thermal constitutive behaviors, considering degradation of the thermal conductivity of fuel particles. The effectiveness of the subroutines is validated through the fem calculation results of an assumed plate-type element in the power water reactor. 2. The material models in an irradiation environment As the material performances of PuO2 are similar as the ones of UO2 (Duyn, 2003), the thermo-mechanical properties of UO2 are used to represent the ones of PuO2 in this study. 2.1. The thermo-mechanical properties for fuel particles 2.1.1. Thermal conductivity The thermal conductivity of UO2 improved by Lucuta et al. (1996) consists of five contributions and can be expressed as

K UO2 ¼ K 0  FD  FP  FM  FR

ð1Þ

where K0 is Harding’s expression for the thermal conductivity of unirradiated UO2; FD quantifies the effect of dissolved fission products; FP describes the effect of precipitated solid fission products; FM is the modified Maxwell factor for the effect of the pore and fission-gas bubbles; FR characterizes the effect of radiation damage.

k0 ¼

1 0:0375 þ 2:165  104 T " #   4:715  109 16361 exp  þ 2 T T 

FD ¼

1:09

þ B3:265

" #  0:0643 pffiffiffi 1 pffiffiffi T ar tan 1:09 0:0643 pffiffiffi þ pffiffiB T B B3:265

# " 0:019B 1   FP ¼ 1 þ 3  0:019B 1 þ exp 1200T 100

ð2Þ

ð3Þ



FM ¼

1P 1 þ ðs  1ÞP

FR ¼ 1 

0:2   1 þ exp T900 80

ð4Þ

ð5Þ

ð6Þ

where T represents the temperature in Kelvin, B is the burnup in at.%, P is the volume fraction of the pores and bubbles; s is the pore shape factor (with a value of 1.5 for spherical bubbles in the absence of other data).

16

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

2.1.2. Thermal expansion The engineering strain of thermal expansion (MATPRO09, 1976) is given as Eq. (7)

where k is the thermal conductivity in W/m K and T depicts temperature in K.

DL ðTÞ ¼ 3:0289  104 þ 8:4217  106 ðT  273Þ L0

2.2.2. The Young’s modulus and Poisson’s ratio (Fisher and Renken, 1964)

þ 2:1481  109 ðT  273Þ2

ð7Þ

where T depicts temperature in K ranged from 300 K to 1530 K. 2.1.3. Elastic modulus and Poisson’s ratio The Young’s modulus and Poisson’ ratio of UO2 fuels can be described as (MATPRO-09, 1976)

h i E ¼ 2:26  1011 1  1:131  104 ðT  273:15Þ ½1  2:62ð1  DÞ;

m ¼ 0:316

ð8Þ

where E is the elastic modulus in Pa, T is temperature in Kelvin and D is the theoretical density (92–98%). In this study, D is set as 95%. 2.1.4. Mises hardening rule The Mises hardening rule (MATPRO-09, 1976) of UO2 66:9  0:0397T þ ð520:0  0:0386TÞep ; T 6 1200  C r¼ ð9Þ 36:6  0:0144T þ ð139:5  0:0688TÞep ; T P 1200  C where r is the Mises equivalent stress in kg f/mm2. the equivalent plastic strain.

ep is

2.1.5. Irradiation swelling The total irradiation swelling SW (t, T) of the fuel pellet is composed of three parts including gas-bubble swelling, solid fission product swelling and the fission densification (MATPRO-09, 1976).

 gs   DV 16450 ¼ 439:6 exp  V T  100

ð10Þ

 ss DV ¼ 0:0025 V

ð11Þ

The above two expressions depict the volume fraction of gas-bubble swelling and solid fission product swelling per 1020 f=cm3 ; T is temperature in K.

 ds DV ¼ f0:0142½1  expð6:7943BUÞ þ 0:00893 V ½1  expð1:1434BUÞg  ADST

ð12Þ

Expression (12) depicts densification of UO2 due to attack of the fission fragments, where BU is burnup in MWD/kg UO2, ADST is the adjusting factor (=0.6). According to Eqs. (10)–(12), the total irradiation swelling SW(t, T) at a certain burnup can be obtained. 2.2. The thermo-mechanical properties for Zircaloy matrix and cladding

ð14Þ

m ¼ 0:3303 þ 8:376  105 ðT  273:15Þ

ð15Þ

where E is the Young’s modulus in Pa, T is the temperature in Kelvin and m is the Poisson’s ratio. Eq. (14) is the Young’s modulus without considering the effects of irradiation hardening, while in this work this expression is adjusted by a factor to predict the effect of the fast neutron fluence. The Young’s modulus (Fisher and Renken, 1964) under irradiation becomes E ¼ ½9:9  105  566:9  ðT  273:15Þ  9:8067  104 =k1

ð16Þ

  /t k1 ¼ 0:88 þ 0:12 exp  25 10

ð17Þ

where k1 is a dimensionless modification factor to allow for the effect of fast neutron fluence; /  t denotes the fast neutron fluence (n/m2) which represents the extent of irradiation damage. 2.2.3. Irradiation hardening plasticity model The strain-hardening curve of unirradiated Zircaloy is described as (Hagrman and Reyman, 1979)

r ¼ K en 



e_

m

103

ð18Þ

where r in Pa is the true stress, e is the true strain and e_ is the true strain rate. If e_ < 105 =s; set e_ ¼ 105 =s. The parameters K; n and m in Eq. (18) separately denote the strength coefficient, strain-hardening exponent and strain rate sensitivity exponent. Their expressions are given as K ¼ 1:17628  109



 þ T 4:54859  105 þ T 3:28185  103 þ 1:72752T

ð19Þ

n ¼ 9:49  102



 þ T 1:165  103 þ T 1:992  106 þ 9:588  1010 T

ð20Þ

m ¼ 0:02

ð21Þ

where T is the temperature in K ranged from 300 K to 730 K. Due to the irradiation hardening effect, the strain hardening exponent is described by multiplying the coefficient k2 given as

k2 ¼ 1:369 þ 0:032  1025 /  t

ð22Þ

The strength coefficient under irradiation is given by adding the value as

2.2.1. Thermal conductivity (MATPRO-09, 1976) k ¼ 7:51 þ 2:09  102 T  1:45  105 T 2 þ 7:67  109 T 3

E ¼ ½9:9  105  566:9  ðT  273:15Þ  9:8067  104

ð13Þ

k3 ¼ 5:54  1018 /  t

ð23Þ

17

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

In Eqs. (22) and (23), /  t is the fast neutron fluence in n/m2; / is the fast neutron flux, which is locationdependent due to the heterogeneous irradiation condition in this study. 2.2.4. Irradiation creep model Under the attack of the fast neutron flux, the metal cladding will experience irradiation creep. The creep strain rate is adopted as (MATPRO-09, 1976)

 e_ c ¼ K/ r þ BeCr expð10000=RTÞt0:5

2.2.5. Irradiation-induced growth model Irradiation growth is defined as different dimensional changes in different directions with a constant volume. As handled in the reference Wang et al. (2011), the irradiation growth effect is not considered in the matrix, and it is mainly taken into account in the cladding. The adopted growth model (Hagrman and Reyman, 1979) gives the form of engineering growth strains in the length and width direction:

DL ðXÞ ¼ A½expð240:8=TÞð/  tÞ1=2 ð1  3f x Þð1 þ 0:02CWÞ L

ð25Þ

DL DL ðYÞ ¼ c  ðXÞ L L

ð26Þ 1=2

where A ¼ 1:407  1016 ðn=m2 Þ ; T is the temperature in K ranged from 310 K to 630 K; / is the fast neutron flux in n/m2 s; and fx is the texture factor in the length direction with the value of 0.13; CW is the cold work induced in manufacturing process, which is ignored in this study. In addition, with the texture factor in the width direction fy adopted as 0.39 a coefficient c = (1  3fy)/ (1  3fx) in Eq. (26) is defined to determine the engineering strain in Y-direction. Thus, the fractional change in the thickness direction due to growth can be obtained based on the prediction of no volume variation induced by irradiation growth. The volume of the full-sized plate at time t and t + Dt can be expressed as

Vðt; TÞ ¼ LðXÞ  LðYÞ  LðZÞ Vðt þ Dt; T þ DTÞ ¼ ðLðXÞ þ DLðXÞÞ  ðLðYÞ þ DLðYÞÞ ð27Þ

such that:     Vðt þ Dt; T þ DTÞ DL DL DL ¼ 1 þ ðXÞ 1 þ ðYÞ 1 þ ðZÞ ¼ 1 Vðt; TÞ L L L

ð29Þ

The corresponding true-strain forms of the three principal growth strains can be transformed as



 DL ðXÞ ; L   DL eðZÞ ¼ ln 1 þ ðZÞ L

eðXÞ ¼ ln 1 þ



eðYÞ ¼ ln 1 þ

 DL ðYÞ ; L ð30Þ

ð24Þ

where e_ c is the creep strain rate and the constants K ¼ 5:129  1029 , B ¼ 7:252  102 , C ¼ 4:967  108 , R ¼ 1:987 ðcal=mol KÞ. T is the temperature in Kelvin and t is the irradiation time in the unit of second. r represents the equivalent stress (Pa) and / is the fast neutron flux (n=m2 s). It can be seen that the creep behavior is time-stress–temperature dependent.

 ðLðZÞ þ DLðZÞÞ

DL 1  1 ðZÞ ¼  L 1 þ DLL ðXÞ 1 þ DLL ðYÞ

ð28Þ

The expression of the fractional change in Z-direction can be calculated as:

3. The three-dimensional constitutive relations and stress update algorithms In order to precisely simulate the microthermo-mechanical coupling behaviors in dispersion nuclear fuel elements, it is critical to involve the special irradiation-induced material constitutive relations for fuel particles, the matrix and cladding in the finite element calculation. The user subroutine UMAT can be used in ABAQUS to define the complex mechanical behaviors of materials in an irradiation environment of a nuclear reactor. For such a nonlinear problem, the whole loading process should be divided into a number of increments. To get the converged values of the mechanical variables at the end of each increment, enough equilibrium iterations should be performed. In the user subroutine UMAT, the stresses and solution-dependent state variables at the end of the increment should be updated according to the incremental constitutive relation and stress update algorithm. Simultaneously, the material consistent tangent DrÞ stiffness moduli @ð should also be provided for the @ðDeÞ equilibrium iteration. In this section, the three-dimensional incremental mechanical constitutive relations and stress update algorithms for the fuel particles, the metal matrix and cladding are developed respectively in a corotational framework suitable for large deformation problems. In the next section, the corresponding consistent stiffness moduli will be derived. 3.1. Incremental constitutive relation for large deformation problems In an irradiation environment, the performances of materials are usually dependent on the irradiation time. The elastic constants are temperature–time-dependent. The constitutive relation should take this special character into account. Considering a typical time increment ½t; t þ Dt with the temperature of an integration point changing from T to T þ DT; the relation between the components of the CaueðtÞ chy stress rtij and the logarithmic elastic strain eij at time t in the corotational framework is given as eðtÞ rtij ¼ 2GðT; tÞeijeðtÞ þ kðt; TÞekk dij

ð31Þ

Similarly, the constitutive relation at time t + Dt can be derived as

18

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

DtÞ rijtþDt ¼ 2GðT þ DT; t þ DtÞeeðtþ ij eðtþDtÞ þ kðt þ Dt; T þ DTÞekk dij

1þDli ðtþDtÞ=li0 anisotropic with Degr dij (with no summaij ¼ ln 1þDl ðtÞ=l i

h i eðtÞ ¼ 2GðT þ DT; t þ DtÞ eij þ Deeij h i eðtÞ þ kðt þ Dt; T þ DTÞ ekk þ Deekk dij eðtÞ

¼ 2GðT þ DT; t þ DtÞeij

þ 2GðT þ DT; t þ DtÞDeeij eðtÞ

þ kðt þ Dt; T þ DTÞekk dij þ kðt þ Dt; T þ DTÞDeekk dij

ð32Þ

eðtþDtÞ ij

Dt where rtþ and e are the components of the Cauchy ij stress and logarithmic elastic tensors at time t + Dt in the current coordinate frame, and Deeij depicts the increments of logarithmic elastic strains; G(T + DT, t + Dt) and kðt þ Dt; T þ DTÞ are the Lame constants at the current temperature and time. Thus, the incremental constitutive relation for an integration point of a material undergoing large deformation can be expressed as

Drij ¼ rijtþDt  rtij ¼ 2½GðT þ DT; t þ DtÞ eðtÞ

 GðT; tÞeij þ 2GðT þ DT; t e ij

þ DtÞDe þ ½kðt þ Dt; T þ DTÞ

strain increments Depij can be expressed as

Depij ¼

þ DTÞDe þ DTÞDe

e kk dij

þ DtÞDe

eðtÞ kk dij

þ Dke

eðtÞ ij

þ 2DGe

¼ 2Gðt þ Dt; T þ DTÞDeeij

þ kðT þ DT; t þ DtÞDeekk dij þ Dr0ij eðtÞ ij

ð33Þ

eðtÞ kk dij

where Dr ¼ 2DGe þ Dke are the known quantities eðtÞ associated with the known logarithmic elastic strains eij . Based on the above developed incremental constitutive relation, the stress update algorithms for the fuel particles, the matrix and cladding will be built in the following. 0 ij

3.2. Stress update algorithm for the cladding For the cladding, the effects of thermal expansion, irradiation hardening plasticity, irradiation growth and irradiation creep are considered. The logarithmic elastic strain increments are denoted as gr p cr Deeij ¼ Detotal  Deth ij ij  Deij  Deij  Deij

ð34Þ

where Detotal are the known total strain increments. ij Deth are the logarithmic thermal expansion strain increij ments as follows

Deth ij ¼ ½lnð1 þ aTþDT ðT þ DT  T 0 ÞÞ  lnð1 þ aT ðT  T 0 ÞÞdij ¼ ln

1 þ aTþDT ðT þ DT  T 0 Þ dij 1 þ aT ðT  T 0 Þ

ð35Þ

where a is the secant thermal expansion coefficient, which can be temperature-dependent. Degr ij are the logarithmic irradiation growth strain increments. Different from the spherical thermal expansion strain tensor, the irradiation growth strain tensor is

ð36Þ

One can see from Eq. (36) that the plastic strain increments are associated with the unknown stresses at time t þ Dt with Depkk ¼ 0. Decr ij are the creep strain increments. According to the constitutive theory of creep, we have

Decr ij ¼

e kk dij ¼ 2Gðt þ Dt; T e ij þ kðT þ DT; t

3sijtþDt De  tþDt p 2r

where Dep is the equivalent plastic strain increment; sij are the components of the deviatoric stress tensor with qffiffiffiffiffiffiffiffiffiffiffi  is the Mises stress with r  ¼ 32 sij sij . sij ¼ rij  13 rkk dij ; r

eðtÞ

 kðt; TÞekk dij þ kðt þ Dt; T

i0

tion for the repeated index i) in the corotational framework. The relative irradiation growth in the i direction with respect to its original length Dli ðt þ DtÞ=li0 varies with the direction i. It should be noted that the irradiation growth effect will not result in the volumetric variation of the cladding, that is, Degr kk ¼ 0. Depij are the plastic strain increments. For the case with plastic strain increments during this time increment, they can be obtained according to the constitutive theory of plasticity. Backward Euler integration of small plastic 3s strain increments depij with depij ¼ 2rij dep , the finite plastic

Dt 3stþ  tþDt  ij  Decr r ; t þ Dt; T þ DT tþ D t  2r

ð37Þ

where Decr is the equivalent creep strain increment, which is related to the equivalent stress, the time and the temperature at the end of this increment; in Section 3.5 the theory will be given how to calculate the equivalent creep strain increment based on the equation of equivalent creep rate. Obviously, the creep strain increments also depend on the unknown stress state at time t þ Dt with Decr kk ¼ 0. As mentioned above that the plastic and creep strain increments are both dependent on the stresses waiting to be solved, it should develop the stress update algorithm to deal with this nonlinear problem. The algorithm is established as follows. Firstly, assuming that no plastic strains occur during this time increment, the elastic strain increments become gr cr Deeij ¼ Detotal  Deth ij ij  Deij  Deij

ð38Þ

Using Eq. (37), substituting Eq. (38) into Eq. (33) yields the trial stresses rpr1 at time t þ Dt as ij pr t rpr1 ij ¼ rij þ Drij ¼ 2Gðt þ Dt; T þ DTÞ

 gr cr  Deth  Detotal ij ij  Deij  Deij   th 0 þ kðT þ DT; t þ DtÞ Detotal kk  Dekk dij þ Drij

 gr ¼ 2Gðt þ Dt; T þ DTÞ Detotal  Deth ij ij  Deij 3Gðt þ Dt; T þ DTÞspr1  pr1  ij  ; t þ Dt; T þ DT Decr r r pr1   th 0 þ kðT þ DT; t þ DtÞ Detotal kk  Dekk dij þ Drij



3Gðt þ Dt; T þ DTÞspr1 ij ¼ rpr2 Decr ij   pr1 r  pr1  r ; t þ Dt; T þ DT

ð39Þ

19

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

where rpr2 are the trial stresses without considering both ij the plastic and creep strain increments with

 gr rpr2 ¼ 2Gðt þ Dt; Tþ DTÞ Deij  Deth þ kðT þ DT; tþ ij  Deij ij

 gr pr1 th 0 DtÞ Dekk  Dekk  Deij dij þ Drij ; sij are the deviatoric part of the trial stresses corresponding to

r

pr1 ij ;

 pr1 is the Mises stress and r

rpr1 ij .

The creep strain tensor being a deviatoric one, it has no contribution to the spherical stress components. The deviatoric stress components can be obtained from Eq. (39) as pr2 spr1 ij ¼ sij 

3Gðt þ Dt; T þ DTÞspr1  pr1  ij  ; t þ Dt; T þ DT Decr r r pr1 ð40Þ

Accordingly, we have 

1þ

  pr1  3Gðt þ Dt; T þ DTÞ pr2  ; t þ Dt; T þ DT spr1 Decr r ij ¼ sij pr1  r

Obviously,

spr1 ij

ð41Þ

can be determined if the Mises stress

r pr1 and the equivalent creep strain increment Decr can be obtained, which are coupled with each other. Multiplying by 1.5 times itself for two sides of Eq. (41) yields

   pr1  2 3 pr1 pr1 3Gðt þ Dt; T þ DTÞ ecr r  ; t þ Dt; T þ DT 1þ D s s 2 ij ij r pr1 3 ¼ spr2 spr2 2 ij ij

 tþDt   3Gðt þ Dt; T þ DTÞsijtþDt Decr r ; t þ Dt; T þ DT r tþDt tþDt  3Gðt þ Dt;T þ DTÞsij Dep  ð46Þ r tþDt

sijtþDt ¼ spr2 ij 

Similarly, after some manipulations the satisfied non tþDt , Decr and Dep can be obtained as linear equation by r



r tþDt þ 3Gðt þ Dt; T þ DTÞDecr r tþDt ; t þ Dt; T þ DT  pr2 þ 3Gðt þ Dt; T þ DTÞDep ¼ r

 ð47Þ

h i Dt  tþDt should locate on the new yield As the point etþ ;r p

Dt  tþDt should obey the following surface, that is, etþ and r p

equation as



r tþDt ¼ ry etp þ Dep ; t þ Dt; T þ DT



ð48Þ

Dt where, etp þ Dep ¼ etþ . Combining the relation between p tþ D t  Decr and r , the nonlinear equation in terms of the equivalent plastic strain increment Dep can be obtained as

 tþDt ðetp þ Dep Þþ3Gðt þ Dt;T þ DTÞDecr gðDep Þ ¼ r

  tþDt ðetp þ Dep Þ þ3Gðt þ Dt;T þ DTÞDep  r  pr2 ¼ 0  r

ð49Þ

ð42Þ

Thus, the nonlinear equation satisfied by the two  pr1 and Decr can be obtained as variables r



r pr1 þ 3Gðt þ Dt; T þ DTÞDecr r pr1 ; t þ Dt; T þ DT  pr2 ¼r





which can be solved by the Newton iteration method. The iterative scheme is



 @g g DepðkÞ þ dDepðkþ1Þ ¼ 0 @ Dep DeðkÞ

ð50Þ

ðkþ1Þ Depðkþ1Þ ¼ DeðkÞ p þ dep

ð51Þ

p

ð43Þ

Simultaneously, the relation between the equivalent  pr1 can be built according creep strain increment Decr and r to the time hardening or strain hardening creep theory  pr1 can be solved shown in Section 3.5. Then, Decr and r with the Newton iteration method with the iteration scheme as  pr1ðkþ1Þ ¼ dr

into consideration in the stress updating. Introducing the contribution of Depij into Eq. (38) and using Eq. (36), the derivation similar to the above leads to the following equation as

r pr2  r pr1ðkÞ  3Gðt þ Dt; T þ DTÞDecr r pr1ðkÞ ;t þ Dt; T þ DT



ecr 1 þ 3G@@rD pr1 r pr1ðkÞ

ð44Þ

@g ep @ D  tþDt @ D ecr @ r  tþDt @ D ep . @r

where

@  tþDt

¼ @rDep þ 3Gðt þ Dt; T þ DTÞ þ 3Gðt þ Dt; T þ DTÞ The convergent Dep can be obtained when

ðkþ1Þ dep ðkþ1Þ 6 1  106 . Next, the Mises stress r  tþDt can be Dep

solved out according to the irradiation hardening function. And the equivalent creep strain increment can be subseDt quently calculated. The deviatoric stress components stþ ij can be obtained on the base of Eq. (46). The stresses can be updated with

r pr1ðkþ1Þ ¼ r pr1ðkÞ þ dr pr1ðkþ1Þ

ð45Þ

ðkþ1Þ When drrðkþ1Þ 6 1  106 ; the converged trial Mises  pr1 can be obtained and Decr can be further stress r

  pr1 > ð1 þ 1  106 Þry etp ; t þ Dt; T þ DT calculated. If r

 with ry etp ; t þ Dt; T þ DT being the calculated yield stress according to the irradiation hardening plastic curve as Eq. (18), it indicates that new plastic strains engender during  pr1 and Decr are the real this time increment. Otherwise, r values at time t þ Dt; with which sijtþDt and

rijtþDt can be

obtained and the stress update process can be completed.

 pr1 > ð1 þ 1  106 Þry etp ; t þ Dt; For the case with r T þ DTÞ; the plastic strain increments Depij should be taken

Dt rtþ ¼ sijtþDt þ 13 r ij

pr2 kk dij

Dt rtþ ij

and the other

solution-dependent state variables can be updated. 3.3. Stress update algorithm for the matrix For the metal matrix, the effects of thermal expansion, irradiation hardening plasticity and irradiation creep are considered. Compared to the metal cladding, the irradiation growth effect is ignored. So the elastic strain increments can be expressed as p cr Deeij ¼ Detotal  Deth ij ij  Deij  Deij

ð52Þ

Since the corresponding stress update algorithm is similar to the one for the cladding material, it is unnecessary to be repeated here.

20

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

3.4. Stress update algorithm for fuel particles

For the case with constant temperature and stress, integration on both sides of Eq. (24) yields

For the fuel particles, the contributions of thermal expansion, plastic strain and irradiation swelling strain are taken into account. So the elastic strain increments can be given as

ecr ¼ 2K/ðr þ BeCr Þ expð10000=RTÞt 0:5

p sw Deeij ¼ Detotal  Deth ij ij  Deij  Deij

De

  1 ¼ ½lnð1 þ SWðt þ DtÞÞ  lnð1 þ SWðtÞÞ dij 3   1 1 þ SWðt þ DtÞ ¼ ln dij 3 1 þ SWðtÞ

ð55Þ

Substitution of Eq. (55) into Eq. (33) yields pr ij

r

 sw ¼r r ¼ 2Gðt þ Dt; T þ DTÞ Detotal  Deth ij ij  Deij

 th sw dij þ Dr0ij ð56Þ þ kðT þ DT;t þ DtÞ Detotal kk  Dekk  Deij t ij þ D

3.5.1. Time hardening Considering a typical time increment for t to t þ Dt, the obtained equivalent creep strain increment based on the time hardening theory is given as

 tþDt  tþDt þ BeC r exp ½10000=RðT þ DTÞ Decr ¼ 2K/tþDt r

 0:5  ðt þ DtÞ  t 0:5 ð58Þ

ð54Þ

One can see that the irradiation swelling only results in the volumetric variations of integration points. The irradiation swelling strain tensor is a spherical one. As is known that the plastic strain increment is related to the unknown stresses at time t þ Dt; the stresses in the fuel particles are updated as follows. Firstly, assuming there are no new plastic strains occurring in the fuel particles, the trial stresses rpr ij are calculated with the elastic strain increments as sw Deeij ¼ Detotal  Deth ij ij  Deij

Based on Eq. (57), two theories are given to obtain the equivalent creep strain increment for cases with temperature and stress variations.

ð53Þ

where Desw ij are the increments of the irradiation swelling strain, which can be obtained from the relative volumetric swelling of the fuel particles SW as sw ij

ð57Þ

pr ij

 pr can be calculated out Then, the trial Mises stress r qffiffiffiffiffiffiffiffiffiffiffiffiffiffi pr pr  pr ¼ 32 spr with r ij sij where sij are the trial deviatoric stress components. Determining whether the yield condition



  pr > 1 þ 1  106 ry etp ; T þ DT is met, described as r

ry ðetp ; T þ DTÞ is obtained with the current temperature T þ DT and the equivalent plastic strain etp at time t. where

If the above yield condition is met, it indicates that new plastic strain increments take place in this time increment. Otherwise, the trial stresses are the real ones, the stress updating is completed. For the case with the plastic strain increments Depij ; the stress update algorithm is similar to the one of the cladding material. Here it is not stated in detail. 3.5. Theory to solve the equivalent creep strain increment In the stress update algorithms for the cladding and matrix materials, the equivalent creep strain increment should be calculated for a certain Mises stress in the Newton iteration. It is critical to compute the equivalent creep strain increment according to the equivalent creep strain rate as Eq. (24). Generally, a model of creep strain rate is obtained with a constant uniaxial stress. For the case with stress variations during a time increment, two theories (Maximov et al., 2014) will be presented in this section.

3.5.2. Strain hardening Differing from the time hardening theory, the equivalent creep strain increment is calculated as the following with the strain hardening theory. For a typical time increment [t, t þ Dt], firstly, an equivalent time t is calculated according to the creep strain ecrðtÞ ¼ EQCR at time t as

 tþDt  tþDt þ BeC r exp½10000=RðT þ DTÞt 0:5 ¼ EQCR 2K/tþDt r 0 12 EQCR A t ¼ @

 tþDt  tþDt þ BeC r exp½10000=RðT þ DTÞ 2K/tþDt r

ð59Þ

Then the equivalent creep strain at t þ Dt is obtained as

 tþDt ecr ðt þ DtÞ ¼ 2K/tþDt r  tþDt þ BeC r exp½10000=RðT þ DTÞ



 tþDt 0:5 2  tþDt þ BeC r EQCR þ 2K/tþDt r  ðt þ DtÞ 0:5  exp½10000=RðT þ DTÞÞ2 Dt

ð60Þ

Thus, the equivalent creep strain increment can be computed as





 tþDt  tþDt þ BeC r Decr ¼ EQCR2 þ 2K/tþDt r 0:5  EQCR  exp½10000=RðT þ DTÞÞ2 Dt

ð61Þ

4. Consistent stiffness modulus As mentioned above, the material consistent tangent rÞ stiffness moduli C ¼ @ð@ðDeDtotal should also be provided for Þ

the equilibrium iteration in the nonlinear fem calculation. In this section, the derivation is presented as follows. According to the stress update algorithm for the cladding in Section 3, for the integration point with both plastic and creep strain increments, we have

Drij ¼ 2GDeeij þ kDeekk dij þ Drð0Þ ij

 gr th ¼ 2G Deij  Depij  Decr ij  Deij  Deij   ð0Þ þ k Dekk  Deth kk dij þ Drij 3spr2 ij De  pr2 p 2r pr2

 3sij ð0Þ th  2G  pr2 Decr  2G Degr  kDeth kk dij þ Drij ij þ Deij  2r ¼ 2GDeij þ kDekk dij  2G 

ð62Þ

21

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

where the Lame constants G and k are the ones at the end of the time increment. Differentiation on both sides of Eq. (62) yields

dðDrij Þ ¼ 2GdðDeij Þ þ kdðDekk Þdij  3Gd  3G

spr2 ij  pr2

r

spr2 ij  pr2

dðDep Þ  3Gd

r

spr2 ij  pr2

r

!

¼ 3G

Dep

pr2

Decr ð63Þ

Define A ¼ 3Gd r pr2 s þ3G rijpr2 dðDecr Þ; then Eq. (63) can be simply expressed as

dðDrij Þ ¼ 2GdðDeij Þ þ kdðDekk Þdij  ðA þ BÞ

ð64Þ

In the above equations, the trial Mises stress is given as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 pr2 s  spr2 ij . Differentiating on both sides of the 2 ij qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pr2  pr2 ¼ 32 spr2 equation r ij  sij , we have

r pr2 ¼

pr2

pr2 3 pr2 s  2dsij 3spr2 ds  pr2 Þ ¼ 2qij ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ¼ kl pr2kl dðr  pr2 2r 2 32 spr2 ij  sij

ð65Þ

Simultaneously, it can be obtained from Eq. (47) as

r pr2  3GDecr  3GDep ¼ r tþDt

ð66Þ

  tþDt ¼ ry etp þ Dep ; t þ Dt; T þ DT , which is the where r consistent condition. Differentiation on both sides of Eq. (66) yields

¼

 tþDt @ Decr @ r dðDep Þ  tþDt @ðDep Þ @r

 tþDt @r dðDep Þ @ðDep Þ

ð67Þ

1  pr2 dðDep Þ ¼ p p cr  dr 3G þ h þ 3Gh h pr2

3spr2 1 kl dskl p cr   pr2 2r 3G þ h þ 3Gh h p

p

cr  tþDt @r @ðD ep Þ ; h

ð68Þ

@ D ecr  tþDt , @r

where h ¼ ¼ which can be calculated out according to the hardening function of plasticity and Eq.  tþDt . (16) and (18) for certain Dep and r Using Eqs. (65) and (66), we have

! pr2

A ¼ 3Gd

3G

sij

r pr2

! pr2

sij

Dep þ 3Gd

 pr2  dðr  pr2 Þ  spr2 dðspr2 ij Þ  r ij  pr2 Þ ðr

r pr2

Decr

!

2

ðDep þ Decr Þ pr2

 pr2  s 3 spr2 dðspr2 r pr2  r tþDt kl kl Þ  r  tþDt  ij  dðspr2  pr2 ij Þ  r   pr2 2 r  pr2 r r pr2 r   r pr2  r tþDt 3 r pr2  r tþDt ¼  dðspr2 g g dðspr2 Þ ij Þ  2 ij kl kl r pr2 r pr2

¼

pr2

sij

where gij ¼ r pr2 . With Eq. (68) we have

pr2

ð69Þ

ð70Þ

Combination of Eq. (69) and (70) yields

AþB ¼

r pr2  r tþDt  dðspr2 ij Þ  pr2 r  tþDt  p r h 3  gij gkl dðspr2  þ p p cr kl Þ ð71Þ 2 r pr2 3G þ h þ 3Gh h

Substituting Eq. (71) into Eq. (64), we have

r pr2  r tþDt dðDrij Þ ¼ 2GdðDeij Þ þ kdðDekk Þdij   dðspr2 ij Þ  pr2 r  tþDt  p  3 r h Þ   gij gkl dðspr2 kl 2 r pr2 3G þ hp þ 3Ghp hcr r tþDt pr2 ¼ 2GdðDeij Þ þ kdðDekk Þdij  dðspr2 ij Þ þ  pr2 dðsij Þ r  tþDt p 3 r h  dðspr2   gij gkl  kl Þ ð72Þ 2 r pr2 3G þ hp þ 3Ghp hcr

  gr þ k Dekk  Deth Since rpr2 ¼ 2G Deij  Deth ij  Deij kk  ij 0 Degr ij Þdij þ Drij ; thus

drpr2 ij ¼ 2GdðDeij Þ þ kdðDekk Þdij

Combination of Eq. (65) and (67) yields

¼

spr2 spr2 @ðDecr Þ @ r  tþDt ij ep Þ þ 3G ij dð D dðDep Þ r pr2 r pr2 @ r tþDt @ðDep Þ

¼

 pr2  pr2 s s Dep þ3Gd rijpr2 Decr ; B ¼ 3G rijpr2 dðDep Þ

 pr2  3GdðDep Þ  3G dr

spr2 spr2 ij ij dðDep Þ þ 3G pr2 dðDecr Þ pr2 r r

p cr sij 3spr2 3G þ 3Gh h kl dskl p p cr  pr2    pr2 2r 3G þ h þ 3Gh h r p cr 3G þ 3Gh h 3 pr2 ¼ p p cr  gij gkl dðskl Þ 3G þ h þ 3Gh h 2

!

spr2 ij  3G pr2 dðDecr Þ   r pr2 sij  pr2

B ¼ 3G

ð73Þ

The deviatoric components of the differential trail stresses can be obtained as

1 pr2 pr2 dsij ¼ drpr2 ij  drkk dij ¼ 2GdðDeij Þ þ kdðDekk Þdij 3 2  GdðDekk Þdij  kdðDekk Þdij ¼ 2GdðDeij Þ 3 2  GdðDekk Þdij 3

ð74Þ

Substituting Eq. (74) into Eq. (72), we have

dðDrij Þ ¼ 2GdðDeij Þ þ kdðDekk Þdij  dðspr2 ij Þ

r tþDt 3 þ pr2 dðspr2 ij Þ  gij gkl 2 r  tþDt  p r h  dðspr2   kl Þ r pr2 3G þ hp þ 3Ghp hcr 2 ¼ kdðDekk Þdij þ GdðDekk Þdij 3   tþDt 2 r þ 2GdðDeij Þ  GdðDekk Þdij 3 r pr2  tþDt  p  3 r h  gij gkl   2 r pr2 3G þ hp þ 3Ghp hcr   2  2GdðDekl Þ  GdðDemm Þdkl 3 Finally, we have

ð75Þ

22

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

Fig. 1. Dispersion fuel plate and the representative volume element.

    tþDt 2 2 r r tþDt dðDrij Þ ¼ k þ G  G pr2 dðDekk Þdij þ 2G pr2 dðDeij Þ  3 3 r r  tþDt  p r h  dðDekl Þgij gkl  3G  r pr2 3G þ hp þ 3Ghp hcr ¼ k  dðDekk Þdij þ 2G  dðDeij Þ   p 3Gh ð76Þ  3G  p p cr  dðDekl Þgij gkl 3G þ h þ 3Gh h tþDt

where k ¼ K  23 G ; G ¼ G rr pr2 ; K ¼ k þ 23 G. And the consistent stiffness tensor considering plasticity and creep can be denoted as  pc   Depc ijkl ¼ Dijkl  Dijkl ¼ Dijkl ðG ; k Þ   p 3Gh  3G  p p cr gij gkl 3G þ h þ 3Gh h

ð77Þ

where Dijkl ðG ; k Þ ¼ G ðdik djl þ dil djk Þ þ k dij dkl . If the consistent stiffness tensor considering only plasticity can be derived out similarly as p  Dep ijkl ¼ Dijkl  Dijkl

 p  3Gh ¼ Dijkl ðG ; k Þ  3G  p gg h þ 3G ij kl ð78Þ

For the case considering only creep, the consistent stiffness tensor can be obtained as

  p   Dec ¼ D  D ¼ D ðG ; k Þ  3G  ijkl ijkl ijkl ijkl

5.1. The finite element model Similar as the used method in our previous studies, the representative volume element (RVE) can be generated, as illustrated in Fig. 1. For an assumed dispersion nuclear plate that the spherical particles are cubically and periodically arranged in the matrix as shown in Fig. 2(a), with the sizes along the length and width directions much larger than the one along the thickness, the RVE is selected as Fig. 2(b). In the fuel meat, the volume fraction of fuel particles is set as 20%, and the diameter of the spherical fuel particle is set as 0.1 mm; the thickness of the cladding is set as 0.4 mm. Finally, 1/8 fraction of the RVE is selected as the finite element model, shown as Fig. 2(c), allowing for the symmetry of RVE. The fission rate in the fuel particles is set as a constant with the value 1  1020 f/(m3 s), and the fast neutron flux in the matrix and cladding is set as 2  1018 n=m2 s which does not vary with burnup. The boundary conditions to determine the temperature and mechanical field are given as: (1) Convection boundary condition is applied on the upside surface Z = H/2: k@T=@n ¼ hðT  T f Þ; where the temperature of the periphery fluid T f is 573 K. And the heat transfer coefficient used is 0:02 W=mm2 K. And the pressure of coolant water 15 MPa is also applied on the upside surface.

 3G cr gij gkl 1 þ 3Gh ð79Þ

Thus, the consistent tangent stiffness moduli for the cladding, the matrix and the fuel particles can be obtained respectively according to Eqs. (77)–(79).

(a)

(b)

5. Validation of UMAT and UMATHT In Sections 3 and 4, the stress update algorithms for the fuel particles, the matrix and cladding, and the consistent tangent stiffness moduli considering irradiation effects are derived out respectively. Accordingly, the user subroutines UMAT are developed to simulate evolution of the thermo-mechanical behaviors in the dispersion nuclear fuel element in an irradiation environment. To ensure that the formed user subroutines are right and effective, validation is carried out through the numerical simulation results of a numerical example-an assumed PuO2/Zircaloy dispersion nuclear fuel plate in this section.

Matrix

X Cladding Z Parcle

(c)

Y

Fig. 2. (a) The sketch map (b) RVE (c) finite element model of the dispersion fuel plate.

23

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27 Table 1 Element and node information for the developed finite element model.

Mesh

Element type

Total number of elements

Total number of nodes

C3D8RT

23,448

27,727

Fig. 3. Distribution of the Mises stresses in the Zircaloy matrix and cladding at 30% FIMA.

Table 2 Comparison of the ABAQUS results and the theoretic ones for the stress components of Node N.

r11 r22 r33 r12 r13 r23

ABAQUS result (MPa)

Theoretical result (MPa)

Error (%)

616.421 617.246 423.584 6.909 56.518 55.585

616.403 617.229 423.59 6.908 56.511 55.578

0.00287 0.00275 0.0013 0.0115 0.0117 0.0118

(2) Except the upside surface Z = H/2, the other surfaces of the finite element models all obey: koT/on = 0. And the symmetric boundary condition is applied for the mechanical fields. (3) The continuous displacement conditions are met on the interfaces between the fuel particles and the matrix and the ones between the meat and the cladding.

where G and k are the Lame constants at the temperature T þ DT and time t þ Dt. According to the logarithmic elastic strain tensor simulated by ABAQUS, the theoretical Cauchy stress tensor corresponding to it can be calculated. By comparing the theoretical Cauchy stresses and the stresses calculated by ABAQUS, the elastic constitutive relation defined in UMAT can be verified. The distribution of Mises stresses in the Zircaloy matrix at 30% FIMA is as shown in Fig. 3. Point N in the matrix is chosen for validation of the elastic constitutive relation of the Zircaloy matrix and cladding. Table 2 shows the theoretical stress components calculated by Eq. (80) and the ABAQUS results. It indicates that the respectively obtained results fit well. So the definition of the elastic constitutive relation of the Zircaloy matrix and cladding in UMAT is right. For the fuel particles, the distribution of Mises stresses at 30% FIMA is as shown in Fig. 4. Point K in Fig. 4 is chosen for validation. Table 3 indicates that the ABAQUS results are consistent with the theoretical results. Thus, the elastic constitutive relation of fuel particles is correctly defined in UMAT.

5.3. Validation of thermal expansion According to the relation between the true strain and the engineering strain, the logarithmic thermal expansion strain induced by temperature variation can be written as:

eth ¼ lnð1 þ aT ðT þ DT  T 0 ÞÞ

ð81Þ

Thus, the theoretical logarithmic thermal expansion corresponding to temperature T þ DT can be calculated according to the temperature T þ DT and thermal

The converged mesh is selected as Fig. 2(c) and the element and node information are shown in Table 1. 5.2. Validation of the elastic constitutive relation The relation between the Cauchy stress and the logarithmic elastic strain at time t þ Dt in the corotational framework is as shown in Eq. (32), which can be written in the matrix form as 8 r11 9 > > > > > > > > > r22 > > > > > > > < =

2

6 6 6 6 r33 frij g ¼ ¼6 6 > > r 12 > > 6 > > > > 6 > > > > 4 r > > 13 > > : ;

r23

2G þ k k k k 2G þ k k k k 2G þ k G G G

38 e e11 9 > > > > > > 7> > > e22 > > 7> > > > 7> 7< e33 = 7 7> c > 7> 12 > > > > 7> > > > 5> > > c13 > > : ;

c23

ð80Þ

Fig. 4. Distribution of Mises stress in the fuel particles at 30% FIMA.

24

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

in Fig. 5(b). It can be seen that the results fit well with each other for each point in the path. So the thermal expansion of the fuel particles in UMAT is right. A path in Fig. 6(a) is chosen for validation of the thermal expansion in the Zircaloy matrix and cladding. Fig. 6(b) indicates that the results fit well for each point in the path. So the thermal expansion of the Zircaloy matrix and cladding in UMAT is right.

Table 3 Comparison of the ABAQUS results and the theoretic ones for the stress components of Node K.

r11 r22 r33 r12 r13 r23

ABAQUS result (MPa)

Theoretical result (MPa)

Error (%)

403.28 401.05 0.557 0.139 0.0890 0.0320

403.28 401.05 0.5579 0.1398 0.089 0.032

4.94e05 7.95e05 0.01359 0.06768 0.01674 0.01078

5.4. Validation of irradiation swelling As the irradiation swelling is an isotropic deformation, the relation between the logarithmic swelling strain and the volumetric swelling is:

expansion coefficient. By comparing the theoretical logarithmic thermal expansion strain with the ABAQUS result, the thermal expansion can be verified. For the fuel particles, a path in Fig. 5(a) is chosen for validation of the thermal expansion strain. Comparison of the theoretical results and the ABAQUS results is shown

1 3

esw ¼ lnð1 þ SWðBUÞÞ

ð82Þ

0.002850

umat results theoretical ones

Thermal Strain

0.002845

0.002840

0.002835

0.002830

0.002825 0.00

0.02

0.04

0.06

0.08

0.10

Z(mm)

(a)

(b)

Fig. 5. (a) Path in the fuel particles (b) comparison of the theoretic result and ABAQUS result.

0.00193

umat results theorotical ones

0.00192

Thermal Strain

0.00191 0.00190 0.00189 0.00188 0.00187 0.00186 0.00185 -0.1

0.0

0.1

0.2

0.3

0.4

Z(mm)

(a)

(b)

Fig. 6. (a) Path in the matrix and cladding (b) comparison of theoretic result and ABAQUS result.

0.5

0.6

0.7

25

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27 0.06 0.05 0.04

εsw

0.03 0.02 0.01

umat results theoretical ones

0.00 -0.01 0

5

10

15

20

25

and t. On the other hand, as in the UMAT back Euler integration is used in calculating the increment of creep strain, the theoretical result can be obtained through the temperature and Mises stress at t þ Dt according to Eq. (58) for time hardening theory and Eq. (61) for strain hardening theory. Through comparison of the two results, the irradiation creep in the matrix and cladding in the UMAT is verified. A path in the matrix is chosen for the comparison as shown in Fig. 8(a). And Fig. 8(b) indicates that the ABAQUS results agree well with the theoretical ones. Thus, the irradiation creep in the matrix and cladding is correctly simulated in UMAT.

30

BU(% FIMA)

5.6. Validation of irradiation growth

Fig. 7. Comparison of the theoretical result and ABAQUS result.

where SWðBUÞ is the relative volumetric swelling of the fuel particles at time t þ Dt. By comparing the theoretical swelling strain calculated by Eq. (82) with the swelling strain calculated by ABAQUS, validity of the irradiation swelling simulated through UMAT can be verified. Fig. 7 shows that the theoretical results are consistent with the ABAQUS results for each burnup. So, the irradiation swelling in the fuel particles in UMAT is right.

5.5. Validation of irradiation creep For a material point, the increment of equivalent creep strain can be obtained through the creep strain at t þ Dt

As mentioned in Section 3, the true irradiation growth strain in the co-rotational framework is

egrij ¼ lnð1 þ Dli ðtÞ=li0 Þdij

ð83Þ

where Dli ðtÞ=li0 is the relative growth in the direction of the principal axis i of the material at t. So, irradiation growth can be verified by comparing the theoretical growth strain calculated by Eq. (88) and that calculated in ABAQUS. A point at the interface of the matrix and cladding is chosen for comparison. Table 4 indicates that the theoretical results are the same as the ABAQUS results. So, irradiation growth in the cladding is correctly simulated in UMAT.

0.00012

umat results theoretical ones

0.00011

Δε cr

0.00010 0.00009 0.00008 0.00007 0.00006 0.00005

0.00

0.02

0.04

0.06

0.08

Z(mm)

(a)

(b)

Fig. 8. (a) path in the matrix (b) comparison of the theoretical result and ABAQUS result. Table 4 Comparison of the ABAQUS results and the theoretic ones for the irradiation growth components.

Length (X) direction Width (Y) direction Thickness (Z) direction

ABAQUS result

Theoretic result

Error (%)

0.00155427 0.000433589 0.00112068

0.001554 0.00043356 0.0011206

0.0046695 0.0049060 0.0044888

26

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

BU=0% BU=5% BU=10% BU=15% BU=20% BU=25% BU=30% fem results

1100

True stress (MPa)

1000 900 800 700 600 500 400 300 200 100 0 -100 0.00

0.05

0.10

0.15

0.20

True strain

(a)

(b)

Fig. 9. (a) Mises stress at 20% FIMA (b) the irradiation hardening curve and results of point M.

5.7. Validation of irradiation-hardening plasticity As there is no plastic strain occurring in the fuel particles and Zircaloy matrix and cladding when irradiation creep is taken into consideration, the case without creep is simulated in order to validate the effect of irradiationhardening plasticity. Fig. 9(a) shows the distribution of Mises stress in the matrix and cladding at 20% FIMA. Point M with a large value of Mises stress is chosen for validation. Fig 9(b) shows the points with the equivalent true strains as the transverse ordinates and the Mises stresses as the longitudinal ordinates together with the irradiation hardening curve of the

matrix material at the corresponding burnup. It can be found from Fig. 9(b) that the obtained stresses and strains obey the irradiation hardening plasticity model. It indicates that the irradiation-induced hardening plasticity in the matrix and cladding is correctly simulated in UMAT. As the plasticity model for the fuel particles is similarly simple. A point J in the fuel particles is chosen as shown in Fig. 10(a). The point formed with the obtained Mises stress and the equivalent plastic strain of node J together with the plasticity curve of the fuel particles at 30% FIMA are plotted in Fig. 10(b). One can observe that the ABAQUS result fits well with the plasticity model. So, the plasticity in the fuel particles in UMAT is correct.

1100

theoretic node J

Mises Stress(MPa)

1000 900 800 700 600 500 0.00

0.02

0.04

0.06

0.08

The equivalent plastic strains

(a) Fig. 10. (a) Mises stress at 20% FIMA, (b) the plastic curve and result of Node J.

(b)

0.10

X. Gong et al. / Mechanics of Materials 77 (2014) 14–27

6. Conclusion In this study, the incremental large-deformation constitutive relations and the stress update algorithms for the fuel particles, the matrix and cladding have been developed; and the corresponding consistent tangent stiffness moduli have been also derived out. Accordingly, the user subroutines UMAT and UMATHT have been programmed to define their thermo-mechanical constitutive relations in fem calculation. According to the numerically obtained thermal–mechanical results, the validity and effectiveness of UMAT are confirmed. This study lays a foundation for evaluation and optimal design of dispersion nuclear elements. Acknowledgements The authors thank for the supports of the Natural Science Foundation of China (11172068, 91226101, 11272092, 11072062, 91026005, 11304324), the Research Fund for the Doctoral Program of Higher Education of China (20110071110013), Strategic Priority Research Program’’ of the Chinese Academy of Sciences (Grant No. XDA03030100). References Andrade-Campos, A., Menezes, L.F., Teixeira-Dias, F., 2006. Numerical analysis of large deformation processes at elevated temperatures. Comput. Methods Appl. Mech. Eng. 195, 3947–3959. Belytschko, T., Liu, W.K., Moran, B., 2000. Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons Ltd., UK. Ding, S., Jiang, X., Huo, Y., Li, L., 2008. Reliability analysis of dispersion nuclear fuel elements. J. Nucl. Mater. 374, 453–460. Ding, S., Huo, Y., Li, L., 2009. Effects of fission heat and fuel swelling on the thermal–mechanical behaviors of dispersion fuel elements. Mech. Adv. Mater. Struct. 16 (7), 552–559. Ding, S., Zhao, Y., Wan, J., et al., 2013. Simulation of the irradiationinduced micro-thermo-mechanical behaviors evolution in ADS nuclear fuel pellets. J. Nucl. Mater. 442 (1–3), 90–99. Duyn, L.V., 2003. Evaluation of the mechanical behavior of a metal–matrix dispersion fuel for plutonium burning, Georgia Institute of Technology: A Thesis for the Degree Master of Science in Mechanical Engineering. Fisher, E.F., Renken, C.J., 1964. Single-crystal elastic moduli and the Hcptransformation in Ti, Zr and Hf. Phys. Rev. 135, 482–494. Gong, X., Ding, S., Zhao, Y., et al., 2013. Effects of irradiation hardening and creep on the thermo-mechanical behaviors in inert matrix fuel elements. Mech. Mater. 65, 110–123. Hagrman, D.L., Reyman, G.A., 1979. MATPRO-version 11, a handbook of materials properties for use in the analysis of light water reactor fuel rod behavior, NUREG/CR-0497, TREE-1280, Review 3. Holden, A.N., 1967. Dispersion Fuel Elements. Gordon and Breach Science Publishers Inc, New York. Jiang, Y., Wang, Q., Cui, Y., Huo, Y., Ding, S., 2011. Simulation of irradiation hardening of Zircaloy within plate-type dispersion nuclear fuel elements. J. Nucl. Mater. 413 (2), 76–89. Kim, Y.S., Hofman, G.L., 2011. Fission product induced swelling of U-Mo alloy fuel. J. Nucl. Mater. 419, 291–301. Lombardi, C., Luzzi, L., Padovani, E., Vettraino, F., 2008. Thoria and inert matrix fuels for a sustainable nuclear power. Prog. Nucl. Energy 50, 944–953. Lucuta, P.G., Matzke, H.S., Hastings, I.J., 1996. A pragmatic approach to modeling thermal conductivity of irradiated UO2 fuel: review and recommendations. J. Nucl. Mater. 232, 166–180.

27

Marelle, V., Huet, F., Lemoine, P., 2004. in: Proceedings of the Eighth International Topical Meeting on Research Reactor Fuel Management, München, Germany, 21–24 March. Mashayekhi, M., Ziaei-Rad, S., Parvizian, J., Niklewicz, J., Hadavinia, H., 2007. Ductile crack growth based on damage criterion: experimental and numerical studies. Mech. Mater. 39 (7), 623–636. MATPRO-09, 1976. A handbook of materials properties for use in the analysis of light water reactor fuel rod behavior, USNRC TREENUREG1005. Maximov, J.T., Duncheva, G.V., Anchev, A.P., Ichkova, M.D., 2014. Modeling of strain hardening and creep behavior of 2024T3 aluminium alloy at room and high temperatures. Comput. Mater. Sci. 83, 381–393. Meyer, M., 2001. Report: The US Program for the Development of Inertmatrix Fuel for Transmutation Systems, Argonne National LaboratoryWest. Nanthikesan, S., Shyam Sunder, S., 1995. Tensile cracks in polycrystalline ice under transient creep: Part I—numerical modeling. Mech. Mater. 21 (4), 265–279. Neeft, E.A.C., Bakker, K., Schram, R.P.C., et al., 2003. The EFTTRA-T3 irradiation experiment on inert matrix fuels. J. Nucl. Mater. 320, 106– 116. Ponthot, J.P., 2002. Unified stress update algorithms for the numerical simulation of large deformation elasto-plastic and elasto-viscoplastic processes. Int. J. Plast. 18, 91–126. Rest, J., Hofman, G., 1999. DART model for irradiation-induced swelling of uranium silicide dispersion fuel elements. Nucl. Technol. 126, 88– 101. Rodriguez, P., Fellow, R.R., 2006. What happens to structural materials in a nuclear reactor: the fascinating radiation enhanced materials science. In: Annual Conference-Indian Nuclear Society, INSAC 2005. Rodriguez, P., Krishnan, R., Sundaram, C.V., 1984. Radiation effects in nuclear reactor materials-correlation with structure. Bull. Mater. Sci. 6 (2), 339–367. Roh, J.H., Bae, J.S., 2010. Thermomechanical behaviors of Ni–Ti shape memory alloy ribbons and their numerical modeling. Mech. Mater. 42 (8), 757–773. Rowcliffe, A.F., Mansur, L.K., Hoelzer, D.T., Nanstad, R.K., 2009. Perspectives on radiation effects in nickel-base alloys for applications in advanced reactors. J. Nucl. Mater. 392, 341– 352. Saliba, R., Taboada, H., Moscarda, M.V., 2003. DART-TM: a thermalmechanical version of DART for LEU VHD dispersed and monolithic fuel analysis. In: 2003 International Meeting on Reduced Enrichment for Research and Test Reactors, Chicago, Illinois, October 5–10. Schram, R.P.C., Klaassen, F.C., 2007. Plutonium management with thorium-based fuels and inert matrix fuels in thermal reactor systems. Prog. Nucl. Energy 49, 617–622. Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. SpringerVerlag. Suzuki, M., Saitou, H., 2005. Light Water Reactor Fuel Analysis Code FEMAXI-6 (Ver. 1), JAEA-Data/Code. Taboada, H., Rest, J., Moscarda, M. V., Markiewicz, M., Estevez, E., 2002. DART-THERMAL: Development and Validation of a New Dart Code Improved Version Proceedings of the 24th intl. Mtg. on Reduced Enrichment for Research and Test Reactors, San Carlos de Bariloche, Argentina, November 3–8. Vatulin, A., Lysenko, V., Kostomarrov, V., Sirotin, V., 1999. Alternative versions of inert matrix fuel for the use of civil and weapons-grade plutonium in reactors. J. Nucl. Mater. 274, 135–138. Verwerft, M., 2007. Report: LWR-Deputy Project Presentation. . Voyiadjis, G.Z., Abu Al-Rub, R.K., Palazotto, A.N., 2006. On the small and finite deformation thermo-elasto-viscoplasticity theory for strain localization problems: algorithmic and computational aspects. Eur. J. Comput. Mech. 15 (7–8), 945–987. Wang, Q., Cui, Y., Huo, Y., Ding, S., 2011. Simulation of the coupling behaviors of particle and matrix irradiation swelling and cladding irradiation growth of plate-type dispersion nuclear fuel elements. Mech. Mater. 43, 222–241. Xu, Z., 2003. Design strategies for optimizing high burnup fuel in pressurized water reactors, a thesis for doctor’s degree in MIT.