...

FEA Heat Conduction - Wilson-Nickell

by user

on
Category: Documents
3

views

Report

Comments

Transcript

FEA Heat Conduction - Wilson-Nickell
NUCLEAR ENGINEERING AND DESIGN 4 (1966) 276-286. NGRTH-HGLLAND PUBLISHING CGMP., AMSTERDAM
APPLICATION OF THE FINITE ELEMENT METHOD
TO HEAT CONDUCTION ANALYSIS *
Edward L. W.iLSON and Robert E. NICKELL
Department of Civi~ Engineering, University of California,
Berkeley, California, USA
Received 20 August 1966
A variational principle is applied to the transient heat conduction analysis of complex solids of arbitrary shape with temperature and heat flux boundary conditions. The finite element discretization
technique is used to reduce the continuous spatial solution into a finite number of time-dependent unknowns. The continuum is divided into subregions (elements) in which the temperature field variable is
approximated by Rayleigh-Ritz polynomial expansions in terms of the values of the temperatures at
prescribed boundary points of the subregions. These temperatures at the node points act as generalized
coordinates of the system and, being common to adjacent subregions, enable appropriate continuity requirements to be satisfied over the entire continuum.
The variational principle yields Euler equations of the Lagrangian form which result in the development of an equivalent set of first order, ordinary differential equations in terms of the nodal temperatures. A unique method of numerical solution of these equations is introduced which is stable and re-"
quires a minimum of computer effort. Elements of various shapes and their associated temperature
fields are discussed for one, two and three-dimensional bodies. The method is developed in detail for
two-dimensional bodies which are idealized by systems of triangular elements. The development of a
digital computer program is discussed and several examples are given to illustrate the validity and
practicality of the method.
1. INTRODUCTION
Several a p p r o x i m a t e m e t h o d s of solution to
heat conduction p r o b l e m s , which a r e based on
v a r i a t i o n a l p r i n c i p l e s , have been introduced in
the p a s t s e v e r a l y e a r s [1-4]. The g e n e r a l i t y of
the v a r i a t i o n a l approach with r e s p e c t to a r b i t r a r y b o u n d a r y conditions and m a t e r i a l p r o p e r t y
v a r i a t i o n s h a s allowed f o r the solution of m a n y
c o m p l e x p r o b l e m s . In t h i s investigation, a v a r i a t i o n a l principle is u s e d in conjunction with the
finite e l e m e n t idealization. This r e s u l t s in a
powerful solution technique f o r the d e t e r m i n a t i o n
of the t e m p e r a t u r e d i s t r i b u t i o n within c o m p l e x
bodies of a r b i t r a r y g e o m e t r y .
In the finite e l e m e n t a p p r o x i m a t i o n of s o l i d s ,
the continuous body is r e p l a c e d by a s y s t e m of
e l e m e n t s . In the c a s e of heat conduction an app r o x i m a t e solution within e a c h e l e m e n t f o r the
t e m p e r a t u r e field is a s s u m e d and heat flux equil i b r i u m equations a r e developed at a d i s c r e t e
n u m b e r of points within the finite e l e m e n t s y s t e m . F o r the c a s e of s t e a d y - s t a t e heat flow the
* Accepted by T. A. Jaeger.
technique is completely described by Zienkiewicz [5]. The purpose of this paper is to extend
the technique to the transient heat conduction
problem.
The advantages of the finite element method,
as compared to other numerical approaches, are
numerous. The method is completely general
with respect to geometry and material properties. Complex bodies composed of many different
anisotropic materials are easily represented.
Temperature or heat flux boundary conditions
may be specified at any point within the finite
element system. Mathematically, it can be shown
that the method converges to the exact solution
as the number of elements is increased; therefore, any desired degree of accuracy may be obtained. In addition, for the steady-state condition
the finite element approach generates heat flow
equilibrium equations which produce a symmetric, positive-definite matrix which may be placed
in a band form and solved with a m i n i m u m of
computer storage and time. Also, for the transient problem, a step-by-step procedure is introduced which takes advantage of the special
characteristics of the matrices.
A P P L I C A T I O N O F THE FINITE E L E M E N T METHOD
2. THE VARIATIONAL EXPRESSION
277
H e n c e the u s e of eq. (1) a s ti.~ g t ~ e r a t i n g f u n c tional is v e r i f i e d .
Adopting a s i m i l a r notation to that of Gurtin
[6] the functional:
3. THE FINITE E L E M E N T IDEALIZATION
~t(o) = ½ f{pco . o + o,i*~O . o j - 2 p . / , . o
V
- 2pCOo* 0} (~,t) d l r - f{Qz, n i , o } (£,t) dS
(1)
is introduced o v e r the volume V and the bound a r y S of a ~ontinuum, w h e r e 0 is the t e m p e r a t u r e field as a function of s p a c e and t i m e , kij is
the t h e r m a l conductivity t e n s o r which m a y be
s p a c e dependent, p(x-') is the m a s s d e n s i t y and
c(x--) is the heat c a p a c i t y p e r unit m a s s . The r a t e
at which heat is g e n e r a t e d p e r unit m a s s is d e noted by the function p(g, t). Initial spatial t e m p e r a t u r e d i s t r i b u t i o n is r e p r e s e n t e d b y 0o(X) and
F o r the p u r p o s e of e x p r e s s i n g the g e n e r a t i n g
functional in t e r m s of a finite n u m b e r of unknowns, the body is subdivided into a finite num o
b e r of r e g i o n s . A typical r e p r e s e n t a t i o n of an
a x i s y m m e t r i c body b y a s y s t e m of t r i a n g u l a r r e gions is shown in fig. 1.
AXISOFSYMMETRY
/~6RAI~IT[
T
t
Oi(~, t) f Qi(g, ~') dT.
=
(2)
0
r e p r e s e n t s p r e s c r i b e d heat flux conditions on the
b o u n d a r y S. The convolution s y m b o l i s defined a s
in ref. [6] by the e x p r e s s i o n :
t
{u • v} (~, t) = f u(~, t - T) v(~; T) tiT.
0
(3)
I
I
I
~,/
/
/I
/
. . . . . . . . .
" -
INSULATOR
, , GLAss ~AeR,C
' STEEL S.EU.
_.T
-ASBESTOS
a. Rocket nozzle.
Let 0 be the set of all functions continuous in
the body and s a t i s f y i n g p r e s c r i b e d t e m p e r a t u r e
b o u n d a r y conditions on S. Then, by taking the
f i r s t v a r i a t i o n of eq. (1) and applying the d i v e r gence t h e o r e m , it is s e e n that 5fit = 0 o v e r 0
(0 ~< t < oo) at a p a r t i c u l a r function 0 if and only
if 0 is a solution of the E u l e r equation:
(kij * O,i),j - ~cO + pcO o + p *p = 0
(4)
with the natural b o u n d a r y conditions:
kff • ed - Qi = o.
(5)
b. Finite element idealization.
Through the u s e of the s i m p l e identity:
Fig. 1. Idealization of rocket nozzle by a system of
finite elements.
t
{pc *
B} (E', t)
= f p(Y~) c(~) 0(.~, 7") dT"
0
= p(~) c(x--) [e(~, t) - eo(X--) ] .
(6)
Eq. (~}) m a y be t r a n s f o r m e d into the m o r e usual
f o r m of the t r a n s i e n t heat conduction equation:
(kij * o,i),j - pc • b - p , p = o .
(7)
The solution f o r the t e m p e r a t u r e d i s t r i b u t i o n
within the body is now d e t e r m i n e d by the s t a n d a r d R a y l e i g h - R i t z p r o c e d u r e in which the gene r a l i z e d c o o r d i n a t e s a r e s e l e c t e d a s the t e m p e r a t u r e s at the nodal points of the finite e l e m e n t
idealization. An a d m i s s i b l e t e m p e r a t u r e field i s
one in which the function is continuous throughout
the d i s c r e t i z e d s y s t e m . It h a s b e e n shown that
278
E.L. WILSON and R. E. NICKEL
continuity of d e r i v a t i v e s of the field between e l e m e n t s i s not a r e q u i r e m e n t of the method [5].
The f o r m of the a s s u m e d t e m p e r a t u r e field within a r e g i o n (element) will depend on the specific
type of e l e m e n t used. F o r a typical element m
the a s s u m e d field is written in the following general form:
em(~, t) = aj(t) ~(.~) ,
j = 1, J,
e¢ = ~ {o(t)} T [c] ,{oCt)} + ~ {o(0} T,[K] ,{o(0}
- {e(t)} T [c] ,{o(o)} - {o(,,)} T ,{Q(t)},
(12)
where
M
tc]=
M
tc"],
tn=
[:],
(8)
(13)
M
w h e r e J is the number of e l e m e n t nodes (e.g. J
e q u a l s t h r e e for a t w o - d i m e n s i o n a l t r i a n g u l a r
element). The spatial functions fj must be s e l e c t e d in such a manner that the t e m p e r a t u r e
field is compatible between e l e m e n t s . The aj a r e
r e l a t e d to the t e m p e r a t u r e s at the element nodes
and a r e c a l c u l a t e d by evaluating eq. (8) at t h e s e
nodes. T h e r e f o r e , the t e m p e r a t u r e at any point
within the e l e m e n t may be e x p r e s s e d in t e r m s of
nodal point t e m p e r a t u r e s by the following m a t r i x
equation:
Om(~,t) = (bin(x-)> {o(t)},
w h e r e {0(t)} i s a column v e c t o r of nodal point
t e m p e r a t u r e s for the c o m p l e t e finite e l e m e n t
s y s t e m and the c o e f f i c i e n t s (bin(x-)) r e p r e s e n t the
spatial approximation. M o s t of the coefficients
a r e z e r o s i n c e the t e m p e r a t u r e field within an
e l e m e n t is c o m p l e t e l y defined b y the nodal point
t e m p e r a t u r e s in or adjacent to the element.
Diffcrentiation of eq. (9) with r e s p e c t to the
spatial c o o r d i n a t e s yields a column vector of
{@(t)} = Yl= {~fn(t)}'
and w h e r e [C] is defined as the heat c a p a c i t y
m a t r i x equal to the sum of the e l e m e n t heat c a pacity m a t r i c e s :
[ Cm] = f P m ( ~ Cm(~ (bm(~) T (bin(x-')) dVm •
Vm
(i4a)
[K] is defined as the conductivity m a t r i x and i s
equal to the sum of the e l e m e n t conductivity m a trices:
[K m] = f
and {Q} is defined a s the t h e r m a l f o r c e v e c t o r ,
equal to the s u m of the e l e m e n t t h e r m a l f o r c e
vectors:
{QmCt)} = f
(10)
The g e n e r a t i n g functional, eq. (1), f o r the
finite e l e m e n t s y s t e m m a y be c a s t in the f o r m of
a s u m m a t i o n o v e r all e l e m e n t s :
M
p'm(~) *p(£, t) (bin(x-')) T dV m
Vm
*
temperature gradients which may be written
symbolically as:
{Om,z{g, t)} = [am(,~)] {oCt)}.
[am(x-')]T [km(~] [am(~] dVm , (14b)
Vm
( Q~Cx, t) n i (bin(x-'))T d S m .
(14c)
The f i r s t v a r i a t i o n of eq. (12) y i e l d s :
[c] {0(0} + [K] ,{oCt)} = [c] {0(o)}, {Q(t)),
C15)
a set of linear equations to be solved for the
nodal point t e m p e r a t u r e s of the finite e l e m e n t
r e p r e s e n t a t i o n as a function of t i m e .
~t = ~
/ ½{p,,,CmOm * 0,,, + Ore,i* k':} * O,,,,j
m=l l, m
- 2Pro *Pro *Om " 2PmCmO~n * Ore} (x, t) d V m
4. SURFACE HEAT TRANSFER
- f { O ~ n i , O m } (~,t) dSm,
The r a t e of heat flux a c r o s s a boundary l a y e r
at the s u r f a c e of a body is:
(11)
Sm
w h e r e M ~s the total n u m b e r of e l e m e n t s of the
s y s t e m . It is important to note that the m a t e r i a l
p r o p e r t i e s p, c, and the conductivity t e n s o r k
may be different for each e l e m e n t of the s y s t e m .
The substitution of eqs. (9) and (10) into ( l l )
produces:
q = h , ( O e - O),
where
layer,
face of
ficient
(16)
0e is the known t e m p e r a t u r e outside the
0 is the unknown t e m p e r a t u r e at the s u r the body, and h is the h e a t - t r a n s f e r c o e f f o r the layer. This s u r f a c e h e a t - t r a n s f e r
A P P L I C A T I O N O F THE FINITE E L E M E N T METHOD
effect is considered by subtracting the following
s u r f a c e integral f r o m eq. (1):
f {h *(Oe -½O) *O} (E,t) dS .
(17)
S
For a finite element idealization the surface
t e m p e r a t u r e field for a typical boundary element
m is expressed in t e r m s of the nodal t e m p e r a t u r e by an approximate relationship:
O,n(X, t) = (an(x)) {o(t)}.
(18)
The approximate surface t e m p e r a t u r e field must
be compatible with the approximate volume t e m p e r a t u r e field, eq. (9).
The inclusion of eq. (17) in the variational
expression leads to the following f o r m of the
element conductivity matrix, eq. (14b), and the
element t h e r m a l force matrix, eq. (14c):
[Kn~ = f [am]T[ kin] [am] d Vm
Vm
+ f hm[dm] T [dm]dSm,
Sm
dSm +
Sm
f Oe*hm
{din} dSm •
Sm
(19b)
The added surface integrals are r e q u i r e d only
for elements which a r e subjected to this heatt r a n s f e r effect.
5. TEMPERATURE B O U N D A R Y
CONDITIONS
A t e m p e r a t u r e distribution as a function of
time may be specified along a boundary of the
body by specifying the t e m p e r a t u r e at a d i s c r e t e
number of nodal points along the boundary. This
is accomplished by rewriting eq. (15) in a matrix
partitioned form:
-----4-
a
I Cbb/
*
L0b!
(21a)
where
{~}a} = {Qa} - [Cab] {0b} - [/tab] * {0b}
+ [Cab ] {Ob(O)} .
(21b)
The effective thermal force matrix, {~a}, is
completely defined and the solution of eq. (21a)
for the unknown t e m p e r a t u r e s may be a c c o m plished by standard techniques. H d e s i r e d , the
unknown thermal forces, {Qb}, which a r e a s s o ciated with the specified t e m p e r a t u r e points, a r e
then calculated d i r e c t l y f r o m the second p a r t i tioned equation, eq. (20b).
The nodal point t e m p e r a t u r e s for a finite element s y s t e m , including t e m p e r a t u r e and heat
flux boundary conditiona, must satisfy eq. (21a).
The solution may be obtained, as suggested by
Blot [2], by the standard mode superposition
method in which the initial step involves the
evaluation of the thermal modes and characteristic roots. Since eq. (21a) may contain several
huntired unknowns for a typical finite element
s y s t e m , the mode superposition approach is not
practical because of the tremendous number of
n u m e r i c a l operations r e q u i r e d to solve a c h a r a c t e r i s t i c value problem of this size.
In o r d e r to develop the step-by-step technique
eq. (21a) is rewritten to reflect the solution at
time t in t e r m s of pseudo-initial values at time
t - At:
[c] {o(t)} ÷ [KI * {o(t)} -- [c] {o(t-,,t)} + {o(z)}.
Over the interval (t, t - At) the nodal point t e m p e r a t u r e vector will be approximated by:
,
k%
[Caa] {0a} + [Kaa] * {0a} = {~a} + [Caa] {0a(0)},
6. STEP-BY-STEP SOLUTION PROCEDURE
Vm
f 0mnk (bin)T
a r e the specified t e m p e r a t u r e s , {Qa} a r e the
known t h e r m a l f o r c e s and {Oh} are the unknown
thermal forces.
The f i r s t partitioned equation, eq. (20a), may
be rewritten:
(19a)
{Qm(t)} = f Pro*Pro (bm~T dVm
+
279
(20)
Cba '1CbbJ ~ ) J
w h e r e {Oa} are the unknown t e m p e r a t u r e s , {0b}
{0(~)} = {F1} ~ [g - (t - At)] { F 2 } .
(23)
The calculation of the thermal force vector r e quires the time integration of the heat flux and
convective boundary conditions as well as the
280
E . L . WILSON and R. E. NICKEL
volumetric heat generated in the body. These
functional values are assumed to vary linearly
over the s a m e time interval (t,t - At).
Evaluating eq. (23) at the end points of the interval, thus insuring continuity of the temperature field with respect to the time variable, the
expression for the nodal point temI~erature vector becomes:
{0(~)} = (t~__~){O}t.A t + [/~ -(t-At)]
At
{o}t.
(24)
Substituting these expressions into eq. (22) r e sults in the following matrix equv~tion for the
finite element system:
([C]
+ ½At
[K]) {O}t = ([C]- ½AI [K]), {Ott_A t
+ ½At {Q}t + ~At {Q}t-At"
(25)
This equation is written in a more convenient
Table 1
Summary of step-by-step solution method
Initial
calculations
1. F o r m [C], [K], and {Q}
form for the purposes of the solution technique
as"
[Xl {oi} -- {oi},
(26a)
where
1
0
t
{Of} = 2 { }t + "~{O}t-At
2
(26b)
[l~] = [K] + ~7 [ C ] ,
(26c)
{o}t.At.
{Qi} = ½{e} + ½IQ}t-At + 2[c]
(,-6d)
7. TWO-DIMENSIONAL TRIANGULAR
ELEMENTS
As an example of the application of the method
to a specific geometry, the development of the
element heat capacity matrix (eq. (14a)), the element conductivity matrix (eq. (14b)) and the element thermal force vector (eq. (14c)) is given
for the plane triangular element (fig. 2). The development of these m a t r i c e s for other types of
elements is analogous. The temperature d i s t r i bution within a typical triangular element m is
assumed to be of the form:
Eqs. (14)
Om(x, y, t) = a 1 (t) + xa2(t ) + ya3(t).
2. Modify for temperature boundary conditions
3. F o r m [K']
Eq. (26c)
2
[~l -- [K] + ~
[c]
(27)
The constants are expressed in t e r m s of the
temperature at the nodal points i, j, and k by:
4. Triangularize [K--]
For
each
time
increment
1. Calculate {Qi}
a2
Eq. (26d)
1
{Oi}
,
(28a)
¢3
2
{qi} = ½{q}t + i {O}t-At + ~ [el {o}t.A t
2. Define
-- [~1 '
OkJ
where
Eq. (26b)
{o~} = ½{o}t + ½{o}t.,,t
3. Evaluate {Oi} by solving
•/33i"
~3j
~3k
ix] {o~} -- {e~}
4. Calculate {O}t
= r4
{o}t-- 2{o~} - {o}t.At
t
Yj'Yk
"
Yk-Yi
Yi-Yj
•
xk-xj
•
xi'x k
xj-x i
•
and
(28b)
A = ½[xfly k -y~) + xi(y j -yk) + xk(y~-~'j)] •
(28c)
5. Repeat for next time increment
APPLICATION OF THE FINITE ELEMENT METHOD
281
[k,.] = [kxx kxy]
xk i
(3])
F r o m eq. (14b) the e l e m e n t conductivity m a t r i x
is defined as:
X j,, ,
xi
[gml -- f [aml T Ik"l [aml d V.,
Vm
Yi
Yk
Yj
Q
.. K m. • . K ~ .
Fig. 2. The t r i a n g u l a r element.
m
• . Kih...
•
The t e m p e r a t u r e can then be written as:
,
(32)
+ ~3rkyy~3s) ,
(33)
• .l~jj...K;~...
•
:
o
1
oi (t) l
•
. . Kh~..
. . . Kt~
. z.
I
Orn(X,y, t) = (brn(X, y))'
(29a)
-" . •
l(hh,
where the typical t e r m :
i:
K r s = A m t m (f32rkxxf~2s +/32r kxyt33s
i %(t)J
where
+ •3rkxyfJ2s
(brn(X, y)> = ( . . . bi(x, y). . . bj(x, y ) . . . bk(x, y))
(29b)
and
in which Am is the a r e a of the element given by
eq. (28c) and t m is the thickness of the element.
The heat capacity m a t r i x for the triangular
element is obtained by evaluating eq. (14a):
b i ( x , Y ) = ~ l i + f32ix + ~3iY ,
b j ( x , y ) = # l j + f32jx + ~3jY ,
(29c)
[C'nl - f Pm Cm [b,n] r [bin] d V m
v,.
bk(x, Y) =/31k + 132kx + (33kY •
The t e m p e r a t u r e gTadients, eq. (10), a r e
evaluated for a triangular element by differentiating eq. (29a):
.
= ...cjm...Cj~...Cj~...,
(34)
~a
...Cki...C
{
O,x 1
O,yj=[am]
,
...Ckk...
(30a)
where the typical t e r m is:
C~s - f pm C,n br bs d V m .
vm
where
Cam]:
~3i
~3j.
~3k.
The t h e r m a l conductivity for two-dimensional
flow is:
The integration over a triangle of the
of the two linear functions, br and b s,
written as a m a t r i x equation in t e r m s
value of the functions at the v e r t i c e s of
angle.
(35)
product
can be
of the
the t r i -
E. L. "4dILSON and R. E. NICKEL
282
lil
=
12
-.--.e, x
a. L i n e a r
--L~r
1
b. Parabolic
c. Cubic
I. ONE-OIMEI'IISIONAL ELEMENTS
(36a)
where the superscripts indicate the nodal point
where the function is evaluated. This e~:pression
simplifies considerabty for the plane ,--lement,
becoming:
,
.
,
.
,
. .
•
,
*
*
.
.
a. Reclan~lle
b. T r , o n g l e
2. TWO-DIMENSIONAL
•
c. S~x-Pomt Triangle
ELEMENTS
...2...1...1...
. ° ,
ECru] P m C m t m A m
=
12
.
. .
o
1.
.
.
.
.
. . , ,
.
.
,
.
2.
.
.
.
.
.
, . ,
.
.
1.
.
.
, . ,
. . . .
...1...1...2...
. , .
,
.
.
.
,
.
,
.
.
,
.
,
a. Rectongulor Rmq
For heat generation within the element the t h e r mal force matrix, eq. (14c), is:
.3. AXISYMMETRIG ELEMENTS
o. Prism
Qk
b
c. T a n - P a i n t Tetrahedron
Tetrahedron
Fig. 3. Summary of finite elements used in the representation of solids.
I
where the typical t e r m Qm is evaluated as:
3
c Six-Point Triangular
RInl
(37)
[Qml = f Pro*Pro [ " m ] T d V m =
=
b Trlanqular Ring
(38)
Equations which include surface heat t r a n s f e r ,
eqs. (25a) and (25b), will not be given here since
their development is straight-forward.
posed elements have been used in the s t r e s s
analysis of solids [7] and their extension to the
heat t r a n s f e r problem involves the introduction
of the scalar t e m p e r a t u r e field in place of the
vector displacement field.
For one-dimensional solids, three are proposed and the t e m p e r a t u r e field associated with
each is as follows:
Two nodal point element:
8. ADDITIONAL ELEMENTS
0 (x, t) = a 1 (t) + xa2(t).
(39a)
Three nodal point element:
Solids of almost any geometric shape may be
represented by a combination of different finite
elements. The only r e q u i r e m e n t is that the a s sumed temperature field is continuous throughout the solid. In addit,on to the plane triangular
element, which is discussed in detail in the p r e vrous section, several different elements for the
representation of various types of solids a r e
summarized in fig. 3. The use of these elements
in heat t r a n s f e r problems has not been c o m pletely investigated; however, most of the p r o -
o(x,t)
=at(t)+xa2(t)
+
x2a3(t).
(39b)
Four nodal point element:
0 (x, t) = a l ( t ) + m2(t ) + x2a3(t) + x3a4(t).
(39c)
It is apparent that one-dimensional elements
with any number of nodal points are admissible.
In addition to the plane t r i a n g u l a r element,
several other two-dimensional e l e m e n t s are
possible. A compatible t e m p e r a t u r e field for the
~-ectangular element is:
APPLICATION OF THE FINITE ELEMENT METHOD
O(x, y, t) = a l ( t ) + Xas(t ) + y a s ( t ) + xYa4(t ) .
283
(40a)
The p r a c t i c a l u s e of the r e c t a n g u l a r e l e m e n t i s
l i m i t e d since it cannot be u s e d in a p p r o x i m a t i n g
c u r v e d boundaries. Another possible t w o - d i m e n sional e l e m e n t is the six-point triangle which is
a s s o c i a t e d with the following t e m p e r a t u r e field:
O(x, y, t) = a l ( t ) + XOts(t) + y a s ( t ) + zSa4(t )
+ xyas(t) + y 2 a 6 ( t ) .
(40b)
Complete t e m p e r a t u r e continuity is m a i n t a i n e d
by the introduction of additional unknown t e m p e r a t u r e s at the m i d - p o i n t of e a c h side.
Ring e l e m e n t s with the s a m e c r o s s section a s
the t w o - d i m e n s i o n a l plane e l e m e n t s m a y be u s e d
in the r e p r e s e n t a t i o n of a x i s y m m e t r i c solids and
f o r the c a s e of a x i s y m m e t r i c t e m p e r a t u r e d i s tribution, only the evaluation of the volume int e g r a l s is different. In the c a s e of non a x i s y m m e t r i c t e m p e r a t u r e distribution the p r o b l e m m a y
be expanded in a F o u r i e r s e r i e s in the c i r c u m f e r e n t i a l d i r e c t i o n ~. F o r the t h r e e - p o i n t t r i angle the a s s u m e d t h r e e - d i m e n s i o n a l t e m p e r a t u r e field is of the following f o r m :
O(x, y,~>, t) = ~ [ ( a ~ l ( t ) + ~ l ( t )
n
+ Ya31(t)] sin n~
Fig. 4. The quadrilateral element.
r i l a t e r a l m a y be e x p r e s s e d in t e r m s of the t e m p e r a t u r e s at the o t h e r f o u r nodes and e l i m i n a t e d
as an unknown in the final set of linear equations.
9. EXAMPLES
To i l l u s t r a t e the solution technique on a t r a n sient heat conduction p r o b l e m for which the exact
solution is known, the o n e - d i m e n s i o n a l e x a m p l e
of constant heat flux applied to a s e m i - i n f i n i t e
solid is c h o s e n [8]. If the origin is chosen at the
s u r f a c e of the solid, the m a t e r i a l p r o p e r t i e s
c h o s e n to be unity, and the applied heat flux p e r
unit t i m e also chosen to be unity, then the solution is:
0 ix, t) = 2
+ ~[af2(tl + x ~ 2 ( t ) + y~2(t)] cos n ~ .
z exp
) - ~x e r f c
.
(411
A d m i s s i b l e t e m p e r a t u r e fields for : h r e e d i m e n sional e l e m e n t s a r e a s follows:
R e c t a n g u l a r right p r i s m
0 ( x , y , z, t) = a l ( t ) . a2(t)x + a3(t)y + ~4(t)z
+ ~5(t)xy + a6(t)xz + ct7(t)yz + a 8 ( t ) x y z •
Tetrahedron
O(x,y, z, t) = a l ( t ) + a2(t)x + a3(tly + a 4 ( t ) z .
Ten-point t e t r a h e d r o n
0(x,y, z, t) = ~l(t) + a2(t)x + a3(t)y + a4(t)z
+ a5 (t)x2 + a6 (t)y2 + a7 (t)z2
+ a8(t)xy + ag(t)xz + a l o ( t ) y z •
In addition to the e l e m e n t s given in fig. 3,
e l e m e n t s m a y be developed by combining t h e s e
basic e l e m e n t s . F o r example, a t w o - d i m e n sional q u a d r i l a t e r a l e l e m e n t m a y be f o r m e d
f r o m f o u r t r i a n g u l a r e l e m e n t s as shown in fig. 4.
The unknown t e m p e r a t u r e at node 5 of the quad-
The spatial t e m p e r a t u r e distributions for v a r i o u s t i m e s a r e plotted in fig. 5. Except at e a r l y
t i m e s a g r e e m e n t with the exact solution is excellent. A c o a r s e m e s h solution at time t = 1.0 is
given in fig. 6. Also, fig. 6 i l l u s t r a t e s the e f f e c t s
of diagonalizing the heat capacity m a t r i x by adding the coefficients of each row and placing the
sum on the diagonal. This lumped heat c a p a c i t y
m a t r i x r e d u c e s by a l m o s t 50~ the n u m b e r of
n u m e r i c a l operations in the s t e p - b y - s t e p s o l u tion; however, the l o s s of a c c u r a c y a p p e a r s to
be s m a l l due to this approximation.
The t e m p e r a t u r e at x = 0 v e r s u s t i m e Js plotted in fig. 7 for two d i f f e r e n t t i m e i n c r e m e n t s .
The approximate solutions oscillate about the
exact solution and c o n v e r g e to the t r u e solution
at l a r g e t i m e s . This type of behavior is typical
even for v e r y l a r g e t i m e i n c r e m e n t s .
A specific application of i n t e r e s t which i n c o r p o r a t e s the convection boundary condition and
the r e d u c t i o n of the t r a n s i e n t solution s c h e m e to
a s t e a d y - s t a t e p r o b l e m is solved by P e a v y [9]
using F o u r i e r a n a l y s i s . An exposed e x t e r i o r
r e c t a n g u l a r c o n c r e t e c o l u m n (see fig. 81 is sub-.
j e c t e d to outside w e a t h e r conditions and an in-
E. L. WILSON and R. E. NICKEL
284
\
o_
:o:.%.* SO-,T,O., .oos,
I
O(~
05,
O£PTH
IO
15
20
Fig. 5. T e m p e r a t u r e distribution at various times v e r s u s depth indicating spatial a c c u r a c y of solution technique.
I¢
0
DISTRIBUTED HEAT CAPACITY
0
LUMPED HEAT CAPACITY
W
EXACT SOLUTION
i-
0.5 m
0
I
10
2,0
OEPTH
3,0
I
I
40
5,0
Fig. 6. T e m p e r a t u r e distribution versus depth at time t = 1.0 indicating effect of lumped v e r s u s distributed heat
capacity for a crude {four element) spatial approximation.
t e r i o r environment s e p a r a t e d by an abutting wall.
The ambient t e m p e r a t u r e and the s u r f a c e heat
t r a n s f e r coefficient on the i n t e r i o r face ~ = 0 a r e
100°F and 0.5 Btu/h ft 2 OF, r e s p e c t i v e l y , and
the c o r r e s p o n d i n g quantities on the e x t e r i o r f a c e
y = l a r e 0OF and 6.0 B t u / h ft 2 OF. The v a r i a t i o n
of the ambient t e m p e r a t u r e and s u r f a c e heat
t r a n s f e r coefficient on the f a c e s x = 0 and x=2a
is indicated in fig. 8. The t h e r m a l conductivity
was taken to be 1.0 B t u / h ft OF.
T h r e e p r o b l e m s , r e p r e s e n t i n g a column having l = 36 in., a = 7 in., and having t h r e e d i f f e r ent p o s i t i o n s f o r the abutting wall, w e r e solved
using a s y s t e m of 296 nodal p o i n t s and 252 quad-
APPLICATIONOFTHEFINITEELEMENTMETHOD
o
0.6 --
285
~
sS
--
E ELEMENTSOLUTION(At • 0 IO1
o, -
/
o~
rj
/
~---Ex~
~,J,00N
FINITE ELEMENT SOLUTION ( A ! , O.OS)
t
"2
/1'I I .."
t
es
oo
I
I
OI
I
I
I
I
O.2
TIME
O3
I
l
I
0.4
I
0.5
Fig. 7. Temperature distribution v e r s u s time at x= 0 indicating stability of solution technique.
IOOi
H:608TU/hr fl 2 °F
ee : IO0* F
IH:O,5 STU/hr fl 2 °F
/
J/
9O
8°: 0 e F
8O
I.y
f - - L I N E OF SYMMETRY
x:2o
7O
6O
y:l
",,~\
,=o~,\\
_
',,,\ ~ . . - - - , : 0
-
'~,
~\ \
-
SOI.N
'~,
•,~\
,~\
4O
Fig. 8. T r a n s v e r s e section of rectangular c o n crete column with abutting wall separating the
inside and outside environment (outside face l o cated at y ffi /).
30 ¸
20
~
/
/
/
~
PEAV-T SOLN.'-"~-_/,~C//I\',
I0
Fig. 9. Steady state temperatures over the width of the
columns against distance measured from the inside
surface of the column (positions of abutting walls ~re
shown).
\
IE
24
OISTANCEFROM INSIDEOF COLUMN INCHES
N
286
E.L. WILSON and R. E. NICKEL
rllateral elements
o v e r the half-width of the
c r o s s s e c t i o n . The r e s u l t i n g t e m p e r a t u r e d i s t r i b u t i o n s at x = 0 and x = a as a function of the d i s t a n c e f r o m the i n t e r i o r face a r e shown in fig. 9.
The r e s u l t s of P e a r y , in t e r m s of a mean t e m p e r a t u r e o v e r the width of the c r o s s section, a r e
shown on the s a m e figure. The total IBM 7094
t i m e f o r all t h r e e solutions was l e s s than one
minute.
10. DISCUSSION
The finite e l e m e n t method has been applied to
the t r a n s i e n t heat conduction a n a l y s i s of c o m p l e x
solids. The method p o s s e s s e s unique advantages
as c o m p a r e d to other n u m e r i c a l approaches with
r e s p e c t to t r e a t i n g v a r i a b l e d i s t r i b u t i o n of t h e r mal ~ r o p e r t i e s , t e m p e r a t u r e and heat flux bound a r y conditions and solids of a r b i t r a r y g e o m e t ric shape. The new s t e p - b y - s t e p solution t e c h nique is stable and p r o v i d e s an efficient digital
c o m p u t e r approach for a l a r g e c l a s s of t i m e dependent p r o b l e m s . An additional advantage,
which has not been d e m o n s t r a t e d in this p a p e r ,
is in the solution of t h e r m a l s t r e s s p r o b l e m s .
The method of t e m p e r a t u r e a n a l y s i s p r e s e n t e d
h e r e is c o m p l e t e l y c o m p a t i b l e with the finite
e l e m e n t solution for the s t r e s s distribution [7,
10,1~]; hence, both p r o b l e m s m a y be solved
s i m u l t a n e o u s l y within the s a m e c o m p u t e r p r o gram.
RE I~ERENCES
[1] M.A. Blot, Thermo-elasticity and irreversible
thermodynamics, J. Appl. Phys. 27 (1956) 240.
[2] M. A. Blot, New methods in heat flux analysis with
application to flight structures, J. Aeron. Sci. 24
(1957) 857.
[3] M.A, Blot, Further developments of new methods
in heat-flow analysis, J. Aero/Space Sci. 26 (1959)
367.
[4] T. J. Lardner, Blot's variational principle in beat
conduction, AIAA J. 1 (1963) 196.
[5] O.C. Zienkiewicz and Y.K. Cheung, Finite elements
in the solution of field problems, The Engineer
220 (1965) 507.
[6] M. E. Gurtin, Variational principles for linear initial-value problems, Quart. Appl. Math. 22 (1964)
252.
[7] E. L. Wilson, Structural analysis of axisymmetric
solids, AIAA J. 3 (1965) 2269.
[8] H. S. Carslaw and J. C. Jaeger, Conduction of heat
in solids, 2nd ed, (Oxford University Press, 1959).
[9] B.A. Peavy, Steady-state heat conduction in an
exposed exterior column ofrectangnlar cross section, J. Res. Natl. Bur. Stand. 69c (1965) 145.
[10] Y. R. Rashid, Analysis of axisymmetric composite
structures by the finite element method, Nucl.
Eng. Design 3 (1966) 163.
[11] R.W. Clough and Y. R. Rashid, Finite element analysis of axisymmetric solids, Proc. ASCE, Eng.
Mech. Die. 91 (1965) 71.
Fly UP