FLUID FLOW AND HEAT TRANSFER
IN A SINGLEPASS, RETURNFLOW HEAT EXCHANGER
BY
SONGLIN YANG
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1985
To My Family
ACKNOWLEDGMENTS
The author wishes to express his sincere appreciation
to his graduate supervisory committee which included
Dr. C. K. Hsieh, chairman; Dr. V. P. Roan; Dr. R. K. Irey;
Dr. C. C. Hsu; and, Dr. H. A. Ingley, III. Special thanks
go to his advisor, Dr. C. K. Hsieh, without whose
suggestion of this research topic, guidance and helpful
advice, this research would not have been possible.
The author also wishes to thank his wife, LiPin. Her
patience and encouragement made his graduate study a
beautiful memory.
Finally, the author wishes to acknowledge the
Interactive Graphics Laboratory (IGL) and the Northeast
Regional Data Center (NERDC) of the State University System
of Florida for allowing him to utilize the computing
facilities on campus.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS....................................... iii
LIST OF TABLES............................................. vi
LIST OF FIGURES........................................... vi
NOMENCLATURE.............................................ix
ABSTRACT ..............................................xiii
CHAPTER
I INTRODUCTION....................................1
II FORMULATION OF PROBLEM..........................6
II.1 Governing Equations......................6
11.2 Boundary Conditions......................9
III GRID GENERATION AND CLUSTERING .................12
III.1 Grid Generation in Domain A............. 14
111.2 Grid Clustering in Domain A.............18
III.3 Grid Generation in Domain B and C.......23
111.4 Evaluation of the Source Terms in
Poisson's Equations...................25
IV NUMERICAL SOLUTION.............................28
IV.1 Transformed Governing Equations and
Boundary Conditions...................28
IV.2 Discretization of Governing Equations
and Boundary Conditions...............31
IV.3 Numerical Solution Procedure.............38
V EVALUATION OF PRESSURE DISTRIBUTION AND
NUSSELT NUMBER .................................44
V.1 Pressure Distribution...................44
V.2 Nusselt Number..........................49
VI RESULTS AND DISCUSSION .........................53
VI.1 Velocity Field.......................... 53
VI.2 Wall Vorticity ..........................59
VI.3 Pressure on Inner Wall..................64
VI.4 Temperature Field.......................67
VI.5 Local Nusselt Number ....................75
VI.6 Mean Nusselt Number.....................75
VI.7 Truncation Error and Convergence........80
VII CONCLUSIONS AND RECOMMENDATIONS................85
APPENDIX
A DERIVATION OF TRANSFORMATION RELATIONS,
GOVERNING EQUATIONS, AND BOUNDARY CONDITIONS
IN THE TRANSFORMED PLANE.......................88
B DERIVATION AND USE OF EQUATION (3.6)..........105
C DERIVATION OF FINITE DIFFERENCE EQUATIONS.....108
D DERIVATION OF PRESSURE GRADIENT EQUATIONS.....119
REFERENCES............................................... 122
BIOGRAPHICAL SKETCH...................................... 125
LIST OF TABLES
Table Page
1 Vortex Sizes for Re=100, 500, and 1000.............60
2 InnerWall Vorticity Data Near the Turning
Point .............................................65
3 Dimensionless Pressure Gradients at Exit.......... 68
4 Comparison of the Computed Nusselt Number at
the Exit with Lundberg's Data .....................78
5 Comparison of Mean Nusselt Number Between SPRF
Heat Exchanger and Heating in a Circular Pipe.....81
6 CPU Time Required for Computation of Stream
Functions Given in Figure 22......................84
LIST OF FIGURES
Figure Page
1 A Schematic Showing the System Under
Investigation..................................... 2
2 Specifications of the SinglePass,
ReturnFlow Heat Exchanger .......................7
5 A Schematic Showing Mapping from the Physical
Plane to the Computational Plane................13
4 Computed Grid Before and After Clustering in
Domain A ........................................19
5 Grid Detail Near the Inner Pipe Wall and the
Turning Point ...................................24
6 Complete Grid System............................26
7 Block Diagram Showing the Numerical Procedure
Used to Solve Stream Function and Vorticity.....42
8 A Control Volume Used to Find the Mean
Temperature by Energy Balance ...................51
9 Velocity Fields for Re=100, 500, and 1000.......54
10 Velocity Field Near Turning Point for Re=100....55
11 Velocity Field Near Turning Point for Re=500....56
12 Velocity Field Near Turning Point for Re=1000...57
13 Vorticity on the Inner Wall.....................61
14 Pressure Distribution on the Inner Wall.........65
15 Isotherms for Pr=0.7 and Re=100, 500, and
1000 ............................................69
16 Isotherms for Pr=20 and Re=100, 500, and
1000.............................................70
17 Isotherms for Pr=100 and Re=100, 500, and
1000 ............................................71
18 Local Nusselt lumberr for Pr=0.7 and Re=100,
500, and 1000 ................................... 74
19 Local Nusselt Number for Pr=20 and Re=100,
500, and 1000................................ ... 76
20 Local Nusselt Number for Pr=100 and Re=100,
500, and 1000................................... 77
21 A Plot of Mean Nusselt Number...................79
22 Comparison of the Computed Stream Functions
Along Constant Clines for Different
Grid Systems..................................... 83
A.1 Schematic Showing Transformation from Physical
to Transformed Plane ............................89
A.2 Unit Tangent and Normal Vectors on Lines of
Constant C and n............................ ......98
A.3 Schematic Showing the Arc Length Along
Contour r ......................................100
B.1 Indexing z Axis for Eq. (B.4)..................106
viii
NOMENCLATURE
al,a2,a3 coefficients defined in Eqs. (B.5) to (B.7)
Ch dimensionless hydraulic diameter defined in Eq.
(5.24), Dh/ri
Cp specific heat at constant pressure
C+ constant defined in Eq. (5.7)
C1C18 metric coefficients defined in Appendix C
C19C28 coefficients defined in Chapter IV
Dh hydraulic diameter
rep,e,ez unit vectors in the cylindrical coordinate
system
F scalar function
h heat transfer coefficient
J Jacobian of transformation
K thermal conductivity
n number of correlated point
n unit normal vector
Nuz local Nusselt number, hDh/k
Num mean Nusselt number defined in Eq. (5.28)
P source term defined in Eq. (3.26)
p dimensionless pressure defined in Eq. (5.3),
p /(Pu )
p pressure
p+ normalized pressure defined in Eq. (5.10),
(Rexp)/(8/C )
Pex dimensionless pressure at exit, Pex/(P~2)
Pfd dimensi less fully developed pressure,
Pfd/(Pum )
pex pressure at exit section
pfd fully developed pressure
Pe Peclet number, PrxRe
Pr Prandtl number, uCp/K
Q source term defined in Eq. (3.27)
Re Reynolds number, Pumri/u
r dimensionless radial coordinate, r/ri
ri dimensionless inner pipe radius at inner pipe
region, 1
rio dimensionless inner pipe radius at annular
region, rio/ri
ro dimensionless outer pipe radius, ro/ri
r radial coordinate in the cylindrical coordinate
system
ri inner pipe radius at inner pipe region
rio inner pipe radius at annular region
r outer pipe radius
r+ ratio of inner to outer radii, rio/ro
S arc length defined in Eq. (3.18)
S arc length defined in Eq. (5.17)
T temperature
Te fluid inlet temperature
TW inner pipe wall temperature
t surface coordinate
t unit tangent vector
u dimensionless axial velocity, u /um
u axial velocity
um mean velocity
v dimensionless radial velocity, v /u*
v radial velocity
W1 relaxation factor for Eq. (4.31)
W2 relaxation factor for Eq. (4.52)
Wr relaxation factor for Eq. (3.16)
W, relaxation factor for Eq. (4.15)
WV relaxation factor for Eq. (4.19)
Wn relaxation factor for Eq. (3.20)
Wg relaxation factor for Eq. (3.21)
z dimensionless axial coordinate, z /ri
z axial coordinate in the cylindrical coordinate
system
E,n transformed coordinates
Vi normalized index defined in Eq. (3.7)
Sdimensionless stream function, V /r2u*
'p stream function
Sdimensionless vorticity, ri/um
S* vorticity
0 dimensionless temperature, (TwT)/(TwTe)
0m dimensionless mean temperature defined in Eq.
(5.25)
p density
a,B ,y
SUBSCRIPTS
i,j
max
n
r
w
wall
z
n
SUPERSCRIPTS
k
n
71
PREFIX
V
A
d
a
dynamic viscosity coefficient
parameter defined in Eq. (3.20)
angular coordinate in the cylindrical
coordinate system
metric coefficients defined in Eqs. (3.20) to
(3.12)
denotes the grid point
maximum
normal derivative
partial differentiation with respect to r
wall
wall value
partial differentiation with respect to z
partial differentiation with respect to 5
partial differentiation with respect to n
iteration count or inner iteration count
iteration count or outer iteration count
quantity that evaluated along constant C line
quantity that evaluated along constant n line
gradient operator
difference operator
differential operator
partial differential operator
S
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
FLUID FLOW AND HEAT TRANSFER
IN A SINGLEPASS, RETURNFLOW HEAT EXCHANGER
BY
SONGLIN YANG
May, 1985
Chairman: C. K. Hsieh
Major Department: Mechanical Engineering
Numerical solutions for flow and heat transfer in a
singlepass, returnflow heat exchanger are given for
Reynolds numbers (Re) of 100, 500, and 1000, and Prandtl
numbers (Pr) of 0.7, 20, and 100. The flow is modeled as
laminar, axisymmetric, incompressible, and in steady state,
and the NavierStokes equations are expressed in terms of
stream function and vorticity.
Prior to the solution of the problem by means of a
finite difference method, the original geometry of the
system in the physical plane is transformed into a
rectangular domain with square grids in a computational
plane. Grid clustering techniques are also employed to
reposition grid lines for good resolution in those regions
where the flow and temperature fields change most
xiii
rapidly. The original partial differential equations are
expressed in finitedifference forms, and in order to
assure convergence in the numerical iteration, the matrix
coefficients in the matrix equations are reconstructed to
make those diagonal terms dominant in the coefficient
matrix.
Numerical data indicate that the flow is accelerated
in the end cap. There is a recirculation zone near the
inner end of the inner pipe, where the main flow is
separated from the wall. The mean Nusselt number for the
exchanger is found to be proportional to the Reynolds
number and the Prandtl number. The present exchanger has a
mean Nusselt number that is about 1.57 to 1.87 times that
for heat transfer in a singular pipe under similar
operating conditions (Re and Pr number). The exchanger is
thus shown to be effective in heat transfer.
CHAPTER I
INTRODUCTION
The singlepass, returnflow (SPRF) heat exchanger
shown in figure 1 consists of a circular pipe inserted
inside an outer pipe. The outer pipe is sealed at one end
so that fluid entering the inner pipe may be deflected into
the annulus formed between the inner and the outer pipe of
the exchanger. Either the inner or the outer pipe wall, or
both, may be heated, and the heat input may be simulated by
different conditions at these walls.
The heat exchanger described above is structurally
similar to the bayonet tube studies by Hurd [1] in 1946;
however, they serve different purposes. In the bayonet
tube, the doublepipe construction is essentially used for
transporting fluid. Hence, the entering and exiting fluids
exchange heat with each other via the heat conducting inner
wall, and whether heat is transferred in or out of the
system depends on the temperature difference between the
fluid in the annulus and the medium surrounding it at the
shell side. This is not so, however, for the heat
exchanger, in which the inner pipe wall may be heated.
Then, the entering and exiting fluids both exchange heat
with this wall, and the system becomes a compact,
2
co
C4
co
0
a
4
C')
r
tro
i
C'
2,
C
cl
p
C,
0
Cl
C
4
c')
C
0)
0:
~ Cl)
0)
2r
"
standalone device. In fact, the exchanger possesses all
the virtues of the bayonet tube, such as the freedom from
stressandstrain problems in the inner pipe and the ease
of replacement of tubes in case of failure, while not
suffering from the lack of versatility that usually
handicaps the use of the bayonet tube.
A review of the literature reveals that there is a
lack of basic research dealing with the heat exchanger.
Hurd [1] studied the mean temperature difference in the
bayonet tube for four different arrangements of fluid flow
directions. Large temperature differences were found when
flow direction in the annulus was counter to the flow
direction of the heat exchanger medium at the shell side.
A more recent study dealt with the use of the heat
exchanger in solar collectors [2]. In this application,
the outside pipe of the exchanger becomes the inner wall of
a dewar, and the outside surface of this pipe is coated
with a spectral selective coating to enhance its heat
absorption. OwensIllinois has used this concept in
developing its SUNPAK collectors with success [3].
It is quite natural for one who is interested in fluid
flow and heat transfer analysis to consider the SPRF
exchanger as a combination of two separate components,
circular pipe and annulus, because their fluid flow and
heat transfer have been thoroughly studied in the
literature, e.g. [4]. However, because of the presence of
the round cap at one end of the heat exchanger, the fluid
flow in the upstream region may be affected by that in the
downstream region. As a remuit, to divide the exchanger
into two subsystems and to patch up their solutions at
points of matching may be an oversimplification of this
complicated problem. A realistic modelling of the system
will require the treatment of the system as a single
unit. This motivates the research of this dissertation.
In this work, the flow is assumed to be laminar,
axisymmetric, incompressible, and in steady state. Stream
function and vorticity are used in the solution of the
problem (Chapter II). The original irregular geometry of
the system in the physical plane is transformed to a
rectangular domain with square grids in the computational
plane. Grid lines are clustered to provide good resolution
in regions where flow and temperature fields change most
rapidly (Chapter III). The problem is solved numerically
with an iterative scheme and, in order to assure a
convergent solution, diagonal matrix elements in the
governing matrix equation are maximized (Chapter IV).
Pressure distribution and Nusselt number are evaluated at
Re=100, 500, and 1000, and Pr=0.7, 20, and 100 (Chapter
V). Thus, laminar flow is stable and the use of the steady
flow form of the governing equation is physically
credible. Numerical results indicate that the SPRF
exchanger is superior to heat transfer in a circular
5
pipe. The present exchanger has a mean Nusselt number that
is about 1.57 to 1.87 times that for heat transfer in a
singular pipe under similar operating conditions (Chapter
VI).
CHAPTER II
FORMULATION OF PROBLEM
The system under investigation is shown in figure 2,
where the dimensional ratios of the pipes are also provided
in the figure. Note that the thickness of the inner pipe
is only onehundredth of its radius, while the inner end of
the pipe that is facing the spherical cap is rounded. The
radius of the inner pipe was chosen such that the cross
section of flow in the pipe is equal to the ring passage
between the pipe and the cap, and this area is also equal
to the cross section of the annulus region; see dashed line
marks in the figure. The length of the pipes is
approximately twelve times the radius of the inner pipe, a
length short enough to show the entrance effect, yet long
enough so that the boundary conditions imposed at the exit
section of the annulus will have little effect on the flow
conditions inside the annulus.
II.1 Governing Equations
For the problem under investigation, the flow is
treated as steady, incompressible, and laminar. The fluid
medium is assumed to have constant thermophysical
properties. Both the flow and the temperature field are
L CM
OD
o cr
C II
5 *M
o o
0,,,
4a
S.. ~
invariant in the 0 direction; the flow is thus
axisymmetrical, which also permits the use of the stream
function in the analysis. A is further assumed that the
Eckert number is small enough so that dissipation is
negligible. Under these conditions, the momentum equations
in r and z directions can be combined to derive the
vorticity transport equation [5]
1 1 2 1 1 (
I 0 + 1 ( + +  n ) (2.1)
r r z r z r 2 +Re zz rr r r 2
r r
where the vorticity a is related to the stream function 9
by the following relationship:
z + rr r = rn (2.2)
zz rr r r
As usual, the stream function is related to the
velocity u and v by
u = r (2.3)
r r
v = (2.4)
The energy equation can be expressed in terms of
stream function as
1 1 1 1
S  e =  (e + e + e ) (2.5)
r r z r z r Pe zz rr r r
where Pe stands for the Peclet number, a product of the
Prandtl number and the Reynolds number, Pe=PrxRe.
In the formulation given above, the spatial variables
r and z have been made dimensionless by dividing the
original r* and z* by ri, the radius of the inner pipe.
The velocities u and v have been normalized with reference
to the mean velocity, um; 0 is a dimensionless temperature,
defined as O=(TwT)/(TwTe), where Tw and Te refer to the
inner wall temperature and the fluid inlet temperature,
respectively. There are five unknowns (u,v,0,v,P) in these
equations, and there are five equations to solve them.
11.2 Boundary Conditions
Because of the large number of conditions to be used
in the analysis, it is convenient to summarize them
according to the locations where they will be applied.
Inlet: The velocity and temperature are both uniform at
the inlet. Here the v component of velocity is
zero.
Exit: Because of the exit section being remote from the
inlet section, u, v, i, and n are considered
invariant in the z direction at exit. This
assumption has been supported by a study in which
it showed that the exact nature of the downstream
condition had little effect on the solution of the
problem [6]. It is further assumed that conduction
in the fluid is small at the exit section, a valid
assumption for Pe > 100 as proposed by Singh [7].
Axis of symmetry: The firstderivative of u and 0 with r
must be zero along the axis of symmetry. Here the
velocity v is also zero.
Walls: Both u and v velocities vanish at the inner and
outer pipe walls. Physically, to satisfy these
conditions, the vorticity must be continuously
generated at these walls. For heat transfer
analysis, a constant temperature condition is
imposed at the inner wall, while the outer wall is
insulated.
The conditions mentioned above can be formulated as
follows.
Inlet:
u = 1 (2.6a)
v = 0 (2.6b)
0 = 1 (2.6c)
2
=  (2.6d)
n = 0 (2.6e)
Exit:
= = =z = = 0
z z z z
(2.7a)
zz = 0 (2.7b)
Axis of symmetry:
u = v = e = = r = 0 (2.8)
r r
Inner pipe wall:
u = v = 0 = O (2.9a)
S= (2.9b)
2
Outer pipe wall:
u = v = 0n = n = 0 (2.10)
where the subscript n for e denotes the normal derivative.
The missing vorticity conditions at the walls may be
derived by using Eq. (2.2), in which a great simplification
can be made by using the facts that (i) u = v = 0 at the
walls, and (ii) the walls may be regarded as one of the
streamlines. For the time being, this condition is
expressed as
S= 1 ( + (2.11)
w = [ r zz rr r r) wall
It will be simplified later in Chapter IV.
CHAPTER III
GRID GENERATION AND CLUSTERING
Because of the irregularity of the geometry under
investigation, the solution of the problem in the original
physical plane is difficult. It is desirable to use grid
generation techniques to transform the grid points in the
physical plane to the computational plane so that
computation can be conveniently carried out in the new
plane.
Refer to the mapping diagrams shown in figure 3. The
physical plane is in (z,r) coordinates, while the
computational plane is in (c,n) coordinates. It is
necessary to transform the domain abcdefghija in the
physical plane to the domain a b c d e f g h i j a in the
computational plane, the mapping being onetoone. Notice
that the fluid enters the heat exchanger through the inner
pipe, and the fluid actions are more intense near the wall
and around the turning point (h) at the tip of the inner
pipe; a fine grid is necessary to study the fluid flow and
heat transfer in these regions. To facilitate grid
generation to meet these requirements, the total system is
subdivided into three domains: domain A in bcdghib, domain
B in abija, and domain C in fgdef; and transformations will
Constant n line
Constant E line
Physical Plane
* *
Sb c d e
I
min
Smax
Computational Plane
Figure 3 A Schematic Showing Mapping from the Physical
Plane to the Computational Plane
nmax
nmin
J '
a
I
be made separately in these regions to map the grid points
from one plane to another. Also note that, after
transformation, the streamlines in the physical plane will
correspond roughly to the constant n lines in the
computational plane.
The transformation will be accomplished in three
steps. In the first step, the domain A will be mapped to
domain A*. The next step is to finetune those grid lines
in this domain so that a dense grid is formed close to the
walls and near the turning point of the inner pipe.
Finally, grids in domains B and C are filled in so that
they match those grids already laid out in domain A. The
details of these transformations will be spelled out in the
three sections that follow.
III.1 Grid Generation in Domain A
Differential equation methods will be used to generate
grids in domain A. Specifically, the elliptic partial
differential equation technique developed by Thompson,
Thames, and Mastin (TTM) will be used [8]. In order to
simplify computation of the velocity and temperature fields
to be made later, the grids in domain A in the
computational plane are taken to be squares, and Laplace
equations
(3.1)
Ezz rr = r
and
nzz + nrr = 0 (3.2)
are used in the physical plane to generate those square
grids in the computational plane; see figure 3. It can be
shown (see Appendix A) that these equations are transformed
to the following two in the computational plane
az 28zn + YZnn = 0 (3.3)
rCE 2Br + yr = 0 (3.4)
where
2 2 2 2
a = Z + r 2 = z z + rr y = z + r (3.5)
n n Tn En 1 E
It is imperative that the physical boundary of the
inner pipe, that is, ihg, be mapped onto i h g and bcd be
mapped to b c d ; both are constant n lines in the
computational plane. Similar relations must be established
for the bi and dg lines, which become constant E lines in
the computational plane. These onetoone mapping
relations are usually referred to in the literature as
boundary conditions for the transformation equations (3.3)
and (3.4). In actuality, they provide the necessary
conditions so that constant n lines conform to the physical
boundary of the actual system.
To lay out z coordinates that enable a dense grid to
be formed near the turning point h, a cubic polynomial
equation
z = Z0 + al i + a2 + a3 (3.6)
proposed by Nietubicz et al. is used [9]. The z0 and zi
refer to the z positions at index i0 and ii, respectively;
i0 refers to the index at the turning point. Ti is a
normalized index, given as
1. i
= 0 (3.7)
1i i i
Equation (3.6) permits the z coordinates to be laid out
with the size of the increment Az increasing with the value
of the index ii. Once the positions of the first and last
point on the z axis are fixed, and correspondingly, values
for the indices i0 and if are assigned, and with the
provision of the sizes of the first and the last z
increment as additional input, this equation can be used to
generate an array of z positions having spacings as
desired. The use of this equation is discussed fully in
Appendix B.
Positions of interior nodal points (z,r) are found by
solving Eqs. (3.3) and (3.4). Since these equations are
nonlinear, a numerical technique must be used to solve
them. In this method, all derivatives in these equations
are represented by central differences, and the resulting
difference equations take the following form:
azi,j2(a+r)zi, +ai+1,j = Y(i,j1zij+1+ i (3.8)
arilj2(a+y)ri,j+ari+l,j = y(ri,j_1+rli+1)+ (n 9)
where
12 2
a = [(z i, zi ) + (ri,j ri, )2] (.10)
4 i,j+1 ij1 1i,j+1 i,1
1 [( i+1,j ilj)( i,j+1zi,j1
+ (r i+ ,r l, )(ri,+ ri,_ )] (3.11)
Y= l(z z )2 2(r1 i,j2] (3.12)
S (Zi+l,jz i1,j + (ri+l,j1rij)
S= (Zi+l,j+li+1,j1 1jll,j+) (3.1)
and
r = (ri+l,j+1ri+1,j1+ri,j1ri1,j1) (3.14)
Iteration is used in conjunction with successive line
overrelaxation to solve these equations. The relaxation
procedure is formulated as follows for each line (i.e.,
for a fixed j index):
k+1 k k
zi, = zi, + w (zi zi ) (3.15)
and
k+1 k k
rij = r i, + wr (r r1 i )
1,]l C^~ 1 ^ ,.) 1* ,^
(3.16)
where 1 < i < imax. In these equations, zi,j and ri,j are
the solutions of Eqs. (3.8) and (3.9); wz and wr are
relaxation parameters for z tnd r. The superscript k for z
and r refers to the iteration count. In the iteration, the
initial guess for the grid locations is made by using
linear interpolation of those grid already located at the
boundaries. It is found that a relaxation parameter of 1.7
yields fast convergence for both wz and wr, and the
iteration is terminated if
k+1 k 4 k+1 k 4
z. z. j < 104 and r. .r., < 10
The grid generated by the TTM solver is shown in
figure 4 (a). The spacing along the z axis is
satisfactory; however, the spacing in the r direction is
not, which will be corrected later. It should be repeated
that these grids are generated by using square grids in the
computational plane. Also notice that, for the
transformation made that is illustrated in this figure, 103
points are used along the z direction, 41 points in the r
direction.
111.2 Grid Clustering in Domain A
The grids generated in the preceding section use
Laplace equations for transformation. If Poisson's
equations are used instead, the source terms in the
19
0' lipi
11111
1~.
iiTT~~~mil IlilllliU II'IIII I '1'1** I1 1111
r 43;
lj niI 14~ i1 lii II
I i l' 'I I I
 ~ ill' I;i I i
11_3
, I OI h
'11 Tj
"4
I'' Ill".N Oil
,'il ittmi~K 11 I I'~i I IT1U!b'; Vc
Poisson's equations may be used to control grid positions,
then a separate gridline clustering effort may be
spared. Sorenson and Stegernfound out that this method of
clustering grid lines may lead to complexities or even
unsatisfactory results in the computation [10]. For this
reason, the grids are not transformed using Poisson's
equations in this work; Laplace equations are used instead,
and the generated grid lines are clustered later by means
of clustering functions. In this effort, Sorenson and
Steger's method will be used as described below.
The grid lines in the z direction have been laid out
satisfactorily; no changes need be made of them. Lines in
the r direction must be clustered so that a dense grid is
where the change in the gradient of the flow variable is
most rapid. In other words, lines of constant n in the
physical plane must be repositioned for each line of
constant 5. This clustering of n lines can be accomplished
in two steps as follows.
(1) The arc length along each constant E line in
figure 4 (a) is computed with the following relation:
i,j+1 = Si,j+[(zi,j+zi,j)2+(ri,j+ri,)21/2 (3.17)
Since the first point is indexed 1, Si,1 = 0 for the sake
of consistency. Note that j refers to the index for the
jth line of constant n; the arc length is computed for a
fixed i (or constant c). Equation (3.17) can be used to
compute the arc length for each grid point in figure 4 (a),
and this yields a table of zi,j and ri,j as functions Si,j
This part of the work is essentially used for locating
grid positions in the physical plane. Once they are found,
the constant n lines in the physical plane may be
discarded.
(2) The new grid lines in the physical plane will be
laid out. An exponential stretching function is used so
that the grid spacing between constant n lines increases
monotonically with the value of j. The following relation
meets this requirement.
S. =S. + AS (+c)1 (5.18)
Si,j+1 Si,j + S 0(1+ ) 1
Once again Si,1 = 0. In this equation, AS0 is the desired
minimum grid spacing near the wall of the outer pipe; e is
chosen such that S. = S. Equation (5.18)
1'max 1'i,max
enables arc lengths, along a constant c line, to be found
that give the desired grid locations for new n lines. The
actual positions for these grids in the physical plane are
found by interpolating the zi,j and ri,j tables constructed
earlier.
There is one parameter E that remains to be found in
Eq. (3.18). As has been mentioned, this parameter controls
the position of the outermost boundary, so
F(s) = S Si = 0 (3.19)
'1max '0max
The NewtonRaphson method may be used to find this e. In
this method, the nth iterated e is expressed as
F(en1
n n1 F(si
= E, 1 (3.20)
1 F1 (, )
where
Fen1 0 n )max 1] (3.21)
F( I Si, n1 [(1+Ein1
1 1,max n i
max Ei
and
F(En1 AS 0 n1 max2
1 (nl 2 1 1+
1
X[n1( 2)1]} (3.22)
Notice that this ei must be found for each C line in the
reconstruction of new grids.
This clustering scheme described above monotonically
increases the grid spacing away from the outer pipe wall.
This creates a slight problem for those grids near the
inner pipe wall because the rate of change of the gradient
of the flow variable is also rapid there. It is necessary
to cluster the grid lines near boundary ihg (figure 3) for
good resolution. To accomplish this, an image method is
employed as follows.
In the image method, the S.. in Eq. (3.18) is
1, max
redefined to be S. m/2, which length corresponds to
1,0max
that length to index (jmax+1)/2, or to that line at its
center. Then lines [(Jmax+1)/2+1] to jmax need not be
regenerated but laid out by regarding them to be the mirror
image of those lines from [(jmax+1)/21] to 1.
The clustered grid lines are shown in figure 4 (b),
where ASO is 0.25% of the inner pipe radius. An enlarged
grid near the turning point is shown in figure 5. The
clustered grid lines appear to be uniformly distributed
near the pipe wall and the turning point of the pipe. The
grid generation task in domain A is thus completed.
III.3 Grid Generation in Domain B and C
The grid generation in domain B and C is rather
straightforward once the grids are laid out in domain A.
Recall that the positions of the grid lines in the z
direction in the physical plane are generated by using Eq.
(3.6), in which values for z0, zf, io, if, AZ1, and Azf
must be provided as input; the grid spacing increases with
the index i. Then, in domain B, the z positions for the
grid lines must be laid out from a to b because close to
the inlet section a dense grid is desired. In domain C,
however, the direction of z must be reversed because, at
the exit section, the flow should approach fully developed;
a dense grid there is unnecessary.
Co

As for the lines of constant n in domains B and C, a
straight ray extension of those in domain A is adequate for
this study. This method of grid generation results in an
orthogonal nonuniform grid net in domains B and C. A plot
of the finished grid is shown in figure 6. Because domains
B and C are now linked to domain A, there are a total of
163 points in the z direction; the number of points across
the flow passage in the r direction is still 41, which is
unchanged.
III.4 Evaluation of the Source Terms
in Poisson's Equations
The dense grids near the walls and the center line of
the exchanger may be considered a result of heat sources
and sinks. Hence, although Poisson's equations were not
used to generate grids, the redistributed grid locations
implicate that these source terms do present now; the
original Laplace transformation relations are no longer
valid, and the forcing functions in the Poisson's equations
must now be determined.
It has been shown in Appendix A that these Poisson's
equations in the computational plane are of the following
form:
az 28En + yzn = J2(Pz +Qz ) (3.23)
ar 26r, + yrn = J2(Pr+Qrn) (3.24)
C(o.C4)
26
~S
a)
4,
m
4
C,
C)
a
a
0
C)
Ct4
where J is the Jacobian of transformation
J = zrn zr
Equations (3.23) and (3.24) may be used to derive the
source expressions as
P = (z Dr r Dz)/J3
Q = (rEDr 
(3.26)
(3.27)
z Dz)/J3
where
Dr = ar 2Brn + yrnn
Dz = az E 28z2n + yznn
Equations (3.26) and (3.27) are to be solved numerically to
find the forcing terms.
(3.25)
CHAPTER IV
NUMERICAL SOLUTION
With the grids generated in the physical plane (z,r),
it is only necessary to solve the problem in the
computational plane (E,n) so that a transformation back to
the physical plane gives the final answers to the
problem. The numerical solution of the problem is provided
in this chapter. The governing equations and boundary
conditions in the transformed plane will be given first in
the section that follows.
IV.1 Transformed Governing Equations
and Boundary Conditions
The governing equations in the computational plane are
derived by using the partial differential identities given
in Appendix A.
Vorticity transport equation:
_2 Jz Jz
r n r2 r n
JReP
=[(Urp ) + (rnr )] (4.1)
r E [( r p~E
Vorticity definition equation:
Jz Jz
aS 28n + +(J2+ ) +(j2Q_ )
En Tn r (r n
= J2r (4.2)
Velocity equations:
1 (z Z P) (4.3)
u = (zJ n znr) (4
1 (r r ) (4.4)
v =  (r n rn) (4.4)
Energy equation:
Jz Jz
00 2B0 + ye +(j2p_ )0 2 +(J Q+ )o
En n) r ( r n
JPe (45)
JPe (ne 90) (4.5)
r nC n1
Boundary conditions can be derived accordingly
(Appendix A). Notice that, in the application of boundary
conditions, all conditions of the first kind will be used
as is; no transformations of them need be made. Other
conditions must be changed to the computational plane, and
they are listed together with those firstkind conditions
below.
Inlet:
u = 1 (4.6a)
v = 0 (4.6b)
0 = 1 (4.6c)
Sr (4.6d)
S= 0 (4.6e)
Exit:
u v= = = = 0 (4.7a)
e (4.7b)
Axis of symmetry:
zu zu = 0 (4.8a)
v = 0 (4.8b)
z e z e =0 (4.8c)
En n
= 0 = 0 (4.8d)
Inner pipe wall:
u = v = e = 0 (4.9a)
1 (4.9b)
Outer pipe wall:
u = v = 0 (4.10a)
en Se = 0 (4.10b)
{ = o (4.10c)
Notice that at both inner and outer walls, the
vorticity boundary conditions can be expressed as
Sw = (4.11)
w 2 nn
Jr
Derivations of these expressions are given in the appendix.
IV.2 Discretization of Governing Equations
and Boundary Conditions
Those governing equations given above will be solved
numerically. In the solution, the grids are indexed as
(i,j), with i pointed in 5 direction and j in n
direction. The governing equations are discretized as
follows.
Vorticity transport equation:
[C7+Re(C17% C18isT7)1j =
(C8ReC16~ )1i+1,j+(C9+ReCl6T )_ i,j
+(C10+ReC16T) i,j+ +(C11ReC16 ) i,ji+C0Tn (4.12)
Vorticity definition equation:
Cli,j = C2 i1+C3 i,+C4ij+1+Ci,1
+C6n. .+COQ0 (4.13)
Velocity equations:
u ij = C12T1 C13IE (4.14)
1,] T1
vij = C14Ti C15 (
(4.15)
Energy equation:
Cli, = (C8PeC16n)9i+1,+(C9+PeCl6) )O ,j
+(C10+PeC16 ) i,j+1+(C11iPeC169 )ij, I
+CO 6 (4.16)
In these equations, CO through C18 represent metric
coefficients listed in Appendix C. Those 9, n, and 0
topped with bars are compact expressions for their
increments in i and j directions; for example,
S i+1,j i,j
n = i,j+1 i,j1
and
b = i + +J
Si+1,j+1 i+1,j1 i1,j1 i1,j+1
It is strictly a matter of computer programming expediency
that these substitutions are made.
The discretized equations given above can be solved by
an iterative method, in which the problems encountered
include the nonlinearity of the vorticity transport
equation and the high convection rate in the energy
equation. Although use of an underrelaxation scheme in
the iterative method may ease these problems, a small
relaxation factor means an excessive computer time. For
these reasons, the matrix coefficients in the matrix
equations are reconstructed to make those diagonal elements
dominant in the coefficient matrix so that an iteration
will yield a convergent solution.
The method used for maximizing the matrix elements is
somewhat similar to the scheme developed by Takemitsu [11]
and is rooted in the upwind scheme that is commonly used
in finitedifferencing flow terms in the solution of fluid
dynamics problems [12]. The vorticity transport equation
is expressed as
C19. = C20n .+C21. +C22.
1,j 1+1,j i1,j 1,j+1
+C23.i,j1 +CO"T (4.17)
This equation is a compact version of the vorticity
equation derived in Appendix C, Eq. (C.35).
The energy equation is derived as
C240ij = C250 +C26e0 +C270
+C280ij_1.COOn (4.18)
In both equations, C19 through C28 represent
C19 = C7+Re[C17*T C18in +2( C16i5j + JC16i j)]
C20 = C8Re(C16!T IC16 1j )
C21 = C9+Re(C16 + IC16I )
C22 = C10+Re(C16 + C16~ )
C23 = C11Re(C16* C16 )
C24 = C1+2Pe( C16Tj + C16Ji )
C25 = C8Pe(C166n C16n )
C26 = C9+Pe(C16Tn + C165 )
C27 = C10+Pe(C16P + IC16~ )
and
C28 = C11Pe(C16, C164l )
Successiveoverrelaxation (SOR) methods will be used
to solve these equations. In this method, the SOR versions
of (4.13), (4.17), and (4.18) are, respectively,
k+l k w( k k+l
+ 1 = (1w ) +  (C2*. + C3i_
iCj *i ), i +1,j il,j
k k+1 cn+ *
+C4P + C5k1 + c6n+ + CO ) (4.19)
i,+1 ij1 1, En
k+l k w k k+1
S+ = (1 )k + + (C20k2 + C21
S1, w i, C19 i+1,j il,j
+C22k k+l *
+C22k + C23k+1i, + CCO ) (4.20)
1 +1 j 1 EnT
k+1 k we k k+l
S = (1 e)0 + (C250 + C26e1
k k+1 *
+C27 + C288k+1 + COe ) (4.21)
2i,j+1 + 10 I n
where w., wn, and we denote the relaxation factors for 4,
n, and e, respectively. Superscript k represents the inner
iteration count in the numerical solution; the superscript
n for ni,j in (4.19) stands for the outer
iteration count discussed later. The crossderivative
_* _x _*
terms n', Q and 09 are evaluated using
* k k+1 k+1 k
n = i+1,j+1 i+1,j1 i1,j1 i1,j+1
S k k+1 k+1 k
Cn i+1,j+1 i+1,j1 + i1,1 i1,j+1
and
* k k+l k+1 kk
n i+l,j+l i+1,j1 i1,j1 ~ l,j+1
Iterations sweeps run from i=1 to imax and j=1 to jmax"
Discretization of the governing equations is now
complete; attention will now be directed to the boundary
conditions which must be processed accordingly. In the
discretization of boundary conditions, it must be noted
that, in the numerical solution to be given later, 9 and n
must be solved prior to the solution of the velocity and
temperature fields. It is thus appropriate to discretize
the and n conditions first. In doing so, conditions of
the first kind will again be used as is. Therefore, only
conditions of the second kind will be addressed. These
secondkind conditions will be discretized in the order of
their appearance in conditions (4.6) to (4.11) below.
The exit conditions Eq. (4.7a) are discretized by
using a secondorder centraldifference scheme. At the
exit section, for example, *, = 0 is discretized to be
= maxi mj 4 (4.22)
i +1,j i 1i,j
max max
This identity will be used to substitute the values of at
index (imax+1,j) when the difference equation (4.13) is
evaluated at the exit section. This method also applies to
n, u, and v.
The wall vorticity conditions (4.11) are approximated
by the relation
"i,w = Y [2(i,w+1 i,w)] (4.25)
Here the bracketed terms are used to approximate 'p and
subscript (w+1) for refers to the grid location being one
point away from the wall. This equation has only
firstorder accuracy but is adequate for use in the present
work because of its conservation property as reported by
Parmentier and Torrance [13].
Conditions for e and u are treated next, and it is
noted that, at the exit, the grids are orthogonal and
P = 0. Equation (4.5) reduces to
Jz
ae + ye + J2 + (Q + Jz Jre = e (4.24)
g nn r n r np
With the use of Eq. (4.7b) and after simplification, Eq.
(4.24) becomes
Jz
Ye + (j2Q + )e =, e (4.25)
Central difference formula is then used to approximate all
nderivative quantities, and 0 is discretized by the
secondorder backward difference, giving
S= 1 i 2 4e + eimx (4.26)
max2,J max max'
Along the axis of symmetry, Eq. (4.8a)is approximated
by
i,1 (z )+1 [ui,2 + (z/z)uil,1 (4.27)
in which the original u and u have been approximated by
(u )i,i = ui,1 ui1,1
(u ) = u ui,
n 1, 1,1
again a firstorder appro:.irmtion.
Temperature condition along the same axis is treated
in a similar manner,
1 + (z /z (4.28)
ei,1 (z /z )+1 i,2 z )i 1] (428)
n S
At the outer pipe wall,
ei,1 (B/y)+1 [ei,2 + (B/Y)Oi_1,1 (4.29)
IV.3 Numerical Solution Procedure
A close examination of Eqs. (4.13) to (4.15), (4.17),
(4.18), and (4.25) reveals that
(1) the p equation, (4.13), contains all p but one a
term;
(2) the a equation, (4.17), contains both 4 and n
terms;
(3) and o are also related in the wall vorticity
condition, Eq. (4.25);
(4) u and v are functions of 4 only; see Eqs. (4.14)
and (4.15);
(5) the e equation (4.18), contains both and e
terms.
These characteristics implicate that
(1) p and Q should be solved for simultaneously first;
(2) 6, u, and v are solved next, with the use of i as
input.
Based on these observations, an algorithm is developed
as follows:
(1) make an initial guess of and n values at all
grid points and assign these values as and
i,j,
(2) update n values at the wall using (4.23)
S2Y (n n
i,w 2 i,w+1 i,w (4.
(3) calculate new values of n from (4.20) using SOR
until convergence is attained; assign these values
as Qi,j and then damp them by using
n+1 n 1 K:
an+1 = an + w ) (4.31)
i,j ,j w i,j ,j
1,3 1, 1j
(4) calculate new values of 1 from (4.19), again using
SOR, until convergence is attained; assign these
values as ai,j and then relax them by using
n+1 n n (4.2)
= + + w ( ) (4.32)
i,j i,j 2 i,j i,j
(5) repeat steps (2) to (4) until both n and 4
converge.
Notice that in steps (3) and (4) above, iteration is
used to find convergent a and i values for each set of
assumed (or revised) and jp.ated 0 and i values; these two
steps are then called the inner iteration in the numerical
solution. Steps (2) and (5), which repeat the total effort
by renewing a and p values, are thus called the outer
iteration.
Based on the experience of numerical work accumulated
in the solution of this problem, w = 1.5 and wn = 0.17
yield convergent solutions; w1 and w2 are both taken to be
0.3, following the work by Rigal [14]. These relaxation
parameters are used consistently for all Reynolds numbers
tested in this work. In the numerical solution, results
for p and a calculated for low Reynoldsnumber cases are
used as the initial guess for high Reynoldsnumber cases.
Thus, for example, the output of t and a for Re=100 is used
as input for Re=500 case in step (1). For Re=100, the
initial guess of the vorticity values is taken to be zero
throughout. As for the stream function, a uniform flow is
assumed.
The convergence criteria adopted are delineated as
follows.
(i) The convergence for outer iteration is set at
n+1 n 4
i,w i,w
which also enables the following conditions to be
met [15]:
n+1 n < 4 a n+1 n < 104
i,j 1 i,j id
(ii) The convergence for inner iteration of Eqs. (4.19)
and (4.20) is set at
k+1 k < 103 and [k+1 k < 10
1,0 i,j i,j i,J
The inner iteration steps are limited to 50 in order to
prevent it from converging to an erroneous solution [5].
It was found that once the outer iteration has converged as
desired, the inner convergence criterion can be met with
just one iteration. A flow chart for the numerical
procedure is given in figure 7.
Once the stream function is found, the velocity field
is calculated from Eqs. (4.14) and (4.15) and the energy
equation can be solved in a straightforward manner. In
the solution of the energy equation, the secondkind
boundary conditions are treated explicitly. Miyakoda's
method [16] is used in this effort, in which the derivative
boundary conditions are incorporated into nodal equations
for interior nodes that are located one point away from the
insulated boundary (that is, the center line or the outer
wall). Hence, one point away from the center line,
Figure 7 Block Diagram Showing the Numerical Procedure
Used to Solve Stream Function and Vorticity
[C23 (z ) 1i,2 C24i+l2+C25i1 2+C26 i3
+C27 (zn/z ) +CO6 (4.33)
(z7/Z )+I 11,1 En
and one point away. from theouterwall 
[C23 ] = C240. +C250 +C260
(/y)+1 .i,2 +1,2 i1,2 i,3
+C27 ( 1 i,1+COen (4.34)
It is also noted that, before applying these equations
to find ei,2, 0i_,1 must be updated by using (4.28) and
(4.29). The relaxation parameter we is taken to be 0.7 for
all cases studied. A fast convergence results if the
iteration proceeds from i=1 to imax' that is following the
flow direction, and J=Jmax to 1, that is, from the inner
wall to the outer wall.
CHAPTER V
EVALUATION OF PRESSURE DISTRIBUTION AND NUSSELT NUMBER
Of the numerous quantities of interest in fluid flow
and heat transfer, the pressure drop and the heat transfer
coefficient are most important. The flow and temperature
fields evaluated in the numerical solution may be used to
compute them as will be shown in this chapter.
V.1 Pressure Distribution
Pressure drop in the heat exchanger can be computed
with the use of the velocity field obtained in the
numerical solution. However, because of the geometric
irregularity of the exchanger, there is a lack of a common
axis along which this pressure drop can be evaluated all
the way from inlet to exit. The centerline of the pipes,
which is commonly used as a basis for evaluating the
pressure drop in pipeflow problems, is not suited for use
in the present investigation because it terminates at the
center of the end cap. The inner pipe wall, fortunately,
has a continuous contour all the way from the inlet to
exit; it is thus used as a basis for evaluating the
pressure drop.
The NavierStokes equations can be used together with
the continuity equation and vorticity to derive the dimen
sionless pressure gradients at the wall as (Appendix D)
1 1 arn
az = " (5.1)
az Re r ar
and
ap (5.2)
ar Re az
where
*2
p = p /(pu ) (5.3)
Since only the total pressure drop and the pressure
gradient along the wall are of concern, pressure at exit is
assigned to be zero,
Pex = (5.4)
This eliminates the need for determination of an arbitrary
constant when the pressure drop is found later.
A refinement is made by defining a new dimensionless
pressure so that the flow characteristics of both
Poiseuille flow and annular flow may be incorporated.
Notice that, for a fully developed flow in an annulus [17]
Re dpfd (55)
(8/C1 (5.5)
(8/C ) dz
In this equation,
C = r2[(1+r2) 12 ] (5.6)
In(1/r+)
*2
and Pfd is dimensionless, fd = Pfd/(pu ); Re is the
Reynolds number; r+ in Eq. (5.6) is a radius ratio,
+ *
r = rio/ro
For a Poiseuille flow under a fully developed
condition, Eq. (5.5) is changed to
Re dPfd
1 (5.7)
8 dz
which can be recast as
Re dpfd + (5.8)
C (5.8)
(8/C+) dz
Comparing Eqs. (5.5) and (5.8) provides a clue for a
new pressure defined as
+ Re
p p (5.9)
(8/C+)
With this definition, the pressure gradients at the inner
wall become
+
L 1 (1 ar ) (5.10)
z 8/C+ r ar
and
=_ 1 a (5.11)
ar (8/C+) az
Correspondingly, the exit pressure condition should satisfy
Px = 0 (5.12)
ex
following Eq. (5.4).
It is also noted that, in the solution of the flow
field, the flow quantities have been assumed to be
invariant in the z direction. It follows that
dp+ dp +
e_ = 1 (5.13)
dz dz
must be met.
The formulation given above completes the preparation
of evaluating the pressure drop whose derivation will now
follow.
Along the inner pipe wall of the exchanger, the
pressure drop can be evaluated using a line integral
+ P+ + dt (5.14)
AP = p2 = t2 dt (5.14)
1 a
where t refers to a surface coordinate. The integrand in
this equation can be expressed by means of a vector
analysis as
+ + +
S= vp* t = ( e + e ) 3 t (5.15)
at az z ar r
where t is a unit vector tangential to the wall. Expres
sions for ap+/az and ap+/ar have been given previously as
Eqs. (5.10) and (5.11). In the transformed plane, they
become (Appendix A)
+ z (rn) z (rn)
1 1 z e(rn) nz n(r) (5.16)
az (8/C+) r J
and
+ r r n
__= 1 n n (5.17)
ar (8/C+) J
Notice that, in the transformed plane, the inner wall
corresponds to lines of constant n. The unit tangential
vector t in Eq. (5.15) thus takes the following form
(Appendix A):
z e + r er
t = (5.18)
The differential ar length along lines of constant is
The differential arc length along lines of constant n is
dt = Y d (
Introducing Eqs. (5.15) through (5.19) into (5.14)
gives
+ + +
p = P2 P
f2 z (rn) z (re)
1 (8/c+) rJ E
+ (1 n n)r} dd (5.20)
(8/C+) J
which is to be integrated by using the trapezoidal rule.
V.2 Nusselt Number
The local Nusselt number along the inner wall can be
evaluated using
hDh Ch aoe
Nuz ( anw
z k 0 m an w
B
Ch n
0 J j w
(5.21)
where Ch is a dimensionless hydraulic diameter, defined as
Ch = Dh/ri
(5.22)
Dh is the hydraulic diameter; eO in Eq. (5.21) denotes a
mixing cup temperature, defined as [18]
(5.19)
r2
J 2 u0r dr
= 1 (5.23)
m r2
2 ur dr
1
in the physical plane. This equation can be used directly
to evaluate the mean temperature in domains B and C, figure
3, where the grid lines in both the physical and
computational planes are orthogonal.
In domain A, the grid lines in the physical plane are
oblique; to use Eq. (5.23) will require a cumbersome
interpolation, which is inconvenient. For this reason, an
alternative expression is used as follows:
em2 = m1 + r dz (5.24)
This expression was derived based on an energy balance of
the control volume shown in figure 8. Since only the inner
wall is heated in the present exchanger, the r value in the
integrand is taken to be unity when em2 is evaluated in the
inner pipe. In the annulus, this value is taken to be 1.01
because of the dimensions of the inner pipe. Also notice
that, in the computational plane, Eq. (5.24) reduces to
z2 a
2Tr K / r dz
9
pA2u2 CpOm2
r
z2
*2*
z, pA1U 1=A2U2=qir(rj) Un
Figure 8 A Control Volume Used to Find the Mean
Temperature by Energy Balance
e2 0 + f2 r de (5.25)
m2 ml +'Pe g 1 J 0n
The scheme developed above works well in the pipe
regions where the flow sections are uniform. At the
turning tip (point h, figure 5), no hydraulic diameter is
defined. The characteristic length, Ch, at this point is
thus arbitrarily taken to be twice of that length from the
tip to the cap wall, that is, the dashed line in figure 2.
Once the local Nusselt number is found, the mean
Nusselt number may be computed using
f Nu dA f Nu r dz
Nu = z = (5.26)
m dA I r dz
Further notice that, at the inlet section, the local
Nusselt number is infinity because of the step change in
temperature. In order to circumvent this problem while
still evaluating a mean Nusselt number useful for the
present investigation, the local Nusselt number at the
inlet is estimated by extrapolating those numbers at two
adjacent points, i=2 and 3. The results are provided in
the next chapter.
CHAPTER VI
RESULTS AND DISCUSSION
Numerical results for the velocity field, temperature
field, wall vorticity, pressure, and Nusselt number will be
presented in this chapter for Reynolds numbers of 100, 500,
and 1000 and Prandtl numbers of 0.7, 20, and 100.
Computation was performed on a Harris 800 computer, and
machine plots were executed with the use of a Gould plot
routine package.
VI.1 Velocity Field
Velocity fields for three Reynolds numbers are
presented in figure 9. Enlarged views of the velocity
fields near the turning point of the inner pipe are
provided in figures 10 through 12. In making these plots,
all arrow heads of the velocity vectors are positioned
along lines that correspond to lines of constant E in the
computational plane. Hence, these arrow heads do not line
up radially when they move closer to the turning point
because they are in domain A (figure 3) and in this domain
these constant E lines are oblique in the physical plane,
also see figure 4. Certainly, the lengths of the arrows
still represent the magnitude of the velocity. It is also
54
M11 I 11141 Ell] 1 111113 1111 11 113
lil A 11 1~ l l I i i I ill11 1 11
I!" Ci 1~ 1111~11!11~
l1111114
If t11111 ,I 11 lr 'ji
oi If 110 1 1I r INrr 1 11 i
CI C
'jU I li f11 l l 0 rr
ItI
0'i~' 0 0
Ii j I II II
I 41I
liii ~ll 0) a
il
I F1
i~i
.Ii, .114.
55
jill L
iil l l 0
I I
L I
lI i '
i i fI I' I 0 :
'ra I r 0)' II
II I ,~ IIb
F r4
I;j I 0
Ii il
/;
/
~~t7/ ); I'I I
/i /
I
. /
,1
T
i
i
i I
noted that, in making these plots, the magnitude of the
velocity has been normalized by the maximum velocity found
for each Reynolds number tesejd. Inasmuch as the maximum
velocity for Re=1000 is larger, the normalized velocities
in this plot are smaller.
A close examination of the computer data shows that,
immediately following the flow entering the inner pipe,
there is a slight bulge of the velocity near the inner pipe
wall (too small to be seen clearly in figure 9). The fluid
velocity near the centerline, however, remains uniform at
its inlet value. Such an observation appears to be in
contradiction to the developing flow analyses in which such
a bulge is not reported. This inconsistency may be
ascribed to the fact that, in the present investigation,
the axial shear term in the momentum equation is retained
in the computation, whereas in the developing flow
analyses, references [18] and [19], such axial shear terms
are dropped. Physically, this bulge of velocity results
from the need for flow to satisfy the continuity equation
and the noslip condition simultaneously. Because of the
closeness of this section to the inlet, the momentum
transfer has yet to reach that far to the center region,
and as a result, only velocity nearby is affected and a
maximum velocity develops locally. Such a phenomenon has
also been discussed in references [20] and [21].
A vortex is found in the annulus right after the flow
makes the 180degree turn. As shown in figures 10 through
12, the size of the recirculation zone increases with
Reynolds number. For a quantification of this observation,
the sizes of the vortex measured in distance from the
turning point to the point of flow reattachment are listed
in Table 1. With the recirculation region, fluid moves in
a clockwise loop. The main flow detaches from the wall,
resulting in a depression of the heat transfer rate as will
be discussed later.
It is also interesting to note that, because of the
blockage of flow by the round cap, the maximum velocity in
the inner pipe is not located on its centerline, but at a
location slightly shifted away from the center. By the
same token, the maximum velocity in the annulus is not
located close to the inner wall, as is normally the case
for an annular flow [17], but is displaced slightly toward
the outer wall. This latter phenomenon is certainly a
result of the accelerated main flow that occurs near the
recirculation region being displaced away from the inner
wall.
VI.2 Wall Vorticity
The vorticity along the inner pipe wall is plotted as
a function of position dimensionlesss) in figure 13. The
turning point of the inner pipe is located at 0.5; the
0
Table 1
Vortex Sizes for Re=100, 500, and 1000
Reynolds Vortex
Number Size
100 0.37
500 1.37
1000 2.3
61
C
o
C
_C
 
o
z ooo C
o
Li iv LU
z r
IZ II C 0 I
a 0 1
cc at )
2 o2
I0

o4 L
Uin 
C)
Ii C J
I 0 
I C
 C C I I
normalized overall length is 1. Again, the vorticity is
dimensionless, normalized by the maximum vorticity value
(786.105) found for the three cases of Reynolds number
tested.
There are three places where the vorticity shows a
deviation: the inlet region, the flow region near the
turning point, and the flow in the recirculation region.
The slight increase in vorticity in the inlet region can be
attributed to the growth of the boundary layer under a no
slip condition. A larger vorticity in this region enables
these conditions to be met. Moving downstream, the
vorticity diffuses via random molecular motion and
convection; the vorticity tends to flatten out as shown in
the figure.
The turning point provides a different perspective.
Here from mathematics, this point is a point of
singularity, and the vorticity approaches infinity.
Physically, the flow accelerates in the 180degree turn,
and correspondingly, the pressure drops as will be
discussed later.
As shown in figure 13, the vorticity changes sign in
the recirculation region, also see Table 2. Finally, as
the fluid moves further downstream toward exit, three
vorticity curves all merge into one.
InnerWall
Table 2
Vorticity Data Near the Turning Point
Re
i: 100 500 1000
92 15.2166 40.7738 62.8808
93 19.4822 54.5176 84.0731
94 25.5662 74.6961 114.966
95 34.3638 105.127 161.401
96 47.5037 152.317 233.401
97 69.3520 227.864 346.933
98 114.130 351.229 513.606
99 210.448 572.942 786.105*
100 362.966 639.621 711.617
101** 278.764 244.890 180.865
102 152.656 1.12102 47.3047
103 26.7717 66.1383 73.1212
104 0.526885 43.7028 42.0327
105 6.66565 29.1544 26.0498
106 7.3296 20.5753 17.5287
107 6.6902 15.4656 12.7805
108 5.88388 12.2948 9.94225
109 5.13125 10.2201 8.10281
110 4.48652 8.83588 6.86040
* Maximum value
** Turning point
VI.3 Pressure on Inner Wall
The pressure distribution on the inner wall is plotted
for the three Reynolds nuirter! in figure 14. Here the
ordinate is a dimensionless pressure defined as Eq. (5.9);
the x axis is again dimensionless as before. Three
additional informations are included in the figure for
contrast. Two of them correspond to the pressure drop for
fully developed flow in a circular pipe and in an annulus,
in which the pressure gradient is linear and is independent
of the Reynolds number; see solid lines for Poiseuille flow
and annular flow. This is certainly a result of the
definition p+. Another set of results, which are plotted
using the stars in the figure, is taken from the numerical
work of Hornbeck [22] who evaluated the total pressure drop
for a hydrodynamic developing flow in a circular pipe;
again the slope of these data are in agreement with that
calculated in the present exchanger. It should be noted
that the level of these curves is arbitrary; only the
pressure gradients are of concern. From the figure, it is
clear that the pressure drop in the inner pipe in the
present exchanger is only affected near the turning
point. To assist this observation, a dashed line is drawn
in the figure to indicate the approximate location where
the pressure is affected in the innerpipe region.
For the heat exchanger, there are three places where
the pressure drop appears unusual. Right after the fluid
z
0
LJ
oc
(
j :
 t
z C
_ 3
2
L1
S
L
.
x
4)
00 0 L
'4
S
(%
CL
c
I)
0 0 0 0
IC 'a ('I + cdf
CD
I
so
0 Oa
C
41
U
0 C
0
ra
0
N
o
C"
0
0
o I
enters the inner pipe, the pressure gradient is positive.
This is not predicted in Hornbeck's results because, in his
analysis, the axial boundaryLayer equation was used and
the momentum equation was linearized locally at any cross
section z1 by means of the velocity at z=z1+Az. However,
in the present investigation, no linearization is made and
the axial diffusion terms are retained in the momentum
equations. As a result, this positive pressure gradient
near the inlet region may be attributed to the inclusion of
the axial diffusion terms in these equations. The
requirement for satisfying the uniform flow boundary
condition at inlet also contributes to this pressure
rise. Notice that this pressure abnormality was also
reported by Morihara and Cheng [20] for flow in a straight
channel, and they used a similar mathematical model as used
in the present study. A point of interest is that both the
pressure rise and the distance over which it occurs are
larger for larger Reynolds numbers, which is not
unexpected.
Before the flow reaches the turning point on the inner
wall, there is a rapid drop in pressure. Obviously, this
drop results from acceleration occurring before the 180
degree turn (see dashed line). There is only partial
recovery of this pressure as the fluid enters the annulus;
as would be expected, the pressure loss is partially
irreversible.
The circulation of flow in the annulus gives rise to
another region of positive pressure gradient. Finally, all
three curves merge toward the exit of the exchanger in
satisfaction of the conditions imposed there, Eqs. (5.4)
and (5.13). The flow becomes fully developed as evidenced
by the solid line drawn in the figure. To test the
accuracy of the pressure computed in the numerical work,
dpx+/dz is evaluated at the exit section. As listed in
Table 3, this pressure gradient has a value ranging from
1.0245 to 1.0781; the error is thus about 8% as compared
with the pressure gradient for annular in a fully developed
condition; also see Eq. (5.13).
VI.4 Temperature Field
Isotherms are plotted in figures 15 through 17.
Recall that, for the problems under investigation, 0=0
along the inner wall, and On=0 along the outer wall.
Notice that, each figure presents results for all three
Reynolds numbers at one Prandtl number. Hence, comparison
between the isotherm distributions of one figure provides
insight into the effect of Reynolds number; a cross
reference of the plots at the same Reynolds number offers
an assessment of the effect of the Prandtl number.
The boundary conditions required that all isotherms
emerge from the junction of the inner wall with the
inlet. When one compares any position downstream of the
Table 3
Dimensionless Pressure Gradients at Exit
+ *
Re dp /dz
ex
100 1.0245
500 1.0781
1000 1.0687
* Correct value: 1
69
0
o
0
ac
Ti
oC
0
o0
rc
0
U3
o
O
0 II
r 4 
0
[3
+
0o
01
7 0
n a0
l 0C
70
o
0
0
o
0
0
0
to
O
1'
O
0
0 r (
bl
a
P
U Cl
71
0
o
0
(d
0
0C
0
0
0
cr.
II
O
01
o
o 
,r
rt
I'
entrance, it is observed that an increase in Reynolds
number at the same Prandtl number produces a steeper
temperature gradient near the (inner) wall.
A cross reference of the plots in the three figures at
fixed Reynolds number reveals that an increase in the
Prandtl number produces the same general effect as an
increase in Reynolds numbers. This is to be expected since
increased Prandtl number implicates a smaller thermal
diffusion as compared to the momentum diffusion.
Therefore, heat entering the fluid at the inner wall is
unable to penetrate as deeply into the fluid. In fact,
figure 16 (a) shows these isotherms are so crowded near the
wall that full presentation of them becomes impossible;
only those isotherms from 0.9 to 0.99999 are thus shown.
For all cases tested, the isotherms are crowded toward
the turning point of the inner pipe. This results from the
high velocity in this region. With respect to the large
wall gradients isotherms observed near the inlet section
and near the turning point, it should be noted that these
imply high convective coefficients in these regions.
A totally different state of affairs occurs in the
recirculation region following the sharp turn. Figure 15
(b) and (c) show a smaller temperature gradient in these
regions, indicating a lower heat transfer coefficient.
Further downstream where the fluid reattaches to the
wall, a crowding of the isotherms; see figure 15 (b) and
(c), indicates slight rise of the convective coefficient.
These observations are summarized in a series of
Nusseltnumber plots presented in the next section.
VI.5 Local Nusselt Number
Local Nusselt numbers for the nine cases are shown in
figures 18 through 20 at three Prandtl numbers. Two pieces
of information are also included in figure 18; one relates
to the Nusselt number for fully developed flow in a
circular pipe under a constantwalltemperature
condition. The Nusselt number for this case is 3.658 [18];
see horizontal line drawn. The other set of data is taken
from the numerical work of Hornbeck [23] who studied the
Nusselt number under a developing flow and heat transfer
condition at the same Prandtl number (0.7). His results
are plotted using stars in the figure.
Referring to figure 18, it is seen that the present
Nuz values are in good agreement with the results of
Hornbeck for Re=100. Moving downstream, the local Nusselt
numbers asymptotically approach those for fully developed
flow. Near the turning point, the turn takes effect,
resulting in an increase in heat transfer and the local
Nusselt numbers increase rapidly. Again, to assist this
observation, a dashed line is drawn in this figure to
c'j
u
4
0
4
O
4
VC
 C, Co
II II I
a:
a X <
a:
z
a:
0 0 0
o o
cl :
0
0
0
O
c
C
0
0
I 0
D cu
C i
o =
I l"
OO
o o
C II
N
0a
0
0 WC
o
C Z
0 )
S0)
a Z
U
C.J 0
C
'L
indicate the approximate location where turn takes
effect. For the higher Reynolds number case, the local
Nusselt numbers near the inlet region of the present study
deviate only slightly from the Hornbeck's results.
For all the cases studied for the heat exchanger, a
large Nusselt number is found near the inlet, the turning
point, and in the reattachment region and low Nusselt
number occurs in the recirculation region.
It is clear that once the flow is deflected into the
annulus, a redevelopment of flow occurs. In the annulus
region, a fully developed condition is reached earlier for
fluid of low Peclet number, see figures 18 through 20. In
order to check on the computations of the present work,
Lundberg and coauthors' results [24] are used. As shown in
Table 4, Lundberg's results for fully developed velocity
and temperature profiles in an annulus appear to be in good
agreement with the present investigation that is computed
at Pr=0.7 and Re=100. If a linear interpolation based on
Lundberg's results at r'=0.5 and 1.0 can serve as a guide,
then the error is only 1.3%.
VI.6 Mean Nusselt Number
Figure 21 presents mean Nusselt number as a function
of Reynolds number with Prandtl number as a parameter.
This figure shows that mean Nusselt number increases with
both the Reynolds and Prandtl numbers. A correlation of
76
0
0
5 0
I 0
ia 0
Z a
0
II '
a *
a( 4
orim a
W Z
t L
0 
Co I. C N
CM N\J
z
77
0 0
o o
I i 
0
= o
t i o
a a
2 g1
z 0
S4 o
,, a)
000 0 C 0
C
S
2 a: X C
L OZ n
a4.' a 1S
'*0 3 bf
;i03 0
10 U 0 0
^o '* 0 cj 
Table 4
Comparison of the Computed Local Nusselt Number
at the Exit with Lundberg's Data
+
r Nusselt Number
0.5* 5.758
0.711** 5.297
1.0* 4.861
Lundberg's data
** Present investigation for Pr=0.7 and
Re=100
102
Nu
01
101
I 1 i I I I 'l
S Pr=100
A Pr=20
SPr=0.7
II I I
Figure 21 A Plot of Mean Nusselt Number
I I I I I 1111
100 L
3x101
3x103
lrrr
the results yields
Nu = 1.15 Re 028Pr 0268 (6.1)
which is valid for Re=100 to 1000, and Pr=0.7 to 100.
Variance evaluated based on (n2) degree of freedom (n
represents number of correlated point) is 2.6 [25].
The present heat exchanger is found to be far more
effective in heat transfer. To verify this point,
Hornbeck's data [23] for air (Pr=0.7) under a developing
velocity and temperature condition in a circular pipe are
used for comparison. As listed in Table 5, the mean
Nusselt number for the SPRF heat exchanger is 1.57 times
that for pipe flow when the Reynolds number is 500.
Doubling this Reynolds number enables this Nusseltnumber
ratio to be raised to 1.87. The superiority of the present
exchanger is thus clearly seen.
VI.7 Truncation Error and Convergence
In order to check the truncation error and the
convergence of the solutions presented, three different
grid systems are used. It is noted that, for the ease of
comparison, the number of grid points along the constant n
line is fixed at 163; the number of grid points along the
constant Sline is tested for 21, 41, and 61. This results
in lines of constant E for these grid systems having the
Table 5
Comparison of Mean Nusselt Number Between
SPRF Heat Exchanger and Heating in a Circular Pipe
Mean Nusselt Number
Mean Nusselt
Pr Re Hornbeck's Data* Present Study Number Ratio
0.7 500 7.4 11.62 1.57
0.7 1000 9.57 17.88 1.87
* Developing velocity and temperature profiles in a
circular pipe
same contours in the physical plane but with different grid
sizes.
Figure 22 shows the competed stream functions along
constant Elines for the three grid systems. In these
figures, i=50 and 152 correspond to the midsection in the
inner pipe and the annular region, respectively; i=101
corresponds to the constant Cline at the turning point
(see figure 5). It is clear that the computed stream
function values based on the use of grids 163x41 converge
to the results of the dense grid (163x61). A close
examination of the difference between these two sets of
data shows that the maximum deviation is only about 3%
relative to the dense grid. However, the CPU time required
for computation with the dense grid is increased by a
factor of about 3.5 (Table 6) using the Harris 800 computer
on campus. For this reason, a denser grid is not used for
the present computation. The results for 163x41 are
considered converged.
0.5
0.4
0.3
0.2
0.1
0.0
0.
0.5
0.4
0.3
0.2
0.1
0.0
0.
0.5
0.4
0.3
0.2
0.1
0.0
0.
Figure 22 Comparison of the Computed Stream Functions
Along Constant Cline for Different Grid Systems
(nTmin)/ (nmaxTmin)
Table 6
CPU Time Required for Computation of
Stream Functions Given in Figure 22
SRe
CPU 100 500 1000
time
Grids
163x21 1 hr 1 hr
30 min 57 min 17 min
163x41 4 hr 4 hr 6 hr
55 min 14 min 19 min
163x61 12 hr 10 hr 18 hr
15 min 3 min 16 min
* CPU time is based
computer on campus
on the use of Harris 800
CHAPTER VII
CONCLUSIONS AND RECOMMENDATIONS
A numerical solution is given for studying the flow
and heat transfer in a singlepass, returnflow heat
exchanger. Stream function and vorticity are used in the
formulation, and the original irregular geometry of the
exchanger in the physical plane is transformed into a
rectangular domain with square grids in a computational
plane. Elliptic partial differential equations are used in
the transformation, and grid lines are clustered with the
use of separate clustering functions which allow fine grids
to be formed in regions where the flow and temperature
fields change most rapidly. The source terms in the
Poisson's equations are also evaluated in the solution.
A finite difference method is used to solve the
problem in the computational plane. In the numerical
algorithm, in order to assure a convergent solution in the
iteration, the diagonal terms in the resulting matrix
equations are reconstructed to make them dominant.
Vorticity boundary conditions at the walls are expressed in
a firstorder approximation. The velocity and temperature
conditions along the centerline of the heat exchanger are
also finite difference, with the firstorder differential
