Heat Transfer: Fundamentals of
Heat Transfer: Fundamentals of
Heat Transfer: Fundamentals of
Heat Transfer
Theory and Applications
Michael C. Wendl
Department of Mechanical Engineering
and School of Medicine
Washington University
c
1999, 2003, 2005 by Michael C. Wendl
Non–lawyer summary: this license allows you to freely reproduce, distribute, and
display verbatim copies of this document, in electronic or print form, for non–
commercial purposes.
This document was implemented entirely with open–source, a.k.a. “free” software:
Typeset using LATEX 2 with graphics by Xfig and plots by Xmgrace, all running
on the Linux operating system (Debian distribution)
Document Version 2.1 (August 2005), see Appendix D for revision history
Wustl–E67–371
Contents
Preface vi
Chapter 1. Introduction 1
1.1. Physical Mechanisms 2
1.2. Concept of Conservation of Energy 3
1.3. The Continuum Assumption 4
1.4. Absolute Versus Relative Temperature 6
Introduction
1
1.1. PHYSICAL MECHANISMS 2
barrel
conduction
gas tube
handguard
radiation
heat shield
convection
Figure 1.1. M–16 cross section showing barrel, gas tube, heat
shield, and outer handguard.
combustion gas in the bore and friction between the bullet and the bore is
conducted through the barrel to its outer surface. The barrel is cooled by
1.2. CONCEPT OF CONSERVATION OF ENERGY 3
E gen E stored
CV
E in CS
E out
Note that the equation integrated over any time period ∆t must also hold
true:
(1.2) ∆Estored = Ein + Egen − Eout .
There will be occasions where the conservation of energy principle will
be required at a surface, for example as a boundary condition at a solid–gas,
liquid–gas, or solid–liquid interface (Fig. 1.3). The thickness of this interface
fluid
energy convected away from surface
values. In other words, they will be continuous functions of space and time. the continuum
This is the so–called continuum assumption. assumption
To illustrate this concept, consider the density ρ of a fluid. A region
of fluid is shown in Fig. 1.4(a). We desire to find the density at a point C
C (x ,o y ,oz ) o
δm ρ
δV defined
here
x δV
(a) (b)
z
Quantity Unit
length meter (m)
mass kilogram (kg)
time second (s)
temperature kelvin (K)
T1
T2
x
x x2
1
heat flow
would not allow energy to flow as shown. Since k is defined as positive, the
leading negative sign is clearly required for consistency. In this sense, the
temperature gradient is very similar to the pressure gradient in pipe flow:
movement proceeds along a negative gradient.
pure metals
alloys
non−metallic solids
liquids
gases
science deals much more comprehensively with such issues. Actual values are
available in many references (e.g. Incropera and Dewitt, 2002, Appendix A).
y
qt
q
f
δy
.
q q, E st q
l r
.
E gen
x
q
n
q δz
b
δx
z
of time: the energy generation and storage terms, Ėgen and Ėst , represented
as volumetric entities in the element, and all heat conduction terms at the
boundaries. Expressing these terms as 1–term Taylor series expansions with
respect to their associated values defined in the center, we find
∂qx δx
(2.3) q l = qx −
∂x 2
∂qx δx
(2.4) q r = qx +
∂x 2
∂qy δy
(2.5) q b = qy −
∂y 2
∂qy δy
(2.6) q t = qy +
∂y 2
∂qz δz
(2.7) q f = qz −
∂z 2
∂qz δz
(2.8) q n = qz +
∂z 2
Energy generation is simply
(2.9) Ėgen = q̇ δx δy δz ,
2.3. THE CONDUCTION EQUATION 11
where q̇ is the per unit volume rate of generation 2.7. Finally, the rate of
change of thermal energy stored in the element can be expressed as
∂T
(2.10) Ėst = ρ cp δx δy δz .
∂t
These terms can be used directly in the conceptual conservation law
in Eq. (1.1) on pp. 3, where the rate of energy entering the element is
Ėin = qb + qf + ql and the rate of energy leaving is Ėout = qn + qr + qt . This
operation yields
∂T
(2.11) ρ cp δx δy δz = q̇ δx δy δz + qb + qf + ql − (qn + qr + qt ) ,
∂t
which becomes
∂T ∂qy δy ∂qz δz
(2.12) ρ cp δx δy δz = q̇ δx δy δz + qy − + qz −
∂t ∂y 2 ∂z 2
∂qx δx ∂qz δz ∂qx δx ∂qy δy
+ qx − − qz + − qx + − qy + ,
∂x 2 ∂z 2 ∂x 2 ∂y 2
after substituting Eqs. (2.3) through (2.6) for the heat conduction terms.
We then cancel and combine terms appropriately to get
∂T ∂qx ∂qy ∂qz
(2.13) ρ cp δx δy δz = q̇ δx δy δz − δx + δy + δz .
∂t ∂x ∂y ∂z
Eq. (2.13) represents the limit in terms of what can be derived strictly
from theory. It cannot be solved, because T and q are both unknowns 2.8.
However, we have additional relationships between T and q in the form of
Fourier’s Law! In light of Eq. (2.2), we can express the heat terms as
∂T
(2.14) qx = −k δy δz
∂x
∂T
(2.15) qy = −k δx δz
∂y
∂T
(2.16) qz = −k δx δy
∂z
which can be substituted into Eq. (2.13) to obtain
∂T ∂ ∂T
(2.17) ρ cp δx δy δz = q̇ δx δy δz − −k δy δz δx
∂t ∂x ∂x
∂ ∂T ∂ ∂T
+ −k δx δz δy + −k δx δy δz .
∂y ∂y ∂z ∂z
2.7This is the unfortunate notation used in Incropera and Dewitt (2002), which can
be easily confused with heat conduction q and heat flux q 00 .
2.8We assume that heat generation q̇ would be prescribed, or would be measurable
for a problem.
2.3. THE CONDUCTION EQUATION 12
All volume terms δx×δy ×δz cancel and all double–negatives cancel, leaving
∂T ∂ ∂T ∂ ∂T ∂ ∂T
(2.18) ρ cp = q̇ + k + k + k
∂t ∂x ∂x ∂y ∂y ∂z ∂z
as the general three–dimensional time–dependent form of the heat conduc-
tion equation. We have explicitly derived this equation in the Cartesian
(rectangular) coordinate system, which will be sufficient for many of the
problems we will study. However, in some instances we will require the
cylindrical and spherical forms of the conduction equation. Most reference
texts address these cases (e.g. Incropera and Dewitt, 2002).
Eq. (2.18) expresses conservation of energy in terms of the single de-
pendent variable of temperature T . We can solve this in principle because
there is one equation and exactly one unknown. In mathematical terms,
this is a linear partial differential equation, since T does not appear in the
form of any products of itself, or its derivatives 2.9. Thus, these problems
are amenable to theoretical treatments. Non–linear systems (which do exist
in many problems) typically require advanced theory and/or computational
treatment.
You should have a conceptual understanding of what the terms in this
equation mean physically. For example, if we multiply ∂/∂x (k ∂T /∂x) in
Eq. (2.18) by dx and compare to terms in Eq. (2.11), we see
∂ ∂T
(2.19) k dx = ql00 − qr00
∂x ∂x
is essentially the net heat flux into the control volume in the x direction.
Similar interpretations hold for the corresponding terms in the y and z
directions. Thus, Eq. (2.18) specifies that the net rate of conduction into
the control volume plus the rate of generation is equal to the rate of change
of energy stored at any point in the domain.
Notice that we have not made any assumptions regarding the conduc-
tivity k in Eq. (2.18). However, if k is a constant, it can be moved outside
of the derivatives, so that
1 ∂T q̇ ∂2T ∂2T ∂2T
(2.20) = + + + ,
α ∂t k ∂x2 ∂y 2 ∂z 2
where q̇ is again the energy generation rate per unit volume (W/m 3 ) and α
is the thermal diffusivity. Mathematically, this is more straightforward to
work with than Eq. (2.18). Further simplifications are possible, for example
2.9Actually, the assumption of linearity of the overall equation also depends on the
heat generation term q̇ as being itself linear. However, if q̇ depended upon T in some
non–linear fashion, e.g. as a power of T , then the equation would no longer be linear.
We shall not discuss such complicated cases here, as the theoretical treatment is quite
difficult.
2.4. BOUNDARY AND INITIAL CONDITIONS 13
2.10Note that we have replaced the regular derivative d with the partial derivative
∂ to indicate that the temperature may be multi–dimensional, whereas we had not yet
introduced this concept in Eq. (2.1).
2.11For the moment, we will assume that h is known for problems involving Eq. (2.29).
Later, we will see that convection heat transfer is largely the study of determining h.
CHAPTER 3
Here, we will examine problems where conduction occurs along only one
primary dimension, say the x direction, and where the problem is steady.
Recall, that Eq. (2.24) on pp. 13 governs this situation, i.e.
d2 T
(3.1) = 0.
dx2
This equation is useful for conceptual understanding of conduction, but is
also a good engineering approximation for cases where there is dominance in
one dimension. For example, consider the two–dimensional domain shown
in Fig. 3.1 and assume that the temperature on the boundaries is given by
the equations
y1
x
(0, 0) x1
x
r = 2.4
400 K 310 K
b = 0.05
R
V1 V2
equation
(3.5) ∆V = I R
Looking more closely at Eqs. (3.3) and (3.4), we see that the linear drop
in temperature gradient coupled with Fourier’s Law allows us to write the
conduction heat transfer in the exact form as
dT T2 − T 1 T1 − T 2 kA
(3.6) q = −k A = −k A = kA = ∆T ,
dx L L L
where the generic problem is defined according to boundary conditions
T (0) = T1 and T (L) = T2 . A little algebra shows the equivalent relationship
L
(3.7) ∆T = q = q Rt ,
kA
whose form is identical to Eq. (3.5) where ∆T is analogous to voltage (driv-
ing potential), q is analogous to current (what “flows” in the circuit), and
Rt is a thermal resistance. More specifically, R t in Eq. (3.7) is the thermal
resistance of conduction, for which we can use the slightly more descriptive
notation
L
(3.8) Rt,cond = .
kA
3.2. CIRCUIT ANALOGY 18
In Eq. (2.29) on pp. 14, we introduced the idea of convection heat trans-
fer as quantified by a convection coefficient h. This equation is a form of
Newton’s Law of Cooling
(3.9) q = h A ( T |x=0 − T∞ ) ,
which shows that this form of convection is also linear. Again, some simple
algebra allows us to arrange this equation as
1
(3.10) T |x=0 − T∞ = q .
hA
Or, put more directly
(3.11) ∆T = q Rt,conv ,
where
1
(3.12) Rt,conv =
hA
is the thermal resistance for convection.
Are such analogs useful for practical problem solving? To answer this
question, refer back to the idea at the beginning of this chapter where we
discussed the pure one–dimensional problem versus the multi–dimensional
problem for which the size of the terms allowed us to use a one–dimensional
approximation. The former case is easily analyzed by solving Eq. (3.1) and
applying boundary conditions. But how can we treat the latter situation,
for example with respect to a problem having composite layers all of differ-
ent thermal conductivities as shown in Fig. 3.4? Here, the top and bottom
R2
2 R1 R4
T1 1 T2
4
T1 R3 T2
3
L1 L 2,3 L4
w/3
cladding
insulation
w/3
glass
w/3
t/2
L
t/2
heat transfer
• heat transfer through the glass pane, whose resistor is also given
by Eq. (3.8) as
t
Rg =
kg L w/3 cladding
R h,f
T oo Rg
T office
Ts,o T s,i
R conv R conv
R h,f
ri To
Ti
ro convection
T oo
rc
be the radius of the bore, rc the radius of the chrome plating, and r o be
the outer radius of the barrel. We assume that the plating is “perfect”, i.e.
there are no interstitial voids, so there is no thermal contact resistance. Also
let’s assume that conditions in the bore are such that the bore temperature
equals the temperature at the inner surface of the barrel. How can this
system be represented by a thermal circuit? We simply construct a series
circuit representing the thermal resistance of the chrome plating, the barrel
itself, and the convection occurring at the outer surface of the barrel, as
depicted in Fig. 3.8. Applying the circuit analogy, we find
3.4. HEAT GENERATION 23
Ti To T oo
R plating R barrel R convection
Ti − T ∞
(3.23) q = ,
Rplating + Rbarrel + Rconvection
where
ln (rc /ri )
(3.24) Rplating = ,
2 π L kchrome
ln (ro /rc )
(3.25) Rbarrel = ,
2 π L kbarrel
1
(3.26) Rconvection = .
2 π ro L h
Here, L is the length of the barrel and h is the convection coefficient at its
outer surface. We notice that for this problem, the thickness of plating is
very thin, therefore rc /ri is just slightly more than unity . The natural log I&D Ex. 3.4
of this is close to zero, therefore, we conclude that chrome–plating the barrel pp 107
does not measurably affect the heat transfer situation.
x
Ts
lower voltage
x=+L
high voltage
Ts x=−L
T oo
ri
rw
w
ire
in
su
la
tio
n
Tw
Ti
Tw Ti T oo
R insul R convec
• T |surf ace : Similar to Tf luid , the device temperature and thus the
resulting surface temperature is usually constrained to a small op-
erating range.
dAs
qx d q conv
Ac
convection
cooling around a x
fin
x dx q x + dx
x + dx
where
hP
(3.42) m2 =
k Ac
is a composite parameter made up of the physical and geometric attributes
of the problem. The general solution for Eq. (3.41) can be written in the
form of exponentials as
(3.43) θ(x) = C1 emx + C2 e−mx .
Likewise, it could also be written in terms of the hyperbolic trigonometric
functions according to the identities
eβ + e−β
(3.44) cosh β =
2
and
eβ − e−β
(3.45) sinh β = .
2
Of course, now we are back in the old position of having to use boundary
conditions for an actual problem to determine the constants of integration.
We once again need 2 boundary conditions, one at the base of the fin x = 0,
and one at the end of the fin at x = L (Fig. 3.13). The first boundary
x=0
x=L
• Fin is very long: Fins convect heat to the surrounding fluid along
their length. If we assume that the fin tip is very long, i.e. L →
∞, eventually all of the heat will be convected away, so that the
temperature at the very end must be equal to T ∞ . That is, there is
no longer any temperature gradient because there is no more heat
to be transferred. According to our original substitution θ(x) =
T (x) − T∞ , this means that θ(∞) → 0. If we substitute this into
Eq. (3.43), it is clear that the first term would be unbounded, so
we must conclude that C1 = 0. We can then use the boundary
condition at the base to find C2 = θ0 , so that the solution is
(3.47) θ(x) = θ0 e−mx .
• Negligible fin tip convection: In this case, we assume the end of the
fin is insulated, so that q = 0 at x = L. Writing this in terms of
Fourier’s Law, we obtain a boundary condition of
dθ
= 0.
dx x=L
We can take the derivative of Eq. (3.43)
dθ
= C1 m emx − C2 m e−mx
dx
and substitute this into the boundary condition to obtain the equa-
tion
C1 emL − C2 e−mL = 0 .
This equation, along with the boundary condition from x = 0 in
Eq. (3.46) allows us to solve for the constants. For example, we
have C1 = θ0 − C2 from Eq. (3.46), which can be substituted into
the above equation
0 = (θ0 − C2 ) emL − C2 e−mL
θ0 emL − C2 emL + e−mL
=
so that C2 can be solved as
emL
C2 = θ 0 .
emL + e−mL
We can substitute this right back into Eq. (3.46) to solve for C 1 as
emL emL
C1 = θ0 − θ0 mL = θ 0 1 −
e + e−mL emL + e−mL
which, with a little algebra, simplifies to
e−mL
C1 = θ 0 .
emL + e−mL
3.5. FIN ANALYSIS 32
for many real–world problems. But how do we do this? As is often the case,
there is an easy way, and a hard way.
Refer to Fig. 3.14. One of our original assumptions was that the problem
temperature
convection
q fin
f
convection
x
is steady state, so that the energy entering the fin via conduction at x = 0
must be equal to the energy convected by the fin to its surroundings over its
total area, i.e. over 0 ≤ x ≤ L. The two possibilities may already be clear
at this point.
The hard way is to integrate Newton’s Law of Cooling about the entire
surface area of the fin
Z Z
qf = h (T (x) − T∞ ) dAs = h θ(x) dAs .
As As
By definition, this gives us the entire energy transferred to the surround-
ings. While this task is somewhat straightforward for the infinite fin using
Eq. (3.47), the other cases are rather more difficult, especially since the fin
tip convection must be considered3.5. Conversely, the easy way is to sim-
ply evaluate the total heat transferred into the fin by way of Fourier’s Law
applied at the base of the fin x = 0
dθ
(3.51) qf = −k Ac .
dx x=0
Since we know the exact temperature distribution, this is a simple matter
of taking a first derivative. For example, in the case of the infinite fin, we
find
dθ
= −θ0 m e−mx ,
dx
3.5The exception is of course the case of the insulated fin tip.
3.6. FIN PERFORMANCE METRICS 34
can compute εf for a variety of cases. For example, we can insert Eq. (3.52)
for the infinite fin into Eq. (3.56) to find
√ r
θ0 h P k A c kP
(3.57) εf = = .
h A c θ0 h Ac
What does this say generally about fins?
• Fin performance is enhanced by using a high conductivity material.
Of course, this is almost obvious.
• An effective fin has a high P/Ac ratio, that is the perimeter is large
compared to the cross–sectional area. For example, assume the
fin has a square cross–section of side length b, then P = 4b and
Ac = b2 , therefore P/Ac = 4/b, which means the ratio, thus the
effectiveness improves as b becomes smaller. Therefore, lots of thin,
closely–spaced fins promote heat transfer, with the provision that
the space between them is still enough to maintain good convection.
3.6Recall from Fig. 3.14 that the gradient along a real fin lessens the heat transfer.
3.6. FIN PERFORMANCE METRICS 36
We can derive a relationship between η f for a single fin and the overall
efficiency for N fins as follows. First, the total area A total is simply the sum
of the areas of all the fins plus the remaining overall area of the unfinned
base Ab . For N fins this is
(3.60) Atotal = N Af + Ab .
According to Eq. (3.58), the heat transfer from one fin is η f hAf θ0 , so that
the heat transfer from N fins is N ηf hAf θ0 . We include the heat transfer
from the remaining unfinned base hAb θ0 , so that the total heat transfer is
(3.61) qtotal = N ηf hAf θ0 + hAb θ0 .
Now, assuming that h is taken as a single constant value for both the finned
surface and the remaining area of the base, we can substitute Eqs. (3.60)
and (3.61) into Eq. (3.59) to find
N ηf hAf θ0 + hAb θ0 N η f Af + A b
(3.62) η0 = = ·
h (N Af + Ab ) θ0 N Af + Ab
This equation is sometimes given in a slightly different form that can be
derived by adding and subtracting 1 from η f as
N (ηf + 1 − 1) Af + Ab [N Af + Ab ] + (ηf − 1) N Af
η0 = = ,
N Af + Ab N Af + Ab
which simplifies to I&D Ex. 3.9
Atotal − (1 − ηf ) N Af N Af pp 144
(3.63) η0 = = 1− (1 − ηf ) .
Atotal Atotal
Eq. (3.63) can be used to calculate the total heat transfer for a fin array.
CHAPTER 4
Transient Conduction
4.1. Generalities
There are of course numerous engineering situations where unsteady ef-
fects are important. For example, the annealed steel quenching problem
that we alluded to in Chapter 1. Consider also solar systems whose input is
a function of the sun’s position, which changes as a function of time, or the
start–up and shut–down of any generic thermal system, for instance inter-
nal combustion engines. Mathematically, these represent fairly complicated
conduction problems, i.e. whereas steady 1–D conduction problems can be
reduced to an ordinary differential equation, transient problems (in general)
result in a partial differential equation, which is typically much more difficult
to solve. Specifically, if we simplify the general equation of heat conduction,
Eq. (2.18) on pp. 12, to model a 1–D transient problem, we get
∂T ∂2T
(4.1) = α ,
∂t ∂x2
where α is the thermal diffusivity. This equation is quite a bit more difficult
to solve than what we’ve seen up to now.
where t0 and θ 0 are variables of integration. Note that on the left hand side,
we are explicitly using limits that correspond to the initial condition (at
t = 0) and to the condition at the time t of interest, i.e.
θ|t0 =0 = θ0 and θ|t0 =t = θ(t) = θ .
Carrying out the integration, we find
θ hAsurf 0 t
ln θ 0 θ0 = − t 0 ,
ρV c
which evaluates to
hAsurf
(ln θ − ln θ0 ) = − t.
ρV c
Using the logarithm identity (ln a − ln b) = ln(a/b), we can write this ex-
pression more conveniently as
θ hAsurf
ln = − t.
θ0 ρV c
4.3. APPLICABILITY OF LUMPED CAPACITANCE 39
Finally, we can exponentiate both sides to obtain the more useful form
hAsurf
(4.3) θ = θ0 exp − t .
ρV c
Substituting back the actual temperatures, we can also write the alternative,
but equivalent form
T (t) − T∞ hAsurf
(4.4) = exp − t ,
T (0) − T∞ ρV c
where T (0) is the initial temperature of the conductor.
According to Eq. (4.4), temperature response is a simple exponential
when spatial gradients are small enough to be neglected. We can write this
equation as
T (t) − T∞
(4.5) = e−t/τ ,
T (0) − T∞
where
ρV c 1
(4.6) τ = = (ρV c)
hAsurf hAsurf
is a “time constant”. Notice that the first term is the standard form for
the convective resistance we examined back in Chapter 3. Increasing this
resistance, i.e. decreasing either h, A surf , or both, increases the response
time of the system. The second term, ρV c, is the lumped thermal capacitance
of the conductor. Likewise, any increase of this term increases response time.
In engineering systems, we would be interested not only in T (t), but also
in how much heat was transferred over a particular time period. Determining
the total energy transferred over some period of time is simply a matter of
integrating Newton’s Law of Cooling over a particular time range
Z t Z t
0
Q = q dt = hAsurf θ dt0 ,
0 0
which gives
(4.7) Q = ρV c θ0 1 − e−t/τ .
Ts2
Ts2 Too
x=0 x=L
2L
x R
H R
sphere: radius R
slab: width 2L
cylinder: radius R
x*
convection
4.4This is the unfortunate nomenclature used in most texts. The Fourier number is
not a dimensionless “number” in the same sense as, for example the Reynolds number in
fluid mechanics. It is a function of time, although it does not have any units. Therefore,
it would be more appropriately referred to as a dimensionless coordinate. However, as
a first approximation, it does also connote the relative effectiveness with which a body
conducts and stores energy (e.g. Incropera and Dewitt, 2002, pp. 256). In that sense, it
could also be loosely interpreted as a “dimensionless number”.
4.5. TRANSIENT ANALYSIS IN ONE DIMENSION 45
Bi ζ1 C1 Bi ζ1 C1 Bi ζ1 C1
0.01 0.0998 1.0017 0.2 0.4328 1.0311 2.0 1.0769 1.1795
0.02 0.1410 1.0033 0.25 0.4801 1.0382 3.0 1.1925 1.2102
0.03 0.1732 1.0049 0.3 0.5218 1.0450 4.0 1.2646 1.2287
0.04 0.1987 1.0066 0.4 0.5932 1.0580 5.0 1.3138 1.2402
0.05 0.2217 1.0082 0.5 0.6533 1.0701 6.0 1.3496 1.2479
0.06 0.2425 1.0098 0.6 0.7051 1.0814 7.0 1.3766 1.2532
0.07 0.2615 1.0114 0.7 0.7506 1.0919 8.0 1.3978 1.2570
0.08 0.2791 1.0130 0.8 0.7910 1.1016 9.0 1.4149 1.2598
0.09 0.2956 1.0145 0.9 0.8274 1.1107 10.0 1.4289 1.2620
0.10 0.3111 1.0160 1.0 0.8603 1.1191 20.0 1.4961 1.2699
∞ 1.5707 1.2733
where t∗ = αt/r02 and J0 and J1 are Bessel functions of the first kind. Here,
eigenvalues are the roots of
J1 (ζn )
ζn = Bi .
J0 (ζn )
So, we find that the problems are more difficult, but their solutions appear
to be more difficult to evaluate, as well. Luckily, we find that the one–term
approximation is again valid for t∗ = F o > 0.2, so that tabulated values
(e.g. Incropera and Dewitt, 2002, pp. 258) can be used directly.
4.5Recall from Eq. (4.8) the form of the Biot number is Bi = hL/k, so that an infinite
Biot number is equivalent to prescribing an infinite convection coefficient h.
4.6. THE SIMILARITY TECHNIQUE 48
by the same equation and initial conditions as the other problems, however,
boundary conditions are special. We still require 2 boundary conditions
because of the second–order spatial derivative, but there is clearly only 1
surface to specify a boundary condition at, i.e. x = 0. There is no boundary
at x → ∞. However, we can use the conceptual boundary condition of
T (x → ∞, t) = T0 ,
which says simply that for some value x that is large enough, the temper-
ature will always be the initial temperature. Physically, this makes sense
because if truly x → ∞, then we can always pick x large enough such that
it is so far away from the boundary at x = 0 that temperature variations
cannot be perceived. For the surface, we can specify any of the 3 standard
types of boundary conditions we’ve studied.
The mathematical significance of this problem is that there is a very
clever trick to simplify the partial differential equation, Eq. (4.1), to an
ordinary differential equation. This trick is known as a similarity transform
and effectively converts the problem of two independent variables, x and t,
into a problem of just one independent variable, the similarity transform
variable η. Such a transform is often seen in fluid mechanics and may be
possible for any problem that doesn’t have naturally–identifiable scales, for
example problems where there’s no clear length scale because of infinite
dimensions. Determining what the similarity transform variable is for a
4.6. THE SIMILARITY TECHNIQUE 49
the process is actually not trivial. If, for example, we prescribe boundary
conditions of the first–kind, e.g. T = T s at x = 0, we obtain the solution
T − Ts
(4.27) = erf η ,
T0 − T s
where erf is the Gaussian Error Function. This function is tabulated in
reference texts (e.g. Incropera and Dewitt, 2002, Appendix B.2, pp. 935).
The surface heat flux can be obtained by applying Fourier’s Law, from which
one obtains
Ts − T 0
(4.28) qs00 = k √ ·
παt
Although somewhat more involved, solutions for boundary conditions of the
second and third kind can be obtained as well (Incropera and Dewitt, 2002).
y=L2 y
x
x = − L1
x=L1
y = − L2
and
∂2θ ∂ 2 θ1 · θ 2
∂ ∂θ1 ∂θ2
2
= = θ2 + θ1 ,
∂x ∂x2 ∂x ∂x ∂x
so that
∂2θ ∂θ2 ∂θ1 ∂ 2 θ1 ∂θ1 ∂θ2 ∂ 2 θ2
= + θ 2 + + θ 1 .
∂x2 ∂x ∂x ∂x2 ∂x ∂x ∂x2
However, recalling that ∂θ2 /∂x = 0 since θ2 is not a function of the x
coordinate, this term simplifies to
∂2θ ∂ 2 θ1
= θ 2 .
∂x2 ∂x2
According to a similar procedure, we see that
∂2θ ∂ 2 θ2
= θ 1 .
∂y 2 ∂y 2
Substituting these terms into Eq. (4.29) and moving everything to the left
hand side, we get
∂ 2 θ1 ∂ 2 θ2
∂θ1 ∂θ2
(4.32) θ2 − + θ1 − = 0.
∂t ∂x2 ∂t ∂y 2
Eq. (4.32) is the original governing equation written in terms of our product
solution hypothesis. If it is valid, then either θ 1 and θ2 must both vanish, or
the terms within the brackets must both vanish. Now, θ 1 and θ2 do not gen-
erally vanish, otherwise the problem is trivial. Examining this equation more
closely we see that the terms in brackets are nothing more than Eqs. (4.30)
and (4.31). Thus, the terms in brackets do in fact vanish, and the original
hypothesis of the product solution is valid. As mentioned above, boundary
conditions can be treated similarly. Note how the presence of a heat gener-
ation term would prevent us from drawing our final conclusion since neither
of the equations for θ1 and θ2 would necessarily be zero. Of course, from
our one–dimensional work, we know various solutions corresponding to θ 1
and θ2 and the rules under which they (and their simplifications) can be
applied. So, once again, we’ve managed to find a significant simplification I&D Ex. 5.7
for a special case, i.e. two–dimensional problem with no heat generation. pp 277
Once it is determined that the product solution is valid, solving a prob-
lem proceeds by treating the one–dimensional problems individually. In
particular
• calculate individual Biot numbers for each coordinate direction us-
ing the appropriate length scale for each coordinate
• likewise, calculate individual Fourier numbers (dimensionless times)
for each coordinate direction
• compute individual one–dimensional solutions, e.g. θ 1 and θ2 , for
the required conditions of the two–dimensional problem
• take the product of the one–dimensional expressions as the final
solution of the problem
CHAPTER 5
Introduction to Convection
u oo Too
Ts x
dx
One of our basic assumptions will be that both T ∞ and Ts are constant, so
that their difference, Ts − T∞ , is also constant. The integration can then be
53
5.1. BOUNDARY LAYER INTRODUCTION 54
simplified as
Z
q = (Ts − T∞ ) h dAs .
As
Let us separately define an average heat transfer coefficient by again apply-
ing Newton’s Law of Cooling
q = h̄ (Ts − T∞ ) As .
From these two equations, it is clear that we can write the average convection
coefficient in terms of its local value as
1
Z
(5.1) h̄ = h dAs .
A s As
So we see that what we were really talking about all along in conduction
was an averaged value of the convection coefficient. Now the task becomes
actually figuring out what this is for specific problems.
For the special case of the flat plate in Fig. 5.1, the average convection
coefficient can be simplified. If we take the length of the plate as L and the
width as W , then As = W × L and dAs = W × dx. Integrating over the
length of the plate from the leading edge 0 → L, we find
Z L
1
(5.2) h̄ = h dx .
L 0
Determining convection coefficients is viewed as the “convection prob-
lem”. However, the problem is not trivial, since now we have the additional
complexity of fluid motion. So, not only is geometry, specific heat, conduc-
tivity, etc. important, we must now also worry about flow conditions, fluid
density, viscosity, etc.. Much of this depends upon the concept of boundary I&D Ex. 6.1
layer theory, which can be used to describe flow regions in the neighborhood pp 329
of a solid surface.
Too
Too
y
δt
x
Ts
FLOW
laminar
turbulent
The Prandtl number is clear in light of the following from the energy equa-
tion
u∞ L u∞ L ν u∞ L ν
= = = Re · P r .
α α ν ν α
Thus, the final non–dimensionalized form of the equations can be written
∂u∗ ∂v ∗
(5.16) + = 0,
∂x∗ ∂y ∗
∂u∗ ∗ ∂u
∗ dP ∗ 1 ∂ 2 u∗
(5.17) u∗ + v = − + ,
∂x∗ ∂y ∗ dx∗ Re ∂y ∗ 2
and
∂T ∗ ∂T ∗ 1 ∂2T ∗
(5.18) u∗ + v∗ = .
∂x ∗ ∂y ∗ Re · P r ∂y ∗ 2
The Reynolds number is familiar from fluid mechanics. It can be interpreted
as the ratio of inertial forces to viscous forces, so it serves as a direct metric
of how well frictional dissipation is able to dampen perturbations in a flow.
This of course bears directly on whether a particular flow is laminar or
turbulent.
Conversely, the Prandtl number is new — it appears only in the energy
equation. Taking a closer look, we see that P r is a direct property of the
fluid itself. It does not contain any non–fluid properties, as for example the
velocity and length scales in the Reynolds number. We can therefore speak
of “the Prandtl number of a fluid”. How do we interpret this parameter
physically? The Prandtl number is the ratio of kinematic viscosity ν to
thermal diffusivity α, the former of which indicates the rate of momentum
diffusion (from the viscous terms in the momentum equations) and the latter
of which is the rate of thermal diffusion
ν momentum diffusion
Pr = = ·
α thermal diffusion
We can say equivalently that P r represents the ratio of the effectiveness
of diffusional energy transport in the velocity boundary layer versus the
thermal boundary layer. In laminar flows, we can thus infer
δ
Pr ∼ ,
δt
that is, P r gives some approximate indication of the thickness of the velocity
boundary layer versus the thermal boundary layer 5.5 (Table 5.1).
Of course, we still have not lost anything in terms of non–linearity of
the system in Eqs. (5.16) through (5.18), however, just as previously, the
5.5This is not the case in turbulent flows where transport is a function of turbulent
mixing, as well as diffusion.
5.3. DIMENSIONLESS PARAMETERS 62
Fluid Pr relationship
gases ∼1 δ ≈ δt
liquid metals 1 δ δt
oils 1 δ δt
number of relevant variables to any problem has been reduced. For example,
the velocity distribution should be of the form
∗ ∗ ∗ ∗ dP ∗
(5.19) u = u x , y , Re, ,
dx∗
where dP ∗ /dx∗ is considered to be determined by the geometry of the prob-
lem and the free–stream conditions and is considered a known parameter as
mentioned previously.
What other dimensionless parameters arise for such configurations? Let
us consider frictional drag on the body surface that arises because of viscous
shear stress. Recall from fluid mechanics that the stress is defined as
∂u µ u∞ ∂u∗
τ = µ = .
∂y
y=0 L ∂y ∗ ∗y =0
This expressions allows us to derive the following form of the skin friction
factor , defined by Cf = τ /(ρu2∞ /2), as
2 ∂u∗
Cf = .
Re ∂y ∗ y∗ =0
Eq. (5.19) indicates that ∂u∗ /∂y ∗ evaluated at y ∗ = 0 is only a function of
(x∗ , Re, dP ∗ /dx∗ ). Therefore, we arrive at the significant conclusion that
the skin friction factor, an important engineering quantity, is only governed
by three “universal” parameters according to the functional form
∗ dP ∗
(5.20) Cf = Cf x , Re, .
dx∗
Likewise, we see that T ∗ is essentially a function of the form
∗ ∗ ∗ ∗ dP ∗
(5.21) T = T x , y , Re, P r, ,
dx∗
so that it has the same dependencies of u ∗ , except with the addition of the
Prandtl number as another parameter. Recall that we cast the convection
coefficient in Eq. (5.3) on pp. 56 as
kf ∂T
h = − ,
Ts − T∞ ∂y y=0
5.4. REYNOLDS–COLBURN ANALOGY 63
derive this solution in the following chapter. The only non–rigorous component of the
process is where we substitute a certain constant derived from the exact solution for its
counterpart from the approximate solution. This will be pointed out.
5.4. REYNOLDS–COLBURN ANALOGY 65
External Convection
N ux = C Rem n
L Pr ,
Ts + T ∞
(6.1) Tf = ,
2
66
6.1. LAMINAR FLOW OVER A FLAT PLATE 67
∂u ∂v
(6.3) + = 0,
∂x ∂y
and
∂T ∂T ∂2T
(6.4) u + v = α ·
∂x ∂y ∂y 2
In Chapter 5, we derived the Reynolds–Colburn analogy using analytical
solutions, both approximate and exact, to this problem. Where do these
solutions come from? The exact solution is not trivial, in fact it is beyond
the mathematical capabilities we have developed thus far. Again, this is a
result of the non–linearity of the system in Eqs. (5.16) through (5.18) on
pp. 61. We can cast the solution procedure, but will be unable to follow it
to completion.
The similarity technique introduced in §4.6 can be used to reduce the
partial differential system to an ordinary differential equation. As we men-
tioned previously, one should suspect that a similarity approach may work
when no inherent scales can be identified in the problem. In this case, the
transformation variables
r
u∞
η = y and
νx
r
ψ νx
f (η) =
u∞ u∞
are appropriate, where ψ = ψ(x, y) is the stream function, which is defined
as
∂ψ ∂ψ
(6.5) u = and v = − ·
∂y ∂x
Transforming the momentum equation appropriately, we find it becomes
2f 000 + f f 00 = 0 ,
6.1. LAMINAR FLOW OVER A FLAT PLATE 68
η f f 0 = u∗ f 00
0 0 0 0.332
0.4 0.027 0.133 0.331
0.8 0.106 0.265 0.327
1.2 0.238 0.394 0.317
1.6 0.420 0.517 0.297
2.0 0.650 0.630 0.267
2.4 0.922 0.729 0.228
2.8 1.231 0.812 0.184
3.2 1.569 0.876 0.139
3.6 1.930 0.923 0.098
4.0 2.306 0.956 0.064
4.4 2.692 0.976 0.039
4.8 3.085 0.988 0.022
5.2 3.482 0.994 0.011
5.6 3.880 0.997 0.005
6.0 4.280 0.999 0.002
6.4 4.679 1.000 0.001
6.8 5.079 1.000 0.000
The first term remains the same, however, to further develop Eq. (6.14), we
make the following observations. The second term can be evaluated using
Eq. (6.12) for v and observing that u = 0 at y = 0 and u = u 0 at y = δ.
In the third term, we utilize the original continuity equation (6.10) to swap
∂v/∂y = − ∂u/∂x. The fourth term is evaluated using the fact that ∂u/∂y
is zero at y = δ. These modifications give
Z δ Z δ Z δ
∂u ∂u ∂u ∂u
(6.15) u dy − u0 dy + u dy = −ν ,
0 ∂x 0 ∂x 0 ∂x ∂y y=0
which can be simplified to
Z δ Z δ
∂u ∂u ∂u
(6.16) u0 dy − 2u dy = ν .
0 ∂x 0 ∂x ∂y
y=0
The u0 can be taken directly under the integral and the u can be “integrated”
in6.3, which yields:
Z δ
∂ ∂u
(6.17) (u0 · u − u · u) dy = ν .
0 ∂x ∂y y=0
Now, the right hand side is basically a term that is evaluated, and the
differential on the left hand side can be taken outside the integral (since the
term is an integral w.r.t. y) and done later, therefore ∂/∂x → d/dx, which
yields:
Z δ
d ∂u
(6.18) (u0 − u) · u dy = ν .
dx 0 ∂y y=0
Eq. (6.18) is an integral form of the boundary layer equations. If we examine
this equation carefully, we see that the right hand side is basically the shear
stress at the surface of the plate. That is,
Z δ
d
(6.19) (u0 − u) · u dy = τw .
dx 0
This quantity is unknown because we do not know the velocity profile u.
Similarly, the left hand side also cannot be evaluated because it contains
u. We apparently have an implicit equation in u that cannot be solved in
closed form.
However, we can apply reasonable approximations for u. For example,
let us assume that u is described by a polynomial 6.4. We know that the flow
is developing, therefore it is a function of both x and y, i.e. u = u(x, y).
However, we would like to use a polynomial that is a function of just a single
6.3i.e. 2u ∂u = ∂u2
∂x ∂x
6.4Here is where the idea of various degrees of accuracy enters into the problem.
Specifically, lower–order polynomials will, in general, yield less accurate results as com-
pared to higher–order polynomials. The order is limited by how many boundary conditions
we can identify because these must be used to evaluate the polynomial coefficients.
6.2. KARMAN–POHLHAUSEN INTEGRAL SOLUTION 71
which simplifies to
39 2
u δ
280 0
• The right hand side of Eq. (6.18) is easily evaluated as
3 νu0
2 δ
• Eq. (6.18) can then be written as
d 39 2 3 νu0
u0 δ =
dx 280 2 δ
• We note that δ is the only term on the left hand side that is a func-
tion of x. Everything else is constant with respect to x, therefore
we can write the equation as
39 2 dδ 3 νu0
u =
280 0 dx 2 δ
which is a simple separable differential equation that can be written
as
140 ν
δ dδ = dx
13 u0
• This equation can be integrated along the flow direction, that is,
in x as
140 ν
Z Z
δ dδ = dx
13 u0
6.2. KARMAN–POHLHAUSEN INTEGRAL SOLUTION 73
somewhat similar, but more tedious and leads to Eq. (5.24). Development
of this result is shown in a number of reference texts (e.g. Burmeister, 1983,
§8.3).
6.3. EMPIRICAL CORRELATIONS 74
Using hlam based on Eq. (6.8) and hturb based on Eq. (6.31), we find I&D Ex. 7.1
4/5
pp 400
(6.33) N uL = 0.037 ReL − 871 P r 1/3 ,
where the limits are 0.6 ≤ P r ≤ 60 and the end of the plate at x = L is
assumed to lie in the turbulent zone, i.e. 5 × 10 5 ≤ Re ≤ 108 . (We assume
that transition occurs at Re = 5 × 105 .)
from fluid mechanics, flow separation occurs at some point where the bound-
ary layer detaches from the surface, forming a low–pressure wake. The flow
is characterized by the Reynolds number
u∞ D
ReD = ,
ν
where u∞ is the velocity upstream of the cylinder and D is the diameter,
and by the coefficient of drag
FD
CD = ,
Af ρ u2∞ /2
where FD is the measured drag force and Af is the projected frontal area
normal to the flow. The character of the boundary layer (laminar or turbu-
lent) strongly influences where separation occurs, and thus the magnitude of
CD (Fig. 6.2). There is a significant decrease as transition and turbulent flow
become important at Reynolds numbers above 200,000 to 500,000. These
6.3. EMPIRICAL CORRELATIONS 76
1000
100
drag coefficient
10
0.1
0.01
0.1 1 10 100 1000 10^4 10^5 10^6
Reynolds number
Internal Convection
boundary layer
inviscid core
boundary layer
x
developing flow (hydrodynamic entrance region) fully developed flow
7.2Again, except for the pressure gradient, which cannot vanish, otherwise there would
be no driving force to sustain the flow.
7.1. LAMINAR PIPE FLOW 79
7.3That is, ρu is the mass flux and c T is the energy per unit mass. Integrating
v
their product, the energy flux, over the cross–section gives the rate at which energy is
transported through the cross–section.
7.2. THE CASE OF CONSTANT HEAT FLUX 81
and by the same argument conclude that this result is not a function of x.
In other words, this expression remains fixed as the flow proceeds in the x
direction. If we substitute surface heat flux from Newton’s Law of Cooling
in Eq. (7.10), qs00 /h = (Tm − Ts ), and from Fourier’s Law of Conduction,
∂T /∂r = −qs00 /k at r = r0 then we find
" #
1 ∂T h qs00 h
= 00 − = ,
Tm (x) − Ts (x) ∂r r=r0
qs k r=r0 k
from which we must conclude that the local convection coefficient h is inde-
pendent of x for fully developed flow 7.4.
Furthermore, Eq. (7.10) indicates that for constant h and constant q s00 ,
the difference Ts − Tm must be a constant that does not vary with x. We
already know that Ts and Tm change with x, so that, if their difference is a
constant, the rates of change of these two temperatures must be the same
in the fully developed regime, i.e.
dTm dTs
(7.12) = ·
dx dx
Let us go one step further. We evaluate ∂θ/∂x using θ as defined in
Eq. (7.11). Our initial assumption of fully developed conditions requires
∂θ/∂x = 0. Evaluating, we find
∂θ ∂ T (r, x) − Ts (x)
=
∂x ∂x Tm (x) − Ts (x)
∂
= [T (r, x) − Ts (x)] [Tm (x) − Ts (x)]−1
∂x
1 ∂ [T (r, x) − Ts (x)]
=
Tm (x) − Ts (x) ∂x
T (r, x) − Ts (x) ∂ [Tm (x) − Ts (x)]
−
[Tm (x) − Ts (x)]2 ∂x
= 0,
The second term vanishes since ∂/∂x of [T m (x) − Ts (x)] was just shown to
be 0 in Eq. (7.12). This leaves
1 ∂ [T (r, x) − Ts (x)]
= 0,
Tm (x) − Ts (x) ∂x
where Tm (x) − Ts (x) is, again, a non–zero constant according to Eq. (7.10).
We can therefore write more concisely
∂ [T (r, x) − Ts (x)]
= 0,
∂x
which leads us to
∂T (r, x) dTs (x)
= ,
∂x dx
where we have changed ∂ to a plain derivative, since T s is only a function
of x. In light of Eq. (7.12), we now conclude
∂T (r, x) dTm
(7.13) = ·
∂x dx
In other words, the streamwise temperature gradient anywhere in the cross– I&D Ex. 8.1
section is equal to the streamwise gradient of the mean temperature. This pp 475
result will be used in solving for the temperature distribution T (r, x) in the
pipe and, ultimately, the Nusselt number for laminar fully developed flow.
We now address the actual problem of laminar fully developed pipe flow
having a constant applied heat flux q s00 at the pipe surface r = r0 . We assume
constant properties and no energy generation. Derivation of the convection
energy equation in cylindrical coordinates is beyond our present scope, but
it is shown in numerous reference texts (e.g. Bejan, 1984, Chapter 1). For
steady flow in (x, θ, r) coordinates (x is in the streamwise axial direction),
we have
2
1 ∂2T
∂T vθ ∂T ∂T ∂ T 1 ∂ ∂T
(7.14) u + + wr =α + 2 + r ,
∂x r ∂θ ∂r ∂x2 r ∂θ 2 r ∂r ∂r
where u is again the axial streamwise velocity component, and v θ and wr are
the azimuthal and radial velocity components. According to our hydrody-
namic solution, u = u (r) is given by Eq. (7.3), while v θ = wr = 0. Also, the
problem is symmetric, so that ∂/∂θ = 0. We will further assume 7.5 that the
mean temperature Tm is a linear function of x, so that ∂T /∂x is a constant
in Eq. (7.13) and therefore ∂ 2 T /∂x2 = 0. Eq. (7.14) then simplifies to
dTm 1 ∂ ∂T
(7.15) u(r) = α r .
dx r ∂r ∂r
If we examine this equation more closely, we notice that since dT m /dx is
constant, the equation can be written as an ordinary differential equation
dTm 1 d dT
(7.16) u(r) = α r .
dx r dr dr
Recasting this equation in dimensionless temperature θ, as defined by Eq. (7.11),
we write
dTm 1 d dθ
(7.17) u(r) = α (Tm − Ts ) r ,
dx r dr dr
which can be re–arranged as
1 d dθ u(r) dTm
(7.18) r = .
r dr dr α (Tm − Ts ) dx
7.5This can be shown according to an energy balance (e.g. Incropera and Dewitt,
2002, §8.3.1).
7.2. THE CASE OF CONSTANT HEAT FLUX 83
Using the solution for u(r) expressed in terms of the average velocity ū
in Eq. (7.6), we can write
" 2 #
d dθ r
(7.19) r = Ar 1− ,
dr dr r0
where A is the constant defined by
2 ū dTm
A = .
α (Tm − Ts ) dx
Eq. (7.19) can be integrated once to get
∂θ A r2 A r4
r = − + C1
∂r 2 4 r02
and then again to obtain
A r2 A r4
θ = − + C1 ln r + C2 ,
4 16 r02
where C1 and C2 are constants of integration.
The appropriate boundary conditions are similar to the hydrodynamic
problem. Temperature has a given value T s at the wall, r = r0 , so that θ = 0
at r = r0 . A boundary condition at r = 0 can be deduced from either of
• the temperature is well–behaved, so that the natural log term must
remain finite at r = 0
• because of symmetry of the overall problem, the temperature profile
must also be symmetric about r = 0
Both these observations yield the same result that C 1 = 0. Applying the
remaining condition, we find
3
C2 = − A r02 ,
16
which enables us to write θ as
A r2 A r4 3
θ = − 2 − A r02 .
4 16 r0 16
A equivalent, but more convenient form is
" 4 2 #
2 3 1 r 1 r
(7.20) θ = − A r0 + − .
16 16 r0 4 r0
The unknown value of A can be determined by employing the definition
of the bulk mean temperature in Eq. (7.9). Let us write this in dimensionless
form by forming θ on both sides (subtract T s , then divide by Tm −Ts ), which
yields
Z r0
Tm − T s 2 T − Ts
(7.21) = 2 u r dr .
Tm − T s ū r0 0 Tm − T s
7.2. THE CASE OF CONSTANT HEAT FLUX 84
Combining this with the expression for Newton’s Law of Cooling in Eq. (7.10),
we can express the convection coefficient as
kf ∂T
h = − .
Tm − Ts ∂r r=r0
Let us now express h in terms of the dimensionless temperature θ from
Eq. (7.11), where T (r, x) = θ [Tm (x) − Ts (x)] + Ts (x). We find
∂θ
(7.25) h = − kf .
∂r r=r0
y T1
u1
T0 x
connecting T0 and T1 (Fig. 7.3). For other cases, viscous dissipation distorts
the straight–line relationship as shown by plotting θ versus y ∗ in Fig. 7.3.
An interesting aspect of this problem is the heat flux at the top wall, defined
0.8
dimensionless "y"
0
=
0.6 Ec 2
* =
Pr Ec
*
0.4 Pr =4
c
*E
Pr
0.2 8
Ec =
Pr *
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
dimensionless temperature
by
dT
q = −k .
dy y=L
From Eq. (7.28), we evaluate the derivative to obtain
µ u21 1
dT T1 − T 0 2 y
= + − 2
dy y=L L 2k L L y=L
T1 − T 0 µ u21
= −
L 2kL
µ u21
T1 − T 0
= 1 −
L 2 k (T1 − T0 )
T1 − T 0 P r Ec
= 1− ,
L 2
which means that the heat flux is
T1 − T 0 P r Ec
(7.30) q = −k 1− .
L 2
There are clearly three cases of interest, depending upon the sign of the
(1 − P r Ec/2) term (again, assuming T 1 > T0 ):
• P r Ec < 2, which makes the right hand side negative and therefore
heat flows in the −y direction, or from the upper wall into the
liquid. This is what we would intuitively expect, given T 1 > T0 .
7.4. EMPIRICAL CORRELATIONS 88
• P r Ec > 2, which makes the right hand side positive and therefore
heat flows in the +y direction, or from the liquid into the wall, even
though the upper wall is at a higher temperature than the lower
wall. Not intuitive.
• P r Ec = 2, which makes the right hand side vanish and thus there
is no heat transfer at the upper wall. The upper wall behaves as
an insulated boundary (q 00 = 0), even though there is no actual
insulation. Also not intuitive.
The opposite case from the above is T 0 = T1 . Here, Eq. (7.28) reduces
to
µ u21 y y
T (y) − T0 = 1− .
2k L L
Non–trivial temperature distributions are due strictly to viscous dissipation
effects. This configuration is symmetric about y = L/2. In fact, the maxi-
mum temperature occurs at the middle of the channel and can be computed
by setting y = L/2 to get
µ u21
Tmax − T0 = .
8k
Though not identical, this case is qualitatively similar to heat generation for
conduction heat transfer in §3.4.
reasonably applied for both the constant temperature and constant heat flux
boundary conditions.
It is emphasized that both Eqs. (7.31) and (7.32) should only be used
for fully turbulent problems. The transition regime between laminar and
fully turbulent flow, 2000 < ReD < 104 , entails additional considerations.
Incropera and Dewitt (2002, §8.5) discuss additional results for these cases.
7.4.2. Non–Circular Tubes. Many engineering applications involve
tubes or ducts having non–circular cross sections. Even for laminar flow,
theoretical approaches are beyond the mathematics we have discussed thus
far. For example, for a rectangular cross section, the momentum equation
cannot be reduced to an ordinary differential equation like what we obtained
for circular pipes in Eq. (7.1). However, results we have obtained for circu-
lar tubes can be applied as an initial approximation to non–circular cross
sections using the hydraulic diameter as a length scale
4 Ac
(7.33) Dh = ,
P
where Ac is the cross–sectional area of the conduit and P is the wetted
perimeter7.6. The value of Dh can be computed directly from geometry, for
example, for a rectangular duct of height a and width b, we obtain
4ab 2ab
(7.34) Dh = = ,
2 (a + b) a+b
Thus, Dh should be used to calculate the non–circular analogs of the dimen-
sionless parameters, e.g. ReDh and N uDh .
For fully turbulent flows, the Dittus–Boelter and Sieder–Tate correla-
tions can be reasonably used in conjunction with the hydraulic diameter.
Results are less accurate for laminar flows, especially in geometries having
sharp corners. Here, Nusselt numbers based on exact solutions should be
used. For example, Table 7.1 shows selected results for channels of rectan-
gular cross section arranged according to the aspect ratio φ = b/a, where
b > a. Kays and Crawford (1980) consider such problems in greater detail.
7.6Eq. (7.33) contains a factor of 4 so that the hydraulic diameter of a circular pipe
is equal to its geometric diameter, i.e. Dh = 4 πD2 /4 / (πD) = D.
` ´
CHAPTER 8
Natural Convection
hot cold
g g
y y
cold hot
tendency for the system to move away from this state of equilibrium and
any heat transfer between the plates will be a result purely of conduction.
However, if the bottom wall is hotter (Fig. 8.1, right), the hot light fluid
will initially be below the cooler heavier fluid. This “top heavy” system will
tend to flow, with the heavy fluid falling by the action of gravity and the
light fluid correspondingly rising. A circulation pattern will arise, so that
heat transfer between the plates is effected primarily by natural convection.
90
8.1. WALL BOUNDED CONVECTION ON A VERTICAL FLAT SURFACE 91
g x
vertical plate
The first term on the right hand side is the buoyancy force and in this
case it is the driving force of the flow. If temperature is the only factor
affecting density8.1, this term can be described by a thermodynamic property
known as the volumetric thermal expansion coefficient
1 ∂ρ
β = − ,
ρ ∂T P
which describes the change of density as a function of temperature at con-
stant pressure. This can be approximated by a simple finite difference ex-
pression
1 ρ∞ − ρ
β ≈ − ,
ρ T∞ − T
which implies
ρ∞ − ρ ≈ ρ β (T − T∞ ) .
This idealization is the so–called Boussinesq approximation for density vari-
ation. Substituting, the streamwise momentum equation becomes
∂u ∂u ∂2u
(8.1) u + v = g β (T − T∞ ) + ν ,
∂x ∂y ∂y 2
which nicely shows how the driving buoyancy force is related to a tempera-
ture difference.
Buoyancy effects are strictly limited to the momentum equation, there-
fore the remaining equations are unaffected. Specifically, for conservation of
mass, we have the usual
∂u ∂v
(8.2) + = 0
∂x ∂y
and for conservation of energy, we likewise have
∂T ∂T ∂2T
(8.3) u + v = α .
∂x ∂y ∂y 2
Eqs. (8.1) through (8.3) are the boundary layer form of the natural convec-
tion equations. Viscous dissipation has been neglected here, primarily due
to the small velocities associated with natural convection.
Examining these equations, we discover a major mathematical compli-
cation as compared to forced convection: T now appears in the momentum
equation in addition to u appearing in the energy equation. Therefore, this
system is fully coupled, as opposed to forced convection, which was uncou-
pled since T only appeared in the energy equation. We can no longer treat
the hydrodynamic problem independently. The equations must be solved
simultaneously.
8.1Density changes can arise via other means, as in for example supersonic flow. We
do not consider such cases here.
8.2. DIMENSIONLESS FORMULATION 93
The Grashof number plays the same role in free convection that the Reynolds
number plays in forced convection, i.e. Re is the ratio of inertial forces to
viscous forces, whereas and Gr is the ratio of buoyancy forces to viscous
forces. We can then recast momentum conservation in Eq. (8.5) as
∂u∗ ∗ ∂u
∗ 1 ∂ 2 u∗
(8.8) u∗ + v = Gr T ∗
+ ,
∂x∗ ∂y ∗ Re ∂y ∗ 2
The three conservation laws, Eqs. (8.4), (8.6), and (8.8), are functions of
the dimensionless parameters Re, P r, and Gr. We would therefore expect
solutions of the natural convection problem, given by the Nusselt number,
to be of the form
N u = N u (Re, P r, Gr) .
There is a subtle aspect to this interpretation. In formulating the Grashof
number, we were able to eliminate the velocity scale u 0 based on the obser-
vation that we could not necessarily characterize it appropriately. Yet, it
could be that an external forcing results in an explicit u 0 , so that a problem
would be a mixture of forced convection and natural convection. The above
presumption about the functional dependence of N u is valid when forced
and natural convection aspects are comparable, i.e.
Gr
∼ 1.
Re2
Looking back at the form of Gr in Eq. (8.5), this condition is equivalent to
g β (Ts − T∞ ) L ∼ u20 .
Consequently, if forced convection dominates a problem,
Gr
1 and N u = N u (Re, P r) .
Re2
That is, free convection effects are so small as to be neglected. Conversely,
if natural convection is dominant, we have
Gr
1 and N u = N u (P r, Gr) .
Re2
Here, the small forced convection effect would be neglected, although strictly
speaking, natural convection occurs exclusively only for Gr/Re 2 → ∞, for
which u0 = 0.
F (P r) k g β (Ts − T∞ ) 1/4 L 1
Z
= dx
L 4 ν2 0 x
1/4
wide range of Ra
2
1/6
0.387 RaD
(8.15) N uL = 0.6 + 8/27 , RaD < 1012 .
9/16
1 + [0.559/P r]
where the Rayleigh number is based on the pipe diameter D.
8.4.3. Rectangular Enclosures. Configurations we have discussed
thus far are primarily of the external flow type, but internal configurations
also arise in applications, for example in the form of various enclosures
(Fig. 8.4). Here, the bottom and top surfaces are adiabatic, while the ver-
cold
hot
H
tical surfaces have specific temperatures, T 1 on the hot side and T2 on the
cold side, and are separated by a distance L. If the Rayleigh number is suffi-
ciently high, a recirculating flow pattern will form, similar to that shown in
the figure. For small RaL , the flow is weak and heat transfer occurs mainly
by conduction, i.e. N uL ∼ 1. For larger RaL , a number of correlations have
been devised, e.g.
P r RaL 0.28 H −0.25
(8.16) N uL = 0.22 , 103 < RaL < 1010 ,
0.2 + P r L
which is valid for aspect ratios 2 < H/L < 10 and P r < 10 5 . Other
correlations are available for larger aspect ratios ( Özişik, 1985; Incropera
and Dewitt, 2002).
CHAPTER 9
Heat Exchangers
Heat exchangers are devices that facilitate heat transfer between two
fluids at different temperatures without allowing them to mix. Usually, the
fluids are separated by a solid boundary, for example as in car radiators,
air conditioners, distillers, etc. These are called indirect contact exchangers.
However, if the fluids do not tend to mix naturally, a direct contact heat
exchanger may be used, e.g. a water chiller.
Heat exchangers are typically classified by flow arrangement and type
of construction. In parallel flow units, fluids enter at the same end and flow
in the same direction, while in counter–flow, the flow of the two fluids is in
opposite directions. Fig. 9.1 illustrates these arrangements for a simple heat
exchanger consisting of concentric tubes.
parallel flow
counter flow
on this basic architecture are shown in reference texts (e.g. Özişik, 1985;
Mills, 1999; Incropera and Dewitt, 2002).
It is also common to classify heat exchangers according to the area avail-
able for heat transfer per unit volume of the unit, i.e.
total area available for heat transfer
ϕ = ·
total volume of heat exchanger
Units having ϕ > 700 m−1 are arbitrarily referred to as being compact heat
exchangers, regardless of their design architecture. Table 9.1 summarizes ϕ
for several classes of heat exchangers.
parallel flow
temperature
H
∆T
C
x=0 x=L
• The outlet temperature of the cool side cannot exceed that of the
hot side. At most, they can be equal, otherwise an un–physical
temperature gradient exists. Therefore, the effectiveness of parallel
flow is limited and not typically used for heat recovery.
• The temperature gradient between the streams, ∆T , is not con-
stant. Therefore the local heat flux q 00 varies with x.
• If we assume that the wall temperature of the material separat-
ing the streams is the average between the two streams, then wall
temperature is roughly constant. This might be an advantageous
design aspect if the wall material is sensitive to stresses induced by
temperature differences along x.
Fig. 9.4 shows thumbnail diagrams for three other cases: counter–flow,
and condensing and boiling heat transfer. For the counter–flow design, the
H
H
temperature
temperature
temperature
C
H
C
C
0 L 0 L 0 L
counter−flow condensing boiling
exit temperature of the cold fluid can in fact be higher than that of the
hot fluid. This makes the counter–flow design more attractive than parallel
flow units if the design requirements permit such a choice. Clearly, the wall
temperature in these units is not constant, which might cause large thermal
stresses that may eliminate certain materials from design consideration. Al-
though not necessarily the case, ∆T might also be constant along x for the
counter–flow arrangement.
Although we do not discuss the physics of condensation and boiling heat
transfer per se, we can examine their application in heat exchangers utilizing
the simple fact that temperature is constant during a phase change. It should
be clear that the flow direction of the phase change stream is irrelevant
because of its constant temperature. For the condenser, heat energy is
liberated from the condensate and absorbed by the cold side, causing its
temperature to rise9.1. For a boiling heat exchanger, the reverse is true: the
hot side transfers energy to boil the fluid on the cold side and cools in the
process.
Thumbnail diagrams can be drawn for more complicated configurations
as well. For example, Fig. 9.5 illustrates a shell–and–tube exchanger having
one shell pass and two tube passes, along with its corresponding thumbnail
diagram.
H
temperature
Figure 9.5. Thumbnail diagram for one shell pass, two tube
pass heat exchanger. We have assumed the tube side contains the
cold fluid.
conductive resistance
9.2By convention, U is the symbol for overall convection coefficient for heat exchang-
ers, rather than some form of h which is what we might have expected. This notation is
the standard (e.g. Özişik, 1985; Mills, 1999; Incropera and Dewitt, 2002).
9.3. LMTD ANALYSIS 104
If the wall thickness is minimal9.4 and the thermal conductivity is high, tube
wall resistance can be neglected to obtain
1
(9.5) U = 1 1 ·
hi + ho
9.3Calculating the overall heat transfer will involve multiplying by A using a form of
t
Newton’s Law of Cooling (see Eq. (9.7)), so that, as long as consistency is preserved with
respect to which area is used, the effect will cancel itself out.
9.4If wall thickness is small, then r ≈ r , which means that ln (r /r ) = ln (1) = 0,
i o o i
resulting in a vanishing wall resistance term.
9.3. LMTD ANALYSIS 105
T Th + d T h hot side
h
Tc dq Tc + d T c cold side
dA
Th,i dq
Th hot side T
h,o
dTh
∆ T0 ∆T ∆T
L
dTc
Tc,i
x=0 x=L
Notice that C0 comes outside the integral since it is a constant. The right
hand side can be multiplied by At /At , i.e.
Z ∆Tx=L Z At
d(∆T ) 1
(9.16) = − C 0 At × U dA .
∆Tx=0 ∆T At 0
Let us now define the average overall heat transfer coefficient for the whole
heat exchanger in the usual fashion as
Z At
1
(9.17) Um = U dA ,
At 0
so that by integrating the left hand side of Eq. (9.16), we get
∆Tx=L ∆Tx=L
(9.18) ln ∆T = ln = − U m At C0 ,
∆Tx=0 ∆Tx=0
or rewriting to remove the minus sign, it becomes
∆Tx=0
(9.19) ln = U m At C0 .
∆Tx=L
Eq. (9.19) is an interesting intermediate result, but is not yet in the form
of Newton’s Law of Cooling in Eq. (9.7). We can go back to Eq. (9.12) and
integrate this over the whole heat exchanger as
Z ∆Tx=L Z q
(9.20) d(∆T ) = − C0 dq 0 ,
∆Tx=0 0
from which we find that the total heat transfer q can be expressed as
∆Tx=L − ∆Tx=0 = − C0 q .
Solving for C0 , we obtain
∆Tx=0 − ∆Tx=L
C0 = ,
q
which can be substituted into Eq. (9.19)
∆Tx=0 ∆Tx=0 − ∆Tx=L
(9.21) ln = U m At ·
∆Tx=L q
9.3. LMTD ANALYSIS 108
so that
(Th,i − Tc,i ) − (Th,o − Tc,o )
∆Tm, parallel =
ln ([Th,i − Tc,i ] / [Th,o − Tc,o ])
(Th,i − Th,o ) − (Tc,i − Tc,o )
=
ln ([Th,i − Tc,i ] / [Th,o − Tc,o ])
∆Thot side + ∆Tcold side
=
ln (large number/small number)
Here, subscript “i” refers to an inlet temperature, while subscript “o” refers
to an outlet temperature. Although the numerator is large, the denominator
is not terribly small, i.e. it is the log of a rather large number.
Were we to duplicate the LMTD analysis for the counter–flow arrange-
ment, we would once again realize Eq. (9.23), except that the bookkeeping
aspect dictates
∆Tx=0 = Th,i − Tc,o and ∆Tx=L = Th,o − Tc,i ,
which can also be deduced from Fig. 9.4. We can then estimate ∆T m as I&D Ex. 11.1
(Th,i − Tc,o ) − (Th,o − Tc,i ) pp 655
∆Tm, counter−f low =
ln ([Th,i − Tc,o ] / [Th,o − Tc,i ])
(Th,i − Th,o ) − (Tc,o − Tc,i )
=
ln ([Th,i − Tc,o ] / [Th,o − Tc,i ])
∆Thot side − ∆Tcold side
=
ln (a number which may be close to unity)
Here, the numerator can be fairly small, but the denominator is the log of
an argument that can be close to unity. Therefore, the denominator is very
small, giving a larger LMTD temperature difference for the same input and
output values (Table 9.3).
Figure 9.8. Cross flow heat exchanger designs for both fluids
un–mixed (left) and one fluid mixed, one fluid un–mixed (right).
While such configurations have been analyzed, the results are usually
considered too complicated for practical use. Instead, these results have
been cast as correction factors to calculate appropriate LMTD temperatures.
The form of the corrected LMTD temperature is
∆Tm, corrected = ∆Tm F ,
where F is the appropriate correction factor and ∆T m is computed on the
basis of counterflow conditions! References (e.g. Özişik, 1985; Incropera
and Dewitt, 2002) typically present F graphically. Parameters in these I&D Ex. 11.2
graphical results depend by convention on the following nomenclature: t is pp 657
the temperature on the tube side (hot or cold), while T is the temperature of
the opposite (shell) stream. If the temperature change of one of the streams
is negligible, as in a phase change, then F = 1. Generally, an exchanger
with both sides being unmixed will give somewhat better performance than
if one side is mixed because the streamwise temperature gradient is more
conserved in the former.
T h, i
T c, o = T h, i
temperature
temperature
T h, o
T c, o
T c, i
T c, i = T h, o
so that the hot side outlet temperature would be cooled to the value of the
cold side inlet temperature, Th, o = Tc, i . Now, maximum energy transfer
would be written
q = Ch (Th, i − Th, o ) = Ch (Th, i − Tc, i ) where Ch = Cmin .
This proves the above concept.
Based on this exercise, it should be clear that that q max is not based on
the maximum ṁcp . If it were, then according to conservation of energy, the
fluid with the minimum ṁcp would undergo a temperature change greater
than the maximum ∆T , which is not possible. According to Eq. (9.24),
computing the heat transfer q depends on knowing . How do we determine
this?
The parameter depends upon the heat exchanger design and flow ar-
rangement. Here, we illustrate the general procedure for deriving assuming
a parallel flow single pass arrangement 9.5. This derivation refers, once again
to Fig. 9.7. From the definition in Eq. (9.24), we get
q
= .
Cmin (Th, i − Tc, i )
Also, based upon the total heat transferred, we have
q = Ch (Th, i − Th, o ) = Cc (Tc, o − Tc, i ) .
By definition, we then have
Ch (Th, i − Th, o ) Cc (Tc, o − Tc, i )
(9.28) = = .
Cmin (Th, i − Tc, i ) Cmin (Th, i − Tc, i )
Now, recall Eq. (9.19) that was used in the derivation of the LMTD analysis
method, again for the parallel flow arrangement
∆Tx=0 Th, i − Tc, i
ln = ln = U m At C0 ,
∆Tx=L Th, o − Tc, o
where we recall that
1 1 1 1
C0 = + = + ·
ṁh ch ṁc cc Ch Cc
Taking the exponential, we can write
Th, i − Tc, i
= e Um A t C0 ,
Th, o − Tc, o
which inverts as
Th, o − Tc, o
= e − Um At C0 .
Th, i − Tc, i
From Eq. (9.28), we can solve for Th, o to get
Cmin (Th, i − Tc, i )
Th, o = Th, i −
Ch
Cc (Tc, o − Tc, i ) Cmin (Th, i − Tc, i )
= Th, i −
Cmin (Th, i − Tc, i ) Ch
Cc (Tc, o − Tc, i )
= Th, i − ,
Ch
which can be used to eliminate Th, o from the exponential equation. This
yields
Cc (Tc, o − Tc, i )
Th, i − Ch − Tc, o
= e − Um At C0 .
Th, i − Tc, i
Add and subtract Tc, i to the numerator of the left hand side and we can
show that the expression becomes
C (T −T )
Th, i − Tc, i Tc, i − Tc, o − c c,Coh c, i
+ = e − Um At C0
Th, i − Tc, i Th, i − Tc, i
Cc Cc
Tc, i 1 + C h
− T c, o 1 + Ch
1 + =
Th, i − Tc, i
Cc Tc, i − Tc, o
1 + 1 + = e − Um At C0
Ch Th, i − Tc, i
Finally, this can be written as
Tc, i − Tc, o e − Um At C0 − 1
=
Th, i − Tc, i 1 + C Cc
h
9.5. –NTU ANALYSIS 114
or
Tc, o − Tc, i 1 − e − Um At C0
=
Th, i − Tc, i 1 + CCc
h
Referring back to Eq. (9.28), we see that the above expression can be directly
substituted to obtain , i.e.
C c 1 − e − Um At C0
=
Cmin 1 + CCc
h
1 − e − Um At C0
= Cmin Cmin
·
Cc + Ch
h i
1 − exp − Um At C1h + C1c
= ·
Cmin C1h + C1c
h i
1 − exp − UCmmin At
Cmin C1h + C1c
= ·
Cmin C1h + C1c
Such expressions are normally written in a short–hand form, where we define
the number of transfer units (NTU) as
Um At
NT U =
Cmin
and the ratio of the minimum to maximum heat capacity rates as
Cmin
C1 = ·
Cmax
It should be clear that, regardless of which side is the minimum side and
which is the maximum,
1 1
Cmin + = 1 + C1 ,
Ch Cc
which allows us to finally define for the parallel flow arrangement as
1 − e−N T U (1 + C1 )
(9.29) = ·
1 + C1
This equation is quite significant in heat exchanger analysis, because one
need not know the outlet temperatures. It suffices only to know the heat
capacity rates9.6, the area, and the overall convection coefficient to calculate
. Moreover, because of the way is defined, the total heat transfer can be
computed from a knowledge of only the inlet temperatures of both streams.
9.6There is, however, one subtle issue: the specific heats may depend upon tempera-
ture. If we were to use the average temperature of a stream to approximate specific heat,
we would need some a priori estimate of the exit temperature.
9.5. –NTU ANALYSIS 115
The procedure for analyzing other heat exchanger designs, e.g. counter–
flow, multiple passes, etc., is similar but more complicated as compared to
what we have shown here. If we simply reverse the design in our simple
tube–type heat exchanger for counter–flow conditions, we find 9.7
1 − e−N T U (1 − C1 ) / 1 − C1 e−N T U (1 − C1 )
: C1 < 1
(9.30) =
N T U/ (N T U + 1) : C1 = 1
Another notable configuration arises when there is a phase change (boil-
ing or condensation), in which case ∆T of the phase change side essentially
vanishes (Fig. 9.4). According to the fact that there is a finite thermal
energy transfer, conservation of energy implies that ṁc p can be taken as in-
finite. Therefore, the phase change side is always the “maximum” side and
the heat capacity ratio can be taken as zero, i.e. C 1 = 0. Moreover, heat
exchanger performance is independent of architecture and flow arrangement!
The analysis gives
(9.31) = 1 − e−N T U .
Reference texts (e.g. Incropera and Dewitt, 2002, Table 11.3, pp. 662–663)
furnish –NTU relationships for a number of configurations and the inverted
relations as well9.8. I&D Ex. 11.3
pp 665
9.7The formula for this case in Incropera and Dewitt (2002, Table 11.3, pp. 662) is
incorrect.
9.8That is, NTU as a function of is provided for the cases that are difficult to invert.
CHAPTER 10
Introduction to Radiation
visible
ultraviolet microwave
X − rays infrared
γ −rays thermal radiation
angular dependence
radiation intensity
wavelength
r y
dAn
d An
d ω = d An / r 2
θ r
(r, θ , φ )
r
dθ = d s / r ds
dω
dA
x
definition of "solid angle" and
φ
comparison to conventional
planar angle
z spherical coordinate system
differential expression10.3
ds
dθ =
·
r
Solid angle is instead defined in the context of a sphere as the ratio of a
subtended area and the square of the radius 10.4. We express this as the
differential
dAn
(10.1) dω = ·
r2
We consider solid angle in the context of an observer anchored at the origin
of a spherical coordinate system whose field of view encompasses a normally–
oriented area dAn a distance r away.
10.2Although we attach a “unit” of radians to the planar angle, the quantity is indeed
dimensionless. Angles expressed in the artificial system of degrees are not dimensionless.
10.3This clear from the observation that we can integrate this equation to obtain the
circumference as
Z θ0 =2π Z s0 =s
˛θ0 =2π
r dθ0 = ds , which gives s = r θ 0 ˛θ0 =0 = 2πr .
θ 0 =0 s0 =0
10.4The “unit” of the solid angle is the steradian, (abbreviated as sr ) which, like its
planar counterpart, is actually dimensionless.
10.1. SOLID ANGLE AND RADIATIVE QUANTITIES 119
Actually, we must take this a step further and define solid angle in terms
of differentials associated with the spherical coordinate system (Fig. 10.4).
Consider emission from a differential element dA along radius r in a specific
r sin θ r dθ
θ r sin θ d φ
r
dθ
dA
dφ
direction defined by angles θ and φ. When viewed from dA, the emitted
radiation passes through the differential element dA n on the boundary of
the hemisphere specified by radius r. By geometric arguments, we see
dAn = (r dφ sin θ) × (r dθ) = r 2 sin θ dφ dθ ,
so that the solid angle relationship is 10.5
dAn r 2 sin θ dφ dθ
(10.2) dω = = = sin θ dφ dθ .
r2 r2
Another geometric consideration in analyzing radiation emitted from dA
is its projected area on dAn . In other words, the question is how much of
dA can an observer anchored on dAn see? We should notice from Fig. 10.4
that when these areas are normal, i.e. θ = 0, our observer will see all of dA,
while at θ = π/2, surface dA cannot be seen at all. A little trigonometry
(Fig. 10.5) should convince us that
(10.3) dAn = dA cos θ .
d An
d An
θ
dA dA θ
power10.6
Z 2π Z π/2
(10.7) Eλ (λ) = Iλ, e (λ, θ, φ) cos θ sin θ dθ dφ ,
0 0
which has units of W/ m2 µm . The total hemispherical emissive power,
i.e. the rate at which radiation is emitted in all directions at all wavelengths
per unit area of the emitter, is obtained by integrating over λ
Z ∞ Z 2π Z π/2
(10.8) E = Iλ, e (λ, θ, φ) cos θ sin θ dθ dφ dλ ,
0 0 0
which has units of 2
W/m .
Although the directional distribution of surface radiation emission varies
according to the properties of specific surfaces, there is an idealization which
closely approximates conditions for many surfaces. This idealized form of
radiation is called diffuse radiation, which simply indicates that there is no
directional dependence. For example, for a diffuse emitter, the radiation
intensity is independent of direction, so that I λ,e (λ, θ, φ) → Iλ,e (λ). Since
Iλ,e (λ) is now constant with respect to directions, we see the emissive power
in Eq. (10.7) can be cast as
Z 2π Z π/2
Eλ (λ) = Iλ, e (λ) cos θ sin θ dθ dφ
0 0
Z 2π Z π/2
= Iλ, e (λ) cos θ sin θ dθ dφ
0 0
Z 2π Z π/2
= Iλ, e (λ) dφ cos θ sin θ dθ
0 0
2π 1 π/2
2
= Iλ, e (λ) φ sin θ
0 2 0
1
= Iλ, e (λ) 2π ×
2
(10.9) Eλ (λ) = π Iλ, e (λ) (diffuse emitter)
It can be similarly shown that E = πIe , where again E is the total emissive
power and Ie is the total radiation intensity. I&D Ex. 12.1
Now let’s look at some additional phenomena of interest in radiation pp 707
heat transfer. First is the concept of irradiation, i.e. incoming or incident
radiation. This is described by its own corresponding intensity, called the
incident spectral intensity, Iλ,i (λ, θ, φ), which, like the corresponding quan-
tity for emitters, is a function of wavelength and direction and is given per
unit area of the intercepting surface normal to this direction per unit solid
angle, per unit dλ. The corresponding flux is the irradiation. Specifically,
there is the spectral irradiation, G λ (λ), which is the rate at which radiation
10.6We include the subscript λ in the notation for E (λ) to indicate that this is a per
λ
unit wavelength dλ about λ quantity, synonymous with the way we have cast Eq. (10.5).
10.2. BLACKBODY RADIATION 122
1e+8
spectral emissive power (W/m**2 micrometer)
T = 50 K
T = 100 K
T = 500 K
1e+6 T = 1000 K
T = 2000 K
T = 5800 K
1e+4
1e+2
1e+0
1e-2
0.1 1 10 100
wavelength (micrometers)
10.7These products yield C = 3.743 × 108 W µm4 /m2 and C = 14, 387 µm K.
1 2
10.2. BLACKBODY RADIATION 124
The third point above has an interesting proof associated with it known as
Wien’s Displacement Law . If the Planck Distribution is differentiated with
respect to λ and set equal to zero then solved, the result is
(10.12) λ T ≈ 2897.8 µm · K .
This is shown in detail in Appendix C. The equation specifies the locus
of peak points for Eλ,b . Because the product λT is constant, the maxi-
mum spectral emissive power drifts toward shorter wavelengths for higher
temperatures.
What we are really interested in, i.e. what will enable us to solve prob-
lems, is the total blackbody emissive power. In principle, this can be derived
by integrating the spectral blackbody emissive power given by Eq. (10.11)
over all wavelengths
Z ∞
C1
(10.13) Eb (T ) = 5 eC2 /(λT ) − 1
dλ .
0 λ
Although this operation is not trivial 10.8 it can be shown that
(10.14) Eb (T ) = σ T 4 ,
where σ = 5.670 × 10−8 W/(m2 · K 4 ) is the Stefan–Boltzmann constant,
which depends upon the radiation constants C 1 and C2 . This equation is
known as the Stefan–Boltzmann Law and is important because it enables
calculation of the amount of radiation emitted in all directions and over all
wavelengths simply by knowing the blackbody temperature. Because the
blackbody is diffuse, we can back out the total blackbody intensity as
Eb σ T4
Ib = = .
π π
Often we need to know the fraction of radiative emission from a black-
body occurring over some finite wavelength range, e.g. λ = 0 → λ 1 . Rigor-
ously, this can be computed by simply integrating Eq. (10.11) appropriately,
e.g.
R λ1 Z λ1
0 Eλ,b dλ 1
f0→λ1 = R ∞ = Eλ,b dλ
0 Eλ,b dλ
σ T4 0
However, like the Stefan–Boltzmann Law, computation of this integral is
not trivial. Introducing the Planck Distribution from above, we see that
this gives
Z λ1
1 C1
f0→λ1 = 4 5 C /(λ T) − 1
dλ .
σT 0 λ e 2
b g(b)
du
Z Z
(10.15) f (u) dx = f (u) du ,
a dx g(a)
d (λ T )
du = d (λ T ) = dλ = T dλ ,
dλ
which leads to the substitution
d (λ T )
dλ → .
T
The upper limit is now λ1 → λ1 T . Therefore, the equation can be rewritten
as
Z λ1 T
1 C1 d (λT )
f0→λ1 =
σ T4 0
5
λ e C 2 /(λ T ) − 1 T
Z λ1 T
C1 1
= 5 d (λ T )
σ 0 (T λ) e 2 /(λ T ) − 1
C
The form and difficulty are similar to Eq. (10.13). Reference texts have tab-
ulated f as a function of the combined variable λT . For example, Incropera
and Dewitt (2002, Table 12.1, pp. 716) and Appendix C (Table C.1)give
abbreviated summaries, while Dunkle (1954) and Siegel and Howell (2001)
contain more comprehensive data.
` ´ `` ´˛˛g(b)
´
Z g(b)
= F g(b) − F g(a) = F u ˛ = f (u) du ,
g(a) g(a)
Eq. (10.16) can be applied more generally for arbitrary ranges of radia-
tion λ = λ1 → λ2 according to
R λ2
λ1 Eλ,b dλ
fλ1 →λ2 =
σ T4
R λ2 R λ1
0 Eλ,b dλ − 0 Eλ,b dλ
=
σ T4
R λ2 R λ1
0 Eλ,b dλ Eλ,b dλ
= 4
− 0
σT σ T4
(10.17) fλ1 →λ2 = f0→λ2 − f0→λ1
Note the special case for a range with a lower bound λ = λ 2 → ∞ corre-
sponds to I&D Ex. 12.3
R∞ pp 717
λ2 Eλ,b dλ
fλ2 →∞ = I&D Ex. 12.4
σ T4
R∞ R λ2 pp 719
0 Eλ,b dλ − 0 Eλ,b dλ
=
σ T4
R∞ R λ2
0 Eλ,b dλ Eλ,b dλ
= 4
− 0
σT σ T4
R λ2
σ T4 Eλ,b dλ
(10.18) fλ2 →∞ = 4
− 0 = 1 − f0→λ2
σT σ T4
10.3. Radiation Characteristics of Real Surfaces
10.3.1. Emissivity. Now that we have developed the concept of an
ideal surface (such that we have a basis of comparison), let us now look at
characteristics of real surfaces. In fact, we will specify properties of real sur-
faces via direct comparison to those of the blackbody. For example, recall
that the blackbody, being the ideal emitter, emits the maximum radiation
possible. Therefore, we can define the emissivity of a body as the ratio of
radiation emitted by its surface to the radiation emitted by a blackbody
at the same temperature. The spectrum of radiation emitted by a real
surface is not necessarily the same as that specified by the Planck Distribu-
tion (Fig. 10.7) Also, the directional distribution may not qualify as diffuse.
Therefore, emissivity can assume various forms according to whether one
talks of emission at a given wavelength or emission in a given direction, or
some overall integrated value.
First, we define a quantity called the spectral directional emissivity,
ελ,θ (λ, θ, φ, T ) of a surface at temperature T as the ratio of the intensity
of the radiation emitted at wavelength λ in the (θ, φ) direction to the inten-
sity of the radiation emitted by a blackbody for the same values of λ and
T . This has the form
Iλ,e (λ, θ, φ, T )
(10.19) ελ,θ (λ, θ, φ, T ) = .
Iλ,b (λ, T )
10.3. RADIATION CHARACTERISTICS OF REAL SURFACES 127
Planck (blackbody)
real surface
Pl
an
ce
ck
rfa
su
(b
la
al
ck
re
b
od
y)
wavelength
Note that the blackbody term in the denominator appears correctly as hav-
ing no dependence upon direction since it is diffuse emitter. Likewise, the
total directional emissivity, εθ , is an average over the wavelength spectrum
and is defined as
Ie (θ, φ, T )
(10.20) εθ (θ, φ, T ) = .
Ib (T )
Note that these two quantities are actually defined in terms of radiation
intensity because of directional dependence.
Conversely, we can also defined an emissivity that has been averaged
over all directions, called the spectral hemispherical emissivity
Eλ (λ, T )
(10.21) ελ (λ, T ) = ,
Eλ,b (λ, T )
which is given in terms of emissive power. We recall from previous definitions
that these quantities represent emissive powers that have been integrated
over the hemispherical area of emission. How does this relate to the direc-
tional emissivity ελ,θ ? The spectral hemispherical emissivity can be written
according to the definitions
R 2π R π/2
Iλ,e (λ, θ, φ, T ) cos θ sin θ dθ dφ
(10.22) ελ (λ, T ) = 0 R 2π0 R π/2 .
0 0 I λ,b (λ, T ) cos θ sin θ dθ dφ
Using Eq. (10.19), substitute Iλ,e = ελ,θ Iλ,b into the top term, take note that
Iλ,b in both numerator and denominator do not depend upon direction (so
they can be taken outside their integral signs), and the relation simplifies to
R 2π R π/2
ελ,θ (λ, θ, φ, T ) cos θ sin θ dθ dφ
(10.23) ελ (λ, T ) = 0 0 R 2π R π/2 .
0 0 cos θ sin θ dθ dφ
If we make the further assumption that ε λ,θ is independent of φ (a rea-
sonable assumption for many surfaces), we can obtain (after evaluating the
10.3. RADIATION CHARACTERISTICS OF REAL SURFACES 128
denominator)
Z π/2
(10.24) ελ (λ, T ) = 2 ελ,θ (λ, θ, T ) cos θ sin θ dθ .
0
The total hemispherical emissivity, ε, accounts for all directions and all
wavelengths and is defined as
E(T ) E(T )
(10.25) ε(T ) = = .
Eb (T ) σ T4
R∞
Now, using our previously established definitions of E = 0 Eλ (λ) dλ and
Eλ (λ, T ) = ελ (λ, T )Eλ,b (λ, T ), we can see that
R∞
ελ (λ, T ) Eλ,b (λ, T ) dλ
(10.26) ε(T ) = 0 .
Eb (T )
If the emissivities of a surface are known, these definitions can be applied
directly to calculate all emission characteristics.
We note (again according to definition) that the directional emissivity
of a diffuse emitter is constant, i.e. independent of direction. For “real”
surfaces, it may be said in general that I&D Ex. 12.5
• emissivity of metallic surfaces is small pp 724
• oxide layers increase the emissivity of metallic surfaces I&D Ex. 12.6
• emissivity of non–conductors is large pp 727
• emissivity of conductors increases with T , but may either decrease
or increase for non–conductors
Moreover, emissivity can depend on other factors, for example local residual
stresses or chemical interactions. Özişik (1985) and Incropera and Dewitt
(2002) contain summary data, while more comprehensive data are available
in reference texts (e.g. Wood et al., 1964).
10.3.2. Absorptivity. It follows that, like for emission, we can quan-
tify the characteristics of irradiation in terms of surface properties. At the
beginning of this chapter, we introduced the idea of a semi–transparent
medium in which a portion of the irradiation is absorbed 10.10, a portion is
reflected10.11, and a portion is transmitted10.12. For the most part, we will
concentrate on situations where transmission vanishes, so that the problem
is one of only absorption and reflection. Recall that such a material is called
opaque.
A simple energy balance for opaque surfaces shows
(10.27) Gλ = Gλ,ref + Gλ,abs + Gλ,trans ,
where Gλ is the irradiation, Gλ,ref and Gλ,abs are the components reflected
and absorbed, respectively, and Gλ,trans is the transmitted component, which
10.10Manifested as heating at the surface.
10.11Thrown back off away from the surface.
10.12Allowed to pass through to lower layers of the material.
10.3. RADIATION CHARACTERISTICS OF REAL SURFACES 129
10.13Although again, we will largely stress absorption and reflection over transmission.
10.4. KIRCHHOFF’S LAW AND GRAY SURFACES 130
Ts
Ts
Eq. (10.46) is applicable to all cases, because ε λ,θ and αλ,θ are indepen-
dent of their respective overall spectral and directional distributions. What
about the constraints on Eq. (10.45)? Are they as restrictive as those for
Eq. (10.44), i.e. that the irradiation is from a blackbody operating at the
same temperature?
Let us write Eq. (10.45) in terms of its respective definitions for ε λ in
Eq. (10.23) and αλ in Eq. (10.30), which yields
R 2π R π/2
0 0 ελ,θ (λ, θ, φ, T ) cos θ sin θ dθ dφ
R 2π R π/2 =
0 0 cos θ sin θ dθ dφ
R 2π R π/2
0 0 αλ,θ Iλ,i (λ, θ, φ) cos θ sin θ dθ dφ
R 2π R π/2 ·
0 0 I λ,i (λ, θ, φ) cos θ sin θ dθ dφ
This relation, and thus Eq. (10.45), will be satisfied for either of the following
two conditions
• If the irradiation Iλ,i (λ, θ, φ) is diffuse then Iλ,i (λ, θ, φ) → Iλ,i (λ) so
this comes outside the integrals and cancels, leaving the two sides
equal since ελ,θ is always equal to αλ,θ
• If the surface is diffuse, so that ελ,θ and αλ,θ are independent of
direction, then they come outside the integrals, in which case the
two sides are again equal
Are there any other conditions in which Eq. (10.44) will be satisfied?
As with the spectral components, let us recast Eq. (10.44) in terms of its
fundamental definitions for ε in Eq. (10.26) and α in Eq. (10.33). We find
R∞ R∞
0 ε λ (λ, T ) E λ,b (λ, T ) dλ αλ (λ) Gλ (λ) dλ
= 0 ·
Eb (T ) G
Assuming ελ = αλ , this expression is satisfied for either of the following two
conditions
• If the irradiation is the result of blackbody emission, then G λ (λ) =
Eλ,b (λ, T ) and G = Eb (T ), in which case the two sides are equal —
this was our original condition
• If the surface is gray, meaning that both ε λ and αλ are indepen-
dent of wavelength λ — they both come outside their respective
integrals, so that the expression is once again satisfied
We see now that there are a number of forms of Kirchhoff’s Law. For I&D Ex. 12.9
convenience, we summarize these forms and their restrictions in Table 10.1. pp 742
I&D Ex. 12.10
pp 744
10.4. KIRCHHOFF’S LAW AND GRAY SURFACES 133
Form Restrictions
ελ,θ = αλ,θ none
ελ = α λ ελ,θ = αλ,θ and irradiation is diffuse or surface is diffuse
ε = α ελ = αλ and thermal equilibrium and black irradiation
or gray surface
CHAPTER 11
Radiation Exchange
θj
d Aj
R
θi Tj
Aj
d Ai Ti
Ai
11.1The terms “configuration factor” and “shape factor” used in some texts are
synonymous.
134
11.1. THE VIEW FACTOR 135
From the form of Eq. (10.5) on pp. 120, we can deduce the rate at which
radiation leaves dAi and strikes dAj as
dqi→j = Ii dAi cos θi dωj,i ,
where Ii is the intensity of the radiation leaving dA i and dωj,i is the solid
angle subtended by dAj when viewed from dAi . Now recall the definition of
solid angle is dAnormal /R2 , which means that, taking the projection of dA j ,
gives
dAj cos θj
dωj,i = .
R2
Substituting, we find
cos θi cos θj
dqi→j = Ii dAi dAj .
R2
The total radiation coming from dAi is that which is emitted plus that which
is reflected, i.e. Ii should be written more exactly as Ii → Ii,e+r , where the
extra subscripts indicate emission and reflection.
Now we make an important assumption about radiation emanating from
dAi . We assume that this surface emits and reflects diffusely, so that the
total radiation is given by the radiosity relationship J i = πIi,e+r . We can
then cast dqi→j as
cos θi cos θj
(11.1) dqi→j = Ji dAi dAj .
πR2
We now make another critical assumption: J i is taken to be uniform over
surface i. Then, the total rate at which radiation leaves surface i and is
intercepted by surface j can be found simply by integrating Eq. (11.1) over
both surfaces, which gives
cos θi cos θj
Z Z
qi→j = Ji dAi dAj
Ai Aj πR2
cos θi cos θj
Z Z
= Ji dAi dAj .
Ai Aj πR2
Uniformity in Ji allows us to take it outside the integral.
Now, from the definition of view factor, the radiation which leaves sur-
face i and is intercepted by surface j is q i→j and the total radiation that
leaves surface i (emitted plus reflected) is A i Ji , therefore, by definition, the
view factor is
radiation from i that hits j qi→j
(11.2) Fij = = ·
total radiation from i Ai Ji
By inspection, we see that view factor can be calculated as I&D Ex. 13.1
1
Z Z
cos θi cos θj pp 798
(11.3) Fij = 2
dAi dAj .
A i Ai Aj πR
Eq. (11.3) is remarkable in the sense that, with the assumptions we have
made, the view factor is dependent strictly on geometry. We do not require
11.2. VIEW FACTOR ALGEBRA 136
concave
convex
flat
In cases where the surface is flat or convex, the surface does not intercept
any of its own emanating radiation. With these observations, we can write
a simple conservation law for the radiation from surface 1 as
A1 J1 = q1→1 + q1→2 + q1→3 + · · · + q1→N .
Dividing by A1 J1 , we find
A1 J1 q1→1 q1→2 q1→3 q1→N
= 1 = + + + ··· + .
A1 J1 A1 J1 A1 J1 A1 J1 A1 J1
Of course, the right hand side simply represents the sum of all the view
factors, i.e.
1 = F11 + F12 + F13 + · · · + F1N .
Note according to the above discussion that F 11 = 0 for flat and convex
surfaces, but F11 6= 0 for concave surfaces. We can write this in compact
form for surface i = 1 as
N
X
(11.6) Fij = 1 ,
j=1
applied to a given surface there are N − 1 reciprocity relations 11.2, and this
can be applied N times, thus giving N (N − 1) relationships. However, only
half of these are independent, giving N (N − 1)/2 relations 11.3. Therefore,
only N 2 − N − N (N − 1)/2 = N (N − 1)/2 view factors actually have to
be calculated directly. For example, for a 3–surface enclosure there are
32 = 9 view factors, only 3 · 2/2 = 3 of which have to be computed directly.
The other 6 factors can be calculated using the algebraic relationships in
Eqs. (11.5) and (11.6).
Example 11.1:
Consider a simple 2–surface enclosure of a sphere within a sphere
(Fig. 11.3). The smaller sphere has an outer surface area of A 1 , while the
larger has an inner surface area of A 2 . Determine the view factors charac-
terizing this problem.
There are N 2 = 22 = 4 relevant view factors, only N (N −1)/2 = 2·1/2 =
1 of which has to actually be calculated directly, i.e. by evaluating the view
factor integral. The trick is to intelligently pick the view factor that must be
calculated! For example, all radiation leaving the inner sphere is intercepted
by the outer sphere, therefore by inspection, we conclude F 12 = 1. In this
case, we skipped having to evaluate Eq. (11.3).
The rest of the problem is trivial. From reciprocity, we see
A1 A1
F21 = F12 = .
A2 A2
From summation, we see
F11 + F12 = 1 → F11 = 0 .
11.2That is, there is a reciprocity relationship to every other surface besides the one
in question. There is no meaningful reciprocity relationship between a surface and itself.
11.3An even more direct proof is to realize that this problem is also quantified by the
number of combinations of N entities taken 2 at a time. From combinatorics, we find
N! N × (N − 1) × (N − 2)! N (N − 1)
CN,2 = = = ·
(N − 2)! 2! (N − 2)! × 2 × 1 2
11.2. VIEW FACTOR ALGEBRA 139
Of course, this also could have been concluded by inspection, since the inner
sphere is completely convex. This would have been the other easy starting
point for the problem. Also from summation, we see
A1
F21 + F22 = 1 → F22 = 1 − .
A2
For more complex geometries, view factors may have to be calculated using
the double–integral definition. Solutions have been obtained for a variety of
configurations (e.g. Incropera and Dewitt, 2002, §13.1.2). ♦♦♦
Another important identity can be derived by considering the case where
a surface is broken up into its component parts (Fig. 11.4). First, it is evi-
2
1
dent that radiation reaching a composite surface is the sum of the radiation
reaching each of its parts
n
X
(11.7) Fi(j) = Fik ,
k=1
n
X
(11.8) A(j) = Ak .
k=1
For example, in Fig. 11.4, we would write F i(j) = Fi(1,2) = Fi1 + Fi2 and
A(j) = A(1,2) = A1 +A2 . Turning this concept around via reciprocity, we can
also determine a component–wise description of the surface that originates
11.2. VIEW FACTOR ALGEBRA 140
2 3
1
1
2
the cylinder, surface 1 is taken as the entire inside surface, so that again
F12 = F13 . Unlike Eqs. (11.5) and (11.6), such relations are problem specific
and will not necessarily exist for every configuration.
In retrospect, we see that N (N − 1)/2 is actually the upper bound for
the number of view factors that have to be calculated directly. Any extra
relations that we can deduce, e.g. by inspection of flat or convex surfaces
and by symmetry observations, decreases this number. Many problems are
such that the view factor integral never actually has to be evaluated.
11.4. EXCHANGE AMONG DIFFUSE GRAY SURFACES 141
“levels”, as the radiation bounces back and forth. This bouncing effect
complicates matters.
In this section, we will relax somewhat the assumption of a blackbody.
Specifically, we will model enclosures having the following properties for each
surface
• isothermal (uniform temperature)
• uniform radiosity and irradiation
• opaque (τ = 0), diffuse, and gray (ε = α)
and will then develop relations for the net rate at which radiation leaves a
surface i, defined as qi . In particular, we can write qi according to conser-
vation of energy at the radiative surface, either in terms of the difference
between radiosity and irradiation, or equivalently, between emission and
absorption (Fig. 11.6). We have
ρ *G*A G*A
J*A E*A
G*A
α*G*A
q q
(11.12) qi = A i Ji − G i
= A i Ei + ρ i Gi − G i
= A i Ei + G i ρi − 1
= A i Ei − G i 1 − ρ i
(11.13) qi = A i Ei − α i Gi .
(1 − ε i ) / ( A i ε i )
Ei Ji
qi
which a real surface differs from the ideal blackbody radiator. In other
words, if we write Eq. (11.14) as
1 − εi
qi = Ebi − Ji ,
Ai εi
it is clear that εi → 1 implies Ji → Ebi . In other words, as the emissivity
goes to unity, the radiosity approaches the blackbody emissive power — the
surface is then equivalent to the ideal blackbody surface, for which ε i = αi =
1 (perfect absorptivity) and ρi = 0 (no reflection). For εi < 1, the surface is
not ideal, so that it is described by Eq. (11.14).
Unfortunately, the results expressed thus far are incomplete because
we do not know the value of the radiosity J i in Eq. (11.14). This must
be calculated based on interaction with the other N − 1 surfaces of the
problem. However, the radiosities for these surfaces then come into play as
11.4. EXCHANGE AMONG DIFFUSE GRAY SURFACES 144
well. Recalling the developments that led to the summation rule in Eq. (11.6)
were based on the conservation of energy of the radiation leaving surface i.
Let us now consider the conservation of energy for the irradiation of surface
i, that is, the total irradiation is the sum of all net radiation transfers to
surface i. We write this as
Ai Gi = q1→i + q2→i + · · · + qN →i
= A1 J1 F1i + A2 J2 F2i + · · · + AN JN FN i (apply Eq. (11.2))
= J1 A1 F1i + J2 A2 F2i + · · · + JN AN FN i
= J1 Ai Fi1 + J2 Ai Fi2 + · · · + JN Ai FiN (reciprocity)
= Ai (J1 Fi1 + J2 Fi2 + · · · + JN FiN ) .
Canceling Ai from both sides, we can write the irradiation as
N
X
Gi = Jj Fij .
j=1
In turn, we can substitute this into Eq. (11.12) for the net rate of radiation
leaving surface i, which yields
XN
qi = Ai Ji − Jj Fij
j=1
N
X
= Ai Ji × 1 − Jj Fij (multiply Ji by 1)
j=1
N
X N
X
= Ai Ji Fij − Jj Fij (Eq. (11.6))
j=1 j=1
XN N
X
= Ai Fij Ji − Jj Fij (take Ji into Σ)
j=1 j=1
N
X
(11.15) qi = Ai Fij (Ji − Jj )
j=1
Eq. (11.15) represents the net radiation leaving surface i written in terms
of the interaction with radiosities of the other surfaces. We again notice a
clear analogy to linear circuits where J i − Jj is the potential, qi is flow,
and (Ai Fij )−1 is a resistance term11.4. Unlike the single resistor implied by
Eq. (11.14), this configuration represents a resistor for each surface interact-
ing with the given surface i (Fig. 11.8). Surface i thus comprises the node
for this circuit. Furthermore, it should be clear that in a given problem of N
11.4This quantity is usually referred to as the spatial resistance, geometric resistance,
or shape resistance.
11.5. THE GRAYBODY MATRIX PROBLEM 145
q
i −> 1
J1
i 1 ) −1
q
i −> 2
J2
i F
−1
(A
)
Fi2
q (A i J
3
i q
i −> 3
Ji
(A i Fi 3 ) −1
(A i
FiN
)
−1
JN
q
i −> N
surfaces, there will be N such circuits of the type shown in Fig. 11.8. That
is, there will be such a circuit for each node, where the nodes go from 1 to
N , each of which represents one surface of the problem.
Let us examine this equation closely. Eq. (11.16) is written for a specific
surface i, so in an actual problem, there would be N such equations for
i : 1 → N . This equation is simply the conservation of energy for each node
(surface) in the equation. In fact, it is really nothing more than Kirchhoff’s
Circuit Law , i.e. the sum of the currents (radiative transfers) for a node
must be zero. This same concept appears in Fig. 11.8.
11.5. THE GRAYBODY MATRIX PROBLEM 146
You have probably further deduced that the circuit analogy for this
problem becomes too unwieldy for large N . For example, the circuit rep-
resentation for N = 3 is shown in Fig. 11.9, but for higher values of N ,
such diagrams quickly become too difficult to reasonably draw. The matrix
Eb3
R
e3
J
R 3
s13
R
e1
R
s23
E J
b1 1
R s12
J2 R e2
Eb2
Here, we know the potential at the two termini from the Stefan–Boltzmann
(1 − ε 1 ) / (Α 1ε 1) (A 1 F1 2 ) −1 (1 − ε 2 ) / (Α 2ε 2)
q J J −q
1 E 1 2 E b2 2
b1
q
12
R
e3
J = E b3
R 3
s13
R
e1
q
13 R
s23
E J
b1 1
q
23
R s12
J2 R e2
Eb2
149
A.1. SEPARATION OF VARIABLES METHOD 150
√
tanh(x2 t + x/t) is not readily separable. We therefore suspect that
p
even this method may be limited in the types of problems that may be
solved. Fortunately, it is known to work for the type of one–dimensional
configurations we are interested in here.
Under the conjecture of Eq. (A.1), partial derivatives have certain forms.
Using the Chain Rule of Calculus, we see
∂T ∂Ψ(x) ∂Γ(t)
(A.2) = Γ(t) + Ψ(x) = Ψ(x) Γ 0 (t) ,
∂t ∂t ∂t
where the prime symbol denotes the derivative of a univariate function.
Notice that the derivative of Ψ(x) with respect to t vanishes because Ψ is
only a function of x, not of t. Similarly, we can apply Chain Rule twice to
find
∂2T
(A.3) = Ψ 00 (x) Γ(t) .
∂x2
Proceeding, we now substitute Eqs. (A.2) and (A.3) into the conduction
equation to obtain
1 Γ 0 (t) Ψ 00 (x)
(A.4) = .
α Γ(t) Ψ(x)
We have now cast the problem in a separated form where the left hand side
is only a function of time t and the right hand side is only a function of
the spatial coordinate x. However, according to principle, Eq. (A.4) must
be valid for all x and t in the problem domain. It follows that each side
must be equal to a constant. Otherwise, either of x or t could be held fixed
while the other could be varied such that Eq. (A.4) would be contradicted.
Therefore,
1 Γ 0 (t)
(A.5) = C0
α Γ(t)
and
Ψ 00 (x)
(A.6) = C0 ,
Ψ(x)
where C0 is a constant. Positive values for C 0 lead to√exponentially increas-
ing behavior and imaginary values (involving i = −1 ) lead to periodic
behavior (Carrier and Pearson, 1976). These responses can be verified by
substitution. We are instead interested in the case where C 0 is negative.
This leads to behavior that decays exponentially in time, a phenomenon
compatible with the type of boundary and initial conditions we will employ.
Therefore, we define
(A.7) C0 = −λ2 ,
where λ > 0.
A.2. SOLUTION PROCEDURE 151
RL
Next, we operate on Eq. (A.16) with 0 sin λn x dx to obtain
Z L Z LX∞
(A.17) F (x) sin λn x dx = Cm sin λm x sin λn x dx .
0 0 m=1
Note, the mode number symbol is arbitrary and we changed n → m in the
summation to avoid ambiguity. We recall that the mode coefficients are
constants rather than functions of x, so Eq. (A.17) can be simplified to
Z L X∞ Z L
(A.18) F (x) sin λn x dx = Cm sin λm x sin λn x dx .
0 m=1 0
The summation sign and coefficients have been taken outside of the integral.
Now we make the observation that the right hand side vanishes except in
the case where m = n due to the orthogonality property. This can be better
visualized if we write out the terms explicitly. We obtain
Z L
C1 sin λ1 x sin λn x dx +
0
Z L
C2 sin λ2 x sin λn x dx +
0
Z L
C3 sin λ3 x sin λn x dx + · · · +
0
Z L
Cn sin λn x sin λn x dx +
0
Z L
Cn+1 sin λn+1 x sin λn x dx +
0
Z L
Cn+2 sin λn+2 x sin λn x dx + · · ·
0
as the explicit representation of the series. Clearly, only mode C n is non–
zero, as governed by orthogonality. The rest of the modes are trivial and
the summation itself vanishes. We can then simplify Eq. (A.18) to
Z L Z L
(A.19) F (x) sin λn x dx = Cn sin λn sin λn x dx ,
0 0
which can be written using the norm as
Z L
(A.20) F (x) sin λn x dx = Cn N (λn ) .
0
Solving Eq. (A.20), we find that the mode coefficients are given by
Z L
1
(A.21) Cn = F (x) sin λn x dx .
N (λn ) 0
It is left as an exercise to the reader to show that N (λ n ) for this case is
simply L/2.
A.4. EXAMPLE: THE UNIT INITIAL CONDITION 154
Eq. (A.15) satisfies the governing equation and the boundary conditions
while Eq. (A.21) determines coefficients such that the initial condition F (x)
is satisfied. If F (x) is very complicated, evaluation of Eq. (A.21) may not
be straightforward. Otherwise, it can usually be handled with the aid of
standard integral tables. These two equations give the final solution to the
problem specified above. Physically, this means that we know the temper-
ature distribution T (x, t) for the entire problem domain 0 ≤ x ≤ L and
t ≥ 0. We can therefore compute quantities of engineering interest such as
heat flux at any position and any time.
n 1 2 3 4 5 6 7 ···
value 2 0 2 0 2 0 2 · · ·
Blackbody Radiation
Iλ,b / σT 5 Iλ,b / σT 5
λT f0→λ λT f0→λ
1100 2.7224e-06 0.000911 6400 2.3088e-05 0.76920
1200 5.2404e-06 0.00213 6500 2.2210e-05 0.77632
1300 8.8324e-06 0.00432 6600 2.1369e-05 0.78317
1400 1.3442e-05 0.00779 6700 2.0563e-05 0.78976
1500 1.8889e-05 0.01285 6800 1.9792e-05 0.79610
1600 2.4913e-05 0.01972 6900 1.9053e-05 0.80230
1700 3.1227e-05 0.02853 7000 1.8345e-05 0.80808
1800 3.7554e-05 0.03934 7100 1.7667e-05 0.81373
1900 4.3654e-05 0.05211 7200 1.7018e-05 0.81918
2000 4.9337e-05 0.06673 7300 1.6396e-05 0.82443
2100 5.4466e-05 0.08305 7400 1.5800e-05 0.82949
2200 5.8955e-05 0.10089 7500 1.5229e-05 0.83437
2300 6.2760e-05 0.12003 7600 1.4682e-05 0.83910
2400 6.5873e-05 0.14026 7800 1.3656e-05 0.84801
2500 6.8312e-05 0.16136 8000 1.2713e-05 0.85625
2600 7.0113e-05 0.18312 8500 1.0672e-05 0.87457
2700 7.1325e-05 0.20536 9000 9.0103e-06 0.88999
2800 7.2005e-05 0.22789 9500 7.6496e-06 0.90304
2898 7.2212e-05 0.25011 10000 6.5296e-06 0.91416
2900 7.2212e-05 0.25056 10500 5.6024e-06 0.92367
3000 7.2005e-05 0.27323 11000 4.8308e-06 0.93185
3100 7.1441e-05 0.29578 11500 4.1852e-06 0.93892
3200 7.0576e-05 0.31810 12000 3.6421e-06 0.94505
3300 6.9458e-05 0.34011 12500 3.1831e-06 0.95401
3400 6.8132e-05 0.36173 13000 2.7932e-06 0.95509
3500 6.6640e-05 0.38291 13500 2.4605e-06 0.95921
3600 6.5017e-05 0.40360 14000 2.1753e-06 0.96285
3700 6.3295e-05 0.42376 14500 1.9299e-06 0.96607
3800 6.1501e-05 0.44337 15000 1.7178e-06 0.96893
3900 5.9657e-05 0.46241 15500 1.5338e-06 0.97149
4000 5.7785e-05 0.48087 16000 1.3736e-06 0.97377
4100 5.5901e-05 0.49873 16500 1.2336e-06 0.97581
4200 5.4019e-05 0.51600 17000 1.1109e-06 0.97765
4300 5.2151e-05 0.53268 18000 9.0777e-07 0.98081
4400 5.0306e-05 0.54878 19000 7.4879e-07 0.98341
4500 4.8492e-05 0.56430 20000 6.2298e-07 0.98555
4600 4.6716e-05 0.57926 25000 2.7633e-07 0.99217
4700 4.4982e-05 0.59367 30000 1.4039e-07 0.99529
4800 4.3293e-05 0.60754 35000 7.8621e-08 0.99695
4900 4.1654e-05 0.62089 40000 4.7364e-08 0.99792
5000 4.0065e-05 0.63373 45000 3.0200e-08 0.99852
5100 3.8528e-05 0.64608 50000 2.0150e-08 0.99890
5200 3.7043e-05 0.65795 55000 1.3952e-08 0.99917
5300 3.5610e-05 0.66936 75000 4.1837e-09 0.99971
5400 3.4230e-05 0.68034 100000 1.3568e-09 0.99991
5500 3.2902e-05 0.69088 ∞ 0 1.0
APPENDIX D
Document History
161
About the Author
162
Creative Commons Public License
Attribution–NoDerivs–NonCommercial 1.0
THE WORK (AS DEFINED BELOW) IS PROVIDED UNDER THE TERMS OF THIS CRE-
ATIVE COMMONS PUBLIC LICENSE (“CCPL” OR “LICENSE”). THE WORK IS PROTECTED
BY COPYRIGHT AND/OR OTHER APPLICABLE LAW. ANY USE OF THE WORK OTHER THAN
AS AUTHORIZED UNDER THIS LICENSE IS PROHIBITED.
BY EXERCISING ANY RIGHTS TO THE WORK PROVIDED HERE, YOU ACCEPT AND
AGREE TO BE BOUND BY THE TERMS OF THIS LICENSE. THE LICENSOR GRANTS YOU
THE RIGHTS CONTAINED HERE IN CONSIDERATION OF YOUR ACCEPTANCE OF SUCH
TERMS AND CONDITIONS.
1. Definitions.
a. “Collective Work” means a work, such as a periodical issue, anthology or encyclopedia,
in which the Work in its entirety in unmodified form, along with a number of other con-
tributions, constituting separate and independent works in themselves, are assembled into
a collective whole. A work that constitutes a Collective Work will not be considered a
Derivative Work (as defined below) for the purposes of this License.
b. “Derivative Work” means a work based upon the Work or upon the Work and other pre-
existing works, such as a translation, musical arrangement, dramatization, fictionalization,
motion picture version, sound recording, art reproduction, abridgment, condensation, or
any other form in which the Work may be recast, transformed, or adapted, except that a
work that constitutes a Collective Work will not be considered a Derivative Work for the
purpose of this License.
c. “Licensor” means the individual or entity that offers the Work under the terms of this
License.
d. “Original Author” means the individual or entity who created the Work.
e. “Work” means the copyrightable work of authorship offered under the terms of this License.
f. “You” means an individual or entity exercising rights under this License who has not pre-
viously violated the terms of this License with respect to the Work, or who has received
express permission from the Licensor to exercise rights under this License despite a previous
violation.
2. Fair Use Rights. Nothing in this license is intended to reduce, limit, or restrict any rights
arising from fair use, first sale or other limitations on the exclusive rights of the copyright owner under
copyright law or other applicable laws.
3. License Grant. Subject to the terms and conditions of this License, Licensor hereby grants
You a worldwide, royalty-free, non-exclusive, perpetual (for the duration of the applicable copyright)
license to exercise the rights in the Work as stated below:
a. to reproduce the Work, to incorporate the Work into one or more Collective Works, and to
reproduce the Work as incorporated in the Collective Works;
b. to distribute copies or phonorecords of, display publicly, perform publicly, and perform
publicly by means of a digital audio transmission the Work including as incorporated in
Collective Works;
The above rights may be exercised in all media and formats whether now known or hereafter devised.
The above rights include the right to make such modifications as are technically necessary to exercise
the rights in other media and formats. All rights not expressly granted by Licensor are hereby reserved.
4. Restrictions. The license granted in Section 3 above is expressly made subject to and limited
by the following restrictions:
a. You may distribute, publicly display, publicly perform, or publicly digitally perform the
Work only under the terms of this License, and You must include a copy of, or the Uni-
form Resource Identifier for, this License with every copy or phonorecord of the Work You
distribute, publicly display, publicly perform, or publicly digitally perform. You may not
offer or impose any terms on the Work that alter or restrict the terms of this License or the
recipients’ exercise of the rights granted hereunder. You may not sublicense the Work. You
must keep intact all notices that refer to this License and to the disclaimer of warranties.
163
CREATIVE COMMONS PUBLIC LICENSE 164
You may not distribute, publicly display, publicly perform, or publicly digitally perform the
Work with any technological measures that control access or use of the Work in a manner
inconsistent with the terms of this License Agreement. The above applies to the Work as
incorporated in a Collective Work, but this does not require the Collective Work apart from
the Work itself to be made subject to the terms of this License. If You create a Collective
Work, upon notice from any Licensor You must, to the extent practicable, remove from the
Collective Work any reference to such Licensor or the Original Author, as requested.
b. You may not exercise any of the rights granted to You in Section 3 above in any manner
that is primarily intended for or directed toward commercial advantage or private monetary
compensation. The exchange of the Work for other copyrighted works by means of digital
file-sharing or otherwise shall not be considered to be intended for or directed toward com-
mercial advantage or private monetary compensation, provided there is no payment of any
monetary compensation in connection with the exchange of copyrighted works.
c. If you distribute, publicly display, publicly perform, or publicly digitally perform the Work
or any Collective Works, You must keep intact all copyright notices for the Work and give the
Original Author credit reasonable to the medium or means You are utilizing by conveying
the name (or pseudonym if applicable) of the Original Author if supplied; the title of the
Work if supplied. Such credit may be implemented in any reasonable manner; provided,
however, that in the case of a Collective Work, at a minimum such credit will appear where
any other comparable authorship credit appears and in a manner at least as prominent as
such other comparable authorship credit.
7. Termination.
a. This License and the rights granted hereunder will terminate automatically upon any breach
by You of the terms of this License. Individuals or entities who have received Collective
Works from You under this License, however, will not have their licenses terminated provided
such individuals or entities remain in full compliance with those licenses. Sections 1, 2, 5,
6, 7, and 8 will survive any termination of this License.
b. Subject to the above terms and conditions, the license granted here is perpetual (for the
duration of the applicable copyright in the Work). Notwithstanding the above, Licensor
reserves the right to release the Work under different license terms or to stop distributing
the Work at any time; provided, however that any such election will not serve to withdraw
this License (or any other license that has been, or is required to be, granted under the terms
of this License), and this License will continue in full force and effect unless terminated as
stated above.
8. Miscellaneous.
a. Each time You distribute or publicly digitally perform the Work or a Collective Work, the
Licensor offers to the recipient a license to the Work on the same terms and conditions as
the license granted to You under this License.
b. If any provision of this License is invalid or unenforceable under applicable law, it shall
not affect the validity or enforceability of the remainder of the terms of this License, and
without further action by the parties to this agreement, such provision shall be reformed to
the minimum extent necessary to make such provision valid and enforceable.
c. No term or provision of this License shall be deemed waived and no breach consented to
unless such waiver or consent shall be in writing and signed by the party to be charged with
such waiver or consent.
CREATIVE COMMONS PUBLIC LICENSE 165
d. This License constitutes the entire agreement between the parties with respect to the Work
licensed here. There are no understandings, agreements or representations with respect to
the Work not specified here. Licensor shall not be bound by any additional provisions that
may appear in any communication from You. This License may not be modified without
the mutual written agreement of the Licensor and You.
Bibliography
Bejan, A., Convection Heat Transfer (John Wiley & Sons, New York NY,
1984).
Burmeister, L. C., Convective Heat Transfer (John Wiley & Sons, New
York NY, 1983).
Lee, T. K., Zhong, X. L., Gong, L., and Quinn, R. (2003). Hypersonic
aerodynamic heating prediction using weighted essentially nonoscillatory
schemes. Journal of Spacecraft and Rockets, 40, 294–298.
Mills, A. F., Heat Transfer (Prentice Hall, Upper Saddle River NJ, 1999),
2nd edition.
Özişik, M. N., Heat Conduction (John Wiley & Sons, New York NY, 1980).
Panton, R. L., Incompressible Flow (John Wiley & Sons, New York NY,
1984).
Taymaz, I., Cakir, K., Gur, M., and Mimaroglu, A. (2003). Ex-
perimental investigation of heat losses in a ceramic coated diesel engine.
Surface and Coatings Technology, 169, 168–170.
Peclet number, 74
Prandtl number, 60, 86
product solution, 50
radiation, 2, 116–117
blackbody, 122
Planck’s Distribution, 123
Stefan–Boltzmann Law, 124
Wien’s Displacement Law, 124
diffuse, 121
emissivity, 126
gray body, 132
intensity, 119
Kirchhoff’s Law, 130
opaque body, 117, 128
re–radiating surface, 147
semi–transparent body, 117, 128
view factor, 134
reciprocity rule, 136
summation rule, 137
radiosity, 122
Rayleigh number, 96
resistors
in parallel, 19
in series, 19
Reynolds number, 56, 60, 78
Reynolds–Colburn analogy, 63
Sieder–Tate correlation, 88
similarity transform, 48, 67, 94
skin friction factor, 62, 64, 68, 73
solid angle, 117
stream function, 67, 95
Taylor series, 9, 28
temperature, 1
absolute and relative, 6
film, 66
gradient, 7–8
thermal
conductivity, 7, 8
contact resistance, 19
diffusivity, 8, 37