Modeling of the transient heat transfer in a nuclear reactor fuel rod using a variational procedure

Nuclear Engineering and Design 51 (1979) 133-142 © North-Holland Publishing Company

133

MODELING OF THE TRANSIENT HEAT TRANSFER IN A NUCLEAR REACTOR FUEL ROD USING A VARIATIONAL PROCEDURE G. LEBON Department of Mechanics, University of Liege, Building B6, Sart Tilman, B-4000 Liege, Belgium Ph. MATHIEU Department of Nuclear Engineering, University of Liege, Val Benoft, B-4000 Liege, Belgium and J. VAN VLIET S.A. Belgonuclbaire, 25, rue du Champ de Mars, B-1050 Brussels, Belgium Received 4 March 1978

The transient temperature field in a nuclear fuel rod consisting of fuel cylindrical pellets surrounded by a metallic sheath is calculated. This is done by using a variational technique based on Lebon-Lambermont's variational principle. The thermal properties of the system are assumed to be temperature dependent, while the boundary conditions are given as functions of time.

1. Introduction The study of transient thermomechanical phenomena in nuclear reactor fuel rods is of great importance because it is directly related to the problem of reactor safety. Transitory effects are not only met in normal operations (for instance, during load-following) but also in accident phases [for instance, during a loss of coolant accident (LOCA)]. It is of crucial importance to be sure that, during these operations, the fuel rod keeps its integrity, which means that neither cladding failures nor gross fuel melting occur. In that respect, it is necessary to know the transient temperature profile in the fuel rod. Once the temperature field is determined, one can compute the stress field as well as the transient structural changes and geometry of both pellets and cladding. Of course, these effects influence, in turn, the temperature distribution which has to be re-evaluated [1]. This paper is concerned with the calculation of the transient temperature field only; mechanical and structural changes are not investigated here. Although the latter effects are ignored, the problem remains intricate because of its highly non-linear character. Non-linearities arise from the dependence of the boundary conditions and the thermal properties of the system with the temperature. The so~rch for an analytic solution is excluded. The purpose of this work is to provide an approximate solution by using a variational technique based on Lebon-Lambermont's variational principle [2,3]. This criterion appears to be particularly well suited for handling problems of heat transfer [4] and fluid mechanics [5]. Sect. 2 is devoted to the description of the system. The basic hypotheses underlying the analysis are stated and the set of equations describing the transient heat transfer problem is established. In sect. 3 the corresponding expression for Lebon-Lambermont's variational principle is formulated. In sect. 4 an approximate solution is ob-

G. Lebon et al. / Transient heat transfer in a nuclear reactor f~tel rod

134

tained by using Kantorovitch's partial integration method. In sect. 5 the formalism is particularized to calculate the temperature profile in a nuclear reactor fuel rod.

2. Basic assumptions and equations A nuclear reactor fuel rod consists of a metallic zircaloy or stainless steel cladding surrounding a stack of cylindrical uranium or uranium-plutonium oxide pellets. Under high temperature and neutron irradiation, these materials undergo irreversible thermal, mechanical and structural changes. These phenomena interact deeply with each other. Moreover, since the treatment of the whole coupled problem is a complicated and costly task, the thermal and the mechanical problems are generally uncoupled and solved by an iterative procedure until convergence, i.e. full consistency between the thermal and mechanical fields, is obtained. Here, we concentrate on the heat transfer problem exclusively. The foregoing analysis rests on the following hypotheses. (1) Any axial heat transfer is neglected; this is justified because the length of the rod is much greater than its diameter and also because the pellets' interfaces offer a high thermal resistance. (2) The only active heat transfer process is conduction: convection due to a gas flow through the pellet cracks is neglected. This is a good assertion because neither the gas amount, nor its flow speed reach the level required to modify appreciably the temperature field. (3) The strain effects on the temperature field in the fuel is taken into account by allowing the fuel thermal conductivity to depend on the porosity. (4) The gap heat transfer coefficient h F which generally depends on the gap width, the temperature at the fuel outer surface and at the cladding inner surface, the inner gas pressure, and the mean temperature, is modeled by a given function of time. The latter can also include effects arising from radiation. (5) In the same way, the heat transfer coefficient h G in the film between cladding and coolant is also approximated by a given function of time. (6) The fuel pellets and the cladding receive separate descriptions so that their temperature distributions are obtained separately and independently. The reasons for adopting an uncoupled description are the following. First, it is believed that the computing time will be appreciably reduced [1 ]. Moreover, this uncoupling allows the use of the same subprogramme in both the fuel and the cladding. Finally, the thermal properties of the fuel rod are given within a certain error percentage which is reflected in the numerical results; therefore, it would be illusory to search for a highly accurate description. The temperature distribution inside the fuel and the cladding is found by solving the set of transient non-linear equations

p i ( r , T ) c i ( r , T ) - ~0T i ( r , t ) = r ~

10r Er k i ( r , T ) ~ 0r T i ( r , t

+W(r,t),

i=F,G,

(2.1)

where the superscripts F and G refer respectively to the fuel and the cladding. In (2.1)p represents the density, c the heat capacity, k the heat conductivity, and W the volumic internal heat source term. The above parabolic partial differential equations must be completedby the initial conditions

Ti(r, O) = T~o(r),

i = F, G ,

(2.2)

and the following boundary conditions; at r =R[n :

qi. n i = _qi(t),

(2.3)

at r-- Rou ti :

qi. ni=hi(t)[Ti(r, t ) - Tie(r, t)],

(2.4)

where

Rin and Rou t denote respectively the inner and outer radii of the pin components, n is the unit normal vec-

G. Lebon et al. / Transient heat transfer in a nuclear reactor fuel rod

135

tor, positively oriented outwards, q is the heat flux vector which according to Fourier's law, reads as 0T q = -k ~rer,

(2.5)

in which e r is the unit vector in the radial direction. Eq. (2.3) shows that the radial heat flux is a prescribed function of time qi(t) Since the scalar product is nonpositive at the inner face, qi. ni has been affected by a minus sign in order to deal only with positive heat flux densities. Relation (2.4) is the well-known Newton law of cooling, wherein h i is the heat transfer coefficient between the fuel and the cladding (i = F) and between the cladding and the coolant (i = G), respectively. In this work, h i is an assigned function of time. The quantity Tie(t) denotes an "external" temperature; when i refers to the fuel, Tff represents the temperature at the inner face of the cladding; when i refers to the cladding, TG is the temperature of the coolant. Since the quantities qi(t) and Tie(t) are prescribed functions of the time, it is clear that the problem described by the above set of equations (2.1)-(2.4)is an uncoupled one, according to hypothesis (6).

3. Lebon-Lambermont's variational principle The problem is to find a variational equation of the form

L T, Or Ot ;r,t

8I(T)-8 t

dtdr=O,

(3.1)

r

such that one recovers, in Euler-Lagrange form, the heat equation (2.1) together with the appropriate boundary and initial conditions. It is well known [6] that there exists no "true" variational principle equivalent to eqs. (2.1)(2.4). A true variational criterion is such that all the quantities appearing in the function L are submitted to variation. Only restricted principles, wherein some quantities are frozen during the variational procedure, are able to restore the above set of equations. Many principles, like the Hamilton principle in analytical mechanics, the Glansdorff-Prigogine local potential [7], the Biot [8] or Lebon-Lambermont [2,3] variational criteria, must be classified in this category. In spite of their restricted character, they have been applied successfully to several problems in heat transfer, solid and fluid mechanics. For a prescribed surface temperature and in the absence of internal heat supply, Lebon-Lambermont's principle [4] can be written as 6t ff V

T 2(gv[T-1]+½o) d V d t = O .

(3.2)

t

An upper dot denotes derivation with respect to time. The subscript t affecting the variational operator 8 means that the time derivative of the variable submitted to variation, namely the temperature, has to be held fixed during variation. The quantity Sv [ T-1 ] is the Legendre transform of the entropy sv, referred to per unit volume. Denoting the internal energy per unit volume by Uv, the Legendre transform SvIT -a ] is expressed as sv[T] =Sv - T-auv = - ( 1 / T ) f v .

(3.3)

Its time derivative is

Sv [T-1 ] = -Uv ~-1 = (1/T 2) uvJ" .

(3.4)

In (3.2) o represents the local entropy production per unit volume, given by [9]

° = T-~-k(T) (grad T)2 .

(3.5)

G. Lebon et al. / Transient heat transfer in a nuclear reactor fuel rod

136

Substituting (3.4) and (3.5) into (3.2) and using the state equation

Uv =f pc dT, the variational principle (3.2) takes the following form: T

fit f f [ J ' . f V t

p(a)c(a)da+ lk(T)(grad T)2]dVdt=O.

(3.6)

0

If the heat conductivity is temperature dependent, it must also be frozen during the variational procedure. With this in mind, it is checked that the necessary condition for (3.6) to hold true is

peT"- div [k(T)grad T] = 0 , which is the equation for heat conduction in the absence of internal heat supply. In the presence of a temperature dependent heat source, a supplementary term must be added to (3.6) which takes the form T

T

f f LTf V t

a. +½k(T)(grad

T)2 -

0

f w(.)d.lav d t = 0 .

(3.7)

0

The above form is only valid for a prescribed surface temperature. For other boundary conditions, e.g. a prescribed flux or a cooling law at the surface, expression (3.7) must be supplemented by surface integrals. In particular, the variational criterion appropriate to the description of the heat transfer in the fuel rod elements takes the form t

fit

'f Rin

oi(r, o0 ci(r, a) da + lki(r, T i)

0

× h i [Ti(r, t) - Tie(r, t)] 2 IRout d t - Rin

Ti(r, t)

- f 0

t } f qi(t) Ti(r, t)lRin dt : 0 .

W(r, or)da r dr dt + IRou t

f 0

(3.8)

0

The surface integrals have to be calculated at the outer and inner surfaces, respectively; the upper limit in the time integral may be chosen arbitrarily. It is easy to verify that the Euler-Lagrange equation corresponding to (3.8) is the heat conduction equation (2.1) and that one recovers also the correct boundary conditions (2.3) and (2.4).

4. The variational Kantorovitch method

To obtain approximate solutions of the temperature distribution in both the fuel core and the sheath, Kantorovitch's partial integration method will be used. It consists in separating the variables r and t in the trial function, namely N

T(r, t)= ~ fn(r)gn(t),

(4.1)

n=l

where the fn (r) are a priori given functions of the radial variable while gn (t) are unknown functions of the time to

G. Lebon et al. / Transientheat transfer in a nuclearreactorfuel rod

137

be determined by the variational technique. After substitution of N

8r =

(4.2)

fn (r) 8gn(t)

n=l

in the expression of the variational principle (3.8), one is led to a relation of the form t

N

N

f dt ~

m=l

0

~out

8gm{ ~ 8n n=l

N

;

Pcfnfm r dr + ~ gn n=l

Rin

9

out Wfmrdr+R°uth[fm( ~ fngn n=l Rin

out k f n f m r dr Rin

Te)]Rout

RinqOem)Rin} = 0 ,

(4.3)

where fn = dfn (r)/dr. Since (4.3) must hold for arbitrary variations 6gm, it follows that the set of differential equations between accolades must vanish identically. Putting Rout

Anm = f

Pcfnfm r dr,

(4.4)

kl',cf~,rdr + RouthOenfm)Rout ,

(4.5)

Wfmr dr + RouthTe(fm)Rout

(4.6)

Rin Rout

Bnm = f Rin

tom =

~

out + RingOem)Rin ,

Rin one obtains the following set of first order differential equations: N n=l

(.4 nm[~n + B n m g n ) = ¢°m ,

m = 1,2 . . . . .

N

(4.7)

or, in a more compact notation, A'g+B.g=oJ,

(4.8)

where A and B are symmetric N X N matrices while to is an N order vector. Solving (4.8) with respect to 8 yields g = - D .g + d ,

(4.9)

with D=-A -a'B

and

d = A -1. to.

(4.10)

Generally (4.9) is a non-linear vectorial differential equation because D and d depend themselves on g. However, when all the thermal properties as well as I¢ are independent of the temperature, eq. (4.9) becomes a linear differential equation; its solution is given by t g = e x p ( - D t ) [ f exp(Dt) d(t) dt + g ( 0 ) ] . o

(4.11)

138

G. Lebon et al. / Transient heat transfer in a nuclear reactor fuel rod

If the source term W is constant with respect to time, (4.11) reduces to g = [I - exp(-Dt)] D-1 • d +g(0) e x p ( - D t ) .

(4.12)

The initial values of g(0) are determined by ensuring that the trial function (4.1), calculated at t = 0, fits best the initial value, i.e. N

fn(r) qn(O) "~ T(r, 0).

(4.13)

n=l

The approximated solution is derived by substitution of the solutions gn into expression (4.1) of the trial functions. The procedure must be repeated for both the fuel and the cladding.

5. Numerical application As an illustration of the method proposed in sect. 4, one calculates the temperature profile in both the fuel and the cladding of a pin, whose geometrical and thermal properties are reported in appendix A. The initial state is obtained by solving a steady state heat transfer problem with a prescribed source and external temperature profile, as reported in appendix B. The transitory conditions for both the fuel and the cladding are specified in appendix C. They roughly simulate a reactor shutdown with perturbed cooling. F In order to simplify somewhat the boundary conditions at Rout, the external temperature TF, which is nothing but the temperature at the inner face of the cladding, is taken as constant and equal to the initial value TG (R G, 0). Similarly, the heat flux qG (t) through the inner face of the cladding is assumed to remain constant and equal to its initial value qG(0). This is the simplest way to uncouple the problem as required by hypothesis (6) of sect. 2. As a trial function, one takes a simple two-parameter parabolic function (5.1)

r(r, t) = y ~ ( r ) g ~ q ) + f 2 ( r ) g 2 ( t ) ,

where

fl(r)

-4\Rou Rou:12 t _Rin]

(5.2)

while f2 (r) = 1 - f , (r).

(5.3)

The motivations for choosing this expression are twofold. First, it has been shown [10] to fit accurately the exact solution of some linear transient problems solved in Carslaw and Jaeger [ 11 ]. Secondly, it provides a simple and meaningful trial function since g~ and g2 represents the temperature values at the inner and outer radii of the fuel pin components. With (5.1) as the trial function, the differential system takes the form

-L2 +L4

- L o + 2L2 - L4

~'2

4Q

- 4 Q - Routh

g2

- M o +M2 +RouthTc

G. Lebon et al. / Transient heat transfer in a nuclear reactor fuel rod

139

where the matrix elements are given by 1 Rout L / - (Rou t - Rin) / J pc(r - Rin~r dr,

j = 0, 2, 4,

(5.5)

j = 0, 2,

(5.6)

Rin 1

M j - (Rout - Rin) /

Rout

J

w(r- Rin) / r dr,

Rin

1 Q = (Rou t - ,JOin~4

.fl°ut

(5.7)

k(r - Rin) 2 r dr.

Rin

To solve the differential system (5.4) with a minimum of computer time, we shall replace the non-linear heat conduction problem by a succession of linearized problems. This is done by dividing the time interval 0 - t into several micro-steps the duration of which is determined by the inverse of the greatest eigenvalue of matrix D. During each of these micro-steps all the thermal properties are assumed temperature independent and equal to their value at the end of the previous step. In order that the solution (4.12) holds in every micro-step, we shall take as the initial value for T the value calculated at the end of the previous micro-step. The schema is explicit and hence noniterative, and can be written as

g(t + At) = { I -- exp [-D(t) At]} D -1 (t)" d(t) + g(t) exp [-O(t) At].

(5.8)

The results concerning the fuel are presented in fig. 1, where the temperature at the inner and outer faces are reported as a function of time. The stationary profile is attained after about 3.5 s, i.e. about 2.5 s after both boundary conditions and volumic heat source remain constant. This is in agreement with the value of the temperature relaxation time in a homogeneous hollow cylinder with an identical mean thermal diffusivity. Similarly, in fig. 2 are presented the temperature at the inner and outer faces of the cladding. It is seen that the stationary state is reached much faster than in the fuel; this is due to the smaller thickness of the cladding and its

T(K)

T(K)

9OO

2250

2000 8O0

175075O

ql 7OO

1250-

o

q2

~

~

65O

~ t (S)

Fig. 1. Temperature variation at the inner (gl) and outer (g2) faces of the fuel as a function of time.

o.o

o.~s

o.~

o.~

~.b

&s tCs)

Fig. 2. Temperature variation at the inner (gl) and outer (g2) faces of the cladding as a function of time.

G. Lebon et al. / Transient heat transfer in a nuclear reactor fuel rod

140

much better conductivity. The stationary state is attained almost instantaneously after that the boundary conditions are frozen. It follows that during the greatest part of the process the cladding is in steady state.

6. Conclusions A rather general procedure for solving a non-linear problem of transient heat conduction in a fuel-cladding nuclear pin has been proposed. It rests on L e b o n - L a m b e r m o n t ' s variational principle and uses Kantorovitch's partial integration method: it consists in separating in the trial function the contributions arising from the radial and the time variables, respectively. The trial function has been given the general form N

T(r, t) = ~

fn (r) gn (t),

(6.1)

n=l

where fn (r) are given functions of the radial variable and gn (t) are a set of unknowns to be determined by the variational technique. The analysis in sect. 4 has not been restricted to any particular choice offn or to a restricted number of terms in the expansion (6.1). The choice (6.1) yields a set of first order non-linear coupled differential equations with the time as the only independent variable. Of course, by using more and more terms in the expansion (6.1), it is reasonable to expect that the quality of the solution is improved. However, such improvement is at the expense of a more and more complicated and costly numerical treatment. In the numerical illustration of the method, we selected a rather simple trial function involving only two parameters, gl and g2, and wherein fn (r) is a parabolic function. This choice has been justified by the good results which it furnishes when used to solve analogous linear problems whose analytic solution is known. In this paper we have considered a system submitted to rather restrictive hypotheses. However, there is no difficulty in extending the above treatment to more realistic conditions corresponding to a more accurate simulation of the fuel rod behaviour. This is the case if the coupling between fuel and cladding is taken into account at the end of each time-step. The generalization to systems wherein either fuel melting or cladding phase change (or melting) occur is also straightforward.

Appendix A: Geometrical and thermal properties

A. 1. Fuel pellet (uranium oxide) Inner radius (cm): RiVn = 0.0229. Outer radius (cm): RoFut -- 0.2121. Cold theoretical density (g/cm3): p F = 10.983. Heat capacity [12](J/gK): c v = 0.194415 + 2.632 10 -4 T - 1.811 10 -7 T 2 + 4.75 10 -11 T 3. Heat conductivity [13] (W/cmK): k r = ko(T)f(~), where ~ is the porosity which is assumed to be time independent, but which varies radially inside the pellet (see table 1), k0 is the thermal conductivity of a 100% dense fuel: 40.05 ko (T) - 7. + 129.4 + 0.6416 10 -12 T 3 - 0.002;f(~) = 1 - 1.029~ - 0.320 10-3 ~2 _ 0.401 10-4 ~3 + 0.158 10-5 ~4.

A.2. Cladding (AIS1304 stainless steel) Inner radius (cm): Ri~n = 0.2121. Outer radius (cm): RoGut = 0.2527. Density (g/cma): p G = 7.95(1 + 1.7 × 10 -s T~

G. Lebon et al. / Transient heat transfer in a nuclear reactor fuel rod

141

Table 1 Radial distribution in the fuel of the initial temperature distribution TF(r, 0), the porosity ~, and the radial part of the internal

heat source ~ (r).

r (cm)

TF(r, O) (K)

~ (%)

~ (r) (W/cm3)

0.2121 0.2015 0.1902 0.1781 0.1650 0.1496 0.1322 0.1152 0.0951 0.0692 0.0229

1290.03 1413.27 1535.18 1658.1 1784.78 1946.42 2102.96 2192.52 2274.46 2347.80 2402.55

7.3 6.9 7.2 8.6 11.7 20.1 2 2 2 2 2

4465 4030 3587 3131 2602 2304 2369 2976 1808 1560

Heat capacity: c t~ = 0.527 35 + 9.56 10-S(T - 273) + 7.425 10-1°(T - 743) 3. Heat conductivity: k G = 0.146(0.640 + 9.65 10-4T).

Appendix B: Initial conditions B. 1. Fuel pellet Initial temperature profile: displayed in table 1. Internal heat source (W/cma): WF = O(r) ¢(t). if(r) is given in table 1, while ¢(0) = 1. Gap heat transfer coefficient (W/cm2K): hF(0) --- 0.564. No inner heat flux.

B.2. Cladding Initial temperature profile: (K): TG(r, 0) : 57(r - RG)/(RCout - RiGn) + 753. Film heat transfer coefficient: hG(0) --- 15. Inner heat flux (W/cm2): qG(0) = 302.87. Coolant temperature: 679.

Appendix C: Transistory conditions C.1. Fuel pellet Internal heat source:¢(t) = 1 - 0.85 t, t < 1 s, ¢(t) =0.15, t > 1 s. Gap heat transfer coefficient: hF(t) = 0.564 -- 0.464 t, hF(t) =0.1, Clad inner temperature: TC(R~n, t) = TG(R~, 0) = 735. No inner heat flux.

tls.

142

G. Lebon et al. / Transient heat transfer in a nuclear reactor fuel rod

C 2. Cladding

Film heat transfer coefficient: h 6 (t) = 15 - 3t, t~ls. Inner heat flux: q G ( t ) = qG(0) = 302.87. Coolant temperature: T¢ = 679 + 17t, t ~< 1 s, =696 t > 1 s. No internal heat source.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

P. Verbeeck and N. Hoppe, Belgonucl6aire Non Proprietary Report, BN 7609-01 (1976). G. Lebon and J. Lambermont, J. Chem. Phys. 59 (1973) 2929. G. Lebon and J. Lambermont, Ann. Phys. 32 (1974) 425. G. Lebon and J. Casas, J. Eng. Math. 8 (1974) 31. G. Lebon, J. Casas and D. Jou, Appl. Sci. Res. 32 (1976) 371. B. Finlayson, The Method of Weighted Residuals (Wiley, New York, 1971). P. Glansdorff and I. Prigogine, Structure, Stability and Fluctuations (Wiley, New York, 1971). M. Blot, Variational Principles in Heat Transfer (Oxford Math. Mono., Oxford University Press, 1970). S.R. de Groot and P. Mazur, Non Equilibrium Thermodynamics (North-Holland, Amsterdam, 1962). G. Lebon and Ph. Mathieu, to be published. H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids (Clarendon Press, Oxford, 1959). J.F. Kerrisk and D.G. Clifton, Nucl. Technol. 16 (1972) 531. S.Y. Ogawa, E.A. Lees and M.F. Lyons, GEAP-5591 Internal Report (1968).