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.
ChiaShun Yih, Dr. Edward K. Walsh and Dr. Tom IP. 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 CBT8611254 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 VELOCITYPRESSURE CORRECTION .............. 110
REFERENCES ................................................ 116
LIST OF FIGURES
Figure M
Fig. 11 Heat transfer in a parallel plate channel: (a) pure conduction;
(b) enhancement by oscillating flows. .... 3
Fig. 12 A schematic of the Thermal Pump (U.S. Patent No. 4590993).
5
Fig. 13 Geometric configuration for the analytic solution of flow
oscillations and the associated heat transfer in parallel plate
channels. .... 10
Fig. 14 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. 15 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. 16 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. 17 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 liquidsolid
combinations. .... 19
Fig. 18 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. 19 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. 21 Configurations of the thermal pump used in the numerical
investigation: (a) conduction model; (b) cross flow model. .... 23
Fig. 31 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. 32 The movable and nonuniform grid in the physical domain: (a)
at=0; (b) at=90; (c) wt=180. .... 33
Fig. 33 The staggered grid in the computational domain. .... 36
Fig. 34 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. 35 Flowchart for determining the time periodic velocity
components U, V and pressure P. .... 42
Fig. 36 Flowchart for solving the time periodic temperature
distribution based on the calculated velocity. .... 45
Fig. 37 Flowchart for the LineByLine iteration method for solving the
associated system of algebraic equations. Dashed lines
represent the innermost DO loops. ... 48
Fig. 38 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. 41 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. 42 Nondimensional 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. 43 Nondimensional 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. 44 Nondimensional 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. 45 Nondimensional temperature T for Case No.5 at the middle of
the connecting channel as a function of nondimensional
position y at several oscillation phases: (a) analytic solution;
(b) numerical result. .... 61
Fig. 46 Nondimensional 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. 47 Nondimensional 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. 48 Stream line pattern for Case No.3 at oscillation phases of
interval A(ot) =150 in the lower half cycle 0tot<1800. .... 66
Fig. 49 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. 410 Stream line pattern for Case No.4 at oscillation intervals of
A(ot) =15 in the lower half cycle Oscot <180. .... 70
Fig. 411 Stream line pattern for Case No.6 at oscillation intervals of
A(ct) =150 in the lower half cycle 00 t<1800. .... 71
Fig. 412 Nondimensional 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. 413 Stream line pattern for Case No.6 at the oscillation intervals of
A(at) =450 over the whole period 00Oat<3600. .... 74
Fig. 414 Stream line pattern for Case No.1 at the oscillation intervals of
A(ct) =150 in the lower half cycle 00t<1800. .... 76
Fig. 415 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. 416 Stream line pattern for Case No.9 at oscillation intervals of
A (ot) =15 over the lower half period 00
Fig. 417 Stream line pattern for Case No.10 at oscillation intervals of
A (t) = 150 in the lower half cycle 00s t <1800. .... 80
Fig. 418 Temperature contours for Case No.3 at oscillation intervals of
A (at) =150 in the lower half cycle 00 t <180. .... 83
Fig. 419 Nondimensional 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. 420 Temperature contours for Case No.4 at oscillation intervals of
A(act) =150 in the lower half cycle 00ot<1800. .... 86
Fig. 421 Nondimensional 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. 422 Temperature contours for Case No.6 at oscillation intervals of
A(ot) =150 in the lower half cycle 0o
Fig. 423 Nondimensional 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. 424 Temperature contours for Case No.9 at oscillation intervals of
A(at) =150 in the lower half cycle 0t <1800. .... 91
Fig. 425 Nondimensional 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. 426 Temperature contours for Case No.10 at oscillation intervals of
A (t) = 150 in the lower half cycle 00ot <1800. .... 93
Fig. 427 Nondimensional 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 41
Table 42
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 nonuniform 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 VelocityPressure 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 fluidtosolid (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 fluidtosolid (wall) thermal diffusivities in analytic solutions or
general convergence indicator in numerical calculations
Ca convergence indicator for VelocityPressure 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 SIMPLETP
(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 noslip 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.
11(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. 11
4
Prandtl number 1, for which the thicknesses of momentum and thermal boundary
layers are equal, are shown in Fig. 11(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. 12. 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. 12
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. 13. 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(ba) 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(va)
where y is the coordinate perpendicular to the channel axis as shown in Fig. 13, 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(ba)
LA
\10\1~x
ic
ie
t
Fig. 13
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 nondimensional velocity, normalized by CoAx,
as a function of nondimensional transverse position in the y direction at different
oscillation phases for a=1 and a=10 is shown in Fig. 14(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 crosssection 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. 14 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(ViPa y/a)
gf (y) _[ cosh(V a) cosh(ia) (1,7)
2aa2(Pr1) 1 tanh( iJa)
hx 4ia
g Y) 1 H cosh[V7 Pr a(ey/a)]
2a2(Pr1) 1 tanh(Ta) cosh[VTaPr a (e1)] (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 fluidtosolid (wall) thermal conductivities, o= r/Kr
is the ratio of thermal diffusivities, and e=b/a.
The nondimensional temperature profiles T/yAx versusy/a for p = a=1 and
e=2 but different Womersley number a are shown in Fig. 15(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. 15 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 + (1H)1] + (Im H) + (1mH)
eff2 (1,12)
16 a2 (Pr2 1) 1 tanh(tfa)
d0x2 VXa
where
I = F/ca tanh(vTa) (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 nondimensional rff/a x2, as a function of Womersley
number for different Prandtl numbers, are presented in Fig. 16(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. 16(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 5x10s 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
103 10 a 10 ' 1 10 10 z 10 S
a
(b) 1
(b ) i n1  i iiiiil 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. 16 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. 17. 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 waterglass 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. 18, where a plot of ./foAx2e versus e for several values of a2Pr
for the liquidsolid combination of liquid lithium and stainless steel are presented.
10 a
101 C
102
Lx
3
10 '
10
10 6
10 3
Fig. 17
0.025
0.020
0.015
3 0.010
0.005
0.000
Fig. 18
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 liquidsolid 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. 19(a) and (b) show the behavior of
the two partially integrated, nondimensional products of U(y,t) and T(0,y,t) for the
case of a waterglass 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. 19(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. 19
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
VelocityPressure 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. 12 were chosen for the present numerical study. These are shown in Fig. 21.
Some simplifications were introduced for the purpose of computational convenience.
Note also that, to more clearly show details, the axialtransverse 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 2bZ
b
 tttttttt" tttttttt
Th Tc
Fig. 21 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. 21.
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. 21(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. 21(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 nondimensional 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 nondimensiona
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
SemiImplicit 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
SIMPLETP (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 SIMPLETP 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 nonuniform 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. 21, 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. 31.
The grid in the physical domain at three different phases of piston movement
are shown in Fig. 32. 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. 31 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. 32 The movable and nonuniform grid in the physical domain:
(a) =t=00; (b) at=900; (c) ot=1800
36
In all the three subregions, the grid points are nonuniformly 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 nonuniform grid spacing, the stretching function found in
the book by Anderson21 was used. It has the form
(<+1) (p1)
(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)
U7JU + (UU),+(VU), = Pg+4! U]+iUj
x 1 + (3,6)
x x y
TXTT+(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. 33,
was also used in the SIMPLETP. 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 VelocityPressure 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 SIMPLETP 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 SIMPLETP.
The notation used here is the same as that of Patankar19 and is shown in Fig.
34. 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. 33
9
O
 O0
L
* O
n
P
1 i i+1
(b)
i1 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
jl
a
r0
" O 
Fig. 34
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)(Ueu) + 6y(PePW)
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 NavierStokes 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. 34(a)). The, socalled, 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(XPx 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. 34(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. 34(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. 35. 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 nonconverged 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 VelocityPressure correction was unable to achieve conver
gence and the value of a, oscillated around a constant value larger than e, after
Fig. 35 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
6r= 6rB.
An almost identical procedure was followed for determining the temperature
distribution based on the velocity obtained (see Fig. 36). The major difference here
from that for calculating U, V and P was that there were no iterations corresponding
to the VelocityPressure 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 multidimensional flow
and heat transfer problems. A LineByLine 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. 36 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
AIM2j AIM2j AIM2
w P
.... AIM_ij AM1j.
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*P2j1+A2 92j+S2
N +AS ,0+ 0 ,+SC
3j *(P3jtl 3j* 3j1 + 3j *3j +3j
(3,22)
RHSM2j = AIM2j* IM2j++AIM2j* qIM2j +AIM2j*(IM2j+SIM 2j
RHSIM J (PIMj l*' M2 N
RHSIMlj = AAIM,. pQMj+A Itlj* (IM1j.+
AS j ,0 0 o C
IMlj ffIMljI'llM1j +A lMl +jM1j
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
(PIM2j
(PIMljj
RHS2j
RHSzj
RHSIMj
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 LineByLine 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 VelocityPressure 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 LineByLine method with lines in Z direction is shown in Fig.
37. 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. 37 Flowchart for the LineByLine 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 LineByLine 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 LineByLine method.
The flowchart of the Vectorized Line Group method is shown in Fig. 38.
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 LineByLine method
shown in Fig.37. The dependencies between the pairs ofA's, RHS's and (p's are still
on the subscripts "i", the same as that in the LineByLine 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 LineByLine 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. 38 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. 31, 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 41.
Table 41. 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 SIMPLETP 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 nondimensionalized 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 nonuniform grid points were used in this
numerical investigation. Fig. 41 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. 41 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 nondimensional 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
VelocityPressure correction, is set at 3%. A typical number of VelocityPressure
corrections required at each time step in the SIMPLETP 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 VelocityPressure 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 41. 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 nonzero, 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. 42 and Fig. 43, 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. 21. 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
NavierStokes 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. 42 Nondimensional 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
S0.05
0.10
0.15 l""" i l tii l lll ll n r
0 60 120 180 240 300 360
wt
Fig. 43 Nondimensional 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 (42 and 43) 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. 19(a)).
However, for a=3, which is near the maximum point in Kff curve (see Fig. 16), 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. 44 and Fig. 45 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. 46.
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. 47. 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. 42 and 43, 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. 44
50 40 30 20 10 0 10
U(O,y,wt)
20 30 40 50
Nondimensional velocity component U for Case No.5 at the middle of
the connecting channel as a function of nondimensional 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 tt=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. 45 Nondimensional temperature T for Case No.5 at the middle of the
connecting channel as a function of nondimensional 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. 46 Nc
the
sev
40 30 20 10 0 10
U(O,y,wt)
20 30 40 50
ndimensional velocity component U for Case No.6 at the middle of
connecting channel as a function of nondimensional 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. 47 Nondimensional 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 simplification 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. 33).
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. 48. 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 nondimensional 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. 48 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 separation at the
entrance would not occur (see Case No.1 in Fig. 414). 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. 48. 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 nondimensional 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. 49 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. 48 since they exist in the left
reservoir during the upper half period 1800 Mt <3600. The time delay in the velocity
peak between the nondimensional 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. 410 and 4
11, respectively, using the same format as that in Fig. 48 but employing 20 non
dimensional units of interval in the isovalues of r. Similar to the results shown in
Fig. 48, 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. 410), 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
S30
40
50
0 720 1440 2160
wt
C) 10
0
3 10
0 20 
'
30
40
50
0 720 1440 2160
wt
Fig. 49 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
oJ13<75}
^^^r'L ^ t
Fig. 410 Stream line pattern for Case No.4 at oscillation intervals of A (et) =15"
in the lower half cycle 00a
Fig. 411 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. 412 presents the pressure distribution for Case No.6 along the x axis at
several oscillation phases over the whole period O0sot<3600, and Fig. 413 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. 414. 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. 412 Nondimensional 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. 413 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 9W, 
1.9
ot= 120(
1.8<11<0
fl450
31.411<
()t=600
(ol=750
81 9
(ft=1500
1 .3
Fig. 414 Stream line pattern for Case No.1 at the oscillation intervals of
A (t) = 150 in the lower half cycle Ostot <1800.
, rr 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;~prpiin~PEPm~liirr~ ~PTPPLV~
