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.