Discontinuous finite elements in thermal analysis

NUCLEAR ENGINEERING AND DESIGN 14 (1970) 211-222. NORTH-HOLLANDPUBLISHINGCOMPANY

DISCONTINUOUS FINITE ELEMENTS IN THERMAL ANALYSIS J. J. THOMPSON and P. Y. P. CHEN University of New South Wales, School of Nuclear Engineering, Box 1, P. O. Kensington, N. S. W. 2033, Australia

Received 16 June 1970

Finite element methods using discontinuous, piecewise linear or quadratic, approximating functions, are developed for transient heat flow analysis. The methods are shown to have definite advantages over the ordinary continuous t-mite element methods, particularly as regards the prediction of the spatial temperature moments which are needed for the ordinary finite element methods of thermal stress analysis.

1. Introduction The ordinary finite element method [ 1 - 4 ] for the calculation of temperatures and/or displacements in solids is well established. The essential feature is continuity of the trial functions, with continuous linear (CLFE) and continuous quadratic (CQFE) approximations being most common. However, if continuous finite element methods are contemplated for both the temperature and displacement problems of reactor thermoelasticity, certain aspects of the thermal problem suggest-the need for an alternative finite element approach. In the CLFE method, thermal contact resistances between some elements, necessary for the typical canned fuel element problem, requires the specification of two temperatures at each of the nodal points along the relevant interface. If a large number of interfaces with thermal contact resistances are involved, as in the study of cracked fuel regions, the size of the problem can be significantly increased. Although there is no difficulty in including the appropriate temperature functions associated with contact resistances in the overall functional to be made stationary in the CLFE method, normal heat fluxes across interfaces remain secondary quantities, estimated in some way from the constant fluxes in the adjacent elements, implied by the linear temperature approximation in each element. In coupled dynamic problems, e.g. fuel element

and coolant, surface heat fluxes are of primary importance. It may also be argued that in transient problems, the maintenance of correct heat balances in elements and for the overall finite region is necessary for good accuracy. If finite element methods are used for the diffusion theory equations of reactor analysis, identical in principle to thermal analysis, overall neutron balance becomes an important criterion for good accuracy in the eigenvalue, e.g. effective multiplication constant. A very common way of establishing the finite difference equations in heat flow and neutron diffusion using orthogonal mesh lines, is to integrate over a cell and obtain equations for mean values rather than nodal values. In this process the overall conservation laws are automatically satisfied. All this suggests that a finite element method which aims at computing mean values and satisfying conservation laws, would be useful for a wide class of problems, particularly those involving discontinuities at interfaces, since the duplication of nodal points would be unnecessary. This conclusion is strengthened by the fact that the CLFE method for stress analysis uses linear displacement approximations and hence, in rectangular cartesian coordinates, constant stress and strain approximations over elements. The working equations of CLFE can therefore only include integrals of temperatures over elements, i.e. terms dependent on average element temperatures or zero order moments. However, the CLFE method for temperature leads to

212

J.J.THOMPSON and P.Y.P.CHEN

these mean values being approximated by weighted sums of temperatures at element nodal points, which are not necessarily the best estimates of the true volume averaged values. Thus in simple steady state slab heat flow with a uniform internal source, the arithmetic average of any two adjacent nodal values underestimates the true mean in the region between the nodes. Similarly, if a continuous quadratic finite element (CQFE) approach is used for displacements, then stresses vary linearly in each element in a rectangular cartesian coordinate system, and the equations of thermoelasticity involve both the zero'th and first order spatial moments of the element temperatures. These are not necessarily well predicted by either the CLFE or CQFE methods of thermal analysis. The object of this paper is to develop methods which are capable of predicting the mean and first moments of element temperatures more accurately than by the ordinary finite element methods. To do so it is necessary to admit discontinuous trial functions and discontinuous variations. In the literature, the emphasis on the need to use continuous approximations arises solely from an implicit acceptance of a particularly simple variational functional as mandatory. A more general approach leads to a class of discontinuous finite element methods, of which two particular cases, the discontinuous linear (DLFE) and discontinuous quadratic (DQFE) methods will be discussed.

2. General discontinuous slab elements Consider the simple transient slab heat flow problem defined by the following equations for t>0:

k a2T aT ~x2+Q(x,t)=pC-~-;

O
(1)

kaTax+hL [ T - T L*(t)] =0", x=L

(2)

kaT ~xx-

(3)

h0 [ T - T * 0(t)] = 0 ;

x =0.

Taking Laplace transforms w.r.t, time, a variational principle for this problem can be written:

"The correct solution, which satisfies (1), (2) and (3), makes the functional

.;,--i

Z

2

;-

o

+(½ hL ~2

_ ~,L~,~) + (½ ho ~2 _ ~ 0 ~ ) )

(4)

stationary with respect to arbitrary small variations 67', which are continuous with at least piecewise continuous first derivatives". In (4), T denotes a Laplace transform, s is the Laplace transform variable, and T O and TL refer to slab boundary temperatures. Calculating the first variation of (4) gives:

d2T_pCsTld x 61=OfL 67' I Q+k~-~ +

x=l

+

-k

+h

x--0

which establishes the result. It should be noted that the same result is obtained if Laplace transforms are not taken, the contribution I(T) from the time derivative term is written -pCTT", and 7" is treated as invariant when T is changed. In the following this approach will be used. In the ordinary finite element methods, the real time functional corresponding to (4) is expressed in terms of an approximating, continuous, temperature distribution, and the working equations derived by requiring that the derivatives o f / w i t h respect to the nodal temperatures vanish. The variations 6T are therefore continuous functions. Discontinuous trial functions and discontinuous variations require a modified functional. The starting point is the restatement of the governing equations in the form aT =0, q+k-~x

Q - ~ xaq = pCT '

(5)

and to regard the 2-vector (T, q) as the unknown. The

DISCONTINUOUS FINITE ELEMENTS

O~(x) .an.. T~.,

213

The first variation in I 0 due to arbitrary continuous variations ~iTn and/iqn, in the nth element only, is

hi.

now

Tn(x)

Xn+l[6Tn(Q n .-

Xn

to

qo f

+6q, xo:O

xn

Xr~

Oqn

(_~

aTn\7 +-~-x )Jdx +/iTCh [~ (q+ - q,-+l)]

XN. | = L

Fig. 1. Discontinuous slab elements.

+ 5q+n [T1 (T~n+ 1

appropriate functional, with the restrictions that /iT and/iq still be continuous and that/iT vanish on the boundaries, is now

L[ q2 3T ] I 0 (T. q) = fO "~'~+q ~x + QT - pC1~i" dx

(6)

+aq;

_ _

(T;. -

T~n)]+ 5y~n [!2 (q+- 1 - qn)]

-0l

which establishes the correspondence. I 0 is aaow easily modified to account for the correct boamdary conditions, in terms of heat transfer coefficiem~ h 0 and h L , and generalised to allow for cor~ct resistances betweem elements, such that

since (AT)n+I = T+ - Tn+l =Rn+l qn+l •

L{E

l

The result, which is the basis of discontinuous finite element approach, is as follows:

l(T, q)= ~ d

-j~ + qn "~x + QnTn

n = 0 Xn

Replacing the slab by a number of finite (slab) elements adds to the problem specification, i.e. eqs. (5) in each element, the requirement that the true (T, q) vector be continuous at element boundaries. With the notation shown in fig. 1, the basic functional 10 (T, q), which now permits approximating functions and variations discontinuous at element boundaries, becomes of the form n+l qn

I 0 (T, q) =

N-1

+1 ~

I(q+ +qn+l)(Tn+l_ T +)

n=0

+RN+ 1 1

]

~2

+

7"o) .

1

+ 2

+ 2~O(qO) + q N ( T L - T~V) +~-L(qN) • (8)

"~+ qn -~x + QnT, rl

-ocr, .Jdx

x n

- o C T , T, dx + ~, ~ (q+ + qn+ 1) [ T~n+1 - T~n]. 1"1

(7)

Note that if a boundary temperature is specified, e.g. T(n=0) = TO, it is sufficient to put h 0 ~ 0% TO = TO

in eq. (8).

214

J.J.THOMPSON and P.Y.P.CHEN

It should be noted that arbitrary variations + fiT+ and 6qn at the (n + 1)th boundary give the first variation f~

1

+

Rn+l

[

Rr~!

:Tn

J a

[-~ (qn - qn+l)] + ½ (Sq+ r n + l -

Rn

~Rn~

+

+ - - - ~ (qn + qn+ 1) 1

in I(T, q). Thus continuity of q and the correct temperature discontinuity are necessary for stationary L

Fig. 2. Slab DLFE functions.

~n

qn(~n) = 2 (qn + qn+l) +An (qn+l - qn)

3. Discontinuous linear and quadratic finite slab

-An~2 < ~n < + An~2

elements

From eq. (8), various methods can be developed by choosing different sets of trial functions for T and q. In the following, the requirement that q be continuous at boundaries is imposed, as this defines unique values for contact resistance problems, implies that integral conservation laws can be satisfied exactly, and keeps the formulation simple.

OTn O,

n

--

n=O,...N+l.

~qn

3.1. DLFE formulation (fig. 2) In each element, q varies linearly between nodal points and T is constant. Because q is continuous, the functional I(T,q) given by eq. (8) can be simplified to N

and the equations for the boundary fluxes and element temperatures obtained by writing

It is immediately obvious that +An/2

~I_ 3T n=

which is the correct element heat balance equation in terms of the average temperature. For the general internal point (n ~ O, N + 1)

N - PCTn (~n) l"n (~n)1 d~n + ½ ~ Rn q2 n=l 1

2

•

+aN+t TL +

(10)

--qn+l +qn = 0 ,

Oqn -Tn ( ~ n ) ~ + Tn (~n) Qn (~n) Vgn

-q0T0 +

f Qn(~n)d~n - pCAn 7"n -An~2

+5n/2

I(T'q)=~n=O _ Afn/2 [~k(qn(~n))2

•

=0,

1, . . .N;

1

JAn+An-l]

31

Oqn=Vn--Vn-1 +qn Rn+ qn+lAn 6k

+ - -

aL,

(9)

+

qn-lAn_l 6k

-o.

~"L

Writing +An/2

In the DLFE method,

o ° =A-1 _ f ~'.(~.) = rn ,

- An~2 < ~n < + An~2

an(~n)d~n,

3k (11)

DISCONTINUOUS FINITE ELEMENTS

Rn

Rn°l

i

~'~n~

215

The ~n coefficient in Tn(~n ) has been chosen to make T'n essentially the first spatial moment, i.e.

~rl*l

+an~2

----a<

t

Tn

i

J

A2 _ZXn/2 Tn(~n) ~nd~n.

Substituting (14) into (9) and requiring

Xn-i

Xn

M aI aI ~n=~n-aqn

Xn.i

M 3sn

Fig. 3. Slab DQFE function. i.e. the zero'th order spatial moment, eqs. (10) and (11) together produce the result, needed later, that

Tn - l n - 1 + qn Rn + Q0 A.2 +

6k

--

An +An-11 2k

QO_ 1An2- 1 6k

0 - pCAn'i n -qn+l + qn = 0. QnAn (12)

=0.

for all elements gives the DQFE simultaneous equations. In particular aI/aTn = 0 gives the heat balance equation as before, i.e.

When n = 0, (and similarly when n = N + 1), the equations corresponding to (I0) and (11) are

Also, aI/aT'n : 0 produces the relation •, , 8Sn p C T n = Q n - 12An '

which is equivalent to the use of the quadratic q trial function in the correct first moment equation

Q°Ao - pCA0 J'0 - ql + q0 = 0 , 2 0

•

T o _ T o + q o ! 1 +2-k- + - ~ - - : 0 .

+An/2

(13)

, _!

f

3.2. DQFE formulation (fig. 3) From eq. (9), the functional to be made stationary when q is continuous, it is clear that if qn(~n) is linear, no additional gain is obtained by attempting to make Tn(~n) linear rather than constant. For Tn(~n ) to be linear, at least a quadratic qn(~n) is necessary. This is the basis of the discontinuous, quadratic finite element method. The DQFE trial functions for slab elements are:

4. Two dimensional discontinuous finite elements Only the DLFE formulation will be developed in

vl

12Tn~n An

I¢

7". ( ~ , ) : L , + - -

qn(~n)

=

qn +qn+l + ~n 2 ~ (q"+ 1

\A ~W -

qn)

INTERNAL (jtl) I' 4~n2

1) -

n

.

,c c ~ SIc' hlc' Tie

(14)

AjP

~•

BOUNDARY (I(B)

Fig. 4. Two dimensional DLFE.

216

J.J.THOMPSON and P.Y.P.CHEN

detail as no new principles are involved in the extension to quadratic and three dimensional elements. With the requirement that the heat flux vector be continuous at element boundaries, and that all functions be at least piecewise continuous over volumes or boundaries, the functional I(T, q) to be made stationary can be written

l(T, 7) =

-

+

~ /I~k] ]e(l+B) r ]

qot,] = ( dP], q]ot) = ' "~'jqlot (16)

q [3,] = ( dp] ,

=

@ *;4

The spatial functions ~/.,(F-) etc., are the usual functions appearing in the continuous linear finite element method, relating internal values to nodal values. Similarly, on the boundary ]]',

ql" qj - T/A " qj + Q/T] n..," ]1 q.., ]] = ('I'e~,i/' q]a) + (qs ]],, 4 )

ncjr/~/]dr,.

(17)

where

f

]e(l+B) jew] S]]'

leB Slc

["

•I,a,//

1

ne~,]/,

gZ

eej],

~,]/

11

(~5) The notation is shown in fig. 4. Each element, index ], is either an internal (]el) or a boundary (]eB) element, and has a number of adjacent elements whose indices form the set cop Index c refers to coolant, and T* is bulk coolant temperature, which may vary in space and time. Assuming an orthogonal (c~,/3) coordinate system, the temperature and the components qc~, q~ of the heat flux become, in the DLFE method,

,

11

The functions ff~ etc., are linear functions of distance along tli'~if' boundary, since q..' varies 1 linearly between the Q and R nodes. T~e coefficients n ., and n , are the a and/3 components of the a,# ~j/, . surtace normal vector from I to ]'. Finally, if £a and £# are the two components of the divergence operator appropriate to the coordinate system, .

~

T]A . 7i = TI (£ot~j, q~) + TI

(£~/., ~).

(l 8)

r/(~) -- rJ, In rectangular cartesian coordinates, the vectors

qo~,](Y) = AoU + Botjo~ + Col,i{3,

£ot~j and £#cbj have elements constant over V].. The functional I(T, 7) may now be written in

q#,/(Y) = A~,/ + B#,/v~ + C~,/13 .

terms of a number of basic matrices and vectors, determined by the geometry and thermal properties of elements and surfaces, and defined as follows:

Using the vectors q/'otand ~ of flux components at the nodal points of the ]th element,

aJ: q~ } etc.,)

f ~ t~;~;ld%,

vj

f

vj

xo--fxd , i.e.,

b

DISCONTINUOUSFINITE ELEMENTS

rig'= f R a, t%a, *'.,//1%,,

In the DQFE method, the approximations are of the form

S.3

11

5. (~-) =a/+8:~+ qe = (×/, r/),

..t

.... ,v'or,t1 ..,l R//,[ . t~,11

rig,. = f St

as..,11 '

qa,] (F) =Aaj +Bada + Cad~ +Da,la2

11

r./C= f

T;

217

koo~,icdSlc

+ E,~:2 + F,~,:o¢

Slc 1

=(*.4).

t

"'.¢ = f < I~'.,,c*.,,: ~,o,

The components of the 3-vector Ti are conveniently chosen as the mean, and the first order a and 13 moments about the element centroid. With the reinterpretation of q/a and 4 as 6-vectors of nodal and mid-nodal flux components, and a corresponding increase in the size of the basic matrices, the only terms in eq. (19) to be changed are those associated with the element temperature. These now become

Slc

H~,ol '~

1

t

= f g [%,,c*.,,cl ~,c. Slc

The functional I(T, q) becomes

;(r,-4)= /e(l+B) E {~(4,a%)+~(4,a%)

- (4, ::::) - (4, 5,':) - (e, #) -

E {(¢.,.g'4)+(4,n{¢~)

+¼ E

ie(I+B) 1"eto/

where

A/= fr(e.¢;)xll d~.,

..I

vj

+ 2 (4, .Ii.q'~)} + ~ { ( r ~ ' q t ) + ( r # 'Icq # )I- ¥ t q1 w: leB

(r:,coo/i./)

~= f o.jxjdV/,

I /_/{c l-

c~qa)

_ i1 (q~, ~Cqo)l _ (qfl, l H~x~qo/)} lc 1 .

vi

(19)

The overall coefficient matrix and input vectors (from Q and T*) follow by summation of the individual element components. Thus an internal element with zero contact resistances gives

cocV=f coc)j t×j×)Jd~. v, If the material is anisotropic, i.e.

3I .= O0 _ (pC)ft. i _ (X/, 4 ) _ (~/, q / ) = 0 3T /

/

a, = ~ 4 -

~4

8I=

r;x,.,

or q = KAT in matrix form, with K symmetric, then AT = K- lq _~Eq, and the first term (1/2k) -4 • -4 in the functional I(T, -4) must be replaced by

' ~ ~qa°oqsq~ o~ :3

218

J.J.THOMPSON and P.Y.P.CHEN

5. Assessment of discontinuous slab elements Steady state heat conduction in a slab, using elements of equal width A, provides a simple application of the DLFE method to a problem for which the exact relations can be obtained analytically. With piecewise constant Q in each element, eqs. (10) and (12) lead to the relation k A--g [Tn+ 1 - 2Tn + Tn+ 1 ] + Qn

+ ~ (Qn+ 1 - 2Qn + Q n - 1) = 0

(20)

for the element average temperatures. As expected, since a linear q variation is correct under these conditions, eq. (20) agrees with the exact relation. This is also true when constant resistances and boundary conditions are included. For the same problem, CLFE predicts the exact relation for the nodal temperatures. The CLFE estimate of mean element temperatures, i.e., the arithmetic average of nodal temperatures, will be in error by an amount e n in the nth element. It is easily shown that the error equation at an internal element is

the grounds that (Q - oCT) plays the role of an effective source with rapid spatial variations. It should also be noted that eq. (20) verifies that the DLFE method leads to the correct differential equations as A -+ 0. More generally, an arbitrary Q can always be expressed as a power series in x, so it is of interest to consider heat sources of the form Q = x m, m = 0, 1 , 2 , . . . , with x = 0 at the centre of the nth element. It can be shown that the exact relation between mean element temperatures is then k 42 [Tn+l - 2 T n + T n + l ) = 0 ;

= am ;

A 2 (en+l -- 2e n + e n _ l ) + 1~ [Qn+l - 2Qn + Q n - l ] = 0

where a m = 3Am (3 m+2 -- 1)/2 m+2 (m :+ I) × (m + 2) (m + 3). The DLFE method produces the relation

A2 [Tn+l - 2Tn + T n - 1 ] = 0

= bm ;

m

0 2 4

2hga ] 2k+hLZ ~ eN

F(2k+3hLA.~

m odd

;

m even

(21)

and the error satisfies boundary conditions of the form

+ ~ L \ 2 k + hLA ]ON -

m even

where b m = A m (3 m + 1)/2 m+l (m + 1). A comparison of the numerical values gives the following results:

!

k l ,52 ( e N - e N - 1 ) +

m odd

QN_I]= 0

(22)

From these error equations, it is clear that the greatest difference between the DLFE and CLFE mean temperature estimates occurs when the heat source is rapidly varying. This occurs at an interface between the fuel and can regions of a reactor fuel element. It might also be anticipated that significant differences will occur in transient problems, on

Exact

DLFE

(am)

(bm)

1 A2/4 13A4/80

l 5A2/12 31A4/80

Thus the DLFE method as formulated, with a constant T and linear q in each element, produces small errors when Q has significant quadratic or higher order (even) components. It is interesting to note that some improvement in accuracy in the static problem can be achieved by using a modified DLFE method, in which the temperature variation in each cell is assumed linear, of the form

DISCONTINUOUS FINITE ELEMENTS (23)

Tn (~n)= Tn - ~ - ~ ( q n + q n + l )

i.e. the heat fluxes from T and q are equated at the centre• In two dimensions the centroid could be used. The modified DLFE method, with Q = x m leads to the relation k A 2 [ r n + 1 - 2T n 'I- T n _ l ] = 0 ;

m odd

= Cm ,

m even

where C m = Am [_3 m (m - 7) + m + 1] / 2 m+2 (m + 1) (m + 2). For the first few even values of m, CO = 1 ,

C2 = A2/4,

C 4 = 31A4/240

and hence correct mean temperatures are predicted for quadratic heat sources• The CLFE method can be shown to predict the correct nodal temperatures for any value of m, provided that the source term in the functional is treated exactly• This is also true for nonuniform element widths, as verified by numerical studies of a slab of width L = 1.2, with Q = k x m , zero heat flux at x = 0 and h L / k = 1.0. The slab was represented by 12 elements, with widths defined by An/An+ 1 = 1.5. Although the correct nodal temperatures were obtained, the CLFE mean element temperatures, i.e. the arithmetic mean o f adjacent nodal temperatures, were in error• A comparison between these CLFE errors, and those obtained by using the modified DLFE method for the first and widest element is as follows:

% Error in mean temperature

219

In applying the modified DLFE to transient problems, it is important to note that it is not valid to represent T by the time derivative of the approximation given by (23), as this leads to instability• The correct approximation is 7"n (~n) = l"n"

(24)

To assess analytically the relative virtues o f the discontinuous and continuous element methods in transient heat flow, it is necessary to choose a realistic slab problem with known analytical solution and to minimise the number of elements• Such a problem, and a severe test of the finite element method, is the symmetric slab quench - surfaces cooled instantaneously from the initial uniform temperature T O to zero - with only one element for the half slab. The governing equation can be written

~2T(x'O)-aT(x'O)'" 3x2 30 '

O~
0>0

with q(0, 0) and T ( 1 , 0 ) = 0. The exact solution is 0°

1

T(x, 0 ) = T O ~ q~i(x) f qbi(x ) dx e -Ki° i=0 0 with

~i (x) = vTcos (xy), Xi = (2i + 1) 7r/2.

•i = X2 ;

The approximating functions for the continuous and discontinuous, linear and quadratic finite element methods are as follows: CLFE:

T(x, 0 ) = TI(0)(1 - x )

DLFE:

T(x, O) = T(O) ,

CQFE:

T(x, O)= TI(O)(1 - 3x + 2 x 2)

m

CLFE

DLFE

1

2 X 10 -1

-

2

4 X 10-2

0.5 X 10-2

3

1 X 10 -2

0.3 X 10-2

4

0.3 X 10-2

0.1 X 10-2

+

r 2 ( o ) (4x

-

q(x, O) = q(O) x

4x 2)

220

J.J.THOMPSON and P.Y.P.CHEN Table 1 Symmetric slab quench analyses.

Method

Differential operator

Mean temperature Tave/T 0

First moment _ T' /T 0

CLFE

D + col col = 3.0

0.5

0.083 e -c° 10

DLFE

D+col col =3.0

1.0 e -colO

0

CQFE

(D + col) (D + co2) wl = 2.486 co2 = 32.181

0.8058 e -col0 + 0.0275 e -co20

0.1053 e -col0 - 0.0220 e w20

DQFE

(D + CO1) (D + co2) co1 = 2.486 co2 = 32.181

0.8143 e -col0 + 0 . 1 8 5 7 e -w20

0.1121 (e -col0

0.81057 e- 2.4674 0

0.1 1074

Exact-Asymptotic

f T(x, O) = T(O) - T'(O) (12x - 6)

DQFE: q(x, O) = q(O) (2x 2 - x ) + q l ( 0 ) (4x - 4 x 2 ) .

The working equations o f the four methods are sufficiently simple to be solved analytically, and the main results are shown in table 1, using zero and first order moments defined b y 1

Tave=f 0

1

T':f0

It is assumed that thermal stress analysis would use a continuous method, so Tave is significant for a CLFE displacement model, while both Tave and T' are significant for a CQFE approach. Table 1 shows that the same order approximation, i.e. linear or quadratic, leads to exactly the same characteristic equation whether the continuous or discontinuous approach is adopted. In the linear models, the mean temperature history is overestimated for most of the time by the discontinuous method, and underestimated consistently by the usual continuous method. Of the two, DLFE is more accurate. The quadratic methods overestimate the fundamental mode decay constant by 0.75%. The

e -colO

e -2"4674

e -t~20)

-

0

amplitude (0.8106) o f the fundamental mode is in error b y + 0.45% (DQFE) and - 0 . 5 9 % ( D L F E ) respectively. However the predicted time histories are quite different in the continuous and discontinuous methods, for small values o f 0. Table 2 gives a comparison between CQFE, DQFE and the exact solution. It is clear that the discontinuous method represents a significant improvement over the usual continuous quadratic finite element method. The results of this slab quench application support the contention that discontinuous methods result in greater accuracy in the prediction of the zero and first order moments o f element temperatures in transient problems.

6. Discussion By treating temperature and heat flux components as separate functions, a general discontinuous finite element formulation has been achieved. Two particular methods, the discontinuous linear and discontinuous quadratic finite element methods, in which only the temperature is discontinuous, have been shown to have definite advantages over the usual continuous methods, where overall conservation laws and temperature spatial moments are important. The proposed methods appear to be the

221

DISCONTINUOUS FINITE ELEMENTS Table 2 Symmetric slab quench results.

0 0 0.02 0.04 0.06 0.08 0.10 0.15 0.20

Tave/T0

-T'/To

Exact

CQFE

DQFE

1.000 0.840 0.775 0.724 0.681 0.643 0.563 0.496

0.833 (-16.7%) 0.781 ( - 7.0%) 0.737 ( - 4.9%) 0.698 ( - 3.6%) 0.663 ( - 3.6%) 0.630 ( - 2.0%) 0.555 ( - 1.4%) 0.490 ( - 1.2%)

1.000 0.873 0.789 0.728 0.682 0.643 0.563 0.496

(+ 3.9%) (+ 1.8%) (+ 0.5%) (+ 0.1%) (-) (-) (-)

"natural" methods for calculating temperatures to be used in continuous finite element analysis o f thermal stress, particularly in transient problems with internal heat sources and contact resistances, and therefore in problems of nuclear structural analysis. In the variational formulation of the static problem, the appropriate functional is no longer an extremum, as it is when temperature only is used in the continuous finite element methods. However, although a variational principle has been invoked in developing the methods, it is clear from the form o f the first variation that the interpretation o f this discontinuous function approach as a particular application o f Galerkin's method is more appropriate, On this basis, considerable scope exists for the development of discontinuous function methods specifically designed for different coordinate systems and particular problems. For this reason extensive numerical studies o f the two proposed DLFE and DQFE methods, as formulated in this paper, have not been included, The application to realistic fuel element analysis will involve temperature dependent conductivity (k) and spatially dependent heat transfer coefficient (h). A characteristic feature o f the proposed discontinuous approximation is the occurrence of these functions as k - 1 and h - 1, not as the simple k and h multipliers of the continuous method. Only numerical studies of specific problems will reveal the effect of this aspect of discontinuous finite element methods. Also further work is required to find the optimum numerical techniques for solving the sets o f equations generated in two or three dimensional problems.

Exact

CQFE

0 0.0598 0.0728 0.0783 0.0800 0.0795 0.0742 0.0669

0.0833 0.0886 0.0893 0.0875 0.0846 0.0813 0.0723 0.0643

DQFE (oo) (+ 48.1%) (+ 22.7%) (+ 11.7%) (+ 5.8%) (+ 2.3%) (+ 2.6%) ( - 3.9%)

0 0.0478 (-20.0%) 0.0705 ( - 3.2%) 0.0803 (+ 2.6%) 0.0834 (+ 4.3%) 0.0830 (+ 4.4%) 0.0763 (+ 2.8%) 0.0680 (+ 1.6%)

7. Notation C h

I(T, q) a,/3 k L q, Q R s S t T T T* V x a,/3 A e 0

O

= heat capacity = heat transfer coefficient = functional of variational principle = component o f divergence operator in two dimensions = thermal conductivity = slab thickness = unit outward normal from a surface = one dimensional and vector heat flux = heat source density = contact resistance at interface = Laplace transform variable = surface = time = temperature = time derivative of T(= aT/at) = bulk coolant temperature = volume of finite element = spatial coordinate for slab = orthogonal coordinates -- thickness o f slab element = error = dimensionless time (= kt/pCL 2) = slab element spatial coordinate, origin at centre = material density

References [1] E.L.Wilson and R.E.Nickell, Application ofthe finite element method to heat conduction analysis, Nucl. Eng. Desgn 4 (1966) 276-286.

222

J.J.THOMPSON and P.Y.P.CHEN

[2] O.C.Zienkiewicz, The finite element method in structural and continuum mechanics (McGraw-Hill, New York, 1967). [3] P.Launay, G.Charpenet and C.Vouillon, Prestressed concrete pressure vessels: Two dimensional thermo-

elasticity computer program, Nucl. Eng. Design 8 (1968) 443-459. [4] A.F.Emery and W.W.Carson, An evaluation of the use of the finite element method in the computation of temperature, SCL-RR-69-83 (August 1969).