Copyright
by
Zhiping Yang
2009
The Dissertation Committee for Zhiping Yang
certifies that this is the approved version of the following dissertation:
Fracture, Friction and Granular Simulation
Committee:
Michael P. Marder, Supervisor
Harry L. Swinney
William D. McCormick
Jack B. Swift
K. Ravi-Chandar
Fracture, Friction and Granular Simulation
by
Zhiping Yang, B.Sc.(Honors)
Dissertation
Presented to the Faculty of the Graduate School of
The University of Texas at Austin
in Partial Fulfillment
of the Requirements
for the Degree of
Doctor of Philosophy
The University of Texas at Austin
May 2009
Dedicated to my wife Hongwei and son Steven.
Acknowledgments
I would like to say thanks to many people for their continuous support and encouragement.
First and foremost, I want to thank my lovely wife Hongwei Shi for taking great care
of me and also bringing me an adorable son Steven. Her patience and endurance contribute
to the making of this thesis significantly.
I also want to thank Dr. Marder for being my supervisor and for his continuous
support both academically and financially. His ever smiling face and humorous comments
encourage me to take on any challenges ahead of me and have fun doing it.
I'd like to thank Dr. Swinney for overseeing Center for NonLinear Dynamics (CNLD).
This is a great research center. I'm honored to be part of it. I have learned a lot from those
weekly seminars and group meetings. What maybe more important is the spirit of coopera-
tion, presentation and inter-change of ideas. It certainly has widened my eyes and deepened
my understanding of science.
I'd like to thank Dr. Robert Deegan for introducing the low temperature fracture
experiment to me. Though his design didn't work out as desired, but I learned many valuable
experimental skills from him.
I'd like to thank Dr. Hepeng for bringing me to the numerical study of granular
systems. His experience and ideas were the fuels that kept me driven.
I'd like to thank Dr. McCormick, Dr. Swift and Dr. Ravi-Chandar for serving as
members of my thesis committee. Their careful inspection and valuable comments made
this thesis much more valuable.
I'd like to thank all other members of CNLD for many of their creative ideas and
suggestions to my research.
Last but not least, I'd like to thank all the technical personnel on the 3rd floor
iv
of RLM. Especially to Allen for machine grounding my steel frames, to Jack for providing
instructions of how to operate machines in the student shop and to Jesse for doing Molecular
Beam Epitaxy (MBE) on all the silicon samples that I broke.
Zhiping Yang
The University of Texas at Austin
May 2009
v
Fracture, Friction and Granular Simulation
Publication No.
Zhiping Yang, Ph.D.
The University of Texas at Austin, 2009
Supervisor: Michael P. Marder
This thesis contains three separate yet closely related topics: fracture, friction and
simulation of mechanical response of a confined granular medium. The first two are experi-
mental investigations and the last one is a numerical study.
In the fracture part, I will describe how to break a piece of silicon in a controlled
way such that the atomic nature of the fracture process can be revealed in a macroscopic
experiment. In the friction part, I will present another experiment using almost exactly the
same setup as for the low temperature fracture experiment to study the properties of static
friction and explore ideas concerning the origin of friction. In the last part, I will construct
a confined granular packing and study how pulses and continuous waves propagate through
it. All these three topics are relevant to geophysical science. I sincerely hope that my study
can ignite some fresh thinking in that area and help other researchers to design models that
can make more precise earthquake predictions.
vi
Contents
Acknowledgments iv
Abstract vi
List of Figures x
Chapter 1 Introduction 1
1.1 Fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.1 Basics of fracture mechanics . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.1.1 Three modes of fracture . . . . . . . . . . . . . . . . . . . . . 6
1.1.2 A brief review of linear elastic fracture mechanics . . . . . . . . . . . . 6
1.1.2.1 Pioneer work by Inglis . . . . . . . . . . . . . . . . . . . . . . 6
1.1.2.2 Energy balance approach by Griffith . . . . . . . . . . . . . . 9
1.1.2.3 Critical stress intensity factor approach by Irwin . . . . . . . 10
1.1.2.4 Beyond the scope of linear elastic fracture mechanics . . . . 11
1.2 Friction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2.1 A brief history of friction . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2.2 Dynamic aspect of static friction and its origin . . . . . . . . . . . . . 16
1.2.3 Purpose of the friction experiment . . . . . . . . . . . . . . . . . . . . 17
1.3 Granular Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.3.1 Introduction to LAMMPS . . . . . . . . . . . . . . . . . . . . . . . . 22
1.3.2 Force-displacement relation and Hertz contact . . . . . . . . . . . . . 22
1.3.3 The attenuation problem and the energy dissipation mechanism . . . . 23
vii
1.3.4 The resonant frequency shift problem and earthquake triggering . . . 24
Chapter 2 Fracture 25
2.1 Discrete nature of fracture process . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.1 Theoretical study of fracture process . . . . . . . . . . . . . . . . . . . 26
2.1.2 Numerical study of fracture process . . . . . . . . . . . . . . . . . . . 28
2.1.3 Previous experimental study of fracture process at CNLD (Center for
NonLinear Dynamics) . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2 Experimental study of crystal silicon fracture at liquid Nitrogen Temperature 41
2.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.2.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.2.3 Potential drop technique . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.2.4 The data acquisition circuit . . . . . . . . . . . . . . . . . . . . . . . . 44
2.2.5 Some major difficulties of this experiment . . . . . . . . . . . . . . . . 46
2.2.5.1 Wafer not breaking straight . . . . . . . . . . . . . . . . . . . 46
2.2.5.2 Wafer sliding . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.2.6 Differential loading to catch velocity gap . . . . . . . . . . . . . . . . . 52
2.3 Data Analysis and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.3.1 Data of partial cracks . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.3.2 Data of a good room temperature crack . . . . . . . . . . . . . . . . . 55
2.3.3 Data of a good low temperature crack . . . . . . . . . . . . . . . . . . 57
2.3.4 Velocity gap confirmed . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Chapter 3 Friction 60
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.2 The inclined plane experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.3 The friction experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.4 Modified rate- and state- friction law . . . . . . . . . . . . . . . . . . . . . . . 67
3.4.1 Rate- and state- equation and its modification . . . . . . . . . . . . . 67
3.5 Using the modified rate- and state- equation to fit our data . . . . . . . . . . 73
3.5.1 Step Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
viii
3.5.2 Continuous data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
Chapter 4 Granular Simulation 80
4.1 Preparation of a confined granular system . . . . . . . . . . . . . . . . . . . . 80
4.1.1 LAMMPS' implementation of the force-displacement relation . . . . . 80
4.1.2 packing generation and boundary conditions . . . . . . . . . . . . . . 81
4.2 Static properties of the prepared system . . . . . . . . . . . . . . . . . . . . . 83
4.3 Dynamic properties of the prepared system . . . . . . . . . . . . . . . . . . . 83
4.3.1 Pulse mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.3.1.1 speed of a pulse . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.3.1.2 attenuation and total absorbing layer . . . . . . . . . . . . . 88
4.3.2 Resonance mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.4 summary and brief report of ongoing studies . . . . . . . . . . . . . . . . . . . 102
Chapter 5 Conclusions 103
Chapter 6 Appendix 104
6.1 Wafer preparation for the fracture experiment . . . . . . . . . . . . . . . . . . 104
6.2 Bug report to LAMMPS admins . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.2.1 The input scripts are as following: . . . . . . . . . . . . . . . . . . . . 109
6.3 Graphs of Relative Increase of Oscillation Intensity . . . . . . . . . . . . . . . 112
Bibliography 117
Vita 125
ix
List of Figures
1.1 Glass cups drop on a concrete floor (photo from the World Wide Web) . . . . 3
1.2 Creep of a crack at the face of a dam (photo from the World Wide Web) . . . 3
1.3 Surface of a straight crack (photo courtesy of Robert Deegan) . . . . . . . . . 4
1.4 A wavy crack of silicon (photo courtesy of Robert Deegan) . . . . . . . . . . 4
1.5 Thermal energy can destroy the velocity gap and make it possible for a crack
to run at any speed. Only when it is largely reduced will the velocity gap be
detectable. (Graph adapted from Holland and Marder's paper[35].) . . . . . . 5
1.6 Schematic pictures of a bond-breaking process at zero temperature provided
an intuitive argument for the existence of a velocity gap. . . . . . . . . . . . . 7
1.7 Basic modes of fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.8 Inglis infinite plate problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.9 Schematic drawing of a crack propagating from x0 to x1 . . . . . . . . . . . . 9
1.10 ancient people using logs to reduce friction, This picture is adapted from
Duncan Dowson's classic book History of tribology[17]. . . . . . . . . . . . . 13
1.11 Drawing of Leonardo da Vinci, This picture is adapted from Duncan Dowson's
classic book History of tribology[17]. . . . . . . . . . . . . . . . . . . . . . . 13
1.12 inter locking mechanism of friction, This picture is adapted from Duncan
Dowson's classic book History of tribology[17]. . . . . . . . . . . . . . . . . 14
x
1.13 Schematic drawing of Bowden and Tabor's experiment of putting two cylin-
drical rods into a cross configuration and measuring the contact resistance,
which was then compared with the real area of contact from the analytical
solution given by Hertz[33, 44] (This graph is adapted from Bowden and
Tabor's paper[86].) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.14 Colored areas were the real area of contact, which was much smaller than the
apparent contact area (The whole picture was the apparent contact area.).
(This photo is adapted from Dieterich and Kilgore's paper[16] .) . . . . . . . 16
1.15 Part of a 3-D reconstruction of a static granular packing realized by the X-ray
micro-tomography method (This picture is adapted from Patrick's paper[21,
72, 84] .). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.16 Birefringent pattern of a photo-elastic disk was used to infer the internal stress
state of every individual disks, which was then assembled to find the contact
force network of the whole packing (This photo is adapted from Martin van
Hecke's paper[6, 89].). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1 One-dimensional lattice model used to simulate a dynamic fracture pro-
cess.(This graph is adapted from Marder and Gross's paper[60].) . . . . . . . 27
2.2 Solution of the 1-d lattice model gave the first indication of the existence of
a velocity gap. A crack would not propagate even the energy stored in the
system was more than the minimum energy required to sustain a running
crack based on the Griffith's energy balance criterion, which was known as
lattice trapping. For a fast running crack (Started at the up right corner of the
graph.), if the energy available to the crack tip gradually decreased towards
the Griffith's threshold, a stable solution did not exist for v smaller than
0.4, hence the crack would abruptly come to a stop. This range of velocity
where no stable solution existed was then called a velocity gap. (This graph
is adapted from Marder and Gross's paper[60].) . . . . . . . . . . . . . . . . . 29
2.3 Silicon fracture simulated with MEAM, which was then compared with ex-
perimental results. (This graph is adapted from Swadener's paper[39].) . . . . 32
xi
2.4 Comparison of numerical simulation results[35, 7] and experimental results[42]
for silicon fracture. All classical Molecular Dynamics simulations disagree
quantitatively with experiment (MEAM, IMSW). The tight binding results
of Bernstein and Hess appear to be compatible with experiment. . . . . . . . 33
2.5 Snapshot of MD simulation with EDIP (Adapted from Bernstein's paper) . . 35
2.6 Comparison of the range of SW potential and EDIP (Adapted from Bern-
stein's paper) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.7 Gross's tensile testing machine, Whomper (This graph is adapted from
Hauch's PhD thesis[32].) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.8 Hauch's sandwich design (This graph is adapted from Hauch's paper.) . . . . 39
2.9 Comparison of experimental results to MD simulation results (This graph is
adapted from Hauch's paper[42]. ) . . . . . . . . . . . . . . . . . . . . . . . . 40
2.10 MD simulation of silicon crack including thermal energy (This picture is
adapted from Holland's paper[35].) . . . . . . . . . . . . . . . . . . . . . . . . 40
2.11 3-D block diagram and 2-D side view of the friction holding system. A silicon
sample sits on top of a steel frame whose middle portion is milled out. Two
pressure blocks each have two feet; one sits on top of the steel frame and the
other on top of the silicon sample right along the edge of the middle window
of the frame. Two C-shape clamps each having a pair of bellows complete
the force loop. This works much the same way as we use our hand to grab a
burger; the C-shape clamp acts as a hand and the steel frame, silicon sample
and the pressure block are the contents of the special burger. . . . . . . . . . 43
2.12 Potential map and current flow through a partially cracked silicon wafer (Nu-
merical solution of a boundary condition partial differential equation given
by Matlab.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.13 The relationship between the potential difference of the inner two tabs and
crack length under constant current assumption. Each dot represents a solu-
tion for a particular crack length and the line is just an interpolation to show
the general trend of the data. . . . . . . . . . . . . . . . . . . . . . . . . . . 45
xii
2.14 Data acquisition circuit for the fracture experiment. The resistance of the
wafer is represented by the resistor between A and B. With a constant current
source on the order of 1 mA, the potential difference between A and B is
measured. Channel 1 picks up the voltage reading directly while channel 2
takes the same voltage data but makes an analog differentiation before the
data is recorded. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.15 It is demonstrated through calculations within software package franc2d that
if only a few percent of the elastic energy goes into Mode II, the crack will
deviate from its original path, which is not desirable in our experiment. . . . 47
2.16 Robert Deegan's three-hole frame design. This design would work if every-
thing functioned ideally, but the reality is that the world we live in is far
from ideal. So, the two piezoelectric transducers can be misaligned with re-
spect to each other and not perfectly perpendicular to the major axis of the
frame. The frame is very vulnerable to any tiny amount of shear force. So in
practice, this design does not work. . . . . . . . . . . . . . . . . . . . . . . . 49
2.17 Boxed low temperature environment surrounding Whomper, The upper sub-
figure is a real photo of the setup; the lower sub-figure is a schematic drawing
of the essential components. A stepper motor is connected to the left-most
steel rod (orange), a steel bar redistributes the load to two smaller rods (blue)
which in turn act upon the steel frame. The steel frame is still the platform
for conducting the fracture experiment. A symmetrical design refocuses all
the load back to the right most steel bar (orange), which is then fixed to
a stationary point on the wall of the outer aluminum structure. Since the
stepper motor is also fixed to the outer aluminum structure on a symmetrical
point of the opposite wall, the force loop is completed. This design has the
advantage of aligning itself and displays much stronger resistance to shear
forces than the design in Figure 2.16. . . . . . . . . . . . . . . . . . . . . . . . 50
2.18 Schematic drawing of mode I crack and quantities used for calculating energy
release rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.19 Bellow for differential loading . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
xiii
2.20 Length and velocity records of cracks moving in jerky fashion (partial cracks).
Note that the velocity never exceeds 1 km/sec and never reaches a steady
state. These experimental signals are characteristic of misaligned samples
that leave behind rough fracture surfaces. . . . . . . . . . . . . . . . . . . . . 55
2.21 A good room temperature break; the crack reaches a velocity over 1 km/sec
and holds it for several microseconds. The distance is calibrated by noting
the final position of the crack tip after the experiment is completed. . . . . . 56
2.22 A good low temperature break. . . . . . . . . . . . . . . . . . . . . . . . . . . 57
2.23 Comparison of the arrests of cracks at room temperature and liquid nitrogen
temperature under comparable energy release rates reveals the existence of a
velocity gap. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.1 Common understanding of friction (The lower portion of this graph is from
World Wide Web.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.2 Slip-stick motion (This graph is adapted from Cochard's paper[1].) . . . . . . 62
3.3 Rough surfaces in contact as a conventional picture of friction . . . . . . . . . 63
3.4 Inclined plane experiment between a polished silicon wafer and a machine-
ground steel surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.5 Schematic diagram of the experimental setup of friction. The four rods in
blue color are the aluminum holders for the two pairs of inductive sensors,
which are machined to exactly the same specification. . . . . . . . . . . . . . 66
3.6 The red and gray thin curves were individual readings from each inductive
sensor and they varied a lot during the process of inflating or deflating the
bellows, but their sum stayed very close to zero as more clearly indicated by
the sub-figure, which is a magnified view of the sum of these two curves. . . . 68
3.7 One set of continuous data under the normal force of 120 N. The green curve
is the prediction of conventional friction, the black curve is the prediction of
the modified rate- and state- friction law treating Dc as a material constant
while the red curve is the prediction of the same modified rate- and state-
friction law but allowing Dc to vary as a function of normal force. . . . . . . 70
xiv
3.8 I sat the silicon sample on top of the steel frame and waited for four different
lengths of time before conducting exactly the same experiment. The results
bear no indication of aging except for a slight difference between data obtained
after waiting for 30 minutes and the rest of them which were obtained after
waiting for much longer time. All the rest of my friction experiment waited
for more than 24 hours, so aging should not present itself in any of my other
data sets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.9 Step data; The normal force for this data set is 180 N and it was held fixed
throughout the experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.10 Slipping of the wafer as the boundary condition was fixed . . . . . . . . . . . 75
3.11 Illustration of the quantities measured for the friction experiment . . . . . . . 76
3.12 Continuous data from experiments where the normal force was constant and
shear force increases slowly and linearly. The green curve shows the sam-
ple response that would be expected from the conventional theory of static
friction. The sample would stretch with the frame up to the limit of static
friction and then slide freely. Instead, as shown by the experimental data
in blue, the sample always moves less than predicted. I provide two fits to
the data based on the modified rate and state friction law. The first of them
(black) treats the asperity length Dc as a constant. The second of them (red)
allows Dc to vary linearly with normal force. All parameters of the model
are held constant across all six panels. . . . . . . . . . . . . . . . . . . . . . . 78
4.1 A snapshot of the system and its contact force network . . . . . . . . . . . . . 84
4.2 Probability distribution curve . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3 Shape of the cosine pulse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.4 A snapshot of a pulse propagating through the system, What point in time? . 86
4.5 Horizontally averaged y-component velocities at equal time interval (colors
to indicate different instants of time.) . . . . . . . . . . . . . . . . . . . . . . 87
4.6 Hysteresis loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.7 sine wave transforms to triangular wave . . . . . . . . . . . . . . . . . . . . . 89
4.8 With total absorbing layer (colors to indicate different instants of time.) . . . 91
4.9 Without total absorbing layer (colors to indicate different instants of time.) . 92
xv
4.10 Horizontal average of y-velocity at resonant condition . . . . . . . . . . . . . 94
4.11 Fit the force on the bottom wall . . . . . . . . . . . . . . . . . . . . . . . . . 95
4.12 resonance curves at 300 kPa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.13 resonance curves at 3.1 MPa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.14 resonance curve at 11.9 MPa . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.15 Frequency shift as a function of amplitudes . . . . . . . . . . . . . . . . . . . 99
4.16 Relative increase of oscillation intensity of the whole contact force network . 100
4.17 Oscillation intensity comparison of driving amplitude 3, 10, 30, 100 to 1 . . . 101
6.1 Two identical balls collide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.2 y-force on ball 2 with radius 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.3 y-force on ball 2 with radius 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
6.4 collision simulated in LAMMPS (2-d and 3-d) . . . . . . . . . . . . . . . . . . 110
6.5 A3/A1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.6 A10/A1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
6.7 A30/A1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6.8 A100/A1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
xvi
Chapter 1
Introduction
To start my thesis, I would like to quote Richard Feynman's comment on what might be
the most important scientific discovery of human being[20].
If, in some cataclysm, all of scientific knowledge were to be destroyed, and only one
sentence passed on to the next generations of creatures, what statement would contain the
most information in the fewest words? I believe it is the atomic hypothesis (or the atomic
fact, or whatever you wish to call it) that all things are made of atomslittle particles
that move around in perpetual motion, attracting each other when they are a little distance
apart, but repelling upon being squeezed into one another. In that one sentence, you will
see, there is an enormous amount of information about the world, if just a little imagination
and thinking are applied.
There is little doubt that the discovery of atoms and the laws governing the dynamics
in the microscopic world - quantum mechanics are the most significant achievement of human
intelligence. But the microscopic world is so small that most of the practical engineering
problems start by making the continuum assumption. To many, the atomic discontinuity is
just theoretical. Out of sight, out of mind.
In this thesis, I will present three separate yet closely related topics that are more or
less concerned the effects of atomic discreteness reflected on a macroscopic scale. The three
topics are: fracture, friction and granular simulation. In the fracture part, I will describe
how to break a piece of silicon in a controlled way such that the atomic nature of the fracture
1
process can be revealed in a macroscopic experiment. In the friction part, I will present
another experiment using almost exactly the same setup as for the low temperature fracture
experiment to study the properties of static friction and explore ideas concerning the origin
of friction. In the last part, I will present a numerical study of mechanical properties of a
confined granular medium.
1.1 Fracture
Fracture is a fascinating topic. Everything breaks, yet things can break in dramatically
different ways: some break within a fraction of a second as when glass cups drop on a
concrete floor. (figure 1.1), while others can take years, such as the creep of a crack at the
face of a dam. (figure 1.2). Some break straight and leave mirror-like surface. One nice
example is illustrated in figure 1.3, which is a high resolution image of the crack surface from
a straight crack sample (courtesy of Robert Deegan). Others can break in a wavy manner
with rich surface features. Dipping a pre-notched hot silicon strip into a bath of cold water
can do the trick. Figure 1.4 shows that silicon can be broken in a wavy manner even it was
exactly the same material as for the straight crack.[76] (courtesy of Robert Deegan).
In the first part of this thesis, I will describe the searching for an interesting dynamic
feature of brittle fracture, namely a velocity gap, which was predicted by Marder in theo-
retical arguments [61] as well as numerical simulations [35, 34]. The theoretical argument
started with an analytical solution of a mass-spring lattice model for Mode I fracture first
proposed by Slepyan [81]. The solution of the lattice model indicated that the crack speed
will be unstable within a certain range (ranging from zero to about 20% of Rayleigh wave
speed). That is, you can not have a crack running steadily at that speed. The physical rea-
son for that has to do with the very nature of the fracturing process, which is the successive
breaking of inter-atomic bonds. Referring to figure 1.6, suppose at some moment the system
is in configuration (a), right after that the bond between the green atom and the atom to
its lower left breaks as in (b). Then, the system evolves from (c) to (e) while it prepares
for the breaking of the next bond, which is indicated in (f). The intuitive argument is that
the process from (c) to (e) can not take too long, otherwise the energy concentrated at the
crack tip will be dissipated through vibrations (sound) and the breaking of the next bond
2
Figure 1.1: Glass cups drop on a concrete floor (photo from the World Wide Web)
Figure 1.2: Creep of a crack at the face of a dam (photo from the World Wide Web)
3
Figure 1.3: Surface of a straight crack (photo courtesy of Robert Deegan)
Figure 1.4: A wavy crack of silicon (photo courtesy of Robert Deegan)
4
Figure 1.5: Thermal energy can destroy the velocity gap and make it possible for a crack
to run at any speed. Only when it is largely reduced will the velocity gap be detectable.
(Graph adapted from Holland and Marder's paper[35].)
could be delayed forever. That also means the crack speed can not be lower than a critical
value otherwise the crack will be unstable. The critical value is then called the velocity gap.
The intuitive argument above did not take thermal energy into account, but random
thermal fluctuations are able to break bonds. A fracture process due to thermal agitation is
also known as creep. For any none-zero temperature (absolute scale), thermal fluctuations
will always be present, which is also why I have to conduct the silicon fracture experiment
at temperatures less than 100 Kelvin. This conclusion was based on a state-of-art MD
simulation of silicon fracture at non-zero temperatures by Holland and Marder[35]. Figure
1.5 represents their main result.
Effort was put to the search for this elusive velocity gap[29, 32, 23] since we physicists
5
know that any theoretical prediction, no matter how intuitively sound it is, has limited value
if it is not backed up by solid experimental evidence.
But, before I go to the details of the experiment, let me first review some of the
basics of fracture mechanics.
1.1.1 Basics of fracture mechanics
1.1.1.1 Three modes of fracture
The first thing to know about fracture is that there are three modes of fracture, linearly
independent, and any fracture can be a linear superposition of them.[3, 78] They are named
Mode I, Mode II and Mode III. Mode I is also known as the opening mode, where the
direction of the applied external force is in the same plane as the material and perpendicular
to the direction of crack. My experimental study of silicon fracture is done exclusively in
this mode. The other two modes are: sliding mode (mode II) and tearing mode (mode III).
A graphical illustration of the basic modes of fracture can be found in figure 1.7.
The reason that I chose mode I for my study is that it is the most clean mode.
Friction between the two crack surfaces does not concern me in that mode of fracture as
compared to the other two modes. Also, a numerical study with a matching condition was
done by Holland[34, 35] and Marder and can be used for direct comparison.
1.1.2 A brief review of linear elastic fracture mechanics
1.1.2.1 Pioneer work by Inglis
The first attempt to quantify the effect of stress concentration due to a crack-like defect
within an elastic body was carried out by C. E. Inglis[37] in 1913. He studied the stress
distribution around an elliptical hole in an infinitely large linear elastic plate loaded at its
outer boundaries (illustrated in figure 1.8).
By demanding the two semi-axis to satisfy b a, the elliptical hole mimics a crack.
He showed that the stress maximum at point A (point A is on the edge of the ellipse and
directly ahead of the major axis) is given by
σA = σ(1 +
2a
b
) (1.1)
6
Figure 1.6: Schematic pictures of a bond-breaking process at zero temperature provided an
intuitive argument for the existence of a velocity gap.
7
(a) Mode I: opening mode (b) Mode II: sliding mode (c) Mode III: tearing mode
Figure 1.7: Basic modes of fracture
Figure 1.8: Inglis infinite plate problem
8
Figure 1.9: Schematic drawing of a crack propagating from x0 to x1
Obviously, there is a singularity for stress at point A if b approaches zero and
a remains finite. That is of course based on the assumption that the plate is infinitely
continuous, which is incorrect for any real material. But, the continuum assumption is so
fundamental for any mathematical analysis of this kind, it became a basic assumption for
the whole fracture mechanics discipline.
1.1.2.2 Energy balance approach by Griffith
When a crack propagates, new surfaces are created. Surface energy are always positive.
Otherwise, the material will spontaneously break into pieces to minimize its total energy.
This is also why bubbles are always spherical. They minimize their energy by minimizing
the surface. Thus, external energy must be supplied to compensate for the energy cost of
these new fracture surfaces.
Based on the energy balance argument, Griffith [27, 28] provided a criterion for the
minimum energy release rate required to make a crack propagate in the 1920s.
Referring to a schematic drawing of a propagating crack in figure 1.9, suppose there
is a crack front originally at some position x0, as it moves to x1, two new surfaces are
created.
The energy cost of creating these two surfaces is:
E = 2Γ(x1− x0)t (1.2)
where Γ is the surface energy density, t is the thickness of the material.
So, the energy release rate is given by eq. (1.2):
9
dE
dl
=
lim
x1→ x0
δE
δ(x1− x0) = 2Γt (1.3)
Of course, the energy release rate calculated this way only provides the absolute
minimum energy required to kick off a crack. It certainly does not guarantee that once
this minimum energy requirement is met, the crack will have to propagate, and usually
it won't. This also leads to a fairly interesting puzzle mentioned in Marder's condensed
matter physics textbook ([64], page 344). There he argues that, if only surface energy
was concerned (which he estimated Γ = 1J/m2 for rocks), all rocks higher than 1.4 cm are
unstable due to gravitational potential energy. In reality, mountains several kilometers high
can stand for millions of years. So, what was missing in this simple argument?
Griffith resolved this puzzle. He was actually the first one to perform a series of
fracture experiments on glass fibers and metal bars of different lengths and cross-sections
in order to find out the relation between atomic bonding strength and fracture strength
of materials. He eventually came to the conclusion that it was not the atomic bonding
strength but rather the maximum size of defects that determined the fracture strength of
the material.
The key to the puzzle is that there is an energy barrier, a very high energy mountain
to overcome in order to get a crack initiated. As I will explain in detail later on, the same
crack tip initiation problem caused some technical difficulties for my experimental study of
silicon fracture, too.
1.1.2.3 Critical stress intensity factor approach by Irwin
Some 30 years after Griffith's energy balance criteria for crack propagation, Irwin[38] in the
Naval Research Lab studied the stress distribution around a crack tip. He found that all
stress distributions around a crack tip share a universal functional form as given by the
formula below
K√
2pir
f(θ) (1.4)
That is, for any given static load, the leading term in the expansion of stress dis-
10
tribution in the vicinity of a crack tip will always be a constant (K) divided by the square
root of the distance to the crack tip times a function of angle dependence. The value of K
is determined by the geometry as well as the loading. What really makes this formula so
valuable is the fact that there is a load and geometry independent material constant KIC ,
which can be experimentally determined under well controlled external conditions, I indi-
cates that it is for mode I and the subscript c stands for critical. For K > Kc the crack tip
is supposed to propagate forward. So, researchers can measure the values of those critical
stress intensity factors for many standard materials and publish them. Engineers can then
calculate the stress intensity factor for the particular shape and size they want to construct
and compare it to the critical stress intensity factor of the material desired. If the value
of the stress intensity factor is greater than the critical one, the structure is bound to fail,
otherwise, it is safe to use (with proper safety margin).
Irwin then proved that this propagation criteria was equivalent to the energy release
rate criteria given by Griffith in equation 1.5, which is close to the relation between force
and energy in Newtonian mechanics.
G = β
K2
E
(1.5)
Where β=1 for plane stress and β = 1 − ν2 for plane strain, ν is the Poisson ratio
and E is the Young's modulus. 1
1.1.2.4 Beyond the scope of linear elastic fracture mechanics
Since 1950s, Linear Elastic Fracture Mechanics (LEFM) become a mature discipline and
research went beyond the linear assumption. Fatigue crack growth, elastic-plastic fracture
mechanics, dynamic fracture mechanics and creep and visco-elastic fracture[22] are just a
few of the major branches in fracture mechanics and the discussion of them is beyond the
scope of this thesis.
1Plane stress is a condition when the third dimension is much smaller than the other two and plane strain
is the opposite extreme.
11
1.2 Friction
Like fracture, friction is another ubiquitous phenomenon. We can not live without it, yet
we lose billions of dollars because of it. We have been studying it since the dawn of modern
science, yet some of its fundamental issues are still actively debated in today's scientific
literature[26, 47, 24, 67, 79]. In the second part of the thesis I will describe an experiment
of dry static friction between two hard surfaces. The purpose of it was originally to trou-
bleshoot the setup for the low temperature fracture experiment, but the result turned out
to be valuable on its own. Basically, I found that for shear force much less than the static
threshold µsN (where µs is the static friction coefficient and N is the normal force), two
surfaces will experience a reproducible relative sliding before they became locked into each
other. This finding inspired us to re-exam the existing phenomenological theory regarding
the origin of friction, namely the rate- and state- friction law.
1.2.1 A brief history of friction
Ancient people accumulated enough empirical knowledge about friction that they were able
to apply it in constructing many of the monumental constructions that still puzzle many of
today's engineers. The Great Pyramid in Egypt and the Great Wall in China are just two
outstanding examples of them. Figure 1.10 shows an example of how ancient people used
logs to reduce friction while moving huge stone monument.
The first scientific work on friction is usually attributed to Leonardo da Vinci (1452-
1519). He probably did some experiments to study the relationship between the magnitude
of a frictional force and the apparent contact area. A good indication of that is some of
his remaining hand-drawings, one of which is shown in figure 1.11. This is truly genius
considering that it was nearly 200 years before Newton wrote down his equation of motion.
Leonardo da Vinci might well be puzzled to find out that the contact area has very
little influence on the amount of frictional force. An intuitive picture of friction he likely
had in his mind was more or less like that shown in figure 1.12. Based on this picture, it is
easy to believe that the frictional force should be linearly proportional to the contact area,
while nature clearly disagrees with our intuitive.
The puzzle of the contact area remains a mystery until Bowden and Tabor[86] pro-
12
Figure 1.10: ancient people using logs to reduce friction, This picture is adapted from
Duncan Dowson's classic book History of tribology[17].
Figure 1.11: Drawing of Leonardo da Vinci, This picture is adapted from Duncan Dowson's
classic book History of tribology[17].
13
Figure 1.12: inter locking mechanism of friction, This picture is adapted from Duncan
Dowson's classic book History of tribology[17].
posed in the 1930s-1940s that the real contact area is only a small fraction of the apparent
contact area and the size of the real contact area is proportion to the normal load. The
experimental method they used was to measure the conductance between two metal parts in
contact. First, they put two cylindrical metal bars in cross contact as shown in figure 1.13.
The contact area inferred from this setup agreed well with the theoretical calculation based
on Hertz contact theory [33, 44] . Then, they put two flat surfaces in contact, the area of
contact thus inferred was much less than the apparent contact area. So, they concluded as
I quote The flat surfaces are held apart by small surface irregularities which form bridges
of an essentially metallic nature. These bridges flow under the applied pressure or their
number increases until their total cross-section is sufficient to enable them to support the
applied load.
Bowden and Tabor's idea about surface irregularities and bridges supporting load
are further confirmed by direct measurements of the real contact area. An optical means
of measuring the real contact area developed by Dieterich and Kilgore[16] works like this:
Two transparent flat objects are put in contact and a beam of monochromatic light passes
through them. Light is scattered by the rough non-contacting parts of the surface, but
transmitted through the real contacts. So, if Bowden and Tabor were right, then most of
the film should be dark while only a few scattered spots are visible. Figure 1.14 was one of
the resultant photos. Obviously, it looked just as predicted.
Bowden and Tabor's explanation is now widely accepted as the canonical picture
14
Figure 1.13: Schematic drawing of Bowden and Tabor's experiment of putting two cylin-
drical rods into a cross configuration and measuring the contact resistance, which was then
compared with the real area of contact from the analytical solution given by Hertz[33, 44]
(This graph is adapted from Bowden and Tabor's paper[86].)
15
Figure 1.14: Colored areas were the real area of contact, which was much smaller than the
apparent contact area (The whole picture was the apparent contact area.). (This photo is
adapted from Dieterich and Kilgore's paper[16] .)
of friction, but it remains a macroscopic theory, a statistical average of multi-asperities;
hence it is unable to account for dynamics on the single asperity level and address how fric-
tion originated from the atomic length-scale. There are basically three different theoretical
approaches to this problem[65, 52, 88, 67]:
1. Phenomenological rate- and state- models[77, 9, 16, 73];
2. Large scale molecular dynamics simulations[87, 8, 43, 70];
3. Minimalistic models[85, 57, 91, 50].
In this thesis I will concentrate on the first approach.
1.2.2 Dynamic aspect of static friction and its origin
Another puzzling aspect of friction is that a frictional force slowly grows with time, also
known as aging[14, 46]. To fit that into Bowden and Tabor's picture, you have to imagine
that the plastic flow is a very slow process and keeps going even the load has long been
stabilized. What is also known is that velocity will destroy the aging history and rejuvenate
16
the contact[9]. This led Dieterich, Rice and Ruina to propose a phenomenological model
that contained a state variable to describe the population of asperities and a rate depen-
dent dynamic equation to describe the dependence of the friction coefficient on the relative
velocity. Equation (1.6) is the original form of the rate- and state- equation they proposed.
µ = µ0 +A ln(1 +
v
v∗
) +B ln(1 +
θ
θ∗
) (1.6)
In its original form, Dieterich, Rice and Ruina kept a static coefficient of friction µ0
without raising too much doubt about its physical meaning. Basically, it means that for a
brand new contact at zero velocity, there is a residual shear resistance, the origin of which
is believed to be the cohesive force between atoms of the real contacts. But, what I found in
my experiment was that the value of µ0 was negligibly small and the origin of an apparently
nonzero µ0 could always be traced to a tiny sliding between the two surfaces before they
locked into each other and became motionless. The operational regime of my experiment
was well within the conventional range of static friction, yet the static friction only appeared
by growing from zero to µ0N by means of sliding the two interfaces a tiny distance. We
interpreted this as a dynamical origin of static friction (µ0N) and we incorporated it into
the rate- and state- equation by removing the constant µ0 and proposing a new functional
form of the state variable. The physical picture we had was fairly close to the interlocking
process imagined by Coulomb, except that the density of asperities was much lower (to
reflect the fact that it is the real contact area instead of the apparent one that supports the
load) and it was rate dependent.
1.2.3 Purpose of the friction experiment
The original purpose of the friction experiment was to find out how much the wafer slipped
before it cracked in low temperature fracture experiments.
In room temperature fracture experiments, Jens Hauch[42] used super glue to fix
the wafer onto the steel frame, and he was able to compensate for the yield of the glue as
the frame was stretched. This mechanism wouldn't work for the low temperature fracture
experiment because many glues would eventually lose their adhesive power when cooled
down. Even if a working glue could be found, it was still virtually impossible to apply the
17
glue after everything was cooled down. And if you glue them together at room temperature
first, the thermal contraction between silicon and steel were so much different that that alone
could break the wafer way before the temperature was cold enough to do the experiment.
So, Robert Deegan came up with the idea of using friction as the gripping mechanism.
He proposed to have all necessary parts cooled down first and then apply a strong enough
normal force through a pair of steel blocks to hold the wafer fixed onto the frame. The
trick was that the wafer sat loosely on top of the frame while cooling down, which avoided
the thermal contraction problem completely. There was a problem in this design however
in assuming that the wafer would not slide if a small shear force were to be applied. Based
on the conventional knowledge of friction, it was indeed correct to believe that no motion
was possible if a shear force less than a certain threshold was applied. But, after breaking a
few pieces of wafer, it became puzzling that all of them broke at much higher energy release
rate than those of room temperature experiments, while they were expected to be lower.
After careful inspection, it was found that the silicon wafer was sliding on top of
the steel frame while the frame was stretched. But, this was not supposed to happen when
a strong enough normal was applied based on conventional understanding of static friction.
So, a review of the current understanding of static friction was presented in the introduction
of the second part of this thesis. Armed with results from my experiment, Marder and I
proposed a modification to the standard rate- and state- friction law, which unified the
understanding of static and dynamic friction and provided a more realistic physical origin
of friction.
1.3 Granular Simulation
Granular materials can be found everywhere around us, such as a pile of sand, a jar of
marbles or just a bag of candies. Most of these granular systems are static due to the fact
that the interaction between each pair of their constitutive particles are non-conservative.
So, either you continuously pump energy into the system (for example, by shaking or rotating
the containers) or the system will quickly dissipate all of its kinetic energy and become static.
Most of these static granular systems are more or less confined due to the boundaries of
their containers or just the weight of the grains on top of them.
18
To study such a confined granular system, field and laboratory experiments are
conducted. Quantities that are commonly measured are usually macroscopic, such as the
force acting upon the surface of a piezo-electric transducer or the time it takes for a pulse
to travel certain distance (or reflect from a boundary).
Macroscopic quantities are limited in the sense that they are the results of averaging
in spacial or temporal directions, but sometime it is the inhomogeneities that are the subjects
of interest. So, different experimental methods are implemented to measure microscopic
quantities.
X-ray micro-tomography[21, 72, 84] can be used to measure exact position of every
particle in a static granular packing, Once that information is obtained, many different
characteristic parameters can be calculated, such as the pair correlation function, radial
distribution function and contact numbers. Figure 1.15 is a partial reconstruction of the
initial packing by the X-ray micro-tomography method.
Another commonly used experimental method to visualize the micro-structure of
a static granular packing is the photo-elastic disk method[6, 89]. It is a 2-D system with
confined transparent disks which would display unique birefringent patterns depending on
the contact forces when imaged through a pair of crossed circular polarizers. By solving an
inverse problem of those patterns, Majmudar and Behringer were able to find the normal
as well as tangential forces between each pair of contacting disks. Once that information
was obtained, it would then be possible to study the contact force distribution or contact
force network in different loading conditions, such as bi-axial compression or shear. Figure
1.16 illustrates a typical birefringent pattern of a photo-elastic disk in contact with three
neighbors.
An obvious disadvantage of those experimental methods is that they can only deal
with static packings, but many interesting features of a granular system are dynamic. For-
tunately, there is now a completely different approach to the study of a confined granular
system especially for studying dynamic quantities and quantities defined on the microscopic
length scale, namely computer simulations.
In the last part of my thesis, I will describe how to apply a general purpose MD
(molecular dynamic) code (LAMMPS) to study the motions of particles in a confined gran-
ular system. Two major subjects are extensively studied:
19
Figure 1.15: Part of a 3-D reconstruction of a static granular packing realized by the X-ray
micro-tomography method (This picture is adapted from Patrick's paper[21, 72, 84] .).
20
Figure 1.16: Birefringent pattern of a photo-elastic disk was used to infer the internal stress
state of every individual disks, which was then assembled to find the contact force network
of the whole packing (This photo is adapted from Martin van Hecke's paper[6, 89].).
21
1. Pulse mode - A single pulse of longitudinal vibration propagates through the system,
measuring the rate of attenuation as well as the speed as a function of confining
pressures and driving amplitudes.
2. Resonance mode - Based on the speed measured in the first part, the base resonant
frequency is estimated and a set of frequency response curves is determined. The
response frequency peak moves towards lower value as the driving amplitude increases,
which is known as the resonant frequency shift.
1.3.1 Introduction to LAMMPS
LAMMPS[75] is the acronym for Large-scale Atomic/Molecular Massively Parallel Simula-
tor. It is a general purpose Molecular Dynamics (MD) simulation package designed and
maintained by Sandia National Lab. It is free software under GNU General Public License
(GPL). It was originally written in Fortran in the mid 1990s and was re-written in C++ (re-
leased to public in 2004). It uses Message Passing Interface (MPI) for co-ordinating memory
and data across variable number of CPUs. For more information on LAMMPS, including all
versions of its source codes, please visit the following website: http://lammps.sandia.gov/
1.3.2 Force-displacement relation and Hertz contact
All MD codes integrate Newton's second law and update particle positions recursively. The
essential differences between them are the force-displacement relations. Those relations were
proposed by researchers to model specific systems and their success usually depends upon
how well the simulation reproduces a corresponding experiment or how well a prediction of
a future experiment it can make.
To model a granular system, it is essential to know how one particle interacts with
another, that is, their force-displacement relation. In the dawn of the last century, Hertz[33,
56, 69] solved the problem of calculating the repelling force of two elastic equal spheres in
contact. The result of his calculation can be summarized in the following formula:
Fij =
4
3
√
RijE
∗
ij(Ri +Rj − rij)
3
2 (1.7)
22
WhereRi Rjare the radius of two interacting balls and Rij is the reduced radius
defined as
Rij =
RiRj
Ri +Rj
(1.8)
Here E∗i , E
∗
j are the modified Young's modulus and E
∗
ij is the reduced Young's
modulus defined as
E∗ij =
E∗i E
∗
j
E∗i + E
∗
j
(1.9)
The definition of the modified Young's modulus is given by
E∗ =
E
1− ν2 (1.10)
E and ν are the Young's modulus and Poisson's ratio respectively.
1.3.3 The attenuation problem and the energy dissipation mecha-
nism
Macroscopic pulse transmission experiments usually only reveal how much energy was dissi-
pated but are incapable of pointing a finger to a specific mechanism that is responsible for the
energy lost. Phenomenological models have been proposed based on empirical knowledge in
specific fields. Two widely cited models for energy dissipation in nonlinear elastic materials
are the Preisach-Mayergoyz hysteretic elastic unit model (PM model[30, 31, 71, 59]) and the
array of energy sinks model[13]. Marder[62] (unpublished) devised a method to differentiate
them. He proposed that if a sine wave was to transmit through a medium described by the
hysteresis model, the sine wave will soon transform to a triangular wave, while if the same
wave is put through a medium described by the array of energy sinks model no change of
wave shape should be expected.
Since LAMMPS had proved to generate sensible result in modeling granular systems,
it is natural to wonder if it would favor one of the model against another. Also, by looking
into the microscopic details of the dynamics, more physical quantities could be revealed and
23
probably led to more physical models.
1.3.4 The resonant frequency shift problem and earthquake trig-
gering
Another problem that drove me to this field was the earthquake triggering mechanism
proposed by Jia and Johnson[45] etc. They argued that based on laboratory results and
field data, dynamic strain on the order of 10−6 could trigger an earthquake. The reason
they put forward was the well known but not well understood phenomenon called resonant
frequency shift, or dynamic strain softening. What they assumed was that there exists a
fault in a critical (near failure) state. If a fairly weak wave of oscillation hits the fault
causing temporary softening of the fault gouge, the softened fault gouge would no longer be
able to bear the load and would eventually generate an earthquake.
Since resonant frequency shift was at the heart of their argument, it was reasonable
to believe that a clear picture of the microscopic detail of particle motions at resonance and
a better understanding of underlining mechanism for the frequency shift could help to build
a better model that makes more precise earthquake predictions.
24
Chapter 2
Fracture
2.1 Discrete nature of fracture process
As mentioned briefly in the introduction, fracture at the smallest scale is a discrete process,
which consists successively snapping of inter-atomic bonds. Obviously, linear elastic frac-
ture mechanics (LEFM) and all other fracture theories that are based on the continuum
assumption are inadequate for capturing features that are rooted in the discrete nature of
fracture process, such as the velocity gap. To follow the dynamics of each and every atom in
a dynamic process seems forbiddingly hard, yet Slepyan and other pioneer researchers uti-
lized a simple mass-spring model, cleverly manipulated the dynamic equations and provided
meaningful solutions as well as intriguing predictions that inspired many young scientists to
join this field[80, 60]. Fortunately, with the incredible increasing power of modern comput-
ers, scientists are now able to study the fracture process by looking into the very details of
those discrete processes without making over-simplified assumptions and make quantitative
predictions that can be directly verified by experiments. In this chapter, I will try to present
a brief survey of the study of dynamic fracture process from three different points of view:
1. the mass-spring model (theoretical) approach[80, 60]
2. the numerical approach (mainly Molecular Dynamic simulations)[19, 60, 34, 39]
3. review of some previous experiments[42, 36]
25
2.1.1 Theoretical study of fracture process
I have no intention of presenting a comprehensive review of the theoretical study of fracture
process. Instead, I will concentrate on a one-dimensional lattice model that absorbed many
of the previous results and had some of its own. The model was proposed by Marder and
Gross[60] in 1994 and was built upon the work of Slepyan[80] and others.
The lower portion of figure 2.1 is one solution of the model, and the dynamics of
the crack tip are magnified to indicate the setup of the equations. The equation concerning
the dynamics of the atom (m,+) was given by the following equation:
u¨m,+ =

um+1,+ −2um,+ + um−1,+
+ 1N (UN − um,+)
+(um,− − um,+) θ(2uf − |um,− − um,+|)
−bu˙m,+
(2.1)
The trick is to find only the steady state solution and apply the following symmetries
um,+ = −um,− (2.2)
um,+(t) = u0,+(t−m/v) (2.3)
Applying the Wiener-Hopf methods, an analytical solution of the equations above
can be found.
The end product of the lattice model was a prediction of the crack speed as a function
of the driving force ∆.1
Figure 2.2 illustrates the prediction of the one-dimensional lattice model for a situ-
ation of N = 100 and v = 0.5.2 Notice that the horizontal line at velocity zero extends far
beyond ∆ = 1, meaning that excessive energy can be stored in the lattice model without
1The driving force ∆ was defined as a ratio of the displacement of the outer boundary UN and the
minimum value of UN (indicated as U
c
N ) at which enough energy was stored to the right of the crack so as
to break the bonds, ∆ ≡ UN
Uc
N
.
2N is the number of the atoms above the crack line, and v is the steady state crack speed.
26
Figure 2.1: One-dimensional lattice model used to simulate a dynamic fracture process.(This
graph is adapted from Marder and Gross's paper[60].)
27
triggering a crack propagating forward. This is known as lattice trapping, first proposed
by Thomson etc. in 1971. Another interesting feature of the graph is the zig-zag pattern
between v = 0 to v = 0.4. It turned out that those solutions were unstable and hence forbid-
den in steady state crack propagation (this forbidden speed range is then call the velocity
gap). Motivated by this finding, Marder et al were interested to see if the velocity gap was
an artifact due to the over-simplified model or it was a generic feature that exists in more
realistic settings. So, they devised a molecular dynamics simulation with much more sophis-
ticated inter-atomic potentials to simulate a propagating crack. Meanwhile, they started to
conduct fracture experiment using real silicon wafers.
2.1.2 Numerical study of fracture process
There is no doubt that the laws of physics in the microscopic scale is quantum mechanics
and the governing equation for the dynamics of electrons and nucleus is the Schrödinger's
equation. But this equation is analytically solvable only for the hydrogen atom. Various
approximation methods are developed to calculate properties of more complex atoms and
molecules. One popular example is the density functional theory (DFT) with the local
density approximation (LDA), which is believed to be the most reliable techniques available
for numerical treatment of materials. However, even that is restricted to perfect crystals
and unit cell sizes around 1000 atoms[51]. For problems involving larger numbers of atoms
(problems requiring a minimum of about 106 atoms to capture just some of the underlying
physics, such as fracture or dislocation loops), approximation methods with much greater
computational efficiency are certainly needed.
Molecular dynamics (MD) simulation with empirical inter-atomic potentials pro-
vides an alternative approach to these problems. Calculations within molecular dynamics
are completely classical, with all the information of the interactions between neighboring
(not necessarily nearest) atoms contained in the effective inter atomic potentials. Different
functional forms of the potentials with different numbers of fitting parameters and fitting
strategies give a variety of predictions concerning certain properties of materials. The valid-
ity and transferability of these potentials are then checked by comparison to experiments.
The advantage of molecular dynamics simulation is its speed and the ability to handle a
large number of particles. For the standard of about half a decade ago, ~ 107atoms followed
28
Figure 2.2: Solution of the 1-d lattice model gave the first indication of the existence of a
velocity gap. A crack would not propagate even the energy stored in the system was more
than the minimum energy required to sustain a running crack based on the Griffith's energy
balance criterion, which was known as lattice trapping. For a fast running crack (Started at
the up right corner of the graph.), if the energy available to the crack tip gradually decreased
towards the Griffith's threshold, a stable solution did not exist for v smaller than 0.4, hence
the crack would abruptly come to a stop. This range of velocity where no stable solution
existed was then called a velocity gap. (This graph is adapted from Marder and Gross's
paper[60].)
29
for a few tens of nanoseconds was achievable by molecular dynamics with fairly sophisticated
inter atomic potential. Today, simulations of billions of atoms for nanoseconds are possible.
Rapid brittle fracture of silicon is one typical example of the problems that requires a large
number of particles. Experimentally, it involves a macroscopically stressed piece of silicon
containing a number of atoms on the order of 1022, within which an atomically sharp crack
tip propagates. Although dynamics of the crack tip is determined microscopically by the
successive failure of the covalent bonds between atoms, the energy needed to break these
bonds is supplied macroscopically. Numerically, many attempts have been performed to
simulate brittle fracture of silicon with different inter atomic potentials, but no consensus
conclusion are reached yet, partly because none of the existing potentials is particularly
superior in predicting the material properties of silicon and many even failed rather badly
in reproducing the basic cohesive energy curve. There have been more than 30 empirical
inter atomic potentials invented for the description of different aspects of properties of sili-
con since 1980s. Among them, the Stillinger-Weber potential (SW)[82], the modified embed
atom method (MEAM)[12] and the environment-dependent inter-atomic potential (EDIP)[5]
are the most frequently used for the description of rapid brittle fracture of silicon.
Though the SW potential[82] has been very successful in predicting many material
properties of silicon, such as the elastic constants and thermal coefficients, it failed to predict
brittle fracture along the experimentally preferred cleavage planes (111) and (110)[35]. At
low or moderate strains a crack will not propagate at all. At high strains, the crack tip
region becomes very disordered. The SW potential does give a type of fracture along the
(100) plane[2], which is quite rough on the atomic scale. However it is known experimentally
that fracture surfaces along the (111) plane can be atomically flat[55].
A simple estimation of the fracture energies along different crystal planes by Grif-
fith's critical energy analysis shows that the minimum energy density required for a crack
to propagate along the (111) plane is around 2.78 Jm2 and along the (110) plane is around
3.28 Jm2 , but the required energy density for a crack to propagate along the (100) plane
is significantly higher, around 4.54 Jm2 based on Hauch's experiments[32]. So, a crack in
silicon is energetically favored to propagate along the (111) plane and (110) plane rather
than (100) plane, which is in agreement with experimental observations. Whether a crack
can propagate along the (100) plane in a brittle manner is still experimentally inconclusive
30
and numerically questionable. But the failure of prediction of brittle fracture along (111)
plane and (110) plane clearly indicates the shortcoming of the SW potential in predicting
the mechanical properties of silicon, especially when the system is far from equilibrium. In
order to get a brittle fracture of silicon, Marder etc. inadvertently proposed a modified SW
potential, in which the coefficient of the three body term had been multiplied by a factor of
2. This modification did predict a brittle fracture along the (111) and (110) planes, but it
worsens the likeness to real silicon in other respects such as the melting temperature been
raised to approximately 3500K while the experimental value is around 1683K and elastic
constants also been shifted some significant amounts[34]. Thus, they concluded that nei-
ther the Stillinger-Weber potential nor their modification can be regarded as a satisfactory
account for simulation of brittle fracture of silicon.
The embedded atom method (EAM) was invented in early 1980s by M. S. Daw
and M. I. Baskes[12] as a semi-empirical model originally designed for metals. The EAM
is based on local electron density theory, which has been shown to accurately describe a
large number of properties in metals. In the later 1980s, M. I. Baskes modified the original
EAM to describe the covalent bonding for materials such as silicon, and called the result
the modified embedded atom method (MEAM)[4].
The basic idea behind the EAM is that each atom in a solid can be viewed as an
impurity embedded in a host comprising all the other atoms. Because the energy of an
impurity is a function of the electron density of the unperturbed host, the cohesive energy
of a solid can be calculated from the embedding energy.
The MEAM was applied to molecular dynamic simulations of brittle silicon fracture
by J. G. Swadener etc. recently[39]. The simulations produce propagating crack speeds
that are arguably in agreement with experimental results. The results of the simulation and
comparisons with experimental results are summarized in figure 2.3.
The experimental results are indicated as open squares; The simulation results are
open cycles (model height=239A), full cycles ( 147A) and full triangles (112A); The solid
line is continuum prediction.
Actually, Holland and Marder had done a very similar study of brittle fracture
simulation of silicon by MEAM a year earlier[35]. Their results, as summarized in figure 2.4
, are quite different however. The predicted minimum energy density for a crack to propagate
31
Figure 2.3: Silicon fracture simulated with MEAM, which was then compared with experi-
mental results. (This graph is adapted from Swadener's paper[39].)
is greater than 5 Jm2 , while the corresponding experimental value is around 2.5
J
m2 .
The difference between Swadener's results and those of Marder's is not due to the
potential, since they were using exactly the same potential. So, Swadener[39] etc. claimed
that all previous simulations failed to agree with experiments due to the ignorance of the
corresponding variation of the effective Young's modulus and Rayleigh wave speed as the
strain within the silicon strip changed. If all of these changes were evaluated properly with
MD simulations adopting the MEAM potential, the simulation would be in agreement with
experiments. But, Marder strongly doubted this conclusion. Marder pointed out that some
of the formulas used in Swadener etc.'s paper are actually borrowed from continuum theory
and not applicable to the special geometry of the relevant problem. The energy stored
in the system could be evaluated directly from the simulation rather than applying some
approximate formulas as done by Swadener et al. The most serious defect of Swadener
et al.'s simulation might be that their system was too small in size and the crack tip had
not reached a steady state as the evaluations of the crack speed and fracture energy were
performed. On the contrary, Marder[60, 63] et al.'s simulation scheme didn't have such a
problem since they applied a method that continuously cut a portion of the system that had
32
Figure 2.4: Comparison of numerical simulation results[35, 7] and experimental results[42]
for silicon fracture. All classical Molecular Dynamics simulations disagree quantitatively
with experiment (MEAM, IMSW). The tight binding results of Bernstein and Hess appear
to be compatible with experiment.
33
fully relaxed after the crack propagated through and pasted a portion of stressed material
in front of the crack tip. By using this method, they could follow the dynamics of the crack
tip for a much longer period of time and found a truly steady state of the crack.
The environment-dependent inter-atomic potential (EDIP) was developed in later
1990s by Martin Z. Bazant[5] etc. The EDIP is similar in form to SW potential but has
an environmental dependence. The interaction model includes two-body and three-body
terms which depend on the local atomic environment through an effective coordination
number. There are 13 adjustable fitting parameters for this potential, half of which can be
theoretically estimated. The theoretical base underlying the model include an analysis of
elastic properties of the diamond and graphitic structures and inversion of ab initio cohesive
energy curves.
The reason that M. Z. Bazant etc. developed a new potential even when there were
already dozens of inter-atomic potentials for silicon in the literature was that they were
unsatisfied with the transferability of all the existing potentials. Most of the potentials give
reasonably good predictions only for limited aspects of silicon properties and it is often
unclear whether the good agreement of simulation and experiment is due to the functional
form of the potential, or the sophisticated fitting strategy.
EDIP had been used to simulate fracture in silicon by Bernstein and his coworkers
recently, but their simulation failed to agree with experiments[7]. The fracture energy they
predicted is almost 10 times the energy observed in experiments. Also, the crack tip in their
simulation proceeds in a very ductile manner accompanied by significant plastic deformation
and disorder, while brittle crack are seen experimentally. A snapshot of their simulation is
given in figure 2.5. It is worth mentioning that Holland and Marder had achieved similar
results two years before Bernstein etc. during their search for the best potential to describe
the brittle fracture of silicon[35].
The failure of all these empirical potentials mentioned above could be due to their
short-range cut-off, which produces unrealistically large cohesive force and prevents bonds
from breaking[35]. Restricted by the value of the surface energy, the depths of different po-
tentials are pretty much the same no matter how different their shapes or ranges might be.
Most of the empirical potentials of silicon such as SW and EDIP include only nearest neigh-
bor interactions, while density functional theory predicts an interaction range considerably
34
Figure 2.5: Snapshot of MD simulation with EDIP (Adapted from Bernstein's paper)
longer. A comparison of the ranges of different potentials is illustrated in figure 2.6.
It is obvious that a shorter cut-off will force a potential to rise more steeply, which
corresponds to a short-range but stronger attractive force. Therefore, even when it is ener-
getically favorable for the crack to propagate, the stress at the crack tip may not be high
enough to break the bond. This is one of the reasons why these empirical potentials failed
to simulate the brittle fracture of silicon.
The UER curve on the graph is the universal energy relation curve, which comes
from an empirical fit to the cohesive energy curve from density function theory[41].
Though the guess of a proper functional form for the potential is essential for the
success of the simulation, it is not the only problem worth concern. The size of the simulation
is in no way directly comparable to experiments even for the most powerful computer built
to date. In fact, a simulation involving ten billion atoms describes only a cube of matter
no more than half a micron along each side. And even for computers powerful enough
to handle a system with so many degrees of freedom, they will not be able to follow the
behavior of all the particles for much longer than a nanosecond. So, comparison of the
simulations with experiments requires a careful scaling analysis. There are two alternative
solutions for this problem. One is to couple molecular dynamics with a continuum simulation
dynamically[7]. The other one is through a clever scaling argument that reduces the system
size significantly and traces the crack tip for much longer time. A cut and paste procedure
35
Figure 2.6: Comparison of the range of SW potential and EDIP (Adapted from Bernstein's
paper)
was carefully manipulated to host the crack tip within the central region of the system[34].
There is yet another problem worth mentioning. It is how to determine whether the
system has reached the steady state or not. For a system of millions of atoms and driven
far from equilibrium, it is a nontrivial task to set a correct criteria for a steady state. If a
system is too small in size or the evolution of the system is not traced long enough, it is
quite possible that the system is not in a steady state but rather in a transient state.
Despite all the difficulties, numerical simulations of dynamic cracks become increas-
ingly realistic and provide a complementary means to study fracture process on the atomic
scale.
2.1.3 Previous experimental study of fracture process at CNLD
(Center for NonLinear Dynamics)
It hardly takes any effort to break a silicon wafer in a casual way. As a matter of fact, caution
has to be taken to prevent that from happening. Dropping it to the ground it shatters into
pieces since it is regarded as one of the perfect brittle materials at room temperature. But,
it is rather difficult to break a silicon wafer in a controlled fashion mostly due to the fact that
precisely controlled boundary condition is hard to achieve and is even harder to maintain
as the crack propagates.
Jay Fineberg and Steven Gross[29] designed a tensile testing machine as schemat-
ically illustrated in figure 2.7 and named it Whomper. It had a stiff outer skirt made of
36
Figure 2.7: Gross's tensile testing machine, Whomper (This graph is adapted from Hauch's
PhD thesis[32].)
bulk Aluminum. The integrity of the frame was further enforced by three Aluminum bars.
A pair of jaws each sitting on top of two pillow blocks were free to move along the long
axis but fixed on the other direction to ensure uni-axial loading. To conduct a fracture
experiment, a piece of sample was fixed to the jaws with a small crack initiated on one side.
A computer controlled stepper motor provided the driving force to tear the sample apart.
The crack speed and energy release rate were measured independently during the fracture
process and compared to theoretical predictions as well as computer simulations.
Jens A Hauch performed a series of fracture experiments on Whomper in later
1990s. The materials he used included Holimmate-100, PMMA, glass and silicon[32]. For
the purpose of searching for the velocity gap, I will only describe the dynamic fracture
experiments done in single crystal silicon.
37
In order to make a direct comparison of the experimental results with the prediction
of MD simulation possible, it was necessary to have a fixed boundary as the crack propagates.
So, Hauch[42] altered the original design by having the silicon sample sandwiched between
two identical frames of steel whose center were milled off as shown in figure 2.8. The two
side beams on each steel frame will act as very stiff extension elements when the whole
structure was subjected to uni-axial loading along y axis. Only a small fraction of the total
load (< 2 %) was transmitted through the silicon sample, so to a good approximation the
boundary can be regarded as fixed when the crack happens.
The results of Hauch's experiment and comparisons to MD simulations are shown
in the figure 2.9. There is an obvious discrepancy between experiment and simulation,
especially at the lower limit of energy release rate at which steady fracture could take
place. The MD simulation predicted almost twice as much energy release rate as the lowest
limit seen experimentally. As stated earlier, this discrepancy led Bernstein and Hess etc.
to propose that quantum mechanical calculation is inevitable for a precise account of the
fracture dynamics of single crystal silicon.
Unfortunately, Hauch's experiment, which was done at room temperature, wasn't
decisive in determining the existence of the velocity gap. Shortly after the propose of the
velocity gap, Marder and Holland realized that thermal fluctuations at room temperature
are energetic enough to kick off an arrested crack to propagate again. That is, thermal
agitation destroyed the visibility of velocity gap and the crack can run at any speed. So,
they modified their MD simulation of silicon fracture to include thermal energy effect. Figure
2.10 gives a snapshot of their simulation.
The result of their simulation indicated that it is necessary to conduct the fracture
experiment under 100 K in order to see the velocity gap without the disturbance of thermal
fluctuation. Figure 1.5 summarized their finding on how thermal energy destroyed the
existence of the velocity gap. The basic idea was that for a crack otherwise would have
stopped, a random thermal fluctuation snapped the bond right ahead of the crack tip and
re-initiated the propagation. Depending on how much thermal energy was available, this
process could happen frequently enough as if there was no velocity gap, that was, the crack
could run at any speed within the velocity gap. Until temperature was reduced under 100
Kelvin, a random thermal agitation became rear enough such that the velocity gap began
38
Figure 2.8: Hauch's sandwich design (This graph is adapted from Hauch's paper.)
39
Figure 2.9: Comparison of experimental results to MD simulation results (This graph is
adapted from Hauch's paper[42]. )
Figure 2.10: MD simulation of silicon crack including thermal energy (This picture is
adapted from Holland's paper[35].)
40
to be detectable, which was the reason that I had to conduct the silicon fracture experiment
in liquid nitrogen environment in order to find the velocity gap.
2.2 Experimental study of crystal silicon fracture at liq-
uid Nitrogen Temperature
2.2.1 Motivation
The primary objective of the low temperature silicon fracture experiment is to search for
the velocity gap.
2.2.2 Experimental setup
As mentioned earlier, Hauch had done some room temperature fracture experiments on
silicon. But, the apparatus he used is incompatible with the requirement of low temperature
fracture experiment. The main issue is the great difference of thermal contraction between
silicon and steel. It is estimated that if silicon is glued onto the steel frame and then
cooled down with it, a few tens of degree Kelvin is enough to break the silicon by thermal
contraction alone. Besides, the glue Hauch used would lose its adhesive power before the
desired temperature was reached and the Whomper is simply too big to cool down entirely.
To overcome these difficulties, new experimental designs were prepared by Robert
Deegan. The basic idea was to cool down first and stretch later. That is, silicon wafer would
be put in position and cooled down with the steel frame loosely, some mechanism would
then fix it upon the steel frame after the desired low temperature (less than 100 Kelvin)
had been reached and stabilized. Most glues don't work in that temperature range, and
even if they do, it is nearly impossible to apply them in situ. So, friction was chosen as the
mechanism to hold up the silicon wafer to the steel frame in low temperature.
Using friction as the mechanism for holding the wafer to the steel frame was believed
to work because we have a method to control the normal force exerted on a pair of steel
blocks through silicon wafer onto the steel frame. We can first put the silicon wafer in its
position and cool it down. While cooling, nothing but a small residual normal force due to
the weight of the pressure blocks is on the silicon wafer. When the desired temperature is
41
reached and it is time to break the sample, a much larger normal force is applied to the
top of the pressure block by inflating a set of four brass bellows, which push the pressure
blocks against a pair of C-shape steel clamp fixed to the frame. See figure 2.11 for a graphic
illustration of the friction holding mechanism.3
2.2.3 Potential drop technique
It is not a trivial task to measure the speed of a propagating crack in silicon since it is on
the order of kilometers per second while the size of the silicon wafer is only a few inches.
That is, it all happens within one micro-second if it happened at all. A very fast electrical
method must be applied to trigger at exactly the right moment and a fast and large enough
data acquisition device is needed to quickly capture, store and transfer data.
The principle of the potential drop technique is fairly simple: as the crack starts to
propagate, the resistance of the silicon wafer will change. Using the same principle as the
four point probe, a semi-constant current is fed through the outer two tabs of the wafer and
the voltage of the inner two tabs is measured. The advantage of this technique is that it is
insensitive to the contact resistance of the tabs to the wafer. By monitoring and recording
the potential drop across the inner two tabs at a high enough frequency, it is possible to
retrieve the information of the position of the crack tip as a function of time and hence the
crack velocity.
One necessary step of this technique is to invert voltage readings back to crack tip
position on the wafer. It is done by solving a two dimensional PDE with some specific
boundary conditions by the finite element method in Matlab.
Figure 2.12 represents a typical snapshot of the solution to the PDE. The geometry
exactly mirrors from the physical dimensions of the silicon strip under test. Currents are
not allow to flow through any of the outer boundaries. (The crack interface can be viewed
as part of the outer boundaries as well since it is implemented by cutting a very thin layer
of material off the main body of the strip.) A constant current flows from the up-most tab
through the body of the silicon strip and sinks at the bottom-most tab. Potential difference
is then measured between the inner two tabs and plotted as a function of crack length in
3Note that for purpose of clarification, only half of the system is shown, which includes one piece of
C-shape clamp one pressure block and two brass bellows. Complete system configuration is shown on the
side view.
42
Figure 2.11: 3-D block diagram and 2-D side view of the friction holding system. A silicon
sample sits on top of a steel frame whose middle portion is milled out. Two pressure blocks
each have two feet; one sits on top of the steel frame and the other on top of the silicon
sample right along the edge of the middle window of the frame. Two C-shape clamps each
having a pair of bellows complete the force loop. This works much the same way as we use
our hand to grab a burger; the C-shape clamp acts as a hand and the steel frame, silicon
sample and the pressure block are the contents of the special burger.43
Figure 2.12: Potential map and current flow through a partially cracked silicon wafer (Nu-
merical solution of a boundary condition partial differential equation given by Matlab.)
figure 2.13.
2.2.4 The data acquisition circuit
The electronic data acquisition circuit is shown in figure 2.14. A semi-constant current
source connected to the outer two tabs of the wafer and the voltage change across the inner
two tabs was monitored by the data acquisition circuit. As the crack starts to propagate,
the voltage on the inner two tabs begin to rise. Once it goes over certain tunable threshold,
the triggering mechanism flips and a time frame of data is captured and transferred to a PC
interface card and recorded to the computer's hard disk. The voltage of the inner two tabs is
buffered and directly fed to the computer through channel 1, while an in-line differentiator
performs an analog differentiation before the data is fed to the computer on channel 2.
44
Figure 2.13: The relationship between the potential difference of the inner two tabs and crack
length under constant current assumption. Each dot represents a solution for a particular
crack length and the line is just an interpolation to show the general trend of the data.
45
Figure 2.14: Data acquisition circuit for the fracture experiment. The resistance of the
wafer is represented by the resistor between A and B. With a constant current source on
the order of 1 mA, the potential difference between A and B is measured. Channel 1 picks
up the voltage reading directly while channel 2 takes the same voltage data but makes an
analog differentiation before the data is recorded.
2.2.5 Some major difficulties of this experiment
2.2.5.1 Wafer not breaking straight
There were two major problems that prevented me from getting the experiment done as it
was designed and I spent over half of my PhD life to track them down. The first problem
was that the wafer did not break straight. This may sound unlikely for the wafers were
pre-cracked in the most preferred fracture plane and subject to a uni-axial loading. Why
would it want to go any other directions other than straight down the road? But, that is
exactly what it did. Obviously, for that to happen, the opening of the middle window of
the steel frame was not purely Mode I. Where were these other modes coming from?
To answer that, I checked three possible sources of imperfection.
First, I checked the crack initiation process. The previous practice was to heat a
notched wafer in a boiling water bath and then quickly dip it into icy cold water. This
method didn't work very well. Some times it failed to generate an initial crack, while
more likely, it would produce an initial crack that is too long (such as half way through).
Besides, this method can cause the metal tabs to lose contact with the wafer. So, I chose to
use another widely used method for introducing initial crack to a wafer called three-point-
46
Figure 2.15: It is demonstrated through calculations within software package franc2d that
if only a few percent of the elastic energy goes into Mode II, the crack will deviate from its
original path, which is not desirable in our experiment.
bending. Basically, one uses a diamond pen to make a tiny scratch at some desired location
on the edge, then puts a small cylindrical object (I used a straighted paper clip) directly
underneath the scratch. The final step is to gently push down on both sides of the scratch.
This method is much more reliable and the initial crack length was partially controllable.
Second, the wafer orientation was suspected. Previously, we used the primary flat
indicated by manufactures to do our alignment. But, due to processes of polishing or
something else, the flat wasn't straight like a razor blade, especially around the corner. In
addition, the uncertainty of the flat was specified as ±1o by the manufacturer. One degree
may sound negligibly small, but after putting a shear force due to 1 degree misalignment
into a finite element fracture code called franc2d4, the code produced a badly deviated crack
(see figure 2.15 for details).
In order to solve the wafer orientation problem, I decided to create a reference plane
by myself. It turned out that all I need to do was to apply three-point-bending a few more
times. Instead of just generating a small initial crack, now I broke it all the way through.
The flat so created was amazingly straight and flat.
4Franc2d is a two dimensional, finite element based program for simulating crack propagation in
planar structures developed by researchers in Cornell University. For more information, please visit
http://www.cfg.cornell.edu/software/franc2d_casca.htm
47
Finishing the first two steps did not help too much in getting a straight break, but
did help narrowing down the problem.
Finally, I started to suspect the design of the experimental loading device. By
looking at Hauch's and Deegan's setup, I eventually realized that there is a stability issue
concerning the deformation of the steel frame between these two devices. Hauch's setup is
intrinsically stable because the forcing is external to the steel frame and self-correcting to
its symmetrical axis. In order to minimize the size of the apparatus to fit it into a small
container for the low temperature experimental requirement, Deegan's setup imposes forcing
from inside of the steel frame (by galvanizing a pair of piezo-electric transducers or inflating
a pair of bellows). That is, the self-correcting mechanism was missing, and to make things
worse, he reduced the thickness of the extension element by a factor of 5 to match the
much weaker forces provided by the piezo-electric transducer or bellows. That makes his
design very vulnerable to shear. The frame is actually unfavorable for expansion; instead
the frame will prefer to deform in the shear direction instead of the expansion direction.
(Hauch's setup is illustrated in figure 2.7 and 2.8, and Deegan's setup is illustrated in figure
2.16).
I tried many different ways to control the shear deformation of the frame. I tried
to rotate the bellows in the side windows hoping to find a balanced configuration such that
any mis-alignment of them can be canceled in the shear direction. That didn't work because
the tunable angle is very limited due to available space. Besides, for real experimental runs,
we would want to create differential loading between these two bellows such that the energy
release rate decreases as the crack propagates, so one configuration wouldn't fit them all. I
also tried to compensate the shear by measuring exactly how much shear deformation there
was and applied a side push to correct it. That worked if the setup was on the optical table,
but it is very hard to think of anyway to do this side pushing correction when the steel
frame was sealed in a container and immersed in a liquid nitrogen bath.
So, eventually I decided to abandon Deegan's design and shifted back to Hauch's
setup. Of course, it is impossible to cool the whole apparatus down to liquid nitrogen
temperature, but that is not really necessary. All that is required is to have the silicon
sample to crack under 100 Kelvin environment, so I decided to enclose the silicon sample
and the steel frame with a very soft Aluminum box and left all the forcing mechanism
48
Figure 2.16: Robert Deegan's three-hole frame design. This design would work if everything
functioned ideally, but the reality is that the world we live in is far from ideal. So, the two
piezoelectric transducers can be misaligned with respect to each other and not perfectly
perpendicular to the major axis of the frame. The frame is very vulnerable to any tiny
amount of shear force. So in practice, this design does not work.
out. Then I filled the aluminum box with liquid nitrogen and created a low temperature
environment just around the steel frame. (a) of Figure 2.17 is a real photo of my setup (red
lines outlined the aluminum box) and (b) is a schematic top view of the setup while the box
was outlined by red dashed line. Four threaded steel rods went through small holes cut on
the side of the box and connected the steel frame to the forcing bars. The holes were sealed
with plastic epoxy commonly used on fish tanks, the coupling of the box to the rest of the
system was negligible because the epoxy is very flexible and the spring constant of the box
is very small compared to that of the steel frame.
The box setup proved to work. I got some decent straight breaks with clean crack
surfaces. But, immediately I ran into another difficulty.
2.2.5.2 Wafer sliding
It seemed that the energy release rate of the modified Whomper was much higher than
Hauch used to get. The only difference was that he used glue while I was using static
49
Figure 2.17: Boxed low temperature environment surrounding Whomper, The upper sub-
figure is a real photo of the setup; the lower sub-figure is a schematic drawing of the essential
components. A stepper motor is connected to the left-most steel rod (orange), a steel bar
redistributes the load to two smaller rods (blue) which in turn act upon the steel frame.
The steel frame is still the platform for conducting the fracture experiment. A symmetrical
design refocuses all the load back to the right most steel bar (orange), which is then fixed
to a stationary point on the wall of the outer aluminum structure. Since the stepper motor
is also fixed to the outer aluminum structure on a symmetrical point of the opposite wall,
the force loop is completed. This design has the advantage of aligning itself and displays
much stronger resistance to shear forces than the design in Figure 2.16.
50
Figure 2.18: Schematic drawing of mode I crack and quantities used for calculating energy
release rate
friction to hold the wafer fixed to the frame. So, what was going wrong?
First, let me show how to calculate the energy release rate, which is exactly the
same way as Hauch did. Figure 2.18illustrates a schematic mode I break. The red arrow
points to the direction of the crack propagation and all relevant quantities for calculating
energy release rate (G) are shown.
The energy release rate for such a system can be calculated by equation 2.4.
G =
1
2
Eδ2
(1− ν2)W (2.4)
Where E is Young's modulus and ν is Poisson's ratio. W is the width of the silicon strip
that was stretched, also corresponds to the width of the middle window of the steel frame.
Finally, δ is the amount of stretch, which is directly measured by a set of inductive sensor.
For Hauch, since he glued the wafer to the frame, it was just necessary to measure
the displacement of the two edges of the middle window of the steel frame for δ.5
5Hauch actually calibrated the yield of the glue and compensated it in his calculation for δ.
51
While assuming that static friction with sufficiently strong normal force could hold
the wafer fixed to the frame as effectively as glues. I was also measuring the displacement
of the edges of the middle window of the steel frame for δ.
An unrealistically high energy release rate could only mean two possibilities:
1. The inductive sensor was not working properly.
2. The wafer was sliding while the frame was stretched. That is, the amount of extension
of wafer is less than the relative displacement of the edges.
To rule out the first possibility, I sent the inductive sensor back to the manufacturer and
asked for a complete calibration. Also, I did consistent tests between the two (three later)
sets of inductive sensors by having them measure exactly the same displacement. I even did
a laser interference experiment to calibrate the inductive sensor by myself. The inductive
sensor was working perfectly.
So, I had to conclude that the wafer was sliding even when it was well within the
static friction threshold. Eventually, I decided that I had to know how much it slides, how
fast it was sliding and what was the relation between the amount of sliding and normal
force. All of that will be detailed in part III of this thesis.
2.2.6 Differential loading to catch velocity gap
In order to determine whether there is a velocity gap, it is not enough just to have a clean
straight break, because that only gives one the crack speed as a function of energy release
rate. What is more vital is to see if a steady propagation of crack is possible as the energy
release rate gradually reduces below the Griffith's threshold. That is, one has to get a
running crack (usually at a speed of 2-3 km/s) arrested before it breaks the wafer through
(the total length of a wafer sample is 3 inch, and a seed crack is around 1 inch.). All that
matters is the manner the crack stops. If it stops gracefully (reducing speed continuously
back to zero), then there is no velocity gap. Otherwise, if the crack speed suddenly drops
to zero from some substantially non-zero value, then there is a forbidden range of velocity,
hence the existence of a velocity gap.
In principal, the velocity gap should also be observable by examining the way the
crack starts. But, the existence of lattice trapping and the fact that the seed crack is not
52
long enough to maximize the stress intensity factor makes it improper to argue that a rapid
rising of crack speed at the beginning of the propagation is due to the existence of a velocity
gap. That is, a crack will always start to propagate only when the energy release rate right
ahead of it is substantially more than the minimum to support a crack to propagate.
The experimental method of creating a gradually decreasing energy release rate is
to have a differential loading, that is, to stretch the wafer harder on one end (the end where
the seed crack was implanted) and less on the other end. It was relatively easy for Deegan
to create and manage a differential loading (refer to his design in figure 2.16) because a
gas regulating system was built to control the gas pressure of individual bellows, hence the
force exerted by each bellow can be separately controlled. In order to realize differential
loading in my setup, I inserted a weak bellow into one end of the rectangular hole of the
steel frame (figure 2.19 provides a close up look of the bellow in action). While conducting
the experiment, I first drove the motor to stretch the wafer to about the Griffith's threshold,
and then slowly inflated the small bellow until the wafer cracked.
2.3 Data Analysis and Discussion
After all the efforts of troubleshooting and perfecting every step of the experiment, I got a
success rate of 1 out of 10. (Not that bad considering that the previous overall record was
less than 1 out of 100.) So, before I show the only sets of good data, let me show how a
typical bad set of data looks.
2.3.1 Data of partial cracks
For both room temperature and low temperature fracture experiment, most commonly seen
are those of partial cracks as shown in figure 2.20. Partial crack is a term used in this
thesis to stand for a particular fracture experiment where the crack propagates in a jerky
fashion that it starts and stops multiple times and usually never reaches speed higher than
1 km/s. For those partial cracks, there are features on the crack surfaces resemble those
reported by Ilan Beery et al[36], and indeed those cracks ran at a speed of a few tens to
a few hundreds of meters per second as they had discovered, but that doesn't mean that
there isn't a velocity gap, because those cracks never ran at a steady state. Whether it is
53
Figure 2.19: Bellow for differential loading
54
Figure 2.20: Length and velocity records of cracks moving in jerky fashion (partial cracks).
Note that the velocity never exceeds 1 km/sec and never reaches a steady state. These ex-
perimental signals are characteristic of misaligned samples that leave behind rough fracture
surfaces.
actually another manifestation of velocity gap remains to be carefully examined.
2.3.2 Data of a good room temperature crack
Figure 2.21 shows a typical data set for a good room temperature crack of silicon, which
contains only one propagation process and has stabilized voltage reading before and after
it. The upper panel is the crack tip position as a function of time and the lower panel
is the corresponding velocity profile of the crack. In order to produce this, the resistance
of the crack is measured before and after the crack. That information combined with the
resistivity of the wafer are used to construct a similar curve as shown in figure 2.13, which
serves as the conversion table to translate the voltage readings to crack tip positions. The
velocity profile is obtained by doing an analog differentiation on the same voltage data and
applying the same conversion rules.
55
Figure 2.21: A good room temperature break; the crack reaches a velocity over 1 km/sec
and holds it for several microseconds. The distance is calibrated by noting the final position
of the crack tip after the experiment is completed.
56
Figure 2.22: A good low temperature break.
2.3.3 Data of a good low temperature crack
Figure 2.22 shows a good data set of low temperature crack. Even though the silicon sample
is exactly 7.62 cm long and it seems that this crack runs dangerously close to that limit, it
is obvious that this crack is arrested because otherwise the voltage reading will shoot over
limit and the data will present a huge vertical jump instead of the flat plateau seen in this
case.
57
Figure 2.23: Comparison of the arrests of cracks at room temperature and liquid nitrogen
temperature under comparable energy release rates reveals the existence of a velocity gap.
2.3.4 Velocity gap confirmed
Comparing the manner the crack was arrested at room temperature and liquid nitrogen
temperature (figure 2.23) provides the best indication yet obtained that a velocity gap does
exist. For the arrest at room temperature, it takes about 4 µs for the crack to descend
from 0.8 km/s to zero. While arrest at low temperature takes only about 1 µs and it jumps
from above 3 km/s back to zero. The gradual arrest at room temperature is attributed
to thermal fluctuation effects and the abrupt arrest at low temperature is believed to be a
strong evidence for the existence of a velocity gap.
2.4 Conclusion
After several years of struggle, I managed to obtain two sets of good data that indicate the
existence of a velocity gap. In combination with a small number of runs previously obtained
58
by Deegan, they constitute the best evidence to date that low-temperature fracture of silicon
differs from room-temperature fracture. All low-temperature breaks arrested more quickly
than all room-temperature breaks. Unfortunately, the small number of cases where the
phenomenon has been seen and doubts cast by the difficulty of the experiment have left
us reluctant to publish. In addition, sliding between frame and sample, which will be
documented in the next Chapter, left us without clear knowledge of the precise value of the
energy release rate in those few good runs we obtained. Thus despite the effort put into this
experiment, I have been forced to leave it still incomplete.
59
Chapter 3
Friction
3.1 Introduction
The purpose of this Chapter is to discuss a set of experiments concerning frictional mo-
tion over very small distances. As mentioned in the introduction, the original purpose of
the friction measurements was to calibrate the energy release rate for our low temperature
experiment which used friction as the gripping mechanism. But, as the attempts to char-
acterize friction proceeded, they became valuable findings on their own. A conventional
picture of friction that most people have in mind is illustrated in figure 3.1. For a small
enough applied shear force, static friction force is always equal and opposite of it and hence
no motion will occur. As the shear force increases to a value that exceeds the limit of static
friction, motion begins. The maximum static friction is determined by the static coefficient
of friction multiplied by the normal force (µsFn). Due to the fact that the dynamic co-
efficient of friction is usually smaller than the static coefficient of friction, the motion will
accelerate a bit at the beginning. But, then the dynamic coefficient of friction increase as
the relative velocity increases. The interplay of these two factors can produce a slip-stick
motion. One typical example was shown in figure 3.2[1].
Few people have ever asked why there is always an equal and opposite static friction
force when a small shear force is applied. Or, what is the origin of a static friction force?
Some who did ask found themselves facing a dilemma that there shouldn't be any static
60
Figure 3.1: Common understanding of friction (The lower portion of this graph is from
World Wide Web.)
61
Figure 3.2: Slip-stick motion (This graph is adapted from Cochard's paper[1].)
friction. The argument goes like this:
Most macroscopically flat surfaces look like endless hills and valleys on a microscopic
length-scale. When two such surfaces come in contact, only a small portion of the surface
is touching the other side as illustrated in figure 3.3. Based on this picture, it is temping to
think of the origin of the static friction as purely geometrical. That is, one needs to provide
energy to push one hill to climb over the other. But, the energy is recovered once it passes
the highest point of the other hill and starts going downward. Statistically speaking, at any
time there will be equal number of hills to climb over one another compared to the hills
going down the valley. So, in the macroscopic limit there shouldn't be any energy required
to move one surface against another, because there is no energy lost in this picture.
So, does geometry matter at all? I believe it does. In reality, there aren't just hills
and valleys, there are also rocks and rivers. That is, there are all kinds of contamination on
most surfaces. Robbins[91] et al. argued that if there are particles in between two surfaces
moving against each other, then there will be a net energy loss and hence friction. But,
62
Figure 3.3: Rough surfaces in contact as a conventional picture of friction
again you will need to have the surfaces in relative motion to have friction, which can not
explain the origin of the static friction.
A more physical argument for the origin of static friction will have to consider the
dynamics on the microscopic level. Once two hills are pressed into each other such as the
one circled in red line in figure 3.3, atoms on the surface of one object get close enough
to feel the Van der Waals force due to the atoms on the surface of the other object. The
detachment of the Van der Waals bonding could cause atoms to vibrate more violently and
dissipate energy through mechanical (sound), thermal (heat) and even electrical (ripping
electrons from one material to another) means. This still sounds like a physical argument
for dynamic friction, and indeed it is. But, now I will argue that it also helps explain the
origin of static friction (or just friction in a more unified sense).
When two surfaces just come into contact, a population of asperities has been born.
There will be a distribution in sizes of asperities as well as a distribution of shear strengths.
And that population will grow due to slow plastic flow, which is known as aging.1
But, as soon as a shear force is applied, the population of asperities will start to
reflect its existence. The two surfaces will move against each other a little bit to create more
asperities and certainly some asperities begin to break. It really is a dynamic process which
happens on a microscopic scale. If the net gain of asperities is enough to counter-balance the
shear force, then the system will be stable and will not show visible motion, which is then
still categorized as static friction. The system will be unstable only when the shear force
is not counter-balanced by asperities growing. Actually, when the relative motion is quick
enough the creation of asperities will be less than the breaking of asperities, which is also
why dynamic coefficient of friction is usually smaller than the static coefficient of friction.
Once that happens, macroscopic motion starts. Depending on the surface conditions such
1Maybe just enlarging some of the real contact areas, but in the sense of an average size of asperities it
is equivalent. Basically, it is the increasing of the real contact area.
63
Figure 3.4: Inclined plane experiment between a polished silicon wafer and a machine-ground
steel surface.
as flatness, lubrication, and contamination, it could end up smooth sliding motion or a
slip-stick motion.
To make a unified approach to the origin of friction, friction can be understood as
a dynamic process of continuously forming and breaking of tiny inter-facial cracks on the
microscopic length-scale.
3.2 The inclined plane experiment
As the starting point to understand why the silicon samples of our low-temperature ex-
periment were sliding against the steel grips intended to hold it in place, , I performed an
inclined plane experiment to measure the conventional coefficient of static friction.
The setup of this experiment is very simple (figure 3.4). I put a wafer on the steel
frame and slowly increased the incline angle α. The surfaces in contact are exactly the same
pair as those I put to measure how much the silicon sample slides. The silicon sample was
doubly polished (industrial standard) and the steel frame was machine ground to a surface
flatness of less than a micron.2
It is clear that the static friction coefficient can be found by the following equation
2I also conducted the inclined plane experiment for a rough steel surface (sand papered, used type-400.),
which only raises the static coefficient of friction by about 20%.
64
µs = tan(α) (3.1)
where α is the maximum incline angle before the wafer begin to show macroscopic
slipping. 3
The static coefficient of friction so found is 0.18 ±0.014 (The experiment was re-
peated 12 (N) times, and the reported value is the mathematical average and the uncertainty
of it which is the standard deviation divided by N
1
2 ).
3.3 The friction experiment
Additional friction experiments were performed in exactly the same manner as for the frac-
ture experiment except that the wafer was not pre-cracked and two pairs of inductive sensors
were glued to the wafer as shown in figure 3.5 to provide careful measurements of sliding
between the sample and the grips.
During the experiment, a pair of normal forces N was applied through the pressure
blocks to press the silicon sample onto the steel frame. The steel frame was then stretched
by a stepper motor. Sext represented how much the stepper motor retracted and was also
directly measured as the control parameter of the experiment. Sf represented how much
the edges of the middle window of the steel frame was displaced. Additionally, the sample
extension S was directly measured by the inductive sensors. If conventional understanding
of static friction were correct, then S would equal Sf as Sext began to increase until it
reached a point that the internal stress of the silicon sample provided a shear force greater
than 2µN .4 Once that point was reached S would stay a constant value while Sext and Sf
could keep increasing. However, my experimental results indicated that this simple picture
was incorrect. Especially since slipping started at the very beginning as Sext started to
increase I concluded that a static coefficient of friction µ0 might not exist at all; instead
static friction has a dynamic origin.
3In order to get a reasonably distributed result, I glued the silicon wafer to a metal block. Otherwise,
the weight of the wafer is so small and its thickness so thin compare to the other two dimensions that the
result scatters over a very wide range.
4The reason for a factor of 2 was that there were two interfaces (both sides of the silicon sample were in
touch with the same type of steel) contributing static friction.
65
Figure 3.5: Schematic diagram of the experimental setup of friction. The four rods in blue
color are the aluminum holders for the two pairs of inductive sensors, which are machined
to exactly the same specification.
66
The reason for using two pairs of inductive sensors instead of one was that the
wafer bends slightly when a normal force is applied. This bending will then be magnified
(due to the long arm of the holder of the aluminum target) and picked up by the inductive
sensor. Obviously, the symmetrical setting of the inductive sensors can be used to eliminate
the bending effect (by taking the average of the reading from both sensors). Figure 3.6
shows how much the bending was picked up by individual inductive sensor and how good
the cancellation was. The reason for such a good cancellation is that the two sensors have
exactly the same height as measured from the surface of the wafer (machining accuracy
limit 1/1000 of an inch) and close (see footnote) to exactly the same distance apart. This
geometry guarantees that a small deflection due to the bending of the wafer creates equal
amount of contraction or elongation of the distance between the inductive coils and their
perspective targets. So, the summation of the reading of the two sensors is very close to
zero as indicated by the blow-up of figure 3.6. (The two humps in the data was due to the
inflation and deflation of the gripping bellows.) 5
3.4 Modified rate- and state- friction law
3.4.1 Rate- and state- equation and its modification
The rate- and state- friction law was first proposed by Dieterich[15] in the 1970s based on
experimental findings. The dynamics between two contacting surfaces in the microscopic
length-scale must be very complicated as it involves forming and breaking of asperities
(micro inter-facial cracks) and unavoidably some surface contamination, such as finger oils,
water molecules, dust and debris from the broken asperities etc. But a surprisingly simple
model with only two macroscopic parameters can describe the motion rather well. The two
parameters are the relative velocity and a state variable that evolves as the surfaces slide.
In general, the friction coefficient µ can be understood as the ratio of horizontal
force F to normal force N,
µ =
F
N
= A ln(1 +
v
v∗
) +B ln(1 +
θ
θ∗
) (3.2)
5An extra factor of 1.0945 was multiply to the blue curve to reflect a small distance mismatch between
holders.
67
Figure 3.6: The red and gray thin curves were individual readings from each inductive sensor
and they varied a lot during the process of inflating or deflating the bellows, but their sum
stayed very close to zero as more clearly indicated by the sub-figure, which is a magnified
view of the sum of these two curves.
68
Here A, B, v∗and θ∗are constants, v is the relative velocity of sliding, and θ is
a state variable that carries all the information about the population of asperities. This
expression differs slightly from the original version because a static coefficient of friction µ0
is missing. Based on experimental findings , the static friction is an ill defined term and
actually comes from the dynamics of the relative sliding. That is, if the state variable and
relative velocity vanish, the friction force vanishes as well. This choice is motivated by our
inability in experiment to observe any lower threshold of horizontal force below which there
is not at first some sliding, which can be most clearly seen in the continuous data (figure
3.7) where the slope of the data (blue curve) is never as steep as the green curve (were there
a non-zero static coefficient of friction). The green curve is what one will expect assuming
the conventional understanding of static friction is correct. That is, the sample extension S
will follow exactly as how much the middle gap of the steel frame is displaced until it reaches
a point that the internal stress of the silicon sample is greater than the threshold of static
friction. Once that point is reached, S will become a constant even though the steel frame is
stretched further. The deviation of the real data (blue curve) to the conventional prediction
(green curve) at the very beginning of the curve (Data close to the origin is when internal
stress of the sample would be very small.) suggests the absence of the static coefficient of
friction. Etsion and Lee and Polycarpou[11] have reported the same phenomenon.
To complete the theory, it is necessary to propose a dynamic equation for the state
variable . There are two common expressions for it.
A first form is from Dieterich[15], and is
dθ
dt
= 1− v
Dc
θ (3.3)
Where, Dc is a constant that describes a sliding length over which the asperity
population reaches a statistical steady state. Another was discussed by Ruina[77] and is
dθ
dt
= − v
Dc
θln(
vθ
Dc
) (3.4)
These two expressions behave differently in the static case where v = 0. The first,
Eq. 3.3, predicts aging, which means that the state variable increases linearly in time if the
system is stationary, and the friction coefficient increases logarithmically. The second, Eq.
69
Figure 3.7: One set of continuous data under the normal force of 120 N. The green curve
is the prediction of conventional friction, the black curve is the prediction of the modified
rate- and state- friction law treating Dc as a material constant while the red curve is the
prediction of the same modified rate- and state- friction law but allowing Dc to vary as a
function of normal force.
70
Figure 3.8: I sat the silicon sample on top of the steel frame and waited for four differ-
ent lengths of time before conducting exactly the same experiment. The results bear no
indication of aging except for a slight difference between data obtained after waiting for
30 minutes and the rest of them which were obtained after waiting for much longer time.
All the rest of my friction experiment waited for more than 24 hours, so aging should not
present itself in any of my other data sets.
3.4, predicts that when the sample is not sliding the state variable does not evolve.
There is compelling experimental evidence dating back to Coulomb that the strength
of static frictional contacts increases over time[14]. However, this evidence was obtained in
specific experimental systems such as wood on wood or rock on rock. Berthoud et al. (17)
show that aging rates of polymeric systems drop rapidly as a function of temperature below
the glass transition temperature. Thus for a system such as ours, with large yield stresses
and far from any melting temperature, one might expect aging effects to be very small. We
show in Fig. 3.8 that the experimental results of silicon on steel we have studied give no
evidence for aging on a time scale ranging from half an hour to 2 weeks.
Because of the absence of aging in our experiments, we have looked for a generaliza-
tion of both Dieterich and Ruina's equations. We suppose that there are two state variables
θ and φ , where φ is proportional to the aging time of contacts, and θ describes how the
71
contacts evolve during sliding.
Equations implementing this idea are
dθ
dt
=
v
Dc
(
1 + φ
1 + v/vθ
− θ) (3.5)
dφ
dt
= αN/τ − vφ
Dc
(3.6)
where N is the normal force, τ is a yield stress, and αis dimensionless.
Ruina similarly presents equations for a pair of state variables. To understand the
behavior of these equations, it is useful to consider a situation where the horizontal force
increases slowly from zero and then holds at a constant final value F. For almost all initial
conditions, the sample initially begins to slide with nonzero velocity v and then heads
toward solutions where v and θ are time-independent. The time-independent solutions are
of two forms. If F is below a critical threshold, the final solution has v = 0 and is given
by F = Bln(1 + θ/θ∗). This solution corresponds to a system gripped by static friction.
However, these solutions are only stable so long as θ is less than a critical value of 1 + φ.
When F/N is larger than Bln(1 + (1 + φ)/θ∗), the steady solutions instead have nonzero v
and θ = (1 + φ)/(1 + v/vs). This situation corresponds to sliding friction.
The critical shear force that divides these two regimes is
µs =
F
N
= Bln(1 + (1 + φ)/θ∗) (3.7)
which we identify with the static coefficient of friction. If instead one waits for a
long time (T ) before beginning to ramp the force up from zero, φ increases by αNT/τ , the
critical shear force increases by the log of this value, and thus the equations can describe
the phenomenon of aging.
In view of the fact that our experiments show no aging, we simplify to
dθ
dt
=
v
Dc
(
1
1 + v/vθ
− θ) (3.8)
and the expression for the static coefficient of friction simplifies to
72
µs =
F
N
= Bln(1 + 1/θ∗) (3.9)
The transition to sliding at this critical force value can be either sub-critical or
super-critical depending on the values of A and B.
3.5 Using the modified rate- and state- equation to fit
our data
We now apply the modified rate- and state- friction equation to fit our data.
First, we write down equations correspond to the experimental geometry shown
in lower portion of Figure 3.5. We assume that the experiment is symmetrical about the
horizontal midpoint and focus on the right-hand side. The amount the sample has slipped
isx ≡ Sf − S. Denote the force of friction by µN . In our experiments, s > 0, and friction
pulls the right-hand side of the sample to the right. Assuming that none of these quantities
changes sign, and neglecting inertia, mechanical equilibrium requires
ksS = 2µN (3.10)
(Sext − Sf )kext = 2µN + Sfkf (3.11)
To simplify the equations further, define
Cext =
kskext
ks + kext + kf
(3.12)
Cx =
ks(kf + kext)
ks + kext + kf
(3.13)
Thus, the two mechanical equilibrium equation can be transformed as
2µN = CextSext − Cxx (3.14)
73
Figure 3.9: Step data; The normal force for this data set is 180 N and it was held fixed
throughout the experiment.
coupling this to the dynamic equation of θ (eq. 3.8), proposed functional form of µ
(eq. 3.2) and recognizing the definition of relative sliding velocity as dxdt =
d(Sf−S)
dt = v give
us a complete set of equation to describe the sliding motion between two surfaces.
In the limit when one increases Sext very slowly and measures S, which was the
condition for the continuous data, a closed form solution for Sext can be obtained as
Sext =
ksS
Cext
− CxDc
Cext
ln(1 + θ∗(1− ekss/2BN )) (3.15)
So, we first fit the continuous data with the quasi-static solution given by eq. 3.15.
The result was shown in figure 3.12. with the spring constants measured independently (The
spring constant of silicon was calculated from its Young's modulus and Poisson's ratio). The
best fits to these experiments set the values of Dc and θ∗ , and B is determined from Eq.
3.9.
To obtain the rest of the fitting parameters, we have relied upon direct integration
74
Figure 3.10: Slipping of the wafer as the boundary condition was fixed
of the dynamic equations. Figure 3.9 compares experimental results in which one repeat-
edly stretches the sample and waits a fixed period of time. The parameters (A and v)are
determined by searching for the best fit. When slipping velocities v are much larger than
the cutoff velocity v*, the dynamic equation has solutions where the slip x is logarithmic
in time, which was confirmed experimentally as shown in Figure 3.10 (which is a magnified
step in figure 3.9), this logarithmic slip is observed until displacements become too small
to measure accurately. Thus, we cannot really determine the cutoff velocity v* but only say
that it is less than 10−5 µm/s.
Figure 3.11 indicates the physical quantities that I measured for the friction ex-
periment. ks is the spring constant of the sample (calculated from Young's modulus and
Poisson ratio), kf is the spring constant of the frame (measured independently). and kext
is the overall compliance of the rest of the system (measured independently). Sample ex-
tension S is how much the sample was stretched, Sext is how much the overall system was
stretched, and Sf is how much the middle window was displaced. Sext can be viewed as the
control variable and S can be taken as the response. The normal force N became a tunable
75
Figure 3.11: Illustration of the quantities measured for the friction experiment
parameter. If N is zero, the wafer should not be stretched at all, because the only force that
stretches the wafer was the friction force. On the other hand, if N is very large S should
approach Sf .
I used two slightly different ways to conduct the experiment:
3.5.1 Step Data
In the first, I ramped up Sext as quickly as the machine was capable of responding and then
stopped for a fixed period of time. The distance between the edges of the middle window
was fixed during that time as confirmed independently while measuring the compliance of
the frame.
A typical step data set is shown in figure 3.9. The normal force for this set of data
was 467 N. The theoretical fit was based on the modified rate- and state- friction law.
A simple calculation can demonstrate that the conventional understanding of friction
will fail to predict the behavior shown in figure 3.9.
76
It is well known that the Young's modulus of single crystal silicon in the [111]
direction is 168.9 GPa. So, the spring constant of the wafer sample is
ks = E[111] ∗ A
W
= 168.9GPa ∗ 3.0inch ∗ 380µm
1.3inch
= 1.48e8
N
m
(3.16)
For an extension of 0.2µm, the shear force will be
ks ∗ S = 1.48e8 N
m
∗ 0.2µm ≈ 29.6 N (3.17)
The static threshold is supposed to be
2 ∗ µ ∗ FN = 2 ∗ 0.18 ∗ 467 ≈ 168 N (3.18)
That is, the conventional understanding of static friction will predict that the wafer
shouldn't slide at all (even at the maximum stretch of the experiment), while the experiment
showed quite the opposite, and the continuous data will further illustrate that the sliding
happened almost immediately after the frame was stretched at which moment the shear
force due to the extension of the wafer was nearly negligible.
3.5.2 Continuous data
To further validate the proposed rate- and state- friction law, I also conducted a series of
experiment where the change of Sext was continuous and fairly slow (5 times slower than the
jump of the step data). Figure 3.12 shows the data and three different ways of fitting. The
green curve gives the prediction of the conventional understanding of friction. The black
and red curve were based on the proposed friction law. The difference is that the black
curve held Dc as a constant while the red curve allowed Dc to be proportional to the normal
force. Dc is one of the fitting parameters in the dynamic equation for the state variable.
(Please refer to equation 3.4.) The physical meaning of Dc is believed to be a characteristic
length-scale of the asperities. It is not clear why allowing Dc to be proportional to normal
force gives a slightly better fit.
77
Figure 3.12: Continuous data from experiments where the normal force was constant and
shear force increases slowly and linearly. The green curve shows the sample response that
would be expected from the conventional theory of static friction. The sample would stretch
with the frame up to the limit of static friction and then slide freely. Instead, as shown by
the experimental data in blue, the sample always moves less than predicted. I provide two
fits to the data based on the modified rate and state friction law. The first of them (black)
treats the asperity length Dc as a constant. The second of them (red) allows Dc to vary
linearly with normal force. All parameters of the model are held constant across all six
panels.
78
3.6 Conclusion
For hundreds of years, there have been certain basic facts about friction that have not been
questioned. One of them is that for any given level of normal force, there is threshold in
shear force below which no motion is possible. Once this threshold is reached, the solid will
begin to slide, and dynamic friction takes over. Even specialists in the modern theories of
friction have retained the concept of static friction in their models.
With the experiments reported in this Chapter, I have shown that the distinction
between static and dynamic friction is not well defined, and only seems natural when mea-
surements of motion are not sufficiently accurate. When motions are measured down to
the nanometer scale, objects supposedly held still by static friction are not really static
at all. They move in predictable and reproducible ways. Furthermore, the mathematical
laws that govern their motion appear to be identical to those now thought to describe dy-
namics friction. Some further work here is called for, since we have not yet shown that
the parameters in the mathematical models we use simultaneously describe the very slow
motions investigated here and the more rapid motions traditionally used to obtain rate and
state theories of friction. However it seems we have obtained a unified picture of static and
dynamic friction, one where the different types of motion are due to different force levels
and loading histories. This work has been published in PNAS[93].
79
Chapter 4
Granular Simulation
Experimental study of the internal structures of a granular system is difficult. Nevertheless,
Behringer[6, 89] etc. used the birefringent pattern of photo-elastic disks to study the force
chain distribution in a 2-D system. Seidler[21, 72, 84] etc. applied X-ray micro-tomography
to measure the position of every granular particle in a 3-D packing. Those experimental
techniques are great for studying static properties, but are difficult to implement for studying
dynamic properties. There is an experimental setup that is designed specifically for dynamic
purposes, which is a huge 2-D air-table run by Troadec[40] etc., but their main focus is to
study the motions of particles in the gas phase.
MD simulations of granular system have become a complementary tool for scien-
tists to study detailed particle motions at the microscopic length-scale in order to better
understand some of the intriguing properties of a granular system.
4.1 Preparation of a confined granular system
4.1.1 LAMMPS' implementation of the force-displacement relation
The version of LAMMPS I used was last updated on Sep-28-2008. In that version of
LAMMPS, the force-displacement relation between two granular particles was not correctly
implemented in the sense that the calculated result did not match the equations given by
the manual. So, I modified the source to implement the Hertz contact force correctly and
80
ran all my simulations on it. (Please refer to appendix 6.2 for details of the modification.)
After the modification, the interaction between two balls are described by the fol-
lowing formula:
−−−−→
FHertz =
√
d− r
√
RiRj
Ri +Rj
[(kn(d− r)−→nij −meffγn−→Vn)− (kt(d− r)∆−→st +meffγt−→Vt)] (4.1)
Where d = R1 +R2 is the contact distance between two balls of radius R1and R2
r is the distance between the centers of contacting balls, r =
√
(x1 − x2)2 + (y1 − y2)2
kn is the elastic constant for normal contact, kt = 27kn (LAMMPS convention)
γn is the viscoelastic constants for normal contact, γt = 12γn (LAMMPS convention)
meff = mimj/(mi +mj) is the effective mass of two balls of mass mi and mj
Vn is the normal component of the relative velocity, and Vt is the tangential com-
ponent of the relative velocity
−→nij is the unit vector between the centers of two balls
4−→st is the tangential displacement vector (truncated to satisfy a frictional yield
criterion)
Obviously, the first term in the square bracket is the normal force and the second
term is the tangential force. The implementation of the tangential force is still not very
physical. It simply relates the tangential spring constant and tangential viscoelastic con-
stant to their respective normal direction constants by two arbitrary ratios 27 and
1
2 (This
limitation has been lifted in the newer version). This convention has been commonly used
in literature, but a more physical implementation, such as that proposed by Mindlin[69] and
Dereciewicz[68] should be implemented to generate more accurate results. Due to the time
constraint, I didn't implement Mindlin tangential force to my modified code, instead I stuck
to LAMMPS convention for the tangential force calculation.
4.1.2 packing generation and boundary conditions
To minimize the complexity and maximize efficiency LAMMPS left out all the pre- and
post- processing to supplemental codes and software. So, I used another program called
81
PFC2D.1
The initial configuration of the system was a 2-D granular gas. A loose granular
packing containing 5,000 balls in a 10 cm by 10 cm 2-D simulation box was ported from
PFC2D with volume fraction about 0.78 (The packing is not confined yet.). The horizontal
direction is defined as the X-axis and the veritical direction is defined as the Y-axis. It
was a poly-disperse granular system with ball diameters evenly distributed between 0.06
cm to 0.08 cm. The material properties of the constitutive grans were assigned based on
the properties of common glasses, which has Young's modulus E = 72 GPa and Poisson's
ratio γ = 0.1. The initial temperature of the system was zero (all velocity components
of all balls were set to zero). The coefficients of friction in ball-ball interaction and ball-
wall interaction were all set to 0.1 from the beginning of the simulation. To confine the
system, the top wall was moved to the negative Y direction gently using a half period
cosine profile, y = Acos( 2piT t); A = −0.1cm T = 0.08s while the bottom wall was fixed
throughout the simulation. A periodic boundary condition was applied to the horizontal
direction. That is, if a ball moves out from the left boundary of the simulation box, it
will re-enter the simulation box from the right at the same height. The integration time-
step for the underlying Newtonian dynamics is set to 2.0e-8 second, which is smaller than
the time for a pulse to travel through the minimum ball radius. Tests using smaller time-
steps show no significant difference in dynamics. That is, 2.0e-8 second is a small enough
time-step to have a converging numerical result. This downward movement of the top wall
was repeated several times to reach a desired final location. Then the system was relaxed
while the boundaries were all held fixed. A stable reading of force acting on the bottom
wall was reached as the system quickly dissipate its kinetic energy. After that, the system
was subject to 100 cycles of large amplitude high frequency (amplitude as large as 1e − 4
cm and frequency as high as 1 MHz) vibration from the top wall to further ensure that the
micro-configuration was stable and wouldn't collapse if pulses or sine waves were to transmit
through the system later. Finally, the system was relaxed2 again for much longer time such
that virtually all kinetic energy was damped out. The system was then ready to be studied.
1PFC2D is itself a fully functional commercial MD code. I was able to obtain a demo version of it and
ran many of the prove-of-concept test script on it. For more information on PFC2D, please visit the vendor's
website : http://www.winternet.com/~icg/pfc.html.
2All relaxation was done while holding all boundary conditions fixed.
82
4.2 Static properties of the prepared system
A snapshot of the prepared system and its contact force network3 is illustrated in figure
4.1. One obvious characteristic of the contact force network is the spacial heterogeneity.
To further quantify the contact forces, a common practice was to plot the probability den-
sity function (PDF). Figure 4.2 was the corresponding probability density function for the
contact forces shown in figure 4.1. It is generally agreed that the PDF should have a peak
or plateau for forces less than the average and an exponential tail for forces greater than
average[66, 74]. Figure 4.2 indicated that my system fit that character rather well especially
for the inset, which was a semi-log plot of the same distribution curve with a linear fit for
forces that were greater than average.
The contact force network has been extensively studied. Hecke[89] argued that there
is a subtle change in the tail of the distribution curve that determined the onset of jamming.
4.3 Dynamic properties of the prepared system
4.3.1 Pulse mode
One of the easiest way to probe a granular system was to send a pulse through it and measure
quantities such as the pulse speed, the attenuation rate etc. The pulse I sent through the
confined granular system was a one-period cosine wave as illustrated in figure 4.3. The
method of generating such a pulse was to move the top wall exactly as the plotted cosine
pulse and held it fixed before and after the pulse.
Figure 4.4 was a snapshot of a pulse propagating through the system. The horizontal
direction was defined as x-axis and the vertical direction was defined as y-axis. The color
code corresponding to the magnitude of the y-component of the particle velocity. Floaters
were usually moving at much higher speed comparing to those jammed particles. That is
why their color was out of the spectrum.
3Only normal forces were plotted and discussed here.
83
Figure 4.1: A snapshot of the system and its contact force network
84
Figure 4.2: Probability distribution curve
Figure 4.3: Shape of the cosine pulse
85
Figure 4.4: A snapshot of a pulse propagating through the system, What point in time?
86
Figure 4.5: Horizontally averaged y-component velocities at equal time interval (colors to
indicate different instants of time.)
4.3.1.1 speed of a pulse
To measure the speed of a pulse, it is necessary to define the arrival of the pulse. Some
common choices were the first raise of 10% of peak magnitude, the peak itself or the first
zero after the peak[18]. Depending on the choices, the measured pulse speed could differ by
10% or more.
I used the peak as the arrival of the pulse and estimated the speed of the pulse in
figure 4.5 to be 1.0e5 cm/s. Because the peak moved about 5 cm between the first peak and
the second last peak and the time spacing between them was 5.0e-5 second (1.0e-5 second
87
interval). This speed was then used to make an estimation of the base resonant frequency
(the first harmonic) utilizing the following formula
f =
v
2L
(4.2)
where v is the speed of the pulse and L is the length of the system.
So, this would predict a resonance peak at 1.0e5/(2 ∗ 10) = 5 kHz, which is a bit
lower than it actually is (~4.96 kHz).
The behavior of pulses propagating in a granular system have been extensively
studied[92, 25, 18]. It was generally accepted that the speed of a pulse in a granular media
is proportional to the power of confining pressure:
v ∝ pα (4.3)
and the power coefficient α is somewhere between 14 and
1
6 .
4.3.1.2 attenuation and total absorbing layer
A puzzling experimental finding in geophysics community is that the attenuation rate is
linearly proportional to frequency over a wide range (as much as 9 decades)[58, 83, 49, 10, 62].
Commonly used dissipation models are unable to explain it. In the quest for a satisfactory
dissipative model, Marder finds two possible candidates. One is a hysteresis model and the
other is a model of array of energy sinks. He believes that both models work in some degree
depending on the physical situation to be modeled. And, he also proposes a method to
distinguish one model from another.
A simple discrete one-dimensional version of a hysteresis model can be described by
the following equation[53]:
u¨i = F (Ei, E˙i)− F (Ei−1, E˙i−1) (4.4)
Here
88
Figure 4.6: Hysteresis loop
Figure 4.7: sine wave transforms to triangular wave
Ei = ui+1 − ui (4.5)
is the extension of the bond between particle i and i+ 1, and F (Ei, E˙i) is the force between
them.
To create hysteresis, the forces between neighboring particles assume the following
form
F (E, E˙) = E + b|E| E˙|E˙| (4.6)
A hysteresis loop of this form of forces is shown in figure 4.6.
A numerical solution of this model (figure ) shows a propagating sine wave transforms
to triangular wave within 20 wavelengths.
89
Another competing model is the arrays of energy sinks model[13]. It describes the
interaction between any two mass points by a large number of over-damped dissipative
modes. The solution of this model shows no sign of change of waveforms.
So in order to find out which one of them gives a better description of the energy
dissipation mechanism in the granular packing that I constructed numerically, all I need to
do is to send a sine wave through the packing and look for changes in the waveform (if any).
Preliminary results indicate that the wave doesn't change its shape even after hun-
dreds of periods, which is clearly in favor of the array of energy sinks model.
One prerequisite of this study is to construct a total absorbing layer in order to
mimic an infinite medium that an outgoing wave never reflects back. Walker and Wax[90]
showed in 1940s that an exponential impedance could provided a total damping to transverse
acoustic waves in isotropic non-homogeneous medium. Khilla[48] showed in 1980s that an
exponential impedance could provided a total damping in a transmission line. But, a total
absorbing layer has never been attempted in a granular simulation.
The current release of LAMMPS is unable to assign a spatially varying damping
coefficient, but it is not that hard to insert a few blocks of code to make that happen.
The modified LAMMPS evaluates a viscous force based on where the contact locates and
by setting an exponentially increasing viscous coefficient in certain region one gets a total
absorbing layer.
Figure 4.8 and figure 4.9 demonstrates how well the total absorbing layer works.
The absorbing layer is between 0 to 50 mm with exponentially increase viscous coefficient
from 50 mm downwards. A pulse is propagating from the top to the bottom (from above
100 mm to 0 mm). In the system with a total absorbing layer, there is almost no reflected
wave. While in the system without a total absorbing layer, a reflected wave exists. (Dashed
line indicates reflected wave.)
4.3.2 Resonance mode
Paul A. Johnson and his coworkers argued that there was a new class of materials called
Nonlinear Mesoscopic Elastic (NME) , whose behavior could not be fully appreciated by
classical nonlinear theory (Landau nonlinearity)[31]. The range of this new class was fairly
wide, including sand, rocks, soils, cement, concrete, artificial granular packings, and dam-
90
Figure 4.8: With total absorbing layer (colors to indicate different instants of time.)
91
Figure 4.9: Without total absorbing layer (colors to indicate different instants of time.)
92
aged materials. The unique characteristics of them included hysteresis, end-point memory,
harmonic generation and resonant frequency shift etc.
Resonant frequency shift was a well known phenomena in the geophysics community,
which is the downward shift of resonant peak as driving amplitudes increase also known
as dynamic softening. But, it was not very clear what was the underlying mechanism.
Experimental investigation was difficult due to the fact that it was a differential effect. In
the case of 2-D photo-elastic disks, the contact force network might be relatively easy to
figure out compared to finding out how some tiny extra load was re-distributed. Besides, it
is dynamical and only significant in statistical measure.
MD simulation of a confined granular system provided a powerful alternative and
much more detailed information regarding the microscopic dynamics, whose collective effect
made the resonance curve shift on the macroscopic scale.
Using the estimated resonance frequency from time-of-flight measurement as a refer-
ence point and driving the top wall with a single harmonic motion at frequencies close to the
resonance, I was able to find a steady state for the confined granular system at frequencies
around the first harmonic .
To drive the system at a particular frequency f , I set the motion of the top wall to
follow exactly as
y = y0 +A sin(2pift) (4.7)
where y0 is the equilibrium position of the top wall maintaining a confining pressure,
A is the dynamic amplitude, amplitude of 1 corresponds to a dynamics strain of 5.0e-8.
To ensure that the system was driven at the first resonant peak, I plotted the
horizontal average of the y-component of particle velocities at different instants. Figure
4.10 shows one of such plots. The fit is a sine wave (half period). This demonstrates clearly
that the system is indeed oscillating at its first harmonic.
The response amplitude was measured by fitting the force on the bottom wall with
a single frequency sine curve. It is possible to have higher harmonics to be generated, but
their magnitude was ignored for the current study. Figure 4.11 was a typical least square
curve fit to the force on the bottom wall. The fit was a search for the best amplitude to
93
Figure 4.10: Horizontal average of y-velocity at resonant condition
minimize the sum of the squared distances between the data points and the curve.
I studied three different confining pressures, which were estimated to correspond to
300 kPa, 3.1 MPa and 11.9 MPa respectively in 3-D by assuming the disks to have unit thick-
ness even though they had different diameters. For each confining pressure, five resonance
curves were measured with driving amplitude to be 1, 3, 10, 30 and 100 respectively.
To find the position of resonant peak, I fitted each individual resonance curve with
the equation 4.8, which was the frequency response function of amplitude for a forced
damped oscillator.4
A(f) =
√
A(f0)
1 + [2Q(f − f0)/f0]2 (4.8)
To summarize those resonance curves, I plotted the resonant peaks (normalized to
the resonant frequency of amplitude of 1) as a function of driving amplitude as shown in
figure 4.15. The blue, red and green curves were the normalized5 peak positions as a function
of driving amplitudes for confining pressures of 300 kPa, 3.1 MPa and 11.9 MPa respectively.
To understand the mechanism of resonant frequency shift, I studied the contact force
network of different driving amplitudes for the packing under 300 kPa confining pressure
4The resonance curves deviated slightly from the resonance curves of a classical damped oscillator, but
the fits were good enough for estimating the peak positions.
5Divide the peak positions of higher amplitudes by the peak position of the smallest amplitude A = 1
for each confining pressures.
94
Figure 4.11: Fit the force on the bottom wall
95
Figure 4.12: resonance curves at 300 kPa
96
Figure 4.13: resonance curves at 3.1 MPa
97
Figure 4.14: resonance curve at 11.9 MPa
98
Figure 4.15: Frequency shift as a function of amplitudes
99
Figure 4.16: Relative increase of oscillation intensity of the whole contact force network
driven at 2.96 kHz (resonant frequency for amplitude 1). One can think of the contact force
network as a stable structure of channels and a small displacement of the top wall creates a
burst of perturbations, which act like a fluid flowing through those channels. For resonance
condition, continuous disturbance from the top wall is just like a continuous source of flow.
An interesting discovery is that statistically the wider the channel, the slower the flow. That
is, the stronger the contact, the less oscillation of the force. Figure 4.16 indicates the ratios
of dynamic intensities on every contact forces between amplitude A = 10 and A = 1. For
each contact force, the width of the line is proportional to the magnitude of the force and
the color code indicates the ratio of dynamic intensity (oscillation amplitudes squared and
averaged over more than 20 periods).
First noticeable feature of figure 4.16 is that the relative increase of oscillation in-
tensity is spatially non-homogeneous. Roughly speaking, the middle portion of the packing
increases above average and the portions close to top and bottom walls increase below
average.
100
Figure 4.17: Oscillation intensity comparison of driving amplitude 3, 10, 30, 100 to 1
To study this feature more quantitatively, the ratios of the oscillation intensities
between amplitude A = 3 and A = 1 are plotted as a function of the magnitudes of the
contact forces, which is then repeated for all larger driving amplitudes (A = 10, A = 30
and A = 100). Figure 4.17 shows the distribution of the relative oscillation intensities6 as
driving amplitudes increase.7.
A clear trend of the relative increase of oscillation intensity is that the greater the
driving amplitude the wider the scattering of the relative increase. This finding can provide
an explanation for the resonant frequency shift if one can somehow relate the magnitude of
the contact forces to the resonant mode of the system.
6Divided by the dynamic intensity on that contact force when amplitude A = 1.
7Each individual plot is given in appendix 6.3.
101
4.4 summary and brief report of ongoing studies
More data and analysis are needed to check if there is any change of waveform in order to
determine which model provides a better energy dissipation mechanism for LAMMPS.
More analysis and theoretical work are needed to provide complete explanation for
resonant frequency shift.
102
Chapter 5
Conclusions
In this thesis I studied three separate but closely related subjects.
In the low temperature fracture experiment, tremendous effort was spent to search
for the velocity gap, a range of speed that no steady crack propagation is possible when the
thermal energy was low enough (< 100 Kelvin). This is a feature solely due to the discrete
nature of the fracture process. The data I gathered provided a good indication that the
velocity gap does exist.
In the friction experiment, the data clearly overturned some commonly held under-
standing of static friction, which eventually led us to propose a modified rate- and state-
equation and unified the mathematical description of friction as well as the physical under-
standing to it.
In the granular simulation, two different modes (pulse mode and resonance mode) of
probing a confined granular system were conducted with a parametric study of the influence
of the confining pressure as well as the driving amplitude. The pulse mode provided infor-
mation on attenuation, which was used to distinguish between two widely accepted energy
dissipation models. Due to the lack of change in the waveform of a single propagating sine
wave, LAMMPS clearly favor the array of energy sinks model over the hysterestic model
as the major energy dissipation mechanism in a confined granular medium. And, the reso-
nance mode provided an intuitive explanation for the resonant frequency shift in terms of
how dynamic loads were re-distributed through the contact force network.
103
Chapter 6
Appendix
6.1 Wafer preparation for the fracture experiment
Here is the procedure to prepare a wafer before it can be used in the fracture experiment.
A set of raw wafer (usually 5 pieces) is firstly baked in an oven under 1000 degree C
for 16~18 hours. This is to grow an insulating layer of silicon dioxide. The purpose of that
is to electrically isolate the wafer from its environment when it is held to the metal frame
lately.
The next step is to put electrical leads on the wafer, which is bridged by a layer of
aluminum. In order to properly position and shape the leads, a set of aluminum masks is
designed. The baked wafers are approximately align with the mark on the mask and fixed
to position by a small piece of adhesive tape. Then, the wafers and masks are carefully
stacked together and bundled by two paper clippers. The stack should be so arranged that
all the parts of wafers that need to be deposited by aluminum are facing the same side and
as close to each other as possible. Then, the stack together with a vaporization bowl and
a piece of crystal are delivered to the MBE (Molecular Beam Epitaxy) chamber on the 3rd
floor to get a layer of 4000 A aluminum deposition.
After the aluminum deposition, a second bake under 550 degree C for 3 hours is
required to drive the aluminum atoms through the silicon dioxide layer and make a good
ohmic contact to the single crystal silicon underneath it.
104
Figure 6.1: Two identical balls collide
Afterward, a set of wires and tabs (small pieces of copper sheet 2~3 mm X 2~3
mm) are prepared. The tabs are made from a piece of thin copper sheet, cut into small
rectangular with desired dimensions and flatted by milling them between flat surfaces of
two big pieces of bulk iron (one at a time to ensure the flatness of the tab). Then, wires are
soldered to the tabs.
The last step is to wire the wafer, that is to connect the tabs to the leads. The
agent that stick the tabs to the aluminum leads is some special purpose silver glue. The
silver glue should be kept refrigerated till the moment it is used. A tiny amount of silver
glue should be applied to the flat side of the copper tab and then small wire campers are
used to hold the tabs firmly on the wafer. After all four tabs are carefully glued onto the
leads of the wafer and firmly held by the campers, the wafer is briefly bake again under 175
degree C for about 4~5 minutes. The baking quickly solidified the silver glue, which makes
the connection permanent (the campers can then be detached).
6.2 Bug report to LAMMPS admins
Dear LAMMPS admins,
I'm trying to illustrate here a bug within the granular package of LAMMPS.
I performed collisions between two identical balls, as indicated in figure 6.1.
In all the simulation, there is no gravity. Equal but opposite velocities are given to
the balls as initial condition.
105
Figure 6.2: y-force on ball 2 with radius 2
The normal force on ball 2 is calculated in LAMMPS and compared to independent
calculations by Matlab. For the first two simulations, γn is set to zero, that is, visco-elastic
damping is turned off for the purpose of pinning down the problem.
The force calculation by Matlab is based on equation 6.1, which can be found on
page 426, LAMMPS user manual:
F =
√
d− r
d
(Kn(d− r)− γnmeffVn) (6.1)
d = R1 +R2 is the contact distance between two balls of radius R1and R2
r is the distance between the centers of contacting balls
Kn is the elastic constant for normal contact
106
Figure 6.3: y-force on ball 2 with radius 3
107
γn is the viscoelastic constants for normal contact
meff = mimj/(mi +mj) is the effective mass of two balls of mass mi and mj
Vn is the normal component of the relative velocity of the two balls
To my surprise, the force calculated this way (the red curve) does NOT match the
force calculated by LAMMPS (the blue circles). However, the forces calculated by LAMMPS
can be matched by calculation based on equation 6.2(the green curve).
F =
√
d− r(Kn(d− r)− γnmeffVn) (6.2)
It is believed that equation 6.1 is correct. It is listed in LAMMPS user manual, and
it is also seen in literature, for example, the experimental paper by Laurent Labous[54], etc.
This finding is further confirmed by looking into the source code (pair_gran_hertzian.cpp,
line 145-154):
// normal damping term
// this definition of DAMP includes the extra 1/r term
xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) xmeff = rmass[j];
if (mask[j] & freeze_group_bit) xmeff = rmass[i];
damp = xmeff*gamman_dl*vnnr/rsq;
ccel = xkk*(radsum-r)/r - damp;
rhertz = sqrt(radsum - r);
ccel = rhertz * ccel;
To me, it is a bit surprised to find out that LAMMPS implemented normal contact force (for
granular particles) by a single parameter Kn instead of Young's modulous E and Possion's
ratio ν. As this will impose an unnecessary restriction of mono-dispersity to the systems
that LAMMPS can simulate, because Kn is actually a function of d, (Kn = Ed). Even for
mono-dispersity, the forces so calculated are off by a factor of
√
d. So, please look into this
matter and fix the bug.
108
Along this line of investigation, I also find the following paragraph misleading (page
198, LAMMPS user manual):
IMPORTANTNOTE: Some models in LAMMPS treat particles as extended spheres
or ellipsoids, as opposed to point particles. In 2d, the particles will still be spheres or ellip-
soids, not circular disks or ellipses, meaning their moment of inertia will be the same as in
3d.
What I found is that the mass of the balls is calculated differently for 2d and 3d
simulations.
M = pir2ρ for 2d
and
M = 43pir
3ρfor 3d.
That is, for 2d simulations, mass of the balls IS calculated as circular disks with
unit thickness.
To illustrate that, I did another 2 simulations with γn equal to 25, such that different
ways of mass calculation could affect the force curve.
6.2.1 The input scripts are as following:
The collision.src file (used as ./lmp_serial < collision.src)
# collision test
dimension 3
atom_style granular
boundary p f p
newton off
read_data collision.in
# gamma_n set to zero, no visco-elastic damping
pair_style gran/hertzian 200000.0 0.0 0.5 0
timestep 0.000001
fix 1 all nve/sphere
109
Figure 6.4: collision simulated in LAMMPS (2-d and 3-d)
110
fix ywalls all wall/gran yplane 0 100 50 0
dump 1 all custom 100 collision_gamma.out tag x y vx vy radius fx fy
run 1000000
The collision.in file
# first line will always be ignored.
# collision test of granular package of LAMMPs
# header section
2 atoms
1 atom types
0 40 xlo xhi
0 40 ylo yhi
-10 10 zlo zhi
# body section
Atoms
1 1 6.0 1.0 14.0 9.0 0.0
2 1 6.0 1.0 14.0 21.0 0.0
Velocities
1 0.0 5.0 0.0 0.0 0.0 0.0
2 0.0 -5.0 0.0 0.0 0.0 0.0
Matlab script
%collision between two identical balls
%parameters as in LAMMPS
kn = 200000.0;
dens=1.0;
gamma_n=0.0;
111
d=6.0; % twice the radius, R1+R2
r=ball2_y-ball1_y; % balls are lined up in x-direction
vn=ball2_vy-ball1_vy; % relative velocity
meff=0.5*(4.0/3.0*pi*(d/2.0)^3)*dens; %effective mass
fy=sqrt((d-r)/d).*(kn*(d-r)-meff*gamma_n*vn); % equation for calculate nor-
mal force as in page 426, LAMMPS manual
fy_m = sqrt(d)*fy;
plot(ball2_fy(5600:7200),'b')
hold on
plot(real(fy(5600:7200)),'r')
plot(real(fy_m(5600:7200)),'g')
xlabel('Time (1.0e-4s)')
ylabel('y-force on ball 2')
6.3 Graphs of Relative Increase of Oscillation Intensity
112
Figure 6.5: A3/A1
113
Figure 6.6: A10/A1
114
Figure 6.7: A30/A1
115
Figure 6.8: A100/A1
116
Bibliography
[1] T. Baumberger A. Cochard, L. Bureau. Stabilization of frictional sliding by normal
load modulation. Trans. ASME, 70:220226, 2003.
[2] F. F. Abraham. Some new Directions in science on computers. World Scientific, 1997.
[3] T. L. Anderson. Fracture mechanics: Fundamentals and applications. CRC Press, 1995.
[4] M. I. Baskes. Application of the embedded-atom method to covalent materials: A
semiempirical potential for silicon. PHYSICAL REVIEW LETTERS, 59:26662669,
1987.
[5] Martin Z. Bazant, Efthimios Kaxiras, and J. F. Justo. Environment-dependent inter-
atomic potential for bulk silicon. PHYSICAL REVIEW B, 56:85428552, 1997.
[6] T. S. Majmudar & R. P. Behringer. Contact force measurements and stress-induced
anisotropy in granular materials. Nature, 435:10791082, 2005.
[7] N. Bernstein and D. Hess. Multiscale simulations of brittle fracture and the quantum-
mechanical nature of bonding in silicon. Mat. Res. Soc. Symp. Proc., 653:Z2.7.1Z2.7.6,
2001.
[8] B. Bhushan, editor. Modern Tribology Handbook. CRC Press, 2001.
[9] Michael L. Blanpied Brian D. Kilgore and James H. Dieterich. Velocity dependent
friction of granite over a wide range of conditions. GEOPHYSICAL RESEARCH LET-
TERS, 20:903906, 1993.
117
[10] Michael J. Buckingham. Theory of acoustic attenuation, dispersion, and pulse propa-
gation in unconsolidated granular material including marine sediments. J. Acoust. Soc.
Am., 102:25792596, 1997.
[11] Andreas A. Polycarpou Chul-Hee Lee. Static friction experiments and verification of an
improved elastic-plastic model including roughness effects. Trans. ASME, 129:754760,
2007.
[12] Murray S. Daw and M. I. Baskes. Semiempirical, quantum mechanical calculation of
hydrogen embrittlement in metals. PHYSICAL REVIEW LETTERS, 50:12851288,
1983.
[13] S. M. Day. Efficient simulation of constant q using coarse-grained memory variables.
Bulletin of the Seismological Society of America, 88:10511062, 1998.
[14] James H. Dieterich. Time-dependent friction in rocks. JOURNAL OF GEOPHYSICAL
RESEARCH, 77:36903697, 1972.
[15] James H. Dieterich. Modeling of rock friction 1 experimental results and constitutive
equations. JOURNAL OF GEOPHYSICAL RESEARCH, 85:21612175, 1979.
[16] JAMES H. DIETERICH and BRIAN D. KILGORE. Direct observation of frictional
contacts: New insights for state-dependent properties. PAGEOPH, 143:283302, 1994.
[17] Duncan Dowson. History of Tribology, 2nd Edition. Wiley, 1998.
[18] Jacco H. Snoeijer Martin van Hecke Ellak Somfai, Jean-Noel Roux and Wim van Saar-
loos. Elastic wave propagation in confined granular systems. PHYSICAL REVIEW E,
72:021301, 2005.
[19] R A Rafey Farid F Abraham, D Brodbeck and W E Rudge. Instability dynamics
of fracture: A computer simulation investigation. PHYSICAL REVIEW LETTERS,
73:272275, 1994.
[20] Richard P. Feynman. The Feynman Lectures on Physics. Addison Wesley Longman,
1970.
118
[21] L. H. Seeley K. H. Kim E. A. Behne S. Zaranek B. D. Chapman S. M. Heald G. T. Sei-
dler, G. Martinez and D. L. Brewe. Granule-by-granule reconstruction of a sandpile
from x-ray microtomography data. PHYSICAL REVIEW E, 62:81758181, 2000.
[22] E. E. Gdoutos. Fracture Mechanics, an introduction. Springer, 2005.
[23] Eri Stendahl Gerde. Fracture and Friction. PhD thesis, The University of Texas at
Austin, 2001.
[24] Eric Gerde and M. Marder. Friction and fracture. Nature, 413:285288, 2001.
[25] Bruno Gilles and Christophe Coste. Low-frequency behavior of beads constrained on a
lattice. PHYSICAL REVIEW LETTERS, 90:174302, 2003.
[26] David L. Goldsby Terry E. Tullis Giulio Di Toro. Friction falls towards zero in quartz
rock as slip velocity approaches seismic rates. Nature, 427:436439, 2004.
[27] A. A. Griffith. The phenomena of rupture and flow in solids. Philosophical Transactions
of the Royal Society of London A, 221:163198, 1921.
[28] A. A. Griffith. The theory of rupture. Proceedings of First International Congress of
Applied Mechanics, UNKNOWN:5563, 1924.
[29] S. P. Gross. Dynamics of Fast Fracture. PhD thesis, University of Texas at Austin,
1995.
[30] R. A. Guyer, K. R. McCall, and G. N. Boitnott. Hysteresis, discrete memory, and
nonlinear wave propagation in rock a new paradigm. PHYSICAL REVIEW LETTERS,
74:34913494, 1995.
[31] Robert A. Guyer and Paul A. Johnson. Nonlinear mesoscopic elasticity: Evidence for
a new class of materials. Physics Today, April:3036, 1999.
[32] Jen Andreas Hauch. Dynamic Fracture in Brittle Materials. PhD thesis, University of
Texas at Austin, 1998.
[33] H. Hertz. unknown. J. reine angew. Math., 92:156171, 1881.
119
[34] Dominic Holland and M. Marder. Ideal brittle fracture of silicon studied with molecular
dynamics. PHYSICAL REVIEW LETTERS, 80:746749, 1998.
[35] Dominic Holland and Michael Marder. Cracks and atoms. ADVANCED MATERIALS,
11:793806, 1999.
[36] Uri Lev Ilan Beery and Dov Sherman. On the lower limiting velocity of a dynamic
crack in brittle solids. JOURNAL OF APPLIED PHYSICS, 93:24292434, 2003.
[37] Inglis. Stresses in a plate due to the presence of cracks and sharp corners. Transactions
of the Institution of Naval Architects, 55:219230, 1913.
[38] G. R. Irwin. Analysis of stresses and strains near the end of a crack traversing a plate.
Journal of Applied Mechanics, 24:361364, 1957.
[39] M. I. Baskes J. G. Swadener and M. Nastasi. Molecular dynamics simulation of brittle
fracture in silicon. PHYSICAL REVIEW LETTERS, 89:085503, 2002.
[40] H Peerhossaini D. Bideau J. Lemaitre, A. Gervois and J. P. Troadec. An air table
designed to study two-dimensional disc packings: preliminary tests and first results. J.
Phys. D: Appl. Phys., 23:13961404, 1990.
[41] Francisco Guinea John Ferrante James H. Rose, John R. Smith. Universal features of
the equation of state of metals. PHYSICAL REVIEW B, 29:29632969, 1984.
[42] M. P. Marder Jens A. Hauch, Dominic Holland and Harry L. Swinney. Dynamic fracture
in single crystal silicon. PHYSICAL REVIEW LETTERS, 82:38233826, 1999.
[43] D. Gourdon M. Ruths J. N. Israelachvili Jianping Gao, W. D. Luedtke and Uzi Land-
man. Frictional forces and amontons law: From the molecular to the macroscopic scale.
J. Phys. Chem. B, UNKNOWN:15, 2003.
[44] K. L. Johnson. Contact Mechanics. Cambridge University Press, 1987.
[45] Paul A. Johnson and Xiaoping Jia. Nonlinear dynamics, granular media and dynamic
earthquake triggering. Nature, 437:871874, 2005.
[46] T. Johnson. Time dependent friction of granite: Implications for pre-cursory slip on
faults. JOURNAL OF GEOPHYSICAL RESEARCH, 86:60176028, 1981.
120
[47] David A. Kessler. A new crack at friction. Nature, 413:260261, 20 SEPTEMBER
2001.
[48] A. M. Khilla. Optimum continuous microstrip tapers are amenable to computer-aided
design. Microwave Journal, unknown:221224, 1983.
[49] Alick C. Kibblewhite. Attenuation of sound in marine sediments: A review with em-
phasis on new low-frequency data. J. Acoust. Soc. Am., 86:716738, 1989.
[50] Lior Kogut and Izhak Etsion. A static friction model for elastic-plastic contacting rough
surfaces. Journal of Tribology, 126:3440, 2004.
[51] J. Kohanoff. Electronic Structure Calculations for Solids and Molecules: Theory and
Computational Methods. Cambridge University Press, 2006.
[52] J. Krim. Friction at macroscopic and microscopic length scales. Am. J. Phys., 70:890
897, 2002.
[53] G. MacDonald L. Knopoff. Attenuation of small amplitude stress waves in solids.
Reviews of Modern Physics, 30:11781192, 1958.
[54] Anthony D. Rosato Laurent Labous and Rajesh N. Dave. Measurements of colli-
sional properties of spheres using high-speed video analysis. PHYSICAL REVIEW
E, 56:57175725, 1997.
[55] B. Lawn. Fracture in Brittle Solids. Cambridge University Press, 1993.
[56] A. E. H. Love. Treatise on the Mathematical Theory of Elasticity. Cambridge University
Press, 1927.
[57] M. Urbakh M. G. Rozman and J. Klafter. Stick-slip motion and force fluctuations in a
driven two-wave potential. PHYSICAL REVIEW LETTERS, 77:683686, 1996.
[58] D. H. Johnston M. N. Toksoz and A. Timur. Attenuation of seismic waves in dry and
saturated rocks: I. laboratory measurements. GEOPHYSICS, 44:681690, 1979.
[59] Valentina Agostini Macro Scalerandi and Paul A. Johnson Pier Paolo Delsanto, Koen
Van Den Abeele. Local interaction simulation approach to modelling nonclassical non-
linear elastic behavior in solids. J. Acoust. Soc. Am., 113:30493059, 2003.
121
[60] M. Marder and Steve Gross. Origin of crack tip instabilities. J. Mech. Phys. Solids,
43:148, 1995.
[61] M. Marder and Xiangming Liu. Instability in lattice fracture. PHYSICAL REVIEW
LETTERS, 71:24172420, 1993.
[62] Michael Marder. Energies of a kinked crack line.
[63] MICHAEL MARDER. Moleculardynamics ofcracks. COMPUTING IN SCIENCE &
ENGINEERING, SEPTEMBER/OCTOBER:4855, 1999.
[64] Michael P. Marder. Condensed Matter Physics. Wiley-interscience, 2000.
[65] Chris Marone. Laboratory-derived friction laws and their application to seismic fault-
ing. Annu. Rev. Earth Planet. Sci., 26:643696, 1998.
[66] Philip T. Metzger. Granular contact force density of states and entropy in a modified
edwards ensemble. PHYSICAL REVIEW E, 70:051303, 2004.
[67] Delphine Gourdon Jacob Israelachvili Michael Urbakh, Joseph Klafter. The nonlinear
nature of friction. Nature, 430:29, 2004.
[68] Mason W. P. Osmer T. F. Mindlin, R. D. and H. Dereciewicz. Effect of an oscillatory
tangential force on contact surfaces of elastic spheres. Proceedings of the First National
Congress of Applied Mechanics, unknown:203228, 1952.
[69] R. D. Mindlin. Compliance of elastic bodies in contact. Journal of Applied Mechanics,
pages 259268, 1949.
[70] JULIA K. MORGAN. Particle dynamics simulations of rate- and state-dependent
frictional sliding of granular fault gouge. Pure appl. geophys., 161:18771891, 2004.
[71] L. A. Ostrovsky and P. A. Johnson. Dynamic nonlinear elasticity in geomaterials.
RIVISTA DEL NUOVO CIMENTO, 24:46, 2001.
[72] Fabrice Barbe Stephane Bourles Xavier Thibault Patrick Richard, Pierre Philippe and
Daniel Bideau. Analysis by x-ray microtomography of a granular packing undergoing
compaction. PHYSICAL REVIEW E, 68:020301, 2003.
122
[73] B. N. J. Persson. Sliding Friction, Physical Properties and Applications. Springer, 2000.
[74] Carly M. Donahue Philip T. Metzger. Elegance of disordered granular packings: A
validation of edward's hypothesis. PHYSICAL REVIEW LETTERS, 94:148001, 2005.
[75] S. J. Plimpton. Fast parallel algorithms for short-range molecular dynamics. Journal
of Computational Physics, 117:119, 1995.
[76] Lisa Patel M. Marder Robert D. Deegan, Shilpa Chheda and Harry L. Swinney. Wavy
and rough cracks in silicon. PHYSICAL REVIEW E, 67:066209, 2003.
[77] Andy Ruina. Slip instability and state variable friction laws. JOURNAL OF GEO-
PHYSICAL RESEARCH, 88:1035910370, 1983.
[78] R. J. Sanford. Principles of Fracture Mechanics. Pearson Education, Inc., 2003.
[79] Gil Cohen Jay Fineberg Shmuel M. Rubinstein. Detachment fronts and the onset of
dynamic friction. Nature, 430:10051009, 2004.
[80] Leonid I. Slepyan. Dynamics of a crack in a lattice. Sov. Phys. Dokl., 26:538540, 1981.
[81] Leonid I. Slepyan. Models and Phenomena in Fracture Mechanics. Springer, 2002.
[82] Frank H. Stillinger and Thomas A. Weber. Computer simulation of local order in
condensed phases of silicon. PHYSICAL REVIEW B, 31:52625271, 1985.
[83] Robert D. Stoll and Robert E. Houtz. Attenuation measurement with sonobuoys. J.
Acoust. Soc. Am., 73:163172, 1983.
[84] M. Saadatfar T. Aste and T. J. Senden. Geometrical structure of disordered sphere
packings. PHYSICAL REVIEW E, 71:061302, 2005.
[85] H.J. Herrmann T. PÄoschel. A simple geometrical model for solid friction. Physica A,
198:441448, 1993.
[86] F. P. Bowden; D. Tabor. The area of contact between stationary and between moving
surfaces. Proceedings of the Royal Society of London. Series A, Mathematical and
Physical Sciences, 169:391413, 1939.
123
[87] Peter A. Thompson and Mark O. Robbins. Origin of stick-slip motion in boundary
lubrication. Science, 250:792794, 1990.
[88] Christiane Caroli Tristan Baumberger and Olivier Ronsin. Self-healing slip pulses along
a gel/glass interface. PHYSICAL REVIEW LETTERS, 88:075509, 2002.
[89] Martin van Hecke. A tale of tails. Nature, 435:10411042, 2005.
[90] L. R. Walker and N. Wax. Non-uniform transmission lines and reflection coefficients.
Journal of Applied Physics, 17:10431045, 1946.
[91] Martin H. Müser Ludgar Wenning and Mark O. Robbins. Simple microscopic theory
of amontons laws for static friction. PHYSICAL REVIEW LETTERS, 86:12951298,
2001.
[92] B. Velicky X. Jia, C. Caroli. Ultrasound propagation in externally stressed granular
media. PHYSICAL REVIEW LETTERS, 82:18631866, 1999.
[93] Hepeng Zhang Zhiping Yang and Michael P. Marder. Dynamics of static friction
between steel and silicon. Proceedings of National Academy of Sciences (PNAS),
105:1326413268, 2008.
124
Vita
Zhiping Yang was born on October 26th, 1979 in Shanghai, China. He is the first son of
Jinming Yang and Weiqiu Wang. He received his high school diploma from SongJiang No.
2 Senior High School in 1998. He spent his freshman year in Fudan University. He went to
City University of Hongkong with full scholarship support from Jockey Club of Hongkong
and acquired his Bachelors of Science in Applied Physics with first class honor in 2002.
From fall 2002, he began graduate school at the University of Texas at Austin pursuing a
Ph. D. in Physics. On June 8th of 2002, he engaged with Hongwei Shi before he left China.
He married her on August 1st, 2004. Steven, his first son was born on May 26th, 2006.
Permanent Address: No. 1086, Jinwei Group 2, Jinshanwei Zheng, Jinshan Qu,
Shanghai, China, 201512
This dissertation was typeset with LATEX2ε
1 by the author.
1LATEX2ε is an extension of LATEX. LATEX is a collection of macros for TEX. TEX is a trademark of the
American Mathematical Society. The macros used in formatting this dissertation were written by Dinesh
Das, Department of Computer Sciences, The University of Texas at Austin, and extended by Bert Kay,
James A. Bednar, and Ayman El-Khashab.
125