LECTURE NOTES ON

INTERMEDIATE HEAT TRANSFER

Mihir Sen

Department of Aerospace and Mechanical Engineering

University of Notre Dame
Notre Dame, IN 46556

April 2, 2000
2
Contents

Preface 9

I Preliminaries 11
1 Mathematical review 13
1.1 Fractals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.1.1 Cantor set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.1.2 Koch curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.1.3 Knopp function . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.1.4 Weierstrass function . . . . . . . . . . . . . . . . . . . . . . . 16
1.1.5 Julia set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.1.6 Mandelbrot set . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.2 Dynamical systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.2.1 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.2.2 Bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.2.3 One-dimensional systems . . . . . . . . . . . . . . . . . . . . . 19
1.2.4 Examples of bifurcations . . . . . . . . . . . . . . . . . . . . . 20
1.2.5 Unfolding and structural instability . . . . . . . . . . . . . . . 23
1.2.6 Two-dimensional systems . . . . . . . . . . . . . . . . . . . . . 23
1.2.7 Three-dimensional systems . . . . . . . . . . . . . . . . . . . . 26
1.2.8 Nonlinear analysis . . . . . . . . . . . . . . . . . . . . . . . . 29
1.3 Singularity theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

II Conduction 33
2 No spatial dimension 35

3
 2.1 Justification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.1.1 Steady state . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.1.2 Transient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2 Convective cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.1 Variable h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.2.2 Radiative cooling . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.2.3 Convective with weak radiation . . . . . . . . . . . . . . . . . 40
2.3 Radiation in an enclosure . . . . . . . . . . . . . . . . . . . . . . . . 41
2.4 Long time behavior . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.4.1 Linear analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.4.2 Nonlinear analysis . . . . . . . . . . . . . . . . . . . . . . . . 42
2.5 Time-dependent T∞ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.5.1 Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.5.2 Oscillatory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.6 Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.6.1 PID control . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.6.2 On-off control . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.7 Two-fluid problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.8 Two-body problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

3 One spatial dimensional 61

3.1 Fin theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.1.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.1.2 One-dimensional approximation . . . . . . . . . . . . . . . . . 61
3.1.3 Fin equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.2 Long time solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.3 Steady state solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.3.1 Uniform cross section . . . . . . . . . . . . . . . . . . . . . . . 65
3.3.2 Shape optimization . . . . . . . . . . . . . . . . . . . . . . . . 67
3.3.3 Annular fin . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.4 Two-dimensional fin analysis . . . . . . . . . . . . . . . . . . . . . . . 69
3.4.1 Eccentric annulus . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.5 Transient conduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4 Multiple spatial dimensions 71

4.1 Steady-state problems . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.2 Transient problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3 Stefan problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4
 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

5 Phase change 73
5.1 Stefan problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.1.1 Neumann’s solution . . . . . . . . . . . . . . . . . . . . . . . . 73
5.1.2 Goodman’s integral . . . . . . . . . . . . . . . . . . . . . . . . 74

III Convection 75
6 One-dimensional forced convection 77
6.1 Hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.1.1 Mass conservation . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.1.2 Momentum equation . . . . . . . . . . . . . . . . . . . . . . . 77
6.1.3 Long time behavior . . . . . . . . . . . . . . . . . . . . . . . . 80
6.2 Energy equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.2.1 Known heat rate . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.2.2 Known wall temperature . . . . . . . . . . . . . . . . . . . . . 82
6.3 Single duct . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
6.3.1 Steady state . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.3.2 General solution . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.3.3 Perfectly insulated duct . . . . . . . . . . . . . . . . . . . . . 86
6.3.4 Constant ambient temperature . . . . . . . . . . . . . . . . . 86
6.3.5 Periodic inlet and ambient temperature . . . . . . . . . . . . . 86
6.3.6 Effect of wall . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.4 Two-fluid configuration . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.5 Regenerator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.6 Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.6.1 Hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.6.2 Thermal networks . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.7 Thermal control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.7.1 Multiple room temperatures . . . . . . . . . . . . . . . . . . . 94
6.7.2 Two rooms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

7 One-dimensional natural convection 97

7.1 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
7.1.1 Mass conservation . . . . . . . . . . . . . . . . . . . . . . . . . 98
7.1.2 Momentum equation . . . . . . . . . . . . . . . . . . . . . . . 98
7.1.3 Energy equation . . . . . . . . . . . . . . . . . . . . . . . . . . 100

5
 7.2 Known heat rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.2.1 Steady state, no axial conduction . . . . . . . . . . . . . . . . 101
7.2.2 Axial conduction effects . . . . . . . . . . . . . . . . . . . . . 106
7.2.3 Toroidal geometry . . . . . . . . . . . . . . . . . . . . . . . . 111
7.2.4 Dynamic analysis . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.2.5 Nonlinear analysis . . . . . . . . . . . . . . . . . . . . . . . . 127
7.3 Known wall temperature . . . . . . . . . . . . . . . . . . . . . . . . . 142
7.4 Mixed condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
7.4.1 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
7.4.2 Steady State . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
7.4.3 Dynamic Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 147
7.4.4 Nonlinear analysis . . . . . . . . . . . . . . . . . . . . . . . . 152
7.5 Thermal control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173

8 Convection in porous media 175

8.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
8.1.1 Darcy’s equation . . . . . . . . . . . . . . . . . . . . . . . . . 175
8.1.2 Forchheimer’s equation . . . . . . . . . . . . . . . . . . . . . . 175
8.1.3 Brinkman’s equation . . . . . . . . . . . . . . . . . . . . . . . 176
8.1.4 Energy equation . . . . . . . . . . . . . . . . . . . . . . . . . . 176
8.2 Forced convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
8.2.1 Plane wall at constant temperature . . . . . . . . . . . . . . . 176
8.2.2 Stagnation-point flow . . . . . . . . . . . . . . . . . . . . . . . 179
8.2.3 Thermal wakes . . . . . . . . . . . . . . . . . . . . . . . . . . 179
8.3 Natural convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
8.3.1 Linear stability . . . . . . . . . . . . . . . . . . . . . . . . . . 181
8.3.2 Steady-state inclined layer solutions . . . . . . . . . . . . . . . 185
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191

9 Multidimensional forced convection 193

9.1 Low Reynolds numbers . . . . . . . . . . . . . . . . . . . . . . . . . . 193
9.2 Potential flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
9.3 Multiple solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
9.4 Plate heat exchangers . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
9.4.1 Cross flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194

6
10 Multi-dimensional natural convection 195
10.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
10.2 Cavities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
10.3 Marangoni convection . . . . . . . . . . . . . . . . . . . . . . . . . . 195
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195

11 Heat exchangers 197

11.1 Basic theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197
11.1.1 Heat transfer coefficients . . . . . . . . . . . . . . . . . . . . . 197
11.1.2 Nondimensional groups . . . . . . . . . . . . . . . . . . . . . . 197
11.1.3 Duct flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197
11.1.4 Parallel flow and counterflow . . . . . . . . . . . . . . . . . . . 197
11.1.5 Crossflow plate heat exchanger . . . . . . . . . . . . . . . . . 199
11.2 HX equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
11.2.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
11.2.2 Effectiveness-NT U relations . . . . . . . . . . . . . . . . . . . 204
11.3 Design methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 204
11.3.1 Mean temperature-difference method . . . . . . . . . . . . . . 204
11.3.2 Effectiveness-NTU method . . . . . . . . . . . . . . . . . . . . 205
11.4 Pressure drop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
11.5 Correlations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
11.6 Extended surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
11.6.1 Fin analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
11.7 Porous medium analogy . . . . . . . . . . . . . . . . . . . . . . . . . 205
11.8 Heat transfer augmentation . . . . . . . . . . . . . . . . . . . . . . . 205
11.9 Maldistribution effects . . . . . . . . . . . . . . . . . . . . . . . . . . 205
11.10Microchannel heat exchangers . . . . . . . . . . . . . . . . . . . . . . 206
11.11Radiation effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
11.12Transient behavior . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206

12 Heat transfer correlations 207

12.1 Least squares method . . . . . . . . . . . . . . . . . . . . . . . . . . . 207
12.2 Genetic algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207
12.2.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
12.2.2 Applications to compact heat exchangers . . . . . . . . . . . . 211
12.2.3 Additional applications in thermal engineering . . . . . . . . . 217
12.2.4 General discussion . . . . . . . . . . . . . . . . . . . . . . . . 220

7
 12.3 Artificial neural networks . . . . . . . . . . . . . . . . . . . . . . . . . 220
12.4 Artificial neural networks . . . . . . . . . . . . . . . . . . . . . . . . . 221
12.4.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 222
12.4.2 Application to compact heat exchangers . . . . . . . . . . . . 227
12.5 Compressible flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 236
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 236
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 236

13 Boiling 237
13.1 Boiling curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 237
13.2 Homogeneous nucleation . . . . . . . . . . . . . . . . . . . . . . . . . 237
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 237

IV Radiation 239
14 Fundamentals of radiation 241
14.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241
14.2 View factors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241

15 Computational methods 243

15.1 Monte Carlo methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 243
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243

V Appendices 245
A Routh-Hurwitz criteria 247

Bibliography 249

8
Preface

These are lecture notes for ME445/AME545: Intermediate Heat Transfer, a second
course on heat transfer for undergraduate seniors and beginning graduate students.
In addition to some undergraduate knowledge of heat transfer, students taking this
course are expected to be familiar with vector algebra, linear algebra, ordinary dif-
ferential equations, particle and rigid-body dynamics, thermodynamics, and integral
and differential analysis in fluid mechanics. The use of computers is essential both
for the purpose of computation as well as for display and visualization of results.
At present these notes are in the process of being written; the student is encour-
aged to make extensive use of the literature listed in the bibliography. The students
are also expected to attempt the problems at the end of each chapter to reinforce
their learning.
I will be glad to receive comments on these notes, and have mistakes brought to
my attention.

Mihir Sen
Department of Aerospace and Mechanical Engineering
University of Notre Dame

Copyright c by M. Sen, 1999

9
10
 Part I

Preliminaries

11
Chapter 1

Mathematical review

1.1 Fractals
Fractals are objects that are not smooth; they are geometrical shapes in which the
parts are in some way similar to the whole. This self-similarity may be exact, i.e. a
piece of the fractal, if magnified, may look exactly like the whole fractal.
A function f (x) is invariant under change of scale if there exists constants a and
b, such that
f (ax) = bf (x) (1.1)
A fractal curve must be nowhere rectifiable (i.e. any part of it cannot be of finite
length) and homogeneous (i.e. any par6 is similar to the whole).
Before discussing examples we need to put forward a working definition of di-
mension. Though there are many definitions in current use, we present here the
Hausdorff-Besicovitch dimension D. If N is the number of ‘boxes’ of side  needed
to cover an object, then
ln N
D = lim (1.2)
→0 ln(1/)

We can check that this definition corresponds to the common geometrical shapes.
1. Point: N = 1, D = 0

2. Line of length l: N = l/, D = 1

3. Surface of size l × l: N = (l/)2 , D = 2

4. Volume of size l × l × l: N = (l/)3 , D = 3

A fractal has a dimension that is not an integer. Many physical objects are
fractal-like, in that they are fractal within a range of length scales. Coastlines are

13
 k=0
k=1
k=2
k=3

Figure 1.1: The Cantor set.

among the geographical features that are of this shape. If there are N units of
a measuring stick of length , the measured length of the coastline will be of the
power-law form N = 1−D , where D is the dimension.

1.1.1 Cantor set

Consider the line [0,1] corresponding to k = 0 in Figure 1.1. Take away the middle
third to leave the two portions; this is shown as k = 1. Repeat the process to get
k = 2, 3, . . .. If k → ∞, what is left is called the Cantor set. Since N = 2k and
 = 1/3k ,its dimension is D = ln 2/ ln 3 = 0.63 . . ..
If we define a function Ck (t) at the kth level so that Ck (t) = (3/2)k if t belongs
to the set and zero otherwise, its integral over the interval is unity. In terms of this
function we can also define
Z t
D(t) = lim Ck (t0 ) dt0 (1.3)
0 k→∞

that is called the devil’s staircase.

1.1.2 Koch curve

Here we start with an equilateral triangle shown in Figure 1.2 as k = 0. The middle
third of each side of the triangle is removed, and two sides of a triangle drawn on
that. This is shown as k = 1. The process is continued, and in the limit gives a
continuous, closed curve that is nowhere smooth. Since N = 3 × 4k and  = 1/3k ,
the dimension of the Koch curve is D = ln 4/ ln 3 = 1.26 . . ..

1.1.3 Knopp function

This is the function
∞
X
K(t) = 2−nH g(2n t) (1.4)
n=0

14
k=0 k=1

k=2

Figure 1.2: The Koch curve.

15
where 0 < H < 1, and g(t) is the periodic triangular function
(
2t for 0 ≤ t ≤ 1/2
g(t) = (1.5)
2(1 − t) for 1/2 < t ≤ 1

defined on [0,1].

1.1.4 Weierstrass function

This is the function
∞
X
W (t) = ω −nH cos (ω n t + φn ) (1.6)
n=0

where a is real, b is odd, and ab > 1 + 3π/2. It is everywhere continuous, but nowhere
differentiable. The related Weierstrass-Mandelbrot function
∞
X
Wm (t) = ω −nH (1 − cos ω n t) (1.7)
n=−∞

satisfies the inveriance relation 1.1.

1.1.5 Julia set

An example of this comes from the application of Newton’s method to find the com-
plex root of the equation z 3 = 1. In this method the following iterative scheme is set
up:
z3 − 1
zk+1 = zk − k 2 (1.8)
3zk
Each one of the three roots has a basin of attraction, the boundaries of which are
fractal.

1.1.6 Mandelbrot set

This is the set of complex numbers c for which

zk+1 = zk2 + c (1.9)

stays bounded as k → ∞. The boundaries of this set shown in Figure 1.3 are again
fractal.

16
 Figure 1.3: Mandelbrot set

1.2 Dynamical systems

A dynamical system is a set of differential equations such as

dxi
= fi (x1 , x2 , . . . , t; λ1 , λ2 , . . . , λp ) for i = 1, . . . , n (1.10)
dt
The x1 , . . . , xn s are state variables and the λ1 , . . . , λp are bifurcation parameters. The
mapping f : X × Rp → Y is a vector field. If f1 , . . . , fn do not depend on time t, the
system is autonomous. A nonautonomous system can be converted into autonomous
one by the change in variable xn+1 = t, from which we can get the additional equation

dxn+1
=1 (1.11)
dt
We will assume that the solutions of the system are always bounded. The system is
P
conservative if the divergence of the vector field ∂fi /∂xi is zero, and dissipative if
it is negative. An attractor of a dissipative dynamical system is the set of {xi } as
t → ∞. The critical (or singular, equilibrium or fixed) points, xi , of equation (1.10)
are those for which
fi (x1 , x2 , . . . , t; λ1 , λ2 , . . . , λp ) = 0 (1.12)
There may be multiple solutions to this algebraic or transcendental equation.
Defining a new coordinate x0i = xi − xi that is centered at the critical point, we
get the local form
dx0i
= fi (x1 − x1 , . . . , xn − xn ) (1.13)
dt
Sometimes we will use the notation
dxi
= g(x1 , . . . , xn ) (1.14)
dt
to indicate the local form, the origin being one critical point of this system.

1.2.1 Stability
The stability of the critical points is of major interest. A critical point is stable if,
given an initial perturbation, the solutions tends to it as t → ∞.

17
Linear stability

The vector field in equation (1.14) can be expanded in a Taylor series to give

dxi X ∂gi
= xj + . . . (1.15)
dt j ∂xj 0

The eigenvalues of the Jacobian matrix

∂gi
A= (1.16)
∂xj

determine the linear stability of the critical point. The critical point is stable if all
eigenvalues have negative real parts, and unstable if one or more eigenvalues have
positive real parts.

Global stability

Consider the dynamical system in local form, equation (1.14). If there exists a func-
tion V (x1 , . . . , xn ) such that V ≥ 0 and dV /dt ≤ 0, then the origin is globally stable,
that is, it is stable to all perturbations, large or small. V is called a Liapunov function.

1.2.2 Bifurcations
The critical point is one possible attractor. There are other time-dependent solutions
which can also be attractors in phase space, as indicated in the list below.

• Point (steady, time-independent)

• Closed curve (limit cycle, periodic)

• Torus (periodic or quasi-periodic)

• Strange (chaotic)

For a given dynamical system, several attractors may co-exist. In this case each
attractor has a basin of attraction, i.e. the set of initial conditions that lead to this
attractor. A bifurcation is a qualitative change in the solution as the bifurcation
parameters λi are changed.

18
1.2.3 One-dimensional systems
A one-dimensional dynamical system is of the type
dx
= f (x) (1.17)
dt

Example 1.1
Consider the linear equation
dx
= ax + b (1.18)
dt
The critical point is
b
x=− (1.19)
a
On defining x0 = x − x, the local form is obtained as
dx0
= ax0 (1.20)
dt
We can take one of two approaches.

(a) Solving
x0 = x00 eat (1.21)
The critical point is a repeller if a > 0, and an attractor if a < 0.

(b) Alternatively, we can multiply equation (1.20) by 2x0 to get

dV
= 2aV (1.22)
dt
where V = x02 . Since V is always nonnegative, the sign of dV /dt is the sign of a. Thus V will
increase with time if a > 0, and decrease if a < 0.
In either case we find that the critical point is unstable if a > 0 and stable if a < 0.

Example 1.2
The nonlinear equation
dx  
= −x x2 − (λ − λ0 ) (1.23)
dt
has critical point which are solutions of the cubic equation
 
x x2 − (λ − λ0 ) = 0 (1.24)

19
 Thus
x(1) = 0 (1.25)
p
x(2) = λ − λ0 (1.26)
p
x(3) = − λ − λ0 (1.27)
where x(i) (i = 1, 2, 3) are the three critical points. The bifurcation diagram is shown in Fig.
1.4.

(i) Critical point x(1) = 0

To analyze the local stability of x(1) = 0, we obtain the local form
dx0
= x0 (λ − λ0 ) − x03 (1.28)
dt
Neglecting the cubic term x03 , this becomes
dx0
= x0 (λ − λ0 ) (1.29)
dt
Thus x(1) = 0 is locally stable if λ < λ0 , and unstable if λ > λ0 .
To analyze the global stability, equation (1.28) can be written as
1 dV
= V (λ − λ0 ) − V 2 (1.30)
2 dt
where V = x02 . For λ < λ0 , V ≥ 0, dV /dt ≤ 0, so that x(1) = 0 is globally stable.
√
(ii) Critical point x(2) = λ − λ0
The local form of the equation around this critical point is
dx0 p   p 2 
=− λ − λ0 + x0 λ − λ0 + x0 − (λ − λ0 ) (1.31)
dt
Linearizing, we get
dx0
= −2(λ − λ0 )x0 (1.32)
dt
so that this critical point is linearly stable.
√
(iii) Critical point x(3) = − λ − λ0
This is similar to the above.

1.2.4 Examples of bifurcations

• Supercritical: Fig. 1.4.
• Subcritical: Fig. 1.5.
• Transcritical: Fig. 1.6.
• Saddle-node: Fig. 1.7.

20
x
stable

stable unstable
λ λ
0

stable

Figure 1.4: Supercritical pitchfork bifurcations

x
stable stable

S-N unstable
unstable unstable
stable
λ

linearly unstable
stable
stable

Figure 1.5: Subcritical version of Fig. 1.4

21
x

stable

stable unstable
λ

unstable

Figure 1.6: Transcritical bifurcation.

x
stable

S-N
unstable
λ

Figure 1.7: Saddle-node bifurcation.

22
1.2.5 Unfolding and structural instability
Adding a small constant to the vector field in equation (1.23), we get
h i
f (x) = −x x2 − (λ − λ0 ) +  (1.33)

The critical points are solutions of

x3 − (λ − λ0 )x −  = 0 (1.34)

To see the nature of the curve, we make an expansion around x = 0, λ = λ0 and

write

x = x0 (1.35)
λ = λ0 + λ0 (1.36)

For small x0 and λ, we get

λ0 x0 = − (1.37)
Figure 1.8 shows the result of adding the imperfection . The dynamical system
without  is thus structurally unstable.

1.2.6 Two-dimensional systems

Consider
dx
= fx (x, y) (1.38)
dt
dy
= fy (x, x) (1.39)
dt
The linearized equation in local form has a Jacobian matrix
!
axx axy
A= (1.40)
ayx ayy

The eigenvalues satisfy a quadratic equation

λ2 + P λ + Q = 0 (1.41)

from which q 
1
λ= −P ± P 2 − 4Q (1.42)
2
The sign of the discriminant
D = P 2 − 4Q (1.43)

23
x
stable

stable unstable

λ0 λ

stable

x
stable

stable λ
unstable
λ0

stable

Figure 1.8: Imperfect bifurcations; (a)  > 0, (b)  < 0.

24
determines the nature of the solution. If D < 0, the eigenvalues are complex and the
solution in phase space is a spiral; if in addition P > 0, the spiral is stable, and if
P < 0, it is unstable. If, on the other hand, D > 0, the eigenvalues are real; the
solutions do not oscillate in time but move exponentially towards (if all eigenvalues
are negative) or away (if at least one eigenvalue is positive) from the critical point.

Example 1.3

dx
= y (1.44)
dt
dy
= −(λ − λ0 )x (1.45)
dt
This is a conservative system which is equivalent to

d2 x
+ (λ − λ0 )x = 0 (1.46)
dt2
The solutions are exponential if λ < λ0 , and periodic if λ > λ0 .

Example 1.4

dx
= y (1.47)
dt
dy
= −ω 2 x − σy (1.48)
dt
For σ > 0, the system is dissipative, and the solutions are damped oscillations.

The occurrence of periodic solutions in two-dimensional systems, permits a Hopf

bifurcation, which is a transition from a time-independent to a periodic behavior
through a pair of complex conjugate imaginary eigenvalues.

Example 1.5
The dynamical system
dx
= (λ − λ0 )x − y − (x2 + y 2 )x (1.49)
dt
dy
= x + (λ − λ0 )y − (x2 + y 2 )y (1.50)
dt

25
 can be converted to polar coordinates. Substituting x = r cos θ and y = r sin θ, we get
dθ dr
−r sin + cos θ = (λ − λ0 )r cos θ − r sin θ − r3 cos θ (1.51)
dt dt
dθ dr
r cos + sin θ = r cos θ + (λ − λ0 )r sin θ − r3 sin θ (1.52)
dt dt
which simplifies to
dr

= r λ − λ0 − r2 (1.53)
dt
dθ
= 1 (1.54)
dt
√
There are two values of r, i..e r = 0 and r = λ − λ0 , at which dr/dt = 0. The first is a critical
point at the origin, and the second a circular periodic orbit that exists only for λ > λ0 . A
linear analysis of equations (1.49) and (AHopftwo) shows that the origin
√ is stable for λ < λ0 .
For λ > λ0 , a similar analysis of equation (1.53) indicates that r = λ − λ0 is a stable orbit.
There is thus a Hopf bifurcation at λ = λ0 .

Example 1.6

dx
= y (1.55)
dt
dy

= −ω 2 x − σ λ − x2 − y 2 (1.56)
dt
This is a Hopf bifurcation at λ = λ0 , at which point the solution goes from time-independent
to periodic.

1.2.7 Three-dimensional systems

Forced Duffing equation
Since a non-autonomous equation can be converted into a three-dimensional au-
tonomous system, we will include the case of the forced Duffing equation here1 .
dx
= y (1.57)
dt
dy
= x − x3 + f (t) (1.58)
dt
1
Sometimes defined with a negative sign in front of x in the second equation.

26
Lorenz equations
An important example is the Lorenz equations:
dx
= σ(y − x) (1.59)
dt
dx
= λx − y − xz (1.60)
dt
dx
= −bz + xy (1.61)
dt
where σ and b are taken to be positive constants, with σ > b + 1. The bifurcation
parameter will be λ.
The critical points are obtained from

y−x = 0
λx − y − xz = 0
−bz + xy = 0

which give
     q   q 
x 0 b(λ − 1) − b(λ − 1)
     q   q 
 y  =  0   b(λ − 1)   − b(λ − 1) 
,   , 
 (1.62)
z 0 λ−1 λ−1

A linear stability analysis of each critical point follows.

(a) x = y = z = 0
Small perturbations around this point give
    
x0 −σ σ 0 x0
d  0    0 
 y  =  λ −1 0   y  (1.63)
dt
z0 0 0 −b z0

The characteristic equation is

(λ + b)[λ2 + λ(σ + 1) − σ(λ − 1)] = 0 (1.64)

q
from which we get the eigenvalues −b, 12 [−(1 + σ) ± (1 + σ)2 − 4σ(1 − λ)]. For
0 < λ < 1, the eigenvalues are real and negative, since (1 + σ)2 > 4σ(1 − λ). At
λ = λ1 , where λ1 = 1, there is a pitchfork bifurcation with one zero eigenvalue. For
λ > λ1 , the origin becomes unstable.

27
 q
(b) x = y = b(λ − 1), z = λ − 1
Small perturbations give
    
x0 −σ σ q
0 x0
d  0      0 
 y = q 1 −1 − b(λ − 1) 
  y  (1.65)
dt q
z0 b(λ − 1) b(λ − 1) −b z0

The characteristic equation is

λ3 + (σ + b + 1)λ2 + (σ + λ)bλ + 2σb(λ − 1) = 0 (1.66)

Using the Hurwitz criteria we can determine the sign of the real parts of the solutions
of this cubic equation without actually solving it. The Hurwitz determinants are

D1 = σ + b + 1
σ + b + 1 2σb(λ + 1)
D2 =
1 (σ + λ)b
= σb(σ + b + 3) − λb(σ − b − 1)
σ + b + 1 2σb(λ − 1) 0
D3 = 1 (σ + λ)b 0
0 σ + b + 1 2σb(λ − 1)
= 2σb(λ − 1)[σb(σ + b + 3) − λb(σ − b − 1)]

Thus the real parts of the eigenvalues are negative if λ < λ3 , where
σ(σ + b + 3)
λ3 = (1.67)
σ−b−1
At λ = λ3 the characteristic equation (1.66) can be factorized to give the eigenvalues
−(σ + b + 1), and ±i 2σ(σ + 1)/(σ − b − 1), corresponding to a Hopf bifurcation. The
periodic solution which is created at this value of λ can be shown to be unstable so
that the bifurcation is subcritical.
Here is a summary of the series of bifurcation with respect to the parameter λ:
• Origin is a stable critical point for λ < λ1 ; becomes unstable at λ = λ1 .
• Two other critical points are created for λ > λ1 ; these are linearly stable in the
range λ1 < λ < λ3 .
• Just below λ3 , i.e. in the range λ2 < λ < λ3 , the two critical points are stable
to small perturbations, but for large enough perturbations produce chaos.
• For λ > λ3 , all initial conditions produce chaos (except for periodic windows).

28
1.2.8 Nonlinear analysis
Center manifold theorem
Consider a vector field fi (x) with fi (0) = 0. The eigenvalues λ of ∂fi /∂xj at the
origin are of three kinds:
(a) Re(λ) > 0 with the generalized eigenspace E u .
(b) Re(λ) < 0 with the generalized eigenspace E s .
(c) Re(λ) = 0 with the generalized eigenspace E c .
There exist manifolds W u , W s , and W c to which E u , E s and E c , respectively, are
tangents. W u , W s , and W c are the unstable, stable and center manifolds, respectively.

1.3 Singularity theory

We have seen that the critical points of a dynamical system, equation (1.10), are
found by solving an equation of the type (1.12), i.e. by finding the singularities of the
function fi . The type of bifurcations that occur with a single bifurcation parameter λ
has been discussed in the previous section. Here we increase the number of parameters
to λ1 , . . . , λp and look at the consequent changes in the bifurcations. A bifurcation (or
catastrophe) set is the set of points in parametric space λ1 , . . . , λp at which equation
(1.12) is satisfied.

Example 1.7
Examine the one-dimensional vector field

f (x) = −(x3 + px + q) (1.68)

Figure 1.9 shows the surface x = x(λ, µ). There are three real solutions if the discriminant
D = 4λ3 + 27µ2 > 0, and only one otherwise. A section of this surface at µ = 0 will give Fig.
1.4 and µ =  will give Fig. 1.8.
The bifurcation set is shown in Fig. 1.10 in (λ, µ) coordinates. It is projection of the
x = x(λ, µ) surface in the (λ, µ) plane.

The quadratic surface

a1 x2 + a2 y 2 + a3 z 3 + a4 xy + a5 xz + a6 yz + a7 x + a8 y + a9 z + a10 = 0 (1.69)

can be classified in terms of eleven canonical surfaces. For gradient systems, i.e.
systems in which fI = ∂φ/∂xi , there are only seven.

29
x

Figure 1.9: Surface x = x(λ, µ).

30
µ

One solution

Three solutions

One solution

Figure 1.10: Bifurcation set.

31
References
1. Tricot, C., Curves and Fractal Dimension, Springer-Verlag, New York, 1995..

Problems
1. Show that
dx
= −x [x − (λ − λ0 )]
dt
has a transcritical bifurcation.
2. Carry out an imperfection analysis on the transcritical bifurcaiton above.
3. Show that
dx
= −x2 + (λ − λ0 )
dt
has a saddle-node bifurcation.
4. Investigate the bifurcations in
dx

= −x x4 − 2x2 + 2 + λ
dt

32
 Part II

Conduction

33
Chapter 2

No spatial dimension

2.1 Justification

T∞,1

Tw,1

Tw,2

T∞,2

Figure 2.1: Wall with fluids on either side.

Consider a wall with fluid on both sides as shown in Fig. 2.1. The fluid temper-
atures are T∞,1 and T∞,2 and the wall temperatures are Tw,1 and Tw,2 . The initial
temperature in the wall is T (x, 0) = f (x).

35
2.1.1 Steady state
In the steady state, we have
Tw,1 − Tw,2
h1 (T∞,1 − Tw,1 ) = ks = h2 (Tw,2 − T∞,2 ) (2.1)
L
from which
h1 L T∞,1 − Tw,1 Tw,1 − Tw,2 h2 L Tw,2 − T∞,2
= = (2.2)
ks T∞,1 − T∞,2 T∞,1 − T∞,2 ks T∞,1 − T∞,2
Thus we have
Tw,1 − Tw,2 T∞,1 − Tw,1 h1 L
 if 1 (2.3)
T∞,1 − T∞,2 T∞,1 − T∞,2 ks
Tw,1 − Tw,2 Tw,2 − T∞,2 h2 L
 if 1 (2.4)
T∞,1 − T∞,2 T∞,1 − T∞,2 ks

The Biot number is defined as

hL
Bi = (2.5)
ks

2.1.2 Transient
∂T ks ∂ 2 T
= (2.6)
∂t ρc ∂x2
There are two time scales: the short (conductive) tk0 = L2 ρc/ks and the long (con-
vective) th0 = Lρc/h. In the short time scale conduction within the slab is important,
and convection from the sides is not. In the long scale, the temperature within the
slab is uniform, and changes due to convection. The ratio of the two tk0 /th0 = Bi. In
the long time scale it is possible to show that
dT
Lρs c + h1 (T − T∞,1) + h2 (T − T∞,2) = 0 (2.7)
dt
where T = Tw,1 = Tw,2.

2.2 Convective cooling

A body at temperature T , such as that shown in Fig. 2.2, is placed in an environ-
ment of different temperature, T∞ , and is being convectively cooled. The governing
equation is
dT
Mc + hA(T − T∞ ) = 0 (2.8)
dt
36
 T∞
T

Figure 2.2: Convective cooling.

with T (0) = Ti . We nondimensionalize using

T − T∞
θ = (2.9)
Ti − T∞
hAt
τ = (2.10)
Mc
The nondimensional form of the governing equation (2.8) is
dθ
+θ =0 (2.11)
dτ
the solution to which is
θ = e−τ (2.12)
This is shown in Fig. 2.3 where the nondimensional temperature goes from θ = 1 to
θ = 0. The dimensional time constant is Mc/hA.

2.2.1 Variable h
If, however, the h is slightly temperature-dependent, then we have
dθ
+ (1 + θ)θ = 0 (2.13)
dτ
which can be solved by the method of perturbations. We assume that
θ(τ ) = θ0 (τ ) + θ1 (τ ) + 2 θ2 (τ ) + . . . (2.14)
To order 0 , we have
dθ0
+ θ0 = 0 (2.15)
dτ
θ0 (0) = 1 (2.16)

37
 θ

Figure 2.3: Convective cooling.

which has the solution

θ0 = e−τ (2.17)
To the next order 1 , we get
dθ1
+ θ1 = −θ02 (2.18)
dτ
= −e−2τ (2.19)
θ1 (0) = 0 (2.20)
the solution to which is
θ1 = −e−t + e−2τ (2.21)
Taking the expansion to order 2
dθ2
+ θ2 = −2θ0 θ1 (2.22)
dτ
= −2e−2t − 2e−3τ (2.23)
θ2 (0) = 1 (2.24)
with the solution
θ2 = e−τ − 2e−2τ + e−3τ (2.25)
And so on. Combining, we get
θ = e−τ − (e−τ − e−2τ ) + 2 (e−τ − 2e−2τ + e−3τ ) + . . . (2.26)

38
 Alternatively, we can find an exact solution to equation (2.13). Separating vari-
ables, we get
dθ
= dτ (2.27)
(1 + θ)θ
Integrating
θ
ln = −τ + C (2.28)
θ + 1
The condition θ(0) = 1 gives C = − ln(1 + ), so that

θ(1 + )
ln = −τ (2.29)
1 + θ
This can be rearranged to give

e−τ
θ= (2.30)
1 + (1 − e−τ )

A Taylor-series expansion of the exact solution gives

h i−1
θ = e−τ 1 + (1 − e−τ (2.31)
h i
= e−τ 1 − (1 − e−τ ) + 2 (1 − e−τ )2 + . . . (2.32)
= e−τ − (e−τ − e−2τ ) + 2 (e−τ − 2e−2τ + e−3τ ) + . . . (2.33)

2.2.2 Radiative cooling

If the heat loss is due to radiation, we can write

dT
Mc + σA(T 4 − T∞
4
)=0 (2.34)
dt
Taking the dimensionless temperature to be defined in equation (2.9), and time to be

σA(Ti − T∞ )3 t
τ= (2.35)
Mc
and introducing the parameter
T∞
β= (2.36)
Ti − T∞
we get
dφ
+ φ4 = β 4 (2.37)
dτ
39
where
φ=θ+β (2.38)
Writing the equation as
dφ
= −dτ (2.39)
φ4 − β4
the integral is ! !
1 φ−β 1 φ
ln − 3 tan−1 = −τ + C (2.40)
4β 3 φ+β 2β β
Using the initial condition θ(0) = 1, we get (?)
" #
1 1 (β + T )(β − 1) T −1
τ= 3 ln + tan−1 (2.41)
2β 2 (β − T )(β + 1) β + (T /β)

2.2.3 Convective with weak radiation

The governing equationis
dT
Mc + hA(T − T∞ ) + σA(T 4 − T∞
4
)=0 (2.42)
dt
with T (0) = Ti . Using the variables defined by equations (2.9) and (2.10), we get
dθ h i
+ θ +  (θ + β 4 )4 − β 4 = 0 (2.43)
dτ
where β is defined in equation (2.36), and
σ(Ti − T∞ )3
= (2.44)
h
If radiative effects are small compared to the convective (for Ti − T∞ = 100 K and
h = 10 W/m2 K we get  = 5.67 × 10−3 ), we can take   1. Substituting the
perturbation series, equation (2.14), in equation (2.43), we get
d    
θ0 + θ1 + 2 θ2 + . . . + θ0 + θ1 + 2 θ2 + . . .
dτ
 4  3
+[ θ0 + θ1 + 2 θ2 + . . . + 4β θ0 + θ1 + 2 θ2 + . . .
 2  
+6β 2 θ0 + θ1 + 2 θ2 + . . . 4β 3 θ0 + θ1 + 2 θ2 + . . . ] = 0 (2.45)
In this case
dθ
+ (θ − θ0 ) + (θ4 − θs4 ) = 0 (2.46)
dτ
θ(0) = 1 (2.47)

40
 As a special case, of we take β = 0, i.e. T∞ = 0, we get
dθ
+ θ + θ4 = 0 (2.48)
dτ
which has an exact solution
1 1 + θ3
τ= ln (2.49)
3 (1 + )θ3

2.3 Radiation in an enclosure

Consider a closed enclosure with N walls radiating to each other and with a central
heater H. The walls have no other heat loss and have different masses and specific
heats. The governing equations are

dTi XN
Mi ci +σ Ai Fij (Ti4 − Tj4 ) + σAi FiH (Ti4 − TH4 ) = 0 (2.50)
dt j=1

where the view factor Fij is the fraction of radiation leaving surface i that falls on j.
The steady state is
T i = TH (i = 1, . . . , N) (2.51)
Linear stability is determined by a small perturbation of the type

Ti = TH + Ti0 (2.52)

from which
dTi0 3
XN
Mi ci + 4σTH Ai Fij (Ti0 − Tj0 ) + 4σTH3 Ai FiH Ti0 = 0 (2.53)
dt j=1

This can be written as

dT0
M = −4σTH3 AT0 (2.54)
dt

2.4 Long time behavior

The general form of the equation for heat loss from a body is
dθ
+ f (θ) = a (2.55)
dτ
with θ(0) = 1. Let
f (θ) = a (2.56)

41
Then we would like to show that θ → θ as t → ∞. Writing

θ = θ + θ0 (2.57)

we have
dθ0
+ f (θ + θ0 ) = a (2.58)
dτ

2.4.1 Linear analysis

If we assume that θ0 is small, then a Taylor series gives

f (θ + θ0 ) = f (θ) + θ0 f 0 (θ) + . . . (2.59)

from which
dθ0
+ bθ0 = 0 (2.60)
dτ
where b = f 0 (θ). The solution is
θ0 = Ce−bτ (2.61)
so that θ0 → 0 as t → ∞ if b > 0.

2.4.2 Nonlinear analysis

Multiplying equation (2.58) by θ0 , we get

1 d 0 2 h i
(θ ) = θ0 f (θ + θ0 ) − f (θ) (2.62)
2 dτ
Thus
d 0 2
(θ ) ≤ 0 (2.63)
dτ
if θ0 and [f (θ + θ0 ) − f (θ)], as shown in Fig. 2.4, are both of the same sign or zero.

2.5 Time-dependent T∞
Let
dT
Mc + hA(T − T∞ (t)) = 0 (2.64)
dt
with
T (0) = Ti (2.65)

42
 f(θ )

f(θ )

θ
Figure 2.4: Convective cooling.

2.5.1 Linear
Let
T∞ = T∞,0 + at (2.66)
Defining the nondimensional temperature as
T − T∞,0
θ= (2.67)
Ti − T∞,0

and time as in equation (2.10), we get

dθ
+ θ = Aτ (2.68)
dτ
where
aMc
A= (2.69)
hA(Ti − T∞,0 )
The nondimensional ambient temperature is
aMc
θ∞ = τ (2.70)
hA(Ti − T∞,0 )

43
The solution to equation (2.68) is

θ = Ce−τ + Aτ − A (2.71)

The condition θ(0) = 1 gives C = 1 + A, so that

θ = (1 + A)e−τ + A(τ − 1) (2.72)

The time shown in Fig. 2.5 at crossover is

1+A
τc = ln (2.73)
A
and the offset is
δθ = A (2.74)
as τ → ∞.

∆θ

θ
∞

τc

Figure 2.5: Response to linear ambient temperature.

44
2.5.2 Oscillatory
Let
T∞ = T ∞ + δT sin ωt (2.75)
where T (0) = Ti . Defining
T − T∞
θ= (2.76)
Ti − T ∞
and using equation (2.10), the nondimensional equation is

dθ
+ θ = δθ sin Ωτ (2.77)
dτ
where
δT
δθ = (2.78)
Ti − T ∞
ωMc
Ω = (2.79)
hA
The solution is
δθ
θ = Ce−τ + sin(Ωτ − φ) (2.80)
(1 + Ω2 ) cos φ
where
φ = tan−1 Ω (2.81)
From the condition θ(0) = 1, we get C = 1 + δΩ/(1 + Ω2 ), so that
 
Ω δθ
θ = 1 + δθ 2
e−τ + sin(Ωτ − φ) (2.82)
1+Ω (1 + Ω2 ) cos φ

2.6 Control
A lumped control system can be represented by
dx
= f (x, u, d) (2.83)
dt
y = h(x, u) (2.84)

where x(t) is the state vector, y(t) is the output vector, d(t) is the disturbance vector,
and u(t) is the control input. Usually u(t) is related to the error e(t) = xs − x(t),
where xs is the desired state. If xs is a function of time, the problem is one of tracking,
but if it is a constant, it is regulation.

45
 For a linear system
dx
= Ax + Bu + Γd (2.85)
dt
y = Cx + Du (2.86)

where A, B, C and D are matrices of appropriate sizes. For a distributed system,

equation (2.83) is replaced by
∂=
(x, u, ∇ · x, ∇2 x, . . .) (2.87)
∂f
We can use internal heat generation, Q(t), to control the temperature of a body
losing heat to the environment which is at a temperature T∞ . The governing equation
is
dT
Mc + hA(T − T∞ ) = Q(t) (2.88)
dt
where T (0) = Ti . Nondimensionalizing the temperature and time as in equations
(2.9) and (2.10), we get
dθ
+ θ = q(τ ) (2.89)
dτ
with θ(0) = 1, where
Q
q= (2.90)
hA(Ti − T∞ )
If we take q = q0 to be a constant, the solution is

θ = (1 − q0 )e−τ + q0 (2.91)

As τ → ∞, θ → q0 asymptotically.
We can design a feedback controller for the lumped system as shown in Fig. 2.6
to attempt to maintain the temperature of the mass at a set temperature θs .

2.6.1 PID control

We will first investigate proportional-integral-derivative (PID) control. The control
input is taken to be Z τ
de
q = K p e + Ki e(τ 0 ) dτ 0 + Kd (2.92)
0 dτ
where the error is
e(τ ) = θs − θ (2.93)

46
 +
Controller Plant
θs − θ

Figure 2.6: Feedback control system.

Proportional control
Let
q = Kp (θs − θ) (2.94)
from which we get
dθ
+ θ = Kp (θs − θ) (2.95)
dτ
The solution is
Kp θs
θ = Ce−(1+Kp )τ + (2.96)
1 + Kp
From the initial condition θ(0) = 1, we have C = 1 − Kp θs /(1 + Kp ), so that
!
Kp θs Kp θs
θ = 1− e−(1+Kp )τ + (2.97)
1 + Kp 1 + Kp
If 1 + Kp > 0, we find that as τ → ∞
Kp
θ(∞) = θs (2.98)
1 + Kp
The set point θs is thus never achieved, and there is an offset θoffset defined by
θ(∞) = θs − θoffset (2.99)
where " #
Kp
θoffset = θs 1 − (2.100)
1 + Kp
the offset can be reduced by choosing a large Kp . For large Kp , the offset is approxi-
mated by
θs
θoffset = (2.101)
Kp

47
Integral control
In this case Z τ
q = Ki (θs − θ(τ 0 )) dτ 0 (2.102)
0
so that Z
dθ τ
+ θ = Ki (θs − θ(τ 0 )) dτ 0 (2.103)
dτ 0
Differentiating with respect to τ , we get
d2 θ dθ
+ + Ki θ = Ki θs (2.104)
dτ 2 dτ
The complementary function is emτ , where

m2 + m + Ki = 0 (2.105)

so that q  
1
m= −1 ± 1 − 4Ki (2.106)
2
The complete solution is
θ = C1 em1 τ + C2 em2 τ + θs (2.107)
where m1 and m2 are the roots from equation (2.106).
There is no offset in integral control. There are damped oscillations for Ki > 1/4,
and no oscillations if 0 < Ki < 1/4.

Integral-derivative control
Derivative control is usually not used alone since any constant θ would trigger no
response from the controller. Combining with integral control, we have
Z τ dθ
q = Ki (θs − θ(τ 0 )) dτ 0 − Kd (2.108)
0 dτ
from which Z
dθ τ dθ
+ θ = Ki (θs − θ(τ 0 )) dτ 0 − Kd (2.109)
dτ 0 dτ
Differentiating with respect to τ , we get
d2 θ 1 dθ Ki Ki
2
+ + θ= θs (2.110)
dτ 1 + Kd dτ 1 + Kd 1 + Kd
This can show overdamped or underdamped oscillatory behavior, depending on the
values of the constants Ki and Kd .

48
PID
If all three control actions are included, we have
Z τ dθ
q = Kp (θs − θ) + Ki (θs − θ(τ 0 )) dτ 0 + Kd (2.111)
0 dτ
from which
Z τ
dθ dθ
+ θ = Kp (θs − θ) + Ki (θs − θ(τ 0 )) dτ 0 + Kd (2.112)
dτ 0 dτ
This can be solved as in the above examples.

Variable ambient temperature

Let the ambient temperature vary sinusoidally as in
T∞ = T ∞ + δT sin ωt (2.113)
The governing equation (2.88) can be nondimensionalized using equations (2.9), (2.10)
and (2.90), where we use T ∞ instead of T∞ . The nondimensional equation is
dθ
+ θ = δT sin Ωτ + q(τ ) (2.114)
dτ
where
δT
δθ = (2.115)
Ti − T ∞
ωMc
Ω = (2.116)
hA
Using PID control, we take the heat source to be given by equation (2.111).
Equation (2.114) becomes
Z τ
dθ dθ
+ θ = δθ sin Ωτ + Kp (θs − θ) + Ki (θs − θ(τ 0 )) dτ 0 − Kd (2.117)
dτ 0 dτ
Differentiating
d2 θ 1 + Kp dθ Ki δθ Ω Ki
2
+ + θ= cos Ωτ + θs (2.118)
dτ 1 + Kd dτ 1 + Kd 1 + Kd 1 + Kd
The complementary functions are em1 τ and em2 τ , where m1 and m2 are obtained from
 v 
u 2
1 1 + Kp u 1 + Kp 4Ki 
m= − ±t −  (2.119)
2 1 + Kd 1 + Kd 1 + Kd

49
Appropriate values of Kp , Ki and Kd will produce m1 and m2 with negative real parts
for which the transient disappears with time and the control system is stable. The
particular solution is
θ = θs + a cos(Ωτ + φ) (2.120)
where
δθ Ω
a = − (2.121)
(1 + Kp )Ω sin φ + [(1 + Kd )Ω2 − Ki ] cos φ
(1 + Kp )Ω
tan φ = (2.122)
(1 + Kp )Ω2 − Ki
For optimal control the values of the parameters can be found to satisfy certain a
priori criteria for optimality.

Delay
Let us apply proportional control with a time delay of δτ , so that

q(τ ) = Kp (θs − θ(τ − δτ )) (2.123)

The governing equation is

dθ
+ θ + Kp θ(τ − δτ ) = Kp θs (2.124)
dτ
with the solution
Kp θs
θ = Cemτ + (2.125)
1 + Kp
where m satisfies the transcendental equation

m + 1 + Kp e−mδτ = 0 (2.126)

For small δτ , we find the approximate value

1 + Kp
m=− (2.127)
1 − Kp δτ
The control system is unstable if m > 0. If stable, the long time solution is
Kp
θ= θs (2.128)
1 + Kp
which is the same as equation (2.98). In Fig. 2.7 the region of instability is where
m > 0.

50
 Kp

m>0
Unstable

m<0
Stable

δτ

Figure 2.7: Convective cooling.

2.6.2 On-off control

This is a common form of thermal control in which the heating or cooling is turned
off at a predetermined temperature and turned on at another. Equation (2.88) for
this case is (
dT Q0 on
Mc + hA(T − T∞ ) = (2.129)
dt 0 off
where for the moment we take the ambient temperature, T∞ to be constant. The heat
rate is Q = Q0 when the system is on and Q = 0 when it is off. With the system in its
off mode, as t → ∞, T → Tmin = T∞ and in its on mode T → Tmax = T∞ + Q0 /hA.
Note that with the nomenclature Tmin and Tmax , we have implicitly assumed that we
are dealing with a heater, though the analysis is also applicable to cooling; in this
case Q0 would be negative and Tmin would be higher than Tmax .
Taking the nondimensional temperature to be

T − Tmin
θ = (2.130)
Tmax − Tmin
(2.131)

and the time as in equation (2.10), the governing equation is

(
dθ 1 on
+θ = (2.132)
dτ 0 off

51
The solution is (
1 + C1 e−τ on
θ= (2.133)
C2 e−τ off
We will assume that the heat source comes on when temperature falls below
a value TL , and goes off when it is above TU . These lower and upper bounds are
nondimensionally
TL − Tmin
θL = (2.134)
Tmax − Tmin
TU − Tmin
θU = (2.135)
Tmax − Tmin
After the initial transients have died disappeared, the oscillatory temperature
looks like that in Fig. 2.8. The heat source is on for a time interval τon during which
time the system temperature goes from θL to θU . The heater source is then switched
off and the system goes from θU back to θL in time τoff . These temperature conditions
can be applied to the solution, equation (2.133).

θ
θ=1

θU

θL
τ τ
on off

Figure 2.8: Response to on-off control.

During the on interval

θL = 1 + C1 (2.136)
θU = 1 + C1 e−τon (2.137)

52
which gives
1 − θL
τon = ln (2.138)
1 − θU
Similarly, in the off interval

θL = C2 (2.139)
θU = C2 e−τoff (2.140)

from which
θU
τoff = ln (2.141)
θL
The total period of the oscillation is

τp = τon + τoff (2.142)

θU (1 − θL )
= ln (2.143)
θL (1 − θU )

If we make a small dead-band assumption, we can write

θL = θs − δτ (2.144)
θU = θs + δτ (2.145)

where δ  1. A Taylor-series expansion gives

1 + δτ /(1 − θs )
τon = ln (2.146)
1 − δτ /(1 − θs )
2δτ
= + ... (2.147)
1 − θs

and
1 + δτ /θs
τoff = ln (2.148)
1 − δτ /θ
2δτ
= + ... (2.149)
θs

so that  
1 1
τp = 2δτ + + ... (2.150)
θs 1 − θs
The period is proportional to the width of the dead band.

53
 Alternatively, the small dead-band period can be calculated by assuming a saw-
toothed (piecewise linear) shape of the temperature-time curve. In equation (2.132),
if we take θ = θs , the equation approximates to
2δ
+ θs = 1 (2.151)
τon
2δ
− + θs = 0 (2.152)
τoff

for the on and off periods, respectively. From this, we get

τp = τon + τoff (2.153)

 
1 1
= 2δ + (2.154)
1 − θs θs
The complete problem can also be looked at as one with variables coefficients.
Writing (
0 τ /τon on
τ = (2.155)
1 + (τ − τon )/τoff off
equation (2.132) becomes
( ) (
dθ τon τon 0 ≤ τ 0 < 1
+ θ= (2.156)
dτ 0 τoff 0 1 ≤ τ0 < 2

The solution is ( 0
1 − (1 − θL )e−τon τ 0 ≤ τ0 < 1
θ= 0 (2.157)
θU e−τoff (τ −1) 1 ≤ τ0 < 2

Variable ambient temperature

Delay
It is possible in real systems that the on and off operations are delayed by a certain
time interval δτ . This is shown in Fig. 2.9 where the temperature is seen to vary
between the limits θ1 and θ2 instead of between θL and θU . The frequency and
amplitude of the oscillation is thus affected by the delay.
To determine the time period of the oscillation, let us begin with the last part of
the off period when the system goes from temperature θL to θ1 in time δτ . From the
exponential solution for this time interval, equation (2.133), we get

θL = C2 (2.158)
θ1 = C2 e−δτ (2.159)

54
 θ

θ2

θU
δτ δτ
2 τ
off

τ
on
θL
θ1
δτ
δτ1

Figure 2.9: On-off control with time delay δτ .

from which
θ1 = θL e−δτ (2.160)
When the system is on, the change in temperature from θ1 to θL occurs in time δ1 ,
so that

θ1 = 1 + C1 (2.161)
θL = 1 + C1 e−δτ1 (2.162)

From this
1 − θ1
δ1 = ln (2.163)
1 − θL
1 − θL e−δτ
= ln (2.164)
1 − θL
Again with the system on, the temperature goes from θL to θU in time τon that is
given by equation (2.138). Next, θ goes from θU to θ2 in time δτ . This gives

θU = 1 + C1 (2.165)
θ2 = 1 + C1 e−δτ (2.166)

55
from which
θ2 = 1 − (1 − θU )e−δτ (2.167)
The next time interval to consider is that when θ changes from θ2 to θU in time δτ2 .
This gives

θ2 = C2 (2.168)
θU = C2 e−δτ2 (2.169)

so that
θ2
δτ2 = ln (2.170)
θU
1 − (1 − θU )e−δτ
= ln (2.171)
θU
Finally, θ changes from θU to θL with the system off. This happens in time τoff that
is given by equation (2.141).
Summing up all these time intervals gives us the total period

τpδ τ = τon + τoff + 2δτ + δτ1 + δτ2 (2.172)

" #
θU (1 − θL ) 1 − θL e−δτ 1 − (1 − θU )e−δτ
= ln + 2δτ (2.173)
θL (1 − θU 1 − θL θU

which is larger due to the delay.

On the other hand the temperature excursion is given by

δθd = θ2 − θ1 (2.174)
= 1 − (1 − θU + θL )e−δτ (2.175)

which has also increased.

For a small delay, i.e. with δτ  1, we get the approximations
" #
θL 1 − θU
τPd = τp + δτ 2 + + + ... (2.176)
1 − θL θU
δθd = (θU − θL ) + δτ (1 − θU + θL ) + . . . (2.177)

2.7 Two-fluid problem

Suppose there is a body in contact with two fluids at different temperatures T∞,1 and
T∞,2 , like in the two examples shown in Fig. 2.10. The governing equation is

56
 T∞,2

T
T∞,2 T T∞,1
T∞,1

(a) (b)

Figure 2.10: Two-fluid problemss.

dT
Mc + h1 A1 (T − T∞,1 ) + h2 A2 (T − T∞,2 ) = 0 (2.178)
dt
where T (0) = Ti . If T∞,1 and T∞,2 are constants, we can nondimensionalize the
equation using the parameters for one of them, fluid 1 for instance. Thus we have
T − T∞,1
θ = (2.179)
Ti − T∞,1
h1 A1 t
τ = (2.180)
Mc
from which
dθ
+ θ + α(θ + β) = 0 (2.181)
dτ
with θ(0) = 1, where
h2 A2
α = (2.182)
h1 A1
T∞,1 − T∞,2
β = (2.183)
Ti − T∞,1
The equation can be written as
dθ
+ (1 + α)θ = −αβ (2.184)
dτ
with the solution
αβ
θ = Ce−(1+α)τ − (2.185)
1+α

57
The condition θ(0) = 1 gives C = 1 + αβ/(1 + α), from which
!
αβ αβ
θ = 1+ e−(1+α)τ − (2.186)
1+α 1+α

For α = 0, the solution reduces to the single-fluid case, equation (2.12). Otherwise
the time constant of the general system is
Mc
t0 = (2.187)
h1 A1 + h2 A2

2.8 Two-body problem

Suppose now that there are two bodies at temperatures T1 and T2 in thermal contact
with each other and exchanging heat with a single fluid at temperature T∞ as shown
in Fig. 2.11.

1 2

Figure 2.11: Two bodies in thermal contact.

The mathematical model of the thermal process is

dT1 ks Ac
M1 c1 + (T1 − T2 ) + hA(T1 − T∞ ) = 0 (2.188)
dt L
dT2 ks Ac
M2 c2 + (T2 − T1 ) + hA(T2 − T∞ ) = 0 (2.189)
dt L

Problems
1. Show that the temperature distribution in a sphere subject to convective cooling tends to
become uniform as Bi → 0.
2. Check one of the perturbation solutions against a numerical solution.
3. Plot all real θ(β, ) surfaces for the convection with radiation problem, and comment on the
existence of solutions.
4. Complete the problem of radiation in an enclosure (linear stability, numerical solutions).

58
5. Lumped system with convective-radiative cooling with nonzero θ0 and θs .
6. Find the steady-state temperatures for the two-body problem and explore the stability of the
system for constant ambient temperature.
7. Consider the change in temperature of a lumped system with convective heat transfer where
the ambient temperature, T∞ (t), varies with time in the form shown. Find (a) the long-time
solution of the system temperature, T (t), and (b) the amplitude of oscillation of the system
temperature, T (t), for a small period δt.

T∞
δt

Tmax

Tmin

t
Figure 2.12: Ambient temperature variation.

8. Study the PID control system for the two-body problem, each of which has a heat source that
can be independently controlled. Show numerical results for the different types of responses
possible (damped, oscillatory, unstable, stable, etc.). Take the ambient temperature, T∞ , to
be (a) constant and (b) oscillatory.
9. Consider on-off control for the two-body problem. Show analytical or numerical results for
the temperature responses of the two bodies. If you do the problem analytically, take the
ambient temperature, T∞ , to be constant, but if you do it numerically, then you can take it
to be (a) constant, and (b) oscillatory.

59
60
Chapter 3

One spatial dimensional

3.1 Fin theory

3.1.1 Definitions
Fin effectiveness f : This is the ratio of the fin heat transfer rate to the rate that
would be if the fin were not there.
Fin efficiency ηf : This is the ratio of the fin heat transfer rate to the rate that would
be if the entire fin were at the base temperature.

3.1.2 One-dimensional approximation

Longitudinal heat flux
Tb − T∞
qx00 = O(ks ) (3.1)
L
Transverse heat flux
qt00 = O(h(Tb − T∞ )) (3.2)

The transverse heat flux can be neglected compared to the longitudinal if

qx00  qt00 (3.3)

which gives a condition on the Biot number

hL
Bi = 1 (3.4)
k

61
 T
∞

T
b
x

Figure 3.1: Schematic of a fin.

3.1.3 Fin equation

Consider the fin shown shown in Fig. 3.1. The energy flows are indicated in Fig. 3.2.
The conductive heat flow along the fin, the convective heat loss from the side, and
the radiative loss from the side are
dT
qk = −ks A (3.5)
dx
qh = hdAs (T − T∞ ) (3.6)
qr = σdAs (T 4 − T∞
4
) (3.7)

where P (x) = dAs /dx is the perimeter. Heat balance gives

∂T ∂qk
ρAc + dx + qh + qr = 00 (3.8)
∂t ∂x
from which
∂T ∂ ∂T
ρAc − ks (A ) + P h(T − T∞ ) + σP (T 4 − T∞
4
)=0 (3.9)
∂t ∂x ∂x
where ks is taken to be a constant.

62
 q (x)+q (x)
h r

T
∞

q (x) q (x+dx)
k k

Figure 3.2: Energy balance.

Boundary conditions
The initial temperature is T (x, 0) = Ti (x). Usually the base temperature Tb is known.
The different types of boundary conditions for the tip are:

• Convective: ∂T /∂x = a at x = L

• Adiabatic: ∂T /∂x = 0 at x = L

• Known tip temperature: T = TL at x = L

• Long fin: T = T∞ as x → ∞

Nondimensionalization
Taking
T − T∞
θ = (3.10)
Tb − T∞
ks t
τ = Fourier modulus (3.11)
L2 ρc
x
ξ = (3.12)
L
63
 A
a(ξ) = (3.13)
Ab
P
p(ξ) = (3.14)
Pb
where the subscript indicates quantities at the base, the fin equation becomes
!
∂θ ∂ ∂θ h i
a − a + m2 pθ + p (θ + β)4 − β 4 = 0 (3.15)
∂τ ∂ξ ∂ξ
where
Pb hL2
m2 = (3.16)
ks Ab
σPb L2 (Tb − T∞ )3
 = (3.17)
ks Ab
T∞
β = (3.18)
Tb − T∞

3.2 Long time solution

The general fin equation is
!
∂θ ∂ ∂θ
a − a + f (θ) = 0 (3.19)
∂τ ∂ξ ∂ξ
where f (θ) includes heat transfer from the sides due to convection and radiation. The
boundary conditions are either Dirichlet or Neumannn type at ξ = 0 and ξ = 1. The
steady state is determined from
!
d dθ
− a + f (θ) = 0 (3.20)
dξ dξ

with the same boundary conditions. Substituting θ = θ + θ0 in equation (??) and

subtracting (3.20), we have
!
∂θ0 ∂ ∂θ0 h i
a − a + f (θ + θ0 ) − f (θ) = 0 (3.21)
∂τ ∂ξ ∂ξ
where θ0 is the perturbation from the steady state. The boundary conditions for θ0
are homogeneous. Multiplying by θ0 and integrating from ξ = 0 to ξ = 1, we have
dE
= I1 + I2 (3.22)
dτ
64
where
Z 1
1
E = a(θ0 )2 dξ (3.23)
2 0
!
Z 1 ∂ ∂θ0
I1 = θ0 a dξ (3.24)
0 ∂ξ ∂ξ
Z 1 h i
I2 = − θ0 f (θ + θ0 ) − f (θ) dξ (3.25)
0

Integrating by parts we can show that

1 Z !2
∂θ0 0
1 dθ0
I1 = θa − a dξ (3.26)
∂ξ 0 0 dξ
Z !2
1 dθ0
= − a dξ (3.27)
0 dξ

since the first term on the right side of equation (3.26) is zero due to boundary
conditions. Thus we know from the above that I1 is nonpositive and from equation
(3.23) that E is nonnegative. If we also assume that

I2 ≤ 0 (3.28)

then equation (3.22) tells us that E must decrease with time until reaching zero.
Thus the steady state is globally stable. Condition (3.28) holds if [θ0 and f (θ + θ0 ) −
f (θ)] are of the same sign or both zero; this is a consequence of the Second Law of
Thermodynamics.

3.3 Steady state solutions

Equation (3.19) reduces to
!
d dθ h i
− a + m2 pθ + p (θ + β)4 − β 4 = 0 (3.29)
dξ dξ

3.3.1 Uniform cross section

For this case a = p = 1, so that

d2 θ h i
2 4 4
− + m θ +  (θ + β) − β =0 (3.30)
dξ 2

65
Convective
With only convective heat transfer, we have

d2 θ
− + m2 θ = 0 (3.31)
dξ 2
the solution to whiich is

θ = C1 sinh mξ + C2 cosh mξ (3.32)

the constants are determined from the boundary conditions. For example, if

θ(0) = 1 (3.33)
dθ
(1) = 0 (3.34)
dξ
we get
θ = − tanh m sinh mξ + cosh mξ (3.35)

Radiative
The fin equation is
d2 θ h i
4 4
− +  (θ + β) − β =0 (3.36)
dξ 2
Let
φ=θ+β (3.37)
so that
d2 φ
− 2 + φ4 = −β 4 (3.38)
dξ
As an example, we will find a perturbation solution with the boundary conditions

φ(0) = 1 + β (3.39)
dφ
(1) = 0 (3.40)
dξ
We write
φ = φ0 + φ1 + 2 φ2 + . . . (3.41)
The lowest order equation is

dφ0 dφ0
= 0, φ0 (0) = 1 + β, (1) = 0 (3.42)
dξ 2 dξ

66
which gives
φ0 = 1 + β (3.43)
To the next order

dφ1 dφ1
2
= φ40 − β 4 , φ1 (0) = 0, (1) = 0 (3.44)
dξ dξ
with the solution
ξ2
φ1 = [(1 + β)4 − β 4 ] − [(1 + β)4 − β 4 ]ξ (3.45)
2
The complete solution is
( )
ξ2
φ = (1 + β) +  [(1 + β) − β ] − [(1 + β)4 − β 4 ]ξ + . . .
4 4
(3.46)
2

so that ( )
ξ2
θ = 1 +  [(1 + β)4 − β 4 ] − [(1 + β)4 − β 4 ]ξ + . . . (3.47)
2

Convective and radiative

(3.48)

3.3.2 Shape optimization

Consider a rectangular fin of length L and thickness δ a shown in Fig. 3.3. The
dimensional equation is
d2 T
− m2 (T − T∞ ) = 0 (3.49)
dx2
where !1/2
2h
m= (3.50)
ks δ
We will take the boundary conditions

T (0) = Tb (3.51)
dT
(L) = 0 (3.52)
dx
The solution is

T = T∞ − (T − T∞ ) [tanh mL sinh mx − cosh mx] (3.53)

67
 T∞
T
b

Figure 3.3: Rectangular fin.

The heat rate through the base per unit width is

dT
q = −ks δ (3.54)
dx x=0
= ks δ(Tb − T∞ )m tanh mL (3.55)
Writing L = Ap /δ, we get
!1/2  !1/2 
2h Ap 2h
q = ks δ (Tb − T∞ ) tanh   (3.56)
ks δ δ ks δ

Keeping Ap constant, i.e. constant fin volume, the heat rate can be maximized if
" !#  !1/2 
1/2 Ap 2h 2h 1/2 3 −5/2 1 −1/2 Ap 2h
δopt sech2 Ap (− )δopt + δopt tanh  =0
δ ks δopt ks 2 2 δ ks δopt
(3.57)
This is equivalent to
3βopt sech2 βopt = tanh βopt (3.58)
where !
Ap 2h
βopt = (3.59)
δopt ks δopt
Numerically, we find that βopt = 1.4192. Thus
 !1/2 2/3
Ap ks Ap
δopt =   (3.60)
βopt 2h

68
  !1/2 2/3
ks Ap
Lopt = βopt  (3.61)
2h

3.3.3 Annular fin

3.4 Two-dimensional fin analysis

3.4.1 Eccentric annulus

3.5 Transient conduction

Let us propose a similarity solution of the transient conduction equation
∂2T 1 ∂T
2
− =0 (3.62)
∂x κ ∂t
as !
x
T = A erf √ (3.63)
2 κt
Taking derivatives we find
!
∂T 1 x2
= A√ exp − (3.64)
∂x πκt 4κt
!
∂2T x x2
= −A √ exp − (3.65)
∂x2 2 πκ3 t3 4κt
!
∂T x x2
= −A √ exp − (3.66)
∂t 2 πκt3 4κt
so that substitution verifies that equation (3.63) is a solution to equation (3.62).

Problems
1. Consider a rectangular fin with convection, radiation and Dirichlet boundary conditions.
Calculate numerically the evolution of an initial temperature distribution at different instants
of time. Graph the results for several values of the parameters.
2. Consider a longitudinal fin of concave parabolic profile as shown in the figure, where δ =
[1 − (x/L)]2 δb . δb is the thickness of the fin at the base. Assume that the base temperature is
known. Neglect convection from the thin sides. Find (a) the temperature distribution in the
fin, and (b) the heat flow at the base of the fin. Optimize the fin assuming the fin volume to
be constant and maximizing the heat rate at the base. Find (c) the optimum base thickness
δb , and (d) the optimum fin height L.

69
 δ(x)
w
δb x

Figure 3.4: Longitudinal fin of concave parabolic profile.

70
Chapter 4

Multiple spatial dimensions

4.1 Steady-state problems

4.2 Transient problems
4.3 Stefan problems
Problems
1. This is a problem

71
72
Chapter 5

Phase change

5.1 Stefan problems

The two phases, indicated by subscripts 1 and 2, are separated by an interface at
x = X(t). In each phase, the conduction equations is
∂ 2 T1 1 ∂T1
2
− = 0 (5.1)
∂x κ1 ∂t
∂ 2 T2 1 ∂T2
2
− = 0 (5.2)
∂x κ2 ∂t
At the interface the temperature should be continuous, so that
T1 (X, t) = T2 (X, t) (5.3)
Furthermore the difference in heat rate into the interface provides the energy required
for phase change. Thus
∂T1 ∂T2 dX
k1 − k2 = Lρ (5.4)
∂x ∂x dt

5.1.1 Neumann’s solution

The material is initially liquid at T = T0 . The temperature at the x = 0 end is
reduced to zero for t > 0. Thus
T1 = 0 at x = 0 (5.5)
T2 → T0 as x → ∞ (5.6)
Assume T1 (x, t) to be
x
T1 = A erf √ (5.7)
2 κ1 t

73
so that it satisfies equations (5.1) and (5.5). Similarly
x x
T1 = A erf √ T2 = T0 − B erf √ (5.8)
2 κ1 t 2 κ2 t

satisfies equation (5.2) and (5.6). The, condition (5.3) requires that
x x
A erf √ = T0 − B erf √ = T1 (5.9)
2 κ1 t 2 κ2 t

This shows that √

X = 2λ κ1 t (5.10)
where λ is a constant. Using the remaining condition (5.4), we get
s
−λ2 κ1 −κ1 λ2 /κ2 √
k1 Ae − k2 B e = λLκ1 ρ π (5.11)
κ2

This can be written as

q
2 √
e−λ2 k2 κ21 (T0 − T1 )e−κ2 λ /κ2 λL π
− √ q = (5.12)
erf λ k1 κ2 T1 erfc(λ κ1 /κ2 ) c1 T1 s

The temperatures are

T1 x
T1 = erf ( √ ) (5.13)
erf λ 2 κ1 t
T0 − T1 x
T2 = T0 − q erfc( √ ) (5.14)
erfc(λ κ1 /κ2 ) 2 κ2 t

5.1.2 Goodman’s integral

74
 Part III

Convection

75
Chapter 6

One-dimensional forced convection

In this chapter we will considering the heat transfer in pipe flows. We will take
a one-dimensional approach and neglect transverse variations in the velocity and
temperature. In addition, for simplicity, we will assume that fluid properties are
constant and that the area of the pipe is also constant.

6.1 Hydrodynamics
6.1.1 Mass conservation
For a duct of constant cross-sectional area and a fluid of constant density, the mean
velocity of the fluid, u, is also constant.

6.1.2 Momentum equation

The forces on an element of length ds, shown in Fig. 6.1, in the positive s direction
are: fv , the viscous force and fp , the pressure force. We can write
fv = −τw P ds (6.1)
∂p
fp = −A ds (6.2)
∂s
where τw is the wall shear stress, and p is the pressure in the fluid. Since the mass of
the element is ρA ds, we can write the momentum equation as
du
ρA ds = fv + fp (6.3)
dt
from which we get
du τw P 1 ∂p
+ u=− (6.4)
dt ρA ρ ∂s

77
 p p+dp

s fp

f
v

ds

Figure 6.1: Forces on an element of fluid.

Integrating over the length L of a pipe, we have

du 4τw p1 − p2
+ u= (6.5)
dt ρDh ρL
where p1 and p2 are the pressures at the inlet and outlet respectively, and the hydraulic
diameter is defined by Dh = 4A/P .
For fully developed flow we can assume that τw is a function of u that depends
on the mean velocity profile, so that we can write
du
+ T (u)u = β ∆p (6.6)
dt
where T (u) is always positive, β = 1/ρL, and ∆p is the pressure difference that is
driving the flow. The wall shear stress is estimated below for laminar and turbulent
flows.

Laminar
The fully developed laminar velocity profile in a circular duct is given by the Poiseuille
flow result !
4r 2
u(r) = um 1 − 2 (6.7)
D
where u is the local velocity, r is the radial corrdinate, um is the maximum velocity
at the centerline, and D is the diameter of the duct. The mean velocity is given by
4 Z D/2
U= u(r) 2πr dr (6.8)
πD 2 0
Substituting the velocity profile, we get
um
U= (6.9)
2
78
The shear stress at the wall τw is given by

∂u
τw = µ (6.10)
∂r r=D/2
4
= −µum (6.11)
D
8µU
= − (6.12)
D
The wall shear stress is linear relationship

τw = αU (6.13)

where
8µ
α= (6.14)
D

Turbulent
For turbulent flow the expression for shear stress at the wall of a duct that is usually
used is  
f 1 2
τw = ρU (6.15)
4 2
Here f is the Darcy-Weisbach friction factor1 . The friction factor may be calculated
from the Blasius equation for smooth pipes

0.3164
f= (6.16)
Re1/4

where the Reynolds number is Re = UD/ν, or the Colebrook equation for rough
pipes
!
1 e/Dh 2.51
= −2.0 log + (6.17)
f 1/2 3.7 Re f 1/2
where e is the roughness at the wall, or similar expressions.
In the flow in a length of duct, L, without acceleration, the pressure drop is given
by
∆p A = τw P L (6.18)
1
Sometimes, confusingly, the Fanning friction factor, which is one-fourth the Darcy-Weisbach
value, is used in the literature.

79
where A is the cross-sectional area, and P is the inner perimeter. Thus
!  
L 1 2
∆p = f ρU (6.19)
4A/P 2
   
L 1 2
= f ρU (6.20)
Dh 2

Example 6.1
Consider a long, thin pipe with pressures p1 and p2 ate either end. For t ≤ 0, p1 − p2 = 0
and there is no flow. For t > 0, p1 −p2 is a nonzero constant. Find the resulting time-dependent
flow. Make the assumption that the axial velocity is only a function of radial position and time.

6.1.3 Long time behavior

Consider the flow in a single duct of finite length with a constant driving pressure
drop. The governing equation for the flow velocity is equation (6.6). The flow velocity
in the steady state is a solution of
T (u)u = β ∆p (6.21)

where ∆p and u are both of the same sign, say nonnegative. We can show that under
certain condtions the steady state is globally stable. Writing u = u + u0 , equation
(6.6) becomes
du0
+ T (u + u0 )(u + u0 ) = β ∆p (6.22)
dt
Subtracting equation (6.21), we get
du0
= −T (u + u0 )(u + u0 ) + T (u)(u) (6.23)
dt
Defining
1
E = u02 (6.24)
2
so that E ≥ 0, we find that
dE du0
= u0 (6.25)
dt dt
= −u0 [T (u + u0 )(u + u0 ) − T (u)u] (6.26)
= −u0 u [T (u + u0 ) − T (u)] − u02 T (u + u0 ) (6.27)

80
If we assume that T (u) is a non-decreasing function of |u|, we see that
u0 u [T (u + u0) − T (u)] ≥ 0 (6.28)
regardless of the sign of either u0 or u, so that
dE
≤0 (6.29)
dt
Thus, E(u) is a Lyapunov function, and u = u is globally stable to all perturbations.

6.2 Energy equation

Consider a section of a duct shown in Fig. 6.2, where an elemental control volume is
shown. The heat rate going in is given by
∂T
Q− = ρAucT − kA (6.30)
∂s
where the first term on the right is due to the advective and second the conductive
transports. c is the specific heat at constant pressure and k is the coefficient of
thermal conductivity. The heat rate going out is
+ ∂Q−
−
Q =Q + ds (6.31)
∂s
The difference between the two is
∂Q−
Q+ − Q− = ds
" ∂s #
∂T ∂2T
= ρAuc − kA 2 ds (6.32)
∂s ∂s
Furthermore, heat is gained from the side at a rate Q, which can be written as
Q = q ds (6.33)
where q is the rate of gain of heat per unit length of the duct.
An energy balance for the elemental control volume gives
∂T
Q− + Q = Q+ + ρA ds c (6.34)
∂t
where the last term is the rate of accumulation of energy within the control volume.
Substituting equations (6.32) and (6.33) in (6.34) we get the energy equation
∂T ∂T q k ∂2T
+u = + (6.35)
∂t ∂s ρAc ρc ∂s2
The two different types of heating conditions to consider are:

81
 Q

_ +
Q s Q

ds

Figure 6.2: Forces on an element of fluid.

6.2.1 Known heat rate

The heat rate per unit length, q(s), is known all along the duct. Defining
x
ξ = (6.36)
L
(T − Ti )ρV AC
θ = (6.37)
Lq
tV
τ = (6.38)
L
gives
∂θ ∂θ d2 θ
+ −λ 2 =1 (6.39)
∂τ ∂ξ dξ
where
k
λ= (6.40)
LV ρc
Boundary conditions may be θ = 0 at ξ = 0, θ = θ1 at ξ = 1.

6.2.2 Known wall temperature

The heating is now convective with a heat transfer coefficient h, and an external
temperature of Tw (s). Thus,
q = P h(Tw − T ) (6.41)
Defining
x
ξ = (6.42)
L
82
 T∞ (t)

Tin (t) u - T (s, t) Tout (t)

Figure 6.3: Fluid duct with heat loss.

T − Ti
θ = (6.43)
Tw − Ti
tV
τ = (6.44)
L
gives
∂θ ∂θ d2 θ
+ − λ 2 + Hθ = H (6.45)
∂τ ∂ξ dξ
where
k
λ = (6.46)
LV ρc
hρL
H = (6.47)
ρV Ac

6.3 Single duct

Consider the duct that is schematically shown in Fig. 6.3. The inlet temperature is
Tin (t), and the outlet temperature is Tout (t), and the fluid velocity is u. The duct
is subject to heat loss through its surface of the form UP (T − T∞ ) per unit length,
where the local fluid temperature is T (s, t) and the ambient temperature is T∞ (t). U
is the overall heat transfer coefficient and P the cross-sectional perimeter of the duct.
We assume that the flow is one-dimensional, and neglect axial conduction through
the fluid and the duct. Using the same variables to represent non-dimensional quan-
tities, the governing non-dimensional equation is
∂T ∂T  
+ + γ T − Te∞ = 0 (6.48)
∂s ∂t
where the nondimensional variables are
s
s = (6.49)
L
tu
t = (6.50)
L
T − T∞
T = (6.51)
∆T
83
The characteristic time is the time taken to traverse the length of the duct, i.e. the
residence time. The ambient temperature is
T∞ (t) = T ∞ + Te∞ (t) (6.52)
where the time-averaged and fluctuating parts have been separated. Notice that the
nondimensional mean ambient temperature is, by definition, zero. The characteristic
temperature difference ∆T will be chosen later. The parameter γ = UP L/ρAuc
represents the heat loss to the ambient.

6.3.1 Steady state

No axial conduction
The solution of the equation
dθ
+ H(θ − 1) = 0 (6.53)
dξ
with boundary condition θ(0) = 0 is
θ(ξ) = 1 − e−Hξ (6.54)

With small axial conduction

We have
d2 θ dθ
λ− − H(θ − 1) = 0 (6.55)
dξ 2 dξ
where λ  1, and with the boundary conditions θ(0) = 0 and θ(1) = θ1 .
We can use a boundary layer analysis for this singular perturbation problem.
The outer solution is
θout = 1 − e−Hξ (6.56)
The boundary layer is near ξ = 1, where we make the transformation
ξ−1
X= (6.57)
λ
This gives the equation
d2 θin dθin
2
− − λH(θin − 1) = 0 (6.58)
dX dX
To lowest order, we have
d2 θin dθin
λ 2
− =0 (6.59)
dX dX
84
with the solution
θin = A + BeX (6.60)
The boundary condition θin (X = 0) = θ1 gives θ1 = A + B, so that

θin = A + (θ1 − A)eX (6.61)

The matching conditions is

θouter (ξ = 1) = θin (X → −∞) (6.62)

so that
A = 1 − e−H (6.63)
The composite solution is then

θ = 1 − e−H + (θ1 − 1 + e−H )e(ξ−1)/λ + . . . (6.64)

6.3.2 General solution

The general solution of this equation is
 Z t 
γt0 e
T (s, t) = f (s − t) + γ e T∞ (t ) dt e−γt
0 0
(6.65)
0

The boundary conditions T (0, t) = Tin (t) and T (s, 0) = T0 (s) are shown in Fig. 6.4.
The solution becomes
( R
eγt Te∞ (t0 ) dt0
0
Tin (t − s)e−γs + γe−γt t−s
t
for t ≥ s
T (s, t) = R (6.66)
T0 (s − t)e−γt + γe−γt 0t eγt Te∞ (t0 ) dt0
0
for t < s

The t < s part of the solution is applicable to the brief, transient period of time in
which the fluid at time t = 0 has still not left the duct. The later t > s part depends
on the temperature of the fluid entering at s = 0. The temperature, Tout (t), at the
outlet section, s = 1, is given by
( R
eγt Te∞ (t0 ) dt0
0
Tin (t − 1)e−γ + γe−γt t−1
t
for t ≥ 1
Tout (t) = R (6.67)
T0 (1 − t)e−γt + γe−γt 0t eγt Te∞ (t0 ) dt0
0
for t < 1

It can be observed that, after an initial transient, the inlet and outlet tempera-
tures are related by a unit delay. The outlet temperature is also affected by the heat
loss parameter, γ, and the ambient temperature fluctuation, Te∞ . The following are
some special cases of equation (6.67).

85
 t
t=s
t>s
T (0, t) = Tin (t)
@R t<s
T (s, 0) = T0 (s)
s

Figure 6.4: Solution in s-t space.

6.3.3 Perfectly insulated duct

If γ = 0 the outlet temperature simplifies to
(
Tin (t − 1) for t ≥ 1
Tout (t) = (6.68)
T0 (1 − t) for t < 1
The outlet temperature is the same as the inlet temperature, but at a previous instant
in time.

6.3.4 Constant ambient temperature

For this Te∞ = 0, and equation (6.67) becomes
(
Tin (t − 1)e−γ for t ≥ 1
Tout (t) = (6.69)
T0 (1 − t)e−γt for t < 1
This is similar to the above, but with an exponential drop due to heat transfer.

6.3.5 Periodic inlet and ambient temperature

We take
Tin (t) = T in + Tbin sin ωt (6.70)
Te∞ (t) = Tb∞ sin Ωt (6.71)
so that equation (6.67) becomes
 h i

 T + b sin ω(t − 1) e−γ
T

 in in
q
+Tb∞ γ −2e γ 2cos
−γ 1+e−2γ
Tout (t) =  2 sin(Ωt + φ) for t ≥ 1 (6.72)

 √
+Ω
 T (1 − t)e
0
−γt
+ γ
Tb γ 2 + Ω2 sin(Ωt + φ0 )
γ 2 +Ω2 ∞ for t < 1

86
 Figure 6.5: Effect of wall.

where
γ(1 − e−γ cos 1) + e−γ Ω sin 1
tan φ = − (6.73)
Ω(1 − e−γ cos 1) − γe−γ sin 1
Ω
tan φ0 = − (6.74)
γ
The outlet temperature has frequencies which come from oscillations in the inlet as
well as the ambient temperatures. A properly-designed control system that senses the
outlet temperature must take the frequency dependence of its amplitude and phase
into account. There are several complexities that must be considered in practical
applications to heating or cooling networks, some of which are analyzed below.

6.3.6 Effect of wall

The governing equations are

∂T ∂T ∂2T
ρAc + ρV Ac − kA 2 + hi Pi (T − Tw ) = 0 (6.75)
∂t ∂x ∂x
2
∂Tw ∂ Tw
ρw Aw cw − kw Aw + hi Pi (Tw − T ) + ho Po (Tw − T∞ ) = 0 (6.76)
∂t ∂x2
(6.77)

Nondimensionalize, using
x
ξ = (6.78)
L
tV
τ = (6.79)
L
87
 T − T∞
θ = (6.80)
Ti − T∞
Tw − T∞
θw = (6.81)
Ti − T∞

we get

∂θ ∂θ ∂2θ
+ − λ 2 + Hin (θ − θw ) = 0 (6.82)
∂τ ∂ξ ∂ξ
2
∂θw ∂ θw
− λw 2 + Hin (θw − θ) + Hout θw = 0 (6.83)
∂τ ∂ξ

where

kw
λw = (6.84)
ρw V Aw cw L
hin Pin L
Hin = (6.85)
ρw Aw cw V
hout Pout L
Hout = (6.86)
ρw Aw cw V

In the steady state and with no axial conduction in the fluid

dθ
+ Hin (θ − θw ) = 0 (6.87)
dξ
d2 θw
−λw + Hin (θw − θ) + Hout θw = 0 (6.88)
dξ 2

If we assume λw = 0 also, we get

Hin
θw = θ (6.89)
Hin + Hout

The governing equation is

dT
+ Hw θ = 0 (6.90)
dξ
where
w w
Hin Hout
Hw = w w
(6.91)
Hin + Hout

88
 1

Figure 6.6: Two-fluids with wall.

6.4 Two-fluid configuration

Consider the heat balance in Fig. 6.6. Neglecting axial conduction, we have

∂Tw
ρw Aw cw + h1 (Tw − T1 ) + h2 (Tw − T2 ) = 0 (6.92)
∂t
∂T1 ∂T1
ρ1 A1 c1 + ρ1 V1 c1 + h1 (T1 − Tw ) = 0 (6.93)
∂t ∂x
∂T2 ∂T2
ρ2 A2 c2 + ρ2 V2 c2 + h2 (T2 − Tw ) = 0 (6.94)
∂t ∂x

6.5 Regenerator
A regenerator is schematically shown in Fig. 6.7.

dT
Mc + ṁc(Tin − Tout ) = 0 (6.95)
dt

89
 Figure 6.7: Schematic of regenerator.

6.6 Networks
A network consists of a number of ducts that are united at certain points. At each
junction, we must have X
Ai ui = 0 (6.96)
i

where Ai are the areas and ui the fluid velocities in the ducts coming in, the sum being
over all the ducts entering the junction. Furthermore, for each duct, the momentum
equation is
dui h i
+ T (ui)ui = β pin i − pi
out
+ ∆p (6.97)
dt
where ∆p is the pressure developed by a pump, if there happens to be one on that
line. We must distinguish between two possible geometries.

(a) Two-dimensional networks: A planar or two-dimensional network is one that is

topologically equivalent to one on a plane in which every intersection of pipes indicates
fluid mixing. For such a graph, we know that

E =V +F −1 (6.98)

where E, V and F are the number of edges, vertices and faces, respectively. In
the present context, these are better referred to as branches, junctions and circuits,
respectively.
The unknowns are the E velocities in the ducts and the V pressures at the
junctions, except for one pressure that must be known. The number of unknowns
thus are E + V − 1. The momentum equation in the branches produce E independent
differential equations, while mass conservation at the juntions give V − 1 independent
algebraic relations. Thus the number of

(b) Three-dimensional networks: For a three-dimensioanl network, we have

E =V +F −2 (6.99)

90
If there are n junctions, they can have a maximum of n(n − 1)/2 lines connecting
them. The number of circuits is then (n2 − 3n + 4)/2. The number of equations to
be solved is thus quite large if n is large.

6.6.1 Hydrodynamics
The global stability of flow in a network can be demonstrated in a manner similar to
that in a finite-length duct. In a general network, assume that there are n junctions,
and each is connected to all the rest. Also, pi is the pressure at junction i, and uij
is the flow velocity from junction i to j defined to be positive in that direction. The
flow velocity matrix uij is anti-symmetric, so that uii which has no physical meaning
is considered zero.
The momentum equation for uij is
duij
+ Tij (uij )uij = βij (pi − pj ) (6.100)
dt
The network properties are represented by the symmetric matrix βij . The resistance
Tij may or may not be symmetric. To simplify the analysis the network is considered
fully connected, but Tij is infinite for those junctions that are not physically connected
so that the flow velocity in the corresponding branch is zero. We take the diagonal
terms in Tij to be also infinite, so as to have uii = 0.
The mass conservation equation at junction j for all flows arriving there is
X
n
Aij uij = 0 for j = 1, . . . , n (6.101)
i=1

where Aij is a symmetric matrix. The symmetry of Aij and uij gives the equivalent
form n X
Ajiuji = 0 for j = 1, . . . , n (6.102)
i=1
which is simply the mass conservation considering all the flows leaving junction j.
The steady states are solutions of
duij
+ Tij (uij )uij = βij (pi − pj ) (6.103)
dt
X
n X
n
Aij uij = 0 or Aji uji = 0 (6.104)
i=1 i=1

We write
uij = uij + u0ij (6.105)
pi = p + p0i (6.106)

91
Substituting in equations (6.100)–(6.102), and subtracting equations (6.103) and
(6.104) we get

du0ij h i
= − Tij (uij + u0ij )(uij + u0ij ) + Tij (uij )uij + βij (p0i − p0j(6.107)
)
dt
X
n X
n
Aij u0ij = 0 or Ajiu0ji = 0 (6.108)
i=1 i=1

Defining
1X n X n
Aij 0 2
E= u (6.109)
2 j=1 i=1 βij ij
we get

dE Xn X n
Aij 0 du0ij
= uij (6.110)
dt j=1 i=1 βij dt
X
n X
n
Aij h i
= − u0ij Tij (uij + u0ij )(uij + u0ij ) − Tij (uij )uij
j=1 i=1 βij
Xn X n
+ Aij u0ij (p0i − p0j ) (6.111)
j=1 i=1

The pressure terms vanish since

X
n X
n X
n X
n
Aij u0ij p0i = Aij u0ij p0i (6.112)
j=1 i=1 i=1 j=1
 
X
n X
n
=  p0 Aij u0ij  (6.113)
i
i=1 j=1

= 0 (6.114)

and
!
X
n X
n X
n X
n
Aij u0ij p0j = p0j Aij u0ij (6.115)
j=1 i=1 j=1 i=1
= 0 (6.116)

The terms that are left in equation (6.111) are similar to those in equation (6.28) and
satisfy the same inequality. Since E ≥ 0 and dE/dt ≤ 0, the steady state is globally
stable. For this reason the steady state is also unique.

92
 p
2
u
20

u
10

p p
1 0
u
30

p
3
Figure 6.8: Star network.

Example 6.2
Show that the flow in the the star network shown in Fig. 6.8 is globally stable. The
pressures p1 , p2 and p3 are known while the pressure p0 and velocities u10 , u20 and u30 are the
unknowns.
For branches i = 1, 2, 3, equation (6.6) is

dui0
+ Ti0 (ui0 )ui0 = βi0 (pi − p0 ) (6.117)
dt
Equation (6.96) at the junction gives
3
X
Ai0 ui0 = 0 (6.118)
i=1

In the steady state

Ti0 (ui0 )ui0 = βi0 (pi − p0 ) (6.119)

3
X
Ai0 ui0 = 0 (6.120)
i=1

Substituting ui0 = ui0 + u0i0 and p0 = p0 + p00 in equations (6.117) and (6.118) and subtracting

93
 equations (6.119) and (6.120), we find that

du0i0
= − [Ti0 (ui0 + u0i0 )(ui0 + u0i0 ) − Ti0 (ui0 )ui0 ] − βi0 p00 (6.121)
dt
3
X
Ai0 u0i0 = 0 (6.122)
i=1

If we define
3
1 X Ai0 0 2
E= u (6.123)
2 i=1 βi0 i
we find that
3
X
dE Ai0 du0i0
= u0i0 (6.124)
dt i=1
βi0 dt
3
X 3
X
Ai0
= − u0i0 [Ti0 (ui0 + u0i0 )(ui0 + u0i0 ) − Ti0 (ui0 )ui0 ] − p00 Ai0 u0i0 (6.125)
i=1
βi0 i=1

The last term vanishes beacuse of equation (6.122). Thus

X3 X3
dE Ai0 0 Ai0 0 2
=− ui0 ui0 [Ti0 (ui0 + u0i0 ) − Ti0 (ui0 )] − u i0 Ti0 (ui0 + u0i0 ) (6.126)
dt i=1
β i0 i=1
β i0

Since E ≥ 0 and dE/dt ≤ 0, E is a Lyapunov function and the steady state is globally stable.

6.6.2 Thermal networks

6.7 Thermal control

6.7.1 Multiple room temperatures
Let there be n interconnected rooms. The wall temperature of room i is Tiw and the
air temperature is Tia . The heat balance equation for this room is

dTiw
Mia cw = hAi (Tia − Tiw ) + UAei (T e − Tiw ) (6.127)
dt
dT a
1 X
Mia ca i = hAi (Tiw − Tia ) + ca (maji + |maji|)Tja
dt 2 j
1 X
− ca (maij + |maij |)Tia + qi (6.128)
2 j

94
where T e is the exterior temperature, mij is the mass flow rate of air from room i to
room j. By definition mij = −mji . Since mii has no meaning and can be arbitrarily
taken to be zero, mij is an anti-symmetric matrix. Also, from mass conservation for
a single room, we know that
X
maji = 0 (6.129)
j

Analysis
The unknowns in equations (6.127) and (6.128) are the 2n temperatures Tiw and Tia .
(i) Steady state with U = 0
(a) The equality
XX XX
(maji + |maji|)Tja − (maij + |maij |)Tia = 0 (6.130)
i j i j

can be shown by interchanging i and j in the second term. Using this result, the sum
of equations (6.127) and (6.128) for all rooms gives
X
qi = 0 (6.131)
i

which is a necessary condition for a steady state.

(b) Because the sum of equations (6.127) and (6.128) for all rooms gives an identity,
the set of equations is not linearly independent. Thus the steady solution is not
unique unless one of the room temperatures is known.

Control
The various proportional control schemes possible are:

• Control of individual room heating

qi = −Ki (Tia − Tiset ) (6.132)

• Control of mass flow rates

maji = fij (Tja , Tia , Tiset ) (6.133)

Similar on-off control schemes can also be proposed.

95
6.7.2 Two rooms
Consider two interconnected rooms 1 and 2 with mass flow m from 1 to 2. Also there
is leakage of air into room 1 from the exterior at rate m, and leakage out of room 2
to the exterior at the same rate. The energy balances for the two rooms give
dT1 1
M1 ca = U1 A1 (T e − T1 ) + (m + |m|)(T e − T1 )
dt 2
1
− (m − |m|)(T2 − T1 ) + q1 (6.134)
2
dT 2 1
M2 ca = U2 A2 (T e − T2 ) − (m − |m|)(T e − T2 )
dt 2
1
+ (m + |m|)(T1 − T2 ) + q2 (6.135)
2
The overall mass balance can be given by the sum of the two equations to give
dT1 dT2
M1 ca + M2 ca = U1 A1 (T e − T1 ) + U2 A2 (T e − T2 ) + |m|T e
dt dt
1 1
+ (m − |m|)T1 − (m + |m|)T2
2 2
+q1 + q2 (6.136)

One example of a control problem would be to change m to keep the temperatures

of the two rooms equal. Delay can be introduced by writing T2 = T2 (t − τ ) and
T1 = T1 (t−τ ) in the second to last terms of equations (6.134) and (6.135), respectively,
where τ is the time taken for the fluid to get from one room to the other.

Problems
1. This is a problem

96
Chapter 7

One-dimensional natural
convection

7.1 Modeling
Let us consider a closed loop, shown in Fig. 7.1, of length L and constant cross-
sectional area A filled with a fluid. The loop is heated in some parts and cooled
in others. The temperature differences within the fluid leads to a chenge in density
and hence a buoyancy force that creates a natural circulation. The spatial coordi-
nate is s, measured from some arbitrary origin and going around the loop in the
counterclockwise direction.
We will make the Boussinesq approximation by which the fluid density is constant
except in the buoyancy term. We will also approximate the behavior of the fluid using
one spatial dimensions. Thus, we will assume that the velocity u and temperature T
are constant across a section of the loop. In general both u and T are functions of

Figure 7.1: A general natural convective loop.

97
 +
m

_
m
s
ds

Figure 7.2: Mass flows in an elemental control volume.

space s and time t, though we will find that u = u(t).

7.1.1 Mass conservation

Consider an elemental control volume as shown in Fig. 7.2. The mass fluxes in and
out are

m− = ρ0 uA (7.1)
−
∂m
m+ = m− + ds (7.2)
∂s
For a fluid of constant density, there is no accumulation of mass within an elemental
control volume, so that the mass flow rate into and out of the control volume must
be the same, i.e. m− = m+ . For a loop of constant cross-sectional area, this implies
that u is the same into and out of the control volume. Thus u is independent of s,
and must be a function of t alone.

7.1.2 Momentum equation

The forces on an element of length ds, shown in Fig. 7.3, in the positive s direction
are: fv , the viscous force, fp , the pressure force, and fg , the component of the gravity
force. We can write

fv = −τw P ds (7.3)

98
 ∂p
fp = −A ds (7.4)
∂s
fg = −ρA ds g̃ (7.5)

where τw is the wall shear stress, and p is the pressure in the fluid. It is impossible to
determine the viscous force fv through a one-dimensional model, since it is a velocity
profile in the tube that is responsible for the shear streass at the wall. For simplicity,
however, we will assume a linear relationship between the wall shear stress and the
mean fluid velocity, i.e. τw = αu. For Poiseuille flow in a duct, which is strictly not
the case here but gives an order of magnitude value for the coefficient, this would be
4µ
α= (7.6)
R
The local component of the acceleration due to gravity has been written in terms of

g̃(s) = g cos θ (7.7)

dz
= g (7.8)
ds
where g is the usual acceleration in the vertical direction, g̃ is its component in the
negative s direction, and dz is the difference in height at the two ends of the element,
with z being measured upwards. The integral around a closed loop should vanish, so
that Z L
g̃(s) ds = 0 (7.9)
0
The density in the gravity force term will be taken to decrease linearly with temper-
ature, so that
ρ = ρ0 [1 − β(T − T0 )] (7.10)
Since the mass of the element is ρ0 A ds, we can write the momentum equation as
du
ρ0 A ds = fv + fp + fg (7.11)
dt
from which we get
du Pα 1 ∂p
+ u=− − [1 − β(T − T0 )] g̃ (7.12)
dt ρ0 A ρ0 ∂s
Integrating around the loop, we find that the pressure term disappears, and
Z L
du Pα β
+ u= T g̃(s) ds (7.13)
dt ρ0 A L 0

where u = u(t) and T = T (s, t).

99
 fv
p+dp fg
fp
p
θ
s
ds

gravity

Figure 7.3: Forces on an element of fluid.

7.1.3 Energy equation

Fig. 7.4 shows the heat rates going into and out of an elemental control volume. The
heat rate going in is given by
∂T
Q− = ρ0 Aucp T − kA (7.14)
∂s
where the first term on the right is due to the advective and second the conductive
transports. cp is the specific heat at constant pressure and k is the coefficient of
thermal conductivity. The heat rate going out is

∂Q−
Q+ = Q− + ds (7.15)
∂s
The difference between the two is

+ − ∂Q−
Q −Q = ds
" ∂s #
∂T ∂2T
= ρ0 Aucp − kA 2 ds (7.16)
∂s ∂s

Furthermore, heat is gained from the side at a rate Q, which can be written as

Q = q ds (7.17)

where q is the rate of gain of heat per unit length of the duct.

100
 Q
+
Q

_
Q
s
ds

Figure 7.4: Heat rates on an elemental control volume.

An energy balance for the elemental control volume gives

∂T
Q− + Q = Q+ + ρ0 A ds cp (7.18)
∂t
where the last term is the rate of accumulation of energy within the control volume.
Substituting equations (7.16) and (7.17) in (7.18) we get the energy equation
∂T ∂T q k ∂2T
+u = + (7.19)
∂t ∂s ρ0 Acp ρ0 cp ∂s2

7.2 Known heat rate

The simplest heating condition is when the heat rate per unit length, q(s), is known
all along the loop. For zero mean heating, we have
Z L
q(s) ds = 0 (7.20)
0

q(s) > 0 indicates heating, and q(s) < 0 cooling.

7.2.1 Steady state, no axial conduction

Neglecting axial conduction, the steady-state governing equations are
Z L
Pα β
u = T (s)g̃(s) ds (7.21)
ρ0 A L 0

101
 dT q(s)
u = (7.22)
ds ρ0 Acp

The solution of equation (7.22) gives us the temperature field

Z
1 s
T (s) = q(s0 ) ds0 + T0 (7.23)
ρ0 Acp u 0

where T (0) = T0 . Using equation (7.9) it can be checked that T (L) = T0 also.
Substituting in equation (7.21), we get
Z L Z s 
Pα β 0 0
u= q(s ) ds g̃(s) ds (7.24)
ρ0 A ρ0 Acp Lu 0 0

from which v
u Z Z 
u β L s
u = ±t q(s0 ) ds0 g̃(s) ds (7.25)
P αLcp 0 0

Two real solutions exist for

Z L Z s 
0 0
q(s ) ds g̃(s) ds ≥ 0 (7.26)
0 0

and none otherwise. Thus there is a bifurcation from no solution to two as the
parameter H passes through zero, where
Z L Z s 
H= q(s0 ) ds0 g̃(s) ds (7.27)
0 0

The pressure distribution can be found from equation (7.12)

dp P αu h i
= − − ρ0 1 − β(T − T0 ) g̃ (7.28)
ds A Z s 
P αu β 0 0
= − − ρ0 g̃ + q(s ) ds g̃ (7.29)
A Acp u 0

from which
Z s " #
P αu 0 0 β Z s Z s00
p(s) = p0 − s − ρ0 g̃(s ) ds + q(s ) ds g̃(s00 ) ds00
0 0
(7.30)
A 0 Acp u 0 0

where p(0) = p0 . Using equations (7.9) and (7.24), it can be shown that p(L) = p0
also.

102
 u

Figure 7.5: Bifurcation with respect to parameter H.

Example 7.1
Find the temperature distributions and velocities in the three heating and cooling distri-
butions corresponding to Fig. 7.6. (a) Constant heating between points c and d, and constant
cooling between h and a. (b) Constant heating between points c and d, and constant cooling
between g and h. (c) Constant heating between points d and e, and constant cooling between
h and a. (d) Constant heating between points a and c, and constant cooling between e and g.
The constant value is q̂, and the total length of the loop is L.

Let us write
Z s
F (s) = q(s0 ) ds0 (7.31)
0
G(s) = F (s)g(s) (7.32)
Z L
H = G(s) ds (7.33)
0

The functions F (s) and G(s) are shown in Fig. 7.7. The origin is at point a, and the coordinate
s runs counterclockwise. The integral H in the four cases is: (a) H = 0, (b) H = q̂L/8, (c)
H = −q̂L/8, (d) H = q̂L/4.
The fluid velocity is s
βH
u=± (7.34)
P αLcp

103
 g f e

h d

a b c

Figure 7.6: Geometry of a square loop.

No real solution exists for case (c); the velocity is zero for (a); the other two cases have two
solutions each, one positive and the other negative. The temperature distribution is given by
F (s)
T − T0 = (7.35)
ρ0 Acp u
The function F (s) is shown in Fig. 7.7. There is no real; solution for case (c); for (a), the
temperature is unbounded since the fluid is not moving; for the other two cases there are two
temperature fields, one the negative of the other.
The pressure distribution can be found from equation (7.29).

Example 7.2
What is the physical interpretation of condition (7.26)?

Let us write
Z L Z s 
H = q(s0 ) ds0 g̃(s) ds (7.36)
0 0
Z L Z s  Z s 
= q(s0 ) ds0 d g̃(s0 ) ds0 (7.37)
0 0 0
Z s L Z s L Z L  Z s 
0 0 0 0 0 0
= q(s ) ds g̃(s ) ds − q(s) g̃(s ) ds ds (7.38)
0 0 0 0 0 0

The first term on the right vanishes due to equations (7.20) and (7.9). Using equation (7.8),
we find that Z L
H = −g q(s)z(s) ds (7.39)
0

104
F F

s s

G G

s s
(a) (b)

F F

s s

G G

s s
(c) (d)

Figure 7.7: Functions F (s) and G(s) for the four cases.

105
 The function z(s) is another way of describing the geometry of the loop. We introduce the
notation
q(s) = q + (s) − q − (s) (7.40)
where 
q(s) for q(s) > 0
q+ = (7.41)
0 for q(s) ≤ 0
and 
− 0 for q(s) ≥ 0
q = (7.42)
−q(s) for q(s) < 0
Equations (7.20) and (7.39) thus becomes
Z L Z L
q + (s) ds = q − (s) ds (7.43)
0 0
"Z Z #
L L
+ −
H = −g q (s)z(s) ds − q (s)z(s) ds (7.44)
0 0

From these, condition (7.26) which is H ≥ 0 can be found to be equivalent to

RL RL
0
q + (s)z(s) ds 0
q − (s)z(s) ds
RL < RL (7.45)
0
q + (s) ds 0
q − (s) ds

This implies that the height of the centroid of the heating rate distribution should be above
that of the cooling.

7.2.2 Axial conduction effects

To nondimensionalize and normalize equations (7.13) and (7.19), we take

t
t∗ = (7.46)
τ
s
s∗ = (7.47)
L
u
u∗ = (7.48)
V G1/2
T − T0
T∗ = (7.49)
∆T G1/2
g̃
g̃ ∗ = (7.50)
g
q
q∗ = (7.51)
qm

106
where
P αL
V = (7.52)
ρ0 A
P 2 α2 L
∆T = (7.53)
βgρ20 A2
ρ0 A
τ = (7.54)
Pα
qm βgρ20 A2
G = (7.55)
P 3 α3 Lcp
Substituting, we get
Z
du∗ ∗
1
+ u = T ∗ g̃ ∗ ds∗ (7.56)
dt∗ 0
∂T ∗ 1/2 ∗ ∂T ∗
1/2 ∗ ∂2T ∗
+ G u = G q + K (7.57)
∂t∗ ∂s∗ ∂s∗2
where
kA
K= (7.58)
P αL2 cp
The two nondimensional parameters which govern the problem are G and K.
Under steady-state conditions, and neglecting axial conduction, the temperature
and velocity are
∗ 1 Z s∗ ∗ ∗
T (s) = q (s1 ) ds∗1 (7.59)
u∗s 0
Z 1 Z s∗ 
∗
u = ± q ∗ (s∗1 ) ds∗1 g̃ ∗ (s∗ ) ds∗ (7.60)
0 0

All variables are of unit order indicating that the variables have been appropriately
normalized.
For α = 8µ/D, A = πD 2 /4, and P = πD, we get
 4
1 Gr D
G = (7.61)
8192π P r L
 2
1 D
K = (7.62)
32 P r L
where the Prandtl and Grashof numbers are
µcp
Pr = (7.63)
k
qm gβL3
Gr = (7.64)
ν2k
107
respectively. Often the Rayleigh number defined by

Ra = Gr P r (7.65)

is used instead of the Grashof number.

Since u∗ is of O(1), the dimensional velocity is of order (8νL/D 2 )Gr 1/2 . The
ratio of axial conduction to the advective transport term is

K
 = (7.66)
G1/2
 
8π 1/2
= (7.67)
Ra

Taking typical numerical values for a loop with water to be: ρ = 998 kg/m3 , µ =
1.003 × 10−3 kg/m s, k = 0.6 W/m K, qm = 100 W/m, g = 9.91 m/s2 , β = 0.207 ×
10−3 K−1 , D = 0.01 m, L = 1 m, cp = 4.18 × 103 J/kgK, we get the velocity and
temperature scales to be

V G1/2 = (7.68)
∆T G1/2 = (7.69)

and the nondimensional numbers as

G = 1.86 × 10−2 (7.70)

K = 4.47 × 10−7 (7.71)
Gr = 3.35 × 1011 (7.72)
Ra = 2.34 × 1012 (7.73)
 = 3.28 × 10−6 (7.74)

Axial conduction is clearly negligible in this context.

For a steady state, equations (7.56) and (7.57) are
Z 1 ∗
u∗ = T g̃ ∗ ds∗ (7.75)
0
∗ ∗
d2 T dT
 ∗2 − u∗ ∗ = −q ∗ (s∗ ) (7.76)
ds ds
∗
Integrating over the loop from s∗ = 0 to s∗ = 1, we find that continuity of T and
∗
equation (7.20) imply continuity of dT /ds∗ also.

108
Conduction-dominated flow
If λ = G1/2 /  1, axial conduction dominates. We can write

u = u0 + λu1 + λ2 u2 + . . . (7.77)
T (s) = T 0 (s) + λT 1 (s) + λ2 T 2 (s) + . . . (7.78)

where, for convenience, the asterisks have been dropped. Substituting into the gov-
erning equations, and collecting terms of O(λ0), we have
Z 1
u0 = T 0 g̃ ds (7.79)
0
d2 T 0
= 0 (7.80)
ds2

The second equation, along with conditions that T 0 and dT 0 /ds have the same value
at s = 0 and s = 1, gives T 0 = an arbitrary constant. The first equation gives u0 = 0.
The terms of O(λ) give
Z 1
u1 = T 1 g̃ ds (7.81)
0
d2 T 1 dT 0
2
= −q(s) + u0 (7.82)
ds ds
The second equation can be integrated once to give
Z
dT 1 s
=− q(s0 ) ds0 + A (7.83)
ds 0

and again "Z #

Z s s00
T1 = − q(s0 ) ds0 ds00 + As + B (7.84)
0 0

Continuity of T 1 (s) and dT 1 /ds at s = 0 and s = 1 give

Z "Z #
1 s00
0 0
B = − q(s ) ds ds00 + A + B (7.85)
0 o

A = A (7.86)

respectively, from which "Z #

Z 1 s00
0 0
A= q(s ) ds ds00 (7.87)
0 0

109
and that B can be arbitrary. Thus
Z "Z # Z "Z #
s s00 1 s00
T1 = q(s0 ) ds0 ds00 + s q(s0 ) ds0 ds00 + T1 (0) (7.88)
0 0 0 0

where T (0) is an arbitrary constant. Substituting in equation (7.81), gives

Z (Z "Z # ) Z "Z # Z
1 s s00 1 s00 1
0 0 00 0 0
u1 = − q(s ) ds ds g̃ ds + q(s ) ds ds00 sg̃ ds (7.89)
0 0 0 0 0 0

The temperature distribution is determined by axial conduction, rather than by the

advective velocity, so that the resulting solution is unique.

Advection-dominated flow
The governing equations are
Z 1
u = T g̃ ds (7.90)
0
dT d2 T
u = q+ (7.91)
ds ds2
where   1. Expanding in terms of , we have

u = u0 + u1 + 2 u2 + . . . (7.92)
T = T 0 + T 1 + 2 T 2 + . . . (7.93)
(7.94)

To O(0 ), we get
Z 1
u0 = T 0 g̃ ds (7.95)
0
dT 0
u0 = q (7.96)
ds
from which
1 Zs
T0 = q(s0 ) ds0 (7.97)
u0 0
Z 1 Z s 
0 0
u0 = ± q(s ) ds g̃ ds (7.98)
0 0

Axial conduction. therefore, slightly modifies the two solutions obtained without it.

110
7.2.3 Toroidal geometry
The dimensional gravity function can be expanded in a Fourier series in s, to give
∞ 
X 
2πns 2πns
g̃(s) = gnc cos + g1s sin (7.99)
n=1 L L

The simplest loop geometry is one for which we have just the terms
2πs 2πs
g̃(s) = g1c cos + g1s sin (7.100)
L L
corresponds to a toroidal geometry. Using

g 2 = (g1c )2 + (g1s )2 (7.101)

gc
φ0 = tan−1 1s (7.102)
g1

equation (7.100) becomes  

2πs
g̃(s) = g cos − φ0 (7.103)
L
Without loss of generality, we can measure the angle from the horizontal, i.e. from
three o’clock point, and take φ0 = 0 so that

g̃(s) = g cos (2πs/L) (7.104)

The nondimensional gravity component is

g̃ = cos(2πs) (7.105)

where the * has been dropped.

Assuming also a sinusoidal distribution of heating

q(s) = − sin(2πs − φ) (7.106)

the momentum and energy equations are

Z 1
u = T (s)g̃(s) ds (7.107)
0
dT d2 T
u = q+ (7.108)
ds ds2
The homogeneous solution is
T h = Beus/ + A (7.109)

111
The particular integral satisfies
d2 Tp u dT p 1
− = sin(2πs − φ) (7.110)
ds2  ds 
Integrating, we have
dTp u 1
− Tp = − cos(2πs − φ) (7.111)
ds  2π
1
= − [cos(2πs) cos φ + sin(2πs) sin φ] (7.112)
2π
Take
T p = a cos(2πs) + b sin(2πs) (7.113)
from which
dT p
= −2πa sin(2πs) + 2πb cos(2πs) (7.114)
ds
Substituting and collecting the coefficients of cos(2πs) and sin(2πs), we get
u cos φ
− a + 2πb = − (7.115)
 2π
u sin φ
−2πa − b = − (7.116)
 2π
The constants are
(u/2π2 ) cos φ + (1/) sin φ
a = (7.117)
4π 2 + u2 /2
−(1/) cos φ + (u/2π2 ) sin φ
b = (7.118)
4π 2 + u2 /2
The temperature field is given by

T = Th + Tp (7.119)

Since T (0) = T (1), we must have B = 0. Taking the other arbitrary constant A to
be zero, we have
    
1 u 1 1 u
T = 2 cos φ + sin φ cos(2πs) + − cos φ + sin φ sin(2πs)
4π + u2 /2 2π2   2π2
(7.120)
The momentum equation gives
(u/2π2 ) cos φ + (1/) sin φ
u= (7.121)
2(4π 2 + u2 /2 )

112
which can be written as
 
3 2 2 1 
u + u 4π  − cos φ − sin φ = 0 (7.122)
4π 2

Special cases are:

• =0
Equations (7.107) and (7.108) can be solved to give

1
T = [cos(2πs) cos φ + sin(2πs) sin φ] (7.123)
2πu
s
cos φ
u = ± (7.124)
4π

On the other hand substituting  = 0 in equation (7.122) gives an additional

spurious solution u = 0.

• →∞
We get that u → 0.

• φ=0
We get 

 0
 q
1
u= 4π
− 4π 2 2 (7.125)

 q
 − 1 − 4π 2 2
4π

The last two solutions exist only when  < (16π 3 )−1/2 .

• φ = π/2
The velocity is a solution of

u3 + u4π 2 2 − =0 (7.126)
2

Figure 7.8 shows u-φ curves for three different values of . Figure 7.9 and 7.10
show u- curves for different values of φ. It is also instructive to see the curve u-
Ra, shown in Figure 7.11, since the Rayleigh number is directly proportional to the
strength of the heating.

113
 1 ε = 0.001
u
0.5
0
-0.5
-1
-3 -2 -1 0 1 2 3

1
u ε = 0.01
0.5
0
-0.5
-1
-3 -2 -1 0 1 2 3
φ

1
u ε = 0.1
0.5
0
-0.5
-1
-3 -2 -1 0 1 2 3

Figure 7.8: u-φ curves.

114
 0.2 φ=0
u 0.1

0.02 0.04 0.06 0.08 0.1

-0.1 ε

-0.2

φ=0.01
0.2
u
0.1

0.02 0.04 0.06 0.08 0.1

-0.1 ε

-0.2

φ=−0.01
0.2
u
0.1

0.02 0.04 0.06 0.08 0.1

-0.1 ε

-0.2

Figure 7.9: u- curves.

115
 φ=π/4
u 0.2
0.1

0.02 0.04 0.06 0.08

-0.1
-0.2

0.2

u
0.15

0.1 φ=π/4

0.05

0.02 0.04 0.06 0.08 0.1

0.1
u

0.08

0.06

0.04 φ=π/4

0.02

0.02 0.04 0.06 0.08 0.1

Figure 7.10: More u- curves.

116
 u 0.2

0.1

10000 20000 30000 40000

Ra
-0.1

-0.2

Figure 7.11: u-Ra for φ = 0.01 radians.

The bifurcation set is the line dividing the regions with only one real solutions
and that with three real solutions. A cubic equation
x3 + px + q = 0 (7.127)
has a discriminant
p3 q 2
D= + (7.128)
27 4
For D < 0, there are three real solutions, and for D > 0, there is only one. The
discriminant for the cubic equation (7.122) is
 3
1 1 1
D= 4π 2 2 − cos φ + ( sin φ)2 (7.129)
27 4π 4
The result is shown in Fig. 7.12.

7.2.4 Dynamic analysis

We rescale the nondimensional governing equations (7.56) and (7.57) by
1
u∗ = ub (7.130)
2πG1/2
117
 0.045

0.04

ε 0.035

0.03

0.025

0.02

0.015

0.01

0.005

0
−100 −80 −60 −40 −20 0 20 40 60 80 100

φ (degrees)
Figure 7.12: Region with three solutions.

118
 1
T∗ = Tb (7.131)
2πG1/2
to get
1 Z
du
+u = T g̃ ds (7.132)
dt 0
∂T 1 ∂T G1/2 K ∂ 2 T
+ u = Gq + (7.133)
∂t 2π ∂s 2π ∂s2
where the hats and stars have been dropped.
We take g̃ = cos(2πs) and q = − sin(2πs − φ). Expanding the temperature in a
Fourier series, we get
∞
X
T (s, t) = T0 (t) + [Tnc (t) cos(2πns) + Tns (t) sin(2πns)] (7.134)
n=1

Substituting, we have
du 1
+ u = T1c (7.135)
dt 2
and
∞
" #
dT0 X dTnc dTnc
+ cos(2πns) + sin(2πns)
dt n=1 dt dt
∞
X
+u [−nTnc sin(2πns) + nTns cos(2πns)]
n=1
= −G [sin(2πs) cos φ − cos(2πs) sin φ]
∞
X
2 1/2
−2πn G K [Tnc cos(2πns) + Tns sin(2πns)] (7.136)
n=1

Integrating, we get
dT0
=0 (7.137)
dt
Multiplying by cos(2πms) and integrating

1 dTmc m 1
+ uTms = G sin φ − πm2 G1/2 KTmc (7.138)
2 dt 2 2
Now multiplying by sin(2πms) and integrating

1 dTms m 1
− uTmc = − G cos φ − πm2 G1/2 KTms (7.139)
2 dt 2 2
119
Choosing the variables
x = u (7.140)
1 c
y = T (7.141)
2 1
1 s
z = T (7.142)
2 1
and the parameters
G
a = sin φ (7.143)
2
G
b = cos φ (7.144)
2
c = 2πG1/2 K (7.145)
we get the dynamical system
dx
= y−x (7.146)
dt
dy
= a − xz − cy (7.147)
dt
dz
= −b + xy − cz (7.148)
dt
The physical significance of the variables are: x is the fluid velocity, y is the horizontal
temperature difference, and z is the vertical temperature difference. The parameter
c is positive, while a and b can have any sign.
The critical points are found by equating the vector field to zero, so that
y−x = 0 (7.149)
a − xz − cy = 0 (7.150)
−b + xy − cz = 0 (7.151)
From equation (7.149), we have y = x, and from equation (7.151), we get z =
(−b + x2 )/c. Substituting these in equation (7.150), we get
x3 + x(c2 − b) − ac = 0 (7.152)
This corresponds to equation (7.122), except in different variables.
To analyze the stability of a critical point (x, y, z) we add perturbations of the
form
x = x + x0 (7.153)
y = y + y0 (7.154)
z = z + z0 (7.155)

120
x

b
2
c

Figure 7.13: Bifurcation diagram for x.

121
Substituting in equation (7.146)-(7.148), we get the local form
      
x0 −1 1 0 x0 0
d  0    0   
 y  =  −z −c −x   y  +  −x0 z 0  (7.156)
dt
z0 y x −c z0 x0 y 0

The linearized version is

    
x0 −1 1 0 x0
d  0    0 
 y  =  −z −c −x   y  (7.157)
dt
z0 y x −c z0

No tilt, with axial conduction (a = 0, c 6= 0)

From equation (7.152), for a = 0 we get

x3 + x(c2 − b) = 0 (7.158)

from which 

 0
√
x=y= b − c2
√ (7.159)


− b − c2
The z coordinate is 
 −b/c

z= −c (7.160)


−c
The bifurcation diagram is shown in Figure 7.13.

Stability of conductive solution

The critical point is (0, 0, −b/c). To examine its linear stability, we look at the
linearized equation (7.157) to get
    
x0 −1 1 0 x0
d  0    0 
 y  =  b/c −c 0   y  (7.161)
dt
z0 0 0 −c z0

The eigenvalues of the matrix are obtained from the equation

−(1 + λ) 1 0
b/c −(c + λ) 0 =0 (7.162)
0 0 −(c + λ)

122
which simplifies to
" #
b
(c + λ) (1 + λ)(c + λ) − =0 (7.163)
c
One eigenvalue is
λ1 = −c (7.164)
Since c ≥ 0 this eigenvalue indicates stability. The other two are solutions of

b
λ2 + (c + 1)λ + (c − ) = 0 (7.165)
c
which are
 s 
1 b
λ2 = −(c + 1) − (c + 1)2 − 4(c − ) (7.166)
2 c
 s 
1 b
λ3 = −(c + 1) + (c + 1)2 − 4(c − ) (7.167)
2 c

λ2 is also negative and hence stable. λ3 is negative as long as

s
b
−(c + 1) + (c + 1)2 − 4(c − ) < 0 (7.168)
c

which gives
b < c2 (7.169)
This is the condition for stability.
In fact, one can also prove global stability of the conductive solution. Restoring
the nonlinear terms in equation (7.156) to equation (7.161), we have

dx0
= y 0 − x0 (7.170)
dt
dy 0 b
= x0 − cy 0 − x0 z 0 (7.171)
dt c
dz 0
= −cz 0 + x0 y 0 (7.172)
dt
Let
b
V (x, y, z) = x02 + y 02 + z 02 (7.173)
c
123
Thus
1 dV b 0 dx0 dy 0 dz 0
= x + y0 + z0 (7.174)
2 dt c dt dt dt
b 02 2b 0 0
= − x + x y − cy 02 − cz 02 (7.175)
c c
b b
= − (x0 − y 0 )2 − (c − )y 02 − cz 02 (7.176)
c c
Since
V ≥ 0 (7.177)
dV
≤ 0 (7.178)
dt
for 0 ≤ b ≤ c2 , V is a Liapunov function, and the critical point is stable to all
perturbations in this region. The bifurcation at b = c2 is thus supercritical.
Stability of convective solution √ √
For b > c2 , only one critical point ( b − c2 , b − c2 , −c) will be considered, the
other being similar. We use the linearized equations (7.157). Its eigenvalues are
solutions of
−(1 + λ) 1 √0
√ c −(c
√ + λ) − b − c = 0
2 (7.179)
b−c 2 b−c 2 −(c + λ)
This can be expanded to give
λ3 + λ2 (1 + 2c) + λ(b + c) + 2(b − c2 ) = 0 (7.180)
The Hurwitz criteria for stability require that all coefficients be positive, which they
are. Also the determinants
D1 = 1 + 2c (7.181)
1 + 2c 2(b − c2 )
D2 = (7.182)
1 b+c
1 + 2c 2(b − c2 ) 0
D3 = 1 b+c 0 (7.183)
0 1 + 2c 2(b − c2 )
should be positive. This requires that
c(1 + 4c)
b< if c < 1/2 (7.184)
1 − 2c
c(1 + 4c)
b> if c > 1/2 (7.185)
1 − 2c
(7.186)

124
With tilt, no axial conduction (a 6= 0, c = 0)
The dynamical system (7.146)-(7.148) simplifies to
dx
= y−x (7.187)
dt
dy
= a − xz (7.188)
dt
dz
= −b + xy (7.189)
dt
√ √ √ +
The √
critical √ are ±( b, b, a/ b). The linear stability of the point P given
√ points
by ( b, b, a/ b) will be analyzed. From equation (7.157), the solutions of
−(1 +√λ) 1 0√
−a/
√ b −λ √ − b =0 (7.190)
b b −λ
are the eigenvalues. This simplifies to
a
λ3 + λ2 + λ(b + √ ) + 2b = 0 (7.191)
b
For stability the Hurwitz criteria require all coefficients to be positive, which they
are. The determinants

D1 = 1 (7.192)
1 2b √
D2 = (7.193)
1 b + a/ b
1 2b √ 0
D3 = 1 b + a/ b 0 (7.194)
0 1 2b
√
should also be positive. This gives the condition (b + a/ b) − 2b > 0, from which, we
have
a > b3/2 (7.195)
for stability. The stable and unstable region for P + is shown in Figure
√ √ 7.14. √
Also
−
shown is the stability of the critical point P with coordinates −( b, b, a/ b).
The dashed circles are of radius G/2, and the angleof tile φ is also indicated. Using
equations (7.143) and (7.144), the stability condition (7.195) can be written as
 1/2
sin φ G
> (7.196)
cos3/2 φ 2

125
 a 3/2
a=b

+
P stable _
Both P + and P
G/2 unstable
φ
b

_
P stable

3/2
a=-b
Figure 7.14: Stability of critical points P + and P − in (b, a) space.

As a numerical example, for the value of G in equation (7.70), P + is stable for the
tilt angle range φ > 7.7◦ , and P − is stable for φ < −7.7◦ . In fact, for G  1, the
stability condition for P + can be approximated as
 1/2
G
φ> (7.197)
2
The same information can be shown in slightly different coordinates. Using
√
x = b for P + and equation (7.144), we get

G x2
= (7.198)
2 cos φ

126
The stability condition (7.195) thus becomes

tan φ < x (7.199)

The stability regions for both P + and P − are shown in Figure 7.15.
The loss of stability is through imaginary eigenvalues. In fact, for P + , substi-
tuting a = b3/2 in √equation (7.191), the equation can be factorized to give the three
eigenvalues −1, ±i 2b. Thus the nondimensional
√ radian frequency of the oscillations
in the unstable range is approximately 2b.
The effect of a small nonzero axial conduction parameter c is to alter the Figure
7.195 in the zone 0 < b < c.

7.2.5 Nonlinear analysis

Numerical

Let us choose b = 1, and reduce a. Figures 7.16 and 7.17 show the x-t and phase
space representation for a = 0.9, Figures 7.18 and 7.19 for a = 0.55, and Figures 7.20
and 7.21 for a = 0.53.
The strange attractor is shown in Figures 7.22 and 7.23.
Comparison of the three figures in Figures 7.24 shows that vestiges of the shape
of the closed curves for a = −0.9 and a = 0.9 can be seen in the trajectories in a = 0.

Analytical

The following analysis is by W. Franco.

We start with the dynamical system which models a toroidal thermosyphon loop
with known heat flux

dx
= y−x
dt
dy
= a − zx (7.200)
dt
dz
= xy − b
dt

For b > 0 two critical points P + and P − appear

!
√ √ a
(x̄, ȳ, z̄) = ± b, b, √ (7.201)
b

127
 x

P +exists
but unstable

P+stable

−π −π π π
φ

_
P stable

_
P exists
but unstable

Figure 7.15: Stability of critical points P + and P − in (φ, x) space.

128
 1.8

1.6

1.4

1.2

1
x

0.8

0.6

0.4

0.2
0 5 10 15 20 25 30
t

Figure 7.16: x-t for a = 0.9, b = 1.

129
 2.5

1.5

1
z

0.5

−0.5
2.5
2
2
1.5
1.5
1
0.5 1
0 0.5
−0.5 0
y
x

Figure 7.17: Phase-space trajectory for a = 0.9, b = 1.

130
 2.5

1.5

1
x

0.5

−0.5

−1
0 5 10 15 20 25 30
t

Figure 7.18: x-t for a = 0.55, b = 1.

131
 4

1
z

−1

−2
4
3
2.5
2 2
1 1.5
1
0 0.5
−1 0
−0.5
−2 −1
y
x

Figure 7.19: Phase-space trajectory for a = 0.55, b = 1.

132
 2.5

1.5

1
x

0.5

−0.5

−1
0 5 10 15 20 25 30
t

Figure 7.20: x-t for a = 0.53, b = 1.

133
 4

1
z

−1

−2

−3
4
3
2.5
2 2
1 1.5
1
0 0.5
−1 0
−0.5
−2 −1
y
x

Figure 7.21: Phase-space trajectory for a = 0.53, b = 1.

134
 3

0
x

−1

−2

−3

−4
0 50 100 150 200 250 300
t

Figure 7.22: x-t for a = 0, b = 1.

135
 6

0
z

−2

−4

−6
5

−5 3
2
1
0
−1
−2
y −10 −3
−4
x

Figure 7.23: Phase-space trajectory for a = 0, b = 1.

136
 6

z
−2 a=0
−4

−6
4

2 4
2
0
0
−2
−2
−4 −4
y
x

0
z

−2
a = 0.9
−4

−6
4

2 4
2
0
0
−2
−2
−4 −4
y
x

0
z

−2 a = - 0.9
−4

−6
4

2 4
2
0
0
−2
−2
−4 −4
y
x

Figure 7.24: Phase-space trajectories for b = 1.

137
The local form respect to P + is

dx0
= y 0 − x0
dt
dy 0 a √
= − √ x0 − bz 0 (7.202)
dt b
dz 0 √ √
= bx0 + by 0
dt
3 3 √
For stability a > b 2 . At a = b 2 the eigenvalues are −1,± 2bi, thus a nonlinear analy-
sis through the center manifold projection is possible. Let’s introduce a perturbation
3
of the form a = b 2 +  and the following change of variables
a
α = √
b
√
β = b

rewriting the local form, dropping the primes and regarding the perturbation the
system becomes

dx
= y−x
dt !
dy 2 
= β + x − βz (7.203)
dt β
dz
= βx − βy
dt
for stability α > β 2 .
Let’s apply the following transformation:
√
2 2β 2
x = w1 + 2 w2 + 2
2β + 1 2β + 1
y = 2w2 (7.204)
√
2β 2 2 (β 2 + 1)
z = −βw1 − 2 w2 + w3
2β + 1 2β 2 + 1

in the new variables

 
−1 0 √0

ẇ =  0 √0 − 2β 
 w + P̂w + l(w) (7.205)
0 2β 0

138
 The center manifold projection is convenient to use if the large-time dynamic
behavior is of interest. In many dimensional systems, the system often settles into
the same large-time dynamics irrespective of the initial condition; this is usually
less complex than the initial dynamics and can be described by far simple evolution
equations.
We first state the definition of an invariant manifold for the equation

ẋ = N(x) (7.206)

where x ∈ Rn . A set S ⊂ Rn is a local invariant manifold for (7.206) if for x0 ⊂ S,

the solution x(t) of (7.206) is in S for | t |< T where T > 0. If we can always choose
T = ∞, then S is an invariant manifold. Consider the system

ẋ = Ax + f (x, y)
ẏ = By + g(x, y) (7.207)

where x ∈ Rn , y ∈ Rm and A and B are constant matrices such that all the eigenvalues
of A have zero real parts while all the eigenvalues of B have negative real parts. If
y = h(x) is an invariant manifold for (7.207) and h is smooth, then it is called a
center manifold if h(0) = 0,h0 (0) = 0. The flow on the center manifold is governed by
the n-dimensional system
ẋ = Ax + f (x, h(x)) (7.208)
The last equation contains all the necessary information needed to determine the
asymptotic behavior of small solutions of (7.207).
Now we calculate, or at least approximate the center manifold h(w). Substituting
w1 = h(w2 , w3 ) in the first component of (7.205) and using the chain rule, we obtain
 
! ẇ2
∂h ∂h  
ẇ1 = ,   = −h + l1 (w2 , w3 , h) (7.209)
∂w2 ∂w3
ẇ3

We seek a center manifold

h = aw22 + bw2 w3 + cw32 + O(3) (7.210)

substituting in (7.209)
 √ 
− 2βw3
 
(2aw2 + bw3 , 2cw3 + bw2 )  √  =
2βw2
   
− aw22 + bw2 w3 + cw32 + k1 w22 + k2 w2 w3 + k3 w32 + O(3)

139
Equating powers of x2 ,xy and y 2 , we find that
√
a = k1 − b 2β
√
c = k3 + b 2β
√
k2 + 2 2β (k1 − k3 )
b =
8β 2 + 1

The reduced system is therefore given by

√
ẇ2 = − 2βw3 + s2 (w2 , w3 )
√
ẇ3 = 2βw2 + s3 (w2 , w3 ) (7.211)

Normal form: Now we carry out a smooth nonlinear coordinate transform of the
type
w = v + ψ(v) (7.212)
to simplify (7.211) by transforming away many nonlinear terms. The system in the
new coordinates is
 
√ ! (νv1 − γv2 ) (v12 + v22 )
√0 − 2β  
v̇ = +  (7.213)
2β 0 2 2
(νv2 + γv1 ) (v1 + v2 )

where ν and γ depend on the nonlinear part of (7.211). This is the unfolding of the
Hopf bifurcation.
Although the normal form theory presented in class pertains to a Jacobian whose
eigenvalues all lie on the imaginary axis, one can also present a perturbed version.
The eigenvalues are then close to the imaginary axis but not quite on it. Consider
the system
v̇ = Av + Âv + f(v) (7.214)
where the Jacobian A has been evaluated at a point in the parameter space where all
its eigenvalues are on the imaginary axis, Â represents a linear expansion of order µ in
the parameters above that point; a perturbed Jacobian. The perturbation parameter
represents the size of the neighborhood in the parameter space. We stipulate the
order of µ such that the real part of the eigenvalues of A + Â is such that, to leading
order, Â does not change the coefficients of the leading order nonlinear terms of
the transformed equation. The linear part A + Â of perturbed Hopf can always be
transformed to !
µ −ω
(7.215)
ω µ

140
The required transformation is a near identity linear transformation

v = u + Bu (7.216)

such that the linear part of (7.214) is transformed to

 
u̇ = A + AB − BA + Â z (7.217)

For the Hopf bifurcation if !

a1 a2
 = (7.218)
a3 a4
then
a1 + a4
µ= (7.219)
ω
Therefore for  small we can write (7.213) as
 
√ ! ! f1 (v)
0 − 2β p22 p23
v̇ = √ v+ v+


 (7.220)
2β 0 p32 p33
f2 (v)

where the perturbation matrix comes from (7.205). Applying a near identity trans-
formation of the form v = u + Bu the system becomes
 
√ ! (νu1 − γu2 ) (u21 + u22)
µ − 2β
u̇ = √ u+


 (7.221)
2β µ 2 2
(νu2 + γu1 ) (u1 + u2 )

which is the unfolding for the perturbed Hopf bifurcation. In polar coordinates we
have

ṙ = µr + νr 3
√
θ̇ = 2β (7.222)

where √
 2
µ=− 2 (7.223)
2β (2β 2 + 1)
40β 6 + 40β 4 + 12β 3 + 10β 2 + 12β + 3
ν=− (7.224)
4 (8β 2 + 1) (2β 2 + 1)4
Appendix

141
 8β (β 2 + 1)
k1 = −
(2β 2 + 1)3
 √ √ 
(β 2 + 1) 4 2 − 8 2β 2
k2 = (7.225)
(2β 2 + 1)3
8β (β 2 + 1)
k3 =
(2β 2 + 1)3


p22 = −
β (2β 2 + 1)
√
 2
p23 = − 2 (7.226)
2β + 1
From the literature
1
ν = (fxxx + fxyy + gxxy + gyyy )
16
1
+ (fxy (fxx + fyy ) − gxy (gxx − gyy ) − fxx gxx − fyy gyy ) (7.227)
16ω
√
in our problem f = f1 , g = f2 , x = v1 , y = v2 and ω = 2β.

7.3 Known wall temperature

The heating is now convective with a heat transfer coefficient U, and an external
temperature of Tw (s). Thus,
q = P U(T − Tw ) (7.228)
Neglecting axial conduction
Z
Pα β L
u = T (s)g̃(s) ds (7.229)
ρ0 A L 0
dT h i
u = γ T − Tw (s) (7.230)
ds
where γ = UP/ρ0 Acp 1 . Multiplying the second equation by e−γs/u /u, we get
d  −γs/u  γ
e T = − e−γs/u Tw (7.231)
ds u
1
The sign of γ appears to be wrong.

142
Integrating, we get
 Z s 
γ −γs0 /u 0 0
T =e γs/u
− e Tw (s ) ds + T0 (7.232)
u 0

Since T (L) = T (0), we get

RL 0
γ eγL/u 0 e−γs /u Tw (s0 ) ds0
T0 = (7.233)
u eγL/u − 1

The velocity is obtained from

Z L   Z s 
Pα β γ −γs0 /u 0 0
u= γs/u
e − e Tw (s ) ds + T0 g̃(s) ds (7.234)
ρ0 A L 0 u 0

This is a transcendental equation that may have more than one real solution.

Example 7.3
Show that there is no motion if the wall temperature is uniform.

Take Tw to be a constant. Then equation (7.230) can be written as

d(T − Tw ) γ
= ds (7.235)
T − Tw u

The solution to this is

T = Tw + K eγs/u (7.236)

where K is a constant. Continuity of T at s = 0 and s = L gives K = 0. Hence T = Tw , and,

from equation (7.229), u = 0.

Assume the wall temperature to be

Tw (s) = − sin(2πs − φ) (7.237)

The temperature field is

! !
b cos(2πs − φ) cos(2πs − φ)
T = − (7.238)
r2 − r1 2π(1 + r22 /4π 2 ) 2π(1 + r1 62/4π 2)

143
7.4 Mixed condition
The following has been written by A. Pacheco-Vega.
It is common, especially in experiments, to have one part of the loop heated with
a known heat rate and the rest with known wall temperature. Thus for part of the
loop the wall temperature is known so that q = P U(T − Tw (s)), while q(s) is known
for the rest. As an example, consider
(
P U(T − T0 ) for 2π
φ
≤ s ≤ π + 2π
φ
q= φ φ (7.239)
q0 for π + 2π < s < 2π + 2π

where T0 and q0 are constants.

7.4.1 Modeling
If we consider a one-dimensional incompressible flow, the equation of continuity indi-
cates that the velocity v is a function of time alone. Thus,
v = v(t). (7.240)
Taking an infinitesimal cylindrical control volume of fluid in the loop πr 2 dθ, see Figure
(7.25), the momentum equation in the θ-direction can be written as
dv dp
ρπr 2 Rdθ = −πr 2 dθ − ρgπr 2 Rdθ cos(θ + α) − τw 2πrRdθ (7.241)
dt dθ
Integrating Eq. (7.241) around the loop using the Boussinesq approximation ρ =
ρw [1 − β(T − Tw )], with the shear stress at the wall being approximated by that
corresponding to Poiseuille flow in a straight pipe τw = 8µv/ρw r 2 , the expression of
the balance in Eq. (7.241) modifies to
Z 2π
dv 8µ βg
+ 2
v= (T − Tw ) cos(θ + α) dθ (7.242)
dt ρw r 2π 0

Neglecting axial heat conduction, the temperature of the fluid satisfies the following
energy balance equation

! 2h
∂T v ∂T  − r (T − Tw ), 0 ≤ θ ≤ π

ρw cp + = (7.243)
∂t R ∂θ 
 2
r
q, π < θ < 2π
Following the notation used by Greif et al. (1979), the nondimensional time, velocity
and temperature are defined as
t v T − Tw
τ= , w= , φ= (7.244)
2πR/V V q/h

144
respectively, where
!1/2
gβRrq
V = . (7.245)
2πcp µ
Accordingly, Eqs. (7.242) and (7.243) become
Z 2π
dw πΓ
+ Γw = φ cos(θ + α) dθ (7.246)
dτ 4D 0

and (
∂φ ∂φ −2Dφ, 0 ≤ θ ≤ π
+ 2πw = (7.247)
∂τ ∂θ 2D, π < θ < 2π
where the parameters D and Γ are defined by
2πRh 16πµR
D= Γ= (7.248)
ρw cp rV ρw r 2 V

7.4.2 Steady State

The steady-state governing equations without axial conduction are
Z 2π
π
w= φ cos(θ + α) dθ (7.249)
4D 0

and 
 − πw φ, 0 ≤ θ ≤ π
D

dφ
= (7.250)
dθ   D
πw
, π < θ < 2π
where w and φ are the steady-state values of velocity and temperature respectively.
Eq. (7.250) can be integrated to give

−(Dθ/πw)

 Ae , 0≤θ≤π
φ(θ) = (7.251)

 D
πw
θ + B, π < θ < 2π
Applying the condition of continuity in the temperature, such that φ(0) = φ(2π) and
φ(π − ) = φ(π + ) the constants A and B can be determined. These are
" #
D 1 D 2 e(−D/w) − 1
A= B= (7.252)
w 1 − e(−D/w) w 1 − e(−D/w)
The resulting temperature filed is

 D e−(Dθ/πw)

 w 1−e−(D/w)
, 0≤θ≤π
φ(θ) =  h i (7.253)

 D θ 2e−(D/w) −1
w π
+ 1−e−(D/w)
, π < θ < 2π

145
Substituion of Eq. (7.250) in Eq. (7.249), followed by an expansion of cos(θ + α),
leads to
(Z
π π D e−(Dθ/πw)
w= cos α cos θ dθ
4D 0 w 1 − e−(D/w)
Z " # )
2π D θ 2e−(D/w) − 1
+ + cos θ dθ
π w π 1 − e−(D/w)
(Z
π π D e−(Dθ/πw)
− sin α sin θ dθ
4D 0 w 1 − e−(D/w)
Z " # )
2π D θ 2e−(D/w) − 1
+ + sin θ dθ (7.254)
π w π 1 − e−(D/w)

and integration around the loop, gives the steady-state velocity as

!
2 cos α (D/w) cos α + π(D/w)2 sin α 1 + e−(D/w)
w = +   2  . (7.255)
2 4 1 + πwD 1 − e−(D/w)

As a final step, multiplying the numerator and denominator by e(D/2w) and rearrang-
ing terms leads to the expresion for the function of the steady-state velocity

cos α (D/w) cos α + π(D/w)2 sin α

G(w, α, D) = w2 − −   2  coth(D/2 w) = 0
2 4 1 + πwD

(7.256)

For α = 0, symmetric steady-state solutions for the fluid velocity are possible since
G(w, 0, D) is an even function of w. In this case Eq.(7.256) reduces to
h i
1 (D/w) 1 + e−(D/w)
w2 = +     . (7.257)
2 4 1 + D 2 [1 − e−(D/w) ]
πw

The steady-state solutions of the velocity field and temperature are shown next.
Figure 7.26 shows the w −α curves for different values of the parameter D. Regions of
zero, one, two and three solutions can be identified. The regions of no possible steady-
state velocity are: −180◦ < α < −147.5◦ and 147.5◦ < α < 180◦ . There is only one
velocity for the ranges −147.5◦ < α < −α0 and α0 < α < 147.5◦ where α0 varies from
90◦ at a value of D = 0.001 to α0 = 32.5◦ when D = 100. Three velocities are obtained
for −α0 < α < −32.5◦ and 32.5◦ < α < α0 , except for the zero-inclination case which
has two possible steady-state velocities. The temperature distribution in the loop,

146
for three values of the parameter D and α = 0 is presented in Figure 7.27. From the
φ−θ curves it can be seen the dependence of the temperature with D. As D increases
the variation in temperature between two opposit points also increases. When has
a value D = 0.1 the heating and cooling curves are almost straight lines, while at
a value of D = 1.0 the temperature decays exponentially and rises linearly. Similar
but more drastic change in temperature is seen when D = 2.5. Figure 7.28 shows
the φ − θ curves for three different inclination angles with D = 2.5. It can be seen
the increase in the temperature as α takes values of α = 0◦ , α = 90◦ and α = 135◦ .
This behaviour is somewhat expected since the steady-state velocity is decreasing in
value such that the fluid stays longer in both parts of the loop. Figure 7.29 shows
the steady-state velocity as a function of D for different angles of inclination α. For
α = 0 we have two branches of the velocity-curve which are symmetric. The positive
and negative values of the velocity are equal in magnitude for any value of D. For
α = 45◦ , the two branches are not symmetric while for α = 90◦ and α = 135◦, only
the positive branch exist.

7.4.3 Dynamic Analysis

The temperature can be expanded in Fourier series, such that
∞
X
φ = φ0 + [φcn (t) cos(nθ) + φsn (t) sin(nθ)] (7.258)
n=1

Substituiting into Eqs. (7.246) and (7.247), we have

dw π2Γ π2Γ
+ Γw = cos αφc1 − sin αφs1 (7.259)
dτ 4D 4D
and
∞
" #
dφc0 X dφcn dφs
+ cos (nθ) + n sin (nθ)
dτ n=1 dτ dτ
∞
X
+2πw [−nφcn sin (nθ) + nφsn cos (nθ)]
n=1
 P∞
 −2D {φ0 + n=1 [φn (t) cos(nθ) + φn (t) sin (nθ)]} , 0 ≤ θ ≤ π
 c c s

= (7.260)


2D, π < θ < 2π

Integrating Eq. (7.260) from θ = 0 to θ = 2π we get

" ∞
#
dφc0 1X φsn
= −D φc0 − [(−1)n − 1] − 1 (7.261)
dτ π n=1 n

147
 Multiplying by cos (mθ) and integrating from θ = 0 to θ = 2π

dφcm D X∞ h i 2n
+ 2πm w φsm = −Dφcm + φsn (−1)m+n − 1 2 (7.262)
dτ π n=1 n − m2
n6=m

Now multiplying by sin (mθ) and integrating from θ = 0 to θ = 2π

dφsm D X∞ h i 2m
− 2πm w φcm = −Dφsm + φcn (−1)m+n − 1
dτ π n=0 m − n2
2
n6=m

2D
− [1 − (−1)m ] (7.263)
πm
for m ≥ 1.
Choosing the variables

w = w (7.264)
C0 = φc0 (7.265)
Cm = φcm (7.266)
Sm = φsm (7.267)

we get an infinte-dimensional dynamical system

dw π2Γ π2Γ
= −Γw + cos α C1 − sin α S1 (7.268)
dτ 4D 4D
∞
dC0 DX Sn
= −D C0 + [(−1)n − 1] + D (7.269)
dτ π n=1 n
dCm D X∞ h i 2n
= −2πm w Sm − D Cm + Sn (−1)m+n − 1 2 (7.270)
dτ π n=1 n − m2
n6=m

dSm D h ∞
X i 2m
= 2πm w Cm − D Sm + Cn (−1)m+n − 1
dτ π n=0 m − n2
2
n6=m

2D
+ [(−1)m − 1] (7.271)
πm
for m ≥ 1. The physical significance of the variables are: w is the fluid velocity, C
is the horizontal temperature difference, and S is the vertical temperature difference.
The parameters of the system are D, Γ and α. D and Γ are positive, while α can
have any sign.

148
 The critical points are found by equating the vector filed to zero, so that
π2 π2
w− cos α C 1 + sin α S 1 = 0 (7.272)
4D 4D
∞
1X Sn
(C 0 − 1) − [(−1)n − 1] = 0 (7.273)
π n=1 n
D X∞ h i 2n
2πm w S m + D C m − S n (−1)m+n − 1 2 = 0 (7.274)
π n=1 n − m2
n6=m

D X∞ h i 2m
2πm w C m − D S m + C n (−1)m+n − 1
π n=0 m − n2
2
n6=m

2D
[(−1)m − 1] = 0
+ (7.275)
πm
However, a convenient alternative way to determine the critical points is by using
a Fourier series expansion of the steady-state temperature field solution given in Eq.
(7.253). The Fourier series expansion is
∞ h
X i
φ= C n cos(nθ) + S n sin(nθ) (7.276)
n=0

Performing the inner product between Eq. (7.253) and cos(mθ) we have
Z
D 1 π
e−(Dθ/πw) cos(mθ) dθ
w 1 − e"−(D/w) 0 #
D Z 2π θ 2e−(D/w) − 1
+ + cos(mθ) dθ
w π π 1 − e−(D/w)
∞
X Z 2π ∞
X Z 2π
= Cn cos(nθ) cos(mθ) dθ + Sn sin(nθ) cos(mθ)dθ (7.277)
n=0 0 n=0 0

Now the inner product between Eq. (7.253) and sin(mθ) gives
Z
D 1 π
e−(Dθ/πw) sin(mθ) dθ
w 1 − e"−(D/w) 0 #
D Z 2π θ 2e−(D/w) − 1
+ + sin(mθ) dθ
w π π 1 − e−(D/w)
∞
X Z 2π ∞
X Z 2π
= Cn cos(nθ) sin(mθ) dθ + Sn sin(nθ) sin(mθ)dθ (7.278)
n=0 0 n=0 0

from which we get

" !#
1 D 3 2e−(D/w) − 1
C0 = 1+ + (7.279)
2 w 2 1 − e−(D/w)

149
  2
D
πw 1 − e−(D/w) cos(mπ) D
Cm =   2  + 2 2 [1 − cos(mπ)] (7.280)
m2 + D 1−e −(D/w) π mw
πw
  3 

 

 D
1 − e−(D/w) cos(mπ) 
πw
Sm = −   2  (7.281)

 1 − e−(D/w) 

 m m2 + D 
πw

To analyze the stability of a critical point (w, C 0 , C 1 , · · · , C m , S 1 , · · · , S m ) we add

perturbations of the form

w = w + w0 (7.282)
C0 = C 0 + C00 (7.283)
C1 = C 1 + C10 (7.284)
..
. (7.285)
0
Cm = C m + Cm (7.286)
S1 = S 1 + S10 (7.287)
..
. (7.288)
0
Sm = S m + S m (7.289)

Substituiting in Eqs. (7.268) to (7.271), we obtain the local form

dw 0 π2Γ π2Γ
= −Γw 0 + cos α C10 − sin α S10 (7.290)
dτ 4D 4D
∞
dC00 0 DX (−1)n − 1 0
= −D C0 + Sn (7.291)
dτ π n=1 n
0
dCm 0
= −2πm w Sm − 2πm S m w 0 − D Cm 0
dτ
D X∞
2n h i
+ (−1) m+n
− 1 Sn0 − 2πm w 0 Sm
0
m≥1
π n=1 n − m
2 2
n6=m

(7.292)
0
dSm 0
= 2πm w Cm + 2πm C m w 0 − D Sm
0
dτ
D X∞
2m h i
+ (−1) m+n
− 1 Cn0 + 2πm w 0 Cm
0
m≥1
π n=0 m − n
2 2
n6=m

(7.293)

150
The linearized version is
dw 0 π2Γ π2Γ
= −Γw 0 + cos α C10 − sin α S10 (7.294)
dτ 4D 4D
∞
dC00 DX (−1)n − 1 0
= −D C00 + Sn (7.295)
dτ π n=1 n
0
dCm 0
= −2πm w Sm − 2πm S m w 0 − D Cm 0
dτ
D X∞
2n h i
+ (−1) m+n
− 1 Sn0 m≥1 (7.296)
π n=1 n2 − m2
n6=m
0
dSm 0
= 2πm w Cm + 2πm C m w 0 − D Sm
0
dτ
D X∞
2m h i
+ (−1) m+n
− 1 Cn0 m≥1 (7.297)
π n=0 m − n
2 2
n6=m

In general, the system given by Eqs. (7.268) to (7.271) can be written as

dx
= f(x) (7.298)
dt
The eigenvalues of the linearized system given by Eqs. (7.294) to (7.297), and in
general form as
dx
= Ax (7.299)
dt
are obtained numerically, such that

|A − λI| = 0 (7.300)

where A is the Jacobian matrix corresponding to the vector field of the linearized
system, I is the identity matrix, and λ are the eigenvalues. The neutral stability
curve is obtained numerically from the condition that <(λ) = 0. A schematic of
the neutral curve is presented in Figure 7.30 for α = 0. In this figure, the stable
and unstable regions can be identified. Along the line of neutral stability, a Hopf-
type of bifurcation occurs. Figure 7.31 shows the plot of w − α curve for a value
of the parameters D = 0.1 and Γ = 0.20029. When α = 0, a Hopf bifurcation
for both the positive and negative branches of the curve can be observed, where
stable and unstable regions can be identified. It is clear that the natural branches
which correspond to the first and third quadrants are stable, whereas the antinatural
branches, second and fourth quadrants are unstable. The symmetry between the
first and third quadrants, and, between the second and fourth quadrants can be

151
notice as well. The corresponding eigenvalues of the bifurcation point are shown in
Figure 7.32. The number of eigenvalues in this figure is 42 which are obtained from a
dynamical system of dimension 42. This system results from truncating the infinite
dimensional system at a number for which the value of the leading eigenvalues does
not change when increasing its dimension. When we increase the size of the system,
new eigenvalues appear in such a way that they are placed symmetrically farther
from the real axis and aligned to the previous set of slave complex eigenmodes. This
behaviour seems to be a characteristic of the dynamical system itself. Figure 7.33
illustrates a view of several stability curves, each for a different value of the tilt angle
α in a D − Γ plane at α = 0. In this plot, the neutral curves appear to unfold when
decreasing the tilt angle from 75.5◦ to −32.5◦ increasing the region of instability.
On the other hand, Figure 7.34 illustrates the linear stability characteristics of the
dynamical system in a w − α plot for a fixed Γ and three values of the parameter
D. The stable and unstable regions can be observed. Hopf bifurcations occur for
each branch of ecah particular curve. However, it is to be notice that the bifurcation
occurs at a higher value of the tilt angle when D is smaller.

7.4.4 Nonlinear analysis

Let us select D = 1.5, and increase Γ. Figures 7.35 and 7.36 show the w − τ time
series and phase-space curves for Γ = 0.95Γcr , whereas Figures 7.37 and 7.38 show
the results for Γ = 1.01Γcr . The plots suggest the appeareance of a subcritical Hopf
bifurcation. Two attractors coexist for Γ ≤ Γcr , these being a critical point and a
strange attractor of fractional dimension. For Γ > Γcr , the only presence is of a
strange attractor.
Now we choose D = 0.1, and increase Γ. Figure 7.39 presents a plot of the w − τ
curve for Γ = 0.99Γcr . Figures 7.40 and 7.41 show the time series plots and phase
space representation for Γ = 1.01Γcr , and Figures 7.42 and 7.43 for Γ = 20Γcr . In
this case, for Γ < Γcr we have stable solutions. For Γ = 1.01Γcr the figures show a
possible limit cycle undergoes a period doubling. This implies a supercritical Hopf
bifurcation. The strange attractor is shown in Figures 7.44 and 7.45.

7.5 Thermal control

Consider the control of temperature at a given point in the loop by modification of
the heating. Both known heat flux and known wall temperatures may be looked at.
In terms of control algorithms, one may use PID or on-off control.

152
 Constant wall
temperature Tw

θ=0, 2π
Cooling
g water out
d=2r ~

θ
o α

R
θ=π

Cooling
water in

Uniform
heat flux q

Figure 7.25: Schematic of a convection loop heated with constant heat flux in one
half and cooled at constant temperature in the other half.

153
 1
D=10
0.8
D=2.5
0.6
0.4
D=1
0.2
D=0.1
-147.5° 32.5°
0
w

-32.5° 147.5°
−0.2
−0.4
−0.6
−0.8
−1

−150 −100 −50 0 50 100 150

α

Figure 7.26: Steady-State velocity field.

154
 3

2.5 D=2.5

2
φ(θ,D,α=0)

1.5 D=1.0

D=0.1
1

0.5

0
0 1 2 3 4 5 6
θ

Figure 7.27: Nondimensional temperature distribution as a function of the parameter

D for α = 0.

155
 6
α= °

4
φ(θ,D=2.5,α)

α= °
3
α= °

0
0 1 2 3 4 5 6 7
θ

Figure 7.28: Nondimensional temperature distribution as a function of thermosyphon

inclination α for D = 2.5.

156
 α=0° α=45°
1
0.8 α=90°

0.6
α=135°
0.4
0.2
0
w

−0.2
−0.4 α=45°

−0.6
−0.8
α=0°
−1

−1 0 1 2
10 10 10 10
D

Figure 7.29: Velocity w as a function of D for different thermosyphon inclinations α.

157
25

20

15
Γ

Unstable
10

Stable
5

0
0 0.5 1 1.5 2 2.5
D

Figure 7.30: Stability curve D and Γ for a thermosyphon inclination α = 0.

158
 1

0.8 Unstable Stable

0.6
Hopf Bifurcation
0.4

0.2

0
w

−0.2

−0.4
Hopf Bifurcation
−0.6
Stable Unstable
−0.8
−1
−150 −100 −50 0 50 100 150
α [°]
Figure 7.31: Stability curve w vs. α for D = 0.1, and Γ = 0.20029.

159
 150

100

50
ℑ(λ)

−50

−100

−150
−0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2
ℜ(λ)
Figure 7.32: Eigenvalues at the neutral curve for D = 0.1, Γ = 0.20029 and α = 0.

160
 25
Unstable

20

15
Γ

10
α=3.5 °

α=21.5 °
5 α=-14.5 °
α=39.5 °
α=57.5 ° α=-32.5 °
α=75.5 ° Stable
0
0 0.5 1 1.5 2 2.5
D
Figure 7.33: Neutral stability curve for different values of the tilt angle α.

161
 1
U
0.8
D=0.1 S
0.6

0.4
D=1.0 D=2.5
0.2
D=1.0
0
w

−0.2
D=0.1
−0.4
D=2.5 U
−0.6

−0.8
S Γ=4.0
−1
−150 −100 −50 0 50 100 150
α°
Figure 7.34: Curve w vs. α Γ = 4.0 and D = 0.1,D = 1.0,D = 2.5.

162
 2.5

1.5

0.5
w

−0.5

−1

−1.5

−2

−2.5
450 460 470 480 490 500 510
τ

Figure 7.35: Curve w vs. τ for D = 1.5, Γ = 0.95Γcr .

163
 2

0
1
φs

−1

−2

−3
4
2 4
0 2
0
−2 −2
φc −4 −4
1 w

Figure 7.36: Phase-space trayectory for D = 1.5, Γ = 0.95Γcr .

164
 2.5

1.5

0.5
w

−0.5

−1

−1.5

−2

−2.5
1900 1920 1940 1960 1980 2000
τ

Figure 7.37: Curve w vs. τ for D = 1.5, Γ = 1.01Γcr .

165
 2

0
1
φs

−1

−2

−3
4
2 4
0 2
0
−2 −2
φc −4 −4
1 w

Figure 7.38: Phase-space trayectory for D = 1.5, Γ = 1.01Γcr .

166
 1.5

1.4

1.3

1.2
w

1.1

0.9

0.8
0 500 1000 1500 2000 2500
τ

Figure 7.39: Curve w vs. τ for D = 0.1, Γ = 0.99Γcr .

167
 2

1.5

0.5
w

−0.5

−1

−1.5

−2
1980 1985 1990 1995 2000
τ

Figure 7.40: Curve w vs. τ for D = 0.1, Γ = 1.1Γcr .

168
 1

0.8

0.6

0.4

0.2
φc1

−0.2

−0.4

−0.6

−0.8

−1
1980 1985 1990 1995 2000
τ

Figure 7.41: Phase-space trayectory for D = 0.1, Γ = 1.1Γcr .

169
 1

0.8

0.6

0.4

0.2
φs1

−0.2

−0.4

−0.6

−0.8

−1
1980 1985 1990 1995 2000
τ

Figure 7.42: Curve w vs. τ for D = 0.1, Γ = 1.1Γcr .

170
 1

0.5

0
φ1
s

−0.5

−1
1
0.5 2
0 1
0
−0.5 −1
φc −1 −2
1 w

Figure 7.43: Phase-space trayectory for D = 0.1, Γ = 1.1Γcr .

171
 0.6

0.4

0.2

0
φc1

−0.2

−0.4

−0.6

−0.8
1088 1090 1092 1094 1096 1098 1100 1102
τ

Figure 7.44: Phase-space trayectory for D = 0.1, Γ = 20Γcr .

172
References
1. Sen, M., Ramos, E. and Treviño, C., On the steady-state velocity of the inclined toroidal ther-
mosyphon, ASME Journal of Heat Transfer, Vol. 107, No. 4, pp. 974–977, 1985.

2. Sen, M., Ramos, E. and Treviño, C., The toroidal thermosyphon with known heat flux, Interna-
tional Journal of Heat and Mass Transfer, Vol. 28, No. 1, pp. 219–233, 1985.

Problems
1. Find the pressure distributions for the different cases of the square loop problem.
2. Consider the same square loop but tilted through an angle θ where 0 ≤ θ < 2π. There is
constant heating between points a and c, and constant cooling between e and g. For the
steady-state problem, determine the temperature distribution and the velocity as a function
of θ. Plot (a) typical temperature distributions for different tilt angles, and (b) the velocity
as a function of tilt angle.
3. Find the steady-state temperature field and velocity for known heating if the loop has a
variable cross-sectional area A(s).
4. Find the temperature field and velocity for known heating if the total heating is not zero.
5. Find the velocity and temperature fields for known heating if the heating and cooling takes
place at two different points. What the condition for the existence of a solution?
6. What is the effect on the known heat rate solution of taking a power-law relationship between
the frictional force and the fluid velocity?
7. For known wall temperature heating, show that if the wall temperature is constant, the
temperature field is uniform and the velocity is zero.
8. Study the steady states of the toroidal loop with known wall temperature including nondi-
mensionalization of the governing equations, axial conduction and tilting effects, multiplicity
of solutions and bifurcation diagrams. Illustrate typical cases with appropriate graphs.
9. Consider a long, thin, vertical tube that is open at both ends. The air in the tube is heated
with an electrical resistance running down the center of the tube. Find the flow rate of the
air due to natural convection. Make any assumptions you need to.

173
 0.2

0.1

0
φ1
s

−0.1

−0.2
0.6
0.4 4
0.2 2
0
0 −2
φc −0.2 −4
1 w

Figure 7.45: Phase-space trayectory for D = 0.1, Γ = 20Γcr .

174
Chapter 8

Convection in porous media

8.1 Governing equations

The continuity equation for incompressible flow in a porous medium is
∇·V = 0 (8.1)

8.1.1 Darcy’s equation

For the momentum equation, the simplest model is that due to Darcy
µ
∇p = − V + ρf f (8.2)
K
where f is the body force per unit mass. Here K is called the permeability of the
medium and has units of inverse area. It is similar to the incompressible Navier-
Stokes equation with constant properties where the inertia terms are dropped and
the viscous force per unit volume is represented by −(µ/K)V. Sometimes a term
cρ0 ∂V/∂t is added to the left side for transient problems, but it is normally left out
because it is very small. The condition on the velocity is that of zero normal velocity
at a boundary, allowing for slip in the tangential direction.
From equations (8.1) and (8.2), for f = 0 we get
∇2 p = 0 (8.3)
from which the pressure distribution can be determined.

8.1.2 Forchheimer’s equation

Forchheimer’s equation which is often used instead of Darcy’s equation is
µ
∇p = − V − cf K −1/2 ρf |V|V + ρf (8.4)
K
175
where cf is a dimensionless constant. There is still slip at a boundary.

8.1.3 Brinkman’s equation

Another alternative is Brinkman’s equation
µ
∇p = − V + µ̃∇2 V + ρf f (8.5)
K
where µ̃ is another viscous coefficient. In this model there is no slip at a solid bound-
ary.

8.1.4 Energy equation

The energy equation is
∂T
(ρc)m + (ρcp )f V · ∇T = km ∇2 T (8.6)
∂t
where km is the effective thermal conductivity, and
(ρc)m = φ(ρcp )f + (1 − φ)(ρc)m (8.7)
is the average heat capacity. Subscripts f and s refer to the fluid and solid respectively,
and φ is the porosity of the material. An equivalent form is
∂T
σ + V · ∇T = αm ∇2 T (8.8)
∂t
where
km
αm = (8.9)
(ρcp )f
(ρc)m
σ = (8.10)
(ρcp )f
See [1].

8.2 Forced convection

8.2.1 Plane wall at constant temperature
The solution to
∂u ∂v
+ = 0 (8.11)
∂x ∂y

176
 K ∂p
u = − (8.12)
µ ∂x
K ∂p
v = − (8.13)
µ ∂y
is
u = U (8.14)
v = 0 (8.15)
For
Ux
= P ex  1 (8.16)
αm
the energy equation is
∂T ∂T ∂2T
u +v = αm 2 (8.17)
∂x ∂y ∂y
or
∂T ∂2T
U = αm 2 (8.18)
∂x ∂y
The boundary conditions are
T (0) = Tw (8.19)
T (∞) = T∞ (8.20)
Writing
s
U
η = y (8.21)
αm x
T − Tw
θ(η) = (8.22)
T∞ − Tw
we get
∂T dθ ∂η
= (T∞ − Tw ) (8.23)
∂x dη ∂x
s !
dθ U 1 −3/2
= (T∞ − Tw ) −y x (8.24)
dη αm 2
∂T dθ ∂η
= (T∞ − Tw ) (8.25)
∂y dη ∂y
s
dθ U
= (T∞ − Tw ) (8.26)
dη αm x
∂2T d2 θ U
= (T∞ − Tw ) 2 (8.27)
∂y 2 dη αm x

177
so that the equation becomes
1
θ00 + η θ0 = 0 (8.28)
2
with

θ(0) = 0 (8.29)
θ(∞) = 1 (8.30)
2 /4
We multiply by the integrating factor eη to get

d  η2 /4 0 
e θ =0 (8.31)
dη

The first integral is

2 /4
θ0 = C1 e−η (8.32)

Integrating again we have

Z η 02 /4
θ = C1 e−η dη 0 + C2 (8.33)
0

With the change in variables x = η 0 /2, the solution becomes

Z η/2 2
θ = 2C1 e−x dx (8.34)
0

√
Applying the boundary conditions, we find that C1 = 1/ π and C2 = 0. Thus
Z η/2
2 2
θ = √ e−x dx (8.35)
π 0
η
= erf (8.36)
2

The heat transfer coefficient is defined as

q 00
h = (8.37)
Tw − T∞
km ∂T
= − (8.38)
Tw − T∞ ∂y
∂θ
= km (8.39)
∂y

178
The local Nusselt number is given by
hx
Nux = (8.40)
km
∂θ
= x (8.41)
∂y y=0
s
1 Ux
= √ (8.42)
π αm
1
= √ P e−1/2 (8.43)
π x

Example 8.1
Find the temperature distribution for flow in a porous medium parallel to a flat plate
with uniform heat flux.

8.2.2 Stagnation-point flow

For flow in a porous medium normal to an infinite flat plate, the velocity field is
u = Cx (8.44)
v = −Cy (8.45)
The energy equation is
∂T ∂T ∂2T
Cx − Cy = αm 2 (8.46)
∂x ∂y ∂y

8.2.3 Thermal wakes

Line source
For P ex  1, the governing equation is
∂T ∂2T
U = αm 2 (8.47)
∂x ∂y
where the boundary conditions are
∂T
= 0 at y = 0 (8.48)
∂y
Z ∞
0
q = (ρcp )f U (T − T∞ ) dy (8.49)
−∞

179
Writing
s
U
η = y (8.50)
αm x
s
T − T∞ Ux
θ(η) = (8.51)
q 0 /km αm
we find that
r   s r !
∂T q0 αm 1 dθ U 1 q0 αm
= θ − x−3/2 + y − x−3/2 (8.52)
∂x km U 2 dη αm 2 km Ux
r s
∂T q 0 αm dθ U
= (8.53)
∂y km Ux dη αm x
r
∂2T q 0 αm dθ U
= (8.54)
∂y 2 km Ux dη αm x
(8.55)
Substituting in the equation, we get
1
θ00 = − (θ + ηθ0 ) (8.56)
2
The conditions (8.48)-(8.49) become
∂θ
= 0 at η = 0 (8.57)
∂η
Z ∞
θ dy = 1 (8.58)
−∞

The equation (8.56) can be written as

1 d
θ00 = − (ηθ) (8.59)
2 dη
which integrates to
1
θ0 = − ηθ + C1 (8.60)
2
0
Since θ = 0 at η = 0, we find that C1 = 0. Integrating again, we get
2 /4
θ = C2 e−η (8.61)
Substituting in the other boundary condition
Z ∞ Z ∞ √
2 /4
1= θ dη = C2 e−η dη = 2 πC2 (8.62)
−∞ −∞

180
from which
1
C2 = √ (8.63)
2 π
Thus the solution is
1 2
θ = √ e−η /4 (8.64)
2 π
or
1 q 0 q αm −Uy 2
T − T∞ = √ ( ) exp( ) (8.65)
2 π km Ux 4αm x

Example 8.2
Show that for a point source

q −U r2
T − T∞ = exp( ) (8.66)
4πkx 4αm x

8.3 Natural convection

8.3.1 Linear stability
This is often called the Horton-Rogers-Lapwood problem, and consists of finding the
stability of a horizontal layer of fluid in a porous medium heated from below. The
geometry is shown in Fig. 8.1.

Figure 8.1: Stability of horizontal porous layer.

The governing equations are

∇·V = 0 (8.67)
µ
−∇p − V + ρf g = 0 (8.68)
K
∂T
σ + V · ∇T = αm ∇2 T (8.69)
∂t
ρf = ρ0 [1 − β (T − T0 )] (8.70)

181
The basic steady solution is

V = 0 (8.71)
 
z
T = T0 + ∆T 1 − (8.72)
"
H !#
1 z2
p = p0 − ρ0 g z + β∆T − 2z (8.73)
2 H

We apply a perturbation to each variable as

V = V + V0 (8.74)
T = T + T0 (8.75)
p = p + p0 (8.76)

Substituting and linearizing

∇ · V0 = 0 (8.77)
µ
−∇p0 − V0 + βρ0 T 0 g = 0 (8.78)
K
∂T 0 ∆T 0
− w = αm ∇2 T 0 (8.79)
∂t H
Using the nondimensional variables
x
x∗ = (8.80)
H
αm t
t∗ = (8.81)
σH 2
HV0
V∗ = (8.82)
αm
T0
T∗ = (8.83)
∆T
Kp0
p∗ = (8.84)
µαm

the equations become, on dropping *s

∇·V = 0 (8.85)
−∇p − V + Ra T k = 0 (8.86)
∂T
− w = ∇2 T (8.87)
∂t
182
where
ρ0 gβKH∆T
Ra = (8.88)
µαm
From these equations we get
∇2 w = Ra∇2H T (8.89)
where
∂2 ∂2
∇H = + (8.90)
∂x2 ∂y 2
Using separation of variables
w(x, y, z, t) = W (z) exp (st + ikx x + iky y) (8.91)
T (x, y, z, t) = Θ(z) exp (st + ikx x + iky y) (8.92)
Substituting into the equations we get
!
d2
− k 2 − s Θ = −W (8.93)
dz 2
!
d2
2
− k W = −k 2 Ra Θ
2
(8.94)
dz
where
k 2 = kx2 + ky2 (8.95)

Isothermal boundary conditions

The boundary conditions are W = Θ = 0 at either wall. For the solutions to remain
bounded as x, y → ∞, the wavenumbers kx and ky must be real. Furthermore, since
the eigenvalue problem is self-adjoint, as shown below, it can be shown that s is also
real.
For a self-adjoint operator L, we must have
(u, Lv) = (Lu, v) (8.96)
If
∂P ∂Q
vL(u) − uL(v) = + (8.97)
∂x ∂y
then
Z Z !
∂P ∂Q
[vL(u) − uL(v)] dV = + dV (8.98)
V ∂x ∂y
Z
= ∇ · (P i + Qj) dV (8.99)
ZV
= n · (P i + Qj) (8.100)
S

183
If n · (P i + Qj) = 0 at the boundaries (i.e. impermeable), which is the case here, then
L is self-adjoint.
Thus, marginal stability occurs when s = 0, for which
!
d2
− k 2 Θ = −W (8.101)
dz 2
!
d2
2
− k W = −k 2 Ra Θ
2
(8.102)
dz

from which !2
d2
− k2 W = k 2 Ra W (8.103)
dz 2
The eigenfunctions are
W = sin nπz (8.104)
where n = 1, 2, 3, . . ., as long as
!2
n2 π 2
Ra = +k (8.105)
k

For each n there is a minimum value of the critical Rayleigh number determined by
!" #
dRa n2 π 2 n2 π 2
=2 +k − 2 +1 (8.106)
dk k k

The lowest critical Ra is with k = π and n = 1, which gives

Rac = 4π 2 (8.107)

for the onset of instability.

Constant heat flux conditions

Here W = dΘ/dz = 0 at the walls. We write

W = W0 + α2 W1 + . . . (8.108)
Θ = Θ0 + α 2 Θ1 + . . . (8.109)
Ra = Ra0 + α2 Ra1 + . . . (8.110)

For the zeroth order system

d2 W0
=0 (8.111)
dz 2
184
 y
x

Figure 8.2: Inclined porous layer.

with W0 = dΘ0 /dz = 0 at the walls. The solutions is W0 = 0, Θ0 = 1. To the next

order
d2 W1
= W0 − Ra0 Θ0 (8.112)
dz 2
d 2 Θ1
+ W1 = Θ0 (8.113)
dz 2
with W1 = dΘ1 /dz = 0 at the walls. The

8.3.2 Steady-state inclined layer solutions

Consider an inclined porous layer of thickness H at angle φ with respect to the
horizontal shown in Fig. 8.2
Introducing the streamfunction ψ(x, y), where
∂ψ
u = (8.114)
∂y
∂ψ
v = − (8.115)
∂x
Darcy’s equation becomes
∂p µ ∂ψ
− − = ρ0 g [1 − β(T − T0 )] sin φ (8.116)
∂x K ∂y

185
 ∂p µ ∂ψ
− + = ρ0 g [1 − β(T − T0 )] cos φ (8.117)
∂y K ∂x
Taking ∂/∂y of the first and ∂/∂x of the second and subtracting, we have
!
∂2ψ ∂2ψ ρ0 gβK ∂T ∂T
+ = − cos φ − sin φ (8.118)
∂x2 ∂y 2 µ ∂x ∂y
The energy equation is
!
∂ψ ∂T ∂ψ ∂T ∂2T ∂2T
αm − = + (8.119)
∂y ∂x ∂x ∂y ∂x2 ∂y 2

Side-wall heating
The non-dimensional equations are
!
∂2ψ ∂2ψ ∂T ∂T
+ 2 = −Ra cos φ − sin φ (8.120)
∂x2 ∂y ∂x ∂y
∂ψ ∂T ∂ψ ∂T ∂2T ∂2T
− = + (8.121)
∂y ∂x ∂x ∂y ∂x2 ∂y 2
where the Rayleigh number is
Ra =? (8.122)
The boundary conditions are
∂T A
ψ = 0, = 0 at x = ± (8.123)
∂x 2
∂T 1
ψ = 0, = −1 at x = ± (8.124)
∂y 2
(8.125)

With a parallel-flow approximation, we assume

ψ = ψ(y) (8.126)
T = Cx + θ(y) (8.127)

The governing equations become

d2 θ dψ
− C = 0 (8.128)
dy 2 dy
∂2ψ dθ
− Ra sin φ + RC cos φ = 0 (8.129)
∂y 2 dy

186
An additional constraint is the heat transported across a transversal section should
be zero. Thus Z 1/2 !
∂T
uT − dy = 0 (8.130)
−1/2 ∂x
Let us look at three cases.
(a) Horizontal layer
For φ = 0◦ the temperature and streamfunction are
" #
RaC 2  2 
T = Cx − y 1 + 4y − 3 (8.131)
24
RaC  
ψ = − 4y 2 − 1 (8.132)
8
Substituting in condition (8.130), we get
 
C 10R − Ra2 C 2 − 120 = 0 (8.133)

the solutions of which are

C = 0 (8.134)
1 q
C = 10(Ra − 12) (8.135)
Ra
1 q
C = − 10(Ra − 12) (8.136)
Ra
The only real solution that exists for Ra ≤ 12 is the conductive solution C = 0. For
C > 12, there are two nonzero values of C which lead to convective solutions, for
which
RaC
ψc = (8.137)
8
12
Nu = (8.138)
12 − RaC 2
For φ = 180◦ , the only real value of C is zero, so that only the conductive solution
exists.
(b) Natural circulation
Let us take C sin φ > 0, for which we get
 
B α
ψc = 1 − cosh (8.139)
C 2
α
Nu = − (8.140)
2B sinh α2 + αC cot φ

187
where

α2 = RC sin φ (8.141)
1 + C cot φ
B = − (8.142)
cosh α2

and the constant C is determined from

!  
B2 sinh α α 2 α
C− − 1 − B cot φ cosh − sinh =0 (8.143)
2C α 2 α 2

(c) Antinatural circulation

For C sin φ < 0, for which we get
!
B β
ψc = 1 − cosh (8.144)
C 2
β
Nu = − β (8.145)
2B sinh 2 + βC cot φ

where

β 2 = −RC sin φ (8.146)

1 + C cot φ
B = − (8.147)
cosh β2

and the constant C is determined from

! !
B2 sin β β 2 β
C− − 1 − B cot φ cosh − sinh =0 (8.148)
2C β 2 β 2

End-wall heating
Darcy’s law is
!
2 ∂T ∂T
∇ψ = R sin φ + cos φ (8.149)
∂x ∂y

The boundary conditions are

∂T A
ψ = 0, = −1 at x = ± (8.150)
∂x 2
∂T 1
ψ = 0, = 0 at x = ± (8.151)
∂y 2

188
 With a parallel-flow approximation, we assume
ψ = ψ(y) (8.152)
T = Cx + θ(y) (8.153)
The governing equations become
dθ dψ
2
−C = 0 (8.154)
dy dy
∂2ψ dθ
2
− R cos φ − RC sin φ = 0 (8.155)
∂y dy
An additional constraint is the heat transported across a transversal section. Thus
Z !
1/2 ∂T
uT − dy = 1 (8.156)
−1/2 ∂x
Let us look at three cases.
(a) Vertical layer
For φ = 0◦ the temperature and streamfunction are
B1 B2
T = Cx + sin(αy) − cos(αy) (8.157)
α α
B1 B2
ψ = cos(αy) + sin(αy) + B3 (8.158)
C C
where
α2 = −RC (8.159)
Substituting in condition (8.130), we get
 
C 10R − R2 C 2 − 120 = 0 (8.160)
the solutions of which are
C = 0 (8.161)
1 q
C = 10(R − 12) (8.162)
R
1 q
C = − 10(R − 12) (8.163)
R
The only real solution that exists for R ≤ 12 is the conductive solution C = 0. For
C > 12, there are two nonzero values of C which lead to convective solutions, for
which
RC
ψc = (8.164)
8
12
Nu = (8.165)
12 − RC

189
 For φ = 180◦ , the only real value of C is zero, so that only the conductive solution
exists.
(b) Natural circulation
Let us take C sin φ > 0, for which we get
 
B α
ψc = 1 − cosh (8.166)
C 2
α
Nu = − α (8.167)
2B sinh 2 + αC cot φ
where
α2 = RC sin φ (8.168)
1 + C cot φ
B = − (8.169)
cosh α2
and the constant C is determined from
!  
B 2 sinh α α 2 α
C− − 1 − B cot φ cosh − sinh =0 (8.170)
2C α 2 α 2
(c) Antinatural circulation
For C sin φ > 0, for which we get
!
B β
ψc = 1 − cosh (8.171)
C 2
β
Nu = − (8.172)
2B sinh β2 + βC cot φ
where
β 2 = −RC sin φ (8.173)
1 + C cot φ
B = − (8.174)
cosh β2
and the constant C is determined from
! !
B 2 sin β β 2 β
C− − 1 − B cot φ cosh − sinh =0 (8.175)
2C β 2 β 2

References
1. Sen, M., Vasseur, P. and Robillard, L., Multiple steady states for unicellular natural convection
in an inclined porous layer, International Journal of Heat and Mass Transfer, Vol. 30, No. 10,
pp. 2097–2113, 1987.

2.

190
Problems
1. This is a problem

191
192
Chapter 9

Multidimensional forced convection

9.1 Low Reynolds numbers

9.2 Potential flow
See [1].

9.3 Multiple solutions

See [2].

9.4 Plate heat exchangers

9.4.1 Cross flow
There is flow on the two sides of a plate, 1 and 2, with an overall heat transfer
coefficient of U.

References
1. Sen, M. and Yang, K.T., Convective heat transfer in two-dimensional potential flows, to be
published.

2. Sen, M. and Vasseur, P., Analysis of multiple solutions in plane Poiseuille flow with viscous heating
and temperature dependent viscosity, Proceedings of the National Heat Transfer Conference,
HTD-Vol. 107, Heat Transfer in Convective Flows, pp. 267–272, 1989.

193
Problems
1. This is a problem

194
Chapter 10

Multi-dimensional natural
convection

10.1 Governing equations

10.2 Cavities
10.3 Marangoni convection
See [1].

References
1. Wang, C.H., Sen, M. and Vasseur, P., Analytical investigation of Bénard-Marangoni convection
heat transfer in a shallow cavity filled with two immiscible fluids, Applied Scientific Research,
Vol. 48, pp. 35–53, 1991.

Problems
1. This is a problem

195
196
Chapter 11

Heat exchangers

11.1 Basic theory

11.1.1 Heat transfer coefficients
Overall heat transfer coefficient
Fouling
Bulk temperature

11.1.2 Nondimensional groups

UL
Reynolds number Re = (11.1)
ν
ν
Prandtl number = (11.2)
κ
hL
Nusselt number Nu = (11.3)
k
Stanton number St = Nu/P r Re (11.4)
Colburn j-factor j = St P r 2/3 (11.5)
2τw
Friction factor f = (11.6)
ρU 2

11.1.3 Duct flow

11.1.4 Parallel flow and counterflow
We define the subscripts h and c to mean hot and cold fluids, i and o for inlet and
outlet, 1 the end where the hot fluids enters, and 2 the other end. Energy balances

197
 cold
(a) parallel flow
hot

cold
(a) counter flow
hot

Figure 11.1: Parallel and counter flow.

give
dq = U(Th − Tc ) dA (11.7)
dq = ṁc Cc dTc (11.8)
dq = −ṁh Ch dTh (11.9)
From equations (11.8) and (11.9), we get
 
1 1
−dq + = d(Th − Tc ) (11.10)
ṁh Ch ṁc Cc
Using (11.7), we find that
 
1 1 d(Th − Tc )
−U dA + = (11.11)
ṁh Ch ṁc Cc Th − Tc
which can be integrated from 1 to 2 to give
 
1 1 (Th − Tc )1
−UA + = ln (11.12)
ṁh Ch ṁc Cc (Th − Tc )2
From equation (11.10), we get
 
1 1
−q + = (Th − Tc )2 − (Th − Tc )1 (11.13)
ṁh Ch ṁc Cc

198
The last two equations can be combined to give

q = UA∆Tlmtd (11.14)

where
(Th − Tc )1 − (Th − Tc )2
∆Tlmtd = (11.15)
ln[(Th − Tc )1 /(Th − Tc )2 ]
is the logarithmic mean temperature difference.
For parallel flow, we have
(Th,i − Tc,i) − (Th,o − Tc,o )
∆lmtd = (11.16)
ln[(Th,i − Tc,i)/(Th,o − Tc,o )]
while for counterflow it is
(Th,i − Tc,o) − (Th,o − Tc,i )
∆lmtd = (11.17)
ln[(Th,i − Tc,o)/(Th,o − Tc,i )]

11.1.5 Crossflow plate heat exchanger

Consider a rectangular plate of size Lx × Ly in the x- and y-directions, respectively,
as shown in Fig. 11.2. The flow on one side of the plate is in the x-direction with a
temperature field Tx (x, y). The mass flow rate of the flow is mx per unit transverse
length. The flow in the other side of the plate is in the y-direction with the corre-
sponding quantities Ty (x, y) and my . The overall heat transfer coefficient between
the two fluids is U, which we will take to be a constant.
For the flow in the x-direction, the steady heat balance on an elemental rectangle
of size dx × dy gives
∂Tx
cx mx dy dx = U dx dy (Ty − Tx ) (11.18)
∂x
where cx is the specific heat of that fluid. Simplifying, we get
∂Tx
2Cx R = Ty − Tx (11.19)
∂x
where R = 1/2U is proportional to the thermal resistance between the two fluids,
and Cx = cx mx . For the other fluid
∂Ty
2Cy R = Tx − Ty (11.20)
∂y
These equations have to be solved with suitable boundary conditions to obtain the
temperature fields Tx (x, y) and Ty (x, y).

199
 x flow

y flow

Figure 11.2: Schematic of crossflow plate HX.

From equation (11.20), we get

∂Ty
Tx = Ty + Cy R (11.21)
∂y

Substituting in equation (11.19), we have

2
1 ∂Ty 1 ∂Ty Ty
+ + 2R =0 (11.22)
Cy ∂x Cx ∂y ∂x∂y

Nusselt (Jakob, 1957) gives an interesting solution in the following manner. Let
the plate be of dimensions L and W in the x- and y-directions. Nondimensional
variables are
x
ξ = (11.23)
L
y
η = (11.24)
W
Tx − Ty,i
θx = (11.25)
Tx,i − Ty,i
Ty − Ty,i
θy = (11.26)
Tx,i − Ty,i
UW L
a = (11.27)
Cx
UW L
b = (11.28)
Cy

200
 The governing equations are then

∂θx
a(θx − θy ) = − (11.29)
∂ξ
∂θy
b(θx − θy ) = (11.30)
∂η

with boundary conditions

θx = 1 at ξ = 0 (11.31)
θy = 0 at η = 0 (11.32)

Equation (11.30) can be written as

∂θy
+ bθy = bθx (11.33)
∂η

Solving for θy we get

 Z η 
0
θy = e−bη C(ξ) + b θx (ξ, η 0)ebη dη 0 (11.34)
0

From the boundary condition (11.32), we get

C=0 (11.35)

so that Z η 0
θy (ξ, η) = be−bη θx (ξ, η 0)ebη dη 0 (11.36)
0

Using the same procedure, from equation (11.29) we get

Z ξ
−aξ −aξ 0
θx (ξ, η) = e + ae θy (ξ 0, η)eaξ dξ 0 (11.37)
0

Substituting for θy , we find the Volterra integral equation

Z ξ Z η
−aξ −(aξ+bη) 0 0
θx (ξ, η) = e + abe θx (ξ 0 , η 0)eaξ +bη dξ 0 dη 0 (11.38)
0 0

for the unknown θx .

We will first solve the Volterra equation for an arbitrary λ, where
Z ξ Z η
−aξ −(aξ+bη) 0 0
θx (ξ, η) = e + abλe θx (ξ 0 , η 0 )eaξ +bη dξ 0 dη 0 (11.39)
0 0

201
Let us express the solution in terms of a finite power series

θx (ξ, η) = φ0 (ξ, η) + λφ1 (ξ, η) + λ2 φ2 (ξ, η) + . . . + λn φn (ξ, η) (11.40)

This can be substituted in the integral equation. Since λ is arbitrary, the coefficient
of each order of λ must vanish. Thus

φ0 (ξ, η) = e−aξ (11.41)

Z ξ Z η 0 0
φ1 (ξ, η) = abe−(aξ+bη) φ0 (ξ 0 , η 0 )eaξ +bη dξ 0 dη 0 (11.42)
0 0
Z ξZ η 0 0
φ2 (ξ, η) = abe−(aξ+bη) φ1 (ξ 0 , η 0 )eaξ +bη dξ 0 dη 0 (11.43)
0 0
..
. (11.44)
Z ξ Z η 0 0
φn (ξ, η) = abe−(aξ+bη) φn−1 (ξ 0 , η 0)eaξ +bη dξ 0 dη 0 (11.45)
0 0
(11.46)

The solutions are

φ0 = e−aξ (11.47)
φ1 = aξe−aξ (1 − e−bη ) (11.48)
1 2 2 −aξ
φ2 = a ξ e (1 − e−bη − bηe−bη ) (11.49)
2
1 3 3 −aξ 1
φ3 = a ξ e (1 − e−bη − bηe−bη − b2 η 2 e−bη ) (11.50)
2×3 2
..
. (11.51)
1 n n −aξ 1
φn = a ξ e (1 − e−bη − bηe−bη − . . . − bn−1 η n−1 e−bη ) (11.52)
n! (n − 1)!

Substituting into the expansion, equation (11.40), and taking λ = 1, we get

θx (ξ, η) = φ0 (ξ, η) + φ1 (ξ, η) + φ2 (ξ, η) + . . . + φn (ξ, η) (11.53)

where the φs are given above.

Example 11.1
Find a solution of the same problem by separation of variables.
Taking
Ty (x, y) = X(x)Y (y) (11.54)

202
 Substituting and dividing by XY , we get
1 ∂Tx 1 dY
+ + 2R = 0 (11.55)
Cx ∂x Cy dy
Since the first term is a function only of x, and the second only of y, each must be a constant.
Thus we can write
dX 1
+ X = 0 (11.56)
dx cx mx (a + R)
dY 1
+ Y = 0 (11.57)
dy cy my (a − R)
where a is a constant. Solving the two equations and taking their product, we have
 
c x y
Ty = exp − + (11.58)
a+R cx mx (a + R) cy my (a − R)

where c is a constant. Substituting in equation (11.21), we get

 
c x y
Tx = exp − + (11.59)
a−R cx mx (a + R) cy my (a − R)
The rate of heat transfer over the entire plate, Q, is given by
Z L−y
Q = Cx [Tx (Lx , y) − Tx (0, y)] dy (11.60)
0
 
lx Ly
= cCx Cy exp − − −2 (11.61)
Cx (a + R) C2 (a − R)
The heat rate can be maximized by varying either of the variables Cx or Cy .

11.2 HX equation
11.2.1 Definitions
The HX effectiveness is

Q
 = (11.62)
Qmax
Ch (Th,i − Th,o )
= (11.63)
Cmin (Th,i − Tc,i )
Cc (Tc,o − Tc,i)
= (11.64)
Cmin (Th,i − Tc,i )

203
where
Cmin = min(Ch , Cc ) (11.65)
Assuming U to be a constant, the number of transfer units is
AU
NT U = (11.66)
Cmin
The capacity ratio is CR = Cmin /Cmax .

11.2.2 Effectiveness-N T U relations

In general, the effectiveness is a function of the HX configuration, its NT U and the
CR of the fluids.
(a) Counterflow
1 − exp[−NT U(1 − CR )]
= (11.67)
1 − CR exp[−NT U(1 − CR )]
so that  → 1 as NT U → ∞.
(b) Parallel flow
1 − exp[−NT U(1 − CR )]
= (11.68)
1 + CR
(c) Crossflow, both fluids unmixed
Series solution (Mason, 1954)
(d) Crossflow, one fluid mixed, the other unmixed
If the unmixed fluid has C = Cmin , then
 = 1 − exp[−CR (1 − exp{−NT UCr })] (11.69)
But if the mixed fluid has C = Cmin
 = CR (1 − exp{−CR (1 − e−N T U )}) (11.70)
(e) Crossflow, both fluids mixed
(f) Tube with wall temperature constant

 = 1 − exp(−NT U) (11.71)

11.3 Design methodology

11.3.1 Mean temperature-difference method
Given the inlet temperatures and flow rates, this method enables one to find the
outlet temperatures, the mean temperature difference, and then the heat rate.

204
11.3.2 Effectiveness-NTU method
The order of calculation is NT U, , qmax and q.

11.4 Pressure drop

It is important to determine the pressure drop through a heat exchanger. This is
given by
" #
∆p G2 ρ1 Aρ1 ρ1
= (Kc + 1 − σ 2 ) + 2( − 1) + f − (1 − σ 2 − Ke ) (11.72)
p1 2ρ1 p1 ρ2 Ac ρm ρ2

where Kc and Ke are the entrance and exit loss coefficients, and σ is the ratio of
free-flow area to frontal area.

11.5 Correlations
11.6 Extended surfaces
Af
η0 = 1 − (1 − ηf ) (11.73)
A
where η0 is the total surface temperature effectiveness, ηf is the fin temperature
effectiveness, Af is the HX total fin area, and A is the HX total heat transfer area.

11.6.1 Fin analysis

Analysis of Kraus (1990) for variable heat transfer coefficients.

11.7 Porous medium analogy

See Nield and Bejan, p. 87.

11.8 Heat transfer augmentation

11.9 Maldistribution effects
Rohsenow (1981)

205
11.10 Microchannel heat exchangers
Phillips (1990).

11.11 Radiation effects

See Ozisik (1981).
Shah (1981)

11.12 Transient behavior

Ontko and Harris (1990)
For both fluids mixed
dT
Mc + ṁ1 c1 (T1in − T1out ) + ṁc2 (T2in − T2out ) = 0 (11.74)
dt
For one fluid mixed and the other unmixed, we have
∂T2 ∂T2
ρAc2 + ρV2 Ac2 + hP (T − T1 ) = 0 (11.75)
∂t ∂x

References

1. Shah, R.K. and London, A.L., 1978, Laminar flow forced convection in ducts, a source book for
compact heat exchanger analytical data, New York, Academic Press, 1978.

2. Sen, M. and Vasseur, P., Analysis of multiple solutions in plane Poiseuille flow with viscous heating
and temperature dependent viscosity, Proceedings of the National Heat Transfer Conference,
HTD-Vol. 107, Heat Transfer in Convective Flows, pp. 267–272, 1989.

Problems
1. This is a problem

206
Chapter 12

Heat transfer correlations

12.1 Least squares method

Possible correlations are
X
n
y(x) = ak xk (12.1)
k=0
y(x) = a0 + a1/2 x1/2 + a1 x + a3/2 x3/2 + a2 x2 + . . . (12.2)
y(x) = axm + bxn + . . . (12.3)
Z n
y(x) = a(k)xk dk where − 1 < x < 1 (12.4)
k=0

A power-law correlation of the form

y = cxn (12.5)

satisfies the invariance condition given by equation (1.1).

12.2 Genetic algorithms

See [1].
Evolutionary programming, of which genetic algorithms and programming are
examples, allow programs to change or evolve as they compute. GAs, specifically,
are based on the principle of Darwinian selection. One of their most important
applications in the thermal sciences is in the area of optimization of various kinds.
Optimization by itself is fundamental to many applications. In engineering, for
example, it is important to the design of systems; analysis permits the prediction of
the behavior of a given system, but optimization is the technique that searches among

207
all possible designs of the system to find the one that is the best for the application.
The importance of this problem has given rise to a wide variety of techniques which
help search for the optimum. There are searches that are gradient-based and those
that are not. In the former the search for the optimum solution, as for example
the maximum of a function of many variables, starts from some point and directs
itself in an incremental fashion towards the optimum; at each stage the gradient of
the function surface determines the direction of the search. Local optima can be
found in this way, the search for global optimum being more difficult. Again, if one
visualizes a multi-variable function, it can have many peaks, any one of which can be
approached by a hill-climbing algorithm. To find the highest of these peaks, the entire
domain has to be searched; the narrower this peak the finer the searching “comb”
must be. For many applications this brute-force approach is too expensive in terms
of computational time. Alternatives, like simulated annealing, are techniques that
have been proposed, and the GA is one of them.
In what follows we will provide an overview of the genetic algorithm and pro-
gramming. A numerical example will be explained in some detail. The methodology
will be applied to one of the heat exchangers discussed before. There will a discussion
on other applications in thermal engineering and comments will be made on potential
uses in the future.

12.2.1 Methodology
GAs are discussed in detail by Holland (1975, 1992), Mitchell (1997), Goldberg (1989),
Michalewicz, (1992) and Chipperfield (1997). One of the principal advantages of this
method is its ability to pick out a global extremum in a problem with multiple local
extrema. For example, we can discuss finding the maximum of a function f (x) in a
given domain a ≤ x ≤ b. In outline the steps of the procedure are the following.

• First, an initial population of n members x1 , x2 , . . . , xn ∈ [a, b] is randomly

generated.

• Then, for each x a fitness is evaluated. The fitness or effectiveness is the pa-
rameter that determines how good the current x is in terms of being close to an
optimum. Clearly, in this case the fitness is the function f (x) itself, since the
higher the value of f (x) the closer we are to the maximum.

• The probability distribution for the next generation is found based on the fitness
values of each member of the population. Pairs of parents are then selected on
the basis of this distribution.

208
 • The offsprings of these parents are found by crossover and mutation. In crossover
two numbers in binary representation, for example, produce two others by inter-
changing part of their bits. After this, and based on a preselected probability,
some bits are randomly changed from 0 to 1 or vice versa. Crossover and mu-
tation create a new generation with a population that is more likely to be fitter
than the previous generation.
• The process is continued as long as desired or until the largest fitness in a
generation does not change much any more.
The procedure can be easily generalized to a function of many variables.
Let us consider a numerical example that is shown in detail in Table 12.1. Sup-
pose that one has to find the x at which f (x) = x(1 − x) is globally a maximum
between 0 and 1. We have taken n = 6, meaning that each generation will have six
numbers. Thus, for a start 6 random numbers are selected between 0 and 1. Now
we choose nb which is the number of bits used to represent a number in binary form.
Taking nb = 5, we can write the numbers in binary form normalized between 0 and
the largest number possible for nb bits, which is 2nb − 1 = 31. In one run the numbers
chosen, and written down in the first column of the table labeled G = 0, are 25, 30,
28, 19, 3, and 1, respectively. The fitnesses of each one of the numbers, i.e. f (x), are
computed and shown in column two. These values are normalized by their sum and
shown in the third column as s(x). The normalized fitnesses are drawn on a roulette
wheel in Figure 12.1. The probability of crossover is taken to be 100%, meaning
that crossover will always occur. Pairs of numbers are chosen by spinning the wheel,
the numbers having a bigger piece of the wheel having a larger probability of being
selected. This produces column four marked G = 1/4, and shuffling to producing
random pairing gives column five marked G = 1/2. The numbers are now split up
in pairs, and crossover applied to each pair. The first pair [0 0 0 1 1] and [1 1 1 0
0] produces [0 0 0 1 0] and [1 1 1 0 1]. This is illustrated in Figure 12.2(a) where
the crossover position is between the fourth and fifth bit; the bits to the right of this
line are interchanged. Crossover positions in the other pairs are randomly selected.
Crossover produces column six marked as G = 3/4. Finally, one of the numbers, in
this case the last number in the list [0 0 1 1 0], is mutated to [0 0 1 0 0] by changing
one randomly selected bit from 1 to 0 as shown in Figure 12.2(b). From the numbers
in generation G = 0, these steps have now produced a new generation G = 1. The
process is repeated until the largest fitness in each generation increases no more. In
this particular case, values within 3.22% of the exact value of x for maximum f (x),
which is the best that can be done using 5 bits, were usually obtained within 10
generations.
The genetic programming technique (Koza, 1992; Koza, 1994) is an extension
of this procedure in which computer codes take the place of numbers. It can be

209
G=0 f(x) s(x) G = 1/4 G = 1/2 G = 3/4 G = 1
11001 0.1561 0.2475 00011 00011 00010 00010
11110 0.0312 0.0495 00011 11100 11101 11101
11100 0.0874 0.1386 11110 00011 10011 10011
10011 0.2373 0.3762 10011 10011 00011 00011
00011 0.0874 0.1386 00011 11110 11011 11011
00001 0.0312 0.0495 11100 00011 00110 00100

Table 12.1: Example of use of the genetic algorithm.

4.95%
24.75%
13.86%

4.95%

13.86% 37.62%

Figure 12.1: Distribution of fitnesses.

210
 0 0 0 1 1 1 1 1 0 0

(a)

0 0 0 1 0 1 1 1 0 1

0 0 1 1 0
(b)

0 0 1 0 0

Figure 12.2: (a) Crossover and (b) mutation in a genetic algorithm.

used in symbolic regression to search within a set of functions for the one which
best fits experimental data. The procedure is similar to that for the GA, except
for the crossover operation. If each function is represented in tree form, though not
necessarily of the same length, crossover can be achieved by cutting and grafting.
As an example, Figure 12.3 shows the result of the operation on the two functions
3x(x + 1) and x(3x + 1) to give 3x(3x + 1) and x(x + 1). The crossover points may
be different for each parent.

12.2.2 Applications to compact heat exchangers

The following analysis is on the basis of data collected on a single-row heat exchanger
referred to as heat exchanger 1 in Section 2.2. In the following a set of N = 214
experimental runs provided the data base. The heat rate is determined by
Q̇ = ṁa cp,a (Taout − Tain ) (12.6)
= ṁw cw (Twin − Twout ) (12.7)
For prediction purposes we will use functions of the type
Q̇ = q̇(Twin , Tain , ṁa , ṁw ) (12.8)
The conventional way of correlating data is to determine correlations for inner and
outer heat transfer coefficients. For example, power laws of the following form
1/3
εNua = a Rem
a P ra (12.9)

211
 * *

* + + x
Parents
3 x xx 1 * 1

3 x

* *

* + + x

Offspring
x3 x * 1 xx 1

3 x

Figure 12.3: Crossover in genetic programming. Parents are 3x(x + 1) and x(3x + 1);
offspring are 3x(3x + 1) and x(x + 1).

Nuw = b Renw P rw0.3 (12.10)

are common. The two Nusselt numbers provide the heat transfer coefficients on each
side and the overall heat transfer coefficient, U, is related to ha and hw by

1 1 1
= + (12.11)
UAa hw Aw εha Aa

To find the constants a, b, m, n, the mean square error

 
1 X 1 1 2
SU = − (12.12)
N Up Ue
must be minimized, where N is the number of experimental data sets, U p is the
prediction made by the power-law correlation, and U e is the experimental value for
that run. The sum is over all N runs.
This procedure was carried out for the data collected. It was found that the SU
had local minima for many different sets of the constants, the following two being
examples.

212
 Correlation a b m n
A 0.1018 0.0299 0.591 0.787
B 0.0910 0.0916 0.626 0.631

Figure 12.4 shows a section of the SU surface that passes though the two minima
A and B. The coordinate z is a linear combination of the constants a, b, m and n
such that it is zero and unity at the two minima. Though the values of SU for the
two correlations are very similar and the heat rate predictions for the two correlations
are also almost equally accurate, the predictions on the thermal resistances on either
side are different. Figure 12.5 shows the ratio of the predicted air- and water-side
Nusselt numbers using these two correlations. Ra is the ratio of the Nusselt number
on the air side predicted by Correlation A divided by that predicted by Correlation
B. Rw is the same value for the water side. The predictions, particularly the one on
the water side, are very different.
There are several reasons for this multiplicity of minima of SU . Experimentally,
it is very difficult to measure the temperature at the wall separating the two fluids,
or even to specify where it should be measured, and mathematically, it is due to the
nonlinearity of the function to be minimized. This raises the question as to which of
the local minima is the “correct” one. A possible conclusion is that the one which
gives the smallest value of the function should be used. This leads to the search for
the global minimum which can be done using the GA.
For this data, Pacheco-Vega et al. (1998) conducted a global search among a
proposed set of heat transfer correlations using the GA. The experimentally deter-
mined heat rate of the heat exchanger was correlated with the flow rates and input
temperatures, with all values being normalized. To reduce the number of possibilities
the total thermal resistance was correlated with the mass flow rates in the form
Twin − Tain
= f (ṁa , ṁw ) (12.13)
Q̇
The functions f (ṁa , ṁw ) that were used are indicated in Table 12.2. The GA was
used to seek the values of the constants associated with each correlation, the objective
being to minimize the variance
1 X p 2
SQ = Q̇ − Q̇e (12.14)
N
where the sum is over all N runs, between the predictions of a correlation, Q̇p , and
the actual experimental values, Q̇e . Since the unknowns are the set of constants a, b,
c and sometimes d, a single binary string represents them; the first part of the string
is a, the next is b, and so on. The rest of the GA is as in the numerical example

213
 −5
3 x 10 THIS FIGURE WILL BE PASTED IN
SU (m2K/W)2

2.5

1.5

0.5 A B

0
−0.2 0 0.2 0.4 0.6 0.8 1 1.2
z
Figure 12.4: Section of SU (a, b, m, n) surface.

214
 Rea 2
0 2 4 6 8 10 x 10
1.5
1.4
1.3
1.2
Ra
1.1
1
0.9
0.8
0.7 Rw
0.6
0.5 4
0 1 2 3 4 5 x 10
Rew

Figure 12.5: Ratio of the predicted air- and water-side Nusselt numbers.

215
 Correlation f a b c d σ

Power aṁ−b −d
w + cṁa 0.1875 0.9997 0.5722 0.5847 0.0252
law
Inverse (a + bṁw )−1 −0.0171 5.3946 0.4414 1.3666 0.0326
linear +(c + dṁa )−1
Inverse (a + ebṁw )−1 −0.9276 3.8522 −0.4476 0.6097 0.0575
exponential +(c + edṁa )−1
Exponential ae−bṁw + ce−dṁa 3.4367 6.8201 1.7347 0.8398 0.0894

Inverse (a + bṁ2w )−1 0.2891 20.3781 0.7159 0.7578 0.0859

quadratic +(c + dṁ2a )−1
Inverse (a + b ln ṁw )−1 0.4050 0.0625 −0.5603 0.2048 0.1165
logarithmic +(c + d ln ṁa )−1
Logarithmic a − b ln ṁw 0.6875 0.4714 0.4902 − 0.1664
−c ln ṁa
Linear a − bṁw − cṁa 2.3087 0.8533 0.8218 − 0.2118

Quadratic a − bṁ2w − cṁ2a 1.8229 0.6156 0.5937 − 0.2468

Table 12.2: Comparison of best fits for different correlations.

given before. The results obtained for each correlation are also summarized in the
table in descending order of SQ . The last column shows the mean square error σ
defined in a manner similar to equations (12.24)-(12.25). The parameters used for
the computations are: population size 20, number of generations 1000, bits for each
variable 30, probability of crossover 1, and probability of mutation 0.03.
Some correlations are clearly seen to be superior to others. However, the differ-
ence in SQ between the first- and second-place correlations, the power-law and inverse
logarithmic which have mean errors of 2.5% and 3.3% respectively, is only about 8%,
indicating that either could do just as well in predictions even though their func-
tional forms are very different. In fact, the mean error in many of the correlations
is quite acceptable. Figures 12.6 shows the predictions of the power-law correlation
versus the experimental values, all in normalized variables. The prediction is seen
to be very good. The quadratic correlation, on the other hand, is the worst in the
set of correlations considered, and Figure 12.7 shows its predictions. It must also be
remarked that, because of the random numbers used in the procedure, the computer

216
 1.4

1.2

+10%
1

0.8
⋅
Qp

−10%
0.6

0.4

0.2

0
0 0.2 0.4 0.6 0.8 1
⋅
Qe

Figure 12.6: Experimental vs. predicted normalized heat flow rates for a power-
law correlation. The straight line is the line of equality between prediction and
experiment, and the broken lines are ±10%.

program gives slightly different results each time it is run, changing the lineup of the
less appropriate correlations somewhat.

12.2.3 Additional applications in thermal engineering

Though the GA is a relatively new technique in relation to its application to thermal
engineering, there are a number of different applications that have already been suc-
cessful. Davalos and Rubinsky (1996) adopted an evolutionary-genetic approach for
numerical heat-transfer computations. Shape optimization is another area that has
been developed. Fabbri (1997) used a GA to determine the optimum shape of a fin.
The two-dimensional temperature distribution for a given fin shape was found using
a finite-element method. The fin shape was proposed as a polynomial, the coefficients

217
 1.4

1.2

+10%
1

0.8
⋅
Qp

−10%
0.6

0.4

0.2

0
0 0.2 0.4 ⋅e
Q
0.6 0.8 1

Figure 12.7: Experimental vs. predicted normalized heat flow rates for a quadratic
correlation. The straight line is the line of equality between prediction and experi-
ment, and the broken lines are ±10%.

218
of which have to be calculated. The fin was optimized for polynomials of degree 1
through 5. Von Wolfersdorf et al. (1997) did shape optimization of cooling channels
using GAs. The design procedure is inherently an optimization process. Androulakis
and Venkatasubramanian (1991) developed a methodology for design and optimiza-
tion that was applied to heat exchanger networks; the proposed algorithm was able
to locate solutions where gradient-based methods failed. Abdel-Magid and Dawoud
(1995) optimized the parameters of an integral and a proportional-plus-integral con-
troller of a reheat thermal system with GAs. The fact that the GA can be used to
optimize in the presence of variables that take on discrete values was put to advan-
tage by Schmit et al. (1996) who used it for the design of a compact high intensity
cooler. The placing of electronic components as heat sources is a problem that has
become very important recently from the point of view of computers. Queipo et al.
(1994) applied GAs to the optimized cooling of electronic components. Tang and
Carothers (1996) showed that the GA worked better than some other methods for
the optimum placement of chips. Queipo and Gil (1997) worked on the multiob-
jective optimization of component placement and presented a solution methodology
for the collocation of convectively and conductively air-cooled electronic components
on planar printed wiring boards. Meysenc et al. (1997) studied the optimization of
microchannels for the cooling of high-power transistors. Inverse problems may also
involve the optimization of the solution. Allred and Kelly (1992) modified the GA
for extracting thermal profiles from infrared image data which can be useful for the
detection of malfunctioning electronic components. Jones et al. (1995) used thermal
tomographic methods for the detection of inhomogeneities in materials by finding lo-
cal variations in the thermal conductivity. Raudensky et al. (1995) used the GA in the
solution of inverse heat conduction problems. Okamoto et al. (1996) reconstructed a
three-dimensional density distribution from limited projection images with the GA.
Wood (1996) studied an inverse thermal field problem based on noisy measurements
and compared a GA and the sequential function specification method. Li and Yang
(1997) used a GA for inverse radiation problems. Castrogiovanni and Sforza (1996,
1997) studied high heat flux flow boiling systems using a numerical method in which
the boiling-induced turbulent eddy diffusivity term was used with an adaptive GA
closure scheme to predict the partial nucleate boiling regime.
Applications involving genetic programming are rarer. Lee et al. (1997) studied
the problem of correlating the CHF for upward water flow in vertical round tubes
under low pressure and low-flow conditions. Two sets of independent parameters
were tested. Both sets included the tube diameter, fluid pressure and mass flux. The
inlet condition type had, in addition, the heated length and the subcooling enthalpy;
the local condition type had the critical quality. Genetic programming was used as a
symbolic regression tool. The parameters were non-dimensionalized; logarithms were
taken of the parameters that were very small. The fitness function was defined as

219
the mean square difference between the predicted and experimental values. The four
arithmetical operations addition, subtraction, multiplication and division were used
to generate the proposed correlations. The programs ran up to 50 generations and
produced 20 populations in each generation. In a first intent, 90% of the data sets was
randomly selected for training and the rest for testing. Since no significant difference
was found in the error for each of the sets, the entire data set was finally used both
for training and testing. The final correlations that were found had predictions better
than those in the literature. The advantage of the genetic programming method in
seeking an optimum functional form was exploited in this application.

12.2.4 General discussion

The evolutionary programming method has the advantage that, unlike the ANN, a
functional form of the relationship is obtained. Genetic algorithms, genetic program-
ming and symbolic regression are relatively new techniques from the perspective of
thermal engineering, and we can only expect the applications to grow. There are a
number of areas in prediction, control and design that these techniques can be effec-
tively used. One of these, in which progress can be expected, is in thermal-hydronic
networks. Networks are complex systems built up from a large number of simple
components; though the behavior of each component may be well understood, the
behavior of the network requires massive computations that may not be practical.
Optimization of networks is an important issue from the perspective of design, since
it is not obvious what the most energy-efficient network, given certain constraints,
should be. The constraints are usually in the form of the locations that must be
served and the range of thermal loads that are needed at each position. A search
methodology based on the calculation of every possible network configuration would
be very expensive in terms of computational time. An alternative based on evolution-
ary techniques would be much more practical. Under this procedure a set of networks
that satisfy the constraints would be proposed as candidates for the optimum. From
this set a new and more fit generation would evolve and the process repeated until
the design does not change much. The definition of fitness, for this purpose, would
be based on the energy requirements of the network.

12.3 Artificial neural networks

See [2].

220
12.4 Artificial neural networks
In this section we will discuss the ANN technique, which is generally considered to
be a sub-class of AI, and its application to the analysis of complex thermal systems.
Applications of ANNs have been found in such diverse fields as philosophy, psychology,
business and economics, sociology, science, a well as in engineering. The common
denominator is the complexity of the field.
The technique is rooted in and inspired by the biological network of neurons in
the human brain that learns from external experience, handles imprecise information,
stores the essential characteristics of the external input, and generalizes previous
experience (Eeckman, 1992). In the biological network of interconnecting neurons,
each receives many input signals from other neurons and gives only one output signal
which is sent to other neurons as part of their inputs. If the sum of the inputs to a
given neuron exceeds a set threshold, normally determined by the electric potential of
the receiver neuron which may be modified under different circumstances, the neuron
fires and sends a signal to all the connected receiver neurons. If not, the signal is
not transmitted. The firing decision represents the key to the learning and memory
ability of the neural network.
The ANN attempts to mimic the biological neural network: the processing unit
is the artificial neuron; it has synapses or inter-neuron connections characterized by
synaptic weights; an operator performs a summation of the input signals weighted by
the respective synapses; an activation function limits the permissible amplitude range
of the output signal. It is also important to realize the essential difference between a
biological neural network and an ANN. Biological neurons function much slower than
the computer calculations associated with an artificial neuron in an ANN. On the
other hand, the delivery of information across the biological neural network is much
faster. The biological one compensates for the relatively slow chemical reactions in
a neuron by having an enormous number of interconnected neurons doing massively
parallel processing, while the number of artificial neurons must necessarily be limited
by the available hardware.
In this section we will briefly discuss the basic principles and characteristics of
the multilayer ANN, along with the details of the computations made in the feedfor-
ward mode and the associated backpropagation algorithm which is used for training.
Issues related to the actual implementation of the algorithm will also be noted and
discussed. Specific examples on the performance of two different compact heat ex-
changers analyzed by the ANN approach will then be shown, followed by a discussion
on how the technique can also be applied to the dynamic performance of heat ex-
changers as well as to their control in real thermal systems. Finally, the potential of
applying similar ANN techniques to other thermal-system problems and their specific
advantages will be delineated.

221
12.4.1 Methodology
The interested reader is referred to the text by Haykin (1994) for an account of
the history of ANN and its mathematical background. Many different definitions of
ANNs are possible; the one proposed by Schalkoff (1997) is that an ANN is a network
composed of a number of artificial neurons. Each neuron has an input/output char-
acteristic and implements a local computation or function. The output of any neuron
is determined by this function, its interconnection with other neurons, and external
inputs. The network usually develops an overall functionality through one or more
forms of training; this is the learning process. Many different network structures and
configurations have been proposed, along with their own methodologies of training
(Warwick et al., 1992).

Feedforward network
There are many different types of ANNs, but one of the most appropriate for engi-
neering applications is the supervised fully-connected multilayer configuration (Zeng,
1998) in which learning is accomplished by comparing the output of the network with
the data used for training. The feedforward or multilayer perceptron is the only con-
figuration that will be described in some detail here. Figure 12.8 shows such an ANN
consisting of a series of layers, each with a number of nodes. The first and last layers
are for input and output, respectively, while the others are the hidden layers. The
network is said to be fully-connected when any node in a given layer is connected to
all the nodes in the adjacent layers.
We introduce the following notation: (i, j) is the jth node in the ith layer. The
line connecting a node (i, j) to another node in the next layer i + 1 represents the
synapse between the two nodes. xi,j is the input of the node (i, j), yi,j is its output,
i,j
θi,j is its bias, and wi−1,k is the synaptic weight between nodes (i−1, k) and (i, j). The
total number of layers, including those for input and output, is I, and the number of
nodes in the ith layer is Ji . The input information is propagated forward through the
network; J1 values enter the network and JI leave. The flow of information through
the layers is a function of the computational processing occurring at every internal
node in the network. The relation between the output of node (i − 1, k) in one layer
and the input of node (i, j) in the following layer is
Ji−1
X i,j
xi,j = θi,j + wi−1,k yi−1,k (12.15)
k=1

Thus the input xi,j of node (i, j) consists of a sum of all the outputs from the previous
i,j
nodes modified by the respective inter-node synaptic weights wi−1,k and a bias θi,j .
The weights are characteristic of the connection between the nodes, and the bias of

222
 node number
↓
- Hg 2,1
w1,1 - g
H - -* g -
j=1
A@AH H @AAH@HH 
@ HHw1,1 HHH 

2,2

- g AA@@ wH2,3Hj A @
g A @ j   
g -
j=2
A @ 1,1 A @


AA @@ AA @@

-g A R g AA R


g -
j=3
AA A

AA AA




.. .. ..
A A
. . .
AU AU

-g g


g -
j = Ji

layer number → i=1 i=2 i=I

Figure 12.8: Schematic of a fully-connected multilayer ANN.

the node itself. The bias represents the propensity for the combined incoming input
to trigger a response from the node and presents a degree of freedom which gives
additional flexibility in the training process. Similarly, the synaptic weights are the
weighting functions which determine the relative importance of the signals originated
from the previous nodes.
The input and output of the node (i, j) are related by

yi,j = φi,j (xi,j ) (12.16)

where φi,j (x), called the activation or threshold function, plays the role of the biological
neuron determining whether it should fire or not on the basis of the input to that
neuron. A schematic of the nodal operation is shown in Figure 12.9. It is obvious
that the activation function plays a central role in the processing of information
through the ANN. Keeping in mind the analogy with the biological neuron, when
the input signal is small, the neuron suppresses the signal altogether, resulting in a
vanishing output, and when the input exceeds a certain threshold, the neuron fires
and sends a signal to all the neurons in the next layer. This behavior is determined by
the activation function. Several appropriate activation functions have been studied
(Haykin, 1994; Schalkoff, 1997). For instance, a simple step function can be used, but
the presence of non-continuous derivatives causes computing difficulties. The most

223
 Σ
x y
i,j i,j

Figure 12.9: Nodal operation in an ANN.

popular one is the logistic sigmoid function

1
φi,j (ξ) = (12.17)
1 + e−ξ/c
for i > 1, where c determines the steepness of the function. For i = 1, φi,j (ξ) = ξ is
used instead. The sigmoid function is an approximation to the step function, but with
continuous derivatives. The nonlinear nature of the sigmoid function is particularly
beneficial in the simulation of practical problems. For any input xi,j , the output of a
node yi,j always lies between 0 and 1. Thus, from a computational point of view, it
is desirable to normalize all the input and output data with the largest and smallest
values of each of the data sets.

Training
For a given network, the weights and biases must be adjusted for known input-output
values through a process known as training. The back-propagation method is a
widely-used deterministic training algorithm for this type of ANN (Rumelhart et al.,
1986). The central idea of this method is to minimize an error function by the method
of steepest descent to add small changes in the direction of minimization. This algo-
rithm may be found in many recent texts on ANN (for instance, Rzempoluck, 1998),
and only a brief outline will be given here.
In usual complex thermal-system applications where no physical models are avail-
able, the appropriate training data come from experiments. The first step in the
training algorithm is to assign initial values to the synaptic weights and biases in the
network based on the chosen ANN configuration. The values may be either positive
or negative and, in general, are taken to be less than unity in absolute value. The
second step is to initiate the feedforward of information starting from the input layer.
In this manner, successive input and output of each node in each layer can all be
computed. When finally i = I, the value of yI,j will be the output of the network.
Training of the network consists of modifying the synaptic weights and biases until

224
the output values differ little from the experimental data which are the targets. This
is done by means of the back propagation method. First an error δI,j is quantified by

δI,j = (tI,j − yI,j )yI,j (1 − yI,j ) (12.18)

where tI,j is the target output for the j-node of the last layer. The above equation
is simply a finite-difference approximation of the derivative of the sigmoid function.
After calculating all the δI,j , the computation then moves back to the layer I − 1.
Since the target outputs for this layer do not exist, a surrogate error is used instead
for this layer defined as

X
JI
δI−1,k = yI−1,k (1 − yI−1,k ) I,j
δI,j wI−1,k (12.19)
j=1

A similar error δi,j is used for all the rest of the inner layers. These calculations are
then continued layer by layer backward until layer 2. It is seen that the nodes of the
first layer 1 have neither δ nor θ values assigned, since the input values are all known
and invariant. After all the errors δi,j are known, the changes in the synaptic weights
and biases can then be calculated by the generalized delta rule (Rumelhart et al.,
1986):
i,j
∆wi−1,k = λδi,j yi−1,k (12.20)
∆θi,j = λδi,j (12.21)

for i < I, from which all the new weights and biases can be determined. The quantity
λ is known as the learning rate that is used to scale down the degree of change made
to the nodes and connections. The larger the training rate, the faster the network will
learn, but the chances of the ANN to reach the desired outcome may become smaller
as a result of possible oscillating error behaviors. Small training rates would normally
imply the need for longer training to achieve the same accuracy. Its value, usually
around 0.4, is determined by numerical experimentation for any given problem.
A cycle of training consists of computing a new set of synaptic weights and biases
successively for all the experimental runs in the training data. The calculations are
then repeated over many cycles while recording an error quantity E for a given run
within each cycle, where
1X JI
E= (tI,j − yI,j )2 (12.22)
2 j=1
The output error of the ANN at the end of each cycle can be based on either a
maximum or averaged value for a given cycle. Note that the weights and biases
are continuously updated throughout the training runs and cycles. The training is

225
terminated when the error of the last cycle, barring the existence of local minima,
falls below a prescribed threshold. The final set of weights and biases can then be
used for prediction purposes, and the corresponding ANN becomes a model of the
input-output relation of the thermal-system problem.

Implementation issues
In the implementation of a supervised fully-connected multilayered ANN, the user
is faced with several uncertain choices which include the number of hidden layers,
the number of nodes in each layer, the initial assignment of weights and biases, the
training rate, the minimum number of training data sets and runs, the learning rate
and the range within which the input-output data are normalized. Such choices are
by no means trivial, and yet are rather important in achieving good ANN results.
Since there is no general sound theoretical basis for specific choices, past experience
and numerical experimentation are still the best guides, despite the fact that much
research is now going on to provide a rational basis (Zeng, 1998).
On the issue of number of hidden layers, there is a sufficient, but certainly not
necessary, theoretical basis known as the Kolmogorov’s mapping neural network ex-
istence theorem as presented by Hecht-Nielsen (1987), which essentially stipulates
that only one hidden layer of artificial neurons is sufficient to model the input-output
relations as long as the hidden layer has 2J1 + 1 nodes. Since in realistic problems
involving a large set of input parameters, the nodes in the hidden layer would be
excessive to satisfy this requirement, the general practice is to use two hidden layers
as a starting point, and then to add more layers as the need arises, while keeping a
reasonable number of nodes in each layer (Flood and Kartam, 1994).
A slightly better situation is in the choice of the number of nodes in each layer
and in the entire network. Increasing the number of internal nodes provides a greater
capacity to fit the training data. In practice, however, too many nodes suffer the same
fate as the polynomial curve-fitting routine by collocation at specific data points, in
which the interpolations between data points may lead to large errors. In addition, a
large number of internal nodes slows down the ANN both in training and in prediction.
One interesting suggestion given by Rogers (1994) and Jenkins (1995) is that
J1 + JI + 1
Nt = 1 + Nn (12.23)
JI
where Nt is the number of training data sets, and Nn is the total number of internal
nodes in the network. If Nt , J1 and JI are known in a given problem, the above
equation determines the suggested minimum number of internal nodes. Also, if Nn ,
J1 and JI are known, it gives the minimum value of Nt . The number of data sets used
should be larger than that given by this equation to insure the adequate determination

226
of the weights and biases in the training process. Other suggested procedures for
choosing the parameters of the network include the one proposed by Karmin (1990)
by first training a relatively large network that is then reduced in size by removing
nodes which do not significantly affect the results, and the so-called Radial-Gaussian
system which adds hidden neurons to the network in an automatic sequential and
systematic way during the training process (Gagarin et al., 1994). Also available
is the use of evolutionary programming approaches to optimize ANN configurations
(Angeline et al., 1994). Some authors (see, for example, Thibault and Grandjean,
1991) present studies of the effect of varying these parameters.
The issue of assigning the initial synaptic weights and biases is less uncertain.
Despite the fact that better initial guesses would require less training efforts, or even
less training data, such initial guesses are generally unavailable in applying the ANN
analysis to a new problem. The initial assignment then normally comes from a random
number generator of bounded numbers. Unfortunately, this does not guarantee that
the training will converge to the final weights and biases for which the error is a
global minimum. Also, the ANN may take a large number of training cycles to reach
the desired level of error. Wessels and Barnard (1992), Drago and Ridella (1992)
and Lehtokangas et al. (1995) suggested other methods for determining the initial
assignment so that the network converges faster and avoids local minima. On the
other hand, when the ANN needs upgrading by additional or new experimental data
sets, the initial weights and biases are simply the existing ones.
During the training process, the weights and biases continuously change as train-
ing proceeds in accordance with equations (12.20) and (12.21), which are the simplest
correction formulae to use. Other possibilities, however, are also available (Kamarthi,
1992). The choice of the training rate λ is largely by trials. It should be selected to
be as large as possible, but not too large to lead to non-convergent oscillatory error
behaviors. Finally, since the sigmoid function has the asymptotic limits of [0,1] and
may thus cause computational problems in these limits, it is desirable to normalize
all physical variables into a more restricted range such as [0.15, 0.85]. The choice
is somewhat arbitrary. However, pushing the limits closer to [0,1] does commonly
produce more accurate training results at the expense of larger computational efforts.

12.4.2 Application to compact heat exchangers

In this section the ANN analysis will be applied to the prediction of the performance
of two different types of compact heat exchangers, one being a single-row fin-tube
heat exchanger (called heat exchanger 1), and the other a much more complicated
multi-row multi-column fin-tube heat exchanger (heat exchanger 2). In both cases,
air is either heated or cooled on the fin side by water flowing inside the serpentine
tubes. Except at the tube ends, the air is in a cross-flow configuration. Details of the

227
analyses are available in the literature (Diaz et al., 1996, 1998, 1999; Pacheco-Vega
et al., 1999). For either heat exchanger, the normal practice is to predict the heat
transfer rates by using separate dimensionless correlations for the air- and water-side
coefficients of heat transfer based on the experimental data and definitions of specific
temperature differences.

Heat exchanger 1

The simpler single-row heat exchanger, a typical example being shown in Figure
12.10, is treated first. It is a nominal 18 in.×24 in. plate-fin-tube type manufactured
by the Trane Company with a single circuit of 12 tubes connected by bends. The
experimental data were obtained in a variable-speed open wind-tunnel facility shown
schematically in Figure 12.11. A PID-controlled electrical resistance heater provides
hot water and its flow rate is measured by a turbine flow meter. All temperatures are
measured by Type T thermocouples. Additional experimental details can be found
in the thesis by Zhao (1995). A total of N = 259 test runs were made, of which only
the data for Nt = 197 runs were used for training, while the rest were used for testing
the predictions. It is advisable to include the extreme cases in the training data sets
so that the predictions will be within the same range.
For the ANN analysis, there are four input nodes, each corresponding to the
normalized quantities: air flow rate ṁa , water flow rate ṁw , inlet air temperature Tain ,
and inlet water temperature Twin . There is a single output node for the normalized heat
transfer rate Q̇. Normalization of the variables was done by limiting them within the
range [0.15, 0.85]. Coefficients of heat transfer have not been used, since that would
imply making some assumptions about the similarity of the temperature fields.
Fourteen different ANN configurations were studied as shown in Table 12.3. As
an example, the training results of the 4-5-2-1-1 configuration, with three hidden
layers with 5, 2 and 1 nodes respectively, are considered in detail. The input and
output layers have 4 nodes and one node, respectively, corresponding to the four
input variables and a single output. Training was carried out to 200,000 cycles to
show how the errors change along the way. The average and maximum values of the
errors for all the runs can be found, where the error for each run is defined in equation
(12.22). These errors are shown in Figure 12.12. It is seen that the the maximum
error asymptotes at about 150,000 cycles, while the corresponding level of the average
error is reached at about 100,000. In either case, the error levels are sufficiently small.
After training, the ANNs were used to predict the Np = 62 testing data which
were not used in the training process; the mean and standard deviations of the error
for each configuration, R and σ respectively, are shown in Table 12.3. R and σ are

228
 Wa
ter
in

Air

229

Wa
ter
ou
t
 2 3 4 5
1
A

A ∆P

7 6

A-A View
8

Figure 12.11: Schematic arrangement of test facility; (1) centrifugal fan, (2) flow
straightener, (3) heat exchanger, (4) Pitot-static tube, (5) screen, (6) thermocouple,
(7) differential pressure gage, (8) motor. View A-A shows the placement of five
thermocouples.

230
 -3
x 10
3

2.5

2
Errors

Maximum error

1.5

0.5
Average
Global error
Error

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
5
Number of cycles x 10

Figure 12.12: Training error results for configuration 4-5-2-1-1 ANN.

231
 Configuration R σ
4-1-1 1.02373 0.266
4-2-1 0.98732 0.084
4-5-1 0.99796 0.018
4-1-1-1 1.00065 0.265
4-2-1-1 0.96579 0.089
4-5-1-1 1.00075 0.035
4-5-2-1 1.00400 0.018
4-5-5-1 1.00288 0.015
4-1-1-1-1 0.95743 0.258
4-5-1-1-1 0.99481 0.032
4-5-2-1-1 1.00212 0.018
4-5-5-1-1 1.00214 0.016
4-5-5-2-1 1.00397 0.019
4-5-5-5-1 1.00147 0.022

Table 12.3: Comparison of heat transfer rates predicted by different ANN configura-
tions for heat exchanger 1.

defined by
Np
1 X
R = Rr (12.24)
Np r=1
v
u Np
uX (Rr − R)2
t
σ = (12.25)
r=1 Np

where Rr is the ratio Q̇e /Q̇pAN N for run number r, Q̇e is the experimental heat-transfer
rate, and Q̇pAN N is the corresponding prediction of the ANN. R is an indication of
the average accuracy of the prediction, while σ is that of the scatter, both quantities
being important for an assessment of the relative success of the ANN analysis. The
network configuration with R closest to unity is 4-1-1-1, while 4-5-5-1 is the one with
the smallest σ. If both factors are taken into account, it seems that 4-5-1-1 would be
the best, even though the exact criterion is of the user’s choice. It is also of interest to
note that adding more hidden layers may not improve the ANN results. Comparisons
of the values of Rr for all test cases are shown in Figure 12.13 for two configurations.
It is seen, that although the 4-5-1-1 configuration is the second best in R, there are
still several points at which the predictions differ from the experiments by more than
14%. The 4-5-5-1 network, on the other hand, has errors confined to 3.7%.

232
 1.15

1.1

1.05
i
Rr
R

0.95

0.9

0.85
0 10 20 30 40 50 60 70
i
r
Figure 12.13: Ratio of heat transfer rates Rr for all testing runs (× 4-5-5-1; + 4-5-1-1)
for heat exchanger 1.

233
 The effect of the normalization range for the physical variables was also stud-
ied. Additional trainings were carried out for the 4-5-5-1 network using the different
normalization range of [0.05,0.95]. For 100,000 training cycles, the results show that
R = 1.00063 and σ = 0.016. Thus, in this case, more accurate averaged results can
be obtained with the range closer to [0,1].
We also compare the heat-transfer rates obtained by the ANN analysis based
on the 4-5-5-1 configuration, Q̇pAN N , and those determined from the dimensionless
correlations of the coefficients of heat transfer, Q̇pcor . For the experimental data used,
the least-square correlation equations have been given by Zhao (1995) and Zhao et
al. (1995) to be

εNua = 0.1368Re0.585
a P ra1/3 (12.26)
Nuw = 0.01854Re0.752
w P rw0.3 (12.27)

applicable for 200 < Rea < 700 and 800 < Rew < 4.5 × 104 , where ε is the fin
effectiveness. The Reynolds, Nusselt, and Prandtl numbers are defined as follows,

Va δ ha δ νa
Rea = ; Nua = ; P ra = (12.28)
νa ka αa
Vw D hw D νw
Rew = ; Nuw = ; P rw = (12.29)
νw kw αw

where the superscripts a and w refer to the air- and water-side, respectively, V is
the average flow velocity, δ is the fin spacing, D is the tube inside diameter, and ν
and k are the kinematic viscosity and thermal conductivity of the fluids, respectively.
The correlations are based on the maximum temperature differences between the two
fluids. The results are shown in Figure 12.14, where the superscript e is used for the
experimental values and p for the predicted. For most of the data the ANN error is
within 0.7%, while the predictions of the correlation are of the order of ±10%. The
superiority of the ANN is evident.
These results suggest that the ANNs have the ability of recognizing all the con-
sistent patterns in the training data including the relevant physics as well as random
and biased measurement errors. It can perhaps be said that it catches the underlying
physics much better than the correlations do, since the error level is consistent with
the uncertainty in the experimental data (Zhao, 1995a). However, the ANN does
not know and does not have to know what the physics is. It completely bypasses
simplifying assumptions such as the use of coefficients of heat transfer. On the other
hand, any unintended and biased errors in the training data set are also picked up by
the ANN. The trained ANN, therefore, is not better than the training data, but not
worse either.

234
 10000
9000
[W]

8000
Q pcor

7000
6000
5000
Q pANN

4000
3000
2000
1000
0
0 1000 2000 3000 4000 5000 6000 7000 8000 9000
e
Q [W]
Figure 12.14: Comparison of 4-5-5-1 ANN (+) and correlation (◦) predictions for heat
exchanger 1.

235
12.5 Compressible flow
References
1. Pacheco-Vega, A., Sen, M., Yang, K.T. and McClain, R.L., genetic-algorithm-based predictions
of fin-tube heat exchanger performance, Heat Transfer 1998, Vol. 6, pp. 137–142, 1998.

2. Dı́az, G., Sen, M., Yang, K.T. and McClain, R.L., Simulation of heat exchanger performance
by artificial neural networks, to be published in International Journal of HVAC&R Research,
1999.

Problems
1. This is a problem

236
Chapter 13

Boiling

13.1 Boiling curve

13.2 Homogeneous nucleation
Problems
1. This is a problem

237
238
Part IV

Radiation

239
Chapter 14

Fundamentals of radiation

14.1 Definitions
14.2 View factors
Problems
1. Consider an unsteady n-body radiative problem. The temperature of the ith body is given
by
n
∂Ti X
= Fij (Tj4 − Ti4 ) + Qi
∂t j=1

What kind of dynamic solutions are possible?

2. The steady-state temperature distribution in a one-dimensional radiative fin is given by
dT
+ hT 4 = 0
dx
Is the solution unique and always possible?

241
242
Chapter 15

Computational methods

15.1 Monte Carlo methods

Problems
1. This is a problem

243
244
 Part V

Appendices

245
Appendix A

Routh-Hurwitz criteria

The polynomial equation

a0 sn + a1 sn−1 + . . . + an−1 s + an = 0

has roots with negative real parts if and only if the following conditions are satisfied:
(i) a1 /a0 , a2 /a0 , . . . , an /a0 > 0
(ii) Di > 0, i = 1, . . . , n
The Hurwitz determinants Di are defined by

D1 = a1
a1 a3
D2 =
a0 a2
a1 a3 a5
D3 = a0 a2 a4
0 a1 a3
a1 a3 a5 ... a2n−1
a0 a2 a4 ... a2n−2
0 a1 a3 ... a2n−3
Dn = 0 a0 a2 ... a2n−4
.. .. .. .. ..
. . . . .
0 0 0 ... an

with ai = 0, if i > n.

247
248
Bibliography

General

Alifanov, O.M., Inverse Heat Transfer Problems, Springer-Verlag, New York, 1994.
Arpaci, V.S., Kao, S.-H. and Selamet, A., Introduction to Heat Transfer, Prentice-
Hall, Upper Saddle River, NJ, 1999.
Aziz, A., Perturbation Methods in Heat Transfer, Hemisphere Pub. Corp., Washing-
ton, 1984.
Baehr, H. D. and Stephan, K., Heat and Mass Transfer, Springer, New York , 1998.
Becker, M., Heat Transfer, A Modern Approach, Plenum Press, New York, 1986.
Bejan, A., Heat Transfer, John Wiley, New York, 1993.
Bennett, C.O., Momentum, Heat, and Mass Transfer, 3rd ed., McGraw-Hill, New
York, 1982.
Burmeister, L.C., Heat Transfer, Dover, New York, 1962.
Ganapathy, V., Applied Heat Transfer, PennWell Pub. Co., Tulsa, OK, 1982.
Hewitt, G.F., Process Heat Transfer, CRC Press, Boca Raton, 1994.
Holman, J.P., Heat Transfer, 7th ed., McGraw-Hill, New York, 1990.
Incropera and DeWitt, Heat and Mass Transfer, John Wiley, 2nd Ed., 1985.
Jakob, Heat Transfer, John Wiley, 1949, Vol. 2, 1957.
Kreith, F. and Bohn, M.S., Principles of Heat Transfer, 5th ed., West Pub. Co., St.
Paul, 1993.
Lienhard, J.H., A Heat Transfer Textbook, Prentice-Hall, Englewood Cliffs, NJ, 1981.

Ozisik, M.N., Heat Transfer, A Basic Approach, McGraw-Hill, New York, 1985.

249
Suryanarayana, N.V., Engineering Heat Transfer, West Pub. Co., Minneapolis/St.
Paul, 1995.
Taine, J., Heat Transfer, Prentice Hall, Englewood Cliffs, N.J., 1993.
White, F.M., Heat and Mass Transfer, Addison-Wesley, Reading, MA, 1988.
Winterton, R.H.S., Heat Transfer, Oxford University Press, New York, 1997.
Wolf, H., Heat Transfer, Harper & Row, New York, 1983.

Conduction

Carslaw, H.S. and Jaeger, J.C., Condution of Heat in Solids, Clarendon Press, Oxford,
UK, 1959.
Grigull, U., Heat Conduction, Springer-Verlag, New York; Hemisphere Pub. Corp.,
Washington, D.C., 1984.
Kakac, S. and Yener, Y., Heat Conduction, Taylor & Francis, Washington, DC, 1993.

Ozisik, M.N., Boundary Value Problems of Heat Conduction, Dover Publications,

New York, 1968.
Poulikakos, D., Conduction Heat Transfer, Prentice Hall, Englewood Cliffs, NJ, 1994.

Convection

Arpaci and Larson, Convection Heat Transfer, Prentice Hall, 1984.

Bejan, A., Convection Heat Transfer, John Wiley, New York, 1984.
Burmeister, L.C., Convective Heat Transfer, 2nd ed., Wiley, New York, 1993.
Gebhart, B., Jaluria, Y., Mahajan, R.L. and Sammakia, B., Buoyancy-Induced Flows
and Transport, Hemisphere Publ. Corp., New York, 1988.
Jaluria, Y., Natural Convection Heat and Mass Transfer, Pergamon Press, New York,
1980.
Kakac, S., Convective Heat Transfer, 2nd ed., CRC Press, Boca Raton, 1995.
Kaviany, M., Principles of Convective Heat Transfer, Springer-Verlag, New York,
1994.

250
Kays, W.M. and Crawford, M.E., Convective Heat and Mass Transfer, 3rd ed.,
McGraw-Hill, New York, 1993.
Oosthuizen, P.H. and Naylor, D., Introduction to Convective Heat Transfer Analysis,
McGraw-Hill, 1998.
Straughan, B., The Energy Method, Stability, and Nonlinear Convection, Springer-
Verlag, New York, 1992.

Radiation

Brewster, M.Q., Thermal Radiative Transfer and Properties, John Wiley, New York,
1992.
Edwards, D.K., Radiation Heat Transfer Notes, Hemisphere Pub. Corp., Washington,
DC, 1981.
Modest, M.F., Radiative Heat Transfer, McGraw-Hill, New York, 1993.
Siegel, R. and Howell, J.R., Thermal Radiation Heat Transfer, 3rd ed., Hemisphere
Pub. Corp., Washington, D.C., 1992.

Computational methods

Bradshaw, Cebeci and Whitelaw, Engineering Calculation Methods of Turbulent Flow,

Academic Press, 1981.
Cebeci and Bradshaw, Physical and Computational Aspects of Convective Heat Trans-
fer, Springer verlag, 1984.
Jaluria, Y., Computational Heat Transfer, Hemisphere Pub. Corp., Washington,
D.C., 1986.
Jaluria, J. and Torrance, K.E., Computational Heat Transfer, Hemisphere Publ.
Corp., Washington, DC, 1986.
Patankar, S. and Spalding, B., Heat and Mass Transfer in Boundary Layers, 2nd Ed.,
International, 1970.
Patankar, S., Numerical Heat Transfer and Fluid Flow, McGraw-Hill/Hemisphere,
1980.
Roache, Computational Fluid Dynamics, Hermosa, 1976.

251
Shih, T.M., Numerical Heat Transfer, Hemisphere Pub. Corp., Washington, D.C.,
1984.
Tannehill, J.C., Anderson, D.A. and Pletcher, R.H., Computational Fluid Mechanics
and Heat Transfer, 2nd Ed., Taylor & Francis, Washington, DC, 1997.

Porous media

Kaviany, M., Principles of heat transfer in porous media, 2nd ed., Springer-Verlag,
New York, 1995.
Nield, D.A. and Bejan, A., Convection in Porous Media, 2nd Ed., Springer-Verlag,
New York, 1999.
Tseng, J.W.C., Radiant Heat Transfer in Porous Media, ?, 1990.

Phase change and two-phase flows

Carey, V. P., Liquid-Vapor Phase-Change Phenomena, Washington, DC, Hemisphere

Publ. Corp., 1992.
Collier, Convective Boiling and Condensation, McGraw-Hill, 1984.
Han and Graham, Transport Processes in Boiling and Two-Phase Systems Including
Near Critical Fluids, McGraw-Hill, 1974.
Tong, L.S. and Tang, Y.S., Boiling Heat Transfer and Two-Phase Flow, Taylor &
Francis, Washington, DC, 1997.

Experimental methods

Measurements in Heat Transfer, Eds. Eckert and Goldstein, 2nd Ed., McGraw-Hill,
1976.

Handbooks

252
CRC Handbook of Thermal Engineering, (Ed.) F. Kreith, CRC Press, Boca Raton,
FL, 2000.
Handbook of Heat Transfer Fundamentals, 2nd Ed., Eds. Rohsenow, Hartnett and
Ganic, 1985.
Numerical Methods in Heat Transfer, John Wiley, New York, NY, 1981.
Handbook of Heat and Mass Transfer, Gulf Pub. Co., Houston, 1986.
Handbook of Numerical Heat Transfer, (Eds.) W.J. Minkowycz, E.M. Sparrow, G.E.
Schneider and R.H. Pletcher, Wiley, New York, NY, 1988.
Handbook of Heat Transfer Applications, 2nd ed., McGraw-Hill, New York, NY, 1985.

Handbook of Heat Transfer Applications, 2nd Ed., Eds. Rohsenow, Hartnett and
Ganic, 1985.
Heat Exchanger Design Handbook, 2nd Ed., Eds. Rohsenow, Hartnett and Ganic,
1983.
Handbook of Single-Phase Convective Heat Transfer, (Eds.) S. Kakac, R.K. Shah, W.
Aung, John Wiley, 1983.
Handbook of Numerical Heat Transfer, John Wiley, 1988.

Serials, journals and periodicals

Advances in Heat Transfer, Eds. Hartnett and Irvine, Academic Press, 1964–.
Progress in Heat and Mass Transfer, Eds. Hartnett and Irvine, Pergamon Press,
1970–.
Proceedings of the International Heat Transfer Conference, 1958–.
ASME Journal of Heat Transfer, 19??–.
International Journal of Heat and Mass Transfer, 19??–.
AIAA Journal of Thermophysics and Heat Transfer, 19??–.
Letters in Heat and Mass Transfer, 19??–.
International Journal in Experimental Thermal and Fluid Sciences, 19??–.
Experimental Heat Transfer, 19??–.
International Journal of Heat and Fluid Flow, 19??–.
Heat Transfer Engineering, 19??–.

253
Heat Transfer–Recent Contents, 19??–.
Annual Review of Heat Transfer, Hemisphere Pub. Corp., New York, 1990-, Annual.

Numerical Heat Transfer, Part A, Applications, Hemisphere Pub. Corp., New York,
NY, ¡1989-.
Numerical heat transfer, Part B, Fundamentals, Hemisphere Pub. Corp., New York,
NY.
Annual Review of Numerical Fluid Mechanics and Heat Transfer, Hemisphere Pub.
Corp., Washington. DC., 1987-, Annual
Experimental Heat Transfer, Hemisphere Pub., Washington, DC, 1987-, Four no. a
year.
Journal of Thermophysics and Heat Transfer, American Institute of Aeronautics and
Astronautics, New York, NY, 1986-, Quarterly
International Communications in Heat and Mass Transfer, Pergamon Press, New
York, 1983-, Bimonthly

254