Title Page
 Table of Contents
 List of Tables
 List of Figures
 Formulation of problem
 Grid generation and clustering
 Numerical solution
 Evaluation of pressure distribution...
 Results and discussion
 Conclusions and recommendation...
 Biographical sketch

Title: Fluid flow and heat transfer in a single-pass, return-flow heat exchanger /
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00097409/00001
 Material Information
Title: Fluid flow and heat transfer in a single-pass, return-flow heat exchanger /
Physical Description: xiv, 125 leaves : ill. ; 28 cm.
Language: English
Creator: Yang, Song-Lin, 1952-
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 1985
Copyright Date: 1985
Subject: Heat exchangers -- Fluid dynamics   ( lcsh )
Heat -- Transmission   ( lcsh )
Mechanical Engineering thesis Ph. D
Dissertations, Academic -- Mechanical Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis (Ph. D.)--University of Florida, 1985.
Bibliography: Bibliography: leaves 122-124.
Additional Physical Form: Also available on World Wide Web
General Note: Typescript.
General Note: Vita.
Statement of Responsibility: by Song-Lin Yang.
 Record Information
Bibliographic ID: UF00097409
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000872413
notis - AEG9664
oclc - 014514574


This item has the following downloads:

PDF ( 4 MBs ) ( PDF )

Table of Contents
    Title Page
        Page i
        Page i-a
        Page ii
        Page iii
    Table of Contents
        Page iv
        Page v
    List of Tables
        Page vi
    List of Figures
        Page vii
        Page viii
        Page ix
        Page x
        Page xi
        Page xii
        Page xiii
        Page xiv
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
    Formulation of problem
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
    Grid generation and clustering
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
    Numerical solution
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
    Evaluation of pressure distribution and Nusselt number
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
    Results and discussion
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
    Conclusions and recommendations
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
    Biographical sketch
        Page 125
        Page 126
        Page 127
Full Text







To My Family


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, Li-Pin. 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.



ACKNOWLEDGMENTS....................................... iii

LIST OF TABLES............................................. vi

LIST OF FIGURES........................................... vi


ABSTRACT ..............................................xiii


I INTRODUCTION....................................1

II FORMULATION OF PROBLEM..........................6

II.1 Governing Equations......................6
11.2 Boundary Conditions......................9


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

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



IN THE TRANSFORMED PLANE.......................88




REFERENCES............................................... 122

BIOGRAPHICAL SKETCH...................................... 125


Table Page

1 Vortex Sizes for Re=100, 500, and 1000.............60

2 Inner-Wall 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


Figure Page

1 A Schematic Showing the System Under
Investigation..................................... 2

2 Specifications of the Single-Pass,
Return-Flow 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

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 C-lines 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



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)

C1-C18 metric coefficients defined in Appendix C

C19-C28 coefficients defined in Chapter IV

Dh hydraulic diameter

rep,e,ez unit vectors in the cylindrical coordinate

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

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

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, (Tw-T)/(Tw-Te)

0m dimensionless mean temperature defined in Eq.

p density

a,B ,y



















dynamic viscosity coefficient

parameter defined in Eq. (3.20)

angular coordinate in the cylindrical
coordinate system

metric coefficients defined in Eqs. (3.20) to

denotes the grid point


normal derivative

partial differentiation with respect to r


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


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




May, 1985

Chairman: C. K. Hsieh
Major Department: Mechanical Engineering

Numerical solutions for flow and heat transfer in a

single-pass, return-flow 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 Navier-Stokes 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


rapidly. The original partial differential equations are

expressed in finite-difference 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


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.


The single-pass, return-flow (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 double-pipe 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,















~ Cl)



stand-alone device. In fact, the exchanger possesses all

the virtues of the bayonet tube, such as the freedom from

stress-and-strain 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. Owens-Illinois 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


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



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 one-hundredth 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


o cr-


5- *M

o o

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=(Tw-T)/(Tw-Te), 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


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


The conditions mentioned above can be formulated as



u = 1 (2.6a)

v = 0 (2.6b)

0 = 1 (2.6c)

= -- (2.6d)

n = 0 (2.6e)


= = =z = = 0
z z z z


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)

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.


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


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 one-to-one. 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




Computational Plane

Figure 3 A Schematic Showing Mapping from the Physical
Plane to the Computational Plane



-J '



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 fine-tune 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



Ezz rr = r


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)


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 one-to-one 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


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-,j-2(a+r)zi, +ai+1,j -= -Y(i,j-1zij+1+ i (3.8)

ari-lj-2(a+y)ri,j+ari+l,j = -y(ri,j_-1+rli+1)+ (n 9)
12 2
a = [(z i, -zi ) + (ri,j -ri, )2] (.10)
4 i,j+1 ij-1 1i,j+1 i,-1

1 [( i+1,j- i-lj)( i,j+1-zi,j-1

+ (r i+ ,-r l, )(ri,+ -ri,_ )] (3.11)

Y= l(z -z )2 2(r1 -i-,j2] (3.12)
S (Zi+l,j-z i-1,j + (ri+l,j-1ri-j)

S= (Zi+l,j+l-i+1,j-1 -1j-l--l,j+) (3.1)


r = (ri+l,j+1-ri+1,j-1+ri-,j-1-ri1,j1) (3.14)

Iteration is used in conjunction with successive line
over-relaxation 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)


k+1 k k
rij = r i, + wr (r -r1 i )
1,]l C^~ 1 ^ ,.) 1* ,^


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 < 10-4 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


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


0' lipi

iiTT~~~mil IlilllliU II'IIII I '1'1** I1 1111

r 4-3;
lj niI 1---4~ i1 lii II

I i l' 'I I I

-- ~ ill' I;i I i


-, I OI h
'11 Tj

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 grid-line 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


(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


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 Newton-Raphson method may be used to find this e. In

this method, the nth iterated e is expressed as

n n-1 F(s-i
= E, 1- (3.20)
1 F1 (, -)


Fen-1 0 n- )max -1] (3.21)
F( -I Si, n-1 [(1+Ein-1
1 1,max n- i
max Ei


F(En-1 -AS 0 n-1 max-2
1 (n-l 2 1 -1+

X[n-1( -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
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)/2-1] 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

straight-forward 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.



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 non-uniform 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


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


az 28En + yzn = J2(Pz +Qz ) (3.23)

ar 26r, + yrn = J2(Pr+Qrn) (3.24)









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 -



z Dz)/J3


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.



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
=-[(Ur-p ) + -(rn-r )] (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 first-kind conditions



u = 1 (4.6a)

v = 0 (4.6b)
0 = 1 (4.6c)

Sr (4.6d)

S= 0 (4.6e)


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

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


Vorticity transport equation:

[C7+Re(C17% -C18isT7)1j =

(C8-ReC16~ )1i+1,j+(C9+ReCl6T )_ i-,j

+(C10+ReC16T) i,j+ +(C11-ReC16 ) i,j-i+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 (


Energy equation:

Cli, = (C8-PeC16n)9i+1,+(C9+PeCl6) )O ,j

+(C10+PeC16 ) i,j+1+(C11i-PeC169 )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,j-1

b = i + +J
Si+1,j+1 i+1,j-1 i-1,j-1 i-1,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 under-relaxation 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 up-wind scheme that is commonly used

in finite-differencing 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 i-1,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 = C8-Re(C16!T IC16 1j )

C21 = C9+Re(C16 + IC16I )

C22 = C10+Re(C16- + C16~ )

C23 = C11-Re(C16* C16 )

C24 = C1+2Pe( C16Tj + C16Ji )

C25 = C8-Pe(C166n C16n )

C26 = C9+Pe(C16Tn + C165 )

C27 = C10+Pe(C16P + IC16~ )


C28 = C11-Pe(C16, C164l )

Successive-over-relaxation (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 = (1-w ) + --- (C2*. + C3i_
iCj *i ), i +1,j i-l,j
k k+1 cn+ -*
+C4P + C5k1- + c6n+ + CO ) (4.19)
i,+1 ij-1 1, En

k+l k w k k+1
S+ = (1- )k + + (C20k2 + C21
S1, -w i, C19 i+1,j i-l,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 cross-derivative
_* _x _*
terms n', Q and 09 are evaluated using

-* k k+1 k+1 k
n = i+1,j+1 i+1,j-1 i-1,j-1 -i-1,j+1

S k k+1 k+1 k
Cn i+1,j+1 i+1,j-1 + i-1,-1 i-1,j+1


-* k k+l k+1 kk
n i+l,j+l i+1,j-1 i-1,j-1 ~ -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

second-kind 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 second-order central-difference 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

first-order 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

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

Ye + (j2Q + )e =, e (4.25)

Central difference formula is then used to approximate all

n-derivative quantities, and 0 is discretized by the

second-order backward difference, giving

S= 1 -i 2 4e + eimx (4.26)
max-2,J max max'

Along the axis of symmetry, Eq. (4.8a)is approximated


i,1 (z )+1 [ui,2 + (z/z)ui-l,1 (4.27)

in which the original u and u have been approximated by

(u )i,i = ui,1 ui-1,1

(u ) = u ui,
n 1, 1,1

again a first-order 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


(2) the a equation, (4.17), contains both 4 and n


(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

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

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

(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


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


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 Reynolds-number cases are

used as the initial guess for high Reynolds-number 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


The convergence criteria adopted are delineated as


(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 < 10-4
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 < 10-3 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 straight-forward manner. In

the solution of the energy equation, the second-kind

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+C25i-1 2+C26 i3

+C27 (zn/z ) +CO6 (4.33)
(z7/Z )+I 1-1,1 En

and one point aw-ay. from the-outer--wall -

[C23 ] = C240. +C250 +C260
(-/y)+1 .i,2 +1,2 i-1,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.


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 pipe-flow 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 Navier-Stokes 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


ap (5.2)
ar Re az


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/C-1 (5.5)
(8/C ) dz

In this equation,

C = r2[(1+r2) 1-2 ] (5.6)

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)

With this definition, the pressure gradients at the inner

wall become

L 1 (1 ar ) (5.10)
z 8/C+ r ar


-=_ 1 a (5.11)
ar (8/C+) az

Correspondingly, the exit pressure condition should satisfy

Px = 0 (5.12)

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


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) n-z n(r) (5.16)
az (8/C+) r J


+ 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)


+ + +
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

Ch n
0 J j w


where Ch is a dimensionless hydraulic diameter, defined as

Ch = Dh/ri


Dh is the hydraulic diameter; eO in Eq. (5.21) denotes a

mixing cup temperature, defined as [18]


J 2 u0r dr
= 1 (5.23)
m r2
2 ur dr

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

pA2u2 CpOm2


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.


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


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~


If t11111 ,I 11 lr 'ji

oi If 110 1 1I r INrr 1 11 i

'jU I li f11 l l 0 rr

0'i~' 0 0
Ii j I II II

I 41I
liii ~ll 0) a


I F1

.Ii, .114.


jill L
iil l l 0

lI i '

i i fI I' I 0 :
'ra I r 0)' II

II -I -,~ IIb

F r-4
I;j I 0

Ii il

~~t7/ ); I'I 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 no-slip 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 180-degree 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


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

Table 1

Vortex Sizes for Re=100, 500, and 1000

Reynolds Vortex
Number Size

100 0.37

500 1.37

1000 2.3




-- -


z ooo C
Li iv LU

z r
IZ II C 0 I-

-a 0- 1
cc at )-

2 o2


o4 L

Uin -


Ii C -J
I 0 -


- 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


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 180-degree 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.


Table 2

Vorticity Data Near the Turning Point

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 inner-pipe region.

For the heat exchanger, there are three places where

the pressure drop appears unusual. Right after the fluid





j :

- t

z C

_- 3





00 0 L





0 0 0 0
IC 'a- ('I + cdf


0 Oa


0 C







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 boundary-Layer 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


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


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

100 -1.0245

500 -1.0781

1000 -1.0687

* Correct value: -1










0 II

r 4 -



7 0

n a0

l 0C










0 r (



U Cl








o -




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

Nusselt-number 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 constant-wall-temperature

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, Co


a X <


0 0 0

o o
cl :




I 0

-D cu
C i

o =
I- l"

o o


0 WC
0 )


a Z

C.J 0


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



5 0

I 0
ia 0

Z a

II '

a- *

a-( 4-

orim a

t L
0 -

Co I. C N


0 0

o o
I i -


= o
-t i o

a a

2 g-1

z 0
S4 o

-,,- a)
000 0 C 0

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




I 1 i I I I 'l

S- Pr=100

A Pr=20



Figure 21 A Plot of Mean Nusselt Number

I I I I I 1111

100 L



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 (n-2) 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 Nusselt-number

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 S-line 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


Figure 22 shows the competed stream functions along

constant E-lines for the three grid systems. In these

figures, i=50 and 152 correspond to the mid-section in the

inner pipe and the annular region, respectively; i=101

corresponds to the constant C-line 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.



















Figure 22 Comparison of the Computed Stream Functions
Along Constant C-line for Different Grid Systems

(n-Tmin)/ (nmax-Tmin)

Table 6

CPU Time Required for Computation of
Stream Functions Given in Figure 22


CPU 100 500 1000


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


A numerical solution is given for studying the flow

and heat transfer in a single-pass, return-flow 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 first-order approximation. The velocity and temperature

conditions along the centerline of the heat exchanger are

also finite difference, with the first-order differential

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs