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