Fluid oscillations and enhanced heat transfer in a parallel plate channel bounded by end reservoirs maintained at differ...

MISSING IMAGE

Material Information

Title:
Fluid oscillations and enhanced heat transfer in a parallel plate channel bounded by end reservoirs maintained at different temperature
Physical Description:
xviii, 120 leaves : ill. ; 28 cm.
Language:
English
Creator:
Zhao, Alex X., 1949-
Publication Date:

Subjects

Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1990.
Bibliography:
Includes bibliographical references (leaves 116-118).
Statement of Responsibility:
by Alex X. Zhao.
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 - 001683420
notis - AHZ5391
oclc - 25034727
System ID:
AA00003339:00001

Full Text









FLUID OSCILLATIONS AND ENHANCED HEAT TRANSFER
IN A PARALLEL PLATE CHANNEL BOUNDED BY END RESERVOIRS
MAINTAINED AT DIFFERENT TEMPERATURE








BY

ALEX X. ZHAO


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


1990























To my parents and my beloved motherland












ACKNOWLEDGEMENTS


The successful completion of this dissertation would have been impossible if

it had not been for the assistance and help I received from numerous people. First

mention must be made for the academic guidance and financial support made

possible by my graduate adviser Dr. Ulrich H. Kurzweg. He, with his strong

background and extensive knowledge in fluid mechanics, heat transfer and applied

mathematics, gave me clear guidance in maintaining the correct research direction.

What I learned from him are not only the relevant knowledge in those fields but also

the skills needed to conduct scientific research. Two other persons who helped me

significantly in the dissertation were Dr. Wei Shyy and Dr. Ranga Narayanan.

Several fruitful discussions with Dr. Shyy helped me to construct the numerical

procedure used successfully and more efficiently for the time periodic problems in

this dissertation. The experience of working with Dr. Narayanan put me in contact

with many interesting topics in fluid mechanics, heat transfer and numerical analysis,

which made my understanding of the enhancement of heat transfer problem

considered here much more clear. I am thankful to the other and former members

of my supervisory committee, Dr. E Rune. Lindgren, Dr. Joseph L. Hammack, Dr.

Chia-Shun Yih, Dr. Edward K. Walsh and Dr. Tom I-P. Shih, for their time and

effort.







Thanks are also expressed here to the Department of Aerospace Engineering,

Mechanics and Engineering Science and the Department of Chemical Engineering

of the University of Florida for their financial assistance in allowing me to pursue my

educational goals within this lovely country. I am also grateful to the early support

of this work by National Science Foundation through the grant CBT-8611254 and the

allocation of 95 Service Units for this study on the CRAY/YMP supercomputer in

the Pittsburgh Supercomputing Center.

I would also like to thank my wife Wei Zhu for being here with me in a

"foreign" country and working hard to help support the family. Without her help, I

would not have been able to successfully complete my education. I also wish to

express gratitude to my wonderful office mate Bill Harter. Through many discussions

with him, I learned a lot about PC computers, about the English language, about

American customs and etc. Finally I wish to thank myself for putting in four hard

working years on this dissertation effort.












TABLE OF CONTENTS


page


ACKNOWLEDGEMENTS

LIST OF FIGURES

LIST OF TABLES

NOMENCLATURE

ABSTRACT


CHAPTERS


I INTRODUCTION ...................................... 1
1.1 Background ...................................... 1
1.2 Literature Review .................................... 5
1.3 Analytical Solutions Within the Connecting Channel ........ 9

II MATHEMATICAL MODEL .............................. 24
2.1 Geometric Configurations ........... .................. 24
2.2 Governing Equations and Boundary Conditions ............. 28

III NUMERICAL TECHNIQUE EMPLOYED ..................... 32
3.1 Grid and Transformation of Governing Equations ......... 33
3.2 A Brief Derivation of the Associated Finite Difference
Equations and Boundary Conditions ..................... 38
3.3 Calculation Procedure .............................. 43
3.4 A Vectorized Line Group Method for Solving
the System of Algebraic Equations ...................... 47

IV RESULTS AND DISCUSSION ............................. 55
4.1 Validation of the Results ............................. 56
4.2 Velocity Field ...................................... 70
4.3 Temperature and Heat Transfer Coefficient ............... 88


...................................... iii

..................... ................ vii

................... ................... xii


..................................... xvii
...................................... xvii







V CONCLUDING REMARKS


.............................. 106


APPENDIX THE VELOCITY-PRESSURE CORRECTION .............. 110

REFERENCES ................................................ 116











LIST OF FIGURES


Figure M

Fig. 1-1 Heat transfer in a parallel plate channel: (a) pure conduction;
(b) enhancement by oscillating flows. .... 3

Fig. 1-2 A schematic of the Thermal Pump (U.S. Patent No. 4590993).
5

Fig. 1-3 Geometric configuration for the analytic solution of flow
oscillations and the associated heat transfer in parallel plate
channels. .... 10

Fig. 1-4 Analytic solution of the dimensionless velocity of oscillating
flow in a parallel plate channel as a function of dimensionless
transverse position at different oscillation phases ot: (a) a=1;
(b) a=10. .... 12

Fig. 1-5 Analytic solution of the dimensionless temperature distribution
in oscillating flow in a parallel plate channel and its walls as a
function of dimensionless transverse position at different
oscillation phases wt: (a) a=1; (b) a=10. .... 14

Fig. 1-6 Analytic solution of the dimensionless time averaged effective
thermal diffusivity enhanced by oscillating flow in a parallel
plate channel as a function of Womersley number a at
different Prandlt numbers: (a) e=2 (conducting walls with
a = =1); (b) e=l (insulating walls). .... 16

Fig. 1-7 Analytic solution of dimensionless time averaged effective
thermal diffusivity enhanced by oscillating flow in a parallel
plate channel at e=2 (conducting walls) as a function of
Womersley number a for several different liquid-solid
combinations. .... 19






Fig. 1-8 Analytic solution of dimensionless time averaged effective
thermal diffusivity enhanced by oscillating flow as a function of
cladding thickness for liquid lithium in a parallel plate channel
with stainless steel walls. .... 19

Fig. 1-9 Analytic solution of heat transferred by an oscillating flow of
water in a parallel plate channel with glass walls for different
a: (a) time averaged value as a function of dimensionless
transverse position; (b) cross section integrated value as a
function of oscillation phase. .... 20

Fig. 2-1 Configurations of the thermal pump used in the numerical
investigation: (a) conduction model; (b) cross flow model. .... 23

Fig. 3-1 The correspondence of blocks between the physical domain
and the computational domain: (a) grid in the physical domain
at ot=0; (b) time independent and uniform grid in the
computational domain; (c) grid in the physical domain at
wt=180o. .... 32

Fig. 3-2 The movable and non-uniform grid in the physical domain: (a)
at=0; (b) at=90; (c) wt=180. .... 33

Fig. 3-3 The staggered grid in the computational domain. .... 36

Fig. 3-4 Typical cells for the velocity component U: (a) an interior cell;
(b) a boundary cell with its south wall at the domain boundary. .... 38

Fig. 3-5 Flowchart for determining the time periodic velocity
components U, V and pressure P. .... 42

Fig. 3-6 Flowchart for solving the time periodic temperature
distribution based on the calculated velocity. .... 45

Fig. 3-7 Flowchart for the Line-By-Line iteration method for solving the
associated system of algebraic equations. Dashed lines
represent the innermost DO loops. ... 48

Fig. 3-8 Flowchart for the Vectorized Line Group iteration method for
solving the associated system of algebraic equations. Dashed
lines represent the innermost DO loops. .... 50







Fig. 4-1 Stream line patterns for Case No.6 calculated with different
number of grid points: (a) grid 91x71; (b) grid 121x101.
Increment between neighboring lines is A r=20 in non-
dimensional units. .... 54

Fig. 4-2 Non-dimensional velocity component U for Case No.5 at the
middle of the connecting channel as a function of oscillation
phase wt: solid line numerical result; dashed line analytic
solution. .... 57

Fig. 4-3 Non-dimensional temperature T for Case No.5 at the middle of
the connecting channel as a function of oscillation phase ft:
solid line numerical result; dashed line analytic solution. .... 57

Fig. 4-4 Non-dimensional velocity component U for Case No.5 at the
middle of the connecting channel as a function of non-
dimensional position y at several oscillation phases: (a)
analytic solution; (b) numerical result. .... 60

Fig. 4-5 Non-dimensional temperature T for Case No.5 at the middle of
the connecting channel as a function of non-dimensional
position y at several oscillation phases: (a) analytic solution;
(b) numerical result. .... 61

Fig. 4-6 Non-dimensional velocity component U for Case No.6 at the
middle of the connecting channel as a function of non-
dimensional position y at several oscillation phases: (a)
analytic solution; (b) numerical result. .... 62

Fig. 4-7 Non-dimensional x direction pressure gradient for Case No.5 at
the middle of connecting channel as a function of oscillation
phase ft: solid line numerical result; dashed line analytic
solution. .... 64

Fig. 4-8 Stream line pattern for Case No.3 at oscillation phases of
interval A(ot) =150 in the lower half cycle 0tot<1800. .... 66

Fig. 4-9 Temporal approach to the final time periodic state of the non-
dimensional velocity component U(x,O, ot) for Case No.3 at
three different locations: (a) x=-152; (b) x=-175; (c) x=-190. .... 69

Fig. 4-10 Stream line pattern for Case No.4 at oscillation intervals of
A(ot) =15 in the lower half cycle Oscot <180. .... 70







Fig. 4-11 Stream line pattern for Case No.6 at oscillation intervals of
A(ct) =150 in the lower half cycle 00 t<1800. .... 71

Fig. 4-12 Non-dimensional pressure distribution along the x axis in Case
No.6 at oscillation intervals of A(wt) =450 over the whole cycle
0" cot<3600. .... 73

Fig. 4-13 Stream line pattern for Case No.6 at the oscillation intervals of
A(at) =450 over the whole period 00Oat<3600. .... 74

Fig. 4-14 Stream line pattern for Case No.1 at the oscillation intervals of
A(ct) =150 in the lower half cycle 00-t<1800. .... 76

Fig. 4-15 Stream line pattern for Case No.8 in the right reservoir at
oscillation intervals of A(at)=900 over the whole period
0at <360. .... 78

Fig. 4-16 Stream line pattern for Case No.9 at oscillation intervals of
A (ot) =15 over the lower half period 00
Fig. 4-17 Stream line pattern for Case No.10 at oscillation intervals of
A (t) = 150 in the lower half cycle 00s t <1800. .... 80

Fig. 4-18 Temperature contours for Case No.3 at oscillation intervals of
A (at) =150 in the lower half cycle 00 t <180. .... 83

Fig. 4-19 Non-dimensional temperature profile for Case No.3 along the
vertical line x=-153 in the left reservoir: (a) at several
oscillation phases; (b) time averaged. .... 85

Fig. 4-20 Temperature contours for Case No.4 at oscillation intervals of
A(act) =150 in the lower half cycle 00ot<1800. .... 86

Fig. 4-21 Non-dimensional temperature profile for Case No.4 along the
vertical line x=-153 in the left reservoir: (a) at several
oscillation phases; (b) time averaged. .... 87

Fig. 4-22 Temperature contours for Case No.6 at oscillation intervals of
A(ot) =150 in the lower half cycle 0o
Fig. 4-23 Non-dimensional temperature profile for Case No.6 along the
vertical line x=-153 in the left reservoir: (a) at several
oscillation phases; (b) time averaged. .... 89







Fig. 4-24 Temperature contours for Case No.9 at oscillation intervals of
A(at) =150 in the lower half cycle 0t <1800. .... 91

Fig. 4-25 Non-dimensional temperature profile for Case No.9 along the
vertical line x=-153 in the left reservoir: (a) at several
oscillation phases; (b) time averaged. .... 92

Fig. 4-26 Temperature contours for Case No.10 at oscillation intervals of
A (t) = 150 in the lower half cycle 00ot <1800. .... 93

Fig. 4-27 Non-dimensional temperature profile for Case No.10 along the
vertical line x=-153 in the left reservoir: (a) at several
oscillation phases; (b) time averaged. .... 94













LIST OF TABLES


The list of ten cases studied in the numerical investigation.

Heat transfer coefficients for several cases in the study.


... 51

... 95


Table 4-1

Table 4-2












NOMENCLATURE


a the half width of the connecting channel

A coefficients in the Finite Difference Equations

Ao the cross section area of the connecting channel

b the distance from the channel center line to the wall center line in the

analytic solution or the half width of the end reservoir in the numerical study

c time averaged length of the end reservoir

Cp specific heat of working fluid

D molecular mass diffusivity of working fluid

h, heat transfer coefficient from the analytic solution

h, heat transfer coefficient from the numerical results

L channel length

P pressure

Pr Prandtl number

q local instantaneous axial heat flux in the channel

Q total time averaged heat rate in the channel

r, the ratio of tidal displacement to channel length

Re Reynolds number based on the maximum of piston velocity and reservoir half

width







RHS

S

Sc

T

T,

Tc

T,



u

U

Up

v

V

Vc

x

6x

y

6y


terms on the Right Hand Side of the associated system of algebraic equations

source term in the Finite Difference Equations

Schmidt number

temperature

high temperature specified at hot end for supplying heat

low temperature specified at cold end for removing heat

the approximate time averaged temperature in right reservoir from the

numerical results

velocity component in axial (x) direction at the walls of U cells

velocity component in axial (x) direction

maximum piston velocity

velocity component in transverse (y) direction at the walls of U cells

velocity component in transverse (y) direction

cross flow velocity

axial coordinate in the physical domain

x direction size of a grid cell in the physical domain

transverse coordinate in the physical domain

y direction size of a grid cell in the physical domain







a Womersley number

p stretching parameter for generating non-uniform grid

y time averaged axial temperature gradient (constant)

6i variable time step size in calculation

6 B basic time step size in calculation

Ax tidal displacement

e the ratio of total channel half width (including the thickness of solid wall) to

the channel half width filled by fluid in the analytic solutions or general

convergence criterion in the numerical calculations

qe convergence criterion for the Velocity-Pressure correction

ep convergence criterion for the final time periodic states

r1 transverse coordinate in the computational domain

0 magnitude of the sinusoidal pressure gradient applied in the analytic solutions

K molecular thermal diffusivity of working fluid

Keff effective thermal diffusivity by oscillating flow

p the ratio of fluid-to-solid (wall) thermal conductivities

v kinematic viscosity of working fluid

t axial coordinate in the computational domain

p density of working fluid

a the ratio of fluid-to-solid (wall) thermal diffusivities in analytic solutions or

general convergence indicator in numerical calculations

Ca convergence indicator for Velocity-Pressure correction







a, convergence indicator for the final time periodic state

r time in the computational domain


ir stream function of the oscillating flow field

o angular frequency of the oscillation






Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment
of the Requirements for the Degree of Doctor of Philosophy


FLUID OSCILLATIONS AND HEAT TRANSFER
IN A PARALLEL PLATE CHANNEL BOUNDED BY END RESERVOIRS
MAINTAINED AT DIFFERENT TEMPERATURE


by

Alex X. Zhao

December 1990

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


Existing analytic and experimental studies have shown that sinusoidal

oscillations of a viscous fluid within open ended conduits connected to reservoirs can

enhance the thermal diffusivity between the hot fluid and cold fluid in the opposite

end reservoirs by some four orders of magnitude in excess of that possible in the

absence of oscillation. The heat transfer coefficients achieved in this process can be

very high and can readily exceed those possible via heat pipes, yet involve no net

convective mass exchange.

A numerical investigation of laminar oscillating flows of incompressible

viscous fluid and the associated enhanced heat transfer in a 2D parallel plate channel

bounded by rectangular end reservoirs which have sinusoidally oscillating piston

boundaries and are maintained at different temperatures forms the topic of the

present dissertation.


xvii






A modification of the widely used SIMPLE algorithm, termed SIMPLE-TP

(for Time Periodic), has been developed to take time periodicity and boundary

movements into account. The method has been successfully applied to the present

time periodic flow and heat transfer problem.

Two different models of heat supply and removal, namely, conduction through

the reservoir walls and convection by cross flow, were considered in this study.

Several case studies involving different combinations of Womersley number, tidal

displacement, length of reservoirs and cross flow velocities were carried out. The

numerical results obtained were compared with analytic solutions where such

comparisons were possible.

It was found that the velocity field in each reservoir is characterized by a high

velocity jet emanating from the channel end during a portion of the oscillation cycle

and that one or more counterrotating vortex(es) or vortex pair(s) exist during the

whole oscillation cycle as long as the Reynolds number, based on the maximum

piston velocity and reservoir width, is large enough. The tidal displacement and cross

flow velocity have a strong effect on the flow patterns in the reservoir and the

reservoir temperature field is mainly determined by the temperature boundary

conditions and less so by the velocity field. High oscillation frequencies coupled with

a large tidal displacement and a relatively large cross flow velocity are found to yield

large enhancements of heat transfer. Heat transfer coefficient as high as 1.8x106

W/m2K can readily be achieved when using water as the working fluid.


xviii












CHAPTER I

INTRODUCTION



1.1 Background

For a sinusoidal oscillatory fluid flow over a flat solid surface, the influence

of the no-slip condition on the solid surface penetrates an approximate distance of

2v_/&c into the fluid, where v is the kinematic viscosity of the fluid and a is the

angular frequency of the flow oscillation. For an oscillatory flow inside a flat plate

channel, this is still true as long as the channel half width, a, is greater than the

penetration distance. Under this condition, the fluid flow can be considered to

consist of two parts, the boundary layers and the core region, which are shown in Fig.

1-1(b).

If a property of the fluid, namely, concentration or temperature, is maintained

at different values at the two ends of the channel, a very strong transport from the

high side to low side will occur. Some theoretical and experimental studies have

shown that the effective coefficient for the transport can be much larger than the

values which characterize the transport without oscillation. The reason for the

enhancement of the transport is that both the area over which the exchange takes

place and the property gradient across the area are increased significantly over the

values existing in the absence of oscillations.

1






2
Consider the case where the temperature of the liquid is different at the two

ends of the channel. Here the heat transfer from the hot end to the cool end can

be enhanced by a flow oscillation in the axial direction. The time averaged heat

transferred axially, without a accompanying mass exchange, can be expressed as

= -AopCpKg (1,1)

where Kef, p and Cp stand for the effective thermal diffusivity, fluid density and

specific heat, respectively, Ao is the channel cross section, and y is the time averaged

axial temperature gradient which is maintained as a constant.

For pure molecular conduction, Kfe simply is the molecular thermal diffusivity

, Ao is equal to 2a when the channel depth is equal to one, as is shown in Fig. 1-

l(a), and the amount of heat transferred per unit depth will be

Qc = -2apCpKY (1,2)

With oscillations, it is necessary to examine the heat transfer process in some

detail since it becomes a time dependent problem. During the axial fluid oscillation

towards the colder side, hot liquid is brought into the core region and moved in

contact with the cooled boundary layers, which in the previous half stroke are in

contact with cooler liquid in the core region. The colder boundary layers will absorb

heat from the hot liquid core and heat up. During the next half stroke the oscillation

is towards the hot side, and cold liquid is brought into contact with the heated

boundary layers. The boundary layers then release heat to the cool core region and

thus become cool again. These two steps of heat transfer in oscillating flows at








(a)











(b)


Heat transfer in a parallel plate channel: (a) pure conduction; (b)
enhancement by oscillating flows.


I
2a Th Tc







2v/w





T
/ c

boundary layers



2a core region Tc

Th


Fig. 1-1







4

Prandtl number 1, for which the thicknesses of momentum and thermal boundary

layers are equal, are shown in Fig. 1-1(b) with the dotted areas representing the cold

fluid. The exchange process is repeated during subsequent oscillation cycles. This

leads to a time periodic problem. Because the heat exchange happens between the

boundary layers and core region, the effective heat transfer surface become much

larger than Ao and can be expressed in terms of a quantity Ax named the tidal

displacement (equal to twice the cross section averaged axial amplitude of the

sinusoidal oscillations). The temperature gradient in the normal of interface between

the boundary layers and core region is also much larger than y and is a function of

time and position. By averaging over the whole oscillation period and channel cross

section, the heat transferred per unit depth from the hot end to the cool end now can

be written as

S= -2apCp, fy (1,3)

where now the effective thermal diffusivity Krg is a function of a, v, <, Ax and some

other factors. Krf can be much larger than K provided a good choice for those

factors is made. Note that there will be no net convection between the two ends

when Ax is kept at a fraction of the channel length, and the process thus can be

considered to be free of mass exchange between the channel ends under proper

conditions (explained in Section 1.3).

If the concentration is maintained at different values at the ends of the

channel, an enhanced axial mass diffusion will occur for appropriate conditions, thus

producing an effective diffusivity much larger than the molecular diffusivity. In







5
recent years, these phenomena have received considerable research attention from

engineers and scientists because of their potential practical applications such as for

nuclear reactor cooling and for gas separation, etc.

A typical device for enhancing heat transfer by oscillating flows, termed the

Thermal Pump, is shown in Fig. 1-2. This represents essentially US Patent No.

4590993 held by U. Kurzweg1 and the University of Florida. This dissertation will

be concerned with a numerical investigation of the laminar flows of incompressible

fluids and the associated heat transfer in a 2D version of the thermal pump when

sinusoidally moving piston boundaries are located in the end reservoirs. The

connecting channel will consist of a single flat plate channel formed by two parallel

plates separated by a fixed distance of 2a.



1.2 Literature Review

The phenomenon of transport enhancement of contaminants was initially

studied by Taylor2(1953) and Aris3(1956) in their investigations of axial dispersion

in steady flows within capillary tubes. Their results indicated that a significant

increase in axial contaminant dispersion will occur in steady laminar flows. Similar

effects were found to exist for oscillating flows by Bowden4 in 1965. Harris and

Goren5(1967) studied, both analytically and experimentally, the mass transfer

through a long tube connecting two reservoirs maintained at different but constant

concentrations by oscillating the flow in the tube and found that the increase of mass

transfer rates is a function of the Womersley number a =a o/v the tidal























Cooler Fluid Out


Heat Transport


Channels


Hot Fluid In Hotter Fluid Out


A schematic of the Thermal Pump (U.S. Patent No. 4590993)1.


Left End
Chamber


Right End
Chamber


Oscillator


Cold Fluid In


Fig. 1-2






7

displacement Ax and the Schmidt number v/D. Rice and Eagleton6(1970) conducted

a similar study to that of Harris5 and found that the mass transfer increase is pro-

portional to the square of tidal displacement Ax. In 1975, Chatwin7 showed theo-

retically that the effective diffusion coefficient in oscillating flows is a harmonic

function of time with a period equal to one half that of the imposed time dependent

pressure gradient. Further theoretical studies by Watsons(1983) and experiments by

Joshi et al.9(1983) gave similar results. Kurzweg and Jaeger"0(1986) first showed

the existence of a tuning effect in enhanced gas dispersion in oscillating flows. All

these studies prove that a contaminant will spread axially in oscillating laminar pipe

flow at rates as much as five orders of magnitude higher than those achievable by

pure molecular diffusion in the absence of oscillations.

The same enhancement in oscillating flow for heat transfer was discovered by

Kurzweg in 1983 in view of the mathematical similarity between the governing

equations for heat conduction and mass diffusion. The equations for contaminant

diffusion with superimposed oscillating flow are mathematically nearly identical to

the equation of heat conduction, with the only major differences being the boundary

conditions and the use of the Prandtl number Pr= v/I instead of the Schmidt number

Sc=D/c. He solved the governing equations for enhanced axial heat transfer in

oscillating flows in a tube of infinite length for low value of CftPr and found the

effective conductivity can readily reach values three orders of magnitude greater than

normal heat conductivity. This information was published in 198511. Kurzweg and

Zhao'2(1984) conducted some experiments to measure the effective axial heat con-







8
duction rates through a capillary bundle connecting two fluid reservoirs maintained

at different temperatures and found very good agreement with theoretical

predictions. Two other papers by Kurzweg(1985)13(1986)14 explained the tuning

effect in the enhancement of heat transfer in detail. In May of 1986, a U.S. patent

was granted to the University of Florida and Professor U. Kurzweg1 on this process.

Zhang15 performed some numerical studies on tubes of finite length with some

portions of their walls or ends maintained at different temperatures and concluded

that the highest rates of heat transfer predicted earlier can only be reached when the

amount of heat supplied and removed from the tube ends is large enough.

Otherwise a heat "bottle neck" will occur and the temperature gradient in the axial

direction of the tube will be reduced considerably. Kaviany's16 recent study on heat

transfer by oscillating flows in tubes bounded by two end reservoirs maintained at

different temperatures verified, both theoretically and experimentally, the results

obtained by Kurzweg and Zhao12.

It should be noted that most of theoretical studies to this point, both analytic

and numerical, have considered only the tube or connecting channel without any

attention to the flow behavior and heat transfer process occurring in the two end

reservoirs which generally have much larger cross section than the connecting

conduit. Kaviany16 has made some calculations on this flow using a very simplified

model of the end chambers. Durst et al.17(1989) did some numerical studies on the

transient laminar flow over a sharp pipe expansion driven by a piston in the wide end

moving as a rectangular function of time starting from an initial rest condition and







9
compared the numerical result with their experimental data18. Neither of them

considered the temperature distribution and their results for velocity can not be

directly used in the analysis of the thermal pump process because of the

simplification made in former and the different velocity conditions used in the latter.

The motivation for this dissertation research was to obtain a clearer under-

standing of the flow and temperature patterns inside of the end reservoirs of a

thermal pump in order to lend further guidance for the practical construction and the

usage of thermal pumps. Also, associated with this study, a proper numerical

procedure based on the SIMPLE algorithm was developed for the simulation of heat

transfer in time periodic flows of viscous incompressible fluids with moving

boundaries.





1.3 Analytic Solutions Within the Connecting Channel

The numerical investigations presented in the chapters II through V can be

considered as an extension of the analytic work done by Kurzweg13, which gave

solutions of the momentum and energy equations in a parallel plate channel of

infinite length. The numerical results in the connecting channel for the more

complicated geometry with end reservoirs should approximate those analytic solutions

when the connecting channel length becomes much longer than the tidal displace-

ment. Thus, it is essential to first briefly evaluate the existing analytic solutions in







10

order to later make a comparison with the numerical results to be obtained wherever

such a comparison is possible.

The configuration dealt with in Ref. 13 is shown in Fig. 1-3. It consists of a

set of channels filled with a viscous fluid and connected to chambers at x = The

individual channels of width 2a each are confined between thermally conducting walls
aP
of width 2(b-a) each and a time periodic pressure gradient a = Oe'"' and a
aT
constant axial temperature gradient = Y are superimposed. The velocity of the

oscillating flow in the central channel generated by the pressure gradient is given by


1 cosh(-Ta y/a)
cosh(ri a)
U(y,t) = REAL Ii Ax cosh(ITa) e (1,4)

2 1 tanh(v-a)



where y is the coordinate perpendicular to the channel axis as shown in Fig. 1-3, W

is the angular frequency of the time periodic pressure gradient, a = a / v the

Womersley number, with v representing fluid kinematic viscosity, and Ax the tidal

displacement for the maximum cross stream averaged axial excursion of fluid

elements during one oscillation cycle. The vertical bars indicate the absolute value

of the quantity shown.

The relation between the magnitude of the axial pressure gradient and the

tidal displacement can be expressed as



















K>


I I


f
Periodic
velocity profile


b


Thermally conducting
fluid


T.
x


Geometric configuration for the analytic solution of flow oscillations
and the associated heat transfer in parallel plate channels.


2(b-a)

LA
\10\1~x


ic
ie
t


Fig. 1-3


r





K~



a~a


solid thermal
conductor



Period
pressure
gradier



--


i










S: P 02,&x

2 1 tanh(-a) (1,5)



In practice, the channel half width a, the tidal displacement Ax, and the

oscillation frequency o as well as the properties of the working fluid are easy to

control. Hence, either they or a combination of some of them, in particular the

Womersley number, are used as the control parameter for this heat transfer problem.

A numerical evaluation of the non-dimensional velocity, normalized by CoAx,

as a function of non-dimensional transverse position in the y direction at different

oscillation phases for a=1 and a=10 is shown in Fig. 1-4(a) and (b), respectively.

It is clear that the velocity profile varies considerably with change in Womersley

number a, from an essentially parabolic shape at very low a to one having a nearly

constant velocity core connected to thin boundary layers of thickness 6 = ~/2 v /

at the walls for large a.

By assuming


T(x,y,t) = y(x + REAL[ag(y)ei"t] (1,6)


,the cross-section dependent portion of the temperature perturbation gly) in the

fluid and g,(y) in solid walls have the forms13










(a)


0.5


0.0


-0.5


-1.0 1 1 I 1 '
-0.8 -0.6 -0.4 -0.2 0.0
U/cAx


(b)


0.5


N 0.0 -


-0.5


-1.0 L
-0.8



Fig. 1-4 An
par
pos


0.2 0.4 0.6 0.8


-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8
U/WAx


alytic solution of the dimensionless velocity of oscillating flow in a
rallel plate channel as a function of dimensionless transverse
sition at different oscillation phases ct: (a) a=1; (b) a=10.










1 a2pr 1 cosh(-a y/a) + H cosh(Vi-Pa y/a)
gf (y) _-[ cosh(V a) cosh(ia) (1,7)

2aa2(Pr-1) 1 tanh( iJa)
hx 4i-a



g Y) 1 H cosh[V7 -Pr a(e-y/a)]
2a2(Pr-1) 1 tanh(Ta) cosh[VTaPr a (e-1)] (1,8)
2a a2Pr 1) tanh(fa) v '0/
Ax Fa


with


H = 11"Pr tanh(F- a) + -Ytanh[VaPr a(e -1)] (1,9)
p tanh(Ti Pr a) + fV- tanh[V/i a Pr a(e -1)]

where i =kf/k, is the ratio of fluid-to-solid (wall) thermal conductivities, o= r/Kr

is the ratio of thermal diffusivities, and e=b/a.

The non-dimensional temperature profiles T/yAx versusy/a for p = a=1 and

e=2 but different Womersley number a are shown in Fig. 1-5(a) and (b). One

observes an appreciable penetration of the temperature variation into the solid wall

at a=l, but almost no penetration occurs for a=10. This is to be expected in view

of the fact that the typical penetration distance into a solid of an oscillating

temperature field is proportional to the inverse square root of the oscillation

frequency (i.e. skin effect). Note that the frequency ratio at a=1 to that at a=10 is

one hundred for the same channel width and the same fluid.











(a) 2.0
a= 1
wall = 1 -
a=1
Pr=1
1.0 channel 10 -
150 60
a I / 300 / t=0 120 90-
S0.0 -
> 210 180 30
270

-1.0 240 330
wall

-2.0
-0.4 -0.2 0.0 0.2 0.4
T/yAx



(b) 2.0 i i I '

-wall I= 1 -
a=1
Pr= I -
-channel
30 60 90 120
0- ut=O 150
S- 330A 180
30 270 240 210

-1.0
wall

-2.0
-0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8
T/yAx

Fig. 1-5 Analytic solution of the dimensionless temperature distribution in
oscillating flow in a parallel plate channel and its walls as a function
of dimensionless transverse position at different oscillation phases (t:
(a) a=I; (b) a=10.







16
Neglecting axial conduction in both the wall and the fluid, the local and

instantaneous axial heat flux can be expressed as

4(y,t) = p C U(y,t) T(0,y,t) (1,10)

The time averaged total heat flow rate, then, will be obtained by integrating

this quantity over the channel width and the oscillation period, i.e.

C 2 n +a
_= p f U(y,t) T(0,y,t)d dd(t) (1,11)
2 o -a

By means of Equations of (1,3), (1,4) and (1,7), and after some rather lengthy

but straight forward manipulations, the effective thermal diffusivity of the oscillation

flow is found to be



Pr[(1 -H)I + (1-H)1] + (I-m H) + (1-mH)
eff2 (1,12)
16 a2 (Pr2 1) 1 tanh(tfa)
d0x2 VXa


where

I = F/ca tanh(v-Ta) (1,13)

m = VPr a tanh(Vi7Pr a) (1,14)

with the bar indicating the complex conjugate of the function shown and H being

defined by equation (1,9).







17
Two graphs of the non-dimensional rff/a x2, as a function of Womersley

number for different Prandtl numbers, are presented in Fig. 1-6(a) and (b). The first

one, (a), for e=2 representing the case of wall thickness being the same as channel

width and the second one corresponding to insulating walls which is mathematically

equivalent to the case of conducting walls with zero thickness. The curves shown are

seen to have a single maximum point which shifts to lower values of a as Pr is

increased. This is termed the tuning effect because the peak indicates the value of

channel width at which the heat transfer is optimized. These tuning points corre-

spond to the condition where the thermal diffusion time from channel axis to the

walls is approximately equal to half of the oscillation period. The observed drop off

in the peak value of r!C/alx2 for very small Prandtl number in the insulating case

is attributed to the inability to store much heat in the boundary layers which become

quite thin for fluids with small kinematic viscosity v such as in liquid metals.

Note that the second set of curves, Fig. 1-6(b), also present information for

the enhanced axial mass diffusion by oscillating flows since insulating walls give the

same boundary conditions as those for diffusional mass transfer and both of them

come from the same type of equations with the only difference being the use of

Prandtl number for heat transfer but the Schmidt number for mass transfer. One

finds that the peak value of rgff/Ax2 for Pr=l is about 102 in the vicinity of a =3

while the corresponding value for the effective mass diffusivity for Schmidt number

1000 at the same Womersley number falls down to 5x10-s by considering the curve

for Pr=103 as a curve for mass transfer for a Sc=10 fluid. Since most liquids have











(a)
e=b/a=2 (conducting walls)
10 -

10 -a Pr=10 /



< 10 -4
10

10 IT 3
10


10-3 10 -a 10 -' 1 10 10 z 10 S
a



(b) 1
(b ) i n1 -- i iiiii|-l i ,i,,. ii 1'" 1 ii i i l "
Se=b/a=l (insulating walls)
10 -1

10 -

< 10- Pr= 103
3 10
100-

10 -4

10

10 10 10 1 10 10 2 10 3
O(


Fig. 1-6 Analytic solution of the dimensionless time averaged effective thermal
diffusivity enhanced by oscillating flow in a parallel plate channel as
a function of Womersley number a at different Prandit numbers: (a)
e=2 (conducting walls with a=p ==1); (b) e=1 (insulating walls).







19
a large value for the ratio of Schmidt number to Prandtl number (about 1000 to 1

for water, for example), a really large increase in heat transfer with a negligible

enhancement of diffusional mass transfer can be achieved, and thus, the process of

heat transfer by thermal pumping in liquids can be considered as essentially free of

mass transfer provided the connecting channel length is longer than the tidal

displacement. Note that this same observations does not follow for gases where

Schmidt number and Prandtl number are approximately the same.

The tuning curves of effective thermal diffusivity for several different liquid-

solid combinations in configurations of equal wall and channel thickness (e=2) are

shown in Fig. 1-7. There is no major difference in the maximum Ki, among these

different material combinations and all curves give a maximum ratio of Kff to oAx2

of about 0.05. It is clear that the effective maximum thermal diffusivity at the peak

of the tuning curves is directly proportional to the frequency of the fluid oscillation

and to the square of the tidal displacement. This means that in the optimization pro-

cedure for a thermal pump, o and Ax should be made as large as possible. The

channel width a, then, should be determined from the valve of a corresponding to

the maximum effective thermal diffusivity and the kinematic viscosity v of the

working fluid. For example, at wo=2xr rad/s in the water-glass combination, the

optimum charmel half width will be a =2mm for optimum heat transfer.

The effect of wall thickness on the heat transfer process in the thermal pump

is shown in Fig. 1-8, where a plot of ./foAx2e versus e for several values of a2Pr

for the liquid-solid combination of liquid lithium and stainless steel are presented.












10 -a
10-1 C

10-2
Lx
3

10 -'

10-

10 -6
10 -3


Fig. 1-7





0.025


0.020


0.015


3 0.010


0.005


0.000



Fig. 1-8


10 -a 10 -1 1 10 10 2


10


a
Analytic solution of dimensionless time averaged effective thermal
diffusivity enhanced by oscillating flow in a parallel plate channel at
e=2 (conducting walls) as a function of Womersley number a for
several different liquid-solid combinations.


2 3 4 5


Analytic solution of dimensionless time averaged effective thermal
diffusivity enhanced by oscillating flow as a function of cladding
thickness for liquid lithium in a parallel plate channel with stainless
steel walls.






21

These curves show that an e=2 geometry, corresponding to equal wall thickness and

channel width, is about optimum when an appreciable fraction of the heat can be

stored in the bounding walls. The reason for normalizing Kef by the extra factor e

is to take into account the fact that the walls take away a fraction 1/e of the cross

sectional area of the thermal pump available for the fluid oscillations.

Integration of the local insantaneous heat flux over the oscillation period

yields the time averaged axial heat flux distribution over the cross section, while an

integration over the channel width yields the instantaneous axial heat flow through

the whole channel as a function of time. Fig. 1-9(a) and (b) show the behavior of

the two partially integrated, non-dimensional products of U(y,t) and T(0,y,t) for the

case of a water-glass combination with e=2. From (a) one observes that the largest

value of UT/yoAx2 occurs for a=1 which is near the tuning point for this particular

condition. Note the much smaller values at a=0.1 and a=10 which occur at de-

tuned conditions. The flux transport becomes confined to the vicinity of the walls as

the Womersley number becomes large. Fig. 1-9(b) shows that the thermal flux is

time periodic at a frequency twice that of the fluid oscillations and that only near

tuning conditions is the positive flux not essentially canceled by the negative heat

flux.

In the following chapter, Chapter II, the geometry for which this numerical

study was performed will be explained in detail. The assumptions and the mathemat-

ical models used to present two different physical situations will also be shown there.

Chapter II will focus on the numerical method employed including a brief derivation











(a) 1.0 -
0.8
0.6
0.4
0.2
"' 0.0
-0.2
-0.4
-0.6
-0.8
-1.0
-0.02


(b)


0.5


-0.5


-1.0


Fig. 1-9


0.00 0.02 0.04 0.06
UT/7yAx2


60 120 180 240 300


0.08


360


Ct


Analytic solution of heat transferred by an oscillating flow of water in
a parallel plate channel with glass walls for different a: (a) time
averaged value as a function of dimensionless transverse position; (b)
cross section integrated value as a function of oscillation phase.







23
of the associated Finite Difference Equations. The results obtained from the numeri-

cal investigation will be presented in Chapter IV, with some comparisons of the

results with analytic solutions discussed in the present section. Concluding remarks

will be found in Chapter V. Some materials useful for a clearer understanding of the

Velocity-Pressure correction in the SIMPLE algorithm are summarized in Appendix.












CHAPTER II

MATHEMATICAL MODEL



The mathematical model for the numerical study presented in this dissertation

is explained in this chapter. First, the geometric configurations for which the

numerical investigation were performed are described. Then, the governing

equations to be solved and boundary conditions assigned are given.



2.1 Geometric Configurations

Two different configurations based upon the thermal pump geometry shown

in Fig. 1-2 were chosen for the present numerical study. These are shown in Fig. 2-1.

Some simplifications were introduced for the purpose of computational convenience.

Note also that, to more clearly show details, the axial-transverse ratio of the

configurations shown is different from their true value used in the calculations as

discussed in Chapter IV. Both configurations essentially consist of a 2D parallel

plate channel of width 2a and length L connected to two 2D rectangular end reser-

voirs of width 2b and time averaged length c. There is a piston moving back and

forth sinusoidally in the axial direction in each of the reservoirs. Choosing a rectan-

gular shape for the reservoirs allows calculations in Cartesian coordinates, although

a coordinate transformation is still necessary because of the moving boundaries.

24









(a) x=-(c+-L). fcosW1


Th insulating Tc
I 1 I I


X= (c+-L+ ACOS(A
2-\r 2b-Z


b


----- tttttttt" tttttttt-
Th Tc



Fig. 2-1 Configurations of the thermal pump used in the numerical
investigation: (a) conduction model; (b) cross flow model.


(b) x=-(C+L)- 1,coswtt






26

The coordinate system used in the physical domain is also shown in Fig. 2-1.

The origin is set at the center point in both axial and transverse directions. The axis

of the connecting channel and reservoirs was chosen as the x axis and then they axis

is in transverse direction. Based on this coordinate system, the positions of the

moving boundaries (piston surfaces) can be expressed as


XL, = -(c + ) Acos(Ot) (2,1)
2 2b


x, = (c + ) cos(ot) (2,2)
2 2b

for the left and right pistons, respectively. Note that the piston movements are in

phase in order to satisfy the law of mass conservation for the assumed incompressible

working fluid entirely filling the channel and reservoirs.

The first configuration, shown in Fig. 2-1(a) and named the conduction model

because of the way that heat is supplied and removed from the end reservoirs, is a

closed system for mass transfer. All walls except the piston surfaces are stationary.

The channel side walls are thermally insulating as indicated by the oblique shading

lines, while the reservoir side walls including the piston surfaces are conducting and

maintained at different temperature at the opposite sides displayed by the blank (for

cold at T,) and black (for hot at T,) regions. The configuration is symmetrical about

the horizontal center line corresponding to the x axis. This symmetry allows the

calculation to be undertaken in only the upper half of the domain as long as suitable

boundary conditions are applied at the center line.






27
The second configuration is shown in Fig. 2-1(b). It uses fluid cross flow to

add and remove heat from the end reservoirs and is termed the cross flow model.

It is assumed that the cross flow enters and exits the reservoirs at a constant but

small velocity V, perpendicular to the reservoir side walls. All walls including the

piston surfaces are insulating in this model. The cross flow is maintained at different

temperatures on the opposite sides as they come into the system. The thermal

conduction at the cross flow exits (i.e. transpiration) is excluded in the numerical

study since it is negligible compared with the heat carried by the flow. This configu-

ration no longer produces symmetric flows about the center horizontal line and the

computation thus has to be done over the entire flow domain.

A very useful characteristic for both configurations is that all information

needed to describe the flow and heat transfer in one oscillation phase step can be

readily obtained from the solution at another phase step one half period earlier as

soon as the flow pattern and temperature distribution in whole physical domain have

reached a final time periodic state. Time periodic conditions existing in the right end

reservoir at a fixed phase point in the oscillation cycle will correspond to those found

in the left reservoir half an oscillation cycle later. Furthermore, if temperature is

normalized so that the given temperatures Th and T, in opposite reservoirs are the

negative of each other, the temperature distribution at wt=1800 are the same as that

in the reflected image of the geometry at wt= 0 except for a sign change. The same

thing happens for every pair of solutions at wt and ot+1800. It is because of this time

symmetry that the computations need be taken over only the range of O"stsl800,







28

which significantly reduced the CPU time and storage required in the numerical

investigation.



2.2 Governing Equations and Boundary Conditions

The general assumptions in this study are following.

a) The working fluid is a homogeneous, Newtonian fluid;

b) The working fluid is incompressible

-constant density;

c) The working fluid has constant material properties

-constant viscosity,

-constant thermal diffusivity,

-constant thermal conductivity;

d) All body forces are negligible;

e) The viscous heating term in the energy equation is negligible;

f) The flow is laminar.

As shown in Section 1.3, liquid is a suitable working fluid in the thermal pump

for heat transfer without mass exchange. The incompressible fluid assumption is

valid for most liquids. The main advantage of the constant material properties

assumption is the ability to separate velocity and pressure solution from the tempera-

ture solution. It was found that in order to reach the time periodic state, the number

of oscillation periods needed for temperature is much more than that for velocity,

while the CPU time used for solving velocity and pressure distribution in each period







29
is much longer than that for temperature. Calculating velocity and temperature

separately did save considerable CPU time.

The dependent variables to be solved for are the x component of velocity, U,

the y component of velocity, V, the pressure P and the temperature T. They and the

independent variables x, y and time t are nondimensionalized as follows.

x* x y t* =t (2,3)
-= .. y- =_ = (2,3)
a a


U* U V* VP' (2,4)
wa wa p o2a2


ST- (Th + T)/2
T* (2,5)
(T T,)/2

where the quantities with are dimensionless and

a -half width of the channel;

W -angular frequency of the piston movement;

p -fluid density;

Th -temperature at the hot end;

Tc -temperature at the cold end.

After nondimensionalization, the governing equations read


V + V = 0 (2,6)


U, + (UU)x + (VU)y = -Px + (1/a2) (U + U )


(2,7)










V, + (UV) + (WV) = -P + (1/a2) (V, + V) (2,8)



T, + (UT)~ + (VT), = (1/Pra2) (T,, + T) (2,9)


where the subscripts indicate partial derivatives and

Pr= v/K Prandlt number;

a= a / v Womersley number;

Note that a conservation form of the governing equations is used here. All

quantities in the equations are dimensionless and the non-dimensional mark has

been omitted for simplicity. This notation, in which dimensionless quantities are

termed without superscript *, will be used from this point through the end of the

dissertation except where the quantity mentioned needs to be expressed in dimen-

sional form.

In view of the geometry explained in Section 2.1 and the non-dimensiona-

lization shown above, the boundary conditions can be mathematically expressed as

follows:

Velocity:

1) U 0 at the stationary walls; (2,10)
= V 0

2) {U = -sin t
2b at the piston surfaces; (2,11)

V= 0










3) u=0

V
a


at the cross flow entrance;


au
ay


4)


along the symmetry line y = 0.


(2,12)





(2,13)


V= 0


Temperature:


1) aT
an


2) T =+1


3) T =-1



4) 0

ere in Eq. (2,12) is
where VF in Eq. (2,12) is


at the insulating walls and symmetry line y =0; (2,14)


at the hot walls and hot cross flow entrance; (2,15)


at the cold walls and cold cross flow entrance; (2,16)



at the cross flow exits. (2,17)


dimensional and in Eq. (2,14) indicates that the partial
iT


derivative normal to the walls.












CHAPTER III

NUMERICAL TECHNIQUE EMPLOYED



The governing equations (2,6)-(2,9) with the boundary conditions (2,10)-(2,17)

specified in the physical domain, which have been presented in Chapter II, have to

be solved numerically and thus require a suitable numerical method. The existing

Semi-Implicit Method for Pressure Linked Equations (SIMPLE) of Patankar19 and

of Patankar and Spalding20 is one such method available for dealing with heat

transfer problems in incompressible viscous flows, but it needs to be modified here

by taking into account time periodicity and boundary movement, which are present

in the fluid oscillation and associated heat transfer problems investigated in this

dissertation. Such a consideration leads to the extended SIMPLE algorithm, termed

SIMPLE-TP (for Time Periodic), developed in this investigation and to be described

in some detail in this chapter. First, the coordinate systems in both the physical

domain and computational domain and the transformation of the governing equations

from the former to latter are described. Second, a brief derivation of the Finite

Difference Equations and the treatment of boundary conditions are presented. In

Section 3.3, an overview of the numerical iteration procedure SIMPLE-TP is given.

Finally, a Vectorized Line Group algorithm for solving the associated system of

algebraic equations also developed in this investigation is presented.







33
3.1 Grid and Transformation of Governing Equations

It is most convenient to transfer the time periodic physical domain to a time

independent computational domain. In other words, the grid system in the physical

domain which is movable and non-uniform needs to be mapped to a grid system in

the computational domain which is fixed and has even spacing in order to allow the

numerical calculations to be carry out properly. For the physical domain geometry

shown in Fig. 2-1, it is found convenient to break both of the domains into three

subregions when making the mapping. The two grid systems and their block corre-

spondence are shown in Fig. 3-1.

The grid in the physical domain at three different phases of piston movement

are shown in Fig. 3-2. One can observe from this figure that the grid is fixed for the

middle blocks corresponding to the connecting channel, but movable in the other two

blocks representing the end reservoirs of the thermal pump. The grid spacing in the

area of each reservoir near the connecting channel, where a cross flow can also enter

or leave the reservoir in they direction, is also independent of time. This subdivision

is of advantage for specifying the boundary conditions correctly and easily when using

cross flow. Note that the same grid structure was used for all cases in both con-

duction and cross flow models in order to allow a good comparison, even if it is not

necessary for the former. The grid points in the rest of the area of the reservoirs

move like an accordion in the x direction because of the piston oscillation. The y

components of the grids were established at the beginning of each run, while the x

components of the grids had to be generated at each time step.









(a)













(b)













(c)









Fig. 3-1 The correspondence of blocks between the physical domain and the
computational domain: (a) grid in the physical domain at ot=0; (b)
time independent and uniform grid in the computational domain; (c)
grid in the physical domain at ot=180.










(a)










(b)











(c)











Fig. 3-2 The movable and non-uniform grid in the physical domain:
(a) =t=00; (b) at=900; (c) ot=1800






36
In all the three subregions, the grid points are non-uniformly spaced such that

those regions near the walls and at the edge of the jet region at y=1 where large

velocity and/or temperature gradients are expected to exist have higher grid density.

In order to generate the non-uniform grid spacing, the stretching function found in

the book by Anderson21 was used. It has the form



(<+1) -(p-1)

(p +12 .-1 (3,1)




,where X and x are the position of grid points before and after stretching; 0 is the

stretching parameter which controls the stretching strength. P is different in diffe-

rent subregions and in different directions. The values of P's in the same direction

for different subregions were chosen so that the grid size nearest the dividing line

between subregions were equal.

The relationship between the coordinates (x,y,t) in physical domain and

(,, rl, r) in computational domain can be expressed as

= ((x,t)
n =1 (y) (3,2)
r=t

Under the condition that the product of the Jacobian and its inverse are equal to

unity, as is required for one to one mapping, the following expressions hold.











=1 1 x x (3,3)
x y7 xt


By means of these, the governing equations (2,6)-(2,9) in the physical domain

can, after some chain rule manipulations, be rewritten in terms of 4, r7 and r as


1- + -V, = 0 (3,4)
x Y7



x 1 1 t P 1 1 U (3,5)
U7-JU + (UU),+(VU), = Pg+-4-! U]+i-Uj



x 1 + (3,6)
x x y



T-XTT+(U), () 1 1 1+- (3,7)
x X y j. Pra2 x t x Ex P "y


Note that there is a new source term appearing in each of the momentum and

energy equations which represents the effect of the grid movement in the physical

domain.






38
3.2 A Brief Derivation of the Associated
Finite Difference Equations and Boundary Conditions

As in the SIMPLE algorithm, the same staggered grid, as shown in Fig. 3-3,

was also used in the SIMPLE-TP. The pressure and temperature are calculated at

the grid cell center designated by o's while the velocity components U and V are

determined at the center of the cell walls indicated by arrows. By means of this

staggered grid configuration, the checkerboard pressure effect introduced for non-

staggered grids can be eliminated and thus allows the Velocity-Pressure correction

to converge easily. In SIMPLE, a control volume approach was used to derive the

Finite Difference Equations (FDE's) in this grid configuration, and a derivation,

based on a stationary grid in the physical domain can, be found in the book by

Patankar19. The derivation of the FDE's in the present SIMPLE-TP followed the

same procedure, but some modifications had to be made since it was carried out in

the computational domain. A brief derivation of the x momentum equation is

presented below as an example of the derivation applied in the SIMPLE-TP.

The notation used here is the same as that of Patankar19 and is shown in Fig.

3-4. The location of U where the differential equation is going to be discretized is

designated by "P" while the location of its four neighboring Us are indicated by

"E"(East), "W'(West), "N"(North) and "S"(South), respectively. The four centers of

cell walls are marked by "e", "w", "n", and "s" also, but in lower case, according to

their directions away from the cell's center. Note that the position indices are placed

on the position of superscripts in order to differ from the partial derivatives which

are expressed by subscripts in the governing equations.
















j+l -1


0


I I -
-- 0 -- 0 0 --0 -

- t t t t -

o- 0 -- o -0



o 0 0
L ---


-I b


The staggered grid in the computational domain.


I I I I I


Fig. 3-3

















-9-


O

- O0


L--


-* O


n


P


1 i i+1


(b)


i-1 i i+1


Typical cells for the velocity component U: (a) an interior cell; (b) a
boundary cell with its south wall at the domain boundary.


0 -t-


j+1










j-
j-




j-l-


-a-


-r0-
" O -


Fig. 3-4


I I 1


I ) I 1 1 1 1


t






41
The equation (3,5) can be rewritten as


xhy Ur yrU + y PE + Y7 UU + x, V o- (3,8)


An integration of this equation over the grid cell S f6r and the time step Sr yields

6X6Y(UP'UO) L6Y(xPXO)(Ue-u) + 6y(Pe-PW)
6r ST
\ r (3,9)
U r U U U,
+ Sy uU- E uU- L] + Sx6vU- vU- = 0
a 2X: C 2X a } 2y


,where all U's have their values at the current time step except Uo which takes the

value of UP at the previous time step. The u and v here represent the velocity

components in the current time step also, but have their values at the previous itera-

tion and thus will go into the coefficients A's or source term S. It is in this way that

the nonlinear Navier-Stokes equations can be solved as linear equations in each step

of the iteration.

Note that the IUs are never located at the centers of the cell walls, thus the

U's with lowercase superscripts have to be replaced by those with uppercase super-

scripts (see Fig. 3-4(a)). The, so-called, hybrid scheme, which is a combination of

central differencing and upwind schemes described in detail by Patankar19, was used

for the convective and diffusion terms. Linear interpolation was used for the

replacement in other terms. Because of the staggered grid, PC and P" are available

without any further action.

The finally obtained FDE for U looks like









APUP =AEUE +AWUW +ANUN +ASUS +AOUO +Sc (3,10)

Where


AE = by 0, 1 uel + MAX(-u e,0) (3,11)
S a 2(XE -X P) 2


A by 0, 1 lul + MAX( u',) (3,12)
AW= y[ (, a2(XP-x W) 2



AN = 6x M 1 lvn + MAX(-v n,0) (3,13)
A a2( N y ) 2 )



Sa2y P -y S) 2


Ao Sxby (3,15)
8'

AP =AE +Aw N A N +As +A (3,16)


Sc = Y x o)(ue UW) (3,17)

Linear interpolation was also needed here to obtain the values of u and v at the

centers of the cell walls.

At the cells, which have one wall and sometimes two coinciding with the

domain boundary, the neighboring velocity component will locate at the center of the






43

cell wall. In Fig. 3-4(b), one of the boundary cells having its south wall coinciding

with the domain boundary is shown as an example. For all of the boundary cells,

proper adjustments to the right hand side in some of the equations (3,11)-(3,17) were

taken according to the geometry and boundary conditions specified. Particularly for

the boundary cell shown in Fig. 3-4(b), one has


As = 0, 1 -vI + MAX( v,0) (3,18)
x 2( P s) 2 1



The changes made in the superscript of yS and v' come from the change of the

location of Us. Equation (3,10) for this boundary cell should be changed also since

US here is known from the boundary conditions.

The derivation of all other FDE's and the treatment of other boundary cells

are straight forward and follow the same procedure explained above. More detail

and some discussion about the discretization scheme can be found in Ref. 19.





3.3 Calculation Procedure

Because of the assumptions of constant density and temperature independent

viscosity, the energy equation does not couple with the continuity and momentum

equations. This simplifies the problem. In this investigation, the velocity and

pressure field was first obtained without involving any temperature consideration.

The energy equation was then solved and the temperature distribution determined







44

based on the calculated velocity field. It was found that the temperature converges

much slower than the velocity does, which showed that a big saving in CPU time was

accomplished by being able to solve for velocity and temperature separately.

A flowchart for solving for U, V and P is shown in Fig. 3-5. The calculation

procedure was begun from an initial condition which was for all cases in this investi-

gation as that of the fluid at rest and the pressure uniform everywhere in the entire

physical domain. In the meantime, the dimensionless time rk=1 was assigned as zero

and a basic time step size rB was defined. Starting from this point, U, V and P at

tk+1= rk+ r were determined based on the grid generated in the physical domain

for this time step by means of the standard SIMPLE algorithm. In the iteration

process, the values for the dependent variables either at the previous time step or

those existing at the same rk but in the previous oscillation period, depending on

which was closer, were chosen as the initial guesses. Typically after stepping off

several oscillation cycles, the values at the same t in the previous cycle gave the most

rapid convergence.

The maximum ratio of the residue at each grid cell which came from substitu-

ting the non-converged velocity into the continuity equation (3,4) to the product of

the local velocity and the cell boundary size was used as a convergence indicator cc.

Whenever the value of a, was less than a prescribed small value e, the solution was

considered to have converged and one could proceed to the next time step. It was

found that sometimes the Velocity-Pressure correction was unable to achieve conver-

gence and the value of a, oscillated around a constant value larger than e, after


















































Fig. 3-5 Flowchart for determining the time periodic velocity components U, V
and pressure P.







46

several time steps, from the initial conditions or from a converged solution, were

taken. Under this situation, convergence could be reestablished by a continual

halving of the time step increments bt until they were sufficiently small. Once

convergence was again obtained, an attempt of doubling tb was taken in order to

use as large a time step size as possible.

After the time had been stepped off to r = t, a check for the convergence to

the final expected time periodic state was undertaken. This consisted of essentially

determining the quantities


1 imax ]max

a = imamax '=1 i=l (3,19)
(U j

i mamax i=1 j=1

and



mx max max

SYmaax i=1 j=1 (v(2=o)


ima max i=1 j=1

which represent the difference between the values of the velocity components at fixed

spatial position between period n and n+1. Note that the symmetric conditions

about the vertical center line have to be used in these formulas too. It was

considered that the periodic solution had been reached when both of them were less







47

than a given small value ep. If the solution was not periodic under the criteria

mentioned above, another cycle of calculation would be started by setting z = 0 and

6-r= 6rB.

An almost identical procedure was followed for determining the temperature

distribution based on the velocity obtained (see Fig. 3-6). The major difference here

from that for calculating U, V and P was that there were no iterations corresponding

to the Velocity-Pressure correction in each time step, thus the variable time step size

was not needed for T. Since the calculation for temperature was associated with the

final time periodic velocity even at the beginning of each run, there were no initial

conditions existing in this situation. A linear temperature distribution in the t

direction and uniform in rn direction was used as the initial guess at r = 0 of the first

period.





3.4 A Vectorized Line Group Method
for Solving the System of Algebraic Equations


The method for solving the associated system of algebraic equations always

play an important role in any overall numerical method for multi-dimensional flow

and heat transfer problems. A Line-By-Line iteration method (also called Successive

Line Relaxation method) is commonly used in the SIMPLE algorithm for 2D or 3D

problems. In this method, one picks a grid line, say in the t direction, and assumes

that the dependent variable at its neighboring lines (i.e. the '1 direction neighbors















































Fig. 3-6 Flowchart for solving the time periodic temperature distribution based
on the calculated velocity.







49

of the points on the chosen line) are known from the last iteration. By doing this,

a system of simultaneous equations with a tridiagnal coefficient matrix is formed,

which reads


A -A 0

-A A -A
3j 3j 3j




0 ... ...
O -- .. ...
o


0
-*- *" 0
0



W P E
-A A1M -AJ
-AIM-2j AIM-2j -AIM-2
w P
.... AIM-_ij AM-1j.


where A's are coefficients,
right hand side which can be expressed as


RHS2,

RHS3,


=W + N +A S +AO 0 -OC
= A2 (P+A2 2jP1+l 2j*P2j-1+A2 92j+S2

-N +AS ,0+ 0 ,+SC
3j *(P3jtl 3j* 3j-1 + 3j *3j +3j


(3,22)


RHSM-2j = AIM-2j* IM-2j++AIM-2j* qIM-2j- +AIM-2j*(IM-2j+SIM -2j

RHSIM J (PIMj l*' M2 N
RHSIM-lj = AAIM-,. pQMj+A It-lj* (IM-1j.+
AS -j ,0 0 o C
IM-lj ffIM-lj-I'llM-1j +A lM-l +jM-1j


This kind of system of algebraic equations can be solved by means of the Thomas

algorithm, a direct elimination method (see Anderson21). The forming and solving


(p2j

(p3j





(PIM-2j

(PIM-ljj


RHS2j

RHSzj





RHSIM-j

RHSIM-


(3,21)


> =







50

procedure, then, was used for all lines in the same direction (the E direction in our

example) and repeated until there was no significant changes in the dependent vari-

ables. The Line-By-Line method is widely preferred because it takes the advantages

of direct and iteration methods together to gain a fast convergence and good accu-

racy for the Velocity-Pressure correction iterations. More discussion concerning the

details of this method can be found in Ref. 19.

The above described method, however, is not suitable for vectorized calcula-

tions possible with supercomputers such as the CRAY and other newer machines.

Roughly speaking, vectorized calculation requires that there is no dependency

between the calculation operations in any pair of DO loop indices within the inner-

most DO loops. Once this requirement is satisfied, the calculations are performed

almost in parallel and the CPU time used is reduced accordingly (see Grenzsch22).

A flowchart for the Line-By-Line method with lines in Z direction is shown in Fig.

3-7. There are two innermost DO loops which are marked by the dashed lines in the

flowchart and both have "i" as their indices according to the line direction. It is easy

to see that Aj depends on A j, RHSi, is calculated from RHSi., in the first inner-

most DO loop and that (ij is obtained from (pi,j in the second innermost DO loop,

and that this procedure does not meet the requirement for a vectorized calculation.

A Vectorized Line Group method was developed in this study to overcome

the difficulties mentioned above especially since it was clear that the computations

for time periodic problems involve a great number of time steps, and thus some

attention have to be paid to the reduction of the CPU time used. Also a super-
















































Fig. 3-7 Flowchart for the Line-By-Line iteration method for solving the
associated system of algebraic equations. Dashed lines represent the
innermost DO loops.






52

computer resource (CRAY/YMP), employing vectorized processor, was available for

this study through the Pittsburgh Supercomputing Center.

By use of this new vectorized method, all lines in one direction, the f

direction for example, were divided into two groups (one for the even and the other

for the odd numbered lines) which lead to two groups of simultaneous algebraic

equations and a corresponding two sets of tridiagnal matrixes in the same way as in

the Line-By-Line method. The Thomas algorithm could then be applied to each line

in one group at the same time and the calculation could be vectorized. By altering

to update the dependent variable for the two groups, a converged solution could be

obtained some three times faster than that obtained by the Line-By-Line method.

The flowchart of the Vectorized Line Group method is shown in Fig. 3-8.

There are three innermost DO loops marked by dashed lines, in which the first and

third are corresponding to the innermost DO loops in the Line-By-Line method

shown in Fig.3-7. The dependencies between the pairs ofA's, RHS's and (p's are still

on the subscripts "i", the same as that in the Line-By-Line method, because the same

line direction is chosen. However, the indices of the innermost DO loops are "j",

which is not the same as that in the Line-By-Line method. It is because of this

difference that the requirement of vectorized calculation is satisfied. The only

drawback of the Vectorized Line Group method is that it requires a rectangular

computational domain. Fortunately this is achievable for most anticipated problems,

when using domain subdivisions.















































Fig. 3-8 Flowchart for the Vectorized Line Group iteration method for solving
the associated system of algebraic equations. Dashed lines represent
the innermost DO loops.







54
When the Vectorized Line Group method was applied to one of the three

subdomains shown in Fig. 3-1, the values of the dependent variables from previous

iteration in the neighboring subdomain were used as the boundary condition along

the interface boundary between the two subdomains. The same treatment was

undertaken for the all three subdomains. Since the solving procedure is iterative, the

continuity condition at such boundaries is insured by updating the values of depen-

dent variables there alteratively.











CHAPTER IV


RESULTS AND DISCUSSION



Ten cases with different values on the control parameters were investigated

in this numerical study. Those parameter values are listed in Table 4-1.


Table 4-1. The list of 10 cases studied in the numerical investigation.
Case b c Ax rs a Up Re V(
No. (cm) (cm) (cm) (cm/s) (cm/s)
1 2 5 0.4 0.013 0.5 0.003 1 0
2 2 5 5 0.167 1 0.125 50 0
3 2 5 6 0.2 3 1.35 540 0
4 2 5 20 0.667 3 4.5 1,800 0
5 2 5 6 0.2 10 15. 6,000 0
6 2 5 20 0.667 10 50. 20,000 0
7 2 5 40 1.333 10 100. 40,000 0
8 0.6 10 6 0.2 3 4.5 540 0
9 2 5 20 0.667 10 50. 20,000 1
10 2 5 20 0.667 10 50. 20,000 10


Note that the connecting channel has half width a = 0.1cm and length L = 30cm

for all cases. The working fluid in this study is assumed to be water with a kinematic







56

viscosity of v=0.01cm2/s, density p=lg/cm3, specific heat C,= 4.18J/gK, and Prandtl

number Pr=1 (corresponding to water at high temperature). The zero cross flow

(Vc=0) represents the conduction model, and Re= a2Ax/a is the corresponding

Reynolds number based on the maximum piston velocity Up= oaAx/2b and the

reservoir half width b, but expressed in terms of Womersley number a, tidal

displacement Ax and the channel half width a through the requirement of mass

conservation and the definition of Womersley number.

The numerical solutions of the governing equations described in Chapter II

obtained by use of the SIMPLE-TP algorithm explained in Chapter III for the ten

cases and some discussion about several interesting phenomena found from those

results will be presented in this chapter with some work done to validate the

numerical results at the beginning. Note that all dependent variables, the velocity

components U, V, the stream function 0i, which is calculated from U, V, the pressure

P, the Temperature T, as well as the independent variables x, y shown in the figures

of this chapter are dimensionless quantities which are non-dimensionalized by the

method shown in equations (2,1)-(2,3), although the marker is omitted as

mentioned in section 2.2.



4.1 Validation of Results

In any numerical process involving calculation by computers, errors always

arise. It is because of this that some estimate of the magnitude of these quantities

is included in all numerical studies.






57

Since the purpose of this dissertation is to obtain a qualitative guidance to the

usage of the thermal pump technique and no detailed experimental data for the time

periodic velocity patterns and temperature profiles in such configuration exist with

which to compare the numerical data, only a qualitative error analysis is carried out.

That is, to demonstrate that the numerical results obtained are qualitatively close to

the true solution of the governing equations for the specified boundary conditions.

The number of grid points chosen is very important for numerical studies of

heat transfer in fluid flows. If a few grid points are taken, the results will be unable

to display the correct structure of the velocity pattern and temperature profile. A

large number of grid points, on the other hand, will need too much computing time,

which may either waste the computer resources or require such long calculations that

they can not be handled by present computers. For time periodic problems, it is

more desirable to use a small number grid points or, in other words, a large grid cell

size because generally a smaller time step size is required for a smaller grid cell size

in order to achieve convergence. If the time step size is very small, the number of

necessary time steps per oscillating cycle will become very large, and thus an

extremely large number of the total time steps will be needed since several

oscillation cycles are required in order to reach the final time periodic state.

A relatively low number of 91x71 non-uniform grid points were used in this

numerical investigation. Fig. 4-1 compares the stream line patterns for Case No.6 at

several different oscillation phases for 91x71 grid points (left part of the figure) and

121 x01 grid points (right part). The procedure used to obtain the stream lines from













































grid: 91x71 grid: 121x101


Fig. 4-1 Stream line patterns for Case No.6 calculated with different number of
grid points: (a) grid 91x71; (b) grid 121x101. Increment between
neighboring lines is A or=20 in non-dimensional units.






59
the numerically calculated velocity components will be explained in the next section.

There is some difference between the two flow patterns but it is not significant. The

changes from one to the other are only in detail and the basic structure of the flows

does not change. The dimensionless basic time step size 6 r is ir/360 for the grid

system 91x71 with several halvings occurring from the basic size to Cr/1440 only in

the first oscillation cycle (see Section 3.3), while the smaller basic time increment

6 B= ir/1440 was needed for the 121x101 grid in order to get a converged solution.

The CPU time used for the latter is more than 4 times of that for the former for the

same convergence criterion. The slight changes in flow structure and the big saving

in calculation time lead to the decision to choose the sparser grid structure of 91x71

in the calculations for other cases.

The value of e, which is the criterion on the convergence indicator oc of the

Velocity-Pressure correction, is set at 3%. A typical number of Velocity-Pressure

corrections required at each time step in the SIMPLE-TP technique for Case No.6

is 11 for this setting in most oscillation cycles except the first one. In the first cycle,

especially in the early part of the cycle, it needs more Velocity-Pressure iterations to

obtain a converged solution for e =3%. The criteria for reaching the final time

periodic state was taken as 1% for the velocity and as 0.01% for the temperature.

Twelve oscillation cycles were typically needed for velocity convergence, while 200

cycles were required for temperature. There were three other criteria needed for

solving the system of algebraic equations associated with velocity, pressure and

temperature. The domain averaged ratio of the difference between calculated







60

dependent variables from two successive iterations to the values of the dependent

variables obtained from the current iteration was chosen as a convergence indicator.

The criteria for velocity, pressure and temperature were 1xlO4, 1x10l, and 1x10"

respectively. For velocity and temperature, less than 10 iterations were needed to

meet the criteria but usually more than 300 iterations were needed to have the

pressure converged. The numbers of iterations needed to satisfy these criteria

settings varied slightly from case to case for the different values of parameters listed

in Table 4-1. 3.5 hours of CPU time on a CRAY/YMP at the Pittsburgh Super-

computing Center was typically required for obtaining the velocity and pressure

values for Case No.6 using the 91x71 grid under the above converge criteria. An

extra 0.5 hour was needed to find the corresponding temperature field.

As shown in Section 1.3, there are closed form solutions for the velocity and

temperature profiles in oscillating flow in a parallel plate channel connected to end

reservoirs at x= + The flow is 1D in the axial direction and the variations of both

velocity and temperature along the transverse direction are independent of x. Note

that the ratio of the tidal displacement to the channel length, r,, in this situation is

always zero for whatever value of Ax except infinity. When the end reservoirs are at

a finite distance x= L/2, and thus r, is non-zero, a distortion of the velocity and

temperature profiles from what they are for infinite channel length will occur and the

above x independence will disappear. The distortion and the x dependence increase

when rL becomes large and under these conditions the solutions have to be obtained

numerically. On the other hand, when r, is small, the distortion should be small and






61
the numerical results for velocity profile and temperature distribution at the middle

of channel should be close to the analytic solution for the channel of infinite length

as long as other conditions are the same and, of course, the numerical results have

converged. Based on this point, a comparison of numerical results obtained in the

middle of connecting channel for the case of small r, with the analytic solution

under the same a and same insulating wall conditions for temperature will show

whether the numerical solutions are correct or not.

The first comparison involves the curves, shown in Fig. 4-2 and Fig. 4-3, for

the velocity component U and temperature T in Case No.5 versus the oscillation

phase ot for an entire oscillation cycle after the time periodic state has been reached.

The numerical results were taken at the position x= y=0.819 in the dimensionless

coordinate system shown in Fig. 2-1. A 5 phase shift was introduced in comparing

the analytic solutions with the numerical results in order to allow the best com-

parison. This phase shift comes from the difference of phase accounting methods in

the two solutions: the analytic solutions are based on the oscillation phase of the

pressure gradient, while the numerical results are based on the phase of the piston

oscillation.

The analytic solution is valid only for a channel of infinite length. A uniform

in space but time sinusoidal pressure gradient exists for this very simple geometry.

The analytic solving procedure starts from introducing this pressure gradient into the

Navier-Stokes equation (see Section 1.3), from which all time periodic solutions are

obtained. The phases of all other quantities, velocity, temperature etc., are delimited























40

30

20

i 10

CO 0
o
T -10

S -20

-30

-40 '" 1
0 60 120 180 240 300 360
ut


Fig. 4-2 Non-dimensional velocity component U for Case No.5 at the middle of
the connecting channel as a function of oscillation phase solid line numerical result; dashed line analytic solution.






63















0.15

0.10

S0.05

'0 0.00
Ci
S-0.05

-0.10

-0.15 l""" i l tii l lll ll n r
0 60 120 180 240 300 360
wt


Fig. 4-3 Non-dimensional temperature T for Case No.5 at the middle of the
connecting channel as a function of oscillation phase solid line numerical result; dashed line analytic solution.






64

from this pressure gradient template. For a geometry of finite channel length, the

pressure gradient is no longer uniform in the entire domain, but uniform in the

region far from the ends when rL is small enough, and thus can not be introduced

a prior at the beginning of the calculation. Instead, its behavior comes part of the

numerical results. The numerical solving procedure begins by giving the movement

of pistons, their position and velocity, in end reservoirs and thus the oscillation phase

of the piston oscillation becomes the basic template. The two templates do not

coincide with each other. In other words, the nearly uniform and time periodic

pressure gradient in the middle of connecting channel generated by oscillation of the

pistons in end reservoirs does not have the same phase as that of piston movement.

Introducing a phase shift is necessary for the best comparison.

One observes from the two figures (4-2 and 4-3) that the numerical results

(solid lines) match the analytic solutions (dashed lines) very well over most parts of

oscillation cycle. By comparing the corresponding zero points of the curves in these

two figures, one finds a 87, not a 90, phase delay of the temperature curves from

those of velocity. It is this difference of 30 in phase delay that makes the net heat

transfer in axial direction possible without an accompanying net mass transfer along

the axis. In the region near the x axis (center line), the phase delay is 90 so that

there is no net heat flux along the axis there (see the curve for a =10 in Fig. 1-9(a)).

However, for a=3, which is near the maximum point in Kff curve (see Fig. 1-6), the

phase delay between temperature and velocity is 750 at the x axis where the

difference from 900 is as much as 15.






65

The second comparison was made for the curves of velocity component U and

temperature T versus y. Fig. 4-4 and Fig. 4-5 show these curves at several oscillation

phases for Case No.5 with rL=0.2. The same phase shift was introduced when

comparing the analytic curves with numerical results. The curves in (a) come from

the analytic solutions (Eq. (1,4) and Eq. (1,6)), while those in (b) are from the

numerical results at x=0. One sees that the shape of the curves are essentially

identical except for the kinks produced by the finite grid point structure used in the

numerical calculations. The velocity curves for Case No.6 are presented in Fig. 4-6.

The curves are still almost the same in the region near the channel walls, but there

are larger differences between the numerical and analytic solutions near the center

line. Since rL=0.67 in this case, it makes sense that large differences should exist

in this case as explained above.

A comparison of curves of pressure gradient for Case No.5 versus oscillation

phase between numerical results atx=0 and that used for analytic solutions is shown

in Fig. 4-7. The dashed line again represents the analytic pressure gradient and the

same phase shift as above was introduced. The behavior of the numerical pressure

gradient does match the analytic one but there is an error larger than that for

velocity and temperature shown in Fig. 4-2 and 4-3, especially near the points at

which the pistons change their movement direction. It was found that this error

could not be reduced by reducing the grid size and/or time step size. Some

simplification is made in the SIMPLE algorithm when the pressure correction

equation is established from the continuity equation (explained in the Appendix, see











(a)


0.5 -


S0.0 -


-0.5


-1.0
-50 -40 -30 -20 -10 0
U(y,At)


(b)


10 20 30 40 50


0.5


0.0


-0.5


-1.0




Fig. 4-4


-50 -40 -30 -20 -10 0 10
U(O,y,wt)


20 30 40 50


Non-dimensional velocity component U for Case No.5 at the middle of
the connecting channel as a function of non-dimensional position y at
several oscillation phases: (a) analytic solution; (b) numerical result.











(a) 1.0
-analytic

0.5
ctt=0 30 60 90 120 150 180

> 0.0


-0.5


-1.0
-0.2 -0.1 -0.0 0.1 0.2
T(y,wt)



(b) 1.o
numerical

0.5 t-t=0 30 60 90 120 150 180


> 0.0 -


-0.5


-1.0
-0.2 -0.1 -0.0 0.1 0.2
T(O,y,wt)


Fig. 4-5 Non-dimensional temperature T for Case No.5 at the middle of the
connecting channel as a function of non-dimensional position y at
several oscillation phases: (a) analytic solution; (b) numerical result.











(a)


0.5


> 0.0


-0.5


- 1.0 ... .. ..
-50 -40 -30 -20 -10 0 10 20 30 40 50
U(y,wt)


(b)


0.5


> 0.0


-0.5


-1.0 -
-50



Fig. 4-6 Nc
the
sev


-40 -30 -20 -10 0 10
U(O,y,wt)


20 30 40 50


n-dimensional velocity component U for Case No.6 at the middle of
connecting channel as a function of non-dimensional position y at
reral oscillation phases: (a) analytic solution; (b) numerical result.





69
















5 0 I i i i i I i i i i l


25 _


x 0



-25


5 0 I I I I I I I I I I I I I I I I I I
0 90 180 270 360
ot


Fig. 4-7 Non-dimensional x direction pressure gradient for Case No.5 at the
middle of connecting channel as a function of oscillation phase solid line numerical result; dashed line analytic solution.






70
also reference 19). This simplifi-cation does not affect obtaining a good pressure

field for steady flows but may introduce more errors for unsteady and especially for

oscillatory flows.





4.2 Velocity Field

It was found that the velocity field under the conditions considered in this

study can not be well presented by the normal ways, i.e. by drawing arrows at several

points in the physical domain with length proportional to the speed and in the

velocity direction at that particular location, because of the complexity of the

oscillation flows and the large differences in speed at different locations within the

end reservoirs. Therefore, plotting the isovalue lines of the time varying stream

function, the stream lines, was chosen as the mode of velocity field presentation.

From the definition of stream function r(xy,t)


U(x,y,t) = a(x,y,t) V(x,y,t) = (x,y,t) (4,1)
ay x

,the value of the stream function (xy,t) at a particular point (xy,t) can be calculated

from r(xOyOt), the value of stream function at its neighboring point (xyet). The

following formula will do this.

(xo,y,)Y (xyt)
*(x,y,t) = *(xo,yo,t) + f U(x, Y y,t)dy f V(XO,y,t)dx (4,2)
(xoyor) (x,,yt)







71
Because of the staggered grid, the values of stream function are easily obtained at

the covers of grid cells without any interpolation (see Fig. 3-3).

The spacing between the streamline isovalues is equal to a constant which may

vary from case to case. This method of representing a time dependent flow field

allows one to take advantage of the fact that the absolute value of the instantaneous

streamfunction gradient is equal to flow speed and hence that regions of high stream-

line concentration (i.e. darkened regions) corresponding to parts of the flow field

having high velocity. The velocity is always in the tangential direction of the stream

lines, and the solid curves represent counterclockwise rotation, the dotted curves zero

rotation and the dashed curves clockwise rotation.

The stream lines in the end reservoirs at the oscillation phases of 150 intervals

for Case No.3 are presented in Fig. 4-8. Only the flow patterns at the phases in

lower half period Oast < 1800 are shown here. The flow profile at any point in upper

half period 180awt <3600 can be found by reflecting the picture at the corresponding

1800 earlier phase about the vertical center line (i.e. y axis). The constant interval

between two neighboring lines is 5 non-dimensional units.

It is observed that as the piston in the left reservoir moves toward the

connecting channel starting from its extreme position of x = -c- at at=0, there
2b
is essentially one large vortex pair present with the flow along the horizontal center

line (i.e. x axis) toward the piston surface. A weaker and smaller second vortex pair

of opposite rotation sense sits in the covers nearest the channel end. Flow is seen

leaving the left reservoir starting at about wt=300 and continuing on to about
















































Fig. 4-8 Stream line pattern for Case No.3 at oscillation phases of interval
A(wt) =15" in the lower half cycle *00ot<180*.






73

ot=1500 mainly via a side flow into the connecting channel. In the right reservoir,

one can see a distinct high velocity jet emanating from the connecting channel and

moving across the entire reservoir along the x axis. It begins at about at =30, the

same time as flow starts leaving the left reservoir, continues flowing along the

channel exit until at least ot=165. The jet collides with the right piston surface at

around ot =210 and disappears at ct =285. The images of these pictures are present

in the left reservoir at ot=300 and wt=1050, respectively. The origin for this high

velocity pulse is clearly the inability for the liquid entering the reservoir during the

suction stage to accommodate itself to the wider cross section of the reservoir as it

would if the Reynolds number were low enough so that flow separa-tion at the

entrance would not occur (see Case No.1 in Fig. 4-14). The appearance of a velocity

jet during the intake portion of the piston stroke has also been noted in the related

experimental work by Durst23 as well as by Kurzweg, Lindgren and Lothrop2. Note

that the large counter rotating vortex pair filling most of the reservoirs remains

during the entire oscillation cycle and its rotating sense is always such that fluid

motion along the x axis is always from the end of connecting channel and toward the

moving pistons, no matter what oscillation phase one is in. The smaller second

vortex pair also exists during the entire oscillation cycle, although it is much weaker

than the primary one. There is third vortex pair showing up in Fig. 4-8. It is

generated near the flow separation point near the channel end in the right reservoir

at about ot=450, shortly after the velocity jet starts, pushed toward the receding

piston by the velocity jet and finally joins the large vortex pair at about ot=135C.







74

The numerical development of the non-dimensional velocity component U at

three different locations along the x axis within the left reservoir under the conditions

corresponding to Case No. 3 are plotted in Fig. 4-9 with (a) at a position of 2a from

the channel end, (b) at 25a, and (c) at 40a. It shows that except at the position

nearest the connecting channel end, the velocity remains negative throughout the

oscillation cycle and is thus always moving the liquid toward the piston surface. The

sharp velocity spikes seen correspond to the periodic velocity jets crossing the

reservoir which are shown in the right reservoirs of Fig. 4-8 since they exist in the left

reservoir during the upper half period 1800 Mt <3600. The time delay in the velocity

peak between the non-dimensional position x =-2 and x =-25 is 0.0875 seconds and

thus indicates a velocity pulse traveling at 26.3cm/s which is only slightly less than

the rms axial velocity of 27cm/s in the channel but considerably larger than the

piston face velocity at the same instant.

The stream lines for Case No.4 and Case No.6 are shown in Fig. 4-10 and 4-

11, respectively, using the same format as that in Fig. 4-8 but employing 20 non-

dimensional units of interval in the isovalues of r. Similar to the results shown in

Fig. 4-8, there are three vortex pairs, two remaining during the entire oscillation cycle

and one existing in a portion of the cycle, and also a velocity jet in both cases.

For Case No.4(Fig. 4-10), the jet within the right reservoir begins at wt=30,

the same as in Case No.3, but collides with the piston surface at about ot=900, much

earlier than the wt =210 value of Case No.3, and disappears at about at=180, also

earlier than the wt=2850 result in Case No.3. The third vortex pair starts at about









a) 10

0

.-10 -

-20

-30

-40

-50 I
0 720 1440 2160



b) ,o



3 -10
0 -
n -20

S-30-

-40

-50
0 720 1440 2160
wt

C) 10

0

3 -10

0 -20 -
'
-30

-40

-50
0 720 1440 2160
wt

Fig. 4-9 Temporal approach to the final time periodic state of the non-
dimensional velocity component U(x,Ot) for Case No.3 at three
different locations: (a) x=-152; (b) x=-175; (c) x=-190.






















.ot=30






I!,-I i ? i



1 aOt=5600




oJ-13<75}
^^^-r'L ^ t


Fig. 4-10 Stream line pattern for Case No.4 at oscillation intervals of A (et) =15"
in the lower half cycle 00a












































Fig. 4-11 Stream line pattern of Case No.6 at oscillation intervals of A(wt) =15
in the lower half cycle 0*ot<180.







78

tot=30, almost at the same time as when the jet pulse begins and a little earlier than

at ot= 45 in Case No.3, and joins the big vortex pair at tf= 750, which is also much

earlier than wt=1500 phase point in Case No.3. The only difference between the

conditions in Case No.3 and Case No.4 is the different values of the tidal displace-

ment Ax. Case No.4 has Ax=20cm, or 10/3 times more than that in Case No.3, which

means there is a larger rms velocity in the channel for Case No.4. It is the faster

fluid movement in channel which produces the faster jet and associated shorter

appearance of the third vortex pair.

Comparing the stream lines of Case No.6 with those in Case No.4, one can not

find as large a difference as that between the stream lines of Case No.3 and Case

No.4. Note that Case No.6 has a=10 and Case No.4 has a=3, which is the only

difference between the two cases and means that the oscillation in Case No.6 is 10/3

times faster than that in Case No.4 since the fluid viscosity and channel width do not

change. The above observation reveals the fact that Womersley number no longer

has a large effect on the flow profile in the end reservoirs for large b/a ratio like its

influence on the axial velocity distributions within the connecting channel (see

Section 1.3). Instead, the tidal displacement Ax has more influence on the flow

profile in reservoirs than Womersley number does. Actually, significant changes in

the reservoir flow patterns are characterized by the Reynolds number, as will be

shown below.

Fig. 4-12 presents the pressure distribution for Case No.6 along the x axis at

several oscillation phases over the whole period O0sot<3600, and Fig. 4-13 shows the






79

corresponding stream lines at the same phase points. The pressure in most parts of

the channel is seen to vary linearly with x for most of the oscillation phases except

those near ot = 900 and wt = 270 when the pistons' movement changes from accelera-

tion into deceleration. The linear pressure distribution in the channel corresponds

to the expected constant pressure gradient there, and as used in the analytic solutions

(see Section 1.3). Examining these last two figures, one finds that the largest

pressure drop occurs at the channel end during that time when the fluid enters into

the channel through that end and when the velocity jet exists from the other end.

It is apparent that the pressure drop existing at one end of the connecting channel

provides the necessary power to form the jet flow at another end. The two ends alter

their roles as the pistons change their directions of motion. The magnitude of the

pressure gradient found in the channel agrees well with the analytic solution given

by Eq. (1,5) for the same values of parameters.

Stream lines with a spacing interval of 0.05 units for Case No.1 are presented

in Fig. 4-14. The corresponding Reynolds number in this case has the low value of

1, which is small enough to avoid flow separations in the case of steady flow over a

step. For a part of the oscillation cycle (about half), when the pistons move essen-

tially in one direction, no separations exists also in this oscillating flow situation. The

fact that flow separation does occur over rest of the oscillation cycle even for such

low Reynolds number is easily understood by noticing that this part of oscillation

period is around the points at which the pistons' motion changes its direction. The

flow pattern in Case No.l, in which both the existing vortex pairs and flow jet have

















0

X


-1


0
a0
X

0--


-250 -150 -0 50 150 250
-250 150 -50 50 1,50 250


0

x


-1


-250 -150 -60 50 150 250


0

-1 0

a -I


-z
-2 L I 1 1 i I
-250 -150 -50 50 160 250

X


2


1

0

x


-1


-250 -150 -50 50 150 250

X


2


1

0

x 0

a -1


I I


-250 -10 -60 50 10
-250 -150 -60 50 150 250


cat=135"

1-









-2-
-250 -150 -60 50 150 250


-50 -0 -0 50 150 50
-250 -150 -60 50 150 250


t=315

I


0
x

-14



-2 -
-250 -150 -60 50 150 250


Fig. 4-12 Non-dimensional pressure distribution along the x axis in Case No.6 at

oscillation intervals of A(ot) =45* over the whole cycle 0 art<3600.


r
ot=45' I
I
1\ I
I
I \ I
I \ I
I \ I
I \I
I


i I
I ot=225'
I
I
I
I / I
I / I
I I
I
I/ (


I ot=9O"



I I


I I

II I


nI I


I co270" I




I



I i


I I


^I I


[


e I I I


i
ot=O" I
I
i
i

I \ I
I
f
f I
I \I


ot=180'
I /,
I /(
I / (
I / I
I / 1
I / I
I / I
r/ I
I/ I


J


I


J


I























ii i:r- x=2250


----NI;,--


Fig. 4-13 Stream line pattern for Case No.6 at the oscillation intervals of
d4(t) =45* over the whole period 0*at<360*




























Ca=t=105 IAIII 9-W, -


-1.9

ot= 120(



-1.8<11<0


-fl450

-31.411<


()t=600


(ol=750

81 9

(ft=1500



-1 .3

Fig. 4-14 Stream line pattern for Case No.1 at the oscillation intervals of

A (t) = 150 in the lower half cycle Ostot <1800.


, -r--r c= 15 0
-'--- c2>


-0.6<%VO


01=900

r;

-r?~y

I


wt=1 350

-1H "as0.


--- I


- -


....._ __..~.i.~~-~I ~I~CPTPTTS~i~FCSLE


~S~IP~ i~DCIE~EIC~PZPPPPZ~


-L


Il;~prpi-in~-PEPm~liirr~ ~PTPPLV~----------