MODELING OF FRICTIONAL GAS FLOW IN A PIEZOELECTRICALLY ACTUATED HIGH-PRESSURE MICROVALVE FOR FLOWRATE CONTROL Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classified information. ____________________________________________ Christopher Alan Johnson Certificate of Approval: ______________________________ ______________________________ Daniel W. Mackowski Jay M. Khodadadi, Chair Associate Professor Professor Mechanical Engineering Mechanical Engineering ______________________________ ______________________________ Robert L. Jackson Stephen L. McFarland Assistant Professor Acting Dean Mechanical Engineering Graduate School MODELING OF FRICTIONAL GAS FLOW IN A PIEZOELECTRICALLY ACTUATED HIGH-PRESSURE MICROVALVE FOR FLOWRATE CONTROL Christopher Alan Johnson A Thesis Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Master of Science Auburn, Alabama December 16, 2005 iii MODELING OF FRICTIONAL GAS FLOW IN A PIEZOELECTRICALLY ACTUATED HIGH-PRESSURE MICROVALVE FOR FLOWRATE CONTROL Christopher Alan Johnson Permission is granted to Auburn University to make copies of this thesis 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 Christopher Alan Johnson, son of John Milton and Janet (Fourroux) Johnson, was born June 4, 1980, in Huntsville, Alabama. He graduated from Virgil I. Grissom High School in 1998. He attended Auburn University for four years and graduated magna cum laude with a Bachelor of Science degree in Mechanical Engineering in December, 2002. He then entered the Graduate School of Auburn University in January, 2003. After working as a graduate teaching assistant for two semesters, he was awarded an Alabama Space Grant Consortium Fellowship. v THESIS ABSTRACT MODELING OF FRICTIONAL GAS FLOW IN A PIEZOELECTRICALLY ACTUATED HIGH-PRESSURE MICROVALVE FOR FLOWRATE CONTROL Christopher Alan Johnson Master of Science, December 16, 2005 (B.S.M.E., Auburn University, 2002) 197 Typed Pages Directed by Jay M. Khodadadi One-dimensional modeling of frictional compressible gaseous flow through a high-pressure piezoelectrically-actuated microvalve is studied. Focusing on the micro- scale gap between the boss and seat plates, variations of flow properties were predicted using two 1-D models: 1. Channel Flow Model, 2. Radial Flow Model. Both models utilized a 4th order Runge-Kutta algorithm to integrate a respective system of nonlinear ordinary differential equations. A channel flow model was developed for steady, compressible flow of a perfect gas between two infinite insulated parallel plates. This model served the purpose of benchmarking the numerical code against analytical expressions for the properties of flow through a constant-area channel. Additionally, utilizing this model, the total vi pressure loss through the outlet tube was found to be negligible in comparison to that of total pressure loss across the seat rings. The radial flow model was developed for steady, axisymmetric, compressible flow of a perfect gas between two insulated, parallel disks flowing radially toward an outlet hole at the center of the bottom disk. This model was implemented to determine the variation of properties of flow between the boss and seat plates of a JPL microvalve. The most notable conclusion from the flow property trends is that of a drastic increase in density and static pressure in contrast to a rather small increase in the Mach number. Also of importance, the total pressure drop was shown to be significant across the seat rings. A 2-D Stokes flow model was derived for incompressible, axisymmetric, radial flow between two concentric parallel disks. The results of this model were used to verify the flow property variations from the radial flow model. In particular, for the Stokes flow model, relations for radial velocity, average velocity, Darcy friction factor, volumetric flowrate, static pressure rise, and total pressure drop were derived. The Stokes flow model trends for both static and total pressure were in accord with the radial compressible flow model trends. In addition, a comparison of Stokes flow values for both the static pressure rise and the total pressure drop to that of the numerical results demonstrates the necessity of accounting for compressibility effects. vii ACKNOWLEDGMENTS The author would like to express a sincere appreciation of his major professor, Dr. Jay M. Khodadadi, for his excellent guidance, enthusiasm, sound advice, encouragement, and patience throughout this study. Dr. Khodadadi was always accessible when difficulties occurred, and often gave his time at odd hours in order to offer advice and assistance. It is difficult to fully express my gratitude in words. Gratitude is also due to my committee members, Drs. Jackson, and Mackowski. I would also like to thank the professors of the Mechanical Engineering Department of Auburn University for their inspiration and inducing in me an abundant interest in the subject matter. I would also like to thank Dr. J-M Wersinger and the Alabama Space Grant Consortium (ASGC) for their support during my research. I am deeply grateful for the ASGC fellowship (NASA Training Grant NGT5-40077). A great debt of gratitude goes to Drs. E-H Yang, Choonsup Lee, Thomas George, and the remaining members of the MEMS Technology Group at JPL. I truly appreciate your patience, assistance with all my questions, and helpfulness throughout this research. Lastly, I would like to thank my fianc?e, Colleen Moynihan, who provided the greatest support and loving encouragement to help complete this research. I also wish to thank my parents for always inspiring me to do my best at everything. viii Style manual or journal used: Guide to Preparation and Submission of Theses and Dissertations 2005 Computer software used: MS Word 2003 ix TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . xv NOMENCLATURE . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Microvalves for Micropropulsion System Applications . . . . . . . 2 1.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . 1.4 Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . 4 CHAPTER 2 LITERATURE REVIEW OF MICROVALVES . . . . . . . . 2.1 Microvalve Categorization . . . . . . . . . . . . . . . . . . . 8 2.1.1 Thermopneumatic Actuation . . . . . . . . . . . . 2.1.2 Shape Memory Alloy (SMA) Actuation . . . . . . . . . . 10 2.1.3 Bi-morph Actuation . . . . . . . . . . . . . . . . . . 12 2.1.4 Electromagnetic Actuation . . . . . . . . . . . . . . . 14 . xxi . . 1 . . 3 . . 8 . . . 9 x 2.1.5 Electrostatic Actuation . . . . . . . . . . . . . . . . . 15 2.1.6 Piezoelectric Actuation . . . . . . . . . . . . . . . . . 16 2.2 Closure . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 3 MICROVALVE MODELING APPROACH . . . . . . . . . 3.1 Background and Microvalve Design Requirements . . . . . . . 3.1.1 Design of the Microvalve . . . . . . . . . . . . . . . . 35 3.2 2-Dimensional Axisymmetic Geometry Model . . . . . . . . . . 37 3.2.1 Sample Calculation for Incompressible Total Pressure Drop of the Inlet Piping Network . . . . . . . . . . . . . . . 39 3.2.2 Outline of the 2-Dimensional Axisymmetric Geometry Model . . . . . . . . . . . . . . . . . . . . . 3.3 1-Dimensional Axisymmetric Geometry Model . . . . . . . . . . 43 3.4 Closure . . . . . . . . . . . . . . . . . . . . . . . . . . 43 CHAPTER 4 MATHEMATICAL FORMULATION . . . . . . . . . . . . . 57 4.1 Channel Flow Model . . . . . . . . . . . . . . . . . . . . 57 4.1.1 Governing Equations . . . . . . . . . . . . . . . 4.2 Derivation of the Differential Equations for Flow Properties . . . . . 60 4.2.1 Derivation of the Differential Equation for the Mach number (M) . . . . . . . . . . . . . . . . . . . . . 17 . . 33 . . 33 . . 42 . . 58 . . 61 xi 4.2.2 Derivation of the Differential Equation for Pressure (P) . . 4.2.3 Derivation of the Differential Equation for Velocity (V) . . . . 65 4.2.4 Derivation of the Differential Equation for Density (?) . . 4.2.5 Derivation of the Differential Equation for Stagnation Pressure (P t ) . . . . . . . . . . . . . . . . . . . . . 66 4.2.6 Derivation of the Differential Equation for Entropy Change (ds) . . . . . . . . . . . . . . . . . . . . . . 4.3 Dimensionless Forms of the Differential Equations . . . . . . . . . 70 4.4 Radial Disk Flow Model . . . . . . . . . . . . . . . . . . . . 73 4.4.1 Governing Equations . . . . . . . . . . . . . . . . . 73 4.5 Derivation of Differential Expressions for Flow Properties . . . . . . 74 4.5.1 Derivation of the Differential Equation for the Mach number (M) . . . . . . . . . . . . . . . . . . . . . 75 4.5.2 Differential Equations for Temperature (T), Stagnation Pressure (P t ) and Entropy Change (ds) . . . . . . . . . . . 77 4.5.3 Derivation of the Differential Equation for Density (?) . . . . 78 4.5.4 Derivation of the Differential Equation for Velocity (V) . . . . 78 4.5.5 Derivation of the Differential Equation for Pressure (P) . . . . 79 4.6 Conversion of Differential Equations in Terms of Radial Coordinate (r) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.7 Dimensionless Form of the Differential Equations . . . . . . . . . . 81 4.8 Closure . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 . . 66 . . 69 . . 83 xii CHAPTER 5 PREDICTIONS OF THE FLOW PROPERTIES UTILIZING THE MODELS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.1 Numerical Analysis . . . . . . . . . . . . . . . . . . . . . 92 5.1.1 Initial Conditions . . . . . . . . . . . . . . . . . . . 92 5.2 Benchmarking of the Computer Code for the Channel Flow Model . . 94 5.2.1 Variable Step Size . . . . . . . . . . . . . . . . . . 96 5.3 Radial Flow Model Results and Discussions . . . . . . . . . . . 97 5.3.1 Variation of the Reynolds Number (Re Dh ) . . . . . . . . . 97 5.3.2 Variation of the Change in Entropy (ds) . . . . . . . . . . 98 5.3.3 Variation of the Mach Number (M + ) and Velocity (V + ) . . . . 98 5.3.4 Variation of the Temperature (T + ) . . . . . . . . . . 5.3.5 Variation of the Total Pressure (P + t ) . . . . . . . . . . . 100 5.3.6 Variation of the Static Pressure (P + ) . . . . . . . . . . . 102 5.3.7 Variation of the Density (? + ) . . . . . . . . . . . . 5.3.8 Variation of the Knudsen Number (Kn) . . . . . . . . . . 104 5.4 Analysis of the Total Pressure Losses Beyond the Seat Rings . . . . 105 5.4.1 Sample Calculation for Incompressible Total Pressure Drop of the Outlet Tube . . . . . . . . . . . . . . . . . . 106 5.5 Closure . . . . . . . . . . . . . . . . . . . . . . . . . 108 . . 99 . . 103 xiii CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS . . . . . . . . . 139 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . 139 6.2 Recommendations for Future Work . . . . . . . . . . . . . . 141 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 A. Correlation for the Modified Friction Coefficient . . . . . . . . . . 147 B. Radial Stokes Flow Between Two Concentric Parallel Disks . . . . . . 155 C. Sample Calculation for the Initial Flow Conditions . . . . . . . . . . 161 D. Reynolds Number Relation in a Radial Disk . . . . . . . . . . . . 165 E. Derivation of the Analytical Expression for L max . . . . . . . . . . 168 F. Derivation of Static Pressure Relations . . . . . . . . . . . . . . 171 xiv LIST OF TABLES Table 2.1 Characteristics of Microvalve Actuation Methods (Ayhan, 2002) . . . 18 Table 2.2 Evaluation of Microvalve Actuation for Microspacecraft Applications (Mueller, 2000) . . . . . . . . . . . . . . . . . 18 Table 3.1 Microvalve Requirements for NASA?s Miniature Spacecraft Propulsion Needs, Compared with Reported Performance to Date (Yang et al., 2004) . . . . . . . . . . . . . . . . . . . Table 3.2 Total Pressure Drop (?P t ) for the Inlet Piping Network . . . . . . . 39 Table 5.1 Measured (Yang et al., 2004) and extrapolated initial flow conditions . . . . . . . . . . . . . . . . . . . . . . . . . 94 Table 5.2 Total Pressure Drop (?P t ) for the Radial Flow Model Compared to the Incompressible Stokes Model . . . . . . . . . . . . . . . 101 Table 5.3 Static Pressure Rise (?P = P o -P i ) for the Radial Flow Model Compared to the Incompressible Stokes Model . . . . . . . . . . 103 Table 5.4 Total Pressure Drop (?P t ) Through the Outlet Tube (Incompressible and Compressible Fanno Analyses) . . . . . . . . 106 Table A.1 Summary of Coefficients for the Modified Friction Correlation . . . 149 . . 34 xv LIST OF FIGURES Figure 1.1 A diagram displaying a range of satellite sizes (source: http://www.sstl.co.uk/) . . . . . . . . . . . . . . . . . . Figure 1.2 Concept drawing of a microsatellite (Janson et al., 1998) . . . . Figure 1.3 A microsatellite developed by the SpaceQuest Ltd. Company in a low Earth orbit (source: http://www.spacequest.com) . . . . . . . . 7 Figure 2.1 Cross-section of a thermopneumatic microvalve (Rich & Wise, 2003) . . . . . . . . . . . . . . . . . . . . . . . . . . 19 Figure 2.2 Cross-section of a Thermopneumatic microvalve (Yang et al., 1997) . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.3 Cross-section of a SMA microvalve , a) closed position, b) open position (Geon et al., 2000) . . . . . . . . . . . . . . . . Figure 2.4 Schematic cross section of the SMA microvalve (A, B, C denote the valve ports) (Skrobanek et al., 1997) . . . . . . . . . . . . . 22 Figure 2.5 Schematic the SMA microvalve cross-section during open operation (Skrobanek et al., 1997) . . . . . . . . . . . . . . . 23 Figure 2.6 Cross-section of a bimetallically actuated microvalve (Jerman et al., 1990) . . . . . . . . . . . . . . . . . . . . . . . . . 24 . . 5 . . 6 . 20 . . 21 xvi Figure 2.7 Cross-section of a 3-way bimetallically actuated microvalve (Messner et al., 1998) . . . . . . . . . . . . . . . . . . . . 25 Figure 2.8 Cross-section of a bimetallic actuator of a 3-way microvalve (Messner et al., 1998) . . . . . . . . . . . . . . . . . . Figure 2.9 Cross-section of a electromagnetically-actuated microvalve (Shinozawa et al., 1997) . . . . . . . . . . . . . . . . . . . 27 Figure 2.10 Schematic of the electromagnetic microvalve (Bintoro & Hesketh, 2005) . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.11 Schematic of the electromagnetic microvalve?s operation (Bintoro & Hesketh, 2005) . . . . . . . . . . . . . . . . . . . . . . 29 Figure 2.12 Schematic of a electrostatically-actuated microvalve (Ohnstein et al., 1990) . . . . . . . . . . . . . . . . . . . . . . . . . 30 Figure 2.13 Schematic of a electrostatically-actuated microvalve (Huff et al., 1990) . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.14 Schematic of a piezoelectrically-actuated microvalve (Roberts et al., 2003) . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 3.1 Schematic diagram of the leak-tight piezoelectric microvalve (Yang et al., 2004) . . . . . . . . . . . . . . . . . . . Figure 3.2 Cross-sectional microvalve configuration, showing the suspension of the boss plate by tensile stressed silicon tethers over the extended valve seat (Yang et al., 2004) . . . . . . . . . . . . . 26 . 28 . 31 . . 45 . . 46 xvii Figure 3.3 Operation principle of the microvalve: (a) Microvalve closed (cross-section A-A), and (b) Microvalve open (cross-section B-B) (Yang et al., 2004) . . . . . . . . . . . . . . . . . . . Figure 3.4 Scanning electron micrographs (SEM) of the seat plate and the boss plate (Yang et al., 2004) . . . . . . . . . . . . . . . Figure 3.5 Packaged high-pressure piezoelectric microvalve (Yang et al., 2004) . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 Figure 3.6 Isometric view of the microvalve enclosed in the housing . . . . . . 50 Figure 3.7 Exploded isometric view of the microvalve enclosed in the housing . . 51 Figure 3.8 Schematic diagram of the microvalve housing and geometry showing the main features (not drawn to scale) . . . . . . . . . . 52 Figure 3.9 Schematic diagram of the piping network leading to the microvalve cavity (not drawn to scale) . . . . . . . . . . . . . . 53 Figure 3.10 Isometric view of the radial flow problem with arrows indicating the directions of the incoming and outgoing flow streams . . . . . . 54 Figure 3.11 Top view of the radial flow problem with the boss plate not shown . . 55 Figure 3.12 Geometry of the 1-D Axisymmetric Model . . . . . . . . . . . . 56 Figure 4.1 The 2-D channel model of the microvalve geometry . . . . . . . . 84 Figure 4.2 The 1-D channel model of the microvalve geometry . . . . . . . . 85 Figure 4.3 Control volume for the changes in flow properties of the 1-D channel model . . . . . . . . . . . . . . . . . . . . . . . 86 Figure 4.4 Forces acting on the control volume of the 1-D channel model . . . . 47 . . 48 . . 87 xviii Figure 4.5 Schematic diagram of the radial flow model displaying radial locations, and the directions of x and r . . . . . . . . . . . . . 88 Figure 4.6 Control volume for the changes in flow properties of the 1-D radial model for (a) side view and (b) top view. . . . . . . . . . . . . 89 Figure 4.7 Forces acting on the control volume of the 1-D radial model . . . . . 90 Figure 5.1 Measured flow rates for an actuated microvalve (Yang et al., 2004) Figure 5.2 Variation of ReDh for an inlet total pressure of 100 psig . . . . . . 110 Figure 5.3 Variation of ReDh for an inlet total pressure of 200 psig . . . . . . 111 Figure 5.4 Variation of ReDh for an inlet total pressure of 300 psig . . . . . . 112 Figure 5.5 Variation of the change in entropy for an inlet total pressure of 100 psig . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 Figure 5.6 Variation of the change in entropy for an inlet total pressure of 200 psig . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Figure 5.7 Variation of the change in entropy for an inlet total pressure of 300 psig . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Figure 5.8 Variation of the Mach number for an inlet total pressure of 100 psig . . . . . . . . . . . . . . . . . . . . . . . . . Figure 5.9 Variation of the Mach number for an inlet total pressure of 200 psig . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 Figure 5.10 Variation of the Mach number for an inlet total pressure of 300 psig . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Figure 5.11 Variation of the velocity for an inlet total pressure of 100 psig . . . . 119 Figure 5.12 Variation of the velocity for an inlet total pressure of 200 psig . . . . 109 . . 116 . . 120 xix Figure 5.13 Variation of the velocity for an inlet total pressure of 300 psig . . . . 121 Figure 5.14 Variation of the temperature for an inlet total pressure of 100 psig . . 122 Figure 5.15 Variation of the temperature for an inlet total pressure of 200 psig . Figure 5.16 Variation of the temperature for an inlet total pressure of 300 psig . Figure 5.17 Variation of the total pressure for an inlet total pressure of 100 psig . . 125 Figure 5.18 Variation of the total pressure for an inlet total pressure of 200 psig . . 126 Figure 5.19 Variation of the total pressure for an inlet total pressure of 300 psig . . 127 Figure 5.20 Comparison of total pressure drop (numerical results and Stokes relation) vs. mass flowrate . . . . . . . . . . . . . . . . Figure 5.21 Variation of the static pressure for an inlet total pressure of 100 psig . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 Figure 5.22 Variation of the static pressure for an inlet total pressure of 200 psig . . . . . . . . . . . . . . . . . . . . . . . . . Figure 5.23 Variation of the static pressure for an inlet total pressure of 300 psig . . . . . . . . . . . . . . . . . . . . . . . . . Figure 5.24 Variations of the static pressure for the Stokes flow and ideal compressible flow relations for an inlet total pressure of 300 psig . Figure 5.25 Variation of the density for an inlet total pressure of 100 psig . . . . 133 Figure 5.26 Variation of the density for an inlet total pressure of 200 psig . . . . 134 Figure 5.27 Variation of the density for an inlet total pressure of 300 psig . . . . 135 Figure 5.28 Variation of the Knudsen number for an inlet total pressure of 100 psig . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 . 123 . 124 . . 128 . . 130 . . 131 . 132 xx Figure 5.29 Variation of the Knudsen number for an inlet total pressure of 200 psig . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 Figure 5.30 Variation of the Knudsen number for an inlet total pressure of 300 psig . . . . . . . . . . . . . . . . . . . . . . . . . Figure A.1 Geometry for the study of Luy et al. (1991) . . . . . . . . . . . 151 Figure A.2 Sample computational data from Luy et al. (1991) . . . . . . . . 152 Figure A.3 LSF for the five-data-point set with Re Dh = 10 and D = 2 . . . . . . 153 Figure A.4 Correlated approximation in comparison with data of Luy et al. (1991) . . . . . . . . . . . . . . . . . . . . . . . . . . 154 Figure B.1 Geometry of the radial flow between two parallel concentric disks . . 160 Figure C.1 The actual and extrapolated experimental data of Yang et al. (2004) . . . . . . . . . . . . . . . . . . . . . . . . . . 164 Figure D.1 Variation of Re + Dh with dimensionless radial coordinate, r + . . . . . 167 . . 138 xxi NOMENCLATURE English Symbols a modified friction coefficient exponential constant for dimensionless obstacle height, E A cross-sectional area, ?m 2 A s wall surface area, ?m 2 b modified friction coefficient exponential constant for the Reynolds number, Re Dh c modified friction coefficient exponential constant for dimensionless pitch length, D C p constant pressure specific heat, J/(kg K) d pitch length of the channel obstacles, ?m ds change in entropy, J/K D dimensionless pitch length of the channel obstacles, d/H D h hydraulic diameter, ?m e obstacle height, ?m E dimensionless obstacle height, e/H f Darcy friction factor F modified friction coefficient h half of the channel height, or half of the spacing between the disks, H/2, ?m xxii h t total enthalpy, J/kg H the channel height or spacing between the disks, ?m K modified friction coefficient multiplication constant Kn Knudsen number, defined as M/(Re Dh ) 1/2 L max maximum length a gas will travel before reaching Mach 1, ?m M Mach number, defined as V/(?RT) 1/2 N total number of points for integration P static pressure, kPa abs. P t total pressure, kPa abs. Q volumetric flowrate, SCCM, ACCM r radial coordinate, ?m r i radial coordinate of the inlet station r o radial coordinate of the outlet station r + dimensionless radial location, r/r i R gas constant, J/(kg K) Re Dh Reynolds number, defined as ?VD h /? T temperature, K T t stagnation temperature, K v applied voltage, volts V velocity or average velocity, m/s V r velocity component in the radial direction, m/s w width into the paper, ?m x channel distance coordinate, ?m xxiii x + dimensionless channel coordinate, defined as x/D h x lim end point of integration z vertical coordinate, ?m Greek Symbols ? ratio of the specific heats, defined as C p /C v ? deflection height of the boss plate above the seat rings, ?m ? viscosity, (kg m)/s ? density, kg/m 3 ? f wall shear stress, N/m 2 ? step size shrinkage percentage Subscripts i denotes values at the inlet o denotes values at the outlet Superscripts + nondimensionalized by the initial value, with the exception of x + 1 CHAPTER 1 INTRODUCTION 1.1 Background The technological applications of microfluidic devices seem to be limitless. Currently, a great number of applicable concepts are being developed in relation to cooling of microelectronics, drug delivery systems, poisonous chemical sensors, lab-on- chip systems, DNA analysis, and micropropulsion systems (Ouellette, 2003). Micropropulsion systems are of great interest to the aerospace field through the concept of spacecraft miniaturization. Currently, the launch cost range from $10,000 to $100,000 per kilogram of the spacecraft (Janson et al., 1998). A reduction in size (and thus mass) to the micron level would greatly reduce the costs of manufacturing and launching of future space missions. These cost reductions could also allow a greater number of launches for a given budget. Furthermore, the size reduction often assists in attaining a necessary power reduction as well. Let us focus on the proposed microsatellites. Generally, ?standard? satellite size is down to a mass of 1000 kg, ?small? size is considered down to 100 kg, a ?microspacecraft? is down to 10 kg, and a ?nanospacecraft? is 1 kg or less (Janson et al., 1998). A diagram displaying a range of satellite sizes is shown in Figure 1.1. The implementation of microsatellites technology could lead to better performance over current satellite technology. In many proposed applications of microspacecrafts, it is 2 suggested to use dense constellations of theses devices to perform the function of older larger satellites. This could provide many advantages over existing larger satellites such as higher data transmission resolution from a distributed antenna array, and a decreased chance of functional failure (Mueller et al., 2000). That is the loss of a single microsatellite may slightly degrade constellation performance but would not cripple operation. A concept drawing of microsatellite is shown in Figure 1.2, and a photo of a microsatellite designed by the SpaceQuest Ltd. Company is shown in Figure 1.3. Miniaturization of the microspacecraft naturally necessitates marked size reduction of all the various components. Among these, the fluid handling and control issues must be addressed within strict design constraints. Successful design, optimization and fabrication of these microfluidic devices depends on the ability to understand and control the fluid flow through the particular device. Among the typical microfluidic devices, micropumps, microvalves, and microthrusters have been studied extensively. In this thesis, results of modeling of a piezoelectrically-actuated high-pressure gas microvalve is presented. It is important to note that inert helium was used as the working fluid for the microvalve operation experiments because of its low molecular weight, that makes it a good choice for leak testing. Actual micropropulsion systems will handle propellants. 1.2 Microvalves for Micropropulsion System Applications The current study focuses on modeling of gaseous compressible flow through a piezoelectrically-actuated JPL (Jet Propulsion Laboratory) microvalve component of envisioned micropropulsion systems. The size restriction of the system requires that the 3 propulsion tank be highly-pressured in order to provide adequate thrust capability and propellant storage. These two coupled design aspects are expressed in a term called an impulse bit, which is thrust integrated over thruster on-time (Newtons-Seconds) (Mueller et al., 2000). Once in orbit, the propulsion system would be used for operations such as altitude control, turning antenna for data transmission, and maneuvering for optical observations (Yang et al., 2004). These operations would not be frequent. Therefore, the high-pressure micropropulsion system would be normally idle. Thus, the microvalve component would have to be normally-closed, and leak-tight under high pressure. In addition, it must also have a fast enough actuation in order to maneuver the spacecraft accurately. 1.3 Objectives The main objective of this thesis is to determine the variations of the flow properties from an inlet to an outlet port for gaseous compressible flow through the JPL microvalve during static operation. This is achieved through the implementation of two 1-D models, namely: 1. Channel Flow Model, 2. Radial Flow Model. A 4 th order Runge-Kutta algorithm was utilized to integrate a system of coupled nonlinear ordinary differential equations for both models. The channel flow model was developed for the purpose of benchmarking the numerical code. The differential equations for the flow properties are similar to those of 4 compressible flow in a constant-area insulated duct with friction. Analytical relations for the flow properties were derived as the benchmarking tools. The radial flow model involves a more accurate approximation of the microvalve geometry. This model accounts for compressibility flow of a perfect gas between two insulated, parallel disks flowing radially toward an outlet hole at the center of the bottom disk. Additionally, a 2-D Stokes flow model was derived for incompressible, axisymmetric, radial flow between two concentric parallel disks. The results of this model were used to verify the flow property variations from the radial flow model. 1.4 Outline of the Thesis A literature review of microvalves is given in Chapter 2. In Chapter 3, the JPL microvalve is discussed in detail, and a modeling strategy is developed. Chapter 4 focuses on the mathematical formulation and derivation of the differential equations for the flow properties for both models. In Chapter 5, a presentation and discussion of the results are given. Lastly, conclusions and recommendations for future work are discussed in Chapter 6. Figure 1.1: Diagram displaying a range of satellite sizes (source: http://www.sstl.co.uk/) 5 Figure 1.2: Concept drawing of a microsatellite (Janson et al., 1998) 6 Figure 1.3: A microsatellite developed by the SpaceQuest Ltd. Company in a low Earth orbit (source: http://www.spacequest.com/) 7 8 CHAPTER 2 LITERATURE REVIEW OF MICROVALVES Microfluidic devices and systems have a number of advantages over traditional technology associated with their size, including lower power consumption, lower production costs, and high reproducibility in manufacturing with batched fabrication. As with any fluid system, fluid control will almost always involve valves. As a vital component to microfluidic systems, an abundant amount of research has gone into developing microvalves. 2.1 Microvalve Categorization Microvalves can be categorized in a number of different ways, such as the working fluid involved (although most are designed for gaseous flow), application, normally-closed or normally-open, and actuation. This review will focus on microvalve designs in terms of the type of actuation utilized. Specifically, the most popularly employed forms of actuation are as follows: 1. Thermopneumatic, 2. Shape Memory Alloy (SMA), 3. Bi-morph, 4. Electromagnetic, 5. Electrostatic, 6. Piezoelectric. 9 2.1.1 Thermopneumatic Actuation The operating principle behind thermopneumatic actuation is the heating of a fluid in a sealed cavity to its boiling point, which causes the bulging of a silicon membrane covering the cavity to close a flow port. Although these microvalves can be designed to be either normally-closed or normally-open, they are by majority designed as normally-open microvalves. A normally-open thermopneumatic microvalve for high flowrates was developed by Rich and Wise (2003). A schematic diagram of the microvalve cross-section is illustrated in Figure 2.1. A small quantity of a volatile fluid (for instance pentane or methanol) at saturated liguid-vapor phase is contained within the actuator cavity below silicon membranes attached to the valve plate. Upon heating the fluid with resistance heaters on the bottom of the cavity, the pressure rises, deflects the valve plate and closes the flow port. This particular valve works well for an application that requires a high flowrate, but does not require a fast actuation (< 1sec.) or cycling rate. Since the actuation is executed by pressure forces, it allows greater deflection over other forms of actuation. This allows for higher flowrates through the valve. This valve was reported to close with 350 mW at a 1000 torr inlet pressure and maintain closure with a 30 mW input. Under a differential pressure of 1500 torr, the valve had a flowrate of 400 SCCM (Standard Cubic Cm per Minute) with a leak rate of 10 -3 SCCM. Also, the response time of the valve was 1 second, which is relatively high when compared to the response time of other forms of actuation. A normally-open thermopneumatic microvalve with a silicone rubber membrane was fabricated by Yang et al. (1997). The valve has three component parts: the heater 10 wafer, the membrane/cavity wafer, and the valve seat. These components are illustrated in Figure 2.2. This particular microvalve is designed to allow for high flowrates through the use of a rubber membrane as opposed to a silicone membrane. The silicone rubber used is MRTV1 manufactured by American Safety Technologies. The rubber has a very low modulus (~ 1 MPa), high elongation, and provides good sealing on rough surfaces. This rubber is well-suited for microfabrication techniques, since it is resistant to certain chemicals used in the process such as: hydrofluoric acid, positive photoresist developer, and alcohol. The working fluid used (Fluoriner) and the actuation fluid (isopropanol) were non-corrosive to the rubber membrane. The maximum deflection of the membrane was 70 ?m. Under an inlet pressure of 20 psig, the valve was shown to shut of a flow off 1340 ccm/min. with a power input of 280 mW. Thermopneumatic microvalves are useful for high flowrates under specific application guidelines. In general, these guideline limitations include slow actuation, and cycling time (i.e. time to cool down). Additionally, in order to keep power consumption low, volatile actuation fluids are used that require only a slight increase in temperature to create a sufficient increase in vapor pressure for actuation. This requires that the microvalve by impervious to environmental temperature change, since it could cause unwanted actuation or malfunction. 2.1.2 Shape Memory Alloy (SMA) Actuation The operating principle behind SMA actuation is the use of SMA (most common is Ti-Ni) for control of a deflection membrane that either opens or closes a flow port. Shape memory alloys are deformed by heating the material to its transition temperature 11 (Mueller, 2000). When the SMA membrane reaches the transition temperature, it deforms to its parent state. This parent state is established during the fabrication using a high-temperature anneal (Mueller, 2000). These valves can be designed for either normally-open or normally-closed operational purposes. However, due to the lack of controllability of the SMA membrane deflection they are limited to use as an on-off valve. A normally closed microvalve using a SMA actuator was designed by Geon et al. (2000). A schematic cross-sectional view of the microvalve operation is illustrated in Figure 2.3. The microvalve consists of three stacked layers: a flat silicon spring (spring constant is 330 N/m), a patterned TiNi SMA actuator, and an orifice die. The silicon spring applies an initial closing force against the orifice, thus making the valve-normally closed. The opening time was 50 msec and the closing time was 18 msec. The maximum flowrate achieved was 0.17 lpm with an inlet air pressure of 34.5 kPa and a leakage of 0.005 lpm. A normally open microvalve using a SMA actuator was designed by Skrobanek et al. (1997). In contrast to a normally-closed valve, instead of overcoming the seating pressure of a spring the SMA actuator acts against the pressure of the fluid. Schematic cross-sections of the microvalve are given in Figures 2.4 and 2.5. For improvement of the microvalve operation the SMA membrane was stress optimized. Skrobanek et al. (1997) describes stress optimization as designing the SMA membrane so that ?homogeneous spatial stress profiles? are created for a given load pattern. For an operational differential pressure limit of 1200 Pa, air flow of 1600 SCCM and a work output of 35 pNm were observed. Response times for closing the valve varies between 12 0.5 and 1.2 seconds. The cooling time (i.e. the opening time) varies between 1 and 2 seconds. Generally, shape memory alloy microvalves have developed for inlet pressures of 100 psi to 400 psi with corresponding maximum flowrates of 6000 SCCM. The fastest corresponding response times are about 1 ms to open and 20 ms to close. However, cycle times are longer, because of the necessary cooling time to reverse the opening or closing operation. Power requirements range from 0.3 to 2 watts. In addition, leak rates as low as 0.01 SCCM have been measured (Mueller, 2000). 2.1.3 Bi-Morph Actuation Bi-morph (or Bimetal) actuation works on the principle of having a membrane composed of two materials with different coefficients of thermal expansion bonded together. Typically, a thin metal (normally Ni of Al) membrane bonded on top of a thin silicon membrane to form a single membrane is utilized. An electrical resistance heater is embedded or diffused into the silicon membrane. Upon heating, the membrane deflects to open the valve due to the metal?s higher coefficient of thermal expansion to that of silicon (Mueller, 2000). By varying the electrical power input of the heater (the temperature of the bimetal membrane), the deflection can be controlled. A normally-open, bimetallically actuated microvalve was developed by Jerman et al. (1990) and tested in the interest of investigating its accuracy of flow control. A cross- sectional view of the microvalve is shown in Figure 2.6. The valve consists of a aluminum-silicon membrane actuator with diffused resistors attached to a central boss plate that is bonded to the etched valve body. For an inlet pressure of 30 psi, the 13 microvalve provided accurate pressure regulation over a flow range from less than 1 cc/min to around 35 cc/min. A normally-closed, bimetallically actuated, 3-way microvalve was developed by Messner et al. (1998). This valve was designed to control differential pressures of 1000 kPa. A cross-sectional view of the microvalve is illustrated in Figure 2.7. Additionally, an isolated, cross-sectional view of the bimetallic actuator is shown in Figure 2.8. The valve has 3 flow ports and two actuation states. In the off-state, fluid port 1 is sealed by the pre-stressed actuator membrane, and fluid port 2 is connected to fluid port 3. In the on-state, fluid port 1 is opened and fluid port 3 is sealed. Therefore, fluid port 1 is connected to fluid port 2. For an inlet pressure of 600 kPa and a temperature range of 0?C to 50?C, flowrates up to 800 ml/min were achieved. The power consumption was about 1 W. Bimetallically-actuated microvalves are typically easier to fabricate as opposed to other actuation forms, because of its simple design and lack of a need for a complicated external actuator. Generally, this form of actuation is capable of producing large force generation and large displacements. However, as with the thermopneumatically-actuated and SMA-actuated microvalves, the relaxation time has a retardation effect on response time. Also, if the silicon substrate is not isolated from the bimetallic membrane, a heat will dissipate into the substrate. This can cause a large increase in power consumption (Tomonari et al., 2003). 14 2.1.4 Electromagnetic Actuation Electromagnetic actuation utilizes electromagnetic forces produced by coil turns to induce a deflection of a magnetic membrane that seals the flow port. Due to the complications of micromachining a sufficient number of coil turns, most electromagnetic valves use external macro-sized coils or permanent magnets (Mueller, 2000). Ideally, the entire microvalve would be of micromachined components. An electromagnetic microvalve was fabricated using a combination of miniature and micromachined components by Shinozawa et al. (1997). This particular microvalve is a good example of the difficulties and inefficiencies of a ?hybrid? electromagnetic microvalve. A schematic of the microvalve is illustrated in Figure 2.9. A small permanent magnet is attached to the valve plug, and is actuated by an electromagnetic force from an external micromachined coil. The vertical displacement of the valve plug opens and closes the valve. The microvalve is about 5x5x5 mm in size and has a 0.5 mm maximum displacement of the valve plug. The smallest controllable amount of water was 0.7 ?l/min and the maximum controllable amount was 900 ?l/min. The authors make note of the difficulties associated with combing components of different scales. These problems include clogging of the flow passage and insufficient actuation force for a millimeter-sized region. A fully-micromachined, normally-open electromagnetically actuated microvalve was designed and fabricated by Bintoro and Hesketh (2005). This is of particular interest, not only because is it fully-micromachined, but also because it was fabricated from using a single silicon wafer. There is no wafer bonding to create layers, which amounts to a lower fabrication cost. A schematic of the microvalve is shown in Figure 15 2.10. In Figure 2.11, there is a schematic displaying the operation of the microvalve. The membrane is separated from the microcoil by a 12 micron gap. To close the microvalve, an electric current is applied to the microcoil to produce an electromagnetic force that displaces the soft magnetic dome over the outlet port. The working fluid used was a 50% diluted methanol with water. This microvalve was able to close for a maximum inlet pressure of 26.5 kPa. The leak rate was in the range of 0.024-0.033 ?l/min over a differential pressure range of 1.5-17 kPa. Also, the range of operational power of the microvalve was 0.3-1.2 W. In general, electromagnetic microvalves have been reported to have response times less than 1 ms, with valve strokes in the range of 10-15 microns. Additionally, electromagnetic valves have low power consumption. However, the electrostatic forces do not fair well against high pressures. Also, for a normally-open microvalve, a loss of power would cause the valve to fail open (Mueller, 2000). 2.1.5 Electrostatic Actuation The operating principle of electrostatic actuation utilizes the electrostatic force between two electrode plates when a voltage is applied between them. The actuator?s closure membrane is a cantilever structure over the inlet port. Typically, these valves are designed for normally-open operation. A normally-open electrostatically actuated silicon microvalve for gaseous flows was developed by Ohnstein et al. (1990). A schematic of the microvalve is shown in Figure 2.12. The microvalve is approximately 600 x 600 microns, with the closure plate being about 350 x 390 microns. The microvalve operational maximum inlet pressure is 16 114 mmHg, with flows of up to 150 SCCM. The leak rate of the closed valve was approximately 0.1 SCCM. A normally-closed, electrostatically-actuated microvalve was developed by Huff et al. (1990). This is one of the few examples of a normally-closed, electrostatic valve. A schematic of the valve is shown in Figure 2.13. This paper presents a unique concept in using the pressure from the incoming fluid to assist the electrostatic actuation used to open the valve. This paper discussed the design and fabrication of the microvalve, but not flow testing. One foreseeable problem with this design is that it limits the maximum allowable inlet pressure to retain closure. Generally, electrostatic valves provided the benefits of both low power consumption and fast response times. However, their greatest drawback is the electrostatic forces are relatively weak. Therefore, these microvalves are not useful under high pressures. Additionally, the electrostatic forces are limited by the deflection, because the force generated between two electrode plates is inversely related to the spacing between them (Roberts et al., 2003). 2.1.6 Piezoelectric Actuation Piezoelectric actuation utilizes piezoelectric materials that deform slightly with an applied voltage. Since a single piezoelectric element provides only a small deflection, the elements are stacked on top of each other, because with this arrangement the deflections are additive. Additionally, this configuration has the advantages of larger deflections for smaller applied voltages, and higher actuation forces at smaller deflections (Mueller, 2000). 17 A normally-open, piezoelectrically-actuated microvalve for micropumping systems was developed by Roberts et al. (2003). A schematic of the microvalve is shown in Figure 2.14. The microvalve has three major components: a piezoelectric actuator element, an enclosed ?hydraulic amplification chamber? (HAC), and a membrane with an attached valve sealing cap. With an applied voltage to the piezoelectric element, the element strains and deflects the piston, thus generating pressure in the ?HAC? which pushes the valve sealing cap against the inlet port. The working fluid was a degassed silicon oil. The microvalve can operate for pressures greater than 300 kPa, and produces strokes of 20-30 microns. The maximum average flowrate through the microvalve was 0.21 ml/s for a valve opening of 17 microns and a differential pressure of 260 kPa. In general, piezoelectric microvalves have an operational applied voltage range of 25-100 volts. Additionally, they have demonstrated leak rates of 0.1 SCCM or less. The advantages of piezoelectric microvalves include fast response times, and substantial actuation forces. The major drawbacks of these microvalves are that they require larger operating voltages than other actuation forms and they require a stacked piezoelectric element arrangement to achieve substantial actuation deflections. The stacked arrangement also makes the fabrication more complicated (Mueller, 2000). 2.2 Closure Six basic actuation methods including examples have been described in the previous sections. The main criteria one must consider when choosing a microvalve actuation method for an application are actuation cycle time (response time), maximum allowable inlet pressure, actuation force, displacement, power consumption, necessary applied voltage for actuation, leakage, and reliability. A general summary of some of these criteria for the different actuation methods is given in Table 2.1 (Ayhan, 2000). Additionally, an evaluation of the microvalve actuation methods as they apply to the application of micropropulsion is given in Table 2.2 (Mueller, 2000). Table 2.1 Characteristics of Microvalve Actuation Methods (Ayhan, 2002) Actuators Force Displacement Response Time Reliability Solenoid Plunger Small Large Medium Good Piezoelectric Very large Medium Fast Good Pneumatic Large Very Small Slow Good SMA Large Large Slow Poor Electrostatic Small Very small Very Fast Very Good Thermopneumatic Large Medium Medium Good Electromagnetic Small Large Fast Good Bimetallic Large Small Medium Poor Table 2.2 Evaluation of Microvalve Actuation for Microspacecraft Applications (Mueller, 2000) Thermopneumatic Bi-Morph Shape Memory-Alloy Electrostatic Piezoelectric Electromagnetic Size & Weight Excellent Excellent Excellent Excellent Excellent Excellent Power Good Good Good Excellent Excellent Excellent Voltage Acceptable Unknown Unknown Poor Poor Acceptable Cycle Time Poor Poor Poor Excellent Excellent Excellent Pressure Marginal Marginal Marginal Poor Unknown Unknown Leakage Poor Poor Poor Poor Unknown Unknown Seating Pressures Acceptable Acceptable Acceptable Poor Good Good 18 Figure 2.1: Cross-section of a thermopneumatic microvalve (Rich and Wise, 2003) 19 Figure 2.2: Cross-section of a thermopneumatic microvalve (Yang et al., 1997) 20 Figure 2.3: Cross-section of a SMA microvalve , a) closed position, b) open position (Geon et al., 2000) 21 Figure 2.4: Schematic cross section of the SMA microvalve (A, B, and C denote the valve ports) (Skrobanek et al., 1997) 22 Figure 2.5: Schematic the SMA microvalve cross-section during open operation (Skrobanek et al., 1997) 23 Figure 2.6: Cross-section of a bimetallically actuated microvalve (Jerman, 1990) 24 Figure 2.7: Cross-section of a 3-way bimetallically actuated microvalve (Messner et al., 1998) 25 Figure 2.8: Cross-section of a bimetallic actuator of a 3-way microvalve (Messner et al., 1998) 26 Figure 2.9: Cross-section of an electromagnetically-actuated microvalve (Shinozawa et al., 1997) 27 Figure 2.10: Schematic of an electromagnetic microvalve (all dimensions are in ?m). 1 = inlet orifice; 2 = microvalve?s base; 3 = Au microcoil; 4 = outlet orifice; 4 ? = gasket (defined by the center of Au microcoil); 5 = circular support (NiFe); 6 = center soft magnetic dome (NiFe); 7 = membrane?s supported legs (NiFe). The component no.?s 6 and 7 form the microvalve?s membrane. (Bintoro and Hesketh, 2005) 28 Figure 2.11: Schematic of an electromagnetic microvalve?s operation. (a) The fluid flows from the inlet orifice to the outlet orifice and beneath the membrane. (b) Current (Icoil) is drawn to the microcoil and it produces an electromagnetic force (FEM) that displaces the membrane downward. When the bottom of the membrane touches the gasket, the fluid flow is choked, i.e. the microvalve is closed. (Bintoro and Hesketh, 2005) 29 Figure 2.12: Schematic of an electrostatically-actuated microvalve (Ohnstein et al., 1990) 30 Figure 2.13: Schematic of an electrostatically-actuated microvalve (Huff et al., 1990) 31 Figure 2.14: Schematic of a piezoelectrically-actuated microvalve (Roberts et al., 2003) 32 33 CHAPTER 3 MICROVALVE MODELING APPROACH 3.1 Background and JPL Microvalve Design Requirements In this Chapter, details of a JPL high-pressure piezoelectrically-actuated microvalve will be given. This will be followed by the steps leading to simplification of the complex microvalve geometry to a 1-D model. The microvalve modeled in this study was designed by Drs. E-H Yang, C. Lee, J. Mueller, and T. George of the MEMS Technology Group at NASA?s Jet Propulsion Laboratory. The following description of the purpose, design, and operating conditions of the microvalve was taken from Yang et al. (2004): Constellations of microspacecrafts (10 kg total mass) are being envisioned to study magnetic fields or radiation belts around the Earth. By using large constellations, tensor mapping of fields and particles may be conducted. The Magnetic Constellation mission, which has recently been approved by NASA, seeks to map Earth?s magnetic field with 50 ? 100 spacecrafts, each equipped with its own magnetometer. Such large constellations of spacecraft are only feasible if very small spacecraft are used in order to keep the total launch mass reasonable and the launch cost affordable. Constellation spacecraft may have propulsive requirements, either to maintain a formation, or to turn (slew) the spacecraft to point an antenna to Earth for data transmission, or to aim a camera for observation. In case of such small spacecraft, significant propulsion system size and mass reduction over current state-of-the-art is required for these subsystems to fit within the greatly reduced mass and size envelope. Thrust and impulse bit capabilities may also be required to be very small depending on spacecraft mass and required pointing accuracy. Required impulse bits may range from the mNs-range for larger craft having relatively coarse attitude requirements, into the ?Ns-range and possible even nNs-range for very tight pointing requirements and very small spacecraft. Thus, there exists a need for very low impulse bit, micro-Newton thrust level propulsion systems in order to provide the required pointing accuracies for the attitude control of the microspacecraft. Such micropropulsion systems require precisely controlled, extremely small propellant flow from a pressurized propellant tank. A fast actuation, leak-tight microvalve at high propellant pressures is required for micropropulsion applications as described in Table 3.1. Table 3.1 Microvalve Requirements for NASA?s Miniature Spacecraft Propulsion Needs, Compared with Reported Performance to Date (Yang et al., 2004) Requirements Target Demonstrated Leak Rate < 0.005 sccm/He 10 -4 sccm/He at 800 psi Response time < 10 ms < 10 ms Inlet Pressure 300 ~ 3000 psi 0 ~ 1000 psi Power < 1 W ~ 4 mW (static) Temperature 0 ~ 75 ?C Not tested yet. Ideally, this valve will be tightly integrated with thruster components to allow for a compact and lightweight overall thruster module design. Solenoid-based miniaturized valves have been developed. Some of these valves meet most of micropropulsion requirements. However, they still consume several Watts 34 35 to operate the valves. Microspacecraft systems are anticipated to have severely limited power budgets. It is therefore desirable to incorporate ?low power? valves meeting all the requirements for micropropulsion. Other previously reported microvalves do not meet the requirements for pressure range and/or leak rate needed for micropropulsion. Recently, leak-tight microvalves operating at 10 atm have been reported, which still fall short of pressure range. Thermally actuated microvalves usually have slow response time (~ 100 ms to complete a cycle), which is unacceptable for micropropulsion applications. This is because the slow valve actuation time causes long thruster on-times and wide impulse bits. Thermally actuated valves also suffer from the risk of random valve opening if ambient heating or cooling occurs, resulting in uncontrolled initiation of the actuation mechanism. Most microvalves reported previously have shown marginal valve seating at high pressures. Microvalves without adequate seating are exposed to severe problems in leakage and pressure handling capability. Therefore, significant efforts are required for the development of high-pressure microvalves to meet the micropropulsion requirements. In response to such needs, a fully characterized leak- tight piezoelectric microvalve, operating under extremely high inlet pressures for micropropulsion applications has been developed. 3.1.1 Design of the Microvalve The piezoelectric microvalve described in this paper consists of a seat plate, a boss plate and an actuator as shown in Figure 3.1. The microvalve components do not contain fragile membranes in order to allow high-pressure operation. Major elements of the microvalve design include its seating configuration with narrow seat rings. The 36 seating configuration is provided by an initial opening pressure attributable to the tensile stress in the silicon tether extended by the valve seat as shown in Figure 3.2. A series of narrow rings on the seat plate is designed to reduce potential leakage due to scratches over a seat ring. The narrow rings reduce contact area, increasing the seating pressure and consequently reducing internal leaks. An additional advantage of the narrow/hard- seat design is that contact pressures may be high enough to crush contaminant particles, thereby also reducing the leakage attributable to contaminants in the flow. The boss plate has a 2 ?m thick silicon dioxide layer as a hard seat material in the boss-center plate. The outer part of the boss plate is a metal-to-metal compression bonded to the seat plate. The boss-center plate covered by the silicon dioxide is slightly thicker than the outer part. This causes the boss-center plate to be pressed toward the seat plate by the stretched tether, enhancing a leak-tight valve operation. The piezoelectric stack actuator exhibits a very high block-pressure (50 MPa in this case), providing enough pressure to overcome the high differential pressure in addition to the downward bending stress from the tethers. The custom designed stack of piezoelectric actuators consists of active zones and an inactive central zone. The piezoelectric stack with mechanically separated (by deep U-grooves) active zones is bonded to the boss plate within a rigid housing. Application of a potential (~40 V) to the stack makes the active zones vertically expand by 5 ?m, lifting the boss center plate, which is bonded to the inactive zone of the stack, away from the seat plate. This actuation creates a channel between the two openings, allowing for the passage of fluids as shown in Figure 3.3. Since the piezoelectric element is essentially a stacked capacitor, the actuator consumes an extremely low power 37 when it is not moving, thus making it possible to achieve a nearly zero-power, normally- closed valve system (Yang et al., 2004). Several figures have been included to assist in visualization and understanding of the microvalve geometry. Scanning electron micrographs (SEM) of the seat plate and the boss plate are shown in Figure 3.4. These SEMs clearly show the geometry and scale of these components. The 1.5 ?m by 10 ?m cross-section of a seat ring is depicted well. In addition, a detailed perspective of the boss plate is shown. The flexible tethers and the center plate that deflects with actuation are easily observed. A photograph depicting the size of the packaged microvalve is shown in Figure 3.5. Additionally, isometric drawings in Figures 3.6 and 3.7 illustrate various components and their configuration within the housing of the packaged microvalve. Lastly, a cross-sectional schematic diagram of the dimensions within the housing is shown in Figure 3.8. 3.2 2-Dimensional Axisymmetic Geometry Model There are several difficulties in computational modeling the 3-dimensional flow field within the complex geometry of the microvalve with all of its details. For instance, the overall complexity and corresponding problem of grid generation with such drastic differences between the millimeter-scale housing of the cavity and the micron-scale gap poses a great challenge. Furthermore, interaction of the fluid with a moving boundary during dynamic operation is another complication that might require dynamic re- meshing. Therefore, instead of solving the 3-D compressible flow equations for the gas flow in the actual geometry, simplifications of geometry are considered for the static operation. Observing Figures 3.3b and 3.8, one notices the relatively large cavity space as compared to the 11 to 15 ?m gap created during open static operation. Therefore, it can safely be assumed that within most of this space the velocity of the gas is negligible ( ) and the gas is stagnant. Thus, the cavity represents the stagnation chamber and the dynamic component of the total (stagnation) pressure within the cavity is assumed to be negligibly small. This leads to the assumption that the static pressure within the cavity is equal to the total pressure (Figure 3.8). Furthermore, it can be hypothesized that a major part of the total pressure drop through the valve occurs within the micron-size gap between the boss plate and seat plate (11 to 15 ?m). This hypothesis is reasonable due to the small spacing between the boss and seat plates, in addition to the excessive friction imposed on the fluid by the rings. In support of this hypothesis, analysis was done to determine the total pressure drop through the inlet piping to the housing cavity (Figures 3.3b, 3.5-3.8) of the microvalve. First, equivalent lengths (L 0?V e /D) and loss coefficients (K) were determined for the piping network shown in Figure 3.9 (White, 1999). Then, the total pressure drop for internal incompressible viscous flow (White, 1999) was found. The results for the incompressible analysis are given in Table 3.2, and a sample calculation for the total pressure drop is given in Sections 3.2.1. 38 Table 3.2 Total Pressure Drop (?P t ) for the Inlet Piping Network (the flowrates were extrapolated from measurements of Yang et al. (2004); details are in Appendix C) 11 1.844E-06 3.522E-08 3.488E-01 0.0441 12 7.350E-06 1.403E-07 1.392E+00 0.1761 13 1.636E-05 3.124E-07 3.107E+00 0.3929 14 2.891E-05 5.521E-07 5.510E+00 0.6967 15 4.499E-05 8.591E-07 8.612E+00 1.0890 11 1.583E-06 5.937E-08 2.996E-01 0.0202 12 8.390E-06 3.146E-07 1.593E+00 0.1076 13 2.047E-05 7.677E-07 3.914E+00 0.2644 14 3.783E-05 1.419E-06 7.301E+00 0.4932 15 6.047E-05 2.268E-06 1.181E+01 0.7979 11 1.225E-06 6.642E-08 2.318E-01 0.0107 12 6.797E-06 3.685E-07 1.292E+00 0.0595 13 1.677E-05 9.094E-07 3.213E+00 0.1481 14 3.115E-05 1.689E-06 6.035E+00 0.2781 15 4.994E-05 2.707E-06 9.815E+00 0.4523 300 100 200 Percentage of total ?P t Inlet Incomp. ?P t [kPa] P ti [psig] H [?m] Q actual [m 3 /min] Mass Flowrate [kg/s] The total pressure drop for the inlet piping network is miniscule compared to the total pressure that is measured upstream of the microvalve system. The percentage of the total pressure drop of the inlet piping network as compared to that of the complete system is less than 1.1%. 3.2.1 Sample Calculation for Incompressible Total Pressure Drop of the Inlet Piping Network The experimental data of Yang et al. (2004) were extrapolated thus leading to 15 flow cases dependent on the initial total pressure (P ti ) and the gap (H) between the boss and seat plates. A detailed description of these 15 cases is given in Section 5.1.1 and the initial conditions are shown in Table 5.1. For this sample calculation the case of P ti = 200 39 psig and a gap of H = 15 ?m is considered. First, a dimensionless equivalent length (L e /D) for the piping network is needed. For the three rounded 90? bends, a dimensionless equivalent length of 12 diameters was determined from Figure 8.16a of Fox et al. (2004). A dimensionless equivalent length of 1000 diameters was used for the needle valve, and 8 diameters for the toggle valve (http://www.carf-engineering.com, Warring, 1982) . The dimensionless equivalent lengths for the different diameter sections of the piping network are then calculated as follows: () ,116481000)12(3 7625.4 1272.203)2.76(29.88 1 =+++ +++ = ? ? ? ? ? ? ? ? mm mm D L e (3.1a) . ,7 1 7 2 == ? ? ? ? ? ? ? ? mm mm D L e (3.1b) ,18 33.0 6 3 == ? ? ? ? ? ? ? ? mm mm D L e (3.1c) For the sharp 90? bend in the last section of the inlet tube, a dimensionless equivalent length of 60 diameters was determined from Figure 8.16b of Fox et al. (2004): () .11060 2.0 5.65.3 4 =+ + = ? ? ? ? ? ? ? ? mm mm D L e (3.1d) The sudden contraction loss coefficients for the three abrupt changes in diameter within the piping network were determined from Figure 6.22 of White (1999). They are listed as follows: From D = 3/16 in. to 1 mm: K 1 = 0.39, (3.2a) 40 41 From D = 1 mm to 1/3 mm: K 2 = 0.37, (3.2b) From D = 1/3 mm to 0.2 mm: K 3 = 0.28. (3.2c) The volumetric flowrate extrapolated from the experimental results of Yang et al. (2004) had to be converted from standard form to actual form (see Appendix C). The mass flowrate is calculated by multiplying the initial actual volumetric flowrate (Q actual ) and the initial density (? i ). The average velocities (V avg ) for the four different diameter sections were found by dividing volumetric flowrate by the respective cross-sectional areas. The Reynolds numbers (Re D ) for the respective piping sections were then calculated. The average velocities and corresponding Re D are listed as follows: For D = 3/16 inch tube: V avg1 = 0.057 m/s, Re D1 = 31, (3.3a) For D = 1 mm tube: V avg2 = 1.28 m/s, Re D2 = 148, (3.3b) For D = 1/3 mm tube: V avg3 = 11.55 m/s, Re D3 = 444, (3.3c) For D = 0.2 mm tube: V avg4 = 32.08 m/s, Re D4 = 739. (3.3d) Finally, the total pressure drop (?P t ) is given by: () .8.11 2 08.32 )110( 2 08.32 28.0 2 55.11 )18( 2 55.11 37.0 2 28.1 )7( 2 28.1 39.0 2 057.0 )1164( 25.2 22Re 64 2 739 64 22 444 64 22 148 64 22 31 64 4 1 22 3 1 kPa V K V D L P m kg avg j avg e D it jj j = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ++ ? ? ? ? ? ? ? ? ++ ? ? ? ? ? ? ? ? + = + ? ? ? ? ? ? ? ? =? ? + ? (3.4) 3.2.2 Outline of the 2-Dimensional Axisymmetric Geometry Model Based on the aforementioned hypothesis and supported by the analysis given in the previous section, a simplified 2-D axisymmetric geometry of the gap region is considered. The flow problem between the boss and seat plates is then considered to be steady, axisymmetric and compressible flow between two parallel disks flowing radially toward an outlet hole at the center of the bottom disk (Figures 3.10 & 3.11). The proposed 2-D flow model with compressibility effects is still a complicated problem to solve. The 2-dimensionality of the problem due to the presence of the rings necessitates utilization of CFD tools. However, the geometry can be simplified further to a 1-D axisymmetric model. 42 43 3.3 1-Dimensional Axisymmetric Geometry Model Employing the Fanno frictional analysis, the 2-D axisymmetric geometry is simplified to a 1-D axisymmetric geometry through the representation of the seat rings with an apparent increase in flow friction. The effective increase in flow friction due to the presence of the seat rings is determined from the experimental data of Luy et al. (1991). In addition, the disks are assumed to be insulated, and thus the flow can be considered adiabatic. This is reasonable because of the short length of the radial geometry, and the large magnitude of the flow velocity relative to that length, which makes heat exchange negligible. The 1-D axisymmetric geometry model is then posed as steady, axisymmetric, compressible frictional flow of a perfect gas between two insulated, parallel disks flowing radially toward an outlet hole at the center of the bottom disk (Figure 3.12). 3.4 Closure The complex microvalve geometry coupled with complicated dynamic operation has been simplified to a 1-D axisymmetric problem under static operating conditions. The total pressure drop from the pressure sensor to the microvalve cavity has been shown to be negligible. This allows for the assumption that the initial conditions at the edge of the seat rings is at the same total pressure measured by the upstream pressure gauge. The major portion of the total pressure drop is hypothesized to occur as the gas flows between the boss and seat plates and over the seat rings. Substitution of increased flow friction due to the presence of the seat rings was proposed. A correlation for experimental data (Luy et al., 1991) is to be derived for a modification in the friction coefficient of the flow. 44 In Chapter 4, a system of coupled simultaneous nonlinear ordinary differential equations are derived for predicting the flow properties within the framework of the proposed 1-D model. Figure 3.1: Schematic diagram of the leak-tight piezoelectric microvalve (Yang et al., 2004) 45 Figure 3.2: Cross-sectional microvalve configuration, showing the suspension of the boss plate by tensile stressed silicon tethers over the extended valve seat (Yang et al., 2004) 46 Figure 3.3: Operation principle of the microvalve: (a) Microvalve closed (cross- section A-A), and (b) Microvalve open (cross-section B-B) (Yang et al., 2004) 47 Figure 3.4: Scanning electron micrographs (SEM) of the seat plate and the boss plate (Yang et al., 2004) 48 Figure 3.5: Packaged high-pressure piezoelectric microvalve (Yang et al., 2004) 49 Figure 3.6: Isometric view of the microvalve enclosed in the housing 50 Figure 3.7: Exploded isometric view of the microvalve enclosed in the housing 51 Figure 3.8: Schematic diagram of the microvalve housing and geometry showing the main features (not drawn to scale) 52 Figure 3.9: Schematic diagram of the piping network leading to the microvalve cavity (not drawn to scale) 53 Figure 3.10: Isometric view of the radial flow problem with arrows indicating the directions of the incoming and outgoing flow streams 54 Figure 3.11: Top view of the radial flow problem with the boss plate not shown 55 Figure 3.12: Geometry of the 1-D Axisymmetric Model 56 57 CHAPTER 4 MATHEMATICAL FORMULATION The flow in the actual microvalve geometry is compressible, 3-dimensional and can be affected by possible non-continuum effects. In addition, interaction of the fluid with a moving boundary during dynamic operation is another complication. Instead of solving the 3-D compressible flow equations, two simpler models were adopted. For both models, a 1-D Fanno analysis is used to determine the variations of flow properties from the inlet to the outlet for the conditions of static operation of the microvalve. For both models, the effect of the presence of seat rings can be incorporated via a modified correlation for the friction coefficient. This Chapter deals with the development of the two models, namely: 1. Channel Flow Model, 2. Radial Flow Model. The channel model was developed for the sake of the benchmarking of a computer code, whereas the radial model is used to uncover the flow features in the microvalve. 4.1 Channel Flow Model A 1-D channel flow model was employed as a preliminary groundwork and a simplified starting point of the study. The intent was to keep the model similar to 58 compressible flow in a constant-area insulated duct with friction. Therefore, the effects of the area change and heat addition are ignored. In doing so, analytical solutions for this model are tabulated in the Fanno flow tables that are available in compressible flow texts (e.g. John, 1984). The results of numerical integration of the nonlinear differential equations derived in this section will then be compared to the analytic solutions. In relation to the microvalve application, the major drawback of the channel flow model is that it disregards the area change, but this restriction will be removed upon development of the radial disk flow model. The 2-D version of the channel model geometry is idealized as steady, compressible flow of a perfect gas between two infinite insulated parallel plates with seat rings on the lower plate (Figure 4.1). The presence of the seat rings are accounted for by a modified friction coefficient in the 1-D version of the channel model (Figure 4.2). Adiabatic conditions prevail because of insulated plates, short length of the geometry, and the large magnitude of the flow velocity relative to that length, which makes heat exchange negligible. In this situation, the flow properties are only affected by the retarding action of friction on the two walls. Therefore, only the significance of friction is to be considered. 4.1.1 Governing Equations In view of the simple geometry under consideration, the Cartesian coordinate system is an obvious choice (Fig. 4.2). The following assumptions are made: 1. Steady flow, 2. One-dimensional flow (properties only depend on one spatial coordinate, i.e. x-coordinate in Fig. 4.2), 3. Perfect gas with constant specific heats, 4. Adiabatic flow with no external work. The governing equations, i.e., the continuity, momentum, and energy equations are then written in the following forms: Since the cross-sectional area does not change, the continuity equation is: V? = constant, (4.1a) or in difference form: 0=+ V dVd ? ? . (4.1b) The momentum equation is: ()AdVVF S xx vv ?= ?? ? ? , (4.2) and it will be discussed in greater detail below. The energy equation is: =+ 2 2 V h constant , (4.3a) or in difference form: 0=+VdVdh . (4.3b) 59 4.2 Derivation of the Differential Equations for Flow Properties The derivation of the expressions for flow properties is similar to the derivation of the Fanno flow expressions. Consider the differential control volume shown in Figure 4.3 that extends from x to x + dx. The flow properties change (identified by the symbol ?d? for ?differential?) as the flow moves from one end of the control volume to the other end. Listed below are some appropriate definitions for the infinite-parallel-plates geometry under consideration. Perimeter =2(h+w), (4.a) with w being the width into the paper. 60 ) Wall Surface Area, . (4.b) ( dxperimeterA s = The hydraulic diameter is defined as: () () h h wh hw perimeter A D w h 4 2 8 22 24 lim 4 == ? ? ? ? ? ? ? ? + == ?? , (4.c) where A = 2hw is the cross-sectional area of the channel. Assuming fully-developed laminar flow of an incompressible fluid between two infinite parallel plates, a parabolic velocity distribution is established between the two plates. The appropriate relation for the wall shear stress (? f ) is (White, 1999): h V dy du hy f ? ?? 3 == += , (4.5a) with ? being the fluid viscosity. The friction Darcy factor ( ) is then: f h D f Vh V f Re 96 )4( 96 4 2 2 1 = ? ? ? ? ? ? ? ? == ? ? ? ? , (4.5b) with V being the average velocity between the two plates. The difference equations for the definitions of the Mach number and the constitutive relation of a perfect gas are: Mach number for a perfect gas is defined as: RT V M ? = , (4.6a) or in two difference forms: T dT V dV M dM ?= 2 2 2 2 , (4.6b) T dT V dV M dM 2 1 ?= . (4.6c) Constitutive relation for a perfect gas is: 2 2 M V RTP ? ? ? == , (4.7a) or in difference form: T dTd P dP += ? ? . (4.7b) 4.2.1 Derivation of the Differential Equation for the Mach number (M) Accounting for all the forces acting on the 1-D differential control volume (Figure 4.4), the momentum equation in the x- direction is: () ()( ) ( )VAVdVVAVAAdPPPA sf ??? ?+=?+? , (4.8a) which simplifies to: AVdVAAdP sf ?? =?? . (4.8b) Substituting for the wall shear stress (? f ) (Eq. 4.5a) in the momentum equation and dividing both sides by the cross-sectional area (A = 2hw) gives: 61 VdV D dx fVdP h ?? =?? 2 2 1 . (4.9) Dividing the momentum equation by 2 2 M V P ? ? = , one gets: 0 22 2 1 =++ V dV M D dx fM P dP h ?? . (4.10) Combining the perfect gas equation (4.7b), the Mach number equation (4.6c) and recalling the continuity relation (Eq. 4.1b) yields: M dM T dT P dP ?= 2 1 . (4.1) Substituting for P dP from Eq. (4.11) and V dV from the Mach number equation (4.6c) into Eq. 4.10, one gets: 0 2 2 22 2 1 2 1 =+++? T dTM M dM M D dx fM M dM T dT h ? ?? . (4.12) In order to obtain an expression for the Mach number as a function of the axial location in the channel, the temperature term must be eliminated. This is done by substituting an expression for temperature as a function of the Mach number. Assuming a perfect gas with constant specific heats, the energy equation (Eq. 4.3a) takes on a difference form: 2 2 dV dTC P ?= . (4.13a) Dividing by and substituting TC P ()1? = ? ?R C P and RT V M ? 2 2 = , this becomes: 2 2 2 2 1 V dV M T dT ? ?= ? . (4.13b) 62 Substituting for 2 2 V dV from the definition of the Mach number (Eq. 4.6b) gives: () 2 2 1 1 1 M MdM T dT ? + ? ?= ? ? , (4.13c) which is the same as the difference form for temperature given in compressible flow texts (Saad, 1993). Equation (4.13c) can also be written as a differential equation if both sides are divided by dx: () ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ?= 2 2 1 1 1 M MT dx dM dx dT ? ? . (4.13d) Substituting the expression for T dT into the most current form of the momentum relation (Eq. 4.12) and upon simplification, one gets: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? = 22 2 2 1 2 1 1 1 MM M M dM D fdx h ? ? , (4.14) that then leads to the differential equation for the Mach number of the fluid in the channel: )1(2 2 1 1 2 23 M M D f M dx dM h ? ? ? ? ? ? ? ? + = ? ? . (4.15) This constitutes as a nonlinear ordinary differential equation that can be integrated to determine the dependence of M on axial distance x. 63 It should be noted that Eq. (4.14) can be integrated analytically provided that f is constant (this is the case in this model). If the distance that a fluid at Mach number M travels to reach M = 1 is denoted by L max , an expression for h D fL max will be formed that depends on ? and M (Appendix E). On the other hand, an explicit expression for M as a function of x does not exist which necessitates the numerical integration of Eq. (4.15). 4.2.2 Derivation of the Differential Equation for Pressure (P) Substituting for V dV from the Mach number equation (4.6c) into equation (4.10) and utilizing Eq. (4.13c), one gets: () 0 2 1 1 1 2 2 2 122 2 1 = ? + ? ?++ M MdMM M dM M D dx fM P dP h ? ?? ?? . (4.16a) Substituting for h D dx f 2 1 and simplifying this relation, one gets: () ? ? ? ? ? ? ? ? ? ? ? ? ? + ?+ ?= 2 2 2 1 1 11 M M M dM P dP ? ? , (4.16b) that is the same as the difference form for pressure given in compressible flow texts (e.g. Saad, 1993). Equation (4.16b) can be written in its differential form: () ? ? ? ? ? ? ? ? ? ? ? ? ? + ?+ ? ? ? ? ? ? ?= 2 2 2 1 1 11 M M M P dx dM dx dP ? ? . (4.16c) 64 4.2.3 Derivation of the Differential Equation for Velocity (V) Substituting for P dP from the perfect gas equation (4.7b) into equation (4.10), one gets: 0 22 2 1 =+++ V dV M D dx fM T dTd h ?? ? ? . (4.17a) Substituting for ? ?d from the Continuity equation (4.1b), h D dx f 2 1 and T dT , one gets: () 0 2 1 1 1 2 1 1 1 2 22 2 2 2 =+ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? + ? + ? ?? V dV M MM M M dM M M MdM V dV ? ? ? ? ? ? . (4.17b) Simplifying this leads to: ? ? ? ? ? ? ? ? ? ? ? ? ? + = 2 2 1 1 1 M M dM V dV ? , (4.17c) which is the same as the difference form for velocity given in compressible flow texts (Saad, 1993). Equation (4.17c) can be rewritten as: ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? = 2 2 1 1 1 M M V dx dM dx dV ? . (4.17d) 65 4.2.4 Derivation of the Differential Equation for Density (?) From the continuity relation (Eq. 4.1b): V dVd ?= ? ? , (4.18a) or: ? ? ? ? ? ? ? ? ? ? ? ? ?= dx dV Vdx d ?? . (4.18b) In addition, differential equations for stagnation or total pressure, P t , and entropy change, ds can be derived. 4.2.5 Derivation of the Differential Equation for Total Pressure (P t ) The differential expression for P t is derived from the isentropic relation, but first an expression for total temperature must be derived beginning with the energy equation. Define the total enthalpy as: 2 2 V hh t += . (4.19a) Assuming a perfect gas with constant specific heats, )( TTChh tpt ?=? , (4.19b) thus leading to: ? ? ? ? ? ? ? ? +=+= TC V TT C V T pp t 2 1 2 22 . (4.19c) 66 Introducing 1? = ? ?R C p into (4.19c), one gets: ? ? ? ? ? ? ? ? ? += RT V TT t ? ? 2 )1( 1 2 . (4.19d) Substituting the definition of the Mach number from equation (4.6a), one gets: ? ? ? ? ? ? ? ? ? += 2 2 )1( 1 MTT t ? . (4.20) This expression is then substituted into the isentropic relation: )1( ? ? ? ? ? ? ? ? ? = ? ? T T P P tt , (4.21a) )1( 2 2 )1( 1 ? ? ? ? ? ? ? ? ? ? += ? ? ? M P P t , (4.21b) )1( 2 2 )1( 1 ? ? ? ? ? ? ? ? ? ? += ? ? ? MPP t . (4.21c) Letting ? ? ? ? ? ? ? += 2 2 )1( 1 MB ? and differentiating with respect to the Mach number, )( )1( 1 )1( MPB dM dP B dM dM M P dM dP P P dM dP ttt ? ?? ? ?? += ? ? + ? ? = (4.22a) Multiplying by dM, and dividing by P t , one gets: dMM P P BdP P B P dP ttt t )( )1( 1 )1( ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? = ? ? . (4.22b) 67 Substituting , )1( ? = ? ? B P P t and , )1( P P B t = ?? ? one gets: dMM M P dP P dP t t )( 2 1 1 1 2 ? ? ? ? ? ? ? ? ? + += . (4.2c) Substituting , we get: 2 2 dMMdM = ? ? ? ? ? ? ? + ? ? ? ? ? ? += 2 2 2 1 1 2 M dM P dP P dP t t ? ? . (4.22d) Plugging in equation (4.16) and simplifying leads to: M dM M M P dP t t ? ? ? ? ? ? ? + ? ?= 2 2 2 1 1 )1( ? , (4.23a) which is the same as the difference form for total pressure given in compressible flow texts (Saad, 1993). Equation (4.23a) can be written in its differential form: dx dM M P M M dx dP tt ? ? ? ? ? ? ? + ? ?= 2 2 2 1 1 )1( ? . (4.23b) 68 4.2.6 Derivation of the Differential Equation for Entropy Change (ds) The derivation of the differential equation for the entropy change (ds) is related to the equation for the change of entropy for a perfect gas: p dp R T dT Cds p ?= . (4.24a) Substituting for C p in terms of R and ? , and plugging equations (4.13c) and (4.16b) for T dT and p dp , respectively, one gets: ( ) M dM M MR ds ? ? ? ? ? ? ? + ? = 2 2 2 1 1 1 ? . (4.24b) By comparing equation (4.24b) to equation (4.23a), one observes that: t t P dP R ds ?= , (4.24c) or: () t t p P dP C ds ? ? 1? ?= . (4.24d) This relation is the same as the difference form for entropy change given in compressible flow texts (Saad, 1993). 69 4.3 Dimensionless Forms of the Differential Equations Due to the small length scale of this particular geometry, it is practical to nondimensionalize the equations using the inlet conditions rather than the Mach 1 state. Denoting the inlet conditions with subscript ?i?, the dimensionless variables are given by: h D x X = + , i M M M = + , i P P P = + , i T T T = + , (4.25) i V V V = + , i ? ? ? = + , i t t t P P P = + . The dimensionless forms of the differential equations are then: )1(2 ) 2 )1( 1( Re 96 2 23 2 22 + ++ + + ? ? + = MM MMMM dx dM i i D i h ? ? , (4.26) + + + + + + + + ? ? ? ? ? ? ? ? ? ? ? ? ? + ?+ ?= M P MM MM dx dM dx dP i i ) 2 )1( 1( ))1(1( 2 2 2 2 ? ? , (4.27) ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ?= + ++ + + + + ) 2 )1( 1( )1( 2 2 2 MM TMM dx dM dx dT i i ? ? , (4.28) 70 ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? = + + + + + + + 2 2 2 1 1 1 MM M V dx dM dx dV i ? , (4.29) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= + + + + + + dx dV Vdx d ?? , (4.30) + + + + + + + + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ?= dx dM M P MM MM dx dP t i it 2 2 2 2 2 1 1 )1( ? , (4.31) () + + ? ?= t t p P dP C ds ? ? 1 . (4.32) The above equations are similar to the equations for Fanno flow in a constant-area channel except that they are nondimensionalized using the inlet conditions instead of the Mach 1 state. Given a set of initial conditions at a specific x + location, these relations can be integrated in the flow direction (x + increasing). The results of the numerical integration can then be conveniently compared to the analytical relations for this specific flow. In order to verify the accuracy of the numerical code and corresponding numerical results, analytical solutions to the above equations (4.26 through 4.32) are determined following a procedure described by Saad (1993). The analytical solutions are: 2 2 2 )1(2 )1(21 + + + ?+ ?+ = MM M M P i i ? ? , (4.3a) 71 2 2 2 )1(2 )1(2 + + ?+ ?+ = MM M T i i ? ? , (4.33b) 2 2 2 )1(2 )1(2 + ++ ?+ ?+ = MM M MV i i ? ? , (4.3c) 2 2 )1(2 )1(211 2 i i M MM MV ?+ ?+ == + ++ + ? ? ? , (4.33d) ( ) ())12 1 2 2 2 )1(2 )1(21 ? + + + + ? ? ? ? ? ? ? ? ? ? ?+ ?+ ? ? ? ? ? ? ? ? = ? ? ? ? i i t M MM M P , (4.3e) () () ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?+ ?+ ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? + ++ + ? ? ? ? 2 1 2 2 2 2 2 12 121 ln MM M M M C ss i i p i . (4.33f) Having derived the differential equations for the channel flow properties, the modified friction coefficient (Appendix A), F, is included in the differential equation for the Mach number (Eq. 4.26) to account for the presence of the seat rings. The differential equation for the Mach number then becomes: )1(2 ) 2 )1( 1( Re 96 2 23 2 22 + ++ + + ? ? + ? ? ? ? ? ? ? ? = MM MMMFM dx dM i i D i h ? ? . (4.33g) 72 4.4 Radial Disk Flow Model The radial disk flow model is posed as steady, axisymmetric, compressible flow of a perfect gas between two insulated, parallel disks flowing radially toward an outlet hole at the center of the bottom disk (Figure 4.5). This geometry is a more realistic approximation of the actual microvalve geometry. The presence of the seat rings is accounted for by employing a modified friction correlation to determine an apparent friction coefficient. Although the data of Luy et al. (1991) were developed for a channel flow geometry, it is reasonable to also apply the correlation of Appendix A to the radial flow model. This is considered reasonable in this case because in taking infinitesimal steps in the radial direction the cross-sectional area is not changing drastically. 4.4.1 Governing Equations Similar to the channel flow model, the following assumptions are made: 1. Steady flow, 2. One-dimensional flow (properties depend on the radial coordinate, r), 3. Axisymmetric geometry, 4. Perfect gas with constant specific heats, 5. Adiabatic flow with no external work. In order to maintain harmony with the derivations of section 4.2, the continuity, momentum, and energy equations are left in terms of the distance in the flow direction (x) and are later converted to the radial coordinate, r (Figure 4.5). Since the cross-sectional area is changing, the continuity equation is: AV? =constant, (4.34a) 73 or in difference form: 0=++ A dA V dVd ? ? . (4.34b) The momentum equation is: ()AdVVF S xx vv ?= ?? ? ? , (4.35) and it will be discussed in greater detail below. The energy equation is: =+ 2 2 V h constant, (4.36a) or in difference form: 0=+VdVdh . (4.36b) 4.5 Derivation of Differential Expressions for Flow Properties The derivation of the expressions for flow properties is similar to the derivation of the Fanno flow through a channel. Consider the differential control volume shown in Figure 4.6 that extends from x to x + dx. The coordinates x and x + dx correspond to r and r ? dr, respectively. Listed below are some appropriate definitions for the concentric radial disks geometry under consideration. Perimeter = 2(2?r). (4.37a) Wall Surface Area, ()dx D A dxperimeterA h s 4 == . (4.37b) The hydraulic diameter is defined as: 74 ( ) () h r hr perimeter A D h 4 22 2244 = ? ? ? ? ? ? ? ? == ? ? , (4.37c) where A = 2?r(2h) is the cross-sectional area and h is half the distance between the disks. Assuming laminar Stokes flow of an incompressible fluid flowing radially between two concentric disks, a parabolic velocity distribution will be observed (Appendix B). The appropriate relation for the wall shear stress (? f ) is (Appendix B): h V dr dp h dz dV hz r f ? ?? 3 === += , (4.38a) with ? being the fluid viscosity. The friction factor ( ) is then: f h Dh f Vh V hV V f Re 96 Re 2424)3(8 4 22 2 1 ===== ? ? ? ? ? ? , (4.38b) with V being the average velocity between the two disks. 4.5.1 Derivation of the Differential Equation for the Mach number (M) Accounting for all the forces acting on the 1-D control volume (Figure 4.7), the momentum equation in the x- direction is: ()() AVdVAdAAdPPPA sf ?? =?++? . (4.39a) Eliminating the 2 nd order terms, this simplifies to: AVdVAAdPPdA sf ?? =??? . (4.39b) Plugging in for surface area (A s ) from equation (4.37b) and for ? f from the definition of friction factor, gives: 75 AVdVA D dx fVAdPPdA h ?? =??? 2 2 1 . (4.0) Bringing all the terms to one side and dividing by 2 2 M V P ? ? = and cross-sectional area, (A), one gets: 0 22 2 1 =+++ V dV M D dx fM P dP A dA h ?? . (4.1) Substituting for P dP from the perfect gas relation (4.7b), one gets: 0 22 2 1 =++++ V dV M D dx fM T dTd A dA h ?? ? ? . (4.2) Using the continuity equation (4.34b) to eliminate ? ?d and A dA gives: 0 22 2 1 =+++? V dV M D dx fM T dT V dV h ?? . (4.3) Substituting for V dV from the Mach number equation (4.6c) yields: 0 2 2 22 2 1 2 1 =++++? T dTM M dM M D dx fM T dT M dM h ? ?? . (4.44) Assuming adiabatic flow, the term T dT is borrowed from equation (4.13c) and upon simplifying, one gets: 76 22 2 2 1 2 1 1 1 MM M M dM D dx f h ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? = . (4.5) This then leads to the differential equation for the Mach number of the fluid: )1(2 2 1 1 2 23 M M D f M dx dM h ? ? ? ? ? ? ? ? + = ? ? , (4.6) that is identical to equation (4.15). 4.5.2 Differential Equations for Temperature (T), Stagnation Pressure (P t ) and Entropy Change (ds) The derivations of dx dT , dx dP t and p C ds are identical to those derived for the channel flow model. These are: () ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ?= 2 2 1 1 1 M MT dx dM dx dT ? ? , (4.7) dx dM M P M M dx dP tt ? ? ? ? ? ? ? + ? ?= 2 2 2 1 1 )1( ? , (4.8) () t t p P dP C ds ? ? 1? ?= . (4.9) 77 4.5.3 Derivation of the Differential Equation for Density (?) From the continuity equation (4.34b): A dA V dVd ??= ? ? , (4.50a) or: dx dA Adx dV Vdx d ? ? ? ? ? ? ?? ? ? ? ? ? ?= ??? . (4.50b) 4.5.4 Derivation of the Differential Equation for Velocity (V) The derivation of dx dV starts with the form of the momentum equation in (4.43), that is: 0 22 2 1 =+++? V dV M D dx fM T dT V dV h ?? . Plugging in for T dT and h D dx f 2 1 from equations (4.13c) and (4.45), respectively, and simplifying, one gets, ? ? ? ? ? ? ? ? ? ? ? ? ? + = 2 2 1 1 1 M M dM V dV ? , (4.51a) or: ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? = 2 2 1 1 1 M M V dx dM dx dV ? . (4.51b) 78 This relation for the radial flow model is identical to equation (4.17d) for the channel flow model. 4.5.5 Derivation of the Differential Equation for Pressure (P) The derivation of dx dP starts with the form of the momentum equation in (4.41): 0 22 2 1 =+++ V dV M D dx fM P dP A dA h ?? . Substituting for V dV and T dT from the Mach number equation (4.6c), and equation (4.13c), respectively, one gets: () 0 2 1 1 1 2 2 2 122 2 1 = ? + ? ?+++ M MdMM M dM M D dx fM P dP A dA h ? ?? ?? . (4.52) Substituting for h D dx f 2 1 from equation (4.45) and simplifying, gives: () 2 2 2 1 1 11 M M M dM A dA P dP ? + ?+ ? ? ? ? ? ? ? ? ??= ? ? , (4.53a) or in the differential equation form: () ? ? ? ? ? ? ? ? ? ? ? ? ? + ?+ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= 2 2 2 1 1 11 M M dx dM M P dx dA A P dx dP ? ? . (4.53b) 79 Note the similarity of this relation to equation (4.16c), except for the first term on the right hand side that accounts for the area change in the radial model. 4.6 Conversion of Differential Equations in Terms of Radial Coordinate (r) Differential equations for flow properties have been derived in terms of distance in the flow direction, i.e. x. For the radial model, conversion of the independent variable x into radial coordinate, r is more convenient. The relations for conversion are shown below. Considering Figure 4.5: rhhrA ?? 4)2(2 == . (4.54a) Differentiating with respect to r, yields: h dr dA ?4= . (4.54b) The two independent variables r and x are related through: r = (constant ? x), (4.55a) with the corresponding difference equation: dr = -dx . (4.55b) Using the above relation, equations (4.46), (4.47), (4.48), (4.49), (4.50b), (4.51b), and (4.53b) become: )1(2 2 1 1 2 23 M M D f M dr dM h ? ? ? ? ? ? ? ? + ?= ? ? , (4.56) 80 () ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ?= 2 2 1 1 1 M MT dr dM dr dT ? ? , (4.57) dr dM M P M M dr dP tt ? ? ? ? ? ? ? + ? ?= 2 2 2 1 1 )1( ? , (4.58) () t t p P dP C ds ? ? 1? ?= , (4.59) ? ? ? ? ? ? ?? ? ? ? ? ? ?= rdr dV Vdr d ??? , (4.60) ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? = 2 2 1 1 1 M M V dr dM dr dV ? , (4.61) () ? ? ? ? ? ? ? ? ? ? ? ? ? + ?+ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= 2 2 2 1 1 11 M M dr dM M P r P dr dP ? ? . (4.62) 4.7 Dimensionless Form of the Differential Equations Similar to the channel flow model, the equations are nondimensionalized using the inlet conditions at the outer edge of the disks rather than the Mach 1 state. Denoting the inlet conditions with subscript ?i?, the dimensionless variables are identical to those for the channel flow model (Eq. 4.25), except for r + that is defined as: 81 i r r r = + , (4.63) where r i is the radial coordinate of the incoming fluid station (Figure 4.5). Introducing the modified friction coefficient, F, the dimensionless forms of the differential equations are then: )1(2 ) 2 )1( 1( Re 96 2 23 2 22 + ++ + + ? ? + ? ? ? ? ? ? ? ? ?= MM MMMF D r M dr dM i i Dh i i h ? ? , (4.64) + + + + + + + + + + ? ? ? ? ? ? ? ? ? ? ? ? ? + ?+ ??= M P MM MM dr dM r P dr dP i i ) 2 )1( 1( ))1(1( 2 2 2 2 ? ? , (4.65) ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ?= + ++ + + + + ) 2 )1( 1( )1( 2 2 2 MM TMM dr dM dr dT i i ? ? , (4.6) ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? = + + + + + + + 2 2 2 1 1 1 MM M V dr dM dr dV i ? , (4.67) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= + + + + + + + + rdr dV Vdr d ??? , (4.68) + + + + + + + + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ?= dr dM M P MM MM dr dP t i it 2 2 2 2 2 1 1 )1( ? . (4.69) () t t p P dP C ds ? ? 1? ?= (4.70) 82 83 The differential equations for the flow properties for the radial flow model are similar to the equations for the channel flow model with only three exceptions. These are: 1. The area-change term in the pressure equation (1/r + ), 2. The area-change term in the density equation (1/r+), 3. Geometric ratio term in the Mach number relation (r i /D h ). It should be noted that the Reynolds number will also change in the flow direction due to the variation of area. 4.8 Closure Two simplified 1-D flow models have been developed. The channel flow model was developed first to serve as a benchmarking tool for a computer code. The more relevant radial flow model is the main contribution of this research study. Through numerical solutions of the coupled nonlinear ordinary differential equations, the variation of the flow properties between the inlet and outlet ports can be elucidated. The solutions of the equations for the two models are presented in Chapter 5. Figure 4.1: The 2-D channel model of the microvalve geometry 84 Figure 4.2: The 1-D channel model of the microvalve geometry 85 Figure 4.3: Control volume for the changes in flow properties of the 1-D channel model 86 Figure 4.4: Forces acting on the control volume of the 1-D channel model 87 Figure 4.5: Schematic diagram of the radial flow model displaying radial locations, and the directions of x and r 88 Figure 4.6: Control volume for the changes in flow properties of the 1-D radial model for (a) side view and (b) top view 89 Figure 4.7: Forces acting on the control volume of the 1-D radial model 90 91 CHAPTER 5 PREDICTIONS OF THE FLOW PROPERTIES UTILIZING THE MODELS In Chapter 4, using the radial and channel flow models, systems of simultaneous nonlinear ordinary differential equations were derived for various flow properties. These flow properties are: - Mach number - velocity - temperature - density - entropy - total pressure - static pressure. In addition to these variables, related quantities such as the Reynolds number and the Knudsen number can change throughout the flow field and expressions for these were also reported. The channel flow model was developed for the sole purpose of benchmarking the accuracy of the computer code through comparison with analytical expressions for the Fanno flow through a constant-area duct (e.g. Saad, 1993). On the other hand, the more realistic radial flow model can be used to determine the variations of the flow properties within the JPL microvalve. In this Chapter, details of the numerical methods, benchmarking information, and comprehensive results of the predictions based on the radial flow model are presented. 92 5.1 Numerical Analysis A 4 th order Runge-Kutta method was employed for numerical integration of the system of simultaneous nonlinear ordinary differential equations. This method was utilized because of its high accuracy and convenience involving initial value problems. The Runge-Kutta algorithm can be applied to an arbitrary number of coupled ordinary differential equations and only requires the appropriate initial conditions, beginning and end points, and step size of the independent variable. Examples of the Runge-Kutta algorithm can be found in Burden and Faires (1997) and White (1991). 5.1.1 Initial Conditions Since all the flow properties except entropy change were nondimensionalized by the initial (i.e. inlet) values (e.g. ? + = ?/? i ), all the associated initial conditions are set equal to unity. Due to the relative nature of the entropy change, this quantity was set to zero at the inlet station. Additionally, in order to integrate the differential equations (Equations 4.26 ? 4.32 and 4.64 ? 4.70), initial condition quantities for the following are required: - initial Mach number (M i ) - ratio of specific heats (?) - modified friction coefficient (F i ) - Reynolds number (Re Dhi ). The working fluid is helium and its properties were taken from Tsederburg et al. (1971). The initial Reynolds number and the initial Mach number are calculated using the following known quantities (Appendix C): - initial velocity (V i ) - initial density (? i ) - hydraulic diameter (D h ) - viscosity of the fluid (?). The following terms were needed to calculate density from the perfect gas equation and hydraulic diameter: - initial total pressure (P ti ) - gas constant (R) - initial temperature (T i ) - spacing between the disks (H = 2h). The initial modified friction coefficient (F i ) is calculated from the following quantities (Appendix C): - dimensionless ring height (E=e/H) - initial Reynolds number (Re Dhi ). In addition, the initial Mach number and the Reynolds number values are used to initialize the Knudsen number ( 2 1 Re ? = iD ii h MKn ). The initial flow conditions were extrapolated from the experimental results reported by Yang et al. (2004) using digitization software (DigXY, 2003) and a Least- Squares Fit (LSF) method of curve-fitting (see Figure 5.1, Section 1 of Appendix C). The JPL microvalve successfully operates at pressures up to 1000 psig, however the flowmeter used could only measure flowrates with an inlet pressure range up to 300 psig. The microvalve has a deflection (shown as displacement in Figure 5.1) range of 1 to 5 ?m from the top of the 10 ?m seat rings. With the effect of the seat rings accommodated in the model by the modified friction factor, the gap between the two disks (H) varies from 11 to 15 ?m. Therefore, 15 sets of initial values were determined for each of the deflections at three inlet total pressures of 100, 200 and 300 psig. The initial flow conditions are tabulated in Table 5.1. In addition, see Appendix C for an illustration of the calculation of the initial flow conditions. 93 Table 5.1 Measured (Yang et al., 2004) and extrapolated initial flow conditions P ti [psig] H [?m] Q [SCCM] Q [ACMM] V i [m/s] M i Re Dhi EDF i 11 12.55 1.84E-06 0.1483 0.00015 0.1892 0.9091 24.7934 23.9907 12 50.00 7.35E-06 0.5416 0.00054 0.7540 0.8333 22.7273 17.9116 13 111.30 1.64E-05 1.1128 0.00110 1.6785 0.7692 20.9790 13.9951 14 196.68 2.89E-05 1.8260 0.00181 2.9660 0.7143 19.4805 11.3339 15 306.06 4.5E-05 2.6520 0.00263 4.6154 0.6667 18.1818 9.4464 11 21.54 1.58E-06 0.1273 0.00013 0.3249 0.9091 24.7934 23.9946 12 114.14 8.39E-06 0.6182 0.00061 1.7213 0.8333 22.7273 17.9319 13 278.54 2.05E-05 1.3924 0.00138 4.2004 0.7692 20.9790 14.0365 14 514.73 3.78E-05 2.3894 0.00237 7.7623 0.7143 19.4805 11.3977 15 822.72 6.05E-05 3.5645 0.00354 12.4068 0.6667 18.1818 9.5329 11 25.00 1.23E-06 0.0985 0.00010 0.3770 0.9091 24.7934 23.9960 12 138.71 6.8E-06 0.5008 0.00050 2.0917 0.8333 22.7273 17.9397 13 342.29 1.68E-05 1.1408 0.00113 5.1618 0.7692 20.9790 14.0523 14 635.75 3.12E-05 1.9674 0.00195 9.5872 0.7143 19.4805 11.4221 15 1019.09 4.99E-05 2.9435 0.00292 15.3681 0.6667 18.1818 9.5660 300 100 200 It is noted that the Mach numbers are extremely low and the initial values of Re Dhi are of the order of 15 or less. The enhanced effect of friction due to the presence of the rings at the inlet station is estimated to be as high as 24 times the case without the rings. 5.2 Benchmarking of the Computer Code for the Channel Flow Model The channel model was used to benchmark the computer code against analytical expressions for the Fanno flow through a constant-area channel (e.g. Saad, 1993). An analytical expression for the maximum channel length (L max ) was derived following a procedure described by Saad (1993). L max is defined as the distance a gas at Mach number M must travel to reach M = 1. Equation 4.14 was integrated with the help of the method of partial fractions to get the relation (Appendix E): 94 () ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + + +== ? +? + 2 2 1 2 2 11 96 Re max 12 1 ln 2 2 max i i M M D L M M L i ih D h ? ? ? ? ? . (5.1) + max L was solved for the given initial Mach number (subsonic) and this quantity was assigned in the code as the end point of the integration. Therefore, the predicted Mach number at can be checked to determine whether it is equal to the expected value of Mach 1. In order for this particular check to be valid, a shrinking step size was necessary and this feature is discussed in detail in the next section. + max L In addition to the end point check, at each point during the updating of the independent variable (x + ), the predicted values for the flow properties were compared to the respective analytical values given by the expressions of Chapter 4 (Equations 4.33a ? f). The error percentages between the numerical and analytical values for all the flow properties were extremely small. Another realizability condition related to the continuity relation for the channel model that must be satisfied at every point is: .1= ++ V? (5.2) This relation was strictly satisfied in the channel flow simulations. Having satisfied these three benchmarking criteria, the computer code was confirmed to be accurate and viable. 95 5.2.1 Variable Step Size It was determined during benchmarking of the computer code that a uniform step size will not be sufficient in producing accurate results. As part of verifying the predictions of the code against analytical expressions, the code was run for various cases for a maximum channel length (L max ) that corresponds to the end point at which Mach 1 is reached. Using a uniform step size, the L max point was corresponding to a Mach number slightly less (2%) to drastically less (15%) than Mach 1, depending on the particular case. This was determined to be due to the exponential-like trend of the Mach number increase in that close to L max the change in Mach number was increasingly sensitive to slight changes in position. This situation was remedied by employing a shrinking step size. The step size (?x + ) was derived from a geometric series and it shrinks at a user-specified percentage of the previous step size, i.e: ++ + ?=? i i xx ? 1 , (5.3) with constant ? chosen to be less than 1. The first step size (?x + 1 ) necessary for the series of steps to end at the user- specified end point is: ( ) )1( )1( 11 N xxx N ? ? ? ? ?=? +++ , (5.4) where x + N is the ending point of integration, x + 1 is the beginning point, ? is the percentage of previous step size and N is the total number of points. 96 The utilization of a shrinking step size resulted in a Mach number accurately (<<1%) corresponding to 1 at the L max point of integration. Through trial and error, the most effective ? value was 0.96 for all 15 cases. 5.3 Radial Flow Model Results and Discussions The numerical predictions of the radial flow model are reported in graphical form for the following flow properties in terms of radial location, r + : - Reynolds number (Re Dh ) - change in entropy (ds) - Mach number (M + ) - velocity (V + ) - temperature (T + ) - total pressure ( ) + t P - static pressure (P + ) - density (? + ) - Knudsen number (Kn). The beginning and end points of integration were r i = 3000 ?m and r o = 100 ?m, respectively. The number of steps used was 1,000, whereas ? was set to 0.96. Increasing the number of points more than 1,000 had no appreciable effect on the results, so using 1,000 points was deemed satisfactory. 5.3.1 Variation of the Reynolds Number (Re Dh ) 97 The dependence of the local Reynolds number on the radial position for initial total pressures of 100 psig, 200 psig, and 300 psig are shown in Figures 5.2, 5.3, and 5.4, respectively. The Reynolds number is based on the hydraulic diameter and varies with changes in density and velocity at each location. Starting at the respective initial values at r + = 1 tabulated in Table 5.1, the Reynolds number is consistently shown to rise as the gas nears the outlet at r + o = r o /r i = 1/30 = 0.033. Through simple algebraic manipulation it can be shown that the Reynolds number varies inversely with the radial position, r + (see Appendix D). Through close agreement of the predicted Reynolds number with the analytical Reynolds number relation, the simulated Reynolds number results are justified. 5.3.2 Variation of the Change in Entropy (ds) For the radial flow model, flow was assumed to be adiabatic and is irreversible due to friction. Therefore, according to the second law of thermodynamics there must be an entropy rise in this system. The variations of the change in entropy with the radial position in Figures 5.5, 5.6 and 5.7 clearly show this trend. The observed entropy rise is directly related to the incurred loss of the total pressure that will be discussed in detail below. 5.3.3 Variation of the Mach Number (M + ) and Velocity (V + ) The Mach number is defined as the fluid velocity divided by the local speed of sound, so variations of the Mach number (M + ) and velocity (V + ) are considered together in this section. In interpreting the predicted Mach number or velocity variations (Figures 5.8-5.13), consider Equation (4.56). Note that the direction of the radial coordinate is opposite to the flow direction. The hydraulic diameter, gas constant (?), friction factor (f), and M are always positive. Therefore, the (1 ? M 2 ) term determines the sign of dr dM . If the Mach number is less than one, the slope is positive and if the Mach number is greater than one the slope is negative. This shows that the presence of friction causes the 98 Mach number tend toward one. Once a Mach number of one is reached, changes downstream no longer affect upstream conditions and flow is considered ?choked? (Oosthuizen and Carscallen, 1997). The operating Mach numbers of the current radial flow cases were subsonic. Therefore, a gradual increase of the Mach number in the flow direction is expected (Figures 5.8-5.10). Due to the negligible change in the gas temperature (section 5.3.4), dimensionless velocity shows similar trends (Figures 5.11- 5.13). The gas is observed to accelerate in the flow direction, however this velocity rise is not significant. In each of the Figures 5.11-5.13, the variation of the dimensionless velocity for an incompressible fluid is also shown. Using continuity, one can easily show that: + + = = r V 1 const.? , (5.) suggesting that an incompressible fluid will dramatically accelerate in the radial direction. The marked disparity of the predictions of the compressible model in comparison to the incompressible flow relation points to the need for accounting for the compressibility effects. 5.3.4 Variation of the Temperature (T + ) In interpreting the predictions for the variation of the temperature, consider the energy equation (4.3a). For an ideal gas with constant specific heats (which is an assumption in the radial model), changes in enthalpy are proportional to changes in temperature. Therefore, the energy equation indicates that for an accelerating flow there is a corresponding decrease in enthalpy or temperature. Conversely, for a decelerating 99 flow there is a corresponding increase in enthalpy. For the radial flow cases, the flow is subsonic and accelerating. Therefore, as shown graphically, there is an extremely mild decrease in temperature (Figures 5.14-5.16). 5.3.5 Variation of the Total Pressure ( ) + t P The total or stagnation pressure of a moving fluid is defined as the static pressure (P) plus its dynamic pressure (?V 2 /2). In interpreting the predictions for the variation of the total pressure (Figures 5.17, 5.18, and 5.19), consider the difference equation 4.24c (same as Eq. 4.49). Since the change in entropy is always positive, this equation indicates that there is a decrease in the total pressure regardless of whether flow is subsonic or supersonic. The numerical predictions of the total pressure drop through the passageway are compared to the total pressure drop relation (Eq. B.23) for Stokes flow from Appendix B in Table 5.2. 100 Table 5.2 Total Pressure Drop ( ) for the Radial Flow Model Compared to the Incompressible Stokes Model oi ttt PPP ??? + 11 1.844E-06 2.1131 11.9311 2.3356 12 7.350E-06 8.4209 27.6364 7.1983 13 1.636E-05 18.7456 38.1025 12.6963 14 2.891E-05 33.1248 43.8836 18.1660 15 4.499E-05 51.5460 46.4365 23.3444 11 1.583E-06 3.5625 9.5476 2.0067 12 8.390E-06 18.8764 29.3631 8.2697 13 2.047E-05 46.0634 44.4654 16.1761 14 3.783E-05 85.1237 53.8081 24.6497 15 6.047E-05 136.0571 58.8977 33.3647 11 1.225E-06 3.9854 7.2085 1.5531 12 6.797E-06 22.1099 23.1575 6.7128 13 1.677E-05 54.5612 35.4285 13.3260 14 3.115E-05 101.3395 43.1155 20.5223 15 4.994E-05 162.4447 47.4225 28.0650 Numerical ?P t [kPa] Stokes ?P t [kPa] P ti [psig] H [?m] Q actual [m 3 /min] Mass Flowrate [mg/min] 200 300 100 It is clear that the total pressure drop for Stokes incompressible flow is consistently lower than the predictions of the radial model that accounts for the compressibility effect. The numerical predictions of the total pressure drop in comparison to the Stokes relation are shown in Figure 5.20. For an imposed total pressure at the inlet station (P ti = constant), both approaches show that as the spacing is changed from 11 to 15 ?m, the incurred total pressure drop is observed to rise with the mass flowrate. Note that the percentage difference between the predictions and Stokes relation diminishes as the spacing is raised. This is expected since the effect of the ring becomes more negligible as the two disks move apart and the Stokes model becomes more realistic. 101 102 5.3.6 Variation of the Static Pressure (P + ) Figures 5.21, 5.22, and 5.23 show the static pressure increasing with the decrease of the radial coordinate, r + . The pressure rise in the flow direction is observed to be weakly affected by the initial conditions or the gap between the two disks. Also shown on these figures are the pressure variation of an ideal incompressible fluid (Bernoulli?s equation) moving radially inward (Section F.1 of Appendix F). In contrast to the pressure rise trends predicted by the compressible model, the Bernoulli?s equation predicts that the pressure will decrease in the flow direction. In order to address these contradicting behaviors, two other simplified models for the static pressure variation were derived in Appendix F. These are: 1. Ideal Compressible Flow (Section F.2), 2. Stokes Flow (Section F.3). For the present simulated cases for which M i ?s are extremely small, P + ? 1. The predictions of the fluid static pressure within the passageway using the Stokes flow model are shown in Figure 5.24 for the cases with P ti = 300 psig and H = 11 and 15 ?m. The more realistic Stokes model that accounts for fluid viscosity predicts the static pressure rise that is in concert with the numerical integrations. The Stokes flow relation for the static pressure rise (Eq. B.13) through the passageway is compared to the numerical pressure rise predictions in Table 5.3. Table 5.3 Static Pressure Rise (?P = P o -P i ) for the Radial Flow Model Compared to the Incompressible Stokes Model 11 2.1131 2.26E+04 2.3328 12 8.4209 2.21E+04 7.1607 13 18.7456 2.18E+04 12.5374 14 33.1248 2.16E+04 17.7381 15 51.5460 2.15E+04 22.4419 11 3.5625 4.26E+04 2.0026 12 18.8764 4.20E+04 8.1734 13 46.0634 4.16E+04 15.6875 14 85.1237 4.13E+04 23.2109 15 136.0571 4.12E+04 30.1630 11 3.9854 6.27E+04 1.5495 12 22.1099 6.22E+04 6.6215 13 54.5612 6.19E+04 12.8519 14 101.3395 6.16E+04 19.1120 15 162.4447 6.15E+04 24.9083 Numerical [kPa] 100 Mass Flowrate [mg/min] 200 300 Static ?P Stokes [kPa] P ti [psig] H [?m] The compressible model predictions for static pressure rise through the passageway are consistently much higher than that of the Stokes flow relation. This marked difference between the two approaches further reinforces the need for accounting for the compressibility effects. 5.3.7 Variation of the Density (? + ) Considering the ideal gas relation, density varies linearly with pressure. However, density varies inversely with temperature, T. Due to the predicted negligible temperature change, density exhibits a similar trend to that of static pressure (Figures 103 5.25-5.27). The drastic increase in density as compared to the relatively small increase in the Mach number suggests that the fluid is compressed markedly more than it is being accelerated. Using the continuity relation, a simple realizability condition must be satisfied by the predicted dimensionless density and velocity values at every radial location. This relation is: + ++ = r V 1 ? . (5.6) This relation was shown to hold very well. 5.3.8 Variation of the Knudsen Number (Kn) The Knudsen number is the ratio of the mean free path of a gas molecule to the characteristic length scale of a given flow problem. The Knudsen number can be defined in terms of the Mach number and Re. The Knudsen number provides a method for determining the validity of the continuum flow assumption or whether a rarefied flow must be considered. Continuum flow is valid when gas molecules are close enough together that they collide with one another frequently and the gas acts as a continuous fluid. The criterion for the continuum flow assumption for large Re is (John, 1984): 01.0 Re <= M Kn . (5.7) Variation of the local Knudsen number throughout the flow is given in Figures 5.28-5.30. It is observed that the Kn is decreasing nonlinearly. This decrease represents a shortening of the mean free path. This is consistent with an increase of fluid density in the flow direction. Density rise forces the gas molecules closer together, thus reducing 104 105 the mean free path. Therefore, the continuum assumption remains valid, and it is not necessary to consider a slip flow condition on the walls. 5.4 Analysis of the Total Pressure Losses Beyond the Seat Rings The radial model is capable of predicting various flow properties as the gas flows between the two disks. Further analysis was done to determine the total pressure drop through the outlet tube of the microvalve for the conditions of incompressibility and compressibility. First, an equivalent length was calculated for the outlet tube shown in Figure 3.4a (Fox et al., 2004). Then, the total pressure drop for internal incompressible viscous flow (Fox et al., 2004) was found. In addition, the channel model code was used to perform the Fanno analysis through the circular outlet tube. The results for both incompressible and compressible analyses are shown in Table 5.4, and a sample calculation for the incompressible total pressure drop is given in Sections 5.4.1. Table 5.4 Total Pressure Drop (?P t ) Through the Outlet Tube (Incompressible and Compressible Fanno Analyses) 11 1.844E-06 2.113062 3.385E+01 3.311E-02 3.285E-05 1.15E+01 1.345E-02 5.086E-04 12 7.350E-06 8.420895 3.317E+01 1.347E-01 1.336E-04 4.57E+01 5.471E-02 2.068E-03 13 1.636E-05 18.74555 3.271E+01 3.040E-01 3.015E-04 1.02E+02 1.235E-01 4.669E-03 14 2.891E-05 33.12479 3.246E+01 5.413E-01 5.370E-04 1.80E+02 2.199E-01 8.314E-03 15 4.499E-05 51.54602 3.235E+01 8.452E-01 8.384E-04 2.80E+02 3.434E-01 1.298E-02 11 1.583E-06 3.562472 6.706E+01 2.818E-02 2.795E-05 1.94E+01 1.145E-02 4.125E-04 12 8.390E-06 18.87638 6.616E+01 1.514E-01 1.501E-04 1.03E+02 6.149E-02 2.216E-03 13 2.047E-05 46.06344 6.547E+01 3.732E-01 3.702E-04 2.50E+02 1.516E-01 5.464E-03 14 3.783E-05 85.12367 6.505E+01 6.943E-01 6.887E-04 4.62E+02 2.820E-01 1.016E-02 15 6.047E-05 136.0571 6.481E+01 1.114E+00 1.105E-03 7.39E+02 4.524E-01 1.630E-02 11 1.225E-06 3.9854 9.727E+01 2.174E-02 2.156E-05 2.17E+01 8.830E-03 3.226E-04 12 6.797E-06 22.10986 9.655E+01 1.215E-01 1.205E-04 1.20E+02 4.935E-02 1.803E-03 13 1.677E-05 54.56123 9.600E+01 3.015E-01 2.991E-04 2.96E+02 1.225E-01 4.475E-03 14 3.115E-05 101.3395 9.565E+01 5.621E-01 5.575E-04 5.51E+02 2.283E-01 8.341E-03 15 4.994E-05 162.4447 9.546E+01 9.028E-01 8.955E-04 8.83E+02 3.667E-01 1.340E-02 100 200 300 H [?m] Q actual [m 3 /min] Mass Flowrate [mg/min] P ti [psig] Outlet Incomp. ?P t [kPa] Outlet Comp. ?P t [kPa] ? o [kg/m 3 ] V avg [m/s] Mach Number Re D 5.4.1 Sample Calculation for Incompressible Total Pressure Drop of the Outlet Tube For this sample calculation the case of P ti = 100 psig and H = 12 ?m is considered. First, a dimensionless equivalent length (L e /D) for the outlet tube is needed. For the sharp 90? bend, a dimensionless equivalent length of 60 diameters was determined from Figure 8.16 of Fox et al. (2004). The dimensionless equivalent length for the 200 micron diameter outlet tube is then calculated as follows: 13060 2.0 5.105.3 =+ + = mm mmmm D L e . (5.8) The mass flowrate is calculated by multiplying the initial actual volumetric flowrate (Q actual ) and the initial density (? i ). The average velocity (V avg ) for the outlet tube is then 106 calculated using the area of the outlet tube and the density (? o ) at the opening of the outlet given by the value at the endpoint of the radial flow model: () s m m kg mg o avg smg kg A m V 1347.0 60 min1 1 10 0001.017.33 421.8 6 2 min 3 = ? ? ? ? ? ? == ? ? & . (5.9) The initial Mach number in the outlet tube is then given by: .1034.1 )293(~)2077(67.1 1347.0 4? ?=== s m i o avg o TR V M i ? (5.10) The Reynolds number (Re D ) in the tube will remain constant and is calculated as follows: ( ) 75.45 10953.1 0002.01347.017.33 Re 5 3 = ? ? ? ? ? ? ? == ? ms kg s m m kg avgo D m DV ? ? . (5.11) Finally, the total pressure drop (?P t ) that is equal to the static pressure drop is given by: () () .1047.517.33 2 1347.0 130 75.45 64 2 2 2 2 3 2 2 kPa V D L fP m kgs m o avg e t ? ?=? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? =? ? (5.12) In determining the compressible Fanno estimation of ?P t , the equivalent length (L e /D), initial Mach number (M oi ), and Reynolds number (Re D ) were used as inputs to the numerical code. 107 108 The total pressure drops for both the incompressible and compressible results are miniscule compared to the total pressure drop predicted by the radial flow model. This supports the initial hypothesis that most of total pressure drop through the valve occurs as the gas flows over the seat rings. 5.5 Closure Variations of flow properties through the microvalve for the radial flow model were presented. The trends for the derived Stokes flow relation for both static and total pressures were in agreement with the predicted trends based on the radial flow model. Additionally, a comparison of the numerically-predicted total pressure drops to that of the Stokes flow calculations, supports the necessity of accounting for compressibility effects. The total pressure drop through the outlet tube was calculated and determined to be negligible compared to the total pressure drop over the rings. Conclusions and recommendations for future work are discussed in Chapter 6. Figure 5.1: Measured flow rates for an actuated microvalve (Yang et al., 2004) 109 r + Re Dh 00.20.40.60.81 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 H=11?m(M i =0.000147, Re i =0.19) H=12?m(M i =0.000537, Re i =0.75) H=13?m(M i =0.001104, Re i =1.68) H=14?m(M i =0.001811, Re i =2.97) H=15?m(M i =0.002631, Re i =4.62) P ti =100psig Figure 5.2: Variation of Re Dh for an inlet total pressure of 100 psig 110 r + Re Dh 00.20.40.60.81 0 50 100 150 200 250 300 350 H=11?m(M i =0.0001262, Re i =0.33) H=12?m(M i =0.0006132, Re i =1.72) H=13?m(M i =0.001381, Re i =4.20) H=14?m(M i =0.00237, Re i =7.76) H=15?m(M i =0.003536, Re i =12.40) P ti =200psig Figure 5.3: Variation of Re Dh for an inlet total pressure of 200 psig 111 r + Re Dh 00.20.40.60.81 0 50 100 150 200 250 300 350 400 450 H=11?m(M i =0.000098, Re i =0.38) H=12?m(M i =0.000497, Re i =2.09) H=13?m(M i =0.001132, Re i =5.16) H=14?m(M i =0.00195, Re i =9.59) H=15?m(M i =0.00292, Re i =15.37) P ti =300psig Figure 5.4: Variation of Re Dh for an inlet total pressure of 300 psig 112 r + ds / C p 00.20.40.60.81 0 0.005 0.01 0.015 0.02 0.025 H=11?m(M i =0.000147, Re i =0.19) H=12?m(M i =0.000537, Re i =0.75) H=13?m(M i =0.001104, Re i =1.68) H=14?m(M i =0.001811, Re i =2.97) H=15?m(M i =0.002631, Re i =4.62) P ti =100psig Figure 5.5: Variation of the change in entropy for an inlet total pressure of 100 psig 113 r + ds / C p 00.20.40.60.81 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 H=11?m(M i =0.0001262, Re i =0.33) H=12?m(M i =0.0006132, Re i =1.72) H=13?m(M i =0.001381, Re i =4.20) H=14?m(M i =0.00237, Re i =7.76) H=15?m(M i =0.003536, Re i =12.40) P ti =200psig Figure 5.6: Variation of the change in entropy for an inlet total pressure of 200 psig 114 r + ds/ C p 00.20.40.60.81 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 H=11?m(M i =0.000098, Re i =0.38) H=12?m(M i =0.000497, Re i =2.09) H=13?m(M i =0.001132, Re i =5.16) H=14?m(M i =0.00195, Re i =9.59) H=15?m(M i =0.00292, Re i =15.37) P ti =300psig Figure 5.7: Variation of the change in entropy for an inlet total pressure of 300 psig 115 r + M + 00.20.40.60.81 1 1.01 1.02 1.03 1.04 1.05 1.06 H=11?m(M i =0.000147, Re i =0.19) H=12?m(M i =0.000537, Re i =0.75) H=13?m(M i =0.001104, Re i =1.68) H=14?m(M i =0.001811, Re i =2.97) H=15?m(M i =0.002631, Re i =4.62) P ti = 100 psig Figure 5.8: Variation of the Mach number for an inlet total pressure of 100 psig 116 r + M + 00.20.40.60.81 1 1.005 1.01 1.015 1.02 1.025 1.03 1.035 1.04 H=11?m(M i =0.0001262, Re i =0.33) H=12?m(M i =0.0006132, Re i =1.72) H=13?m(M i =0.001381, Re i =4.20) H=14?m(M i =0.00237, Re i =7.76) H=15?m(M i =0.003536, Re i =12.40) P ti =200psig Figure 5.9: Variation of the Mach number for an inlet total pressure of 200 psig 117 r + M + 00.20.40.60.81 1 1.005 1.01 1.015 1.02 H=11?m(M i =0.000098, Re i =0.38) H=12?m(M i =0.000497, Re i =2.09) H=13?m(M i =0.001132, Re i =5.16) H=14?m(M i =0.00195, Re i =9.59) H=15?m(M i =0.00292, Re i =15.37) P ti =300psig Figure 5.10: Variation of the Mach number for an inlet total pressure of 300 psig 118 r + V + 00.20.40.60.81 1 1.01 1.02 1.03 1.04 1.05 1.06 H=11?m(M i =0.000147, Re i =0.19) H=12?m(M i =0.000537, Re i =0.75) H=13?m(M i =0.001104, Re i =1.68) H=14?m(M i =0.001811, Re i =2.97) H=15?m(M i =0.002631, Re i =4.62) Incompressible P ti =100 psig Figure 5.11: Variation of the velocity for an inlet total pressure of 100 psig 119 r + V + 00.20.40.60.81 1 1.005 1.01 1.015 1.02 1.025 1.03 1.035 1.04 H=11?m(M i =0.0001262, Re i =0.33) H=12?m(M i =0.0006132, Re i =1.72) H=13?m(M i =0.001381, Re i =4.20) H=14?m(M i =0.00237, Re i =7.76) H=15?m(M i =0.003536, Re i =12.40) Incompressible P ti =200 psig Figure 5.12: Variation of the velocity for an inlet total pressure of 200 psig 120 r + V + 00.20.40.60.81 1 1.005 1.01 1.015 1.02 H=11?m(M i =0.000098, Re i =0.38) H=12?m(M i =0.000497, Re i =2.09) H=13?m(M i =0.001132, Re i =5.16) H=14?m(M i =0.00195, Re i =9.59) H=15?m(M i =0.00292, Re i =15.37) Incompressible P ti =300 psig Figure 5.13: Variation of the velocity for an inlet total pressure of 300 psig 121 r + T + -T + i 00.250.50.751 -80 -70 -60 -50 -40 -30 -20 -10 0 x10 11 ,H=11?m(M i =0.000147, Re i =0.19) x10 10 ,H=12?m(M i =0.000537, Re i =0.75) x10 9 ,H=13?m(M i =0.001104, Re i =1.68) x10 8 ,H=14?m(M i =0.001811, Re i =2.97) x10 8 ,H=15?m(M i =0.002631, Re i =4.62) P ti =100psig Figure 5.14: Variation of the temperature for an inlet total pressure of 100 psig 122 r + T + -T + i 00.250.50.751 -8 -7 -6 -5 -4 -3 -2 -1 0 x10 11 ,H=11?m(M i =0.0001262, Re i =0.33) x10 9 ,H=12?m(M i =0.0006132, Re i =1.72) x10 8 ,H=13?m(M i =0.001381, Re i =4.20) x10 7 ,H=14?m(M i =0.00237, Re i =7.76) x10 7 ,H=15?m(M i =0.003536, Re i =12.40) P ti =200psig Figure 5.15: Variation of the temperature for an inlet total pressure of 200 psig 123 r + T + -T + i 00.250.50.751 -6 -5 -4 -3 -2 -1 0 x10 11 ,H=11?m(M i =0.000098, Re i =0.38) x10 9 ,H=12?m(M i =0.000497, Re i =2.09) x10 8 ,H=13?m(M i =0.001132, Re i =5.16) x10 8 ,H=14?m(M i =0.00195, Re i =9.59) x10 7 ,H=15?m(M i =0.00292, Re i =15.37) P ti =300psig Figure 5.16: Variation of the temperature for an inlet total pressure of 300 psig 124 r + P t+ 00.20.40.60.81 0.94 0.95 0.96 0.97 0.98 0.99 1 H=11?m(M i =0.000147, Re i =0.19) H=12?m(M i =0.000537, Re i =0.75) H=13?m(M i =0.001104, Re i =1.68) H=14?m(M i =0.001811, Re i =2.97) H=15?m(M i =0.002631, Re i =4.62) P ti =100psig Figure 5.17: Variation of the total pressure for an inlet total pressure of 100 psig 125 r + P t+ 00.20.40.60.81 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 H=11?m(M i =0.0001262, Re i =0.33) H=12?m(M i =0.0006132, Re i =1.72) H=13?m(M i =0.001381, Re i =4.20) H=14?m(M i =0.00237, Re i =7.76) H=15?m(M i =0.003536, Re i =12.40) P ti =200psig Figure 5.18: Variation of the total pressure for an inlet total pressure of 200 psig 126 r + P t+ 00.20.40.60.81 0.98 0.985 0.99 0.995 1 H=11?m(M i =0.000098, Re i =0.38) H=12?m(M i =0.000497, Re i =2.09) H=13?m(M i =0.001132, Re i =5.16) H=14?m(M i =0.00195, Re i =9.59) H=15?m(M i =0.00292, Re i =15.37) P ti =300psig Figure 5.19: Variation of the total pressure for an inlet total pressure of 300 psig 127 Mass Flow Rate [mg/min] ? P t [k P a ] 20 40 60 80 100 120 140 160 5 10 15 20 25 30 35 40 45 50 55 P ti = 100 psig P ti = 200 psig P ti = 300 psig Stokes (P ti =100psig) Stokes (P ti =200psig) Stokes (P ti =300psig) Figure 5.20: Comparison of total pressure drop (numerical results and Stokes relation) vs. mass flowrate 128 r + P + 00.20.40.60.81 1 6 11 16 21 26 H=11?m(M i =0.000147, Re i =0.19) H=12?m(M i =0.000537, Re i =0.75) H=13?m(M i =0.001104, Re i =1.68) H=14?m(M i =0.001811, Re i =2.97 H=15?m(M i =0.002631, Re i =4.62) Ideal Incompressible (H = 15 ?m) P ti =100psig <1 Figure 5.21: Variation of the static pressure for an inlet total pressure of 100 psig 129 r + P + 00.20.40.60.81 1 6 11 16 21 26 H=11?m(M i =0.0001262, Re i =0.33) H=12?m(M i =0.0006132, Re i =1.72) H=13?m(M i =0.001381, Re i =4.20) H=14?m(M i =0.00237, Re i =7.76) H=15?m(M i =0.003536, Re i =12.40) Ideal Incompressible (H = 15?m) P ti =200psig <1 Figure 5.22: Variation of the static pressure for an inlet total pressure of 200 psig 130 r + P + 00.20.40.60.81 1 6 11 16 21 26 H=11?m(M i =0.000098, Re i =0.38) H=12?m(M i =0.000497, Re i =2.09) H=13?m(M i =0.001132, Re i =5.16) H=14?m(M i =0.00195, Re i =9.59) H=15?m(M i =0.00292, Re i =15.37) Ideal Incompressible (H = 15 ?m) P ti = 300 psig <1 Figure 5.23: Variation of the static pressure for an inlet total pressure of 300 psig 131 r + P + 00.250.50.751 1 1.0025 1.005 1.0075 1.01 1.0125 1.015 Stokes [Incompressible] 11 ?m Stokes [Incompressible] 15 ?m P ti = 300 psig, H = 11, 15 ?m Figure 5.24: Variations of the static pressure for the Stokes flow relation for an inlet total pressure of 300 psig 132 r + ? + 00.20.40.60.81 5 10 15 20 25 30 H=11?m(M i =0.000147, Re i =0.19) H=12?m(M i =0.000537, Re i =0.75) H=13?m(M i =0.001104, Re i =1.68) H=14?m(M i =0.001811, Re i =2.97) H=15?m(M i =0.002631, Re i =4.62) P ti =100psig Figure 5.25: Variation of the density for an inlet total pressure of 100 psig 133 r + ? + 00.20.40.60.81 5 10 15 20 25 30 H=11?m(M i =0.0001262, Re i =0.33) H=12?m(M i =0.0006132, Re i =1.72) H=13?m(M i =0.001381, Re i =4.20) H=14?m(M i =0.00237, Re i =7.76) H=15?m(M i =0.003536, Re i =12.40) P ti =200psig Figure 5.26: Variation of the density for an inlet total pressure of 200 psig 134 r + ? + 00.20.40.60.81 5 10 15 20 25 30 H=11?m(M i =0.000098, Re i =0.38) H=12?m(M i =0.000497, Re i =2.09) H=13?m(M i =0.001132, Re i =5.16) H=14?m(M i =0.00195, Re i =9.59) H=15?m(M i =0.00292, Re i =15.37) P ti =300psig Figure 5.27: Variation of the density for an inlet total pressure of 300 psig 135 r + Kn = M R e -1 / 2 00.20.40.60.81 0 0.0002 0.0004 0.0006 0.0008 0.001 0.0012 H=11?m(M i =0.000147, Re i =0.19) H=12?m(M i =0.000537, Re i =0.75) H=13?m(M i =0.001104, Re i =1.68) H=14?m(M i =0.001811, Re i =2.97) H=15?m(M i =0.002631, Re i =4.62) P ti =100psig Figure 5.28: Variation of the Knudsen number for an inlet total pressure of 100 psig 136 r + Kn = M R e -1 / 2 00.20.40.60.81 0 0.0002 0.0004 0.0006 0.0008 0.001 H=11?m(M i =0.0001262, Re i =0.33) H=12?m(M i =0.0006132, Re i =1.72) H=13?m(M i =0.001381, Re i =4.20) H=14?m(M i =0.00237, Re i =7.76) H=15?m(M i =0.003536, Re i =12.40) P ti =200psig Figure 5.29: Variation of the Knudsen number for an inlet total pressure of 200 psig 137 r + Kn = M R e -1 / 2 00.20.40.60.81 0 0.0002 0.0004 0.0006 H=11?m(M i =0.000098, Re i =0.38) H=12?m(M i =0.000497, Re i =2.09) H=13?m(M i =0.001132, Re i =5.16) H=14?m(M i =0.00195, Re i =9.59) H=15?m(M i =0.00292, Re i =15.37) P ti =300psig Figure 5.30: Variation of the Knudsen number for an inlet total pressure of 300 psig 138 139 CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS 6.1 Conclusions In order to determine the variations of flow properties over the rings for the conditions of static operation of the microvalve, two 1-D flow models were developed and the results analyzed. The main conclusions and achievements of this study are summarized as follows: 1. A channel flow model was developed for steady, compressible frictional flow of a perfect gas between two infinite insulated parallel plates. This code was implemented in benchmarking the numerical code against analytical expressions for the properties of flow through a constant-area channel. In addition, this model was used to perform a 1-D Fanno analysis of flow through the outlet tube of the microvalve in order to determine the total pressure drop incurred. 2. Based on the data for flow through a channel with rectangular fins mounted on the bottom wall (Luy et al., 1991), a correlation for a modified friction coefficient was derived that accounts for the fins? presence. The modified friction coefficient allowed for a simpler 1-D analysis of flow through the microvalve to be accomplished. 3. A radial flow model was developed for steady, axisymmetric, compressible frictional flow of a perfect gas between two insulated, parallel disks flowing radially 140 toward an outlet hole at the center of the bottom disk. This model was implemented to determine the variation of properties of flow between the boss and seat plates of a JPL microvalve. The most notable conclusion from the flow property trends is that of a drastic increase in density and static pressure in contrast to a rather small increase in the Mach number. One can infer from these findings that the gas is being compressed more than it is being accelerated. Also of importance, the total pressure drop was shown to be significant across the seat rings. 4. A 2-D Stokes flow model was derived for incompressible, axisymmetric, radial flow between two concentric parallel disks. The results of this model were used to verify the flow property variations from the radial flow model. In particular, for the Stokes flow model, relations for radial velocity, average velocity, Darcy friction factor, volumetric flowrate, static pressure rise, and total pressure drop were derived. The Stokes flow model trends for both static and total pressure were in accordance with the radial compressible flow model trends. In addition, a comparison of Stokes flow values for both the static pressure rise and the total pressure drop to that of the numerical results demonstrates the necessity of accounting for compressibility effects. 5. Implementing the channel flow model for the outlet tube, the total pressure loss through the outlet tube was found to be negligible in comparison to that of total pressure loss across the seat rings. Therefore, the hypothesis that most of total pressure drop occurs across the rings between the boss and seat plates is upheld. 141 6.2 Recommendations for Future Work The suggestions for future work are summarized as follows: 1. The variation of the flow properties for the microvalve models involve conditions at relatively low Mach numbers. However, the models are incapable of functioning with Mach numbers approaching 1. An investigation into higher Mach number conditions would be of interest. This could be made possible through the implementation of a perturbation method within the existing numerical code. 2. An investigation into non-continuum flow effects would be of interest. The current two models that are implemented in the computer code can easily handle slip conditions. However, the Knudsen numbers for the microvalve flow conditions did not merit the use of a slip condition on the walls. 3. This thesis focused on modeling the flow through the microvalve for static operating conditions and ignored the complications of the microvalve?s dynamic operation. It would be of value to investigate the dynamic operation through a moving- boundary microvalve model. 4. A 2-D or 3-D CFD analysis of the microvalve would be of interest. Despite the anticipated difficulties of grid generation and compressibility effects involved with this approach, the results would be of great consequence to the understanding of the flow through the microvalve. 142 REFERENCES Ayhan, A. F., ?Design of a Piezoelectrically Actuated Microvalve for Flow Control in Fuel Cells,? (Master?s Thesis) University of Pittsburgh, 2002. Bintoro. J. S., and P.J. Hesketh, ?An Electromagnetic Actuated On/Off Microvalve Fabricated on Top of a Single Wafer,? Journal of Micromechanics and Microengineering, Vol. 15, No. 6, (June 1, 2005): pp. 1157-1173. Burden, R. L. and J. Faires, Numerical Analysis, 6 th ed., Brooks/Cole Pub. Co., Pacific Grove, CA, 1997. DigXY 1.2.1., Computer Digitization Software, Thunderhead Engineering Consultants, Manhattan, KS, 2003. Fox, R. W., A. T. McDonald, and P. J. Pritchard, Introduction to Fluid Mechanics, 6 th ed., John Wiley & Sons, Inc., New York, NY, 2004. Geon, H., H. Kahn, S. M. Phillips, A. H. Heuer, ?Fully-Microfabricated, Silicon Spring Biased, Shape Memory Actuated Microvalve,? Solid-State Sensor and Actuator Workshop, Hilton Head Island, SC, (June 4-8, 2000): pp. 230-233. Huff, M. A., M. S. Mettner, T. A. Lober, and M. A. Schmidt, ?A Pressure-Balanced Electrostatically-Actuated Microvalve,? Technical Digest, 1990 Solid-State Sensor and Actuator Workshop, (1990): pp. 123-127. 143 Janson, S., H. Helvajian, S. Amimoto, G. Smith, D. Mayer, and S. Feuerstein, ?Microtechnology for Space Systems,? IEEE Aerospace Applications Conference Proceedings, Vol. 1, (March 21-28, 1998): pp. 409-418. Jerman, H., ?Electrically-Activated, Micromachined Diaphragm Valves,? Technical Digest, IEEE Solid-State Sensor and Actuator Workshop, (June 4-7, 1990): pp. 65-69. John, J. E. A., Gas Dynamics, 2 nd ed., Prentice-Hall Inc., Upper Saddle River, NJ, 1984. Luy, C-D, C-H Cheng, and W-H Huang, ?Forced Convection in Parallel-Plate Channels with a Series of Fins Mounted on the Wall,? Applied Energy, Vol. 39, No. 2, (1991): pp. 127-144. Messner, S., M. Mueller, V. Burger, J. Schaible, H. Sandmaier, and R. Zengerle, ?Normally-Closed, Bimetallically Actuated 3-Way Microvalve for Pneumatic Applications,? Proceedings of the IEEE MicroElectroMechanical Systems (MEMS), (January 25-29, 1998): pp. 40-44. Mueller, J., ?A Review and Applicability Assessment of MEMS-Based Microvalve Technologies for Microspacecraft Propulsion,? in Micropropulsion for Small Spacecraft, M. Micci and A. Ketsdever, Eds., Reston, VA: AIAA, 2000, Vol. 187, Progress in Astronautics and Aeronautics, Ch. 19. Mueller, J., I. Chakraborty, S. Vargo, C. Marrese, V. White, D. Bame, R. Reinicke, and J. Holzinger, ?Towards Micropropulsion Systems On-A-Chip: Initial Results of Component Feasibility Studies,? IEEE Aerospace Conference Proceedings, Vol. 4, (March 18-25, 2000): pp. 149-168. 144 Ohnstein, T., T. Fukiura, J. Ridley, and U. Bonne, ?Micromachined Silicon Microvalve,? Proceedings of the IEEE MicroElectroMechanical Systems ? An Investigation of Micro Structures, Sensors, Actuators, Machines, (February 11-14, 1990): pp. 95- 98. Oosthuizen, P. H., and W. E. Carscallen, Compressible Fluid Flow, McGraw-Hill, New York, NY, 1997. Ouellette, J., ?A New Wave of Microfluidic Devices,? Industrial Physicist, Vol. 9, No. 4, (August/September 2003): pp. 14-17. Rich, C. A., and K. D. Wise, ?A High-Flow Thermopneumatic Microvalve with Improved Efficiency and Integrated State Sensing,? Journal of MicroElectroMechanical Systems, Vol. 12, No. 2, (April 2003): pp. 201-208. Roberts, D. C., H. Li, J. L. Steyn, O. Yaglioglu, S. M. Spearing, M. A. Schmidt, and N. W. Hagood, ?A Piezoelectric Microvalve for Compact High-Frequency, High- Differential Pressure Hydraulic Micropumping Systems,? Journal of MicroElectroMechanical Systems, Vol. 12, No. 1, (February 2003): pp. 81-92. Saad, M. A., Compressible Fluid Flow, 2 nd ed., Prentice-Hall, Englewood Cliffs, NJ, 1993. Shinozawa, Y., T. Abe, and T. Kondo, ?Proportional Microvalve Using a Bi-Stable Magnetic Actuator,? Proceedings of the IEEE MicroElectroMechanical Systems (MEMS), (January 26-30, 1997): pp. 233-237. Skrobanek, K. D., M. Kohl, and S. Miyazaki, ?Stress-Optimized Shape Memory Microvalves,? Proceedings of the IEEE MicroElectroMechanical Systems (MEMS), (January 26-30 1997): pp. 256-261. 145 Tomonari, S., H. Yoshida, M. Kamakura, K. Yoshida, K. Kawahito, M. Saitoh, H. Kawada, S. Juodkazis, and H. Misawa, ?Efficient Microvalve Driven by a Si-Ni Bimorph,? Japanese Journal of Applied Physics, Part 1: Regular Papers and Short Notes and Review Papers, Vol. 42, No. 7A, (July 2003): pp. 4593-4597. Tsederburg, N. V., V. N. Popov, and N. A. Morozova, Thermodynamic and Thermophyical Properties of Helium, Translated from Russian by IPST staff, Ed. A. F. Alyab?ev, Israel Program for Scientific Translations, Jerusalem, 1971. Warring, R. H., Handbook of Valves, Piping and Pipelines, Gulf Publishing Company, Houston, TX, 1982. White, F. M., Viscous Fluid Flow, 2 nd ed., McGraw-Hill, New York, NY, 1991. White, F. M., Fluid Mechanics, 4 th ed., WCB/McGraw-Hill, New York, NY, 1999. Yang, E-H, Choonsup Lee, Juergen Mueller, and Thomas George, ?Leak-Tight Piezoelectric Microvalve for High-Pressure Gas Micropropulsion,? Journal of MicroElectroMechanical Systems, Vol. 13, No. 5, (Oct. 2004): pp. 799-807. Yang, X., C. Grosjean, Y-C Tai, C-M Ho, ?A MEMS Thermopneumatic Silicon Membrane Valve,? Proceedings of the IEEE MicroElectroMechanical Systems (MEMS), (January 26-30, 1997): pp. 114-118. Websites: http://www.carf-engineering.com/ http://www.spacequest.com/ http://www.sstl.co.uk/ 146 APPENDICES 147 Appendix A Correlation for the Modified Friction Coefficient A.1 Channel Flow Model The simplest model for the microvalve system is posed as steady compressible flow of a perfect gas between two infinite parallel plates (Figure. 4.1). The gas emerging from the inlet port flows over the repeating obstacles that represent the actual seat rings in the microvalve. This represents a 2-dimensional flow problem that requires the solution of the governing partial differential equations. Instead of following such a complicated path, a 1-D channel flow model is proposed using a modified friction factor that accounts for the repeating rectangular obstacles (Figure 4.2). In essence, this constitutes substituting the friction coefficient, i.e. 96/Re Dh , with a friction factor that accounts for the height of the obstacles (E = e/H), the spacing between the obstacles (D = d/H), and the Reynolds number (Figure A.1). Using computational data of Luy et al. (1991), a correlation was generated to calculate the modified friction factor as a function of the Re Dh and geometrical parameters (D and E). The original computational data of Luy et al. (1991) has parameter ranges for Re Dh , E, and D of 10-300, 0.1-0.5, and 1-4, respectively. A sample of the CFD-based results of Luy et al. (1991) are presented in Figure A.2. Based on the trends of the data of Luy et al. (1991), the suggested correlation is of the following separated form: cD b aE eeKeF f f h D Re channel modified == , (A.1) with E = e/H and D = d/H being the dimensionless obstacle height and spacing, respectively. The coefficients K, a, b, and c are to be determined using the procedure outlined next. Naturally, the proposed ratio is always greater than unity due to the presence of the obstacles. Using a graph digitizing software (DigXY, 2003), a total of 68 points were gathered. By organizing the data into sets holding two of the three parameters (Re Dh , D, and E) constant, the Least-Squares Fit (LSF) method was used to determine each of the three coefficients a, b, and c. For example, for five data points with Re Dh = 10 and D = 2, the LSF method leads to the following relation: 2.610205E channel modified e f f ? . (A.2) The five data points along with the relation (A.2) are shown in Figure A.3. The exponential coefficient 2.610205 and all other constants for the sets of f vs. E data are categorized as ?a? constants. Similarly, the constants for f vs. Re Dh data are grouped as ?b? constants, and f vs. D data as ?c? constants. Using a weighted averaging scheme, the correlation?s exponential constants (a, b and c) were determined from the groups of constants from each set. These exponential constants are summarized in Table A.1 where the minimum, maximum, weighted average, and standard deviation for each constant are given. These average values were used exclusively from this point on to evaluate the exponential relations. 148 For each of the 68 data points, the value of f was divided by the value calculated from the exponential relations. For example, for the data point with E = 0.1, Re Dh = 10 and D = 2: .5471.0 1.03971 )2(13361.0)10(00117.0)1.0(866.3 Re 2,10Re,1.0 == = ? === eee eee F K cD b aE data DE h D h D (A.3) The universal constant K was determined by averaging all the K values for each point. The maximum, minimum, weighted average and standard deviation of the coefficient K are also listed in Table A.1. Table A.1 Summary of Coefficients for the Modified Friction Correlation a b c K Max 4.230523 0.003419 -0.04646 1.20007 Min 2.610205 0.00016 -0.26996 0.329265 W. Avg. 3.86602 0.00117 -0.13361 0.713826 St. Dev. 0.496064 0.000927 0.066154 0.227608 The final form of the correlated data of Luy et al. (1991) is: DE eeeF f f h D 13361.0 Re00117.0 86602.3 channel modified 713826.0 ? == . (A.4) When comparing the actual friction coefficient to that predicted by equation (A.4), an average error of 14.2 % and maximum error of 40.6% were determined. The comparison between the actual data and those calculated from the correlation (A.4) is shown in the Figure A.4. 149 A.2 Application of the Modified Friction Correlation to the Radial Flow Model In applying this correlation to the geometry of the radial flow model, the ranges of the governing parameters must be comparable. The Re Dh for the radial flow model varies with radial distance, since density times velocity is not constant due to area change. The Re Dh and E are contained in the ranges of 0.19-450 and 0.67-0.91, respectively, which are close to the ranges of the correlation?s parent data. However, the range of D for the model is 18-25 which is well beyond the range of the correlation of 0-4. With the limited data and range used to develop the D relation in the correlation, the trend can not be accurately captured for use with the radial flow model?s geometric parameters. Furthermore, the use of the correlation in its current form results in erroneous friction factors less than 96/Re Dh , which would give the impossible result of friction conditions less than that of Fanno flow. To correct for this problem, the knowledge that as ??D results in a friction factor of 96/Re Dh is utilized. Therefore, the D relation part of the correlation is set to 1 for use with the radial flow model geometry. Therefore, the following relation is used: h D eeF f f E Re00117.0 86602.3 channel modified 713826.0== . (A.5) 150 Figure A.1: Geometry for the study of Luy et al. (1991) 151 Figure A.2: Sample computational data from Luy et al. (1991) 152 E f mo d i f i e d /f ch a n n e l 0.1 0.2 0.3 0.4 0.5 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3 Luy et al. (1991) Least Squares Fit Re Dh =10,D=2 Figure A.3: LSF for the five-data-point set with Re Dh =10 and D=2 153 Re F/ ( e aE e cD ) 0 50 100 150 200 250 300 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 f/(e aE e cD ) Ke bRe Figure A.4: Correlated approximation in comparison with data of Luy et al. (1991) 154 Appendix B Radial Stokes Flow Between Two Concentric Parallel Disks 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 B.1, with the distance between the disks being 2h. It is assumed that for this axisymmetric flow, the velocity components in the z and ? directions are negligible ( ). The flow is laminar, incompressible and steady. 0? r V B.1 Variation of the Radial Velocity Component By applying these assumptions to the governing equations, the continuity and momentum equations in the r and z directions are simplified to: ,0)( = ? ? r rV r (B.1) , 2 2 z V r P r V V rr r ? ? + ? ? ?= ? ? ?? (B.2) ,0= ? ? z P (B.3) 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 ( ) that automatically satisfies the continuity equation (Eq. B.1). The new form of the momentum equation in the r direction will be: r rVzf =)( 155 3 2'' r f r f dr dP ?? += . (B.4) Note that pressure is just a function of r (Eq. B.3), so the partial differentiation of pressure in Eq. B.2 is changed to ordinary differentiation in Eq. B.4. Equation B.4 is not solvable analytically because of the nonlinear term ( 3 2 r f ). This term represents the convective effect in the Navier-Stokes equation. This term can be neglected by assuming very low Reynolds number (very low fluid velocity, i.e. Stokes flow, which is valid for Re << 1 (White, 1991)). Equation B.4 then becomes: r f dr dP '' ?= . (B.5a) Upon separating variables, each side depends on the independent variables r and z, leading to: constant '' == f dr dP r ? . (B.5b) Upon integrating twice, one gets: 21 2 2 )( CzCz dr dPr zf ++= ? . (B.6) The physical boundary conditions are: At : , (B.7a) hz ?= 0for 0 == r V At : 0=z 0for 0 =?= ? ? z V r . (B.7b) 156 Upon applying the boundary conditions, one gets: .1 2 )( 2 2 ? ? ? ? ? ? ? ? ?? ? ? ? ? ? == h z dr dPh r zf V r ? (B.8) Denoting the constant in Eq. (B.5b) by K, we have: , r K dr dP = (B.9) integration of which leads to: ,ln ? ? ? ? ? ? ? ? =?=? i o io r r KPPP (B.10) with subscripts i and o denoting the inlet and outlet positions, respectively. Then: , )ln( i o r r r P dr dP ? = (B.1) and finally: .1 )ln( 1 2 )( 2 2 ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? == h z r r P r h r zf V i o r ? (B.12) It is observed that under the Stokes flow conditions, at a given radial position, the radial velocity component varies parbolically. 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 dPrh dzrVQ ? ? ? ? ? ? === ? + ? . (B.13) 157 B.2 Friction Factor The friction factor for the radial Stokes flow between two concentric disks can now be determined. The wall shear stress, ? f , is computed from the wall velocity gradient: dr dP h dz dV hz r f == += ?? . (B.14) An expression for dr dP in terms of average velocity is needed. Dividing equation B.13 by cross-sectional area, A = 2?r(2h), one gets the average velocity: ? ? ? ? ? ? ? === o i r r r Ph dr dPh hr Q V ln 33)2(2 22 ??? . (B.15) Solving Eq. B.15 for dr dP , one gets: 2 3 h V dr dP ? = . (B.16) Substituting for dr dP in Eq. B.14, one gets: h V f ? ? 3 = . (B.17) Nondimensionalizing by the dynamic pressure, one gets the Darcy friction factor, f: Vh V f f ? ? ? ? 24 4 2 2 1 == . (B.18) Taking note of the hydraulic diameter, D h , 158 h r hr perimeter A D h 4 )2(2 )22(44 -sectionalcross === ? ? . (B.19) Finally, Equation B.18 becomes, h D f Re 96 = . (B.20) B.3 Total Pressure Drop Relation The total 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 outlet position, denoted ?o?and the inlet 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 ?+?= ?+?= +?+=? ? ? ?? (B.21) Substituting for the average velocities from equation B.15, one gets: ? ? ? ? ? ? ? ? ?+?=? 22 io 11 22 2 32 tt PPP io rr h Q ? ? . (B.2) Substituting for ?P from equation B.13 gives: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?+ ? ? ? ? ? ? =? 22io 11 824 1 tt ln 3 PP io o i rr r r Q h h Q ? ? ? ? . (B.23) 159 Figure B.1: Geometry of the radial flow between two parallel concentric disks 160 Appendix C Sample Calculation for the Initial Flow Conditions C.1 Extrapolation of the Experimental Data Figure 5.1 shows the limited experimental results of Yang et al. (2004) for both microvalve displacement (?) and volumetric flowrate (Q) vs. applied voltage (v). It is a common practice in industry to specify mass flowrate of a gaseous medium in terms of volumetric flowrate at standard conditions. The volumetric flowrate of gaseous flows can be expressed in Standard Cubic Centimeters per Minute (SCCM). The standard form is defined as flowrates expressed at standard reference conditions (20?C and 1 atm). The units of the volumetric flow rate must be converted from SCCM to actual flowrate (ACMM) before calculating the velocity of the flowing gas. Using digitization software (DigXY, 2003) and a Least-Squares Fit (LSF) method of curve-fitting, equations for the experimental data were extrapolated from Figure 5.1. The LSF equations are listed below. Microvalve displacement (?) in ?m vs. applied voltage (v) in volts: 02.0)1291.0( += ?? . (C.1) Volumetric flowrate in SCCM (Q) vs. v in volts for P ti = 100 psig: 8.024.02.0 2 ?+= ??Q . (C.2) Volumetric flowrate in SCCM (Q) vs. v in volts for P ti = 200 psig: 161 4486.0763.15983.0 2 +?= ??Q . (C.3) Volumetric flowrate (Q) in SCCM vs. v in volts for P ti = 300 psig: 7742.0494.2749.0 2 +?= ??Q . (C.4) The actual and extrapolated experimental data of Yang et al. (2004) are summarized in Figure C.1. C.2 Conversion of the Units for Flowrate For illustration purposes, the details involved in the conversion of the units for the volumetric flowrate are given for one case. Considering an initial total pressure of 200 psig and a deflection of 4 ?m, the volumetric flowrate for this case is determined. Using equation C.1, the voltage corresponding to a displacement of 4 ?m is calculated: volts83.30 1291.0 02.04 = ? =? . Substituting this voltage into equation C.3 to determine the standard flowrate, one gets: .SCCM73.5144486.0)83.30)(763.1()83.30)(5983.0( 2 standard =+?=Q This standard volumetric flowrate can then be converted to an actual value in cubic meters per minute (ACMM) via the following relation: .10783.3 10 1 kPa 1378.95 kPa 101.35 293 K 293 )73.514( cm10 m1 min. m5 6 36 3 actual STP STP actual standardactual 3 ? ?= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = K P P T T QQ (C.5) 162 C.3 Calculation of the Initial Flow Properties Continuing with the current example, the initial velocity can be calculated and used to determine the initial Mach number and Reynolds number. Dividing the volumetric flowrate by cross-sectional area, A = 2?r i (H), one gets, s m m i s mm V 39.2 60 .min1 )1014)(003.0(2 10783.3 6 min 5 3 = ? ? ? ? ? ? ? ? ? ? = ? ? ? . The gas being considered is helium. Using the initial velocity to compute the initial Mach number, M i , one gets, 00237.0 293)2077(67.1 39.2 === s m i i RT V M ? . Substituting the initial velocity into the Reynolds number equation gives: ( )( ) 71.7 10953.1 102839.225.2 Re 5 6 3 = ? ? ? ? ? ? ? ? == ? ? ? s mkg s m m kg hi iD m DV h ? ? . The properties of helium were taken from Tsederburg et al. (1971). For the current example, the dimensionless ring height, E = e/H, and the Reynolds number can be used to determine the initial modified friction coefficient. Substituting the initial parameters for the current example into equation A.5 (see Appendix A), the initial modified friction coefficient (F i ) is calculated: .397.11713826.0 713826.0 )71.7(00117.00.714286)(86602.3 Re00117.0 86602.3 channel modified == == ee ee f f F h DE i 163 applied voltage (v) [volts] Fl o w r a t e ( Q ) [ S CCM ] Di s p l a c e m e n t ( ? )[ ? m] 0 1020304 0 100 200 300 400 500 600 700 800 900 1000 1 2 3 4 5 Yang et al. (2004) (P ti = 100 psig) Yang et al. (2004) (P ti =200psig) Yang et al. (2004) (P ti =300psig) Displacement vs. Applied Voltage Simulation Cases Curve-fit (P ti =100psig) Curve-fit (P ti =200psig) Curve-fit (P ti =300psig) 0 Figure C.1: The actual and extrapolated experimental data of Yang et al. (2004) 164 Appendix D Reynolds Number Relation in a Radial Disk A relation for the Reynolds number (Re Dh ) dependence on the radial distance, r + , is determined. This relation can be used to verify the predictions of the computer code for the Reynolds number. The Reynolds number at a radial position is defined as: ? ? h D VD h =Re , (D.1) where ? is density, V is average radial velocity, D h is hydraulic diameter, and ? is the fluid viscosity. Consider the expression for mass flowrate, i.e: VAm ?= & . (D.2) Substituting equation D.2 into equation D.1, one gets: ?A Dm h D h & =Re . (D.3) Substituting for cross-sectional area (A) and hydraulic diameter (D h ) from equation 4.37c gives: ??r m h D & =Re , (D.4) 165 showing that the Reynolds number varies inversely with the radial distance. Note that for the initial Re Dh at the inlet port (Re i ): ?? i i r m & = Re . (D.5) Dividing equation D.4 by D.5, one gets: + + === r r r i D i D h h 1 Re Re Re . (D.6) Equation D.6 shows that Re + Dh is inversely proportional to dimensionless radial distance, r + . A graph of the relation is shown in Figure D.1. 166 r + Re Dh / Re i 00.250.50.751 1 4 7 10 13 16 19 22 25 28 Figure D.1: Variation of Re + Dh with dimensionless radial coordinate, r + 167 Appendix E Derivation of the Analytical Expression for L max E.1 Analytical L max for the Channel Flow Model Using a procedure described by Saad (1993), an analytical expression for L max for a channel (constant f) can be derived through integration of Equation 4.14. Setting the limits of integration to M at the entrance of the channel and M = 1 at L max , the integral of Equation 4.14 becomes: ( ) dM MM M D fdx M L h ?? ? ? ? ? ? ? ? + ? = 1 32 2 0 2 1 1 12 max ? ? . (E.1) Noting that dM 2 = 2MdM, Equation E.1 becomes: 2 1 42 2 max 2 2 1 1 1 dM MM M D fL M h ? ? ? ? ? ? ? ? + ? = ? ? . (E.2) Before integrating, the method of partial fractions is used to put the right hand side of Equation E.2 into an easier form for integration. Letting x = M 2 , and denoting 2 )1( ? = ? a the integrand of Equation E.2 becomes: ax C x B x A axx x MM M + ++= + ? = ? ? ? ? ? ? + ? ? 1 )1( 1 1 1 22 42 2 1 2 ? . (E.3) 168 Multiplying both sides by x 2 (1 + ax), one gets: 2 )1()1(1 CxaxBaxAxx ++++=? . (E.4) Combining terms will lead to: BxaBAxCaAx ++++=? )()(1 2 . (E.5) Equating terms having the same power of x gives: A = - (1 + a), B = 1, C = a(1 + a). Thus, Equation E.2 becomes: 2 1 242 1 max 2 1 )1(11 dM Ma aa MM a D fL M h ? ? ? ? ? ? ? ? ? + + ++ ?? = ? . (E.6) Integrating Equation E.6, one gets: 1 2 2 2 1max 2 )1ln()1( 1 )ln()1( M h Maa M Ma D fL ? ? ? ? ? ? +++?+?= ? . (E.7) Applying the limits and combining like terms gives: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + + +++?= 2 2 2 1max 1 )1( ln)1( 1 1 Ma Ma a M D fL h ? . (E.8) Subsituting ?a? into the equation, one gets: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ++ + ? = 2 2 2 2 max 2 1 12 )1( ln 2 11 M M M M D fL h ? ? ? ? ? . (E.9) 169 Recalling that for channel flow, f = 96/Re Dh , the dimensionless L max is given by: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ++ + ? == + 2 2 2 2 max max 2 1 12 )1( ln 2 11 96 Re M M M M D L L h D h ? ? ? ? ? . (E.10) 170 Appendix F Derivation of Static Pressure Relations for the Radial Disk Geometry F.1 Ideal Incompressible Flow Relation for Static Pressure The derivation for this relation is started with the Bernoulli?s equation. The Bernoulli?s equation is: 2 2 12 2 1 ii VPVP ?? +=+ . (F.1) Rearranging Equation F.1 for P, one gets: ( ) 22 2 1 VVPP ii ?+= ? . (F.2) An expression for V can be found using the continuity equation as follows: ii VhrVhr )2(2)2(2 ?? = . (F.3) Rearranging Equation F.3 for V, one gets: r rV V ii = . (F.4) Substituting V from Equation F.4 into Equation F.2 gives: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?+= 2 2 2 1 1 r r VPP i ii ? . (F.5) Equation F.5 is nondimensionalized to give: 171 ( ) ? ? ? ? ? ? ? ? ?+= ? ++ 2 2 2 1 11 r P V P i i ? . (F.6) Substituting for i P ? from the perfect gas relation, one gets: ( ) ? ? ? ? ? ? ? ? ?+= ? ++ 2 2 2 1 11 r RT V P i i . (F.7) Introducing the Mach number (Eq. 4.6a) into equation F.7 gives: ( ) ? ? ? ? ? ? ? ? ?+= ? ++ 2 2 2 1 11 rMP i ? . (F.8) F.2 Ideal Compressible Flow Relation for Static Pressure The derivation for the ideal compressible flow case is started with the difference form of the Bernoulli?s equation. Neglecting gravity effects, the difference form of the Bernoulli?s equation (White, 1999) is: 0=+VdV dP ? . (F.9) Substituting for ? from ideal gas relation assuming constant temperature, and integrating, one gets: constantln 2 2 1 =+ VPRT . (F.10) Setting the right hand side to conditions at the inlet port gives: 172 2 2 12 2 1 lnln ii VPRTVPRT +=+ . (F.1) An expression for V can be found from the continuity equation: VhrVhr iii ???? )2(2)2(2 = , (F.12) which upon introducing the perfect gas relation gives, P P r rV V iii = . (F.13) Rearranging equation F.11 and substituting for V from equation F.13, one gets: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= 2 2 1 2 ln P P r r RT V P P iii i . (F.14) Introducing the Mach number (Eq. 4.6a) into equation F.14 gives: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= 2 2 1 2 ln P P r r M P P ii i i ? . (F.15) Equation F.15 can be nondimensionalized to give: ( ) ? ? ? ? ? ? ? ? ?= ? +++ 2 2 1 2 ln PrMP i ? . (F.16) This equation can not be solved explicitly for P + , however it can be evaluated numerically for values of ?, M i and r + < 1. For the present simulated cases for which M i ?s are extremely small, it can be shown that P + ? 1. 173 F.3 Stokes Flow Relation for Static Pressure (Incompressible Fluid) The Stokes flow relation for the variation of the static pressure is derived by integrating Equation B.11: ?? ? = r r i o io P P ii r dr r r PP dP )ln( )( . (F.17) Equation F.17 becomes: ? ? ? ? ? ? ? ?? += i i o io i r r r r PP PP ln )ln( )( . (F.18) Equation B.13 yields an expression for )ln( )( i o io r r PP ? , that is: 3 4 3 )ln( h Q r r P i o ? ? ?= ? . (F.19) Substituting Equation F.19 into Equation F.18, one gets: ? ? ? ? ? ? ? ? ?= i i r r h Q PP ln 4 3 3 ? ? , which finally becomes the dimensionless expression for static pressure: ++ ?= r Ph Q P i ln 4 3 1 3 ? ? . (F.20) 174