SEC. 3 .1
C HAPTER
3
A PPROXIMATE M ETHODS
FOR DIRECT H EAT
C ONDUCTION P ROBLEMS
3.1
I NTRODUCTION
In this chapter some approximate methods are given for solving direct transient
heat conduction problems. These direct techniques are the first stage of solution
procedures for the inverse heat conduction problems. The second stage, which
involves specific algorithms for the I HCP, is developed in Chapter 4. The
approximate heat conduction solution procedures o f this chapter are combined
in Chapters 5 and 6 with the I HCP algorithms of Chapter 4.
3.1.1
78
heat conduction problems are those using Duhamel's theorem and Green's
functions both of which yield convolutiontype integrals. Since these methods
are very similar, only Duhamel's theorem is given. A numerical approximation
of the integral is also given. The heatconducting body can be of arbitrary shape
and of nonhomogeneous material. Although the heat transfer may be two or
threedimensional, in this chapter there is a single forcing input which is only
time dependent.
A new technique which utilizes Duhamel's theorem for solving connected
basic geometries is called the unsteady surface element method.!  ~ I t also
involves solving an inverse problem but there are no measurement errors and
the matching conditions are a t surfaces rather than at interior points o f the
bodies.
The solution of a nonlinear transient heat conduction problem requires an
approach that discretizes the partial differential equation. Two methods for
solving the nonlinear transient heat conduction equation are the finite difference
and finite element (FE) methods. An alternative method is used in this chapter
in which conservation o f energy is applied directly to finite control volumes;
it is called the finite control volume procedure.
A nonlinear problem is one for which the heat conduction equation o r a
boundary condition is a nonlinear function o f temperature. O ne example of a
nonlinear partial differential equation is:
_~
ax
(k aT)=pc aatT
ax
(3.1.1 )
where k a nd/or pc are functions of temperature T. I f k and pc are functions of
position, x, only, then Eq. (3 .1.1) is linear. An example o f a nonlinear boundary
condition is the radiation condition,
a
 k  TI
=u£[T!,T4(O, t)]
a x .. = 0
V arious N umerical A pproaches
Partial differential equations including the transient heat conduction equation
can be solved using a variety of methods including exact and numerical procedures. The exact methods include the classical methods of separation of
variables and Laplace transforms.
Two types o f numerical procedures are discussed in this chapter. O ne is
based on an integral formulation of the mathematical model and the other on a
differential form o f the model.
The transient heat conduction equation can be either linear o r nonlinear.
F or the linear case the partial differential equation formulation can be equivalently presented by an integral equation. The advantage of the integral form
is that the space dependence can be represented exactly leaving only the time
dependence to be approximated.
Two different approaches that employ integral equations for solving transient
79
I NTRODUCTION
3.1.2
(3.1.2)
S cope o f C hapter
This chapter gives a derivation of Duhamel's integral (Section 3.2), and a
numerical approximation thereof. Section 3.3 discusses finite control volume
solutions of the transient heat conduction equation. Both linear and nonlinear
cases are included.
Section 3.3.5 demonstrates that the linear finite control volume approach
can yield expressions that are analogous to those obtained from the Duhamel's
theorem method. This is important because the same linear I HCP algorithms
can be employed using many different methods of approximating the heat
conduction equation and boundary conditions; such methods include a
numerical form o f Duhamel's theorem, finite differences, and finite elements.
CHAP.3
80
3 .2
3.2.1
M ETHOOS FOR DIRECT HEAT C ONDUCTION PROBLEMS
D UHAMEL'S T HEOREM
SEC. 3 .2
D UHAMEL'S T HEOREM
81
are used to represent the heat fluxes between times of, respectively,
OtOAt, AI t OAz,···, AMI to 'M
D erivation o f D uhamel's T heorem
The corresponding heat fluxes are denoted q l' qz, . .. , qM o r in general
Duhamel's theorem can be considered to be a result of the principle of superposition and thus it is only valid for linear cases. There are several ways of
deriving it, one of which uses Laplace transforms. See Ozisik (Reference 6, p. 197)
and Luikov (Reference 7, p. 344). Another derivation uses the concept o f superposition; that is the method used in Arpaci (Reference 8, p. 307) and Myers
(Reference 9, p. 155). T he following derivation employs the principle of superposition and it also produces a convenient numerical approximation.
Duhamel's theorem employs a "building block" solution which is used with
the suPerposition principle to obtain the temperature a t any point and time.
One such solution is </J(r, t) which is for the temperature rise at a point r in a
heatconducting body due to the heat flux,
O,
t<O
q(t)= { 1, t >O
(3.2.1)
This is sometimes called a unit step heat flux o r even a constant unit heat flux.
The thermal properties of the body are independent of temperature but can be
functions of position, and the temperature distribution need not be onedimensional. F or the derivation herein the only requirement is t hat the prescribed surface heat flux is given by the product,
q(s, t) = f(t)g(s)
where s is a coordinate measured along the surface. For convenience, the heat
flux is written as a function o f time only.
The surface heat flux is approximated in the manner shown in Figure 3.1.
The heat fluxes a t times
(3.2.2)
The initial temperature distribution in the body is the constant value of To.
Then using the principle o f superposition, the temperature a t location r, at
time t M , is composed o f contributions due to the heat flux components ql
through qM inclusive:
T(r, t M)= TO +ql[cp(r, tM Ao)</J(r, tM  AI)]
+qz[</J(r, t MAdcp(r, t MAZ)]
(3.2.3)
+qM[</J(r, tM  AM d</J(r, tM  AM)]
where </J(r, tM  AM) = </J(r, 0 )=0.
Using equal time steps and times o f
Aj=i.:1A, ; =0,1, . .. , M
(3.2.4)
A j Ai = j.:1A i.:1A=(j  ;).:1A=Aji
(3 .2.5)
yields:
Then Eq. (3.2.3) can be written as
(3.2.6)
For the limiting case o f .:1A+O a nd for a fixed time tM = t, Eq. (3.2.6) becomes
the integral
(3.2.7)
From the relation
q (t)
a</J(r, t  A) a</J(r, 1  A)
aA
at
(3.2.8)
Duhamel's theorem (sometimes called Duhamel's integral) becomes
T(r, 1)=To+ f~ q(A) a</J(r~:A) dA
~~...4:'':'_~£.1
1M
o
Al
A2
A M_1
FIGURE 3.1 Approximate representation
o f q (t) by M steps in heat Hux .
(3.2.9)
Equation (3.2.9) is a heat flux form o f Duhamel's theorem; it is a convolution
because there is a product o f two functions, one of A and the other o f t A;
that is, there is folding o r convoluting o f one function with respect to the other.
For a temperaturebased form of Duhamel's theorem, see Eq. (2.4.1) a nd
References 6  10.
C HAP.3
82
3 .2.2
M ETHODS FOR DIRECT H EAT C ONDUCTION PROBLEMS
SEC. 3.2
D UHAMEL'S THEOREM
Then the T+ values at , ; = 0.1, 0.2, a nd 0.3 a re used to obtain the ~tPl values,
N umerical A pproximation o f D uhamel's T heorem
x
0.0\
+
~l/Io=l/I, = k 1/11 = 40 (0.003943)=9.8575E7 Cm1/W
A numerical approximation for Eq. (3.2.9) is o btained from Eq. (3.2.6). Note that
<p(r, tM  A.._ d  <p(r, tM  A..) = <p(r, t M.+ d  <p(r, t M.)=d<p(r, t M.)
83
X
(3.2.10)
~I/I, =1/121/1, = k (tP; 1/Ii>=6.69725E6 Cm 2/W
a nd thus Eq. (3.2.6) c an be written as
~1/I2 =tPl1/I2 = i (tPj  tP;) = 1.029025E5Cm z/W
M
L
T (r, t M )= T o+
q.d<p{r, t M .)
(3.2.11)
n= 1
T he q values are evaluated at 1= 0.5, 1.5, a nd 2.5 s so that
T he notation can be simplified by omitting the r dependence and writing
q l=0.5E6, qz=1.5E6, q l=2 .5E6W/m z
and thus the temperatures found from Eqs. (a), (b), a nd (c) a re
M
T M=To+
L
q .d<PM.,
d<Pi=<Pi+l<P,
(3.2.12)
TI = 30 + 0.5E6(9 .8575E7) = 3O.493°C
n =1
where q. is best evaluated a t time ( nt)dt as indicated by Eq. (3.2.2). I f t he
actual heat flux is constant over each time step, Eq. (3 .2.12) gives an exact
expression for the temperature, TM ; o n the other hand, Eq. (3.2.12) yields
approximate results i fthe t rue heat flux varies over the time steps. This equation
is quite important in investigation o f the I HCP because it gives a convenient
expression for the temperature in terms o f t he heat flux components. The
observation that the sum o f t he subscripts on the right side o f Eq. (3.2.12) is
M can aid in memorizing this equation.
Tz = 30 +0.5E6(6.69725E6) + 1.5E6(9.8575E7)
= 34.827°C
Tl = To + 0.5E6( 1.029025E5) + 1.5E6(6.69725E6)
+ 2.5E6(9.8575E7) = 47.655°C
T he exact temperatures are obtained from Carslaw and Jaeger (Reference 10, p. 77),
T = To+(:~)(4a1)3/2 i l erfc[(41;)I/
Z
]
(3.2.13)
where
EXAMPLE 3.1. Calculate the temperature at 1 cm inside a thick steel plate ( k=
40 W/mC, a = 1 0 5 m1/s) that is initially a t 30°C with an applied heat flux o f
(3.2.14)
The temperatures T" T2 , a nd Tl a re obtained from Eq. (3.2.13) at times 1 = 1 ,1=2, a nd
1= 3 s o r equivalently for I ; = 0.1, 0.2, a nd 0.3. T he resulting exact values are
q =4 ot, 4 o=106 W/m 2s, t ins
starting a t time zero. Find values of the temperature a t times 1,2, and 3 s. Use Eq.
(3.2.12) to approximate each temperature and compare it with the exact value.
Solution. T he temperatures arc denoted TJ> T2 , and Tl a nd correspond t o 1 = 1, 2,
a nd 3 s; from Eq. (3 .2.12) they are given by
TI = To+ql~</>o
(a)
Tl = To+ql~</>I +q2~</>O
(b)
Tl=To+q'~</>2+q2~</>' +ql~tPO
(c)
Hence the errors in the approximate temperatures are + 0.305,0.753, a nd 0.930, respectively. These errors are increasing, but relative to the exact temperature rises they are
0
decreasing; namely, 163%, 18.5%, and 5.6%, respectively.
T he tPi values are for a semiinfinite body with a unit step increase in surface heat flux;
tPi values can be found using entries in Table 1.2 for T+ defined by Eq. (1.6.18c) and
letting q c= 1, T(x, t)=</>, and To=O; hence
t P=(:) T+
=e~l) T+
Also the time 1 is related to , ; by Eq. (1.6.18d) which gives
10 5
1 = ==0.11
2
•x
10 4
+
al
3.2.3
M atrix F orm o f D uhamel's T heorem
I t is frequently advantageous to perform algebraic manipulations utilizing a
matrix form o f t he model. This subsection displays a matrix form for the
numerical statement o f D uhamel's theorem.
An expansion o f t he constant heat flux numerical approximation o f
D uhamel's integral, Eq. (3.2.12), with M replaced by 1 to M + r1, is
1'. =
To +q1d<po
T l = TO+q1d<Pl +q1d<po
84
C HAP.3
M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS
T",+ 1 = To+qIAtP",+q2AtP"'1 + . .. +q",AtPl +q",+ ,AtPo
T2
TJ
AtPo
AtPl
AtP2
AtPo
AtP,
AtP",_,
AtP",
AtP"'2
AtP",_,
AtP"'J
AtP"'l
T"'+'_I
AtP"'+,2 AtP"'+,J AtP",+,4
q",
(3.2.16)
q"'+1
where 1 is a vector o f ones. I n a more compact form Eq. (3.2.16) is
(3.2.17)
T =Xq+Tol
where T and q are the appropriate vectors displayed in Eq. (3.2.16) a nd X is
the lower triangular matrix.
X=
AtPo
AtP"'1
AtP",
AtP"'2
AtP",_,
AtP", 3
AtP"'2
(3.2.19)
1
X=[~::
AtP,1
[q",]
q = ~"'+I
AtPo
AtP,2
.0.
_
1 1 ;Oq
. .J
AtPo
AtP,
AtPo
AtP'1
AtPr2" . AtPo
(3.2.18a)
The X matrix is called the pulse sensitivity coefficient matrix for q. Its structure
is lower triangular with AtPo's along the main diagonal. AtP,'s along the diagonal
t.",+d''';''''I;o
[
(3 .2.20a, b)
qU+,1
t"'I''';O
AtP",+,2 AtP"'+,  J AtP",+,4
(3.2.18b)
This is called the s tandard;{orm for the linear I HCP . This equation is n ot
restricted to representing a numerical form o f Duhamel's Theorem. I t can be
derived using Taylor series for the linear inverse heat conduction problem.
The various components o f Eq. (3.2.19) are
T'"
T = : "'+1
•
[
1''''+'1
q"'+,I
AtPo
AtP,
i <j
T =Xq+1'lq;o
+ Tol
AtPo
AtP,
AtPl
O.
oqi
f e yv1 iJ <2 rO+GI .,.e.
ql
ql
q3
X
. I}
which shows that X i) depends only o n the difference between i a ndj for linear
problems.
In the sequential I HCP methods investigated in C hapter 4. the heat flux
components. q l. q2 • ... , q",  I. are considered previously estimated and are
denoted iiI>' ..• ii",  " Estimates of q"" q"'+I • ...• q"'+, _1 are needed. Moreover. there may be a known heat flux at x = L , a known convection condition a t
x =L. o r known internal volumetric heating. F or these cases it is convenient
to use the modification of Eq. (3.2.17).
AtPo
T",
T",+,
x .. =01i={AtPii. i~j
(3.2.15)
which can be written in expanded matrix form as
'Ii
85
D UHAMEL'S THEOREM
just below the main one. and so on. The components of the X matrix are called
sensitivity coefficients because they are equal to the first derivatives of the
temperature with respect t o h eat flux components. From Eq. (3.2.15). one can
verify that the ijth component o f X is
T",=To+qIAtP"'1 +q2 AtP"'  2+ '" +q"'  IAtPl +q",AtPo
T",+,_, = To+q IAtP",+,2 + ... +Q",+,2AtPl +q",+,_,AtPo
S EC.3.2
(3.2.21)
]
(3.2.22)
f",+r  d,.,;,.,.I ; '" ; ,., •• 1 = 0
Equations (3.2.19)  (3.2.22) apply for the numerical convolution as well as for
the finite difference methods. F or the case of known heating o r energy sources
other than where q is to be estimated. X and Tlq ; o must be properly interpreted.
The quantity tPi needed in X is the temperature rise at the sensor location
for a unit step change in the surface heat flux a t t = 0 for the same partial differential equation and boundary conditions as the original problem except the
differential equation and boundary conditions (other than at x=O) are homo
86
C HAP.3
M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS
geneous. An example is the problem for a hollow cylinder,
~~
r ar
(r
O T)=PC a T _ g(r, t ),
or
at
k
aTI
ar
k 
=q(t),
a <r<b
unknown heat flux
(3 .2.23)
(3.2.24a)
r=G
aTI
kar
=h[T(b,t)T,., (t)]
(3.2.24b)
r =b
(3 .2.25)
T (r, 0) = F(r)
S EC.3.3
DIFFERENCE M ETHODS
87
I t is a ppropriate to compare Eqs. (3 .2. 17) a nd (3.2 . 19). Although the symbolism used in these two equations is very similar, the two equations are really
quite different. Equation (3 .2.17) is displayed more clearly by Eq. (3 . 1.16) for
which temperatures are given starting at time t I a nd continuing up to time
tM + .1' Moreover, there is only one heat input t o the body; namely, the unknown surface heat flux, and also the initial temperature is t he uniform value of
To . Equation (3 .2.19) is much more general and can include the conditions for
Eq . (3.2.16) o r Eq. (3.2.17 ), ~15' seKiilg 111 :;±wla:1fl:,_o +o . E quation
(3.2.19) can also cover the case of an arbitrary temperature distribution at
time t M  1 a nd a k nown heat flux at x =L [ or o ther nonhomogeneous conditions
as shown in Eqs. (3.2.23) a nd (3 .2.24)].
T he eI>(r, t) problem is the solution of the problem of aT/aqc where q (t)=qc'
See Section 1.6. T he eI>(r, t) problem is
~~
r ar
(r
ael» = pc ael>
or
at
_ k ael>
or
k
I
=1
(3 .2.26)
(3.2.27a)
r =G
~~Lb =hel>(b, t)
cf>(r, 0) = 0
(3.2.27b)
(3 .2.28)
Note that the only nonhomogeneous term is a t the heated surface, r =a . The
el>i value is the solution of Eqs. (3.2.26) (3.2.28) for the sensor location and
time t i '
T he symbol t MI, .. =o represents the calculated temperature for the model
at time tM for the estimated heat flux c omponents q l, q2, ·· ·' q MI ' b ut qM is
set equal to zero. F or example, for the model given by Eqs. (3.2.23)(3.2.25),
known values o f g(r, t) a re used in Eq. (3.2.23), known values of T,.,(t) are used
in Eq. (3.2.24b), and the known initial temperature distribution, F(r), is used; the
temperature at the sensor location is calculated a t time tM with these conditions
and with the heat flux, q(t), equal to q I for 0 < t < I lt, q2 for I lt t o 21lt, . . . , qM  I
for t M 2 < t < t M I , a nd q = 0 for t M 1< t < tM' T he temperature
t M+ d, .. : , . .. , =0 is similarly found a t the sensor location a nd a t t M + I with the
known qi values for i =1 , 2, . .. , Ml a nd also q M=qM+I=O. F or the case
g(r, t) = 0 in Eq. (3 .2.23), T,.,(t) = To in Eq. (3.2.24b) and T (r, 0) = To in Eq. (3.2.25),
the t lq=o vector is simply given by
M I
L
qillel>M  i+ To
3.3
DIFFERENCE M ETHODS
Duhamel's theorem is a powerful technique for solving a wide variety of linear
heat conduction problems. Unfortunately, many situations exist in which either
(I) this technique becomes cumbersome o r (2) it is n ot applicable. The most
severe limitation is for nonlinear problems such as for bodies with temperaturedependent thermal properties. In many heat transfer problems, the temperature
change of a body exposed to a heat flux is sufficiently large that the changes in
thermal properties are appreciable. Fortunately, numerical methods can be
used t o convert the nonlinear partial differential equation of heat conduction
into a system of linear algebraic equations involving the temperature a t discrete
locations.
3.3.1 F inite C ontrol V olume P rocedure f or C onstant
P roperty P lanar G eometries
Many methods exist for discretizing partial differential equations. T wo p opular
ones are (1) finite differences a nd (2) finite elements. An alternative approach
in which conservation o f energy is applied directly to control volumes is a dopted
here. The finite control volume (FCV) procedure uses a control volume o f
arbitrary and finite size as in Figure 3.2 a nd the integral form of the energy
equation is applied to this control volume. Conservation of energy for a stationary solid requires that the sum o fthe rate of heat flow across the control volume
boundary and the time rate of change o f energy within the control volume is
equal to the rate a t which energy is produced within the control volume. In
equation form, the conservation o f energy can be written as
i =1
ff q.dA+ fr fff pedY fff e'"dY
M I
Tlq =o=
L
qillel>Mi+1
+ To
i =1
M I
L
i= 1
q illel>Mi+rI+To
=
(3.2.29)
A
r
(3.3.1)
r
where A a nd  Yare the closed surface area and volume respectively o f the control
volume. q is the heat flux vector leaving the control volume surface, e is t he
88
C HAP.3
M ETHODS F OR DIRECT HEAT C ONDUCTION PROBLEMS
Region of
interest
FIGURE 3 .2 Typical finite control volume.
specific energy (energy per unit mass), and em is the energy production rate per
unit volume. In the control volume shown in Figure 3.2, a p art o f the control
volume boundary coincides with the boundary o f t he body, a nd t he remaining
control volume boundary is completely within the body. Consequently, the
surface integral o fthe h eat flux can involve both heat flux from external sources
(boundary conditions) and heat conducted t o t he given control volume from
all neighboring control volumes.
I n its present form, Eq. (3.3.1) is too general to demonstrate the details o f the
FCV procedure. Let us restrict o ur a ttention to onedimensional p lanar p roblems and divide the body o f interest into increments o f e qual size (called elements) as shown in Figure 3.3. Each element has its own thermal properties
a nd a n element boundary can also correspond t o a material interface. The
control volume boundaries a re chosen t o lie a t t he midplane o f each element;
this selection o f c ontrol volume boundary is somewhat arbitrary. T he choice
for control volume boundaries permits two different materials t o be present
within a single control volume. Each element contains two nodes, o ne a t each
boundary o f the element which are identified by a (. ) i n Figure 3.3.
First, the general conservation o f energy Eq. (3.3.1) is applied t o t he control
volume surrounding node 1; Figure 3.40 provides an expanded view o f this
control volume. T he result is
x , p osition_
~~jr r  I
I
q (t)
 •
Element 2
:
I
,=~~':''::'''';=::J' Control
,
volume
I
boundary
N ode: Node
I
I
Element boundary
2
3
I
~~~~~~~~~~~~~~
3
Lx
FIGURE 3 .3
boundaries.
q~t)
'======~=:::~
N l
N
2
Onedimensional model for planar geometry showing element and control volume
r  ,
1
q(t)+
I
+k :!' I
1
I
1
""
liT
" , _ __ J
1
( a)
_ k a TI
 q(t)+ !!..
ax Axl2
dt
( Ax/2
Jo
p cTdx=O
(3.3.2)
where the energy content is given by e =cT a nd e"'=O. F or t he arbitrary interior
node j a nd the backface node N (see Figure 3.4b a nd c, respectively), application
o f the conservation o f energy equation yields
r '
• +. +.
k:;1 L _____ lk :!'I . .
"" I
X jT
j l
a TI
aTI
d
k+k+ax XJ+ A x/2
ax x r A x/2 dt
f ,XJ+Ax/2
pcTdx=O
(3.3.3)
1
I
.J
j
% J+T
j +l
(b)
x r A x/2
a TI
d
qN(t)+k  a
+d f ,XN
pcTdx=O
x X NAx/2
' xNAx/2
(3.3 .4)
T his control volume approach insures that the heat conducted o ut o f o ne
control volume is c onducted into the adjacent control volume. In order t o
evaluate the integrals and derivatives in Eqs. (3.3.2)  (3.3.4) a nd convert them
(e)
FIGURE 3 .4 Expanded view o f various con·
trol volumes. (a), Control volume surrounding
surface node I ; (b). control volume surrounding
arbitrary interior node j ; (c). control volume
surrounding backface surface node N .
89
90
C HAP.3
M ETHODS FOR DIRECT H EAT C ONDUCTION P ROBLEMS
into useful computational algorithms, it is necessary to make some assumptions
about the temperature profile within each element. Many assumptions can be
made concerning the element temperature profile; however, the more physics
exercised, the better. As an example, for a steadystate, planar geometry with
specified temperatures as boundary conditions and constant k, the temperature
profile is linear with position. Thus for the transient case, a linear temperature
profile is assumed for the element. F or an element bounded by X j  l ~X~Xj,
the assumed temperature profile is
T ( x )=T. l
r
( xXj)
T. ( xxj_d
+ J ( Xj Xj_ d '
Xj)
( Xjl 
X jl
~X~Xj
(3.3.5)
where 1 j= T (xj) is a node point (or nodal) temperature. Equation (3.3.5) is
written in the form of a Lagrangian interpolation polynomial. The local heat
flux is determined by differentiating Eq. (3.3.5) and using Fourier's law
k aT
( 1j1j d
q (x)= ;=k
, X j_l<X<Xj
(IX
X jXjl
(3.3.6)
The linear element temperature profile produces a constant heat flux within
elements; however, the heat flux changes from element to element.
Assuming that the volumetric heat capacity (pc) is a constant for each element,
a numerical expression for the energy storage term requires the evaluation of
integrals of the form T(x)dx. I t is convenient to perform these integrals over
portions of elements and then add the appropriate contribution to the control
volume energy balance. The two integrals that occur are evaluated using Eq.
(3.3.5) t o obtain
J
(3.3.7)
(3.3.8)
Note that Eqs. (3.3.7) and (3.3.8) are proportional to the average element
temperature over h alf of the element.
From the foregoing linear element temperature profiles, the finite control
volume energy balance for nodes 1, j and N, are
k
( T1 T2 )
Ax d
Ax
 q(t)+pc 2 "dt ( iTl + iT2 )=0
k (1j1j+d _ k(1jl1j)+ Ax !!.(1T,
Ax
Ax
pc 2 dt 4 r
+ pc
Ax d
2dt
q N(t)k
(i1j+i1j+d=0,
( TN1TN)
Ax d I
Ax
+PC
(4 TNl+i TN)=0
2dt
These three equations represent a system of N firstorder ordinary differential
equations in time. By assuming an element temperature profile and integrating
out the spatial dependence, the F ey procedure yields a system of ordinary
differential equations (as opposed to a single partial differential equation).
Note that each of the energy storage terms in Eqs. (3.3.9) (3.3.11) contains a
contribution from a node outside the control volume; this occurs because the
linear element temperature profile depends on two nodal temperatures, one
of which lies outside the control volume of interest. The foregoing results are
sometimes referred to as distributed capacitance effects because the thermal
capacitance of a control volume may be distributed over as many as three nodes.
The finite difference method is another. technique for arriving at a system
of ordinary differential equations representing heat conduction within a solid.
This procedure is also called lumped capacitance because all of the control
volume thermal capacitance is lumped at the center node point for all interior
control volumes. This is accomplished by replacing (i 1j _ I + i 1j) by 1j and
( i1j+i1j+d with 1j in Eq. (3.3.10); for surface control volumes, ( iT1 + !T2 )
becomes TI and ( iTN 1 + iTN) becomes TN' Additional details on difference
methods are available in Myers,9 Smith, I S Patankar,21 Anderson et al.,22 and
S hihP
Another method for the numerical solution of heat conduction problems is
the finite element (FE) method. I1  14 Many methods can be used to develop
finite element type systems of differential equations. The Galerkintype method
of weighted residuals is one of the simpler approaches and requires that the
integral weighted average of the partial differential equation be zero over some
domain. F or example, for a slab
i
L
(aT aT)
2
w(x)  ex dx=O
o
at
O X2
(3.3.12)
where L is the thickness of the slab and w(x) is a weighting function to be
specified. There is an infinite number of possible weighting functions; the trick
is to choose a weighting function that gives good results. A popular weighting
function is the socalled tent function which is shown in Figure 3.5 :
X j+lX
X j+l
J
j =2,3, . .. , N 1
91
DIFFERENCE M ETHODS
(3.3.9)
"JT,)
1 +4
SEC. 3.3
= 0,
(3.3.10)
(3 .3.11)
 x/
all other x
(3.3.13)
Note the similarity between the weighting function Wj(x) and the Lagrangian
interpolation polynomial in Eq. (3 .3.5), and that the weighting function is
different for each node j . Since the weighting function Wj(x) is zero over a large
92
C HAP.3
M ETHOOS FOR DIRECT H EAT C ONDUCnON PROBLEMS
w (x)
SEC. 3.3
DIFFERENCE M ETHODS
93
All o fthe F CV, LC, a nd F E results can be p ut in a similar form :
d
(X
q(t)
d (fJT, + 1'Tz)=  A2 (T,  T2 )+ :.t
uX
p cux
d
 (1'TN'
dt
(X
+ P TN)=Al ( TN, uX
qN(t)
TN)  
p cAx
(3.3.19)
(3 .3.21)
where fJ a nd l' h ave t he following values :
FIGURE 3.5 Weighting runction ror Galerkintypc method o r weighted
residuals.
P
LC
F CV
FE
p ortion o f t he d omain O~x~L, Eq. (3.3.12) c an be written as
"}
i
+1
" }I
(aT aT)
2
Wj(x) 
at
(X
2
ax
d x=O
(3.3.14)
Utilizing t he l inear element temperature profile assumption a nd E q. (3.3.13) for
Wj(x), t he F E e quation for a n a rbitrary interior node j c an be written as
k ( 1)1)+1)  k ( 1),  1)
A x!!.,
z
Ax
Ax
+ pc 2 dt h 1), +"3"1)
Ax d
+ pc T di ( i1)+l1)+d=O, j =2, 3, . .. , N l
T he o~ly difference between the F CV a nd F E results is the different weighting
coefficients o n t he t hermal capacitance terms; the F CV gives a greater weighting
t o t he center temperature 1). T he F E e quations for t he surface nodes are
l'
P+1'
1/8
1/6
1/2
1/2
1/2
o
(3.3.22)
E quations (3.3.19)  (3.3.21) a re restricted t o t emperatureindependent thermal
properties b ut t he c oncepts c an readily be extended t o Tvariable cases.
3.3_2
(3.3.15)
1/2
3/8
2/6
O ther B oundary C onditions a nd M aterial I nterfaces
In Section 3.3.1, a p lanar b ody with c onstant p roperties subjected t o a specified
heat flux b oundary c ondition was considered. A heat flux given by a convective
heat transfer coefficient a nd a t emperature difference driving potential is also a
common b oundary c ondition. F or t his' l atter case t he energy balance o n t he
control volume s urrounding n ode I c an b e obtained directly from Eq. (3.3.19).
Ax
Ax d
q(t)+pc Tdi<1 T, +1 Tl)=O
(3.3.16)
d
(X
h
 (PT, + 1'Tl )=   2 ( T,T2 )+(T",,T,)
d
t
Ax
p cAx
(3.3.17)
k
( T1 T2 )
Note t hat convective b oundary c onditions involve the surface temperature.
T he F CV a pproach c an be readily extended to elements o f n onuniform size
and differing material properties; such a control volume is shown in Figure 3.6.
L emmon a nd H eaton'6 d emonstrated t hat if the weighting function w(x) is
chosen t o be t he D irac delta function
Wj
~
\

()
~
x = u(xXj)=
{I,
0,
X =X·
}
all o ther x
Material
interface
I'
( e+ 1)
(3.3.18)
t hen t he F E results for interior nodes a re identical t o t he lumped thermal
capacitance (LC). This method is also called the collocation method. Additional
details o f the F E m ethod c an be found in Zienkiewicz, 1 2 N orrie a nd deVries ' 3
Becker et aI., ' 4 a nd B aker. 24
'
jl
j +1
j
FIGURE 3.6 C ontrol volume for material interface.
(3.3.23)
94
C HAP.3
M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS
T he two adjacent elements are. denoted (e) a nd ( e+ 1) with each having its own
temperature dependent properties (pc)le) a nd kle ). Applying an energy balance
kle + l ) 1 )1)+1 _ k le ) 1)I1)+(pc)e(XrXj_I)~(1'1)_1 + /11))
X j+lx j
X jXj_1
dt
+ (pc)le+ I)(Xj+ 1  Xj) ~ (/11)+1'1)+ d=O
Expanding U(l + I\l) in a Taylor series about time l, gives
dU(l)
u (t+l\t)=U(l)+ """dt  I\l+ . ..
(3.3.24)
(3.3.25)
where the subscript R indicates right half o f element (e). Similarly, for the left
half o f element ( e+ 1), the average temperature is
f r+ I ) = 11)+*1)+ 1
96
DIFFERENCE METHODS
(3.3.28)
If the terms 0 (l\t 2 ) are ignored, then u(t + I\t) is given by
T he element thermal conductivity is evaluated a t the midplane of each element;
for a linear element temperature profile, the appropriate temperature for element (e) is f le)=(1)_1 +1))/2. There are several choices for the appropriate
temperature a t which t o evaluate the volumetric heat capacity. The simplest
procedure is to use the average element temperature to evaluate (pc)le) and
assign this value to both halves o f the element. A slightly more complicated
approach is to compute the average temperature of each half of an element and
evaluate the volumetric heat capacity a t this temperature. F or example, the
right half o f element (e) contributes t o control volume j a nd the average temperature of this half of element (e) is
f~) = l[fle ) + 1)] = 1[t(1)1 + 1)) + 1)]
= i1)+*1)1
SEC. 3 .3
(3.3.26)
At this point, a word of caution is a ppropriate for the evaluation of temperaturedependent thermal properties. Most thermal properties are known only within
±5% at best. Consequently, it is not particularly fruitful to develop highly
sophisticated techniques to precisely handle temperaturedependent thermal
properties that are known imprecisely. Many times, the simplest approach is
more than adequate.
u(t + I\t) z u(t)+ f [u(t), l ]l\t
(3.3.29)
I f u(t.) is denoted as u ., then Eq. (3.3.28) becomes
U .+ 1 = u.+
f (u., t .)l\t=u.+
~~I.l\l
(FO)
(3.3.30)
The algorithm given by Eq. (3.3.30) is shown in Figure 3.7; it is simply the
extrapolation of the slope a t p oint (un' tn) t o (un+ I , tn+ d. As an example of the
Euler method, it is applied to Eq. (3.3.19) with /1=1/2, 1 '=0 (LC). Since the
subscript o n T represents node number, a superscript is used t o indicate time,
2
nl\l
Ti)+~
p cux
T j+ 1 = T j  2p(Tj 
(3.3.31)
where p =a.l\t/l\x 2 is a F ourier number based o n the timespac~ .grid. The
Euler method for the LC formulation is also known as the explIcit method
since T j+ 1 in Eq. (3.3.31) can be written as an explicit function of known temperatures only at time n ; a nother name is the forward difference method because
the time derivative is replaced by a forward time difference. Explicit results are
also obtained from Eqs. (3.3.20) and (3.3.21) for the LC formulation.
As can be seen from Figure 3.7, the Euler method underpredicts the temperature if the true solution is concave upward. In an attempt to improve stability, a
backward time difference (BO) is often used; analogous to Eq. (3.3.28), the
I
/
3 .3.3 N umerical T echniques f or S olving S ystems o f
F irstOrder O rdinary D ifferential E quations
u (t)
'/
/
/ '
I n principle, an analytical solution of the system of energy balance Eqs. (3.3.19) (3.3.21) is possible. However, numerical techniques offer solutions t o more
general problems and can be readily extended to the temperaturedependent
thermal properties. Numerous numerical techniques are available.
O ne o f the simpler techniques is the Euler (or forward difference, FO)
method which can be developed from a Taylor series expansion. Let u(t) be a
function of time that satisfies the firstorder ordinary differential equation
du
d t = f (u, t)
sIOpeU~+l
I
(3.3.27)
FIGURE 3 .7
methods.
Extrapolation of
slope u~ = { (un,tn)
Schematic o f the Euler (forward difference) and implicit (backward difference)
96
C HAP.3
M ETHODS FOR DIRECT H EAT C ONDUCTION PROBLEMS
2
(
3.3.32)
Neglecting terms of 0(llt 2), the backward difference integrator gives
(3.3.33)
The backward time difference integrator is also known as the fully implicit
method because the righthand side contains the unknown Uo + I . A central
difference (CD) time derivative is known to be more accurate (under certain
conditions) than either the forward o r backward difference and can be derived
from the following Taylor series expansions:
Ilt u" ( Ilt)2
uO+I/2=uo+u~2+2~ 2 + ...
(3.3.34)
,
(3.3.35)
U +I/1=UO +IUO +1
O
97
DIFFERENCE M ETHODS
d
a.
0+ 1
0  (RT, + T )0+1 =  0 _ (To+l_ T~+I)+O q
dt I ' 1 Y 2
Ilx1 1
pcllx
equation is
, ..
" Ilt
U o=Uo+luo+lut+Uo+1 2 1···
SEC. 3.3
Ilt U:+ 1 ( Ilt)2
2 +T! 2 _ . ..
O' d (fJ TI+yT1)0=8' ..a. 2(T~T~)+0' q:
d
uX
p cux
t
where 8' = 1  0 a nd the superscript n means evaluated a t time to· Adding the
foregoing two equations gives
d
d
a.
0  (fJT1+yT1)0+ 1 + 8' d (fJT1+ yT1)0=  0 A 2 (T~+ 1  T~+ I )
dt
t
uX
 8'  ; (T~ 
Ilx
T~)+~ (Oq"+ 1 + 8'q")
p cux
The general integrator can be applied directly to the lefthand side ofEq. (3.3.39)
with the result
Tft+ 1 _ TO)
( T o+1  TO)
fJ ( 1
1+
2
2
Ilt
Y
Ilt
(3.3.40)
Equating these two expressions and solving for Uo + 1 gives the central difference
integrator
(3.3.36)
which is also known as the CrankNicolson method. von Rosenberg 1 8 describes
the CrankNicolson procedure as successive applications of the forward and
backward integrators. The forward, backward, and central difference integrators can all be combined into a general integrator (GI) as follows:
U +1 = uo+ [OU~+I +(10)u~]llt
o
(3.3.37)
1
0 = 1,
q "+e= OqO+ 1 +O'q"
(3.3 .41)
Equation (3.3.40) can be simplified to yield
Ilt
p cux
(pO+fJ)T~+I_(pOy)T~+ 1 =  (pO' fJ)T~ +(pO' +y)T~ + "A q0+e
(3.3.42)
Applying the G I t o Eqs. (3.3.20) and (3.3.21) gives
(3.3.43)
+2(fJ  pO')Tj+(p8' +y)Tj+ I , j =2, 3, . .. , N l
forward difference o r Euler
0 =2' central difference o r CrankNicolson
where for convenience the following definition is employed
(pOy)Tj~: +2(pO+ f J)Tr 1 _(pOy)TH: =(p8' + y)Tj_1
The three special cases are:
0=0,
(3.3.39)
(3.3.38)
Ilt
(pOy)TN~11 +(pO+fJ)TN+I =(pO' + y)TN 1 +(fJ  pO')TN ; qN+6
p cux
backward difference o r fully implicit
3 .3.4 G eneral F orm o f D ifference E quations f or H eat
C onduction i n P lanar B ody
The general integrator given by Eq. (3.3.37) for firstorder equations converts
an ordinary differential equation into a single difference equation. The next
step is to convert the system o f ordinary differential equations, Eqs. (3.3.19) (3.3 .21), into a system o f algebraic difference equations. The procedure for
applying the GI, Eq. (3.3 .37), is to evaluate the differential equation at times
n + 1 and n, multiply these equations by 0 a nd ( 1 0), respectively, and add the
two equations. Using Eq. (3.3.19) as an example gives
(3.3.44)
Equations (3.3.42)  (3.3.44) form a linear system of N algebraic equations with
N unknowns. These algebraic equations have a very special structure that is
known as tridiagonal and is best visualized by writing them in matrix form for
constant thermal properties,
pO+fJ
 (pOy)
 (pOy) 2(pO+:P)
 (pOy)
TN~\
 (pOy) pO+fJ
TN 1
+
C HAP.3
98
fJ  pO'
y+p(:J'
y+p(:J'
2(fJ  p(:J')
M ETHOOS FOR DIRECT H EAT CONDUCTION P ROBLEMS'
y + p(:J'
y+p(:J'
[
y+p(:J'
fJ  p(:J'
][
T~]
T~
At
S EC.3.3
Comments on Accuracy. I n order t o o btain the desired accuracy it is necessary
to use a sufficient number o f nodes N and t o t ake a "small" time step. In many
cases a value of N equal t o 20 is sufficient. The time steps must be a t least as
small as the time steps in the measured temperatures but can be considerably
smaller; it is convenient t o choose the computational time step, At, so that the
time step between measurements, At m• as , is equal t o
[qn+8]
0
6
i'NI + pcAx
T'N
q'N+ 8
(3.3.45)
T he tridiagonal coefficient matrix has nonzero terms along only three diagonals.
This tridiagonal structure allows for a very efficient specialization of the Gauss
elimination procedure for solving linear systems o f algebraic equations. The
details of the algorithm, known as the Thomas algorithm, are presented in
Section_ .3.2.
6
T he form o f Eq. (3.3.45) shows the explicit dependence o f the new temperatures on the old temperatures and the boundary conditions. From a computational point o f view, it is more convenient t o lump the righthand side of
Eq. (3.3.45) into a constant vector.
The Euler integrator (forward difference, (:J = 0) for the LC (fJ = 1/2, y = 0) is
worthy o f special note. F or this case, the coefficient matrix for the unknown
temperatures become diagonal; consequently, each unknown temperature
T7+ I can be expressed explicitly in terms o f known temperatures ( TiI' Ti,
Ti+d. F or either the FCY (fJ=3/8, y =I/8) o r F E (fJ=2/6, y =I/6) method in
conjunction with the Euler integrator «(:J = 0), all three diagonals o f the coefficient
matrix are nonzero; hence, the term explicit integrator is inappropriate for
these cases.
EXAMPLE 3.2. Give the difference equations for the CrankNicolson approximation
for a plate that is heated a t x =O a nd insulated a t x =L. Let Ax = L/3 a nd use a lumped
thermal capacitance. The thermal properties are independent of temperature.
Solution. The difference equations can be obtained using Eq. (3.3.45) for N = 4 (see
Figure 3.3); there are two boundary nodes and two interior nodes. From Eq. (3.3.22)
for the lumped capacitance approximation, fJ equals 1/2 and y equals O. F rom Eq. (3.3.38)
for the CrankNicolson approximation, 8 equals 1/2. F or the insulated boundary,
q~+ 1 /2 = 0.
Using these values in Eq. (3.4.45) and multiplying the first and last equations by 2
gives the four difference equations in matrix form as,
[ ,+,
t
~.
[''
p/2
0
0
p
p +l
 p/2
0
p
0
l p p/2
l p
p/2
p
0
or·]
oJrJ [~."'J
0
 p/2
p +l
p
o
T j+1
o
Ti
p/2
l p
T:
T4
A t= Atm. ..
i
(3.3.46)
where i is a positive integer. In order t o choose the space step, Ax, a s well as
the time step, At, experience is needed for each particular case because their
values are affected by many factors including the time period of interest, variation of thermal properties with temperature, and the time variation o f the
surface heat flux. F or more discussion, see books o n numerical methods such
as References 6, 8, 9 a nd 15.
3.3.5
S tandard F orm o f T emperature E quations f or I HCP
The purpose of this section is t o demonstrate that the difference equations
produced by the FD, FE, and FCY methods can be used t o o btain the same
standard equation, Eq. (3.2.19), that was obtained from Duhamel's theorem.
The discussion in this section is confined t o the linear I HCP b ut nearly the
same approach can be used for the quasilinear analysis which is discussed in
Section 6.2.
T he analysis can be performed using either algebraic o r matrix notation
but the latter analysis is preferred because it is more compact and general.
I f n + 1 is replaced by the index M, Eq. (3.3.45) can be written as
(3.3.47)
where A is the square matrix o n the left of Eq. (3.3.45), B is the square matrix
on the right of Eq. (3.3.45) and C = A t/ pcAx. The vectors TM a nd qM are
.
~~[::l ~[f"')
(n~)
T i+ 1
 p/2
99
DIFFERENCE M ETHODS
p +l
w herelhe subscripts relate to the space nodes and the superscripts t o time. The
=
gM vector represents any known energy sources due to an applied heat flux
T~+I
2 At
+ pcAx
at x =L, volumetric heating, o r convective heat transfer. I f the CrankNicolson
method of solution is selected, (:J = 1/2 and
0
0
0
q M+81 = qMI/2
0
which is the surface heat flux evaluated a t time tM 
1/ 2
= (M  1/2)At. Observe
100
C HAP.3
M ETHODS FOR DIRECT H EAT C ONDUCTION PROBLEMS
t hat q MI/2 m eans exactly the same heat ftux component as denoted q M in the
Duhamel procedure [see Eq. (3.2.2)].
Equation (3.3.47) is solved for TM by mUltiplying by A  I t o o btain
TM = DTM I + EqM + A  lgM
(3.3.49a)
where
D =A  'B,
E =CA '
SEC. 3.3
DIFFERENCE M ETHODS
E quations (3.3.53a,b,c) are vector equations for the temperatures a t each
node in the body. Temperature sensors are located a t only a few o fthese nodes, in
general. F or simplicity, a single sensor is considered a nd t he space dependence
is indicated by a subscript K. T hen Eqs. (3.3.53a,b,c) c an be written for the K th
node as
(3.3.54a)
(3.3.49b,c)
(Though it is n ot a g ood numerical p rocedure t o find the inverse o f A when
~lc~lating T M it is convenient here t o d o so symbolically.) Notice t hat TM
,.
IS a h near funchon o f qM . Replacing M in Eq. (3.3.49a) with M + 1 gives
(3.3 .50)
T M+ 1 ; "DTM+ EqM+l + A  lgM+l
a nd t hen using Eq. (3.3.49a) for TM yields
T M+ I = D 2T M 1 + DEqM + EqM+l + DA  lgM+ A  lgM+1
(3.3.51)
(3 .3.54b)
M+2
T"
=tM+21,.,=,.,+1 =,".>=0+ X 31q M+ X 21q M+I + X I lqM +2
"
( 33 .54c)
.
where (J in q M + 9 I is c hosen t o be 1 t o simplify the notation. T he s ymbols
X I I' X 21, a nd X 31 d enote sensitivity coefficients a t l ocation K a nd a re n ot
functions o f t he time index M a s indicated by the corresponding matrix terms
of E, D E a nd D 2E. T he sensitivity components X I I' X 21, X 31 a re given by
(3.3.54d)
This expression is l inear in b oth qM a nd qM + I. R epeating the process for M
replaced by M + 2 in Eq. (3.3.49a) gives
T M+l=D3TMl + D1EqM + DEqM+l + EqM+l
+ D2A  lgM + DA  lgM+I + A  lgM+l
101
(3.3.54e)
(3.3.52)
which is linear in qM, qM + 1 a nd qM + 2.
D ue t o t he linearity o f Eqs. (3.3.49a), (3.3.51), a nd (3.3.52), they c an be written
as
(3.3.54f)
TM = TMlq.,=o+EqM
T M+ I = TM+ Ilq"=q.,+'=o+DEqM +EqM+ 1
(3.3.53a)
(3.3.53b)
Equations (3.3.54) can be written in the standard matrix form with the K
subscript o n T a nd Tlq=o o mitted as
(3.3.55a)
T M+l=tM+llq"=q"" _q.,+1=0+D1EqM + DEqM+I + EqM+l
(3.3.53c)
where
where the temperatures evaluated a t qM, . .. , equal t o 0 a re defined by
TMlq., =o=DTMI+A  lgM
T=[~:+I] ,
T
(3.3.53d)
M
T M+ Ilq"=q.,+1 =0=D1T M 1 + DA  lgM + A  lgM+ I
(3.3.53e)
TM+llq"=q.,+I=q.,+> =0=D3TMl + D1A  lgM + DA  lgM+1 + A  lgM+l
(3.3.55b,c)
+2
•
(3.3.55d)
(3.3.53f)
~otice t hat. the temperature components defined by Eqs. (3.3.53d,e,f) are
hnear functIOns o f only the temperature distribution a t t ime t M  1 a nd t he
known heat sources a t times t M , tM + I , a nd tM + 1 • T hey are equal t o t he corresponding temperatures given by Eqs. (3.3.35a,b,c) if the future heat ftux vectors
are equal t o zero. Equations (3.3.35a,b,c) display the important characteristic
o f linearity in t he t emperatures defined by Eqs. (3 .3.35d,e,f) a nd t he h eat fluxes
qM, qM+I, a nd qM+2
Whereas the temperature vectors in Eqs. (3.3.47) (3.3.53) are for the nodes a t
time t , t he temperature vectors in Eqs. (3.3.55) are for different times. E quation
M
(3 .3.55a) is called the s tandard form for the temperature because the same
equation was derived using Duhamel's theorem [see Eq. (3 .2.19)]. In Eqs.
(3 .3.55b,c,d) the superscript M d enotes the same time dependence as the M
subscripts in Eqs. (3 .2.20a,b) a nd (3 .2.22a). T he X 1 1, X 2 " a nd X 31 values
1 02
C HAP.3
M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS
correspond to
17.
used in the Duhamel theorem method and the numerical values of X I I and
X 2 1 and llcp I , and so on, will be nearly identical if the finite difference
(or F E or FCV) calculations are accurate. The symbol cpj represents the temperature rise a t time t l and node X K for the numerical solution of Eq. (3.3.47) for
q M+91 = 1 for i =M = 1 ,2,3, . .. a nd gM = 0.
Since the same standard form, Eq. (3 .2.19) o r (3 .3.55a), is obtained for T,
I HCP algorithms can be developed that do not depend on the type o f numerical
approximation of the transient heat conduction equation. Advantage is taken
of this fact in Chapter 4 to derive general algorithms. In Chapter 5 the algorithms
of Chapter 4 are implemented using Duhamel's integral, and Chapter 6 uses
the finite control volume method.
llcpo,
Keltner, N. R. a nd Beck, J. V., Unsteady Surface Element Method, J . Heat Transfer 103, 759 764 (1981).
2. Beck, J. V. a nd Keltner, N. R., Transient Thermal Contact o f T wo SemiInfinite Bodies O ver a
C ircular Area, in Spacecraft Radiative Transfer and Temperature Control, T. E. H orton, ed.,
Vol. 83 o f Progress ill Astronautics and Aeronautics (1982), AIAA, 61  82.
3. Beck, J . V., Schisler, l . P. a nd Keltner, N. R., Simplified Laplace Transform Inversion for
Unsteady Surface Element Method for Transient Conduction, A I A A Journal, 22 , 1328 1333,
(1984).
4.
Myers, G. E ., T he C ritica l Time Step for Finite Element Solutions t o T woDimensional Heat
Conduction Transients, A SME J . Heat Transfer 1 00,120 127, (1978).
18. von Rosenberg, D. U., Methods for the Numerical Solulion o f Partial DijJerential Equation .•,
Americal Elsevier, New York, 1969.
19. Richtmyer, R. D. a nd M orton, K. W ., Difference Methods for Inirial Value Problems, 2nd ed.,
Interscience Publishers, New York, 1967.
20 . Beck, J. V., G reen's Function Solution for Transient Heat Conduction Problems, Int . J . Heat
Mass Transfer, 2 7,1235  1244 (1984).
21. P atankar, S. V., Numerical Heat Transfer and Flllid Flow, Hemisphere, W ashington, 1980.
22 . Anderson, D. A., T annehill, J. C., a nd Pletcher, R. H ., Computational Fluid Mechanics and
Heat Transfer, Hemisphere, Washingto n, 1984.
23 . Shih, T. M., Numerical Heat Transfer, Hemisphere, Washington, 1984.
24. Baker, A. J., Finite Element Compulational Fluid Mechanics, Hemisphere, 1983.
PROBLEMS
3.1. O btain mathematical expressions for the surface temperature o f a
semiinfinite body that is initially at the uniform temperature of To
and then is exposed to the heat flux history given by,
REFERENCES
1.
1 03
PROBLEMS
Keltner, N. R. a nd Beck, J . V., Surface Temperature Measurement Errors, J . Heat Transfer,
lOS, 312  318 (1983).
qo, O <t<t l
q(t)=  qo, tl < t<2t l
0,
t <O a nd
l
t >2tl
A different expression is needed for each time interval. Plot the group o f
T  To ( k PC)1/2
qo
tl
versus t it 1 for t it I from 0 to 4. Why d o the temperatures change after
t = 2t 1 even though the net energy added is zero?
S.
Litkouhi, B. a nd Beck, J. V., Intrinsic Thermocouple Analysis Using M ultinode Unsteady
Surface Element Method, AIAA Paper N o . 83  1937 . (To be published in AIAA Journal.)
6. Ozisik, M. N. Heat Conduction, Wiley, New York, 1980.
7. Luikov, A. V., Analyt ical Hear Diffusion Theory, Academic Press, New York, 1968.
8. Arpaci, V. S., Conduction Hear Transfer, AddisonWesley, Reading, MA, 1966.
9. Myers, G . E., Analytical Methods in Conduction Heat Transfer, M cGrawHiII,"New York, 1971.
10.
11.
~
\ .)
C arslaw, H. S. a nd J aeger, J. C , Conduclion o f Heat in Solids, 2nd ed ., O xford, L ondon,
19~ .
Bathe, KlausJurgen, Finite Element Procedures in Engineering Analysis, PrenticeHall,
Englewood Cliffs, NJ, 1982.
12. Zienkiewicz, O. C , T he Finite Element M ethod, 3rd ed., McGrawHili, New York, 1977 .
13 . Norrie, D. H. a nd deVries, G., An Introduclion 10 Finite Element Analysis, Academic P re ss,
New York, 1978.
14. Becker, E. B., C arey, G . F., a nd O den, J . T., Finile Elements an Introduction, Vol. I , PrenticeHall, Englewood Cliffs, NJ, 1981.
IS. Smith, G . D., Numerical SolUlion o f Parlial Differential Equations, O xford University Press,
New York, 1965.
16. Lemmon, E. C. a nd H eaton, H. S ., Accuracy, Stability, and Oscillation C haracteristics o f
F inite Element Method for Solving Heat Conduction Equation, ASME P aper N o. 69WA/HT35 .
3 .2. T he heat flux history a t the surface of a body is
10 W/m2,
q(t) = { 0
2 s<t<4s
otherwise
The temperature a t X l is 400 K until t =2s, 401 K a t t =4s, and 406 K
a t t = 6s. The heat conduction problem is linear. Find the temperature
rise a t X 1 a t times 2 and 4s for a unit step rise a t t = 0 in the surface
heat flux .
3 .3. Find exact expressions for the surface temperature of a semiinfinite
body that is initially a t the uniform temperature of To a nd that is subjected to a heat flux that varies in a triangular fashion with time as
given by
40 t ,
l
q(t)= ~0(2tlt),
where 40 is a constant.
O <t<t l
tl < t<t 2 =2t l
otherwise
1 04
C HAP.3
M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS
105
PROBLEMS
20
w here G(x, l )=(kla)ocp(x, t)/Ol. (G(x, t) is a Green's function .)
a. Is it possible t o use a numerical approximation o f Eq. (a) when
G(O, t) is needed a t t = O? G ive a reason for your answer.
b. D erive the numerical approximation o f Eq. (a),
G ive expressions for the three time regions o f (a) 0 < 1< 1I , (b) 1I < 1<
211 a nd (c) 2 tl < t.
P lot t he surface temperature in a dimensionless form versus t it I
for t he r ange o f 0 t o 3.
M
3.4. A h eat flux is given by
T M=To+
L
qft H Mft+IL11
(b)
n =1
j
qG'
q (l) =
,,qG
qG + ql> _ tG ( 11G ),
t
a
H i == G (x, l i1/2) Ie
F or a semiinfinite body initially a t To c alculate the surface temperature
a t time tl> using Eq. (3.2.12).
. Use L11=0.5, IG= 1, qG= 1, q ,,=2, a nd k pc= 1. C ompare t he answer
with the exact value.
c.
U se Eq. (b) for t he l inear h eat flux case o f E xample 3.1 a nd c ompare
t he values.
3.9. S how t hat if CPi is given by
3 .5. Derive t he a pproximation o f t he Duhamel theorem given by
M
T M=To+
L
qftaM  ft+1
.. .. 1
where q i=q a t I i a nd
;
t hat t he s um
((I)
4>; ,
:=
0
M ",l
for c"'S 0
j=(cpl')2cpl~1 +CPj~2)~ ; =2, 3, .~
..
,....
a nd ' f'j IS e t emperature n se a t time t i for the heat flux q = I . T he
h eat flux is composed o f c onnected linear segments starting with
q o=O . ( The superscript 1 denotes a linear h eat flux " building block.")
T M=
1= cp\l)/M;
. .1.(1)'
T M + I = e  I TM + B(e<  l)qM",e("+ 1)<
3 .10.
3 .6. W rite a c omputer o r p rogrammable calculator p rogram for the numerical convolution as given by Eq. (3.2.12). Allow for a t least 10 qft a nd
10 L1CPi c omponents. T he i nput is t o b e To a nd t he qft'S a nd L1CP/s. T he
o utput is
Assume t hat t he temperature rise for a n u nknown surfa~e h eat flux is
L1 T (r, t) which is a k nown function. By t aking t he L aplace transform
o f Eq. (3.2.9), d erive
q
(t)==~I ~[L1T(r, I )]
s 9'[q,(r,
In
w here s is t he L aplace transform parameter. F or
TI , T2 ,···, 7 ;0'
1
Verify y our p rogram by checking Example 3.1.
c p=2 ( _ k
n pc
3 .7. A semiinfinite body initially a t t he uniform temperature o f z ero is
exposed a t t = 0 t o t he heat flux.
)1/2
,L1T
440(al)3/2
3a ( )1/2
k
n
find q(I).
q = (nt)1/2
3 .11.
Let a = I , k = I, a nd L1t = 1.
C alculate the temperature a t x =O using Eq. (3.2.12). Use a t least
10 time steps. T he e xact temperature is 1. C omment o n t he accuracy
o f t he procedure.
P artial Answers: 0.900, 0.784, 0.841, 0.865, a nd 0.880
q (l)G(x, 1 l)dl
U sing the Laplace transform show t hat t he t emperature rise in a body
exposed t o t he linearly varying heat flux,
q (l) == 401
is e qual t o
L1T(r, 1)==40
3 .8. An alternate way to write the convolution integral is
i f~
q~CPM  ft
c an b e written in t he s equential form
th
T (x, I) = To +
L
ft=1
(a)
f~ cp(r,I')dl'
where q,(r, I) is for a unit step heat flux starting a t 1 =0.
106
C HAP.3
METHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS
3.12. Using the divergence theorem of Gauss, demonstrate that Eq. (3.3.1)
becomes
iJe
iJt
3 .13. F or onedimensional (radial) cylindrical geometries, use the steadystate temperature profile as the element temperature profile and develop
the equations analogous t o Eqs. (3.3.9)  (3.3.11). Assume the control
volume boundaries are a t the midpoint of the elements.
3 .14. Repeat Problem 3.13 for onedimensional (radial) spherical geometries.
3 .15. Verify the finite element equation, Eq. (3.3.15).
A x=L/4
3 .17. a.
F or two nodes in a plate, one at x = 0 a nd t he other at x = L , give
the two CrankNicolson, lumped capacitance equations. The
surface at x = 0 is exposed to a heat flux a nd the surface at L is
insulated.
b. F or the temperature a t n ode 2, derive a n expression analogous to
Eq. (3.3.54a). Give T 1'l f "=0
a s a function of T~  1, T 1'  1, p, a nd A t/p cAx .
c. Give X 1 1 as a function o f p a nd A t/p cAx.
3 .18. T he onedimensional Green's function equation analogous to
Duhamel's theorem is given in Problem 3.8. Write the equation in the
form
tM)=To+~
E"
ex M
q(A.)G(x, tMA.)dA.
J.';
=To+kj~1 "_1
q(A.)G(x, tuA.)dA.
Approximate q(A.) for t j _ 1 < A. < t l by
t jA.
A .t j _ 1
q(A.) = qjI M
+ qj  ;;t"
where qj = q (tJ, t j = i At .
Derive the approximate relation,
M
T (x, t M)=To +
I
i)==~
12 (x, M, i ) == ex
k
L,
J.".
'
[ iql1It(X, M, i )
i =l
(a)
_I
G(X,IMA.)dA.
A
A.t G (x, tM  A.)d..1.
Show t hat E q. (a) c an be written as
M
T (x , t M)=To+
I
j =O
W hat are the a ju(x) expressions?
which is exposed t o a timevariable heat flux at x =O a nd a convective
boundary condition at x = L . T he properties are temperature independent. Use t he backward difference time approximation a nd lumped
capacitance.
T (x,
where
I I(x, M ,
V ·q+p  =e'"
3 .16. Give the difference equations in matrix form for a plate with
1 07
PROBLEMS
a jM(x)qMj