FLUID FLOW AND HEAT TRANSFER IN CAVITIES WITH INLET AND OUTLET
PORTS: EFFECT OF FLOW OSCILLATION AND APPLICATION TO DESIGN OF
MICROVALVES
Except where reference is made to the work of others, the work described in this
dissertation is my own or was done in collaboration with my advisory committee.
This dissertation does not include proprietary or classified information.
_____________________________________
Seyed Mahdi Saeidi
Certificate of Approval:
________________________ ________________________
James A. Guin Jay M. Khodadadi, Chair
Profesor Profesor
Chemical Engineering Mechanical Engineering
________________________ ________________________
Amnon J. Meir Robert L. Jackson
Professor Assistant Professor
Mathematics and Statistics Mechanical Engineering
______________________________
Stephen L. McFarland
Acting Dean
Graduate School
FLUID FLOW AND HEAT TRANSFER IN CAVITIES WITH INLET AND OUTLET
PORTS: EFFECT OF FLOW OSCILLATION AND APPLICATION TO DESIGN OF
MICROVALVES
Seyed Mahdi Saeidi
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of
the Requirements for the
Degree of
Doctor of Philosophy
Auburn, Alabama
December 16, 2005
iii
FLUID FLOW AND HEAT TRANSFER IN CAVITIES WITH INLET AND OUTLET
PORTS: EFFECT OF FLOW OSCILLATION AND APPLICATION TO DESIGN OF
MICROVALVES
Seyed Mahdi Saeidi
Permission is granted to Auburn University to make copies of this dissertation at its
discretion, upon the request of individuals or institutions and at their expense.
The author reserves all publication rights.
_________________
Signature of Author
________________
Date of Graduation
iv
VITA
Seyed Mahdi Saeidi, son of Mohammad and Zohreh Saeidi, was born on August
17, 1976, in Yazd, Iran. He graduated from Nemooneh High School, Yazd, in 1993. He
attended Sharif University of Technology, Tehran, Iran for four years and graduated with
a Bachelor of Science degree in Mechanical Engineering, specialized in Heat Transfer
and Fluid Mechanics, in September 1997. He then entered the Department of Mechanical
Engineering, Shiraz University, Shiraz and graduated with a Master of Science degree in
October 2000. After working as a CFD analyst in Iran-Khodro Powertrain Company (IP-
CO), Tehran for one and a half years, He attended Auburn University to study Doctor of
Philosophy in the Mechanical Engineering Department.
v
DISSERTATION ABSTRACT
FLUID FLOW AND HEAT TRANSFER IN CAVITIES WITH INLET AND OUTLET
PORTS: EFFECT OF FLOW OSCILLATION AND APPLICATION TO DESIGN OF
MICROVALVES
Seyed Mahdi Saeidi
Doctor of Philosophy, December 16, 2005
(Master of Science, Shiraz University, 2000)
(Bachelor of Science, Sharif University of Technology, 1997)
203 Typed Pages
Directed by Dr. Jay M. Khodadadi
A Computational Fluid Dynamics (CFD) code is developed and tested to solve the
two-dimensional form of the Navier-Stokes and energy equations. The flow field and
heat transfer in a square cavity with inlet and outlet ports is investigated thoroughly. The
effect of introducing a constant and oscillating inlet velocity on the forced convection in a
square cavity is studied. As a practical application of the generic problems studied, a
three-dimensional model of a piezoelectrically-actuated microvalve is presented.
The effects of the Reynolds number, the ports width, and position of the outlet
port on the steady-state fluid flow and heat transfer in a sealed cavity have been studied
vi
in great detail. The trends of the pressure drop, local and mean Nusselt numbers are
elucidated for a wide range of operating conditions. It is concluded that placing the outlet
port in parallel to the inlet port reduces the pressure drop significantly, whereas placing
the outlet port at the corners of the cavity increases the heat transfer rate in the cavity.
The transient evolution and periodic characteristic of the flow in a cavity with an
oscillating inlet velocity have been studied. To maximize heat transfer and minimize
pressure drop, the location of the inlet and outlet ports are chosen based on the result of
the previous study. A detailed discussion of the effects of the Reynolds number and the
frequency of the oscillating velocity is presented. The results indicate that both the flow
and thermal fields reach their periodic states after a certain period of time. The results
show that heat transfer rate greatly increases when the period of the oscillating inlet
velocity is close to convective time scale of the cavity (St ? 1).
Computational modeling of liquid flow in a NASA JPL (Jet Propulsion Library)
piezoelectrically-actuated microvalve is discussed. The three-dimensional velocity and
pressure fields are obtained for different deflections. By changing the mass flow rate at
the inlet, the pressure drop between the inlet and outlet ports is found and a loss
coefficient is determined for every deflection. The predicted pressure drop values are
compared to the experimental data for water flow within the microvalve. A correlation is
presented for the mass flow rate versus the pressure drop coefficient. The results show
that the pressure drop increases exponentially by decrease of the deflection.
vii
ACKNOWLEDGMENTS
The path to completing this PhD has included the discovery of new friends and
colleagues. I have many people to thank. First, I would like to express my deepest
appreciation to my advisor, Dr. Jay M. Khodadadi, for his patient, friendly, and unfailing
support over the past three and a half years.
My other dissertation committee members Drs. Meir, Guin, and Jackson have
provided excellent guidance with respect to my dissertation work. Conversations with
Dr. Meir have been beneficial and enjoyable; he has often inspired renewed interest in
my own research, and has helped me to look at old problems in new ways. I appreciate
Dr. Chris Roy, Aerospace Engineering Department, for his careful attention to my work.
Many thanks also go to the Alabama Supercomputer Center and Auburn
University College of Engineering Linux Cluster for the CPU time and technical support.
This research work was supported by Sam Ginn College of Engineering, Department of
Mechanical Engineering.
I would like to thank Dr. Yang and his colleges at NASA Jet Propulsion
Laboratory for providing the information and assistance with the liquid microvalve study.
Finally, gratitude is due to my parents for their inspiration and support during
these years. I am indebted to my family for encouraging me to pursue this degree, and
supporting me every step along the way.
viii
Style manual or journal used:
Guide to Preparation and Submission of Theses and Dissertations
Computer software used:
MS Word 2003
xiv
TABLE OF CONTENTS
LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . xix
LIST OF FIGURES. . . . . . . . . . . . . . . . . . . . . . . . . . . xx
NOMENCLATURE . . . . . . . . . . . . . . . . . . . . . . . . . . xxix
CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW. . . . . . . . . 1
1.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . 4
1.1.1 Review of Thermal Energy Storage Literature. . . . . . . . 4
1.1.2 Review of Microvalve Literature . . . . . . . . . . . . . 8
1.1.3 Review of Isolated Applications Literature . . . . . . . . .11
1.2 Closure and Structure of the Dissertation . . . . . . . . . . . . .15
CHAPTER 2 HIGH-ACCURENCY COMPUTATIONAL METHODS FOR
ELLIPTIC PROBLEMS AND VERIFICATION OF THE
CFD CODE . . . . . . . . . . . . . . . . . . . . . . . 40
2.1 Mathematical Description of Physical Phenomena . . . . . . . . 40
2.1.1 General Transport Equations . . . . . . . . . . . . . . 40
2.2 Staggered Grid System . . . . . . . . . . . . . . . . . . 41
xv
2.3 Discretized Equations . . . . . . . . . . . . . . . . . . . . 43
2.4 The Hybrid Differencing Scheme . . . . . . . . . . . . . . . 46
2.5 Quadratic Upwind Differencing Scheme: The QUICK Schemes . . . 47
2.5.1 QUICK Scheme of Leonard (1979). . . . . . . . . . . . 48
2.5.2 A Stable and Fast-Converging Variation of the QUICK Scheme:
The QUICK of Hayase et al. (1992) . . . . . . . . . . . 51
2.5.3 High-Order Approximation for Boundary Conditions . . . . 52
2.6 Benchmarking of the Code . . . . . . . . . . . . . . . . . . 53
2.6.1 2-Dimensional Duct Flow . . . . . . . . . . . . . . . 53
2.6.1.1 HYBRID Scheme. . . . . . . . . . . . . . . . 54
2.6.1.2 QUICK Scheme. . . . . . . . . . . . . . . . .54
2.6.1.3 Results for the 2-D Channel . . . . . . . . . . . 55
2.6.2 The Backward-Facing Step Flow. . . . . . . . . . . . . 57
2.6.2.1 Results for the 2-D Back-Step Flow . . . . . . . . 58
2.6.3 The Forward-Facing Step Flow. . . . . . . . . . . . . .62
2.6.3.1 Results for the 2-D Forward-Step Flow . . . . . . . 63
2.7 Closure . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
CHAPTER 3 STEADY-STATE FORCED CONVECTION IN A SQUARE
CAVITY WITH INLET AND OUTLET PORTS . . . . . . . . . 77
3.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . 78
3.1.1 Governing Equations . . . . . . . . . . . . . . . . . . . .79
xvi
3.1.2 Dimensionless Form of the Governing Equations . . . . . . . . . 80
3.2 Computational Details . . . . . . . . . . . . . . . . . . . . . . 81
3.2.1 Grid Independence Study and Benchmarking of the Code. . . . . . 82
3.2.2 Parameters for Numerical Simulation. . . . . . . . . . . . . . 83
3.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . 83
3.3.1 Flow Fields in the Cavity. . . . . . . . . . . . . . . . 83
3.3.2 Pressure Drop in the Cavity . . . . . . . . . . . . . . 86
3.3.3 Temperature Fields in the Cavity. . . . . . . . . . . . . 88
3.3.4 Variation of the Local Nusselt Number on the Walls of the
Cavity. . . . . . . . . . . . . . . . . . . . . . . .89
3.3.5 Variation of the Total Nusselt Number of the Cavity . . . . 90
3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
3.5 Closure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
CHAPTER 4 FLOW FIELD AND HEAT TRANSFER IN A CAVITY WITH
INLET AND OUTLET PORTS DUE TO INCOMING
FLOW OSCILLATION. . . . . . . . . . . . . . . . . . . 105
4.1 Problem Formulation . . . . . . . . . . . . . . . . . . . 106
4.1.1 Governing Equations . . . . . . . . . . . . . . . . 107
4.1.2 Dimensionless Form of the Governing Equations . . . . . 108
4.2 Computational Details. . . . . . . . . . . . . . . . . . . .110
4.3 Results and Discussions . . . . . . . . . . . . . . . . . . .111
xvii
4.3.1 Time Required to Reach the Periodic State . . . . . . . . 113
4.3.2 Periodic Evolution of the Temperature Field . . . . . . . .115
4.3.3 Periodic Evolution of the Temperature Field . . . . .. . . 116
4.3.4 Transient Evolution of the Mean Nusselt Number on Four
Walls. . . . . . . . . . . . . . . . . . . . . . . 117
4.3.5 Variation of the Total Nusselt Number of the Cavity . . . . 119
4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 120
4.5 Closure. . . . . . . . . . . . . . . . . . . . . . . . . . . 121
CHAPTER 5 COMPUTTIONAL MODELING OF A
PIEZOELECTRICALLY-ACTUATED MICROVALVE
FOR THE CONTROL OF LIQUID FLOWRATE. . . . . . . . .132
5.1 Geometry of the Microvalve Model . . . . . . . . . . . . . .133
5.2. Problem Formulation . . . . . . . . . . . . . . . . . . . 134
5.2.1 Possible Non-Continuum Effects. . . . . . . . . . . . 136
5.3 Radial Stokes Flow Between Two Parallel Disks . . . . . . . . .136
5.3.1 Determination of the Radial Velocity. . . . . . . . . . . 137
5.3.2 Determination of the Friction Factor. . . . . . . . . . . 140
5.3.3 Determination of the Total Pressure Loss. . . . . . . . . 141
5.4 Computational Details. . . . . . . . . . . . . . . . . . . .142
5.5 Results and Discussions . . . . . . . . . . . . . . . . . . .142
xviii
5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . .146
5.7 Closure . . . . . . . . . . . . . . . . . . . . . . . . . 146
CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS . . . . . . . . 167
REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
xix
LIST OF TABLES
Table 1.1 Total entropy generated in the cavity for various cases . . . . . . . 14
Table 2.1 The location of the eye of the vortices for the back-step problem . . . 59
Table 2.2 Comparison of the locations of computed and measured
separation and reattachment points for the back-step problem . . . . .60
Table 2.3 The location of the eye of the vortices for the forward-step
problem . . . . . . . . . . . . . . . . . . . . . . . . . 65
Table 2.4 The location of separation and reattachment points for
the forward-step problem. . . . . . . . . . . . . . . . . . . 66
Table 4.1 Number of cycles needed to reach the periodic states for
fluid flow (N
f
) and temperature fields (N
t
). . . . . . . . . . . . 115
Table 4.2 Mean Nusselt numbers on the four walls under steady
state conditions. . . . . . . . . . . . . . . . . . . . . . 118
Table 5.1 The Reynolds number at the outlet port of the microvalve
for different deflections at different total pressures . . . . . . . . 136
Table 5.2 Comparison of the mass flow rates (mg/s) obtained from different
analyses for deflection 1 ?m . . . . . . . . . . . . . . . . .143
xx
LIST OF FIGURES
Figure 1.1 Schematic diagram of a generic region with inlet/outlet ports and
bounded by fixed and flexible boundaries . . . . . . . . . . . . 16
Figure 1.2 (a) Schematic diagram of FOUP-nitrogen charge system with
plenum injector and (b) top view of the FOUP nitrogen
charge system showing the nitrogen inlet/outlet locations
and dimensions . . . . . . . . . . . . . . . . . . . . . . 17
Figure 1.3 Flow in the mid plane and heat transfer for inlet Re = 715
, Tin = 31
o
C and back wall stratified from 20 to 30
o
C
From top to bottom: (a) streaklines obtained by average
of 60 video frames, (b) velocity vectors obtained by
using PIV, (c) velocity vectors obtained from numerical
model, (d) heat flux (kW / m
2
) into the back wall obtained
from numerical model. . . . . . . . . . . . . . . . . . . . 18
Figure 1.4 Heat transfer as a function of inlet Re, stratified 20-30 C back
wall and 30
o
C inlet flow (Rosengarten et al., 1999) . . . . . . . . 19
Figure 1.5 Air-cooled enclosure with two different inlet-outlet port
locations ( d / L = 1 / 15) (Omri and Nasrallah, 1999) . . . . . . . 20
Figure 1.6 Cooling efficiency (a) when Re is varied and (b) when Ri is varied:
xxi
(solid symbols) for configuration A and (open symbols)
for configuration B. . . . . . . . . . . . . . . . . . . . . 20
Figure 1.7 Configuration B: (a) time sequences of instantaneous
temperature profile at the horizontal centerline, and
(b) the cooling efficiency over time . . . . . . . . . . . . . . 21
Figure 1.8 Schematic diagram of the cylindrical cavity studied . . . . . . . . 22
Figure 1.9 Strong efficiency using different fluids for each configuration . . . . 23
Figure 1.10 Schematic diagrams of various configurations . . . . . . . . . . 24
Figure 1.11 Average temperature, cooling efficiency, and Nusselt
number for Re = 50 . . . . . . . . . . . . . . . . . . . . . . 25
Figure 1.12 Average temperature, cooling efficiency, and Nusselt
number for Re = 500 . . . . . . . . . . . . . . . . . . . . 26
Figure 1.13 Stack piezoelectric type normally closed microvalve . . . . . . . . 27
Figure 1.14 Stack piezoelectric type three way microvalve . . . . . . . . . . .28
Figure 1.15 Microvalve array conceptual design. To modulate flow,
a voltage is applied between the microvalve diaphragm
and substrate, causing the diaphragm to collapse unstably
into contact with the substrate, thereby closing the
inlet port and restricting flow . . . . . . . . . . . . . . . . . 28
Figure 1.16 (a) Scanning electron microscope (SEM) photograph of the
fabricated microvalve array (15X) and (b) an anisotropically
etched inlet port (90X) . . . . . . . . . . . . . . . . . . . 29
Figure 1.17 Flow rate as a function of pressure differential across the
xxii
microvalve array for fully open arrays for the three diaphragm size 30
Figure 1.18 Flow rate as a function of number of microvalve open for
400 and 500 ?m arrays for a constant pressure differential of 10 kPa. 30
Figure 1.19 Predicted array flow rate as a function of pressure differential
compared with measured flow rate for the three microvalve
diaphragm sizes. Dashed lines are predicted flow rates
assuming a constant gap of 5?m . . . . . . . . . . . . . . . 31
Figure 1.20 Conceptual schematic diagram of the JPL piezoelectric microvalve . . 32
Figure 1.21 Measured flowrates for an actuated microvalve. Piezoelectric
actuation has been successfully demonstrated for inlet
pressures in the range of 0 ~ 1000 psi. Flow rates at inlet
pressures above 300 psi are above the measurement
range of the flow meter used in the tests . . . . . . . . . . . . . 33
Figure 1.22 Schematic diagram of the cavity and the solid body showing the
grid and boundary conditions for aspect ratio of 0.25 whith L = 0.05 . 34
Figure 1.23 Variation of Gr/Re
2
with exit port locations with aspect ratios
of (a) 0.25 and (b) 4. . . . . . . . . . . . . . . . . . . . . 35
Figure 1.24 Physical model for one heat source. . . . . . . . . . . . . . . 36
Figure 1.25 Physical model for two heat sources (case 2) . . . . . . . . . . . 37
Figure 1.26 A schematic diagram of the calculation domain shown the
grid and the boundary conditions(Shuja et al., 2000b) . . . . . . . 38
xxiii
Figure 2.1 2-D Staggered Grid System. . . . . . . . . . . . . . . . . . 67
Figure 2.2 1-D, Uniform, Staggered Grid Showing the Nodes
Involved in the Evaluation of ? on Cell Faces of a Control-Volume . 68
Figure 2.3 Non-uniform grid 100?80 . . . . . . . . . . . . . . . . . . 69
Figure 2.4 Streamlines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
Figure 2.5 Velocity Vectors . . . . . . . . . . . . . . . . . . . . . . 69
Figure 2.6 Velocity profiles in the developing zone of a channel
(HYBRID, Re = 100) . . . . . . . . . . . . . . . . . . . . 70
Figure 2.7 Growth and decay of the overshoot parameter along the X- axis
(HYBRID scheme, 100
?
80 grid) . . . . . . . . . . . . . . . 71
Figure 2.8 Developing of the centerline velocity along the X axis (100
?
80 grid). . 72
Figure 2.9 The normalized wall shear stress versus mesh refinement. . . . . . 72
Figure 2.10 Schematic diagrams of the backward-facing flow geometries . . . . 73
Figure 2.11 Streamlines for Re = 200 (Case A) with uniform inlet flow. . . . . .73
Figure 2.12 Streamlines for Re = 200 (Case A) with parabolic inlet flow . . . . 73
Figure 2.13 Streamlines for Re = 200 (Case B). . . . . . . . . . . . . . . 73
Figure 2.14 Streamlines for Re = 200 (Case C). . . . . . . . . . . . . . . 73
Figure 2.15 Streamlines for Re = 450 (Case A) with uniform inlet flow . . . . . 74
Figure 2.16 Streamlines for Re = 450 (Case A) with parabolic inlet flow . . . . 74
Figure 2.17 Streamlines for Re = 450 (Case B). . . . . . . . . . . . . . . 74
Figure 2.18 Streamlines for Re = 450 (Case C) . . . . . . . . . . . . . . . 74
xxiv
Figure 2.19 Schematic diagrams of the forward-facing flow geometries . . . .75
Figure 2.20 Streamlines for Re = 200 (Case A) with uniform inlet flow . . . . . 75
Figure 2.21 Streamlines for Re = 200 (Case A) with parabolic inlet flow . . . . 75
Figure 2.22 Streamlines for Re = 200 (Case B) . . . . . . . . . . . . . . .75
Figure 2.23 Streamlines for Re = 450 (Case A) with uniform inlet flow. . . . . .76
Figure 2.24 Streamlines for Re = 450 (Case A) with parabolic inlet flow . . . . 76
Figure 2.25 Streamlines for Re = 450 (Case B) . . . . . . . . . . . . . . .76
Figure 2.26 Streamlines for Re = 1000 (Case A) with uniform inlet flow . . . . .76
Figure 2.27 Streamlines for Re = 1000 (Case A) with parabolic inlet flow. . . . .76
Figure 2.28 Streamlines for Re = 1000 (Case B) . . . . . . . . . . . . . . 76
Figure 3.1 Schematic diagram of a cavity with inlet and outlet ports . . . . . . 94
Figure 3.2 The absolute value of the stream function at the
center of primary vortex,?
min
, as a function of the
grid size (Re = 1000, S
o
= 2.1 and W = 0.2). . . . . . . . . . . .95
Figure 3.3 Streamlines for Re = 500 and W = 0.25 for 9 different
positions of the outlet port . . . . . . . . . . . . . . . . . .96
Figure 3.4 Streamlines for Re = 40 and W = 0.15 for 9 different
positions of the outlet port . . . . . . . . . . . . . . . . . .97
Figure 3.5 1 Streamlines for Re = 100 with three outlet ports of widths
W = 0.05, 0.15 and 0.25 positioned at three positions . . . . . . . 98
Figure 3.6 Streamlines for (a) W = 0.15 and S
o
= 0.925, and (b)
W = 0.25 and S
o
= 0.15, for Reynolds numbers of 10,
xxv
40, 100 and 500 . . . . . . . . . . . . . . . . . . . . . . 99
Figure 3.7 Variations of the pressure drop coefficient of the cavity
for different positions of the outlet port with different Re numbers:
(a) W = 0.05, (b) W = 0.15 and (c) W = 0.25 . . . . . . . . . . 100
Figure 3.8 Temperature fields for Re = 500 and W = 0.25 for 9 different
positions of the outlet port [contour level increment of 0.05] . . . .101
Figure 3.9 Temperature fields for (a) W = 0.15 and S
o
= 0.925, and (b)
W = 0.25 and S
o
= 0.15, for Reynolds numbers of 10, 40,
100 and 500 [contour level increment of 0.05]. . . . . . . . . . 102
Figure 3.10 Variations of the local Nusselt number on the four walls
of the cavity for different positions of the outlet port with
Re = 500 and W = 0.25 . . . . . . . . . . . . . . . . . . 103
Figure 3.11 Variations of the total Nusselt number of the cavity for different
positions of the outlet port with different Re numbers:
(a) W = 0.05, (b) W = 0.15 and (c) W = 0.25 . . . . . . . . . . 104
Figure 4.1 Schematic diagram of a cavity with inlet and outlet ports. . . . . . 122
Figure 4.2 Velocity oscillation digram . . . . . . . . . . . . . . . . . 123
Figure 4.3 Convection and diffusion time scale in relation to the period
of the oscillation . . . . . . . . . . . . . . . . . . . . . .123
xxvi
Figure 4.4 Comparison of the mean Nusselt number of the system
for various time steps (St = 10) . . . . . . . . . . . . . . . . 124
Figure 4.5 Dependence of the (a) streamlines and (b) temperature
contours for the steady case with Re = 100, 300 and 500
[contour level increment of 0.05 for the temperature field] . . . . . 125
Figure 4.6 Cycle-to-cycle changes of the streamline values over the
entire domain for different Strouhal numbers . . . . . . . . . . 126
Figure 4.7 Cycle-to-cycle changes of the temperature values
over the entire domain for different Strouhal numbers. . . . . . . 127
Figure 4.8 Periodic evolution of the streamlines for St = 0.5 during the 13
th
cycle for a phase angle increment . . . . . . . . . . . . 128 6/?
Figure 4.9 Periodic evolution of the temperature field for St = 0.5 during
the 21
th
cycle for a 6/? phase angle increment . . . . . . . . . 129
Figure 4.10 Transient evolution of the instantaneous mean Nusselt
numbers on the four walls for St = 0.1. . . . . . . . . . . . . 130
Figure 4.11 Transient evolution of the instantaneous mean Nusselt
numbers on the four walls for St = 10 . . . . . . . . . . . . . 130
Figure 4.12 Transient evolution of the instantaneous total Nusselt
number of the system: (a) St = 0.1 and 0.5, and (b) St = 1, 2 and 10 . 131
Figure 5.1 Structure of the liquid-compatible microvalve (Yang et al., 2005) . . 147
Figure 5.2 Operating principle of the microvalve in actuated ?on? state
(Yang et al., 2005) . . . . . . . . . . . . . . . . . . . . 148
Figure 5.3 images of the (a) valve seat, (b) upper boss plate (c) valve
xxvii
sealing surface of the lower boss plate and (d) upper boss
bonding side of the lower boss plate (Yang et al., 2005) . . . . . . 149
Figure 5.4 Fully assembled and packaged microvalve (Yang et al., 2005). . . . 151
Figure 5.5 One quarter of the microvalve structure . . . . . . . . . . . . 152
Figure 5.6 Isometric view of the microvalve geometry used for
the computer modeling. . . . . . . . . . . . . . . . . . . 153
Figure 5.7 Isometric view of the 3-dimensional computational Mesh . . . . . 154
Figure 5.8 Isometric view of the 3-dimensional computational mesh . . . . . 155
Figure 5.9 Isometric view of the 3-dimensional computational mesh . . . . . 156
Figure 5.10 Isometric view of one of the rings . . . . . . . . . . . . . . .157
Figure 5.11 Geometry of the radial flow between two parallel disks . . . . . . 158
Figure 5.12 Isometric views showing pathlines of selected particles
released at the inlet plane for a deflection 2.7 ?m and
pressure drop of 2.4 psi . . . . . . . . . . . . . . . . . . . 159
Figure 5.13 Streamlines colored by stagnation pressure in the inlet pipe
and 10-micron box . . . . . . . . . . . . . . . . . . . . 160
Figure 5.14 Streamlines colored by stagnation pressure over the rings
and the outlet hole. . . . . . . . . . . . . . . . . . . . . 161
Figure 5.15 The comparison between two and three dimensional numerical
results and experimental data (Yang et al. 2005). . . . . . . . . 162
Figure 5.16 Pressure drop versus mass flowrate for a 3.5 ?m deflection
(Yang et al., 2005) . . . . . . . . . . . . . . . . . . . . 163
xxviii
Figure 5.17 The comparison between three-dimensional numerical
results with different deflections . . . . . . . . . . . . . . . 164
Figure 5.18 Pressure drop coefficient for different deflections. . . . . . . . . 165
Figure 5.19 The measured deflections as function of applied
voltage(Yang et al., 2005) . . . . . . . . . . . . . . . . . 166
xxix
NOMENCLATURE
English Symbols
C
p
pressure drop coefficient
D
h
hydraulic diameter,m
ER aspect ratio of step
f frequency of the velocity?s oscillation, 1/s or Darcy friction factor, defined by
equation 5.17
h half of the distance between the disks, m
H height of the cavity, m
K pressure drop coefficient of microvalve, defined by equation 5.22
L length of the cavity, m
Nu average or mean Nusselt number
p pressure, N/m
2
P Dimensionless pressure, defined as P = p / ?U
2
P
t
Total pressure, N / m
2
Pr Prandtl number, defined as Pr = ?/?
Q volumetric flow rate, m
3
/s
Re Reynolds number
s coordinate adopted for distance along the walls, m
S dimensionless coordinate, i.e., S = s/H
xxx
St Strouhal number, St = H?f/U
m
t time, s
t* dimensionless time, t* = U?t/H
t
conv
convection time scale, s
t
diff
diffusion time scale, s
T temperature, K
u x-direction velocity, m/s
u
in
inlet velocity m/s
u
m
mean velocity
u
o
amplitude of oscillating velocity
U Dimensionless velocity in the X-direction
v y-direction velocity, m/s
V Dimensionless velocity in the Y-direction
w width of the inlet/outlet ports, m
W dimensionless width of the inlet/outlet ports, W = w / H
Greek Symbols
? thermal diffusivity, m
2
/s or phase angle
? diffusion coefficient
?
?
cycle-to-cycle variation of temperature field
?
?
cycle-to-cycle variation of the stream function
? dimensionless temperature, i.e., ? = (T-T
in
)/(T
w
-T
in
)
xxxi
? dynamic viscosity, kg/s?m
? kinematic viscosity, m
2
/s
? dimensionless coordinate
? density, kg/m
3
? period, s
?
rz
wall shear stress, N/m
2
? generic variable
? stream function value
?
c
stream function value at the center of the primary vortex
Subscripts
0 reference value
b related to the bottom wall
i related to the inlet or inner position
in related to the inlet
l related to the left wall
m mean value
o related to the outlet or outer position
r related to the right wall
t related to the top wall
tot total value
w related to the wall
1
CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW
Forced convection within a lid-driven sealed cavity and buoyancy-driven
convection within a differentially heated cavity stand out as simply-stated academic
problems. For many decades, these two problems have provided excellent opportunities
for researchers to test and benchmark computational tools for solution of complex flows
with recirculation. In addition to their geometrical simplicity, both of these flows offer
simultaneous existence of diverse flow regimes involving boundary layers and multiple
separated regions. To this end, great attention has been focused on these flows and
variations of them by Dr. Khodadadi?s research team at the Fluid Mechanics Research
Laboratory in the Mechanical Engineering Department, Auburn University. These
studies have focused on:
1. Steady fluid flow and heat transfer in a lid-driven square cavity due to a thin fin
(Shi and Khodadadi, 2002),
2. Transient evolution to periodic fluid flow and heat transfer in a lid-driven
square cavity due to an oscillating thin fin (Shi and Khodadadi, 2004, and 2005),
3. Laminar natural convection heat transfer in a differentially-heated square cavity
due to a thin fin (Shi and Khodadadi, 2003).
The logical extension of the previous work is to concentrate on more complicated
flow situations in cavities through introduction of multiple inlet and outlet ports. These
2
flow systems offer a wider range of interesting possibilities in comparison to the lid-
driven and differentially heated cavities and have a great number of engineering
applications.
From an academic viewpoint, any fluid system that is enclosed by solid/flexible
enclosure surfaces and has inlet/outlet ports is a generic model for a variety of practical
systems (Figure 1.1). These engineering systems can be utilized for mixing, separation,
compression/expansion of various fluid constituents and might entail more complicated
physical processes involving phase change, chemical reactions, etc. For instance, mixed
convection studies within cavities with inlet and outlet ports that contain heat sources are
of great interest in electronic cooling, ventilation of buildings and design of solar
collectors. A related topic of interest is the thermal energy storage inside enclosures with
inlet and outlet ports. Simulations of unsteady mixed convection in thermal storage
applications have received a great deal of attention for many decades. The relevant
literature review of thermal energy storage within enclosures having inlet and outlet ports
will be presented in Section 1.1.1.
More recently, transport phenomena at micron-size scales has attracted the
attention of researchers. Many microfluidic devices are characterized by enclosures that
are supplied by various entering fluid streams and in turn have many outlet fluid streams.
In this connection, design of a great number of microfluidic devices (microfabricated
fluid filters, piezoelectrically-actuated micropumps and microvalves) can also benefit
from the knowledge that is gained from the results of the present dissertation. A review
of literature on microvalves will be given in Section 1.1.2.
3
Other isolated applications related to purge of enclosures in wafer manufacturing,
mass transfer in cavities, relaxation tanks employed in petroleum and transformer
industries, transient removal of contaminated fluid from cavities and food processing are
also noteworthy. For instance, control of Airborne Molecular Contamination (AMC)
plays an increasingly important role in semiconductor manufacturing processes. A useful
method for reducing AMC is purging of the wafers with nitrogen gas. A typical
configuration of purging a Front Opening Unified Pod (FOUP), as part of a wafer box for
handling 300-mm wafers is shown in Figure 1.2 (Hu and Taso, 2005). Although the
geometrical details of this application are complicated, the salient features of the results
of the current study can be used for design of such practical systems. A review of the
isolated applications will be given in Section 1.1.3.
It is clear that a basic understanding of the complicated flow regimes and the
associated thermal fields can be of a great interest to a wide range of practical
engineering systems. In this dissertation, results of studies of:
1. Steady forced convection in a cavity with one inlet and one outlet ports,
2. Transient flow field and heat transfer in a cavity due to an oscillating inlet
velocity,
3. Steady flow in a liquid microvalve,
will be presented. In Chapter 2 a stable fast-converging variation of the QUICK scheme
is discussed. Three classic benchmark problems are also solved to verify the finite-
volume-based code that is developed for this study. The first two studies are fundamental
in nature and will be covered in Chapters 3 and 4. The more applied study of flow in a
4
microvalve is covered in Chapter 5. The relevant literature review for the above
problems are given in the remainder of this Chapter.
1.1 Literature Review
The literature review is presented in three sections. The first section is dedicated to
the previous work that is related to the thermal energy storage improvement with focus
on problems that include cavities with inlet and outlet ports. In the second part, some of
the studies on microvalve fabrication and modeling are reviewed. Finally, the isolated
applications will be reviewed.
1.1.1 Review of Thermal Energy Storage Literature
With rising energy costs and an increasing demand for renewable energy sources,
Thermal Energy Storage (TES) systems are becoming an interesting option. TES is a key
component of any successful thermal system and a reliable TES should allow minimal
thermal energy losses, thus leading to reduced energy costs. Thermal energy storage is
considered an advanced energy technology, and there has been an increasing interest in
using this essential technique for applications such as heating, hot water utilities, air
conditioning, etc. A very common application of TES is in solar-energy-based systems.
Solar radiation is a time?dependent energy source with an intermittent character. The
heating demands of a residential house are also time-dependent. However, the energy
source and the demands of a building, in general, do not match each other, especially in
solar heating applications.
5
In the following sections, some of the previous work that have been carried out to
improve heat transfer within enclosures are reviewed. Most of the reviewed articles are
related to the problems with inlet and outlet ports.
Chan et al. (1983) studied transient, two-dimensional, mixed convection flows in
a thermal storage tank using the Marker And Cell (MAC) method to solve the
Boussinesq-based governing equations. Their results showed that the device can store
energy at a faster rate when hot water is charged into the tank from the top and cold water
is extracted from the bottom. However, the directions of the inlet and outlet streams,
which can be either vertical or horizontal, were found to have negligible effect on the
thermal storage efficiency.
Homan and Soo (1998) studied the steady flow of a wall jet issuing into a large-
width cavity for which the primary axis is normal to the direction of the jet inflow.
Numerical solutions of the two-dimensional Navier-Stokes equations were obtained for
the inlet Reynolds numbers of 10 to 50 and tank width to inlet height ratios of 16 to 128.
The length and velocity scales of the wall jet boundary layer exhibit close agreement with
the classic wall jet similarity solution and published experimental data. However, the
width of the region for which the comparison proves to be favorable has a limited range.
This departure from a self-similar evolution of the wall jet is shown to result from the
finite domain width and its influence on the large recirculation cell located immediately
above the wall jet boundary layer.
An experimental and numerical investigation was undertaken to study the heat
transfer process in horizontal mantle heat exchangers used in solar water heaters by
Rosengarten et al. (1999). Three-dimensional flow simulation was performed as well as
an experimental study by using the Particle Image Velocimetry (PIV) technique. The
experimental and numerical results for the inlet Reynolds number of 715 are shown in
Figure 1.3. Note that the inlet port is placed at the left bottom corner of the cavity,
whereas the outlet port is at the right bottom corner. The results of this study suggest that
most of the heat transferred was in the bottom half of the cavity. In Figure 1.4 the
variation of the heat transfer as a function of the inlet Reynolds number is presented. The
heat transfer data show proportionally with Re
0.5
that is the typical behavior in forced
boundary layer convective heat transfer.
Omri and Nasrallah (1999) presented a numerical study of transient mixed
convection of laminar flows in an air-cooled cavity in which air was injected at a
temperature lower than the initial temperature of the cavity. The relevant geometries are
shown in Figure 1.5. A Control Volume Finite Element Method (CVFEM) using
triangular elements was employed to discretize the governing equations. In order to
investigate the quality of the ventilation and cooling, the dynamic and thermal fields were
numerically determined for a large range of the Reynolds and Richardson numbers and
for different inlet-outlet locations. The effect of the Reynolds and Richardson numbers
on the thermal efficiency of the system are shown in Figure 1.6. The thermal efficiency
is defined as:
.1
in
out
T
T
?=? (1.1)
These authors also examined the cooling efficiency for both transient and steady regimes
of operation. The results of this part of their study are presented in Figure 1.7. The
6
7
results showed that air injected from the bottom of the hot wall was more effective in
removing of the heat from the cavity.
Sensible heat storage in fluids generates thermal stratification. In order to
improve thermodynamic system efficiency, stratification should be promoted
substantially. To this end, Bouhdjar and Harhad (2002) presented a numerical study of
transient mixed convection using different fluids as heat storage media in cylindrical
cavities with different aspect ratios. The effect of the working fluid is investigated
through the variation of physical properties represented through the Prandtl number. The
system consists of a cavity with fluid injection at the top and fluid discharge at the bottom
(Figure 1.8). Three representative fluids were investigated, i.e. Torada oil, ethylene
glycol and water. These authors studied cavities with aspect ratios varying from 3 to 1/3.
Flow analysis was made through typical transient temperature distributions for the three
fluids and for different configurations. The performance of thermal energy storage using
these fluids are analyzed through the transient thermal storage efficiency. In Figure 1.9,
the thermal efficiency of the system is shown at different times for Prandtl numbers, 204
and 3.1 and five aspect ratios (0.33, 0.5, 1, 2 and 3.3).
A numerical study was conducted by Singh and Sharif (2003) to investigate
mixed convective cooling of a two-dimensional rectangular cavity with differentially
heated side walls. The geometries of Singh and Sharif (2003) are shown in Figure 1.10.
The horizontal walls were assumed to be adiabatic. The cold fluid was blown into the
cavity from an inlet on the side wall of the cavity and leaves through an outlet on the
opposite side wall. This configuration of mixed convective heat transfer has application
in building energy systems, cooling of electronic circuit boards, and solar collectors,
8
among others. The objective of the research was to optimize the relative locations of
inlet and outlet in order to have the most effective cooling in the core of the cavity by
maximizing the heat-removal rate and reducing the overall temperature in the cavity.
Various placement configurations of the inlet and outlet were examined for a range of the
Reynolds and Richardson numbers. For a given Reynolds number, the Richardson
number was varied from 0, which represents pure forced convection, to 10, which implies
a dominant buoyancy effect. The cooling efficiency, average temperature, and local and
average Nusselt numbers on the hot wall versus the Richardson number for the Reynolds
numbers 50 and 500 are presented in Figures 1.11 and 1.12, respectively. It is observed
that maximum cooling effectiveness is achieved if the inlet is kept near the bottom of the
cold wall, while the outlet is placed near the top of the hot wall.
1.1.2 Review of Microvalve Literature
The impetus for control of fluid flow on a small scale has motivated the recent
development of microvalves and micropumps using well-established microelectronic
fabrication techniques. The resulting fluidic MicroElectroMechanical System (MEMS)
devices are potentially important in applications that include biomedical dosing,
biochemical reaction systems, microfluidic mixing and regulation and the design/control
of microspacecrafts/microsatellites.
Microvalves developed so far can be classified in two categories:
1. Active microvalves (with an actuator),
2. Passive microvalves (without an actuator).
9
The active microvalves can be classified based on the type of their actuators. The
actuator in a microvalve is a component that applies the force to the moving part of the
valve to control the volumetric flow rate. There are normally-closed and normally-open
microvalves when the actuator is at rest. There are many different forms of actuation
with the most common being electric, magnetic, pneumatic, thermopnuematic and
piezoelectric actuations. In this review, the focus is on piezoelectrically-actuated
microvalves.
A normally closed microvalve using a stack-type piezoelectric actuator was
designed by Esashi (1990). The schematic geometry is illustrated in Figure 1.13. It has a
glass-silicon structure with a movable diaphragm. To avoid leakage, a free standing Ni
gasket was used. Nitrogen gas flow rate of 0.1 to 85 ml/min at a gas pressure of 7.5 kPa
was measured. A three-way microvalve was fabricated modifying this structure as shown
in Figure 1.14 by Nakagawa et al. (1990). A sample injector was fabricated using a pair
of these valves. Shoji and Esashi (1998) fabricated a microvalve having a large stroke for
liquid flow control using a bimorph-type piezoelectric actuator. Very small dead volume
microvalves were demonstrated using stack-type piezoelectric actuators. In order to have
disposable flow channel, it was separated from the actuator part. Since the actuator part
can be re-used many times, it reduces the cost and is useful for medical applications
(Shoji et al., 1991). A flow rate of 24 ?l/min at 10.1 kPa (0.1 atm) was achieved. The
on/off ratio of this microvalve was about 14. The total size of the fabricated microvalve
is about 8 mm ? 10 mm ? 10 mm, whereas the size of the stack of the piezoactuator is 2
mm ? 3 mm ? 9 mm. These authors also proposed another type of microvalve which can
10
be driven by a low power actuator such as piezo-disk. The valve was designed to reduce
the influence of the input pressure and to fit the targeted medical application.
Vandelli et al. (1998) developed a microvalve and investigated the pressure drop
versus mass flow rate of the system both numerically and experimentally. The device
consists of a parallel array of surface-micromachined binary microvalves working
cooperatively to achieve precision flow control on a macroscopic level. The schematic
diagram of each microvalve is shown in Figure 1.15, where the structure of the
microvalve arrays is shown. The anisotropically-etched inlet port is shown in Figures
1.16a and 1.16b, respectively. Flow rate across the microvalve array is proportional to
the number of microvalves open, yielding a scalable high-precision fluidic control
system. Device design and fabrication, using a one-level polycrystalline silicon surface-
micromachining process combined with a single anisotropic bulk etching process were
detailed. Performance measurements on fabricated devices confirm feasibility of the
fluidic control concept and robustness of the electromechanical design. Air-flow rates of
150 ml/min for a pressure differential of 10 kPa were obtained for diaphragm size of 300
?m (Figure 1.17). As it is seen in Figure 1.18, a linear flow control was achieved over a
wide range of operating flow rates. In this figure the mass flow rate versus the number of
open valves is sketched for 400 and 500 ?m arrays. A continuum fluidic model based on
the incompressible low Reynolds number flow theory was implemented using a finite-
difference approximation. The model accurately predicted the effect of microvalve
diaphragm compliance on the flow rate. In Figure 1.19, a comparison between the
numerical and experimental results is demonstrated. This graph shows a good agreement
11
between the theoretical predictions and experimental data over the entire range of flow
conditions tested experimentally.
Yang et al. (2004) fabricated and tested a normally-closed leak-tight
piezoelectrically-actuated microvalve for high-pressure gas micropropulsion (Figure
1.20). The microvalve operates at extremely high upstream pressures (0 ~ 1000 psi) for
microspacecraft applications. In order to control the mass flow rate, the boss plate moves
by the force applied by a piezoelectric actuator. This distance between the boss and seat
plate is called displacement or deflection and varies from 10 to 14 microns. The
measured flow rate of the microvalve for different applied displacement and inlet total
pressures is shown in Figure 1.21. The applied voltage determines the displacement
(solid line). The measured flow rate at an actuation voltage of 10 V, is approximately 52
SCCM at an inlet total pressure of 300 psig. A 1-D computational/theoretical analysis
has been conducted by Johnson (2005) to find the total pressure drop versus mass flow
rate in the gap between the seat and boss plates.
1.1.3 Review of Isolated Applications Literature
Mixed convection in a cavity with an internal heat source has received
considerable attention by researchers because of its importance in many engineering
applications. Cooling of electronic equipment is an application of this type of research.
The electronic parts are cooled down by air that blows through the equipment by a fan.
The heat transfer is a combination of free convection (buoyancy-driven flow) and forced
convection (pressure driven flow). The ratio of Gr/Re is introduced to determine the
relative importance of free convection and forced convection. Another application of the
2
12
same flow is in solar collectors. In solar collectors the heat source is solar energy and the
objective is to maximize the heat transfer rate from the sun radiance to the working fluid.
The working fluid in these systems is chosen based on the type of the application. Oil is
the most common fluid in the solar power plants, whereas air is commonly used in solar
air conditioning systems. Water is also chosen as the working fluid to provide the hot
water in houses by using solar energy. In this section some of studies that were carried
out to investigate heat transfer in enclosures with inside heat sources and inlet and outlet
ports are reviewed.
Combined free convection and forced convection in a partially divided enclosure
with a finite-size heat source was studied by Hsu et al. (1997). The enclosure was
partially divided by a conductive vertical divider protruding from the floor or the ceiling
of the enclosure. The emphasis in this research was placed on the influence of the
location and the height of the divider, as well as the locations of the source and the
outflow opening. The computation was carried out for wide range of the Reynolds
number (50 ~2000), whereas Gr/Re
2
varied from 10
-3
to 10. The results indicated that the
average Nusselt number and the dimensionless surface temperature on the heat source
depend on the location and the height of the divider. The heat transfer rate increased
when the heat source was close to the inlet port.
Shuja et al. (2000a) investigated flow in a cavity with a heat generating body.
The heat transfer characteristics and the irreversibility generated in the cavity depend on
mainly the cavity size, the aspect ratio of the heat generating body, and inlet/outlet port
locations (Figure 1.22). The effect of exit port locations on the heat transfer
characteristics and irreversibility generation in a square cavity with heat generating body
13
was investigated. A numerical simulation was carried out to predict the velocity and
temperature fields in the cavity. To examine the effect of the solid body aspect ratio on
the heat transfer characteristics, two extreme aspect ratios (0.25 and 4.0) were considered
in the analysis. Fifteen different locations of exit port are introduced while air was used
as an environment in the cavity. The results showed that non-uniform cooling of the
solid body occurs where the exit port location is on the left side of the bottom wall of the
cavity, whereas the inlet is placed at the bottom of the left wall (the outlet port number 13
and above in Figure 1.23). In this case, heat transfer reduces while irreversibility
increases in the cavity. These findings are valid for both aspect ratios of the solid body.
Hsu and Wang (2000) carried out a numerical study on mixed convective heat
transfer in an enclosure with discrete heat sources. The physical models with one and
two heat sources were studied (Figures 1.24 and 1.25). The discrete heat sources were
embedded on a vertical board, which was positioned on the bottom wall of an enclosure.
An external airflow enters the enclosure through an opening in one vertical wall and exits
from another opening on the opposite wall. This study simulated a practical system, such
as air-cooled electronic devices with heated elements. The emphasis was placed on the
influence of the governing parameters, such as the Reynolds number, Re, the buoyancy
parameter, Gr/Re
2
, location of the heat sources, and the conductivity ratio, r
k
, on the
thermal phenomenon in the enclosure. The coupled equations of the simulated model
were solved numerically using the Cubic Spline Collocation (CSC) method. The
computational results indicated that both the thermal field and the average Nusselt
number (Nu) depend strongly on the governing parameters, position of the heat sources,
as well as the property of the heat-source-embedded board.
Conjugate heat transfer in a cavity is an important consideration with regard to
cooling of micro-electronic equipment. Shuja et al. (2000b) carried out a heat transfer
analysis of conditions taking place in a square cavity with a heat source located in it.
Natural convection within the fluid in combination with conduction heat transfer inside
the heat generating solid body was examined. Air or water were considered as the fluid
in the cavity while a steel substrate was considered as the heat generating solid body.
The location of the solid was changed in the cavity to examine the cooling conditions
(Figure 1.26). The entropy analysis of the system was carried out to determine the
irreversibility ratio for each location of the solid body in the cavity. It was found that the
heat transfer from the solid body surfaces increases where the surfaces are facing the inlet
and the exit of the cavity. In Table 1.1 the total entropy generated by the cavity for
various cases is presented. The irreversibility ratio is defined as the ratio of volumetric
entropy ( ) generated by fluid friction and heat transfer:
'''
S
.
'''
'''
transferheat
frictionfluid
S
S
=? (1.2)
The entropy generated attains the maximum value for air when the solid body is located
at the center of the cavity. In this case, the irreversibility ratio reduces to a minimum
value.
14
Table 1.1 Total entropy generated in the cavity for various cases, Shuja et al. (2000b)
1.2 Closure and Structure of the Dissertation
In this research, attention was focused on the flow field and heat transfer in
simple enclosures with inlet and outlet ports. A highly-accurate CFD software (Chapter
2) is employed to investigate the effects of the width of the inlet and outlet ports, the
position of the ports and the period of the inlet velocity oscillation on the forced
convection in an enclosure. A complete parametric study is performed in order to gain a
full understanding of the convection processes. In light of the academic nature of this
research, the results of this study are significant contribution to the understanding of the
fluid flow and heat transfer in an enclosure with inlet and outlet ports. The results for
steady and oscillating convection regimes are presented in Chapters 3 and 4, respectively.
In Chapter 3, a steady problem is solved for a cavity with one inlet and one outlet port.
The effects of the Reynolds number, different outlet positions and the inlet and outlet
ports width on the heat transfer and pressure drop of the cavity are investigated. Upon
determining the best position of the outlet port to achieve maximum heat transfer rate and
minimum pressure drop from Chapter 3, the effect of oscillation of the inlet flow on the
total rate of the heat transfer on the walls of the cavity was investigated in Chapter 4.
The inlet velocity oscillates sinusoidally with different frequencies of oscillation. The
15
16
task of this study was improving the heat transfer in the cavity by tuning a proper
frequency of oscillation. In Chapter 5 a practical problem is solved. In this chapter the
flow field within a NASA JPL (Jet Propulsion Laboratory) piezoelectrically-actuated
microvalve is studied. Although the geometry of the microvalve is very complicated, the
basic concepts are similar to the generic cavity problem in Chapter 3. In Chapter 5 great
attention is paid to find an empirical equation for the required pressure of the system as a
function of mass flow rate. Chapter 6 is dedicated to the conclusions of this research
work and some suggestions for future effort in continuation of the present study.
17
Flexible Boundary
Outlet
Inlet
Fixed Boundary
Inlet/Outlet
Inlet/Outlet
Figure 1.1 Schematic diagram of a generic region with inlet/outlet ports and bounded by
fixed and flexible boundaries
Wafer
Outlet
Support
structure
Pl enum
Inlet(slit)
(a)
FOUP
Hole B
(1.5cm)
Hole A
(1.5cm)
Hole D
(1.5cm)
Hole C
(1.5cm)
90?X
60?X
Wafer
Support
Structure
Wafer
(b)
Figure 1.2 Schematic diagrams of (a) FOUP-nitrogen charge system with plenum injector
and (b) top view of the FOUP nitrogen charge system showing the nitrogen inlet/outlet
locations and dimensions (Hu and Tsao, 2005)
18
Figure 1.3 Flow in the mid plane and heat transfer for inlet Re = 715, T
in
= 31
o
C and
back wall stratified from 20 to 30
o
C. From top to bottom: (a) streaklines obtained by
average of 60 video frames, (b) velocity vectors obtained by using PIV, (c) velocity
vectors obtained from numerical model, (d) heat flux (kW / m
2
) into the back wall
obtained from numerical model (Rosengarten et al., 1999)
19
Figure 1.4 Heat transfer as a function of inlet Re, stratified 20-30
o
C back wall and 30
o
C
inlet flow (Rosengarten et al., 1999)
20
Figure 1.5 Air-cooled enclosures with two different inlet-outlet port locations (d / L = 1 /
15) (Omri and Nasrallah, 1999)
Figure 1.6 Cooling efficiency (a) when Re is varied and (b) when Ri is varied: (solid
symbols) for configuration A and (open symbols) for configuration B (Omri and
Nasrallah, 1999)
21
Figure 1.7 Configuration B: (a) time sequences of instantaneous temperature profile at
the horizontal centerline, and (b) the cooling efficiency over time (Omri and Nasrallah,
1999)
22
Figure 1.8 Schematic diagram of the cylindrical cavity studied (Bouhdjar and Harhad,
2002)
23
Figure 1.9 Thermal efficiency using different fluids for each configuration (Bouhdjar and
Harhad, 2002)
24
Figure 1.10 Schematic diagrams of various configurations (Singh and Sharif, 2003)
25
Figure 1.11 Average temperature, cooling efficiency, and Nusselt number for Re = 50
(Singh and Sharif, 2003)
26
Figure 1.12 Average temperature, cooling efficiency, and Nusselt number for Re = 500
(Singh and Sharif, 2003)
27
Figure 1.13 Stack piezoelectric type normally closed microvalve (Esashi et al., 1990)
28
Figure 1.14 Stack piezoelectric type three-way microvalve (Nakagawa et al., 1990)
Figure 1.15 Microvalve array conceptual design. To modulate flow, a voltage is applied
between the microvalve diaphragm and substrate, causing the diaphragm to collapse
unstably into contact with the substrate, thereby closing the inlet port and restricting flow
(Vandelli et al., 1998)
29
Figure 1.16 (a) Scanning electron microscope (SEM) photograph of the fabricated
microvalve array (15X) and (b) an anisotropically etched inlet port (90X) (Vandelli et al.,
1998)
30
Figure 1.17 Flow rate as a function of pressure differential across the microvalve array
for fully open arrays for the three diaphragm size (Vandelli et al., 1998)
Figure 1.18 Flow rate as a function of number of microvalve open for 400 and 500 ?m
arrays for a constant pressure differential of 10 kPa (Vandelli et al., 1998)
31
Figure 1.19 Predicted array flow rate as a function of pressure differential compared with
measured flow rate for the three microvalve diaphragm sizes. Dashed lines are predicted
flow rates assuming a constant gap of 5?m (Vandelli et al., 1998)
32
Figure 1.20 Conceptual schematic diagram of the JPL piezoelectric microvalve (Yang et
al., 2004)
33
Figure 1.21 Measured flowrates for an actuated microvalve. Piezoelectric actuation has
been successfully demonstrated for inlet pressures in the range of 0 ~ 1000 psi. Flow
rates at inlet pressures above 300 psi are above the measurement range of the flow meter
used in the tests (Yang et al., 2004)
34
Figure 1.22 Schematic diagram of the cavity and the solid body showing the grid and
boundary conditions for aspect ratio of 0.25 with L = 0.05 (Shuja et al., 2000a)
35
Figure 1.23 Variation of
2
Re
Gr
with exit port locations with aspect ratios of (a) 0.25 and
(b) 4 (Shuja et al., 2000a)
36
Figure 1.24 Physical models for one heat source (Hsu and Wang, 2000)
37
Figure 1.25 Physical model for two heat sources (case 2) (Hsu and Wang, 2000)
38
Figure 1.26 A schematic diagram of the calculation domain showing the grid and the
boundary conditions (Shuja et al., 2000b)
39
CHAPTER 2 HIGH-ACCURACY COMPUTATIONAL METHODS
FOR ELLIPTIC PROBLEMS AND VALIDATION OF THE CODE
The numerical solutions of elliptic problems involving fluid flow, heat and mass
transfer, and other related processes are based on the laws governing the transport
phenomena, i.e. conservation of mass, momentum, energy and chemical species, etc. All
these laws can be expressed in terms of differential equations. All these equations possess
a common form, thus making it possible to construct a general solution procedure for
elliptic problems. This chapter deals with the development of high-accuracy schemes
within the finite-volume methodology.
2.1 Mathematical Description of Physical Phenomena
2.1.1 General Transport Equations
There are significant commonalities among all equations governing the transport
phenomena. Upon introduction of a general variable ?, the conservative form of the
governing equations, including equations for scalar quantities such as temperature and
pollutant concentration etc., can be written in the following form:
()
().
?
???
??
Sgraddivudiv
t
+?=
?
?
?
?
?
?
+
?
?
(2.1)
40
Equation (2.1) is the so-called transport equation for genetic property, ?. Every term
corresponds to various transport processes. The rate of change term (unsteady term) and
41
the convective term are on the left hand side, whereas the diffusive term (? being the
diffusion coefficient) and source term are found on the right hand side. By setting ? equal
to 1, u, v, w and T or h and selecting appropriate values for ? and source term, one
obtains the continuity, momentum and energy equations. All the terms except unsteady,
convective and diffusive terms can be lumped into the source term. For instance, for the
Navier-Stokes equation, the source terms involve the pressure gradient. The
commonalities of all transport equations and flexibility of source term make it possible to
develop a general numerical procedure for the transport phenomena.
2.2 Staggered Grid System
The finite-volume method starts with the discretization of all the governing
equations within the chosen calculation domain. First, one needs to decide where to store
the velocities, pressure field values and all other related dependent variables. It seems
logical to define theses variables at the same point in space, that is referred to as the
"colocational" approach. One can easily show that if the velocity components and
pressure are stored at the same location, a wavy pressure field can act like a uniform field
in the discretized momentum equations. A remedy for this problem is to use a staggered
grid system. Such a staggered grid system was first used by Harlow and Welch (1965) in
their MAC (Marker And Cell) method.
Figure 2.1 is an illustration of a 2-D staggered grid system. The locations for the
discrete values of the u and v velocity components are shown in Figure 2.1 by arrows,
whereas the grid points (main grid points) are indicated by circles. The dashed lines
indicate the control-volume faces for the main grid points. One can see that the velocity
42
components are calculated at points that lie on the face of the control volume for
pressure. The u velocity component is only staggered in the x-direction, whereas the v
velocity component is only staggered in y-direction with respect to main grid points. A 3-
D staggered grid system can be constructed in a similar approach.
There are several advantages to the use of the staggered grid system.
1. The mass flow rate across the main control-volume faces can be calculated without any
interpolation of the relevant velocity component.
2. Adoption of a staggered grid system can prevent a wavy velocity field. Given a
staggered grid arrangement, the velocity components on the faces of the P control-
volume are used in the discretized continuity equation. Given this, a wavy velocity field
can not satisfy the continuity equation. As a result, this will lead the solution to a
reasonable velocity distribution.
3. Use of a staggered grid system can prevent a wavy pressure field. The pressure
difference between two adjacent main grid points are used for the momentum equation
and the wavy pressure field would no longer be accepted as uniform pressure fields and
could not be a possible solution.
Because of the above reasons, the staggered grid system was used for the research
activities reported here. On the other hand, the staggered grid system has some
disadvantages. Since the discrete values of u, v and p are stored at different locations,
extra work is needed for grid set-up and more space is required to store all the grid
information. Secondly, since u and v are not located at the same point, extensive
interpolation is generally required in the post-processing procedure. For high-accuracy
calculations, the interpolation procedure also needs to be of the same or a higher order
scheme to get reliable computed fields.
2.3 Discretized Equations
The most important step in the utilization of the finite-volume method is the
discretization of the general transport equation. The control volume integration is the key
step of the discretization procedure that distinguishes it from other CFD techniques. The
spatial integration of Equation (2.1) over a control volume and a finite time step
temporally, combined with the Green's Theorem yield the following relation:
() () ()
?? ?????
?
?+?+
?
?+?+
+
?
?
?
?
?
?
?
?
?
?
??=
?
?
?
?
?
?
?
?
?
?
?+
?
?
?
?
?
?
?
?
?
?
?
?
V
tt
t
tt
tV
tt
tAA
tt
t
dVdtSdtdAgradndtdAundVdt
t
?
?????
, (2.2)
where
n
and dV are the normal direction vector and the volume of the control volume,
respectively. Up to this point, no approximation or assumptions are introduced. In order
to evaluate the integrals in Equation (2.2), various differencing schemes can be
employed. In general, utilization of a high-order scheme will lead to a high-accuracy
solution.
For the 2-D case, our attention is focused on the main grid point P and the control
volume for P is shaded using crossed-lines in Figure 2.1. The neighboring nodes are
identified by W, E, N, S and the control volume faces by w, e, n, s. It is convenient to
define two variables F and D (both with dimensions ML
-2
t
-1
) that represent the convective
mass flux per unit area and the diffusion conductance at cell faces, respectively. So we
have:
43
,voruF ??= (2.3a)
.
y
or
x
D
??
??
= (2.3b)
If the ? value at node P is assumed to prevail over the whole control volume in the
unsteady term and the central-difference scheme is used for both convective and diffusive
terms, Equation (2.2) can be written as:
() ()()()()
,
0
??
?
?+?+
?+
?+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?+?
?
?
?
?
?
?
?
???
?
?
?
?
?
?
?
?=
?+?+??
tt
t
tt
t
sn
we
tt
t
snwePP
VdtSdt
y
A
y
A
x
A
x
A
dtvAvAuAuAV
?
????
???????????
(2.4)
where superscript 0 in the first term on the left hand side means the value of ? at the
previous time step. Quantity
?
S
is the source term. As suggested by Patankar (1980), the
average value
?
S
can be expressed as:
PPC
SSS ?
?
+=
, (2.5)
where S
C
stands for the constant part of
?
S
, whereas S
P
is coefficient of ?
P
.
For a uniform grid (i.e. ), the ? value at the cell
face can be written as:
1,...,2,1,
1
?==
+
nixx
ii
??
.
22
22
PS
s
NP
n
PW
w
EP
e
??
?
??
?
??
?
??
?
+
=
+
=
+
=
+
=
(2.6)
44
The diffusive flux on the cell faces can be evaluated as follows:
.
SP
SP
ss
s
PN
PN
nn
n
WP
WP
ww
wPE
PE
ee
e
y
A
y
A
y
A
y
A
x
A
x
A
x
A
x
A
?
???
?
???
?
???
?
???
?
?=
?
?
?
?
?
?
?
?
?
?
?
?
?=
?
?
?
?
?
?
?
?
?
?
?
?
?=
?
?
?
?
?
?
?
?
?
?
?=
?
?
?
?
?
?
?
?
?
(2.7)
Upon substituting Equations (2.6) and (2.7) into Equation (2.4), the fully implicit 2-D
dicretization equation takes the final form:
.
00
CpNNSSEEWWPP
Saaaaaa +++++= ??????
(2.8)
It should be noted that a similar relation can be derived for 3-D problems relating the
value at the P cell to the neighboring 6 points.
It should be noticed that the values of all coefficients depend on the differencing
scheme for the unsteady term, convective term and diffusive term. If the implicit Euler
scheme is used for the unsteady term and the central differencing scheme is used for both
the convective and diffusive terms, the corresponding coefficients are as follows:
()
()
()
()
.
22
22
22
22
0
0
0
nn
n
PN
nn
nN
ss
s
SP
ss
sS
ee
e
PE
ee
eE
ww
w
WP
ww
wW
snwe
P
P
PPSNEWP
Av
A
y
F
Da
Av
A
y
F
Da
Au
A
x
F
Da
Au
A
x
F
Da
FFFFF
t
V
a
SFaaaaaa
?
?
?
?
?
?
?
?
?
?
?
=?=
+
?
=+=
?
?
=?=
+
?
=+=
?+?=?
?
?
=
??+++++=
(2.9)
45
2.4 The Hybrid Differencing Scheme
Non-dimensional grouping denoted as Pe (the Peclet number) is introduced as a
measure of the relative strengths of convection and diffusion. The Peclet number is
defined as:
.
/ x
u
D
F
Pe
?
?
?
== (2.10)
It is very important that the relationship between the magnitude of the Peclet number and
the directionality of influence, known as the transportiveness, is borne out in the
discretization scheme.
Although the central differencing scheme is accurate to 2
nd
order, one of its major
inadequacies is the inability to identify the flow direction. Figure 2.2 provides the basis
for a 1-D control-volume formulation for the transport of the scalar quantity ? using a
uniform staggered grid system (i.e. equal spacing in the horizontal direction between
consecutive grid nodes). The value of the property ? on the west face is influenced by
both ?
P
and ?
W
equally in the central differencing scheme. In a strongly convective flow
(Pe > 2) from west to east, the above treatment is unsuitable because the west cell face
should receive much greater influence from node W in comparison to node P. The
upwind differencing scheme takes into account the flow direction when determining the
value at a cell face. In summary, upon adoption of the upwind scheme the convected
value of ? at a cell face is taken to be equal to the value at the upstream node. For
instance, when the flow is in the positive direction, u
w
> 0 and u
e
> 0, the upwind scheme
requires:
PeWw
and ???? ==
. (2.11)
46
Alternately, when the flow is in the negative direction, u
w
< 0 and u
e
< 0, the upwind
scheme requires:
EePw
and ???? ==
. (2.12)
From the above discussion, the upwind differencing scheme is only accurate to the 1
st
order but accounts for transportiveness.
The hybrid differencing scheme of Spalding (1972) is based on a combination of
the central and upwind differencing schemes. The central differencing scheme is
employed for small Peclet numbers (Pe < 2) and the upwind scheme is employed for
large Peclet numbers (Pe ? 2).
The hybrid differencing scheme exploits the favorable properties of the upwind
and central differencing schemes. It switches between these two schemes depending on
the Peclet number. It satisfies the transportiveness requirements. It is fully conservative
and is unconditionally bounded since the main coefficients are always positive. It
produces physically realistic solutions, and has been widely used in various
computational fluid dynamics procedures. The disadvantage of hybrid scheme is that the
accuracy in terms of the Taylor's series truncation error is only first-order.
2.5 Quadratic Upwind Differencing Scheme: The QUICK Schemes
Since the hybrid scheme switches to the upwind scheme at high Peclet numbers,
the accuracy of the hybrid and upwind schemes is only first-order in terms of the Taylor's
series truncation error. The first-order accuracy makes them prone to numerical errors. In
order to obtain more accurate results, such errors should be minimized by employing
high-order discretization which involves more neighboring points.
47
48
The use of upwind quantities ensures that the schemes are very stable since it
obeys the "transportiveness" requirement. The central differencing scheme involves one
neighbor point on each side and does not take into account the flow direction. So it has
been proved that the central differencing scheme is unstable at high Pe numbers. More
accurate high-order schemes, which preserve upwinding for stability, are needed.
2.5.1 QUICK Scheme of Leonard (1979)
The Quadratic Upstream Interpolation for Convective Kinetics (QUICK) was
outlined by Leonard (1979). It uses a three-point upstream-weighted quadratic
interpolation to approximate the cell face value. Leonard (1979) proved that this scheme
has greater formal accuracy than the central difference scheme and still retains the basic
stable convective sensitivity property with upstream-weighted interpolation procedure.
Soon after Leonard's (1979) introduction of the QUICK scheme for the discretization of
the convective terms, a series of papers appeared in the literature. Leschziner (1980), Han
et al. (1981), Polland and Siu (1982) addressed the practical implementation and testing
of the scheme in flows more complex than those originally studied by Leonard (1979).
Upon reference to Figure 2.2, the values of ? at the grid points are identified with
capital letter subscripts (?
WW
, ?
W
, ?
P
, ?
E
and ?
EE
), whereas the values of ? at the cell
faces of the P control volume are denoted with lower case letter subscripts (?
w
and ?
e
).
The control volume surface values for ? are obtained by a quadratic fit via using the
values of ? at three connective nodes (the two nodes located on either side of the surface
of interest, plus the adjacent node on the upstream side).
By reference to Figure 2.2, when u
w
> 0, a quadratic fit through WW, W and P is
used to evaluate ?
w
:
WWPWw
????
8
1
8
3
8
6
?+=
. (2.13)
When u
w
< 0, a quadratic fit through W, P and E is used to evaluate ?
w
:
EWPw
????
8
1
8
3
8
6
?+=
. (2.14)
When u
e
> 0, a quadratic fit through W, P and E is used to evaluate ?
e
:
WEPe
????
8
1
8
3
8
6
?+=
. (2.15)
When u
e
< 0, a quadratic fit through P, E and EE is used to evaluate ?
e
:
EEPEe
????
8
1
8
3
8
6
?+=
. (2.16)
Upon substitution of the above approximation of the values at the cell faces into
the 1-D convection-diffusion discretization equation (assuming u
e
> 0 and u
w
> 0):
()()????=
?
?
?
?
?
?
?+??
?
?
?
?
?
?+
WPwPEe
WWPWwWEPe
DD
FF
????
??????
8
1
8
3
8
6
8
1
8
3
8
6
(2.17)
The above discretized equation is then written to fit into the standard form:
uEEWWWWWWPP
Saaaa +++= ????
(2.18)
with
( )
PweEWWWP
SFFaaaa ??+++=
(2.19)
49
and
WW
a
W
a
E
a
p
S
u
S
w
F
8
1
?
wew
FFD
8
6
8
1
++
ee
FD
8
3
?
0 0
Patankar (1980) pointed out that all coefficients ( and neighbor coefficients
) must always be positive in order to get stable and bounded solutions. According to
the above analysis, the main coefficients ( and ) are not guaranteed to be
positive and the coefficients ( and ) are negative. For example, if u
P
a
nb
a
E
a
W
a
EE
a
WW
a
w
> 0 and
u
e
> 0, the east coefficient E
a
becomes negative when the cell Peclet number Pe = F
e
/
D
e
> 8/3. Similarly, the west coefficient W
a
can become negative too if the flow is in
the other direction. This gives rise to stability problems and unbounded solutions under
certain flow conditions if Leonard?s QUICK scheme is used.
Another problem using Leonard?s QUICK scheme is the fact that the discretized
equation not only involves immediate-neighbor nodes (W and E) but also nodes farther
away (WW and EE). This gives rise to a penta-diagonal system that has five non-zero
coefficients for a 1-D convection-diffusion problem. So, it is not possible to solve the
final algebraic equations using the tri-diagonal matrix algorithm (TDMA) directly in this
situation. It is also obvious that more memory is needed to store a penta-diagonal system
compared to a tri-diagonal system.
50
2.5.2 A Stable and Fast-Converging Variation of the QUICK Scheme: The QUICK
Scheme of Hayase et al. (1992)
As shown above, Leonard?s QUICK scheme can be unstable due to the negative
main coefficients under certain conditions. Han et al. (1981), Pollard and Siu (1982) and
Hayase et al. (1992) have re-formulated the QUICK scheme in different ways in order to
alleviate stability problems. They lumped the troublesome negative coefficients into the
source term to retain positive main coefficient.
The Hayase et al.'s (1992) QUICK scheme for a uniform grid (Figure 2.2) can be
summarized as follows:
()
()
()
().023
8
1
023
8
1
023
8
1
023
8
1
??+=
>??+=
eEEEPEe
wEPWPw
eWPEPe
wWWWPWw
ufor
ufor
ufor
ufor
?????
?????
?????
?????
(2.20)
This form of the QUICK scheme formulation has a very simple structure, i.e. the value of
? on the cell face is given as the sum of an upwind estimation with correction terms. The
correction terms can conveniently be lumped into the source term.
Hayase et al. (1992) proved that this QUICK scheme satisfies all the following five
constraints:
(1) Consistency of surface flux calculations at the control-volume faces,
(2) All coefficients are positive in the discretization equation,
(3) Negative-slope linearization of the source term,
51
52
(4) Node coefficient equals the sum of neighboring node coefficients, and,
(5) Constraints on the summations of coefficients.
Extensive tests were performed by Hayase et al. (1992) and the results showed
that this QUICK scheme produces physically realistic and stable numerical solutions of
the finite-volume approximated transport equations. It also exhibits more robustness and
converges considerably faster than any other previous formulation.
Another attractive characteristic of this scheme is that the main coefficients are
the same as those using the upwind scheme. No farther away nodes are involved to
calculate the main coefficients, so the TDMA can be used directly to solve the final
algebraic matrix. It is very easy and straightforward to implement this scheme in a CFD
solver that uses TDMA. It also needs the same size memory to accommodate all the main
coefficients as that using the upwind scheme.
2.5.3 High-Order Approximation for Boundary Conditions
It is well known that the solution of partial differential equations not only depends
on the equation itself but also on the boundary conditions. Similarly, the accuracy of the
numerical solution for the transport equations is not only affected by the scheme used for
the points within the calculation domain, but is also dependent on the approximation of
the boundary conditions. Hayase et al. (1992) reported that high-order approximations for
the boundary conditions are needed in order to obtain highly accurate results. All the
approximations are 3
rd
order accurate. With these approximations for the boundaries and
the QUICK scheme for the rest of the computational domain, the entire computational
domain is of the 3
rd
order accuracy to approximate the convective terms.
2.6 Benchmarking of the Code
In order to investigate the accuracy of the control-volume-based code that is used
for solving the problems in the next two chapters, three related benchmark problems are
solved. These are:
1. 2-D duct flow,
2. 2-D backward?facing step flow,
3. 2-D forward?facing step flow.
These problems are very well-established benchmark problems that involve inlet and
outlet ports. Therefore, benchmarking of the code for these problems is directly relevant
to the topic of this thesis.
2.6.1 2-Dimensional Duct Flow
Fluid flow within two parallel plates is solved using both the HYBRID and
QUICK differencing schemes. The distance between the two plates (H) is equal to 0.05
m and the length of plates in the flow direction is 0.6 m (12 H). The hydraulic diameter
is defined as 2 ? H. The Reynolds number that is based on the hydraulic diameter and
the inlet velocity is defined as follows:
.
2
Re
?
Hu
in
=
(2.21)
For all the cases in this section, the Reynolds number is equal to 100. Three grid systems
(50?40, 80?60 and 100?80) are employed for each differencing scheme
53
54
2.6.1.1 HYBRID Scheme
The grid is non-uniform with an expansion ratio of 1.1 in both directions. In the
Y-direction, the grid is denser near the impermeable plates and in the X-direction, the
grid expands along the X-axis toward the outlet. The generated grid (100?80) is shown
in Figure 2.3. The streamlines and velocity vectors are presented in Figures 2.4 and 2.5,
respectively. Due to the simplicity of the geometry and low Reynolds number in the
duct, the streamlines are generally parallel to the duct walls except near the entrance
region. No recirculation zones are observed in the duct and the flow reaches its fully-
developed state after passing the entrance region. The velocity profile is parabolic in the
fully-developed region whereas next-to-wall fluid velocity overshoots near the entrance
of the duct. The overshoot phenomena will be discussed in detail in the following
section.
2.6.1.2 QUICK Scheme
A uniformly-distributed grid system is used in combination with the QUICK
scheme. The QUICK scheme gives better results at the outlet, compared to the HYBRID
scheme, but within the developing zone near the inlet, the HYBRID scheme gives results
that are more accurate. The reason for this is that with the HYBRID scheme, the grid is
very dense in the entrance zone and there are more grids near the inlet plane in
comparison to the QUICK scheme.
2.6.1.3 Results for the 2-D Channel Studies
The velocity profiles at different axial stations along the duct using the HYBRID
scheme with a 100?80 grid are presented in Figure 2.6. Next-to-wall velocity overshoots
are clearly observed within the developing zone ( 45.0/ ?Hx ). The QUICK scheme
(with uniform grid) could not detect the overshoot because overshoots are generated very
close the inlet section. The overshoot profile progresses toward the mid-plane (y / H =
0.5) as the flow moves through the duct. At the outlet plane ( ) the velocity
profile is parabolic. Darbandi and Schneider (1998) studied the laminar flow between
two parallel plates for a wide range of Reynolds numbers using different mesh densities.
Their results also showed an overshoot at the entrance region.
12/ =Hx
To determine the magnitude of the overshoot at every axial station with respect to
the centerline velocity at the same axial location, the strength of maxima is defined as:
,
max
cl
cl
U
UU
Z
?
= (2.22)
where U
max
and U
cl
are the maximum velocity and centerline velocity at a given axial
station, respectively. This quantity is shown in Figure 2.7. The X axis is
nondimensionalized by the entrance length (L
e
). The distance downstream from the
entrance to the location at which the fully-developed state is achieved, is called the
entrance length (L
e
). The entrance length in laminar channel flow is observed to be a
function of the Reynolds number as (White, 1999):
.Re06.0=
H
L
e
(2.23)
For this problem the entrance length is 0.151 m for the Reynolds number of 100. The
nondimensionalized L
e
based on distance between two plates (H) is 3.02. The initial rise
55
of the overshoot parameter reached 14.3 percent at 041.0?
e
L
x
followed by its monotonic
decline (Figure 2.7).
The axial development of the centerline axial velocities using both the QUICK
and HYBRID schemes are shown in Figure 2.8 using the respective 100?80 grids. As it
is expected from the analytical solution, the ratio of the centerline velocity to the inlet
velocity at the fully-developed state is 1.5. Both QUICK and HYBRID schemes show
good agreement with the analytical solution near the outlet.
Under the fully-developed condition, the shear stress value at the end of the duct
can be determined using an analytical formula:
,
2
1
?
?
?
?
?
?
?
?
?
?
?
?
?
=
H
y
dx
dp
H
xy
? (2.24)
with
LH
Q
dx
dp
3
12??
= , (2.25)
where ? is the fluid viscosity, H is the distance between the two plates, L is the duct
length and Q is the volumetric flow rate. For this problem, the wall friction coefficient
under the fully-developed condition is 120.0
2
==
in
w
f
U
C
?
?
. The computed results have
been compared to this value in Figure 2.9. The Y axis values are nondimensionalized
wall shear stresses that are normalized by the analytical solution value (0.12). The values
on the X axis (m) are the total number of nodes that is equal to NI?NJ. It is obvious that
the QUICK scheme predicts the wall shear stress much better than the HYBRID scheme.
56
The HYBRID scheme loses accuracy as the number of nodes increases, but the QUICK
scheme estimation improves with finer mesh refinement.
2.6.2 The Backward-Facing Step Flow
The problem of flow over a 2-D backward-facing step is a very well-established
benchmark problem. There are some available experimental and numerical work on this
problem. Freitas (1995) reported the results of flow simulation over a backward-facing
step that were obtained using various commercial codes. The computational results were
compared to the experimental data of Armaly et al. (1983). It should be noted that Abu-
Mulaweh (2003) reviewed the research on laminar mixed convection flow over backward
and forward-facing steps. For this study, the experimental results of Armaly et al. (1983)
are used for comparison to computational results of our finite-volume based code.
A 2-D channel problem with a sudden asymmetric expansion is solved using the
QUICK scheme. This expansion is called a backward-facing step when the upstream
height (h) is smaller than downstream height. The height of the step is S = 4.9 mm and
the height of the inlet is 2.5=h mm. The expansion ratio is defined
as . The problem is solved in three steps, called Cases A, B and C.
Schematic diagrams of the three geometries are shown in Figure 2.10. For Case A that
has the simplest geometry, the expansion (step) is located at the beginning of the channel.
Both uniform and parabolic velocity profiles as the inlet boundary condition are used. At
the outlet, the fully-developed boundary condition is applied. The Reynolds number is
based on the hydraulic diameter that is equal to . The length of the channel
downstream of the step is chosen to be
94.1/)( =+= hShER
h2
125=Lmm S?? 25 . This length is long enough
57
so that the code can capture all possible phenomena such as the recirculating vortices and
the flow redevelopment that occurs in the channel. For Cases B and C, an extension is
added to the upstream end of the step. For Case B, the length of the extension is 100
and for Case C, its length is . In other words, the sudden expansion for these
cases is placed at stations 100 and 200 downstream of the inlet plane. For
Cases B and C, uniform velocity profiles are used as the inlet boundary condition, thus
allowing the flow to develop in the axial direction and assume a ?near-fully-developed?
profile at the step. The fully-developed outlet boundary condition is applied to the Cases
B and C similar to Case A.
mm
200 mm
mm mm
2.6.2.1 Results for the 2-D Back-Step Flow
A uniform grid distribution is used for all cases. Case A has a 100?50 grid. For
Cases B and C, additional 80?50 and 160?50 grid systems are added for their upstream
extensions, respectively. Note that in all cases there are 100 nodes distributed in X
direction between the step and the outlet of the channel.
The results are presented for two Reynolds numbers 200 and 450. When Re=200,
one vortex forms right after the step. The separating streamline of this vortex originates
at the step and reattaches on the bottom wall of the channel. For Re=450, there is another
vortex occurring underneath the top wall. The locations of the eye of the vortices are
reported in Table 2.1. The coordinate system used is Cartesian and the origin is placed at
the bottom corner of the step. The values are nondimensionalized with the step height
(S). Subscript 1 is used for the vortex neighboring the step and subscript 2 is used for the
58
top vortex. The vortices are clearly observed in Figures 2.11 to 2.18 for which
streamlines are sketched for all cases.
Table 2.1 The location of the eye of the vortices for the back-step problem
Case B
Grid: 180?50
Case A
Grid:
100?50
Uniform
inlet flow
Case A
Grid:
100?50
Parabolic
inlet flow
ER = 1.92
ER =1.96
Case C
Grid:
260?50
Re = 200
Sx
e
/
1
1.04 1.82 1.68 1.80 1.68
Sy
e
/
1
0.49 0.58 0.58 0.58 0.58
Re = 450
Sx
e
/
1
1.82 3.90 3.35 3.87 3.35
Sy
e
/
1
0.49 0.58 0.58 0.58 0.58
Sx
e
/
2
n/a 9.89 9.29 9.80 9.52
Sy
e
/
2
n/a 1.91 1.95 1.91
1.95
x, y : Horizontal and vertical distances from the bottom corner of the step, respectively
Subscripts:
e: Eye of vortex
S: Step height
1: Vortex created close to the bottom corner of the step
2: Vortex created under the top wall of the channel
59
Table 2.2 Comparison of the locations of computed and measured separation and
reattachment points for the back-step problem
Case B
Grid: 180?50
Experimental
(Armaly et al.
1983)
ER = 1.94
Case A
Grid: 100?50
Uniform inlet
flow
Case A
Grid: 100?50
Parabolic inlet
flow
ER = 1.92
ER =1.96
Case C
Grid:
260?50
Re = 200
5.14 4.86 5.02 4.87 Sx
r
/
1
5.0 2.86
?=0.0 ?=0.051 ?=0.047 ?=0.051
Re = 450
Sx
r
/
1
9.5 4.78 8.97 8.69 8.87 8.70
Sx
s
/
2
7.6 n/a 7.76 7.86 7.64 7.85
Sx
r
/
2
11.3 n/a 11.84 10.89 11.82 10.91
4.09 3.03 4.18
3.06
3.7 n/a
?=0.0 ?=0.0168 ?=0.0162 ?=0.0160
Sxx
sr
/)(
22
?
x, y : Horizontal and vertical distances from the bottom corner of the step, respectively
Subscripts:
s, r: Separation and reattachment points, respectively
S: Step height
1: Vortex created close to the bottom corner of the step
2: Vortex created under the top wall of the channel
60
The separation and reattachment points for the vortices are compared with the
experimental results reported by Armaly et al. (1983) in Table 2.2. The variables and
stand for the x-coordinates of the separation point and reattachment point, respectively.
Except Case A with uniform inlet flow, other cases show an acceptable agreement with
the experimental observations of Armaly et al. (1983). It is therefore concluded that the
inlet velocity variation at the step is very critical to the successful simulation of this flow.
Above the step at x = 0, the velocity variation is compared to the parabolic velocity
distribution. To determine the difference between the velocity profiles, the variable ? is
defined as:
s
x
r
x
,
1
2
n
u
uun
j
j
p
j
pj
?
=
?
?
?
?
?
?
?
? ?
=?
(2.6)
where is the total number of nodes in vertical direction (j) from the top corner of the
step to the top wall of the channel, is the velocity at point j, and is the velocity at
the same point if it was a perfect parabolic velocity profile. The computed values of ? are
summarized in Table 2.2.
n
j
u
j
p
u
The expansion ratio ( ER ) of the experimental work is 1.94 that is a bit different
from that is used for the present study. Because of the restriction of the
QUICK scheme, the grid had to be uniform, so it was not possible to duplicate the
desirable
92.1=ER
ER . The closest possible ER values to 1.94 were 1.92 and 1.96. An expansion
ratio of 1.92 was chosen but Case B was run with both 1.92 and 1.96.
61
2.6.3 The Forward-Facing Step Flow
The flow over a 2-dimensional forward-facing step is another interesting problem
for benchmarking of the CFD codes. Mei and Plotkin (1986) studied the flow field in a
forward-facing step using the vorticity-stream function form of the Navier-Stokes
equations. They also developed an analytical solution of the problem by using the
conformal mapping technique. The flow was assumed inviscid in the analytical part.
Other notable studies were carried out by St?er et al. (1999) and Abu-Mulaweh and Chen
(1993).
In the present study, two geometries are considered (Cases A and B) and the
relevant geometries are shown in Figure 2.19. Case A is a forward step that has no
extension downstream of the step. In Case B, an extension has been added downstream
of the step. The height of the step is 9.4=S mm and the height of the inlet is
. The expansion ratio is defined as1.10=h mm 94.1)/( =?= ShhER . For Case A, two
different inlet velocity profiles are applied as the inlet boundary condition, namely
uniform and parabolic flows, whereas a uniform flow is used as the inlet boundary
condition for Case B. A fully-developed outlet boundary condition is used for all cases.
Similar to the back-step flow (section 2.6.2), the length of the upstream section is chosen
to be and the length of the extension for Case B is 100 mm . The
Reynolds number is based on the hydraulic diameter that is equal to . The QUICK
scheme is used for all the cases to solve the flow field.
125=Lmm ?? 25 S
h2
62
2.6.3.1 Results for the 2-D Forward-Step Flow
A uniform grid distribution is used for all cases. Case A has a 100?50 grid. For
Cases B an additional 80?50 grid is added as its downstream extension. In both cases
there are 100 nodes distributed in the X direction from the inlet port to the step (wider
portion of the channel).
The problem is solved for three Reynolds numbers of 200, 450 and 1000. For all
the Reynolds numbers, a recirculation zone forms at the bottom corner of the step. As the
Reynolds number increases, the size and strength of this vortex increases. Flow
separation occurs on the bottom wall. The flow will then reattach on the vertical wall of
the step. The streamlines are shown in Figures 2.20 to 2.28 for different cases. For Case
B, another vortex forms right at the beginning of the extension where the flow has just
passed the step. The flow separates at the edge of the step and reattaches to the bottom
wall of the extension downstream. This vortex is longer and thinner than the vortex that
appeared at the bottom corner of the step. This vortex can be viewed in Figures 2.22,
2.25 and 2.28. The locations of the eye of the vortices are summarized in Table 2.3. The
coordinate system is Cartesian and the origin is placed at the bottom corner of the step.
The values are nondimensionalized with the step height ( ). The subscript 1 is used for
the vortex neighboring the step and the subscript 2 is used for the vortex on the step.
S
In Table 2.4, the separation and reattachment points for the vortices are reported.
The variables and stand for the x-coordinates of the separation point and
reattachment point, respectively. Comparing the values in the first two columns of
Tables 2.3 and 2.4, it is concluded that there is not a significant difference between the
results of different boundary conditions that are applied to Case A. This is because the
s
x
r
x
63
64
channel is long enough that the velocity profile at the inlet does not affect the
characteristics of the vortex. For a constant Reynolds number, the x-coordinate of the
eye of the bottom vortex is almost the same for all cases, whereas the y-coordinate of the
eye in Case B (3
rd
column, Table 2.3) is different than Case A (1
st
and 2
nd
columns, Table
2.3). This difference is more recognizable for the low Reynolds number flows (Re = 200
and 450). In addition, the coordinates of the separation and reattachment points of the
bottom vortex are almost the same for all cases for a constant Reynolds number (Table
2.4).
2.7 Closure
In order to obtain high accuracy computational results for elliptic problems with
inlet and outlet ports, a stable fast-converging variation of the QUICK differencing
scheme within the finite-volume framework is discussed. The software developed for
this study is verified by comparing the results with several well-known benchmark
problems in the literature. Three classic benchmark problems are solved. These were the
two-dimensional channel, backward- and forward-facing step flows. The results are
successfully compared to the available analytical and experimental information.
Table 2.3 The location of the eye of the vortices for the forward-step problem
Case A
Grid:
100?50
Uniform
inlet flow
Case A
Grid:
100?50
Parabolic
inlet flow
Case B
Grid:
180?50
Re = 200
sx
e
/
1
-0.264 -0.264 -0.260
sy
e
/
1
0.109 0.109 0.065
Re = 450
sx
e
/
1
-0.264 -0.264 -0.260
sy
e
/
1
0.109 0.109 0.152
sx
e
/
2
- - 0.262
sy
e
/
2
- - 1.02
Re = 1000
sx
e
/
1
-0.264 -0.264 -0.260
sy
e
/
1
0.196 0.196 0.196
sx
e
/
2
- - 0.26
sy
e
/
2
- - 1.02
x, y : Horizontal and vertical distances from the bottom corner of the step, respectively
Subscripts:
e: Eye of vortex
S: The step height
1: The vortex created close to the bottom corner of the step
2: The vortex created on the bottom wall of the extension
65
Table 2.4 The location of separation and reattachment points for the forward-step
problem
Case A
Grid:
100?50
Uniform
inlet flow
Case A
Grid:
100?50
Parabolic
inlet flow
Case B
Grid:
180?50
Re = 200
sx
s
/
1
-0.442 -0.442 -0.313
Re = 450
sx
s
/
1
-0.569 -0.57 -0.566
sx
r
/
2
- - 0.602
Re = 1000
sx
s
/
1
-0.845 -0.846 -0.860
sx
r
/
2
- - 1.265
x, y : Horizontal and vertical distance from the bottom corner of the step, respectively
Subscripts:
s, r: Separation and reattachment points, respectively
S: Step height
1: Vortex created close to the bottom corner of the step
2: Vortex created on the bottom wall of the extension
66
n
N
V Cell
s
S
P Cell
W
e w
P
E
U Cell
J+1
j+1
J
j
J-1
j-1
J-2
I-2 I-1 I I+1 i-1 i i+1
Figure 2.1 2-D staggered grid system
67
X
?
e
?
w
?
EE
?
E
?
P
?
W
u
w
?
WW
u
e
WW W P E EE
Figure 2.2 1-D uniform staggered grid showing the nodes involved in the evaluation of ?
on cell faces of a control-volume
68
Figure 2.3 Non-uniform grid 100?80 (HYBRID)
Figure 2.4 Streamlines
Figure 2.5 Velocity Vectors
69
u/U
in
y/
H
00.511.5
0
0.1
0.2
0.3
0.4
0.5
x/H = 0
x/H = 0.0045
x/H = 0.02
x/H = 0.1
x/H = 0.25
x/H = 0.45
x/H = 1
x/H = 12
Figure 2.6 Velocity profiles in the developing zone of a channel (HYBRID, Re = 100)
70
x/L
e
Z[
%
]
0 0.1 0.2 0.3 0.4 0.5 0.6
0
3
6
9
12
15
Figure 2.7 Growth and decay of the overshoot parameter along the X- axis (HYBRID
scheme, 100
?
80 grid)
71
Hybrid Scheme
x/H
u
cl
/U
in
024681012
1
1.1
1.2
1.3
1.4
1.5
1.6
QuickScheme
Figure 2.8 Developing of the centerline velocity along the X-axis (100
?
80 grid)
72
Wall Shear Stress
Figure 2.9: The normalized wall shear stress versus mesh refinement
73
Figure 2.10 Schematic diagrams of the backward-facing flow geometries
Figure 2.11 Streamlines for Re = 200 (Case A) with uniform inlet flow
Figure 2.12 Streamlines for Re = 200 (Case A) with parabolic inlet flow
Figure 2.13 Streamlines for Re = 200 (Case B)
Figure 2.14 Streamlines for Re = 200 (Case C)
Figure 2.15 Streamlines for Re = 450 (Case A) with uniform inlet flow
Figure 2.16 Streamlines for Re = 450 (Case A) with parabolic inlet flow
Figure 2.17 Streamlines for Re = 450 (Case B)
Figure 2.18 Streamlines for Re = 450 (Case C)
74
Figure 2.19 Schematic diagrams of the forward-facing flow geometries
Figure 2.20 Streamlines for Re = 200 (Case A) with uniform inlet flow
Figure 2.21 Streamlines for Re = 200 (Case A) with parabolic inlet flow
Figure 2.22 Streamlines for Re = 200 (Case B)
75
Figure 2.23 Streamlines for Re = 450 (Case A) with uniform inlet flow
Figure 2.24 Streamlines for Re = 450 (Case A) with parabolic inlet flow
Figure 2.25 Streamlines for Re = 450 (Case B)
Figure 2.26 Streamlines for Re = 1000 (Case A) with uniform inlet flow
Figure 2.27 Streamlines for Re = 1000 (Case A) with parabolic inlet flow
Figure 2.28 Streamlines for Re = 1000 (Case B)
76
77
CHAPTER 3 STEADY-STATE FORCED CONVECTION IN A SQUARE
CAVITY WITH INLET AND OUTLET PORTS
Forced convection within a lid-driven sealed cavity and buoyancy-driven
convection within differentially heated cavities have provided excellent opportunities for
researchers to test and benchmark computational tools for solution of complex flows with
recirculation. In addition to their geometrical simplicity, both of these flows offer
simultaneous existence of diverse flow regimes involving boundary layers and multiple
separated regions. To this end, great attention has been focused on these flows and
variations of them by Dr. Khodadadi?s research group at Auburn University:
1. Steady fluid flow and heat transfer in a lid-driven square cavity due to a thin fin
(Shi and Khodadadi, 2002),
2. Transient evolution to periodic fluid flow and heat transfer in a lid-driven
square cavity due to an oscillating thin fin (Shi and Khodadadi, 2004, and 2005),
3. Laminar natural convection heat transfer in a differentially-heated square cavity
due to a thin fin (Shi and Khodadadi, 2003).
A logical extension of the earlier work is to concentrate on more complicated flow
situations in cavities upon introduction of multiple inlet and outlet ports. Such complex
systems offer a wide range of interesting flow regimes and have a great number of
engineering applications. For instance, mixed convection studies within cavities that
78
contain heat sources are of great interest in electronic cooling, ventilation of buildings
and design of solar collectors. The relevant literature was reviewed in Section 1.1.3.
Simulations of unsteady mixed convection in thermal storage applications have also
received a great deal of attention for many decades (Section 1.1.1). Other isolated
applications related to mass transfer in cavities (Occhialini and Higdon, 1992), relaxation
tanks employed in petroleum and transformer industries (Jayaram and Thangaraj, 1995),
transient removal of contaminated fluid from cavities (Fang et al., 1999) and food
processing (Verboven et al., 2003) are also noteworthy. Design of microfabricated fluid
filters (Leung Ki et al., 1998) and piezoelectric liquid-compatible microvalves (Yang et
al., 2004) can also benefit from the knowledge that is gained from the present study.
Given the variety of applications sharing the cavity geometry with inlet and outlet
ports, the present parametric study was undertaken in order to gain further knowledge of
the complicated steady-state flow and thermal fields in such a configuration.
3.1 Problem Formulation
The schematic drawing of the system and the coordinates are shown in Figure 3.1.
The height and width of the cavity are denoted by H and L, respectively. In this study,
only a square cavity is considered (H = L). The depth of the enclosure perpendicular to
the plane of the diagram is assumed to be long. Hence, the problem can be considered to
be two-dimensional. The fluid enters the cavity through an inlet port of width w
i
that is
located on the left vertical wall extending from y = (H - w
i
) to y = H. An exit port of
width w
o
can be present on any of the four walls. In this study the widths of the inlet and
outlet ports are identical (w
i
= w
o
= w). A special coordinate system (s) along the walls is
adopted with its origin at x = 0 and y = H, as identified by the dashed lines in Figure 3.1.
The s-coordinates of the symmetry planes of the inlet and outlet ports are 0.5w
i
and s
o
,
respectively. The temperature of the fluid entering the cavity is T
in
, whereas the four
solid walls of the cavity are maintained at temperature T
w
with T
w
> T
in
.
3.1.1 Governing Equations
The fluid within the enclosure is an incompressible fluid and the fluid properties
are constant. The flow within the enclosure is assumed 2-D, steady and laminar. The
gravity effect is negligible. The governing equations of continuity, momentum and
thermal energy are:
,0=
?
?
+
?
?
y
v
x
u
(3.1)
,
1
2
2
2
2
?
?
?
?
?
?
?
?
?
?
+
?
?
+
?
?
?=
?
?
+
?
?
y
u
x
u
x
p
y
u
v
x
u
u ?
?
(3.2)
,
1
2
2
2
2
?
?
?
?
?
?
?
?
?
?
+
?
?
+
?
?
?=
?
?
+
?
?
y
v
x
v
y
p
y
v
v
x
v
u ?
?
(3.3)
,
2
2
2
2
?
?
?
?
?
?
?
?
?
?
+
?
?
=
?
?
+
?
?
y
T
x
T
y
T
v
x
T
u ?
(3.4)
in which u, v, p and T are the velocity components in the x and y directions, pressure and
temperature, respectively. Quantities ?, ? and ? are the density, kinematic viscosity and
thermal diffusivity of the fluid, respectively. In the energy equation, the viscous
dissipation terms are neglected.
79
At the inlet port a constant velocity and temperature are imposed. The no-slip
condition is imposed on the solid surfaces. At the outlet port, the flow and temperature
are assumed fully-developed. So the boundary conditions can be expressed as follows:
At the inlet (x = 0, y = (H - w
i
) to H): u = u
in
, v = 0, T = T
in
, (3.5a)
On the four solid walls: u = v = 0, T = T
w
. (3.5b)
At the outlet port (s = s
o
? 0.5w
o
to s
o
+ 0.5w
o
):
For vertical outlet ports: 0,0 =
?
?
=
?
?
x
T
x
u
. (3.5c)
For horizontal outlet ports: 0,0 =
?
?
=
?
?
y
T
y
v
. (3.5d)
3.1.2 Dimensionless Form of the Governing Equations
Dimensionless form of the governing equations can be obtained via introducing
dimensionless variables. The lengths can be scaled by the cavity length H and the
velocities can be scaled by the inlet fluid velocity, u
in
. As for the temperature, the two
extreme values (T
in
and T
w
) are used. The dimensionless variables are then defined as:
.
2
inw
in
in
inin
TT
TT
u
p
P
u
v
V
u
u
U
H
y
Y
H
x
X
?
?
==
====
?
?
(3.6)
Based on the dimensionless variables above, the dimensionless equations for the
conservation of mass, momentum and thermal energy are:
,0=
?
?
+
?
?
Y
V
X
U
(3.7)
80
,
Re
12
2
2
2
2
?
?
?
?
?
?
?
?
?
?
+
?
?
+
?
?
?=
?
?
+
?
?
Y
U
X
U
H
w
X
P
Y
U
V
X
U
U
(3.8)
,
Re
12
2
2
2
2
?
?
?
?
?
?
?
?
?
?
+
?
?
+
?
?
?=
?
?
+
?
?
Y
V
X
V
H
w
Y
P
Y
V
V
X
V
U
(3.9)
.
RePr
12
2
2
2
2
?
?
?
?
?
?
?
?
?
?
+
?
?
=
?
?
+
?
?
YXH
w
Y
V
X
U
????
(3.10)
The Reynolds number is defined as Re = u
in
(2w) / ? and the Prandtl number is Pr
= ? / ?. The Reynolds number signifies the ratio of inertia to viscous forces, whereas the
Prandtl number indicates the ratio of the momentum to thermal diffusivities.
Denoting W = w / H, the dimensionless form of the boundary conditions are:
At the inlet (X = 0, Y = (1 - W) to 1): U = 1, V = 0, ? = 0, (3.11)
On the four solid walls: U = V = 0, ? = 1. (3.12)
At the outlet port (S = s / H = S
o
? 0.5W to S
o
+ 0.5W), the velocity and temperature fields
were taken to be fully-developed. From the foregoing formulation, it is clear that the
dimensionless parameters governing this problem are Re, Pr, S
o
and W. In this study, the
Prandtl number of the fluid is fixed to 5. The Reynolds numbers considered are 10, 40,
100 and 500.
3.2 Computational Details
The governing equations were iteratively solved by the finite-volume-method
using Patankar's SIMPLE algorithm (Patankar, 1980). A two-dimensional uniformly-
spaced staggered grid system was used. Hayase et al.'s (1992) QUICK scheme was
81
82
utilized for the convective terms, whereas the central difference scheme was used for the
diffusive terms. In order to keep consistent accuracy over the entire computational
domain, a third-order-accurate boundary condition treatment suggested by Hayase et al.
(1992) was adopted.
3.2.1 Grid Independence Study and Benchmarking of the Code
In order to determine the proper grid size for this study, a grid independence test
was conducted with Re = 1000, S
o
= 2.1 and W = 0.2. Seven uniformly-spaced grid
densities were used for the grid independence study. These grid densities were 50
x
? 50
y
,
75
x
? 75
y
, 100
x
? 100
y
and 125
x
? 125
y
, 150
x
? 150
y
, 175
x
? 175
y
and 200
x
? 200
y
. The
minimum value of the stream function of the primary vortex (?
min
) is commonly used as
a sensitivity measure of the accuracy of the solution and was selected as the monitoring
variable for the grid independence study. Figure 3.2 shows the dependence of quantity
?
min
on the grid size. Considering both the accuracy and the computational time, the
following calculations were all performed with a 160
x
? 160
y
grid system that
accommodated resolving all the three widths of the outlet port on the uniform grid
system. The code was recently tested and verified extensively via comparing the results
with established forced and buoyancy-driven convection benchmark problems (Shi and
Khodadadi, 2002 and 2003). Benchmarking results for two-dimensional channel flow,
backward-facing and forward-facing step problems were given in section 2.6 of this
thesis.
83
3.2.2 Parameters for Numerical Simulations
Tolerance of the normalized residuals upon convergence is set to 10
-6
for most
cases, whereas for some high Re cases a less stringent condition (10
-5
) was imposed. For
Re numbers other than 500, the under-relaxation parameters for u, v, and T are all set to
0.6, whereas under-relaxation parameter for pressure correction is set to 0.3. For Re =
500, these were set to 0.6, 0.6, 0.2 and 0.3.
3.3 Results and Discussion
In order to understand the flow field and heat transfer characteristics of this
problem, a total of 108 cases were considered. This involved studying the effect of
placing an outlet port at 9 different positions for a fixed position of the inlet port. The
Reynolds numbers are 10, 40, 100 and 500. The widths of both ports are equal to 5%,
15% and 25% of the width of the cavity. The typical number of iterations needed to
converge was 11,000 ~ 60,000. Convergence difficulties were encountered for seven
cases with Re = 500 and W = 0.05 when the outlet port was not located next to the lower
right corner. Similar difficulties were also repeated by running these cases on the
commercial CFD package FLUENT. For these cases, hybrid differencing scheme was
utilized. All the calculations were performed on a Cray SV1 of the Alabama
Supercomputer Center (Huntsville, Alabama).
3.3.1 Flow Fields in the Cavity
The streamlines corresponding to nine cases with Re = 500 and W = 0.25 are
shown in Figure 3.3 for different positions of the outlet port. Similar streamline patterns
84
for nine cases with Re = 40 and W = 0.15 are shown in Figure 3.4. For each of the
individual cavities shown, two arrows are drawn to help the reader clearly identify the
inlet and outlet ports. Nine different positions of the outlet port are shown corresponding
to variation of S
o
from 0.875 to 3.5 in Figure 3.3 and 0.925 to 3.5 in Figure 3.4. It is
observed that the streamlines representing the throughflow of the incoming fluid
(streamline values of zero and inlet mass flowrate) do not coincide with the two ends of
the outlet ports. The streamlines for the top row correspond to cases where the outlet port
is positioned on the top wall or closest to it. The fresh fluid entering the cavity travels the
shortest distance possible before leaving the cavity. For these three cases of Figure 3.3,
the presence of a primary CW (clockwise) rotating vortex that occupies 75 to 88 percent
of the cavity is clearly observed. Two CCW (counter clockwise) rotating vortices that
are located at the bottom two corners - but are much smaller that the primary vortex ? are
the other permanent features of these three cases. Considering similar three positions of
the outlet ports, the flow fields appeared very similar for the other Re numbers studied
(e.g. the top row of Figure 3.4) except for Re = 10, with the strength the corner vortices
diminishing as Re decreases. For the Re = 10 cases corresponding to the same outlet port
positions, the cavity was generally dominated by two large CW rotating vortices at the
bottom two corners regardless of the width size. For an outlet port located at S
o
= 3.5
with Re = 10, these two vortices merged to form a single CW vortex occupying the
bottom of the cavity and a third small CCW vortex was observed at the top right hand
corner. Going through the middle row of Figure 3.3, one can observe the changes in the
flow field as the position of the outlet port is lowered along the right side wall. The area
occupied by the CW rotating primary vortex diminishes as the outlet port is lowered. A
85
CCW rotating vortex that is created next to the top right corner covers a larger portion of
the cavity as the outlet port is lowered. For the same three positions of the outlet ports,
the features of the flow fields were retained for the other Re numbers studied (such as the
middle row of 3.4). For the Re = 10 cases corresponding to these three outlet port
positions, flow is generally void of vortices except for a small CCW rotating one at the
top right corner and a CW rotating one at the bottom left corner. The flow fields shown
in the bottom row of Figure 3.3 are obtained as the outlet port is moved along the bottom
wall to the left corner, with the bottom right corner case showing a situation where the
two inlet and outlet ports are present on the left wall. As the outlet port is moved to the
left, the primary CW vortex is found to shrink in size and occupies the space adjacent to
the left wall. The CCW rotating vortex at the top right corner gains prominence as the
outlet port is moved to the left side occupying the right half of the cavity in the most
extreme case. An interesting feature of this CCW rotating vortex is its double-extrema
property when the outlet port is located at S
o
= 1.5. Similar trends were observed for
other cases when the Reynolds number was lowered.
The effects of the position and width of the outlet ports on the flow field for Re =
100 are illustrated in Figure 3.5. Going across each row, the width of the ports is varied
whereas its location is fixed. Keeping everything else fixed, as the width of the ports
increases, the size of the CW rotating primary vortex shrinks dramatically. For the top
two rows with S
o
= 2.5 and S
o
= 1.5, as the width of the ports increases, the CCW
vortices of the bottom left corner vanish and a CCW vortex gains prominence at the top
right corner. Going across the bottom row of Figure 3.5 that shows cases with both ports
on the left wall, it is observed that the two CCW rotating vortices at the two corners of
the right wall coalesce as the width of the port is increased.
The effect of the Reynolds number on the flow pattern is elucidated in Figure 3.6
for two position and width combinations of the outlet port, i.e. (a) W = 0.15 and S
o
=
0.925, and (b) W = 0.25 and S
o
= 1.5. When the two ports are positioned on the left wall
with W = 0.15 (Figure 3.6a), the rise in the value of the Reynolds number brings about a
more complicated flow field with the number of CW and CCW vortices rising. As many
as five vortices coexist for the Re = 500 case. The change in the number of vortices for
the cases with the outlet at the middle of the bottom wall (Figure 3.6b) is not observed,
however the strength of the two existing vortices increases as Re is raised.
3.3.2 Pressure Drop in the Cavity
The dimensionless pressure drop coefficient for the cavity is an important
parameter in the design of piezoelectric microvalves (Yang et al., 2004). This quantity is
related to the difference between the average pressures of the inlet and outlet ports by the
following equation:
.
2
1
2
in
outin
p
u
pp
C
?
?
= (3.13)
The average pressures are calculated by integrating the pressure over the inlet and outlet
ports as follows:
i
w
in
in
w
dlp
p
i
?
=
0
and
o
w
out
out
w
dlp
p
o
?
=
0
, (3.14)
86
87
where w
i
and w
o
are the widths of the inlet and outlet ports, respectively. The trapezoidal
rule of numerical integration was used to evaluate the above integrals. In Figure 3.7, the
variation of the pressure drop coefficient versus the position of the outlet port (S
o
) for
different Re numbers and widths of ports is shown. The pressure drop coefficient is a
strong function of the position of the outlet port. For the cases that the outlet port is on
the opposite or the same wall as the inlet port, the pressure drop coefficient is smaller in
comparison to cases where it is located on the adjacent walls (top and bottom). This
drastic change is illustrated clearly where the position of the outlet port changes from the
end of a vertical wall to a position at the beginning of a horizontal wall and vise versa.
But as it is seen in Figure 3.7, the pressure drop coefficient does not alter markedly as the
position of the outlet port changes on the same wall. The pressure drop coefficient attains
its maximum value with the outlet port located at the left side of the bottom wall, whereas
for most cases it gains the minimum value with the outlet port at the middle of the right
wall. For a couple of cases, the minimum is exhibited when the outlet port is at the
bottom of the left wall or the top of the right wall, although the values for these three
positions are very close to each other. It is also concluded that the pressure drop
coefficient is increasing with the increase of port?s width but is not a strong function of it.
For design considerations, a low-leakage valve is realized with the highest pressure drop
coefficient for a given flowrate (or Reynolds number). In the present configuration, such
an objective is obtained with both ports placed on the same wall.
88
3.3.3 Temperature Fields in the Cavity
The contours of dimensionless temperature (?) corresponding to nine cases with
Re = 500 and W = 0.25 are shown in Figure 3.8 for different positions of the outlet port.
There is a one-to-one correspondence between each cavity in this figure to those shown
in Figure 3.3. In effect, the influence of the complicated flow field on the temperature
field can be elucidated. The value of ? on the four solid walls is 1, whereas the value of ?
of the fluid entering the cavity is zero, and the contour levels are incremented by 0.05.
For the cases shown in the top row, the fluid temperature adjacent to the left, bottom and
lower part of the right walls exhibit moderate spatial variations. On the other hand, fluid
temperature gradients are marked next to the top portion of the right wall and next to the
top wall. The temperature gradient is always marked at the interface of the throughflow
stream and the CW rotating primary vortex suggesting intense heat exchange from the
fresh incoming fluid to the CW primary vortex. The core of the primary vortex is
generally isothermal. As the location of the outlet port is varied in the clockwise
direction, it is noted that the regions of intense fluid temperature gradients are found on
both sides of the outlet port. Marked temperature gradients are also noted on both sides
of the throughflow stream indicating intense heat exchange with the neighboring vortices.
The effect of the Reynolds number on the temperature field is highlighted in
Figure 3.9 for two position and width combinations of the outlet port, i.e. (a) W = 0.15
and S
o
= 0.925, and (b) W = 0.25 and S
o
= 1.5. The temperature fields shown in Figure
3.9 have one-to-one correspondence to the flow fields shown in Figure 3.6. When the
two ports are positioned on the left wall with W = 0.15, the rise in the value of the
Reynolds number brings about steeper fluid temperature gradients next to the left and
bottom walls, in addition to greater heat exchange on both sides of the throughflow.
With the outlet at the middle of the bottom wall, the temperature gradient adjacent to the
right wall is also seen to be enhanced as the Reynolds number is raised.
3.3.4 Variation of the Local Nusselt Number on the Walls of the Cavity
In order to evaluate the heat transfer rate along the walls, it is necessary to
observe the variations of the local Nusselt numbers on these walls. These are defined as:
,
walli
n
Nu
?
?
?=
?
(3.15)
with the subscript i representing b, l, r and t that correspond to the bottom, left, right and
top walls, respectively. The normal direction to each wall (?X or ?Y) is symbolized with
n. In order to present the local Nusselt number variations on the four walls
simultaneously, the use of the S coordinate system was adopted. For example, variations
of the Nu
l
, Nu
b
, Nu
r
and Nu
t
are plotted in graphical form for intervals S = 0 ~ 1, S = 1 ~
2, S = 2 ~ 3 and S = 3 ~ 4, respectively. This corresponds to traversing the inner wall of
the cavity in the counterclockwise direction starting at the top left corner. Variations of
the local Nusselt number on the four walls of the cavity for different positions of the
outlet port with Re = 500 and W = 0.25 are shown in Figure 3.10. In this figure, no data
for the Nusselt number is shown for S = 0 ~ 0.25 that corresponds to the location of the
inlet port. Discontinuities for each individual curve correspond to changing location of
the outlet port of the same width (W = 0.25). The Nusselt numbers shown in this figure
correspond to the flow and temperature fields shown in Figures 3.3 and 3.8, respectively,
and therefore are directly affected by those complex variations. The three corners (S = 1,
89
2 and 3) that do not have outlet ports on either side exhibit very poor heat transfer rate
due to their immediate proximity to small non-energetic vortices (Figure 3.3) that are
dominated by slow fluid motion. The two ends of the outlet port enjoy very intense heat
transfer rate due to the steep temperature gradients already discussed in Figure 3.8. In
between these two extreme rates of heat transfer, the Nusselt number exhibits a variety of
monotonic rise, montonic decline, and far more complicated variations involving local
minima and maxima points. These are directly dependent on the local status of the fluid
temperature gradient on the walls that were previously shown in Figure 3.8.
3.3.5 Variation of the Total Nusselt Number of the Cavity
Another variable utilized to evaluate the heat transfer rate is the total Nusselt
number of the cavity. Since all the four walls are active in heat exchange, the total
Nusselt number is the sum of the individual mean Nusselt numbers, iNu (i.e. integrals of
Equation 3.15 on respective walls). Therefore, the total Nusselt number is defined as:
,
4
1
4
1
12
2
1
?
?
?
?
==
i
ii
S
S
i
i
itot
SS
dSNu
NuNu
i
i
(3.16)
with S
i1
and S
i2
being the S-coordinates of the two ends of the i-th wall. The total Nusselt
number of the cavity as a function of the position of the outlet port for different Re
numbers and widths are given in Figure 3.11. The dependence of the total Nusselt
numbers take form of monotonically rising and decaying, piecewise functions that exhibit
local maxima when the outlet ports are positioned with one end at the three corners (S =
90
91
1, 2 and 3). Between the two outlet ports with one end at the corners, the one with a
higher S
o
coordinate consistently gives rise to a higher heat transfer rate. On the other
hand, the local minima total heat transfer of the cavity is achieved with the outlet port
located at the middle of the walls (S
o
= 1.5, 2.5 and 3.5). These important findings can be
used in design of practical heat exchange devices. By considering Figures 3.7 and 3.11,
it is clear that the best heat exchange cavity is realized when the outlet port?s position is
equal to (2+0.5W). This will accommodate both maximum heat transfer and minimum
pressure drop of the cavity.
3.4 Conclusions
1. For high Re and with the shortest distance between the inlet and outlet ports along the
top wall, a primary CW rotating vortex that covers about 75 to 88 percent of the
cavity is observed. Similar cases with smaller Re exhibit identical flow patterns but
with weaker vortices as Re is lowered. As the outlet port is moved along the right
wall toward the bottom wall, the CW primary vortex diminishes its strength, however
a CCW vortex that is present next to the top right corner covers a greater portion of
the cavity. With the outlet port moving left along the bottom wall, the CW primary
vortex is weakened further and the CCW vortex occupies nearly the right half of the
cavity.
2. For design purposes, inlet and outlet ports should be parallel to each other in order to
minimize pressure drop, even if they are designed to be placed at the same wall. In
order to realize excessive pressure drops, it is recommended that they are
92
perpendicular to each other. The distance between inlet and outlet ports and the
width size of ports do not greatly affect the value of pressure drop.
3. The temperature fields are directly related to the presence of the multiple vortices in
the cavity. Regions of high temperature gradient are consistently seen at the interface
of the throughflow with the neighboring vortices. Intense fluid temperature gradients
are observed next to solid walls on both sides of the outlet port.
4. Local Nusselt numbers are low at three corners when no outlet port is present in their
vicinity, whereas intense heat transfer rate is observed on the two sides of the outlet
port. Between these minima and maxima points, the local Nusselt number can vary
drastically depending on the flow and temperature fields adjacent to the respective
wall.
5. By placing the outlet port with one end at three corners, maximum total Nusselt
number of the cavity can be achieved. Minimum total heat transfer of the cavity is
achieved with the outlet port located at the middle of the walls. The most effective
location for the outlet port that accommodates maximum heat transfer and minimum
pressure drop is when it is located at (2+0.5W).
3.5 Closure
In this Chapter the laminar steady-state flow field and heat transfer in a sealed
square cavity with the one inlet and one outlet port was investigated. Simulations for 108
cases were performed for four Reynolds numbers (40, 100, 500 and 1000), 9 different
positions of the outlet port and 3 port widths (5, 15 and 25 percent of the side wall of the
cavity). It was found that in order to maximize the heat transfer rate and minimize the
93
pressure drop with the inlet port placed on top of the left wall, the best location for the
outlet port is on the bottom of the right wall. In Chapter 4, the inlet and outlet ports are
fixed at these positions and the effect of the oscillation of the inlet flow on the heat
transfer rate and pressure drop is studied.
w
i
T
w
94
Figure 3.1 Schematic diagram of a cavity with inlet and outlet ports
w
o
H
u
in
T
in
L
x
y
T
w
T
w
T
w
s
s = 0
?
?
s
o
GridSize
??
mi
n
?
25 50 75 100 125 150 175 200
0.021
0.022
0.023
0.024
0.025
0.026
0.027
Re= 1000
S
o
=2.1
W=0.2
Figure 3.2 The absolute value of the stream function at the center of primary vortex,?
min
,
as a function of the grid size (Re = 1000, S
o
= 2.1 and W = 0.2)
95
96
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
Figure 3.3 Streamlines for Re = 500 and W = 0.25 for 9 different positions of the outlet
port
97
?
?
?
? ??
?
?
?
?
?
?
?
?
?
?
?
?
Figure 3.4 Streamlines for Re = 40 and W = 0.15 for 9 different positions of the outlet
port
?
?
?
?
?
?
?
?
?
?
?
?
5.2=
o
S
5.1=
o
S
?
?
?
?
?
?
9.0?
o
S
05.0=W 15.0=W 25.0=W
Figure 3.5 Streamlines for Re = 100 with three outlet ports of widths W = 0.05, 0.15 and
0.25 positioned at three positions
98
?
?
?
?
?
?
40Re =
100Re =
10Re =
500Re =
?
?
?
?
?
?
?
?
?
?
(a) (b)
Figure 3.6 Streamlines for (a) W = 0.15 and S
o
= 0.925, and (b) W = 0.25 and S
o
= 1.5,
for Reynolds numbers of 10, 40, 100 and 500
99
S
o
C
p
123
10
0
10
1
Re= 10
Re= 40
Re= 100
Re= 500
S
o
C
p
123
10
0
10
1
10
2
Re= 10
Re= 40
Re= 100
Re= 500
S
o
C
p
123
10
0
10
1
10
2
Re= 10
Re= 40
Re= 100
Re= 500
(b)
(a)
(c)
Figure 3.7 Variations of the pressure drop coefficient of the cavity for different positions
of the outlet port with different Re numbers: (a) W = 0.05, (b) W = 0.15 and (c) W = 0.25
100
101
?
?
?
? ??
?
?
?
?
?
?
?
?
?
?
?
?
Figure 3.8 Temperature fields for Re = 500 and W = 0.25 for 9 different positions of the
outlet port [contour level increment of 0.05]
?
?
?
?
?
?
10Re =
40Re =
100Re =
500Re =
?
?
?
?
?
?
?
?
?
?
(b)(a)
Figure 3.9 Temperature fields for (a) W = 0.15 and S
o
= 0.925, and (b) W = 0.25 and S
o
=
1.5, for Reynolds numbers of 10, 40, 100 and 500 [contour level increment of 0.05]
102
s
Loc
a
l
N
u
s
s
e
l
t
N
u
m
b
e
r
01234
10
-2
10
-1
10
0
10
1
10
2
S
o
=0.875
S
o
=1.5
S
o
=1.875
S
o
=2.5
S
o
=2.875
S
o
=3.5
s
Loc
a
l
N
u
s
s
e
l
t
N
u
m
b
e
r
01234
10
-2
10
-1
10
0
10
1
10
2
S
o
=1.125
S
o
=2.125
S
o
=3.125
Figure 3.10 Variations of the local Nusselt number on the four walls of the cavity for
different positions of the outlet port with Re = 500 and W = 0.25
103
104
Figure 3.11 Variations of the total Nusselt number of the cavity for different positions of
the outlet port with different Re numbers: (a) W = 0.05, (b) W = 0.15 and (c) W = 0.25
S
o
To
t
al
N
u
s
s
e
l
t
N
u
m
b
e
r
123
5
10
15
20
25
Re = 10
Re = 40
Re = 100
Re = 500
S
o
To
t
al
N
u
s
s
e
l
t
N
u
m
b
e
r
123
5
10
15
20
Re = 10
Re = 40
Re = 100
Re = 500
(a)
(c)
(b)
S
o
To
t
a
l
N
u
s
s
e
l
t
Nu
m
b
e
r
123
0
5
10
15
20
25
30
35
40
Re = 10
Re = 40
Re = 100
Re = 500
105
CHAPTER 4 FLOW FIELD AND HEAT TRANSFER IN A CAVITY WITH
INLET AND OUTLET PORTS DUE TO INCOMING FLOW OSCILLATION
The most common technique for active control of heat transfer in various
convection-dominated heat exchange systems is oscillation of the fluid stream and/or a
boundary. A simplified model of a heat exchange unit is that of a cavity with one inlet
and one outlet port. Many practical problems can be simplified into this model. For
instance, mixed convection studies within cavities that contain heat sources are of great
interest in electronic cooling, ventilation of buildings and design of solar collectors. The
relevant literature for this class of problems was reviewed in Section 1.1.3. Simulations
of unsteady mixed convection in thermal storage applications have also received a great
deal of attention for many decades (Section 1.1.1). Other applications (mass transfer in
cavities, relaxation tanks employed in petroleum and transformer industries, transient
removal of contaminated fluid from cavities and food processing) were also noted by
Saeidi and Khodadadi (2004), who investigated steady laminar flow and heat transfer
within a cavity with inlet and outlet ports (Chapter 3). Nine different positions of the
outlet port on the right, left, bottom and top wall were studied with inlet Reynolds
numbers equal to 10, 40, 100 and 500.
In extending the work of Chapter 3, the objective of this study is to investigate the
transient and periodic characteristics of a cavity with inlet and outlet ports with an
oscillating velocity variation at the inlet port. The effects of the Strouhal number and the
oscillating Reynolds number are studied. The results will prove very helpful to design of
effective systems that might utilize oscillation of velocity to control heat transfer. The
scope of the current Chapter also has commonality with results of previous work on
resonant heat transfer enhancement in grooved channels that has attracted the attention of
the researches in past two decades. Some of the notable work that have been reported in
this area are Patera and Mikic (1986), Greiner (1991) and Nishimura et al. (1998).
4.1 Problem Formulation
The schematic drawing of the system and the coordinates is shown in Figure 4.1.
The height and width of the cavity are denoted by H and L, respectively. In this study,
only a square cavity is considered (H = L). The depth of the enclosure perpendicular to
the plane of the diagram is assumed to be long. Hence, the problem can be considered to
be two-dimensional. Based on the steady results for this problem (Chapter 3), in order to
maximize heat transfer and minimize the pressure drop, the inlet and outlet ports are
placed at the top of the left wall and the bottom of the right wall, respectively. The inlet
and outlet port widths are equal (w
i
= w
o
= w) and fixed to 15 percent of the length of the
side wall. The velocity at the inlet port is uniform and time-dependent:
tf2sinuuu
omin
?+=
, (4.1)
where f is the frequency of the oscillation. Quantities u
m
and u
o
(u
m
? u
o
) are the mean
and oscillating inlet velocities, respectively (Figure 4.2). The inlet fluid is maintained at
the temperature ( ) different from the remaining walls of cavity ( ), with .
in
T
w
T
inw
TT >
106
4.1.1 Governing Equations
The fluid within the enclosure is an incompressible fluid and the its properties are
constant. The flow within the enclosure is assumed 2-D, unsteady and laminar. The
gravity effect is negligible. The governing equations of continuity, momentum and
thermal energy are:
,0=
?
?
+
?
?
y
v
x
u
(4.2)
,
1
2
2
2
2
?
?
?
?
?
?
?
?
?
?
+
?
?
+
?
?
?=
?
?
+
?
?
+
?
?
y
u
x
u
x
p
y
u
v
x
u
u
t
u
?
?
(4.3)
,
1
2
2
2
2
?
?
?
?
?
?
?
?
?
?
+
?
?
+
?
?
?=
?
?
+
?
?
+
?
?
y
v
x
v
y
p
y
v
v
x
v
u
t
v
?
?
(4.4)
,
2
2
2
2
?
?
?
?
?
?
?
?
?
?
+
?
?
=
?
?
+
?
?
+
?
?
y
T
x
T
y
T
v
x
T
u
t
T
?
(4.5)
in which u, v, p and T are the velocity components in the x and y directions, pressure and
temperature, respectively. Quantities ?, ? and ? are the density, kinematic viscosity and
thermal diffusivity of the fluid, respectively. In the energy equation, the viscous
dissipation terms are neglected.
At the inlet port a time-dependent velocity and a constant temperature are
imposed. The no-slip condition is imposed on the solid surfaces. At the outlet port, the
107
flow and temperature are assumed fully-developed. So the boundary conditions can be
expressed as follows:
At the inlet (x = 0, y = (H - w
i
) to H):
108
)(tfuuu
omin
?2sin+= , v = 0, T = T
in
. (4.6a)
On the four solid walls:
u = v = 0, T = T
w
. (4.6b)
At the outlet port (x = H, y = 0 to w
o
):
0,0 =
?
?
=
?
?
x
T
x
u
. (4.6c)
4.1.2 Dimensionless Form of the Governing Equations
Dimensionless form of the governing equations can be obtained via introducing
dimensionless variables. The lengths can be scaled by the cavity length H and the
velocities can be scaled by the mean inlet fluid velocity, u
m
. As for the temperature, the
two extreme values (T
in
and T
w
) are used. The dimensionless variables are then defined
as:
.,,
2
,,,,
*
H
tu
t
TT
TT
u
p
P
u
v
V
u
u
U
H
y
Y
H
x
X
m
inw
in
m
mm
=
?
?
=?=
====
?
(4.7)
The governing equations for continuity, momentum and thermal energy are then
written in dimensionless form:
,0
Y
V
X
U
=
?
?
+
?
?
(4.8)
?
?
?
?
?
?
?
?
?
?
+
?
?
+
?
?
?=
?
?
+
?
?
+
?
?
2
2
2
2
Re
12
*
Y
U
X
U
H
w
X
P
Y
U
V
X
U
U
t
U
m
, (4.9)
?
?
?
?
?
?
?
?
?
?
+
?
?
+
?
?
?=
?
?
+
?
?
+
?
?
2
2
2
2
2
Re
1
*
Y
V
X
V
Y
P
Y
V
V
X
V
U
t
V
m
H
w
, (4.10)
?
?
?
?
?
?
?
?
?
??
+
?
??
=
?
??
+
?
??
+
?
??
2
2
2
2
RePr
1
2
*
YX
H
w
Y
V
X
U
t
m
. (4.11)
The mean Reynolds number is defined as Re
m
= u
m
2w / and the Prandtl
number is Pr = . For t
?
?? /
*
> 0, the dimensionless form of the boundary conditions can
be expressed as follows:
Nothing that w / H = 0.15, at the inlet (X = 0, 0.85 ? Y ? 1):
),2sin(1
*
tSt
u
u
U
m
o
??+= ?
and 0=V
0=?
, (4.12a)
On the walls:
1,0 =?==VU
, (4.12b)
At the outlet (X = 1, 0 ? Y ? 0.15):
0=
?
?
X
U
and
0=
?
??
X
. (4.12c)
The Strouhal number that is a dimensionless parameter widely used in study of
oscillating fluid flow phenomena is defined as
m
u
fH
St
?
=
.
The range of variation of the
Strouhal number is illustrated schematically in Figure 4.3. Note that the quantity ? is the
period of oscillation that is the inverse of f. Therefore, Re
m
, Pr, St and u
o
/ u
m
are the
109
dimensionless groups that govern this problem. The instantaneous Reynolds number is
defined as
?
=
w2u
Re
in
. Thus, the instantaneous Reynolds number varies sinusoidally as
the inlet velocity oscillates. For this study, the ratio u
o
/ u
m
is set to 2/3. Thus, with Re
m
= 300, the instantaneous Reynolds number varies from 100 to 500 (mean value of 300).
Five cases are studied with Strouhal numbers of 0.1, 0.5, 1, 2 and 10, respectively and the
Prandtl number of the fluid is fixed to 5.
4.2 Computational Details
The governing equations were solved by the finite-volume-method using
Patankar's (1980) SIMPLE algorithm. A two-dimensional uniformly-spaced staggered
grid system was used. Hayase et al.'s (1992) QUICK scheme was utilized for the
convective terms, whereas the central difference scheme was used for the diffusive terms.
In order to keep consistent accuracy over the entire computational domain, a third-order-
accurate boundary condition treatment suggested by Hayase et al. (1992) was adopted.
Details of a grid independence test using six different grid densities (75 ? 75, 100 ? 100,
125 ? 125, 150 ? 150, 175 ? 175 and 200 ? 200) were given in Chapter 3 of this
dissertation for the steady version of this problem. A 160 ? 160 uniform grid was
selected for all cases reported here. The steady-state fluid flow and temperature fields for
a case with Re
m
= 300 is treated as the starting field (t = 0) in this study. For the unsteady
computations, tolerance of the normalized residuals upon convergence is set to 10
-5
for
every time step. The under-relaxation parameters for u, v, and T are all set to 0.6,
whereas under-relaxation parameter for pressure correction is set to 0.3. The average
110
111
number of iterations needed to achieve convergence at every time step were 110, 254,
381, 598 and 1117 for the Strouhal numbers of 0.1, 0.5, 1, 2 and 10, respectively. The
three-time-level-method that is a 2
nd
order implicit scheme is used to approximate the
unsteady terms. For a period of oscillation (?), four different time steps ?t = ? / 72, ? /
90, ? / 120 and ? / 180 (St = 10) were adopted for the time step independence test. After
comparing the variations of the mean Nusselt number of the system shown in Figure 4.4
(St = 10), the time step ?t = ? / 120 was found small enough to be used for all the cases in
this study. All the calculations were performed on a SGI Altix 350 of the Alabama
Supercomputer Center (Huntsville, Alabama).
4.3 Results and Discussions
Before presenting the unsteady flow and thermal fields, the behavior of steady
velocity and temperature fields are briefly described. The streamlines and temperature
contours corresponding to the steady cases with Re = 100, 300 and 500 are shown in
Figure 4.5a-b. The streamlines are shown on the left column (Figure 4.5a), whereas the
temperature contours are presented on the right column (Figure 4.5b), thus allowing one
to clearly observe the influence of the complicated flow field on the temperature field.
The general features of the flow fields are similar for the three Re numbers studied. It is
observed that the streamlines representing the throughflow of the incoming fluid
(streamline values of zero and inlet mass flowrate) do not coincide with the two ends of
the outlet ports. The area occupied by the clockwise (CW) rotating primary vortex
diminishes very little as the Re is increased, whereas a counterclockwise (CCW) rotating
vortex next to the top right corner covers a larger portion of the cavity. Another CCW
rotating vortex is permanently featured at the bottom left corner with its strength growing
with the rise of the Reynolds number. As for the contours of the dimensionless
temperature ( ) shown in Figure 4.5b, the value of ? ? entering the cavity is zero,
whereas on the four solid walls it is 1, and the contour levels are incremented by 0.05.
The temperature gradient is always marked at the interface of the throughflow stream and
the neighboring CW and CCW vortices suggesting intense heat exchange from the fresh
incoming fluid to those rotating structures. The core of the primary CW vortex is
generally isothermal. Regions of intense fluid temperature gradients are found on both
sides of the outlet port. A rise in the value of the Reynolds number brings about steeper
fluid temperature gradients next to the left and bottom walls, in addition to greater heat
exchange on both sides of the throughflow.
Starting with the steady state velocity and temperature fields for Re = 300, the
oscillating incoming velocity of the fluid was imposed at t = 0. After a certain period of
time, the fluid flow and temperature fields will reach their respective periodic states. The
evolution of the periodic flow and temperature fields, the transient response of the mean
Nusselt numbers on the four walls and the time dependence of the Nusselt number of the
system are discussed next.
4.3.1 Time Required to Reach the Periodic State
In this study, the initial field is the steady-state solution for Re = 300. The system
will take a certain period of time to reach the periodic state. Due to uncoupling of the
momentum and energy equations (natural convection is neglected), the fluid flow and
112
temperature fields could reach their respective periodic states at different times. The
cycle-to-cycle variations of the fluid flow and temperature fields were quantitatively
defined and strictly monitored, exhibiting monotonic decay with time. The time required
to reach the periodic state is an overall characteristic parameter for this system. In order
to obtain this parameter, the criterion needed to determine if the system has reached the
periodic state needs to be defined. For a system reaching its periodic state, the cycle-to-
cycle variation should be negligible. Two variables ?
?
and ?
?
are introduced to indicate
the cycle-to-cycle variations of the fluid flow and temperature fields. For the n-th cycle,
the cycle-to-cycle variation of fluid flow can be defined as:
[][]
[][]
?
=
??+?
??
??+???+?
??
?
=
N
k
tkn
ji
NYjNXi
tkn
ji
tkn
ji
NYjNXi
N
1
)1(
,
,1,,1
)2(
,
)1(
,
,1,,1
max
max
1
?
??
?
?
??
?
, (4.13)
where N, NX and NY are the number of the time steps within one period (120 in all the
cases studied), the number of grid points in the x and y direction (160 ? 160),
respectively. The quantity is the stream function value at grid point (i, j) at
time instant t = (n-1)?+k??t.
tkn
ji
??+? ?
?
)1(
,
Similarly, the cycle-to-cycle variation of the temperature field for the n-th cycle
can be defined as:
[][]
?
=
??+???+?
??
?=
N
K
tkn
ji
tkn
ji
NYjNXi
N
1
)2(
,
)1(
,
],1,,1
max
1
??
?
???
. (4.14)
113
The measures of the cycle-to-cycle variation of the stream function field (?
?
) as a
function of t / ? for different Strouhal numbers are shown in Figure 4.6. It takes a greater
number of cycles to achieve the periodic state for cases with high Strouhal numbers in
comparison to the cases with smaller St numbers. Note that the y axis is in logarithmic
scale. Similarly, in Figure 4.7 the variation of ?
?
as a function of t / ? is shown that
exhibits similar trends as ?
?
. Comparing Figures 4.6 and 4.7, the flow and temperature
fields reach their periodic states at different times. For all the cases investigated, it takes
fewer cycles for the flow field to reach the periodic state in comparison to the
temperature field.
In order to define the criteria needed to decide when the periodic solution is
achieved, one case (St = 0.1) was selected as the test case. Qualitative examination of the
streamline and temperature fields (to be discussed in Figures 4.8 and 4.9), showed that it
takes 5 cycles (?
?
= 0.00013) to reach the periodic state for the flow field, whereas 18
cycles (?
?
= 0.00079) are needed for the temperature field. Based on this, the following
criteria for attaining the periodic state are introduced:
0007.0
0001.0
<
<
?
?
?
?
. (4.15)
These values are smaller than those observed for the test case, so it is more
stringent in comparison to the qualitative inspection of the flow and thermal fields.
According to this criteria, the number of cycles (N
f
and N
t
) needed to reach the periodic
state are summarized in Table 4.1. The number of cycles needed for the thermal field to
reach the periodic state is found to be greater than the velocity field. The difference
between N
t
and N
f
increases as St increases.
114
Table 4.1 Number of cycles needed to reach the periodic states for fluid flow (N
f
)
and temperature fields (N
t
)
N
f
, N
t
St = 0.1 St = 0.5 St = 1 St = 2 St = 10
5, 18 13, 21 23, 51 45, 93 57, 107
4.3.2. Periodic Evolution of the Flow Field
The periodic flow fields during the 13
th
cycle for the case with St = 0.5 are
presented for every 30 degrees phase angle in Figure 4.8. During the first quarter of the
cycle (i.e. ), as the Reynolds number increases (i.e. 300 ? Re ? 500), more
throughflow enters the cavity. In effect, the area covered by the throughflow is enlarged
and the CCW vortex of the top right corner (Figure 4.5a) is washed away. The primary
CW rotating vortex does not greatly change, however it interacts with two rotating
vortices on the left wall. During the second and third quarters of the cycle
( ), when the Reynolds number is decreasing (i.e. 100 ? Re ? 500) the
throughflow diminishes and its coverage area is observed to shrink drastically.
Meanwhile the CW rotating primary vortex expands by merging with two small vortices
on the left wall. The small vortices were located at the bottom left corner of the cavity
(CCW) and right below the inlet port (CW), respectively. During the fourth quarter of
the cycle, (i.e.
2/0 ????
2/32/ ?????
????? 22/3 and 100 ? Re ? 300), the primary vortex shrinks with the
raising of the Reynolds number and the two small vortices that had merged with it
115
reappear on the left wall. A very interesting phenomenon that is observed during the
fourth quarter is the penetration of the throughflow into the space between the left wall
and the primary CW vortex. The throughflow is observed to touch the middle of the left
wall near the end of the cycle. Starting with the second quarter of the cycle, a CCW
vortex that is formed next to the top right corner of the cavity strengthens and its region
of influence changes as the Reynolds number reduces. By the end of the fourth quarter,
this vortex diminishes and finally, it is washed away by the energetic throughflow.
4.3.3. Periodic Evolution of the Temperature Field
The contours of dimensionless temperature fields during the 21
st
cycle for the case
with St = 0.5 are presented for every 30 degrees phase angle in Figure 4.9. There is a
one-to-one correspondence between each cavity in this figure to those shown in Figure
4.8. In effect, the influence of the complicated flow field on the temperature field can be
elucidated. Similar to the temperature contours for the steady cases (Figure 4.5b), ?
values are 0 and 1 at the inlet plane and on the four solid walls, respectively, and the
contour levels are incremented by 0.05. The temperature within the core of the CW
rotating primary vortex remains fairly uniform throughout the cycle, whereas steep
temperature gradients are observed on its boundary. Marked heat exchange takes place
on the boundary of the primary vortex where it interacts with the evolving throughflow
and the two rotating vortices on the left wall. The temperature gradients next to the left
wall do not change drastically since the two rotating vortices on this surface are sites of
low speed flow. Similarly, portions of the top and right walls closest to the top right
corner do not experience steep variations in temperature gradient due to the frequent
116
presence of a low-speed CCW rotating vortex next to this corner. The most marked
temperature gradients and therefore the greatest heat exchange between the fluid and the
walls are observed on the left segment of the top wall, the bottom segment of the right
wall and the right segment of the bottom wall.
4.3.4. Transient Evolution of the Mean Nusselt Numbers on the Four Walls
The transient evolution of the instantaneous mean Nusselt numbers on the four
walls are shown in Figures 4.10 and 4.11 for St = 0.1 and 10, respectively. The
definitions of the instantaneous mean Nusselt numbers on the walls are given as surface
integrals of the instantaneous local Nusselt numbers, i.e.:
,ds
n
ds*)t,s(Nu*)t(Nu
1
0
wall
1
0
ii
??
?
??
?==
(4.16)
with the subscript i denoting b, l, r or t and variable s being X or Y. Variable n will take
on X or Y opposite to s. For comparison purposes, the mean Nusselt numbers on the four
walls under steady state conditions for Re = 100, 300 and 500 are summarized in Table
4.2.
Table 4.2 Mean Nusselt numbers on the four walls under steady-state conditions
Re
l
Nu
b
Nu
r
Nu
t
Nu
100 4.13 14.2 10.1 15.8
300 5.56 19.4 16.7 20.9
500 6.79 30.3 28.0 27.9
117
As for the variations of the instantaneous mean Nusselt numbers on the walls in Figures
4.10 and 4.11, it is observed that regardless of the Strouhal number, the Nusselt number
on the left wall is consistently the lowest among the four walls. This was already
discussed above in connection with the presence of two rotating vortices next to the left
wall. For the lowest Strouhal number case (St = 0.1) shown in Figure 4.10, one can
observe that the oscillations of the mean Nusselt numbers on the four walls are easily
recognizable. This is due to the long period of oscillation in relation to the convective
time scale (
?== Stu/Ht
mconv
). In other words, as the inlet velocity oscillates
slowly, convection has enough time to transmit the oscillation of the inlet velocity to
every point within the cavity. The amplitude of the oscillations exhibited on the four
walls are also marked in general except for the left wall. For the highest Strouhal number
case (St = 10) shown in Figure 4.11, the time for the fluid to travel from one side of the
cavity to another side is greater than the period of the velocity?s oscillation. Therefore,
the signature of flow oscillation at the inlet port is only observed on the left and top walls
of the cavity. The top wall is in the closest proximity of the throughflow that brings in
the oscillation information, whereas the left wall also feels the penetration of the
throughflow discussed earlier. The amplitude of oscillations of the Nusselt numbers on
the bottom and right walls are heavily damped compared to the low frequency oscillation
case of Figure 4.10. This can be explained in light of the minimal intraction of the
throughflow with these walls because of the frequent presence of the vortices next to
them.
4.3.5. Variation of the Total Nusselt Number of the Cavity
118
Another variable utilized to evaluate the heat transfer rate is the total Nusselt
number of the cavity. Since all the four walls are active in heat exchange, the total
Nusselt number is the sum of the individual mean Nusselt numbers ( iNu ). Therefore,
the total Nusselt number of the cavity is defined as:
,
)(
4
1
)(
4
1
)(
12
*
**
2
1
?
?
?
?
==
i
ii
S
S
i
i
itot
SS
dStNu
tNutNu
i
i
(4.17)
with S
i1
and S
i2
being the S-coordinates of the two ends of the i-th wall. The
instantaneous total Nusselt number of the cavity for different Strouhal numbers are given
in Figure 4.12. The two low frequency cases with St = 0.5 and 0.1 are shown in Figure
4.12a, whereas the cases with St = 1, 2 and 10 are shown in Figure 4.12b. The periodic
nature of the total Nusselt number of the system is consistently maintained regardless of
the Strouhal number. The amplitude of oscillations consistently decrease as the Strouhal
number is raised covering the low to high frequency cases studied. Except for a short
duration during each cycle for St = 0.1, heat transfer rate from the cavity is always
enhanced in comparison to the steady-state performance of the system for Re = 300. For
the Strouhal numbers of 0.5 and 1, the instantaneous heat transfer attained after 15 cycles
is always above the performance for the steady case with Re = 500. This suggests that
oscillating the inlet velocity about a Strouhal number of unity when the frequency of
oscillation resonates with the convection time scale is the optimum condition for this
system.
119
120
4.4 Conclusions
1. Both the flow patterns and temperature fields in a cavity due to an oscillating
velocity at the inlet port reach their respective periodic states after certain time durations.
It takes more cycles for a temperature field to reach its periodic state in comparison to the
corresponding flow field. For cases with greater Strouhal numbers, it takes more cycles
to reach a periodic state.
2. The oscillation of the velocity at the inlet port causes the cyclic growth and
decay of the throughflow. The throughflow stream is in constant contact with the CW
rotating primary vortex, which in turn interacts with two rotating vortices on the left wall.
A CCW rotating vortex at the top right corner also experiences periodic growth and
decay.
3. The temperature field is directly affected by the flow field to the extent that it
supports minimal heat transfer on the left wall. In contrast, certain segments of the other
three walls and the boundary of the throughflow are zones of active heat exchange.
4. For the low Strouhal number of 0.1, the mean Nusselt numbers on the four
walls clearly exhibit large amplitudes of oscillation and periodicity. On the other
extreme, with St = 10, the amplitudes of oscillation are degraded. These behaviors are
directly linked to the relation between the period of oscillation and the convection time
scale.
5. Regardless of the Strouhal number, heat transfer enhancement in comparison to
the performance of the system for the mean Reynolds number is consistently observed.
The best heat transfer rate is realized when the Strouhal number is close to unity. This is
121
the case when the period of the incoming stream resonates with the convection time
scale.
4.5 Closure
In completing the steady-state parametric study of fluid flow and heat transfer in a
cavity with inlet and outlet ports (Chapter 3), the effect of inlet flow oscillation was
discussed in this Chapter. Even though only one geometrical configuration (fixed inlet
and outlet locations and fixed width) in combination with fixed Re
m
and u
o
/ u
m
values
was studied, the results point to the positive influence of flow oscillation on heat transfer
enhancement.
H
u
in
T
in
x
y
T
w
T
w
T
w
w
o
w
i
T
w
?
?
L
Figure 4.1 Schematic diagram of a cavity with inlet and outlet ports
122
u
in
u
o
?
u
m
t
Figure 4.2 Velocity oscillation diagram
t
diff
= H
2
/ ? = (H / 2 w).Re?t
conv
t
t
conv
= H / u
m
123
Figure 4.3 Convection and diffusion time scales in relation to the
period of the oscillation (Re > 1)
Fast Oscillation Slow Oscillation
? = 1 / f
St < 1 St > 1 St=1
High Frequency Low Frequency
t
*
Nu
tota
l
0.01 0.02 0.03 0.04 0.05
18.3
18.4
18.5
18.6
18.7
18.8
18.9
?t=?/72
?t=?/90
?t=?/120
?t=?/180
Figure 4.4 Comparison of the total mean Nusselt number of the system for various time
steps (St = 10)
124
?
?
?
?
?
?
?
?
?
?
?
?
Re = 300
Re = 500
Re = 100
(a) (b)
Figure 4.5 Dependence of the (a) streamlines and (b) temperature contours for the steady
case with Re = 100, 300 and 500 [contour level increment of 0.05 for the temperature
field]
125
Cycle ( t / ? )
?
?
10 20 30 40 50 60
10
-4
10
-3
10
-2
10
-1
10
0
St= 10
St= 2
St= 1
St= 0.5
St= 0.1
Figure 4.6 Cycle-to-cycle changes of the streamline values over the entire domain for
different Strouhal numbers
126
Cycle ( t / ? )
?
?
10 20 30 40 50 60
10
-4
10
-3
10
-2
10
-1
10
0
St= 10
St= 2
St= 1
St= 0.5
St= 0.1
Figure 4.7 Cycle-to-cycle changes of the temperature values over the entire domain for
different Strouhal numbers
127
? = ? ? = 5?/6
? = 7?/6 ? = 4?/3 ? = 3?/2
? = 2?/3
? = ?/2 ? = ?/3 ? = ?/6
? = 5?/3 ? = 11?/6 ? = 2?
Figure 4.8 Periodic evolution of the streamlines for St = 0.5 during the 13
th
cycle for a
6/? phase angle increment
128
? = ? ? = 5?/6
? = 7?/6 ? = 4?/3 ? = 3?/2
? = 2?/3
? = ?/2 ? = ?/3 ? = ?/6
? = 5?/3 ? = 11?/6 ? = 2?
Figure 4.9 Periodic evolution of the temperature fields for St = 0.5 during the 21
st
cycle
for a phase angle increment [contour level increment of 0.05] 6/?
129
t/?
W
a
l
l
N
u
s
s
e
l
t
N
um
be
r
s
024681012
8
12
16
20
24
28
32
36
40
44
48
52
Nu
l
Nu
b
Nu
r
Nu
t
Figure 4.10 Transient evolution of the instantaneous mean Nusselt numbers on the four
walls for St = 0.1
t/?
W
a
l
l
N
u
s
s
e
l
t
N
um
be
r
s
02468101214161820
0
4
8
12
16
20
24
28
Nu
l
Nu
b
Nu
r
Nu
t
Figure 4.11 Transient evolution of the instantaneous mean Nusselt numbers on the four
walls for St = 10
130
Re = 500
Re = 300
Re = 100
Cycle ( t / ? )
M
e
a
n
Nu
s
s
e
l
t
Nu
m
b
e
r
0 5 10 15 20 25 30 35 40 45 50 55 60
10
15
20
25
30
St= 0.5
St= 0.1
(a)
Re = 500
Re = 300
Re = 100
Cycle ( t / ? )
M
ean
N
u
s
s
e
l
t
N
u
m
b
er
0 5 10 15 20 25 30 35 40 45 50 55 60
10
15
20
25
30
St= 10
St= 2
St= 1
(b)
Figure 4.12 Transient evolution of the instantaneous total Nusselt number of the system:
(a) St = 0.1 and 0.5, and (b) St = 1, 2 and 10
131
132
CHAPTER 5 COMPUTATIONAL MODELING OF A PIEZOELECTRICALLY-
ACTUATED MICROVALVE FOR THE CONTROL OF LIQUID FLOWRATE
The packaging of a number of microfluidic devices such as micropumps,
microvalves and microfilters are characterized by an enclosure with multiple fluid inlet
and outlet ports. This feature of these microfluidic devices is in common with the
objectives of this research directed at studies of flow and heat transfer in enclosures with
inlet and outlet ports. In exploiting this commonality, the results of a research study of
liquid flow in a complex three-dimensional microvalve are given in this Chapter.
Microvalves that are compatible with both liquids and gases are envisioned to be
integral parts of proposed micropropulsion units for microspacecrafts and microsatellites.
The design, operation and performance of such a leak-tight piezoelectrically-actuated
microvalve for high-pressure gas micropropulsion was recently reported by Yang et al.
(2004). In the process of enhancing this class of microvalves, a leak-tight, low-power,
liquid-compatible, piezoelectrically-actuated microvalve has been designed, fabricated,
and characterized by researchers of the MEMS Technology Group at NASA JPL. The
microvalve consists of a custom-designed piezoelectric stack actuator bonded onto silicon
valve components. The microvalve has a silicon membrane called boss plate (Figure 5.1)
for isolating the piezoelectric actuator from the liquid effluents. The valve seat
configurations include 12 narrow-edge seating rings and tensile-stressed silicon tethers
133
that enable the normally closed and leak-tight operation. A concentric series of narrow
rings (Figures 5.2 and 5.3) simulate a ?knife-edge? seal by reducing the contact area,
thereby increasing the seating pressure and consequently reducing internal leaks (Lee and
Yang, 2004).
5.1 Geometry of the Microvalve Model
The flow enters the microvalve from four inlet ports positioned at the corners of
the upper boss plate (Figure 5.3b) and exits the valve from an outlet port fabricated in the
middle of the seat plate (Figure 5.3a). The diameter of the inlet and outlet pipe is 200
?m. The length of the inlet pipes is 7100 ?m, whereas the length of the outlet pipe is
600 ?m. The actual size of the microvalve under investigation is shown in Figure 5.4.
The liquid enters a 10 ?m high box after passing through the inlet pipe and then flows
into a bigger box with 6 mm ? 5.2 mm ? 0.56 mm dimensions. There are twelve (12)
thin seat rings with different radii that surround the outlet hole. The thickness of the
rings is 1 ?m and their height is 10 ?m. The diameter of the innermost ring is 110 ?m
and other rings are evenly spaced concentrically at a distance of 150 ?m from each other.
A square-shaped lower boss plate is shown on top of the outlet port with dimensions of
1.8 mm ? 1.8 mm ? 0.4 mm. The flow is directed to the outlet pipe after passing through
the space between the lower boss plate and the rings. The spacing between the lower
boss plate and the seat plate ranges from 11 to 13.5 ?m, which means the distance
between the top of the valve seat rings and lower boss plate varies between 1 and 3.5 ?m.
Three deflections of 1, 2.7, and 3.5 ?m were studied in great detail.
134
The bulk of this Chapter is devoted to presenting the CFD-based results of
simulation of flow in the complicated 3-dimensional microvalve. Results of a 2-
dimensional CFD analysis are also presented for the flow between the boss and seat
plates. The 2-dimensional model is studied in order to find a more accurate estimation of
the pressure drop over the rings. In addition, the CFD-based results will be compared to
simplified 2-D analysis of Stokes flow and experimental observations.
5.2 Problem Formulation for the 3-D Model
A 3-D structured mesh with 330,000 hexahedral cells that represents the major
features of the microvalve was generated for the numerical model. Due to the symmetry
of the microvalve about two planes, only a one-quarter section of the microvalve was
necessary for the model. An isometric view of a one-quarter section of the microvalve is
shown in Figure 5.5. The solid components are depicted in shades of grey, whereas the
fluid is shown in red. The composite model of the microvalve is shown in Figure 5.6.
The symmetry planes are demonstrated in red color. The mesh is denser in regions that
experience excessive pressure drop (i.e. between the boss and seat plates and where the
flow enters the microvalve). Different views of the computational mesh are shown in
Figures 5.7 to 5.10. The green lines identify where the various blocks meet.
The continuum-based continuity and momentum equations are valid for this
problem (Karniadakis and Beskok, 2001). The fluid is incompressible and Newtonian,
whereas the flow is taken to be laminar and steady. The governing equations of
continuity and momentum in index notation form are:
,0=
?
?
i
i
x
u
(5.1)
() .
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
+
?
?
?=
?
?
i
j
j
i
ji
ji
j
x
u
x
u
xx
p
uu
x
??
(5.2)
In view of the complexity of the geometry in comparison to the problems of
Chapters 3 and 4, version 6.2 of the commercial code FLUENT was utilized for solving
the governing equations. The no-slip boundary condition was applied on the walls,
whereas symmetry boundary conditions were used for the planes of symmetry. A second
order upwind scheme was chosen to discretize the momentum equation, whereas the
SIMPLE algorithm was used to couple the pressure field and velocities. Two different
types of boundary conditions were applied to the inlet and outlet ports. At first, a gage
total pressure was specified at the inlet port and the outlet total pressure was set to zero
(gage). In order to verify the results, the mass flow rate of the system, which was found
by applying the above boundary conditions, was chosen as the inlet boundary condition
and flow at the outlet was assumed fully developed. The problem was solved again with
the assigned boundary condition. The same total pressure was found at the inlet. A total
of 15 cases were run for 3 deflections of 1, 2.7 and 3.5 ?m (Figure 5.10) and 5 inlet total
pressures of 3, 5, 10, 15 and 20 psig with water as the working fluid. The Reynolds
numbers based on the outlet pipe diameter for different cases are reported in Table 5.1.
The experimental mass flow rates are used to calculate the mean velocity at the outlet of
the microvalve. The Reynolds number varies from 0.63 to 6.35
135
136
Table 5.1 The Reynolds number at the outlet port of the microvalve for different
deflections at different total pressures
?P
t
= P
t inlet
- P
t outlet
Applied
Voltage (V)
Deflection
(microns)
3 psig 5 psig 10 psig 15 psig 20 psig
10 1 0.63 1.27 2.12 2.54 2.86
20 2 1.16 2.33 2.96 4.23 4.87
30 2.7 1.48 2.54 3.60 4.97 5.71
40 3.5 1.69 3.17 3.91 5.61 6.35
5.2.1 Possible Non-Continuum Effects
For the macro-scale problems, the no-slip boundary condition for the wall
surfaces is acceptable and widely used in numerical studies. But recent experimental
studies (Derek et al. 2002 and 2005) in micro- and nano-scale media indicate that the
proper boundary condition depends both on the characteristic length scale of the flow and
the chemical and physical properties of the solid surface. Therefore, the importance of
boundary conditions increases as the dimensions of microfluidic devices decrease.
Although the current study is within the micro-scale range, the no-slip boundary
condition is still valid, because the working fluid is a liquid (water). In gaseous systems,
the slip boundary conditions should be considered to obtain reliable results.
5.3 Radial Stokes Flow between Two Concentric Parallel Disks
In order to be able to assess the validity of the results of 3-D CFD simulations,
comparison to limiting cases or simplified models are always warranted. In view of this,
a simplified linearized Stokes analysis was carried out for the radial region near the
outlet. This analysis is deemed to be appropriate since the Reynolds numbers of the flow
is known to be between 0.0008 and 0.008 at the outer surface of the rings.
The flow between the seat and boss plates can be simplified as an axisymmetric
radial flow between two concentric parallel disks. The geometry is shown in Figure 5.11.
The distance between the disks is 2h. The effects of the rings on the flow and pressure
drop are neglected. In addition, it is assumed that for this axisymmetric flow, the velocity
components in the z and ? directions are negligible ( 0?
r
V ). The flow is laminar,
incompressible and steady.
5.3.1 Determination of the Radial Velocity
By applying these assumptions to the governing equations, the continuity and
momentum equations in the r and z directions are simplified to (Bird, 1960):
,0)( =
?
?
r
rV
r
(5.3)
,
2
2
z
V
r
P
r
V
V
rr
r
?
?
+
?
?
?=
?
?
?? (5.4)
,0=
?
?
z
P
(5.5)
where and P are the velocity in the radial direction and the fluid static pressure,
respectively. Quantity ? is the density and ? represents the viscosity of the working fluid.
r
V
Ideally, the above equations should be solved analytically to find the pressure
drop (or rise) as a function of the mass flow rate. In doing so, introduce a new variable f
137
( ) that automatically satisfies the continuity equation (Eq. 5.3). The new form
of the momentum equation in the r direction will be:
r
rVzf =)(
3
2''
r
f
r
f
dr
dP
?? += . (5.6)
Note that pressure is just a function of r (Eq. 5.5), so the partial differentiation of pressure
in Eq. 5.4 is changed to ordinary differentiation in Eq. 5.6.
Equation 5.6 is not solvable analytically because of the nonlinear term (
3
2
r
f
).
This term represents the convective effect in the Navier-Stokes equations. This term can
be neglected by assuming very low Reynolds number (very low fluid velocity, i.e. Stokes
flow). The Reynolds number analysis shows that the maximum Reynolds number
between the boss and seat plates is 0.008. Therefore, the flow can be safely assumed to
be in the Stokes regime and the convective term in the momentum equation is negligible.
Equation 5.6 then becomes:
r
f
dr
dP
''
?= . (5.7a)
Upon separating variables, each side depends on independent variables r and z leading to:
.
''
constf
dr
dP
r ==?
(5.7b)
Upon integrating twice, one gets:
21
2
2
)( CzCz
dr
dPr
zf ++=
?
. (5.8)
The physical boundary conditions are:
At : , or f = 0, (5.9a) hz ?= 0=
r
V
138
At : 0=z 0=
?
?
z
V
r
, or f?= 0. (5.9b)
By applying the boundary conditions:
.1
2
)(
2
2
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
==
h
z
dr
dPh
r
zf
V
r
?
(5.10)
Denoting the constant in Eq. (5.7b) by K, we have:
,
r
K
dr
dP
= (5.11)
integration of which leads to:
,ln
?
?
?
?
?
?
?
?
=?=?
i
o
io
r
r
KPPP (5.12)
with subscripts i and o denoting the inner and outer radial positions, respectively. Then:
,
)ln(
i
o
r
r
r
P
dr
dP ?
= (5.13)
.1
)ln(
1
2
)(
2
2
?
?
?
?
?
?
?
?
??
?
?
?
?
??
==
h
z
r
r
P
r
h
r
zf
V
i
o
r
?
(5.14)
The volumetric flowrate expression is determined at a given radial location as:
)ln(3
4
3
4
)2(
33
o
i
h
h
r
r
r
Ph
dr
dP
r
h
dzrVQ
?
?
?
?
?
?
===
?
+
?
. (5.15)
As expected, the volumetric flow rate varies linearly with the differential static pressure.
Also note that since r
o
> r
i
, the static pressure rises radially inward
139
5.3.2 Determination of the Friction Factor
In addition, the friction factor for radial Stokes flow between two concentric disks
is determined. The wall shear stress, ?
rz
, is computed from the wall velocity gradient:
dr
dP
h
dz
dV
hz
r
rz
==
+=
??
. (5.16)
An expression for
dr
dP
in terms of average velocity is needed. Dividing Eq. 5.15 by
cross-sectional area, A = 2?r(2h), one gets the average velocity:
.
3)2(2
2
dr
dPh
hr
Q
V
??
==
(5.17)
Solving (Eq. 5.16) for
dr
dP
, one gets:
.
3
2
h
V
dr
dP ?
=
(5.18)
Substituting for
dr
dP
in Eq. 5.16, one gets:
.
3
h
V
rz
?
? =
(5.19)
Nondimensionalizing by the dynamic pressure, one gets the Darcy friction factor, f:
.
244
2
2
1
VhV
f
?
?
?
?
==
(5.20)
Taking note of the hydraulic diameter, D
h
,
,4
)2(2
)22(44
sec
h
r
hr
perimeter
A
D
tionalcross
h
===
?
?
?
(5.21)
140
Eq. 5.20 then becomes,
.
Re
96
h
D
f =
(5.22)
5.3.3 Determination of the Total Pressure Loss
The total or stagnation pressure is defined as the pressure attained if the flow was
brought isentropically to rest at that point. Furthermore, the total pressure is the sum of
the static pressure and dynamic pressure. Therefore, the equation for change in total
pressure between the outer radial position, denoted ?o? and the inner radial position,
denoted ?i?, is:
()( )
()()
().VVP
VVPP
VPVPPP
2
i
2
o
2
1
2
i
2
o
2
1
io
2
i
2
1
i
2
o
2
1
ott
io
?+?=?
?+?=
+?+=?
?
?
??
t
P
(5.23)
Substituting for the average velocities from equation 5.17, one gets:
.
11
32
PPP
2222
2
tt
io
?
?
?
?
?
?
?
?
?+?=?
io
rr
h
Q
?
?
(5.24)
Substituting for ?P from equation 5.15 gives:
.
11
8
ln
3
4
1
PPP
222
ttt
io
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?+
?
?
?
?
?
?
?
?
=?=?
ioo
i
rr
Q
r
r
hh
Q
?
??
?
(5.25)
141
142
The differential total pressure varies quadratically with Q. Equation (5.25) provides a
convenient analytic expression for the total pressure loss for a limiting case of parallel
disks with no rings.
5.4 Results and Discussion
Two different views of the pathlines of fluid flow (tracked by buoyant particles
that follow the flow) are shown in Figure 5.12, for the case of 2.7 ?m deflection and a
total pressure drop of 2.4 psi. The monitored buoyant particles are released at the inlet
plane from two rake sets that are placed normal to each other (Figure 5.13). The
monitored pathlines released from these rakes are colored by the local total pressure
values with the scale shown on the left side of both figures. The fluid particles that are
shown in Figure 5.12 prefer to flow over the tether rather than beneath it due to lower
resistance. Total pressure is almost constant within the inlet pipe (Figure 5.13) but drops
suddenly and markedly where the flow impacts the valve seat plate that is located 10 ?m
away from the inlet port. Once the flow emerges into the microvalve cavity, the total
pressure remains almost constant, due to the very small velocity in such a relatively large
volume. This uniformity of total pressure is maintained until the flow passes over the
rings in the gap between the seat and the boss plates. As it is shown in Figure 5.14, the
total pressure monotonically decreases as the flow passes over each ring. Passage of the
flow through this narrow gap causes most of the total pressure drop in the system. There
was also a lesser total pressure drop associated with the small box under the inlet port
(Figure 5.13). Since most of the total pressure drop occurs as the fluid goes over the
rings, especially for small deflections, a 2-D axisymmetric CFD analysis was also
143
performed, taking into account the geometry of the rings. This was deemed necessary in
view of the relatively small number of the grid points used to resolve the spacing between
the boss plate and the rings in the fully 3-D model. The total pressure drop over the rings
determined from the 2-D analysis verified the accuracy of the 3-D microvalve model
results, which are consistently lower than the 2-D predictions.
The mass flow rates calculated from the 2-dimensional Stokes, 2-dimensional
CFD and 3-dimensional CFD analyses are compared to the experimental results for 1
micron deflection (Table 5.2). In the first column, the differences between the inlet and
outlet total pressures are reported. The Reynolds number of the outlet pipe changes from
1.27 to 2.86. The Stokes analysis predicts high mass flow rates compared to others
because the effect of the rings is not considered in this analysis. The experimental results
are greater than 2-dimensional CFD and lower than 3-dimensional CFD predictions.
Table 5.2 Comparison of the mass flow rates (mg/s) obtained from different analyses for
deflection = 1 ?m
Total pressure drop for the
rings
Total pressure drop for
system
?P
t
[psig]
Re
(Outlet)
Stokes Flow 2-D Model 3-D Model Experimental
5 1.27 18.3 0.260 0.0876 0.2
15 2.54 55.1 0.781 0.263 0.4
20 2.86 73.5 1.04 0.352 0.45
Total pressure drop versus mass flowrate results are compared between the
numerical (2-D and 3-D) and experimental data for a 1-?m deflection in Figure 5.15.
The observed trend for the numerical results is a straight line that passes through the
origin. The experimental results for a deflection of 1 micron lie between the 2-D and 3-D
predictions. The experimental data are close to 2-D results for low inlet pressure cases
but for the cases with high inlet pressures, the experimental data tend toward 3-D results.
A similar comparison is demonstrated in Figure 5.16 for a deflection of 3.5 ?m. In
contrast to the 1- ?m deflection case, the experimental results are very different from the
numerical results in this case, although the difference between 2-D and 3-D numerical
simulations is still reasonable. In Figure 5.17, the 3-D numerical predictions for the total
pressure drop between the inlet and outlet ports versus mass flowrate are given for
deflections of 1, 2.7 and 3.5 ?m. For a constant total pressure difference, the mass flow
rate increases for greater values of deflection.
The values of the total pressure drop coefficient, K versus deflection are shown in
Figure 5.18. The total pressure drop coefficient is defined as the slope of each line in
Figure 5.17:
.
sec
?
?
?
?
?
?
?
?
?
?
?
=
?
mg
psi
m
P
K
t
(5.26)
Since the total pressure drop changes linearly with respect to change in mass
flowrate, K is a constant number. The value of K is exponentially dependent on
deflection, and is therefore extremely sensitive to changes in deflection. The empirical
dependence obtained by plotting a trendline through the points in Figure 5.18 is:
144
,29.188
3.1 deft
e
m
P
K
??
?
=
?
=
(5.27)
where def is the value of deflection in microns,
t
P? is the total pressure drop (the
difference between the inlet and outlet total pressures) in psi and is the mass flowrate
of the system in
?
m
sec
mg
. The exponential trend exhibited in Figure 5.18 can explain the
large difference between the experimental and numerical results for high deflection cases.
The other reason could be the error in determining the deflection precisely. The
measured deflections of the PZT stack actuator did not incorporate any changes in
deflection during liquid flow, which could give rise to approximately 1-?m error in the
value of the deflection. The measured values of the deflection are shown in Figure 5.19.
The measurements are done at atmospheric pressure. At zero voltage, the boss plate is
sitting on the rings (zero deflection). The boss plate moves up as the applied voltage
increases. The upper line in Figure 5.19 demonstrates the measured deflection while the
voltage is increasing. The deflection is about 4.5 microns when the voltage is 50 V that
is the maximum designate deflection. The lower line represents the deflection values for
the same voltages as the voltage decreases from 50 V to zero while the boss plate moves
down toward the rings. Deflection uncertainty of the order of 1 micron is observed when
the applied voltage to the piezoelectric actuator is 30 V. This uncertainty in the value of
deflection can lead to marked changes in the K value.
145
146
5.5 Conclusions
1. The predicted mass flow rate in the microvalve varies linearly with the imposed
total pressure difference. This is expected because the flow is in the low range of
Reynolds numbers with the highest value being 6.
2. The results indicate the sensitivity of the total pressure drop coefficient to
deflection. Small change in deflection can greatly modify the value of K. Therefore, the
differences between the numerical and experimental results can be explained in view of
the error tolerance of the measured deflections.
5.6 Closure
In this Chapter, a CFD-based simulation of flow in a complicated 3-dimensional
NASA JPL microvalve was presented. A 2-dimensional CFD analysis was also carried
out for the flow between the boss and seat plates. The CFD-based results are compared
to a simplified 2-D analysis of Stokes flow between to disks and the experimental
observations.
147
Inlet Hole
Seat Ring
Outlet Hole
Lower-Boss
Upper-Boss
PZT Actuator
Seat
Figure 5.1 Structure of the liquid-compatible microvalve (Yang et al., 2005)
148
Outlet
Boss Plate
Inlet Inlet
Rings
Figure 5.2 Operating principle of the microvalve in actuated ?on? state Yang et al.,
2005)
(a)
Inlet Holes
Outlet Hole
(b)
Figure 5.3 SEM images of the (a) valve seat and (b) upper boss plate (Yang et al., 2005)
149
Center Plate
Tethers
(c)
Inlet Holes
(d)
Figure 5.3 SEM images of the (c) valve sealing surface of the lower boss plate
and (d) upper boss bonding side of the lower boss plate (Yang et al., 2005)
150
Figure 5.4 Fully assembled and packaged liquid-compatible microvalve (Yang et
al., 2005)
151
Figure 5.5 One-quarter of the microvalve structure
152
Figure 5.6 Isometric view of the microvalve geometry used for the computer modeling
153
154
Outlet
Inlet
Figure 5.7 Isometric view of the overall 3-dimensional computational mesh
10-micron-high Box
Inlet
Figure 5.8 Isometric view of the 3-dimensional computational mesh near the inlet pipe
and 10-micron-high box
155
Rings
Outlet
Figure 5.9 Isometric view of the 3-dimensional computational mesh showing the detail of
one-quarter of the valve seat rings
156
Boss Plate
Deflection (1 ~ 3.5 microns)
157
Figure 5.10 Isometric view of one of the rings
1 micron
10 microns
Ring
Seat Plate
Figure 5.11 Geometry of the radial flow between two parallel disks
158
(a)
(b)
Figure 5.12 Isometric views showing pathlines of selected particles released at the inlet
plane for a deflection 2.7 ?m and total pressure drop of 2.4 psi
159
Figure 5.13 Streamlines colored by stagnation pressure in the inlet pipe and 10-micron
box for a deflection 2.7 ?m and total pressure drop of 2.4 psi
160
Figure 5.14 Streamlines colored by stagnation pressure over the rings and the outlet hole
for a deflection 2.7 ?m and total pressure drop of 2.4 psi
161
Mass Flow Rate [mg/sec]
?
P
t
[p
s
i
g
]
0 0.1 0.2 0.3 0.4 0.5
0
4
8
12
16
20
24
3-D (Microvalve)
Experimental
2-D (Rings)
Deflection = 1 micron
Figure 5.15 Comparison between two and three dimensional numerical results and
experimental data (Yang et al., 2005) with 1 micron deflection
162
Mass Flo
163
w
rate [mg/sec]
?
P
t
[p
s
i
g
]
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
2
4
6
8
10
12
14
16
18
20
3-D (Micro alve)
Experimental
2-D (Ring
v
s)
Deflection = 3.5 mic onr
Figure 5.16 Comparison between two and three dimensional numerical results and
experimental data (Yang et al,. 2005) with 3.5 micron deflection
Mass Flow Rate [mg/sec]
?
P
t
[p
s
i
g
]
0.25 0.5 0.75 1
4
8
12
16
20
24
1.0
2.7
3.5
Deflection [micron]
Figure 5.17 Comparison among three-dimensional numerical results with different
deflections
164
3.5
1
2.7
1
10
100
Deflect
165
io
n [microns]
K
[
psi
/
(
m
g
/
sec)
]
def
eK
??
=
3.1
29.188
Figure 5.18 Total pressure drop coefficient for different deflections
Voltage [ V ]
D
e
fl
e
c
ti
o
n
[
m
i
c
r
o
n
s
]
0 10203040
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Forward Measurment
Backward Measurment
50
Figure 5.19 The measured deflections as a function of applied voltage (Yang et al., 2005)
166
168
CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS
In order to investigate forced convection and fluid flow in enclosures with inlet
and outlet ports, a high-accuracy CFD software is developed and tested. The steady-state
and transient flow and heat transport phenomena in an enclosure with constant and
oscillating inlet velocities have been studied in great detail. The main conclusions and
achievements are summarized as follows:
1. A general-purpose primitive-variable-based finite-volume program was
developed. It is based on the SIMPLE algorithm. The QUICK scheme was used for the
convective terms and the central difference scheme was used for the diffusive terms. In
order to capture the transient behaviors for the system with an oscillating inlet velocity,
the Three-Time-Level method was used for the unsteady terms. The benchmarking of the
code for the applicable limiting cases was very successful.
2. For steady forced convection in a square cavity with inlet and outlet ports, it
was found that to maximize the heat transfer rate and minimize the pressure drop, the best
location for the outlet port is on the bottom of the right wall when the inlet port is placed
on top of the left wall. For design purposes, inlet and outlet ports should be parallel to
each other in order to minimize the pressure drop, even if they are designed to be placed
on the same wall. In order to realize excessive pressure drops, it is recommended that
they are perpendicular to each other. By placing the outlet port with one end at three
corners, maximum total Nusselt number of the cavity can be achieved. Minimum total
168
heat transfer of the cavity is achieved with the outlet port located at the middle of the
walls.
3. For a square cavity with an oscillating inlet velocity, both the flow and thermal
fields reach their respective periodic states after certain periods of time. The flow fields
reach their periodic state earlier than the thermal fields. Regardless of the Strouhal
number, heat transfer enhancement in comparison to the performance of the system for
the mean Reynolds number is consistently observed. The best heat transfer rate is
realized when the Strouhal number is close to unity. This is the case when the period of
the incoming stream resonates with the convection time scale.
4. In the microvalve flow analysis, it is found that the mass flow rate varies
linearly with the total pressure difference because of the low Reynolds numbers of the
system. The results also indicate the sensitivity of the total pressure drop coefficient to
deflection. Small change in deflection can greatly modify the total pressure drop
coefficient. The differences between the numerical and experimental results can be
explained in view of the error tolerance of the measured deflections.
The results of this study are useful to design of systems with enclosure that have
inlet and outlet ports, such as heat exchangers, thermal energy storage and microfluidic
MEMS devices. The MEMS systems that can benefit from this work are microvalves
and micropumps. The suggestions for future work are summarized as follows:
1. Flexible boundaries: Many micropump systems have flexible walls that move
with a piezoelectric actuator. It will be very helpful to model the flow to investigate the
pressure drop in many novel micropump designs.
169
2. Rotating objects: By placing a rotating object in the cavity, the code will be
capable to model flow and heat transfer in enclosures with a moving object. This has
applications in microinjection systems, industrial mixers and heat exchangers.
3. Multi inlet/outlet ports: This research has focused on cavities with one inlet
and one outlet port. Adding more inlet/outlet ports will increase the capability of the
code to model fluid flow and heat transfer in more complicated geometries.
170
REFERENCES
Abu-Mulaweh, H. I., 2003, ?A Review of Research on Laminar Mixed Convection
Flow Over Backward- and Forward-Facing Steps,? Int. J. Thermal Sciences, Vol. 42, No.
9, pp. 897-909.
Abu-Mulaweh, H. I., B. F. Armaly and T. S. Chen, 1995, ?Laminar Natural
Convection Flow Over a Vertical Backward-Facing Step,? J. Heat Transfer, Vol. 117,
No. 4, pp. 895-901.
Abu-Mulaweh, H. I. and T. S. Chen, 1993, ?Measurements of Laminar Mixed
Convection Flow Over a Horizontal Forward-Facing Step,? J. Thermophysics and Heat
Transfer, Vol. 7, No. 4, pp. 569-573.
Armaly, B. F., F. Durst, J. C. Pereira, and B. Schoenung, 1983, ?Experimental and
Theoretical Investigation of Backward-Facing Step Flow,? J. of Fluid Mechanics, Vol.
127, pp. 473-496.
Bird, R. B., W. E. Stewart and E. N. Lightfoot, 1960, "Transport Phenomena," John
Wiley and Sons, Inc., NY.
Bouhdjar, A. and A. Harhad, 2002, ?Numerical Analysis of Transient Mixed
Convection Flow in Storage Tank: Influence of Fluid Properties and Aspect Ratios on
Stratification,? Renewable Energy, Vol. 25, pp. 555-567.
171
Chan, A. M. C., P. S. Smereka and D. Giutsi, 1983, ?A Numerical Study of Transient
Mixed Convection Flows in a Thermal Storage Tank,? J. Solar Energy Engineering, Vol.
105, No. 3, pp. 246-253.
Darbandi, M., G. E. Schneider, 1998, ?Numerical Study of the Flow Behavior in the
Uniform Velocity Entry Flow Problem,? Numerical Heat Transfer; Part A: Applications
Vol. 34, No. 5, pp. 479-494.
Esashi, M., 1990, ?Integratcd Microflow Control Systems,? Sensors and Actuators,
A21-A23, pp. 161-168.
Fang, L. C., D. Nicolaou and J. W. Cleaver, 1999, ?Transient Removal of a
Contaminated Fluid from a Cavity,? Int. J. Heat Fluid Flow, Vol. 20, No. 6, pp. 605-613.
Ferziger, J. H. and Peri?, M., 1999, ?Computational Methods for Fluid Dynamics,?
2
nd
Edition, Springer-Verlag.
Freitas, C. J., 1995, ?Perspective: Selected Benchmarks from Commercial CFD
Codes,? J. Fluids Engineering, Vol. 117, pp. 208-218.
Ghia, U., Ghia, K. N. and Shin, C. T., 1982, "High-Re Solutions for Incompressible
Flow Using the Navier-Stokes Equations and Multigrid Method," J. Comput. Phys., Vol.
48, pp. 387-411.
Greiner, M., 1991, ?An Experimental Investigation of Resonant Heat Transfer
Enhancement in Grooved Channels,? Int. J. Heat Mass Transfer, Vol. 34, No. 6, pp.
1383?1391.
Han, T., J. A. C. Humphrey and Launder, B. E., 1982, ? A Comparison of Hybrid and
Quadratic-Upstream Differencing in High Reynolds Number Elliptic Flows,? Comp.
Methods Appli. Mech. Eng., Vol. 29, pp. 81-95.
172
Harlow, F. H. and J. E. Welch, 1965, ?Numerical Calculation of Time-Dependent
Viscous Incompressible Flow of Fluid with Free Surface,? Phys. Fluids, Vol. 8, pp. 2182-
2189.
Hayase, T., J. A. C. Humphrey and R. Grief, 1992, ?A Consistently Formulated
QUICK Scheme for Fast and Stable Convergence Using Finite-Volume Iterative
Calculation Procedures,? J. Comp. Phys., Vol. 98, pp. 108-118.
Homan, K. O. and S. L. Soo, 1998, ?The Steady Horizontal Flow of a Wall Jet into a
Large-Width Cavity,? J. Fluids Engineering, Vol. 120, No. 1, pp. 70-75.
Hsu, T-H, P-T Hsu and S-P How, 1997, ?Mixed Convection in a Partially Divided
Rectangular Enclosure,? Numerical Heat Transfer, Part A, Vol. 31, pp. 655-683.
Hsu, T. H. and S. G. Wang, 2000, ?Mixed Convection in a Rectangular Enclosure
with Discrete Heat Sources,? Numerical Heat Transfer, Part A, Vol. 38, pp. 627-652.
Hu S-C and Tsao J-M, National Taipei University of Technology, TAIWAN, private
communication.
Jayaram, S. and D. Thangaraj, 1995, ?Effects of Eddies on Convected Space Charge
Density in Relaxation Tanks,? in Proceedings of the Electrical/Electronics Insulation
Conference, IEEE, pp. 123-127.
Johnson, C. A., 2005, ?Modeling of Gas Flow in a Piezoelectrically Actuated High-
pressure Microvalve for Flowrate Control,? Masters Thesis, Mech. Eng. Dept., Auburn
University, AL, 2005.
Karniadakis, G. E. and A. Beskok, 2001, ?Micro Flows: Fundamentals and
Simulation,? Springer, New York, NY.
173
Lee, C. and E-H Yang, 2004, ?Piezoelectric Liquid-Compatible Microvalve for
Integrated Micropropulsion,? Solid State Sensor, Actuator and Microsystems Workshop,
Hilton Head Island, South Carolina.
Leonard, B. P., 1979, ?A Stable and Accurate Convective Modelling Procedure Based
on Quadratic Upstream Interpolation,? Comput. Methods Appl. Mech. Eng., Vol. 19, pp.
59-98, 1979.
Leschziner, M. A., 1980, ?Practical Evaluation of Three Finite Difference Schemes
for the Computation of Steady-State Recirculating Flows,? Comput. Methods Appl. Mech.
Eng., Vol. 23, pp. 293-312.
Leung Ki, Y-S, D. Maillefer, G. Rey-Mermet, P. A. Monkewitz, P. Renaud and H. T.
G. van Lintel, 1998, ?Microfluidics of a Planar Microfabricated Fluid Filter,? in Micro-
ElectroMechanical Systems (MEMS) ? 1998, DSC-Vol. 66, ASME, pp. 165-170.
Mei, R. W. and A. Plotkin, 1986, ?Navier-Stokes Solution for Laminar
Incompressible Flows in Forward-Facing Step Geometries,? AIAA J., Vol. 24, No. 7, pp.
1106-1111.
Nakagawa, S., S. Shoji and M. Esashi, 1990, ?A Microchemical Analyzing System
Integrated on a Silicon Wafer,? Proc. IEEE-MCMS Workshop, pp. 89-94.
Nishimura, T., K. Kunitsugu and A. M. Morega, 1998, ?Fluid Mixing and Mass
Transfer Enhancement in Grooved Channels for Pulsatile Flow,? Enhanced Heat
Transfer, Vol. 5, pp. 23?37.
Occhialini, J. M. and J. J. L. Higdon, 1992, ?Convective Mass Transport from
Rectangular Cavities in Viscous Flow, ? J. Electrochemical Society, Vol. 139, No. 10, pp.
2845-2855.
174
Omri, A. and S. B. Nasrallah, 1999, ?Control Volume Finite Element Numerical
Simulation of Mixed Convection in an Air-Cooled Cavity,? Numerical Heat Transfer,
Part A, Vol. 36, No. 6, pp. 615-637.
Patankar, S. V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere Pub.
Co., Washington, DC.
Patera, A. T. and B. B. Mikic, 1986, ?Exploiting Hydrodynamic Instabilities:
Resonant Heat Transfer Enhancement,? Int. J. Heat Mass Transfer, Vol. 29, No. 8, pp.
1127?1138.
Pollard, A. and A. l. Siu, 1982, ?The Calculation of Some Laminar Flows Using
Various Discretization Schemes,? Compt. Methods Appl. Mech. Eng., Vol. 35, pp.293-
313.
Roache, P. J., 1998, ?Fundamental of Computational Fluid Dynamics,? Hermosa
Publishers, NM.
Rosengarten, G., M. Behnia and G. Morrison, 1999, ?Some Aspects Concerning
Modelling the Flow and Heat Transfer in Horizontal Mantle Heat Exchangers in Solar
Water Heaters,? Int. J. Energy Research, Vol. 23, pp. 1007-1016.
Saeidi, S. M. and J. M. Khodadadi, 2004, ?Forced Convection in a Square Cavity
with Inlet and Outlet Ports,? ASME Heat Transfer/Fluids Engineering Summer
Conference, Paper HT-FED2004-56246, Charlotte, NC.
Shi, X. and J. M. Khodadadi, 2002, ?Laminar Fluid Flow and Heat Transfer in a Lid-
Driven Cavity Due to a Thin Fin,? J. Heat Transfer, Vol. 124, No. 6, pp. 1056-1063.
175
Shi, X. and J. M. Khodadadi, 2003, ?Laminar Natural Convection Heat Transfer in a
Differentially Heated Square Cavity Due to a Thin Fin on the Hot Wall,? J. Heat
Transfer, Vol. 125, No. 4, pp. 612-623.
Shi, X. and J. M. Khodadadi, 2004, ?Fluid Flow and Heat Transfer in a Lid-driven
Cavity due to an Oscillating Thin Fin: Transient Behavior,? J. of Heat Transfer, Vol. 126,
No. 6, pp. 924-930.
Shi, X. and J. M. Khodadadi, 2005, ?Periodic State of Fluid Flow and Heat Transfer
in a Lid-driven Cavity Due to an Oscillating Thin Fin,? Accepted for Publication in the
Int. J. Heat Mass Transfer.
Shoji, S. and M. Esashi, 1988, ?Micromachining for Chemical Sensors,? Chem.
Sensor Techn., Vol. 1, pp. 179-272.
Shoji, S. and M. Esashi, 1994, ?Microflow Devices and Systems,? J. Micromech. and
Microeng. Vol. 4, pp. 157-171.
Shoji, S., B. H. van der Schoot, N. F. de Rooij, and M. Esashi, 1991, ?Smallest Dead
Volume Microvalves for Integrated Chemical Analyzing Systems,? Technical Digest of
Transducers, Vol. 91, pp. 1052-1057.
Shuja, S. Z., B. S. Yilbas and M. O. Budair, 2000a, ?Natural Convection in a Square
Cavity with a Heat Generating Body: Entropy Considerations,? Heat and Mass Transfer,
Vol. 36, pp. 343-350.
Shuja, S. Z., B. S. Yilbas and M. O. Iqbal, 2000b, ?Mixed Convection in a Square
Cavity due to Heat Generating Rectangular Body. Effect of Cavity Exit Port Locations,?
Int. J. Numerical Methods Heat Fluid Flow, Vol. 10, No. 8, pp. 824-841.
176
Singh, S. and M. A. R. Sharif, 2003, ?Mixed Convective Cooling of a Rectangular
Cavity with Inlet and Exit Openings on Differentially Heated Side Walls,? Numerical
Heat Transfer, Part A, Vol. 44, pp. 233-253.
Spalding, D. B., 1972, ?A Novel Finite-Deference Formulation for Differential
Expressions Involving Both First and Second Derivative,? Int. J. Numer. Methods Eng.,
Vol. 4, pp. 551-559.
St?er, H., A. Gyr and W. Kinzelbach, 1999, ?Laminar Separation on a Forward
Facing Step,? European J. Mechanics, B/Fluids, Vol. 18, No. 4, pp. 675-692.
Vandelli, N., D. Wroblewski, M. Velonis, and T. Bifano, 1998, ?Development of a
MEMS Microvalve Array for Fluid Flow Control,? J. of MicroElectoMechanical
Systems, Vol. 7, No. 4, pp. 395-403.
Tretheway, D. C. and C. D. Meinharta, 2002, ?Apparent Fluid Slip at Hydrophobic
Microchannel Walls,? Physics of Fluids, Vol. 14, No. 3, pp. L9-L12.
Tretheway, D. C. and C. D. Meinharta, 2004, ?A Generating Mechanism for
Apparent Fluid Slip in Hydrophobic Microchannels,? Physics of Fluids, No. 5, Vol. 16,
pp. 1509-1515.
Verboven, P., A. K. Datta, N. T. Anh, N. Scheerlinck and B. M. Nicolai, 2003,
?Computation of Airflow Effects on Heat and Mass Transfer in a Microwave Oven,? J.
Food Engineering, Vol. 59, pp. 181-190.
White, M. F., 1999, ?Fluid Mechanics,? McGraw-Hill, New York, NY.
Yang, E-H, 2004, JPL MEMS Technology Group, private communication.
177
Yang, E-H, C. Lee, J. Mueller and T. George, 2004, ?Leak-Tight Piezoelectric
Microvalve for High-Pressure Gas Micropropulsion,? J. of MicroElectroMechanical
Systems, Vol. 13, No. 5, pp. 799-807.
Yang, E-H, C. Lee, S. M. Saeidi and J. M. Khodadadi, 2005, ?Fabrication,
Characterization and Computational Modeling of a Piezoelectrically Actuated
Microvalve for Liquid Flow Control,? Under Review by J. MEMS.