The transient heat transfer analysis of solids with radiative boundary condition using finite element analysis

INT. CCt~4. H ~ % T M A S S T R A N S F E R 0735-1933/85 $3.00 + .00 VO1. 12, pp. 169-178, 1985 ©Pergamon Press Ltd. Printed in the United States

THE TRANSIENT HEAT TRANSFER ANALYSIS OF SOLIDS WITH RADIATIVE BOUNDARY CONDITION USING FINITE ELEMENT ANALYSIS

B h i m a v a r a p u S. Reddy Graduate Student

Anand M. Sharan Associate Professor

Faculty of Engineering Memorial Unlv~rsity of Newfoundland St. John's, Newfoundland, Canada AIB 3X5

(C~,~nicated by J.P. Hartnett and W.J. Minkowycz)

KBSTRACT I n t h i s work, t h e s y s t e m o f e q u a t i o n s t o o b t a i n t h e t e m p e r a t u r e distribution within a solid are derived using the finite element analysis. These nonlinear set of algebraic equations are solved using the Gauss-Seidel iteration technique. The r e s u l t s indicate t h a t t h i s m e t h o d c a n be a v a l u a b l e t o o l i n t h e s t u d y o f h e a t transfer process in the metallurgical or other related industries.

Introduction The heat transfer process in some of the metallurgical processes is quite involved [1-3]; for example, during the cooling of castings or heating of ingots before forging. complicated shapes.

These castings or ingots can be of very

Therefore, the solution of heat transfer problems by

exact methods is n o t possible.

In such situations, the heat transfer

process is studied either by finite difference or finite element method. The heat transfer process in this problem involves all the three modes of heat transfer which are: the conduction, convection and radiation.

In this

paper, the equations for the heat transfer process of a solid subjected to

169

170

B.S. Reddy and A.M. Sharan

Vol. 12, No.

nonlinear boundary conditions using the finite element analysis have been derived.

Then, these equations are solved using the Gauss-Seldel iteration

technlque. The Mathematical Formulation The governing dlferentlal equation for heat-conduction in solids is

~)2T ~2T ~)2T ~I' kx - ~ + ky ~ + kz --~,2+ Q ~ ~c ~ with

(I)

boundary conditions T ffiT B on Sl,

(2)

and kx ~

~x + ky ~

£y + k z

-~-~- %

h(T-T®~" o e(T4-'l' ~ ) ffi 0

+ q +

on S2

(3) The union of both the surfaces S I and S 2 forms the complete boundary of the solid having a Volume V.

The functional formulation that is

equivalent to Eqn. (i) and its boundary conditions, Eqns.

(2) and (3), can

be written as 1 x(~) ~r 2 +ky (~v)2+kz(~z)2_2QT x- I g[k r

+ 2PcT ~ ] d V

V

+ f[qT+ gh (T_T~2 + oetT5 _ T 4T)]dS

(4)

S Since T is not continuous over the region, rather it is defined over the individual elements, Eqn. (4) can be written as the summation over

individual elements.

Defining the elemental matrices [4]

T e ffi [Ne ] {T }

~r e

{ge}T

(5a)

~{T}

~e

(5 b)

~fe ~re (5 c)

(5d)

0 k ~ Z

Vol. 12, No. 2

[Be]

SOLIDS WITH R A D I A T I V E ~ "~Nie ~x

~Nj_~ e 3x

~Nk e" ~x

8NI e

~Nj e

~Nk e

~y

~y

~Y

~Nie

~Nje

~Nke

~z

3z

=

3z

{ge}

[Be]{T}

=

00NDITION

,

and

,

171

(Se)

(Sf)

One can rewrite Eqn. (4) for each element as X e= f ~ {T}T[Be]T[De][Be]{ T}dV-f (Q-oc[N e] %--~--,[Ne]{T}dV Ve Ve + f

q[Ne]{T} dS + f

Sle

~ {T}T[Ne]T[Ne]{T}dS

s2e

- f hT [Ne]{T}dS + f h T 2 dS - S o T4[ e] Z}dS $2e $2 e s2e +f

~[Ne]{T}[Ne]{T}[Ne]{T}[Ne]{T}[Ne]{T}dS

(6)

S2 e The functional X can be minimized by equating its first derivative with respect to vector {T} to zero, i.e. ~i ~2 ~n (7)

. . . . .

Differentiating Eqn. (6) with respect to {T} and using the following relationships ~--~-f

~ {T}T [Be]T[De][Be]{T}dV- f [Be]T[De][Be]{T}dV

Ve

f Ve

(8a)

Ve

(Q-pc[N e] --~--)[N ]{T}dV = Y (Q-pc[N el -~-~-)[N ]Tdv Ve (8b)

f Sle f s2e

q[Ne]{T}dS = f

q[Ne]Tds

(8c)

SIe ~h [T]T[Ne]T[Ne] {z}dS = f $2 e

h[Ne]T[Ne]{T}dS

(8d)

172

B.S. Reddy and A.M. Sharan d

T{YT

Vol. 12, No.

I h Too [Ne]{T}dS = f h T® [Ne]Tds e e S2 S2

s

"I

fT

e

(Se)

dS = 0

(8f)

S2 0 C Tt

s

[Ne]{T}dS

= I

e

0 e T4

[ N e ] T dS

(8 g)

e

S2

S2

3

o

~-TfY s e

e

[Ne]{T}[Ne]{T}[Ne]{T}[Ne]{T}[Ne]{T}dS

5

S2 o e [Ne]TIN e]{T}[N e]{T}[N e]{T}[N e]{T}dS

=I

(8h)

e S2 one can rewrite Eqn. (7) in the simplified form as

where the expressions for the matrices [ce], He],

{Fe} and {FR e} are

given in the Appendix B. The vector {F~} is the global force vector due to radiative heat transfer,

and is obtained by summing the elemental force vectors {FRe}. This

elemental force vector contains the terms of higher powers of the nodal temperatures; so does the global one.

Consequently, Eqn. (9) represents the

system of equations which mast be solved using nonlinear numerical techniques.

One of the techniques which can be used is the Gauss-Seldel

iteration technique. Numerical Example and Results In castings,

generally an alloy is enclosed by sand mold or other

materials which have poor thermal conductivity.

With this in mind, the heat

transfer process within two types of materials has been analysed in this work.

One of them is an alloy steel (AISI 4140) and the other one is a

ceramic material. The heat transfer process of the body shown in Fig. i, has been studied by using the finite element method.

It is a ceramic body whose sides are

maintained at 300°C, the bottom surface is insulated and the heat is lost to the surrounding fluid media by convection and radiation processes from the top surface.

Initially the whole body is at 300°C.

This is allowed to cool

"2

VOl. 12, NO. 2

SOLIDSWIR'H P A D I A T I V E ~ O O N D I T I C t q

in an environment at 50°C.

173

This problem has been analyzed as a two-

dimensional problem. The finite elements used were the two-dimenslonal linear tralngular elements.

The details about the element are given in Appendix A.

Using the central finite difference scheme, one can rewrite Eqn. (9)

as ([K G] + ~2 -

-

[CG] ){TG}t+At.(~[cG]_[KG]){TG)t

F5 t

+ 2{FG}t

(FG)t+At

(lO)

In Eqn. (10), the h i g h e r powers of nodal temperatures are p r e s e n t i n {FRG}t and {FR~}t+At . -

The temperatures at time, t , are known, whereas

temperatures at time, t+At, are to be c a l c u l a t e d . temperatures Eqn. (10).

(T}t+At are assumed and s u b s t i t u t e d

At f i r s t ,

the

on the r i g h t hand side of

At t h i s s t a g e , Eqn. (10) becomes 2

[[KG] + ~ - [cG]]{TG}t+At " {A} where, the vector {A} is known. {TG}t+At •

(ii)

Eqn. (12) is solved for the unknown,

If the assumed nodal temperatures and those calculated do not

meet the convergence criteria, then {TG}t+At calculated is used as the assumed temperature vector for the next iteration. Fig. 2 shows the locatlons of nodes on the right half of Fig. i. Because of the symmetry about the y-axls, the study of temperature distribution was done for the right half body only. distribution of various nodes is shown in Fig. 3.

The temperature This figure shows the

continuous decline of temperature for all the nodes with time. time, node 21 is at the lowest temperature.

At any given

This Is because it lles on the

symmetrical axis loosing heat by convection as well as radiation.

It is the

farthest point from the insulated as well as the isothermal surfaces.

There

are oscillatlons in the temperatures initially which decay with time. The temperature distribution of the steel (AISI 4140), a solld of the same geometry as before, is shown in Fig. 4. state much faster than the ceramic body.

This body reaches the steady

The other characteristics of the

temperature distribution are the same as that of the ceramic. The Eqn. (9), which has been derived in this work, is falrly general. The sollds studied in the numerical examples have simple shapes but,

174

B.S. Reddy and A.M. Shmran

Vol. 12, No.

problems having intricate geometry and variable thermal properties can be easily analyzed using this technique.

In the latter case, the properties

should be calculated at each time increment. Conclusions In this work, the finite element analysis has been used to obtain the matrix differential equation for the heat transfer of a body which exchanges heat with its surroundings by the convection and radiation processes.

The

resulting system of equations are solved using the Gauss-Seidel iteration technique.

To illustrate the method, two numerical examples were studied.

The results obtained for different types of materials, such as ceramic and an alloy steel show that this method can be readily used to study various metallurgical processes such as cooling of castings within a mold or heating of ingots for forging or the heat transfer process in a rolling process where heat loss by radiation is quite significant. Nomenclature A

area of the triangular element

[B]

gradient matrix

C

specific heat

[D]

material property matrix

[Cc ]

global capacitance matrix

n

total number of elements

{Fc }

global force vector

{FG }

global force vector due to radiation

h

convection heat transfer coefficient

[Kc l

global conduction matrix

k

thermal conductivty in the x-direction X

thermal conductivity in the y-direction

k Y k

thermal conductivity in the z-direction Z

£x' £y' £z

direction cosines normal to the surface

[Nel

element shape function matrix

Q

heat generated within the body

q

heat flux

TB

boundary temperature

Vol. 12, No. 2

SCLIDSWITH N A D I A T I V E ~ C f l ~ g I T I C N

{T}

elemental nodal temperature vector

{TG}

global nodal temperature vector

t

time

At

time increment

175

density of the material o

Stefan-Boltzmann constant

e

emissivity of the body References

I.

Merle, R. A., Finite Element Analysis of Optimal Heating of a Slab with Temperature Dependent Thermal Conductivity, Int. J. of Heat and Mass Transfer, 22, 4347, (1979).

2.

Sun, R. C., Determination of the Forging ~eatlng Schedule for a Large Hastelloy Alloy X Ingot Metallurgical Transactions, l, 1881, (1970).

3.

Finlayson, P. C., and Schofleld, J. S., J. Iron Steel Institute, 193, 238, (1959).

4.

Segerlind, L. J., Applied Finite Element Analysis, pp. 71-77. Wiley and Sons, Inc., New York (1976).

John

Appendix A C o n s i d e r a t h r e e uode t r i a n g u l a r

e l e m e n t as shown i n F i g . A . 1 .

The

nodes of this element are denoted by i, J, k in the counter-clockwise direction.

The shape functions for th~s linear tralngular element are

1 N 8 = ~-f (a 8 + bBx + csy), 8-i,l,k where a i = xjy k - xkY j aj = XkY i - yiYk

a k = xiY j - x j y i b i = YJ - Yk b j = Yk - Yi

bk "

Yi

-

YJ

cl = Yk - x j

cj

-

%

~ xj

Yl -

~k

xi

176

B.S. Reddy and A.M. Sharan

%7ol. 12, No.

The temperature T is given by

T = [Ni Nj Nk]

Tj

Appendix B The following are the expressions for the elemental matrices: [ce] ffi f

pc[Ne]T[N e] dV

Ve [ke] ffif [Be]T[De][Be]dV + Ve

h [Ne]T[N e] dS

$2 e

{F e} ffif

Q[Ne]Tdv +

Ve {FEe} =

f

f

q[Ne] r dS - f

Sle

f

h T

[Ne]T dS

$2 e

o ¢ [Ne]r [Ne]{T}[Ne]{T}[Ne]{T}[Ne](T}

dS

s2e - f

o

e T 4 [Ne]T dS

e

S2 CONVECTION. RROIRTION

BOUNORR¥

>c c~ z o ~D

iT

.J !-

T

C UJ Z I.-

[

U~

INSULRT[O L,

I-

The s y s t e m c o n f i g u r a t i o n

x OOUNORR¥

0.211

J

"!

FIG. 1 with various boundary conditions

Vol. 12, No. 2

SOLIDS WITH RADIATIVE ~

~ITI(I~

T CONVECT I O N . f l R O I R T | O N

22

21

BOUNDARY

23

24

25

16

gD

"

11

6

10"

2

3 INSULATED

4

5

X

BDUMORflT

FIG. 2 The dlscretlzatlon of the system into finite elements

- - NODE NODE

m

oo

$'.8S

2~. I't

2~.so

T;.((.z.uxEs)

The tlme-temperature

S~), 33

*~.17

ss. o0

FIG. 3 plot of the ceramic body

177

178

B.S. Reddy and A.M. Sharan

if

Vol. 12, N o

O -- NODE I

iili , 2S

ml ,4

w~ L

I

I

;:

w

N

o m

~ 00 The

s:Is

zk.s't TIME

2k.so

sk.ss

4b.,7

ss.oo

FIG. 4 plot of AISI

4140

(MINUTES)

time-temperature

:,Ky K) THE £LE/~NT

Node [ (x L . Yi )

The

details

ode j (xj. yj)

FIG. A. i of the t r i a n g u l a r

element

steel