A model for analysis of fuel behaviour in transients

Nuclear Engineering and Design 42 (1977) 319-355 © North-Holland Publishing Company

319

A MODEL FOR ANALYSIS OF FUEL BEHAVIOUR IN TRANSIENTS S. BANERJEE, H.J. BRIDGES, J,J.M. TOO Whiteshell Nuclear Research Establishment, Atomic Energy of Canada, Pinawa, Manitoba ROE 1LO, Canada and T.R. HSU Department o f Mechanical Engineering, University o f Manitoba, Winnipeg, Manitoba, Canada Received 3 January 1977

For many applications, analysis of fuel element behaviour must take non-linear thermal, and elasto-plastic effects into account. This is particularly true if the fuel undergoes large deformations and rapid temperature transients. To meet this need a multi-dimensional fuel model based on finite element stress and thermal analysis has been developed. The model is solved for the transient temperature distribution by a step-by-step time incremental procedure. The temperature is then introduced into the elasto-plastic analysis as a thermal load and stresses and deformations are calculated. A model for treatment of creep and a special element for the gap between fuel pellet and cladding is incorporated together with semi-empirical procedures for calculating fission gas release, fuel pellet to cladding heat transfer coefficients, etc. The fuel model has been compared with both analytical solutions and in-reactor experimental results. The observed and predicted results are in good agree4nent.

PART I. THEORETICAL BASIS 1. Introduction Nuclear fuel elements usually consist of uranium dioxide pellets encased in leak-tight zirconium alloy cladding. The pellets and cladding form the first two lines of containment for radioactive fission products. Their behaviour under normal operating and transient conditions such as load-following, start-up, shutdown, accident conditions, refuelling, must therefore be understood. As our understanding of fuel related phenomena has increased, models to predict fuel behaviour have begun to evolve on sound theoretical foundations. In particular, ther thermo-mechanical behaviour of fuel elements has received increasing attention in the work of Ma [1 ], Gittus [2], Levy and Wilkinson [3], Rashid [4], Royl et al. [5], Hsu et al. [6,7] and Too et al. [8]. The analysis of thermo-mechanical behaviour is complicated by the following factors, amongst others, which must be considered: (1) interactions between fuel pellets, and between pellets and cladding; (2) temperature dependent materials properties (because of the large temperature changes and gradients that can occur); (3) the effect of strain rate on materials properties (for fast transients); (4) loading-unloading-reloading cycles which may occur during transients; and (5) deformation and change in properties of fuel pellets and cladding under irradiation. This paper describes a finite element analysis involving quasi-coupled heat transfer, elasto-plastic stress analysis and creep analysis which forms the basic framework to predict fuel behaviour taking many of these effects into ac-

320

,9. Baner]ee et al. / Analysis o/fuel behaviour in transients

count. As the analysis handles many non-linear thermal and elasto-plastic effects, it is being thoroughly validated by comparison with a range of experiments. Our approach is to check each part first and then check the whole combined analysis. Until now the predictions agree well with experiment and closed form solutions; they will continue to be verified by being compared with more experimental data.

2. T h e o r e t i c a l

basis for fuel model

Thermal and stress analyses comprise the main theoretical structure of the fuel model. Both these analyses are based on the finite element approach. The usual procedure in finite element stress analysis is to divide a solid into a finite number of elements of regular shape, e.g. triangular, quadrilateral. The system of discrete elements linked at the node points then represents the original continuum. The equations of equilibrium and continuity at the node points, together with the appropriate stress-strain relationships, provide a set of equations. With the prescribed boundary conditions these equations can be solved to give the stresses in each element and displacements at every node.

) (a) Triangular torus element Z

~K_ _~K

zl

i

I I

i i

If

Ir

I

I

I

i

i I

Z2

I K

_zi

rO

._l _3_ _- -~L ~ ' z -"~ - , I

s

I

rr5

~r41 r2

I !

[ J

I I

I I

Irl ir3

(b)

,

I

I I

fm

Global coordinates of a trlan~ular element

fr

(c) Global coordinates of s quadrilateral element

Z

(d) L o c a l c o o r d i n a t e s

of a triangular

element

Fig. 1. S h a p e s a n d c o o r d i n a t e s y s t e m o f t o r u s e l e m e n t s ,

,.~ Baner]ee et al. / Analysis of fuel behaviour in transients

321

The same basic approach can be applied to thermal analysis except that the two-dimensional displacement field is now replaced by the temperature and its time derivative fields. Two types of elements (from those shown in fig. 1) were used in the work described here. These are the triangular and quadrilateral toroidal elements shown in fig. l(a). The quadrilateral element may be viewed as four smaller triangular elements as shown in fig. l(c).

2.1. Finite element thermal analysis 2.1.1. Thermal energy balance The present analysis is based on the framework proposed by Wilson and Nickell [9]. The temperature over an element is assumed to be a linear function of the spatial coordinates. The nodal temperatures can then be expressed by the matrix equation

Tin(r, t) = {bin(r)) r ( T ( t ) } ,

(I)

where

Tm(r,t), {T(t) } = temperatures in the element and at nodal points, respectively, r

= spatial coordinates, and

{bin (r) }

= interpolation function.

The curly brackets denote a column matrix. Square brackets will represent a square or rectangular matrix. The temperature gradients are

{Tm, l(r,t) } = {am(r) }T {T(t) }.

(2)

For a solid structure in thermal equilibrium the temperature field must be continuous between elements. This condition can be satisfied by minimizing the functional derived by Gurtin [10] M

= ~ f ~ ( P m C m T m T m + Ym, tg~]Tm, ] - 2 P m q m T m - 2 p r a C m T ~ T m ) (r,t) do m m=l l)m -

f

(Q~nniTm)

(r, t) CLam,

(3)

Sm

where Pro, Cm, k~] = mass density, specific heat and thermal conductivity of element material, M = total number of elements in the system, qm = heat generation in the element Um = volume of element, and nt = outward normal to the element surface Sin. Substituting eq. (1) and (2) into eq. (3) gives = 1 {T(t) }'r [C] {T(t) } - {T(t) }T [C] {T(0) } + ½{T(t) }T [K ] {T(t) } - {T(t) }T {Q(t) }, where M

[C] = ~ m=l

M

[cm],

[K] = ~ m=l

M

[Km],

{Q(t)}= ~ m=l

{Qm(t)}.

(4)

322

X Baner/ee et al.

/ Analysis o f fuel behaviour in transients

The element heat capacity matrix is

[Cm] = f

pro(r) {bin(r)} Cm(r) {bin(r)} T dora •

(s)

I)m

The element conductivity matrix is

[K m] = f

(6)

{am(r)}TKm(r) {am(O} dom •

om

The element heat flow is

(Qm(t)} = f pro(r)q(r, t) (bin(r)} T dora + f Q~n(r, t)ni {brn(r)}TdSm • Om

(7)

8m

The first variation of eq. (4) over T(t) for 0 < t < oo yields

[c] ({T(t)}- {T(0)})+ [gl {T(t)} = (O(t)}.

(8)

The physical interpretation of eq. (8) at any time, t is as follows: Rate at which heat ] Rate at which Rate at which heat is1 stored in elements [ + I flows from elements/ = externalheat ] [adjacent to the nodeI lenters the node I adjacent to the node /

2.1.2. Conductivity matrix The element conductivity matrix [K m ] given in eq. (6) for a triangular element may be expressed in terms of the interpolation function {bm(r, t)} in eq. (1) as

f (a{~} a'[~}T+ ar

[Km] =a K \ ~7

a{~) a{~}T)rdA ' a-7-- az

(9)

A

in which

or

{¢} = (bm(r,t)} ,

+ (bjbkr - akz

} ,

(1 O)

( b/r + a/z where A is the area of the triangular element and aj, ak, bk are the local coordinates of the nodes [see fig. l(d)]. The differentials are a{¢} t

. 1

a-S- -2,4 ((bj - bD b,, -b,.} and

a{~b}T - 1 az

2A

{(ak -- aj) --ak a/ } .

Assume K is constant over the triangular element. Then, substitution of the above expressions into eq. (9) and

S. Baner]ee et al. / Analysis o f fuel behaviour in transients

323

evaluation of the integral for one radian of the ring element yields Kr c

F(bj

- b~)~ +

[Kin] =-~-|Sym./L

- ai)

bk(b] - bk) -- ak(ak - ai)

- b l ( b I - bk + a/(ak - al)

b~ + a~

- b / b k - aka/

J ,

#+ 4 where rc is the radius of the center of gravity of the triangle. As described earlier,quadrilateralelements may be treated as four triangularelements, as shown in fig.I(c).The heat capacity term associated with the nodal point 5 is assumed to be distributedto the four adjacent nodes. The 5 X 5 element conductivity matrix is reduced to a 4 X 4 matrix by the assumption that there is no external heat flow at node 5.

Z 1.3. Heat capacity matrix The heat capacity matrix of the structure is given by the direct summation of the element capacity matrices. The element capacity matrix is given in eq. (5). Since the interpolation function {bin(r, t) } for the axisymmetric structure is defined by eq. (10), the element capacity matrix may be obtained by integrating eq. (5). However, a simpler method may be used for calculating this matrix. The method is based on the concept of lumping the heat capacity of the element at the four nodes of the element, a scheme which is similar to the lumped mass idealization in vibration problems. Assume that for a unit rate of change of temperature at node n(J'n = 1), only the material within half the distance to the next node on each side of the node contributes to storing heat in the body area (hatched) in fig. 2. If the volume of the ith element is I'm and the specific heat Cm and the density Pm of the material are constant, then the total heat capacity of the element is A m = ZrnPrnCm .

(11)

This quantity is lumped at the four nodes of the element such that the magnitude of the capacity at each node is proportional to its contribution to the volume. Fig. 3 shows the cross-sectionof a quadrilateraltoroidal element. The element may be divided into four quadrilateralsby connecting the mid-points of the oi)positesides.The c o m m o n vertex of all the quadrilateralsis located at a point with radius (12a)

r c = ~(r i + r/+ r~ + rt).

r¢L \7

rK

roK ro / " ' . ' " 4

~'." " \ i ,.

3

'~'d I"

Fig. 2. Contribution o f node i to heat capacity of an element.

rj

Ll

/K

r

Fig. 3. l~umping of heat capacity in a typical quadrilateral element.

S~ Banerjee et al. / Analysis or fuel behaviour in transients

324

The radius of the centroid of each quadrilateral is approximately

rci=~[2rl+½(r]+rt) + r c ] ,

rci=~[2r/+½(rk +rt)+rc],

(12b,c)

rck = ~ [2rg + ~(rj + rt) + rc],

rcl = ~ [2rl + ½(rk + ri) + re].

(12d,e)

The total heat capacity is lumped at each node and is proportional to the radius of the centroid of the quadrilateral sub-element which contains that node. The element heat capacity matrix given in eq. (5) has thus become the diagonal matrix rci 0 0 0 ] Pm Cmgm rci 0 0 . (13) [cm ] 4rc sym. rcg

J

rcI

2.1.4. Numerical integration For a solid structure in thermal equilibrium, eq. (8) must be satisfied. This equation, when expressed as

[C] { T } + [K] { T } = { Q } , represents a set of non-linear differential equations which may be solved by numerical integration techniques. For instances t and t + At, this equation may be expressed as follows:

[C] {:it} + [K]

{Tt} = {at},

{Tt+at} = {at+at},

[c'] {iPt+at} + [K']

where At is the time increment and [C'] and [K'] are the respective capacity and conductivity matrices at temperature Tt+at. For simplicity, these two matrices are assumed to be constant within every time increment as shown in fig. 4. The difference between these two expressions is then

[C] ( ( : i t + a t } -

{:it}) + [K]

({Tt+at}- {Tt}) = {Qt+at}- {Qt}.

Again, by assuming a linear variation of temperature within each time increment (fig. 4),

{:it+at} =

{rt+at}-

{rt}

(t + At) -- t

or {:it+at} - {:it} = ((Tt+at} -

{Tt})l At- {:it }.

Eq. (14) can thus be simplified to

At

T

I

I

I

I

I

I

I

I

I

I~'q t

I

~t

t "tAt

Fig. 4. Instantaneous values of [k] and [c].

(14)

S. Baner]ee et al. / Analysis of fuel behaviour in transients

325

where {ATt} = {Tt+At} - {Tt} ,

{AQt} = {at+At}- {Qt}.

Now, with {ikt} = ((Tt - Tt-At)/At }, the above equation may be expressed in the following form: [K*] {ATt} = {Q~'},

(15)

where [K*] = [C]/At + [K],

{Q~'} = ~

{ATt_At} + {AQt} .

(16a,b)

Eq. (15) can now be solved directly for the nodal temperature increments at the end of the time increment. The nodal temperatures at given instant t can thus be calculated. The matrix [K*] is a function of time and may be computed at selected intervals when the change of magnitude of K is significant. Since the solution of eq. (15) requires the temperature condition at the end of the previous time step, a starting procedure is required. The equation used for the first time increment can be obtained from eq. (8)

[Cl

[K*] {Tat} = {Qat} + ~

{To},

(17)

in which (T0) are the initial nodal temperatures.

2. 2. Non-linear elastic-plastic finite element analysis Nuclear fuel subjected to power transients may deform plastically. A description of the relationship between stress and strain under continuous loading beyond the elastic limit of a material will be stress and strain dependent, and is, therefore, termed non-linear. Various theories for such elastic-plastic material behaviour have been proposed by Yamada et al. [ 11, 12] and Zienkiewicz et al. [ 13]. Most of the resulting constitutive relations are expressed in differential form and must be solved step by step to allow for their inherent non-linearity. These equations require a statement of the work-hardening character of materials, i.e. whether the material is capable of sustaining stresses over and above the initial yield stress. This hardening character is commonly expressed in terms of a single scalar point function which is a function of the stresses and strains themselves. This function is usually referred to as the yield function or plastic potential. The following derivations are based on the plasticity theories of Yamada et al. [11] and Ueda and Yamakawa [14], which assume the existence of incremental stationarity of the potential energy of a system of finite elements (incremental force balance). Since the linear elastic, finite element program is also based on the energy principle, it can, therefore, be modified in a relatively straightforward manner to incorporate stepwise linearized plastic behaviour. For a selected class of problems (small to moderately large deformations), the published results are encouraging, so that adoption of the incremental equilibrium approach at this time is warranted. A generalization has been introduced by considering material properties to be dependent on both temperature and strain rate. The basic assumptions are (1) inertial effects are neglected, and (2) the thermal and stress fields are assumed to be quasi-coupled as the thermal and stress fields are both made dependent on the current deformed geometry and are solved incrementally.

2.2.1. Derivation of the elastic-plastic sn'ffness matrix Although a global virtual work formulation can be accommodated in the analysis, at present incremental stationarity of potential energy is applied to give force-displacement equilibrium in each loading step: {AQt) = [Ktl {aut},

(18)

326

S. Banerjee et al. /Analysis o f fuel behaviour in transients

where

(AQI) = (ASI) + ( A L l ) + ( A F t ) + ( A P t ) ,

(19)

with

[Kt] = incremental elastic-plastic stiffness matrix, (Aut) = deformation vector, (AS t) --- incremental concentrated nodal point load vector,

M

(AL l) = ~

m=l

{eLr~) = incremental load vector,

{~Qt) = overall incremental load vector, (AFI) = incremental body force load vector, and (APt) = incremental pressure load vector. System matrices and vectors are obtained by summation of elemental properties such as in

M [Ktl = ~

[Kin] e •

(20)

m=l The total strain at a given stress level must be obtained by integrating small loads, and stress and strain increments along the loading path. For a small thermo-mechanical load increment: (de) = total strain

(dee) + elastic strain at constant T, e-"

(deT) + thermal expansion from ref. temp.

(del-e) + strain due to temp. dependent mat. props.

(de~:e ) + s.train due to C - dependent mat. prop.

(dee). plastic strain

The elastic strain is of the functional form expressed in a generalized Hooke's law

{e) = (e)((o)T, [Ce]-I(T, e ) ) ,

(22)

where (o) = stress vector, T = temperature and [Ce] = elasticity matrix. ~: = rate of effective strain. The temperature and strain rate effect on the material properties is expressed implicity by [Ce]-I(T, e ) ; the ex. plicit dependence on temperature results from thermal expansion. Expansion of eq. (22) in differential form gives

ae

d(e) = ~ o

d(o)+

~

a(e)

(a [Ce] -1

d T + a [ C e ] _ 1~

aT

dT+

a [Ce]-I d e ) ae

'

(23)

From Hooke's Law, (o) = [Ce] (e) and from the expression for thermal strain, (e) = ( a ) T w h e r e (o~)= thermal expansion coefficient vector, the first two terms on the right-hand side of eq. (23) may be calculated and eq. (23) becomes d(e) = [ C e ] - l d ( o ) + (at) d T +

c - / a [Ce] -1

lO)~dT+

a [Ce]- 1.--'~

ae

oe].

(24)

For plastic strain a generalized plastic potential F is used to include strain rate and temperature dependence as well as the hardening parameter K as

F=F((o), (ev), K,T,~-).

(25)

S. Banetjee et al. /Analysis of fuel behaviour in transients

321

To define plastic behaviour uniquely, we use the von Mises yield criterion F = 0, or {aF/h’} = {u’} from the incompressibility condition, and the Prandtl-Reuss stress-strain relation {dep} = dX{a’}, where h is a proportionality constant and {u’} is the deviatoric stress tensor. By substituting u’ in the stress-strain relation, we obtain {den} = dX @F/h’}

.

W-3

Solving eqs. (21), (24) and (26) for the stress increment {do) = [C,] {de} - [C,] [{a} dT+ y

leads to the constitutive

{u} dT+ F

(27)

From eq. (25) and the von Mses yield criterion the following relationship

W=

g 1

equation

exists:

T {do}+~dKt~dTt~df=O, I

or expressing the hardening

in terms of plastic strains,

(28) By substituting

(aF/h’}

= {a’} and eq. (27) into eq. (28), the proportionality

factor dX is obtained

as

(29)

It can be shown that the plasticity [CPI = [Gl

g (

ET )I

terms in eq. (27) will contain the components

[Ccl IS,

(30)

I

with (31) Let the elastic-plasticity G,l

= [Gel - &I

then a new constitutive

{Au) = [c,,] Following

matrix be of the form, 3

equation

{Ae) - [C,,]

the derivation

(32) for the elastic-plastic

({u} AT t q

given in Yamada [ 121:

analysis can be written as

{u} AT+ 9

{a} A$

- [“I

~‘““‘($+

A;-) (33)

328

S. Baner/ee et aL / Analysis o f fuel behaviour in transients -

~2

O n, t

t

t

t

O~Ozz

X =

OrrOzO t

I

I

t2

t

F

I

t

t

t

Ozz O00

OttO00 , i

0020

Ozz OzO Ozz GOt,

OrrOOr

I

t

I

t

t

t

(34)

P~

OzO F

O00 GOt

Ozz Orz

OrrOr

I

000 OzO

t

OzO OOr

t

t

Go00rz

I

OzO Or,z

o'g, ?

t

OOrOrz

I~

Orz

in which G = shear modulus of elasticity and

aF/T ~o

j Ec~] L ~ j

where

(35)

'

1--I) 1 -

1 -

2v

p

2v

1-v 1 - 2v

/)

[Ce] = 2G 1 -

/)

2v

1 -

2v

sym. 1 1-v 1

-

(36)

2v

0

0

0

i-

0

0

0

o

0

0

0

0

0

The elasticity matrix shown above is used for a general isotropic solid in cartesian coordinates. Zienkiewicz et al. [13] have shown that

{aF} ~ = ' T '-~ aF o =31~, 2 ' aep - {a' }r, where H' is the slope of the effective stress (0-)/effective strain (~-) curve. Thus,

Since

,-' . . . . .

O/~OrOzOOLG=l

T 2~2 =

it follows that

Substituting eqs.(35) and (38) into eq. (31), we get

S = ~-(3G + H ' ) ~ -2 . The plasticity matrix [Cp] may be calculated by substituting eqs. (34) and (39) into eq. (30) leading to [Cp] :

2G" (9G/2~ -2) [gz~zm-~,l 3a + tt'

-

2G [[7~rr,,~ ] So '

(39)

(40)

S. Baner]ee et al./ Analysisof fuel behaviour in transients

329

where [ ~ ] is the symmetric deviatoric stress matrix given in eq. (34). The elastic-plasticity matrix in eq. (32) may thus be written as

2p

1 -

So

t t

P 1 -

OrrOzz

2p

S o t

P 1 - 2p #

[Cep] = 2G

#

orroo0 SO

t

OrrOzO

1 -- P

Ozz

2p

So

1 -

#

OrrOOr

t

I

OrrOrz

So

I

OzzOO0 So

1 - 2v t

1- p 1 - 2p #

Ozz OzO

t

t

Ozz OOr

t

GO0 O0 r

So

• (41)

2

So #

t

o So

GO0Ozo

So t

So t

sym.

I

v

So t

t2

So t

1

So

2

So l

t

OzzOrz

So

t

#

OzO O0 r

t

#

Go0 Orz

OzO Orz

So

So

t2

Oor So t ;

OorO~ So

1 2

So

This form is similar to that of Yamada et al. [11 ]. 2. 2. 2. Material constitutive equations

The slope H ' of the effective stress/effective plastic strain curve may be obtained as H'-

do_ 1 dep 1/E T - 1 / E '

(42)

where ep, O = effective plastic strain and stress respectively, E T = tangent modulus, and E = Young's modulus. I-Isu and Bertel [15] have suggested that the experimentally determined constitutive relations be approximated by a generalized family of continuous functions based on the proposal of Richards and Blacklock [16]. This leads to ~- = E~-/{1 + (E~-/[(1 - E ' l E ) o k

+ g,~-l)n }l/n,

where E ' = lime__,: E, crk = auxiliary effective stress close to elastic to plastic transition and n determines the absorptions of the transition. The relation for E T can then be deduced as d°=E/1

ET = d~

t +I(1 -- E'IE)-01, E~" + E'g]

--~

1 + (1 - E'/E)-~ ~ + E'~-

Temperature and strain rate dependencies can be taken into account by specifying E, E', ~k and n as functions of T and ~-. This relationship is convenient from the computational viewpoint.

2.2.3. Thermo-elastic-plastic stress analysis with kinematic hardening When a structure is subjected to repeated loading, the usual finite element elasto-plastic analysis based on isotropic hardening of materials must be modified to allow for the Bauschinger effect. One possible model is Ziegler's modification of Prager's kinematic hardening theory [ 17]. An illustration of kinematic hardening in conjunction with the yon Mises yield criterion for a biaxial stress rate is shown in fig. 5. As can be seen from this figure, the loading curve shifts by an amount atl upon a change of stress state from 1 to 2. This hardening scheme has been

330

S. Baner]ee et aL / Analysis o/fuel behaviour in transients o'22 (x22

~ ~

~ 2

~- Crll

y

/" AT

o.i

o'11 T Fig. 6. Yield surfaces as function of temperature.

Fig. 5. Shift of yield surfaces during kinematic hardening,

used by Armen et al. [18] in conjunction with the finite element method for structures subjected to cyclic loading. For finite element analysis of structures subject to strong thermal loading the temperature dependence of material properties must be included. Therefore the yield function and flow law become temperature dependent. The change of a yon Mises-type yield surface due to a temperature increment AT is illustrated in fig. 6. The yield surface expands (or contracts) with changes in temperature. This has been included in the elasto-plastic analyses described previously. For structures subjected to combined cyclic thermo-mechanical loading, both kinematic hardening and temperature effects must be included. We consider the loading path from stress state 1 to state 2 to be composed of paths from 1 to 1' for a finite temperature increment ATand from 1' to 2 due to kinematic hardening. This is illustrated in fig. 7. The flow rule which includes kinematic hardening and temperature and strain rate dependence then has a form similar to that for isotropic hardening

de/~.= dX (aF/aoi/),

~a

o"22

AT

Z -NI

Fig. 7. General plastic loading path.

S~ Baner/ee et aL / Analysis or fuel behaviour in transients

331

where the plastic potential F is modified to accommodate the shift of the yield surface

F = r(o0, at/, T, e), in which aii is the yield surface translation tensor shown in fig. 5. The plastic multiplier dX during kinematic hardening based on Ziegler's model [17] may be shown to be 0F v + ~ OF OF de dX = Ooi/ [Co~tdekl- Cil~l([3ktdT + DkldT Die de)] + ~ d T + - 3 [C(T, ~b-~//* + C07a~-~t*] (OF/Ooii) where the coefficient C(T, e ) in the above expression is the equivalent hardening coefficient for a multi-axial state of stress as derived in Isakson et al. [18], and the following translated quantities are defined by Hsu and Bertels

[2o]: translated stress tensor:

o//* = o i / - ati,

translated stress deviators:

Oil* = oi/* - -~Okk 6i/,

where 6ij = Kronecker delta. The constitutive equation for the incremental stress and strain components can be shown to be

(Ao) = [Cep ] ( A e ) - [Cep I (~)AT+

(DT)AT+(Dg) ae)

[Cels(~-*)/aFI,.3_fAT+_~aF,,,~_).

(43)

with

[Cep] = [Ce]-

9 [Ce]

,[~'*]~{~-*].T[Ce]

, s = 9(~-*){~-*) x C(T) + 9{~-*)T[ce] {~-*),

The main difference between eqs. (33) and (43) is that the latter involves the translation of the yield surface during the cyclic loads. Both types of stress-strain relationship are similar in form to the usual elastic constitutive relation. Thus, the standard finite element formulation for elastic analysis is equally applicable to the elasto-plastic analysis in the incremental sense, i.e. [K(o)] {Au) = (AQ).

(44)

The stiffness matrix is now related to the elasto-plastic matrix [Cep] which is in fact dependent on the current stress state. The stiffness matrix, of course, will have to be updated at each load increment.

2. 2.4. Numerical solution to elasto-plastic analysis The most commonly used solution scheme for the non-linear system of eqs. (44) is Euler's method. It advances the current state by incrementing the load and updating the stiffness at each load step. The application of this method is straightforward; however, the linear incremental material law intrinsic to the method leads the solution to step out from the true solution. Unless very small increments are used, errors usually accumulate throughout the loading history. Alternative methods have been proposed by various authors (notably initial stress-strain methods by Zienkiewicz et al. [ 13] which use constant stiffness, at the cost of iteration instead of updating the stiffness). These methods are more efficient in some applications. In the present analysis, a mixed variable stiffness method and one-step correction method is used. The algorithm of this method is stated in the following:

332

S. Baner/ee et al. /Analysis or fuel behaviour in transients

(1) for a given set of values of displacement Uo, strain eo, incremental load'AQ0 and the stiffness [Ko(oo)], the displacement u 1 and strain el can be calculated from eqs. (44) and (21)

(AUl) = [go(oo)] -1 (AQo) ; (2) (01} is determined by eq. (33) or eq. (43) and the unbalanced load is evaluated by (AR 1) = {AQo) = (AQ 0) - f [ B ] T (Ao1)do, incrementing load vector (AQ1) where [B] is the matrix relating strain to displacement; (3) current stiffnesses [K 1(o l)] are obtained based on the new (o); (4) the displacements are obtained from (Au 2 ) = [K1 (Ol)]-1 ((AQ 1 ) + (ARI)); and (5) steps (1)-(4) can be repeated throughout the loading process with {Aui) = [Ki-1 (°i-1 )]-1 ((AQi_ 1) + (ARi-'I ) ) .

2.3. Creep analysis Strictly speaking, the thermo-elasto-plastic stress analysis described above is valid only if the fuel element is under thermal-mechanical load for a short time at a low temperature range. The fuel is exposed to high temperature, large thermal gradients, long in-reactor residence time, etc. Under these conditions time dependent material behaviour becomes important, and creep deformation can no longer be ignored. Consideration of creep effects must therefore be included. In recent years creep of nuclear reactor structures has been discussed by Zudans et al. [21], Rashid et al. [22], Gittus [23] and Krempl [24]. The approach used in this paper is patterned after the theory of plasticity. However, creep deformation is analysed separately and treated as a component of the total deformation. Due to the non-linear dependence of creep strains on stresses, exact duplication of an actual creep test curve is often quite difficult, even for a one-dimensionally stressed situation. Multi-dimensional creep tests are usually expensive. Thus, most creep data are obtained from uniaxial tests. For multi-axial creep, it is necessary to generalize these data to include the multi-dimensional stress state. To do this, four basic requirements have to be established: (1) a creep law based on uniaxial creep test data; (2) a creep potential function which describes the dependence of creep on stresses: (3) a definition of equivalent creep strain rate for a multidimensional system; and (4) a constitutive relation of creep strain rate and stresses. The in-reactor creep behaviour of a material depends on stresses o, temperature T, time t, and neutron flux ~b, as discussed by Gittus [25] and Robertson [26]. Thus, a general expression can be written as ~- = if(o, T, t, ~b).

(45)

However, laws for specific materials must come from correlation of creep test data. In this paper a modified Drucker-Prager yield function is adopted for the creep potential function,

G = AI~ + J2.

(46)

Here A is referred to as the compressibility constant, 11 and Jz are the first stress invariant and the second deviatoric stress invariant, respectively. By analogy with the conventional theory of plasticity, the constitutive relation for creep strain rate ~. is given

~j =/3(~G/~oij), where/3 is a scalar quantity which may vary during loading.

(47)

333

S. Baner/ee et al. /Analysis of fuel behaviour in transients

The equivalent creep strain rate is defined as ,2 .c • c',1/2

= [~6ijeij

)

(48)

•

In eqs. (47) and (48) the dot ( ' ) indicates the derivative with respect to time. The total strain rate can now be expressed as the sum of the total strain rate {~'} obtained from eq. (21) and the creep strain rate in eq. (47). Adopting matrix notation, we have {~ } = {~,} + {~ c}.

(49)

A complete constitutive relation for the combined thermal-elasto-plastic and creep analysis can thus be obtained from eqs. (47), (49) and {4}=[Ce]-l{a}+

(

(~}+ a [ C e ] - l ( ° }

aT

)

J'+

a [Ce]-' {o}

a(~}

d{4}+X

aF

a-To}

+3

aG

a(o)

•

(50)

With all the previously mentioned requirements established, a thermal-elasto-plastic and creep problem is now reduced to the solution of eq. (50). The two constants A and 3 will now be discussed. The fuel cladding and UO2 pellet will be treated separately. 2. 3.1. Creep model for the fuel cladding In common with the usual assumption in creep analysis, we assume that creep deformation does not involve volume change and the material is incompressible. The condition of incompressibility can be stated mathematically as

~tt = 0.

(51)

To satisfy the above condition, the compressibility constant A in eq. (46) is set to zero. Thus, the creep potential becomes G =J2 •

(52)

Eq. (47) may be expressed as '¢=

co

(53)

3(aJ2/aoo).

Combining of eqs. (48) and (53) results in (54)

= [3(~ SijSij) 1/2 ,

where Stj is the deviatoric stress tensor as previously defined. From eq. (54) the following expression can be obtained for 3:

~3 ~ ( ~ / ¢ ) ,

(55)

=

where F is the equivalent (or effective) stress as defined in section 2.2. Introducing eq. (55) into eq. (53), the constitutive relation for the creep of the cladding material can be written in matrix form as

~o/

[sym.

b

Ozz

a

oo~

3re

(56)

S. Baner/ee et al. / Analysis or fuel behaviour in transients

334

Where a,b and c have the value of 1, -½ and 3 respectively. A suitable creep law for the cladding material can be introduced for the equivalent creep strain rate e c in eq. (56). Eq. (50) is therefore fully determined. Z3.2. Creep model for the U02 pellet material Because of the porosity of UO 2 pellets, creep deformation is likely to involve a volume change. The assumption of incompressibility immediately becomes questionable and the equation used for fuel cladding is no longer valid. The compressibility constant A in eq. (46) has to be redefined. To relate the compressibility constant to the voids in the pellets, the material is modelled as a composite of voids and solids resembling closely packed hollow spheres. The structural material of the hollow sphere is assumed to be incompressible. The spheres are subjected to an interior pressure Po. From the above assumptions one can derive an expression for A in terms of void ratio 8. The expression is presented here without proof,

A = ~ [3/2m]m8 ,

(57)

where m is the exponent of a power law stress dependence. It should be mentioned here that any relevant model of a composite material can be used to arrive at an expression relating A with 8 provided a solution can be obtained analytically or numerically. Here closely packed hollow spheres are used merely for the mathematical convenience that a solution is readily obtained. Eq. (47) can be rewritten

(58)

4c]= j3(2AllSil + aS2 ~ \ aoq]" The equivalent creep strain rate is now defined as 2 C e*-2-- (3dq ~)l/z

(59)

where ~ci = [J(aJ=/aoq) = deviatoric creep strain rate. Thus, ~ can be determined from eq. (59) = g(~-e/~-).

(60)

Substitution into eq. (58) of A and ~ from eqs. (57) and (60), respectively, results in the constitutive relation for the creep of UOa pellets in matrix form "c

a

egr .c egz

_i ~c

.c G06

0

b

b

0

OfF

a

b

0

Ozz

(61) sym.

a

0

G06

C

orz

where a=

b=~ C=~

v+l

2 \2m/

3

m

~-o.5,

.

A suitable creep law for UO 2 pellets can now be introduced in place of the equivalent creep strain rate in eq. (61). Thus eq. (50) is fully determined and a solution can now be sought specifically for UO2 pellets.

S. Baner]ee et al. /Analysis o f fuel behaviour in transients

335

2.3.3. Solution procedure The solution procedure follows exactly the samle approach as that given in the elasto-plastic analysis. However, the integration must now be carried out with respect to time as well as space. The creep strain rate enters the finite element formulation as a pseudo load vector,

{ARc} = _ f f [~]r [D] {~c} do d t ,

(62)

where {AR c} = the increment pseudo load vector. The rate of stress change can be obtained from eq. (50)

{~'}=([ce]-[ce] {af/a°){af/a°}rs[Ce])({~}_ [{,~]. a[Ce]-' (o}]T)ar [Ce] {c3//aO}S{af/aT} ,j,_ ([i] _ [Ce] {af/ao}~{af/a°}T )LCeJ"" {~e}.

(63)

2.4. The gap element The presence of a gas mixture between the fuel and the sheath introduces difficulties in both the thermal and stress analysis described in the previous sections. A standard quadrilateral element is used in these areas, with a special development of the conductivity and stiffness matrices. The gap is considered to close if the average thickness becomes less than an arbitrary amount and to open if the stress normal to the principal axis of the gap (see fig. 8) changes from compression to tension.

2.4.1. Conductivity The "equivalent conductivity" of the element is calculated as (64)

Ceff = heffAZgap ,

where C e f f = effective conductivity, heft = effective gap heat transfer coefficient, AZgap= mean gap thickness. The effective gap heat transfer coefficient is treated as the sum of solid-solid conduction, gas conduction and

gap

Principal axis

Fig. 8. Gap element.

S. Banerjee et al. / Analysis of fuel behaviour in transients

336 radiation terms,

heft = hr + hg + h s ,

(65)

where h r = radiation term, hg = gas conduction term, h s = solid-solid conduction term. The radiation term, is calculated as hf =

asBe(T ~ + r~) (TF + T s ) ,

(66)

where OsB = Stefan-Boltzmann constant, e = emissivity, TF = absolute fuel temperature, Ts = absolute sheath temperature. The gas conduction term is modified by the temperature drop at each surface due to gas molecule/surface collisions by increasing the gap by a "temperature jump distance" as suggested by Dushman [27]:

1

- Ps ~i

gl + g2

boTg

Yi

(67)

" (gl + g2)i '

where

(gl + g2) P~ rg

= temperature jump distance, = gas pressure, = gas temperature, Y; = mole fraction of gas i, fgl + g2)t = temperature jump distance of gas i at standard pressure and temperature = constant. bo Then the gas conduction contribution is

c~ hg - c(RF + RS ) + (gx + g2) + AZ~ap '

(68)

where C~ = gas mixture conductivity, RF = fuel surface roughness, R s = sheath surface roughness, c = constant. Cg is a function of both pressure and temperature and can be expressed as

Cg(P, T) = X ( T ) Y(P, r ) .

(69)

For a single gas X ( T ) Ubish [28] has suggested the form

X ( T ) = Ko Ts

(70)

where Ko, S are constants for the gas. For a mixture

X(T) =

zt YiXi( T)MI /3 riy~Ml/s

(71)

S. Baner]eeet al. / Analysis or fuel behaviour in transients

337

According to Lenoir et al. [29], Y(P,T) is given by

where -Pc = zlYiPct, Tc = XiY~Tci, Pci = critical pressure of gas i, Tc~ = critical temperature of gas i. The solid-solid conduction term is given by Ross and Stoute's expression [30]

Cm

P

(73)

hs - ao~1/20"y' where Cm

= 2 C F C s [ ( C F + CS)

= harmonic mean fuel/sheath conductivity, = interfacial contact pressure, = root mean squared surface roughness of fuel and sheath, Oy = yield strength of the sheath material, a0 = constant. The error introduced by using Ceff in an isotropic manner when a portion of the term h7 is derived from considering radiation (with is anisotropic), introduces acceptable errors because the aspect ratio of the gap element is large.

P

2.4.2. Pressure loading in contained volumes The pressure of the gases in the fuel sheath depends on the temperature of the gas, the quantity of the gas and the free volume of the container. At each computation step a new value of gas pressure is calculated. The gas in an element is assumed to be at the mean temperature of the element nodal temperatures. The ideal gas law leads to

P =n R / ~ !

Vt/Tt,

(74)

where nR = quantity of gas, Vt = free volume of element l, Tt = free volume of element l, The change in pressure is applied in the next load step as (APe).

2.4.3. Gap snffness The gap is treated as an elastic anisotropic element in the stress analysis. The anisotropy is designed to produce non-zero stresses normal to the principal axis of the element and zero stresses elsewhere. In the present analysis no sliding friction is considered. The general isotropic elastic relationship dot

b'

P

1 -v

1 -v

p

doz

1 g ( 1 - v)

1-v

0

der

0

dez

0

deo

(1 -- vX1 - 2v) do o tim

1

-

2v

A~

338

S. Banerjee et aL / Analysis o f fuel behaviour in transients

is modified by the assumption v = 0 and to have only normal stresses, to give sin 4 ~b

do r

doz doo I

=E'

sin2ff COS2~ 0

sin ~ cosaff

dEr

Cos4 ~

0

sin 3 ff cos

dez

0

0

de0

sin 2 ~ cos 2

derz

sym

dorz

(75)

where ff is the angle between the R axis and the principal axis c~fthe gap element and E ' is the effective modulus of elasticity assigned to the element. When the gap is closed E ' is assigned an arbitrarily large value that will result in little further deformation of the element under the interfacial pressures expected. This value of E' is temperature independent. When the gap is open an estimate of the "stiffness" of the gas in the connected volume is used. This is calculated by reducing the thickness of each gap element by 10% of the gap closing tolerance. Then a new pressure may be calculated from P' = nR

v;/T

,

(76)

l

where V; is the modified volume of element l. The first order estimate of E' is E ' = ( P - P') AZgap/0.1 rig,

(77)

where dg = gap closing tolerance and AZgap = mean gap thickness. The 'fictitious stresses' produced by the use o f E ' ¢ 0 are considered to be first order estimates of the pressure load change. The solution procedure reduces these stresses to zero by evaluating ARgap = - f [B] T (Ogap) dV o

(78)

at each step when the gap is open. Since the load ARgap is applied at the same time as the pressure load A/l, the result is a globally correct pressure load and a first order estimate of the incremental pressure load. The rationale for this approximate treatment is as follows. A 'correct' procedure would require that the displacement of every node on a connected gas volume be related to the displacement of every other node on that volume. The bandwidth of the resulting set of simultaneous equations is large enough to increase the computation time unacceptably. I r E ' is set to zero, an oscillation can arise because the load increment lags the deformation by one solution step. The compromise adopted has produced acceptable results in problems that otherwise displayed tlnacceptable oscillatory behaviour.

PART II. VALIDATION OF MODEL

3. Comparison with analytical solutions The finite element thermal and stress analysis described in section 2 has been developed into subroutines in a general computer program for predicting fuel behaviour. Each subroutine has been tested separately by comparison with several analytical solutions. The whole program has also been compared with experiments as described later.

S~ Baner]ee et al. / Analysis or fuel behaviour in transients

3" R

0,3"

339

IOolement$

23 , , e~ 8 , , o , , , ~ / Pi ~ ' 1 2 [ 3 1 " l s l e l : ' l e l g l '° t " , ,3 ,,;,5,e , r , 8 , 9 ~o~,2~

El I

PI :7650 psi (52.7MPo) E = I0 x 106psi ( 6 8 . 9 M P a }

b' : 0.33 G : Ct:

t/1/////#///////

~.76 x l O 6 p s i

(25.9MP(I)

O/°F

O"

I0 k$i I I

"1':':

IZ" 0.0

Fig. 9. Geometry, finite element idealization and material properties o f a thick-walled cylinder.

3.1. The elastic-plastic behaviour o f a thick-walled cylinder The stress analysis subroutine was applied to a structure subjected to mechanical loads only. The analytical solutions by Hodge and White [31 ] for a thick-walled cylinder loaded by uniform pressure applied at its inner surface have been selected for the comparisons presented here. The dimensions, material properties and finite element idealization of the structure were as shown in fig. 9. The analytical solutions for the tangential stress and radial displacement of the outer surface of the cylinder are compared with the finite element calculations in fig. 10. The two solutions agree well.

I.O

I.oAnalytical

o 0.8

solutions

Finite e l e m e n t

solutions

(TEPSA program)

O.S -

°S 0,6

O

'

0.6

b"

%, Q. 0.4

0.4

PLASTIC

~ ~ ELASTIC --

0.2

0.0 1.0

0.2

1 1.2

I 1.4

I 1.6 r/o

I 1.8

I 2.0

0.0 0.0

/

o/

b=2a

I 0.8

I 1.6

2.4 4Gub

I 3.2

4.0

%0 Fig. 10. Comparison of finite element results with analytical solution of a pressurized thick-walled cylinder.

S. Baner/ee et al. /Analysis or fuel behaviour in transients

340

l

Z

0o

0o

°F

Oo

®

@

@

720(75 ° )

680(69 ° )

5 6 ° ( 5 6 °)

590(58 ° )

8

®

@

,Z

Temperature

'®

'@

162°(170 °) 153°(154 °)

125°(125 ° ) 86°{86 °)

®

r@

@

@

2950(:51301 276°{279 °) i223°(224°) 151°(152 °)

®

r

®

0o

0°

19°(17 ° )

Oo(O 0)

'®

~®420142 ° )

'@ 750(7401

@

0°(0 °)

~ 0 o ( 0 °)

®

o 509 155301 4610148001 380°(371 °) 246°(251°J 125°(123° )

©

®

'®

®

'@

.60Q °

400 °

.200 °

O°(O °)

®

Fig. 11. Temperature distribution in a solid cylinder. (Note: numbers in parentheses are the temperature calculated by the TEPSA program; numbers in closed circles and the nodal numbers in F.E. calculation.)

3. 2. The temperature distribution in a cylinder The thermal analysis subroutine was applied to predict temperature distribution in a cylinder under non-uniform heating conditions. The analytical solution is available in Carslaw and Jaeger [32]. Consider a solid cylinder of homogeneous material with radius a and length 2L. The surface temperature of the cylinder is kept at zero except at its mid-plane where the temperature varies according to f(r). By using the coordinate system shown in fig. 11, the following differential equation and boundary conditions apply: ~2T

I ~ T a2T

8r 2 + r~r + az 2 = 0 T(r,z)lz=L = T(r,z)lr--a = O,

r(r,Z)lz=o =/'(r). The solution of the above equation satisfying the specified boundary conditions may be found in Carslaw and Jaeger [32] and is o~

T(r,z)= 2 P1- J°(c~nr) sinh(L-z)C~n~s~nh(La~ /o rf(r)J°(~nr)dr'

S. Baner/ee et al. / Analysis or fuel behaviour in transients

341

where a n are the positive roots of the transcendental equation Jo(ana) = O .

Results from this analytical solution are compared in fig. 11 with computed results for the following input data L =a = 25 mm;

f ( r ) = 811 (1 - r ) K .

The highest deviation of the results occurs at the central point in the cylinder (i.e. r = z = 0). This maximum deviation is approximately 10 K. Good agreement is obtained for other parts of the cylinder.

4. C o m p a r i s o n w i t h a l u m i n u m t u b e d e f o r m a t i o n tests *

Although the present finite element analysis has been validated to some extent against analytical solutions, experimental verification has been complicated by the following factors: (a) difficulty in controlling and measuring the combined thermal and mechanical loadings; (b) difficulty in measuring large plastic deformation of solids at elevated temperature; and (c) scarcity of temperature-dependent material properties. These difficulties have led in general to progress in thermo-plastic analysis being more rapid than its experimental validation. The most relevant reference in this area is due to Phillips and Kasper [33]. However, their emphasis was on the observation of the shii't of yield surface of a tube subject to various combinations of thermal and mechanical loads. The specimen temperature was limited to 505.37 K, the limiting characteristic of the special strain gages used in the investigation. No attempt was made in their work to correlate experimental results with theoretical analyses. 4.1. E x p e r i m e n t a l setup

The present experimental investigation was designed to verify various aspects of the finite element analysis described previously. The test rig was designed in such a way that it had many characteristics similar to the fuel tested, e.g. internally pressurized thin-walled tubes heated by a solid internal heater. The tubes were subjected to combined thermal and mechanical loads for either steady-state or transient conditions. Repeated (or cyclic) loads could also be programmed to assess the residual stress-strain induced. The test rig used nitrogen as the pressurizing medium and a tubular heater was inserted inside the tubular specimen to provide heat to the specimen by radiation. Resistance heating of the specimen was not suitable because the tubular specimens undergo large plastic deformation resulting in local thinning of tube walls which in turn causes uneven heating of the specimen. 4.1.1. Specimens

Two types of tubular specimens were used in the present experiments; plain tubes and tubes with a circumferential collar. The geometries and dimensions of these tubes are shown in fig. 12. The tubes were made of 6061 T6 aluminum alloy. This material was chosen because its mechanical properties are readily available. Typical stressstrain curves for this material at various temperatures are shown in fig. 13.

* SI units are used throughout the text. However, English units are used to avoid possible inaccuracy due to conversion in cases where the test specimen is specified or the equipment manufactured in these units. For these cases, conversion factors are shown on the corresponding figures.

342

S, Baner/ee et al. / Analysis or fuel behaviour in transients INTEGRAL

x

. . . . .

I / 8 " ROUNDS

z o

~ I

.,,

1

~-~J

.-.I,.,.

!

1

MATERIAL : A L U M I N U M 6 0 6 1 - T 6

Fig. 12. Geometries and dimensions of tubular specimens.

4.1.2. Heater

A special internal heater was developed for the present experiment. Heat is generated by passing electric current through a thin-walled tube made of 0.20 mm stainless steel shim stock. To minimize differential thermal expansion between the stainless heater and the concentric aluminum specimen surrounding it, an expansion joint made from a series of bent steel strips was included in one and as shown in fig. 14. The o.d. of the heater was 38.1 mm.

4.1.3. Specimen grips As the tests were done at high temperature and pressure, specimen grips capable of sealing hot high pressure gas in the specimen and also providing sufficient electric insulation between the specimen and the heater unit had to be designed. The system shown in figs. 15 and 16 proved satisfactory. As seen from fig. 15, the tubular specimen is screwed into the grip and the O-ring situated between the grip and the stepped brass terminals provide adequate sealing for the high pressure gas. The Teflon electric tape wrapped around the brass terminals serves as electric insulator.

70

STRESSPOWER -

4o

~

3o

I

~~?~6 ~

,\0 .b'2~

g

/

Tt. =LO

~/

~9

N

"

0

4

8

12

16

20

24

28

TRUE STRAIN , %

Fig. 13. Temperature-dependent mechanical properties of tube material. NB. lpsi = 6.59 kla; 1 ksi = 6.59 MPa.

S. Banerjee et aL / Analysis of fuel behaviour in transients N.B,

I"=

343

25,4mm

= F L E X BLE CONNECTOR I S.STEEL SHAFT~

WELD I0 TO 3 i ~SPOT ,

//

SPOT WELD 3 TO I--

~,,

~8 ~-9 ~--

--I0

12 MARINITE INSULATING SPACER~

MARINITEINSULATING SPACER

WELD II TO 7-~

'

i I i

S.STEEL CAP

THREADED CONNECTOR

NOTE: ITEM 7 MATERIAL TO BE ROLLED TO I - I / 2 " I . D . HEAT TREATED AND SPOT WELDED ALONG SEAM I=- ~ : = : = ~ ~1

;

± 7

'

1.5" [.D,xO.OOS"WALLx9-3/B"LONG S.STEEL TUBE HEATING ELEMENT (PRIMARY)

t L-6 O.060"THICK FIREBRICK H.T.CEMENT LAYER

L 5 1.375"O,D.xO.O35"WALLxS"LONG S. STEEL SLEEVE - 4 MINERAL WOOL FILTER 3 3/4"O.D.xO.O35"WALL x 8 - i 5 / 1 6 " LONG S.STEEL TUBE HEATING ELEMENT (SECONDARY)

2 MINERAL WOOL FILTER

Fig. 14. Assembly drawing of heater design.

4.1.4. Arrangement o f equipment As shown in figure 17, the assembled specimen grip is fastened to and contained in a test chamber which is evacuated for high temperature work and protects personnel in case of tube rupture. Chrome-alumel thermocouples (0.25 mm) were spot-welded onto the outside surface of the tube. The thermocouple located at mid-level was used for the specimen temperature control. In most cases, two strain gages capable of measuring strain up to 2% at 588.70 K were cemented on selected points at the tube surface. The radial displacement of the tube was measured by a linear variable differential transformer (LVDT) which was mounted on a platform capable of moving freely in the horizontal plane as shown in fig. 18. This was necessary to allow for rigid body motion of the specimen during deformation. Thermal and mechanical loads were controlled by Paramac controllers driven by load curves scratched on the drums of a dual-channel 'Data-trak'. These Paramacs controlled the opening of the high pressure gas valve for specimen pressurization and also the electric current output of the d.c. supply for heating the specimen. The measured data including thermocouples, strain gages, pressure transducer and LVDT, etc. could either be recorded on multi-channel recorders after signal conditioning, or be fed into a data acquisition system and processed in a large computer as illustrated in fig. 17. A video-camera/tape recorder/CRT monitor unit was used to monitor and record variation of tube profiles continuously during the tests. 4. 2. Experimental methods~results Three sets of experimental results will be presented as follows: (1) A 1.19 mm tube was first pressurized by nitrogen gas to 1.95 MPa. It was then heated by controlling the temperature increase at the half-length point at a rate of 0.292 K/sec as shown in fig. 19. Also shown in this figure is the measured radial displacement by a LVDT placed at a mid-level of the tube. Measurements were done till the tube ruptured, usually about 20 min after the heating started. The gas pressure inside the tube was maintained at a constant level of 1.95 MPa by the Paramac controller. As a result, the system became an 'open' system with continuous bleeding of nitrogen gas into the system in order to maintain a constant pressure as the tube heated up.

344

~ Baner/ee et al. / Analysis of fuel behaviour in transients

rced Teflon

Ring

nless Steel - Energizing Inlet

3n Sealant

Fig. 15. Tube grip assembly. Consequently the upper end of the tube; where cold gas entered, became cooler and the measured temperature profiles along the tube length were as shown in fig. 20. The maximum temperature in the tube occurred at 1 in. be low half-length. This was where the tube burst (fig. 21). A 152.4 m m long section of the tube, 76.2 mm from the top grip was used for the finite element model verification. A pseudo-cap was added to the model as shown in fig.

Fig. 16. Photograph of tube grip assembly.

D.C. POWERIN

GABPRESS.FEED-BACK

I

~

II

[~

VIDEO

,I

LVDT ~

FT~7~-7

~

i E.E. CODE~

DE~OW~ROOT~

COMPUTER I

CRTMONITOR

CAMERA

1I

GASOUTLET r"""~~~~L~ili~'P L:t c'- ~ ~R~E

_I U.M. ZDMBTO/IS81

9

B~E~M~N

----IZ .........' ~ - - ~ - - ~ ~

/I

PARAMAC CONTROLLER FOR TEMPERATURE i~ DATA-TRAK

/PREB~ORE

I

EH. 1-TEMP.

CONTROL

MOTOR ; \

CH. 2-PRESS.

CONTROL

~

'/;~"

'x h

~

~

,) . . . . . . . . . . . . ),,

"

~,~

1 THERMOCOUPLEFEED-BACK

Fig. 17. G e n e r a l a r r a n g e m e n t

of apparatus.

_-7."

346

S. Banerjee et aL

/Analysis of fuel behaviour in transients

Fig. 18. Perspective view of mofified LVDT assembly. Key: (1) LVDT; (2) tubular specimen; (3) micrometer head; (4) guide plates; (5) floating platform; (6) damp; (7) support; (9) knife edges.

22. A total of 183 nodes and 227 elements were used in the calculation. The internal pressure to the tube was maintained at 1.95 MPa and the thermal load was applied according to the temperature variation and distribution shown in fig. 20. The radial displacement calculated by the F.E. code is compared to that measured by LVDT in fig. 19. It is seen that agreement is good except very near rupture. (2) A tubular specimen with a distinct step (collar) at the half-length was used to study the effect of geometric discontinuities. A collar introduces to some extent effects similar to that of a bearing pad on a fuel sheath. The LVDT was positioned 25.4 mm below the collar. The tube was first pressurized to 1.38 MPa and then heated up

BOTTOM GRIP

wt

--IOI ......

~400F

~

-~

I.C

1.85"0.D x 0.047 TK

~2oo

i 3ooi--

c O

I

/

I

TOP GRIP

I'~LVDT5" 'FR~OMTOPJ ' ~ l

,~

/ p

,./

i ~

~

I-" z

~ 1 0 I hi laJ

to

(3 a
0[=-" -5

~;~ 0

K/ 5

I0

15

TIME (Minutes)

Fig. 19. Applied loads and correlation of radial deformation of tube under constant pressure.

20

S. Banerjee et al. / Analysis or fuel behaviour in transients

~d

w

.d

.-..I.-

400

Fig. 20. Temperature variation along tube length.

Fig. 21. Photograph of ruptured tubes.

347

348

S. Banerjee et al. / Analysis or fuel behaviour in transients 6.9384"

84

No. of nodes in Region I -

9g

No. o f nodes in Region II =

84

Total No. o f nodes = 183 5.9 " No, o f e l e m e n t s in Region I ~ 132

No. of elements in Region I! ~ 95 Total No. o f elements - Z27

5.894"

O.II8" On

Fig. 22. F i n i t e e l e m e n t m o d e l f o r t h e c o r r e l a t i o n o f e x p e r i m e n t a l m e a s u r e m e n t s . N o . o f n o d e s in r e g i o n I = 9 9 ; n o . o f n o d e s in reg i o n II = 8 4 ; t o t a l n o . o f n o d e s = 1 8 3 . N o . o f e l e m e n t in r e g i o n I = 1 3 2 ; n o o f e l e m e n t in r e g i o n 1I = 9 5 ; t o t a l n o . o f n o d e s = 2 2 7 .

by programmed electric current supplied to the internal heater. One noticeable difference in results is the measured temperature profile along the length of the tube as shown in fig. 23. It is interesting to note that temperature dropped in the collar region. Again, to maintain constant gas pressure in the tube, cold nitrogen was purged into the tube during the test, which made the upper part of the tube colder. The temperature profile shown in fig. 23 was used as the input loading to the finite element model similar to fig. 22 except with the provision for a collar. Fig. 24 compares the LVDT measurements to the radial displacement calculated by the finite element program. Again, good agreement was obtained except very near rupture. The final shape of the ruptured tube is shown in the top of fig. 21. (3) A combination of thermal and mechanical load was applied to an aluminum tube of uniform thickness as described under (1). The variation of the loading conditions is illustrated in the upper part of fig. 25. Two temperature ramps were programmed. They were 0.89 K/sec for the earlier portion of the test followed by a lower rate of 0.46 K/sec Gas pressure was maintained at constant level of 0.31 MPa and started to increase after 5 min at a rate of 0.0043 MPa/sec. The measured temperature profdes along the tube length at various instants are shown in fig. 26. These profiles, along with the pressure loadings, were used as input to the finite element stress analysis com-

S, Banerjee et al.

/Analysis o / f u e l

behaviour in transients

349

\

\

T i t Cpt. No.

X~ I

t

I 6 MI

E '5

7

L,J

I

/

/

I /

200

IO0

/ I

300

400

TEMPERAI~URE,*F Fig. 2~. Temperature profiles along a tube with a collar.

400

" -

Gas Pressure , kPa ( g a u g e )

-

I1 -----

f

Programmed Temperature at the Collar, *C

tO 300

-

x

x

Measured Temperature,*C

9

o~

....

t. _o

o ®

(,~

"6

Measured Radial Displ.,inch

8

Radial Displ Cal'd by Tepsa Code

-

200

\

oJ

\ \

E oJ

\

I.I00

-

\

x~ \ \y..

OI

0

I

2

3

4

5

I0

15

T i m e , min

Fig. 24. Applied loads and correlation o f result o f a tube with a collar.

I

20

S. Baner]ee et al. /Analysis o f fuel behaviour in transients

350

4ool.~

o

°-300

/

'5°I-~

/

I~ 'ooI- /

°°l-x

Th.CpI.Nos

/ ~"

l~2OOF

IZ

/

../

p.,/

.....

oL o , F ; i i i i , , , 0

I

2

~

4

"" 5

6

7

8

I I I I 9 1 0 1 1 1 2

TIME, Min.

Measured by LVDT

4,5,6 e

{

lO

LVDT 'o

"" 7

6

7, 4 i

I

I

/

/

/

i ./

/

/ /./ / /

i Ii I I/ ii

/

z/

//

117

I II 11 /~If/ 11 / ,///

2 I I

0 0

t

I

I

I

I

I

I

I

I

I

I

I

I

2

3

4

6

6

"r

8

9

I I0

I II

Fig. 25. Correlation o f result o f a tube subject to c o m b i n e d variation o f thermal and m e c h a n i c a l loads.

I 12

i fit" I00

,'

~"," I/ 200

/

/"

TEMPERATURE

I 300

I 400

, °C

Fig. 26. Temperature profiles along a tube.

puter code and the correlation of the measured and calculated radial displacement at 19.05 mm below the halftube length is shown in the lower part of fig. 25. The shape of the ruptured tube is shown in fig. 21.

4. 2.1. Discussion of results Although only a limited number of experimental results have been reported here on structures of relatively simple geometries under moderate loads, the agreement between experimental data and predictions has been encouraging. This has been summarized in figs. 19, 24 and 25. The present experimental setup for generating and controlling combined thermal and mechanical loads is also applicable to more complicated experiments which will be done in future.

5. Comparison with an in-reactor experiment To illustrate the applicability and validity of the thermal-elasto-plastic analysis model, a fuel pin is analysed for an in-reactor experiment by Notley et al. [34]. The experiment was run for a 5-hr period between reactor stari. up and shutdown, and is the shortest in-reactor test available in the literature. It is believed that with such a short power duration the creep effect on the thermal-mechanical behaviour of a fuel pin would be negligible. It is there-

U2

I II

i

~-

¢e..

U2=O

! --B

~f

Fig. 27(a) F u e l pin c o n f i g u r a t i o n a n d finite e l e m e n t i d e a l i z a t i o n . A = p e l l e t m i d p o i n t ; B = p e l l e t i n t e r f a c e ; h = p e l l e t l e n g t h = 2 4 . 0 0 5 n m ; g = gap t h i c k n e s s = 0 . 0 7 6 m m ; r = p e l l e t radius -' 9 . 2 4 5 m m ; t = s h e a t h t h i c k n e s s = 0 . 5 5 8 m m . I

I

I

I

i

I

I

i

I

240

2OO

1.0

E

"v'• e'~,,,e0

0.8

160~ 0

Z

)

~- 0 . 6 ¢/)

d20 w

.J Iz

K 0 a.

0

•

0.4

80

e

¢.D Z

•---.-0- TEST •

0.2

STRAIN

COMPUTED STRAIN FUEL

_1 "' n,

40

POWER

0

~ O J o / ° -O.2

I

I

2

4

I

6

I

8

I

I•

I

12

I

14

I

16

I

18

I

t x IO 3 (sec) Fig. 2 7 ( b ) V a r i a t i o n o f t a n g e n t i a l s t r a i n o n the o u t s i d e surface o f t h e s h e a t h at a ridge. F u e l p o w e r h i s t o r y is s h o w n for reference.

352

S. Banerjee et aL /Analysis o f fuel behaviour in transients 2.0

I

I

I

18

I

I

I

I

I

I

I

I

i

I

I

I

i

I

I

54.6 MPa _

450°C 54.6

1.0

78 MPa

,< < 06 -

I

///78 /)

t4 I 2

I

EXPERIMENTAL PREDICTED

1.6

o~o

i

AN~AL ~ T ~

~

/,J

~ A

ff - 450Oc

"

/ / / - ' ~ "

MPa

/

~

I0

/

i

-08

~

04

-

~oo~

MPa

sTRAI~

-o6

~

-o4

m~

02

o 0

I

I

I

I

I

I

I

I

I

I

I

I

I

I

i

I

I

I

I

I

1

50

I00

15@

200

250

300

550

400

450

500

550

600

650

700

750

800

850

900

950

rO00

1050

I100

TIME (HRS.) Fig. 28. Measured and predicted creep strains of a fuel cladding subjected to various temperature and loading conditions.

fore used for the validation o f the thermal-elasto-plastic model involving no creep effects. A segment of the fuel pin including 5 pellets is analysed with constraining condition as shown in fig. 27(a). Because of symmetry, only one-quarter of the segment needs to be considered. The finite element idealization is also shown in fig. 27(a). The fuel pin is subjected to a temperature of 490 K and a constant coolant pressure of 7.76 MPa. The power history in terms of heat generation is shown in fig. 27(b). Both computed and experimental results of hoop strain vs. time, at the location of a ridge * are presented in fig. 27(b). These results demonstrate that the predictions are in agreement with the experimental results.

6. Validation of creep model For the purpose of validation of the creep model, two sample problems, for which creep test results are available, have been analysed. The first problem was taken from an out-reactor uniaxial creep test by Kohn [35] on a Zr-2.5% Nb fuel cladding with outside radius 0.674 cm and thickness 0.043 cm. The test was carried out under various temperature and loading conditions. At the start of the test the specimen was subjected to an axial tensile loading of 78 MPa with temperature at 673 K. After 700 hr the loading was reduced to 54.6 MPa, and the temperature was kept at 673 K. At 870 hr the temperature was increased to 723 K, while the load was kept at 54.6 MPa. The load was then increased back to 78 MPa while the temperature remained at 723 K. The total test time was approximately 1060 hr. The predicted axial creep strain is shown together with the experimental observation in fig. 28. Clearly, the resuits are in close agreement. In the same figure the predicted hoop strain is also given without comparison because no recorded test data are available. However, the expected contraction of the cladding is clearly demonstrated by the negative hoop strains. In order that the creep model may be used with confidence a more critical validation, namely comparison with in-reactor test, is desirable. Thus the second problem was taken from an in-reactor biaxial creep test by Kohn [35]. * A ridge is a type of fuel cladding deformation peculiar to pelletized fuel pins. The maximum deflections coincide with pellet interfaces and extend around the circumference of the fuel pin giving a characteristic bamboo-rod appearance.

35 3

S. Baner]ee et al. / Analysis or fuel behaviour in transients C9

.7620 cm,

[

i

I~"

"F '///~////////////////

~l

'

~

P = .537 M P a

T,

11.1125 cm

I"l

I

1,5

I

I

I

I//////A

~1

I

I t

_

~ I

I~I

I

t

+ Predicted • E~;periment

'~

=-I

/ ///////////////////////////J~77"/77~

4

~

I.O

0.5

I 200

I 400

I 600

I 800

I I000

i 1200

I 1400

1600

Time (hrs.)

Fig. 29(a). Dimensions of a specimen used in the {n-reactor biaxial test; 27(b) correlation of creep strain for an in-reactor test. The test was done using eight identical Zr-2.5% Nb tubes which were subjected to the same internal pressure, 0.537 MPa. The structural dimensions of a typical cladding are shown in fig. 29(a). The total test time was approximately 1060 hr. To obtain the strains at different time intervals the tubes were withdrawn from the reactor test loop one by one, each at a chosen time interval. The hoop strain was measured immediately after withdrawal of each tube. The predicted hoop strains of the cladding are compared with the test results in fig. 29(b). The maximum discrepancy between the predicted and measured strains was found to be approximately 0.35%. This discrepancy may be partly attributed to the possible difficulty of maintaining a uniform temperature and pressure over all the test specimens and throughout the testing process, and partly attributed to elastic recovery when the tubes were withdrawn. In ,dew of these intrinsic difficulties, the comparison may be considered satisfactory.

7. Conclusions

and

recommendations

Though the complete fuel model described in this paper has not been compared with extensive in-reactor tests, each separate part appears to be accurate and satisfactory comparison was obtained with our in-reactor experiment where creep was unimportant. More in-reactor tests to check the model are needed and are being planned. The process of validation by comparison with out-reactor experiments with continuous thermal and mechanical loads is also continuing. Even in its partially validated form the fuel model is of use to designers as three-dimensional axisymmetric and two-dimensional plane strain problems can be calculated. The model is particularly useful in predicting fuel beha,dour in fast power and thermohydraulic transients, and during load-following.

354

S. Banerjee et aL /Analysis of fuel behaviour in transients

The immediate requirement is to increase computational efficiency and decrease the ratio of computational time to real time. The present stress analysis is based on small deformation theory. A stress analysis based on large deformation theory is being developed and will be incorporated in the code. It promises to reduce computer time requirements substantially without loss in accuracy. The fuel model is also being improved to include (1) fuel failure criteria; (2) frictional forces at fuel pellet-cladding and fuel pellet-fuel pellet interfaces; (3) cladding oxidation during temperature transients; and (4) better modules for fission gas release, fuel swelling, and pellet cracking. For certain applications a fully three-dimensional fuel model would be desirable. To be of practical use, major developments to increase computational efficiency are required.

Acknowledgement The authors are indebted to E. Kohn and M.G. Wright for supplying the creep data and material properties, to P.A. O'Connor for his assistance in computation, to W.C. Harrison for help in defining some experiments, and to W.C. Harrison and A.D. Lane for their encouragement and support in the development of this work.

References [1] [2] [3] [4] [5] [6]

B.M. Ma, Nucl. Eng. Des. 7 (1968) 181. J.Ho Gittus, Nucl. Eng. Des. 28 (1974) 252. S. Levy and J.P.D. Wilkinson, 2nd SMiRt Conf., Berlin, Germany, Sept. 1973, Paper D1/5. Y.R. Rashid, Nucl. Eng. Des. 29 (1974) 22. P. Royl, A. Watanabe and A. Agrawal, 2nd SMiRT Conf., Berlin, Germany, Sept. 1973, Paper C3/1". T.R. Hsu, A.W.M. Bertels, B. Arya and S. Banerjee, in: Computational Methods in Non-linear Mechanics, ed. J.T. Oden et al. (TICOM Press, 1974) pp. 521-540. [7] T.R. Hsu, A.W.M. Bertels, S. Banerjee and W.C. Harrison, AECL-5233 (1976). [8] J.J.M. Too, T.R. Hsu and A.W.M. Bertels, in: Fuel Element Analysis (ASME special publication), ed. Y.R. Rashid and F.C. Weiler (1975) pp. 23-36. [9] E.L. Wilson and R.E. Nickell, Nucl. Eng. Des. 4 (1966) 276. [10] M.E. Gurtin, Quart. Appl. Math. 22 (1964) 252. [11] Y. Yamada, N. Yoshimura and T. Sakurai, Int. J. Mech. Sci. 10 (1968) 343. [12] Y. Yamada, in: Recent Advances in Matrix Methods of Structural Analysis and Design, ed. R.H. Gallagher et al. (University of Alabama Press, 1971) p. 283. [13] O.C. Zienkiewicz, S. Valliapan and I.P. King, Int. J. Numer. Methods Eng. 1 (1969) 75. [14] Y. Ueda and T. Yamakawa, in: Advances in Computational Methods in Structural Mechanics and Design, ed. J.T, Oden et al. (University of Alabama Press, 1972) pp. 375-392. [15] T.R. Hsu and A.W.M. Bertels, AIAA J. 12 (1974) 1450. [16] R.M. Richard and J.R. Blacklock, AIAA J. 7 (1969) 432. [17] H. Ziegler, Quart. Appl. Math. 17 (1959) 1. [18] H. Armen Jr., G. Isakson and A. Pifko, Int. J. Numer. Methods 2 (1970) 189. [19] G. lsakson, H. Armen Jr. and A. Pifko, Report NASA Cr-803 (1967). [20] T.R. Hsu and A.W.M. Bertels, J. Pressure Vessel Technol., Trans. ASME, Feb. (1976) 17. [ 21 ] Z. Zudans, M.M. Reddi, H.M. Fishman and H.C. Tsai, 2rid SMiRT Conf., Berlin, Germany, Sept. 1973, Paper L4/1 *. [22] Y.R. Rashid, H.T. Tang and E.B. Johansson, Nucl. Eng. Des. 29 (1974) 1. [23] J.H. Gittus, 3rd SMiRT Conf., London, UK, Sept. 1975, Vol. 1, Paper C2/1". [24] E. Krempl, Nucl. Eng. Des. 29 (1974) 125. [25] J. Gittus, Creep Viscoelasticity and Creep Fracture in Solids (John Wiley, New York, 1975). [26] J.A.L. Robertson, Irradiation Effects in Nuclear Fuels (Gordon & Breach, New York/London/Paris, 1969). [27] S. Dushman, Scientific Foundations of Vacuum Technique (John Wiley, New York, 1959).

S. Baner]ee et al. / Analysis of fuel behaviour in transients

355

[28] H. von Ubisch, S. Hall and R. Srivastav, Proc. Int. Conf. Peaceful Uses of Atomic Energy, Geneva, A/Conf. 15, Vol. 7 (1958) pp. 697-700. [29] J.M. Lenoir, W.A. Junk and E.W. Comings, Chem. Eng. Progr. 49 (1953) 539. [30] A.M. Ross and R.L. Stoute, Report AECL-1552 (1962). [31] P.G. Hodge and G.N. White, J. Appl. Mech., Trans. ASME 17 (1950) 180. [32] H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids (Oxford University Press, 1959) p. 219. [33] A. Phillips and R. Kasper, J. Appl. Mech., Trans. ASME 40 (1973) 891. [34] J.J.F. Notley, M.J. Pettigrew and H. Vidal, Report AECL-4072 (1972). [35] E. Kohn, Private communication, April 1976.