Time-dependent enhanced heat transfer in oscillating pipe flow

MISSING IMAGE

Material Information

Title:
Time-dependent enhanced heat transfer in oscillating pipe flow
Physical Description:
xiii, 162 leaves : ill. ; 28 cm.
Language:
English
Creator:
Zhang, Guo-Jie, 1941-
Publication Date:

Subjects

Subjects / Keywords:
Heat -- Transmission   ( lcsh )
Pipe -- Fluid dynamics   ( lcsh )
Viscous flow   ( lcsh )
Laminar flow   ( lcsh )
Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1988.
Bibliography:
Includes bibliographical references.
Statement of Responsibility:
by Guo-Jie Zhang.
General Note:
Typescript.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 001113215
notis - AFK9835
oclc - 19897631
sobekcm - AA00004834_00001
System ID:
AA00004834:00001

Full Text











TIME-DEPENDENT ENHANCED HEAT TRANSFER
IN OSCILLATING PIPE FLOW











By

GUO-JIE ZHANG


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


1988





























To my beloved motherland














ACKNOWLEDGEMENTS


The author wishes to express his deep gratitude to

the chairman of his committee, Dr. Ulrich H. Kurzweg, for

the valuable assistance and advice in guiding this research

work to its completion and for thoroughly reviewing the

entire manuscript leading to the realization of this

dissertation. Also he wishes to express his appreciation to

Drs. E. Rune Lindgren, Lawrence E. Malvern, Arun K. Varma

and David W. Mikolaitis for the many helpful discussions in

the formulation of the problem and constructive suggestions

for overcoming many difficulties in the solution process.

Thanks are expressed here also to Dean Eugene R. Chenette,

to the department chairman, Dr. Martin A. Eisenberg, and to

Dr. Charles E. Taylor for their support in allowing the

author to pursue his educational goals within this lovely

country. Part of the work presented here was funded by a

grant from the National Science Foundation, under contract

number CBT-8611254. This support is gratefully

acknowledged.


iii

















TABLE OF CONTENTS


Page


ACKNOWLEDGEMENTS .

LIST OF FIGURES .

LIST OF TABLES .

KEY TO SYMBOLS .


ABSTRACT .

CHAPTERS

I INTRODUCTION


II FORMULATION OF THE PROBLEM .
Governing Equations .
Boundary Conditions .
Model 1 . .
Model 2 . .
Model 3 . .
Initial Conditions .
Calculation of Tidal Displacement
Effective Heat Flux .


III NUMERICAL TECHNIQUES EMPLOYED . .
Transformation . .
Crank-Nicolson Method for Momentum Equation .
ADI Method for Axisymmetric Heat Equations. .
Convergence Criteria . .
Grid Generation . .

IV NUMERICAL RESULTS AND DISCUSSION .
Part 1. Oscillatory Pipe Flow Features .
Velocity Profiles . .
Lagrangian Displacements . .
Tidal Displacements . .
Phase Lags . .
Part 2. Enhanced Heat Transfer Investigation
Periodic Temperature Build-up
in Thermal Pumping Process .


. . . iii

. . . vi

. . . ix

. . . x
ix

X


xii


18
18
. 24
24
28
30
. 33
35
. 37


r

r
r
r



rr









Temperature Distribution in Model 1 .
Temperature Distribution in Model 2 .
Temperature Distribution in Model 3 .
Heat Flux versus Tidal Displacement
(Model 2) . .
Influence of Thermodynamic Properties .
Heat Flux versus Tidal Displacement
(Model 1) . .
Influence of Wall Thickness .
Influence of Pipe Diameter .
Variation of Axial Temperature Gradient
In Model 3 . .
Comparison of Enhanced Oscillatory Heat
Transfer and Heat Conduction .
Enhanced Feat Flux as a Function of
Wormersley Number . .
Tuning Curves . .

V CONCLUDING REMARKS . .

APPENDIX: ETP COMPUTER CODE . .

REFERENCES . . .

BIOGRAPHICAL SKETCH . .


95
. 100
. 108

. 111
. 115

. 119
. 124
. 126

. 130

. 133

. 135
* 137

. 143

* 149

. 157

. 161















LIST OF FIGURES


Ficrure Page

1-1 F(a) Curve . . 7

2-1 Thermal Pumping Device . .... 19

2-2 Model 1---Fixed End Temperature Model. 25

2-3 Model 2---Periodic Heat and Cold Sources
on Insulated Wall . 29

2-4 Model 3---Pipe with Extended Conducting Sections 31

3-1 Grid System Used in the Numerical Simulations. 39

3-2 Coordinate Transformation ... 42

4-1 1-D Velocity Profiles in Oscillating Flow for
Wormersley Number a = 1, 10, 100, and 1000 65

4-2 Velocity Profiles (a = 5, after Uchida) 66

4-3 Magnified View of Velocity Profile near Wall 67

4-4 Lagrangian Displacement for a = 0.1, 1.0 70

4-5 Lagrangian Displacement at a = 10 71

4-6 Dimensionless Cross-Section Averaged
Displacement versus Time (a = 0.1 1.0) 73

4-7 Dimensionless Cross-Section Averaged
Displacement versus Time (a = 2 50). 73

4-8 Relationship Between Dimensionless Tidal
Displacement AX and Wormersley Number a 77

4-9 Relationship Between Tidal Displacement Ax and
Exciting Pressure Gradient in Water 78

4-10 Phase Variation Along Radius for Different
Wormersley Numbers . 84








4-11 Temperature Build-up Process in Oscillating
Flow (Model 1, a = 1, Ax = 2cm). .. 89

4-12 Temperature Build-up Process in Oscillating
Flow (Model 1, a = 1, Ax = 5cm) ... 90

4-13 Temperature Build-up Process in Oscillating
Flow (Model 1, a = 1, Ax = 10cm) 91

4-14 Temperature Build-up Process in Steady Flow
(Model 2, Ugve = 1.5 cm/sec) . 93

4-15 Build-up Time versus Tidal Displacement
(Model 2, a = 1). . 94

4-16 Temperature Distribution in Oscillating Pipe
Flow (Model 1, a = 1, Ax = 1cm) 96

4-17 Temperature Distribution in Oscillating Pipe
Flow (Model 1, a = 1, Ax = 2cm) 97

4-18 Temperature Distribution in Oscillating Pipe
Flow (Model 1, a = 1, Ax = 5cm) 98

4-19 Temperature Distribution in Oscillating Pipe
Flow (Model 2, a = 1, Ax = 1cm) 101

4-20 Temperature Distribution in Oscillating Pipe
Flow (Model 2, a = 1, Ax = 5cm). 102

4-21 Temperature Distribution in Oscillating Pipe
Flow (Model 2, a = 1, Ax = 10cm) 103

4-22 Temperature Distribution in Oscillating Pipe
Flow (Model 2, a = 1, Ax = 20cm) 104

4-23 Temperature Distribution in Oscillating Pipe
Flow (Model 1, a = 1, Ax = 30cm) .. 105

4-24 Temperature Distribution in Steady Flow
(Model 2, Uave = 0.5 7.5 cm/sec) 106

4-25 Temperature Distribution in Oscillating Pipe
Flow (Model 3, a = 1, Ax = 10cm) 109

4-26 Magnified View of Temperature in the Central
Pipe Section (Model 3, a = 1, Ax = 10cm) 110

4-27 Heat Flux in Oscillating Flow and Steady Flow
(Model 2, a = 1, Water as Working Fluid). 112


vii








4-28 Influence of Thermodynamic Properties of H20
on the Enhanced Heat Flux (Model 2, a = 1,
Ax = 10 cm). . . 117

4-29 Heat Flux versus Tidal Displacement
(Model 1, a = 1, Pr = 7.03). 120

4-30 Heat Flux versus Tidal Displacement
(Model 1, a = 3, Pr = 7.03). 121

4-31 Influence of Wall Thickness on Axial Heat Flux
(Model 1, Water-Glass, a = 1, Ax = 5cm). 126

4-32 Influence of Pipe Diameter on Heat Flux for
Fixed Frequency (Model 3, Water-glass, a = 3,
Ax = 10cm) . . 128

4-33 Typical Iso-Temperature Contour in Oscillating
Pipe Flow (Model 3, Water-Glass, a = 3,
Ax = 10cm). . 129

4-34 Variation of Temperature T1 and T2 versus Ax
(Model 3, Water-Glass, a = 3) ... 131

4-35 Comparison of Enhanced Heat Transfer and Heat
Conduction in Oscillating Pipe Flow
(Model 3, Water-Glass, a = 3) . 134

4-36 Variation of Axial Heat Flux versus Wormersley
Number (Model 3, H20-Glass, Hg-Steel,
Ax = 10cm) . . 136

4-37 Computed Tuning Curves (Model 3, H20-Steel and
Hg-Steel, Ax = 10cm) . 138

4-38 Tuning Curve versus Wormersley Number
(after Kurzweg). . 140

4-39 Ratio of Heat Flux due to Conduction to
Enhanced Heat Flux versus Wormersley Number
(Model 3, H20-Steel, Hg-Steel, Ax = 10cm) 141


viii














LIST OF TABLES


Tables Pages

4-1 Dimensionless Tidal Displacement at Different
Wormersley Numbers . .... 72

4-2 Phase Lags Along Radius (Working Medium: H20,
Ax = 10cm, a 0.1 20) . 80

4-3 Comparison of Phase Lags with Different
Working Mediums (Ax = 10cm, a 1, 5). 83

4-4 Comparison of Enhanced Heat Flux Using
Numerical Velocity with Heat Flux Using
Analytic Velocity (Model 3, H20-Glass,
R1 = 0.1cm, R2 = 0.15cm, Pr = 7.03, a = 3) 86

4-5 Enhanced Axial Heat Flux via Tidal Displacement. 113

4-6 Enhanced Axial Heat Flux in Steady Flow 114

4-7 The Influence of Properties of Water on the
Enhanced Axial Heat Flux . 116

4-8 Variation of the Axial Temperature Gradient
versus Wormersley Numbers (Water-Glass,
Ax = 10 cm) . . 132














KEY TO SYMBOLS


x x coordinate

r radial coordinate

t time

5 coordinate (x) in transformed plane

1 coordinate (r) in transformed plane

T transformed time

L pipe length

R1 pipe inner radius

R2 pipe outer radius

w oscillating frequency

6 boundary layer thickness

P pressure

p density

c specific heat

v Kinematic viscosity

p dynamic viscosity

K thermal diffusivity

k thermal conductivity

Ke coefficient of enhanced heat diffusivity

Pr Prandtl number

A = l/p ap/axI a measure of the maximum axial
pressure gradient (cm/sec2)








a Wormersley number a = Ju/ov

T Temperature

-I 7 = aT/ax time-averaged axial temperature
gradient

S = r/R1 dimensionless radial distance

g radial temperature distribution function

U velocity

Uo representative velocity

f velocity shape function

X Lagrangian displacement

DX dimensionless tidal displacement

Ax dimensional tidal displacement

Qtotal time averaged total enhanced axial heat flow over
pipe cross-section

axial heat flux




Subscript

f fluid

w wall

h hot

c cold

th thermal

eq equivalent

adj adjacent









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



TIME-DEPENDENT ENHANCED HEAT TRANSFER
IN OSCILLATING PIPE FLOW


By

Guo-jie Zhang

April 1988


Chairman: Ulrich H. Kurzweg
Major Department: Aerospace Engineering, Mechanics,
and Engineering Science

The problem of time-dependent enhanced heat transfer

in an incompressible viscous laminar fluid subjected to

sinusoidal oscillations in circular pipes which are

connected to a hot reservoir at one end and a cold reservoir

at the other end has been examined numerically in detail.

Three models were designed for the investigation of such an

enhanced thermal pumping process and a computer code (ETP)

was developed to implement all the numerical calculations.

To increase the understanding of the mechanism of

thermal pumping, the periodic velocity profiles and

Lagrangian displacements as well as tidal displacements at

various Wormersley numbers (from a = 0.1 to 1000) were

studied. Some transient problems of enhanced axial heat

transfer in oscillating pipe flow such as the periodic final

temperature build-up process in oscillating pipe flow were

also examined. The time-dependent temperature distribution

xii








in the different models was numerically studied in detail.

The enhanced axial heat flux magnitude versus different

tidal displacements with water and mercury as the working

fluids bounded by pipe walls of different material were

observed and the quadratic coefficients found. The

influence of the variation of water properties on the

enhanced axial heat flux was numerically examined and the

results show that the enhanced axial heat flux can vary

about 150 percent even within the temperature range from 0C

to 1000C. The effects of wall thickness and pipe diameter

in enhanced thermal pumping were also studied and the

optimum wall thickness was found to be about 20 percent of

the pipe radius in the water-glass combination. The tuning

effect in the water-steel and the mercury-steel cases was

examined and the results show good agreement with analytic

predictions. A comparison of the enhanced axial heat flux

with the axial heat flow due to heat conduction at various

tidal displacement and Wormersley numbers shows that the

latter is quite small and negligible provided the tuning

condition is satisfied.

This study has shown that the enhanced thermal pump

is indeed a very effective tool for those problems where

large amounts of heat must be transported without an

accompanying convective mass exchange. The investigation

also indicates that turbulent flow in the reservoirs is

preferable to laminar conditions and should receive more

attention in future studies.


xiii














CHAPTER I
INTRODUCTION


Enhanced heat transport in a viscous laminar fluid

subjected to sinusoidal oscillations in a very long pipe

which connects a hot fluid reservoir at one end and a cold

fluid reservoir at the other end (Fig. 2-1) has been

recognized and studied recently by Kurzweg [15, 16, 17, 22].

The results obtained show that with this oscillatory pipe

flow the heat transferred axially from the hot end to the

cold end can be orders of magnitude larger than that

obtained by pure molecular conduction in the absence of

oscillations. In addition, the more important thing of

interest is that this heat transfer process involves no net

convective mass transport. Major assumptions made in the

above cited studies on enhanced heat diffusion are that a

constant time-averaged non-zero axial temperature gradient

is always present in the oscillating flow and that the axial

molecular conduction along the wall and in the oscillating

fluid is negligible.

Discovery of this enhanced heat transport phenomenon

was made possible by earlier studies on axial dispersion of

contaminants within steady laminar flows through capillary

tubes by Taylor in 1953 [32], and Aris in 1956 [4]. These








2

earlier studies show that when a small quantity of a

contaminant is introduced into a circular pipe, the

dispersion of the resultant contaminant cloud is greatly

enhanced by the flow of the fluid. Bowden's 1965 results

show that similar dispersion effects occur in oscillatory

flow [8]. This enhanced axial dispersion of contaminants in

the presence of laminar oscillatory flow within capillary

tubes was studied in 1975 by Chatwin who suggested that the

assumption of constant time-averaged axial contaminant

gradient can be made [10]. Recently, Bohn et al. extended

this work to the study of gas component transfer in binary

gas mixtures when these are confined to single tubes and a

sinusoidal pressure variation is applied [7]. Further

studies in 1983 by Watson [38] show that the effective

diffusion of contaminants is proportional to the square of

the tidal displacement. This has been experimentally

verified by Joshi et al. [13] and by Jaeger [12], both in

1983, and most recently by Kurzweg and Jaeger [19], in 1987.

All these results show that the contaminant would spread

axially in both steady and oscillatory laminar pipe flow at

rates as much as five orders of magnitude higher than in the

absence of fluid motion.

The first significant research work extending these

enhanced axially diffusion studies to the heat transfer

problem in oscillatory flow within very slender pipes or

flat plate channels is due to Kurzweg [15, 16, 22]. In








3
early 1983 Kurzweg suggested that a similar dispersion

process should occur in the heat transfer area because of

the similarity in both the governing diffusion and heat

transfer equations [16], and the first preliminary theory

was formulated in 1985 [15, 16], in which, referring to

Chatwin's idea, a time averaged constant axial temperature

gradient assumption was used. The instantaneous temperature

distribution was taken to be of the form [15, 22].


T = y[x + R1 g(n)e ] (1-1)


where y = aT/ax is the time-averaged axial temperature

gradient, R1 is the tube radius, x is the axial distance

along the capillaries (with x = L/2 and x = -L/2 denoting

the ends), L is the pipe length under consideration, w is

the oscillating frequency of fluid, g(S) is a radial

temperature distribution, and = r/R1 is the dimensionless

radial distance. The theoretical analysis shows that under

certain conditions such enhanced axial heat diffusivity can

indeed be significantly larger than the axial molecular

conduction [15, 16] and this has been verified by some

experimental measurements by Kurzweg and Zhao [22].

In order to better understand the physical mechanisms

of this interesting and potentially useful heat transfer

technique, we shall examine in greater detail the thermal

pumping model shown in Fig. 2-1. It is assumed that a

bundle of very thin and long tubes is connected to a








4

reservoir which supplies unlimited hot liquid at one end to

a second reservoir which supplies unlimited cold liquid at

the other end. The liquid in the pipes is oscillating

axially with an amplitude such that none of the liquid which

is originally in the middle portion of pipes ever runs into

either reservoir. That is, there is no net convective mass

exchange between two reservoirs. The largest axial fluid

dimensionless displacement (when cross-section averaged) is

referred to as the nondimensional tidal displacement and is

denoted by "AX" (it should not be confused with the

dimensional tidal displacement "Ax" frequently used in the

present study). At time t = 0, the fluid within the pipes

is set into axial oscillations at angular frequency w and

tidal displacement AX. After a short transient, this

oscillatory motion will lead to very large axial heat flows

which can be readily made to exceed those possible with heat

pipes.

Before exploring the mechanisms of this enhanced heat

transport, it is necessary to introduce some new concepts

which are commonly used in the study of this type of

oscillatory motion.

As is well known [36], for high frequency viscous

laminar axial oscillations within fluid flow along rigid

pipes, the non-slip boundary condition creates a very thin

Stokes' viscous boundary layer of thickness


6 = 72 v/-


(1-2)








5

where v is the fluid kinematic viscosity. For room

temperature water at a frequency of 10 Hz, this viscous

boundary layer is approximately 1.7.10-2 cm. The

corresponding thermal boundary layer thickness is about


6th = 6//Pr (1-3)

where Pr is the Prandtl number. Note that both 6 and 6th

decrease in thickness as the oscillating frequency w

increases.

In the theoretical analysis, it is always assumed

that a fully developed velocity profile of the oscillating

flow exists within the pipes. At high frequency w, this

flow consists essentially of a slug flow over most of the

fluid core bounded by a thin boundary layer of width S as

discussed by Uchida [36]. Neglecting end effects, the fully

developed oscillating laminar velocity profile in pipes, due

to a periodic axial pressure gradient, is found to be [19]



U(S,t) = Uof(C)eiW (1-4)


where Uo is a representative velocity, r = r/R1 is again the

dimensionless radial distance, f() the velocity shape

function, and w the angular velocity of the oscillatory

flow. The explicit form of f(r) is




af() = 1 (1-iaC ) (1-5)








6

where A = R12 8p/a8x/pvUo is the nondimensional pressure

gradient maximum acting along the capillaries, a = R1J/wlv is

the Wormersley number measuring the ratio of inertia to

viscous forces, v is the fluid kinematic viscosity of the

carrier liquid, and p is the fluid density. This velocity

profile will reduce to the familiar Poiseuille parabolic

shape as the angular frequency w becomes small, while at

moderate frequency, f() has the shape demonstrated by

Uchida [36].

Another new term commonly used when dealing with

oscillating flow is the cross-stream averaged dimensionless

tidal displacement AX which can be mathematically defined as


4Uo 1i
AX =I4 J~ff(s)dr (1-6)



and on integration, yields [18]



|ap/ax| 2
AX = -1/2 1 + F() (1-7)




where the complex function F(a) has the form



F(a) =i J( -ia)


with the prime denoting differentiation. Using the

definition of Kelvin function: Jo(,/-ia) = ber a + i bei a




















































(


u---


OD




00000


+ ++ + 4.


co r-

0d


Y- 1
, I


iCOM 'cT nt)


o d


LL.


-----


30
d1-









8

[1], the complex function F(a) (Fig. 1-1) can be further
written as:



F(a) = Fr(a) + i Fi(a) = i ber a + i bei'a (1-8)
ber a + i bei a


This dimensionless tidal displacement is related to the

maximum of the periodic pressure gradient via [18]



|2 a << 1
AX 1/2 p2 J2 (1-9)
1- ---- a >> 1




Apparently, for any fixed tidal displacement Ax and

oscillating frequency w, the axial pressure gradient

required is proportional to the inverse square of the

Wormersley number when a is small; however, it is

independent of a when a is very large. This also implies

that the exciting axial pressure gradient |ap/axI is

approximately proportional to the fluid kinematic viscosity

v and inversely proportional to the square of the pipe

radius R1 when Wormersley number is very small (it happens

only at low w, small R1, and large v, for example, oil),

while it is almost independent of the fluid kinematic

viscosity and the pipe radius when the Wormersley is very

large (it happens only at very large R1, high w, small v,








9
for instance, a liquid metal). Note, if the tidal

displacement AX is fixed, while allowing the oscillating

frequency to change, the axial pressure gradient |ap/ax| can

become very large when the oscillating frequency becomes

large. This is mainly due to inertial effects and not so

much due to viscous drag forces which dominate the

oscillatory flow at small Wormersley number.

With the definition of the above quantities, we are

now in a position to explore the details of the enhanced

axial heat transfer in oscillatory flow within pipes. It is

assumed that a constant temperature gradient exists along

the pipe in the axial direction and that a very large time-

dependent radial temperature gradient variation is

superimposed. When the fluid moves towards the cooler side

(we term this the positive stroke), the hotter fluid within

the pipe core which is initially brought into the pipe from

the hot reservoir produces a large radial heat flow via

conduction to the cooler portions of the fluid within the

Stokes' boundary layer and to the cooler pipe wall; while

during the negative (or reverse) stroke, i.e., when the

fluid moves towards the hotter side, the higher temperature

in the boundary layer and the pipe wall conducts the heat

back into the cooler fluid core. This coupled radial heat

conduction with an axial convective transport leads to an

enhanced axial heat flux along the entire length of the

pipes.








10

Further, from the system point of view, the heated

fluid near the cold reservoir will eventually be ejected

into the cold reservoir and mixed there with the lower

temperature liquid. Contrarily, near the hot reservoir side

the fluid within the pipes which has been cooled during each

positive stroke is pushed out into the hot reservoir and

mixes with the higher temperature liquid. This process of

thermal pumping is what leads to a time-averaged heat flow

from the hot reservoir into the cold reservoir. It differs

essentially from the working principle of a normal pump.

For a normal pump one can draw the analogy with transport of

a one-way vehicle which transports passengers as well as the

carrier from one point to another. For the thermal pump one

can draw the analogy with a two-way busline which

periodically loads and unloads the passengers (heat) from

the hot reservoir to the cold reservoir, heat can be

continually transferred, and the carrier, in the time-

averaged sense, does not move. This property is

particularly important for those systems where a large

amount of heat transfer is needed while the working fluid is

required to remain in the system (as in nuclear reactors).

Note that the axial heat conduction, in general, is assumed

to be very small in this thermal pumping process compared

with the enhanced axial heat flow [16].

Apparently, the heat transport rate in thermal

pumping is governed by both the thermal properties of the








11
working medium and pipe wall and the characteristics of the

oscillatory motion of the fluid. The enhanced axial heat

flow does increase with increasing oscillating frequency as

this thins out the boundary layer and leads to an

accompanying increase in the radial temperature gradient.

This observation holds only as long as the thermal

properties of the liquid and the wall are compatible. If

the molecular conduction of the fluid in the radial

direction is very small, then even high frequency

oscillatory motion will not produce a large increase in the

rate of enhanced axial heat flow. This is because such a

system fails to supply the "passengers" enough time to be

loaded onto the "bus" and to be unloaded from the "bus".

Obviously, the system just wastes energy. On the other

hand, if the molecular conduction of the fluid in the radial

direction is very large, but the frequency of the

oscillatory motion is low, once again one can not expect

that there will be an efficient enhanced axial heat transfer

between two reservoirs because the "bus" is now moving too

slowly.

The above observation can be confirmed by an analysis

of the performance of a water-glass combination (i.e., water

is the working fluid medium within a glass tube) and of a

mercury-steel combination. For the former, it is necessary

to employ rather small diameter tubes and low frequency

oscillations with large tidal displacement, for it has








12
relatively poor heat diffusive properties as compared with

the mercury-steel case. This small tube diameter-lower

oscillation frequency set-up is necessary in order to ensure

that there is sufficient time to transfer the excessive heat

content of the bulk water core to the tube wall during the

positive stroke of each period and to permit the transfer of

excess heat from wall to the cooler core fluid during the

negative stroke. Otherwise, the water would just carry a

portion of its heat content back and forth in the pipe and

the condition for achieving optimum enhanced axial heat

transfer could not be met. For the mercury-steel case, one

can chose relatively large pipe diameters and higher

oscillation frequencies with smaller tidal displacement

because of the higher thermal conductivity. This assures

that only very short times are needed to exchange the heat

between the core of the fluid and the wall. It should be

pointed out that the suggestion of using smaller tidal

displacement is purely due to the mechanical considerations

and that one always tries to keep Ax as large as possible in

order to produce large axial heat flows.

From the above discussion, it can be concluded that

the process of enhanced heat transfer via oscillatory

pumping requires a precise tuning of parameters governing

the enhanced heat transport. Indeed there is an expected

"tuning effect" as discussed in references [16, 20, 21].








13

The tuning effect is a very important concept in the

study of the presently considered heat transfer process. It

shows that there will be an optimum combination between

thermal properties of the working medium and wall and the

characteristics of the oscillatory motion. The qualitative

aspects of the tuning effect have been observed earlier for

both the case of a flat channel and that of the cylindrical

pipe [15, 16]. From Fig. 4-38, one can see that an optimum

for axial heat transfer occurs only at or near the tuning

point which depends on the oscillating frequency and the

thermodynamic properties of the fluid and wall. As has been

pointed out by Kurzweg [20], in order to obtain the optimum

enhanced heat transfer one has to carefully select suitable

values for the pipe size, the material of pipe and the

working medium as well as the manner of oscillatory motion.

The nondimensionalized enhanced heat diffusivity is

defined as

Ke
A = w~X2 (1-10)


where Ke = 0/ypc is the coefficient of enhanced heat

diffusivity, 0 is the axial heat flux, 7 is the time-

averaged axial temperature gradient, p is the density of the

fluid, and c is the specific heat. One can show that this

nondimensional enhanced heat diffusivity is a function of

both the Wormersley number and the Prandtl number [21] and

hence that the dimensional axial heat diffusivity Ke is a








14
function of the tube radius R1, the oscillating frequency w,

the kinematic viscosity v, and the square of the

dimensionless tidal displacement AX. This can be explained

from the fact that the radial heat flow is proportional to

the product of the representative radial temperature

gradient YAX/Sth and the surface area per unit depth of rAX

available for cross-stream heat transport.

The use of large tidal displacement is always

beneficial in the enhanced axial heat transfer within

oscillating pipe flow. However, in order to avoid the

direct convective net mass exchange between the two

reservoirs, the tidal displacement must be limited to less

than about one half of the pipe length.

As has already been predicted by theoretical studies

and will be confirmed by the present numerical simulations,

the axial heat transfer will be further enhanced if the

rigid surface (part of the rigid wall with finite thickness)

has a non-zero thermal diffusivity and hence heat storage

capability.

Note that the existing considerations are restricted

to laminar flow. Turbulent flow conditions can occur in

oscillating pipe flow at higher values of wAx2/,v [26, 27]

and apparently would destroy the assumptions of the current

analytic model of the thermal pumping process. Fortunately,

the condition for optimum enhanced heat transfer in such

oscillating pipe flow obtained at the tuning point requires









15
very slender pipes, such that the Reynolds number is usually

small enough so that the oscillating motion falls within the

laminar range [22].

The theoretical aspects of the oscillatory enhanced

axial heat transfer process have been developed much further

than its experimental and numerical counterparts. The

theoretical predictions are quite limited and consider only

cases under certain simplifying assumptions [8]. Numerical

work is necessary in order to not only to examine the

correctness of the theoretical analysis but also to further

the development of advanced enhanced thermal pumping

devices. Numerical studies are not only fast, economical

and accurate, but they also offer a handy way to access

complex geometries which can not be handled analytically.

It is the purpose of this study to extend the

analytic work on thermal pumping by a detailed numerical

study. We intend first to examine some transient problems

of axial heat transfer in oscillatory pipe flow, such as the

development of the velocity profile at various Wormersley

numbers in contrast with those of reference [36], where only

several special cases with intermediate Wormersley number a

were discussed, and to examine numerically the relationship

between the tidal displacement and the required

corresponding strength of the periodic pressure gradient as

a function of Wormersley number a and of tidal displacement

Ax.








16

Next, we examine the build-up process of the

temperature distribution in a pipe which connects a hot

reservoir at one end to a cold reservoir at the other end

and see whether there actually exists a constant time-

averaged temperature gradient along the pipe axis when the

final periodic state is eventually reached. Note that a

time-averaged linear temperature distribution along the

axial direction is an essential assumption in the existing

theoretical studies.

The third part of this investigation which forms the

main effort, is a computer-aided numerical simulation of the

thermal pumping technique, including an investigation of the

variation of the enhanced axial heat flux versus the tidal

displacement, the variation of enhanced axial heat flux

versus different Wormersley numbers, and a study of the

variation of heat flux versus different Prandtl numbers. It

also includes a study of the influence of wall thickness and

pipe diameter as well as the change of the fluid properties

on such an enhanced axial heat flux and an examination of

the tuning effect in the conducting wall case. Further, we

compare axial heat transfer in oscillating flow with that in

the steady flow, and compare quantitatively this enhanced

heat flux with that induced by the pure axial molecular

conduction under various Wormersley numbers and different

fluid-solid wall combinations. Finally, the author examines

the phase lag phenomenon in oscillatory flow in order to








17
further reveal the mechanism whereby this new heat transfer

technique functions.

The computational work presented here was completed

on the Vax-11/750 electronic computer in the Department of

Aerospace Engineering, Mechanics, and Engineering Science,

and on Vax-11/780 in the Center for Instructional and

Research Computing Activities, University of Florida,

Gainesville, Florida. The numerical approach used, for

solving the presently considered heat transfer problem,

employed a second-order Crank-Nicolson scheme and an

Alternating Direction Implicit method (ADI) with Thomas

algorithm. A Fortran computer code named ETP (Enhanced

Thermal Pumping) was developed to implement all the

calculations.















CHAPTER II
FORMULATION OF THE PROBLEM


A single tube with inner radius R1, outer radius R2

and length L (such that L >> R1 and R2) connects a large

reservoir of hot fluid at temperature T = Th at one end, to

another large reservoir of the same fluid at cold

temperature T = Tc at the other end (Figs. 2-1 and 2-2).

The tidal displacement Ax and the oscillating frequency w

are adjusted to assure that the turbulent flow does not

occur. The tube is oriented in such a manner so that the

effect of gravity on the oscillatory motion of the fluid in

the pipe is negligible. Water and mercury are employed as

the working mediums. They are understood to be Newtonian

and incompressible fluids. The variations of thermal

properties of the fluid with temperature during the heat

transport process are assumed to be negligible. With these

assumptions the problem of enhanced time-dependent heat

transfer induced by a simple harmonic oscillatory laminar

fluid motion in a very long circular pipe can then readily

be formulated.

Governing Equations

The use of a very slender capillary circular tube

with constant cross-sectional area and neglecting end









































0

-J
LL

LU 3-
I-0
W ~

0
I



















Z

)-8


I-
CL~
I0


Qz


-L



0
0
-o
0U
u,

















o




-J5

O
0
C..


Q


(LV
a \.


\0




1r



>0

C) a
I o-

cn



0r
w -1



4-
EC


01
















4-1
r10








20

effects insures a laminar axisymmetrical one-dimensional

time-dependent motion. It is convenient to employ

cylindrical coordinates for this problem and we denote the

coordinate in the axial direction of the tube by x, and the

radial direction by r. The axial velocity can be taken to

be independent of x and, in order to satisfy the

requirements of continuity the other velocity components

must vanish. We shall further assume that the pressure

gradient induced by moving the piston (Fig. 2-1) is harmonic

and has the form [29]

1 ap
S= A cos wt (2-1)


where A = l/p ap/ax| is a constant which measures the

maximum pressure gradient existing along the x-axis.

Clearly this equation implies that we are now dealing with a

time-dependent sinusoidal pressure gradient which is

constant over the pipe cross-section at any instant and the

pressure varies linearly along the x-axis. The simplified

Navier-Stokes equation for this problem is [36]


aU v 8 aU
t~ Acos(wt) + B [r(a)J
0 < r < R1
(2-2)

where U(r,t) is the time and radially dependent axial

velocity component.








21

The corresponding temperature T(x,r,t) of the fluid

within the pipe is governed by the heat conduction equation

[15]



aT OT 1 a aT a2T
at- -U 8x + Kr r (r + a x2


0 < r < R1

(2-3)

where R1 is the inner radius of the tube and xf the fluid

thermal diffusivity which is related the thermal

conductivity kf by

kf
= pc

Here p is the fluid density, and c is the specific heat.

Note that the viscous heating term has been neglected in

equation (2-3) since it is very small for most experimental

conditions; this is justified provided one does not deal

with very high Prandtl number fluids such as oils. The

temperature in the wall can be determined from the solution

of

8T 1 aa T a2T
at w[I r arr-8ar) + ~ax\


R1 < r < R2

(2-4)
where Kw is the thermal diffusivity of the conducting wall

in R1 < r < R2 and is defined by









22
kw
PwcW

where kw, py, and cw, are the wall conductivity, density,

and the specific heat, respectively. By introducing the

following non-dimensional terms, equation (2-2), (2-3) and

(2-4) can be treated more easily.

t x
t / x 6

r U
r* U* -
6 A/w

T
T =
Th Tc


where 6 = J12i/ is again the fluid viscous boundary layer

thickness. The dimensionless governing system of equations

can then be written as

aU* 1 a2U* 1 au*
= cos(t*) + -2 [ r*2 + -w-

0 < r < R1

(2-5)
aT* OT* 1 82T* 1 aT* 82T*
at = CU* x + 2Pr + + a x


0 < r < R1

(2-6)
and

T* 1 a2T* 1 T* a2T*
a*= 2S + + ax r *

R1 < r < R2


(2-7)











where

Cf
Pr

'UCw
S -
Cw

and
-A
C= ~^
076


This nondimensionalization has some advantages in the

computing process to be carried out below. The

dimensionless velocity and its distribution over the cross-

section found from the momentum equation (2-5) are expected

to be universal for any Wormersley number a = R1Jw p/ and

any associated quantities, such as the tidal displacement

and Lagrangian displacement. Its final periodic form is of

the form given by (1-4). The dimensionless temperature in

the pipe is only related to the Prantdl number and the

dimensionless velocity, while that in the wall is related to

the ratio of wall heat diffusivity to the kinematic

viscosity, as seen from equation (2-6) and (2-7).

The governing dimensionless equations (2-5), (2-6),

and (2-7) are a set of second-order parabolic type of

partial differential equations expressed in cylindrical

coordinates. To solve this set of simultaneous equations, a

corresponding set of boundary conditions and initial

conditions are required. Since the velocity is assumed to

be a function of r and t only, just two boundary conditions








24
are needed for solving the momentum equation, while for the

temperature T(x,r,t) the heat conduction equations require

two boundary conditions in both the r and x directions as

well as compatibility conditions along the interface between

the fluid and the solid wall. It should be pointed out that

the initial conditions are less important than the boundary

conditions if one seeks only a final periodic state.

Boundary Conditions

The boundary conditions for the velocity are the

usual ones for viscous flow, namely, that the velocity is

zero at the inner surface of the wall (r* = Rl*). Also by

symmetry, the normal derivative of velocity at the axis is

zero. That is,



U*(R1 *, t*) = 0 (2-8)


and

8U*(0, t*)
ar = 0 (2-9)


The boundary conditions for the temperature depend on the

particular model investigated.

Model 1

In this model, it is assumed that the temperature at

each end of the pipe is equal to that of the connecting

fluid filling the reservoirs when the fluid moves into the

capillary tube (i.e., during one half of each cycle). The

























r(-

0

E
E


XO




0
3Wr









0
k II










a,





X
*-








26
end boundary temperature at x = 0 and x = L will take the

values of temperature of the adjacent pipe fluid elements.

With this model, we tried to simulate the real experimental

observations where it can be clearly seen that a fluid jet

exits the pipe during the outstroke, while the elements of

the fluid enter the tube in a funnel pattern during the

instroke. This observation shows that there is enough time

during each oscillating cycle to allow the fluid elements

within the exiting jet to fully mix with the fluid particles

which originally are in the reservoir before the next

instroke, provided the oscillating frequency is not too

high. The end temperature boundary conditions are here

taken as


T*h (when fluid enters pipe)
T*(0, r*, t*) =

T adj (when fluid leaves pipe)


(2-10)


[T*c (when fluid enters pipe)
T (L r t ) =

T*adj (when fluid leaves pipe)


(2-11)


where T*c, T*h are the nondimensional temperatures of cold

and hot reservoirs, respectively. While T*adj is the








27

temperature of the fluid elements which are adjacent to the

corresponding ends at a particular instant.

The thermal boundary conditions along the outer

surface r* = R2 of the pipe wall are taken to satisfy the

insulating wall condition, and along the axis of the tube

the temperature is assumed to meet the symmetric boundary

requirement, i.e., the radial derivatives of temperature

along axis are equal to zero. We thus have

9T*
r I R = 0 (2-12)

and

8T*
S* = 0 (2-13)
r =0

The compatibility conditions along the interface

between the fluid and the solid wall are that the radial

heat flux and temperature are constant across the interface.

That is


a8T* aT*
kf r = kw --r* at r = R1 (2-14)
fluid wall

and


T* = T* at r = R1 (2-15)
fluid wall


Since numerically the same nodes are chosen along the

interface, condition (2-15) will be automatically satisfied.










Model 2

Here it is assumed that a heat source rim of width 2b

is directly mounted on the interface between the solid wall

and the fluid at x = L/2, while two small cold rim sources

of width b each are mounted at x = 0 and x = L (see Fig. 2-

3). The wall thickness is assumed to be zero. This model

is intended to simulate the enhanced heat transfer process

of oscillatory fluid in an infinitely long pipe which is

heated and cooled by the alternative evenly distributed heat

and cold sources along the outer surface of a very thin wall

which possesses infinite heat conductivity. Apparently,

this geometry can be well simulated by the present model

with periodic boundary conditions. Mathematically the

boundary conditions are



T*(x*, Rl*, t*) = T*h X2* < X* > X3* (2-14)


0 < x* X1*
T (x*, R1*, t*) = Tc or (2-15)
X4 : x* x< L*


aT* X1* < x* < X2*
= 0 or (2-16)
X3* < X* < X4*

where Xl*, X2* X3*, X4*, X5* and X6* are nondimensional

coordinates of points which correspond to X = 0, b, L/2-b,

L/2+b, L-b, and L, respectively in Fig. 2-3. The periodic

end boundary conditions are given as
























U)





0



V H



0
U ^
:3






o H

O H II

V 0
0)
IWU



a








0C
Oar'
Pf-r








30

T*(0, r*, t*) = T*(L*, r*, t*) (2-17)


and


8T* aT*
ax X* ax *(2-18)
x 0 x =L




Model 3

In model 1, the boundary conditions for the heat

equation at X = 0 and L are somewhat artificial. The

temperature distribution in the vicinity of pipe ends may

not precisely match the constant temperature end conditions.

This is particularly true in the large tidal displacement

and/or very high oscillating frequency cases since a

temperature drop occurs at the hot end and a rise at the

cold end as fluid particles exit the pipe. This leads to a

discontinuity of temperature which will lead to unexpected

dispersion errors if an even order numerical method is used

[3, 24, 30]. To improve the end boundary conditions

presented in model 1, the following extended model has been

considered. It is assumed that an extension pipe which is

some 5 times the length of the central pipe is attached to

the original pipe used in model 1. The heat and cold source

are assumed to be located on the outer surface of the wall

in the extended sections of the pipe as well as at both ends

of the tube as shown in Fig. 2-4. This model is used to























0








IS
0


0






-.I
CH
o
U II



cr
0


u
>N





0













Ir









32

investigate the situation where the heat contained in the

jet is exchanged by pure heat conduction with the

surrounding fluid elements in each reservoir without any

convective mixing. In this case, the boundary conditions

along the outer surface of the wall are


T* = T*h



T = T*c


O s X* < X1*



X2* X3*


8T*
-ir^ 0


Xl X* < X2*


(2-21)


where Xl*, X2*, and X3* are nondimensional coordinates of

points which correspond to X = 5L, 6L, and 11L,

respectively, in Fig. 2-4. At both ends, we have


T *(0, r*, t*) = T*h


(2-22)


T* (L*, r*, t*) = T*c


(2-23)


As in the other models, the symmetry condition along the

axis of pipe requires that


aT*
-ar r=O= 0 (2-24)


and


(2-19)


(2-20)


and








33

and the radial flux and temperature continuity condition

along the interface are


aT* | T*
kf far li = kw w (2-25)
fluid wall

and



T* = T* (2-26)
T fluid wall


Here again kf and kw are the same as described in the

previous models.

Initial Conditions

The initial condition chosen depends on the problem

under consideration. If one's goal is to investigate the

periodic quasi-steady solution only, the initial conditions

chosen should be as close to the quasi-steady state as

possible so that a rapid convergence rate at low CPU time

cost is achievable. If one intends to study the transient

process, then various initial conditions should be supplied

according to purpose of the investigating cases selected.

For the velocity initial condition we choose for all our

studies



U*(x*, r*, 0) = 0 (2-27)



For the temperature, if the purpose of investigation is to

examine the build-up process of the periodic quasi-steady








34

state in the thermal field, the initial condition should be

taken as



T*(x*, r*, 0) = 0 (2-28)



However, for the other cases, in order to gain faster

convergence, the initial temperature can be assumed to have

a linear distribution in the axial direction, and thus have

the form


(T** T*h)xx*
T*(x*, r*, 0) = T*h + LL* (2-29)



where LL* is the dimensionless length within which the axial

linear temperature distribution is assumed to hold, and xx*

is the dimensionless distance measured from the "origin"

which is chosen only for the purpose of this linear

temperature initialization. Both the LL* and the origin

selected depend on the model considered. For Model 1 and

Model 3, LL* is equal to L*, and the origin is taken at the

left end (Model 1), or at the left intersection between the

central pipe and the left extension pipe (Model 3). However,

for Model 2, equation (2-29) is valid only for the right

half portion of the pipe as LL* is taken to be L*/2, and the

origin is chosen at the middle section of the pipe. The

initial temperature of the left half pipe can then be found

from the plane symmetric condition about the origin cross-

section.








35

Calculation of Tidal Displacement

An important quantity encountered in the study of the

enhanced heat transfer process in oscillatory pipe flow is

the tidal displacement, which is usually required to be

smaller than one half of the total pipe length in order to

avoid any convective mass exchange occurring between the two

reservoirs. It has already been defined in the introduction

chapter (1-6) and in the present computation the dimensional

tidal displacement takes the form


2 R1
AX = R12 f x(r,w/w)rdr (2-30)



where R1 is the inner radius of the pipe and x(r,w/w) is the

Lagrangian displacement of the fluid elements located along

a radius within the capillary tube at t = w/w. It is

assumed these elements are initially lined up at axial

position x = L/2 (Model 1 and Model 2) or x = 5.5L (Model

3), half way between the tube ends. Numerically, the

dimensional Lagrangian displacement at time t is computed

via the equation


t
x(r,t) = U(r,t)dt (2-31)
0

It is obvious that the dimensional Lagrangian displacement

x(r,t) is a function of both time t and the radial position

r. Note that, since the existence of phase lags between the









36
stimulating pressure gradient and the displacements vary for

different Wormersley numbers, in actual calculations, the

tidal displacement is equal to the sum of the absolute

maximum cross-section averaged Lagrangian displacement and

the absolute minimum cross-section averaged Lagrangian

displacement within a period. With the same non-

dimensionalizaton procedure used in the previous section for

the governing equations, the dimensionless Lagrangian and

tidal displacement can be written as


R*
AX = 2 f X*(r*,r/w)r*dr* (2-32)



The nondimensional Lagrangian displacement within the period

can then be computed by



X* (r*,t*) = U*(r*,t*)dt* (2-33)
*0

where


=
X (A/w2)

and

AX
AX =


It is pointed out that this dimensionless tidal displacement

Ax* differs from AX defined by equation (1-7) by a constant

1/2.









37

Effective Heat Flux

The time averaged total effective axial heat flow

within the pipe has the form



Qtotal = cpw dt I UTrdr (2-34)



where c is the specific heat, p is the density and again, w

is the oscillatory frequency. The per unit area effective

heat flow (termed heat flux) can then be written as



-= (2-35)
rR1


Equation (2-34) follows from the fact that pcUT is the

convective heat flux. The dimensionless total effective

heat flow can be written as


2n R1*
Q*total = cpw dt* I U*T*r*dr* (2-36)


and the dimensionless heat flux by



Q total
S= *2 (2-37)















CHAPTER III
NUMERICAL TECHNIQUES EMPLOYED


Equations (2-5), (2-6) and (2-7) which were derived

in Chapter II can not generally be solved analytically in

any simple manner. Therefore, it is necessary to seek a

numerical solution approach to the problem.

As is well known, in order to numerically solve a set

of simultaneous governing differential equations more

accurately and efficiently, an optimized grid system plays a

very important role. For the present purpose the grid net-

works were generated in the following way: within the tube,

a non-uniform grid (see Fig. 3-1) is clustered along the

radial direction in the vicinity of the inner surface of the

pipe where the larger gradients of velocity and temperature

are expected to be present, while along the x-axis the grid

is distributed uniformly except for those calculations

dealing with model 3.

The solution process was carried out in the

computational plane rather than directly in the physical

plane. Thus, a transformation which converts the governing

equations as well as the boundary condition plays an

essential role.








39









r












1). For Model 1
1). For Model I


wall


2). for Model 3


Fig 3-1 Grid Systems Used in the Numerical Simulations








40

The second-order implicit unconditionally stable

Crank-Nicolson scheme and the ADI method (2, 31] are

employed to break up the transformed equations into finite

difference form. A computer code designated ETP (for

Enhanced Thermal Pumping) has been developed for obtaining

the desired results. Details of this procedure will be

described in this chapter.

Transformation

The numerical solution of a system of partial

differential equations can be greatly simplified by a well-

constructed grid. It is also true that an improper choice

of grid point locations can lead to an apparent instability

or lack of convergence. Early work using finite difference

methods was restricted to some simple problems which can be

numerically solved in the physical domain. As experience

was gained, general mappings were used to transform the

physical plane into the computational domain as well as the

governing differential equations [31]. Such a grid

transformation technique brought to the numerical simulation

numerous advantages, and the computational work was no

longer restricted to a few simple geometries. A body-fitted

non-uniform grid in the physical domain could be used, which

retains the uniform spaced grid system features in the

computational domain [33, 34, 35]. However, several

requirements must be placed on the transformation: first,

the mapping must be one to one, second, the grid lines must









41
be smooth and have small skewness in the physical domain,

third, grid node point spacing should be small where large

numerical errors are expected (i.e., large solution gradient

regions) in the physical domain (Fig. 3-2).

In present study, the grid system in the physical

space is numerically generated by solving an algebraic

equation which clusters the grid lines in the region where

large gradients of the physical quantities are expected so

as to gain higher resolution of these physical quantities.

Because of the simple pipe geometry, the requirements of

grid smoothness and orthogonality are not a serious problem.

For the present purpose, the transformation is a simple one

which transforms a non-uniform grid network in the physical

plane onto a uniform grid system in the computational plane.


x* = x*(e)

r* = r*(s) (3-1)

t* = t*(r)


The inverse transformation can be found as


= ( (x*)
S= P7 (r*) (3-2)

r = r (t*)

where x*, r* are the dimensionless coordinates and t* the

dimensionless time in the physical domain, ,. 1 are the

transformed coordinates and r is the transformed time in the




























3 C
i-'-
0
e
0



0


-H
o
0
4V

4.





































41 0
r-.
u



































0.





1-I


'I











computational domain. With this transformation, the

transformed governing equations have the form:


au a2U* au*
f(r) + a()-) g-- + b(,q) a (3-3)



aT* BT* a2T* 8T* a2T*
r -TP1 + P2--a- + P3 + P4 (3-4)


and

OT* aT* 82T* 8T* 82T*
r w + w2 2-+ w3 + w4- (3-5)


where

at*
ar
f(r) = Or cos (t*) (3-6)


1 at* 1
a() = 2 Or (ar*/a,)2 (3-7)


1 at* 1 1 a2r*
b() -2 r*(ar*/,) (ar*/a,)3 (3-8)


and

at CU* 1 a2x*
P1 r tx*/a~ 2Pr(9x*7/a)3 + 2 + (3-9)


at* 1
P2 = a7 2Pr(ax*/aO2) (3-10)


at* 1 1 1 a2r*
P3 r 2Pr r(r (r*/8Oj) (ar*/)3 q2 (3-11)

at* 1
P4= r 2Pr (ar /a-)2 (3-12)









44

at* 1 -1 a2x*
W1 -ar 2PrW (xl 0 8 (3-13)
72Pr;(ax*/aa) a3 "


at* 1 1
W2 = 87 2Prw (ax*/a)2 (3-14)


at* 1 1 1 a2r*
W3= r 2PrW (r*(ar*/a8,) (ar*/a?)3 7J (315)


at* 1 1
W4 ~9 2Pr; (ar*/a8)2 (3-16)



As mentioned in chapter II, equations (3-3), (3-4) and (3-5)

are a set of second-order parabolic partial difference

equations in cylindrical coordinates. In addition, since

the oscillating flow is considered as incompressible in the

present studies, the momentum equation (3-3) can then be

independently solved at each time step. As a result, the

time-dependent update velocity found can then be substituted

into the heat conduction equation (3-4) as a known quantity

at the same time step level. Eventually, the coupled heat

conduction equations (3-4) and (3-5) are solved

simultaneously to obtain the temperature distribution both

in the fluid and in the wall at any time. To best solve

this set of equations in terms of accuracy and efficiency,

the proper choice of numerical technique and grid net-work

is dictated by an understanding of the physical aspects of

the problem.









45

The same transformation should be also applied on all

the boundary conditions proposed in the three different

models. For the temperature compatibility conditions along

the interface between the fluid and the wall (2-25) one has

the following transformed forms:


aT* 1 aT* 1
a17 (Or a) fluid a8q (ar /l) 'wall


or


3T* aT*
Slow= Ka wall (3-18)


where


kw ar*/8a pipe
Kg = -r*/a- wall (3-19)
kf ar*/8 1,|ll



To make the subsequent form of the corresponding

finite difference governing equations less cumbersome, the

superscript will be dropped from the variables T*, r*,

t*, and x*, and in addition, U* will be replaced by V.

In the process of deriving the finite difference

governing system equations, the second-order central

differences in the domain and forward or backward

differences along the boundaries or interface of the fluid

and wall have been employed at each nodal point.









46
Crank-Nicolson Method for Momentum Equation

The second-order accurate Crank-Nicolson method is

quite well known and widely used. It is an unconditional

stable, implicit scheme for solving the parabolic types of

partial difference equations. When the Crank-Nicolson

method is applied to equation (3-3), the finite difference

algorithm at a typical node k in the radial direction and at

the time step n assumes the simpler form:



(Bkk-1 + DV + A Vk+ n+ = C-
SkOk-l kk kk+lj k


k = 1, 2, .....km

(3-20)
where



C= f n+l+ fn Bk-1 + EkVk- Akk+ ]n (3-21)



1
Ak 4 ( 2 a + b )k (3-22)



1
Bk = 4 (2 a + b )k (3-23)



Dk = 1 + ak (3-24)


Ek = 1 ak


(3-25)









47
The boundary conditions along the interface (i.e.,

the inner surface of the wall, k = kmid) and the axis (k =

1) then become




V |k = kmid= 0


(3-26)
8V

8ri k = 1= 0



The initial conditions of the velocity at all nodal

points is taken to be zero.



k in = 1

k = 1, 2, 3, .....kmid


(3-27)



Equation (3-20) associated with the boundary conditions (3-

26) can then be written in the matrix form. This yields a

set of linear system algebraic equations which can be solved

in terms of the nodal values of velocity in the capillary

tube by using either an iterative method or a Thomas

algorithm at each time step. The explicit form is










A* n n+l n
D A V C
1 1 1 1
B2 D2 A2 V2 C2

B3 D3 A3 V3 C3





Bkd-1Dkd-1 Akd-1 kd-1 kd-1
Bkd Dkd Vd kd


(3-28)

where

A = A + B


kd = kmid-1

ADI Method for Axisymmetric Heat Equations

The governing PDEs (3-3), (3-4) and (3-5) are all of

the second-order parabolic type. Thus it might be suggested

that the Crank-Nicolson scheme used in solving the momentum

equation (3-3) can also be applied to the axisymmetric heat

equations (3-4) and (3-5), and one can then take advantage

of the tridiagonal matrix form while using this

unconditionally stable technique. However, when attempting

to use such a formulation, one immediately finds that the

resulting system of linear algebraic equations is no longer

of the tridiagonal type (3-20) but rather a non-tridiagonal

matrix system requiring substantial CPU time to solve. This

difficulty can be avoided by applying the unconditionally

stable Alternating Direction Implicit method (ADI),








49

developed by Peaceman and Rachford and Douglas in 1955.

According to this scheme, the entire solution process at

each time step is "split" into two portions, i.e., the first

half of solution processes for kth column (radial

direction), while the other half processes for the jth row

(axial direction).

With the ADI scheme, second-order central differences

are used to approximate the values of derivatives at each

nodal point in equations (3-4), (3-5). The finite

difference algorithm for those equations during the first

half of each time step for the jth column are then



IBP T + DP T + AP T n+1/2 -n
BPjjk jjk-l + DPjjk jjk + APjjk jk+1 n+/2 PXjj,k


k = 1, 2, ....kmid


(3-29)
and



BWjjk jj,k-l + DWjjkTjj,k + AWjjk jjk+1 n+1/2 Xj,k


k = kmid+1,...kmax


(3-30)

where the subscript jj is used to emphasize the specific

column currently to be computed. It can be seen that the

set of difference equations now is in the tridiagonal form








50

since the right-hand-side terms in the equations (3-29), (3-

30) contain only known values from the previous results and

the boundary conditions. These values can be computed by


PXjk = -CPjjkTjj-lk + EPjjkTjjk


(3-31)


and


WX -CW...T + EW. T -FW.. T
j,k = j,k jj-l,k j,k jj,k jjjj+l,J

(3-32)

The computational algorithm is implemented column by

column and the unknown value (Tj,k)n can then be solved by

either an iterative or a direct method. In order to do

this, equation (3-29) needs to be assembled into the

following matrix form.


PD1 PA1

PB2 PD2 PA2

PB3 PD3 PA3


* in


PBkd-lPDkd-1 PAkd-1

PBkd PDkd


n+1/2


PXl

PX2

PX3


Tkd-1 PXkd-1

j kd PXkd

j = 1, 2, 3 .....jmax


(3-33)


jj- FP,kjj+l,k n


T1

T2

T3








51

Similarly, equation (3-30) can be assembled into the matrix

form as seen in (3-34).


WD WA n T n+1/2 n
WDkd+1 WAkd+1 n kd+1 n+1/2 kd+1

WBkd+2 WDkd+2 WAkd+2 Tkd+2 WXkd+2





WBkm WDkm WAkm km WXkm-l

kWBmax WDmax j Tkmax j WX kmax


j = 1, 2, 3 .....jmax


(3-34)


It should be pointed out that the coefficients and/or

the right hand side terms in the row marked with symbol "*"

in those matrixes needs to be properly modified as well as

the limit of either subscript j or subscript k according to

the different boundary conditions in the various models

considered. For example, along the axis, the symmetry

condition requires that aT/ar = 0 and this could be

accomplished numerically by equating the temperature To and

T2 and by employing a new combining coefficient of PA1* =

PA1 + PB1 rather than the original PA1 in the first matrix

for evaluating the temperature distribution within the pipe.

In addition, if the temperature distributions along certain

parts of the outer surface of the wall are given, then the








52

limit of subscript k will be ended with k = kmax-l rather

than k = kmax shown in equation (3-34). The same arguments

are also applicable for another subscript j in the axial

direction.

By using a similar procedures for the second half of

each time step, the finite difference algorithm for the kth

row then becomes


ICP kT + GPj + 4 +FP T+ 1n+1 p Tt1/2
[CPjkkj-lkk k+ GPjkkTjkk + FPjkk+lkk n+1 +kk


j = 1, 2, ........jmax


(3-35)



CWjkk j-lkk + GWj,kkj lkk Wjkkj+ = kk / 2j,kk,


j = 1, 2, ........jmax


(3-36)
where the right hand terms can be computed by



pyn+1/2= -CP T + HP Tm API T n+1/2
P,kk j,kk j,kk-l j,kk j,kk j,kk j,kk+l1


(3-37)


kk (-CBjkkTjkk-1 + HWjkk jkk- Wjkk jkk+l 1


(3-38)











Similarly, to emphasize that the current calculation is in

the kth row, we use the symbol kk rather than k in above

formula and the matrix form of equation (3-35) for the kth

row can then be written as


PG2 PF2

PC3 PG3 PF3

PC4 PG4 PF4






PC
*


*I n+1/2


PG. PF.
1 Pjm-l jm-1
PCjm PG. jm


Tjm-1,

T. jI


y2

PY3

PY4






Pjm-l
PX j
jm .


n+1/2













k


k = 2, 3 .......kmid

(3-39)


and the matrix form of equation (3-36) becomes


WG2 WF2

WC3 WG3 WF3

WC4 WG4 WF 4
4 4 4


* In+1/2f


WC. WG. WF
jm-1 jm-1 Wjm-1
WC.jm DG
3m jm ]


T.
jm-1
T.


N+1 WY2

WY3

WY4






WX
jm-1

k Wjm


k = kmid+l, kmid+2.......km


(3-40)


N+1/2













k








54

As indicated above, the symbol shows that the

coefficients in that row as well as the right-hand-side

terms need to be properly modified corresponding to the

different boundary conditions. The matrix terms in

equations (3-29), (3-30), (3-31) and (3-32), (3-35), (3-36),

(3-37) and (3-38) can be computed by


1
APj -( P3 + P4)jk


13 -
BPjk (- --P3 P4).


CPj,k


1
P1 P2)jk


DP k. = (1 + 2*P4)j k


EPk = (1 2*P2)j k



FPk= -( -P1 + P2)jk



GPjk = (1 + 2*P2)j k



HPjk = (1 2*P4)jk
j,kjk


(3-41)


and


1
AW = -( WS + P4)'
j,k 2(j-W + j,k


1
= (-2-W3 W4) jk


1
CWjk = ( 2 W1 W2)j


DWjk = (1 + 2*W4) jk


EWj = (1 2*P2)j k



FWjk = -( 2 W1 + W2)j k



GWj = (1 + 2*W2)j k



HWjk = (1 -2*W4).k


(3-42)










In order to rewrite the equation in a simple form we define

the following linear operator :


LY [P]j,k = [BP, DP, AP]j k



LY [W]j,k = [BW, DW, AW]j,k

(3-43)


LX [P]j,k = [CP,
= Iicp


GP, FP]j k


LX [W]j,k =[CW, GW, FWk]jk


and column vectors


(TY)j,k = (Tjk-' Tjk' Tjk+l


(3-44)


(TX)j,k = (Tjlk Tk' Tj+l,k )


Equations (3-29) which are in the radial direction can be

simply written as


LY (P]
j,k


LY [W
j,k


n+1/2
(TY}
j,k



n+1/2
(TY}
j,k


j = 1, 2 ..........

k = 1, 2, .........


= PX
j,k


= WX
j,k


(3-45)


(3-46)









56
Equations (3-33) and (3-34), which are in the axial

direction, then assume the following simple forms:


n+1/2
LX [P]
j,k





n+1/2
LX [W]
j,k


n+1
(TX)
j,k





n+1
(TX)
j,k


n+1/2
= PY
j,k





n+1/2
= WY
j,k


(3-47)







(3-48)


j = 1, 2, ..........

k = 1, 2, .......


where the right side terms can be estimated by


PX
j,k


= (-CP, EP, -FP)


j,k


(TX)


j,k


j,k


j,k


= (-CW, EW, -FW)





= (-BP, HP, -AP)


(TX)
j,k j,k


(3-49)


j,k


{TY)


j,k


WY = (-BW, HW, -AW) (TY)
j,k j,k j,k









57
It is seen, as the result of the ADI "splitting"

procedure which has been employed in the algorithm

associated with different boundary conditions for the

various models, that only a tridiagonal system of linear

algebraic equations needs to be solved (i.e., during step 1,

the coupled tridiagonal matrix (3-33), (3-34) are solved for

each jth column of the grid points, while during step 2, the

coupled tridiagonal matrix (3-39) and (3-40) are then solved

for each kth row of grid points).

Once the periodic steady solution has been obtained,

we can calculate both the tidal displacement and the heat

flux at different locations within the pipe. The numerical

technique used here to integrate equations (2-38), (2-39)

and (2-41) for evaluating the tidal displacement and the

heat flux at each specified location can be obtained either

by the Trapezoidal rule (with end correction) or Simpson's

rule. Both numerical approaches are essentially fourth-

order methods.

Convergence Criteria
As is known, once the calculation work has started,

the time-matching process will be in the loop forever unless

a criterion can be derived that indicates when the goal of

the current computing work has been reached and further

solution-matching processes do not produce significant

increases in accuracy. Such a criterion depends on the

purpose of the calculation. If one's goal is to study the









58
transient process, one can indicate a time limit to pause or

quit the current computing work.

However, if one seeks a final periodic state

solution, one has to develop a criterion to test if the

solution can be considered acceptable.

It seems that the temperature is a good measure of

the accuracy of the overall solution process so that the

most efficient way is to apply the convergence test on the

temperature rather than on the velocity. For the present

purpose, two convergence criteria were alternatively used.

The choice of the criterion depends on what was the main

goal in the particular study case. If the main effort is to

observe the distribution of the quasi-steady temperature at

any phase angle wt within a period, the testing is done by

comparing the temperature residual, i.e., by inspecting the

averaged temperature difference of each nodal point at thee

same wt between adjacent periods. This can be written as



SRe (T11 T22 )2 k )/2
Resl = ---- < el
(JMAX)(KMAX)


j = 1, 2 .........JMAX

k = 1, 2, .......KMAX


(3-50)
where Tllj,k is the value of the temperature of node (j,k)

at ot in the current period, T22j,k is the value of the
), iste au o h









59
temperature of node (j,k) at the same wt in the previous

period, and EI is the convergence parameter.

If the goal of the investigation is to examine the

effective heat transfer, the convergence criterion is

established by computing the residual



Z ( 01 2 ) 1/
Res2 = E2 (3-51)
JSEC


j = 1,2,......JSEC


where JSEC is the number of cross-sections where the axial

heat flux was examined, E2 is the convergence parameter, and

41, 02 are the cross-section averaged heat flux in the

current period and previous period, respectively. The

summation is carried out over all sections (JSEC) where the

heat flux is computed.

The value of the convergence parameter can be

determined by a balance of the acceptable solution accuracy

and the cost of CPU time. In the current study both el and

e2 were taken between 0.001-0.01 as the residual of the

temperature or the heat fluxes for acceptable value.

Grid Generation
As is well known, the solution accuracy and efficiency

in a large degree depends on the grid system used in the

numerical calculation. A good grid is characterized by

small skewness, high smoothness and capability of high








60

resolution in the large gradient regions in the physical

plane. It has been shown that rapid changes of grid size

and highly skewed grids can result in undesirable errors

[25]. The success of a numerical simulation of a complex

thermal fluid dynamics problem does depend strongly upon the

grid system used in the computation.

In the present study, the grid mesh was generated with

the emphasis on the high resolution capability of its simple

geometric boundaries. The technique used here is one of the

algebraic schemes which cluster the grid lines near the

region desired [34), namely,



ASk = AS (1 + e) k-1 (3-52)




where AS0 is the minimum specified grid spacing next to the

wall or to some inner interfaces within the tube. The

parameter e is determined by a Newton-Raphson iteration

process so that the sum of the above increments matches the

known arc length between k = 1 and k = kmax. The grid

networks used in Model 1 and Model 2 are uniformly

distributed along the axial direction in both the fluid and

the wall; however, in Model 3, most grid points are

clustered near the central portion of the pipe so as to gain

higher resolution in the region of interest, while along the

radial direction a non-uniform grid was employed in all

three models. The minimum specified grid spacing ASo used









61

depends on the boundary layer thickness 6, namely, the

kinematic viscosity of working fluid and the oscillating

frequency. In the present study, it was found that a good

choice of this value is ASo = 0.056 for laminar flow cases

with a total of about 15-20 nodal points distributed along

the radius.














CHAPTER IV
NUMERICAL RESULTS AND DISCUSSION


The problem of time-dependent enhanced heat

conduction subjected to sinusoidal oscillations can now be

solved numerically for the boundary conditions appropriate

to a long capillary tube according to the various models

described in previous chapters. The computational tasks

fall into two categories: the first part is a numerical

study concerned with the characteristics of the oscillatory

pipe flow, which includes an investigation of the velocity

profiles, the Lagrangian and the tidal displacement

trajectories. The second part of this study includes a

thorough investigation of the final periodic state

temperature build-up process in oscillating pipe flow, the

periodic temperature distribution in the pipe and wall, and

a study of the relationship between the enhanced axial heat

flux and the tidal displacement. It also includes an

investigation of the tuning effect and a comparison of the

enhanced axial heat transfer with the corresponding pure

molecular axial heat conduction as well as the investigation

of the influences of the pipe radius and the wall thickness

on the enhanced axial heat transfer.








63

Properties of oscillating laminar pipe flow have been

analytically discussed by S. Uchida (36] for several

different Wormersley numbers, and these results offer a

useful reference for comparing with the present numerical

studies. However, in the area of thermal fields associated

with enhanced thermal pumping there is little information

available for comparison, except for some recent results of

Kaviany (13].

The present computational work was carried out on the

Vax-11/750 computer in the Department of Aerospace

Engineering, Mechanics, and Engineering Science, and on the

Vax-11/780 in the Center for Instructional and Research

Computing Activities, University of Florida, Gainesville,

Florida. Note that it is very time-consuming to build-up a

final periodic state, for example, if the grid size is

101x22 and the time steps are 2000 per period, it takes 7.5

minutes (CPU) with the VAX-11/780 machine or almost 1 hour

(CPU) with the VAX-11/750 machine to run only one period.

It usually takes 20-30 periods to reach the final periodic

state solution.

The process of selecting model size (i.e., total

nodal points) is a synthetic balance among the storage

requirement, the solution accuracy, and the cost of CPU time

for finding an acceptable solution. Such a selection allows

each case to be solved with a minimum of expense in

computing operation thereby making it possible to do the








64
large number of runs needed to obtain enough data points to

plot curves of the desired parameters.

The results discussed in this chapter are presented

for the purpose of illustrating the effects to be

encountered in working with the enhanced thermal pump.

These results should be effective in gaining insight into

some interesting features of this enhanced heat transfer

process. Since the analytic approximations [15, 16, 18, 22]

are effective in describing the flow and the heat transfer

aspects and to give considerable insight into the problem,

they are often used to compare with the numerical results

and should give good agreement when applicable.

Part 1. Oscillatory Pine Flow Features

In order to better understand the mechanism of the

enhanced thermal pump, it is necessary to examine the

mechanical features of the oscillating pipe flow. As

mentioned in previous chapters, the present interest in the

enhanced thermal pumping is confined to an investigation of

the central part of the slender pipe, thus the flow field

can be well approached by a 1-D time-dependent laminar model

which neglects the ends effects.

Velocity Profiles

Fig. 4-1 shows the numerically computed time-

dependent velocity profiles at different phase angles of the

exciting pressure when Wormersley numbers are equal to 1,

10, 100, and 1000, respectively. It can be clearly seen









65
















0 10












a 100







51




a =1000











Wt = ge s s Is z 20se 1e zIe I Z 0l 3W 33 336e X



Fig 4-1 1-D Velocity Profiles in Oscillating Flow
for Wormersley Number a = 1, 10, 100 and 1000








66

that the velocity profile at a = 1 presents a quasi-

parabolic shape at any instant within a period and is in

phase with the stimulating pressure gradient. However, at

higher frequency, for example, a = 100 and 1000, the

velocity profiles can be clearly separated into two regions:

in the vicinity of the wall the flow shows a typical thin

boundary layer, while in the region far away from the wall

the fluid moves as if it were frictionless slug flow. In

fact, within this core region the velocity distribution is

independent of the distance from the wall [29]. It can be

also seen that the phases of the velocity profiles at higher

frequencies cases are shifted about w/2 with respect to the














7 1 7 .


Fig 4-2 Velocity Profile (a=5) [29]










(R1-r) /6



0.2





0.1


0.5 V


0.5 V


2400 2100


1800 (R1-r)/6


V 0.5


(R1-r) /6
2700 300 3300 360


0.5 0


Fig 4-3 Magnified View of Velocity Profile Near Wall
(Wormersley Number a=10, H20, 6=0.014cm)


(Rl-r) /
T 00


0.1 I








68

stimulating pressure gradient. At the intermediate

frequency case (a = 10 of Fig. 4-1), the slug flow boundary

is not so evident, but one can still note a boundary layer

near the wall. The same pattern of the velocity profiles

associated with moderate Wormersley numbers can also be

found in the reference [36] where the velocity profile at

Wormersley number a less than 10 was presented (Fig. 4-2).

In order to better see the variation of the time-dependent

velocity profile within the boundary layer, a set of closer

view of the velocity profile with respect to Wormersley

number at a = 10 during phase intervals ot = 300 is plotted

in Fig. 4-3.

It should be emphasized that the solutions shown here

are under the assumption that the secondary velocity in the

radial direction is negligible compared to the axial

velocity component, namely, the non-linear inertial terms

are not considered in the governing equation (2-1). This is

a reasonable approximation for moderate oscillatory

frequency w except near the ends of the pipe. However, at

high Wormersley numbers, care should also be taken to avoid

violating another assumption, namely, that of

incompressibility and the concomitant condition that the

oscillating phase does not change between the tube ends.

High Wormersley number can be obtained either by increasing

the oscillating frequency w or by increasing the pipe

diameter for a given fluid. An oscillating flow can be








69

considered incompressible if Axw/2 < 0.05C, where C is the

speed of sound in the fluid. For mercury, C = 1360 m/sec,

one requires a = 27.1, when Ax = 20 cm and R1 = 0.1 cm. To

avoid a appreciable phase difference between the tube ends,

one requires that L/C << 2n/w. Both the restrictions are

met in the examples to be considered below.

The Lagrangian Displacements

An alternative interesting representation to the

oscillatory velocity field are the Lagrangian displacements

of the fluid elements at different radii within the pipe.

They have been plotted in Figs. 4-4 and 4-5 at time

intervals of wt = 300 for Wormersley number a = 0.1, 1.0 and

10. It is noted that since both pipe diameter and the

working fluid were fixed in this test, Figs. 4-4 and 4-5

represent the relationship between the Lagrangian

displacement and the oscillating frequency. The

trajectories plotted in Figs. 4-4 and 4-5 have been

normalized by A/w2, where A = 1/p dp/dxl is the amplitude of

the sinusoidal pressure gradient as defined in equation (2-

1) and w is the angular velocity. Suffice it here to point

out that for the lower frequency case (for example a 0.1

and 1.0) the Lagrangian displacement trajectory shows a

foreseeable parabolic pattern at any moment. Nevertheless,

the essential distinction between low frequency oscillatory

pipe flow and steady Hagen-Poiseuille flow is that in the

former, the Lagrangian displacement trajectories as well as



















a -0. 1



27(fy


a -1.0


= r/R1


.003


if = r/R1


Lagrangian Displacement for a = 0.1 and 1.0


Fig 4-4






























---AX -


= r/R1


cat=2100

X udt


Lagrangian Displacement at a = 10


Fig 4-5









72
the velocity profiles are periodic so that the fluid

particles do not translate axially upon time averaging,

while in the later case they will.

For intermediate Wormersley number (a = 10), the

trajectories of the Lagrangian displacement departs

considerably from the standard parabolic shape. This

phenomenon can be even more clearly seen in the a = 100

case. Evidently, the higher the oscillating frequency, the

thinner the boundary layer (6 = J2T/w).

Tidal Displacements

Figs. 4-6 and 4-7 demonstrate time variation of the

cross-section averaged dimensionless Lagrangian displacement

at Wormersley numbers in the range from a = 0.1 to 50. The

tidal displacement can be obtained by summing the absolute

maxima and the absolute minima of these curves. The

corresponding non-dimensional tidal displacements with

respect to Wormersley number from a = 0.1 to a = 12 are

listed in table 4-1 below


Table 4-1 Dimensionless Tidal Displacement
at Different Wormersley Numbers

a 0.1 1.0 2.0 3.0 4.0 5.0


AX 0.00246 0.14295 0.67303 1.20213 1.41160 1.49041


a 6.0 7.0 8.0 10.0 12.0


AX 1.56294 1.61879 1.66046 1.71900 1.76698

























a 0.


-.05 0.0 .05 .10

Dimensionless Cross-section Averaged
Displacement Versus Time
a = 0.1 1.0


Xave


-2. -1. 0. 1.


Dimensionless Cross-section Averaged
Displacement Versus Time
a 2 50


-.10


Fig 4-6


Xave


.15


Fig 4-7








74

It is noticed that for very small Wormersley number

cases (a < 0.5), the cross-section averaged displacement

varies like a sinusoidal function with respect to time (the

ordinate), and as Wormersley number increases, the cross-

section averaged displacements are no longer symmetric about

the ordinate but rather favor positive values of x. If the

Wormersley number is further increased, eventually, the

cross-section averaged displacement almost entirely lies on

the right half plane. A similar feature can also be seen in

Figs. 4-4 and 4-5. In fact, this just again shows the

existence of the phase shift in the oscillating pipe flow.

When the Wormersley number is small the cross-section

averaged displacements as well as the Lagrangian

displacements show a sinusoidal variation with respect to

time and this has w/2 phase lag with respect to the

stimulating pressure gradient which is assumed to be a

cosine function of wt with zero initial phase angle.

However, for higher Wormersley number, the phase lags

increase to almost r with respect to the phase of the

exciting pressure gradient. Note that the Lagrangian

displacements in Figs. 4-4 and 4-5 were computed by lining

up all the fluid particles on the plane x = 0 at wt = 0,

however, by using non-zero velocity at ut = 00 as shown in

Fig. 4-1. This implies that if the phase of the exciting

pressure gradient is taken as the base of the measurement,

it is generally not possible to assure the same phase to the








75

velocity and Lagrangian displacement as well as the tidal

displacement. For example, at phase of the stimulating

pressure gradient wt = 00, the corresponding velocity phase

may be w/2 and the Lagrangian displacement phase may be w.

In fact, this is just as true for very high Wormersley

numbers; however, the phase difference with respect to the

exciting pressure gradient phase is less than the value

shown above for low Wormersley number. One can well see

that if the phase lags of the Lagrangian displacement or the

tidal displacement is given, one can certainly re-draw the

diagram shown in Figs. 4-6, 4-7 by lining up the phase with

itself, and then an exact symmetric pattern of the tidal

displacement curve similar to that in the very small

Wormersley number case can be obtained. Unfortunately, the

phase lags are a function of the Wormersley number and they

are not known in advance.

In order to verify the numerical method used in this

study and to check the ETP code developed, a comparison of

the computed dimensionless tidal displacement to the one

using the analytic equations (1-7) and (1-8) given by

Kurzweg [18] for Wormersley number varying from 0.1 to 100

has been calculated and plotted in Fig. 4-8. The solid line

shows the analytic solution obtained by use of equations (1-

7) and (1-8), while the dashed line shows the results with

the ETP code developed in this study. The agreement is

quite good, particularly when the Wormersley is less than








76

3.0. However, at high Wormersley number the numerical

solution shows a very slight deviation from the analytic

solution. This deviation is believed due to an inaccurate

numerical integration over the cross-section using

relatively large time steps (our time steps per period in

the calculation were between 1000-2500). A comparison with

using 104 time steps per period for Wormersley number a = 10

was studied and shows some improvement. However, using such

small time steps in the present investigation is beyond the

capacity of the current VAX computer facility used. The

numerical error becomes particularly serious as the

oscillating frequency becomes large where the extremely thin

boundary layer requires more grid nodal points to resolve

the flow variables in the vicinity of the wall.

Fig. 4-8 shows that as the Wormersley number gets

large, the dimensionless tidal displacement tends to the

limit of 2.0, which agrees with the limit of 1.0 in the

analytical solution given by Kurzweg [18] for the reason

that the normalization parameter used in [18) is twice as

large as that in the present numerical simulation.

Fig. 4-9 shows the required stimulating axial

pressure gradient used in the present study for a pipe

radius R1 = 0.1 cm and water (v = 0.01 cm2/sec) taken as the

working medium versus the dimensional tidal displacement Ax

in cm for various Wormersley numbers (namely, oscillating

frequency). It is evident from these results that for fixed










77










--- Numerical Solution

-- --Analytic Solution


1.5


4)


0










r410
0 .
0
-H















10 100 101 1D2

Wormersley Number a



Fig 4-8 Relationship Between Dimensionless Tidal
Displacement AX and Wormersley Number a




































































0 1 1 I 1-I tI I I i
O O O 0


-

I
I



1-
0











ito
0)



0


44
01





0)
10





V
.0)




'4


(uoi) xv 4uaaovTTdsTG IgpTI


x
M
4P4a
0Q)

0C

r-
$4
>:1
0) 0

*c







tP
Q-H

CO








*4
to
S0
. hJ.

.1 (-
)u4














Irr








79

tidal displacement, the required axial pressure gradient in

the large Wormersley number case is orders of magnitude

higher than that in the small Wormersley number case. This

may eventually put some constraint on the use of very high

frequency in the enhanced heat transfer technique.

Fortunately, to meet the tuning condition, the required

Wormersley number in this case is of order 1. As already

discussed in the introduction, in order to gain the benefit

of axial heat transfer, the use of large tidal displacement

is always preferred. However, such an increase must be

limited by the requirement of no convective net mass

transfer occurring between two reservoirs and may also be

constrained by the ability of the device to withstand the

increase in the exciting pressure gradient.

Phase Lags

We are now in the position to study the phase lags of

the Lagrangian displacement in the oscillatory pipe flow.

Fig. 4-10 shows the phase variations (in degrees) along

radius for Wormersley number varying from a = 0.1 to 4.

Some numerical results are also shown in Tables 4-2 and 4-3.

All of the data shown in Fig. 4-10 and the tables have the

phase angle measured relative to the exciting pressure

gradient. Two features can be seen in Fig. 4-10; first, in

the core portion, the phase lags are almost equal to wr/2

when the Wormersley number is small, while the lags are

almost w when the Wormersley number is large, and second,









80

the phase lags vary along radius, especially in the boundary

layer. It is such phase lags that allow the existing

temperature gradient in the very thin boundary layer of the

oscillating pipe flow to act as region of temporary heat

storage. It absorbs heat when the temperature of the core


Table 4-2


Nodal
point
K


1


2


3


4


5


6


7


8


9


10


Phase Lags Along Radius
(Working Medium: H20, Ax = 10 cm)


a = 0.1


q=r/R1


0.0000


0.1667


0.3333


0.5000


0.5833


0.6667


0.7500


0.8333


0.9167


1.0000


phase


90.63


90.63


90.63


90.63


90.63


90.63


90.63


90.63


90.63


0.00


qr=r/R1


0.0000


0.0366


0.0957


0.1913


0.2593


0.3458


0.4559


0.5958


0.7737


1.0000


phase


102.06


102.06


102.06


101.88


101.70


101.52


101.16


100.62


99.72


0.00


a = 2.


q=r/R1


0.0000


0.1095


0.2164


0.2836


0.3620


0.4535


0.5603


0.6849


0.8303


1.0000


phase


134.28


134.10


133.38


132.66


131.76


130.50


128.70


126.00


122.58


0.00









81

Table 4-2 -- continued

Nodal a = 3 a = 4 a = 7
point
K q=r/R1 phase P7=r/R1 phase v,=r/R1 phase


1 0.0000 170.83 0.0000 179.82 0.0000 181.62


2 0.0917 170.47 0.1242 179.28 0.0970 181.62


3 0.1818 169.75 0.2370 178.02 0.1892 181.62


4 0.2703 168.31 0.3395 176.04 0.2769 181.44


5 0.3573 166.51 0.4324 173.70 0.3604 180.72


6 0.4427 164.00 0.5168 171.00 0.4397 179.28


7 0.5267 161.12 0.5935 168.12 0.5152 177.12


8 0.6091 157.88 0.6630 165.24 0.5869 174.06


9 0.6901 154.64 0.7262 162.18 0.6552 170.28


10 0.7696 150.69 0.7835 159.12 0.7201 165.78


11 0.8478 146.73 0.8356 156.06 0.7818 160.74


12 0.9246 142.42 0.8828 153.18 0.8406 155.16


13 1.0000 0.00 0.9257 150.30 0.8964 149.04


14 0.9646 147.60 0.9495 142.56


15 1.0000 0.00 1.0000 0.00











Table 4-2 -- continued


Nodal
Point
K


1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


a = 10


q7=r/R1


0.0000


0.0884


0.1740


0.2566


0.3364


0.4135


0.4880


0.5599


0.6295


0.6966


0.7616


0.8243


0.8849


0.9434


1.0000


phase


180.54


180.72


181.08


181.80


182.52


182.70


181.80


179.64


175.68


170.46


163.80


156.06


147.24


137.52


0.00


a = 14


a = 20


v,=r/Rl


0.0000


0.0834


0.1647


0.2440


0.3215


0.3970


0.4709


0.5428


0.6130


0.6815


0.7484


0.8136


0.8773


0.9394


1.0000


phase


179.46


179.10


180.18


184.50


190.25


192.41


189.53


183.06


174.07


164.00


153.20


141.70


130.19


118.32


0.00


,q=r/R1


0.0000


0.1635


0.3032


0.4227


0.5247


0.6120


0.6865


0.7502


0.8046


0.8512


0.8909


0.9249


0.9540


0.9788


1.0000


phase


179.10


180.54


185.93


182.34


186.65


184.14


184.14


181.98


178.38


173.71


167.95


161.12


153.57


146.01


0.00












83

Table 4-3 Comparison of Phase Lags With Different
Working Mediums (AX = 10)


a = 1.


Phase Lags


Nodal
Point

K



1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


q =r/R1



0.0000


0.0161


0.0366


0.0626


0.0957


0.1378


0.1912


0.2593


0.3458


0.4559


0.5957


0.7737


0.0000


a = 5.


Water


102.06


102.06


102.06


102.06


102.06


101.88


101.88


101.70


101.52


101.16


100.62


99.72


0.00


Phase Lags


Mercury


102.14


102.14


102.14


102.14


102.14


102.14


102.14


102.14


101.78


101.42


101.06


99.98


0.00


qi =r/R1



0.0000


0.0714


0.1429


0.2143


0.2857


0.3571


0.4286


0.5000


0.5714


0.6429


0.7143


0.7857


0.8571


0.9286


1.0000


Water


181.62


181.44


180.90


179.64


178.02


176.04


173.34


170.28


166.50


162.18


157.50


152.10


146.16


139.86


0.00


178.02


175.86


173.37


170.11


166.51


162.20


157.52


152.13


146.37


139.90


0.00







































(Z

Sa)




00




SQ)


4O4


0
4-
> 0





l(-l


snlpea ssaouosuawjju


I I I 0 0
a e as
O O a a a O








85

slug flow is higher than that of the boundary layer, and it

releases heat as the core temperature is relatively lower.

This large temperature gradient enhanced by the existing

velocity phase lags allows a large amount heat to be

conductively transferred radially within a very short time

and subsequently to be transferred axially by a convective

coupling.

Part 2. The Enhanced Heat Transfer Investigation

We have examined some mechanical characteristics in

oscillating pipe flow, and compared the computed solution of

the velocity field with Uchida's solution. The results for

the velocity profiles are in good agreement. Nevertheless,

before a detailed examination of the thermal field, it is

first necessary to test the current developed ETP code when

applying the temperature equations (2-6),(2-7). It is seen

that the energy equations strongly depend on the velocity

distribution and its build-up process, so that one can use

analytic periodic velocity state (Eqs. 1-4 and 1-5, with no

build-up process), and the computed velocity (with build-up

process) to verify the correctness of the resulting thermal

variables. The enhanced heat flux is a function of both

velocity and temperature (Eqs. 2-34 and 2-35) and was chosen

for a comparison of the analytic and numerical results of

the problem. Part A of table 4-4 shows the results of the

computed enhanced axial heat flux as well as the axial

conduction heat flux when using the analytic velocity









86

Table 4-4 The Comparison of Enhanced Heat Flux Using
Numerical Velocity with Heat Flux Using Analytical Velocity
(Model 3, Water-Glass, Pr = 7.03, a = 3)

A. Heat Flux Using Analytic Velocity


(cm) (w/cm2oK)


0.9839 0.0230


2.9519 0.1855


4.9194 0.3864


6.8868 0.5810


7.8847 0.6775


9.8390 0.8638

B. Heat Flux

1.0033 0.0234


3.0101 0.1870


5.0164 0.3899


7.0226 0.5790


8.0402 0.6746


10.0328 0.8597

* Model 3, Water-Glass
* .......enhanced axial


(w/cm 0K)


0.0149


0.0132


0.0098


0.0075


0.0067


0.0055


(w/cm -K)


0.0811


0.0724


0.0541


0.0415


0.0369


0.0302


(w/cm4 K)


0.0238


0.0213


0.0159


0.0122


0.0109


0.0089


Using Computed Velocity


0.0149


0.0133


0.0097


0.0076


0.0068


0.0055


heat flux


0.0811


0.0728


0.0533


0.0420


0.0374


0.0302


0.0232


0.0206


0.0155


0.0117


0.0104


0.0085


* Of .....axial heat flux by conduction in fluid
* Ow.....axial heat flux by conduction in pipe wall
* P = O/Ax2








87

(Eq. 1-4) for a = 1, Pr = 7.03, pipe radius R1 = 0.1 cm and

glass walls of thickness AR = R2- R1 = 0.05 cm (which is

almost equal to the boundary layer thickness 6 = 0.047 cm

for this case) with Model 3. There were 22 points

distributed along the radius (13 points in fluid), and the

smallest distance next to the wall is equal to 0.166, while

along the axial direction 101 point were used (60 points in

the central pipe section). Thirty period runs were

considered in order to insure that the final periodic state

was achieved. An exception was the Ax = 1 cm case where

only 23 period runs were computed. Under the exact same

conditions described above, a calculation using the

numerical velocity with build-up process was tested. Part B

of table 4-4 gives the axial heat flux results using the

corresponding numerical generated velocity profile. A

comparison of results shown in Table 4-4 is quite good. One

sees that the disagreement for most terms is less than

0.001. Note for Ax of 10 cm, the ratio O/of is about 15 and

gets larger as Ax is increased further.

Periodic Temperature Build-up in Thermal Pumping Process

The purpose of this facet of the present study was to

examine the temperature field build-up pattern and to

determine how long before the final periodic temperature

state can be reached. The author also intends to verify

whether or not there indeed exists a time averaged linear

axial temperature distribution within the thermal pump at