
Citation 
 Permanent Link:
 https://ufdc.ufl.edu/UF00097409/00001
Material Information
 Title:
 Fluid flow and heat transfer in a singlepass, returnflow heat exchanger
 Creator:
 Yang, SongLin, 1952 ( Dissertant )
Hsieh, C. K. ( Thesis advisor )
Roan, V. P. ( Reviewer )
Irey, R. K. ( Reviewer )
Hsu, ChenChi ( Reviewer )
Ingley, H. A. ( Reviewer )
 Place of Publication:
 Gainesville, Fla.
 Publisher:
 University of Florida
 Publication Date:
 1985
 Copyright Date:
 1985
 Language:
 English
 Physical Description:
 xiv, 125 leaves : ill. ; 28 cm.
Subjects
 Subjects / Keywords:
 Boundary conditions ( jstor )
Geometric lines ( jstor ) Geometric planes ( jstor ) Heat transfer ( jstor ) Inlets ( jstor ) Nusselt number ( jstor ) Pressure gradients ( jstor ) Reynolds number ( jstor ) Velocity ( jstor ) Vorticity ( jstor ) Dissertations, Academic  Mechanical Engineering  UF Heat  Transmission ( lcsh ) Heat exchangers  Fluid dynamics ( lcsh ) Mechanical Engineering thesis Ph. D
 Genre:
 bibliography ( marcgt )
nonfiction ( marcgt )
Notes
 Abstract:
 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, axisymraetric, 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 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.
 Thesis:
 Thesis (Ph. D.)University of Florida, 1985.
 Bibliography:
 Bibliography: leaves 122124.
 Additional Physical Form:
 Also available on World Wide Web
 General Note:
 Typescript.
 General Note:
 Vita.
 Statement of Responsibility:
 by SongLin Yang.
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for nonprofit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
 Resource Identifier:
 029491645 ( AlephBibNum )
AEG9664 ( NOTIS ) 014514574 ( OCLC )

Downloads 
This item has the following downloads:

Full Text 
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

Full Text 
PAGE 1
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
PAGE 2
UNIVERSITY OF FLORIDA 3 1262 08552 3412
PAGE 3
To My Family
PAGE 4
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 wnose 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. Ill
PAGE 5
TABLE OF CONTENTS Page ACKNOWLEDGMENTS iii LIST OF TABLES vi LIST OF FIGURES vii NOMENCLATURE ix ABSTRACT xiii CHAPTER I INTRODUCTION 1 II FORMULATION OF PROBLEM 6 11. 1 Governing Equations 6 11. 2 Boundary Conditions 9 III GRID GENERATION AND CLUSTERING 12 111.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 III. 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 IV
PAGE 6
VI RESULTS AND DISCUSSION 53 VI. 1 Velocity Field 55 VI. 2 Wall Vorticity 59 VI. 3 Pressure on Inner Wall 64 VI. 4 Temperature Field 6? 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
PAGE 7
LIST OF TABLES Table P^gs 1 Vortex Sizes for Re=100, 500, and 1000 60 2 InnerWall Vorticity Data Near the Turning Point 63 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 VI
PAGE 8
LIST OF FIGURES Figure Page 1 A Schematic Showing the System Under Investigation 2 2 Specifications of the SinglePass, ReturnFlow Heat Exchanger 7 3 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 r'ind 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 Re100, 500, and 1 000 70
PAGE 9
17 Isotherms for Pr=100 and Re=100, 500, and 1 000 71 18 Local Nusselt Number for Pr=0.7 and Re=100, 500 , and 1 000 74 19 Local Nusselt Number for Pr=20 and Re=100, 500 , and 1 000 76 20 Local Nusselt Number for Pr=100 and Re=100, 500 , and 1 000 77 21 A Plot of Mean Nusselt Number 79 22 Comparison of the Computed Stream Functions Along Constant ?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 5 and n 98 A. 3 Schematic Showing the Arc Length Along Contour r 100 B.1 Indexing z Axis for Eq. (B.4) 106 Vlll
PAGE 10
NOMENCLATURE a>, 32,33 coefficients defined in Eqs. (B.5) to (B.7) Ci. diraensionless hydraulic diameter defined in Eq, (5.24), D^/rCp 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 D^ hydraulic diameter e^.e^je^ 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 Nu^ local Nusselt number, hD^/k NUjjj mean Nusselt number defined in Eq. (5.28) P source term defined in Eq. (3.26) p diraensionless pressure defined in Eq. (5.3), p* pressure ix
PAGE 11
p+ normalized pressure defined in Eq. (5.10), (Rexp)/(8/C^) diraensionless pi;pssure at exit, Pex^^'^^m ^ Pex pÂ„^ djmensignless fully developed pressure, p pressure at exit section pt, fully developed pressure Pe Peclet number, PrxRe Pr Prandtl number, uCp/K Q source term defined in Eq. (3.27) Re Reynolds number, Pu^^r^/y r dimensionless radial coordinate, r/rj_ r^ dimensionless inner pipe radius at inner pipe region, 1 10 dimension^ess^inner pipe radius at annular region, r^^/v^ * , * r diraensionless outer pipe radius, r^/rj^ r* radial coordinate in the cylindrical coordinate system rinner pipe radius at inner pipe region r* inner pipe radius at annular region 10 ^ ^ r outer pipe radius r"*" ratio of inner to outer radii, ^iq/^q S arc length defined in Eq. (3.18) S arc length defined in Eq. (3.17) T temperature Tg fluid inlet temperature Ty inner pipe wall temperature
PAGE 12
* , * t surface coordinate t unit tangent vector u diraensionless axial velocity, u /u^ u axial velocity Ujjj mean velocity V diraensionless radial velocity, v /u V radial velocity m W> relaxation factor for Eq. (4.31) Wp relaxation factor for Eq. (4.32) Wj. relaxation factor for Eq. (3.16) Wg relaxation factor for Eq. (4.15) W^ relaxation factor for Eq. (4.19) Wjj relaxation factor for Eq. (3.20) Wq relaxation factor for Eq. (3.21) z diraensionless axial coordinate, z /v^ z axial coordinate in the cylindrical coordinate system C,n transformed coordinates ^j^ normalized index defined in Eq. (3.7) ij) diraensionless stream function, i) /va^u ^ stream tunction Ji diraensionless vorticity, n r^i/Uju ft vorticity Q diraensionless temperature, (T^T)/(T^Tg ) Qjj. diraensionless mean temperature defined in Eq. (5.25) P density XI
PAGE 13
y dynamic viscosity coefficient e parameter defined in Eq. (3.20) (f) angular coordinate in the cylindrical coordinate system a,B,Y metric coefficients defined in Eqs. (3.20) to (3.12) SUBSCRIPTS i,j denotes the grid point max maximum n normal derivative r partial differentiation with respect to r w wall wall wall value z partial differentiation with respect to z 5 partial differentiation with respect to C Ti partial differentiation with respect to n SUPERSCRIPTS k iteration count or inner iteration count n iteration count or outer iteration count C quantity that evaluated along constant C line n quantity that evaluated along constant n line PREFIX V gradient operator A difference operator d differential operator 3 partial differential operator Xll
PAGE 14
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, axisymraetric, incompressible, and in steady state, and the NavierStok.es 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 Xlll
PAGE 15
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. XIV
PAGE 16
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 hent with this wall, and the system becomes a compact.
PAGE 17
c o Â•H Â•P CO bO Â•H m 0) > 0) T3 C e , W (U x: W) c Â•H IS o Si CO +> CO B (U jC O CO Â•H
PAGE 18
standalone device. In fact, the exchanger possesses all tne 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 [2j. 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
PAGE 19
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 result, 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, axisyrametric, 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 l.nes 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
PAGE 20
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).
PAGE 21
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
PAGE 22
* Â•!K O * Â•!* O U Qi CD Xi O X P cd (U > o rH c +9 (U ce: 03 CQ CO PU I fH C o CO c o Â•H P CO o Â•H Â«H Â•H O
PAGE 23
invariant in the ^ direction; the flow is thus axisymmetrical , which also permits the use of the stream function in the analysis. It 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] Â— il> n 1^0. + ^ ^p =4Â— (^ + ^ + Â— a ^) (2.1) r r z r z r 2 z Re zz rr r r 2 r r where the vorticity f2 is related to the stream function ij) by the following relationship: *^^ + 'i'^^ T"^t. = r*" (2.2) zz rr r r As usual, the stream function is related to the velocity u and v by u = 7 'l'^ (2.3) V = f 1^2 (2.4) The energy equation can be expressed in terms of stream function as 2. ip Q _J_^ = 1_ (0 ^.0 +J_0) (2.5) r r z r ^z r Pe ^ zz rr r r'^ \^'JJ
PAGE 24
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 diraensionless by dividing the original r and z by rj_, the radius of the inner pipe. The velocities u and v have been normalized with reference to the mean velocity, Ujj,; is a dimensionless temperature, defined as e=(T^T)/(T^TQ) , where T^ and Tg refer to the inner wall temperature and the fluid inlet temperature, respectively. There are five unknowns {n,v ,Q,^ ,^) in these equations, and there are five equations to solve them. I I. 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, tj , and 9. 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
PAGE 25
10 in the fluid is small at the exit section, a valid assumption for Pe > 100 as proposed by Singh [?]Â• Axis of symmetry: The first* derivative of u and Q 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 = (2.6b) = 1 (2.6c) i> = j(2.6d) n = (2.6e) Exit; u =v = ^ = Q =0 (2.7a)
PAGE 26
11 \z =0 (2.7b) Axis of symmetry: u =v = e = ri, = Q. = (2.8) Inner pipe wall: u=v=0=O (2.9a) 4= J (2.9b) Outer pipe wall: u = v = G=i(;=0 (2.10) n where the subscript n for 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 = at the walls, and (ii) the walls may be regarded as one of the streamlines. For the time being, this condition is expressed as n=[ Â— (i; + ^ _1^)] (2.11 ) w r zz rr r ^r wall It will be simplified later in Chapter IV.
PAGE 27
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 (Â£;,n) coordinates. It is necessary to transform the domain abcdefghija in the physical plane to the domain abcdefghija 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 12
PAGE 28
Constant n line Constant E, line Physical Plane 13
PAGE 29
14 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. I II . 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 [8j. 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 5 + K = (3.1) zz rr
PAGE 30
15 and n + n =0 (3.2) 'zz rr v^"^y 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 "^C 2^^. ^ ^^n = Â° (3.3) ar^^ 26r^ + yr =0 (3.4) 55 5n nri where 2 2 2 2a = z + r , 3 = z z + r r , y = z + v (3.5) Tin i T] E. r\ E, t, It is imperative that the physical boundary of the X^ Xinner pipe, that is, ihg, be mapped onto i h g , and bed be KX)(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 5 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
PAGE 31
16 equation z, = Zq . a^., . a^tj . a^^l (3.6) proposed by Nietubicz et al. is used [9]. The Zq and Zj_ refer to the z positions at index iQ and i^, respectively; Iq refers to the index at the turning point. "V^ is a normalized index, given as 1 . ip, 4. = r^ r^ (3.7) 1 Ir. 1 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 ij_. Once the positions of the first and last point on the z axis are fixed, and correspondingly, values for the indices Iq and i^ 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
PAGE 32
17 difference equations take the following form: az._,^.2(an)^i,jaz,^^^. = "Y ( ^i , j_, ^^i , j,i ) ^ f ^.^ (38) ari_,,j2(a^y)r,^..ar,,,^. = "T ( ^i , j_, ^i , j.^ ) i ?^^ (39) where Â« =Tt^^i, 0^1^1.31^'^ (^i,d.ii,di)'^ (5.10) e =T t(^i.i,j"ii,j^("i,j^i"i,ji^ h^u^,^^i^,j^ ^ ^'u^,3'i^,3^ ^ (3.12) ^n ^ ("i^1,j^1^i1,jl'"i1,jl"^i1,j1^ (^Â•''^^ and ^5n = (^i.1.j.1^i.1,j1^^i1,d1^1,j.l) (^'^^ Iteration is used in conjunction with successive line overrelaxation to solve these equations. Tne relaxation procedure is formulated as follows for each ti line (i.e., for a fixed j index): z. . = z. . + w^ (z. .z. .) (3.15) and r^"''. = r^ , ^ w^ (r. ,r^ ,) (3.16)
PAGE 33
18 where 1 < i < imcv In these equations, zand r,.: are met A '> J 'f o the solutions of Eqs. (3.8) and (3.9); w^ and Wj, are relaxation parameters for z find 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 w^. and Wj, , and the iteration is terminated if k+1 k z . .z . < 10 4 and k+1 k r . .r . . 1,3 i.J < 10 4 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. III. 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
PAGE 34
19 Â«! C Â•H cd a o Q bO C Â•H Li 0) P ra iH U (^ p <: C CO 0) o m t4 o T3
PAGE 35
20 Poisson's equations may be used to control grid positions, then a separate gridline clustering effort may be spared. Sorenson and Steger* found 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 effor", 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 C This clustering of n lines can be accomplished in two steps as follows. (1) The arc length along each constant 5 line in figure 4 (a) is computed with the following relation: S. . , = S. .+[(z. . ,z. .)^+(r. . ,r. .)^]^^^ (3.17) i,J+1 i,J 1,0+1 1,0 1,0+1 1,0 Since the first point is indexed 1 , Sj_ , =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
PAGE 36
21 fixed i (or constant 5). 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 z,.; and r,a as functions S; ^. This part of the work is essentially used for locating grid positions in the physical plane. Once they are found, the constant ri 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. h.i.^ = s^_ . . aSqCUcjJ'' (3.18) Once again S. . = 0. In this equation, aSÂ„ is the desired minimum grid spacing near the wall of the outer pipe; e is chosen such that S. . = S. . . Equation (3.18) l"i in '. ' I '^max '''max enables arc lengths, along a constant 5 line, to be found that give the desired grid locations for new ^ lines. The actual positions for these grids in the physical plane are found oy interpolating the Zj_ ^ and rj_ ^ tables constructed earlier . There is one parameter Â£ that remains to be found in Eq. (3.18). As has been mentioned, this parameter controls the position of the outermost boundary, so
PAGE 37
22 F(e) = S. . S =0 (3.19) ^"^raax '"'max The NewtonRaphson method may be used to find this e. In this method, the nth iterated e is expressed as F(e"~^) ,n ^ ^n1 3^_^ (3.20) ^ ^ F'(sr') 1 where FC."') = S, , . ^ UU^Ih^"'^'' U (3.21 and (e. ) X[ernj._2)1]} (3.22) 1 "max Notice that this e^ must be found for each E. 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.
PAGE 38
23 In the image method, the S. . in Eq. (3.18) is '^max redefined to be S. . /2, which length corresponds to 'Â•^max that length to index ( Jniax"^^ ^Z^' ^"^ ^Â° that line at its center. Then lines [ ( Jmax'^'' ^/^"^^ ^ ^Â° Jmax "^^d not be regenerated but laid out by regarding them to be the mirror image of those lines from [{j^g^^ + '])/2']] to 1. The clustered grid lines are shown in figure 4 (b), where ASq 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 Zq, z^, Iq, i^, Az^ , and Az^ 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.
PAGE 39
24 p c Â•H O ou bO C Â•H C u =1 0) X! Â•P 73 C CO (U +> CO 0) 3 Â•H CO P (U Q XJ Â•H CD bO Â•H
PAGE 40
25 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_. 26z. + TZ^^ = J^(Pz^+Qz ) (3.23) 55 Cn nn S n ar,, 2Br, + yr = J^(Pr,+Qr ) (3.24) 55 5n nn <; n
PAGE 41
26 nrtis
PAGE 42
27 where J is the Jacobian of transformation ^ = ^5^ ^n^? (3.25) Equations (3.2:)) and (3.24) may be used to derive the source expressions as v/here and P = (z^Dr r^Dz)/J= (rpDr z^Dz)/JDr = ar,^ 26r^ + yr CC Cn r\r] Dz = az^^ 26z^ + YZ 5? 5n nn (3.26) (3.27) Equations (3.26) and (3.27) are to be solved numerically to find the forcing terras.
PAGE 43
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 (5,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: t2 5 Jz p Jz. ^^W^,\^ .^(.^,^.^,^)] (4.1) 28
PAGE 44
29 Vorticity definition equation: Jz Jz c^^^^ 23^^^ . Y^^^ .(J^P. ^)*^ .(jV ^)^^ = J VQ. (4.2) Velocity equations 1 ^ = jF ^^5^ ^^^ (4.3) (4.4) Energy equation: Â«^? 230^^ . ye^^ .(j2p4^)0^ .(j2q. _!l)e^ i^(,^e^ ^%) (4.5) 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 V = = 1 (4.6a) (4.6b) (4.6c)
PAGE 45
^ = 2 30 li (4.6d) Exit; a = M (4.6e) u =v = i> = a = (4.7a) 9 = ^ e (4.7b) Axis of symmetry; z u z u, = (4.8a) V = (4.8b) z z =0 (4.8c) ^], = a = (4.8d) (4.9a) (4.9b) (4.10a) (4.10b) (4.10c) Inner pipe
PAGE 46
31 "w = 3^n ^4.11) J r Derivations of these expressions are given in le 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 C direction and j in n direction. The governing equations are discretized as follows. Vorticity transport equation: [cy+ReCciy'F.cisiF^)]^. . = (C8ReCl6^^)^^^^^j + (C9+ReCl6l'^)^._^^j + {C^0+ReC^6i>.)^. . .+(C11ReCl6ij;.)f^. . ^+COn (4.12) Vorticity definition equation: C1t, . . C2*l.1,o*"ti.i_j.C4*,^j,,.C5*._j., +C6fi. .+CO^r^ (4.1J)) Velocity equations u. . = C12^Â„ C13K (4.14) V. . = C14^^ C15^g (4.15)
PAGE 47
52 Energy equation: Cie. . = (C8PeCl6?^)e +(C9+PeCi6;F_)e + (C10HPeCl6;F.)e. . . + (C11PeCl6^.)e +C0 S^^ (4.16) In these equations, CO through C18 represent metric coefficients listed in Appendix C. Those 'I', ^, and Q topped with bars are compact expressions for their increments in i and j directions; for example, h = h.^,j ni,j and *n "^1,0 + 1 "^1,01 '''^n '''i + 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 underralaxation scheme in the iterative method may ease these problems, a small relaxation factor means an excessive computer time. For
PAGE 48
33 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 terras in the solution of fluid dynamics problems [12]. The vorticity transport equation is expressed as C19^. . = C20f2. . .+C21f^. . .+C22J^. . . i,J 1+1 , J 11 ,J i,J+1 +C23^^^j_^+C0n^^ (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 C240. . = C25e. . .+C26e. . .+C27Q. . 1. J 1+1 ,3 11 ,3 i,J+1 +C280, . .^COO^ (4.18) In both equations, C19 through C28 represent C19 = C7+Re[C17K CIS'J^^ +2(C16''^ + ci6;F^)] C20 = C8Re(Cl6iF^ Cl6ii^^ )
PAGE 49
C21 = C9+Re(Cl6^^ + Cl64i J ) C22 = C10+Re(Cl6^^ + Cl64iJ ) 34 C23 = C11Re(Cl6^^ C16^, and C24 = C1+2Pe(Cl6^^ + C16^^ ) C25 = C8Pe(Cl6?^ C16^^ ) C26 = C9+Pe(Cl6<'^ + C16^^ ) C27 = C10+Pe(Cl64'. + C16^^ ) C28 = C11Pe(Cl6ii)^ C164;^ ) 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, 'i:l = ^^V^i,o^cf^^
PAGE 50
35 *=279j,j,1 * C28ej;].^ . COe;^) (4.21) where w^, v^, and Wq denote the relaxation factors for 4), f^, and 0, respectively. Superscript k represents the inner iteration count in the numerical solution; the superscript n for J^,A in (4.19) stands for the outer > J iteration count discussed later. The crossderivative _* _* _* terms ^ c > ^r t ai^d Â©^ are evaluated using T* _ ,k k+1 k+1 _ k '''en '^i + l.j + l " '^i+Ljl " "^il.jl *i1,j + 1 ^^ = ^ Â• .^ Â• ^ ^ Â• .^ Â• , + ^ Â• 1 y, ^ y, y, Cn 1+1, j + 1 1 + 1, j_i i_i,j_i 11, j + 1 and "^n ^i + 1,j + 1 ' ^i + 1,j1 " ^i1,j1 " ^i1,j + 1 Iterations sweeps run from i = 1 to ijjj^^ and j = 1 to jjjjax* 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, ^ and ^ must be solved prior to the solution of the velocity and temperature fields. It is thus appropriate to discretize the ^ and ^ conditions first. In doing so, conditions of
PAGE 51
36 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, i)^ = is discretized to be il;. . . = ^, (4.22) ^^max^^'J '^max^'J This identity will be used to substitute the values of \li at index (imax+'''j^ ^^^^ ^^^ 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 Here the bracketed terms are used to approximate '>Â„^> s^i^d 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 Parraentier and Torrance [13].
PAGE 52
37 Conditions for e and u are treated next, and it is noted that, at the exit, the grids are orthogonal and ij;_ _ Q. Equation (4.5) reduces to Jz "^5 ' '\n ' '^''\ ' ^'^^ ' r^^ = ^ %^ ^4.24) With the use of Eq. (4.7b) and after simplification, Eq. (4.24) becomes Jz ye + (J^Q + Â— r^)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 G =^ (Oi _2 . 40. _^ . + 30. .) (4.26) ^ ^ ^max "^'^ ^max ' '^ ^max'^ Along the axis of symmetry, Eq. (4.8a)is approximated by ^.1 = (zyz^).1 tUi,2 ^ ^%/^)^1.1^ ^4.27) in which the original u^ and u have been approximated by ^^)i,1 = ^.1 ^1,1
PAGE 53
38 (u ). . = u. u. . n 1 , I 1 ,^ 1,1 again a firstorder approximation. Temperature condition along the same axis is treated in a similar manner, Qi.1 = (z /z ).1 tei,2 ^ (%/^)Qi1,lJ (4.28) At the outer pipe wall, ^i,1 =Ti7TnT ^Â®i,2 ^ (e/^)9i1,1^ (4.29) IV. 3 Numerical Solution Procedure A close examination of Eqs. (4.13) to (4.15), (4.17), (4.18), and (4.23) reveals that (1) the ^ equation, (4.13), contains all ^i but one q terra; (2) the Q equation, (4. 17), contains both 4, and q terms; (3) ill and ^ are also related in the wall vorticity condition, Eq. (4.23); (4) u and v are functions of ijj only; see Eqs. (4.14) and (4.15); (5) the equation (4. 18), contains both 4; and terras .
PAGE 54
39 These characteristics implicate that (1) 4) and Q. should be solved for simultaneously first; (2) 6, u, and v are solved next, with the use of ^ as input. Based on these observations, an algorithm is developed as follows: (1) make an initial guess of i^ and Q values at all grid points and assign these values as ^^ a and "Id(2) update n values at the wall using (4.23) n 2y / ,n ^,w =^ ^"^^ , n s j2^ ^^i,w+1 " '^i,w'' (4.30) (3) calculate new values of ^ from (4.20) using SOR until convergence is attained; assign these values as n^ . and then damp them by using iÂ»0 1,0 w, (n . . n . . ) (4.31) (4) calculate new values of ij; from (4.19), again usin^ SOR, until convergence is attained; assign these values as ^^ ^ and then relax them by using ip. . = ^ . .+w^(i^. . ^ . .) (4.32) (5) repeat steps (2) to (4) until both n and ^ converge .
PAGE 55
40 Notice that in steps (3) and (4) above, iteration is used to find convergent a and ^ values for each set of assumed (or revised) and updated n and ^ 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 n and ^^) 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 w^ = 0.17 yield convergent solutions; w^ and Wo 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 ij; and a calculated for low Reynoldsnumber cases are used as the initial guess for high Reynoldsnumber cases. Thus, for example, the output of ^ and fi 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 ^i,w "i,w < 10
PAGE 56
41 which also enables the following conditions to be met [15]: ,n+1 ,n 1,J 1Â»J < 10 and 1,0 i,J < 10 4 (ii) The convergence for inner iteration of Eqs. (4.19) and (4.20) is set at ,k+1 k 1,0 1,0 < 10 3 and 1,0 1,0 < 10 3 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,
PAGE 57
42 Set initial guess fpr ipl j and ^. j Update wall vorticity Q^^ t.; Solve for vorticity fi. .Â• by SOR un yes(^.^j) Damped by .i;^U^J^j+w^(^i,j"?,j) jk Solve for stream function ij^. ^^ by SOR noyes(i..^j) Damped by ^^^5=^;^j+W2(^,.^j4^';,j) I noFigure 7 Block Diagram Showing the Numerical Procedure Used to Solve Stream Function and Vorticity
PAGE 58
43 tC23 (z^/z^).l ^Qi,2 = C24e,^^^2^C25G,_^^2^C260,^3 . Â•^^7 (/;,\, ^1,1CQ^. (4.33) and one point awiay. from theouterwall t^23X67^]ei,2 = C24e.^^ ^2^0250. _^^2^C26e. ^3 ^^27j^f^0i_,,,.COe^^ (4.34) It is also noted that, before applying these equations to find 0j_ 2> Â®i1 1 must. be updated by using (4.28) and (4.29). The relaxation parameter Wq is taken to be 0.7 for all cases studied. A fast convergence results if the iteration proceeds from i = 1 to imax'^^^^ is following the flow direction, and = Jiiiax ^Â° '' ' that is, from the inner wall to the outer wall. Â•
PAGE 59
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 e.aluated all the way from in^et to exit. xhe 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. 44
PAGE 60
45 The NavierStokes equations can be used together with the continuity equation and vorticity to derive the diraensionless pressure gradients at the wall as (Appendix D) 8z Re^r 8r^ ^^' ' ^ and where 8r ~ Re dz ^^ 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, P,^ = (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 diraensionless 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] Â«^ '"" = 1 (5.5) (8/C*) dz
PAGE 61
46 In this equation, C^ = r2[(i+r^2) _.JlzÂ£l_] (5.5) Â° ln(1/r^) and p^, is diraensionless, p^, = p^,/(pu ); Re is the Reynolds number; r in Eq. (5.6) is a radius ratio, ^ = '^io/'^oFor a Poiseuille flow under a fully developed condition, Eq. (5.5) is changed to Re ^Pfd _ f . which can be recast as ""^ *" = C* (5.8) (8/C"^) dz Comparing Eqs. (5.5) and (5.8) provides a clue for a new pressure defined as "^ P (5.9) (S/C"") With this definition, the pressure gradients at the inner wall become
PAGE 62
and 12. 8Z a/c"" r 8 r 47 (5.10) 12. 9a ar (8/C'') 8z (5.11) Correspondingly, the exit pressure condition should satisfy p" = (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 dPex ^Pfd dz dz = 1 (5.13) 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 t + + + f
PAGE 63
48 where t refers to a surface coordinate. The integrand in this equation can be expressed by means of a vector analysis as if=vp^ . t= (1;, .If a^) Â• t (5.15) where t is a unit vector tangential to the wall. Expressions for ap'^/az and 3p"^/9r have been given previously as Eqs. (5.10) and (5.11). In the transformed plane, they become (Appendix A) + . . zivn) z (rf2)^ 12_ ^ ___[__ 1 _l D D L (5.16) 3z (8/C^) r J and + . r Q, r^n IS= 1 nS L_IL (5.17) 3r (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): A A z^e^ + r_e^ t = ^ ^ ^ ^ (5.18) ;[ y The differential arc length along lines of constant n is
PAGE 64
49 dt = ^ Y d ? (5.19) Introducing Eqs. (5.15) through (5.19) into (514) gives . + + + Ap = P2 P] = ^1 (8/C^) "^"^ ' ^ " ^ ^ "Jr 1 dS (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 m m where Cy^ is a dimensionless hydraulic diameter, defined as C, = D, /r. (5.22) h h 1 Dj^ is the hydraulic diameter; '^^ in Eq. (5.21) denotes a .nixing cup temperature, defined as [18]
PAGE 65
50 / u9r dr Q = Â— Â• (5.23) I J, 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: V = %1 ^ p/z^4f ^ d(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 0jjj2 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
PAGE 66
51 P^^lS^ml pA2^2S^ni2 * 2 * pA^u^=pA2U2=pTT(r.) u^ Figure 8 A Control Volume Used to Find the Mean Temperature by Energy Balance
PAGE 67
52 '^u2%^*^k\l<\^' <5'=^' The scheme developed above works well in the pipe regions where the flow sections are uniform. At the turning tip (point h, figure 3), no hydraulic diameter is defined. The characteristic length, C^, 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 / Nu dA / Nu^ r dz Nu = ; I, = r^ (5.26) ra J dA J 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.
PAGE 68
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 plotroutine package. VI. 1 Velocity Field Velocity fields for three Reynolds numbers are presented in figure 9Enlarged 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 g 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 53
PAGE 69
54 .111 I n ..III I It 1^1 Hill 'fl o o o o o o o o ^Â— II 0) cd u o Â«H 03 X) ft 0) Â•H P Â•H O O rÂ— I >
PAGE 70
55 / HI ! f ;ii o Â•p c Â•H o cu W) c Â•rH c tl Ei CO CD rÂ— I 0) Â•H >, p Â•H O O rH 0) > 0)
PAGE 71
56 m I l5?^^t^ / ' ii^'v "^^' ^ I ' ,1 o o Ln II 0) p fe.
PAGE 72
57 OS o Â«H P C Â•H O a. bO c Â•H fl U EH t4 CO 0) u bO
PAGE 73
58 noted that, in making these plots, the magnitude of the velocity has been normalized by the maximum velocity found for each Reynolds number tesifed. 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 terras 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 J.
PAGE 74
59 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 tov.ard 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 ( dimensionless) in figure 13. The turning point of the inner pipe is located at 0.5; the
PAGE 75
Table 1 Vortex Sizes for Re=100, 500, and 1000 60 Reynolds Number
PAGE 76
61
PAGE 77
62 normalized overall length is 1. Again, the vorticity is diraensionless , normalized by the maximum vorticity value (785.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 noslip 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.
PAGE 78
63 Table 2 InnerWall Vorticity Data Near the Turning Point iN. ^^ ^
PAGE 79
64 VI. 3 Pressure on Inner Wall The pressure distribution on the inner wall is plotted for the three Reynolds numbers 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
PAGE 80
65 CO 0) c c =3 U
PAGE 81
66 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 z^ by means of the velocity at z=z+a2However, in the present investigation, no linearization is made and the axial diffusion terras 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 Iso 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 180degree 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 .
PAGE 82
67 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, dpgj,/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 9^=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
PAGE 83
Table 3 Dlniensionless Pressure Gradients at Exit 68 Re
PAGE 84
nil! 69 mmih
PAGE 85
70 3j\
PAGE 86
71 M o o ! I fe
PAGE 87
72 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.
PAGE 88
73 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 13; 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 3658 [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 Nu^ 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
PAGE 89
74 o o in fa
PAGE 90
75 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 studydeviate 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.5^. 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
PAGE 91
76 o o
PAGE 92
77 o o
PAGE 93
78 Table 4 Comparison of the Computed Local Nusselt Number at the Exit with Lundberg's Data + r
PAGE 94
79 10' Nu. 10^ I10' 3x10 3x10^ Re Figure 21 A Plot of Mean Nusselt Number
PAGE 95
80 the results yields Nu = 1.15RBÂ°528pr02S8 (6.1) m 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 r\line is fixed at 163; the number of grid points along the constant E,line is tested for 21, 41, and 61. This results in lines of constant c for these grid systems having the
PAGE 96
81 Table 5 Comparison of Mean Nusselt Number Between SPRF Heat Exchanger and Heating in a Circular Pipe Pr
PAGE 97
82 same contours in the physical plane but with different grid sizes . Figure 22 shows the comprtited stream functions along constant 5lines 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 ^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 (153x61). 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.
PAGE 98
83 0.5 0.4 0.2 0.1 i=50(about half way (163x61)inside the inner pipe) 63x61) ]___i I . I I L. 0. .1 .2 .3 .4 .5 .6 .7 .8 .9 1 0.5 ' 1 ' rÂ— 1 1 1 1 ' 1 Â— Â•(163x61 )n 0.4 [. i = 152(about half way inside the annular region) 0.3 0.2 I 0.1 0.0 (163x21) 0. (163x41) 3 .4 .5 .6 .7 .8 (^^min^/^Vx^nin) Figure 22 Comparison of the Computed Stream Functions Along Constant Cline for Different Grid Systems
PAGE 99
84 Table 6 CPU Time Required for Computation of Stream Functions Given in Figure 22 \^^^^^ Re \ CPU *""^^ \ time Grids \
PAGE 100
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 fiov7 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 differenced, with the firstorder differential 85
PAGE 101
86 terms of velocity and temperature expressed in a firstorder approximation. Solutions are obtained for Pr=0.7, 20, and 100, and Re=100, 500,* and 1000. Numerical data show a marked increase in velocity as flow makes a 180degree turn in the end cap. Here a rapid drop of pressure is also noticed. Once the flow enters annulus, the main flow is separated from the inner wall, and a recirculation is found near the wall. This results in a drop in the heat transfer coefficient in this recirculation region. Moving further downstream, the flow reattaches to the wall, causing a slight rise in the heat transfer. The mean heat transfer coefficient for the heat exchanger is found to be proportional to the Reynolds number and Prandtl number. The numerical results are correlated to develop an empirical equation. For a check of the accuracy in the solution, the pressure gradient and the Nusselt number at the exit of the exchanger are also evaluated and compared with the theory (analytical equation) and practice (literature data). The agreement is satisfactory. The singlepass, returnflow heat exchanger is found to have a mean Nusselt number that is about 1.57 to 1.37 times that for heat transfer in a singular pipe under comparable working conditions (Pr=0.7, and Re=500, and 1000). The exchanger is thus shown to be effective in heat transfer.
PAGE 102
87 Based on the study made in this dissertation, directions of further research are charted as follows. 1. There are a wide range of variations in the conditions used which the SPRF exchanger can be operated. They include (i) change in the flow direction, (ii) change in the system geometry, particularly that gap size from the inner pipe to the end cap, (iii) change in thermal boundary conditions on both walls, and (iv) introducing minor modifications in the construction of the exchanger, such as eccentricity in the pipe layout, flow perturbation devices, etc. They provide a wide range of conditions to be tested in the future. 2. Because of the complexity of the geometry of the exchanger, it provides an opportunity for testing various numerical methods, particularly those approximation methods in the finitedifference formulation of the boundary conditions. 3. A study is not complete without experiment. It is highly desirable that hardware be built to test the results obtained in the numerical solution. A limited experimental program serves well to validate the numerical solution; until then the numerical solution can be used with confidence for acquiring data for various operating conditions.
PAGE 103
APPENDIX A DERIVATION OF TRANSFORMATION RELATIONS, GOVERNING EQUATIONS, AND BOUNDARY CONDITIONS IN THE TRANSFORMED PLANE This appendix gives the derivation of the transformation relations, governing equations, and boundary conditions in the physical and computational planes when elliptic partial differential equations are used to generate the computational grid system. Referring to figure A.1, the independent variables in the physical plane are z(C,n) and r(C,n), while those in the computational plane are ?(z,r) and n(z,r). In the computational plane, the domain of interest is rectangular in shape and grids in this plane are uniformly spaced (i.e., square grid). Since all the computations are to be made in the transformed domain, the independent variables in the governing equations and boundary conditions must be changed from the physical domain to the transformed domain. Consider the vector z(5,n) r(5,n) (A.1) which may be differentiated to give 88
PAGE 104
89 $r _^ A 0) c CO a, Â•a s u o CQ c CO EH o p CO O Â•H CQ >! Â£: as o u o Â•H P CO a L< O Â«H CQ C CO u EH c Â•H o X! CO p CO a Si o CO fa
PAGE 105
dz dr z^ z r^ r d5 dn 90 (A. 2) Cramer's rule may be used to solve for d? and dn as d? dn . n Â— r 5 5 dz dr (A. 3) where J is the transformation Jacobian J = Similarly, 9(2, r) ^ 3C5,n) z_ z C n r _ r 5 n dS dn ?z ?r ^z ^r dz dr (A. 4) The following identities may be established by comparing (A. 3) and (A. 4) C = r /J z n 5r ^Â« "z = 'i'' (A. 5) and ^r = \/^
PAGE 106
91 The original partial differentials in the physical plane may be transformed to the computational plane by using the chain rule for partial derivatives. It follows that dz ^z 35 ^z 8ti = (r ^ r ^)/J (A. 6) 3r ^r 95 ^r 3n = (Z ^ + z^ ^)/J (A. 7) Ti 35 5 9n 2 1 = ^(r 2^^2 3Z ^z 35 ' "z 3ti' _ _3 _3 /Â•_3 \ _3 _3 /_3 s ^ZZ 35 ^z 3Z 35 '^zz 3ri ^z 3Z 3n = 5 ^+ 5 [ ( 5 ^+ n ^)^] '^zz 35 ^z^'^z 35 ^z 3n 8C o 2 Â„ 2 ^zz 35 '^z 2 ^zz 3n ^z 2 35 3ri + 2 5 n ^ (A. 8) ^z"z 358n and
PAGE 107
92 3^ ^ 9 ? 8^ a 2 8^ r r 95 on Combining (A. 8) and (A. 9) gives the Laplace operator ^^ ^^ = (C? + C?) ^ + 2(? n^ + Â£ n ) ^ Â„ 2 ^ 2 ' z r' .^2 '"z z "r r' a^Sn 9z or 35 ^2 2s 3^ / ^ r ^ 3 ^ ^^z ^^^72" ^^zz " ^rr^ T? 3ri + (n + n ) 4(A. 10) zz rr 3n in which the coefficients for partial differentials on the righthand side may be expressed in terms of J as and ?2 + e^ = (z^ + r2)/j2 = a/j2 z r ^ n n 2 2 C n + 5 n = (r^r + zz )/J = b/J z z r r C n 5 n 2 2,2 2s/t2 /,2 riÂ„ + nÂ„ = (z^ + r )/J = Y/J where 2 2 a = z + r n n = r^r^ + z_z n ? n
PAGE 108
2 2 Then, 93 2 2 2 9 ^ i 1 ^ J oÂ„ 9 + ?r = 7(> a 5 2B Â— .2 , 2 ~ J 9z 3r 35 953n 9ti ^ ^zz ^rr' dE, 'zz ^rr^ an (A. 11 ) If Laplace equations C +5 =0 zz rr (A. 12) and (A. 13) are used in the transformation, Eq. (A.11) may be simplified to be 1/8 ,Â„ af ^ 2 2 .2 '" 2 3Z 9r J 95 Y Â— y. 3 5 9ti 9ri (A. 14) If Poisson equations ?zz ^ S^r = P^^'^^ (A. 15) and ^zz ^ ^rr = ^^^'^) (A. 16) are used, Eq. (A.11) reduces to
PAGE 109
94 2 2 2 2 2 \ a Â— 2 e + Y ^) p P ~ P 2 ? 3z 3r J 3C 859n 3ti + P TT + Q T^ (A. 17) Those relations derived above can be used to derive transformation equations listed below. (1) From Eq. (A. 14), the Laplace transformation equations in the computational plane become az^^ 2B2^ + yz =0 (A. 18) C? 5n nn and ar^^ 26r^ + yr =0 (A. 19) Â£,i 5ri nn (2) From Eq. (A. 17), the Poisson transformation equations in the computational plane become az 2bz + vZ = J^(Pz + Qz ) (A. 20) 5C 5n nn 5 n and ar 2Br + ^r = J^(Pr^ + Qr ) (A. 21) CS 5n nn 5 n From these two equations, the source terms P and Q can be derived as P = (z Dr r Dz)/J^ (A. 22) n n and
PAGE 110
95 Q = (r^Dz z^Dr)/J^ (A. 23) where Dr = ar^^ 26r^^ + Yr^^ and Dz = az^^ 26z^^ + Yz^^, (3) With the use of Eqs. (A. 7) and (A. 17), the vorticity definition equation (2.2) in the computational plane takes the following form Jz Jz = J^r^ (A. 24) (4) With the use of Eqs. (A. 6) and (A. 7), the velocity equations (2.3) and (2.4) in the computational plane become ^ = jF ^^i\ ^nh^ ^^25) and ^ JF ^H\ ^^) ^^26) (5) With the use of Eqs. (A. 6), (A. 7), and (A. 17), the vorticity transport equation (2.1) in the computational plane becomes 2 Jz Jz aa^^ 26fi^^ + yn^^ ^ a ^ (j2p _ ^)^^ + (j^q + ^)^^ r = ^ [i\^, ^,\) ^f (r,>^5 r^^^)] (A. 27)
PAGE 111
96 (6) The energy equation (2.5) in the computational plane is also derived as Â«^5 230^, . ye^, . (J2p ^)e^ . (j^q . ^)e^ Derivation of the governing equations in the computational plane is now complete. To derive boundary conditions, one must first derive the expressions for the normal and tangential vectors in the computational plane. For a scalar function F in the physical plane, the gradient of F is VF = F e + F e (A. 29) z z r r This equation can be recast, with the use of (A. 6) and (A. 7), as VF =Â• 7 ^^\h 'i'/'z * (/? ^ ^/n>Â«ri (*30) On a line of constant n, n^ = and n^ = 1 ; Eq. (A. 30) reduces to Â„ . 1. (,^e^ . z^e^) (A.31)
PAGE 112
97 Similarly, on a line of constant 5, 5^ =1, 5^ =0, and '^ =T^\K \K^ (A32) With reference to figure A. 2, the unit normal vector to lines of constant n can be expressed as n^^) ^^ ^ (A.33) The unit normal vector to lines of constant 5 is ^ / \ r e z e The unit tangent vectors can be derived with the use of cross products as follows: (n) :(n)v : _ ^g^z ^ ^K^r t^^^ = n^^^X e, = ^^ ^^ (A. 35) and The directional derivatives of F along normal and tangential directions can be expressed as
PAGE 113
~^^ms^' < QJ P c CO Â•p CO c o u u o CQ c o m u o p o CO s ti o 3 T3 C CO P C
PAGE 114
99 8n 9F ;(Ti) 3F _ (n) :i^ 8F ^(Z n 7T? on and 9F ;(?) VF = 4 ^^ r^ VF rr r^ VF (A. 37) (A. 38) (A. 39) (A. 40) The arc length dt along any contour r in the (z,r^ plane is (see figure A. 3) dt (r ^ = \ (dz)2 + (dr)2 In the transformed plane, dt (r) J^(dC) ^ + BdCdn + a(dn)^ (A. 41) (A. 42) where (A. 2) has been used for dz and dr. If r is a line of constant n, Eq. (A. 42) reduces to ,(n) dt' ' = \^( d? Similarly, if r is a line of constant C, then (A. 43)
PAGE 115
100 Figure A. 3 Schematic Showing the Arc Length Along Contour r
PAGE 116
101 Boundary conditions will now be derived in the transformed plane. Note that all conditions of the first kind can be used as is; no transformation is necessary of them. As for the other conditions, notice the fact that, in both inlet and exit domains B and C, grids are orthogonal in both original and transformed coordinates. So, metrics such as r_, z , r^.., and z all vanish. Furthermore, terms such as r^ and r ^ also vanish because 5n n5 r = ^(2^) = ^(^^) = r (A. 45) and r =0. With these in mind, consider Eq. (2.7a) first. It is clear that this condition is equivalent to u =v =;h =n =0 (A. 45) Next consider Eq. (2.7b). To express it in the transformed plane requires use of Eq. (A. 8) in which the last three partial differentials all vanish because n^ = 0* '^^^ remaining two terms on the righthand side of this equation may be transformed by use of Eqs. (A. 5) and (A. 6). It follows that
PAGE 117
102 Q =^0^ (A. 47) Along the axis of symmetry, that condition Uj, = is changed to z u z u^ = (A. 48) using (A. 7). Similarly, e = becomes z,e z e^ = (A. 49) At the outer pipe wall, e^ = is changed to y9 ee, = (A. 50) because of (A. 37). The vorticity boundary conditions at the inner and outer pipe walls can be derived using a lengthy procedure as follows. At the wall, u = v = 0; Eqs. (43) and (4.4) become z i> z J; =0 (A. 51 ) and '^\ r^*^ = (A. 52) These walls may be regarded as streamlines along which
PAGE 118
103 4^=0 (A. 53) everywhere . Using this relation in (A. 51) yields l*^ = (A. 54) which is also valid everywhere along the wall. Differentiating Eq. (4.3) with respect to 5 gives (>A.55J (Jr)2 where the second term in the numerator vanishes because of (A. 53) and (A. 54). It follows that u ^ r(z ^ + z ij; z i> z J> ) (A. 56) Along the wall, u^. = 0. Also notice that (A. 53) can be used to write ^l>^^ = (A. 57) Then from Eq. (A. 56). 4/, = (A. 58)
PAGE 119
104 Equations (A. 53), (A. 54), (A. 57), and (A. 58) may now be used in (4.2) to write ^ at the wall as n = I^^^ (A. 59) w j2^ nn
PAGE 120
APPENDIX B DERIVATION AND USE OF EQUATION (3.6) This appendix gives the details on how Eq. (3.6) is derived and used to locate grid lines in z direction. It is necessary to locate the grid lines in figure (B.1) so that Zq corresponds to index ig, zj to i^, Zj_ to ij^,ZÂ£_^ to if1 , and finally z^ to i^. In particular, because of the z^ function to be chosen later, it is necessary that the following conditions be met for i. = i^ (B.1) z^ = z^_^ for i^ = i^1 (B.2) z. = zÂ„ for i. = iÂ„ (B.3) if if In order to meet the other conditions mentioned, but not formulated above, z^ is chosen to be a cubic polynomial of the form z. = Zq + a.'?. + a^'V^ + a^'V^ (B.4) where i . i^ 1 1^ 1q 105
PAGE 121
106 4N < 1 Â—JÂ— 44CQ o c Â•H X 0) c in 0) bO Â•H
PAGE 122
107 There are three unknown coefficients (a>, a2, a^) in Eq. (B.4) and there are three conditions ((B.1), (B.2), (B.3)) to solve them. These coefficients are found to be a. = z^ Zq a2 a^ (B.5) a^ = [Az^ h(z^ Zq) a^(h^ h)]/(h^ h) (B.6) a = [Az^ + Az^ 2h(z^ ZQ)]/(h 3h^ + 2h^) (B.7) where AZ^ = Z^ Zq AZ^ = 2f ^t.^ h = (i^ Xq)' In practice, values of Iq, i^, Zq, z^, Az^ and Az^ must be preassigned. With suitable values of Az^ and Az^ (to be tested), Eq. (B.4) yields a monotonous function of increasing grid sizes from Zq to z^.
PAGE 123
APPENDIX C DERIVATION OF FINITE DIFFERENCE EQUATIONS Transformation equations (3.3) and (3.4) and governing equations (4.1) and (4.5) will be changed to finitedifference forms in this appendix. For the ease of formulation of these equations and later computation in the transformed plane, the grids in the computational plane are taken to be squares, i .e . , AC = An = 1 . In the indexing scheme, i refers to 5 and j, n axis. Defining F = F(c,n) as a twice dif ferentiable function of 5 and n, the finite difference representations for the first and second derivatives of F at an interior point (i,j) may be expressed, by means of a secondorder central difference scheme, as follows: (^c'l.a "<^.1, 0^1. a>/2 '^Â•^' (^Â«)i,j = ^1*1,0 2^,0 *^ii,d ( = 3) (F ) = (F F + F ^ CTi^i,j ^ i + 1,j + 1 1 + 1, j1 ^ 11, j1 ^1, + 1^/4 (C.4) and 108
PAGE 124
(^.>l,a = "i,J*i ""i,J * 'i,3i 109 (C.5) The firstorder backward and forward differences of F and F are given respectively as and (F ). . = F. , . F. . (C.8) ^ c'i,a 1+1, J i,j (F ). . = F. . , F. . (C.9) n 1,3 1,0+1 1,0 To facilitate representation of partial differentials in computer programs, the changes of F in both g and ri directions are written compactly as F = F. , . F. , . (CIO) 5 1+1, J 11,0 F = F. . , F. . . (C.11) and ^^Ti ^ ^1 + 1,0 + 1 " ^i + 1,01 " ^i1,d1 Then, the transformation parameters a, 3, y, J, PÂ» 3.nd become
PAGE 125
110 and where and a = (i2 , ?2)/4 e = (5^i^ + ?g?n)/4* Y = (5^ + ?^)/4 J = (^c^n ^ ^n^?)/^ P = (i^Dr ?,Dz)/(2J^) = (r^Dz z^Dr)/(2J^) Dr = iv^^).^. 26(r,,),^. . Y(r^^),^. Dz = <^(z,,)i,j 23(z^^),^. . T(z^^),^. The transformation equations (3.3) and (3.4) take the following forms in the computational plane az. . . 2(a+y)z. . + az . . 11, J i,J 1+1, J Â•"^1,0*1 * ^i,ai* *2^5n 'ciJ) and cir. . . 2(ci + y)r. . + "r. . . 11,0 iÂ»a 1+1, J ''('l.a.l * 'i.3^^ * 2 ?5n (C^t'
PAGE 126
111 These equations are to be used to map the grid points from the computational plane to the physical plane. With the finite difference equations given above, the governing equations in Chapter IV can be changed to finitedifference forms as follows: Vorticity transport equation: [2(a+T) + ^ + Re(Â— ^ ^^ 1ij; )J 4p2 c ^^2 n i,j 1 / .2 Jz 1 / .2, J: J~ Jz R^ 4F ^^5 ' ^^ 4F ^c% " 2 ~\r. (C.15) Vorticity definition equation: U . 1
PAGE 127
112 TT^^ mrr^. CC.18) 1,0 4Jr % 4Jr ^5 t Energy equation: 2(ci+Y)e. . It is convenient to introduce the following expressions for metric coefficients CO = b/2 CI = 2(a+Y) P Jz C2 = a + (J^P + 2^)/2 C3 = a (J^P + 2^)/2 P Jz C4 = y + (J^Q 2f)/2 C5 = Y (J Q 2r^/2 C6 = J^r C7 = 2(a+Y) + J^/r^
PAGE 128
113 C8 = a + (J^P 2F^/2 C9 = a (j2p _ ^)/2 CIO = Y + (J^Q + 2f)/2 ? Jz C11 = Y (J'^Q + 2f)/2 C12 = z^/(4Jr) C13 = i^/(4Jr) C14 = r^/(4Jr) C15 = r^/(4Jr) and C16 J/(4r) C17 = (J?^)/(4r^) CIS = (J?^)/(4r^) so that Eqs. (C.15) to (C.19) can be recast compactly as [C7 + Re(C17K C18^ )]fi. . ReCl6(? n^ lb n ) + COf^, n E, ' E, n Sn (C.20)
PAGE 129
114 + C6j2. . + CO^^ Â• (C.21) u = C12^ C134i, (C.22) V = C14^^ C15;^^ (C.23) and PeClbC^j; 0^ + ij; i ) + COe^ (C.24) Rearranging the righthand sides of Eqs. (C.20) and (C.24) yield [C7 + RaiCM^, C18^ )]q. . = (C8 ReClS^;^)^. . + (C9 + ReCl6i )n. , , + (CIO + ReC^bi,^)^^^^^^ + (C11 ReCIbi )fii ^_i + C05^^ (C.25) and Cle.^ . = (C8 PeCl6;^^)o.,,^j . (C9 . PeCl6^^ ) e^., ^ . + (C10 + PeCl6;i; )0 + (C11 PeClbij; )0. , , + COi^^ (C.26) Equations (C.21) to (C.23) and (C.25) and (C.2b) are the
PAGE 130
115 centraldifference versions of the original governing partial differentialequations. The numerical solution of Eqs. (C.25) and (C.26) maypresent a problem because it does not always converge. In order to rectify this problem, one way to do it is to make the diagonal dominant in the matrix expressions for these equations. Notice that the upwind method cannot be directly used here because the convection terms in those transport equations have been simplified to the extent that their original formats have lost. Use is thus made of Takemitsu's method which is slightly modified as presented below. There are four terms on the righthand side of Eq. (C.25) that can be regrouped to relate convection as CV = Re[ C^6^^{n.. . ^. . ) + Cl6^,(n. . . fi. . .)] (C.27) Introducing U = C16.)^ , V = C16^^ (C.28) permits rewriting Eq. (C.27) more compactly as CV , RetU(Â«.^,_. a..^_j) V(a._.^^ a._._,)] (c.29)
PAGE 131
116 Notice that C16 is always positive, while ^ and ^^ may change sign, depending on the local stream function. The terms in Eq. (C.29) have been written in central difference forms and they must be changed to either forward or backward difference to assure a convergent solution. Experience shows its choice hinges on the sign of U and V; a forward difference must be used if either of them is negative and vice versa. To meet these requirements, Eq. (C.29) is written as i,j + 1 i,j i,j i,jV' (C.30) which reduces to the following four situations (i) U and V are both positive: CV = 2Re[\]{n. . n. ^ .) V(n. . n. . J] (C.31) iij 11, J i,J i,J1 (ii) U is positive and V is negative CV = 2Re[U(n. . n. , .) V(f^. . . n. .)] (C.32) (iii) U is negative and V is positive CV = 2Re[\]{a. . . n. .) V(n. . ^. . ,)] (C.33) 1+1,0 1,0 1,0 1,01
PAGE 132
117 (iv) U and V are both negative: CV = 2Re[U(fi. . . ^. .) V(j^. . . n. .)] (C.34) 1+1, J ijj iÂ»J"'"' 'to The general form of CV, Eq. (C.30), is now substituted in Eq. (C.25) and with some regrouping of terms changes the latter to {C7 + Re[C17?, C184i + 2(u + v)]}q. . ^ 5 n 1 I I I 1 > = [C8 Re(U U )]a. , . + [C9 + Re(U + U ) ] n . . . II 1+1, J II '~i,j + [CIO Re(V V )]n. . . + [C11 + Re(V + I I 1 , J + I + CO^, V hjfi 5ti i,d1 (C.35) Following the similar approach, a general form of the energy equation is derived as [CI + 2Pe( U + V ) ]0. . II II 1 , J = [C8 Pe(U u)]0. , , + [C9 + Pe(U + U hje l 1+1, J II LI, J + [C10 Pe(V I v)]0. . , + [C11 + Pe (V + v)]0. . . II i,j+' i,j~i + CO0 5n (C.36) These equations will lead to a convergent solution. In closing this appendix, a short note is included to show that the scneme developed above is rooted in the upwind scheme mentioned earlier. According to this scheme.
PAGE 133
118 it is necessary to consider the flow retains its past identity when it is migrated to a new location. The backward and forward represei^tations of ^ discussed above indeed come from this origin. In fact, it can be shown that the sign of U and V are related to that of u and v, which can be verified by using those equations in domain B,
PAGE 134
APPENDIX D DERIVATION OF PRESSURE GRADIENT EQUATIONS The expressions for pressure gradients in z and r directions at the wall will be derived in this appendix. In the derivation, the NavierStokes equations will be used together with the continuity equation to derive the pressure gradient expressions. Finally, the vorticity definition is introduced to reduce the differential order in the final expressions. For axisyraraetrical flow of an incompressible medium in a cylindrical pipe or annulus, the nondimensional continuity equation and momentum equations are given as follows: 8z 3r r 3u au 9p 1 , 9^u 9^u 1 9u. ,^ ^. and 9v 9v 3d 1 ,9^v 9^v 1 9v v^ .^ ^. dz 3r 3r Re g^2 g^2 r 3r ^2 where X^o P = P /(P^m ) 119
PAGE 135
120 Equation (D.2) is used to write the pressure gradient In z direction at the wall as L. 1 (!%,!%. 1^) (D.4) The continuity equation is differentiated with respect to z 8^U 8 V 1 3v ^ 2 3z9r r 3z dZ and substituted into (D.4) to give (D.5) 3p 1 ( 9^u a^v J_ in _ 1 iv_^ n^ K ^ 3z ~ Re \ 2 " az8r "^ r 3r r ^z^ Ku.o) dr Recall that the vorticity is defined as = =lf17 (D7) Differentiating this equation with respect to r gives ar 3z3r " , 2 ^'^Â•^^ 3r Introducing these two equations into (D.6) gives the final expression for the pressure gradient in z direction at the wall
PAGE 136
121 3z Re ^r ar^ ^^Â•y'' The pressure gradient in r direction at the wall may be derived by following a similar procedure. Equation (D.3) will be used instead, and this equation is simplified by using the condition that u = v = at the wall. It follows that 3r Re \ 2 2 ^ r ar" ^^* '^^ 8z 3r Equation (D.I) is differentiated with respect to r a^u a^v^ii^.v ^^^^^^ ar9z g^2 r ar ^2 and Eq. (D.7) is differentiated with respect to z to give 2 2 in _ 9 V a u , . az g^2 " araz ^^'^^ and finally introducing these two equations into (D.10) gives the pressure gradient in r direction at the wall as iÂ£. _ _L iiL fn 1^ a r Re a z ^ u . I :)
PAGE 137
REFERENCES 1. Hurd, N. L., "Mean Temperature Difference in the Field or Bayonet Tube," Industrial and Engineering Chemiscry, Volume 38, Number 12, 1946, pp. 12661271. 2. Beekly, D. C. , and Mather, G. R. , "Analysis and Experimental Test of a HighPerforraance , Evacuated Tubular Collector," NASA, CR150874, 1978. 3. Moan, K. L., An Analysis of the Low Loss Evacuated Tubular Collector Using Air as the Heat Transfer Fluid, OwensIllinois, Inc., Toledo, Ohio. Unpublished report. 4. Shah, R. K., and London, A. L., Laminar Flow Forced Convection in Ducts, Advances in Heat Transfer, Supplement 1, Academic Press, Inc., New York, 1978. 5. Ehrlich, L. W., "The Numerical Solution of a NavierStokes Problem in a Stenosed Tube: A Danger in Boundary Approximations of Implicit Marching Schemes," Computer and Fluid, Volume 7, 1979, pp. 247256. 6. Gosraan, A. D., Pun, W. M. , Runchal, A. K., Spalding, D. B. , and Wolfshtein, M., Heat and Mass Transfer in Recirculating Flows, Academic Press, Inc., London, 1969. 7. Singh, S. N., "Heat Transfer by Laminar Flow in a Cylindrical Tube," Applied Scientific Research, Section A, Volume 7, 1958, pp. 325340. 8. Thompson, J. F., Thames, F. C, and Mastin, C. W., "BoundaryFitted Curvilinear Coordinate Systems for Solution of Partial Differential Equation on Field Containing Any Number of Arbitrary TwoDimensional Bodies," NASA, CR2729, 1977. 9. Nietubicz, C. J., Heavey, K. R., and Steger, J. L. , "Grid Generation Techniques for Projectile Configurations," Proceedings of 1982 Army Numerical Analysis and Computer Conference, Vicksburg, Mississippi. ARO Rept. 823, 1983. 122
PAGE 138
123 10. Sorenson, R. L., and Steger, J. L., "Simplified Clustering of Nonorthogonal Grids Generated by Elliptic Partial Differential Equations," NASA, TM73252, 1977. 11. Takeraitsu, N. , "On a FiniteDifference Approximation for the SteadyState NavierStokes Equations," Journal of Computational Physics, Volume 36, 1980, pp. 236248. 12. Roache, P. J., Computational Fluid Dynamics, Revised Printing, Herraosa Publishers, Albuquerque, New Mexico, 1976. 13. Parmentier, E. M. , and Torrance, K. E., "Kinematically Consistent Velocity Fields for Hydrodynamic Calculations in Curvilinear Coordinates," Journal of Computational Physics, Volume 19, 1975, pp. 404417. 14. Rigal, A., "Acceleration of the Convergence in Viscous Flow Computations," Journal of Computational Physics, Volume 43, 1981, pp. 177179. 15. Roache, P. J., "Semidirect Calculation of Steady Twoand ThreeDimensional Flows," in Numerical Methods in Laminar and Turbulent of the First International Conference, John Wiley & Sons, Inc., New York, 1978, pp. 1727. 16. Miyakoda, K., "Contribution to the Numerical Weather PredictionComputation with Finite Difference," Japanese Journal of Geophysics, Volume 3, 1962, pp. 75190. 17. Bird, R. B. , Stewart, W. E., and Lightfoot, E. N., Transport Phenomena, John Wiley & Sons, Inc., New York, I960. 18. Kays, W. M. , and Crawford, M. E., Convective Heat and Mass Transfer, Second Edition, McGrawHill, New York, 1980. 19. White, F. M., Viscous Fluid Flow, McGrawHill, New York, 1974. 20. Morihara, H. , and Cheng, R. T.S., "Numerical Solution of the Viscous Flow in the Entrance Region of Parallel Plates," Journal of Computational Physics, Volume 11, 1973, pp. 550572.
PAGE 139
124 21. Vrentas, J. S., Duda , J. L., and Bargeron, K. G., "Effect of Axial Diffusion of Vorticity on Flow Development in Circular Conducts: Part I. Numerical Solutions," A. I. Ch. E. Journal, Volume 12, Number 5, 1966, pp. 837844. Â• 22. Hornbeck, R.W., "Laminar Flow in the Entrance Region of a Pipe," Applied Scientific Research, Section A, Volume 13, 1964, pp. 224232. 23. Hornbeck, R. W., "An AllNumerical Method for Heat Transfer in the Inlet of a Tube," American Society of Mechanical Engineering, Paper 65WA/HT36, 1965. 24. Lundberg, R. E., Reynold, W. C. , and Kays, W. M. , "Heat Transfer with Laminar Flow in Concentric Annuli with Constant and Variable Wall Temperature and Heat Flux," NASA, TN D1972, 1963. 25. Carnahan, B. , Luther, H. A., and Wilkes, J. 0., Applied Numerical Method, John Wiley & Sons, Inc., New York, 1969.
PAGE 140
*' \'"' r BIOGRAPHICAL SKETCH SongLin Yang, was born on November 6, 1952, in Hsinchuang, Taipei, Taiwan, Republic of China. He received his early childhood education in his hometown; and his high school education in Taipei. After graduation from high school, he entered Chung Yuan University, Chungli, Taiwan, where he received his Bachelor of Science degree in mechanical engineering (1976) and Master of Science degree in applied physics (1980). On March 7, 1981, he married LiPin Hsiao and now they have a lovely son named Albert ChungPu Yang. He is a student member of the American Society of Mechanical Engineers (ASME) and the American Society of Heating, Refrigeration and Air Conditioning Engineering (ASHRAE). He is also a member of Tau Beta Pi, the National Engineering Honor Society. 125
PAGE 141
I certify that I have read this study and in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. e&~ r^ C. K. Hsieh, Chairman Professor of Mechanical Engineering I certify that I have read this study and in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. V. P. Roan Professor of Mechnical Engineering I certify that I have read this study and in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. R. K. Irey Professor of Mechanical Engineering I certify that I have read this study and in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. IX w dL hh u C. C. Hsu Professor of Engineering Science
PAGE 142
I certify that I have read this study and in ray opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Associate Prolessorof Mechanical Engineering This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School, and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. May, 1985 iLua/I L^tAt,; Dean, College of Engineering Dean for Graduate Studies and Research

