Unformatted text preview: cen58933_ch05.qxd 9/4/2002 11:41 AM Page 265 CHAPTER NUMERICAL METHODS
I N H E AT C O N D U C T I O N
o far we have mostly considered relatively simple heat conduction problems involving simple geometries with simple boundary conditions because only such simple problems can be solved analytically. But many
problems encountered in practice involve complicated geometries with complex boundary conditions or variable properties and cannot be solved analytically. In such cases, sufficiently accurate approximate solutions can be
obtained by computers using a numerical method.
Analytical solution methods such as those presented in Chapter 2 are based
on solving the governing differential equation together with the boundary conditions. They result in solution functions for the temperature at every point in
the medium. Numerical methods, on the other hand, are based on replacing
the differential equation by a set of n algebraic equations for the unknown
temperatures at n selected points in the medium, and the simultaneous solution of these equations results in the temperature values at those discrete
points.
There are several ways of obtaining the numerical formulation of a heat
conduction problem, such as the finite difference method, the finite element
method, the boundary element method, and the energy balance (or control
volume) method. Each method has its own advantages and disadvantages, and
each is used in practice. In this chapter we will use primarily the energy balance approach since it is based on the familiar energy balances on control volumes instead of heavy mathematical formulations, and thus it gives a better
physical feel for the problem. Besides, it results in the same set of algebraic
equations as the finite difference method. In this chapter, the numerical formulation and solution of heat conduction problems are demonstrated for both
steady and transient cases in various geometries. S 5
CONTENTS
5–1
Why Numerical Methods 266
5–2
Finite Difference Formulation of
Differential Equations 269
5–3
OneDimensional
Steady Heat Conduction 272
5–4
TwoDimensional
Steady Heat Conduction 282
5–5
Transient Heat Conduction 291
Topic of Special Interest:
Controlling the
Numerical Error 309 265 cen58933_ch05.qxd 9/4/2002 11:41 AM Page 266 266
HEAT TRANSFER 5–1 ·
g0
1d
— — r 2 dT + — = 0
—
r 2 dr
dr
k () r0
Sphere 0
dT(0)
——– = 0
dr I WHY NUMERICAL METHODS? The ready availability of highspeed computers and easytouse powerful software packages has had a major impact on engineering education and practice
in recent years. Engineers in the past had to rely on analytical skills to solve
significant engineering problems, and thus they had to undergo a rigorous
training in mathematics. Today’s engineers, on the other hand, have access to
a tremendous amount of computation power under their fingertips, and they
mostly need to understand the physical nature of the problem and interpret the
results. But they also need to understand how calculations are performed by
the computers to develop an awareness of the processes involved and the limitations, while avoiding any possible pitfalls.
In Chapter 2 we solved various heat conduction problems in various geometries in a systematic but highly mathematical manner by (1) deriving the
governing differential equation by performing an energy balance on a differential volume element, (2) expressing the boundary conditions in the proper
mathematical form, and (3) solving the differential equation and applying the
boundary conditions to determine the integration constants. This resulted in a
solution function for the temperature distribution in the medium, and the solution obtained in this manner is called the analytical solution of the problem.
For example, the mathematical formulation of onedimensional steady heat
conduction in a sphere of radius r0 whose outer surface is maintained at a uni·
form temperature of T1 with uniform heat generation at a rate of g0 was expressed as (Fig. 5–1) T(r0) = T1 1 d 2 dT
r
dr
r 2 dr
dT(0)
0
dr Solution:
·
g0
T(r) = T1 + — (r02 – r 2)
6k
·
4 π g0r 3
·
dT ———
Q(r) = –kA — =
dr
3 FIGURE 5–1
The analytical solution of a problem
requires solving the governing
differential equation and applying
the boundary conditions. ·
g0
k 0 and T(r0) T1 (51) whose (analytical) solution is
T(r) T1 ·
g0 2
(r
6k 0 r 2) (52) This is certainly a very desirable form of solution since the temperature at
any point within the sphere can be determined simply by substituting the
rcoordinate of the point into the analytical solution function above. The analytical solution of a problem is also referred to as the exact solution since it
satisfies the differential equation and the boundary conditions. This can be
verified by substituting the solution function into the differential equation and
the boundary conditions. Further, the rate of heat flow at any location within
the sphere or its surface can be determined by taking the derivative of the solution function T(r) and substituting it into Fourier’s law as
·
Q (r) kA dT
dr k(4 r 2) ·
g0r
3k ·
4 g0r 3
3 (53) The analysis above did not require any mathematical sophistication beyond
the level of simple integration, and you are probably wondering why anyone
would ask for something else. After all, the solutions obtained are exact and cen58933_ch05.qxd 9/4/2002 11:41 AM Page 267 267
CHAPTER 5 easy to use. Besides, they are instructive since they show clearly the functional dependence of temperature and heat transfer on the independent variable r. Well, there are several reasons for searching for alternative solution
methods. T ,h k = constant No
radiation T ,h h, T h = constant
T = constant FIGURE 5–2
Analytical solution methods are
limited to simplified problems
in simple geometries. An
ovalshaped
body 2 Better Modeling
We mentioned earlier that analytical solutions are exact solutions since they
do not involve any approximations. But this statement needs some clarification. Distinction should be made between an actual realworld problem and
the mathematical model that is an idealized representation of it. The solutions
we get are the solutions of mathematical models, and the degree of applicability of these solutions to the actual physical problems depends on the accuracy of the model. An “approximate” solution of a realistic model of a
physical problem is usually more accurate than the “exact” solution of a crude
mathematical model (Fig. 5–3).
When attempting to get an analytical solution to a physical problem, there
is always the tendency to oversimplify the problem to make the mathematical
model sufficiently simple to warrant an analytical solution. Therefore, it is
common practice to ignore any effects that cause mathematical complications
such as nonlinearities in the differential equation or the boundary conditions.
So it comes as no surprise that nonlinearities such as temperature dependence
of thermal conductivity and the radiation boundary conditions are seldom considered in analytical solutions. A mathematical model intended for a numerical solution is likely to represent the actual problem better. Therefore, the
numerical solution of engineering problems has now become the norm rather
than the exception even when analytical solutions are available. No
radiation Long
cylinder 1 Limitations
Analytical solution methods are limited to highly simplified problems in simple geometries (Fig. 5–2). The geometry must be such that its entire surface
can be described mathematically in a coordinate system by setting the variables equal to constants. That is, it must fit into a coordinate system perfectly
with nothing sticking out or in. In the case of onedimensional heat conduction in a solid sphere of radius r0, for example, the entire outer surface can be
described by r r0. Likewise, the surfaces of a finite solid cylinder of radius
r0 and height H can be described by r r0 for the side surface and z 0 and
z H for the bottom and top surfaces, respectively. Even minor complications in geometry can make an analytical solution impossible. For example, a
spherical object with an extrusion like a handle at some location is impossible
to handle analytically since the boundary conditions in this case cannot be expressed in any familiar coordinate system.
Even in simple geometries, heat transfer problems cannot be solved analytically if the thermal conditions are not sufficiently simple. For example, the
consideration of the variation of thermal conductivity with temperature, the
variation of the heat transfer coefficient over the surface, or the radiation heat
transfer on the surfaces can make it impossible to obtain an analytical solution. Therefore, analytical solutions are limited to problems that are simple or
can be simplified with reasonable approximations. h, T Simplified
model Realistic
model A sphere Exact (analytical)
solution of model,
but crude solution
of actual problem Approximate (numerical)
solution of model,
but accurate solution
of actual problem FIGURE 5–3
The approximate numerical solution
of a realworld problem may be more
accurate than the exact (analytical)
solution of an oversimplified
model of that problem. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 268 268
HEAT TRANSFER
z 3 Flexibility
Engineering problems often require extensive parametric studies to understand the influence of some variables on the solution in order to choose the
right set of variables and to answer some “whatif ” questions. This is an iterative process that is extremely tedious and timeconsuming if done by hand.
Computers and numerical methods are ideally suited for such calculations,
and a wide range of related problems can be solved by minor modifications in
the code or input variables. Today it is almost unthinkable to perform any significant optimization studies in engineering without the power and flexibility
of computers and numerical methods. L T
T(r, z) 0
r0
T0 4 Complications
r Analytical solution:
T(r, z) – T
J0(λnr) sinh λn(L – z)
————– = ∑ ———— –—————
λnJ1( λnr0) sinh (λnL)
T0 – T
n=1
where λn’s are roots of J0(λnr0) = 0 FIGURE 5–4
Some analytical solutions are very
complex and difficult to use. Some problems can be solved analytically, but the solution procedure is so
complex and the resulting solution expressions so complicated that it is not
worth all that effort. With the exception of steady onedimensional or transient
lumped system problems, all heat conduction problems result in partial
differential equations. Solving such equations usually requires mathematical
sophistication beyond that acquired at the undergraduate level, such as orthogonality, eigenvalues, Fourier and Laplace transforms, Bessel and Legendre
functions, and infinite series. In such cases, the evaluation of the solution,
which often involves double or triple summations of infinite series at a specified point, is a challenge in itself (Fig. 5–4). Therefore, even when the solutions are available in some handbooks, they are intimidating enough to scare
prospective users away. 5 Human Nature FIGURE 5–5
The ready availability of highpowered computers with sophisticated
software packages has made
numerical solution the norm
rather than the exception. As human beings, we like to sit back and make wishes, and we like our wishes
to come true without much effort. The invention of TV remote controls made
us feel like kings in our homes since the commands we give in our comfortable chairs by pressing buttons are immediately carried out by the obedient
TV sets. After all, what good is cable TV without a remote control. We certainly would love to continue being the king in our little cubicle in the engineering office by solving problems at the press of a button on a computer
(until they invent a remote control for the computers, of course). Well, this
might have been a fantasy yesterday, but it is a reality today. Practically all
engineering offices today are equipped with highpowered computers with
sophisticated software packages, with impressive presentationstyle colorful
output in graphical and tabular form (Fig. 5–5). Besides, the results are as
accurate as the analytical results for all practical purposes. The computers
have certainly changed the way engineering is practiced.
The discussions above should not lead you to believe that analytical solutions are unnecessary and that they should be discarded from the engineering
curriculum. On the contrary, insight to the physical phenomena and engineering wisdom is gained primarily through analysis. The “feel” that engineers
develop during the analysis of simple but fundamental problems serves as
an invaluable tool when interpreting a huge pile of results obtained from a
computer when solving a complex problem. A simple analysis by hand for
a limiting case can be used to check if the results are in the proper range. Also, cen58933_ch05.qxd 9/4/2002 11:41 AM Page 269 269
CHAPTER 5 nothing can take the place of getting “ball park” results on a piece of paper
during preliminary discussions. The calculators made the basic arithmetic
operations by hand a thing of the past, but they did not eliminate the need for
instructing grade school children how to add or multiply.
In this chapter, you will learn how to formulate and solve heat transfer
problems numerically using one or more approaches. In your professional life,
you will probably solve the heat transfer problems you come across using a
professional software package, and you are highly unlikely to write your own
programs to solve such problems. (Besides, people will be highly skeptical
of the results obtained using your own program instead of using a wellestablished commercial software package that has stood the test of time.) The
insight you will gain in this chapter by formulating and solving some heat
transfer problems will help you better understand the available software packages and be an informed and responsible user. 5–2 I FINITE DIFFERENCE FORMULATION
OF DIFFERENTIAL EQUATIONS The numerical methods for solving differential equations are based on replacing the differential equations by algebraic equations. In the case of the popular finite difference method, this is done by replacing the derivatives by
differences. Below we will demonstrate this with both first and secondorder
derivatives. But first we give a motivational example.
Consider a man who deposits his money in the amount of A0 $100 in a
savings account at an annual interest rate of 18 percent, and let us try to determine the amount of money he will have after one year if interest is compounded continuously (or instantaneously). In the case of simple interest, the
money will earn $18 interest, and the man will have 100 100 0.18
$118.00 in his account after one year. But in the case of compounding, the
interest earned during a compounding period will also earn interest for the
remaining part of the year, and the yearend balance will be greater than $118.
For example, if the money is compounded twice a year, the balance will be
100 100 (0.18/2) $109 after six months, and 109 109 (0.18/2)
$118.81 at the end of the year. We could also determine the balance A directly from
A A0(1 i)n ($100)(1 0.09)2 $118.81 (54) where i is the interest rate for the compounding period and n is the number of
periods. Using the same formula, the yearend balance is determined for
monthly, daily, hourly, minutely, and even secondly compounding, and the results are given in Table 5–1.
Note that in the case of daily compounding, the yearend balance will be
$119.72, which is $1.72 more than the simple interest case. (So it is no wonder
that the credit card companies usually charge interest compounded daily when
determining the balance.) Also note that compounding at smaller time intervals, even at the end of each second, does not change the result, and we suspect that instantaneous compounding using “differential” time intervals dt will
give the same result. This suspicion is confirmed by obtaining the differential TABLE 5–1
Yearend balance of a $100 account
earning interest at an annual rate
of 18 percent for various
compounding periods
Compounding
Period Number
of
YearEnd
Periods, n Balance 1 year
1
$118.00
6 months
2
118.81
1 month
12
119.56
1 week
52
119.68
1 day
365
119.72
1 hour
8760
119.72
1 minute
525,600
119.72
1 second
31,536,000 119.72
Instantaneous
119.72 cen58933_ch05.qxd 9/4/2002 11:41 AM Page 270 270
HEAT TRANSFER equation dA/dt
stitution yields iA for the balance A, whose solution is A
A ($100)exp(0.18 1) A0 exp(it). Sub $119.72 which is identical to the result for daily compounding. Therefore, replacing a
differential time interval dt by a finite time interval of t 1 day gave the
same result, which leads us into believing that reasonably accurate results can
be obtained by replacing differential quantities by sufficiently small differences. Next, we develop the finite difference formulation of heat conduction
problems by replacing the derivatives in the differential equations by differences. In the following section we will do it using the energy balance method,
which does not require any knowledge of differential equations.
Derivatives are the building blocks of differential equations, and thus we
first give a brief review of derivatives. Consider a function f that depends on
x, as shown in Figure 5–6. The first derivative of f(x) at a point is equivalent
to the slope of a line tangent to the curve at that point and is defined as f (x) df(x)
dx f (x + ∆ x) lim x→0 f
x lim f(x x→0 x)
x f(x)
(55) ∆f
f (x) which is the ratio of the increment f of the function to the increment x of the
independent variable as x → 0. If we don’t take the indicated limit, we will
have the following approximate relation for the derivative: ∆x Tangent line
x x + ∆x FIGURE 5–6
The derivative of a function at a point
represents the slope of the function
at that point. Plane wall Tm + 1
Tm
Tm – 1 01 2 L
m m +1… M
M –1
m– 1 m+ 1
–
–
2
2 … m –1 x FIGURE 5–7
Schematic of the nodes and the nodal
temperatures used in the development
of the finite difference formulation
of heat transfer in a plane wall. f(x x)
x f(x)
(56) This approximate expression of the derivative in terms of differences is the
finite difference form of the first derivative. The equation above can also be
obtained by writing the Taylor series expansion of the function f about the
point x,
f(x T (x) 0 df(x)
dx x x) f(x) x df(x)
dx 1 2 d 2f(x)
x
2
dx2 ··· (57) and neglecting all the terms in the expansion except the first two. The first
term neglected is proportional to x2, and thus the error involved in each step
of this approximation is also proportional to x2. However, the commutative
error involved after M steps in the direction of length L is proportional to x
since M x2 (L/ x) x2 L x. Therefore, the smaller the x, the smaller
the error, and thus the more accurate the approximation.
Now consider steady onedimensional heat transfer in a plane wall of thickness L with heat generation. The wall is subdivided into M sections of equal
thickness x L/M in the xdirection, separated by planes passing through
M 1 points 0, 1, 2, . . . , m 1, m, m 1, . . . , M called nodes or nodal
points, as shown in Figure 5–7. The xcoordinate of any point m is simply
xm m x, and the temperature at that point is simply T(xm) Tm.
The heat conduction equation involves the second derivatives of temperature with respect to the space variables, such as d 2T/dx2, and the finite difference formulation is based on replacing the second derivatives by appropriate cen58933_ch05.qxd 9/4/2002 11:41 AM Page 271 271
CHAPTER 5 differences. But we need to start the process with first derivatives. Using
1
Eq. 5–6, the first derivative of temperature dT/dx at the midpoints m 2 and
1
m 2 of the sections surrounding the node m can be expressed as
dT
dx Tm Tm
x 1
2 m 1 dT
dx and Tm
m 1
2 Tm 1 (58) x Noting that the second derivative is simply the derivative of the first derivative, the second derivative of temperature at node m can be expressed as
dT
dx 2 dT
dx2 m dT
dx 1
2 Tm 1
2 m x Tm 1 Tm
x 1 x x m Tm Tm 1 2Tm
x2 Tm 1 (59) which is the finite difference representation of the second derivative at a general internal node m. Note that the second derivative of temperature at a node
m is expressed in terms of the temperatures at node m and its two neighboring
nodes. Then the differential equation
·
g
k d 2T
dx 2 0 (510) which is the governing equation for steady onedimensional heat transfer in a
plane wall with heat generation and constant thermal conductivity, can be expressed in the finite difference form as (Fig. 5–8)
Tm 1 2T m
x2 Tm 1 ·
gm
k 0, m 1, 2, 3, . . . , M 1 (511) ·
where gm is the rate of heat generation per unit volume at node m. If the surface temperatures T0 and TM are specified, the application of this equation to
each of the M 1 interior nodes results in M 1 equations for the determination of M 1 unknown temperatures at the interior nodes. Solving these
equations simultaneously gives the temperature values at the nodes. If the
temperatures at the outer surfaces are not known, then we need to obtain two
more equations in a similar manner using the specified boundary conditions.
Then the unknown temperatures at M 1 nodes are determined by solving
the resulting system of M 1 equations in M 1 unknowns simultaneously.
Note that the boundary conditions have no effect on the finite difference
formulation of interior nodes of the medium. This is not surprising since the
control volume used in the development of the formulation does not involve
any part of the boundary. You may recall that the boundary conditions had no
effect on the differential equation of heat conduction in the medium either.
The finite difference formulation above can easily be extended to two or
threedimensional heat transfer problems by replacing each second derivative
by a difference equation in that direction. For example, the finite difference
formulation for steady twodimensional heat conduction in a region with Plane wall
Differential equation:
·
g
d 2T —
—– +
=0
2
k
dx
Valid at every point Finite difference equation:
·
Tm – 1 – 2Tm + Tm + 1 gm
—–—————— + — = 0
—
k
∆x2
Valid at discrete points
∆x FIGURE 5–8
The differential equation is valid at
every point of a medium, whereas the
finite difference equation is valid at
discrete points (the nodes) only. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 272 272
HEAT TRANSFER heat generation and constant thermal conductivity can be expressed in rectangular coordinates as (Fig. 5–9)
Tm
m, n + 1 n+1 ∆y m – 1, n
∆y n
n–1 m, n m–1 m m+1 5–3 ·
Qcond, left Volume
element
of node m ·
gm ·
Qcond, right
A general
interior node 0 01 2 1, n Tm, n 1 2Tm, n
y2 Tm, n 1 ·
gm, n
k 0 (512) for m 1, 2, 3, . . . , M 1 and n 1, 2, 3, . . . , N 1 at any interior node
(m, n). Note that a rectangular region that is divided into M equal subregions
in the xdirection and N equal subregions in the ydirection has a total of
(M 1)(N 1) nodes, and Eq. 5–12 can be used to obtain the finite difference equations at (M 1)(N 1) of these nodes (i.e., all nodes except those
at the boundaries).
The finite difference formulation is given above to demonstrate how difference equations are obtained from differential equations. However, we will use
the energy balance approach in the following sections to obtain the numerical
formulation because it is more intuitive and can handle boundary conditions
more easily. Besides, the energy balance approach does not require having the
differential equation before the analysis. FIGURE 5–9
Finite difference mesh for twodimensional conduction in
rectangular coordinates. Plane wall Tm m + 1, n ∆x ∆x
x 2Tm, n
x2 m, n – 1 y 1, n m–1 m m+1
∆x ∆x
∆x FIGURE 5–10
The nodal points and volume
elements for the finite difference
formulation of onedimensional
conduction in a plane wall. L
M x I ONEDIMENSIONAL
STEADY HEAT CONDUCTION In this section we will develop the finite difference formulation of heat conduction in a plane wall using the energy balance approach and discuss how to
solve the resulting equations. The energy balance method is based on subdividing the medium into a sufficient number of volume elements and then
applying an energy balance on each element. This is done by first selecting
the nodal points (or nodes) at which the temperatures are to be determined and
then forming elements (or control volumes) over the nodes by drawing lines
through the midpoints between the nodes. This way, the interior nodes remain
at the middle of the elements, and the properties at the node such as the
temperature and the rate of heat generation represent the average properties of
the element. Sometimes it is convenient to think of temperature as varying
linearly between the nodes, especially when expressing heat conduction between the elements using Fourier ’s law.
To demonstrate the approach, again consider steady onedimensional heat
·
transfer in a plane wall of thickness L with heat generation g(x) and constant
conductivity k. The wall is now subdivided into M equal regions of thickness
x
L/M in the xdirection, and the divisions between the regions are
selected as the nodes. Therefore, we have M 1 nodes labeled 0, 1, 2, . . . ,
m 1, m, m 1, . . . , M, as shown in Figure 5–10. The xcoordinate of any
node m is simply xm m x, and the temperature at that point is T(xm) Tm.
Elements are formed by drawing vertical lines through the midpoints between
the nodes. Note that all interior elements represented by interior nodes are
fullsize elements (they have a thickness of x), whereas the two elements at
the boundaries are halfsized.
To obtain a general difference equation for the interior nodes, consider the
element represented by node m and the two neighboring nodes m
1 and
m 1. Assuming the heat conduction to be into the element on all surfaces,
an energy balance on the element can be expressed as cen58933_ch05.qxd 9/4/2002 11:41 AM Page 273 273
CHAPTER 5 Rate of heat
conduction
at the left
surface Rate of heat
conduction
at the right
surface Rate of heat
generation
inside the
element Rate of change
of the energy
content of
the element or
·
Q cond, left ·
Q cond, right Eelement
t ·
Gelement 0 (513) since the energy content of a medium (or any part of it) does not change under
steady conditions and thus Eelement 0. The rate of heat generation within
the element can be expressed as
·
Gelement ·
gmVelement ·
gm A x (514) ·
where gm is the rate of heat generation per unit volume in W/m3 evaluated at
node m and treated as a constant for the entire element, and A is heat transfer
area, which is simply the inner (or outer) surface area of the wall.
Recall that when temperature varies linearly, the steady rate of heat conduction across a plane wall of thickness L can be expressed as
·
Q cond kA T
L (515) where T is the temperature change across the wall and the direction of heat
transfer is from the high temperature side to the low temperature. In the case
of a plane wall with heat generation, the variation of temperature is not linear
and thus the relation above is not applicable. However, the variation of temperature between the nodes can be approximated as being linear in the determination of heat conduction across a thin layer of thickness x between two
nodes (Fig. 5–11). Obviously the smaller the distance x between two nodes,
the more accurate is this approximation. (In fact, such approximations are the
reason for classifying the numerical methods as approximate solution methods. In the limiting case of x approaching zero, the formulation becomes exact and we obtain a differential equation.) Noting that the direction of heat
transfer on both surfaces of the element is assumed to be toward the node m,
the rate of heat conduction at the left and right surfaces can be expressed as
·
Q cond, left kA Tm Tm 1 ·
Q cond, right and x kA Tm Tm 1 x (516) Tm Tm 1 x kA Tm Tm 1 x ·
gm A x 0 (517) which simplifies to
Tm 1 2Tm
x2 Tm 1 ·
gm
k 0, m 1, 2, 3, . . . , M 1 (518) Volume
element k Tm + 1 Tm Linear Linear
∆x
m–1 Tm – 1 – Tm
k A ————–
∆x
A Substituting Eqs. 5–14 and 5–16 into Eq. 5–13 gives
kA Tm – 1 ∆x
m m+1 Tm + 1 – Tm
k A ————–
∆x
A FIGURE 5–11
In finite difference formulation, the
temperature is assumed to vary
linearly between the nodes. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 274 274
HEAT TRANSFER
·
g2 A∆ x
T1 – T2
k A ———
∆x T2 – T3
k A ———
∆x
1 2 3
Volume
element
of node 2 T2 – T3 ·
T1 – T2
k A ——— – k A ——— + g2 A∆ x = 0
∆x
∆x
or
·
T1 – 2T2 + T3 + g2 A∆ x 2 / k = 0
(a) Assuming heat transfer to be out of the
volume element at the right surface. ·
g2 A∆ x
T1 – T2
k A ———
∆x T3 – T2
k A ———
∆x
1 2 which is identical to the difference equation (Eq. 5–11) obtained earlier.
Again, this equation is applicable to each of the M 1 interior nodes, and its
application gives M 1 equations for the determination of temperatures at
M 1 nodes. The two additional equations needed to solve for the M 1 unknown nodal temperatures are obtained by applying the energy balance on the
two elements at the boundaries (unless, of course, the boundary temperatures
are specified).
You are probably thinking that if heat is conducted into the element from
both sides, as assumed in the formulation, the temperature of the medium will
have to rise and thus heat conduction cannot be steady. Perhaps a more realistic approach would be to assume the heat conduction to be into the element on
the left side and out of the element on the right side. If you repeat the formulation using this assumption, you will again obtain the same result since the
heat conduction term on the right side in this case will involve Tm Tm 1 instead of Tm 1 Tm, which is subtracted instead of being added. Therefore,
the assumed direction of heat conduction at the surfaces of the volume elements has no effect on the formulation, as shown in Figure 5–12. (Besides, the
actual direction of heat transfer is usually not known.) However, it is convenient to assume heat conduction to be into the element at all surfaces and not
worry about the sign of the conduction terms. Then all temperature differences
in conduction relations are expressed as the temperature of the neighboring
node minus the temperature of the node under consideration, and all conduction terms are added. 3
Volume
element
of node 2 T3 – T2 ·
T1 – T2
k A ——— + k A ——— + g2 A∆ x = 0
∆x
∆x
or
·
T1 – 2T2 + T3 + g2 A∆ x 2 / k = 0
(b) Assuming heat transfer to be into the
volume element at all surfaces. FIGURE 5–12
The assumed direction of heat transfer
at surfaces of a volume element has
no effect on the finite difference
formulation. Boundary Conditions
Above we have developed a general relation for obtaining the finite difference
equation for each interior node of a plane wall. This relation is not applicable
to the nodes on the boundaries, however, since it requires the presence of
nodes on both sides of the node under consideration, and a boundary node
does not have a neighboring node on at least one side. Therefore, we need to
obtain the finite difference equations of boundary nodes separately. This is
best done by applying an energy balance on the volume elements of boundary
nodes.
Boundary conditions most commonly encountered in practice are the specified temperature, specified heat flux, convection, and radiation boundary
conditions, and here we develop the finite difference formulations for them
for the case of steady onedimensional heat conduction in a plane wall of
thickness L as an example. The node number at the left surface at x 0 is 0,
and at the right surface at x L it is M. Note that the width of the volume element for either boundary node is x/2.
The specified temperature boundary condition is the simplest boundary
condition to deal with. For onedimensional heat transfer through a plane wall
of thickness L, the specified temperature boundary conditions on both the left
and right surfaces can be expressed as (Fig. 5–13)
T(0)
T(L) T0
TM Specified value
Specified value (519) where T0 and Tm are the specified temperatures at surfaces at x 0 and x L,
respectively. Therefore, the specified temperature boundary conditions are cen58933_ch05.qxd 9/4/2002 11:41 AM Page 275 275
CHAPTER 5 incorporated by simply assigning the given surface temperatures to the boundary nodes. We do not need to write an energy balance in this case unless we
decide to determine the rate of heat transfer into or out of the medium after the
temperatures at the interior nodes are determined.
When other boundary conditions such as the specified heat flux, convection,
radiation, or combined convection and radiation conditions are specified at a
boundary, the finite difference equation for the node at that boundary is obtained by writing an energy balance on the volume element at that boundary.
The energy balance is again expressed as
·
Gelement ·
Q 0 Plane wall
35°C 0 82°C 0 1 for heat transfer under steady conditions. Again we assume all heat transfer to
be into the volume element from all surfaces for convenience in formulation,
except for specified heat flux since its direction is already specified. Specified
heat flux is taken to be a positive quantity if into the medium and a negative
quantity if out of the medium. Then the finite difference formulation at the
node m 0 (at the left boundary where x 0) of a plane wall of thickness L
during steady onedimensional heat conduction can be expressed as (Fig.
5–14)
kA T0 T1 ·
g0(A x/2) x 0 T1 T0 ·
g0(A x/2) x ·
Special case: Insulated Boundary (q0
kA T1 T0 0 (522) 0) ·
g0(A x/2) x 0 (523) 2. Convection Boundary Condition
hA(T T0) kA T1 T0
x ·
g0(A x/2) FIGURE 5–13
Finite difference formulation of
specified temperature boundary
conditions on both surfaces
of a plane wall. ∆x
—
–
2
·
g0 1. Specified Heat Flux Boundary Condition
kA M T0 = 35°C
TM = 82°C (521) where A x/2 is the volume of the volume element (note that the boundary ele·
ment has half thickness), g0 is the rate of heat generation per unit volume (in
3
W/m ) at x 0, and A is the heat transfer area, which is constant for a plane
wall. Note that we have x in the denominator of the second term instead of
x/2. This is because the ratio in that term involves the temperature difference
between nodes 0 and 1, and thus we must use the distance between those two
nodes, which is x.
The finite difference form of various boundary conditions can be obtained
·
from Eq. 5–21 by replacing Q left surface by a suitable expression. Next this is
done for various boundary conditions at the left boundary. ·
q 0A L … (520) all sides ·
Q left surface 2 0 (524) Volume element
of node 0
T1 – T0
k A ———
∆x ·
Qleft surface 0 0 1
∆x 2 … L
x ∆x T1 – T0 · ∆ x
·
–
Qleft surface + k A ——— + g0 A — = 0
∆x
2 FIGURE 5–14
Schematic for the finite difference
formulation of the left boundary
node of a plane wall. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 276 276
HEAT TRANSFER 3. Radiation Boundary Condition
Tsurr
4
A(Tsurr ∆x
—
–
2 ε T1 – T0
k A ———
∆x
hA(T – T0)
0 0 1 2 ∆x A … hA(T
L
x T0 ·
g0(A x/2) x 0 (525) ·
q 0A Medium B
kB kAA Tm + 1 – Tm
kB A ————–
∆x m–1 m
∆x x m+1
∆x ∆x ∆x
—– —–
22 kA T1 T0
x ·
g0(A x/2) 0 (526) T0) kA T1 T0
x ·
g0(A x/2) 0 (527) hA(T T0) 4
A(Tsurr T04) kA T1 T0
x ·
g0(A x/2) 0 (528) 6. Interface Boundary Condition Two different solid media A and B
are assumed to be in perfect contact, and thus at the same temperature
at the interface at node m (Fig. 5–16). Subscripts A and B indicate
properties of media A and B, respectively. Interface Tm – 1 – Tm
kA A ————–
∆x T04) 5. Combined Convection, Radiation, and Heat Flux Boundary
Condition FIGURE 5–15
Schematic for the finite difference
formulation of combined convection
and radiation on the left boundary
of a plane wall. ·
·
gA,m gB,m 4
A(Tsurr T0) hcombined A(T T1 – T0 · ∆ x
+ kA ——— + g0 A —– = 0
∆x
2 A T1 or ∆x hA(T – T0) + εσ A(T 4 – T 4 )
surr
0 Medium A
kA kA 4. Combined Convection and Radiation Boundary Condition
(Fig. 5–15) ·
g0 εσ A(T 4 – T 4 )
surr
0 T04) A Tm + 1 – Tm
Tm – 1 – Tm
kA A ————– + kB A ————–
∆x
∆x
∆x ·
∆x
· A —– + g A —– = 0
+ gA, m
B, m
2
2 FIGURE 5–16
Schematic for the finite difference
formulation of the interface boundary
condition for two mediums A and B
that are in perfect thermal contact. Tm Tm 1 x kBA Tm Tm 1 x ·
gA, m(A x/2) ·
gB, m(A x/2) 0 (529) ·
In these relations, q0 is the specified heat flux in W/m2, h is the convection
coefficient, hcombined is the combined convection and radiation coefficient, T is
the temperature of the surrounding medium, Tsurr is the temperature of the
surrounding surfaces, is the emissivity of the surface, and is the Stefan–
Boltzman constant. The relations above can also be used for node M on the
right boundary by replacing the subscript “0” by “M” and the subscript “1” by
“M 1”.
Note that absolute temperatures must be used in radiation heat transfer
calculations, and all temperatures should be expressed in K or R when a
boundary condition involves radiation to avoid mistakes. We usually try to
avoid the radiation boundary condition even in numerical solutions since it
causes the finite difference equations to be nonlinear, which are more difficult
to solve. Treating Insulated Boundary Nodes as Interior Nodes:
The Mirror Image Concept
One way of obtaining the finite difference formulation of a node on an insulated boundary is to treat insulation as “zero” heat flux and to write an energy
balance, as done in Eq. 5–23. Another and more practical way is to treat the
node on an insulated boundary as an interior node. Conceptually this is done cen58933_ch05.qxd 9/4/2002 11:41 AM Page 277 277
CHAPTER 5 by replacing the insulation on the boundary by a mirror and considering the
reflection of the medium as its extension (Fig. 5–17). This way the node next
to the boundary node appears on both sides of the boundary node because of
symmetry, converting it into an interior node. Then using the general formula
(Eq. 5–18) for an interior node, which involves the sum of the temperatures of
the adjoining nodes minus twice the node temperature, the finite difference
formulation of a node m 0 on an insulated boundary of a plane wall can be
expressed as
Tm 1 2Tm
x2 Tm 1 ·
gm
k 0 → T1 2T0
x2 T1 ·
g0
k 0 EXAMPLE 5–1 01 x 2 Mirror Equivalent
interior
node Mirror
image (530) which is equivalent to Eq. 5–23 obtained by the energy balance approach.
The mirror image approach can also be used for problems that possess thermal symmetry by replacing the plane of symmetry by a mirror. Alternately, we
can replace the plane of symmetry by insulation and consider only half of the
medium in the solution. The solution in the other half of the medium is simply the mirror image of the solution obtained. Insulated
boundary
node Insulation x 2 1 01 x 2 FIGURE 5–17
A node on an insulated boundary
can be treated as an interior node by
replacing the insulation by a mirror. Steady Heat Conduction in a Large Uranium Plate Consider a large uranium plate of thickness L 4 cm and thermal conductivity
k
28 W/m · °C in which heat is generated uniformly at a constant rate of
·
g
5 106 W/m3. One side of the plate is maintained at 0°C by iced water
while the other side is subjected to convection to an environment at T
30°C
with a heat transfer coefficient of h 45 W/m2 · °C, as shown in Figure 5–18.
Considering a total of three equally spaced nodes in the medium, two at the
boundaries and one at the middle, estimate the exposed surface temperature of
the plate under steady conditions using the finite difference approach. SOLUTION A uranium plate is subjected to specified temperature on one side
and convection on the other. The unknown surface temperature of the plate is
to be determined numerically using three equally spaced nodes.
Assumptions 1 Heat transfer through the wall is steady since there is no indication of any change with time. 2 Heat transfer is onedimensional since
the plate is large relative to its thickness. 3 Thermal conductivity is constant.
4 Radiation heat transfer is negligible.
Properties The thermal conductivity is given to be k 28 W/m · °C.
Analysis The number of nodes is specified to be M 3, and they are chosen
to be at the two surfaces of the plate and the midpoint, as shown in the figure.
Then the nodal spacing x becomes x L
M 1 0.04 m
31 0.02 m We number the nodes 0, 1, and 2. The temperature at node 0 is given to be
T0
0°C, and the temperatures at nodes 1 and 2 are to be determined. This
problem involves only two unknown nodal temperatures, and thus we need to
have only two equations to determine them uniquely. These equations are obtained by applying the finite difference method to nodes 1 and 2. Uranium
plate
0°C 0 h
T k = 28 W/m·°C
·
g = 5 × 106 W/m3
L
0 1 2 x FIGURE 5–18
Schematic for Example 5–1. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 278 278
HEAT TRANSFER Node 1 is an interior node, and the finite difference formulation at that node
is obtained directly from Eq. 5–18 by setting m 1: T0 2T1
x2 T2 ·
g1
k 0 → 0 2T1
x2 T2 ·
g1
k → 2T1 0 T2 ·
g1 x2
k
(1) Node 2 is a boundary node subjected to convection, and the finite difference
formulation at that node is obtained by writing an energy balance on the volume
element of thickness x/2 at that boundary by assuming heat transfer to be into
the medium at all sides: hA(T T2) kA T1 T2
x ·
g2(A x/2) 0 Canceling the heat transfer area A and rearranging give 1 T1 hx
T2
k ·
g2 x2
2k hx
T
k (2) Equations (1) and (2) form a system of two equations in two unknowns T1 and
T2. Substituting the given quantities and simplifying gives T1 2T1 T2
1.032T2 71.43
36.68 (in °C)
(in °C) This is a system of two algebraic equations in two unknowns and can be solved
easily by the elimination method. Solving the first equation for T1 and substituting into the second equation result in an equation in T2 whose solution is T2 136.1°C This is the temperature of the surface exposed to convection, which is the
desired result. Substitution of this result into the first equation gives T1
103.8°C, which is the temperature at the middle of the plate. h
T Plate 0 1
2 cm 2 x 2 cm Finite difference solution:
T2 = 136.1°C
Exact solution:
T2 = 136.0°C FIGURE 5–19
Despite being approximate in nature,
highly accurate results can be
obtained by numerical methods. Discussion The purpose of this example is to demonstrate the use of the finite
difference method with minimal calculations, and the accuracy of the result
was not a major concern. But you might still be wondering how accurate the result obtained above is. After all, we used a mesh of only three nodes for the
entire plate, which seems to be rather crude. This problem can be solved analytically as described in Chapter 2, and the analytical (exact) solution can be
shown to be T(x) ·
0.5ghL2/k
hL ·
gL
k Th x ·
gx2
2k Substituting the given quantities, the temperature of the exposed surface of the
plate at x L 0.04 m is determined to be 136.0°C, which is almost identical to the result obtained here with the approximate finite difference method
(Fig. 5–19). Therefore, highly accurate results can be obtained with numerical
methods by using a limited number of nodes. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 279 279
CHAPTER 5 EXAMPLE 5–2 Heat Transfer from Triangular Fins Consider an aluminum alloy fin (k 180 W/m · °C) of triangular cross section
with length L 5 cm, base thickness b 1 cm, and very large width w in the
direction normal to the plane of paper, as shown in Figure 5–20. The base of
the fin is maintained at a temperature of T0
200°C. The fin is losing heat
to the surrounding medium at T
25°C with a heat transfer coefficient of
h
15 W/m2 · °C. Using the finite difference method with six equally spaced
nodes along the fin in the xdirection, determine (a) the temperatures at the
nodes, (b) the rate of heat transfer from the fin for w
1 m, and (c) the fin
efficiency. h, T T0
b
0 SOLUTION A long triangular fin attached to a surface is considered. The nodal
temperatures, the rate of heat transfer, and the fin efficiency are to be determined numerically using six equally spaced nodes. Triangular fin w 1 2 0.05 m
61 L
M 1 0 → kAleft Tm Tm 1 kAright x all sides Tm 0.01 m Tm 1 x hAconv(T Tm) 0 Note that heat transfer areas are different for each node in this case, and using
geometrical relations, they can be expressed as Aleft
Aright
Aconv (Height Width)@m
(Height Width)@m
2 Length Width 1
2
1
2 2w[L (m
2w[L (m
2w( x/cos ) 1/2) x]tan
1/2) x]tan Substituting, 2kw[L
2kw[L 1)
2 (m
(m x]tan
1)
2 Tm x]tan Tm 1 x
Tm Tm 1 x h 5 2w x
(T
cos Tm) x L ∆x
——–
cos θ [L – (m + 1 )∆ x]tan θ
–
2 θ The temperature at node 0 is given to be T0 200°C, and the temperatures at
the remaining five nodes are to be determined. Therefore, we need to have five
equations to determine them uniquely. Nodes 1, 2, 3, and 4 are interior nodes,
and the finite difference formulation for a general interior node m is obtained
by applying an energy balance on the volume element of this node. Noting that
heat transfer is steady and there is no heat generation in the fin and assuming
heat transfer to be into the medium at all sides, the energy balance can be expressed as ·
Q 4 180 W/m · °C. Analysis (a) The number of nodes in the fin is specified to be M 6, and their
location is as shown in the figure. Then the nodal spacing x becomes x 3 ∆x Assumptions 1 Heat transfer is steady since there is no indication of any
change with time. 2 The temperature along the fin varies in the x direction only.
3 Thermal conductivity is constant. 4 Radiation heat transfer is negligible.
Properties The thermal conductivity is given to be k b/ 2
tan θ = —–
L θ 0 m–1 m m+1 ∆x
L – (m – 1 )∆ x
–
2
1 )∆ x]tan θ
[L – (m – –
2 FIGURE 5–20
Schematic for Example 5–2 and the
volume element of a general
interior node of the fin. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 280 280
HEAT TRANSFER Dividing each term by 2kwL tan / x gives 1 x
(Tm
L 1
)
2 (m 1 Tm) 1 (m 1
)
2 x
(Tm
L 1 h( x)2
(T
kL sin Tm)
Tm) 0 Note that b/2
L tan
Also, sin 5.71° (5.5 0.5 cm
5 cm 0.1 → tan 10.1 5.71° 0.0995. Then the substitution of known quantities gives m)Tm 1 (10.00838 2m)Tm (4.5 m)Tm 1 0.209 Now substituting 1, 2, 3, and 4 for m results in these finite difference equations for the interior nodes: m
m
m
m kAleft T4 5 ∆x
—–
2 FIGURE 5–21
Schematic of the volume element of
node 5 at the tip of a triangular fin. T5
x θ 4 ∆x
—–
2 8.00838T1 3.5T2
900.209
3.5T1 6.00838T2 2.5T3
0.209
2.5T2 4.00838T3 1.5T4
0.209
1.5T3 2.00838T4 0.5T5
0.209 (1)
(2)
(3)
(4) The finite difference equation for the boundary node 5 is obtained by writing an
energy balance on the volume element of length x/2 at that boundary, again by
assuming heat transfer to be into the medium at all sides (Fig. 5–21): ∆x/2
—–—
cos θ
∆x
—– tan θ
2 1:
2:
3:
4: hAconv (T T5) 0 where Aleft 2w x
tan
2 and Aconv 2w x/2
cos Canceling w in all terms and substituting the known quantities gives T4 1.00838T5 0.209 (5) Equations (1) through (5) form a linear system of five algebraic equations in five
unknowns. Solving them simultaneously using an equation solver gives T1
T4 198.6°C,
194.3°C, T2
T5 197.1°C,
192.9°C T3 195.7°C, which is the desired solution for the nodal temperatures.
(b) The total rate of heat transfer from the fin is simply the sum of the heat
transfer from each volume element to the ambient, and for w 1 m it is determined from cen58933_ch05.qxd 9/4/2002 11:41 AM Page 281 281
CHAPTER 5
5 ·
Q fin 5 ·
Q element, m m0 hAconv, m(Tm T) m0 Noting that the heat transfer surface area is w x/cos for the boundary nodes
0 and 5, and twice as large for the interior nodes 1, 2, 3, and 4, we have ·
Q fin h h wx
[(T0 T ) 2(T1 T )
cos
2(T4 T ) (T5 T )]
wx
[T0
cos 2(T1 (15 W/m2 · °C) T2 T3 2(T2 T4) (1 m)(0.01 m)
[200
cos 5.71° T5
2 T) 2(T3 T) 10T ]
785.7 192.9 10 25] 258.4 W
(c) If the entire fin were at the base temperature of T0
of heat transfer from the fin for w 1 m would be ·
Q max 200°C, the total rate hAfin, total (T0 T ) h(2wL/cos )(T0 T )
(15 W/m2 · °C)[2(1 m)(0.05 m)/cos5.71°](200
263.8 W 25)°C Then the fin efficiency is determined from
fin ·
Q fin
·
Q max 258.4 W
263.8 W 0.98 which is less than 1, as expected. We could also determine the fin efficiency in
this case from the proper fin efficiency curve in Chapter 3, which is based on
the analytical solution. We would read 0.98 for the fin efficiency, which is identical to the value determined above numerically. The finite difference formulation of steady heat conduction problems usually results in a system of N algebraic equations in N unknown nodal temperatures that need to be solved simultaneously. When N is small (such as 2 or 3),
we can use the elementary elimination method to eliminate all unknowns except one and then solve for that unknown (see Example 5–1). The other unknowns are then determined by back substitution. When N is large, which is
usually the case, the elimination method is not practical and we need to use a
more systematic approach that can be adapted to computers.
There are numerous systematic approaches available in the literature, and
they are broadly classified as direct and iterative methods. The direct methods are based on a fixed number of welldefined steps that result in the solution in a systematic manner. The iterative methods, on the other hand, are
based on an initial guess for the solution that is refined by iteration until a
specified convergence criterion is satisfied (Fig. 5–22). The direct methods
usually require a large amount of computer memory and computation time, Direct methods:
Solve in a systematic manner following a
series of welldefined steps.
Iterative methods:
Start with an initial guess for the solution,
and iterate until solution converges. FIGURE 5–22
Two general categories of solution
methods for solving systems
of algebraic equations. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 282 282
HEAT TRANSFER and they are more suitable for systems with a relatively small number of equations. The computer memory requirements for iterative methods are minimal,
and thus they are usually preferred for large systems. The convergence of iterative methods to the desired solution, however, may pose a problem. 5–4 y … N ∆y
∆y Node (m, n) … n+1
n
n–1 2
1
0 ∆x ∆x … 012 …
m
m–1 m+1 M x FIGURE 5–23
The nodal network for the finite
difference formulation of twodimensional conduction in
rectangular coordinates. m, n + 1 n+1
∆y Volume
element m – 1, n n ·
gm, n m, n m + 1, n Rate of heat conduction
at the left, top, right,
and bottom surfaces m, n – 1
∆x y m–1 ∆x
m TWODIMENSIONAL
STEADY HEAT CONDUCTION In Section 5–3 we considered onedimensional heat conduction and assumed
heat conduction in other directions to be negligible. Many heat transfer problems encountered in practice can be approximated as being onedimensional,
but this is not always the case. Sometimes we need to consider heat transfer in
other directions as well when the variation of temperature in other directions
is significant. In this section we will consider the numerical formulation and
solution of twodimensional steady heat conduction in rectangular coordinates
using the finite difference method. The approach presented below can be extended to threedimensional cases.
Consider a rectangular region in which heat conduction is significant in the
x and ydirections. Now divide the xy plane of the region into a rectangular
mesh of nodal points spaced x and y apart in the x and ydirections,
respectively, as shown in Figure 5–23, and consider a unit depth of z 1
in the zdirection. Our goal is to determine the temperatures at the nodes,
and it is convenient to number the nodes and describe their position by
the numbers instead of actual coordinates. A logical numbering scheme for
twodimensional problems is the double subscript notation (m, n) where
m 0, 1, 2, . . . , M is the node count in the xdirection and n 0, 1, 2, . . . , N
is the node count in the ydirection. The coordinates of the node (m, n) are
simply x
m x and y
n y, and the temperature at the node (m, n) is
denoted by Tm, n.
Now consider a volume element of size x
y 1 centered about a gen·
eral interior node (m, n) in a region in which heat is generated at a rate of g and
the thermal conductivity k is constant, as shown in Figure 5–24. Again
assuming the direction of heat conduction to be toward the node under
consideration at all surfaces, the energy balance on the volume element can be
expressed as ∆y
n–1 I Rate of heat
generation inside
the element Rate of change of
the energy content
of the element or
m+1 x FIGURE 5–24
The volume element of a general
interior node (m, n) for twodimensional conduction in
rectangular coordinates. ·
Q cond, left ·
Q cond, top ·
Q cond, right ·
Q cond, bottom ·
Gelement Eelement
t 0 (531) for the steady case. Again assuming the temperatures between the adjacent nodes to vary linearly and noting that the heat transfer area is
y1
y in the xdirection and Ay
x1
x in the ydirection,
Ax
the energy balance relation above becomes cen58933_ch05.qxd 9/4/2002 11:41 AM Page 283 283
CHAPTER 5 ky Tm Tm, n 1, n kx x Tm, n kx Dividing each term by x
Tm 1, n 2Tm, n Tm Tm, n 1 y
Tm, n ky
Tm, n 1 y Tm Tm, n 1, n x
·
gm, n xy 0 (532) y and simplifying gives
Tm, n 1, n 2Tm, n 1 2 Tm, n ·
gm, n
k 1 2 x y 0 (533) for m 1, 2, 3, . . . , M 1 and n 1, 2, 3, . . . , N 1. This equation is identical to Eq. 5–12 obtained earlier by replacing the derivatives in the differential equation by differences for an interior node (m, n). Again a rectangular
region M equally spaced nodes in the xdirection and N equally spaced nodes
in the ydirection has a total of (M 1)(N 1) nodes, and Eq. 5–33 can be
used to obtain the finite difference equations at all interior nodes.
In finite difference analysis, usually a square mesh is used for simplicity (except when the magnitudes of temperature gradients in the x and
ydirections are very different), and thus x and y are taken to be the same.
Then x
y l, and the relation above simplifies to
Tm 1, n Tm 1, n Tm, n 1 Tm, n 1 4Tm, n ·
gm, nl 2
k 0 (534) That is, the finite difference formulation of an interior node is obtained by
adding the temperatures of the four nearest neighbors of the node, subtracting
four times the temperature of the node itself, and adding the heat generation
term. It can also be expressed in this form, which is easy to remember:
Tleft Ttop Tright Tbottom 4Tnode ·
g nodel 2
k 0 (535) When there is no heat generation in the medium, the finite difference equation for an interior node further simplifies to Tnode (Tleft Ttop Tright
Tbottom)/4, which has the interesting interpretation that the temperature of each
interior node is the arithmetic average of the temperatures of the four neighboring nodes. This statement is also true for the threedimensional problems
except that the interior nodes in that case will have six neighboring nodes instead of four. Boundary Nodes
The development of finite difference formulation of boundary nodes in two(or three) dimensional problems is similar to the development in the onedimensional case discussed earlier. Again, the region is partitioned between
the nodes by forming volume elements around the nodes, and an energy balance is written for each boundary node. Various boundary conditions can be
handled as discussed for a plane wall, except that the volume elements
in the twodimensional case involve heat transfer in the ydirection as well as
the xdirection. Insulated surfaces can still be viewed as “mirrors, ” and the cen58933_ch05.qxd 9/4/2002 11:41 AM Page 284 284
HEAT TRANSFER
Volume element
of node 2 h, T 1 ·
Qtop 2 3 ·
Qleft ∆y mirror image concept can be used to treat nodes on insulated boundaries as interior nodes.
For heat transfer under steady conditions, the basic equation to keep in mind
when writing an energy balance on a volume element is (Fig. 5–25) Boundary
subjected
to convection ·
Q ·
Qright
·
Qbottom ·
g2V2
·
·
·
·
Qleft + Qtop + Qright + Qbottom + —— = 0
k FIGURE 5–25
The finite difference formulation of
a boundary node is obtained by
writing an energy balance
on its volume element. EXAMPLE 5–3 Convection
h, T y
1 2 3 4 5 6 7 8 10 11 12 13 14 15 ∆ x = ∆y = l ∆y ·
9 qR ∆y ∆x ∆x 90°C
∆x x
∆x h, T h, T
2 1 2 3 5 4
(a) Node 1 (536) Steady TwoDimensional Heat Conduction
in LBars Consider steady heat transfer in an Lshaped solid body whose cross section is
given in Figure 5–26. Heat transfer in the direction normal to the plane of the
paper is negligible, and thus heat transfer in the body is twodimensional. The
thermal conductivity of the body is k 15 W/m · °C, and heat is generated in
·
the body at a rate of g
2 106 W/m3. The left surface of the body is insulated, and the bottom surface is maintained at a uniform temperature of 90°C.
The entire top surface is subjected to convection to ambient air at T
25°C
with a convection coefficient of h 80 W/m2 · °C, and the right surface is sub·
jected to heat flux at a uniform rate of q R 5000 W/m2. The nodal network of
the problem consists of 15 equally spaced nodes with x
y
1.2 cm, as
shown in the figure. Five of the nodes are at the bottom surface, and thus their
temperatures are known. Obtain the finite difference equations at the remaining nine nodes and determine the nodal temperatures by solving them. ∆x FIGURE 5–26
Schematic for Example 5–3 and
the nodal network (the boundaries
of volume elements of the nodes are
indicated by dashed lines). 1 0 whether the problem is one, two, or threedimensional. Again we assume,
for convenience in formulation, all heat transfer to be into the volume element from all surfaces except for specified heat flux, whose direction is already specified. This is demonstrated in Example 5–3 for various boundary
conditions. 4 ∆x ·
gVelement all sides (b) Node 2 FIGURE 5–27
Schematics for energy balances on the
volume elements of nodes 1 and 2. SOLUTION Heat transfer in a long Lshaped solid bar with specified boundary
conditions is considered. The nine unknown nodal temperatures are to be determined with the finite difference method.
Assumptions 1 Heat transfer is steady and twodimensional, as stated. 2 Thermal conductivity is constant. 3 Heat generation is uniform. 4 Radiation heat
transfer is negligible.
Properties The thermal conductivity is given to be k 15 W/m · °C.
Analysis We observe that all nodes are boundary nodes except node 5, which
is an interior node. Therefore, we will have to rely on energy balances to obtain
the finite difference equations. But first we form the volume elements by partitioning the region among the nodes equitably by drawing dashed lines between
the nodes. If we consider the volume element represented by an interior node
to be full size (i.e., x
y 1), then the element represented by a regular
boundary node such as node 2 becomes half size (i.e., x
y/2 1), and
a corner node such as node 1 is quarter size (i.e., x/2
y/2 1). Keeping
Eq. 5–36 in mind for the energy balance, the finite difference equations for
each of the nine nodes are obtained as follows:
(a) Node 1. The volume element of this corner node is insulated on the left and
subjected to convection at the top and to conduction at the right and bottom
surfaces. An energy balance on this element gives [Fig. 5–27a] cen58933_ch05.qxd 9/4/2002 11:41 AM Page 285 285
CHAPTER 5 0 h Taking x x
(T
2
y T1) y T2 T1
2
x k · xy
g1
22 x T4 T1
2
y k 0 l, it simplifies to –2 hl
T1
k T2 hl
T
k T4 ·
g1l 2
2k (b) Node 2. The volume element of this boundary node is subjected to convection at the top and to conduction at the right, bottom, and left surfaces. An
energy balance on this element gives [Fig. 5–27b] h x(T T2) Taking x k y y T3 T2
2
x kx T5 T2
y k y T1 T2
2
x y
·
g2 x
2 0 l, it simplifies to T1 4 2hl
T2
k T3 ·
g2l 2
k 2hl
T
k 2T5 (c) Node 3. The volume element of this corner node is subjected to convection
at the top and right surfaces and to conduction at the bottom and left surfaces.
An energy balance on this element gives [Fig. 5–28a] h x
2 Taking x y
(T
2
y T3) k x T6 T3
2
y k y T2 T3
2
x · xy
g3
22 0 3 (5) 4 5 h, T l, it simplifies to T2 2hl
T3
k 2 2hl
T
k T6 ·
g3l 2
2k (d ) Node 4. This node is on the insulated boundary and can be treated as an
interior node by replacing the insulation by a mirror. This puts a reflected image
of node 5 to the left of node 4. Noting that x
y
l, the general interior
node relation for the steady twodimensional case (Eq. 5–35) gives [Fig. 5–28b] T5
or, noting that T10 2 1 Mirror h, T T1 T5 T10 4T4 ·
g4l 2
k (a) Node 3 0 T2 3 2 4T4 2T5 90 ·
g4l 2
k 4 T6 T11 4T5 ·
g5l 2
k 5 6 h, T
7 6
5 (e) Node 5. This is an interior node, and noting that x
y
l, the finite
difference formulation of this node is obtained directly from Eq. 5–35 to be
[Fig. 5–29a] T4 (b) Node 4 FIGURE 5–28
Schematics for energy balances on the
volume elements of nodes 3 and 4. 90° C, T1 10 6 0 12 11
(a) Node 5 (b) Node 6 FIGURE 5–29
Schematics for energy balances on the
volume elements of nodes 5 and 6. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 286 286
HEAT TRANSFER or, noting that T11 90°C, T2 T4 4T5 T6 90 ·
g5l 2
k (f ) Node 6. The volume element of this inner corner node is subjected to convection at the Lshaped exposed surface and to conduction at other surfaces.
An energy balance on this element gives [Fig. 5–29b] y
(T
2 x
2 h ky
Taking x y T3 6 7 8 T5 T6
x l and noting that T12 2hl
T6
k 6 2T5 9 h x(T ·
qR T7)
k 13 k 90°C, it simplifies to T7 180 ·
3g6l 2
2k 2hl
T
k (g) Node 7. The volume element of this boundary node is subjected to convection at the top and to conduction at the right, bottom, and left surfaces. An energy balance on this element gives [Fig. 5–30a] h, T h, T T12 T6
y T7 T6
kx
2
x
y
x T3 T6
·3xy 0
k
g6
2
4
y T6) 15 FIGURE 5–30
Schematics for energy balances on the
volume elements of nodes 7 and 9. Taking x y k y T6
2 T13 T7
y T8 T7
kx
2
x
y
T7
y
·
g7 x
0
2
x l and noting that T13 T6 2hl
T7
k 4 90°C, it simplifies to T8 180 ·
g7l 2
k 2hl
T
k (h) Node 8. This node is identical to Node 7, and the finite difference formulation of this node can be obtained from that of Node 7 by shifting the node
numbers by 1 (i.e., replacing subscript m by m 1). It gives T7 2hl
T8
k 4 T9 180 ·
g8l 2
k 2hl
T
k (i ) Node 9. The volume element of this corner node is subjected to convection
at the top surface, to heat flux at the right surface, and to conduction at the
bottom and left surfaces. An energy balance on this element gives [Fig. 5–30b] h x
(T
2 Taking x T9)
y ·y
qR
2 k x T15 T9
2
y l and noting that T15 T8 2 hl
T9
k 90 k y T8 T9
2
x · xy
g9
22 90°C, it simplifies to ·
qRl
k hl
T
k ·
g9l 2
2k 0 cen58933_ch05.qxd 9/4/2002 11:41 AM Page 287 287
CHAPTER 5 This completes the development of finite difference formulation for this problem. Substituting the given quantities, the system of nine equations for the
determination of nine unknown nodal temperatures becomes T1 T3 –2.064T1 T2 T4
4.128T2 T3 2T5
T2 2.128T3 T6
T1 4T4 2T5
T2 T4 4T5 T6
2T5 6.128T6 T7
T6 4.128T7 T8
T7 4.128T8 T9
T8 2.064T9 11.2
22.4
12.8
109.2
109.2
212.0
202.4
202.4
105.2 which is a system of nine algebraic equations with nine unknowns. Using an
equation solver, its solution is determined to be T1
T4
T7 112.1°C
109.4°C
97.3°C T2
T5
T8 110.8°C
108.1°C
96.3°C T3
T6
T9 106.6°C
103.2°C
97.6°C Note that the temperature is the highest at node 1 and the lowest at node 8.
This is consistent with our expectations since node 1 is the farthest away from
the bottom surface, which is maintained at 90°C and has one side insulated,
and node 8 has the largest exposed area relative to its volume while being close
to the surface at 90°C. Irregular Boundaries
In problems with simple geometries, we can fill the entire region using simple
volume elements such as strips for a plane wall and rectangular elements for
twodimensional conduction in a rectangular region. We can also use cylindrical or spherical shell elements to cover the cylindrical and spherical bodies
entirely. However, many geometries encountered in practice such as turbine
blades or engine blocks do not have simple shapes, and it is difficult to fill
such geometries having irregular boundaries with simple volume elements.
A practical way of dealing with such geometries is to replace the irregular
geometry by a series of simple volume elements, as shown in Figure 5–31.
This simple approach is often satisfactory for practical purposes, especially
when the nodes are closely spaced near the boundary. More sophisticated approaches are available for handling irregular boundaries, and they are commonly incorporated into the commercial software packages.
EXAMPLE 5–4 Actual boundary Approximation Heat Loss through Chimneys Hot combustion gases of a furnace are flowing through a square chimney made
of concrete (k
1.4 W/m · °C). The flow section of the chimney is 20 cm
20 cm, and the thickness of the wall is 20 cm. The average temperature of the FIGURE 5–31
Approximating an irregular
boundary with a rectangular mesh. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 288 288
HEAT TRANSFER hot gases in the chimney is Ti 300°C, and the average convection heat transfer coefficient inside the chimney is hi 70 W/m2 · °C. The chimney is losing
heat from its outer surface to the ambient air at To 20°C by convection with
a heat transfer coefficient of ho 21 W/m2 · °C and to the sky by radiation. The
emissivity of the outer surface of the wall is
0.9, and the effective sky temperature is estimated to be 260 K. Using the finite difference method with
x
y
10 cm and taking full advantage of symmetry, determine the
temperatures at the nodal points of a cross section and the rate of heat loss for
a 1mlong section of the chimney. SOLUTION Heat transfer through a square chimney is considered. The nodal
temperatures and the rate of heat loss per unit length are to be determined with
the finite difference method.
Assumptions 1 Heat transfer is steady since there is no indication of change
with time. 2 Heat transfer through the chimney is twodimensional since the
height of the chimney is large relative to its cross section, and thus heat conduction through the chimney in the axial direction is negligible. It is tempting
to simplify the problem further by considering heat transfer in each wall to be
onedimensional, which would be the case if the walls were thin and thus the
corner effects were negligible. This assumption cannot be justified in this case
since the walls are very thick and the corner sections constitute a considerable
portion of the chimney structure. 3 Thermal conductivity is constant. Symmetry lines
(Equivalent to insulation) Properties
0.9. h1 T1 1 2 3 4 5 6 7 8 h0, T0
Tsky 9 Representative
section of chimney FIGURE 5–32
Schematic of the chimney discussed in
Example 5–4 and the nodal network
for a representative section.
h, T
1 h, T
2 1 2 (a) Node 1 (a) Node 1. On the inner boundary, subjected to convection, Figure 5–33a hi x
(T
2i T1) k y T2 T1
2
x k x T3 T1
2
y 4
(b) Node 2 FIGURE 5–33
Schematics for energy balances on the
volume elements of nodes 1 and 2. 1.4 W/m · °C and Analysis The cross section of the chimney is given in Figure 5–32. The most
striking aspect of this problem is the apparent symmetry about the horizontal
and vertical lines passing through the midpoint of the chimney as well as the
diagonal axes, as indicated on the figure. Therefore, we need to consider only
oneeighth of the geometry in the solution whose nodal network consists of nine
equally spaced nodes.
No heat can cross a symmetry line, and thus symmetry lines can be treated
as insulated surfaces and thus “mirrors” in the finite difference formulation.
Then the nodes in the middle of the symmetry lines can be treated as interior
nodes by using mirror images. Six of the nodes are boundary nodes, so we will
have to write energy balances to obtain their finite difference formulations. First
we partition the region among the nodes equitably by drawing dashed lines between the nodes through the middle. Then the region around a node surrounded
by the boundary or the dashed lines represents the volume element of the node.
Considering a unit depth and using the energy balance approach for the boundary nodes (again assuming all heat transfer into the volume element for convenience) and the formula for the interior nodes, the finite difference equations
for the nine nodes are determined as follows: 0
3 The properties of chimney are given to be k Taking x y l, it simplifies to –2 hi l
T1
k T2 T3 hi l
T
ki 0 0 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 289 289
CHAPTER 5 (b) Node 2. On the inner boundary, subjected to convection, Figure 5–33b k
Taking x y T1 T2
2
x
y x
(T
2i hi T2) 0 kx T2 (4) hi l
T
ki 2T4 3 4 5 (8) 7 8 9 Mirror T1
T2
T4 T4
T5
T8 T6
T7
T8 4T3
4T4
4T5 Mirror FIGURE 5–34
Converting the boundary
nodes 3 and 5 on symmetry lines to
interior nodes by using mirror images. (c) Nodes 3, 4, and 5. (Interior nodes, Fig. 5–34) Node 3: T4
Node 4: T3
Node 5: T4 Mirror
images Mirror
image hi l
T2
k 3 (4) 0 y l, it simplifies to T1 2 6 T4 1 0
0
0 (d ) Node 6. (On the outer boundary, subjected to convection and radiation) 0
ho
Taking x y T2 x T3 T 6
2
y k x
(T
2o k y T7 T6
2
x x4
(T
2 sky T6) T64) 0 l, it simplifies to 2 T3 ho l
T6
k ho l
T
ko l
k 4
(Tsky T64) (e) Node 7. (On the outer boundary, subjected to convection and radiation,
Fig. 5–35) k Taking x 2T4 y T6 y T6 T7
T4 T 7
kx
2
x
y
ho x(To T7) y T8 T7
2
x
4
x(Tsky T74) 4 k 0
6 l, it simplifies to 4 2ho l
T7
k 2ho l
T
ko T8 2 l
k 4
(Tsky T74) (f ) Node 8. Same as Node 7, except shift the node numbers up by 1 (replace
4 by 5, 6 by 7, 7 by 8, and 8 by 9 in the last relation) 2T5 T7 4 2ho l
T8
k 2ho l
T
ko T9 2 l
k 4
(Tsky T84) (g) Node 9. (On the outer boundary, subjected to convection and radiation,
Fig. 5–35) k y T8 T9
2
x 0 ho Insulation x
(T
2o T9) x4
(T
2 sky T94) 0 7 8 9 h, T
Tsky FIGURE 5–35
Schematics for energy balances on the
volume elements of nodes 7 and 9. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 290 290
HEAT TRANSFER Taking x y l, it simplifies to T8 1 ho l
T9
k ho l
T
ko l
k 4
(Tsky T94) This problem involves radiation, which requires the use of absolute temperature, and thus all temperatures should be expressed in Kelvin. Alternately, we
could use °C for all temperatures provided that the four temperatures in the radiation terms are expressed in the form (T
273)4. Substituting the given
quantities, the system of nine equations for the determination of nine unknown
nodal temperatures in a form suitable for use with the GaussSeidel iteration
method becomes T1
T2
T3
T4
T5
T6
T7
T8
T9 (T2
(T1
(T1
(T2
(2T4
(T2
(2T4
(2T5
(T8 T3 2865)/7
2T4 2865)/8
2T4 T6)/4
T3 T5 T7)/4
2T8)/4
T3 456.2 0.3645 10 9 T64)/3.5
T6 T8 912.4 0.729 10 9 T74)/7
T7 T9 912.4 0.729 10 9 T84)/7
456.2 0.3645 10 9 T94)/2.5 which is a system of nonlinear equations. Using an equation solver, its solution
is determined to be 23 40
55
60
55
40 T1
T4
T7 40 Temperature, °C
55
60
55 40 89 152 89 138 152 256 273 273 256 138 273 152 138 256 273 256 138 89 138 152 138 89 545.7 K
411.2 K
328.1 K 272.6°C
138.0°C
54.9°C T2
T5
T8 529.2 K
362.1 K
313.1 K 40 55 60 55 FIGURE 5–36
The variation of temperature
in the chimney. 40 T3
T6
T9 425.2 K
332.9 K
296.5 K 152.1°C
59.7°C
23.4°C 23 40
55 The variation of temperature in the chimney is shown in Figure 5–36.
Note that the temperatures are highest at the inner wall (but less than
300°C) and lowest at the outer wall (but more that 260 K), as expected.
The average temperature at the outer surface of the chimney weighed by the
surface area is 60 Twall, out 55
40 0.5 332.9 (0.5T6 T7
(0.5 1 328.1 313.1
3 23 256.1°C
89.0°C
39.9°C T8 0.5T9)
1 0.5)
0.5 296.5 318.6 K 23 Then the rate of heat loss through the 1mlong section of the chimney can be
determined approximately from cen58933_ch05.qxd 9/4/2002 11:42 AM Page 291 291
CHAPTER 5 ·
Q chimney 4
4
ho Ao (Twall, out To)
Ao (Twall, out Tsky)
2
(21 W/m · K)[4 (0.6 m)(1 m)](318.6 293)K
0.9(5.67 10 8 W/m2 · K4)
[4 (0.6 m)(1 m)](318.6 K)4 (260 K)4]
1291 702 1993 W We could also determine the heat transfer by finding the average temperature of
the inner wall, which is (272.6 256.1)/2 264.4°C, and applying Newton’s
law of cooling at that surface: ·
Q chimney hi Ai (Ti Twall, in)
(70 W/m2 · K)[4 (0.2 m)(1 m)](300 264.4)°C 1994 W The difference between the two results is due to the approximate nature of the
numerical analysis.
Discussion We used a relatively crude numerical model to solve this problem
to keep the complexities at a manageable level. The accuracy of the solution obtained can be improved by using a finer mesh and thus a greater number of
nodes. Also, when radiation is involved, it is more accurate (but more laborious)
to determine the heat losses for each node and add them up instead of using
the average temperature. 5–5 I TRANSIENT HEAT CONDUCTION So far in this chapter we have applied the finite difference method to steady
heat transfer problems. In this section we extend the method to solve transient
problems.
We applied the finite difference method to steady problems by discretizing
the problem in the space variables and solving for temperatures at discrete
points called the nodes. The solution obtained is valid for any time since under
steady conditions the temperatures do not change with time. In transient problems, however, the temperatures change with time as well as position, and
thus the finite difference solution of transient problems requires discretization
in time in addition to discretization in space, as shown in Figure 5–37. This is
done by selecting a suitable time step t and solving for the unknown nodal
temperatures repeatedly for each t until the solution at the desired time is obtained. For example, consider a hot metal object that is taken out of the oven
at an initial temperature of Ti at time t 0 and is allowed to cool in ambient
air. If a time step of t 5 min is chosen, the determination of the temperature distribution in the metal piece after 3 h requires the determination of the
temperatures 3 60/5 36 times, or in 36 time steps. Therefore, the computation time of this problem will be 36 times that of a steady problem. Choosing a smaller t will increase the accuracy of the solution, but it will also
increase the computation time.
In transient problems, the superscript i is used as the index or counter
of time steps, with i
0 corresponding to the specified initial condition.
In the case of the hot metal piece discussed above, i
1 corresponds to
t1
t 5 min, i 2 corresponds to t 2
t 10 min, and a general t
i
i
i
Tm +1 Tm +1 Tm +1
–1
+1 i+1
∆t i
i
Tm – 1 Tm i 1
0 ∆t
∆x
0 1 ∆x
m–1 i
Tm +1 ∆x
m m+1 x FIGURE 5–37
Finite difference formulation of timedependent problems involves discrete
points in time as well as space. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 292 292
HEAT TRANSFER
i
time step i corresponds to ti i t. The notation Tm is used to represent the
temperature at the node m at time step i.
The formulation of transient heat conduction problems differs from that of
steady ones in that the transient problems involve an additional term representing the change in the energy content of the medium with time. This additional term appears as a first derivative of temperature with respect to time in
the differential equation, and as a change in the internal energy content during
t in the energy balance formulation. The nodes and the volume elements in
transient problems are selected as they are in the steady case, and, again assuming all heat transfer is into the element for convenience, the energy balance on a volume element during a time interval t can be expressed as Heat transferred into
the volume element
from all of its surfaces
during t Heat generated
within the
volume element
during t The change in the
energy content of
the volume element
during t or
·
Q t t ·
G element Eelement (537) All sides ·
where the rate of heat transfer Q normally consists of conduction terms for
interior nodes, but may involve convection, heat flux, and radiation for boundary nodes.
Velement C T, where is density and C is
Noting that Eelement mC T
the specific heat of the element, dividing the earlier relation by t gives
·
Q ·
G element All sides Eelement
t Velement C T
t (538) or, for any node m in the medium and its volume element,
·
Q
All sides Volume element
(can be any shape) Node m ρ = density
V = volume
ρV = mass
C = specific heat
∆ T = temperature change i
i
∆U = ρVC∆T = ρVC(Tm+ 1 – Tm ) FIGURE 5–38
The change in the energy content of
the volume element of a node
during a time interval t. ·
G element Velement C i
Tm 1 i
Tm t (539) i
i
where Tm and Tm 1 are the temperatures of node m at times ti i t and ti 1
i
i
Tm represents the temperature change
(i
1) t, respectively, and Tm 1
of the node during the time interval t between the time steps i and i
1
(Fig. 5–38).
i
i
Note that the ratio (Tm 1 Tm)/ t is simply the finite difference approximation of the partial derivative T/ t that appears in the differential equations
of transient problems. Therefore, we would obtain the same result for the
finite difference formulation if we followed a strict mathematical approach
instead of the energy balance approach used above. Also note that the finite
difference formulations of steady and transient problems differ by the single
term on the right side of the equal sign, and the format of that term remains the
same in all coordinate systems regardless of whether heat transfer is one,
i
i
two, or threedimensional. For the special case of Tm 1 Tm (i.e., no change
in temperature with time), the formulation reduces to that of steady case, as
expected. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 293 293
CHAPTER 5 The nodal temperatures in transient problems normally change during each
time step, and you may be wondering whether to use temperatures at the previous time step i or the new time step i 1 for the terms on the left side of Eq.
5–39. Well, both are reasonable approaches and both are used in practice. The
finite difference approach is called the explicit method in the first case and
the implicit method in the second case, and they are expressed in the general
form as (Fig. 5–39)
·
Qi Explicit method: ·
G ielement Velement C i
Tm 1 i
Tm t All sides ·
Qi Implicit method: 1 ·i 1
G element Velement C i
Tm 1 i
Tm t All sides If expressed at i + 1: Implicit method · · ∑ Q + Gelement All sides i
i
Tm+ 1 – Tm
= ρVelementC ————–
∆t If expressed at i: Explicit method FIGURE 5–39
The formulation of explicit and
implicit methods differs at the time
step (previous or new) at which the
heat transfer and heat generation
terms are expressed. (540) (541) It appears that the time derivative is expressed in forward difference form in
the explicit case and backward difference form in the implicit case. Of course,
it is also possible to mix the two fundamental formulations of Eqs. 5–40 and
5–41 and come up with more elaborate formulations, but such formulations
offer little insight and are beyond the scope of this text. Note that both formulations are simply expressions between the nodal temperatures before and
i
after a time interval and are based on determining the new temperatures Tm 1
i
using the previous temperatures Tm. The explicit and implicit formulations
given here are quite general and can be used in any coordinate system regardless of the dimension of heat transfer. The volume elements in multidimensional cases simply have more surfaces and thus involve more terms in
the summation.
The explicit and implicit methods have their advantages and disadvantages,
and one method is not necessarily better than the other one. Next you will see
that the explicit method is easy to implement but imposes a limit on the allowable time step to avoid instabilities in the solution, and the implicit method
requires the nodal temperatures to be solved simultaneously for each time step
but imposes no limit on the magnitude of the time step. We will limit the discussion to one and twodimensional cases to keep the complexities at a manageable level, but the analysis can readily be extended to threedimensional
cases and other coordinate systems. A Plane wall ·
gm Volume
element
of node m i
Tm+ 1
i
Tm kA Transient Heat Conduction in a Plane Wall
Consider transient onedimensional heat conduction in a plane wall of thick·
ness L with heat generation g(x, t) that may vary with time and position and
constant conductivity k with a mesh size of x L/M and nodes 0, 1, 2, . . . ,
M in the xdirection, as shown in Figure 5–40. Noting that the volume element of a general interior node m involves heat conduction from two sides and
the volume of the element is Velement A x, the transient finite difference formulation for an interior node can be expressed on the basis of Eq. 5–39 as
kA Tm Tm 1 x kA Tm Tm 1 x ·
gm A x A xC i
Tm 1 i
Tm t (542) i
i
Tm – 1 – Tm
————– i
i
Tm + 1 – Tm
kA ————–
∆x ∆x ∆x
01 2 ∆x m –1 m m +1
∆x M –1 Mx FIGURE 5–40
The nodal points and volume elements
for the transient finite difference
formulation of onedimensional
conduction in a plane wall. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 294 294
HEAT TRANSFER Canceling the surface area A and multiplying by x/k, it simplifies to
Tm 2Tm 1 Tm ·
gm x2
k 1 x2 i
(T
tm 1 i
Tm) (543) where
k/ C is the thermal diffusivity of the wall material. We now define
a dimensionless mesh Fourier number as
αt
x2 (544) Then Eq. 5–43 reduces to
Tm 1 2Tm Tm ·
gm x2
k 1 i
Tm 1 i
Tm (545) Note that the left side of this equation is simply the finite difference formulation of the problem for the steady case. This is not surprising since the formui
i
Tm. Also, we are still not
lation must reduce to the steady case for Tm 1
committed to explicit or implicit formulation since we did not indicate the
time step on the left side of the equation. We now obtain the explicit finite difference formulation by expressing the left side at time step i as
i
Tm i
2Tm 1 i
Tm ·i
gm x2
k 1 i
Tm 1 i
Tm (explicit) i
This equation can be solved explicitly for the new temperature Tm
the name explicit method) to give
i
Tm 1 i
(Tm i
Tm 1) 1 1 ·i
gm x2
k i
2 ) Tm (1 (546) (and thus (547) for all interior nodes m 1, 2, 3, . . . , M 1 in a plane wall. Expressing the
left side of Eq. 5–45 at time step i 1 instead of i would give the implicit
finite difference formulation as
∆–
x
—
2
A i
hA(T – T0 ) ∆–
x
ρA — C
2 i
Tm T0i + 1 – T0i
————– i
2Tm 1 i
Tm ·i
gm 1
1 1 x2 i
Tm 1 i
Tm k (implicit) (548) ∆t which can be rearranged as ·
g0
i
T1 – T0i
kA ———
∆x i
Tm ∆x ∆x
0 1
1 1 2 … FIGURE 5–41
Schematic for the explicit finite
difference formulation of the
convection condition at the left
boundary of a plane wall. L x 1
1 (1 i
2 ) Tm 1 i
Tm 1
1 ·i
gm 1 k x2 i
Tm 0 (549) The application of either the explicit or the implicit formulation to each of the
M 1 interior nodes gives M 1 equations. The remaining two equations are
obtained by applying the same method to the two boundary nodes unless, of
course, the boundary temperatures are specified as constants (invariant with
time). For example, the formulation of the convection boundary condition at
the left boundary (node 0) for the explicit case can be expressed as (Fig. 5–41)
hA(T T0i) kA T1i T0i
x x
·
gi0 A
2 i1 A x T0
C
2 T0i
t (550) cen58933_ch05.qxd 9/4/2002 11:42 AM Page 295 295
CHAPTER 5 which simplifies to
T0i 1 1 2 2 hx i
T0
k 2 T1i 2 hx
T
k ·
g i0 x 2
k (551) Note that in the case of no heat generation and
0.5, the explicit
i
finite difference formulation for a general interior node reduces to Tm 1
i
i
(Tm 1 Tm 1)/2, which has the interesting interpretation that the temperature
of an interior node at the new time step is simply the average of the temperatures of its neighboring nodes at the previous time step.
Once the formulation (explicit or implicit) is complete and the initial condition is specified, the solution of a transient problem is obtained by marching
in time using a step size of t as follows: select a suitable time step t and determine the nodal temperatures from the initial condition. Taking the initial
i
i
temperatures as the previous solution Tm at t 0, obtain the new solution Tm 1
at all nodes at time t
t using the transient finite difference relations. Now
i
using the solution just obtained at t
t as the previous solution Tm, obtain
i1
the new solution Tm at t 2 t using the same relations. Repeat the process
until the solution at the desired time is obtained. Stability Criterion for Explicit Method: Limitation on t
The explicit method is easy to use, but it suffers from an undesirable feature
that severely restricts its utility: the explicit method is not unconditionally stable, and the largest permissible value of the time step t is limited by the stability criterion. If the time step t is not sufficiently small, the solutions
obtained by the explicit method may oscillate wildly and diverge from the actual solution. To avoid such divergent oscillations in nodal temperatures, the
value of t must be maintained below a certain upper limit established by the
stability criterion. It can be shown mathematically or by a physical argument
based on the second law of thermodynamics that the stability criterion is sati
i
isfied if the coefficients of all Tm in the Tm 1 expressions (called the primary
coefficients) are greater than or equal to zero for all nodes m (Fig. 5–42). Of
i
course, all the terms involving Tm for a particular node must be grouped together before this criterion is applied.
Different equations for different nodes may result in different restrictions on
the size of the time step t, and the criterion that is most restrictive should be
used in the solution of the problem. A practical approach is to identify the
equation with the smallest primary coefficient since it is the most restrictive
and to determine the allowable values of t by applying the stability criterion
to that equation only. A t value obtained this way will also satisfy the stability criterion for all other equations in the system.
For example, in the case of transient onedimensional heat conduction in a
plane wall with specified surface temperatures, the explicit finite difference
equations for all the nodes (which are interior nodes) are obtained from
i
i
Eq. 5–47. The coefficient of Tm in the Tm 1 expression is 1 2 , which is
independent of the node number m, and thus the stability criterion for all
nodes in this case is 1 2
0 or
t
x2 1
2 interior nodes, onedimensional heat
transfer in rectangular coordinates (552) Explicit formulation:
T0i 1 a0T0i ··· T1i 1 a1T1i ··· M
i
Tm 1 i
amTm ··· M
i
TM 1 i
aMTM ··· Stability criterion:
am 0, m 0, 1, 2, . . . m, . . . M FIGURE 5–42
The stability criterion of the
explicit method requires all primary
coefficients to be positive or zero. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 296 296
HEAT TRANSFER When the material of the medium and thus its thermal diffusivity is known
and the value of the mesh size x is specified, the largest allowable value of
the time step t can be determined from this relation. For example, in the case
of a brick wall (
0.45 10 6 m2/s) with a mesh size of x 0.01 m, the
upper limit of the time step is
t (0.01 m)2
2(0.45 10 6 m2/s) 1 x2
2 111s 1.85 min The boundary nodes involving convection and/or radiation are more restrictive than the interior nodes and thus require smaller time steps. Therefore,
the most restrictive boundary node should be used in the determination of the
maximum allowable time step t when a transient problem is solved with
the explicit method.
To gain a better understanding of the stability criterion, consider the explicit
finite difference formulation for an interior node of a plane wall (Eq. 5–47) for
the case of no heat generation,
i
Tm 80°C
50°C 50°C
20°C m–1 m m–1 m+1 m m+1 Time step: i + 1 Time step: i FIGURE 5–43
The violation of the stability criterion
in the explicit method may result in
the violation of the second law of
thermodynamics and thus
divergence of solution. 0°C k = 28 W/m·°C
·
g = 5 × 106 W/m3 h
T α = 12.5 × 10 – 6 m2/s
∆x ∆x
0 0 1 L
2 Tinitial = 200°C FIGURE 5–44
Schematic for Example 5–5. x i
(Tm 1 i
Tm 1) (1 i
2 )Tm i
i
Assume that at some time step i the temperatures Tm 1 and Tm 1 are equal but
i
i
i
i
less than Tm (say, Tm 1 Tm 1 50°C and Tm 80°C). At the next time
step, we expect the temperature of node m to be between the two values (say,
70°C). However, if the value of exceeds 0.5 (say,
1), the temperature of
node m at the next time step will be less than the temperature of the neighboring nodes (it will be 20°C), which is physically impossible and violates the
second law of thermodynamics (Fig. 5–43). Requiring the new temperature of
node m to remain above the temperature of the neighboring nodes is equivalent to requiring the value of to remain below 0.5.
The implicit method is unconditionally stable, and thus we can use any time
step we please with that method (of course, the smaller the time step, the better the accuracy of the solution). The disadvantage of the implicit method is
that it results in a set of equations that must be solved simultaneously for each
time step. Both methods are used in practice. EXAMPLE 5–5
Uranium plate 1 Transient Heat Conduction in a Large Uranium
Plate Consider a large uranium plate of thickness L 4 cm, thermal conductivity k
28 W/m · °C, and thermal diffusivity
12.5 10 6 m2/s that is initially at
a uniform temperature of 200°C. Heat is generated uniformly in the plate at a
·
constant rate of g
5 106 W/m3. At time t
0, one side of the plate is
brought into contact with iced water and is maintained at 0°C at all times, while
the other side is subjected to convection to an environment at T
30°C with
a heat transfer coefficient of h 45 W/m2 · °C, as shown in Figure 5–44. Considering a total of three equally spaced nodes in the medium, two at the boundaries and one at the middle, estimate the exposed surface temperature of the
plate 2.5 min after the start of cooling using (a) the explicit method and (b) the
implicit method. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 297 297
CHAPTER 5 SOLUTION We have solved this problem in Example 5–1 for the steady case,
and here we repeat it for the transient case to demonstrate the application of
the transient finite difference methods. Again we assume onedimensional heat
transfer in rectangular coordinates and constant thermal conductivity. The number of nodes is specified to be M 3, and they are chosen to be at the two surfaces of the plate and at the middle, as shown in the figure. Then the nodal
spacing x becomes 0.04 m
31 L x M 1 0.02 m We number the nodes as 0, 1, and 2. The temperature at node 0 is given to be
T0 0°C at all times, and the temperatures at nodes 1 and 2 are to be determined. This problem involves only two unknown nodal temperatures, and thus
we need to have only two equations to determine them uniquely. These equations are obtained by applying the finite difference method to nodes 1 and 2.
(a) Node 1 is an interior node, and the explicit finite difference formulation at
that node is obtained directly from Eq. 5–47 by setting m 1: T1i 1 T2i) (T0 ·
g1 x2
k 2 ) T1i (1 (1) Node 2 is a boundary node subjected to convection, and the finite difference
formulation at that node is obtained by writing an energy balance on the volume
element of thickness x/2 at that boundary by assuming heat transfer to be into
the medium at all sides (Fig. 5–45): T2i) hA(T kA T1i T2i
x i1 x
·
g2 A
2 A x T2
C
2 A T2i i
i
T1 – T2
kA ———
∆x x Dividing by kA/2 x and using the definitions of thermal diffusivity
the dimensionless mesh Fourier number
t/( x)2 gives k/ C and T2i) i
which can be solved for T2 T2i 1 1 2 1 2 2(T1i T2i) ·
g2 x2
k T2i 1 T2i to give hx i
T2
k 2T1i 2 ·
g2 x2
k hx
T
k (2) Note that we did not use the superscript i for quantities that do not change with
time. Next we need to determine the upper limit of the time step t from the
i
stability criterion, which requires the coefficient of T1 in Equation 1 and the coi
efficient of T2 in the second equation to be greater than or equal to zero. The
i
coefficient of T2 is smaller in this case, and thus the stability criterion for this
problem can be expressed as 1 2 2 hx
k 0 → 2(1 1
→
h x/k) t 2 (1 x2
h x/k) i
hA(T – T2 ) i
T2 + 1
i
T2 0 2h x
(T
k ·
g2 Volume element
of node 2 0 1 ∆–
x
—
2 2 x FIGURE 5–45
Schematic for the explicit
finite difference formulation of the
convection condition at the right
boundary of a plane wall. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 298 298
HEAT TRANSFER since
t/( x)2. Substituting the given quantities, the maximum allowable
value of the time step is determined to be t 2(12.5 10 6 m2/s)[1 (0.02 m)2
(45 W/m2 · °C)(0.02 m)/28 W/m · °C] 15.5 s Therefore, any time step less than 15.5 s can be used to solve this problem. For
convenience, let us choose the time step to be t
15 s. Then the mesh
Fourier number becomes αt
( x)2 (12.5 10 6 m2/s)(15 s)
(0.02 m)2 0.46875 (for t 15 s) Substituting this value of and other given quantities, the explicit finite difference equations (1) and (2) developed here reduce to T1i
T2i 1
1 0.0625T1i
0.9375T1i 0.46875T2i 33.482
0.032366T2i 34.386 The initial temperature of the medium at t 0 and i 0 is given to be 200°C
0
0
1
throughout, and thus T1
200°C. Then the nodal temperatures at T1
T2
1
and T2 at t
t 15 s are determined from these equations to be T11
T21 2
2
Similarly, the nodal temperatures T1 and T2 at t
determined to be TABLE 5–2
The variation of the nodal
temperatures in Example 5–5 with
time obtained by the explicit
method Time
Step, i
0
1
2
3
4
5
6
7
8
9
10
20
30
40 Time,
s
0
15
30
45
60
75
90
105
120
135
150
300
450
600 0.0625T10 0.46875T20 33.482
0.0625 200 0.46875 200 33.482 139.7°C
0.9375T10 0.032366T20 34.386
0.9375 200 0.032366 200 34.386 228.4°C Node
Temperature, °C
i
T1
200.0
139.7
149.3
123.8
125.6
114.6
114.3
109.5
108.9
106.7
106.3
103.8
103.7
103.7 i
T2
200.0
228.4
172.8
179.9
156.3
157.1
146.9
146.3
141.8
141.1
139.0
136.1
136.0
136.0 T12
T22 2t 2 15 30 s are 0.0625T11 0.46875T21 33.482
0.0625 139.7 0.46875 228.4 33.482 149.3°C
0.9375T11 0.032366T21 34.386
0.9375 139.7 0.032366 228.4 34.386 172.8°C Continuing in the same manner, the temperatures at nodes 1 and 2 are determined for i 1, 2, 3, 4, 5, . . . , 50 and are given in Table 5–2. Therefore,
the temperature at the exposed boundary surface 2.5 min after the start of
cooling is
2.5
TL min T210 139.0°C (b) Node 1 is an interior node, and the implicit finite difference formulation at
that node is obtained directly from Eq. 5–49 by setting m 1: T0 (1 2 ) T1i 1 T2i 1 ·
g0 x2
k T1i 0 (3) Node 2 is a boundary node subjected to convection, and the implicit finite difference formulation at that node can be obtained from this formulation by expressing the left side of the equation at time step i 1 instead of i as cen58933_ch05.qxd 9/4/2002 11:42 AM Page 299 299
CHAPTER 5 2h x
(T
k T2i 1) 2(T1i 1 T2i 1) hx i
T2
k 1 2 ·
g2 x2
k T2i 1 T2i which can be rearranged as 2 T1i 1 1 2 2 ·
g2 x2
k hx
T
k T2i 0 (4) Again we did not use the superscript i or i 1 for quantities that do not change
with time. The implicit method imposes no limit on the time step, and thus we
can choose any value we want. However, we will again choose t
15 s, and
thus
0.46875, to make a comparison with part (a) possible. Substituting
this value of and other given quantities, the two implicit finite difference
equations developed here reduce to –1.9375T1i 1 0.46875T2i
0.9375T1i 1 1.9676T2i
0
Again T1
and for i 1
1 T1i
T2i 33.482
34.386 0
0 0
T2
200°C at t
0 and i
0 because of the initial condition,
0, these two equations reduce to 1.9375T11 0.46875T21
0.9375T11 1.9676T21 200
200 33.482
34.386 1
1
The unknown nodal temperatures T1 and T2 at t
solving these two equations simultaneously to be T11
Similarly, for i 168.8°C T21 and t 0
0 15 s are determined by 199.6°C TABLE 5–3 1, these equations reduce to 1.9375T12 0.46875T22
0.9375T12 1.9676T22 168.8
199.6 33.482
34.386 0
0 2
2
The unknown nodal temperatures T1 and T2 at t
t
2 15
determined by solving these two equations simultaneously to be T12 150.5°C and T22 30 s are 190.6°C Continuing in this manner, the temperatures at nodes 1 and 2 are determined
for i
2, 3, 4, 5, . . . , 40 and are listed in Table 5–3, and the temperature
at the exposed boundary surface (node 2) 2.5 min after the start of cooling is
obtained to be
2.5
TL min T210 143.9°C which is close to the result obtained by the explicit method. Note that either
method could be used to obtain satisfactory results to transient problems, except, perhaps, for the first few time steps. The implicit method is preferred
when it is desirable to use large time steps, and the explicit method is preferred
when one wishes to avoid the simultaneous solution of a system of algebraic
equations. The variation of the nodal
temperatures in Example 5–5 with
time obtained by the implicit
method Time
Step, i Time,
s 0
1
2
3
4
5
6
7
8
9
10
20
30
40 0
15
30
45
60
75
90
105
120
135
150
300
450
600 Node
Temperature, °C
i
T1 i
T2 200.0
168.8
150.5
138.6
130.3
124.1
119.5
115.9
113.2
111.0
109.4
104.2
103.8
103.8 200.0
199.6
190.6
180.4
171.2
163.6
157.6
152.8
149.0
146.1
143.9
136.7
136.1
136.1 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 300 300
HEAT TRANSFER EXAMPLE 5–6 South Warm
air Sun’s
rays Trombe
wall
Heat
loss Heat
gain
Vent
Cool
air Glazing FIGURE 5–46
Schematic of a Trombe wall
(Example 5–6). TABLE 5–4
The hourly variation of monthly
average ambient temperature and
solar heat flux incident on a vertical
surface for January in Reno, Nevada
Time
of
Day
7
10
1
4
7
10
1
4 Ambient
Solar
Temperature, Radiation,
°F
Btu/h · ft2 AM–10 AM
AM–1 PM
PM–4 PM
PM–7 PM
PM–10 PM
PM–1 AM
AM–4 AM
AM–7 AM 33
43
45
37
32
27
26
25 114
242
178
0
0
0
0
0 Solar Energy Storage in Trombe Walls Dark painted thick masonry walls called Trombe walls are commonly used on
south sides of passive solar homes to absorb solar energy, store it during the
day, and release it to the house during the night (Fig. 5–46). The idea was proposed by E. L. Morse of Massachusetts in 1881 and is named after Professor
Felix Trombe of France, who used it extensively in his designs in the 1970s.
Usually a single or double layer of glazing is placed outside the wall and transmits most of the solar energy while blocking heat losses from the exposed surface of the wall to the outside. Also, air vents are commonly installed at the
bottom and top of the Trombe walls so that the house air enters the parallel flow
channel between the Trombe wall and the glazing, rises as it is heated, and enters the room through the top vent.
Consider a house in Reno, Nevada, whose south wall consists of a 1ftthick
Trombe wall whose thermal conductivity is k
0.40 Btu/h · ft · °F and whose
thermal diffusivity is
4.78 10 6 ft2/s. The variation of the ambient tem·
perature Tout and the solar heat flux q solar incident on a southfacing vertical surface throughout the day for a typical day in January is given in Table 5–4 in 3h
intervals. The Trombe wall has single glazing with an absorptivitytransmissivity
product of
0.77 (that is, 77 percent of the solar energy incident is absorbed by the exposed surface of the Trombe wall), and the average combined
heat transfer coefficient for heat loss from the Trombe wall to the ambient is determined to be hout 0.7 Btu/h · ft2 · °F. The interior of the house is maintained
at Tin 70°F at all times, and the heat transfer coefficient at the interior surface of the Trombe wall is hin
1.8 Btu/h · ft2 · °F. Also, the vents on the
Trombe wall are kept closed, and thus the only heat transfer between the air in
the house and the Trombe wall is through the interior surface of the wall. Assuming the temperature of the Trombe wall to vary linearly between 70°F at the
interior surface and 30°F at the exterior surface at 7 AM and using the explicit
finite difference method with a uniform nodal spacing of x 0.2 ft, determine
the temperature distribution along the thickness of the Trombe wall after 12,
24, 36, and 48 h. Also, determine the net amount of heat transferred to the
house from the Trombe wall during the first day and the second day. Assume the
wall is 10 ft high and 25 ft long. SOLUTION The passive solar heating of a house through a Trombe wall is considered. The temperature distribution in the wall in 12h intervals and the
amount of heat transfer during the first and second days are to be determined.
Assumptions 1 Heat transfer is onedimensional since the exposed surface of
the wall is large relative to its thickness. 2 Thermal conductivity is constant.
3 The heat transfer coefficients are constant.
Properties The wall properties are given to be k
0.40 Btu/h · ft · °F,
4.78 10 6 ft2/s, and
0.77.
Analysis The nodal spacing is given to be x 0.2 ft, and thus the total number of nodes along the Trombe wall is M L
x 1 1 ft
0.2 ft 1 6 We number the nodes as 0, 1, 2, 3, 4, and 5, with node 0 on the interior surface of the Trombe wall and node 5 on the exterior surface, as shown in Figure
5–47. Nodes 1 through 4 are interior nodes, and the explicit finite difference
formulations of these nodes are obtained directly from Eq. 5–47 to be cen58933_ch05.qxd 9/4/2002 11:42 AM Page 301 301
CHAPTER 5 Node 1 (m
Node 2 (m
Node 3 (m
Node 4 (m T1i
T2i
T3i
T4i 1):
2):
3):
4): 1 (T0i
(T1i
(T2i
(T3i 1
1
1 T2i)
T3i)
T4i)
T5i) (1
(1
(1
(1 )T1i
)T2i
)T3i
)T4i 2
2
2
2 Trombe wall (2) k = 0.40 Btu/h·ft·°F
α = 4.78 × 10 – 6 ft2/s (3)
(4) The interior surface is subjected to convection, and thus the explicit formulation of node 0 can be obtained directly from Eq. 5–51 to be T0i 1 1 hin x i
T0
2
k 2 Substituting the quantities hin,
into this equation gives T0i 1 ·
Aq isolar T5i) 0 (2T1i 126.0) (5) kA T4i T5i i1 A x x T5
C
2 T5i
t (553) which simplifies to T5i 1 1 2 hout x i
T5
k 2 2 T4i 2 hout x i
Tout
k 2 ·i
qsolar x
(554)
k where
t/ x2 is the dimensionless mesh Fourier number. Note that we
kept the superscript i for quantities that vary with time. Substituting the quantities hout, x, k, and , which do not change with time, into this equation gives T5i 1 2.70 ) T5i (1 (2T4i i
0.70Tout 0.770q·isolar) (6) ·
where the unit of q isolar is Btu/h · ft2.
Next we need to determine the upper limit of the time step t from the stability criterion since we are using the explicit method. This requires the identification of the smallest primary coefficient in the system. We know that the
boundary nodes are more restrictive than the interior nodes, and thus we examine the formulations of the boundary nodes 0 and 5 only. The smallest and thus
i
the most restrictive primary coefficient in this case is the coefficient of T0 in the
formulation of node 0 since 1 3.8
1 2.7 , and thus the stability criterion for this problem can be expressed as 1 3.80 0 → αx
x2 1
3.80 Substituting the given quantities, the maximum allowable value of the time step
is determined to be t x2
3.80α 3.80 30°F ∆ x = 0.2 ft The exterior surface of the Trombe wall is subjected to convection as well as to
heat flux. The explicit finite difference formulation at that boundary is obtained
by writing an energy balance on the volume element represented by node 5,
i
hout A(Tout hout, Tout 0 1 2 3 4 5L x x, k, and Tin, which do not change with time, 3.80 ) T0i (1 2 Initial temperature
distribution at
7 AM (t = 0) 70°F
hin, Tin hin x
2
Tin
k T1i ·
qsolar (1) (0.2 ft)2
(4.78 10 6 ft2/s) 2202 s FIGURE 5–47
The nodal network for the Trombe
wall discussed in Example 5–6. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 302 302
HEAT TRANSFER Therefore, any time step less than 2202 s can be used to solve this problem.
For convenience, let us choose the time step to be t 900 s 15 min. Then
the mesh Fourier number becomes αt
( x)2 (4.78 10 6 ft2/s)(900 s)
(0.2 ft)2 0.10755 (for t 15 min) Initially (at 7 AM or t 0), the temperature of the wall is said to vary linearly between 70°F at node 0 and 30°F at node 5. Noting that there are five nodal
spacings of equal length, the temperature change between two neighboring
nodes is (70 30)°F/5 8°F. Therefore, the initial nodal temperatures are T00
T30 70°F,
46°F, Then the nodal temperatures at t
from these equations to be T01
T11
T21
T31
Temperature T41 °F
170 T51 150 1st day
2nd day 130
7 PM
110 1 AM 90
1 PM
70
50
Initial
temperature
30 0 7 AM 0.2
0.4
0.6
0.8
Distance along the Trombe wall FIGURE 5–48
The variation of temperatures in
the Trombe wall discussed
in Example 5–6. 1 ft T10
T40 62°F,
38°F,
t T20
T50 54°F,
30°F 15 min (at 7:15 AM) are determined (1 3.80 ) T00
(2T10 126.0)
(1 3.80 0.10755) 70 0.10755(2 62 126.0) 68.3° F
(T00 T20) (1 2 ) T10
0.10755(70 54) (1 2 0.10755)62 62°F
(T10 T30) (1 2 ) T20
0.10755(62 46) (1 2 0.10755)54 54°F
(T20 T40) (1 2 ) T30
0.10755(54 38) (1 2 0.10755)46 46°F
(T30 T50) (1 2 ) T40
0.10755(46 30) (1 2 0.10755)38 38°F
0
·
(1 2.70 ) T50
(2T40 0.70Tout 0.770q 0 )
solar
(1 2.70 0.10755)30 0.10755(2 38 0.70 33 0.770
41.4°F 114) Note that the inner surface temperature of the Trombe wall dropped by 1.7°F
and the outer surface temperature rose by 11.4°F during the first time step
while the temperatures at the interior nodes remained the same. This is typical
of transient problems in mediums that involve no heat generation. The nodal
temperatures at the following time steps are determined similarly with the help
of a computer. Note that the data for ambient temperature and the incident
solar radiation change every 3 hours, which corresponds to 12 time steps,
and this must be reflected in the computer program. For example, the value of
·
·
·
q isolar must be taken to be q isolar 75 for i 1–12, q isolar
242 for i 13–24,
·i
·i
q solar
178 for i 25–36, and q solar
0 for i 37–96.
The results after 6, 12, 18, 24, 30, 36, 42, and 48 h are given in Table 5–5
and are plotted in Figure 5–48 for the first day. Note that the interior temperature of the Trombe wall drops in early morning hours, but then rises as the solar
energy absorbed by the exterior surface diffuses through the wall. The exterior
surface temperature of the Trombe wall rises from 30 to 142°F in just 6 h because of the solar energy absorbed, but then drops to 53°F by next morning as
a result of heat loss at night. Therefore, it may be worthwhile to cover the outer
surface at night to minimize the heat losses. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 303 303
CHAPTER 5 TABLE 5–5
The temperatures at the nodes of a Trombe wall at various times Time
0
6
12
18
24
30
36
42
48 h
h
h
h
h
h
h
h
h (7
(1
(7
(1
(7
(1
(7
(1
(7 Nodal Temperatures, °F Time
Step, i
AM)
PM)
PM)
AM)
AM)
PM)
PM)
AM)
AM) T0 T1 T2 T3 T4 T5 0
24
48
72
96
120
144
168
192 70.0
65.3
71.6
73.3
71.2
70.3
75.4
75.8
73.0 62.0
61.7
74.2
75.9
71.9
71.1
81.1
80.7
75.1 54.0
61.5
80.4
77.4
70.9
74.3
89.4
83.5
72.2 46.0
69.7
88.4
76.3
67.7
84.2
98.2
83.0
66.0 38.0
94.1
91.7
71.2
61.7
108.3
101.0
77.4
66.0 30.0
142.0
82.4
61.2
53.0
153.2
89.7
66.2
56.3 The rate of heat transfer from the Trombe wall to the interior of the house during each time step is determined from Newton’s law using the average temperature at the inner surface of the wall (node 0) as QiTrombe wall ·
Q iTrombe wall t hin A(T0i Tin) t hin A[(T0i T0i 1)/2 Therefore, the amount of heat transfer during the first time step (i
during the first 15min period is Q1
Trombe wall hin A[(T01 T00)/2 Tin] t
(1.8 Btu/h · ft2 · °F)(10 25 ft2)[(68.3
95.6 Btu 70)/2 Tin] t
1) or 70°F](0.25 h) The negative sign indicates that heat is transferred to the Trombe wall from the
air in the house, which represents a heat loss. Then the total heat transfer during a specified time period is determined by adding the heat transfer amounts
for each time step as
I QTrombe wall
i1 I ·
Q iTrombe wall hin A[(T0i T0i 1)/2 Tin] t (555) i1 where I is the total number of time intervals in the specified time period. In this
case I 48 for 12 h, 96 for 24 h, and so on. Following the approach described
here using a computer, the amount of heat transfer between the Trombe wall
and the interior of the house is determined to be QTrombe wall
QTrombe wall
QTrombe wall
QTrombe wall 17, 048 Btu after 12 h
2483 Btu after 24 h
5610 Btu after 36 h
34, 400 Btu after 48 h ( 17, 078 Btu during the first 12 h)
(14, 565 Btu during the second 12 h)
(8093 Btu during the third 12 h)
(28, 790 Btu during the fourth 12 h) Therefore, the house loses 2483 Btu through the Trombe wall the first day as a
result of the low startup temperature but delivers a total of 36,883 Btu of heat
to the house the second day. It can be shown that the Trombe wall will deliver
even more heat to the house during the third day since it will start the day at a
higher average temperature. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 304 304
HEAT TRANSFER TwoDimensional Transient Heat Conduction m, n + 1 n+1
∆y Volume
element m – 1, n n ·
gm, n m, n m + 1, n ∆y
m, n – 1 n–1 ∆x y m–1 ∆x
m m+1 x FIGURE 5–49
The volume element of a
general interior node (m, n) for twodimensional transient conduction
in rectangular coordinates. Consider a rectangular region in which heat conduction is significant in the
x and ydirections, and consider a unit depth of z 1 in the zdirection.
·
Heat may be generated in the medium at a rate of g(x, y, t), which may vary
with time and position, with the thermal conductivity k of the medium assumed to be constant. Now divide the xyplane of the region into a rectangular mesh of nodal points spaced x and y apart in the x and ydirections,
respectively, and consider a general interior node (m, n) whose coordinates are
x m x and y n y, as shown in Figure 5–49. Noting that the volume element centered about the general interior node (m, n) involves heat conduction
from four sides (right, left, top, and bottom) and the volume of the element is
x
y1
x y, the transient finite difference formulation for
Velement
a general interior node can be expressed on the basis of Eq. 5–39 as
ky Tm kx Tm, n 1, n kx x
Tm, n Tm, n 1 Tm 1, n Tm 1, n Tm, n y 1 Tm, n 1 ky y ·
gm, n x y y Taking a square mesh ( x
simplifying, Tm, n x yC Tm Tm, n 1, n x
i
Tm 1 i
Tm (556) t l ) and dividing each term by k gives after Tm, n 1 4Tm, n ·
gm, nl 2
k i
Tm 1 i
Tm (557) where again
k/ C is the thermal diffusivity of the material and
t/l 2
is the dimensionless mesh Fourier number. It can also be expressed in terms of
the temperatures at the neighboring nodes in the following easytoremember
form:
Tleft Ttop Tright Tbottom 4Tnode ·
gnodel 2
k i1
Tnode i
Tnode (558) Again the left side of this equation is simply the finite difference formulation
of the problem for the steady case, as expected. Also, we are still not committed to explicit or implicit formulation since we did not indicate the time
step on the left side of the equation. We now obtain the explicit finite difference formulation by expressing the left side at time step i as
i
Tleft i
Ttop i
Tright i
Tbottom i
4Tnode ·i
gnodel 2
k i1
Tnode i
Tnode (559) Expressing the left side at time step i 1 instead of i would give the implicit
formulation. This equation can be solved explicitly for the new temperature
i1
Tnode to give
i1
Tnode i
(Tleft i
Ttop i
Tright i
Tbottom) (1 i
4 ) Tnode ·i
gnodel 2
k for all interior nodes (m, n) where m 1, 2, 3, . . . , M 1 and n
3, . . . , N 1 in the medium. In the case of no heat generation and (560) 1, 2,
the 1
,
4 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 305 305
CHAPTER 5 explicit finite difference formulation for a general interior node reduces to
i1
i
i
i
i
Tnode (Tleft Ttop Tright Tbottom)/4, which has the interpretation that the
temperature of an interior node at the new time step is simply the average
of the temperatures of its neighboring nodes at the previous time step
(Fig. 5–50).
i
i
The stability criterion that requires the coefficient of Tm in the Tm 1 expression to be greater than or equal to zero for all nodes is equally valid for twoor threedimensional cases and severely limits the size of the time step t that
can be used with the explicit method. In the case of transient twodimensional
i
i
heat transfer in rectangular coordinates, the coefficient of Tm in the Tm 1 expression is 1 4 , and thus the stability criterion for all interior nodes in this
case is 1 4
0, or
t
l2 1
4 (interior nodes, twodimensional heat
transfer in rectangular coordinates) Time step i:
30°C i
Tm 20°C 40°C Node m 10°C
Time step i + 1: (561) i
Tm+ 1 25°C Node m where x
y l. When the material of the medium and thus its thermal
diffusivity are known and the value of the mesh size l is specified, the
largest allowable value of the time step t can be determined from the relation
above. Again the boundary nodes involving convection and/or radiation are
more restrictive than the interior nodes and thus require smaller time steps.
Therefore, the most restrictive boundary node should be used in the determination of the maximum allowable time step t when a transient problem is
solved with the explicit method.
The application of Eq. 5–60 to each of the (M 1) (N 1) interior nodes
gives (M 1) (N 1) equations. The remaining equations are obtained by
applying the method to the boundary nodes unless, of course, the boundary
temperatures are specified as being constant. The development of the transient
finite difference formulation of boundary nodes in two (or three) dimensional problems is similar to the development in the onedimensional case discussed earlier. Again the region is partitioned between the nodes by forming
volume elements around the nodes, and an energy balance is written for each
boundary node on the basis of Eq. 5–39. This is illustrated in Example 5–7.
EXAMPLE 5–7 FIGURE 5–50
In the case of no heat generation
1
and
, the temperature of an
4
interior node at the new time step is
the average of the temperatures of
its neighboring nodes at the
previous time step. Transient TwoDimensional Heat Conduction
in LBars Consider twodimensional transient heat transfer in an Lshaped solid body that
is initially at a uniform temperature of 90°C and whose cross section is given
in Figure 5–51. The thermal conductivity and diffusivity of the body are k
15 W/m · °C and
3.2 10 6 m2/s, respectively, and heat is generated in
·
the body at a rate of g
2 106 W/m3. The left surface of the body is insulated, and the bottom surface is maintained at a uniform temperature of 90°C
at all times. At time t 0, the entire top surface is subjected to convection to
ambient air at T
25°C with a convection coefficient of h
80 W/m2 · °C,
·
and the right surface is subjected to heat flux at a uniform rate of q R 5000
2
W/m . The nodal network of the problem consists of 15 equally spaced nodes
with x
y 1.2 cm, as shown in the figure. Five of the nodes are at the bottom surface, and thus their temperatures are known. Using the explicit method,
determine the temperature at the top corner (node 3) of the body after 1, 3, 5,
10, and 60 min. Convection
h, T y
1 2 3 4 5 6 7 8 10 11 12 13 14 15 ∆ x = ∆y = l ∆y ·
9 qR ∆y ∆x ∆x 90°C
∆x x
∆x ∆x FIGURE 5–51
Schematic and nodal network for
Example 5–7. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 306 306
HEAT TRANSFER h, T
1 SOLUTION This is a transient twodimensional heat transfer problem in rectangular coordinates, and it was solved in Example 5–3 for the steady case.
Therefore, the solution of this transient problem should approach the solution
for the steady case when the time is sufficiently large. The thermal conductivity and heat generation rate are given to be constants. We observe that all nodes
are boundary nodes except node 5, which is an interior node. Therefore, we will
have to rely on energy balances to obtain the finite difference equations. The region is partitioned among the nodes equitably as shown in the figure, and the
explicit finite difference equations are determined on the basis of the energy
balance for the transient case expressed as h, T
2 1 2 3 ·
Qi ·
G ielement Velement C i
Tm 1 t All sides 5 4
(a) Node 1 (b) Node 2 FIGURE 5–52
Schematics for energy balances on the
volume elements of nodes 1 and 2. i
Tm ·
·
The quantities h, T , g , and q R do not change with time, and thus we do not
need to use the superscript i for them. Also, the energy balance expressions are
simplified using the definitions of thermal diffusivity
k/ C and the dimensionless mesh Fourier number
t/l2, where x
y l.
(a) Node 1. (Boundary node subjected to convection and insulation, Fig.
5–52a) h x
(T
2 T1i) y T2i T1i
2
x k i k i x T4 T1
2
y
· xy
g1
22 i
x y T1
C
22 1 T1i
t Dividing by k/4 and simplifying, 2hl
(T
k T1i) 2(T2i i
which can be solved for T1 T1i 1 1 4 1 T1i) 2(T4i ·
g1l 2
k T1i) T1i 1 T1i to give 2 hl
T1i
k 2 T2i hl
T
k T4i ·
g1l 2
2k (b) Node 2. (Boundary node subjected to convection, Fig. 5–52b) T2i) h x(T y T3i T2i
2
x
y T1i T2i
k
2
x k kx
·
g2 T5i
y
y
x
2 i
Dividing by k/2, simplifying, and solving for T2 T2i 1 1 4 2 hl
T2i
k T1i T2i 1 T3i x y T2i
C
2 1 T2i
t gives 2T5i 2hl
T
k ·
g2l 2
k cen58933_ch05.qxd 9/4/2002 11:42 AM Page 307 307
CHAPTER 5 (c) Node 3. (Boundary node subjected to convection on two sides, Fig. 5–53a)
i i k x T6 T3
2
y y T2i T3i
2
x · xy
g3
22 y
(T
2 x
2 h T3i) k i
Dividing by k/4, simplifying, and solving for T3 T3i 1 1 4 hl
T3i
k 4 2 3 (5) 4 5 h, T
i
x y T3
22
1 1 T3i
t 10 6 gives T4i 2 1 Mirror h, T (a) Node 3 T6i 2 FIGURE 5–53
Schematics for energy balances on the
volume elements of nodes 3 and 4. ·
g3l 2
2k hl
T
k (b) Node 4 (d ) Node 4. (On the insulated boundary, and can be treated as an interior node,
Fig. 5–53b). Noting that T10 90°C, Eq. 5–60 gives T4i 1 4 ) T4i (1 2T5i T1i ·
g4l 2
k 90 (e) Node 5. (Interior node, Fig. 5–54a). Noting that T11
gives T5i 1 4 ) T5i (1 T2i T4i T6i 90°C, Eq. 5–60 ·
g5l 2
k 90 (f ) Node 6. (Boundary node subjected to convection on two sides, Fig. 5–54b) h y
(T
2 x
2 T6i)
i k
i x T3 T6
2
y y
2 T7i T6i kx x
·3xy
g6
4 i
T12 T6i y
3 x y T6i
C
4 ky
1 T5i T6i
x 3 2 4 5 6 T6i 1 1 4 3 4 2T3i 1 T6i
t gives 12 2T7i 4 90 4 FIGURE 5–54
Schematics for energy balances on the
volume elements of nodes 5 and 6. hl
T
k 3 ·
g6l 2
k (g) Node 7. (Boundary node subjected to convection, Fig. 5–55) T7i) h x(T k i
T13 T7i
y T8i T7i
kx
2
x
y
i
i
y
y T6 T7
·
k
g7 x
2
2
x i
Dividing by k/2, simplifying, and solving for T7 T7i 1 1 4 2 hl
T7i
k T6i T8i 1 h, T h, T
6 7 8 9
·
qR x y T7i
C
2 1 T7i
t gives 2 (b) Node 6 (a) Node 5 hl
Ti
3k 3 4T5i 7 6
5 11
i
Dividing by 3k/4, simplifying, and solving for T6 h, T 90 2hl
T
k ·
g7l 2
k 13 15 FIGURE 5–55
Schematics for energy balances on the
volume elements of nodes 7 and 9. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 308 308
HEAT TRANSFER (h) Node 8. This node is identical to node 7, and the finite difference formulation of this node can be obtained from that of node 7 by shifting the node numbers by 1 (i.e., replacing subscript m by subscript m 1). It gives T8i 1 1 4 hl
T8i
k 2 T7i T9i 2 ·
g8l 2
k 2hl
T
k 90 (i ) Node 9. (Boundary node subjected to convection on two sides, Fig. 5–55) h i
i
y
x T15 T9
q·R
k
2
2
y
k y T8i T9i
· xy
g9
2
22
x x
(T
2 T9i) i
Dividing by k/4, simplifying, and solving for T9 T9i 1 1 4 2 hl
T9i
k 2 T8i i
x y T9
C
22
1 1 T9i
t gives ·
qR l
k 90 hl
T
k ·
g9l 2
2k This completes the finite difference formulation of the problem. Next we need
to determine the upper limit of the time step t from the stability criterion,
i
i
which requires the coefficient of Tm in the Tm 1 expression (the primary coefficient) to be greater than or equal to zero for all nodes. The smallest primary coi
efficient in the nine equations here is the coefficient of T3 in the expression,
and thus the stability criterion for this problem can be expressed as 1 4 hl
k 4 → 0 1
4(1 hl/k) → t l2
4 (1 hl/k) since
t /l 2. Substituting the given quantities, the maximum allowable
value of the time step is determined to be t 4(3.2 10 6 2 m /s)[1 (0.012 m)2
(80 W/m2 · °C)(0.012 m)/(15 W/m · °C)] 10.6 s Therefore, any time step less than 10.6 s can be used to solve this problem. For
convenience, let us choose the time step to be t
10 s. Then the mesh
Fourier number becomes t
l (3.2 2 10 6 m2/s)(10 s)
(0.012 m)2 0.222 (for t 10 s) Substituting this value of and other given quantities, the developed transient
finite difference equations simplify to T1i
T2i
T3i
T4i
T5i 1
1
1
1
1 0.0836T1i 0.444(T2i T4i 11.2)
0.0836T2i 0.222(T1i T3i 2T5i 22.4)
0.0552T3i 0.444(T2i T6i 12.8)
0.112T4i 0.222(T1i 2T5i 109.2)
0.112T1i 0.222(T2i T4i T6i 109.2) cen58933_ch05.qxd 9/4/2002 11:42 AM Page 309 309
CHAPTER 5 T6i
T7i
T8i
T9i 1
1
1
1 0.0931T6i
0.0836T7i
0.0836T8i
0.0836T9i 0.074(2T3i 4T5i 2T7i
0.222(T6i T8i 202.4)
0.222(T7i T9i 202.4)
0.444(T8i 105.2) 424) Using the specified initial condition as the solution at time t
0 (for i
0),
sweeping through these nine equations will give the solution at intervals of
10 s. The solution at the upper corner node (node 3) is determined to be
100.2, 105.9, 106.5, 106.6, and 106.6°C at 1, 3, 5, 10, and 60 min, respectively. Note that the last three solutions are practically identical to the
solution for the steady case obtained in Example 5–3. This indicates that steady
conditions are reached in the medium after about 5 min. TOPIC OF SPECIAL INTEREST Controlling the Numerical Error
A comparison of the numerical results with the exact results for temperature distribution in a cylinder would show that the results obtained by a numerical method are approximate, and they may or may not be sufficiently
close to the exact (true) solution values. The difference between a numerical solution and the exact solution is the error involved in the numerical
solution, and it is primarily due to two sources:
• The discretization error (also called the truncation or formulation
error), which is caused by the approximations used in the formulation
of the numerical method.
• The roundoff error, which is caused by the computer ’s use of a
limited number of significant digits and continuously rounding (or
chopping) off the digits it cannot retain.
Below we discuss both types of errors. T(xm, t) Discretization Error
The discretization error involved in numerical methods is due to replacing
the derivatives by differences in each step, or the actual temperature distribution between two adjacent nodes by a straight line segment.
Consider the variation of the solution of a transient heat transfer problem
with time at a specified nodal point. Both the numerical and actual (exact)
solutions coincide at the beginning of the first time step, as expected, but
the numerical solution deviates from the exact solution as the time t increases. The difference between the two solutions at t
t is due to the approximation at the first time step only and is called the local discretization
error. One would expect the situation to get worse with each step since the
second step uses the erroneous result of the first step as its starting point
and adds a second local discretization error on top of it, as shown in Figure
5–56. The accumulation of the local discretization errors continues with the
increasing number of time steps, and the total discretization error at any Local
error Actual solution
T(x0, t) Global
error
T3 T2
T0 T1 Step 1
t0 Numerical
solution Step 2
t1 Step 3
t2 t3 Time FIGURE 5–56
The local and global discretization
errors of the finite difference
method at the third time step
at a specified nodal point. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 310 310
HEAT TRANSFER step is called the global or accumulated discretization error. Note that the
local and global discretization errors are identical for the first time step.
The global discretization error usually increases with the increasing number of steps, but the opposite may occur when the solution function
changes direction frequently, giving rise to local discretization errors of opposite signs, which tend to cancel each other.
To have an idea about the magnitude of the local discretization error,
consider the Taylor series expansion of the temperature at a specified nodal
point m about time ti,
T(xm, ti t) T(xm, ti) t T(xm, ti)
t 12
t
2 2 T(xm, ti)
t2 ··· (562) The finite difference formulation of the time derivative at the same nodal
point is expressed as
T(xm, ti)
t T(xm, ti t)
t T(xm, ti) i
Tm 1 i
Tm t (563) or
T(xm, ti t) T(xm, ti) t T(xm, ti)
t (564) which resembles the Taylor series expansion terminated after the first two
terms. Therefore, the third and later terms in the Taylor series expansion
represent the error involved in the finite difference approximation. For a
sufficiently small time step, these terms decay rapidly as the order of derivative increases, and their contributions become smaller and smaller. The
first term neglected in the Taylor series expansion is proportional to t 2,
and thus the local discretization error of this approximation, which is the
error involved in each step, is also proportional to t 2.
The local discretization error is the formulation error associated with a
single step and gives an idea about the accuracy of the method used. However, the solution results obtained at every step except the first one involve
the accumulated error up to that point, and the local error alone does not
have much significance. What we really need to know is the global discretization error. At the worst case, the accumulated discretization error
(t0/ t)( t)2
t0 t,
after I time steps during a time period t0 is i( t)2
which is proportional to t. Thus, we conclude that the local discretization
error is proportional to the square of the step size t 2 while the global discretization error is proportional to the step size t itself. Therefore, the
smaller the mesh size (or the size of the time step in transient problems),
the smaller the error, and thus the more accurate is the approximation. For
example, halving the step size will reduce the global discretization error by
half. It should be clear from the discussions above that the discretization error can be minimized by decreasing the step size in space or time as much
as possible. The discretization error approaches zero as the difference
quantities such as x and t approach the differential quantities such as dx
and dt. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 311 311
CHAPTER 5 Roundoff Error
If we had a computer that could retain an infinite number of digits for all
numbers, the difference between the exact solution and the approximate
(numerical) solution at any point would entirely be due to discretization error. But we know that every computer (or calculator) represents numbers
using a finite number of significant digits. The default value of the number
of significant digits for many computers is 7, which is referred to as single
precision. But the user may perform the calculations using 15 significant
digits for the numbers, if he or she wishes, which is referred to as double
precision. Of course, performing calculations in double precision will require more computer memory and a longer execution time.
In single precision mode with seven significant digits, a computer will
register the number 44444.666666 as 44444.67 or 44444.66, depending on
the method of rounding the computer uses. In the first case, the excess digits are said to be rounded to the closest integer, whereas in the second case
they are said to be chopped off. Therefore, the numbers a 44444.12345
and b 44444.12032 are equivalent for a computer that performs calculations using seven significant digits. Such a computer would give a b
0 instead of the true value 0.00313.
The error due to retaining a limited number of digits during calculations
is called the roundoff error. This error is random in nature and there is no
easy and systematic way of predicting it. It depends on the number of calculations, the method of rounding off, the type of computer, and even the
sequence of calculations.
In algebra you learned that a b c a c b, which seems quite
reasonable. But this is not necessarily true for calculations performed with
a computer, as demonstrated in Figure 5–57. Note that changing the sequence of calculations results in an error of 30.8 percent in just two operations. Considering that any significant problem involves thousands or even
millions of such operations performed in sequence, we realize that the accumulated roundoff error has the potential to cause serious error without
giving any warning signs. Experienced programmers are very much aware
of this danger, and they structure their programs to prevent any buildup of
the roundoff error. For example, it is much safer to multiply a number by
10 than to add it 10 times. Also, it is much safer to start any addition
process with the smallest numbers and continue with larger numbers. This
rule is particularly important when evaluating series with a large number of
terms with alternating signs.
The roundoff error is proportional to the number of computations performed during the solution. In the finite difference method, the number of
calculations increases as the mesh size or the time step size decreases.
Halving the mesh or time step size, for example, will double the number of
calculations and thus the accumulated roundoff error. Controlling the Error in Numerical Methods
The total error in any result obtained by a numerical method is the sum of
the discretization error, which decreases with decreasing step size, and the
roundoff error, which increases with decreasing step size, as shown in Figure 5–58. Therefore, decreasing the step size too much in order to get more Given:
a
b
c
D
E Find: 7777777
7777776
0.4444432
a
a b
c c
b Solution:
D 7777777 7777776 0.4444432
1 0.4444432
1.444443 (Correct result) E 7777777 0.4444432 7777776
7777777 7777776
1.000000 (In error by 30.8%) FIGURE 5–57
A simple arithmetic operation
performed with a computer
in single precision using seven
significant digits, which results in
30.8 percent error when the order
of operation is reversed.
Error Total
error Discretization
error
Roundoff error
Optimum
step size Step size FIGURE 5–58
As the mesh or time step size
decreases, the discretization error
decreases but the roundoff
error increases. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 312 312
HEAT TRANSFER accurate results may actually backfire and give less accurate results because of a faster increase in the roundoff error. We should be careful not to
let roundoff error get out of control by avoiding a large number of computations with very small numbers.
In practice, we will not know the exact solution of the problem, and thus
we will not be able to determine the magnitude of the error involved in the
numerical method. Knowing that the global discretization error is proportional to the step size is not much help either since there is no easy way of
determining the value of the proportionality constant. Besides, the global
discretization error alone is meaningless without a true estimate of the
roundoff error. Therefore, we recommend the following practical procedures to assess the accuracy of the results obtained by a numerical method.
• Start the calculations with a reasonable mesh size x (and time step
size t for transient problems) based on experience. Then repeat the
calculations using a mesh size of x/2. If the results obtained by
halving the mesh size do not differ significantly from the results
obtained with the full mesh size, we conclude that the discretization
error is at an acceptable level. But if the difference is larger than we
can accept, then we have to repeat the calculations using a mesh size
x/4 or even a smaller one at regions of high temperature gradients.
We continue in this manner until halving the mesh size does not cause
any significant change in the results, which indicates that the
discretization error is reduced to an acceptable level.
• Repeat the calculations using double precision holding the mesh size
(and the size of the time step in transient problems) constant. If the
changes are not significant, we conclude that the roundoff error is not
a problem. But if the changes are too large to accept, then we may try
reducing the total number of calculations by increasing the mesh size
or changing the order of computations. But if the increased mesh size
gives unacceptable discretization errors, then we may have to find a
reasonable compromise.
It should always be kept in mind that the results obtained by any numerical method may not reflect any trouble spots in certain problems that require special consideration such as hot spots or areas of high temperature
gradients. The results that seem quite reasonable overall may be in considerable error at certain locations. This is another reason for always repeating
the calculations at least twice with different mesh sizes before accepting
them as the solution of the problem. Most commercial software packages
have builtin routines that vary the mesh size as necessary to obtain highly
accurate solutions. But it is a good engineering practice to be aware of any
potential pitfalls of numerical methods and to examine the results obtained
with a critical eye. SUMMARY
Analytical solution methods are limited to highly simplified
problems in simple geometries, and it is often necessary to use a numerical method to solve real world problems with complicated geometries or nonuniform thermal conditions. The cen58933_ch05.qxd 9/4/2002 11:42 AM Page 313 313
CHAPTER 5 numerical finite difference method is based on replacing derivatives by differences, and the finite difference formulation of a
heat transfer problem is obtained by selecting a sufficient number of points in the region, called the nodal points or nodes,
and writing energy balances on the volume elements centered
about the nodes.
For steady heat transfer, the energy balance on a volume element can be expressed in general as
·
Q ·
gVelement 0 also be used to solve a system of equations simultaneously at
the press of a button.
The finite difference formulation of transient heat conduction problems is based on an energy balance that also accounts
for the variation of the energy content of the volume element
during a time interval t. The heat transfer and heat generation
terms are expressed at the previous time step i in the explicit
method, and at the new time step i 1 in the implicit method.
For a general node m, the finite difference formulations are
expressed as All sides whether the problem is one, two, or threedimensional. For
convenience in formulation, we always assume all heat transfer to be into the volume element from all surfaces toward the
node under consideration, except for specified heat flux whose
direction is already specified. The finite difference formulations for a general interior node under steady conditions are expressed for some geometries as follows:
Onedimensional
steady conduction
in a plane wall:
Twodimensional
steady
Ttop
T
conduction left
in rectangular
coordinates: Tm 2Tm
( x)2 1 Tright Tbottom Tm 1 4Tnode ·
gm
k 0 ·
gnodel 2
k 0 where x is the nodal spacing for the plane wall and x
y l is the nodal spacing for the twodimensional case. Insulated boundaries can be viewed as mirrors in formulation, and
thus the nodes on insulated boundaries can be treated as interior nodes by using mirror images.
The finite difference formulation at node 0 at the left boundary of a plane wall for steady onedimensional heat conduction
can be expressed as
·
Q left surface kA T1 T0
x ·
g0(A x/2) 0 ·
where A x/2 is the volume of the volume, g0 is the rate of heat
generation per unit volume at x 0, and A is the heat transfer
area. The form of the first term depends on the boundary condition at x 0 (convection, radiation, specified heat flux, etc.).
The finite difference formulation of heat conduction problems usually results in a system of N algebraic equations in N
unknown nodal temperatures that need to be solved simultaneously. There are numerous systematic approaches available in
the literature. Several widely available equation solvers can Explicit
method:
Implicit
method: ·
G ielement ·
Qi Velement C i
Tm 1 t All sides ·
Qi ·i 1
G element 1 Velement C i
Tm i
Tm 1 All sides i
Tm t i
i
where Tm and Tm 1 are the temperatures of node m at times
i
i
ti i t and ti 1 (i 1) t, respectively, and Tm 1 Tm represents the temperature change of the node during the time interval t between the time steps i and i 1. The explicit and
implicit formulations given here are quite general and can be
used in any coordinate system regardless of heat transfer being
one, two, or threedimensional.
The explicit formulation of a general interior node for oneand twodimensional heat transfer in rectangular coordinates
can be expressed as Onei
Tm 1
dimensional case:
Twoi1
Tnode
dimensional
case: i
(Tm i
(Tleft (1 1 i
Tm 1) (1 i
Tright i
Ttop i
4 ) Tnode i
2 ) Tm ·i
gm x2
k i
Tbottom) ·i
gnodel 2
k where
t
x2
is the dimensionless mesh Fourier number and
k/ C is the
thermal diffusivity of the medium.
The implicit method is inherently stable, and any value of t
can be used with that method as the time step. The largest value
of the time step t in the explicit method is limited by the stai
bility criterion, expressed as: the coefficients of all Tm in the
i1
Tm expressions (called the primary coefficients) must be
greater than or equal to zero for all nodes m. The maximum
value of t is determined by applying the stability criterion to
the equation with the smallest primary coefficient since it is the cen58933_ch05.qxd 9/4/2002 11:42 AM Page 314 314
HEAT TRANSFER most restrictive. For problems with specified temperatures or
heat fluxes at all the boundaries, the stability criterion can be 1
expressed as
for onedimensional problems and
2
the twodimensional problems in rectangular coordinates. 1
4 for REFERENCES AND SUGGESTED READING
1. D. A. Anderson, J. C. Tannehill, and R. H. Pletcher.
Computational Fluid Mechanics and Heat Transfer. New
York: Hemisphere, 1984. 7. W. J. Minkowycz, E. M. Sparrow, G. E. Schneider, and
R. H. Pletcher. Handbook of Numerical Heat Transfer.
New York: John Wiley & Sons, 1988. 2. C. A. Brebbia. The Boundary Element Method for
Engineers. New York: Halsted Press, 1978. 8. G. E. Myers. Analytical Methods in Conduction Heat
Transfer. New York: McGrawHill, 1971. 3. G. E. Forsythe and W. R. Wasow. Finite Difference
Methods for Partial Differential Equations. New York:
John Wiley & Sons, 1960. 9. D. H. Norrie and G. DeVries. An Introduction to Finite
Element Analysis. New York: Academic Press, 1978. 4. B. Gebhart. Heat Conduction and Mass Diffusion. New
York: McGrawHill, 1993.
5. K. H. Huebner and E. A. Thornton. The Finite Element
Method for Engineers. 2nd ed. New York: John Wiley &
Sons, 1982. 10. M. N. Özisik. Finite Difference Methods in Heat Transfer.
¸
Boca Raton, FL: CRC Press, 1994.
11. S. V. Patankhar. Numerical Heat Transfer and Fluid Flow.
New York: Hemisphere, 1980.
12. T. M. Shih. Numerical Heat Transfer. New York:
Hemisphere, 1984. 6. Y. Jaluria and K. E. Torrance. Computational Heat
Transfer. New York: Hemisphere, 1986. PROBLEMS*
Why Numerical Methods?
5–1C What are the limitations of the analytical solution
methods?
5–2C How do numerical solution methods differ from analytical ones? What are the advantages and disadvantages of
numerical and analytical methods?
5–3C What is the basis of the energy balance method? How
does it differ from the formal finite difference method? For a
specified nodal network, will these two methods result in the
same or a different set of equations?
5–4C Consider a heat conduction problem that can be solved
both analytically, by solving the governing differential equation and applying the boundary conditions, and numerically,
by a software package available on your computer. Which *Problems designated by a “C” are concept questions, and
students are encouraged to answer them all. Problems designated
by an “E” are in English units, and the SI users can ignore them.
Problems with an EESCD icon
are solved using EES, and
complete solutions together with parametric studies are included
on the enclosed CD. Problems with a computerEES icon
are
comprehensive in nature, and are intended to be solved with a
computer, preferably using the EES software that accompanies
this text. approach would you use to solve this problem? Explain your
reasoning.
5–5C Two engineers are to solve an actual heat transfer
problem in a manufacturing facility. Engineer A makes the necessary simplifying assumptions and solves the problem
analytically, while engineer B solves it numerically using a
powerful software package. Engineer A claims he solved
the problem exactly and thus his results are better, while engineer B claims that he used a more realistic model and thus his
results are better. To resolve the dispute, you are asked to solve
the problem experimentally in a lab. Which engineer do you
think the experiments will prove right? Explain. Finite Difference Formulation of Differential Equations
5–6C Define these terms used in the finite difference formulation: node, nodal network, volume element, nodal spacing,
and difference equation.
5–7 Consider three consecutive nodes n 1, n, and n 1 in
a plane wall. Using the finite difference form of the first derivative at the midpoints, show that the finite difference form of
the second derivative can be expressed as
Tn 1 2Tn
x2 Tn 1 0 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 315 315
CHAPTER 5
T(x) nodes 0, 1, 2, 3, 4, and 5 with a uniform nodal spacing of x.
Using the finite difference form of the first derivative (not the
energy balance approach), obtain the finite difference formulation of the boundary nodes for the case of insulation at the left
boundary (node 0) and radiation at the right boundary (node 5)
with an emissivity of and surrounding temperature of Tsurr. Tn + 1 Tn OneDimensional Steady Heat Conduction Tn – 1 5–11C Explain how the finite difference form of a heat conduction problem is obtained by the energy balance method.
5–12C In the energy balance formulation of the finite difference method, it is recommended that all heat transfer at the
boundaries of the volume element be assumed to be into the
volume element even for steady heat conduction. Is this a valid
recommendation even though it seems to violate the conservation of energy principle? ∆x ∆x
n–1 n n+1 x FIGURE P5–7
5–8 The finite difference formulation of steady twodimensional heat conduction in a medium with heat generation
and constant thermal conductivity is given by
Tm 1, n 2Tm, n Tm 1, n Tm, n 1 2Tm, n 2 Tm, n 1 2 x y ·
gm, n
k 0 in rectangular coordinates. Modify this relation for the threedimensional case.
5–9 Consider steady onedimensional heat conduction in a
plane wall with variable heat generation and constant thermal
conductivity. The nodal network of the medium consists of
nodes 0, 1, 2, 3, and 4 with a uniform nodal spacing of x.
Using the finite difference form of the first derivative (not the
energy balance approach), obtain the finite difference formulation of the boundary nodes for the case of uniform heat flux q·0
at the left boundary (node 0) and convection at the right boundary (node 4) with a convection coefficient of h and an ambient
temperature of T .
5–10 Consider steady onedimensional heat conduction in a
plane wall with variable heat generation and constant thermal
conductivity. The nodal network of the medium consists of 5–13C How is an insulated boundary handled in the finite
difference formulation of a problem? How does a symmetry
line differ from an insulated boundary in the finite difference
formulation?
5–14C How can a node on an insulated boundary be treated
as an interior node in the finite difference formulation of a
plane wall? Explain.
5–15C Consider a medium in which the finite difference
formulation of a general interior node is given in its simplest
form as
Tm (a)
(b)
(c)
(d)
(e) 2Tm
x2 1 Tm 1 ·
gm
k 0 Is heat transfer in this medium steady or transient?
Is heat transfer one, two, or threedimensional?
Is there heat generation in the medium?
Is the nodal spacing constant or variable?
Is the thermal conductivity of the medium constant or
variable? 5–16 Consider steady heat conduction in a plane wall whose
left surface (node 0) is maintained at 30°C while the right surface (node 8) is subjected to a heat flux of 800 W/m2. Express
the finite difference formulation of the boundary nodes 0 and 8
30°C Insulation
Tsurr
Radiation 1 2 FIGURE P5–10 3 4 ∆x ε ∆x
0 W
–
800 —2
m No heat generation ·
g(x) 01 2 3 5 FIGURE P5–16 4 5 6 78 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 316 316
HEAT TRANSFER for the case of no heat generation. Also obtain the finite difference formulation for the rate of heat transfer at the left
boundary.
5–17 Consider steady heat conduction in a plane wall with
variable heat generation and constant thermal conductivity.
The nodal network of the medium consists of nodes 0, 1, 2, 3,
and 4 with a uniform nodal spacing of x. Using the energy
balance approach, obtain the finite difference formulation of
the boundary nodes for the case of uniform heat flux q·0 at the
left boundary (node 0) and convection at the right boundary
(node 4) with a convection coefficient of h and an ambient
temperature of T .
5–18 Consider steady onedimensional heat conduction in a
plane wall with variable heat generation and constant thermal
conductivity. The nodal network of the medium consists of
nodes 0, 1, 2, 3, 4, and 5 with a uniform nodal spacing of x.
Using the energy balance approach, obtain the finite difference
formulation of the boundary nodes for the case of insulation at
the left boundary (node 0) and radiation at the right boundary
(node 5) with an emissivity of and surrounding temperature
of Tsurr.
5–19 Consider steady onedimensional heat conduction in a
plane wall with variable heat generation and constant thermal
conductivity. The nodal network of the medium consists of
nodes 0, 1, 2, 3, 4, and 5 with a uniform nodal spacing of x.
The temperature at the right boundary (node 5) is specified.
Using the energy balance approach, obtain the finite difference
formulation of the boundary node 0 on the left boundary for the
case of combined convection, radiation, and heat flux at the left
boundary with an emissivity of , convection coefficient of h,
ambient temperature of T , surrounding temperature of Tsurr,
and uniform heat flux of q·0. Also, obtain the finite difference
formulation for the rate of heat transfer at the right boundary.
·
q0 Tsurr ε
Convection T5 ·
g(x) Radiation
∆x
0 1 2 3 45 boundary (node 0) and radiation at the right boundary (node 2)
with an emissivity of and surrounding temperature of Tsurr.
5–21 Consider steady onedimensional heat conduction in a
plane wall with variable heat generation and variable thermal
conductivity. The nodal network of the medium consists of
nodes 0, 1, and 2 with a uniform nodal spacing of x. Using
the energy balance approach, obtain the finite difference formulation of this problem for the case of specified heat flux q·0
to the wall and convection at the left boundary (node 0) with a
convection coefficient of h and ambient temperature of T , and
radiation at the right boundary (node 2) with an emissivity of
and surrounding surface temperature of Tsurr.
Tsurr ·
q0
·
g(x)
k(T) Convection
h
T Radiation ε ∆x
0 1 2 FIGURE P5–21
5–22 Consider steady onedimensional heat conduction in a
pin fin of constant diameter D with constant thermal conductivity. The fin is losing heat by convection to the ambient air at
T with a heat transfer coefficient of h. The nodal network of
the fin consists of nodes 0 (at the base), 1 (in the middle), and 2
(at the fin tip) with a uniform nodal spacing of x. Using the
energy balance approach, obtain the finite difference formulation of this problem to determine T1 and T2 for the case of specified temperature at the fin base and negligible heat transfer at
the fin tip. All temperatures are in °C.
5–23 Consider steady onedimensional heat conduction in a
pin fin of constant diameter D with constant thermal conductivity. The fin is losing heat by convection to the ambient air
at T with a convection coefficient of h, and by radiation to
the surrounding surfaces at an average temperature of Tsurr. h, T
Tsurr FIGURE P5–19
5–20 Consider steady onedimensional heat conduction in a
composite plane wall consisting of two layers A and B in perfect contact at the interface. The wall involves no heat generation. The nodal network of the medium consists of nodes 0, 1
(at the interface), and 2 with a uniform nodal spacing of x.
Using the energy balance approach, obtain the finite difference
formulation of this problem for the case of insulation at the left T0 Radiation ε
0
∆x 1 h, T
Convection FIGURE P5–23 D
2 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 317 317
CHAPTER 5 The nodal network of the fin consists of nodes 0 (at the base),
1 (in the middle), and 2 (at the fin tip) with a uniform nodal
spacing of x. Using the energy balance approach, obtain the
finite difference formulation of this problem to determine
T1 and T2 for the case of specified temperature at the fin base
and negligible heat transfer at the fin tip. All temperatures
are in °C.
5–24 Consider a large uranium plate of thickness 5 cm and
thermal conductivity k 28 W/m · °C in which heat is gener·
ated uniformly at a constant rate of g 6 105 W/m3. One
side of the plate is insulated while the other side is subjected
to convection to an environment at 30°C with a heat transfer
coefficient of h
60 W/m2 · °C. Considering six equally
spaced nodes with a nodal spacing of 1 cm, (a) obtain the finite
difference formulation of this problem and (b) determine the
nodal temperatures under steady conditions by solving those
equations.
5–25 Consider an aluminum alloy fin (k 180 W/m · °C) of
triangular cross section whose length is L 5 cm, base thickness is b 1 cm, and width w in the direction normal to the
plane of paper is very large. The base of the fin is maintained
at a temperature of T0 180°C. The fin is losing heat by con25°C with a heat transfer
vection to the ambient air at T
coefficient of h 25 W/m2 · °C and by radiation to the surrounding surfaces at an average temperature of Tsurr 290 K.
Using the finite difference method with six equally spaced
nodes along the fin in the xdirection, determine (a) the temperatures at the nodes and (b) the rate of heat transfer from the
fin for w 1 m. Take the emissivity of the fin surface to be 0.9
and assume steady onedimensional heat transfer in the fin. from 100°C to 200°C. Plot the fin tip temperature and the rate
of heat transfer as a function of the fin base temperature, and
discuss the results.
5–27 Consider a large plane wall of thickness L
0.4 m,
thermal conductivity k
2.3 W/m · °C, and surface area
A 20 m2. The left side of the wall is maintained at a constant
temperature of 80°C, while the right side loses heat by con15°C with a heat transvection to the surrounding air at T
fer coefficient of h
24 W/m2 · °C. Assuming steady onedimensional heat transfer and taking the nodal spacing to be
10 cm, (a) obtain the finite difference formulation for all nodes,
(b) determine the nodal temperatures by solving those equations, and (c) evaluate the rate of heat transfer through the wall.
5–28 Consider the base plate of a 800W household iron having a thickness of L 0.6 cm, base area of A 160 cm2, and
thermal conductivity of k 20 W/m · °C. The inner surface of
the base plate is subjected to uniform heat flux generated by
the resistance heaters inside. When steady operating conditions
are reached, the outer surface temperature of the plate is measured to be 85°C. Disregarding any heat loss through the upper
part of the iron and taking the nodal spacing to be 0.2 cm,
(a) obtain the finite difference formulation for the nodes and
(b) determine the inner surface temperature of the plate by
Answer: (b) 100°C
solving those equations.
Insulation
Resistance heater, 800 W ∆ x = 0.2 cm Triangular fin w 0 h, T 1 2 3 x 160 cm2 T0 FIGURE P5–28 b
0 b/2
tan θ = —–
L θ
1 2 3 ∆x 4 5 x L FIGURE P5–25
5–26 85°C Base plate Reconsider Problem 5–25. Using EES (or other)
software, investigate the effect of the fin base
temperature on the fin tip temperature and the rate of heat
transfer from the fin. Let the temperature at the fin base vary 5–29 Consider a large plane wall of thickness L
0.3 m,
thermal conductivity k
2.5 W/m · °C, and surface area
A 12 m2. The left side of the wall is subjected to a heat flux
·
of q 0 700 W/m2 while the temperature at that surface is measured to be T0 60°C. Assuming steady onedimensional heat
transfer and taking the nodal spacing to be 6 cm, (a) obtain the
finite difference formulation for the six nodes and (b) determine the temperature of the other surface of the wall by solving those equations.
5–30E A large steel plate having a thickness of L
5 in.,
thermal conductivity of k 7.2 Btu/h · ft · °F, and an emissivity of
0.6 is lying on the ground. The exposed surface of cen58933_ch05.qxd 9/4/2002 11:42 AM Page 318 318
HEAT TRANSFER
Tsurr = 295 K Tsky
Radiation ε 6
T = 25°C
h = 13 W/m2·°C Convection
h, T
0
1
2
Plate
3
4
5 1 in. 5
4
18 cm 3
2
1 3 cm 0 0.6 ft 6
Soil 95° 7
8 FIGURE P5–32 9
10 5–34 x FIGURE P5–30E
the plate exchanges heat by convection with the ambient air
80°F with an average heat transfer coefficient of
at T
h 3.5 Btu/h · ft2 · °F as well as by radiation with the open sky
at an equivalent sky temperature of Tsky 510 R. The ground
temperature below a certain depth (say, 3 ft) is not affected by
the weather conditions outside and remains fairly constant at
50°F at that location. The thermal conductivity of the soil can
be taken to be ksoil 0.49 Btu/h · ft · °F, and the steel plate can
be assumed to be in perfect contact with the ground. Assuming
steady onedimensional heat transfer and taking the nodal
spacings to be 1 in. in the plate and 0.6 ft in the ground, (a) obtain the finite difference formulation for all 11 nodes shown in
Figure P5–30E and (b) determine the top and bottom surface
temperatures of the plate by solving those equations.
5–31E Repeat Problem 5–30E by disregarding radiation heat
transfer from the upper surface.
Answers: (b) 78.7°F, 78.4°F
5–32 Consider a stainless steel spoon (k
15.1 W/m · C,
0.6) that is partially immersed in boiling water at 95°C in
a kitchen at 25°C. The handle of the spoon has a cross section
of about 0.2 cm 1 cm and extends 18 cm in the air from the
free surface of the water. The spoon loses heat by convection
to the ambient air with an average heat transfer coefficient of
h 13 W/m2 · °C as well as by radiation to the surrounding
surfaces at an average temperature of Tsurr 295 K. Assuming
steady onedimensional heat transfer along the spoon and taking the nodal spacing to be 3 cm, (a) obtain the finite difference
formulation for all nodes, (b) determine the temperature of the
tip of the spoon by solving those equations, and (c) determine
the rate of heat transfer from the exposed surfaces of the spoon.
5–33 Repeat Problem 5–32 using a nodal spacing of 1.5 cm. Reconsider Problem 5–33. Using EES (or other)
software, investigate the effects of the thermal
conductivity and the emissivity of the spoon material on the
temperature at the spoon tip and the rate of heat transfer from
the exposed surfaces of the spoon. Let the thermal conductivity vary from 10 W/m · °C to 400 W/m · °C, and the emissivity
from 0.1 to 1.0. Plot the spoon tip temperature and the heat
transfer rate as functions of thermal conductivity and emissivity, and discuss the results.
5–35 One side of a 2mhigh and 3mwide vertical plate
at 130°C is to be cooled by attaching aluminum fins (k
237 W/m · °C) of rectangular profile in an environment at
35°C. The fins are 2 cm long, 0.3 cm thick, and 0.4 cm apart.
The heat transfer coefficient between the fins and the surrounding air for combined convection and radiation is estimated to be 30 W/m2 · °C. Assuming steady onedimensional
heat transfer along the fin and taking the nodal spacing to be
0.5 cm, determine (a) the finite difference formulation of this
problem, (b) the nodal temperatures along the fin by solving
these equations, (c) the rate of heat transfer from a single fin,
130°C T = 35°C
0.4
cm 0.5 cm
0.3
cm 0
1 2 3 x
4 3m
2 cm FIGURE P5–35 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 319 319
CHAPTER 5 and (d) the rate of heat transfer from the entire finned surface
of the plate. 10 cm
9.2 cm 5–36 A hot surface at 100°C is to be cooled by attaching 3cmlong, 0.25cmdiameter aluminum pin fins (k
237 W/m · °C) with a centertocenter distance of 0.6 cm. The
temperature of the surrounding medium is 30°C, and the combined heat transfer coefficient on the surfaces is 35 W/m2 · °C.
Assuming steady onedimensional heat transfer along the fin
and taking the nodal spacing to be 0.5 cm, determine (a) the finite difference formulation of this problem, (b) the nodal temperatures along the fin by solving these equations, (c) the rate
of heat transfer from a single fin, and (d) the rate of heat transfer from a 1m 1m section of the plate. Tair = 8°C 1 cm
1 cm 20 cm
Steam
200°C 3 cm FIGURE P5–38 0.6
cm 5–39 Reconsider Problem 5–38. Using EES (or other)
software, investigate the effects of the steam temperature and the outer heat transfer coefficient on the flange tip
temperature and the rate of heat transfer from the exposed surfaces of the flange. Let the steam temperature vary from 150°C
to 300°C and the heat transfer coefficient from 15 W/m2 · °C to
60 W/m2 · °C. Plot the flange tip temperature and the heat
transfer rate as functions of steam temperature and heat transfer coefficient, and discuss the results. 0.25
cm 100°C 5–40 0.5 cm
0 1 2 3 4 5 6 x Using EES (or other) software, solve these systems of algebraic equations. (a) 3x1 x2 3x3
x1 2x2 x3
2x1 x2 x3 FIGURE P5–36
5–37 Repeat Problem 5–36 using copper fins (k
W/m · °C) instead of aluminum ones. 386 (b) 4x1 Answers: (b) 98.6°C, 97.5°C, 96.7°C, 96.0°C, 95.7°C, 95.5°C 5–38 Two 3mlong and 0.4cmthick cast iron (k
52
W/m · °C,
0.8) steam pipes of outer diameter 10 cm are
connected to each other through two 1cmthick flanges of
outer diameter 20 cm, as shown in the figure. The steam flows
inside the pipe at an average temperature of 200°C with a heat
transfer coefficient of 180 W/m2 · °C. The outer surface of the
pipe is exposed to convection with ambient air at 8°C with a
heat transfer coefficient of 25 W/m2 · °C as well as radiation
with the surrounding surfaces at an average temperature of
Tsurr 290 K. Assuming steady onedimensional heat conduction along the flanges and taking the nodal spacing to be 1 cm
along the flange (a) obtain the finite difference formulation for
all nodes, (b) determine the temperature at the tip of the flange
by solving those equations, and (c) determine the rate of heat
transfer from the exposed surfaces of the flange. x3 Answers: (a) x1
1.62 5–41
(a) 2x2 0.5x3
2
x3 x2 x3
1
x1 x2 x3 2, x2 3, x3 2
11.964
3
1, (b) x1 2.33, x2 2.29, Using EES (or other) software, solve these systems of algebraic equations.
3x1
2x1 (b) 0
3
2 2x2 x3 x4
x1 2x2 x4
x2 3x3 x4
3x2 x3 4x4 6 3x1 x2
2
x2 3x2
1
2x1 x4
2 8 2x3
2x3
4x3 3
2
6 6.293
12 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 320 320
HEAT TRANSFER 5–42 Using EES (or other) software, solve these systems of algebraic equations.
4x1 x2 2x3 x4
x1 3x2 x3 4x4
x1 2x2 5x4
2x2 4x3 3x4 (a) (b) x4
2 2x1
x2
1 4x2
x1
3x1 2x3 x4
2x2 2x4
3
x4 5x3
2
x2 8x4
3 Insulation 6
1 ·
g = 6 × 106 W/ m3
260 5
1
3
10
15 5 cm
5 cm (a)
(b)
(c)
(d)
(e) Tright Tbottom 4Tnode ·
gnodel 2
k 0 Is heat transfer in this medium steady or transient?
Is heat transfer one, two, or threedimensional?
Is there heat generation in the medium?
Is the nodal spacing constant or variable?
Is the thermal conductivity of the medium constant or
variable? 5–44C Consider a medium in which the finite difference
formulation of a general interior node is given in its simplest
form as
Tnode
(a)
(b)
(c)
(d)
(e) (Tleft Ttop Tright 2 325 Convection
T = 20°C, h = 50 W/ m2 · °C 5–43C Consider a medium in which the finite difference
formulation of a general interior node is given in its simplest
form as
Ttop 1 200°C TwoDimensional Steady Heat Conduction Tleft 290 240 350°C 305 3 5 Tbottom)/4 Is heat transfer in this medium steady or transient?
Is heat transfer one, two, or threedimensional?
Is there heat generation in the medium?
Is the nodal spacing constant or variable?
Is the thermal conductivity of the medium constant or
variable? 5–45C What is an irregular boundary? What is a practical
way of handling irregular boundary surfaces with the finite difference method?
5–46 Consider steady twodimensional heat transfer in a long
solid body whose cross section is given in the figure. The temperatures at the selected nodes and the thermal conditions at
the boundaries are as shown. The thermal conductivity of the
body is k 45 W/m · °C, and heat is generated in the body uni·
formly at a rate of g 6 106 W/m3. Using the finite difference method with a mesh size of x
y 5.0 cm, determine
(a) the temperatures at nodes 1, 2, and 3 and (b) the rate of heat
loss from the bottom surface through a 1mlong section of the
body. FIGURE P5–46
5–47 Consider steady twodimensional heat transfer in a long
solid body whose cross section is given in the figure. The measured temperatures at selected points of the outer surfaces are
as shown. The thermal conductivity of the body is k
45
W/m · °C, and there is no heat generation. Using the finite difference method with a mesh size of x
y 2.0 cm, determine the temperatures at the indicated points in the medium.
Hint: Take advantage of symmetry.
150 180 200 150°C 1 180 2 3 180 180
4 5 6 200 200
7 8 9 180 180
2 cm
2 cm 150 180 200 180 150 FIGURE P5–47
5–48 Consider steady twodimensional heat transfer in a long
solid bar whose cross section is given in the figure. The measured temperatures at selected points of the outer surfaces are
as shown. The thermal conductivity of the body is k
20
W/m · °C, and there is no heat generation. Using the finite difference method with a mesh size of x
y 1.0 cm, determine the temperatures at the indicated points in the medium.
Answers: T1 185°C, T2 T3 T4 190°C 5–49 Starting with an energy balance on a volume element,
obtain the steady twodimensional finite difference equation
for a general interior node in rectangular coordinates for T(x, y)
for the case of variable thermal conductivity and uniform heat
generation. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 321 321
CHAPTER 5
200°C W/m3. Plot the temperatures at nodes 1 and 3, and the rate of
heat loss as functions of thermal conductivity and heat generation rate, and discuss the results. 1 2 3 4 180 5–52 Consider steady twodimensional heat transfer in a long
solid bar whose cross section is given in the figure. The measured temperatures at selected points on the outer surfaces are
as shown. The thermal conductivity of the body is k
20
W/m · °C, and there is no heat generation. Using the finite difference method with a mesh size of x
y 1.0 cm, determine the temperatures at the indicated points in the medium.
Hint: Take advantage of symmetry. 200
Insulation
(a)
120 100 120 1 100°C Answers: (b) T1 143°C, T2 T4 100°C 2 120 x 120 y 150
3 140 1 200 FIGURE P5–48 100 4 6 300 (a) Insulation 100 100 1 100 2 100 3 100°C 4 1 cm
1 cm
200 200 200 200 200°C (b) ·
g = 107 W/ m3 FIGURE P5–52 1 2 120 120 3 4 150 150
0.1 m
0.1 m Convection
200 200 1 Reconsider Problem 5–50. Using EES (or other)
software, investigate the effects of the thermal
conductivity and the heat generation rate on the temperatures at
nodes 1 and 3, and the rate of heat loss from the top surface.
Let the thermal conductivity vary from 10 W/m · °C to 400
W/m · °C and the heat generation rate from 105 W/m3 to 108 ·
qL 2 4 h, T 3 5 6 1.5 cm
1.5 cm
120°C FIGURE P5–53 7 8
Insulation 200 5–53 Consider steady twodimensional heat transfer in an
Lshaped solid body whose cross section is given in the figure.
The thermal conductivity of the body is k 45 W/m · °C, and
·
heat is generated in the body at a rate of g 5 106 W/m3. FIGURE P5–50
5–51 250 5
Insulation 100°C 100 3 1 cm 300 5–50 Consider steady twodimensional heat transfer in a long
solid body whose cross section is given in the figure. The temperatures at the selected nodes and the thermal conditions on
the boundaries are as shown. The thermal conductivity of the
body is k 180 W/m · °C, and heat is generated in the body
·
uniformly at a rate of g 107 W/m3. Using the finite difference
method with a mesh size of x
y
10 cm, determine
(a) the temperatures at nodes 1, 2, 3, and 4 and (b) the rate
of heat loss from the top surface through a 1mlong section of
the body. 200 2 250 (b) 200 1 cm Insulation 100 150 4 140 136°C T3 Insulation 180 Insulation 150 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 322 322
HEAT TRANSFER The right surface of the body is insulated, and the bottom surface is maintained at a uniform temperature of 120°C. The
entire top surface is subjected to convection with ambient air
30°C with a heat transfer coefficient of h
55
at T
W/m2 · °C, and the left surface is subjected to heat flux at a uniform rate of q·L 8000 W/m2. The nodal network of the problem consists of 13 equally spaced nodes with x
y
1.5 cm. Five of the nodes are at the bottom surface and thus
their temperatures are known. (a) Obtain the finite difference
equations at the remaining eight nodes and (b) determine the
nodal temperatures by solving those equations.
5–54E Consider steady twodimensional heat transfer in a
long solid bar of square cross section in which heat is gener·
ated uniformly at a rate of g 0.19 105 Btu/h · ft3. The cross
section of the bar is 0.4 ft 0.4 ft in size, and its thermal conductivity is k 16 Btu/h · ft · °F. All four sides of the bar are
70°F with
subjected to convection with the ambient air at T
a heat transfer coefficient of h 7.9 Btu/h · ft2 · °F. Using the
finite difference method with a mesh size of x
y 0.2 ft,
determine (a) the temperatures at the nine nodes and (b) the
rate of heat loss from the bar through a 1ftlong section.
Answer: (b) 3040 Btu/h
h, T
1 2 3 ·
g
4 5 7 h, T 6 8 h, T 9 h, T FIGURE P5–54E
5–55 Hot combustion gases of a furnace are flowing through
a concrete chimney (k 1.4 W/m · °C) of rectangular cross
Tsky section. The flow section of the chimney is 20 cm 40 cm,
and the thickness of the wall is 10 cm. The average temperature
of the hot gases in the chimney is Ti 280°C, and the average
convection heat transfer coefficient inside the chimney is hi
75 W/m2 · °C. The chimney is losing heat from its outer surface
15°C by convection with a heat
to the ambient air at To
18 W/m2 · °C and to the sky by
transfer coefficient of ho
radiation. The emissivity of the outer surface of the wall is
0.9, and the effective sky temperature is estimated to be
250 K. Using the finite difference method with x
y
10 cm and taking full advantage of symmetry, (a) obtain the
finite difference formulation of this problem for steady twodimensional heat transfer, (b) determine the temperatures at the
nodal points of a cross section, and (c) evaluate the rate of heat
loss for a 1mlong section of the chimney.
5–56 Repeat Problem 5–55 by disregarding radiation heat
transfer from the outer surfaces of the chimney.
5–57 Reconsider Problem 5–55. Using EES (or other)
software, investigate the effects of hotgas temperature and the outer surface emissivity on the temperatures at
the outer corner of the wall and the middle of the inner surface
of the right wall, and the rate of heat loss. Let the temperature
of the hot gases vary from 200°C to 400°C and the emissivity
from 0.1 to 1.0. Plot the temperatures and the rate of heat loss
as functions of the temperature of the hot gases and the emissivity, and discuss the results.
5–58 Consider a long concrete dam (k 0.6 W/m · °C,
0.7 m2/s) of triangular cross section whose
s
exposed surface is subjected to solar heat flux of q·s
800 W/m2 and to convection and radiation to the environment at 25°C with a combined heat transfer coefficient of 30
W/m2 · °C. The 2mhigh vertical section of the dam is subjected to convection by water at 15°C with a heat transfer
coefficient of 150 W/m2 · °C, and heat transfer through the
2mlong base is considered to be negligible. Using the finite
difference method with a mesh size of x
y
1 m and
assuming steady twodimensional heat transfer, determine the
temperature of the top, middle, and bottom of the exposed surAnswers: 21.3°C, 43.2°C, 43.6°C
face of the dam. ε To, ho
Chimney 1
·
qs 10 cm
h
T Hot gases
20 cm
Ti, hi 1m
2
3 Water 10 cm 1m
4 10 cm FIGURE P5–55 40 cm 1m 5 1m
6 10 cm FIGURE P5–58 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 323 323
CHAPTER 5 5–59E Consider steady twodimensional heat transfer in a
Vgrooved solid body whose cross section is given in the figure. The top surfaces of the groove are maintained at 32°F
while the bottom surface is maintained at 212°F. The side surfaces of the groove are insulated. Using the finite difference
method with a mesh size of x
y 1 ft and taking advantage of symmetry, determine the temperatures at the middle of
the insulated surfaces. 32°F 2 4 3 5 9 7 6 10 8 Insulation Insulation 1 11 1 ft
1 ft dimensional heat transfer and (b) determine the unknown nodal
temperatures by solving those equations.
Answers: (b) 85.7°C, 86.4°C, 87.6°C 5–62 Consider a 5mlong constantan block (k
23
W/m · °C) 30 cm high and 50 cm wide. The block is completely submerged in iced water at 0°C that is well stirred, and
the heat transfer coefficient is so high that the temperatures on
both sides of the block can be taken to be 0°C. The bottom surface of the bar is covered with a lowconductivity material so
that heat transfer through the bottom surface is negligible. The
top surface of the block is heated uniformly by a 6kW resistance heater. Using the finite difference method with a mesh
size of x
y 10 cm and taking advantage of symmetry,
(a) obtain the finite difference formulation of this problem for
steady twodimensional heat transfer, (b) determine the unknown nodal temperatures by solving those equations, and
(c) determine the rate of heat transfer from the block to the iced
water.
6 kW heater Insulation 212°F FIGURE P5–59E Reconsider Problem 5–59E. Using EES (or
other) software, investigate the effects of the
temperatures at the top and bottom surfaces on the temperature
in the middle of the insulated surface. Let the temperatures at
the top and bottom surfaces vary from 32°F to 212°F. Plot the
temperature in the middle of the insulated surface as functions
of the temperatures at the top and bottom surfaces, and discuss
the results. 0°C 5–60 5–61 Consider a long solid bar whose thermal conductivity is
k 12 W/m · °C and whose cross section is given in the figure.
The top surface of the bar is maintained at 50°C while the bottom surface is maintained at 120°C. The left surface is insulated and the remaining three surfaces are subjected to
25°C with a heat transfer
convection with ambient air at T
coefficient of h
30 W/m2 · °C. Using the finite difference
method with a mesh size of x
y 10 cm, (a) obtain the
finite difference formulation of this problem for steady two50°C Insulation Convection
h, T 10 cm
10 cm 10 cm
10 cm
Insulation FIGURE P5–62
Transient Heat Conduction
5–63C How does the finite difference formulation of
a transient heat conduction problem differ from that of a
steady heat conduction problem? What does the term
i
i
A xC(Tm 1 Tm)/ t represent in the transient finite difference formulation?
5–64C What are the two basic methods of solution of transient problems based on finite differencing? How do heat
transfer terms in the energy balance formulation differ in the
two methods?
5–65C The explicit finite difference formulation of a general
interior node for transient heat conduction in a plane wall is
given by
i
Tm 120°C FIGURE P5–61 0°C 1 i
2Tm i
Tm 1 ·
g im x2
k T im 1 i
Tm Obtain the finite difference formulation for the steady case by
simplifying the relation above. cen58933_ch05.qxd 9/4/2002 11:43 AM Page 324 324
HEAT TRANSFER 5–66C The explicit finite difference formulation of a general
interior node for transient twodimensional heat conduction is
given by
i1
Tnode i
(Tleft (1 i
Ttop 4 i
Tright i
)Tnode ·
g(x, t )
·
q0 i
Tbottom)
·i
gnodel 2
k 5–67C Is there any limitation on the size of the time step t
in the solution of transient heat conduction problems using
(a) the explicit method and (b) the implicit method?
5–68C Express the general stability criterion for the explicit
method of solution of transient heat conduction problems.
5–69C Consider transient onedimensional heat conduction
in a plane wall that is to be solved by the explicit method. If
both sides of the wall are at specified temperatures, express the
stability criterion for this problem in its simplest form.
5–70C Consider transient onedimensional heat conduction
in a plane wall that is to be solved by the explicit method.
If both sides of the wall are subjected to specified heat
flux, express the stability criterion for this problem in its simplest form.
5–71C Consider transient twodimensional heat conduction
in a rectangular region that is to be solved by the explicit
method. If all boundaries of the region are either insulated or at
specified temperatures, express the stability criterion for this
problem in its simplest form.
5–72C The implicit method is unconditionally stable and
thus any value of time step t can be used in the solution of
transient heat conduction problems. To minimize the computation time, someone suggests using a very large value of t
since there is no danger of instability. Do you agree with this
suggestion? Explain.
5–73 Consider transient heat conduction in a plane wall
whose left surface (node 0) is maintained at 50°C while the
right surface (node 6) is subjected to a solar heat flux of 600
W/m2. The wall is initially at a uniform temperature of 50°C.
Express the explicit finite difference formulation of the boundary nodes 0 and 6 for the case of no heat generation. Also,
obtain the finite difference formulation for the total amount
of heat transfer at the left boundary during the first three
time steps.
5–74 Consider transient heat conduction in a plane wall with
variable heat generation and constant thermal conductivity.
The nodal network of the medium consists of nodes 0, 1, 2, 3,
and 4 with a uniform nodal spacing of x. The wall is initially
at a specified temperature. Using the energy balance approach,
obtain the explicit finite difference formulation of the boundary
·
nodes for the case of uniform heat flux q 0 at the left boundary ∆x
0 Obtain the finite difference formulation for the steady case by
simplifying the relation above. h, T 1 2 3 x 4 FIGURE P5–74
(node 0) and convection at the right boundary (node 4) with a
convection coefficient of h and an ambient temperature of T .
Do not simplify.
5–75
tion. Repeat Problem 5–74 for the case of implicit formula 5–76 Consider transient heat conduction in a plane wall with
variable heat generation and constant thermal conductivity.
The nodal network of the medium consists of nodes 0, 1, 2,
3, 4, and 5 with a uniform nodal spacing of x. The wall is initially at a specified temperature. Using the energy balance approach, obtain the explicit finite difference formulation of the
boundary nodes for the case of insulation at the left boundary
(node 0) and radiation at the right boundary (node 5) with an
emissivity of and surrounding temperature of Tsurr.
5–77 Consider transient heat conduction in a plane wall with
variable heat generation and constant thermal conductivity.
The nodal network of the medium consists of nodes 0, 1, 2, 3,
and 4 with a uniform nodal spacing of x. The wall is initially
at a specified temperature. The temperature at the right boundary (node 4) is specified. Using the energy balance approach,
obtain the explicit finite difference formulation of the boundary
node 0 for the case of combined convection, radiation, and heat
flux at the left boundary with an emissivity of , convection coefficient of h, ambient temperature of T , surrounding temper·
ature of Tsurr, and uniform heat flux of q 0 toward the wall. Also,
obtain the finite difference formulation for the total amount of
heat transfer at the right boundary for the first 20 time steps. Tsurr
Radiation TL ·
g(x, t )
∆x ·
q0
0
h, T
Convection FIGURE P5–77 L
1 2 3 4 x cen58933_ch05.qxd 9/4/2002 11:43 AM Page 325 325
CHAPTER 5 5–78 Starting with an energy balance on a volume element,
obtain the twodimensional transient explicit finite difference
equation for a general interior node in rectangular coordinates
for T(x, y, t) for the case of constant thermal conductivity and
no heat generation.
5–79 Starting with an energy balance on a volume element,
obtain the twodimensional transient implicit finite difference
equation for a general interior node in rectangular coordinates
for T(x, y, t) for the case of constant thermal conductivity and
no heat generation.
5–80 Starting with an energy balance on a disk volume element, derive the onedimensional transient explicit finite difference equation for a general interior node for T(z, t) in a
cylinder whose side surface is insulated for the case of constant
thermal conductivity with uniform heat generation.
5–81 Consider onedimensional transient heat conduction in
a composite plane wall that consists of two layers A and B with
perfect contact at the interface. The wall involves no heat generation and initially is at a specified temperature. The nodal
network of the medium consists of nodes 0, 1 (at the interface),
and 2 with a uniform nodal spacing of x. Using the energy
balance approach, obtain the explicit finite difference formulation of this problem for the case of insulation at the left boundary (node 0) and radiation at the right boundary (node 2) with
an emissivity of and surrounding temperature of Tsurr.
Tsurr
Insulation
A ε B Radiation
∆x
0 ∆x
1 2 5–84 Consider a large uranium plate of thickness L 8 cm,
thermal conductivity k 28 W/m · °C, and thermal diffusivity
12.5 10 6 m2/s that is initially at a uniform temperature
of 100°C. Heat is generated uniformly in the plate at a constant
·
rate of g 106 W/m3. At time t 0, the left side of the plate is
insulated while the other side is subjected to convection with
20°C with a heat transfer coefficient of
an environment at T
h 35 W/m2 · °C. Using the explicit finite difference approach
with a uniform nodal spacing of x 2 cm, determine (a) the
temperature distribution in the plate after 5 min and (b) how
long it will take for steady conditions to be reached in the plate.
5–85 Reconsider Problem 5–84. Using EES (or other)
software, investigate the effect of the cooling
time on the temperatures of the left and right sides of the plate.
Let the time vary from 5 min to 60 min. Plot the temperatures
at the left and right surfaces as a function of time, and discuss
the results.
5–86 Consider a house whose south wall consists of a 30cmthick Trombe wall whose thermal conductivity is k
0.70
W/m · °C and whose thermal diffusivity is
0.44 10 6
m2/s. The variations of the ambient temperature Tout and the
·
solar heat flux q solar incident on a southfacing vertical surface
throughout the day for a typical day in February are given in
the table in 3h intervals. The Trombe wall has single glazing
with an absorptivitytransmissivity product of
0.76 (that
is, 76 percent of the solar energy incident is absorbed by the
exposed surface of the Trombe wall), and the average combined heat transfer coefficient for heat loss from the Trombe
wall to the ambient is determined to be hout 3.4 W/m2 · °C.
20°C at all
The interior of the house is maintained at Tin
times, and the heat transfer coefficient at the interior surface of
the Trombe wall is hin 9.1 W/m2 · °C. Also, the vents on the
Trombe wall are kept closed, and thus the only heat transfer between the air in the house and the Trombe wall is through the Interface FIGURE P5–81
5–82 Consider transient onedimensional heat conduction in
a pin fin of constant diameter D with constant thermal conductivity. The fin is losing heat by convection to the ambient air at
T with a heat transfer coefficient of h and by radiation to the
surrounding surfaces at an average temperature of Tsurr. The
nodal network of the fin consists of nodes 0 (at the base), 1 (in
the middle), and 2 (at the fin tip) with a uniform nodal spacing
of x. Using the energy balance approach, obtain the explicit
finite difference formulation of this problem for the case of a
specified temperature at the fin base and negligible heat transfer at the fin tip.
5–83 Repeat Problem 5–82 for the case of implicit
formulation. Sun’s
rays
Trombe
wall
Heat
gain Heat
loss Glazing FIGURE P5–86 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 326 326
HEAT TRANSFER 5–88 TABLE P5–86
The hourly variations of the monthly average ambient
temperature and solar heat flux incident on a
vertical surface
Ambient
Temperature, °C Time of Day
7
10
1
4
7
10
1
4 AM–10 AM Solar
Insolation, W/m2 0
4
6
1
2
3
4
4 AM–1 PM
PM–4 PM
PM–7 PM
PM–10 PM
PM–1 AM
AM–4 AM
AM–7 AM 375
750
580
95
0
0
0
0 Reconsider Problem 5–87. Using EES (or other)
software, plot the temperature at the top corner as
a function of heating time varies from 2 min to 30 min, and discuss the results.
5–89 Consider a long solid bar (k 28 W/m · °C and
12 10 6 m2/s) of square cross section that is initially at a uniform temperature of 20°C. The cross section of the bar is
20 cm 20 cm in size, and heat is generated in it uniformly at
·
a rate of g 8 105 W/m3. All four sides of the bar are sub30°C with
jected to convection to the ambient air at T
a heat transfer coefficient of h
45 W/m2 · °C. Using the
explicit finite difference method with a mesh size of x
y 10 cm, determine the centerline temperature of the bar
(a) after 10 min and (b) after steady conditions are established.
h, T
1 interior surface of the wall. Assuming the temperature of the
Trombe wall to vary linearly between 20°C at the interior surface and 0°C at the exterior surface at 7 AM and using the explicit finite difference method with a uniform nodal spacing of
x 5 cm, determine the temperature distribution along the
thickness of the Trombe wall after 6, 12, 18, 24, 30, 36, 42, and
48 hours and plot the results. Also, determine the net amount of
heat transferred to the house from the Trombe wall during the
first day if the wall is 2.8 m high and 7 m long.
5–87 Consider twodimensional transient heat transfer in an
Lshaped solid bar that is initially at a uniform temperature
of 140°C and whose cross section is given in the figure. The
thermal conductivity and diffusivity of the body are k
15
W/m · °C and
3.2 10 6 m2/s, respectively, and heat is
·
generated in the body at a rate of g 2 107 W/m3. The right
surface of the body is insulated, and the bottom surface is
maintained at a uniform temperature of 140°C at all times. At
time t
0, the entire top surface is subjected to convection
25°C with a heat transfer coefficient
with ambient air at T
of h
80 W/m2 · °C, and the left surface is subjected to
·
uniform heat flux at a rate of q L 8000 W/m2. The nodal network of the problem consists of 13 equally spaced nodes with
x
y 1.5 cm. Using the explicit method, determine the
temperature at the top corner (node 3) of the body after 2, 5,
and 30 min.
Convection ·
qL 2 4 h, T 3 5 6 1.5 cm
1.5 cm 7 5
10 cm
10 cm 7 8 6 h, T 9 h, T FIGURE P5–89
5–90E Consider a house whose windows are made of
0.375in.thick glass (k 0.48 Btu/h · ft · °F and
4.2
10 6 ft2/s). Initially, the entire house, including the walls and
the windows, is at the outdoor temperature of To 35°F. It is
observed that the windows are fogged because the indoor temperature is below the dewpoint temperature of 54°F. Now the
heater is turned on and the air temperature in the house is
raised to Ti 72°F at a rate of 2°F rise per minute. The heat
transfer coefficients at the inner and outer surfaces of the wall
can be taken to be hi 1.2 and ho 2.6 Btu/h · ft2 · °F, respectively, and the outdoor temperature can be assumed to remain
constant. Using the explicit finite difference method with a
mesh size of x 0.125 in., determine how long it will take
Window
glass
∆x
12 House
Fog 140°C FIGURE P5–87 4 h, T hi
8 3 ·
g Ti Insulation 1 2 FIGURE P5–90E 34 0.375 in To
ho Outdoors cen58933_ch05.qxd 9/4/2002 11:43 AM Page 327 327
CHAPTER 5 for the fog on the windows to clear up (i.e., for the inner surface temperature of the window glass to reach 54°F).
5–91 A common annoyance in cars in winter months is the
formation of fog on the glass surfaces that blocks the view.
A practical way of solving this problem is to blow hot air or to
attach electric resistance heaters to the inner surfaces. Consider
the rear window of a car that consists of a 0.4cmthick glass
(k 0.84 W/m · °C and
0.39 10 6 m2/s). Strip heater
wires of negligible thickness are attached to the inner surface
of the glass, 4 cm apart. Each wire generates heat at a rate of
10 W/m length. Initially the entire car, including its windows,
3°C. The heat transfer
is at the outdoor temperature of To
coefficients at the inner and outer surfaces of the glass can be
taken to be hi 6 and ho 20 W/m2 · °C, respectively. Using
the explicit finite difference method with a mesh size of x
0.2 cm along the thickness and y 1 cm in the direction normal to the heater wires, determine the temperature distribution
throughout the glass 15 min after the strip heaters are turned
on. Also, determine the temperature distribution when steady
conditions are reached.
Thermal
symmetry line
Inner
surface Outer
surface
Glass Heater
10 W/ m
0.2 cm
1 cm
Thermal
symmetry line Tsky
Radiation Convection
ho, To ε 15 cm ε
Convection
hi, Ti Radiation
Ti FIGURE P5–93
are maintained at a constant temperature of 20°C during the
night, and the emissivity of both surfaces of the concrete roof
is 0.9. Considering both radiation and convection heat transfers
and using the explicit finite difference method with a time step
of t 5 min and a mesh size of x 3 cm, determine the
temperatures of the inner and outer surfaces of the roof at 6 AM.
Also, determine the average rate of heat transfer through the
roof during that night.
5–94 Consider a refrigerator whose outer dimensions are
1.80 m
0.8 m
0.7 m. The walls of the refrigerator are
constructed of 3cmthick urethane insulation (k
0.026
W/m · ° C and
0.36 10 6 m2/s) sandwiched between
two layers of sheet metal with negligible thickness. The refrigerated space is maintained at 3°C and the average heat transfer
coefficients at the inner and outer surfaces of the wall are
6 W/m2 · °C and 9 W/m2 · °C, respectively. Heat transfer
through the bottom surface of the refrigerator is negligible. The
kitchen temperature remains constant at about 25°C. Initially,
the refrigerator contains 15 kg of food items at an average
specific heat of 3.6 kJ/kg · °C. Now a malfunction occurs and
the refrigerator stops running for 6 h as a result. Assuming the FIGURE P5–91
5–92 Refrigerator
Air Repeat Problem 5–91 using the implicit method
with a time step of 1 min. 5–93 The roof of a house consists of a 15cmthick concrete
slab (k 1.4 W/m · °C and
0.69 10 6 m2/s) that is 20 m
wide and 20 m long. One evening at 6 PM, the slab is observed
to be at a uniform temperature of 18°C. The average ambient
air and the night sky temperatures for the entire night are predicted to be 6°C and 260 K, respectively. The convection heat
transfer coefficients at the inner and outer surfaces of the roof
can be taken to be hi 5 and ho 12 W/m2 · °C, respectively.
The house and the interior surfaces of the walls and the floor Concrete
roof 3°C Heat
25°C
ho hi Food
15 kg FIGURE P5–94 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 328 328
HEAT TRANSFER temperature of the contents of the refrigerator, including the
air inside, rises uniformly during this period, predict the temperature inside the refrigerator after 6 h when the repairman arrives. Use the explicit finite difference method with a
time step of t 1 min and a mesh size of x 1 cm and disregard corner effects (i.e., assume onedimensional heat transfer in the walls).
5–95 Reconsider Problem 5–94. Using EES (or other)
software, plot the temperature inside the refrigerator as a function of heating time as time varies from 1 h to
10 h, and discuss the results. Special Topic: Controlling the Numerical Error coordinates for T(x, y, z, t) for the case of constant thermal conductivity and no heat generation.
5–108 Consider steady onedimensional heat conduction in a
plane wall with variable heat generation and constant thermal
conductivity. The nodal network of the medium consists of
nodes 0, 1, 2, and 3 with a uniform nodal spacing of x. The
temperature at the left boundary (node 0) is specified. Using
the energy balance approach, obtain the finite difference formulation of boundary node 3 at the right boundary for the case
of combined convection and radiation with an emissivity of ,
convection coefficient of h, ambient temperature of T , and
surrounding temperature of Tsurr. Also, obtain the finite difference formulation for the rate of heat transfer at the left
boundary. 5–96C Why do the results obtained using a numerical
method differ from the exact results obtained analytically?
What are the causes of this difference?
5–97C What is the cause of the discretization error? How
does the global discretization error differ from the local
discretization error? Tsurr Radiation 5–98C Can the global (accumulated) discretization error be
less than the local error during a step? Explain.
5–99C How is the finite difference formulation for the first
derivative related to the Taylor series expansion of the solution
function?
5–100C Explain why the local discretization error of the finite difference method is proportional to the square of the step
size. Also explain why the global discretization error is proportional to the step size itself.
5–101C What causes the roundoff error? What kind of
calculations are most susceptible to roundoff error?
5–102C What happens to the discretization and the roundoff errors as the step size is decreased? Suggest some practical ways of reducing the
roundoff error. 5–103C 5–104C What is a practical way of checking if the roundoff
error has been significant in calculations?
5–105C What is a practical way of checking if the discretization error has been significant in calculations? ε ·
g(x) T0 h, T
Convection ∆x
0 1 2 3 FIGURE P5–108
5–109 Consider onedimensional transient heat conduction
in a plane wall with variable heat generation and variable thermal conductivity. The nodal network of the medium consists of
nodes 0, 1, and 2 with a uniform nodal spacing of x. Using
the energy balance approach, obtain the explicit finite difference formulation of this problem for the case of specified heat
·
flux q 0 and convection at the left boundary (node 0) with a convection coefficient of h and ambient temperature of T , and radiation at the right boundary (node 2) with an emissivity of
and surrounding temperature of Tsurr.
5–110 Repeat Problem 5–109 for the case of implicit
formulation.
5–111 Consider steady onedimensional heat conduction in a
pin fin of constant diameter D with constant thermal conductivity. The fin is losing heat by convection with the ambient air Review Problems
5–106 Starting with an energy balance on the volume element, obtain the steady threedimensional finite difference
equation for a general interior node in rectangular coordinates
for T(x, y, z) for the case of constant thermal conductivity and
uniform heat generation.
5–107 Starting with an energy balance on the volume element, obtain the threedimensional transient explicit finite difference equation for a general interior node in rectangular Convection
h, T ε 0
1 ∆x FIGURE P5–111 Tsurr Radiation
2 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 329 329
CHAPTER 5 at T (in °C) with a convection coefficient of h, and by radiation to the surrounding surfaces at an average temperature of
Tsurr (in K). The nodal network of the fin consists of nodes 0 (at
the base), 1 (in the middle), and 2 (at the fin tip) with a uniform
nodal spacing of x. Using the energy balance approach, obtain the finite difference formulation of this problem for the
case of a specified temperature at the fin base and convection
and radiation heat transfer at the fin tip.
5–112 Starting with an energy balance on the volume element, obtain the twodimensional transient explicit finite difference equation for a general interior node in rectangular
coordinates for T(x, y, t) for the case of constant thermal conductivity and uniform heat generation.
5–113 Starting with an energy balance on a disk volume element, derive the onedimensional transient implicit finite difference equation for a general interior node for T(z, t) in a
cylinder whose side surface is subjected to convection with a
convection coefficient of h and an ambient temperature of T
for the case of constant thermal conductivity with uniform heat
generation.
5–114E The roof of a house consists of a 5in.thick concrete
slab (k 0.81 Btu/h · ft · °F and
7.4 10 6 ft2/s) that is
45 ft wide and 55 ft long. One evening at 6 PM, the slab is observed to be at a uniform temperature of 70°F. The ambient air
temperature is predicted to be at about 50°F from 6 PM to
10 PM, 42°F from 10 PM to 2 AM, and 38°F from 2 AM to 6 AM,
while the night sky temperature is expected to be about 445 R
for the entire night. The convection heat transfer coefficients at
the inner and outer surfaces of the roof can be taken to be
hi 0.9 and ho 2.1 Btu/h · ft2 · °F, respectively. The house
and the interior surfaces of the walls and the floor are maintained at a constant temperature of 70°F during the night, and
the emissivity of both surfaces of the concrete roof is 0.9.
Considering both radiation and convection heat transfers and
using the explicit finite difference method with a mesh size of x 1 in. and a time step of t 5 min, determine the temperatures of the inner and outer surfaces of the roof at 6 AM.
Also, determine the average rate of heat transfer through the
roof during that night.
5–115 Solar radiation incident on a large body of clean water
(k 0.61 W/m · °C and
0.15 10 6 m2/s) such as a lake,
a river, or a pond is mostly absorbed by water, and the amount
of absorption varies with depth. For solar radiation incident at
a 45° angle on a 1mdeep large pond whose bottom surface is
black (zero reflectivity), for example, 2.8 percent of the solar
energy is reflected back to the atmosphere, 37.9 percent is absorbed by the bottom surface, and the remaining 59.3 percent
is absorbed by the water body. If the pond is considered to be
four layers of equal thickness (0.25 m in this case), it can be
shown that 47.3 percent of the incident solar energy is absorbed by the top layer, 6.1 percent by the upper mid layer, 3.6
percent by the lower mid layer, and 2.4 percent by the bottom
layer [for more information see Çengel and Özisik, Solar En¸
ergy, 33, no. 6 (1984), pp. 581–591]. The radiation absorbed by
the water can be treated conveniently as heat generation in the
heat transfer analysis of the pond.
Consider a large 1mdeep pond that is initially at a uniform
temperature of 15°C throughout. Solar energy is incident on
the pond surface at 45° at an average rate of 500 W/m2 for a period of 4 h. Assuming no convection currents in the water and
using the explicit finite difference method with a mesh size of
x 0.25 m and a time step of t 15 min, determine the
temperature distribution in the pond under the most favorable
conditions (i.e., no heat losses from the top or bottom surfaces
of the pond). The solar energy absorbed by the bottom surface
of the pond can be treated as a heat flux to the water at that surface in this case. Sun
Tsky
Radiation Convection
ho, To ε Solar
radiation
Concrete
roof 45° 0
1 ε
Radiation Convection
hi, Ti 2
3 Ti
4 Top layer Solar pond Upper mid layer
Lower mid layer
Bottom layer x FIGURE P5–114E ·
qs, W/ m2 FIGURE P5–115 Black cen58933_ch05.qxd 9/4/2002 11:43 AM Page 330 330
HEAT TRANSFER 5–116 Reconsider Problem 5–115. The absorption of solar
radiation in that case can be expressed more accurately as a
fourthdegree polynomial as
·
g(x)
·
q s(0.859 3.415x 6.704x2 6.339x3 2.278x4), W/m3 ·
where q s is the solar flux incident on the surface of the pond in
2
W/m and x is the distance from the free surface of the pond
in m. Solve Problem 5–115 using this relation for the absorption of solar radiation.
5–117 A hot surface at 120°C is to be cooled by attaching
8 cm long, 0.8 cm in diameter aluminum pin fins (k
237
W/m · °C and
97.1 10 6 m2/s) to it with a centertocenter distance of 1.6 cm. The temperature of the surrounding
medium is 15°C, and the heat transfer coefficient on the surfaces is 35 W/m2 · °C. Initially, the fins are at a uniform temperature of 30°C, and at time t 0, the temperature of the hot
surface is raised to 120°C. Assuming onedimensional heat
conduction along the fin and taking the nodal spacing to be
x
2 cm and a time step to be t
0.5 s, determine the
nodal temperatures after 5 min by using the explicit finite difference method. Also, determine how long it will take for
steady conditions to be reached. T0 Convection
h, T
2 cm 0
1 2
3 4 FIGURE P5–117
5–118E Consider a large plane wall of thickness L 0.3 ft
and thermal conductivity k 1.2 Btu/h · ft · °F in space. The
wall is covered with a material having an emissivity of
0.80 and a solar absorptivity of s 0.45. The inner surface of the wall is maintained at 520 R at all times, while the
outer surface is exposed to solar radiation that is incident at a
·
rate of q s 300 Btu/h · ft2. The outer surface is also losing heat
αs
T0 by radiation to deep space at 0 R. Using a uniform nodal spacing of x 0.1 ft, (a) obtain the finite difference formulation
for steady onedimensional heat conduction and (b) determine
the nodal temperatures by solving those equations.
Answers: (b) 522 R, 525 R, 527 R 5–119 Frozen food items can be defrosted by simply leaving
them on the counter, but it takes too long. The process can
be speeded up considerably for flat items such as steaks by
placing them on a large piece of highly conducting metal,
called the defrosting plate, which serves as a fin. The increased
surface area enhances heat transfer and thus reduces the defrosting time.
Consider two 1.5cmthick frozen steaks at 18°C that resemble a 15cmdiameter circular object when placed next to
each other. The steaks are now placed on a 1cmthick blackanodized circular aluminum defrosting plate (k
237
0.90) whose outer
W/m · °C,
97.1 10 6 m2/s, and
diameter is 30 cm. The properties of the frozen steaks are
1.55 kJ/kg · °C, k
1.40 W/m · °C,
970 kg/m3, Cp
0.95, and the heat of fusion is
0.93 10 6 m2/s, and
hif 187 kJ/kg. The steaks can be considered to be defrosted
when their average temperature is 0°C and all of the ice in the
steaks is melted. Initially, the defrosting plate is at the room
temperature of 20°C, and the wooden countertop it is placed on
can be treated as insulation. Also, the surrounding surfaces can
be taken to be at the same temperature as the ambient air, and
the convection heat transfer coefficient for all exposed surfaces
can be taken to be 12 W/m2 · °C. Heat transfer from the lateral
surfaces of the steaks and the defrosting plate can be neglected.
Assuming onedimensional heat conduction in both the steaks
and the defrosting plate and using the explicit finite difference
method, determine how long it will take to defrost the steaks.
Use four nodes with a nodal spacing of x 0.5 cm for the
steaks, and three nodes with a nodal spacing of r 3.75 cm
for the exposed portion of the defrosting plate. Also, use a time
step of t 5 s. Hint: First, determine the total amount of heat
transfer needed to defrost the steaks, and then determine how
long it will take to transfer that much heat.
Symmetry line
20°C
Heat
transfer ·
qs
6
∆x
0 1 2 3 ε 5 1
2
3
4 Tsurr
Defrosting
plate Radiation FIGURE P5–118E Steaks FIGURE P5–119 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 331 331
CHAPTER 5 5–120 Repeat Problem 5–119 for a copper defrosting plate
using a time step of t 3 s. Design and Essay Problems
5–121 Write a twopage essay on the finite element method,
and explain why it is used in most commercial engineering
software packages. Also explain how it compares to the finite
difference method.
5–122 Numerous professional software packages are available in the market for performing heat transfer analysis, and
they are widely advertised in professional magazines such as
the Mechanical Engineering magazine published by the American Society of Mechanical Engineers (ASME). Your company
decides to purchase such a software package and asks you
to prepare a report on the available packages, their costs, capabilities, ease of use, and compatibility with the available
hardware, and other software as well as the reputation of the
software company, their history, financial health, customer
support, training, and future prospects, among other things.
After a preliminary investigation, select the top three packages
and prepare a full report on them.
5–123 Design a defrosting plate to speed up defrosting of flat
food items such as frozen steaks and packaged vegetables and
evaluate its performance using the finite difference method
(see Prob. 5–119). Compare your design to the defrosting plates currently available on the market. The plate must perform well, and it must be suitable for purchase and use as a
household utensil, durable, easy to clean, easy to manufacture,
and affordable. The frozen food is expected to be at an initial
temperature of 18°C at the beginning of the thawing process
and 0°C at the end with all the ice melted. Specify the material,
shape, size, and thickness of the proposed plate. Justify your
recommendations by calculations. Take the ambient and surrounding surface temperatures to be 20°C and the convection
heat transfer coefficient to be 15 W/m2 · °C in your analysis.
For a typical case, determine the defrosting time with and
without the plate.
5–124 Design a fireresistant safety box whose outer dimensions are 0.5 m
0.5 m
0.5 m that will protect its combustible contents from fire which may last up to 2 h. Assume
the box will be exposed to an environment at an average temperature of 700°C with a combined heat transfer coefficient of
70 W/m2 · °C and the temperature inside the box must be below 150°C at the end of 2 h. The cavity of the box must be as
large as possible while meeting the design constraints, and the
insulation material selected must withstand the high temperatures to which it will be exposed. Cost, durability, and strength
are also important considerations in the selection of insulation
materials. cen58933_ch05.qxd 9/4/2002 11:43 AM Page 332 ...
View
Full Document
 Spring '10
 Ghaz
 Numerical Analysis, Heat Transfer, Partial differential equation, finite difference

Click to edit the document details