Experimental and Analytical Investigation of a Dynamic Gas Squeeze Film
Bearing including Asperity Contact Effects
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.
Manoj Deepak Mahajan
Certificate of Approval:
George T. Flowers
Professor
Mechanical Engineering
Robert L. Jackson, Chair
Assistant Professor
Mechanical Engineering
Jay M. Khodadadi
Professor
Mechanical Engineering
Joe F. Pittman
Interim Dean
Graduate School
Experimental and Analytical Investigation of a Dynamic Gas Squeeze Film
Bearing including Asperity Contact Effects
Manoj Deepak Mahajan
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 15, 2006
Experimental and Analytical Investigation of a Dynamic Gas Squeeze Film
Bearing including Asperity Contact Effects
Manoj Deepak Mahajan
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
iii
Vita
Manoj Mahajan, son of Deepak and Neeta Mahajan was born on January 19, 1981,
in Nasik District in India. He attended Government College of Engineering, Pune and
graduated in November 2002 with the degree of Bachelor of Engineering in Mechanical
Engineering. He joined the Masters program in the department of Mechanical Engineering
at Auburn University in August 2004.
iv
Thesis Abstract
Experimental and Analytical Investigation of a Dynamic Gas Squeeze Film
Bearing including Asperity Contact Effects
Manoj Deepak Mahajan
Master of Science, December 15, 2006
(B.E., Government College of Engineering, Pune University, 2002)
123 Typed Pages
Directed by Robert L. Jackson
This thesis presents a theoretical and an experimental investigation of planar gas
squeeze film bearings. The thickness and pressure profile of the gas squeeze film are ob-
tained by simultaneously solving the Reynolds equation and the equation of motion for the
squeeze film bearing. This work also accounts for the force due to surface asperity contact
in the equation of motion. When the surfaces are in contact, the model predicts the contact
force as a function of film thickness. Computational simulations are performed to study
the development of the squeeze film from its initial state to a pseudo-steady state condition
and to evaluate its load carrying capacity. For certain cases, the simulation results correlate
well with the pre-established analytical results. However, corrections must be made to the
analytical equations when they are used out of their effective range. In the experimental
study, a squeeze film is developed due to an applied relative normal motion between two
v
parallel circular plates of which one circular plate is effectively levitated. Theoretical results
for the squeeze film thickness match qualitatively with its experimental counterpart.
On successful testing of macro-scale gas squeeze film bearings, micro-scale bearing
surfaces are fabricated. Experimental investigation of micro-scale bearings suggests that
these bearings have significant potential for a wide range of applications in Micro-Electro
Mechanical Systems (MEMS).
vi
Acknowledgments
I wish to thank my advisor, Dr. Robert L. Jackson, for providing me with an oppor-
tunity to conduct research in the exciting field of tribology. His guidance, support and
encouragement have helped me towards the successful completion of my Master?s Degree
in Engineering. He is a great mentor and working with him was a wonderful experience.
I would like to extend my appreciation and thanks to Dr George Flowers for all the help
with the test setup and being part of my graduate committee. I would also like to thank
Dr. Jay Khodadadi for serving as a graduate committee member.
Many thanks to Klaus Hornig, Dr. Roland Horvath and Alfonso Moreira for their help
and motivation at the Vibration Analysis Laboratory. Special thanks to Mr. Charles Ellis
and Abhishek for all the help at Alabama Microelectronics Science and Technology Center.
I wish to acknowledge my companions here at Auburn, Harish, Rajendra, Harshavardhan,
Ananth and Ravi Shankar for their friendship.
I wish to dedicate this work to my parents and sisters for their enduring love, immense
moral support and encouragement in the journey of life.
vii
Style manual or journal used LATEX: A Document Preparation System by Leslie
Lamport (together with the style known as ?aums?) and Bibliography as per Tribology
Transactions
Computer software used TEX (specifically LATEX), MATLAB 7.0.4, MS Office
PowerPoint 2003 and the departmental style-file aums.sty
viii
Table of Contents
List of Figures xi
List of Tables xiii
Nomenclature xiv
1 Introduction 1
2 Background 3
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Lubrication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.3 Squeeze Film Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Incompressible Squeeze Film Bearings . . . . . . . . . . . . . . . . . . . . . 6
2.5 Compressible Squeeze Film Bearings . . . . . . . . . . . . . . . . . . . . . . 7
2.6 Squeeze Film Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.7 Micro-Scale Squeeze Film Bearings and Effects Due to Molecular Dynamics 18
3 Objectives 24
4 Numerical Investigation 26
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2 Thin Squeeze Film Bearings . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2.1 Formulation of Coupled Dynamics . . . . . . . . . . . . . . . . . . . 26
4.2.2 Computational Scheme . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.3 Ultra-Thin Squeeze Film Bearings . . . . . . . . . . . . . . . . . . . . . . . 47
4.3.1 Formulation of Coupled Dynamics . . . . . . . . . . . . . . . . . . . 47
4.3.2 Solution Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5 Experimental Investigation 50
5.1 Thin Squeeze Film Bearings . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.1.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.1.2 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2 Ultra-thin Squeeze Film Bearings . . . . . . . . . . . . . . . . . . . . . . . . 63
5.2.1 Design of Micro-Scale Bearing Surfaces . . . . . . . . . . . . . . . . 63
5.2.2 Fabrication of Micro-Scale Bearing Surfaces . . . . . . . . . . . . . . 65
5.2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
ix
6 Summary 73
Bibliography 75
Appendices 79
A Contact Force 80
B C-program to solve the coupled dynamics for thin squeeze films 86
C Tables of Experimental and Simulation Results 93
C.1 Configuration 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
C.2 Configuration 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
C.3 Configuration 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
D Calibration of Capacitance Sensor 98
E Computer program to solve the dynamics for ultra-thin squeeze films102
x
List of Figures
2.1 Classification of lubrication . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 The squeeze film thickness, h, as a function of normalized time, T,
(h=hm(1+epsilon1?cos(?t) and epsilon1=?h/hm) . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Analogy of a squeeze film as a spring-damper system . . . . . . . . . . . . . 14
4.1 Schematic of a planar squeeze film bearing . . . . . . . . . . . . . . . . . . . 27
4.2 Scheme for discretization of i-spatial variable and j-time variable . . . . . . 28
4.3 Free body diagram of a squeeze film bearing (one degree of freedom) . . . . 30
4.4 Algorithm for computational simulation . . . . . . . . . . . . . . . . . . . . 34
4.5 Dimensionless P as a function of normalized T at R=0.5, ?=1000 and H=1-
0.5?sin(T) (16 nodes) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.6 Dimensionless P as a function of normalized T at R=0.5, ?=1000 and H=1-
0.5?sin(T) (160 nodes) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.7 Variation in film height for different time steps . . . . . . . . . . . . . . . . 37
4.8 Dynamic behavior of the squeeze film height for given input conditions as a
function of time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.9 Comparison of the dimensionless squeeze film force, F(T)/(pi?R02?patm), as
a function of normalized time, T, from the numerical simulation and from
Langlois squeeze film force model [8] . . . . . . . . . . . . . . . . . . . . . . 40
4.10 Change in percentage error between Fi from the numerical simulation and
Fn as predicted by Salbu?s Eq. as a function of the squeeze number, ? . . . 43
4.11 The results of numerical parametric study for hm as a function of Fi at a
constant frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.12 The results of numerical parametric study for hm as a function of ? . . . . 46
xi
4.13 Algorithm for computational simulation of the ultra-thin squeeze film bearing 49
5.1 Schematic explaining the experimental setup . . . . . . . . . . . . . . . . . 50
5.2 Test stand used for experimental investigation . . . . . . . . . . . . . . . . . 51
5.3 Disks used for levitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.4 Photograph illustrates the experimental setup with the capacitance sensor . 53
5.5 Electrical circuit for capacitance measurement . . . . . . . . . . . . . . . . . 54
5.6 Photograph illustrates the experimental setup with the laser displacement
measurement system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.7 Plot of DC voltage from the laser beam displacement measurement system 56
5.8 Experimental and simulation results for the squeeze film height against the
amplitude of vibration for bearing configuration 1 (Capacitance measurement
system) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.9 Experimental and simulation results for the squeeze film height against the
amplitude of vibration for bearing configuration 2 (Laser displacement mea-
surement system) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.10 Experimental and simulation results for the squeeze film height against the
amplitude of vibration for bearing configuration 3 (Laser displacement mea-
surement system) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.11 LASICAD drawing of micro bearing surfaces (not to scale and only portions
of the textured sections are shown in this enlarged view) . . . . . . . . . . . 64
5.12 Fabrication procedure for micro-scale bearing surfaces . . . . . . . . . . . . 65
5.13 Schematics of single unit cell (not to scale) . . . . . . . . . . . . . . . . . . 68
5.14 Experimental results of micro-scale bearing . . . . . . . . . . . . . . . . . . 70
5.15 Dimensionless mean pressure, Pn, as a function of normalized time, T, for
the ultra-thin squeeze film bearing . . . . . . . . . . . . . . . . . . . . . . . 71
A.1 Surface profile of a rough bearing surface . . . . . . . . . . . . . . . . . . . 83
A.2 Dimensionless contact load as a function of dimensionless mean surface sep-
aration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
xii
List of Tables
2.1 Assumptions and parameters used to estimate the mean free path . . . . . 20
2.2 Types of flows as per Knudsen number [34] . . . . . . . . . . . . . . . . . . 20
4.1 Comparison of load carrying capacity between Fi in numerical simulation
and Fn by Salbu?s Eq. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.1 Bearing configurations used for the experimental purpose . . . . . . . . . . 51
5.2 Silicon Wafer Cleaning Procedure [42] . . . . . . . . . . . . . . . . . . . . . 66
xiii
Nomenclature
English Symbols
A area of contact (m2)
B bearing breadth (m)
bei, ber Kelvin functions of order zero
bei1, ber1 Kelvin functions of order one
C capacitance (F)
Cairgap capacitance between the squeeze film bearing surfaces (F)
d gap due to air (m)
E modulus of elasticity (Pa)
F squeeze film force (N)
Fi applied load, i.e. weight of levitating disk (N)
Fd squeeze film damping force (N)
Fn mean positive film force (N)
f frequency (Hz)
G gas constant (J?K?1?mole?1)
H dimensionless film thickness, h/h0
H? dimensionless film thickness, h/hm
h thickness of squeeze film in z direction (m)
xiv
h0 initial film thickness (m)
hm mean film thickness (m)
Kn Knudsen number
k dielectric constant for air
m mass of the levitated plate (kg)
P dimensionless pressure, p/patm
p pressure (Pa)
patm atmospheric pressure (Pa)
QP rarefaction coefficient
R dimensionless radius, r/R0
R0 radius of area of contact (m)
Rq RMS surface roughness (m)
r radial coordinate in cylindrical polar coordinates
Sy yield strength (Pa)
T normalized time (rad), ?t
Tg gas temperature (K)
t time (s)
UR, U?R physical components of surface motion in cylindrical polar coordinates
Ux, U?x dimensionless surface velocity components, ux/V, u?x/V
Uy, U?y dimensionless surface velocity components, uy/V, u?y/V
ux, u?x surface velocity components of respective bearing surfaces
xv
in x-direction (m?s?1)
uy, u?y surface velocity components of respective bearing surfaces
in y-direction (m?s?1)
V reference velocity (m?s?1)
V in input voltage to the capacitance sensor (V)
V o/p output voltage measured across the capacitor (V)
Wn mean load capacity
X, Y dimensionless right-handed Cartesian coordinates, x/B, y/B
x, y, z right-handed Cartesian coordinates (m)
Z0 amplitude of vibration of base plate (m)
zb displacement of base plate (m)
zt displacement of top disk (m)
Greek Symbols
?R step-size for R
?T step-size for T (rad)
?h amplitude of oscillation of squeeze film at a steady state
epsilon1 excursion ratio
epsilon10 permittivity of free space (F?m?1)
? polytropic coefficient
? angular coordinate in cylindrical polar coordinates (rad)
xvi
? bearing number
? mean free path (m)
? viscosity (Pa?s)
?eff effective viscosity (Pa?s)
? Poisson?s Ratio
? density (kg?m?3)
? squeeze number based on initial film height
?? squeeze number based on mean film height
?XP, ?Y P pressure flow factors in x and y direction
? angular velocity (rad?s?1), 2pif
Subscripts:
cont contact between rough surfaces
i nodal value in radial direction
j nodal value in time
n nominal or apparent value
xvii
Chapter 1
Introduction
Research in the field of tribology, the science and technology of friction, wear and
lubrication, has been continuously improving the performance of mechanical systems ever
since the industrial revolution. A British government report (the ?Jost Report?) in 1966
estimated a potential savings of ? 515 million per annum only for the United Kingdom
by better application of tribological principles and practices [1]. Application of principles
of tribology led to the development of many different types of bearings used for different
purposes. Thrust bearings, hydrodynamic bearings, hydrostatic bearings, rolling element
bearings, etc. are widely used as means to reduce the frictional losses in mechanical systems
based on the selection criteria. The selection criteria include speed, load, life, maintenance,
space requirements and environmental conditions etc.
The largely underutilized squeeze film effect has a potential for use as a means to
lubricate surfaces by creating squeeze film bearings. Applications of such squeeze film
bearings can be used in read/write heads in hard disk drives, manufacturing processes,
vibrating machinery, low-speed applications and hydrodynamic bearing start up and shut
down. In this research, the thin gas squeeze film bearings are studied extensively. Both
numerical and experimental results confirm that squeeze film bearings can be used as means
for lubrication. Micro-scale bearing surfaces are fabricated to study ultra-thin squeeze films
for lubrication of miniaturized mechanical systems.
1
In Chapter 2 of this thesis, a review of the squeeze film effect with the perspective
of lubrication and damping is presented. The various lubrication regimes are described.
Hydrodynamic lubrication as well as squeeze film lubrication are defined. Physics and
fundamentals of the squeeze film effect are explained. Incompressible and compressible
squeeze film bearings, with their governing equations, are reviewed. The latter part of
chapter 2 gives a background of squeeze film damping. It gives a concise overview of the
cut-off frequency, squeeze film damping and spring forces. This is followed by a review of
ultra-thin squeeze film bearings where gas rarefaction effects are significant.
In Chapter 3, all the research objectives as well as the specific goals pertaining to these
objectives are stated.
Chapter 4 is divided into two sections, namely Thin and Ultra-thin squeeze film bear-
ings. Each section covers subsections such as formulation of coupled dynamics (the equation
of motion and the Reynolds equation for the bearing), numerical scheme to solve the dy-
namics, followed by the results of the numerical investigation.
Chapter 5 explains all the observations made during the measurements of thin and
ultra-thin squeeze film bearings. Here, experimental results are compared with the simula-
tion results.
Chapter 6 summarizes the thesis. In summary, a thorough study was conducted on the
squeeze film bearings for their use as a potential means to lubricate surfaces.
2
Chapter 2
Background
2.1 Introduction
This chapter discusses in detail the different types of lubrication regimes. The phe-
nomenon of squeeze film effect with the perspective of lubrication is then described. In-
compressible and compressible squeeze film bearings with their governing equations are
reviewed. Then, the squeeze film damping phenomenon and its cut-off frequency are ex-
plained. Lastly, ultra-thin (nano-scale) squeeze film bearings including effects due to molec-
ular dynamics are studied.
2.2 Lubrication
?Lubrication is an application of a lubricant between two surfaces in relative motion for
the purpose of reducing friction and wear or other forms of surface deterioration? [2]. The
lubricant is usually a fluid, but in some cases it can be a solid such as a powder. Lubrication
is broadly classified into fluid-film, boundary and mixed regimes (see Fig. 2.1). In fluid-
film lubrication, bearing surfaces are completely separated by either a liquid or a gaseous
lubricating film [3]. If the loads are high or speeds are low then the contact between high
or tall asperities is likely to occur. This is a boundary lubrication regime where a suitable
molecular layer of lubricant covers the high asperities. Hence, metal welding due to adhesion
3
is avoided [3]. The lubrication regime between boundary and fluid-film is categorized as
mixed or partial lubrication. In the mixed lubrication regime, effects due to both, boundary
and fluid-film lubrication are observed [3]. In fluid-film lubrication, a thin fluid film between
Lubrication
Fluid-Film Mixed Boundary
Hydrodynamic Hydrostatic Elastohydrodynamic
Sliding
Squeeze film
Figure 2.1: Classification of lubrication
bearing surfaces is obtained by either hydrostatic or hydrodynamic action. Hydrostatic
lubrication is a phenomenon of maintaining a lubricating film by external means; whereas,
hydrodynamic lubrication is self-acting. In hydrodynamic lubrication, positive film pressure
between conformal surfaces is developed due to relative motion and fluid viscosity [3]. The
topic of interest here is squeeze film lubrication, which is a type of hydrodynamic lubrication
4
where a lubricating film is developed due to relative normal motion and fluid viscosity.
Elastohydrodynamic lubrication is a form of hydrodynamic lubrication where lubricating
surfaces are elastically deformable [3].
2.3 Squeeze Film Effect
The term ?squeeze film? defines a fluid film contained between two conformal, moving
surfaces with velocities of the surfaces normal to the planes of the containment [4]. If the
bearing surfaces approach each other then the motion is termed as ?positive squeeze?. Con-
versely, if the bearing surfaces move apart then the motion is termed as ?negative squeeze?
[4]. A relative normal motion between two parallel surfaces can produce a squeeze film
which can completely separate the surfaces and contribute to lubrication. This phenomenon
is known as the ?squeeze film effect? [3].
The load-carrying capacity results from the fact that a viscous flow cannot be squeezed
out of the gap without any delay; therefore, providing a cushioning effect and the film
equilibrium is established through a balance between viscous flow forces and compressibility
effects [5]. Thus, the flow of fluid at the boundary is approximately reduced to zero due
to high viscous forces resulting from alternate compression and decompression of the fluid
[6]. Alternate compression and decompression produces a steady-state film pressure which
oscillates about its mean value [6]. The average steady-state film pressure is greater than
the atmospheric pressure over one cycle (T to T+2pi) and thus provides the squeeze film lift
and lubrication [6]. Applications of the squeeze film effect are clutch packs in automotive
5
transmission, engine piston pin bearings, human hip and knee joints, damper films for jet
engine ball bearings and piston rings [7].
The dynamics of squeeze film is coupled as it includes both the equation of motion for
squeeze film bearing and the Reynolds equation. Here, the Reynolds equation governs the
generated fluid pressure. The general form of the Reynolds equation is given by [3]
?
?x
parenleftBigg
?h3
12?
?p
?x
parenrightBigg
+ ??y
parenleftBigg
?h3
12?
?p
?y
parenrightBigg
= ??x
parenleftbigg?h(u
x+u?x)
2
parenrightbigg
+ ??y
?
??h
parenleftBig
uy+u?y
parenrightBig
2
?
?+?(?h)?t (2.1)
Here, h is the squeeze film height and ? is the viscosity of the fluid. The time dependent
term ?(?h)?t is known as the ?squeeze term? as it represents squeezing motion of the fluid.
Here, ux, uy and u?x, u?y are surface velocity components of bottom and top surfaces in x
and y direction, respectively.
2.4 Incompressible Squeeze Film Bearings
In incompressible squeeze film bearings, the lubricant is a liquid and its density is
assumed to be constant in the operating range. As per Hamrock [3], when two surfaces
approach each other, it takes a finite amount of time to squeeze out the fluid, and this
action provides a lubricating effect. It is also interesting to note that it takes an infinite
amount of time to theoretically squeeze out all the fluid. For an incompressible fluid,
viscosity along with density is assumed to be constant in the Reynolds equation [3]. Also, if
6
only normal motion is considered and sliding velocities are zero then the Reynolds equation
in rectangular coordinates is given as [3]
?
?x
parenleftbigg
h3?p?x
parenrightbigg
+ ??y
parenleftbigg
h3?p?y
parenrightbigg
=12??h?t (2.2)
Eq. (2.2) in cylindrical-polar coordinates is expressed as
?
?r
parenleftbigg
rh3?p?r
parenrightbigg
+1r ???
parenleftbigg
h3?p??
parenrightbigg
=12?r?h?t (2.3)
The Reynolds equation given in the form of Eqs. (2.2) and (2.3) can be analytically
solved for pressure, film thickness and finite squeeze time (i.e. the amount of time for the
film to be squeezed out). Hamrock [3] provides analytical expressions for pressure, film
thickness and finite squeeze time for various geometries such as parallel-surface bearings
with infinite width, journal bearings with no rotation, a parallel circular plate approaching
a plane surface and a long cylinder near plane. As per Hamrock [3], the parallel film shape
produces the largest normal load-carrying capacity.
2.5 Compressible Squeeze Film Bearings
In contrast, for compressible squeeze film bearings, the lubricant used is gaseous (such
as ambient air). Langlois [8] was one of the first to extensively study the isothermal gas
7
squeeze films with the assumption of a thin and continuous gas film with a constant viscos-
ity. Langlois [8] also assumed that the density of the gas is proportional to the pressure,
which means that the gas squeeze film obeys the ideal gas law under isothermal condi-
tions. Langlois [8] derived the equation that governs the pressure variation for a thin, flat
isothermal gas squeeze film
?
?X
parenleftbigg
H3P ?P?X
parenrightbigg
+ ??Y
parenleftbigg
H3P ?P?Y
parenrightbigg
= ?
parenleftbigg ?
?X
bracketleftbigPH parenleftbigU
x +U?x
parenrightbigbracketrightbigparenrightbigg
+?
parenleftbigg ?
?Y
bracketleftBig
PH
parenleftBig
Uy +U?y
parenrightBigbracketrightBigparenrightbigg
+??PH?T (2.4)
Here, ? is the bearing number and ? is the squeeze number.
The definitions of ? and ? are
? = 6?BV/patmh02 (2.5)
? = 12?B2?/patmh02 (2.6)
Here, B is bearing breadth. The reference velocity, V, is used to obtain normalized Ux, U?x,
Uy and U?y (see Nomenclature), where the order of magnitude of these normalized velocities
is unity [8]. Likewise, Eq. (2.4) in cylindrical-polar coordinates is given as
8
?
?R
parenleftbigg
RH3P ?P?R
parenrightbigg
+ 1R ???
parenleftbigg
H3P ?P??
parenrightbigg
= ?
braceleftbigg ?
?R
bracketleftbigRPH parenleftbigU
R +U?R
parenrightbigbracketrightbigbracerightbigg
+?
braceleftbigg ?
??
bracketleftbigPH parenleftbigU
? +U??
parenrightbigbracketrightbigbracerightbigg+R??(PH)
?T (2.7)
Later, Langlois [8] dealt with the exact solution to squeeze film equations (Eqs. (2.4) and
0 2 4 6 8 10 12 14 16
8
8.5
9
9.5
10
10.5
11
11.5
12
T  (rad)
hh
hm
?h
h
h=hm(1+??cos(?t))
Figure 2.2: The squeeze film thickness, h, as a function of normalized time, T,
(h=hm(1+epsilon1?cos(?t) and epsilon1=?h/hm)
9
(2.7)) by assuming the squeeze film thickness as a function of time (see Fig. 2.2)
h = hm (1+epsilon1?cos?t) (2.8)
Based on the assumption of a film thickness given by Eq. (2.8), Langlois [8] introduced
a perturbation parameter for pressure of the order of epsilon1 in the Reynolds equation. Then,
Langlois [8] solved Eq.(2.4) to obtain the squeeze film force (i.e. the force due to the
difference between the pressure of the squeeze film and the ambient pressure for a given
instant of time). This solution is given for the squeeze film between two flat long parallel
plates and two parallel disks. Since the present work considers axisymmetric parallel disks
without lateral surface motions, the governing Reynolds equation as per [8] is
?
?R
parenleftbigg
RH3P ?P?R
parenrightbigg
= R??(PH)?T (2.9)
where ? is as given
? = 12?R02?/patmhm2 (2.10)
As per [8], if the squeeze number is very large then the gas between bearing surfaces
does not leak, and the squeeze film can be considered as incompressible. At low squeeze
10
numbers, the frequency of squeeze motion is low, and therefore gas leaks out [8]. The
squeeze film force between parallel disks given by [8] is
F (T) =
parenleftBig
pipatmR02epsilon1
parenrightBig
[?g1 (?)cosT+g2 (?)sinT] (2.11)
where,
g1 (?)=1?
radicalbigg2
? ?
ber??(bei1???ber1??)?bei??(ber1??+bei1??)
(ber??)2 +(bei??)2 (2.12)
g2 (?)=
radicalbigg2
? ?
ber??(ber1??+bei1??)+bei??(bei1???ber1??)
(ber??)2 +(bei??)2 (2.13)
Eq.(2.11) gives a good estimate of the squeeze film force between parallel disks at a
particular instant of time.
For compressible squeeze film bearings, Salbu [9] showed that the squeeze film effect
can be used to operate bearings in a highly vibrational environment, and also provided
analytical equations for bearing load carrying capacity. Salbu [9] assumes that the squeeze
film thickness is a known sinusoidal function of time. One disk is held stationary and
the other oscillates sinusoidally about a mean film thickness in a direction normal to the
surfaces with an amplitude of oscillation, ?h, and frequency, ?. Salbu [9] uses a simplified
11
form of the Reynolds equation (Eq.(2.9)). Then, the boundary conditions for the planar
radial squeeze film bearing are assumed as
P (R,T = 0)=1 (2.14)
i.e. initial pressure between disks is atmospheric. At the outer periphery,
P (R = 1,T)=1 (2.15)
i.e. the pressure is atmospheric at all times. At equilibrium, the mean positive film force,
is equal to the applied load. When ??? , ?the mass content rule? can be imposed on
Eq. (2.9) and the following equation can be derived to approximate the mean load carrying
capacity [9]
Wn= Fnpip
atmR02
=
??
?
bracketleftBigg1+ 3
2epsilon1
2
1?epsilon12
bracketrightBigg1/2
?1
??
? (2.16)
Salbu [9] numerically modeled Wn as a function of ? for various epsilon1 values and compared
it with the asymptotic values of Wn, as given by Eq. (2.16). For ??10, Salbu [9] observed
little variation between the numerical and analytical Wn predicted from Eq. (2.16) at the
same epsilon1. Thus, Eq. (2.16) may be used to predict Wn for ??10. For ??10, Eq.(2.16) is
12
used to validate the computational simulation results of the current analysis presented in
this thesis and will be referred to as Salbu?s Equation.
Further work on dynamic gas squeeze film bearing theory was performed by Minikes et
al. [10]. In [10], instead of assuming a squeeze film height as a function of time, it is deter-
mined by equilibrium of forces. Minikes et al. [10] used a piezoelectrically vibrating disk to
levitate another disk for experimental and numerical work. A sinusoidal voltage was applied
on the electrodes which gave rise to harmonic deformation of the piezoelectric disk at the
excitation frequency. The piezoelectric disk was also statically and dynamically deforming.
However, in the present work both disks are considered to be rigid. A similar approach of
coupled dynamics can also be found in [11], where the fluid film lubrication forces and seal
dynamics were solved simultaneously for noncontacting gas face seals, although the faces
were annular and not circular.
2.6 Squeeze Film Damping
?Damping is the element, present in all real systems, which dissipates vibrational en-
ergy, usually as heat, and so attenuates the motion? [12]. A squeeze film can also be modeled
as a spring-damper system (Fig. 2.3).
An early study of squeeze film by Griffin et al. [13] suggests that squeeze film between
two parallel plates provides viscous damping action over a certain frequency range. As
per Griffin et al. [13], if the displacements to be damped are small (sub micrometer)
then squeezing of a thin gas film between two parallel flat surfaces can produce substantial
13
Figure 2.3: Analogy of a squeeze film as a spring-damper system
damping forces at very high frequencies. The damping force is proportional to the relative
velocity over certain ranges of operation, and this type of damper is termed a viscous damper
[13]. Squeeze film damping due to relative axial or tilting motion between two closely spaced
plates is analyzed by Griffin et al. [13]. Griffin et al. [13] provides a critical frequency below
which the squeeze film acts as a damper and above which it acts as a spring. At the critical
frequency, both the damping and spring forces are equal. Analytical expressions are also
provided by Griffin et al. [13] for the squeeze film force and critical frequency for special
cases of infinitely wide parallel plates, annular parallel plates and parallel disks. The critical
frequency for parallel disks is given by [13]
?c= ?hm
2pa
2.07?R02 (2.17)
Blech [14] also analyzed squeeze film cut-off frequencies for different geometries and
divided the squeeze film force into damping and spring force components. The spring force
14
is greater than the damping force for all frequencies except the cut-off frequency, and below
the cut-off frequency damping is comparable to spring action [14]. As per Blech [14], the
maximum damping occurs at the cut-off frequency when the spring and damping forces are
equal, and above the cut-off frequency, the spring force increases while the damping force
decreases with increases in the squeeze number, ?.
Blech [14] gives the damping force as a function of the squeeze number, ?, and excursion
ratio, epsilon1, between parallel circular disks
Fd=?
radicalbigg2
? ?
bracketleftbigA
c
parenleftbigber
1
???bei
1
??parenrightbig+B
c
parenleftbigber
1
??+bei
1
??parenrightbigbracketrightbigepsilon1 (2.18)
Similarly, the spring force by [14] is
Fs=1+
radicalbigg2
? ?
bracketleftbigA
c
parenleftbigber
1
??+bei
1
??parenrightbig+B
c
parenleftbigber
1
???bei
1
??parenrightbigbracketrightbigepsilon1 (2.19)
Here,
Ac= bei
??
parenleftBig
ber2??+bei2??
parenrightBig (2.20)
15
and
Bc=? ber
??
parenleftBig
ber2??+bei2??
parenrightBig (2.21)
The critical frequency for parallel disks is given by [14] as
?c= hm
2pa
1.93?R02 (2.22)
For isothermal conditions, Blech [14] (Eq. (2.22)) underestimates the cut-off frequency by
7 % as compared to Griffin et al. [13] (Eq. (2.17)) because Blech [14] assumes one term
approximation to the cut-off frequency. A sample cut-off frequency using Eq. (2.17) [13] and
assumption of isothermal flow for an experimental result is calculated to be 1377 Hz. For
this case (See Apendix C.2), hm is 13.38 ?m and the frequency of oscillation is 800 Hz. Thus,
the frequency at which the bearing is operated is less than the cut-off frequency where the
spring force is more than the damping force [14]. Similarly, a few more calculations based
on the experimental results suggest that in the present work the range of operation of the
squeeze film bearing is either below or above the cut-off frequency, and for such cases the
spring force is always greater than the damping force [14].
The gas squeeze film stiffness and damping torques on a circular disk oscillating about
its diameter were analyzed by Ausman [15]. In [15], the linearized Reynolds equation for
small squeeze film motions is solved for pressure between the disks where one disk oscillates
in a tilting motion about its diameter. This pressure is then integrated over the surface
16
area to obtain the total squeeze film torque. Then Ausman [15] separates the total squeeze
film torque into stiffness and damping components. As per Ausman [15], the torque which
opposes the angular deflection is stiffness torque; whereas, the torque which opposes angular
rate is damping torque. The gas damping is observed as viscous friction of the gas as it flows
in and out between the disks, and gas stiffness is observed as compressibility of a trapped gas
between the disks [15]. Ausman[15] concludes that at higher frequencies, the gas is trapped
and does not leak, resulting in higher compressibility. This means that higher frequencies
produce higher stiffness torque, and at lower frequencies, gas has more than sufficient time
to flow in and out, resulting in higher damping torque[15].
Etsion [16] analyzed squeeze film effects in liquid lubricated radial face seals and ob-
tained damping coefficients. Green et al. [17] calculated dynamic damping and stiffness
coefficients of the fluid films in mechanical face seals, considering squeeze film effects along
with hydrostatic and hydrodynamic effects. Work specific to compressible squeeze film
damping was done by Blech [18] for annular squeeze-film plates in relative motion. Re-
search into the effect of squeeze film damping is also currently prominent in the design of
micro-electro mechanical systems (MEMS) and microstructures ([19]-[26]). Some of this
work is discussed in detail in the next section.
17
2.7 Micro-Scale Squeeze Film Bearings and Effects Due to Molecular Dynam-
ics
Micro-Electro-Mechanical Systems (MEMS) are miniature systems used to combine
electro-mechanical functions with dimensions varying from a micrometer to a millimeter.
Due to the micrometer size, lubrication in MEMS becomes a critical design parameter
and design and selection of MEMS bearings is a challenge for researchers. Usage of liquid
lubricant in MEMS leads to a power dissipation problem and so is not a good choice for
lubrication [27]. An alternative to oil-based lubricant, gas lubricated bearings can be used
in MEMS. Gas bearings can support their loads on pressurized thin gas films. As per
Epstein [28], for micromachines such as turbines, gas bearings have several advantages over
electromagnetic bearings, such as no temperature limits, high load carrying capability, and
relatively simple fabrication. ?The relative load-bearing capability of a gas bearing improves
as size decreases since the volume-to-surface area ratio (and thus the inertial load) scales
inversely with size? [28]. An example of a fully functional gas film bearing is seen in the MIT
Microengine project [29]. Here, a rotor of a micro-gas turbine generator is supported by a
journal air bearing. As per Breuer [30], the gas lubrication system in MEMS should be easy
to fabricate with sufficient performance and robustness. In another example, a self-acting
gas thrust bearing was designed, fabricated and tested on a silicon microturbine [31]. Wong
et al. [31] compared a hydrodynamic gas thrust bearing to an existing hydrostatic one,
and observed that a hydrodynamic approach is much simpler to fabricate and the required
source of pressurized gas can be eliminated.
18
If the high speed relative normal motion is already available in MEMS, then potentially
the squeeze film effect can be used as a means to lubricate surfaces in MEMS, and create
MEMS bearings. In such cases, there is no special geometry needed as two planar surfaces
serve as a bearing, and the system is also maintenance free because ambient air serves as
the lubricant. Breuer [30], in a summary on some of the issues of lubrication in MEMS,
suggests that fundamental issues of fluids and solid physics such as gas surface interactions,
momentum and energy accommodation phenomenon and surface contamination effects are
vital parameters for ultra-thin lubrication. For such a small scale, continuum assumptions
are not always valid, which can be seen in the following calculations. The RMS roughness,
Rq, of MEMS surfaces seems rather smooth. The range of Rq is usually varies from 0.07 to
0.25 nm [32]. For a complete surface separation, film thickness of 10?Rq is a good approx-
imation. Thus, a minimum film thickness of 2.5 nm should provide a full film lubrication.
The above calculations show that ultra-thin fluid films can be used in MEMS as a means
of lubrication.
The mean free path of a gas is given as [33]
?=16?5p
radicalBigg
GTa
2pi (2.23)
The mean free path can be estimated for the normal operating conditions (See Table 2.1)
in a squeeze film bearing. Using Eq. (2.23) and the assumptions in Table 2.1, the mean free
path, ?, is calculated to be 8.7 nm. The mean free path is used to calculate the Knudsen
19
Table 2.1: Assumptions and parameters used to estimate the mean free path
Flow Isothermal
Air Temperature, Ta 293 K
Viscosity of air, ? 1.82?10?5 N?s/m2
Pressure, p 101325 N/m2
Gas Constant, G 8.3144 J?mol?1?K?1
number to further characterize the type of fluid flow. The Knudsen number, Kn is defined
as
Kn= Mean free path (?)film height (h) (2.24)
If the mean free path is very small compared to the characteristic length (i.e. length scale
of the problem), then the fluid flow is considered to be in the continuum range. If the
mean free path is comparable to the characteristic length then continuum theory of fluid
mechanics does not hold well. Typical values of the Knudsen number corresponding to
different types of flow are tabulated in Table 2.2. Substituting values of ?=8.7 nm and
Table 2.2: Types of flows as per Knudsen number [34]
Range of Knudsen Number Type of Flow
Kn? 10?3 Continuum Flow
10?3 <Kn< 0.1 Slip Flow
0.1 <Kn< 10 Transition Flow
10 ?Kn Molecular Flow
h=2.5 nm in Eq. (2.24) gives Kn?3.5. Thus, the mean free path is not negligible when
compared to the film thickness and the flow between bearing surfaces cannot be considered
20
as a continuum any more. For Kn?3.5, the fluid flow falls into the transition flow regime
(See Table 2.2).
If the film height is increased to 8 ?m then the fluid flow will at most change to slip
flow. Hsia et al. [35] experimentally studied gas bearings at ultra low clearances and pointed
out that for ultra-thin films (below 0.25 microns), slip-flow theory (effects due to molecular
dynamics) needs to be considered when modeling the fluid dynamics. There are various
models which provide means to include rarefaction effects due to reduction in the density
of gas (i.e. slip flow or transition flow regime) in the Reynolds equation ([24],[33],[36],[38]).
Although most of this work is done for the squeeze film damping, the governing equa-
tions generally remain the same and can be used for squeeze film lubrication. Li [36] modeled
an ultra-thin gas squeeze film using the modified gas film lubrication (MMGL) equation in-
cluding coupled roughness and rarefaction effects. A rectangular geometry is considered in
this work and a linearized form of MMGL equation [36] is solved.
?
?x
parenleftbiggparenleftBig
?XP
parenrightBig
0
?P
?x
parenrightbigg
+ ??y
parenleftbiggparenleftBig
?Y P
parenrightBig
0
?P
?y
parenrightbigg
=?2
parenleftbigg
??H?t +?P?t
parenrightbigg
(2.25)
Here,
?2 = 12??h
m2patm
(2.26)
and ?XP and ?Y P are pressure flow factors in x and y direction respectively.
21
Then, Li [36] transforms the linearized MMGL problem to the continuum gas film
problem to obtain an analytic solution. Darling et al.[37] provides an analytical equation
for the squeeze film force in the form of a complex number. The real part of the complex
number is the spring force while the imaginary part is the damping force. In an another work
by Pandey et al. [24], surface roughness and rarefaction effects are included in the analysis
of squeeze film damping in MEMS. This work [24] accounts for nonlinear effects rather than
solving the linearized form of the MMGL equation. An explicit finite difference solution of
the nonlinear MMGL equation suggests that rarefaction effects reduce the spring force and
the damping force below the cut-off frequency but above the cut-off frequency the spring
force reduces but the damping force increases when compared with the linearized solution
[24]. Rarefaction effects increase the value of the cut-off frequency and so increase the range
of frequency where damping is dominant [24]. Thus, in order to achieve lubrication from
the squeeze film effect, it should be better to operate the squeeze film bearings above the
cut-off frequency.
Nayfeh et al. [23] modeled flexible microstructures under the effect of squeeze film
damping. The form of Reynolds equation used by [23] is
?
?x
parenleftbigg
H3P ?P?x
parenrightbigg
+ ??y
parenleftbigg
H3P ?P?y
parenrightbigg
=12?eff ?(PH)?t (2.27)
22
Here, the effective viscosity, ?eff, is given by Veijola et al. [38] which accounts for the
gas rarefaction effects and is given as
?eff= ?1+2Kn+0.2?Kn0.788e?Kn/10 (2.28)
Eqs. (2.27) and (2.28) can be easily used for modeling the squeeze film effect for ultra-thin
films and therefore the present work models ultra-thin squeeze film effect using Eqs. (2.27)
and (2.28) (See Chapter 4).
23
Chapter 3
Objectives
The objective of this research is to experimentally and analytically study dynamic
compressible squeeze film bearings. The squeeze film bearings are hydrodynamic bearings
where fluid pressure is generated due to relative normal motion between the bearing surfaces.
The dynamic motion of gas squeeze film bearings is governed by fluid flow forces, forces due
to motion of both the bearing surfaces and the contact forces due to surface asperities. The
lubricant used in studied compressible squeeze film bearings is air or gas. Gas lubricated
bearings have certain advantages over the liquid lubricated bearings such as low frictional
heat and a greater temperature range [39].
In the first part of this thesis, macro-scale dynamic gas squeeze film bearings having
thin squeeze films will be studied. The coupled dynamics of macro-scale gas squeeze film
bearings which includes the Reynolds equation for fluid flow, the equation of motion for the
squeeze film bearing and the contact force due to initial surface asperity interaction will
be formulated. Then, coupled equations will be numerically solved to study the dynamic
behavior of gas squeeze films from its initial state to pseudo-steady state. Based on numer-
ical investigation, the load carrying capacity to the squeeze number, ? and the excursion
ratio, epsilon1 will be correlated. A test setup will be built to experimentally generate, control and
measure the gas squeeze films between disk shaped surfaces. Experimental and numerical
results will be compared.
24
In the latter part of this thesis, micro-scale dynamic gas squeeze film bearings having
ultra-thin squeeze films will be studied. The fluid dynamics for micro-scale dynamic gas
squeeze film bearings (ultra-thin) as it changes from macro-scale dynamic gas squeeze film
bearings (thin) will be modified, and the changes will be incorporated into the coupled
dynamic model. Prototypes of micro-scale squeeze film bearings will be fabricated for the
experimental study. Finally, experimental and simulation results will be compared.
25
Chapter 4
Numerical Investigation
4.1 Introduction
This chapter describes in detail the numerical scheme to solve the coupled dynamic
system. The results obtained from the numerical investigation are then discussed. The
chapter is divided into two parts, thin squeeze film bearings and ultra-thin squeeze film
bearings.
4.2 Thin Squeeze Film Bearings
If the squeeze film is thin enough to operate in the continuum flow regime then the
bearing is considered to have a thin squeeze film (as opposed to an ultra-thin squeeze film).
This section covers the formulation of the coupled dynamics, the computational scheme and
the simulation results.
4.2.1 Formulation of Coupled Dynamics
See Fig. 4.1 for a schematic of the dynamic planar gas squeeze film bearing. The
bearing configuration consists of two metal disks. Ideally, these metal disks are subjected
to only one degree of freedom (permitting motion only in the vertical direction). The
bottom disk is used as an oscillating base and the top disk freely levitates just above the
bottom disk. The oscillating base is achieved by mounting on a vibrating shaker. At rest,
26
the upper disk sits on the top surface of the oscillating base. Once the base is subjected to
adequate oscillatory motion, the air squeeze film is generated between the oscillating base
and the upper disk.
Levitating disk 
Air squeeze film
Oscillating base
z
zt
zb=Z0sin(T)
h
2
R
Figure 4.1: Schematic of a planar squeeze film bearing
Assumingthatthesqueezefilmisaxisymmetricandisothermal, thenormalizedReynolds
equation in polar coordinates for an ideal gas between two parallel circular disks is given as
1
R
?
?R
parenleftbigg
RH?3P ?P?R
parenrightbigg
= ???(PH
?)
?T (4.1)
where the squeeze number, ??, is as given
?? = 12?R02?/patmh02 (4.2)
27
Eq. (4.1) differs from Eq. (2.9) in the manner in which the squeeze film height is
normalized. In Eq. (4.1) the squeeze film height, h, is normalized by the initial film height,
h0, whereas in Eq. (2.9) the squeeze film height, h, is normalized by the mean film height,
hm (see nomenclature). The initial film height, h0, is known; whereas the mean film height,
hm, is unknown before solving the fluid dynamics. Thus, the fluid dynamics is modeled for
Eq. (4.1). Eq. (4.1) is a nonlinear parabolic partial differential equation. The boundary
conditions of Eq. (4.1) are given by Eqs. (2.14) and (2.15). The initial pressure is ambient
and the initial film thickness is determined from a contact force model (see Appendix A).
The normalized Reynolds equation (Eq. 4.1) is solved numerically.
Figure 4.2: Scheme for discretization of i-spatial variable and j-time variable
The left hand side of Eq. (4.1) is discretized using a central-finite difference scheme.
Fig. 4.2 explains the meshing scheme for the squeeze film bearing in space (i) and time (j).
For (0<R<1),
28
1
R
?
?R
bracketleftbigg
PRH?3?P?R
bracketrightbigg
? H
?j+13
Ri
bracketleftbigg
Ri+1
2
Pi+1
2,j+1
parenleftbiggP
i+1,j+1 ?Pi,j+1
?R
parenrightbiggbracketrightbigg
?H
?j+13
Ri
bracketleftbigg
Ri?1
2
Pi?1
2,j+1
parenleftbiggP
i,j+1 ?Pi?1,j+1
?R
parenrightbiggbracketrightbigg
(4.3)
where,
Ri+1
2
= Ri+Ri+12 (4.4)
Ri?1
2
= Ri+Ri?12 (4.5)
Pi?1
2,j+1
= Pi,j+1 +Pi?1,j+12 (4.6)
Pi+1
2,j+1
= Pi,j+1 +Pi+1,j+12 (4.7)
and at the left boundary (R=0), the pressure gradient ?P?R is zero.
P0,j = P0,j+1 (4.8)
At the right boundary (R=1), the pressure is always atmospheric.
P(R = 1,T) = 1 (4.9)
The right hand side of Eq. (4.1) is discretized as
???(PH
?)
?T ? ?
?Pi,j+1H?j+1 ?Pi,jH?j
?T (4.10)
29
In order to investigate the dynamic squeeze film bearing, it is necessary to consider the
coupled dynamics of the system, which includes the Reynolds equation and the equation of
motion using a force balance.
2
2
dt
zdm t
mg( ) ApP atmn 1? contF
Levitating disk
Figure 4.3: Free body diagram of a squeeze film bearing (one degree of freedom)
A free body diagram of the top levitating disk (Fig. 4.3) shows all the forces acting on
the top disk. Thus, the equation of motion can be written as
md
2zt
dt2 = (Pn ?1)?patm ?A+Fcont ?Fi (4.11)
Here, Fi is the weight of the top disk as employed in the simulation i.e.
Fi = mg (4.12)
30
Substituting Eq. (4.12) into Eq. (4.11) gives
md
2zt
dt2 = (Pn ?1)?patm ?A+Fcont ?mg (4.13)
Here, the term Fcont represents the asperity contact force. The contact force model is
presented in the Appendix A. In the equation of motion, contact force between bearing
surfaces needs to be considered since the surfaces can be in contact at start-up, shut-down
and if the squeeze film force is not large enough to achieve continuous separation of the
surfaces. The model detailed in the Appendix A provides an average contact force, Fcont,
as a function of surface separation or film thickness (h), rms surface roughness (Rq), modulus
of elasticity (E), yield strength (Sy) and Poisson?s Ratio (?). Using an exponential fit to
the results of Eqs. (A1-A7), an analytical expression of Fcont as a function of h is obtained.
In the numerical simulation, this relationship is used to predict the contact force between
the bearing surfaces. Initially, at t = 0, the two surfaces are stationary and the contact
force, Fcont, balances the applied load, Fi. But during the development of the squeeze film,
the two surfaces are pushed out of contact by the reaction force between bearing surfaces.
Inclusion of Fcont in the dynamics also predicts the initial (at t = 0) film thickness h0 when
Fcont=mg, and also indicates when a full film of lubrication is achieved (the contact force
reaches zero). As the squeeze film thickness, h, is also the relative displacement between
the two bearing surfaces (See Fig. 4.1)
31
zt = zb +h (4.14)
Then taking the time derivative of Eq. (4.14) twice yields
d2zt
dt2 =
d2zb
dt2 +
d2h
dt2 (4.15)
and since the base is excited by a sine wave
zb = Z0sinT (4.16)
Differentiating Eq. (4.16) twice with respect to t yields
d2zb
dt2 = ??
2Z0sinT (4.17)
Substituting Eqs. (4.15) and (4.17) into Eq. (4.13) yields
md
2h
dt2 = m?
2Z0sinT +(Pn ?1)?patm ?A+Fcont ?mg (4.18)
Eq. (4.18) can be normalized for h and t. Here, h is normalized to H? (h/h0) and t to
T (?t). The normalized form of Eq. (4.18) is
d2H?
dT2 =
1
h0?2
parenleftbigg
?2Z0sinT + (Pn ?1)?patm ?Am + Fcontm ?g
parenrightbigg
(4.19)
32
Eq. (4.19) is a second-order nonlinear partial differential equation and its initial con-
ditions (at t=0) are given as
dH?
dT = 0 (4.20)
Pn = 1 (4.21)
Fcont = mg (4.22)
The fourth-order Runge-Kutta method is used to solve the equation of motion (Eq.
(4.19)). Eq. (4.19) is also coupled to the fluid dynamics of the squeeze film through the
dimensionless mean pressure, Pn, calculated from the discretized Reynolds equation (Eqs.
(4.3) and (4.10)).
4.2.2 Computational Scheme
The coupled dynamics (the Reynolds equation and the equation of motion) is solved
simultaneously. Since the Reynolds equation is parabolic it must be solved at each individual
step of the Runge-Kutta method. Fig. 4.4 illustrates the steps of an algorithm to solve the
coupled dynamics. First, all constants, initial and boundary conditions are defined. Then
the fourth-order Runge-Kutta method is implemented to solve the equation of motion.
For each step of the Runge-Kutta method, the discretized Reynolds equation is solved to
33
Define constants, initial conditions and boundary conditions
Step k of 4th order Runge-Kutta method to solve equation of motion for H
Solve discretized Reynolds equation for P
Is pseudo-steady state achieved?
No
Yes
Stop
k=1
If k<4
k=k+1
Increment T
Update P and H
Runge-Kutta Method
Figure 4.4: Algorithm for computational simulation
34
determine change in pressure due to change in the film height. This change in pressure is
used in the next step of the Runge-Kutta method to account for change in the film height.
Pressure and film thickness are updated at the last step. This process is continued until a
pseudo-steady state is achieved. The simulation code is programmed in C (see Appendix
B).
In order to confirm that the mesh density is adequate in space, only the discretized
Reynolds equation is solved with an assumption that ? = 1000 and the squeeze film height
is known as a function of time, i.e. H=1-0.5?sin(T). P as a function of T is plotted for 16
(Fig. 4.5) and 160 nodes (Fig. 4.6). From the plots (Figs. 4.5 and 4.6), it can be seen that
the 16-node grid is as accurate as the 160-node grid. The deviation of pressure for grid size
of 16 to that of 160 is only 0.75 %. Thus, the 16-node grid is used in numerical simulations
in order to reduce computational effort.
35
?2?3?4?
(rad)
0
Figure 4.5: Dimensionless P as a function of normalized T at R=0.5, ?=1000 and H=1-
0.5?sin(T) (16 nodes)
?2?3?4?
(rad)
0
Figure 4.6: Dimensionless P as a function of normalized T at R=0.5, ?=1000 and H=1-
0.5?sin(T) (160 nodes)
36
After defining the grid size, the value of the time-step ?t is refined by simultaneously
solving the coupled dynamics. Fig. 4.7 shows that ?t=?T/?=10?7 s is small enough
(deviation with step size, ?t of 10?8 s is 0.2214 %) and thus used to characterize squeeze
film dynamics, but some error (deviation with step size, ?t of 10?8 s is 2.3868 %) occurs
for ?t=?T/?=10?6 s.
2.5 2.5001 2.5002 2.5003 2.5004 2.500576
78
80
82
84
86
88
90
t (s)
h(
?m)
?t=?T/?=10    s?6
?t=?T/?=10    s?7
?t=?T/?=10    s?8
Figure 4.7: Variation in film height for different time steps
37
4.2.3 Results
The numerical simulation is performed from t=0 to the time when the squeeze film
thickness achieves a pseudo-steady state (see Fig. 4.8). The simulation input parameters
used are f =5000 Hz, Y=5 ?m, R0=0.0325 m and Fi=1.57 N. At t=0, film thickness is
calculated solely from the contact force model and is used as an initial film thickness in the
simulation. The base oscillation is initiated and the simulation is run to a pseudo-steady
state.
Dynamic behavior of the squeeze film
From Fig. 4.8, it can be observed that between time 0 to 0.01 seconds the surfaces come
in and out of contact and the contact force is dominant in this region. After 0.01 seconds,
a full film is developed and the film height continues to increase until approximately 0.8
seconds. After approximately 0.8 seconds, the squeeze film achieves a pseudo-steady state
where it oscillates at a mean squeeze film height (see Fig. 4.8 enlarged view). Simulation
output parameters at pseudo-steady state are hm=81.52 ?m, epsilon1=0.063 and ?=10.7.
Comparison with the pre-established models
Simulation results are used to verify the predictions of the squeeze film force model
derived by Langlois [8]. The excursion ration, epsilon1, and the squeeze number, ?, are substituted
into Eq. (2.11) [8] and the squeeze film force, F(T), is obtained. The squeeze film force,
F(T), from simulation result is [(P ?1) ? patm ? A]. Fig. 4.9 shows that Langlois [8] and
38
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
t (s)
h (
?m)
Simulation inputs
Fi=1.57 N
R0=0.0325 m
f=5000 Hz
Z0=5 ?m
Simulation outputs
hm=81.5 ?m
?=0.063
?=10.7
Enlarged View
Figure 4.8: Dynamic behavior of the squeeze film height for given input conditions as a
function of time
39
0 2 4 6 8 10 12 14?0.06
?0.04
?0.02
0
0.02
0.04
0.06
0.08
T (rad)
F(T)
/(pi
?R 02
?p at
m)
Langlois Solution
Numerical Simulation
Figure 4.9: Comparison of the dimensionless squeeze film force, F(T)/(pi?R02?patm), as a
function of normalized time, T, from the numerical simulation and from Langlois squeeze
film force model [8]
40
present model compare moderately well and so [8] approximately predicts the periodic
squeeze film force. Since F(T) given by [8] is a function of sine and cosine (Eq. (2.11)),
the net squeeze film force over one period is zero. This does not hold true because the
numerical analysis and Salbu [9] predict a net load carrying force is formed. Thus, F(T)
by [8] is a good estimation but is not entirely accurate. As per the present model, the net
squeeze film force over one period equals to the weight of the levitating disk.
The excursion ratio, epsilon1, and the squeeze number, ?, calculated from the current simu-
lation results (See Fig. 4.8) can also be used to verify and correlate with the predictions
of Salbu?s Eq. Substituting values R0=0.0325, f=5000 Hz and hm=81.5 ?m (See Fig. 4.8)
into Eq. (2.10) gives squeeze number, ?=10.7. As ? is greater than the limiting minimum
value of 10 defined by Salbu [9], Salbu?s Eq. can be used to calculate the load carrying
capacity. A comparison of the applied load, Fi and predicted load, Fn by Salbu?s Eq. is
shown in Table 4.1.
Table 4.1: Comparison of load carrying capacity between Fi in numerical simulation and
Fn by Salbu?s Eq.
Fi 1.57 N
Fn by Salbu?s Eq. 1.67 N
for epsilon1=0.063
Percentage error 6.4
between Fi and Fn
As the percentage error between Fi and Fn by Salbu?s Eq. is only 6.4; it appears that
the coupled model is in agreement with the pre-established model [9] and can be used to
study the dynamic behavior of squeeze film bearings. More extensive comparisons of the
simulation results to Salbu?s predictions are presented in the next section.
41
Development of Semi-analytical Equation
Simulations are performed by changing the simulation parameters such as f, R0, Z0
and Fi. Thus a wide range of values for ? (0.4<?<14.6) are obtained. Also the percentage
error between Fi and Fn by Salbu?s Eq. is calculated for these simulations. For ?>10,
it is observed that the percentage error between Fi and Fn by Salbu?s Eq. is fairly low
(less than 7). However, for ?<10, the percentage error increases and over-predicts the load
carrying capacity (Fig. 4.10). Thus, Salbu?s Eq. can be used as a basis to calculate load
carrying capacity if ? is greater than 10, but if it is less than 10, the percentage error
increases significantly and a modified form of Salbu?s Eq. is needed. In order to obtain
a semi-analytical equation to extend the range of Salbu?s Eq., it is modified by the curve
fitting technique using the percentage error between Fi and Fn as predicted by Salbu?s Eq.
Here, the level of certainty of the fit is 95 % and value of R-square is 0.9985. An exponential
equation f(?) is fit to the percentage error of the data as a function of ? (see Fig. 4.10)
and is given by
f (?) = 460?exp(?4.5??)+19?exp(?0.15??) (4.23)
Wn= Fnpip
atmR02
=
parenleftbigg 1
1+f (?)
parenrightbigg??
?
bracketleftBigg1+ 3
2epsilon1
2
1?epsilon12
bracketrightBigg1/2
?1
??
? (4.24)
A comparison between Fi and Fn using Eq. (4.24) for a large number of points within
0.4<?<10 is performed. It is observed that the average percentage error between Fi and
42
2 4 6 8 10 12 14
0
20
40
60
80
100
120
140
160
?
percentage error between  
F i 
and 
F n
 by Salbu?s Eq.
percentage error between Fi  and  Fn by Salbu?s Eq.
Exponential fit to the data
Figure 4.10: Change in percentage error between Fi from the numerical simulation and Fn
as predicted by Salbu?s Eq. as a function of the squeeze number, ?
43
Fn from Eq. (4.24) for 0.4<?<10 is only 1.04. Thus Eq. (4.24) confirms its usefulness for
lower ? values in this range.
Parametric Study
A parametric study using the analytical simulation is performed to study the behavior
of the mean squeeze film height as a function of Fi (Fig. 4.11). As shown, hm decreases with
increases in Fi, for a given set of simulation parameters. It can also be observed from Fig.
4.11, that if the frequency is doubled and other parameters kept constant, the mean squeeze
film height, hm, changes but relatively less than when load is varied. The variation of hm
with Fi is similar in nature for both the frequencies. Numerical simulations are performed
to study the behavior of the mean squeeze film height as a function of the amplitude of
vibration, Z0. An experimental investigation is also performed to study the behavior of the
mean squeeze film height as a function of Z0. Both numerical and experimental results are
compared, and so the results are presented in the next chapter on the experimental work.
The behavior of the mean squeeze film height, hm, as a function of frequency, ?, is
also studied by keeping all other simulation parameters (m=6.97 gm, R0=9 mm and Z0=
5 ?m) constant. It can be observed from Fig. 4.12, that if the frequency, ? increases, the
mean film height also increases. For lower frequencies, the slope of the graph is higher
which indicates more increase in the mean squeeze film height with increase in frequency.
For higher frequencies, the slope of the graph is lower which indicates less increase in the
mean squeeze film height with increase in frequency.
44
Fi (N)
hm (?m)
Figure 4.11: The results of numerical parametric study for hm as a function of Fi at a
constant frequency
45
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x 105
86
88
90
92
94
96
98
100
102
104
? (rad)
h m
 (?
m)
Figure 4.12: The results of numerical parametric study for hm as a function of ?
46
4.3 Ultra-Thin Squeeze Film Bearings
From the background (Section 2.7), if the squeeze film bearing operates in the slip, tran-
sition or molecular flow regime then the bearing is considered as an ultra-thin squeeze film
bearing. This section covers the formulation of the coupled dynamics and computational
scheme as it differs from the thin squeeze film bearings.
4.3.1 Formulation of Coupled Dynamics
The dynamics of the ultra-thin squeeze film bearings is generally the same as the thin
squeeze film bearings except in the ultra-thin squeeze film bearings, a modified form of the
Reynolds equation is used.
1
R
?
?R
parenleftbigg
RH?3P ?P?R
parenrightbigg
= ??eff ?(PH
?)
?T (4.25)
where ??eff is as given
??eff = 12?effR02?/patmh02 (4.26)
The effective viscosity, ?eff, is given by Eq.(2.28) [38] which accounts for the gas rarefaction
effects. The discretization scheme for the left hand side of Eq. (4.25) and the boundary
conditions are the same as that used for the thin squeeze film bearings (See Eqs (4.3)-(4.9)).
47
The right hand side of Eq. (4.25) is discretized as
?eff??(PH
?)
?T ? ?eff
?Pi,j+1H?j+1 ?Pi,jH?j
?T (4.27)
The equation of motion, Eq.(4.19), is the same for both thin and ultra-thin squeeze film
bearings. The initial conditions are given by Eqs. (4.20)-(4.22).
4.3.2 Solution Methodology
The algorithm (Fig. 4.13) to solve the coupled dynamics (the modified Reynolds equa-
tion Eq. (4.25) and the equation of motion Eq. (4.19)) for the ultra-thin squeeze film
bearing is same as the thin squeeze film bearings (Eqs. (4.1) and (4.19)) except for one
change. Here, the modified squeeze number, ??eff is calculated each time before solving
the discretized Reynolds equation.
4.3.3 Results
A computer program (see Appendix E) as per the algorithm detailed in Fig. 4.13 is
written to solve the discretized Reynolds equation (Eqs. 4.22 and 4.24). Here, the squeeze
film thickness is assumed as a known function of time and only the modified Reynolds
equation (Eqs. 4.22 and 4.24) is solved. Simulation and experimental results are compared
and the results are presented in the next chapter (see Section 5.2.3).
48
Define constants, initial conditions and boundary conditions
Step k of 4th order Runge-Kutta method to solve equation of motion for H
Solve discretized Reynolds equation for P
Is pseudo-steady state achieved?
No
Yes
Stop
k=1
If k<4
k=k+1
Increment T
Update P and H
Find Kn and calculate ? eff and 
1 *eff
Runge-Kutta Method
Figure 4.13: Algorithm for computational simulation of the ultra-thin squeeze film bearing
49
Chapter 5
Experimental Investigation
5.1 Thin Squeeze Film Bearings
This section covers in detail the effort taken towards reaching the goal of generation,
control and measurement of the squeeze film between two bearing surfaces. The different
techniques for the measurement of the squeeze film height are described here followed by a
discussion of the results.
5.1.1 Experimental Setup
A test rig (See Fig. 5.1) has been built in order to generate a squeeze film and measure
the squeeze film thickness. The experimental setup is used to perform a parametric inves-
tigation with the variation of frequency, amplitude of vibration and mass of the levitating
disk. Each subsystem of the experimental setup is described in detail.
Figure 5.1: Schematic explaining the experimental setup
50
Test Stand
strip of tape
base surface
top disk
electrodynamic shaker
Figure 5.2: Test stand used for experimental investigation
A rectangular bracket holding the lower bearing surface is used as an oscillating base
(See Fig. 5.2). On the top of the lower bearing surface a disk is placed. The top disk is
tethered to the rectangular bracket using four small strips of tape. This attachment is done
in such a way that the top disk can only float with one degree of freedom i.e permitting the
motion only in vertical direction. This setup may result in additional load being placed on
the bearing. Different disk sizes used for levitation can be seen in the Fig. 5.3 and their
specifications are tabulated in Table 5.1.
Table 5.1: Bearing configurations used for the experimental purpose
Configuration Radius (mm) Mass (g)
1 32.5 160
2 7.5 2.62
3 9 6.97
51
Config. 2
Config. 3
Config. 1
Figure 5.3: Disks used for levitation
Electrodynamic shaker
The test stand is fastened to an electrodynamic shaker (LDS Model V408) ( See Fig.
5.4 and 5.6). The electrodynamic shaker provides sinusoidal oscillations to the lower bearing
surface.
Gain Control and Dynamic Signal Analyzer
A gain control (LDS PA 500L amplifier) and a dynamic signal analyzer (HP 3566 5A)
are used together for vibration control (Fig. 5.4 and 5.6). The source from the dynamic
signal analyzer is set as a sine wave. The frequency of the source is set to a fixed value
and the level of the source is changed to vary the amplitude of the vibration. The output
from this source is given to the gain control mechanism which sends the signal to the
electrodynamic shaker.
52
A dynamic signal analyzer is also used to read the amplitude of vibration of the oscillat-
ing base as given by a laser beam displacement measurement system (LBDMS). It records
the analog voltage from the LBDMS. The amplitude of this analog voltage is then converted
into the amplitude of oscillation of the bearing surfaces in microns(?m) by multiplying the
scale set on the LBDMS. The LBDMS is discussed in greater detail later.
Capacitance sensor
Signal analyzer
Optical sensor head
Gain control
Shaker
Test stand
Vibrometer
Data acquisition board
Signal generator
Capacitance sensor
Figure 5.4: Photograph illustrates the experimental setup with the capacitance sensor
If there is a gas squeeze film between two parallel plates then the plates form a parallel
plate capacitor with air as the dielectric. The capacitance between two parallel plates is
53
theoretically given as
C = Aepsilon10kd (5.1)
Thus, a capacitance sensor (Fig. 5.4) is built in order to measure the capacitance between
the base and the levitating disk and then calibrated (Appendix D) to provide the squeeze
film height in microns. First, electrical connections are made to both the bearing surfaces.
The electrical circuit in Fig. 5.5 [46] is built to measure the capacitance due to an air-gap
between these surfaces. Vin is the input voltage to the electrical circuit which is provided
Figure 5.5: Electrical circuit for capacitance measurement
by the signal generator. The frequency of AC input voltage is set to 2000 Hz and RMS Vin
is set to 7.07 V. Rst is a standard resistor with a value of 33 k?. Cairgap is the parallel plate
capacitor formed by the two bearing surfaces across which Vo/p, output voltage is measured
by a data acquisition system, NI BNC 2140. The voltage ratio (RMS output voltage/RMS
input voltage) is needed to calculate the capacitance (see Appendix D). This capacitance
54
between the parallel plates is then calibrated using thin film shims to provide squeeze film
height. Thus, using this experimental set-up, the squeeze film lift and height generated
between the plates can be measured.
Laser Beam Displacement Measurement System
Data acquisition board
Optical sensor head
Signal analyzer
Gain control
Test stand
Shaker
Vibrometer
Figure 5.6: Photograph illustrates the experimental setup with the laser displacement mea-
surement system
The LBDMS used comprises of an optical sensor head (Polytec OFV 2610) and a laser
vibrometer (Polytec OFV 2610)(See Fig. 5.6). The laser beam is incident on the top
levitating disk by an optical sensor head. An analog voltage signal proportional to the
displacement of the top disk is then available at the output of the laser vibrometer. The
scale on the laser vibrometer is set to the finest value of 20 ?m/V and the resolution at this
scale is 0.08 ?m. This analog voltage is given as an input to the data acquisition system
55
(NI BNC 2140). National Instrument?s LabView is used to record the DC estimate of the
input analog signal at a given scan rate in a text file. This text file is later imported and
plotted by MATLAB (Fig. 5.7) for further analysis. The change in DC voltage from the
initial to final state multiplied by 20 ?m/V gives the squeeze film height in microns. For
instance, the average measured film thickness of the experimental results shown in Fig 5.7
is 17.15 ?m.
0 1 2 3 4 5 6 7 8 9 10 11 12?0.5
?0.4
?0.3
?0.2
?0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
t (s)
DC estimate of analog voltage (V)
Change in DC voltage
Base oscillation started
Base oscillation stopped
Initial state
Final state
Figure 5.7: Plot of DC voltage from the laser beam displacement measurement system
56
5.1.2 Experimental Results
Configuration 1
A comparison between the experimental and the numerical results is first made using
the capacitance measurement system on configuration 1 (See Appendix C, Section C1). Re-
sults are plotted for the mean squeeze film height as a function of the amplitude of vibration
for several frequencies (See Fig. 5.8). These numerical results show that, for an increase
in the amplitude of vibration, Z0, there is a corresponding increase in the mean squeeze
film height, hm, if all other parameters remain unchanged. The experimental results for
configuration 1 agree qualitatively but the offset between the experimental and the numer-
ical results is large. hm for configuration 1 is measured using the capacitance sensor. It
is observed that small metal particles or dust between two surfaces can significantly affect
the capacitance measurements, which requires that the surfaces must be frequently cleaned
between tests. It should also be noted that the surfaces of this configuration have significant
imperfections (i.e. roughness and waviness). Even after cleaning the surfaces, during the
start-up of the test, initial contact may create small wear particles and the absolute clean-
liness of surfaces cannot be guaranteed. This makes the capacitance measurements difficult
and also not repeatable, LBDMS is used to record the film thicknesses for configurations 2
and 3 which are also smoother.
57
2 3 4 5 6 7 8 90
20
40
60
80
100
120
140
160
180
Z0 (?m)
h m
 (?
m) f=800 Hz Simulation Results
f=800 Hz Experimental Results
f=1050 Hz Simulation Results
f=1050 Hz Experimental Results
f=1300 Hz Simulation Results
f=1300 Hz Experimental Results
Figure 5.8: Experimental and simulation results for the squeeze film height against the
amplitude of vibration for bearing configuration 1 (Capacitance measurement system)
58
Configuration 2
The plot in Fig. 5.9 shows the squeeze film thickness plotted as a function of the
amplitude of vibration for the base. The measurements are made using LBDMS. Each test
is run ten times in order to assure the repeatability of the setup. The average standard
deviation of the test data (See Appendix C, Section C2) is found out to be 0.31 ?m, and
thus it can be said that experiments show good repeatability. It can be seen from Fig. 5.9
that if the frequency is small then it takes a higher amplitude to generate the squeeze film
height. Conversely, if the frequency is increased then the amplitude of vibration is lower.
These qualitative trends are also agreed by the theoretical results as discussed earlier. For
configuration 2, both simulation and experimental results match qualitatively but the offset
between them is very high. Also, there is a convergence problem for configuration 1 (Fig.
5.9) for lower frequencies and amplitudes. For such cases, the solution did not converge to
a pseudo-steady state even after running simulations for a longer duration.
Configuration 3
Measurements are also carried out for configuration 3 (See Appendix C, Section C3)
using LBDMS and the average standard deviation for configuration 3 is found to be 0.18
?m which suggests that experiments are repeatable. Observations for configuration 2 (Fig.
5.9) and configuration 3 (Fig. 5.10) are similar in nature. The average deviation between
numerical and experimental data is calculated to be 44 ?m (? 73 %) from the results of
configurations 2 and 3. For a similar type of comparison with the results of [10], the average
59
1 2 3 4 5 6 7 8 90
10
20
30
40
50
60
70
Z0    (?m)
h m
(?m)
f=800 Hz Simulation Results
f=800 Hz Experimental Results
f=900 Hz Simulation Results
f=900 Hz Experimental Results
f=1000 Hz Simulation Results
f=1000 Hz Experimental Results
f=1500 Hz Simulation Results
f=1500 Hz Experimental Results
Figure 5.9: Experimental and simulation results for the squeeze film height against the
amplitudeofvibrationforbearingconfiguration2(Laserdisplacementmeasurementsystem)
60
deviation is 38 ?m (? 51 %). From this comparison, it should be noted that both Minikes
and Bucher [10] and a current numerical simulations overestimate the film thickness by
similar magnitudes. This is because the assumptions made in the numerical model cannot
be perfectly implemented while conducting experiments. The numerical model assumes that
the disks are flat and rigid, and that the system has a single degree of freedom (motion only
in the vertical direction). However, during the experimental tests, it is observed that there
is a component of motion acting in the horizontal plane. The surfaces are not perfectly flat
due to surface roughness and waviness and consequently are not perfectly leveled. Thus,
the assumption of a planar squeeze film cannot be experimentally achieved. It is possible to
include the effects due to tilting of the squeeze film in the theoretical work to achieve a better
comparison but it will lead to a much more complex theoretical model. The other possible
reasons for differences between the analytical and experimental results can be attributed
to dust or small wear particles present in between the plates, misalignment between the
surfaces, and the sensitive nature of the laser setup to external vibrations. However, the
experimental results do confirm that the squeeze film effect can be used in real applications
to generate gas squeeze film bearings.
61
1 2 3 4 5 6 7 810
20
30
40
50
60
70
Z0 (?m)
h m
(?m)
f=800 Hz Simulation Results
f=800 Hz Experimental Results
f=900 Hz Simulation Results
f=900 Hz Experimental Results
f=1000 Hz Simulation Results
f=1000 Hz Experimental Results
f=1500 Hz Simulation Results
f=1500 Hz Experimental Results
Figure 5.10: Experimental and simulation results for the squeeze film height against the
amplitudeofvibrationforbearingconfiguration3(Laserdisplacementmeasurementsystem)
62
5.2 Ultra-thin Squeeze Film Bearings
This section covers in detail the design and fabrication of micro-scale squeeze bearing
surfaces which is later used for the testing of ultra-thin squeeze film bearings. LASICAD
software is used for the design of micro-scale bearing surfaces and the fabrication is carried
out in the Alabama Microelectronics Science and Technology Center.
5.2.1 Design of Micro-Scale Bearing Surfaces
See Fig. 5.11 for an enlarged view of examples of the LASICAD drawing for each
surface. Here, the relative scale between different geometries is preserved. Three circular
geometries (Diameter 10, 50, 100 ?m) and one square geometry (50 ?m X 50 ?m) are
drawn using LASICAD. Each sample micro bearing surface represents a grid of a particular
geometry on a 3 cm X 3 cm square. For example, a sample bearing surface (Fig. 5.11
top left) is a grid of 50 ?m diameter circular areas where the centers of two consecutive
areas are separated 100 ?m apart. Other samples are also seen in Fig. 5.11 where top right
represents a grid of 50 ?m X 50 ?m squares, bottom left represents a grid of 10 ?m diameter
circles and bottom right represents a grid of 100 ?m diameter circles. As per the LASICAD
drawing, a TLD (Transportable LASI Drawing) file is utilized to externally fabricate the
photomask. The photomask is a square glass plate having the drawing encapsulated on a
chromium layer.
63
Figure 5.11: LASICAD drawing of micro bearing surfaces (not to scale and only portions
of the textured sections are shown in this enlarged view)
64
5.2.2 Fabrication of Micro-Scale Bearing Surfaces
The process flow of fabrication of micro bearing surfaces is diagrammed in Fig. 5.12.
The steps involved in this fabrication are described as below.
Cleaning of the wafer
Hard bake
Hexamethyldisalizane (HMDS)
Soft bake
Photoresist application
Mask Alignment and Developing
Etching
Dicing
3X3 cm4X4 cm
Photoresist removal
Figure 5.12: Fabrication procedure for micro-scale bearing surfaces
Cleaning of the wafer
A new silicon wafer of p-type and surface orientation < 100 > is cleaned as per the B,
C, and D silicon wafer cleaning procedure (See Table 5.2 [42]).
65
Table 5.2: Silicon Wafer Cleaning Procedure [42]
B Removal of Residual Organic/Ionic Contamination
1 Hold a wafer in a (5:1:1) solution of H2O?NH4OH ?H2O2 for 10
min at a temperature of 75-80 0C
2 Quench the solution under running deionized (DI) water for 1 min
3 Clean a wafer in DI water for 5 min
C Hydrous Oxide Removal
1 Immerse a wafer in a (1:50) solution of HF ?H2O for 15 sec
2 Clean a wafer under running DI water for 30 sec
D Heavy Metal Clean
1 Hold a wafer in a (6:1:1) solution of H20?HCl?H2O2 for
10 min at a temperature of 75-80 0C
2 Quench the solution under running deionized (DI) water for 1 min
3 Clean a wafer in running DI water for 20 min
Hard bake
Hard baking of the wafer is done by an Imperial IV microprocessor oven. The wafer is
kept inside the oven for 20 minutes at 120 0C. The hard bake process removes any moisture
content from the wafer.
Hexamethyl disalizane (HMDS)
After the hard bake, the wafer is kept in a HMDS chamber for 20 minutes. Here, the
wafer surface is primed with HMDS to promote better adhesion to the photoresist.
Photoresist application
Photoresist is a light-sensitive material to which patterns are first transferred from
the photomask. A liquid photoresist (AZ 5214) is applied in a liquid form onto the wafer
66
surface. Then the wafer held on a vacuum chuck undergoes rotation at a speed of 3000 rpm
for 40 sec. This process provides a layer of photoresist of even thickness.
Soft baking
The photoresist-coated wafer is then transferred to a hot plate for soft baking or pre-
baking. Soft baking is performed on a hot plate at 105 0C for 1 minute. It improves the
adhesion of the photoresist to the wafer and also drives off solvent from the photoresist
before the wafer is introduced into the exposure system.
Mask Alignment and Exposure
In this step, the photo mask is aligned with the surface of the wafer. The wafer is held
on a vacuum chuck, and moved into position below the photo mask. The spacing between
the photo mask and wafer surface is in the range of 25 to 125 ?m. Following alignment, the
photoresist is exposed for 10 seconds with high-intensity ultraviolet light. After this step,
the wafer is developed in the AZ 514 developer.
Etching
Etching is performed to remove material between the circular and square areas so that
micro-scale bearing areas are formed on the wafer (in the form of posts). These micro-scale
bearing areas can be observed as posts protruding out from the base silicon. Deep reactive
ion etching is used to etch out the wafer. The machine used for this process is Surface
Technology System?s (STS) advanced silicon etcher. An etch depth of approximately 80
67
400 ? m
20 ? m

- 10 ? m
80 ? m
c
400 ? m
200 ? m

- 100 ? m
80 ? m
d
400 ? m
100 ? m

- 50 ? m
80 ? m
a
400 ? m
100 ? m
50 ? m
b
80 ? m
Figure 5.13: Schematics of single unit cell (not to scale)
68
?m is achieved from 120 cycles of Morgan SOS 1 process. The process took 40 minutes to
complete 120 cycles. The depth of the etch is checked using a microscope.
Photoresist removal
Photoresist is stripped off from a wafer using a Matrix machine.
Dicing
Dicing is used to cut the wafer and separate out the four different samples(Fig. 5.12).
The wafer is attached to a plastic film before starting the dicing operation. The plastic film
is then supported on a steel rim. Micro Automation?s dicing saw is used. After dicing, four
samples of micro bearing surfaces with different sizes and geometries are collected in a petri
dish. A schematic of the resulting unit cell can be seen in Fig. 5.13.
5.2.3 Results
Measurements of micro-scale bearing samples using LBDMS encountered many prob-
lems. In the first approach, the bearing sample was tethered to the rectangular bracket
using four small strips of tape. This approach produced the squeeze film lift upon start of
the base oscillation. However, after stopping the base oscillation the squeeze film lift did
not always return to zero as the weight of the micro-scale bearing sample was supported by
the attached four strips. Another method of testing was tried to test the sample. Here, the
sample was constrained to translate in horizontal plane. This was achieved by placing small
circular posts of tape at the midpoints of each side of the sample so that the circular post
69
and sample have in total four point contacts. This test method also did not work as the
friction between the circular posts and the bearing sample was far more and the squeeze film
could not generate. In the last test method, the micro bearing sample was not constrained
0 1 2 3 4 5 6?10
?8
?6
?4
?2
0
2
4
6
8
t  (s)
h  (
 ?m)
Base oscillation started
Base oscillation stopped
Figure 5.14: Experimental results of micro-scale bearing
by any means. The starting and stopping of base oscillation was performed in approxi-
mately 2 seconds. Most of the tests using this method failed because the sample tended
to translate in horizontal direction hampering the readings from LBDMS. One good result
(See Fig. 5.14) was observed for the bearing configuration d (Fig. 5.13) having diameter of
70
100 ?m and height of 80 ?m. Here, once the base oscillation was started, the squeeze film
achieved a steady state, and after stopping the base oscillation, it returned to the initial
state (Fig. 5.14). For this test case, the frequency of the base oscillation, f, was 2000 Hz
0 2 4 6 8 10 12 14 16 18 200.975
0.98
0.985
0.99
0.995
1
1.005
1.01
1.015
1.02
1.025
T (rad)
P n
Figure 5.15: Dimensionless mean pressure, Pn, as a function of normalized time, T, for the
ultra-thin squeeze film bearing
and the amplitude, Z0, was 0.63 ?m. The mean squeeze film thickness, hm, obtained was
2 ?m. As the weight of the sample is 0.01079 N, the load on one post of bearing surface is
71
47.955 ?N. If the amplitude of oscillation of the squeeze film height is assumed to be 0.63
?m, then the excursion ratio is 0.63/2 (i.e. 0.315). In order to compare the analytical and
experimental results, the squeeze film height is assumed as a known function of time (Eq.
(5.2)) in analytical work.
h = hm?(1?epsilon1?sin(T)) (5.2)
The Reynolds equation (Eq. (4.22)) is solved for pressure with the assumptions as made in
Eq. (5.2) (See Appendix E). The pressure trace obtained is shown in Fig. 5.15. The pressure
trace for every cycle is the same. However, the pressure profile is not sinusoidal because
of the nonlinear pressure-squeeze film thickness relationship [6] as seen in the governing
Reynolds equation (Eq. (4.25)). The mean pressure over one cycle results into a load
carrying capacity of 40.16 ?N. The % deviation of the theoretical result (from the pressure
trace) and the actual load (from weighing the sample) is 16.25. This comparison shows
that the experimental results and theoretical results for these micro-scale bearing samples
are in good agreement with each other. In future, more comparisons of experimental and
simulation results are needed to validate the numerical model.
72
Chapter 6
Summary
The tribology of macro-scale systems such as power plants, automobiles, air-craft en-
gines etc has been greatly studied in the last century. As the miniaturization of mechanical
systems is a need of future technology, it is important to address the issues related to
friction, wear and lubrication for such micro-scale mechanical systems. The research under-
taken starts with the experimental and analytical study of macro-scale squeeze film bearings
which is later extended to study micro-scale squeeze film bearings.
The first part of this research is to extensively study the squeeze film bearings where
the squeeze films are characterized as thin films (the film thickness is in the range of 9
to 23 ?m). A coupled dynamic model with asperity contact effects is developed to study
compressible dynamic squeeze films between disk shaped surfaces, in which, one disk is
excited by a sinusoidal displacement. The model presented is general and can be used to
investigate dynamic squeeze films with input parameters, frequency (f) and amplitude of
vibration (Z0), mass (m), area of contact (A) and the surface properties. From the results
of the numerical simulations, a comparison with Fi and Fn from Salbu?s Eq. is made.
Based on these comparisons a new semi-analytical equation is developed to predict the load
carrying capacity for 0.4???10 using an exponential curve to fit the simulation results.
Experimental results disagree quantitatively because of the inability to perfectly model the
73
experimental system. Qualitatively, the experimental results are in fairly good agreement
with the numerical results.
The second part of this research is to study the micro-scale squeeze film bearings, where
the squeeze films are characterized as ultra-thin (using Knudsen number, the film thickness
should be less than 8 ?m). A sample of 3 cm X 3 cm patterned with an array of micro
bearing areas having 100 ?m diameter is used for an experimental purpose. A single unit
cell of this pattern is shown in Fig 5.13. A squeeze film thickness of approximately 2 ?m
is measured experimentally when the micro bearing was operated at a frequency of 2000
Hz and an amplitude of 0.63 ?m. To compare with experimental results, the squeeze film
thickness is assumed as a known function of time and the discretized Reynolds equation is
solved. The deviation of load carrying capacity from the simulation to the actual load is
found to be 16 %. These results are in good agreement with each other, although more
extensive work is needed to confirm the results.
In conclusion, numerical simulations and experiments have been performed to investi-
gate the compressible squeeze film bearings. Both experimental and analytical results have
shown that the squeeze film bearings have a good potential for lubrication in macro and
micro-scale mechanical systems. These results are important in an age where considerable
effort has been made to develop gas bearings for MEMS in research laboratories around the
world.
74
Bibliography
[1] Jost, H.P. (2006), Tribology: How a word was coined 40 years ago, in Tribology and
Lubrication Technology. pp 24-28.
[2] Kajdas, C., Harvey, S.S.K., and Wilusz, E. (1990), Encyclopedia of Tribology, Elsevier
Science Publishing Company Inc., Amsterdam, pp 478.
[3] Hamrock, B.J. (1991), ?Fundamentals of Fluid Film Lubrication,? NASA Office of
Management, Scientific and Technical Information Program: [Supt. of Docs., U.S.
G.P.O.,distributor], Washington, D.C., pp 2-9.
[4] Parkins, D.W. and Stanley, W.T. (1982) ?Characteristics of an Oil Squeeze Film,? J
LUBR TECHNOL TRANS ASME, V 104, 4,pp 497-503.
[5] Minikes, A. and Bucher, I. (2003), ?Noncontacting lateral transportation using gas
squeeze film generated by flexural traveling waves - Numerical analysis? J of the Acous-
tical Society of America, v 113, 5,pp 2464-2473.
[6] Wiesendanger, M. (2001), ?Squeeze film air bearings using piezoelectric bending ele-
ments,? Ecole Polytechnique Fdrale de Lausanne, Lausanne: EPFL.
[7] Khonsari, M.M. and Booser, E.R. (2001) ?Applied tribology : bearing design and
lubrication ,? John Wiley,New York, pp 258.
[8] Langlois, W.E. (1962), ?Isothermal Squeeze Films,? Quarterly of Applied Mathematics,
20,2, pp 131-150.
[9] Salbu, E.O.J. (1964), ?Compressible squeeze films and squeeze bearings,? ASME J.
Basic Eng., 86, pp 355-366.
[10] Minikes, A. and Bucher, I. (2003), ?Coupled dynamics of a squeeze-film levitated mass
and a vibrating piezoelectric disc: Numerical analysis and experimental study,? J.
Sound Vib., 263, 2, pp 241-268.
[11] Miller, B.A. and Green, I. (2001), ?Numerical Formulation for the Dynamic Analysis
of Spiral-Grooved Gas Face Seals,? ASME J. Tribol., 123, 2, pp 395-403.
[12] Pretlove, A.J. (1985), BASIC mechanical vibrations, Butterworth and Co. (Publishers)
Ltd., Kent, England, pp 118.
[13] Griffin, W.S., Richardson, H.H., and Yamanami S. (1966), ?A study of Fluid Squeeze-
Film Damping,? ASME J. Basic Eng., 88, pp 451-456.
75
[14] Blech, J.J. (1983), ?On Isothermal Squeeze Films,? J. Lubr. Technol. Trans. ASME,
105, pp 615-620.
[15] Ausman, J.S. (1967) ?Gas Squeeze Film Stiffness and Damping Torques on a Circular
Disk Oscillating About Its Diameter,? J Lubr Technol Trans ASME, 89, pp 219-221.
[16] Etsion, I. (1980), ?Squeeze Effects in Radial Face Seals,? J Lubr Technol Trans ASME,
102, 2, pp 145-152.
[17] Green, I. and Etsion, I. (1983), ?Fluid Film Dynamic Coefficients in Mechanical Face
Seals,? ASME J. Lubr. Technol., 105, 2, pp 297-302.
[18] Blech, J.J. (1985), ?Annular Compressible Squeeze Films,? ASME J. Tribol., 107, 4,
pp 544-547.
[19] Starr, J.B. (1990), ?Squeeze-film damping in solid-state accelerometers,? Techn. Digest,
Solid State Sensor and Actuator Workshop Hilton Head Island, pp 44-47.
[20] Pan, F., Kubby, J., Peeters, E., Tran, A.T., and Mukherjee, S. (1998), ?Squeeze film
damping effect on the dynamic response of a MEMS torsion mirror,? J. Micromech.
Microeng, 8, 3, pp 20-208.
[21] Bao, M., Yang, H., Yin, H., and Sun, Y. (2002), ?Energy transfer model for squeeze-film
air damping in low vacuum,? J. Micromech. Microeng., 12, 3, pp 341-346.
[22] Bao, M., Yang, H., Sun, Y., and French, P.J. (2003), ?Modified Reynolds equation and
analytical analysis of squeeze-film air damping of perforated structures,? J. Micromech.
Microeng., 13, 6, pp 795-800.
[23] Nayfeh, A.H. and Younis, M.I. (2004), ?A new approach to the modelling and simula-
tion of flexible microstructures under the effect of squeeze film damping,? J. Micromech.
Microeng, 14, 2, pp 170-181.
[24] Pandey, A.K. and Pratap, R. (2004), ?Coupled nonlinear effects of surface roughness
and rarefaction on squeeze film damping in MEMS structures,? J. Micromech. Micro-
eng, 14, 10, pp 1430-1437.
[25] Same, H. and Wenjing, Y. (2004), ?On the squeeze-film damping of micro-resonators
in the free-molecule regime,? J. Micromech. Microeng., 14, 12, pp 1726-1733.
[26] Zhang, C., Xu, G., and Jiang, Q. (2004), ?Characterization of the squeeze film damping
effect on the quality factor of a microbeam resonator,? J. Micromech. Microeng., 14,
10, pp 1302-1306.
[27] Strawhecker, K., Asay, D. B., McKinney, J., Kim, S. H. (2005), ?Reduction of adhesion
and friction of silicon oxide surface in the presence of n-propanol vapor in the gas
phase,? Tribology Letters, 19, 1, pp 17-21.
76
[28] Epstein, A.H. (2004), ?Millimeter-scale, MEMS gas turbine engines,? J. Eng. Gas
Turbines Power, 126, 2, pp 205-226.
[29] Epstein, A.H., et al. (1997), ?Micro-Heat Engines, Gas Turbines, And Rocket En-
gines -The MIT MICROENGINE PROJECT,? American Institute of Aeronautics and
Astronautics, pp 1-12.
[30] Breuer, K. (2001), Lubrication in MEMS, Ed. M. Gad el Hak, CRC Press, pp 1-49.
[31] Wong, C.W., Zhang, X., Jacobson, S.A., and Epstein, A.H. (2004), ?A Self-Acting
Gas Thrust Bearing for High-Speed Microrotors,? J Microelectromech Syst, 13, 2, pp
158-164.
[32] Borionettia, G., Bazzali, A., and Orizio, R. (2004), ?Atomic force microscopy: a pow-
erful tool for surface defect and morphology inspection in semiconductor industry,?
Eur. Phys. J. Appl. Phys., 27, pp 101-106.
[33] Sun, Y., Chan, W.K., and Liu, N. (2003), ?A slip model with molecular dynamics,? J
Micromech Microengineering, 12, 3, pp 316-322.
[34] Bird, G.A. (1994), Molecular Gas Dynamics and the Direct Simulation of Gas Flows,
Oxford: Clarendon.
[35] Hsia, Y.T. and Domoto, G.A. (1983), ?An Experimental Investigation of Molecular
Rarefaction Effects in Gas Lubricated Bearings at Ultra Low Clearances,? J Lubr
Technol Trans ASME, 105, 1, pp 120-130.
[36] Li, W.L. (1999), ?Analytical modelling of ultra-thin gas squeeze film,? Nanotechnology,
10, pp 440-446.
[37] Darling, R.B., C., H., and J., X. (1997), ?Compact analytical modeling of squeeze film
damping with arbitrary venting conditions using a Green?s function approach,? Sens
Actuators A Phys, 70, 1-2, pp 32-41.
[38] Veijola, T. and Turowski, M. (2001), ?Compact damping models for laterally moving
microstructures with gas-rarefaction effects,? J Microelectromech Syst, 10, 2, pp 263-
273.
[39] Constantinescu, V.N. (1969), Gas Lubrication, The American Society of Mechanical
Engineers, New York, USA, pp 621.
[40] Jackson, R.L. and Green, I. (2005), ?A Finite Element Study Of Elasto-Plastic Hemis-
perical Contact Against a Rigid Flat,? ASME J. Tribol., 127, 2, pp 343-354.
[41] Jackson, R.L. and Green, I. (2006), ?A statistical model of elasto-plastic asperity con-
tact between rough surfaces,? Trib. Intl., 39, pp 906-914.
77
[42] Jaeger, R.C. (1998), Introduction to Microelectronic Fabrication, Prentice Hall, Inc.,
pp 316.
[43] Greenwood, J.A. and Williamson, J.B.P. (1966), ?Contact of Nominally Flat Surfaces,?
Proc. R. Soc. London, 295, pp 300-319.
[44] McCool, J.I. (1987), ?Relating Profile Instrument Measurements to the Functional
Performance of Rough Surfaces,? ASME J. Tribol., 109, 2, pp 264-270.
[45] Front, I. (1990), ?The Effects of Closing Force and Surface Roughness on Leakage in
Radial Face Seals,? MS Thesis,Technion, Israel Institute of Technology, Israel.
[46] Horton, B.D. (2004), ?Magnetic Head Flyability on Patterned Media,? MS Thesis, The
George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technol-
ogy, Atlanta, GA.
78
Appendices
79
Appendix A
Contact Force
By fitting equations to finite element results, Jackson and Green [40] and [41] provide
the following equations to predict the elastic perfectly-plastic contact of a sphere and a rigid
flat. P?F is a ratio of contact load to critical contact load, Pc. For 0??????t
P?F = (??)3/2 (A.1)
where, ?? is the ratio of penetration or indentation depth between spherical asperities
(?) to the critical interference (?c) and ??t is the value that defines the effective transition
from elastic to plastic behavior of ??. For Jackson and Green [40], ?? is 1.9.
For ??t???
P?F =
bracketleftbigg
exp
parenleftbigg
?14 (??)5/12
parenrightbiggbracketrightbigg
(??)3/2 +4HGCS
y
bracketleftbigg
exp
parenleftbigg
? 125 (??)5/9
parenrightbiggbracketrightbigg
(??)?AE (A.2)
where,
HG
Sy = 2.84
?
?1?exp
parenleftBigg
?0.82
parenleftBigg
piCey
2
???parenleftbigg ??
??t
parenrightbiggB/2parenrightBiggparenrightBigg0.7?
? (A.3)
In Eq. (A.3), HG is the limiting average contact pressure, Sy is yield strength and ey is
uniaxial yield strain (ratio of yield strength to equivalent elastic modulus).
80
The critical interference to cause initial yielding, ?c, is derived independently of the
hardness, to be
?c =
parenleftbiggpi?C?S
y
2E?
parenrightbigg2
R (A.4)
where, R is the radius of the hemispherical asperity and E? is the equivalent elastic modulus.
B and C are functions of the material properties given as
B = 0.14?exp(23ey) (A.5)
C = 1.295?exp(0.736?) (A.6)
This model then assumes that the individual asperity contact between rough surfaces
can be approximated by hemispherical contact with a rigid flat. Then, statistical relation-
ships from Greenwood and Williamson [43] are used to model an entire surface of asperities
with a range of heights described by a Gaussian distribution, G(z). These statistical equa-
tions are given as
Fcont = ?An?
integraldisplay
d
PF(z ?d)G(z)dz (A.7)
where, the average asperity radius of curvature, R, and the asperity surface density, ?, are
needed to model asperity contact and are obtained from a profilometer produced surface
profile using the methods outlined in McCool [44]. The distance between the surfaces can
81
be described in two ways: (1) the distance between the mean of the surface heights, h, and
(2) the distance between the mean of the surface asperities or peaks, d. These values of h
and d are related by
h = d+ys (A.8)
The value of ys is derived by Front [45] and given as
ys = 0.045944?R (A.9)
where, ? is the area density of the asperities. Eq. (A.7) is then numerically integrated to
predict Fcont as a function of h. The surface profile of one of the rough bearing surfaces is
used in the work shown in Fig. A.1. Dimensionless contact load for the surface is plotted
against dimensionless mean separation (See Fig. A.2). An exponential fit to the data in Fig.
A.2 is obtained. The resulting fit, Eq. (A.10) is then included in the numerical simulation
to predict the contact force.
Fcont =
parenleftbigg
?9.98?exp
parenleftbigg
?3.621? h?
parenrightbigg
+9.824?exp
parenleftbigg
?3.595? h?
parenrightbiggparenrightbigg
?AE? (A.10)
Here, the goodness of fit is
Sum of squares due to error, SSE= 7.586e-008
The ratio of the sum of squares of the regression, R-square= 0.9988
82
0 50 100 150 200 250 300 350 400 450 500?0.3
?0.25
?0.2
?0.15
?0.1
?0.05
0
0.05
0.1
0.15
0.2
Sampling length (?m)
Surface asperity height (
?m)
Figure A.1: Surface profile of a rough bearing surface
83
1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 60
0.5
1
1.5
2
2.5
3 x 10
?3
Dimensionless mean separation, h/Rq
Dimensionless load, 
F cont
/(A
nE
? )
Figure A.2: Dimensionless contact load as a function of dimensionless mean surface sepa-
ration
84
Adjusted R-square= 0.9988
Root mean squared error, RMSE= 1.967e-005
85
Appendix B
C-program to solve the coupled dynamics for thin squeeze films
include <stdio.h>
include <process.h>
include <conio.h>
include <dos.h>
include <math.h>
define MPI 3.1415926535897932385E0
FILE *f1;
FILE *f2;
FILE *f3;
FILE *f4;
double p[17][2]; // 2-D array for pressure
double meanpressure(double h1,double sigma, double p[][2],double R[],double Rst,double
tst,double h[]);
main()
/*********************All Variable Declaration**************************/
double t; // time
double f; // frquency in Hz
double w; // frequency in radians
double Rst; // Step-size for R
double h[2];
double hdot[2]; // Derivative of Height
double Mult=0; //Dummy variable
double tst; // Time step
double Ro=0; // Radius of circular area (i.e. Area of contact)
double visc; // Vissinity of the oil at T=293
double R[17]; //1-D array for Radius
double Rinv[17];
double pmean=0; //Mean pressure
double pa=0; //atmospheric pressure
double g=9.81; //accn due to gravity
double E=0;
double F=0; //Dummy Variable
double Aml=2.0*pow(10,-6); // Amlitude of the Shaker
double m=7*pow(10,-3);; // Mass of the levitating plate
double Ar; // Area of Contact
86
double sig=0.17*pow(10,-6);
double CntFr=0;
double sigma=0;
double YM=0;
double Dum1=0;
double Dum2=0;
double Dum3=0;
double Dum4=0;
double tstDum2=0;
/********Variables Used for Runge-Kutta********/
double K11=0;
double K12=0;
double K21=0;
double K22=0;
double K31=0;
double K32=0;
double K41=0;
double K42=0;
/********Variables Used for Runge-Kutta********/
double K1=0;
double K2=0;
double K=0;
double D1=1;
double D2=0;
double D3=0;
double ho=0;
double h1=0;
double hdt=0;
double htst;
double qtst;
double C1=0;
double C2=0;
double C3=0;
double C4=0;
double C5=0;
double x=0;
/****************Assignment of Values********************/
f=1500;
w=2*MPI*f;
visc=1.76/100000;
ho=30*5.554*sig;
Ro=9.0*pow(10,-3);
87
pa=101325;
sigma=(12*visc*w*Ro*Ro)/(pa*ho*ho);
Ar=MPI*Ro*Ro;
YM=57.67*pow(10,9); /************Dummy Variables ********************/
Dum1=YM*Ar;
Dum2=1/(w*w*ho);
Dum3=w*w*Aml;
Dum4=(pa/m)*Ar;
/**********Step Sizes *************/
Rst=1/16.0;
tst=0.005;
tstDum2=tst*Dum2;
htst=tst/2;
qtst=tst/4; /************ Initial Conditions *****************/
t=0;
h[0]=1;
hdot[0]=0;
pmean=1;
int k=0;
int j=0;
/************ Assignment of Valuesof R @ different node points*******/
R[0]=0;
for(k=0;k<16;k++)
/***************Mean Pressure ***********************/
{
R[k+1]=R[k]+Rst;
}
/***************Assignment of Valuesof R @ Ends Here ***********/
/******Assignment for pressure at time t=0 at all the radius is 1*******/
for(k=0;k<=16;k++)
{
p[k][0]=1;
} /*********Assignment for pressure for initial guess **************/
for(k=0;k<=16;k++)
{
p[k][1]=1;
}
/*******************Assignment of Boundary Conditions ***********/
p[16][1]=1;
// i is for radius //
// and j is for time //
88
x=(h[0]*ho)/sig;
CntFr=(-9.98*exp(-3.621*x)+9.824*exp(-3.595*x))*Dum1;
printf(?
f1 = fopen (?Height.txt?, ?wt?); /****3****/
f2 = fopen (?Hdot?, ?wt?); /****3****/
f3 = fopen (?MeanPr?, ?wt?); /****3****/
f4 = fopen (?Time?, ?wt?); /****3****/
for (int j=0;j<200000000;j++)
{
/******************* RK1 *************************/
if(h[0]<2)
{
x=(h[0]*ho)/sig;
CntFr=(-9.98*exp(-3.621*x)+9.824*exp(-3.595*x))*Dum1;
}
else
{ CntFr=0;
}
K11=tst*hdot[0];
K12=tstDum2*(Dum3*sin(t)-0.1*Dum3*sin(7*t)+
(CntFr/m)+(pmean-1)*Dum4-9.810007193613373);
/******************* RK1 *************************/
/**************H and HDOT After RK1**********************/
h1=h[0]+0.5*K11;
hdt=hdot[0]+0.5*K12;
Contact Force
if(h1<2)
{
x=(h1*ho)/sig;
CntFr=(-9.98*exp(-3.621*x)+9.824*exp(-3.595*x))*Dum1; }
else
{
CntFr=0;
}
/******************Mean Pressure ***********************/
pmean=meanpressure(h1,sigma,p,R,Rst,htst,h);
/*******************RK2*********************************/
t=t+htst;
K21=tst*(hdot[0]+0.5*K12);
K22=tstDum2*(Dum3*sin(t)-0.1*Dum3*sin(7*t)+
(CntFr/m)+(pmean-1)*Dum4- 9.810007193613373);
h1=h[0]+0.5*K21;
89
hdt=hdot[0]+0.5*K22;
if(h1<2)
{
x=(h1*ho)/sig;
CntFr=(-9.98*exp(-3.621*x)+9.824*exp(-3.595*x))*Dum1;
}
else
{
CntFr=0;
}
pmean=meanpressure(h1,sigma,p,R,Rst,htst,h);
/*******************RK3*************************/
K31=tst*(hdot[0]+0.5*K22);
K32=tstDum2*(Dum3*sin(t)-0.1*Dum3*sin(7*t)+
(CntFr/m)+(pmean-1)*Dum4- 9.810007193613373);
/*******************RK3*************************/
h1=h[0]+K31;
hdt=hdot[0]+K32;
if(h1<2)
{
x=(h1*ho)/sig;
CntFr=(-9.98*exp(-3.621*x)+9.824*exp(-3.595*x))*Dum1;
}
else
{
CntFr=0;
}
pmean=meanpressure(h1,sigma,p,R,Rst,tst,h);
/************RK4******************/
t=t+htst;
K41=tst*(hdot[0]+K32);
K42=tstDum2*(Dum3*sin(t)-0.1*Dum3*sin(7*t)+
(CntFr/m)+(pmean-1)*Dum4- 9.810007193613373);
h[1]=h[0]+0.16666667*(K11+2*K21+2*K31+K41);
hdot[1]=hdot[0]+0.16666667*(K12+2*K22+2*K32+K42);
/**************Mean Pressure *******************/
h1=h[1];
pmean=meanpressure(h1,sigma,p,R,Rst,tst,h);
h[0]=h[1];
hdot[0]=hdot[1];
for(int i=0;i<16;i++)
{
90
p[i][0]=p[i][1];
}
if(j
{
fprintf(f1,?
fprintf(f2,?
fprintf(f3,?
fprintf(f4,?
}
}
fclose(f1);
fclose(f2);
fclose(f3);
fclose(f4);
}
fprintf (f1, ?
doublemeanpressure(doubleh1,doublesigma, doublep[][2],doubleR[],doubleRst,double
tst,double h[])
{
double h3sig,A,B,C,D,R1,R2,P1,P2,pMean;
double error[17];
double M[17];
double err=0;
int s;
h3sig=(pow(h1,3))/(sigma);
B=h1/tst;
for(s=0;s<=16;s++)
{
M[s]=p[s][1];
}
/****************Assignment of Dummy Variable ******************/
do
{
for(s=15;s>0;s?)
{ R1=(R[s+1]+R[s])/2;
R2=(R[s-1]+R[s])/2;
A=(h3sig*(R1+R2))/(2*Rst*R[s]);
C=((-p[s][0]*h[0])/tst)-((h3sig)/(2*Rst*R[s]))*
(R1*pow(p[s+1][1],2)+R2*pow(p[s-1][1],2));
p[s][1]=(-B+sqrt(B*B-4*A*C))/(2*A);
if (s==1)
{
91
p[0][1]=p[1][1];
}
}
for(s=15;s>0;s?)
{
D=p[s][1];
p[s][1]=1.1*p[s][1]-0.1*M[s];
error[s]=fabs((p[s][1]-M[s])/p[s][1]);
M[s]=D;
}
p[0][1]=p[1][1];
err=0;
for(s=15;s>0;s?)
{
err=err+error[s];
}
err=err/(15); /// Average error
}
while(err > 0.000001);
P1=0;
P2=0;
for(s=1;s<=15;s=s+2)
{
P1=P1+p[s][1];
}
for(s=2;s<=14;s=s+2)
{
P2=P2+p[s][1];
}
pMean=(p[0][1]+4*P1+2*P2+p[16][1])/48;
return pMean;
}
92
Appendix C
Tables of Experimental and Simulation Results
C.1 Configuration 1
Frequency Amplitude Simulation hm Experimental hm
(Hz) (?m) (?m) (?m)
800 6.1792 143.2729 0.9447
800 6.669 147.0127 1.2371
800 8.9486 162.775 2.5657
1050 2.4304 6.0489 3.1891
1050 3.617 106.1366 12.0375
1050 4.936 123.2343 21.4303
1050 5.9116 134.5224 23.4826
1300 2.6374 76.4061 6.5227
1300 2.9202 81.0329 7.6587
1300 4.598 107.2562 16.8267
93
C.2 Configuration 2
Frequency 800 Hz
Amplitude (?m) 7.4 7.6 7.8 8.0 8.2
Test 1 13.4707 14.269 15.6517 16.274 16.1508
Test 2 13.0979 14.6449 15.0356 15.8266 16.0957
Test 3 13.2905 13.9359 15.6428 15.3388 16.1507
Test 4 13.0951 13.6214 15.0053 15.6878 16.0518
Test 5 13.5354 14.617 14.988 15.7182 16.0818
Test 6 13.9357 15.3383 15.6701 15.6848 16.3287
Test 7 13.445 14.469 15.4086 15.6693 16.0388
Test 8 13.1422 14.3621 15.6795 15.7049 16.5886
Test 9 13.4914 14.8486 15.1449 15.4356 15.8794
Test 10 13.3572 14.8437 15.2333 15.5237 15.508
Average (?m) 13.38611 14.49499 15.34598 15.68637 16.08743
Simulation (?m) DNC1 DNC DNC 57.3054 60.1417
Frequency 900 Hz
Amplitude (?m) 5.8 6 6.2 6.4 6.6
Test 1 11.336 11.7322 12.2815 13.0025 13.9145
Test 2 10.6783 11.7956 11.8943 12.5222 13.2427
Test 3 11.01 12.0665 12.23 12.8047 13.263
Test 4 10.8241 11.8172 11.999 12.7467 13.1827
Test 5 11.3899 11.2529 11.8836 13.0562 13.414
Test 6 11.1745 11.6384 12.3865 12.9316 13.8134
Test 7 11.1602 11.7326 12.3004 12.8417 13.5801
Test 8 10.6472 11.1909 12.4745 12.6746 13.4392
Test 9 11.1084 11.1556 12.2482 12.6892 13.2696
Test 10 10.99 11.2763 12.5702 12.896 13.8474
Average (?m) 11.03186 11.56582 12.22682 12.81654 13.49666
Simulation (?m) DNC 55.5434 58.0158 60.4626 62.4467
1 DNC -Did Not Converge
94
Frequency 1000 Hz
Amplitude (?m) 4.6 4.8 5.0 5.2 5.4
Test 1 11.0867 11.262 11.9811 12.5825 12.6396
Test 2 11.1729 11.1629 12.1185 12.72 12.6326
Test 3 11.264 12.0024 12.2013 12.4823 13.5024
Test 4 10.945 11.9393 12.2688 12.1882 13.1681
Test 5 11.0825 11.485 12.1646 13.0589 12.8527
Test 6 10.958 11.7505 11.9228 12.9179 12.8413
Test 7 10.6522 12.3604 12.3521 12.4905 14.1083
Test 8 10.7868 11.7528 11.8116 12.1864 12.8677
Test 9 11.2518 11.6818 11.8853 12.2711 12.6816
Test 10 10.9668 11.3539 12.1838 11.9896 13.3451
Average (?m) 11.01667 11.6751 12.08899 12.48874 13.06394
Simulation (?m) DNC 56.4912 59.1097 61.2036 62.9875
Frequency 1500 Hz
Amplitude (?m) 1.8 2.0 2.2 2.4 2.6
Test 1 9.9339 10.9988 11.0939 11.6938 12.9752
Test 2 10.0439 10.1946 10.8504 12.5324 12.111
Test 3 8.9149 10.3972 10.8781 11.6999 12.7567
Test 4 8.8607 10.5234 10.7668 11.1444 12.5647
Test 5 8.9149 10.9258 11.3802 11.6934 13.2707
Test 6 10.2578 10.9915 11.0566 11.328 12.8714
Test 7 9.9627 10.8179 11.1325 11.9491 12.9994
Test 8 10.0065 10.7677 11.1432 11.8122 12.9973
Test 9 10.0962 10.6744 11.2106 11.8719 11.9577
Test 10 9.9072 10.7097 11.0775 11.5196 12.8166
Average (?m) 9.68987 10.7001 11.05898 11.72447 12.73207
Simulation (?m) 53.5354 56.0012 58.1313 60.0423 61.7946
95
C.3 Configuration 3
Frequency 800 Hz
Amplitude (?m) 6.4 6.6 6.8 7.0 7.2
Test 1 21.4024 21.4089 22.1387 22.7954 22.7595
Test 2 20.9746 21.5593 21.8336 22.7572 23.0541
Test 3 21.2554 21.6139 21.7966 22.5867 22.9858
Test 4 20.9251 21.739 21.7729 22.6405 23.2238
Test 5 20.8396 21.7471 22.0179 22.6409 23.1163
Test 6 21.0984 21.4937 21.9323 22.8155 22.9272
Test 7 21.0307 21.4902 21.8789 22.5575 23.0514
Test 8 20.8029 21.438 22.096 22.4713 23.0452
Test 9 20.8746 21.4432 21.7509 22.4318 23.0638
Test 10 20.8737 21.5169 22.0775 22.5943 23.1871
Average (?m) 21.00774 21.54502 21.92953 22.62911 23.04142
Simulation (?m) 65.592 66.6991 67.7418 68.7306 69.6732
Frequency 900 Hz
Amplitude (?m) 5.2 5.4 5.6 5.8 6.0
Test 1 18.7737 20.1916 20.8962 20.8628 21.7042
Test 2 18.8144 20.0742 21.0646 20.8695 21.5064
Test 3 19.0361 19.9377 20.9214 21.19 21.1382
Test 4 19.1751 19.528 20.9238 21.3347 21.2448
Test 5 19.0633 20.0621 20.9489 21.0129 21.7445
Test 6 18.9867 20.0116 20.8519 21.1412 21.3917
Test 7 19.1702 19.5744 20.4802 21.506 21.196
Test 8 19.0511 19.1967 21.1435 21.0589 21.05
Test 9 19.3148 19.4528 20.7809 20.8976 21.1931
Test 10 19.2281 19.8355 20.9035 21.222 21.3466
Average (?m) 19.06135 19.78646 20.89149 21.10956 21.35155
Simulation (?m) 65.4042 66.5021 67.5448 68.5408 69.4961
96
Frequency 1000 Hz
Amplitude (?m) 4 4.2 4.4 4.6 4.8
Test 1 20.8198 22.0756 22.7689 23.0927 23.2616
Test 2 20.1229 22.0891 22.4801 22.8036 23.2182
Test 3 20.3533 21.9974 22.4754 22.8771 23.4326
Test 4 20.5557 22.0531 22.5285 22.5587 22.9707
Test 5 20.887 21.7522 22.322 22.7171 23.0604
Test 6 20.6034 22.1179 22.0915 23.1343 23.7415
Test 7 20.6502 22.0388 22.6978 22.7592 23.4185
Test 8 20.498 21.9899 22.8898 22.4465 23.3308
Test 9 20.6819 22.1545 22.6626 22.5107 23.8755
Test 10 20.8721 21.6264 22.6115 22.9328 23.7882
Average (?m) 20.60443 21.98949 22.55281 22.78327 23.4098
Simulation (?m) 62.8183 64.0531 65.2203 66.331 67.3937
Frequency 1500 Hz
Amplitude (?m) 1.8 2.0 2.2 2.4 2.6
Test 1 16.2606 16.4348 16.944 17.1475 17.508
Test 2 16.4536 16.7723 16.9895 17.1686 17.2069
Test 3 16.2021 16.406 16.9542 17.2596 17.7354
Test 4 16.3595 16.9441 17.0921 17.2073 17.3995
Test 5 16.4152 16.5749 16.9935 17.4292 17.4173
Test 6 16.145 16.484 16.9261 17.1501 17.5845
Test 7 16.301 16.3783 16.9782 17.3174 17.7091
Test 8 16.2998 16.3724 16.9907 17.02 17.7404
Test 9 16.114 16.678 17.0715 17.2938 17.7186
Test 10 16.1834 16.5505 16.9701 17.2764 17.7282
Average (?m) 16.27342 16.55953 16.99099 17.22699 17.57479
Simulation (?m) 54.5747 56.4424 58.2134 59.8992 61.5094
97
Appendix D
Calibration of Capacitance Sensor
Equation for calibrating physical capacitance, C to the mean squeeze film thickness,
hm in ?m is
hm = A?epsilon10 ?k1.853?C1.706 +0.0784 (D.1)
where, physical capacitance, C in F is given by
C = ?Rst ?(?1+
1
2.511?(VRexp)4 +6.011?(VRexp)3 ?4.72?(VRexp)2 +2.29?(VRexp)?0.1417) (D.2)
and
VRexp= experimental voltage ratio as recorded from LabView
?= frequency of AC voltage, Vin (2?pi?2000 rad)
Rst= resistor used in the electric circuit (33 k?)
A= area of contact (m2)
epsilon10= permittivity of free space (8.85E-6 F??m?1)
k= dielectric constant of air (1)
98
Calibration Procedure
1. Read the experimental voltage, Vo/p using LabView.
2. Calculate the experimental voltage ratio.
VRexp = Rms Vo/pRms V
in
(D.3)
3. Experimental voltage ratio is calibrated in terms of the physical voltage ratio. Stan-
dard capacitors are utilized for this calibration. For each standard capacitor, the
experimental voltage ratio is calculated. Then, using these standard capacitor values
and Eq. (D.5), physical voltage ratios are calculated.
As,
VRphys = RR+ 1
j?C
(D.4)
Thus,
VRphys2 = ?
2C2R2
1+?2C2R2 +2?RC (D.5)
Using a curve fitting technique, a polyfit is obtained for VRphys in terms of VRexp
and is given by Eq. (D.6). Thus, the experimental voltage ratio is converted to the
99
physical voltage ratio using the 4th order polynomial (Eq.(D.6)).
VRphys = 2.511?(VRexp)4 +6.011?(VRexp)3 ?4.72?(VRexp)2 +
2.29?(VRexp)?0.1417 (D.6)
4. Capacitance due to parallel plate can be calculated using the physical voltage ratio.
A quadratic equation in terms of C is written as
bracketleftBiggparenleftBigg
1? 1VR
phys2
parenrightBigg
?2Rst2
bracketrightBigg
C2 +[2?Rst]C +1 = 0 (D.7)
Here,
?=2pi?2000 rad
Rst=33 k?
and VRphys is calculated as per step 3. Solving Eq. (D.7) and taking positive root,
C is calculated as
C = ?Rst ?(?1+ 1VR
phys
) (D.8)
5. Mylar shims of thickness 38.1 ?m are used to calibrate the squeeze film thickness to
the capacitance. The two bearing surfaces are separated using mylar shims so that
99.7 percent volume between the two plates is air and 0.3 percent volume is mylar.
Thus, two surfaces form two capacitors in parallel, one due to air and the other due
100
to mylar. Capacitance due to each dielectric medium is calculated using
C = Aepsilon10kd (D.9)
Here, k for air =1 and k for mylar = 3.2. The total capacitance is the sum of
capacitance due to air and mylar. This theoretical capacitance is then calibrated in
terms of experimental capacitance (Eq. D.8). This results in a power function which
represents the above calibration.
Ccal = 1.853?C1.706 +0.0784 (D.10)
Here, the goodness of fit is
Variance Reduction= 99.99
S/(N - P) : 0.00002754
RMS (Y - Ycalc) : 0.00262
6. The mean squeeze film height is calculated using Eq. (D.11) as
hm = Aepsilon10kC
cal
(D.11)
101
Appendix E
Computer program to solve the dynamics for ultra-thin squeeze films
import java.io.*;
import java.math.*;
class Reynoldsequation
{
public static void main(String args[])
{
/*************All Variable Declaration****************/
double Rst; //Step-size for R
double[][] p=new double[17][2]; // 2-D array for pressure
double[] h=new double[2];
double tst; // Time step
double[] R=new double[17]; //1-D array for Radius
double[] M=new double[17];
double[] Rinv=new double[17];
double[] error=new double[17];
double err=0; //Mean pressure
double A=0; //Dummy Variable
double B=0; //Dummy Variable
double C=0; //Dummy Variable
double D=0;
double P1=0; //Dummy Variable
double P2=0; //Dummy Variable
double pmean=1;
double sig=0;
double sig1=0;
double R1=0;
double R2=0;
double K1=0;
double K2=0;
double h3sig;
double visc=1.8*Math.pow(10,-5);
double visc1=0;
double t;
102
double N=0;
double Exc=0.315;
int n=0;
double KN=0;
double hm=2*Math.pow(10,-6);
double ME=1.1467*Math.pow(10,-3);
/*************Assignment of Values*********************/
h[0]=1;
sig=929.68;
/**********Step Sizes *************/
tst=0.0001;
/************ Initial Conditions *****************/
t=0;
err=1;
/**************Assignment of Valuesof R @ different node points*********/
R[0]=0;
N=16;
n=16;
Rst=1/N;
for(int k=0;k<N;k++)
{
R[k+1]=R[k]+Rst;
}
/**************Assignment of Valuesof R @ Ends Here ***********/
/*******Assignment for pressure at time t=0 at all the radius is 1*********/
for(int k=0;k<=N;k++)
{
p[k][0]=1;
}
/*********Assignment for pressure for initial guess **************/
for(int k=0;k<=N;k++)
{
p[k][1]=1;
}
/*******************Assignment of Boundary Conditions **********/
p[n][1]=1;
103
// i is for radius //
// and j is for time //
/******* Finite Difference Algorithm Starts Here************/
try
{
BufferedWriter out6 = new BufferedWriter (new FileWriter (?Radial.txt?));
BufferedWriter out7 = new BufferedWriter (new FileWriter (?Error.txt?));
BufferedWriter out1 = new BufferedWriter (new FileWriter (?sigma1000R01.txt?));
BufferedWriter out2 = new BufferedWriter (new FileWriter (?sigma1000R051.txt?));
BufferedWriter out3 = new BufferedWriter (new FileWriter (?sigma1000R781.txt?));
BufferedWriter out4 = new BufferedWriter (new FileWriter (?sigma1000R15161.txt?));
BufferedWriter out5 = new BufferedWriter (new FileWriter (?Pressure.txt?));
for (int j=0;j<1600000;j++)
{
/******************* RK1 *************************/
t=t+tst;
h[1]=1-Exc*Math.sin(t);
KN=ME/(h[1]*hm*pmean*101325);
visc1=visc/(1+2*KN+0.2*(Math.pow(KN,0.780))*Math.exp(-KN/10));
sig1=visc1*sig;
h3sig=(Math.pow(h[1],3))/(sig1);
B=h[1]/tst;
/*************Assignment of Dummy Variable ****************/
for(int i=0;i<n;i++)
{
M[i]=p[i][0];
}
p[n][1]=1;
/**************Assignment of Dummy Variable*****************/
do
{
for(int i=15;i>0;i?)
{
R1=(R[i+1]+R[i])/2;
R2=(R[i-1]+R[i])/2;
A=(h3sig*(R1+R2))/(2*Rst*R[i]);
C=((-p[i][0]*h[0])/tst)-((h3sig)/(2*Rst*R[i]))*
(R1*Math.pow(p[i+1][1],2)+R2*Math.pow(p[i-1][1],2));
104
p[i][1]=(-B+Math.sqrt(B*B-4*A*C))/(2*A);
if (i==1)
{
p[0][1]=p[1][1];
}
}
for(int i=15;i>0;i?)
{
D=p[i][1];
p[i][1]=1.1*p[i][1]-0.1*M[i];
error[i]=Math.abs((p[i][1]-M[i])/p[i][1]);
M[i]=D;
}
p[0][1]=p[1][1];
err=0;
for(int i=15;i>0;i?)
{
err=err+error[i];
}
err=err/(n-1); /// Average error
}
while(err > 0.0000001);
for(int i=0;i<N;i++)
{
p[i][0]=p[i][1];
}
if(j%200000==0)
{
for(int i=0;i<=N;i++)
{
out6.write(Double.toString(p[i][1]));
out6.newLine();
}
}
P1=0;
P2=0;
for(int i=1;i<=15;i=i+2)
{
105
P1=P1+p[i][1];
}
for(int i=2;i<=14;i=i+2)
{
P2=P2+p[i][1];
}
pmean=(p[0][1]+4*P1+2*P2+p[16][1])/48;
if(j%52==0)
{
out5.write(Double.toString(pmean));
out5.newLine();
}
h[0]=h[1];
}
/*******Integration of Pressure using Simpson?s Rule Ends Here *********/
out1.close();
out2.close();
out3.close();
out4.close();
out5.close();
out6.close();
out7.close();
}
catch(IOException E2)}}
/******** Finite Difference Algorithm Ends Here************/
} }
106