Citation
Fluid flow and heat transfer in a single-pass, return-flow heat exchanger

Material Information

Title:
Fluid flow and heat transfer in a single-pass, return-flow heat exchanger
Creator:
Yang, Song-Lin, 1952- ( Dissertant )
Hsieh, C. K. ( Thesis advisor )
Roan, V. P. ( Reviewer )
Irey, R. K. ( Reviewer )
Hsu, Chen-Chi ( Reviewer )
Ingley, H. A. ( Reviewer )
Place of Publication:
Gainesville, Fla.
Publisher:
University of Florida
Publication Date:
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 )
non-fiction ( marcgt )

Notes

Abstract:
Numerical solutions for flow and heat transfer in a single-pass, return-flow heat exchanger are given for Reynolds numbers (Re) of 100, 500, and 1000, and Prandtl numbers (Pr) of 0.7, 20, and 100. The flow is modeled as laminar, axisymraetric, incompressible, and in steady state, and the Navier-Stokes equations are expressed in terms of stream function and vorticity. Prior to the solution of the problem by means of a finite difference method, the original geometry of the system in the physical plane is transformed into a rectangular domain with square grids in a computational plane. Grid clustering techniques are also employed to reposition grid lines for good resolution in those regions where the flow and temperature fields change most rapidly. The original partial differential equations are expressed in finite-difference forms, and in order to assure convergence in the numerical iteration, the matrix coefficients in the matrix equations are reconstructed to make those diagonal terms dominant in the coefficient 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 122-124.
Additional Physical Form:
Also available on World Wide Web
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by Song-Lin Yang.

Record Information

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 non-profit 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:

PDF ( 4 MBs ) ( .pdf )

fluidflowheattra00yangrich_Page_065.txt

fluidflowheattra00yangrich_Page_094.txt

fluidflowheattra00yangrich_Page_090.txt

fluidflowheattra00yangrich_Page_127.txt

fluidflowheattra00yangrich_Page_141.txt

fluidflowheattra00yangrich_Page_087.txt

fluidflowheattra00yangrich_Page_060.txt

fluidflowheattra00yangrich_Page_106.txt

fluidflowheattra00yangrich_Page_042.txt

fluidflowheattra00yangrich_Page_061.txt

fluidflowheattra00yangrich_Page_027.txt

fluidflowheattra00yangrich_Page_091.txt

fluidflowheattra00yangrich_Page_073.txt

fluidflowheattra00yangrich_Page_081.txt

fluidflowheattra00yangrich_Page_110.txt

fluidflowheattra00yangrich_Page_008.txt

fluidflowheattra00yangrich_Page_113.txt

fluidflowheattra00yangrich_Page_131.txt

fluidflowheattra00yangrich_Page_112.txt

fluidflowheattra00yangrich_Page_007.txt

fluidflowheattra00yangrich_Page_026.txt

fluidflowheattra00yangrich_Page_037.txt

fluidflowheattra00yangrich_Page_074.txt

fluidflowheattra00yangrich_Page_132.txt

fluidflowheattra00yangrich_Page_059.txt

fluidflowheattra00yangrich_Page_092.txt

fluidflowheattra00yangrich_Page_062.txt

fluidflowheattra00yangrich_Page_040.txt

fluidflowheattra00yangrich_Page_057.txt

fluidflowheattra00yangrich_Page_064.txt

fluidflowheattra00yangrich_Page_136.txt

fluidflowheattra00yangrich_Page_034.txt

fluidflowheattra00yangrich_Page_050.txt

fluidflowheattra00yangrich_Page_088.txt

fluidflowheattra00yangrich_Page_142.txt

fluidflowheattra00yangrich_Page_025.txt

fluidflowheattra00yangrich_Page_044.txt

fluidflowheattra00yangrich_Page_076.txt

fluidflowheattra00yangrich_Page_055.txt

fluidflowheattra00yangrich_Page_086.txt

fluidflowheattra00yangrich_Page_011.txt

fluidflowheattra00yangrich_Page_068.txt

fluidflowheattra00yangrich_Page_134.txt

fluidflowheattra00yangrich_Page_022.txt

fluidflowheattra00yangrich_Page_006.txt

fluidflowheattra00yangrich_Page_036.txt

fluidflowheattra00yangrich_Page_024.txt

fluidflowheattra00yangrich_Page_108.txt

fluidflowheattra00yangrich_Page_140.txt

fluidflowheattra00yangrich_Page_102.txt

fluidflowheattra00yangrich_Page_100.txt

fluidflowheattra00yangrich_Page_049.txt

fluidflowheattra00yangrich_Page_058.txt

fluidflowheattra00yangrich_Page_056.txt

fluidflowheattra00yangrich_Page_119.txt

fluidflowheattra00yangrich_Page_067.txt

fluidflowheattra00yangrich_Page_097.txt

fluidflowheattra00yangrich_Page_001.txt

fluidflowheattra00yangrich_Page_013.txt

fluidflowheattra00yangrich_Page_115.txt

fluidflowheattra00yangrich_Page_004.txt

fluidflowheattra00yangrich_Page_096.txt

fluidflowheattra00yangrich_Page_111.txt

fluidflowheattra00yangrich_Page_082.txt

fluidflowheattra00yangrich_Page_016.txt

fluidflowheattra00yangrich_Page_019.txt

fluidflowheattra00yangrich_Page_046.txt

fluidflowheattra00yangrich_Page_031.txt

fluidflowheattra00yangrich_Page_122.txt

fluidflowheattra00yangrich_Page_133.txt

fluidflowheattra00yangrich_Page_129.txt

fluidflowheattra00yangrich_Page_085.txt

fluidflowheattra00yangrich_Page_104.txt

fluidflowheattra00yangrich_Page_029.txt

fluidflowheattra00yangrich_Page_137.txt

fluidflowheattra00yangrich_Page_099.txt

fluidflowheattra00yangrich_Page_063.txt

fluidflowheattra00yangrich_Page_109.txt

fluidflowheattra00yangrich_Page_089.txt

fluidflowheattra00yangrich_Page_069.txt

fluidflowheattra00yangrich_Page_039.txt

fluidflowheattra00yangrich_Page_126.txt

fluidflowheattra00yangrich_Page_035.txt

fluidflowheattra00yangrich_Page_033.txt

fluidflowheattra00yangrich_Page_051.txt

fluidflowheattra00yangrich_Page_010.txt

fluidflowheattra00yangrich_Page_054.txt

fluidflowheattra00yangrich_Page_095.txt

fluidflowheattra00yangrich_Page_015.txt

fluidflowheattra00yangrich_Page_030.txt

fluidflowheattra00yangrich_Page_139.txt

fluidflowheattra00yangrich_Page_101.txt

fluidflowheattra00yangrich_Page_018.txt

fluidflowheattra00yangrich_pdf.txt

fluidflowheattra00yangrich_Page_028.txt

fluidflowheattra00yangrich_Page_098.txt

fluidflowheattra00yangrich_Page_075.txt

fluidflowheattra00yangrich_Page_012.txt

fluidflowheattra00yangrich_Page_043.txt

fluidflowheattra00yangrich_Page_070.txt

fluidflowheattra00yangrich_Page_118.txt

fluidflowheattra00yangrich_Page_020.txt

fluidflowheattra00yangrich_Page_121.txt

fluidflowheattra00yangrich_Page_080.txt

fluidflowheattra00yangrich_Page_093.txt

fluidflowheattra00yangrich_Page_038.txt

fluidflowheattra00yangrich_Page_045.txt

fluidflowheattra00yangrich_Page_032.txt

fluidflowheattra00yangrich_Page_003.txt

fluidflowheattra00yangrich_Page_123.txt

fluidflowheattra00yangrich_Page_135.txt

fluidflowheattra00yangrich_Page_083.txt

fluidflowheattra00yangrich_Page_053.txt

fluidflowheattra00yangrich_Page_078.txt

fluidflowheattra00yangrich_Page_116.txt

fluidflowheattra00yangrich_Page_107.txt

fluidflowheattra00yangrich_Page_084.txt

fluidflowheattra00yangrich_Page_052.txt

fluidflowheattra00yangrich_Page_048.txt

fluidflowheattra00yangrich_Page_072.txt

fluidflowheattra00yangrich_Page_041.txt

fluidflowheattra00yangrich_Page_125.txt

fluidflowheattra00yangrich_Page_130.txt

fluidflowheattra00yangrich_Page_128.txt

fluidflowheattra00yangrich_Page_117.txt

fluidflowheattra00yangrich_Page_105.txt

fluidflowheattra00yangrich_Page_017.txt

fluidflowheattra00yangrich_Page_021.txt

fluidflowheattra00yangrich_Page_014.txt

fluidflowheattra00yangrich_Page_138.txt

fluidflowheattra00yangrich_Page_120.txt

fluidflowheattra00yangrich_Page_077.txt

fluidflowheattra00yangrich_Page_023.txt

fluidflowheattra00yangrich_Page_047.txt

fluidflowheattra00yangrich_Page_114.txt

fluidflowheattra00yangrich_Page_079.txt

fluidflowheattra00yangrich_Page_009.txt

fluidflowheattra00yangrich_Page_124.txt

fluidflowheattra00yangrich_Page_103.txt

fluidflowheattra00yangrich_Page_005.txt

fluidflowheattra00yangrich_Page_066.txt

fluidflowheattra00yangrich_Page_071.txt


Full Text











FLUID FLOW AND HEAT TRANSFER
IN A SINGLE-PASS, RETURN-FLOW HEAT EXCHANGER






BY






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

patience and encouragement made his graduate study a

beautiful memory.

Finally, the author wishes to acknowledge the

Interactive Graphics Laboratory (IGL) and the Northeast

Regional Data Center (NERDC) of the State University System

of Florida for allowing him to utilize the computing

facilities on campus.













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 Inner-Wall Vorticity Data Near the Turning
Point .............................................65

3 Dimensionless Pressure Gradients at Exit.......... 68

4 Comparison of the Computed Nusselt Number at
the Exit with Lundberg's Data .....................78

5 Comparison of Mean Nusselt Number Between SPRF
Heat Exchanger and Heating in a Circular Pipe.....81

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













LIST OF FIGURES


Figure Page

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

2 Specifications of the Single-Pass,
Return-Flow Heat Exchanger .......................7

5 A Schematic Showing Mapping from the Physical
Plane to the Computational Plane................13

4 Computed Grid Before and After Clustering in
Domain A ........................................19

5 Grid Detail Near the Inner Pipe Wall and the
Turning Point ...................................24

6 Complete Grid System............................26

7 Block Diagram Showing the Numerical Procedure
Used to Solve Stream Function and Vorticity.....42

8 A Control Volume Used to Find the Mean
Temperature by Energy Balance ...................51

9 Velocity Fields for Re=100, 500, and 1000.......54

10 Velocity Field Near Turning Point for Re=100....55

11 Velocity Field Near Turning Point for Re=500....56

12 Velocity Field Near Turning Point for Re=1000...57

13 Vorticity on the Inner Wall.....................61

14 Pressure Distribution on the Inner Wall.........65

15 Isotherms for Pr=0.7 and Re=100, 500, and
1000 ............................................69

16 Isotherms for Pr=20 and Re=100, 500, and
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 C-lines for Different
Grid Systems..................................... 83

A.1 Schematic Showing Transformation from Physical
to Transformed Plane ............................89

A.2 Unit Tangent and Normal Vectors on Lines of
Constant C and n............................ ......98

A.3 Schematic Showing the Arc Length Along
Contour r ......................................100

B.1 Indexing z Axis for Eq. (B.4)..................106


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)

C1-C18 metric coefficients defined in Appendix C

C19-C28 coefficients defined in Chapter IV

Dh hydraulic diameter

rep,e,ez unit vectors in the cylindrical coordinate
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, (Tw-T)/(Tw-Te)

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 SINGLE-PASS, RETURN-FLOW HEAT EXCHANGER

BY

SONG-LIN YANG

May, 1985


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

Numerical solutions for flow and heat transfer in a

single-pass, return-flow heat exchanger are given for

Reynolds numbers (Re) of 100, 500, and 1000, and Prandtl

numbers (Pr) of 0.7, 20, and 100. The flow is modeled as

laminar, axisymmetric, incompressible, and in steady state,

and the Navier-Stokes equations are expressed in terms of

stream function and vorticity.

Prior to the solution of the problem by means of a

finite difference method, the original geometry of the

system in the physical plane is transformed into a

rectangular domain with square grids in a computational

plane. Grid clustering techniques are also employed to

reposition grid lines for good resolution in those regions

where the flow and temperature fields change most


xiii








rapidly. The original partial differential equations are

expressed in finite-difference forms, and in order to

assure convergence in the numerical iteration, the matrix

coefficients in the matrix equations are reconstructed to

make those diagonal terms dominant in the coefficient

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 single-pass, return-flow (SPRF) heat exchanger

shown in figure 1 consists of a circular pipe inserted

inside an outer pipe. The outer pipe is sealed at one end

so that fluid entering the inner pipe may be deflected into

the annulus formed between the inner and the outer pipe of

the exchanger. Either the inner or the outer pipe wall, or

both, may be heated, and the heat input may be simulated by

different conditions at these walls.

The heat exchanger described above is structurally

similar to the bayonet tube studies by Hurd [1] in 1946;

however, they serve different purposes. In the bayonet

tube, the double-pipe construction is essentially used for

transporting fluid. Hence, the entering and exiting fluids

exchange heat with each other via the heat conducting inner

wall, and whether heat is transferred in or out of the

system depends on the temperature difference between the

fluid in the annulus and the medium surrounding it at the

shell side. This is not so, however, for the heat

exchanger, in which the inner pipe wall may be heated.

Then, the entering and exiting fluids both exchange heat

with this wall, and the system becomes a compact,







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

"-









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

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

stress-and-strain problems in the inner pipe and the ease

of replacement of tubes in case of failure, while not

suffering from the lack of versatility that usually

handicaps the use of the bayonet tube.

A review of the literature reveals that there is a

lack of basic research dealing with the heat exchanger.

Hurd [1] studied the mean temperature difference in the

bayonet tube for four different arrangements of fluid flow

directions. Large temperature differences were found when

flow direction in the annulus was counter to the flow

direction of the heat exchanger medium at the shell side.

A more recent study dealt with the use of the heat

exchanger in solar collectors [2]. In this application,

the outside pipe of the exchanger becomes the inner wall of

a dewar, and the outside surface of this pipe is coated

with a spectral selective coating to enhance its heat

absorption. Owens-Illinois has used this concept in

developing its SUNPAK collectors with success [3].

It is quite natural for one who is interested in fluid

flow and heat transfer analysis to consider the SPRF

exchanger as a combination of two separate components,

circular pipe and annulus, because their fluid flow and

heat transfer have been thoroughly studied in the

literature, e.g. [4]. However, because of the presence of









the round cap at one end of the heat exchanger, the fluid

flow in the upstream region may be affected by that in the

downstream region. As a remuit, to divide the exchanger

into two subsystems and to patch up their solutions at

points of matching may be an oversimplification of this

complicated problem. A realistic modelling of the system

will require the treatment of the system as a single

unit. This motivates the research of this dissertation.

In this work, the flow is assumed to be laminar,

axisymmetric, incompressible, and in steady state. Stream

function and vorticity are used in the solution of the

problem (Chapter II). The original irregular geometry of

the system in the physical plane is transformed to a

rectangular domain with square grids in the computational

plane. Grid lines are clustered to provide good resolution

in regions where flow and temperature fields change most

rapidly (Chapter III). The problem is solved numerically

with an iterative scheme and, in order to assure a

convergent solution, diagonal matrix elements in the

governing matrix equation are maximized (Chapter IV).

Pressure distribution and Nusselt number are evaluated at

Re=100, 500, and 1000, and Pr=0.7, 20, and 100 (Chapter

V). Thus, laminar flow is stable and the use of the steady

flow form of the governing equation is physically

credible. Numerical results indicate that the SPRF

exchanger is superior to heat transfer in a circular




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 one-hundredth of its radius, while the inner end of

the pipe that is facing the spherical cap is rounded. The

radius of the inner pipe was chosen such that the cross

section of flow in the pipe is equal to the ring passage

between the pipe and the cap, and this area is also equal

to the cross section of the annulus region; see dashed line

marks in the figure. The length of the pipes is

approximately twelve times the radius of the inner pipe, a

length short enough to show the entrance effect, yet long

enough so that the boundary conditions imposed at the exit

section of the annulus will have little effect on the flow

conditions inside the annulus.



II.1 Governing Equations

For the problem under investigation, the flow is

treated as steady, incompressible, and laminar. The fluid

medium is assumed to have constant thermophysical

properties. Both the flow and the temperature field are








































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=(Tw-T)/(Tw-Te), where Tw and Te refer to the

inner wall temperature and the fluid inlet temperature,

respectively. There are five unknowns (u,v,0,v,P) in these

equations, and there are five equations to solve them.



11.2 Boundary Conditions

Because of the large number of conditions to be used

in the analysis, it is convenient to summarize them

according to the locations where they will be applied.

Inlet: The velocity and temperature are both uniform at

the inlet. Here the v component of velocity is

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

that the fluid enters the heat exchanger through the inner

pipe, and the fluid actions are more intense near the wall

and around the turning point (h) at the tip of the inner

pipe; a fine grid is necessary to study the fluid flow and

heat transfer in these regions. To facilitate grid

generation to meet these requirements, the total system is

subdivided into three domains: domain A in bcdghib, domain

B in abija, and domain C in fgdef; and transformations will












Constant n line


Constant E line


Physical Plane


* *
Sb c d e

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 fine-tune those grid lines

in this domain so that a dense grid is formed close to the

walls and near the turning point of the inner pipe.

Finally, grids in domains B and C are filled in so that

they match those grids already laid out in domain A. The

details of these transformations will be spelled out in the

three sections that follow.



III.1 Grid Generation in Domain A

Differential equation methods will be used to generate

grids in domain A. Specifically, the elliptic partial

differential equation technique developed by Thompson,

Thames, and Mastin (TTM) will be used [8]. In order to

simplify computation of the velocity and temperature fields

to be made later, the grids in domain A in the

computational plane are taken to be squares, and Laplace

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

relations are usually referred to in the literature as

boundary conditions for the transformation equations (3.3)

and (3.4). In actuality, they provide the necessary

conditions so that constant n lines conform to the physical

boundary of the actual system.

To lay out z coordinates that enable a dense grid to

be formed near the turning point h, a cubic polynomial









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


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

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

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


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


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

and

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


Iteration is used in conjunction with successive line
over-relaxation to solve these equations. The relaxation
procedure is formulated as follows for each line (i.e.,

for a fixed j index):


k+1 k k
zi, = zi, + w (zi -zi ) (3.15)

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 < 10-4 and r. .-r., < 10-


The grid generated by the TTM solver is shown in

figure 4 (a). The spacing along the z axis is

satisfactory; however, the spacing in the r direction is

not, which will be corrected later. It should be repeated

that these grids are generated by using square grids in the

computational plane. Also notice that, for the

transformation made that is illustrated in this figure, 103

points are used along the z direction, 41 points in the r

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 4-3;
lj niI 1---4~ 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 grid-line clustering effort may be

spared. Sorenson and Stegernfound out that this method of

clustering grid lines may lead to complexities or even

unsatisfactory results in the computation [10]. For this

reason, the grids are not transformed using Poisson's

equations in this work; Laplace equations are used instead,

and the generated grid lines are clustered later by means

of clustering functions. In this effort, Sorenson and

Steger's method will be used as described below.

The grid lines in the z direction have been laid out

satisfactorily; no changes need be made of them. Lines in

the r direction must be clustered so that a dense grid is

where the change in the gradient of the flow variable is

most rapid. In other words, lines of constant n in the

physical plane must be repositioned for each line of

constant 5. This clustering of n lines can be accomplished

in two steps as follows.

(1) The arc length along each constant E line in

figure 4 (a) is computed with the following relation:



i,j+1 = Si,j+[(zi,j+-zi,j)2+(ri,j+-ri,)21/2 (3.17)


Since the first point is indexed 1, Si,1 = 0 for the sake

of consistency. Note that j refers to the index for the

jth line of constant n; the arc length is computed for a









fixed i (or constant c). Equation (3.17) can be used to

compute the arc length for each grid point in figure 4 (a),

and this yields a table of zi,j and ri,j as functions Si,j-

This part of the work is essentially used for locating

grid positions in the physical plane. Once they are found,

the constant n lines in the physical plane may be

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

this method, the nth iterated e is expressed as


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


where

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

and


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

X[n-1( -2)-1]} (3.22)


Notice that this ei must be found for each C line in the

reconstruction of new grids.

This clustering scheme described above monotonically

increases the grid spacing away from the outer pipe wall.

This creates a slight problem for those grids near the

inner pipe wall because the rate of change of the gradient

of the flow variable is also rapid there. It is necessary

to cluster the grid lines near boundary ihg (figure 3) for

good resolution. To accomplish this, an image method is

employed as follows.









In the image method, the S.. in Eq. (3.18) is
1, max
redefined to be S. m/2, which length corresponds to
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)/2-1] to 1.

The clustered grid lines are shown in figure 4 (b),

where ASO is 0.25% of the inner pipe radius. An enlarged

grid near the turning point is shown in figure 5. The

clustered grid lines appear to be uniformly distributed

near the pipe wall and the turning point of the pipe. The

grid generation task in domain A is thus completed.



III.3 Grid Generation in Domain B and C

The grid generation in domain B and C is rather

straight-forward once the grids are laid out in domain A.

Recall that the positions of the grid lines in the z

direction in the physical plane are generated by using Eq.

(3.6), in which values for z0, zf, io, if, AZ1, and Azf

must be provided as input; the grid spacing increases with

the index i. Then, in domain B, the z positions for the

grid lines must be laid out from a to b because close to

the inlet section a dense grid is desired. In domain C,

however, the direction of z must be reversed because, at

the exit section, the flow should approach fully developed;

a dense grid there is unnecessary.























































































































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 non-uniform grid net in domains B and C. A plot

of the finished grid is shown in figure 6. Because domains

B and C are now linked to domain A, there are a total of

163 points in the z direction; the number of points across

the flow passage in the r direction is still 41, which is

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-
=-[(Ur-p ) + -(rn-r )] (4.1)
r E [( r p~E








Vorticity definition equation:

Jz Jz
aS 28n + +(J2+ ) +(j2Q_ )
En Tn r (r n

= J2r (4.2)

Velocity equations:

1 (z Z P) (4.3)
u = (zJ n znr) (4

1 (r r ) (4.4)
v = -- (r n rn) (4.4)


Energy equation:

Jz Jz
00 2B0 + ye +(j2p_ )0 2 +(J Q+ -)o
En n) r ( r n

JPe (45)
JPe (ne 90) (4.5)
r nC n1

Boundary conditions can be derived accordingly

(Appendix A). Notice that, in the application of boundary

conditions, all conditions of the first kind will be used

as is; no transformations of them need be made. Other

conditions must be changed to the computational plane, and

they are listed together with those first-kind conditions

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 =

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


+(C10+ReC16T) i,j+ +(C11-ReC16 ) i,j-i+C0Tn (4.12)


Vorticity definition equation:

Cli,j = C2 i1+C3 i-,+C4ij+1+Ci,-1

+C6n. .+COQ0 (4.13)


Velocity equations:

u ij = C12T1 C13IE (4.14)
1,] T1


vij = C14Ti C15 (


(4.15)









Energy equation:

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


+(C10+PeC16 ) i,j+1+(C11i-PeC169 )ij, I


+CO 6 (4.16)


In these equations, CO through C18 represent metric

coefficients listed in Appendix C. Those 9, n, and 0

topped with bars are compact expressions for their

increments in i and j directions; for example,



S i+1,j i-,j


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

and
b = i + +J
Si+1,j+1 i+1,j-1 i-1,j-1 i-1,j+1

It is strictly a matter of computer programming expediency

that these substitutions are made.

The discretized equations given above can be solved by

an iterative method, in which the problems encountered

include the nonlinearity of the vorticity transport

equation and the high convection rate in the energy

equation. Although use of an under-relaxation scheme in

the iterative method may ease these problems, a small

relaxation factor means an excessive computer time. For









these reasons, the matrix coefficients in the matrix

equations are reconstructed to make those diagonal elements

dominant in the coefficient matrix so that an iteration

will yield a convergent solution.

The method used for maximizing the matrix elements is

somewhat similar to the scheme developed by Takemitsu [11]

and is rooted in the up-wind scheme that is commonly used

in finite-differencing flow terms in the solution of fluid

dynamics problems [12]. The vorticity transport equation

is expressed as


C19. = C20n .+C21. +C22.
1,j 1+1,j i-1,j 1,j+1

+C23.i,j1 +CO"T (4.17)


This equation is a compact version of the vorticity

equation derived in Appendix C, Eq. (C.35).

The energy equation is derived as


C240ij = C250 +C26e0 +C270


+C280ij_-1.COOn (4.18)


In both equations, C19 through C28 represent



C19 = C7+Re[C17*T C18in +2( C16i5j + JC16i j)]


C20 = C8-Re(C16!T IC16 1j )









C21 = C9+Re(C16 + IC16I )

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

C23 = C11-Re(C16* C16 )

C24 = C1+2Pe( C16Tj + C16Ji )

C25 = C8-Pe(C166n C16n )

C26 = C9+Pe(C16Tn + C165 )

C27 = C10+Pe(C16P + IC16~ )

and

C28 = C11-Pe(C16, C164l )


Successive-over-relaxation (SOR) methods will be used
to solve these equations. In this method, the SOR versions
of (4.13), (4.17), and (4.18) are, respectively,


k+l k w( k k+l
+ 1 = (1-w ) + --- (C2*. + C3i_
iCj *i ), i +1,j i-l,j
k k+1 cn+ -*
+C4P + C5k1- + c6n+ + CO ) (4.19)
i,+1 ij-1 1, En


k+l k w k k+1
S+ = (1- )k + + (C20k2 + C21
S1, -w i, C19 i+1,j i-l,j
+C22k k+l -*
+C22k + C23k+1i,- + CCO ) (4.20)
1 +1 j- 1 EnT










k+1 k we k k+l
S = (1- e)0 + (C250 + C26e1

k k+1 -*
+C27 + C288k+1 + COe ) (4.21)
2i,j+1 + 10- I n


where w., wn, and we denote the relaxation factors for 4,

n, and e, respectively. Superscript k represents the inner

iteration count in the numerical solution; the superscript

n for ni,j in (4.19) stands for the outer

iteration count discussed later. The cross-derivative
_* _x _*
terms n', Q and 09 are evaluated using



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


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

and

-* k k+l k+1 kk
n i+l,j+l i+1,j-1 i-1,j-1 ~ -l,j+1


Iterations sweeps run from i=1 to imax and j=1 to jmax"

Discretization of the governing equations is now

complete; attention will now be directed to the boundary

conditions which must be processed accordingly. In the

discretization of boundary conditions, it must be noted

that, in the numerical solution to be given later, 9 and n

must be solved prior to the solution of the velocity and

temperature fields. It is thus appropriate to discretize

the and n conditions first. In doing so, conditions of









the first kind will again be used as is. Therefore, only

conditions of the second kind will be addressed. These

second-kind conditions will be discretized in the order of

their appearance in conditions (4.6) to (4.11) below.

The exit conditions Eq. (4.7a) are discretized by

using a second-order central-difference scheme. At the

exit section, for example, *, = 0 is discretized to be


= maxi mj 4 (4.22)
i +1,j i -1i,j
max max-


This identity will be used to substitute the values of at

index (imax+1,j) when the difference equation (4.13) is

evaluated at the exit section. This method also applies to

n, u, and v.

The wall vorticity conditions (4.11) are approximated

by the relation



"i,w = Y [2(i,w+1 i,w)] (4.25)


Here the bracketed terms are used to approximate -'p and

subscript (w+1) for refers to the grid location being one

point away from the wall. This equation has only

first-order accuracy but is adequate for use in the present

work because of its conservation property as reported by

Parmentier and Torrance [13].








Conditions for e and u are treated next, and it is

noted that, at the exit, the grids are orthogonal and

P = 0. Equation (4.5) reduces to


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

n-derivative quantities, and 0 is discretized by the

second-order backward difference, giving


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


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

by



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


in which the original u and u have been approximated by


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









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

again a first-order appro:.irmtion.

Temperature condition along the same axis is treated

in a similar manner,


1 + (z /z (4.28)
ei,1 (z /z )+1 i,2 z )i 1] (428)
n S

At the outer pipe wall,



ei,1 (B/y)+1 [ei,2 + (B/Y)Oi_-1,1 (4.29)



IV.3 Numerical Solution Procedure

A close examination of Eqs. (4.13) to (4.15), (4.17),

(4.18), and (4.25) reveals that

(1) the p equation, (4.13), contains all p but one a

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

used as the initial guess for high Reynolds-number cases.

Thus, for example, the output of t and a for Re=100 is used

as input for Re=500 case in step (1). For Re=100, the

initial guess of the vorticity values is taken to be zero

throughout. As for the stream function, a uniform flow is

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 < 10-4
i,j 1- i,j id


(ii) The convergence for inner iteration of Eqs. (4.19)

and (4.20) is set at


k+1 k < 10-3 and [k+1 k < 10-
1,0 i,j i,j i,J


The inner iteration steps are limited to 50 in order to

prevent it from converging to an erroneous solution [5].

It was found that once the outer iteration has converged as

desired, the inner convergence criterion can be met with

just one iteration. A flow chart for the numerical

procedure is given in figure 7.

Once the stream function is found, the velocity field

is calculated from Eqs. (4.14) and (4.15) and the energy

equation can be solved in a straight-forward manner. In

the solution of the energy equation, the second-kind

boundary conditions are treated explicitly. Miyakoda's

method [16] is used in this effort, in which the derivative

boundary conditions are incorporated into nodal equations

for interior nodes that are located one point away from the

insulated boundary (that is, the center line or the outer

wall). Hence, one point away from the center line,





















































Figure 7 Block Diagram Showing the Numerical Procedure
Used to Solve Stream Function and Vorticity










[C23 (z ) 1i,2 C24i+l2+C25i-1 2+C26 i3


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


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


[C23 ] = C240. +C250 +C260
(-/y)+1 .i,2 +1,2 i-1,2 i,3

+C27 ( 1 i,1+COen (4.34)


It is also noted that, before applying these equations
to find ei,2, 0i-_,1 must be updated by using (4.28) and

(4.29). The relaxation parameter we is taken to be 0.7 for
all cases studied. A fast convergence results if the

iteration proceeds from i=1 to imax' that is following the
flow direction, and J=Jmax to 1, that is, from the inner

wall to the outer wall.












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 pipe-flow problems, is not suited for use

in the present investigation because it terminates at the

center of the end cap. The inner pipe wall, fortunately,

has a continuous contour all the way from the inlet to

exit; it is thus used as a basis for evaluating the

pressure drop.









The Navier-Stokes equations can be used together with

the continuity equation and vorticity to derive the dimen-

sionless pressure gradients at the wall as (Appendix D)



1 1 arn
az = --" (5.1)
az Re r ar

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









In this equation,



C = r2[(1+r2) 1-2 ] (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) n-z 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 r-4
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 no-slip condition simultaneously. Because of the

closeness of this section to the inlet, the momentum

transfer has yet to reach that far to the center region,

and as a result, only velocity nearby is affected and a

maximum velocity develops locally. Such a phenomenon has

also been discussed in references [20] and [21].









A vortex is found in the annulus right after the flow

makes the 180-degree turn. As shown in figures 10 through

12, the size of the recirculation zone increases with

Reynolds number. For a quantification of this observation,

the sizes of the vortex measured in distance from the

turning point to the point of flow reattachment are listed

in Table 1. With the recirculation region, fluid moves in

a clockwise loop. The main flow detaches from the wall,

resulting in a depression of the heat transfer rate as will

be discussed later.

It is also interesting to note that, because of the

blockage of flow by the round cap, the maximum velocity in

the inner pipe is not located on its centerline, but at a

location slightly shifted away from the center. By the

same token, the maximum velocity in the annulus is not

located close to the inner wall, as is normally the case

for an annular flow [17], but is displaced slightly toward

the outer wall. This latter phenomenon is certainly a

result of the accelerated main flow that occurs near the

recirculation region being displaced away from the inner

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
-I-0

-




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 180-degree turn,

and correspondingly, the pressure drops as will be

discussed later.

As shown in figure 13, the vorticity changes sign in

the recirculation region, also see Table 2. Finally, as

the fluid moves further downstream toward exit, three

vorticity curves all merge into one.











Inner-Wall


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




4-1






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 boundary-Layer equation was used and

the momentum equation was linearized locally at any cross

section z1 by means of the velocity at z=z1+Az. However,

in the present investigation, no linearization is made and

the axial diffusion terms are retained in the momentum

equations. As a result, this positive pressure gradient

near the inlet region may be attributed to the inclusion of

the axial diffusion terms in these equations. The

requirement for satisfying the uniform flow boundary

condition at inlet also contributes to this pressure

rise. Notice that this pressure abnormality was also

reported by Morihara and Cheng [20] for flow in a straight

channel, and they used a similar mathematical model as used

in the present study. A point of interest is that both the

pressure rise and the distance over which it occurs are

larger for larger Reynolds numbers, which is not

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-




r-t


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

Nusselt-number plots presented in the next section.



VI.5 Local Nusselt Number

Local Nusselt numbers for the nine cases are shown in

figures 18 through 20 at three Prandtl numbers. Two pieces

of information are also included in figure 18; one relates

to the Nusselt number for fully developed flow in a

circular pipe under a constant-wall-temperature

condition. The Nusselt number for this case is 3.658 [18];

see horizontal line drawn. The other set of data is taken

from the numerical work of Hornbeck [23] who studied the

Nusselt number under a developing flow and heat transfer

condition at the same Prandtl number (0.7). His results

are plotted using stars in the figure.

Referring to figure 18, it is seen that the present

Nuz values are in good agreement with the results of

Hornbeck for Re=100. Moving downstream, the local Nusselt

numbers asymptotically approach those for fully developed

flow. Near the turning point, the turn takes effect,

resulting in an increase in heat transfer and the local

Nusselt numbers increase rapidly. Again, to assist this

observation, a dashed line is drawn in this figure to
















c'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 g-1



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


S-Pr=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 (n-2) degree of freedom (n

represents number of correlated point) is 2.6 [25].

The present heat exchanger is found to be far more

effective in heat transfer. To verify this point,

Hornbeck's data [23] for air (Pr=0.7) under a developing

velocity and temperature condition in a circular pipe are

used for comparison. As listed in Table 5, the mean

Nusselt number for the SPRF heat exchanger is 1.57 times

that for pipe flow when the Reynolds number is 500.

Doubling this Reynolds number enables this Nusselt-number

ratio to be raised to 1.87. The superiority of the present

exchanger is thus clearly seen.



VI.7 Truncation Error and Convergence

In order to check the truncation error and the

convergence of the solutions presented, three different

grid systems are used. It is noted that, for the ease of

comparison, the number of grid points along the constant n-

line is fixed at 163; the number of grid points along the

constant S-line is tested for 21, 41, and 61. This results

in lines of constant E for these grid systems having the















Table 5

Comparison of Mean Nusselt Number Between
SPRF Heat Exchanger and Heating in a Circular Pipe




Mean Nusselt Number
Mean Nusselt-
Pr Re Hornbeck's Data* Present Study Number Ratio


0.7 500 7.4 11.62 1.57



0.7 1000 9.57 17.88 1.87



* Developing velocity and temperature profiles in a
circular pipe









same contours in the physical plane but with different grid

sizes.

Figure 22 shows the competed stream functions along

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

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

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

corresponds to the constant C-line at the turning point

(see figure 5). It is clear that the computed stream

function values based on the use of grids 163x41 converge

to the results of the dense grid (163x61). A close

examination of the difference between these two sets of

data shows that the maximum deviation is only about 3%

relative to the dense grid. However, the CPU time required

for computation with the dense grid is increased by a

factor of about 3.5 (Table 6) using the Harris 800 computer

on campus. For this reason, a denser grid is not used for

the present computation. The results for 163x41 are

considered converged.









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 C-line for Different Grid Systems


(n-Tmin)/ (nmax-Tmin)













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

exchanger. Stream function and vorticity are used in the

formulation, and the original irregular geometry of the

exchanger in the physical plane is transformed into a

rectangular domain with square grids in a computational

plane. Elliptic partial differential equations are used in

the transformation, and grid lines are clustered with the

use of separate clustering functions which allow fine grids

to be formed in regions where the flow and temperature

fields change most rapidly. The source terms in the

Poisson's equations are also evaluated in the solution.

A finite difference method is used to solve the

problem in the computational plane. In the numerical

algorithm, in order to assure a convergent solution in the

iteration, the diagonal terms in the resulting matrix

equations are reconstructed to make them dominant.

Vorticity boundary conditions at the walls are expressed in

a first-order approximation. The velocity and temperature

conditions along the centerline of the heat exchanger are

also finite difference, with the first-order differential




Full Text

PAGE 1

FLUID FLOW AND HEAT TRANSFER IN A SINGLE-PASS, RETURN-FLOW HEAT EXCHANGER BY SONG-LIN 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, Li-Pin. Her patience and encouragement made his graduate study a beautiful memory. Finally, the author wishes to acknowledge the Interactive Graphics Laboratory (IGL) and the Northeast Regional Data Center (NERDC) of the State University System of Florida for allowing him to utilize the computing facilities on campus. 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 Inner-Wall 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 Single-Pass, Return-Flow 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 Re-100, 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) C1-C18 metric coefficients defined in Appendix C C19-C28 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 SINGLE-PASS, RETURN-FLOW HEAT EXCHANGER BY SONG-LIN YANG May, 1985 Chairman: C. K. Hsieh Major Department: Mechanical Engineering Numerical solutions for flow and heat transfer in a single-pass, return-flow heat exchanger are given for Reynolds numbers (Re) of 100, 500, and 1000, and Prandtl numbers (Pr) of 0.7, 20, and 100. The flow is modeled as laminar, axisymraetric, incompressible, and in steady state, and the Navier-Stok.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 finite-difference forms, and in order to assure convergence in the numerical iteration, the matrix coefficients in the matrix equations are reconstructed to make those diagonal terms dominant in the coefficient 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 single-pass, return-flow (SPRF) heat exchanger shown in figure 1 consists of a circular pipe inserted inside an outer pipe. The outer pipe is sealed at one end so that fluid entering the inner pipe may be deflected into the annulus formed between the inner and the outer pipe of the exchanger. Either the inner or the outer pipe wall, or both, may be heated, and the heat input may be simulated by different conditions at these walls. The heat exchanger described above is structurally similar to the bayonet tube studies by Hurd [1] in 1946; however, they serve different purposes. In the bayonet tube, the double-pipe construction is essentially used for transporting fluid. Hence, the entering and exiting fluids exchange heat with each other via the heat conducting inner wall, and whether heat is transferred in or out of the system depends on the temperature difference between the fluid in the annulus and the medium surrounding it at the shell side. This is not so, however, for the heat exchanger, in which the inner pipe wall may be heated. Then, the entering and exiting fluids both exchange 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

stand-alone device. In fact, the exchanger possesses all tne virtues of the bayonet tube, such as the freedom from stress-and-strain problems in the inner pipe and the ease of replacement of tubes in case of failure, while not suffering from the lack of versatility that usually handicaps the use of the bayonet tube. A review of the literature reveals that there is a lack of basic research dealing with the heat exchanger. Hurd [1] studied the mean temperature difference in the bayonet tube for four different arrangements of fluid flow directions. Large temperature differences were found when flow direction in the annulus was counter to the flow direction of the heat exchanger medium at the shell side. A more recent study dealt with the use of the heat exchanger in solar collectors [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. Owens-Illinois has used this concept in developing its SUNPAK collectors with success [3]. It is quite natural for one who is interested in fluid flow and heat transfer analysis to consider the SPRF exchanger as a combination of two separate components, circular pipe and annulus, because their fluid flow and heat transfer have been thoroughly studied in the literature, e.g. [4]. However, because of the presence of

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 one-hundredth of its radius, while the inner end of the pipe that is facing the spherical cap is rounded. The radius of the inner pipe was chosen such that the cross section of flow in the pipe is equal to the ring passage between the pipe and the cap, and this area is also equal to the cross section of the annulus region; see dashed line marks in the figure. The length of the pipes is approximately twelve times the radius of the inner pipe, a length short enough to show the entrance effect, yet long enough so that the boundary conditions imposed at the exit section of the annulus will have little effect on the flow conditions inside the annulus. II . 1 Governing Equations For the problem under investigation, the flow is treated as steady, incompressible, and laminar. The fluid medium is assumed to have constant thermophysical properties. Both the flow and the temperature field are

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, t|j , 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 one-to-one. Notice that the fluid enters the heat exchanger through the inner pipe, and the fluid actions are more intense near the wall and around the turning point (h) at the tip of the inner pipe; a fine grid is necessary to study the fluid flow and heat transfer in these regions. To facilitate grid generation to meet these requirements, the total system is subdivided into three domains: domain A in bcdghib, domain B in abija, and domain C in fgdef; and transformations will 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 fine-tune those grid lines in this domain so that a dense grid is formed close to the walls and near the turning point of the inner pipe. Finally, grids in domains B and C are filled in so that they match those grids already laid out in domain A. The details of these transformations will be spelled out in the three sections that follow. 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 -K-X-)(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 one-to-one mapping relations are usually referred to in the literature as boundary conditions for the transformation equations (3.3) and (3.4). In actuality, they provide the necessary conditions so that constant n lines conform to the physical boundary of the actual system. To lay out z coordinates that enable a dense grid to be formed near the turning point h, a cubic polynomial

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,j-az,^^^. = "Y ( ^i , j_, ^^i , j,i ) ^ f ^.^ (3-8) ari_,,j-2(a^y)r,^..ar,,,^. = "T ( ^i , -j_, ^i , j.^ ) i ?^^ (3-9) where « =Tt^^i, 0^1-^1.3-1^'^ (^i,d.i-i,d-i)'^ (5.10) e =T t(^i.i,j-"i-i,j^("i,j^i-"i,j-i^ -h^-u^,^^i-^,j^ ^ ^'u^,3-'i-^,3^ ^ (3.12) ^n ^ ("i^1,j^1-^i-1,j-l'"i-1,j-l"^i-1,j-1^ (^•''^^ and ^5n = (^i.1.j.1-^i.1,j-1^^i-1,d-1-^-1,j.l) (^-'^^ Iteration is used in conjunction with successive line over-relaxation 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 i-H 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 grid-line 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 Newton-Raphson method may be used to find this e. In this method, the nth iterated e is expressed as F(e"~^) ,n ^ ^n-1 3^_^ (3.20) ^ ^ F'(sr') 1 where FC."--') = S, , . ^ UU^I-h^"'^'' -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 straight-forward once the grids are laid out in domain A. Recall that the positions of the grid lines in the z direction in the physical plane are generated by using Eq. (3.6), in which values for 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 non-uniform grid net in domains B and C. A plot of the finished grid is shown in figure 6. Because domains B and C are now linked to domain A, there are a total of 163 points in the z direction; the number of points across the flow passage in the r direction is still 41, which is 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 first-kind 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^)]^. . = (C8-ReCl6^^)^^^^^j + (C9+ReCl6l'^)^._^^j + {C^0+ReC^6i>.)^. . .+(C11-ReCl6ij;.)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. . = (C8-PeCl6?^)e +(C9+PeCi6;F_)e + (C10H-PeCl6;F.)e. . . + (C11-PeCl6^.)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 n-i,j and *n "^1,0 + 1 "^1,0-1 '''^n '''i + 1,j + 1 " *i + 1,j-1 " *i-1,j-1 '^i-1,j + 1 It is strictly a matter of computer programming expediency that these substitutions are made. The discretized equations given above can be solved by an iterative method, in which the problems encountered include the nonlinearity of the vorticity transport equation and the high convection rate in the energy equation. Although use of an under-ralaxation 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 up-wind scheme that is commonly used in finite-differencing 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 1-1 ,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 1-1 ,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 = C8-Re(Cl6iF^ Cl6ii^^| )

PAGE 49

C21 = C9+Re(Cl6^^ + Cl64i J ) C22 = C10+Re(Cl6^^ + Cl64iJ ) 34 C23 = C11-Re(Cl6^^ C16^, and C24 = C1+2Pe(Cl6^^ + |C16^^ ) C25 = C8-Pe(Cl6?^ |C16^^ ) C26 = C9+Pe(Cl6<'^ + C16^^| ) C27 = C10+Pe(Cl64'. + C16^^ ) C28 = C11-Pe(Cl6ii)^ C164;^| ) Successive-over-relaxation (SOR) methods will be used to solve these equations. In this method, the SOR versions of (4.13), (4.17), and (4.18) are, respectively, '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 cross-derivative _* _* _* terms ^ c > ^r t ai^d ©^ are evaluated using T* _ ,k k+1 k+1 _ k '''en '^i + l.j + l " '^i+Lj-l " "^i-l.j-l *i-1,j + 1 ^^ = ^ • .^ • ^ ^ • .^ • -, + ^ • -1 y, ^ y, y, Cn 1+1, j + 1 1 + 1, j_i i_i,j_i 1-1, j + 1 and "^n ^i + 1,j + 1 ' ^i + 1,j-1 " ^i-1,j-1 " ^i-1,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 second-kind conditions will be discretized in the order of their appearance in conditions (4.6) to (4.11) below. The exit conditions Eq. (4.7a) are discretized by using a second-order central-difference scheme. At the exit section, for example, 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 first-order 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 n-derivative quantities, and 0^. is discretized by the second-order 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 first-order approximation. Temperature condition along the same axis is treated in a similar manner, Qi.1 = (z /z ).1 tei,2 ^ (%/^)Qi-1,lJ (4.28) At the outer pipe wall, ^i,1 =Ti7TnT ^®i,2 ^ (e/^)9i-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.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 Reynolds-number cases are used as the initial guess for high Reynolds-number 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 straight-forward manner. In the solution of the energy equation, the second-kind boundary conditions are treated explicitly. Miyakoda's method [16] is used in this effort, in which the derivative boundary conditions are incorporated into nodal equations for interior nodes that are located one point away from the insulated boundary (that is, the center line or the outer wall). Hence, one point away from the center line,

PAGE 57

42 Set initial guess fpr ipl j and ^. j Update wall vorticity Q^^ t.; Solve for vorticity fi. .• by SOR u-n yes(^.^j) Damped by .i;^U^J^j+w^(^i,j-"?,j) jk Solve for stream function i|j^. ^^ by SOR -noyes(i|..^j) Damped by ^^^5=^;^j+W2(^,.^j-4^';,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,1-CQ^. (4.33) and one point awiay. from theouterwall t^23--X67^]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> ®i-1 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 pipe-flow problems, is not suited for use in the present investigation because it terminates at the center of the end cap. The inner pipe wall, fortunately, has a continuous contour all the way from the inlet to exit; it is thus used as a basis for evaluating the pressure drop. 44

PAGE 60

45 The Navier-Stokes 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 n-S 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 (5-14) 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) f-t 0) •H -P •H O O r— I >

PAGE 70

55 / HI ! f ;ii o •p c •H o cu W) c •rH c t-l E-i 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 no-slip condition simultaneously. Because of the closeness of this section to the inlet, the momentum transfer has yet to reach that far to the center region, and as a result, only velocity nearby is affected and a maximum velocity develops locally. Such a phenomenon has also been discussed in references [20] and [21 J.

PAGE 74

59 A vortex is found in the annulus right after the flow makes the 180-degree turn. As shown in figures 10 through 12, the size of the recirculation zone increases with Reynolds number. For a quantification of this observation, the sizes of the vortex measured in distance from the turning point to the point of flow reattachment are listed in Table 1. With the recirculation region, fluid moves in a clockwise loop. The main flow detaches from the wall, resulting in a depression of the heat transfer rate as will be discussed later. It is also interesting to note that, because of the blockage of flow by the round cap, the maximum velocity in the inner pipe is not located on its centerline, but at a location slightly shifted away from the center. By the same token, the maximum velocity in the annulus is not located close to the inner wall, as is normally the case for an annular flow [17], but is displaced slightly 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 180-degree turn, and correspondingly, the pressure drops as will be discussed later. As shown in figure 13> the vorticity changes sign in the recirculation region, also see Table 2. Finally, as the fluid moves further downstream toward exit, three vorticity curves all merge into one.

PAGE 78

63 Table 2 Inner-Wall 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 inner-pipe 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 boundary-layer 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 Nusselt-number plots presented in the next section. VI. 5 Local Nusselt Number Local Nusselt numbers for the nine cases are shown in figures 18 through 20 at three Prandtl numbers. Two pieces of information are also included in figure 13; one relates to the Nusselt number for fully developed flow in a circular pipe under a constant-wall-temperature condition. The Nusselt number for this case is 3-658 [18]; see horizontal line drawn. The other set of data is taken from the numerical work of Hornbeck [23] who studied the Nusselt number under a developing flow and heat transfer condition at the same Prandtl number (0.7). His results are plotted using stars in the figure. Referring to figure 18, it is seen that the present 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°-528pr0-2S8 (6.1) m which is valid for Re=100 to 1000, and Pr=0.7 to 100. Variance evaluated based on (n-2) degree of freedom (n represents number of correlated point) is 2.6 [25]. The present heat exchanger is found to be far more effective in heat transfer. To verify this point, Hornbeck's data [23] for air (Pr=0.7) under a developing velocity and temperature condition in a circular pipe are used for comparison. As listed in Table 5, the mean Nusselt number for the SPRF heat exchanger is 1.57 times that for pipe flow when the Reynolds number is 500. Doubling this Reynolds number enables this Nusselt-number ratio to be raised to 1.87. The superiority of the present exchanger is thus clearly seen. VI. 7 Truncation Error and Convergence In order to check the truncation error and the convergence of the solutions presented, three different grid systems are used. It is noted that, for the ease of comparison, the number of grid points along the constant 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 5-lines for the three grid systems. In these figures, i=50 and 152 correspond to the mid-section in the inner pipe and the annular region, respectively; i=101 corresponds to the constant ^-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 C-line 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 single-pass, return-flow heat exchanger. Stream function and vorticity are used in the formulation, and the original irregular geometry of the exchanger in the physical plane is transformed into a rectangular domain with square grids in a computational plane. Elliptic partial differential equations are used in the transformation, and grid lines are clustered with the use of separate clustering functions which allow fine grids to be formed in regions where the 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 first-order approximation. The velocity and temperature conditions along the centerline of the heat exchanger are also finite differenced, with the first-order 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 180-degree 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 single-pass, return-flow 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 finite-difference 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 ^^^7-2" ^^zz " ^rr^ T? 3ri + (n + n ) 4(A. 10) zz rr 3n in which the coefficients for partial differentials on the right-hand 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 + z-z )/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^ (A-32) 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 right-hand 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. (4-3) 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, z-j to i-^, Zj_ to ij^,Z£_^ to if-1 , 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 second-order central difference scheme, as follows: (^c'l.a "<^.1, 0-^-1. a>/2 '^•^' (^«)i,j = ^1*1,0 -2^,0 *^i-i,d ( = -3) (F ) = (F F + F ^ CTi^i,j ^ i + 1,j + 1 1 + 1, j-1 ^ 1-1, j-1 ^-1, + 1^/4 (C.4) and 108

PAGE 124

(^.>l,a = "i,J*i ""i,J * 'i,3-i 109 (C.5) The first-order 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 1-1,0 F = F. . , F. . . (C.11) and ^^Ti ^ ^1 + 1,0 + 1 " ^i + 1,0-1 " ^i-1,d-1 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 . . 1-1, J i,J 1+1, J -•"^1,0*1 * ^i,a-i* *2^5n 'c-iJ) and cir. . . 2(ci + y)r. . + "r. . . 1-1,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 right-hand 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 central-difference versions of the original governing partial -differential-equations. 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 up-wind 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 right-hand 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,j-V-' (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 1-1, J i,J i,J-1 (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,0-1

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,d-1 (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 Navier-Stokes 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 = =lf-17 (D-7) 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 Chemis-cry, Volume 38, Number 12, 1946, pp. 1266-1271. 2. Beekly, D. C. , and Mather, G. R. , "Analysis and Experimental Test of a High-Perforraance , Evacuated Tubular Collector," NASA, CR-150874, 1978. 3. Moan, K. L., An Analysis of the Low Loss Evacuated Tubular Collector Using Air as the Heat Transfer Fluid, Owens-Illinois, 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. 247-256. 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. 325-340. 8. Thompson, J. F., Thames, F. C, and Mastin, C. W., "Boundary-Fitted Curvilinear Coordinate Systems for Solution of Partial Differential Equation on Field Containing Any Number of Arbitrary Two-Dimensional Bodies," NASA, CR-2729, 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. 82-3, 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 Finite-Difference Approximation for the Steady-State Navier-Stokes 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. 404-417. 14. Rigal, A., "Acceleration of the Convergence in Viscous Flow Computations," Journal of Computational Physics, Volume 43, 1981, pp. 177-179. 15. Roache, P. J., "Semidirect Calculation of Steady Twoand Three-Dimensional Flows," in Numerical Methods in Laminar and Turbulent of the First International Conference, John Wiley & Sons, Inc., New York, 1978, pp. 17-27. 16. Miyakoda, K., "Contribution to the Numerical Weather Prediction--Computation with Finite Difference," Japanese Journal of Geophysics, Volume 3, 1962, pp. 75-190. 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, McGraw-Hill, New York, 1980. 19. White, F. M., Viscous Fluid Flow, McGraw-Hill, 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. 550-572.

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. 837-844. • 22. Hornbeck, R.W., "Laminar Flow in the Entrance Region of a Pipe," Applied Scientific Research, Section A, Volume 13, 1964, pp. 224-232. 23. Hornbeck, R. W., "An All-Numerical Method for Heat Transfer in the Inlet of a Tube," American Society of Mechanical Engineering, Paper 65-WA/HT-36, 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 D-1972, 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 Song-Lin 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 Li-Pin Hsiao and now they have a lovely son named Albert Chung-Pu 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 Prolessor-of 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