Impacts of Kinematic Links with a Granular Material
Except where reference is made to the work of others, the work described in this
dissertation is my own or was done in collaboration with my advisory committee. This
dissertation does not include proprietary or classifled information.
Seunghun Lee
Certiflcate of Approval:
George Flowers
Professor
Mechanical Engineering
Dan B. Marghitu, Chair
Professor
Mechanical Engineering
Nels Madsen
Professor
Mechanical Engineering
John Y. Hung
Professor
Electrical and Computer Engineering
George Flowers
Dean
Graduate School
Impacts of Kinematic Links with a Granular Material
Seunghun Lee
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulflllment of the
Requirements for the
Degree of
Doctor of Philosophy
Auburn, Alabama
August 10, 2009
Impacts of Kinematic Links with a Granular Material
Seunghun Lee
Permission is granted to Auburn University to make copies of this dissertation at its
discretion, upon the request of individuals or institutions and at
their expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Seunghun Lee, son of Kyuwan Lee and Junghye Ahn was born on March 25, 1971,
in Seoul, South Korea. He had completed the degree of B.S. in 1994 from Korea Military
Academy and the degree of M.S. in 2000 from University of Tsukuba, Japan. He joined the
Mechanical Engineering Department at Auburn University in September 2005.
iv
Dissertation Abstract
Impacts of Kinematic Links with a Granular Material
Seunghun Lee
Doctor of Philosophy, August 10, 2009
(M.S., University of Tsukuba, 2000)
186 Typed Pages
Directed by Dan B. Marghitu
The impact of kinematic links into a granular material is studied. In the impact
process through the granular material deflned as a conglomeration of discrete solids, the
most important force governing the motion of a body is the resistance force of the granular
matter. The resistance force acting on the body is studied as a linear superposition of
a static (depth-dependent) resistance force and a dynamic (velocity-dependent) frictional
force by considering the characteristics of the granular material acting like solid and liquid.
The impact models considered are spheres, mathematical pendulum, compound pendulum,
open kinematic chains, and elastic links. The impact is observed from the contact moment
until the rest using various initial conditions. In this study, the impact angle and the initial
velocity are the dominant conditions deciding the whole impact process. The results of the
simulations and the experiments are analyzed based on the initial impact conditions. For
most of the analyzed impact systems including elastic objects, the system under high force
impact (higher initial velocity) comes to rest faster in the granular matter than the same
system under low-force impacts (lower initial velocity).
v
Acknowledgments
I would like to thank Professor Marghitu for providing me the opportunity to conduct
research in this exciting fleld and guidance during the course of this work. I would also
like to thank the government of South Korea for ofiering me the chance to study abroad
at government expense. I would also like to thank all my advisory committee members,
Professors George Flowers, Nels Madsen, and John Y. Hung for their suggestions. Finally,
I wish thank all the professors, stafi, and students I have interacted with throughout this
study and Auburn University for a peaceful study environment.
I wish to dedicate this work to my wife, Hye Yon Cho, whose endless love and support
in my life has enabled what I have accomplished including this work and has always been
the origin of my desire to move forward.
vi
Style manual or journal used Journal of Approximation Theory (together with the style
known as \auphd"). Bibliograpy follows van Leunen?s A Handbook for Scholars.
Computer software used The document preparation package TEX (speciflcally
LATEX) together with the departmental style-flle auphd.sty, Adobe Illustrator CS, and
Mathematica 6.0.
vii
Table of Contents
List of Figures x
1 Introduction 1
2 Description of the problem 5
2.1 Coordinate system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Description of the models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.1 Impact of a rigid sphere . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.2 Single impact of a planar kinematic chain . . . . . . . . . . . . . . . 9
2.2.3 Multiple impacts of a planar kinematic chain with n links . . . . . . 12
2.3 Generalized coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3 Force analysis 16
3.1 Dynamic frictional force Fd . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2 Static resistance force Fs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.1 Horizontal static resistance force Fsh . . . . . . . . . . . . . . . . . . 20
3.2.2 Vertical static resistance force Fsv . . . . . . . . . . . . . . . . . . . 22
3.3 Air drag force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.4 Application to the impact models . . . . . . . . . . . . . . . . . . . . . . . . 24
4 Impact of rigid objects 28
4.1 Impact of a rigid sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.1.1 Impact of a mathematical pendulum . . . . . . . . . . . . . . . . . . 29
4.1.2 Oblique impact of a rigid sphere . . . . . . . . . . . . . . . . . . . . 36
4.2 Single impact with a granular medium . . . . . . . . . . . . . . . . . . . . . 42
4.2.1 Impact of a free rigid link . . . . . . . . . . . . . . . . . . . . . . . . 42
4.2.2 Impact of a planar kinematic chain . . . . . . . . . . . . . . . . . . . 52
4.3 Multiple impacts of a planar kinematic chain . . . . . . . . . . . . . . . . . 70
5 Impact of flexible links 83
5.1 Free elastic link . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.1.1 Transverse vibration of a free link . . . . . . . . . . . . . . . . . . . 85
5.1.2 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.1.3 Equations of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.1.4 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.2 Elastic compound pendulum . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
5.2.1 Transverse vibration of an elastic compound pendulum . . . . . . . . 99
5.2.2 Modeling of an articulated elastic compound pendulum . . . . . . . 102
viii
5.2.3 Simulation results of the articulated elastic compound pendulum . . 106
5.2.4 Modeling of a cantilevered elastic compound pendulum . . . . . . . 114
5.2.5 Simulation results of the cantilevered elastic compound pendulum . 118
6 Experiments 129
6.1 Experimental system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
6.2 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
6.2.1 Impact of a free link . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
6.2.2 Impact of a rigid compound pendulum . . . . . . . . . . . . . . . . . 137
7 Conclusion 144
Bibliography 147
Appendices: Mathematica Programs 152
A Mathematical pendulum 153
B Oblique impact of a sphere 155
C Rigid free link 157
D Rigid compound pendulum 159
E Rigid double pendulum 162
F Two link kinematic chain with two impact points 165
G Elastic free link 168
H Elastic compound pendulum 171
ix
List of Figures
2.1 Coordinate system for rigid impact bodies . . . . . . . . . . . . . . . . . . 8
2.2 Coordinate system for exible impact bodies . . . . . . . . . . . . . . . . . 8
2.3 Oblique impact of a sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Impact of a mathematical pendulum . . . . . . . . . . . . . . . . . . . . . . 10
2.5 Impact of a free link . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.6 Single impact of a kinematic chain . . . . . . . . . . . . . . . . . . . . . . . 12
2.7 Multiple impacts of a kinematic chain . . . . . . . . . . . . . . . . . . . . . 13
3.1 Air drag force impeding the motion of the compound pendulum . . . . . . 24
3.2 Moving angle, qm, of a link in planar motion . . . . . . . . . . . . . . . . . 26
4.1 Free body diagram of a mathematical pendulum in impact . . . . . . . . . 29
4.2 Impact results of the mathematical pendulum into a granular medium for
q(0) = 15o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.3 Impact results of the mathematical pendulum into a granular medium for
q(0) = 45o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.4 Impact results of the mathematical pendulum into a granular medium for
q(0) = 75o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.5 Free body diagram of a rigid sphere in oblique impact . . . . . . . . . . . . 36
4.6 Impact results of the rigid sphere into a granular medium for q(0) = 15o . . 39
4.7 Impact results of the rigid sphere into a granular medium for q(0) = 45o . . 40
4.8 Impact results of the rigid sphere into a granular medium for q(0) = 75o . . 41
4.9 Free body diagram of a free rigid link in impact . . . . . . . . . . . . . . . 43
x
4.10 Impact results of the rigid cylinder type link for q(0) = 0o . . . . . . . . . . 47
4.11 Resistance force Fd and Fs for q(0) = 0o . . . . . . . . . . . . . . . . . . . . 48
4.12 Total resistance force FR for q(0) = 0o . . . . . . . . . . . . . . . . . . . . . 49
4.13 Impact results of the rigid cylinder type link for q(0) = 32o . . . . . . . . . 50
4.14 Impact results of the rigid cylinder type link for q(0) = 55o . . . . . . . . . 51
4.15 Free body diagram of a single impact of a planar kinematic chain . . . . . 52
4.16 Compound pendulum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.17 Impact results of the rigid compound pendulum for q(0) = 22o . . . . . . . 59
4.18 Impact results of the rigid compound pendulum for q(0) = 31o . . . . . . . 60
4.19 Impact results of the rigid compound pendulum for q(0) = 45o . . . . . . . 61
4.20 Impact results of the rigid compound pendulum for q(0) = 61:5o . . . . . . 62
4.21 Double pendulum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.22 Angularvelocities _q1 and _q2 resultsofthedoublependulumimpactforq1(0) =
30o, q2(0) = 75o, and _q1(0) = _q2(0) = -1, -3, -5 rad/s . . . . . . . . . . . . . 65
4.23 Velocities vTx and vTz results of the double pendulum impact for q1(0) = 30o,
q2(0) = 75o, and _q1(0) = _q2(0) = -1, -3, -5 rad/s . . . . . . . . . . . . . . . 66
4.24 Angularvelocities _q1 and _q2 resultsofthedoublependulumimpactforq1(0) =
75o, q2(0) = 30o, _q1(0) = -5 rad/s, and _q2(0) = -1, -3, -5 rad/s . . . . . . . 67
4.25 Velocities vTx and vTz results of the double pendulum impact for q1(0) = 75o,
q2(0) = 30o, _q1(0) = -5 rad/s, and _q2(0) = -1, -3, -5 rad/s . . . . . . . . . . 68
4.26 Free body diagram of multiple impacts of a planar kinematic chain . . . . 71
4.27 Two link kinematic chain with two impact points model . . . . . . . . . . . 71
4.28 Velocities vT1z and vT2z results of the two link kinematic chain for q1(0) =
330o, q2(0) = 45o, and vAz(0) = 1, 3, 5 m/s . . . . . . . . . . . . . . . . . . 77
4.29 Velocities vT1z and vT2z results of the two link kinematic chain for q1(0) =
330o, q2(0) = 45o, and _q1(0) = ?_q2(0) = 1, 3, 5 rad/s . . . . . . . . . . . . 78
4.30 Velocities vT1z and vT2z results of the two link kinematic chain for q1(0) =
330o, q2(0) = 45o, and _q1(0) = 1, 3, 5 rad/s . . . . . . . . . . . . . . . . . . 79
xi
4.31 Velocities vT1z and vT2z results of the two link kinematic chain for q1(0) =
330o, q2(0) = 45o, and _q2(0) = -1, -3, -5 rad/s . . . . . . . . . . . . . . . . 80
5.1 Transverse vibrations of a exible link . . . . . . . . . . . . . . . . . . . . . 83
5.2 Free elastic link . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.3 Displacement zT and velocity vTz results of the exible link for q1(0) = 32o
and _qz(0) = 1.26, 1.87, 2.33 m/s . . . . . . . . . . . . . . . . . . . . . . . . 95
5.4 Displacement q2 and velocity _q2 results of the exible link for q1(0) = 32o
and _qz(0) = 1.26, 1.87, 2.33 m/s . . . . . . . . . . . . . . . . . . . . . . . . 96
5.5 Displacement zT and velocity vTz results of the exible link for q1(0) = 55o
and _qz(0) = 1.45, 1.98, 2.43 m/s . . . . . . . . . . . . . . . . . . . . . . . . 97
5.6 Displacement q2 and velocity _q2 results of the exible link for q1(0) = 55o
and _qz(0) = 1.45, 1.98, 2.43 m/s . . . . . . . . . . . . . . . . . . . . . . . . 98
5.7 Articulated elastic compound pendulum . . . . . . . . . . . . . . . . . . . . 100
5.8 Cantilevered elastic compound pendulum . . . . . . . . . . . . . . . . . . . 100
5.9 Angle q1 and angular velocity _q1 results of the articulated elastic compound
pendulum for q1(0) = 30o and _q1(0) = -1, -3, -5 rad/s . . . . . . . . . . . . 107
5.10 Displacement q2 and velocity _q2 results of the articulated elastic compound
pendulum for q1(0) = 30o and _q1(0) = -1, -3, -5 rad/s . . . . . . . . . . . . 108
5.11 Angle q1 and angular velocity _q1 results of the articulated elastic compound
pendulum for q1(0) = 60o and _q1(0) = -1, -3, -5 rad/s . . . . . . . . . . . . 109
5.12 Displacement q2 and velocity _q2 results of the articulated elastic compound
pendulum for q1(0) = 60o and _q1(0) = -1, -3, -5 rad/s . . . . . . . . . . . . 110
5.13 Angle q1 and angular velocity _q1 results of the articulated elastic compound
pendulum for q1(0) = 75o and _q1(0) = -1, -3, -5 rad/s . . . . . . . . . . . . 111
5.14 Displacement q2 and velocity _q2 results of the articulated elastic compound
pendulum for q1(0) = 75o and _q1(0) = -1, -3, -5 rad/s . . . . . . . . . . . . 112
5.15 Angle q1 and angular velocity _q1 results of the cantilevered elastic compound
pendulum for q1(0) = 22o and _q1(0) = -1.75, -3.38, -4.66 rad/s . . . . . . . 119
5.16 Displacement q2 and velocity _q2 results of the cantilevered elastic compound
pendulum for q1(0) = 22o and _q1(0) = -1.75, -3.38, -4.66 rad/s . . . . . . . 120
xii
5.17 Angle q1 and angular velocity _q1 results of the cantilevered elastic compound
pendulum for q1(0) = 31o and _q1(0) = -3.31, -6.24, -8.41 rad/s . . . . . . . 121
5.18 Displacement q2 and velocity _q2 results of the cantilevered elastic compound
pendulum for q1(0) = 31o and _q1(0) = -3.31, -6.24, -8.41 rad/s . . . . . . . 122
5.19 Angle q1 and angular velocity _q1 results of the cantilevered elastic compound
pendulum for q1(0) = 45o and _q1(0) = -2.66, -6.47, -9.06 rad/s . . . . . . . 123
5.20 Displacement q2 and velocity _q2 results of the cantilevered elastic compound
pendulum for q1(0) = 45o and _q1(0) = -2.66, -6.47, -9.06 rad/s . . . . . . . 124
5.21 Angle q1 and angular velocity _q1 results of the cantilevered elastic compound
pendulum for q1(0) = 61:5o and _q1(0) = -2.70, -6.54, -9.17 rad/s . . . . . . 125
5.22 Displacement q2 and velocity _q2 results of the cantilevered elastic compound
pendulum for q1(0) = 61:5o and _q1(0) = -2.70, -6.54, -9.17 rad/s . . . . . . 126
6.1 Free link applied for experiments . . . . . . . . . . . . . . . . . . . . . . . . 130
6.2 Compound pendulum applied for experiments . . . . . . . . . . . . . . . . 130
6.3 Motion measurement system Optotrak 3020 . . . . . . . . . . . . . . . . . 131
6.4 Experimental and simulation results of the cylinder type rigid link for q(0) = 0o 134
6.5 Experimental and simulation results of the cylinder type rigid link for q(0) =
32o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
6.6 Experimental and simulation results of the cylinder type rigid link for q(0) =
55o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
6.7 Experimental and simulation results of the compound pendulum for q(0) = 22o 139
6.8 Experimental and simulation results of the compound pendulum for q(0) = 31o 140
6.9 Experimental and simulation results of the compound pendulum for q(0) ? 45o 141
6.10 Experimental and simulation results of the compound pendulum for q(0) =
61:5o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
6.11 Stopping time results of the compound pendulum . . . . . . . . . . . . . . 143
xiii
Chapter 1
Introduction
A granular material is a simple mixture of solid particles. Grains such as rice, soils
including sand, and artiflcial granules such as fertilizer, glass beads, ball bearings, and pills
can be classifled as granular materials. Even snow is considered as a granular medium in
the study of avalanche phenomenon. Powder is a special class of granular material. It is
composed of much flne particles compared with the usual granular material. Due to its
characteristics, inhomogeneous, disordered, and anisotropic on a microscopic scale [1], the
powder can ow much easily when agitated.
An individual granule is solid and shares the physical properties of solid matter. A
granular in its conglomerated state acts like a solid in the case when the kinetic energy
of individual grains is low and the grains keep a stationary state. However individual
grains begin to loose the stationary state and the granular material begins to uidize acting
liquid-like in case when the external driving forces including tilting and shaking exceed the
stationary conditions. Roughly speaking the granular medium may ow such as ordinary
uids. According to circumstances, the granular material may act even like a gas when the
external driving forces act intensively and the contacts between the grains become highly
sporadic.
The flrst approach of the granular medium was based on the contact and the collision of
solid particles [2, 3, 4, 5, 6, 7, 8, 9, 10]. The other approach was based on uid mechanics and
hydrodynamics approach taking the uid-like characteristics of the materials into account
[11, 12, 13, 14, 15].
Even though granular materials show similarities, they exhibit unusual behaviors com-
pared with solids, liquids, or gases. For example, when granular materials are oscillated by
shaking or vibrating containers, they become inhomogeneous in space and time. They are
1
segregated and reorganized [16, 17, 18, 19, 20, 21, 22]. And the compaction happens also by
the vibration of the package container [23, 24, 25, 26, 27, 28]. In some cases, the pressure of
a granular material located in tall containers is non-linear with the depth of the material due
to the contact force generated between granules or between granules and the container walls
[29]. These contact forces constitute a network of forces at large scale [30, 31, 32, 33, 34, 35].
Granular materials show also complicated behaviors such as formation of arches, convection
rolls, pattern formation, and dynamical instabilities [36].
In industry, the granular materials are the second-most manipulated material (the flrst
one is water) [36]. Not only most raw materials exist as the form of granules but also many
flnal products are manufactured as granular materials. Therefore the study of this material
can be applied to industries such as food, pharmaceutical, mining, chemical, semiconductor,
agriculture, and construction. Disasters in nature such as landslides and avalanches show
the importance of the study of the granular medium.
The impact with a granular material is also an interesting subject in the fleld of engi-
neering. The impact with a granular material is not a simple problem because the phenom-
ena of contact, collision, and ow simultaneously happen when a body impacts the material.
Of course, the major phase state of granular material is decided by the penetrating velocity
of an impacting body in the granular material. In the case of high speed impact, the char-
acteristic of the granular material in the vicinity of a rigid body is similar to a uid. The
granular material acts like a solid when the penetrating speed of the body is slow. However
in usual impact cases the behaviors of the granular materials exhibit the combined form of
these characteristics.
Most early studies of the impact penetrating problems started from high speed impact
and were motivated by military applications. Benjamin Robins carried out experiments in
gunnery and wrote \New Principles in Gunnery" (1742), which contains a description of
his ballistic pendulum. His work was translated into German by Leonhard Euler (1745)
[37] and Euler added annotations of his own. Jean-Victor Poncelet developed Poncelet?s
model in his book \Cours de Mecanique Industrielle" (1829). This model is considered as
2
the most venerable and classical empirical formula to calculate the penetration resistance
of a target stricken by a at-nosed projectile [38].
The study having begun from high speed impact area attracts attention recently again
at relatively low speed regime in granular physics. Penetrating granular medium with a
rigid body at very low speed and the impact cratering including the impact penetration
are the representative areas. In the area of penetrating granular medium at low speed,
the horizontal resistance force by granular materials [39], jamming and uctuations of the
resistance force [40, 41], shape efiects on the resistance force [42], and the vertical resistance
force [43, 44] were studied. In the case of the impact cratering, the size, the depth, and
the form of the crater with the initial impact conditions of a rigid body were studied
[45, 46, 47, 48, 49, 50, 51, 52]. To develop a force law model for the granular impacts and
to flnd a mathematical formula in order to measure the impact force of objects dropped
represent a new interest in this fleld of the study [51].
The resistance force models, linear to the depth of a body [50, 53], linear to the velocity
of a body [46], and linear to the square of velocity of a body [49] have been studied to explain
the motion in the granular materials. Tsimring and Volfson studied the impact cratering
by penetration of large projectiles into a dry granular medium [48]. In their study they
proposed for the resistance force model the sum of a velocity dependent drag force and
a depth dependent resistance force. A depth dependent resistance force (static resistance
force) of a rigid body in a granular matter has been studied in granular physics prior to
their study. The static resistance force model has been developed for the horizontal motion
in [39, 40, 41, 42] and for the vertical motion in [43, 44].
The models of the resistance force impeding the motion of a rigid body in a granular
materials and high speed motion capture cameras enabled to simulate and to conflrm the
velocity of the bodies. Ambroso, Kamien, and Durian [49] studied the time dependance of
a rigid sphere impacting and concluded the stopping time taking for a rigid body to stop
from the moment of impact is a function of the geometry and the initial impact velocities.
Hou, Peng, Liu, Lu, and Chan [50] studied the deceleration of projectiles impacting with
3
the medium and they concluded the stopping time is not a linear function of initial impact
velocity. The paper of Katsuragi and Durian has sparkled new interest in the fleld of impact
with a granular matter regarding the stopping time of a rigid body. They introduced the
resistance force model proposed by [48] to study the impact of a rigid sphere and verifled
with a line-scan digital CCD camera [51]. They analyzed an interesting phenomenon, how
rapidly a sphere impacting a granular medium slows upon collision and they clarifled the
relation between the stopping time of a rigid body and the initial impact velocity in a
vertical impact situation. Analysis shows that as the speed at which the spheres impact
the medium increases, the sooner it will come to rest for the vertical impact. Lee and
Marghitu extended the study to the model of a rigid body obliquely impacting the medium
and reported similar results even in the case of oblique impact [54].
As mentioned, granular materials are ubiquitous. Hence the impact with a granular can
take place in various granular material environments. This means its study can be applied
in numerous areas such as robotic locomotion, tracked vehicles, and heavy duty construction
equipments. Especially military operation environments are closely related with this study
because in general the military operations happens in outflelds such as soil, sand, mud,
or mixed environments. Therefore developing or planning to develop multi-legged military
robots for surveillance or carrying cannot avoid the continuous impact with the granular
materials. From this view point, the study of the impact with granular medium can be
applied to development including the design, the manufacture, and the optimization of the
operations of these multi-body systems.
In this study we focused on modeling and simulation of kinematic chains impacting a
granular medium using the resistance force model verifled by Katsuragi [51] as the sum of a
velocity dependent drag force and a depth dependent resistance force. We also analyzed the
relation among initial impact velocities, stopping time, and penetrating depth based on the
experimental and the simulation results. To our best knowledge this is the flrst time when
a mathematical model is proposed, analyzed, and experimentally verifled for the impact of
kinematic links impacting into a granular matter.
4
Chapter 2
Description of the problem
The purpose of this study is to analyze the impact and penetrating motion of a particle
and kinematic chains into a granular material. As introduced, the impact of spheres and
cylinders with a granular material have been studied by various researches from high to low
speed impact taking into account the depth and the width of the crater. However, most
studies were restricted in one vertical direction. There had not been studies regrading the
general impact conditions such as oblique impact and penetration with low initial impact
velocity although almost all real impact problems happen under this oblique condition.
In this study, the motions of impact bodies are restricted to planar motion. In order
to analyze the collision efiect, the models of a single impact point and multiple impact
points are chosen as analytical models. The single impact point model includes the oblique
impact of a sphere, a mathematical pendulum, and single impact of a planar kinematic
chain. For multiple impact points model, a planar kinematic chain having multi contacts
with a granular medium is studied. The impact of a planar kinematic chain is separately
studied as a rigid and as a exible model.
Basically the most important component in modeling the impact with a granular mate-
rial is to understand the resistance forces originated from the granular material. There had
been studies about granular physics in order to explain the unusual characteristics of gran-
ular materials at high and low speeds. However there are not enough researches studying
the problems of resistance forces impeding the motion of a rigid body in a granular mate-
rial. In this study, the resistance forces afiecting the motion of a body are studied based
on the characteristics of a granular material acting uid-like and solid-like. Based on these
5
characteristics of a granular material, the resistance forces are studied as the dynamic fric-
tional force (velocity-dependent force) which is modeled from uid dynamics and the static
resistance force (depth-dependent force) which is modeled from solid particle contacts.
There are various mathematical modeling methods such as the equations of Newton-
Euler, Lagrange, and Kane. In this study, the descriptions of the dynamics of the system of
particles and kinematic chainsare formulatedusing the Newton-Euler equations. Lagrange?s
and Kane?s equations can be also applied and were tested but there are no difierences in
the numerical results using difierent methods.
Initial impact conditions are decisive factors for the motion of the particle and the
kinematic chain during penetration phase in the case that there are no drive force and no
moment to move and to rotate the links. Of course, the initial states of a granular material
including density and volume fraction are also a very important components because the re-
sistance forces strongly depend on them. However the initial states of the granular material
utilized in experiments and simulations are assumed not to change during the experiments
in order to investigate the efiects of the initial impact conditions. An actual granular ma-
terial test bed was prepared as equally as possible at every impact experiment. Under this
assumption, the most important initial conditions for impact problem are the angle and the
velocity at the impact moment. The experiments and the simulations are performed under
these conditions of various initial impact angles and velocities and the results are compared.
The simulations of the models are basically performed utilizing NDSolve function of
Mathematica and this software is considered as su?cient to solve the O.D.E.s formulated
by Newton-Euler equations. Due to the characteristics of the resistance forces acting always
opposite to the direction of motion of an impact object, the simulation should solve repeat-
edly, should start and stop when the motion of the rigid body changes direction until the
rest of the rigid body. The simulation are performed with variable and flxed step methods.
The experiments are performed several times for a certain initial impact angle and
velocity. At every experiment, the granular material of test bed is plowed and atted
to provide almost the same state of the material and to avoid ageing efiects arising in
6
repeated experiments. The data attained from the repeated experiments are averaged to
minimize the efiects of abnormalities such as uctuation phenomenon which can occur
during experiments. Experimental data are applied to flnd the resistance force coe?cients
and the experimental results will be compared with the simulation results to verify the
impact models. As an experimental equipment, NDI (Northern Digital Inc.) optotrak
system is utilized for the motion capturing of impact objects and the digitizing of the
captured data. The system can capture the position of markers within the RMS accuracy
of 0.1mm and can track up to 256 markers. The position of the impact body was captured
with 3D data cartesian coordinate system with 500 frames/s.
The impact and penetrating motions of a particle and kinematic chains obtained from
simulations and experiments will be analyzed and compared based on the initial impact
conditions. We analyze how the initial impact conditions afiect the penetration including
the depth and transition of the velocity after impact. Especially, the efiects of the initial
impact conditions, deceleration and stopping time from the moment of impact until the rest
of an impact body will be evaluated. Until now only the vertical impact of a sphere type
was reported [51] and our study will evaluate the efiects of initial impact angle including
oblique impact for rigid and elastic links.
2.1 Coordinate system
The motion of the considered systems are in 2D (planar motion). For describing the
motion of the impact object composed of rigid bodies, the system utilizes a flxed frame
because it is su?cient to describe the planar motion of the models with a flxed cartesian
coordinate system. The x-axis is horizontal, with the positive sense directed to the right,
the z-axis is vertical, with the positive sense directed downward, and y-axis is perpendicular
to the plane of motion as shown in Fig. 2.1. The planar motion of an impact rigid body
will be in the xz plane. The origin of the coordinate system will be flxed. The unit vectors
for the x-axis is ?0, the y-axis is ?0, and the z-axis is k0. The angle between the z-axis
and the link i of kinematic chain consisting of n interconnected rigid links is denoted by qi
7
q1
q2
qn
Granular medium
qz
qx
O
00
0 x
zy
k
?
?
Figure 2.1: Coordinate system for rigid impact bodies
k0
0?
0?
(0) O
q1
k
?
(1)
Granular medium
Figure 2.2: Coordinate system for exible impact bodies
8
(i = 1;2;3;:::;n). The displacement from the origin to a certain point is denoted by qx
in x-axis and by qz in z-axis. For describing the motion of kinematic chain composed of
exible bodies, the system utilizes a flxed frame and a rotating reference frame as shown in
Fig. 2.1. The motion of the exible bodies are described in a flxed reference frame (0) of
unit vectors [?0;?0;k0] and a mobile (rotating) reference frame (1) of unit vectors [?;?;k].
2.2 Description of the models
2.2.1 Impact of a rigid sphere
A general representations of a sphere and a mathematical pendulum impacting a gran-
ular material are shown in Figs. 2.3 and 2.4. The impact is initiated when the sphere (mass
ms and diameter ds) strikes the surface of a granular material. In this study, its diameter ds
is considered relatively small and then its motion is assumed to be the motion of a particle.
The initial impact conditions are represented as initial linear velocity and initial impact
angle. The motion of the sphere will be observed based on the initial impact conditions.
In particular, the impact by the sphere does not change the direction of the velocity after
the impact moment but the mathematical pendulum can changes the direction of the ve-
locity on z direction when the sphere passes through vertical axis. Therefore the impact by
the mathematical pendulum may lead to several outcomes depending on initial impact the
angles and the velocity.
2.2.2 Single impact of a planar kinematic chain
Impact of a free link
The link utilized for the impact model has the length L, the mass mc, and the diameter
dc as shown in Fig. 2.5. The impact is initiated when the end T strikes the surface of a gran-
ular material. This problem is studied based on the continuum modeling and additionally a
exible model is simulated in order to compare the exible elastic body and the rigid body.
The initial conditions are given by the linear vertical impact velocity and angular position.
9
Granular medium
ds
v(0)
q(0)
O
x
zy
00
0
k
?
?
Figure 2.3: Oblique impact of a sphere
q(0)
ds
?q(0)
x
zy
Granular
medium
L
O
massless rod
00
0
k
?
?
Figure 2.4: Impact of a mathematical pendulum
10
L
T
A
q
0
k0
0?
?O
Granular medium
C
dc
qz
qx
Figure 2.5: Impact of a free link
Even though there is no tangential force at the initial impact moment, after impact the
horizontal resistance force is generated by the rotating link.
Single impact of a planar kinematic chain with n links
A schematic representation of a kinematic chain consisting of n interconnected rigid
links 1;2;3;:::;n is shown in Fig. 2.6. The kinematic chain has the rotational joints at
Ai (i = 0;1;2;:::;n ? 1) and each link i has the length Li, the mass mi, and the same
diameter as dc. The impact is initiated when the end last T of the link n strikes the surface
of a granular material. The initial combination of the linear and rotational impact motion
are considered as initial impact conditions. The impact may lead to several consequences
depending on initial conditions such as the initial impact angles of the links and the velocity
of the tip at the impact moment as well as the characteristics of a granular medium such as
resistance force coe?cients. The case when only the last link n is penetrating the granular
medium is considered and modeled in this study. This problem is also studied as an elastic
model and a rigid model for the special case n = 1 (compound pendulum).
11
A1
A2
A0
1
2
T
Ai?1 A
i
-1An
i
n
Granular medium
zy
q1(0)
qi(0)
qn(0)
q2(0)
q1? (0)
iq? (0)
nq? (0)
q2? (0)
x
00
0
k
?
?
Figure 2.6: Single impact of a kinematic chain
2.2.3 Multiple impacts of a planar kinematic chain with n links
A schematic representation of a kinematic chain consisting of n interconnected rigid
links 1;2;3;:::;n and l joints Ai (i = 1;2;3;:::;l) is shown in Fig. 2.7. Each link i has
the length Li, the mass mi, and the same diameter, dc. Immediately before impact the
end links of the kinematic chain can be resting on the surface surfaces, impacting the
surface, or not interacting with the surface. The impact is initiated when a speciflc end
Tj (j = 1;2;3;:::;m) or all ends of the chain strike the surface Sj. The angle of the
surface with the horizontal is j. The combination of the linear and the rotational impact
displacements are considered as initial impact conditions. The impact initiated at Tj may
also lead to several outcomes depending on the initial impact conditions as well as the
characteristics of a granular medium.
12
T1
A1
A2
T2
T3
Tm
Al
?m
?3
?1 ?2
1
2
qz
qx O
n
q1(0)
qn(0)
q2(0)
q3(0)
q1? (0)
nq? (0)
q2? (0)
x
zy
Ai
i
qi(0) iq? (0)
00
0
k
?
?
Granular
medium
Figure 2.7: Multiple impacts of a kinematic chain
2.3 Generalized coordinates
In the case of the oblique impact of a sphere, the impact and the penetrating motion of
the rigid body is restricted to xz plane. The numbers of degrees of freedom for the oblique
impact of the sphere is 2. Hence the generalized coordinates of the model can be expressed
by displacements in the xz plane as
q = fqx;qzgT ; (2.1)
where qx and qz represents the displacements from the origin to the mass center of the
sphere in x-axis and in z-axis respectively.
In order to describe the motion of a rigid free link impacting and penetrating a granular
medium, at least three position coordinates are required even though the motion is restricted
to xz plane. Two displacement coordinates qx and qz and one angular position coordinates
13
q can be considered as generalized coordinates of the model as
q = fqx;qz;qgT : (2.2)
For a exible free link, the generalized coordinates include the coordinates of the rigid free
link and additionally the elastic coordinates because the position of each point on the link
depends on the deformation of the link which is represented by inflnity superposition of
deformation. The generalized coordinates can be given as
q = fqx;qz;q1;q2 :::q1gT : (2.3)
In the model of a single impact of kinematic chains, the impact and the penetrating
motion of each link is translational and rotational simultaneously (Fig. 2.6). Each link i
has one relative degree of freedom with respect to the link i?1, and therefore the number
of D.O.F. is n. Angles, qis, are changing with time at the instant of interest and even
translational motion of each link can be described by the angles. Therefore the angles qis
are appropriate generalized coordinates describing the motion of the kinematic chain and
the generalized coordinates can be expressed by n?1-dimensional vector as
q = fq1;q2;q3;:::;qngT: (2.4)
The mathematical and the compound pendulum are particular cases with n = 1. However,
for the exible compound pendulum, the generalized coordinates can be represented as the
rigid and elastic coordinates as
q = fq1;q2;q3;:::q1gT : (2.5)
For multiple impacts of a planar kinematic chain, the number of D.O.F. is n + 2
(Fig. 2.7). The generalized coordinates, qx and qz, are the displacements form the origin to a
14
certain particular point. The angles, qi, and the displacements, qx and qz, are changing with
time at the instant of interest and can be considered as generalized coordinates describing
the motion of the kinematic chain. The generalized coordinates can be expressed by (n +
2)?1-dimensional vector as
q = fqx;qz;q1;q2;q3;:::;qngT: (2.6)
15
Chapter 3
Force analysis
For the impact with penetration into a granular medium, the most important inter-
action forces between the body and the medium is the resistance force of the medium on
the body from the moment of impact until the body stops. Even after the body stops,
the reaction forces between the kinematic chain and the granular medium exist due to the
gravity force and the static resistance force.
There are various models of resistance forces that act upon the body during penetration
process. The Bingham?s model is deflned as FR = F0+bv and has been utilized in the case
when the drag is viscous. The model is originated for a viscoplastic material that behaves
as a solid at low stresses but ows as a viscous uid at high stress. Bruyn and Walsh
applied this model and calculated the penetration of spheres into loose granular material
[46]. The Poncelet?s model FR = F0 + cv2 is originated from high speed impact and has
been applied in the situation when the drag is inertia, where v is velocity of the body
relative to a medium, b and c are drag coe?cients. The Bingham?s model has recently been
advocated for granular impact, while the Poncelet?s model has long been used for ballistics
applications [49].
Allen, Mayfleld, and Morrison [55, 56] suggested the resistance force model including
the Bingham?s and the Poncelet?s model as FR = Av2 +Bv +C in the study of projectile
penetrating sand, where A is a classical uid dynamic drag parameter, B is a deceleration
parameter due to kinetic flction on the surface of the projectile, and C is the deceleration
force due to the inherent structural characteristics of the target material.
However, some recent research results [39, 40, 41, 42, 43, 44] indicate that the static
resistance force, deflned as a constant for the Bingham?s and the Poncelet?s models, depends
on the depth of penetration linearly or non-linearly. Lohse, Rauhe, Bergmann, and van der
16
Meer formulated the resistance force model made up of the only static resistance force
[53]. They suggested the resistance force as FR = kz, where k is a constant and z is the
penetrated depth. Tsimring and Volfson [48] proposed a resistance force model based on
the Poncelet?s model and the static resistance force model studied in [39] as a generalized
Poncelet?s model. The total resistance forces acting on a moving body into a granular
matter can be generalized as the sum of the static resistance force characterized by a
depth-dependent friction force as well as the dynamic frictional force characterized by a
velocity-dependent drag force as
FR = Fs(z)+Fd(v); (3.1)
where Fs is the static resistance force and Fd is the dynamic frictional force. With this
approach the resistance force model can utilize the solid-like and uid-like characteristics
simultaneously.
Anther possible force impeding the body penetrating a granular medium is generated
by air resistance. In some case, air drag force plays an important role adding the disturbing
force to an object moving such as the motion of a car. The air drag force depends on the
density of the medium and on the velocity of the moving object. However the efiect of this
resistance force in the problem of the impacting with a granular is not important because
the simulation results show ignorable difierences.
3.1 Dynamic frictional force Fd
The dynamic frictional force Fd is conceptually the same force as \drag" widely used
in the fleld of uid dynamics. Even though the state of a granular matter is not a uid but
a solid, individual grains begin to loose stationary state in their contacts and a granular
material begins to uidize acting liquid-like when external driving forces including tilting
and shaking exceed the condition of stationary. From this uid-like behavior of a granular
17
matter, the resistance force is assumed to be as a drag force impeding the body motion
[57, 58, 59, 60, 61].
In case of the impact with a granular material, once a body penetrates through the
granular matter after impact, the matter contacting with the moving body partly acts like a
uid and this resistance force begins to be generated by the relative movement of the body
with the granular material and acts in a direction opposite to the instantaneous relative
velocity of the body. Hence, in the case of impact with the granular material, the force is
generated from the moment of impact until the rest of the body. The dynamic frictional
force acts no more as the resistance force impeding the motion of a body when the body
and the granular material keep stationary state. The force is a velocity-dependent force
unlike ordinary resistive friction force such as the static resistance force which depends on
depth of body in a granular material.
Like the drag force model in uid dynamics, the dynamic frictional force Fd also had
been modeled as a linear drag force or a quadratic drag force. Therefore the Bingham?s
model utilizing a linear drag force model had been advocated for low speed granular impact,
while the Poncelet?s model utilizing a quadratic drag force had long been used for ballistics
applications [49]. However the research results using a dilute granular ow condition which
is not afiected by the static resistance force [57, 59, 61] and the research results at low speed
impact [48, 51] show that the quadratic drag force model is more suitable for the dynamic
frictional force Fd than the linear drag force model even at relatively low speed regime. The
dynamic frictional force can be modeled as
Fd = ? vjvj ?d ?g Ar v?v; (3.2)
where ?v=jvj term is applied because the dynamic frictional force acting on the opposed
direction of the velocity vector v, ?d is a drag coe?cient determined from experimental
results, ?g is the density of the granular medium, and Ar is the reference area of the body.
18
3.2 Static resistance force Fs
The static resistance force is an internal resisting force and appears when an external
force is applied to a body. This resistance force does not depend on the velocity or the
acceleration of a body and acts as an internal stress. Therefore this force acts as a main
resistance force when a body penetrates a granular medium at low speed.
The existence of this force not existing in a uid can be observed easily. When a body is
put on the granular material, the body starts to penetrate vertically through the granular
material but does not sink beyond a certain depth. The body stops to penetrate in the
material due to this resistance force unlike an usual uid in which a body sinks unlimitedly
without a lift force such as buoyancy.
Even though this resistance force play an important role in a body penetrating a gran-
ular material, the force had not been highlighted comparatively in the early studies. Most
models including the Bingham?s and the Poncelet?s models have considered this force as a
mere constant. Because the dynamic frictional force originated from a velocity-dependent
drag force acts as a main resistance force in the study at relatively high speed impact, the
static resistance force was modeled using a simple form.
However, the theoretical modeling of this force originated from the characteristics of
a granular material even in stationary state is not that easy to be determined. Various
research results explain the origin of this resistance force as the network of forces generated
from contacts between granules within the material pile [30, 33, 31, 34, 32, 35].
In the case of impact with a granular material problem, the force acting on a body is not
generated only by interacting between granules in contacts. Describing interacting forces
of a certain moment during the continuous penetrating process, the force interacts between
a penetrating body and granules located on the body circumference and the interacts of
the contacting forces among the granules in the vicinity of the body are expanded to other
granules of a granular container. As a consequence, the interacts of the contacts between
granules become to form a force chain in a granular containing bed as bulk scale.
19
The reason this resistance force is not easy explained in the impact penetrating problem
is that the force is transmitted following consecutive granular medium movement propaga-
tion such as chain reaction beginning with serial displacements of granules placed at the path
of the moving body and ending with large scale reorganization of granules in the container.
That means the forces and their spatial correlations, speciflcally in response to forces at the
granular system boundaries do not represent flxed structure. Force chains are connected
newly or disconnected continuously until the rest of penetration. This process is strongly
in uenced by properties of a granular matter, the external shape of a body, packaging
state of medium, even the form of a vessel containing medium [41, 62]. Any insigniflcant
changes during the continuous penetration can increase uncertainties of the propagation
chain. These uncertainties in a real impact problem make the forces not to propagate
uniformly through the granular matter but localized along directional force chains. The
inhomogeneous force propagation and the irregular grain reorganization cause the uctua-
tion of the static resistance force in direction and magnitude [41, 42, 63, 64, 65]. Due to
this complicated process of transmitting the force, there are few general continuum theories
completely describing the static resistance force [44]. In this study, we considered horizontal
and vertical static resistance forces developed by theoretical and empirical approaches.
3.2.1 Horizontal static resistance force Fsh
Horizontal static resistance force is deflned as an internal impeding force acting on the
horizontal direction. When a body penetrates a granular material, it is not easy to separate
and to measure this force individually without the efiect of dynamic frictional force Fd. The
experiments related with the horizontal static resistance forces were performed at very slow
speed such as 0.04 - 1.4 mm/s [39, 40, 41].
Albert, Pfeifer, Barab?asi, and Schifier [39] applied a probability approach to model
this resistance force. They used the force chain model developed by Coppersmith, Liu,
Majumdar, Narayan, and Witten [30] in order to model the mean horizontal static resistance
force. According to the chain model of Coppersmith et al., the force chains transmitting the
20
propagation of contact forces are built up based on the force fractions that act on a given
particle of the matter. The force fraction is unequally split and transmitted to the adjacent
particles. The traces of these force fractions are built as the form of chain in the granular
material pile. For a particle to move horizontally relative to a granular matter which is
restricted by the force chain, the probability of the particle immersed at a certain depth
position, exceeding a critical force necessary to make the granule slip relative to another,
should be 1 when a certain external force is applied horizontally on this particle. Expanding
this concept of small particle to a body, the probabilities of all granules contacting with a
body, exceeding a critical force necessary to make the granules slip relative to others, should
simultaneously be 1 at any depth position when a certain external force F is applied on a
body horizontally for a body immersed in the granular matter to penetrate horizontally.
The calculated minimum force F satisfying this probability condition of a vertically
immersed cylinder type body is F = ?h g?g H2 dc, where g is the gravitational acceleration,
H is the immersed length of the cylinder, dc is the diameter of the cylinder, and ?h is a
constant depending on the medium properties such as surface friction, morphology, packing
of the grains, etc. [39]. Using the same process, the horizontal static resistance force of the
cylinder including the state of immersion, at any slope, can be generalized as
Fsh = ? vhjv
hj
?h g?g z2T dc ; (3.3)
where vh is the horizontal velocity vector of a body and zT is the depth of the immersed
cylinder tip. For a sphere, the horizontal static resistance force is
Fsh = ? vhjv
hj
?h g?g ds2 zT ; (3.4)
where ds is the diameter of the sphere. Equations (3.3) and (3.4) show that the horizontal
static resistance force is a function of the granular properties, the immersed surface area,
and the depth.
21
Even though Eqs. (3.3) and (3.4) do not re ect the periodical force uctuation such
as stick and slip phenomena which were conflrmed in some real cases, experimental data
show that the equation can be applied for the static resistance force model not only in the
case of little uctuation phenomenon but also for large uctuation because the calculated
force represents the average static force of uctuation [41, 42]. In addition, this modeling
approach reveals that inertia of the body and the friction force between granules and the
body surface in contact have negligible contribution to the horizontal static resistance force.
In other words, any body only with the same shape can cause almost the same results
regardless of the body?s properties such as density and surface friction.
3.2.2 Vertical static resistance force Fsv
The vertical static resistance force is deflned as an internal impeding resistance acting
on the vertical axis. Although most researches regarding the impact with agranularmaterial
focus on the vertical direction, the studies about the static resistance force in this direction
are not enough or are limited. Simple models of this force applied as a mere constant [46, 49]
and linear function of the immersed depth of a body in a restricted range [51, 53]. However,
basically this force is known as a nonlinear function of the immersion depth [43, 44].
The reason that there are very few continuum models for this force in a granular
system is originated from complex force distribution. Even simple experimental result
shows the contacts of granules increase exponentially by the force loaded externally [31, 34]
and especially the efiects by the container bottom boundary increase the non-linearity of
this resistance force [43, 62]. In addition, the problem of the direction loaded by external
force makes the modeling of this force even more di?cult.
Hill, Yeung, and Koehler [44] suggested an empirical equation with coe?cients calcu-
lated from the experimental data as
Fsv = ? vvjv
vj
?v (zT=l)? g?g V ; (3.5)
22
where vv is the vertical velocity vector of a body, V is the immersed volume of the body,
and l is the lateral dimension. The coe?cients ?v and ? are depending on the shape of
the body, the properties of the granular matter, the shape of medium container, and the
moving direction such as plunging and withdrawing.
When Eq. (3.5) is applied with the coe?cients ?v and ? to a dynamic model, the
noticeable conditions for modeling are the inclination and the moving direction of the body.
The coe?cients based on the empirical data [44] reveal that the inclination of a body has
little efiect on the vertical static resistance force however this force changes drastically by
the moving directions. For the cylinder type body, whether the axis is vertical or horizontal,
the coe?cients ?v and ? are approximately 10 and 1.4 for plunging motion and 0.5 and 1.7
for withdrawing motion in their experiments.
3.3 Air drag force
Air can also act as the medium impeding the motion of an impacting body besides
the granular medium. For some cases, the air drag force acts as a main resistance force
disturbing the motion. However, the air drag force that can act as an additional resistance
force acting on the impact object is not considered in this study because the force is too
weak compared with the resistance force of the granular medium. The simulation of a
moving compound pendulum shows the difierence between the air drag force applied to the
model and not applied to the model is as small as possibly to be ignored. The calculated
time of the force applied to the compound pendulum model for the pendulum released at
45o to pass through vertical plane is 0.164954 s and that of the model not applying the force
is 0.164932 s. The difierence between these times is only 0.000022 s. The air drag force is
calculated as F = 1=2Cd Ar ?air v2C where Cd is the drag coe?cient, ?air is the density of
air, and vC is the velocity of the center of the mass [66]. As the condition of air for this
simulation, the air density ?air = 1.2045 kg/m3 (20?C and 1 ATM) and the drag coe?cients
Cd = 0.93 for cylinder [66] are utilized. The applied dimensions of the compound pendulum
are the length 0.15 m, the diameter 0.00635 m, and the density of the steel 7.7?103 kg/m3.
23
Tim e (s)
0.00 0.05 0.10 0.15
0.00000
0.00005
0.00010
0.00015
Dra
g
forc
e
(N)
Figure 3.1: Air drag force impeding the motion of the compound pendulum
Figure 3.1 shows the air drag force acting at the compound pendulum. The reason that the
air drag resistance is small even compared with the dynamic frictional force is the difierence
between the density of the air and the granular medium (sand). In this study, the density
of the sand applied as a granular medium is 2500 kg/m3. Considering the weakness of the
air drag resistance force, only resistance force by the granular medium is assumed to act
against an impact in this study.
3.4 Application to the impact models
In the case when the impact object is considered as a particle, there is no restriction
to apply Eqs. (3.2), (3.3), (3.4), and (3.5) because the acting points of the dynamic fric-
tional force and the static resistance force are consistent with each other. In addition the
magnitude and the direction of the dynamic frictional force depend only on that acting
point. However, in the case when the impact object is regarded as a continuum, there can
be difierences between the application points of the resistance forces (between the dynamic
frictional force and the static resistance force). It is also di?cult to decide the magnitude
24
and the direction of the dynamic frictional force because the magnitude and the direction
of the velocities of points on the impact object are not the same.
In this study, the centroid point of the immersed part of the body is assumed to be
application point for the resistance force. Under these assumptions, Eqs. (3.2), (3.3), (3.4),
and (3.5) are newly rewritten as
Fd = ? vEjv
Ej
?d ?g Ar vE ?vE = ?vE ?d ?g Ar jvEj; (3.6)
Fsh = ? vExjv
Exj
?h g?g z2T dc = ??sign(vEx)?h g?g z2T dc? ?0; (3.7)
Fsh = ? vExjv
Exj
?h g?g d2s zT = ??sign(vEx)?h g?g ds2 zT? ?0; (3.8)
Fsv = ? vEzjv
Ezj
?v (zT=l)? g?g V =
h
?sign(vEz)?v (zT=l)? g?g V
i
k0; (3.9)
where vE is the velocity vector of the centroid point of the immersed part of a body,
vEx represents the velocity vector composed of the horizontal component of vE, and vEz
represents the velocity vector of the vertical component of vE. Equations (3.7) and (3.8)
represent the horizontal static resistance force for a cylinder type body and for a sphere
respectively. With this approach, the resistance force FR is applied as FR(vE;zT) in the
equations of motion.
When the impact object is characterized as a particle, the terms of the resistance force
are calculated easily. The velocity vector of the resistance force application point E, vE,
is the same as the velocity vector of the mass center C, vC. The reference area Ar is
calculated as ?d2s=4 for sphere type objects. The immersed depth zT is the same as the
vertical coordinate of the mass center C and the immersed volume V is calculated as ?d3s=6
for the sphere type objects and d3s for rectangular typed objects. In this case, the term of
lateral dimension l for calculating the vertical static resistance force is ds. Therefore the
25
E
qm
q
Granular medium vEx
Ar
vEzvE
zTz
T /2
dc dc
Figure 3.2: Moving angle, qm, of a link in planar motion
resistance force of a small sphere type object is calculated as
FR = Fd +Fsh +Fsv
=
"
?vCx ?d ?g ?d
2s
4
q
v2Cx +v2Cz ?sign(vCx)?h g?g q2z ds
#
?0
"
?vCz ?d ?g ?d
2s
4
q
v2Cx +v2Cz ?sign(vCz)?v (qz=ds)? g?g ?ds
3
6
#
k0; (3.10)
where qz is the scalar value of the vertical position vector of the mass center C.
For a link in planar motion (see Fig. 3.2), the position and the velocity vector of the
resistance force application point E should be deflned in order to calculate the resistance
force of link type objects. The reference area Ar is calculated as
Ar = dc zTcosq jsin(q ?qm)j; (3.11)
qm = tan?1
v
Ex
vEz
?
; (3.12)
where q is the angle between vertical plane and the object and qm represents the moving
angle of the link penetrating the granular medium as shown in Fig. 3.2. The term, lateral
26
dimension l, for calculating the vertical static resistance force for this case is dc and the
immersed volume V is calculated with
V = ?dc
2
4
zT
cosq: (3.13)
Therefore the resistance force of a cylinder type link is calculated as
FR = Fd +Fsh +Fsv
=
"
?vEx ?d ?g dc zTcosq
flfl
flflsin
q ?tan?1
v
Ex
vEz
??flfl
flfl
q
v2Ex +v2Ez
?sign(vEx)?h g?g z2T ds
#
?0
"
?vEz ?d ?g dc zTcosq
flfl
flflsin
q ?tan?1
v
Ex
vEz
??flfl
flfl
q
v2Ex +v2Ez
?sign(vEz)?v (zT=ds)? g?g ?dc
2
4
zT
cosq
#
k0: (3.14)
27
Chapter 4
Impact of rigid objects
In this chapter, rigid bodies impacting a granular medium are studied. Their equations
of motion (E.O.M.) are formulated and simulated. The resistance forces and the gravity
forces are considered in the E.O.M.s of the system. The E.O.M.s are formulated using
Newton-Euler equations. Other equations such as Lagrange?s equation and Kane?s dynami-
cal equation also can be applied to flnd the E.O.M.s. The simulation results of the E.O.M.s
are the same as Newton-Euler equations.
First, a one degree of freedom (D.O.F.) mathematical pendulum is studied. An oblique
impact of a sphere also studied as the motion of a particle has two D.O.F. After modeling of
the motion of a particle type sphere, the impacts of rigid bodies are studied. Single impact
of a planar kinematic chain including a compound pendulum and a double pendulum are
studied. Additionally the impact of a free rigid link is modeled as the special case of the
single impact of a planar kinematic chain. In order to conflrm the efiect of resistance forces
acting on multiple links simultaneously, multiple impacts of a planar kinematic chain are
modeled.
As mentioned previously the planar motions of the models are described in a flxed
cartesian coordinate system. The x-axis is horizontal, with the positive sense directed to
the right, the z-axis is vertical, with the positive sense directed downward, and y-axis is
perpendicular to the plane of motion. The unit vectors for the flxed frame are ?0, ?0, and k0.
The angle between the z-axis and the link i is denoted by qi, (i = 1;2;3;:::;n). The x and
z-axis distances from the origin to a reference point are denoted by qx and qz respectively.
The simulation is performed using Mathematica and the equations of motion is solved
by NDSolve function. When the flxed step is required to solve E.O.M.s, the applied step
size is 10?6 s. Various initial impact angles are applied as initial conditions in simulations
28
Fsv
x
zy
Granular
medium Fsh
Fd
mg
L
C
q
?q
v
0k
0?
0?
ds
Figure 4.1: Free body diagram of a mathematical pendulum in impact
from low to high impact angle in order to conflrm the efiects of initial impact angles to
the penetrating process from the impact moment until the kinematic chain stops. The
simulation results are basically compared using difierent initial impact velocities from low
to relatively high for the same initial impact angle.
In this study, the penetrating velocity and the penetrated distance are considered as
the basic data for analyzing the impact with a granular medium. The penetrating velocity
is analyzed by comparing the stopping times meaning the time intervals from the impact
moment until the kinematic chains stop.
4.1 Impact of a rigid sphere
4.1.1 Impact of a mathematical pendulum
Modeling
First, the impact of a single mathematical pendulum into the granular material is
modeled. For the mathematical pendulum the net forces acting on the sphere are the
gravity G, the static resistance force Fs, and the dynamic frictional force Fd as shown in
Fig. 4.1.
29
All the forces act at the center of the mass of the sphere of the mathematical pendulum.
In the vertical z-axis, the forces acting on the mathematical pendulum are composed of
the gravity force G, the vertical static resistance force Fsv, and the vertical component of
dynamic frictional force, Fdv. The forces acting in x-axis are the horizontal static resistance
force Fsh and the horizontal component of dynamic frictional force, Fdh. The Newton-Euler
equation is given by
ms L2 ?q??0 =
h
rC ?(G+Fs +Fd)
i
??0 ; (4.1)
where L is the length of the pendulum, q is the angle of pendulum with the vertical, and ?q
is the angular acceleration. The vector rC is the position vector from the origin O to the
mass center C and given as
rC = Lsinq?0 +Lcosqk0: (4.2)
The vector of gravity force G is also represented as
G = ms gk0: (4.3)
All the forces acting on the sphere are assumed to have their application point at the mass
center of the sphere, C. Therefore the velocity vector of the resistance force application
point E and the dynamic frictional force Fd calculated by Eq. (3.6) are represented as
vE = vC = L _q cosq?0 ?L _q sinqk0 (4.4)
Fd = ?vE ?d ?g Ar jvE j = ?vC ?d ?g Ar jvC j
= ?d ?g Ar L2 j _qj
h
? _q cosq?0 + _q sinqk0
i
; (4.5)
30
where the reference area is calculated as Ar = ?d2s=4. The horizontal and vertical static
resistance forces, Fsh and Fvh, are
Fsh = ?sign(_qcosq)?h g?g z2T ds ?0; (4.6)
Fsv = sign(_qsinq)?v (zT=ds)? g?g V k0; (4.7)
where the immersed volume is calculated as V = ?ds3=6 and the immersed depth of the
mathematical pendulum, zT, is calculated as
zT = Lcosq ?Lcosq(0): (4.8)
The initial condition q(0) is the angle of the mathematical pendulum at the impact moment.
The resistance force FR is represented by Eqs. (4.5), (4.6), (4.7), and (4.8) as
FR = Fd +Fsh +Fsv
= ?
"
?d ?g ?ds
2
4 L
2 j _qj _q cosq +sign(_qcosq)?hg?g
Lcosq ?Lcosq(0)
?2
ds
#
?0
+
"
?d ?g ?ds
2
4 L
2 j _qj _q sinq
+sign(_qsinq)?v
Lcosq ?Lcosq(0)
ds
??
g?g?ds
3
6
#
k0: (4.9)
Simulation results
Figures 4.2, 4.3, and 4.4 represent the simulation results of the impact of the math-
ematical pendulum, depicted in Figs. 2.4 and 4.1. The length from the joint to the mass
center of the sphere is L = 0:5 m, the diameter of the sphere is ds = 0:0254 m, its density
is ?s = 7:7?103 kg/m3, and the density of the granular medium (sand) is ?g = 2:5?103
kg/m3. The dynamic frictional force coe?cient ?d = 6.5, the horizontal static resistance
force coe?cient ?h = 8, and the vertical static resistance force coe?cients ?v = 22 and ?
= 1.1 are used in the simulation. The gravitational acceleration g is utilized as 9.81 m/s2
31
Table 4.1: Stopping time of the mathematical pendulum into a granular medium
Impact angle (o) _q(0) (rad/s) ts (s)
-1 0.0642752
15 -3 0.0513544
-5 0.0471527
-1 0.0381696
45 -3 0.0295588
-5 0.0267668
-1 0.0354039
75 -3 0.0267909
-5 0.0240349
in this study. The simulation is performed for the difierent impact angles (q(0) = 15, 45,
and 75o) and the difierent impact angular velocities (_q(0) = -1, -3, and -5 rad/s) from the
impact moment until the angular velocity _q becomes zero.
As shown in Figs. 4.2, 4.3, and 4.4, the penetrating angle increases when the initial
impact angular velocity is increasing. However, the stopping time deflned as the time
interval from the impact moment until the angular velocity of the pendulum becomes zero,
ts, decreases as shown in table 4.1. The increasing of the initial velocity causes the stopping
time into the granular medium to decrease. The faster the mass of the mathematical
pendulum impacts the surface of the granular medium, the sooner it will come to a stop.
This interesting phenomenon involving how rapidly a particle type rigid body vertically
strikes the granular medium slowing down upon contact was reported in [51]. The results
observed only in vertical impact are kept in the impact of the mathematical pendulum
which has both horizontal and vertical components of penetrating motion as the initial
impact velocity v(0) at which a particle vertically impacts the medium increases the sooner
it will come to a stop.
32
0.00 0.01 0.02 0.03 0.04 0.05 0.06
13. 0
13. 5
14. 0
14. 5
15. 0
Tim e (s)
q
(o
)
0.00 0.01 0.02 0.03 0.04 0.05 0.06
-5
-4
-3
-2
-1
0
Tim e (s)
?q(rad/s)
?q (0) = -1 rad/s
?q (0) = -3 rad/s
?q (0) = -5 rad/s
Figure 4.2: Impact results of the mathematical pendulum into a granular medium for
q(0) = 15o
33
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
43. 5
44. 0
44. 5
45. 0
Tim e (s)
q
(o
)
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
-5
-4
-3
-2
-1
0
Tim e (s)
?q(
rad/
s)
?q (0) = -1 rad/s
?q (0) = -3 rad/s
?q (0) = -5 rad/s
Figure 4.3: Impact results of the mathematical pendulum into a granular medium for
q(0) = 45o
34
q
(o
)
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
73. 5
74. 0
74. 5
75. 0
Tim e (s)
Tim e (s)
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
-5
-4
-3
-2
-1
0
?q(rad/
s)
?q (0) = -1 rad/s
?q (0) = -3 rad/s
?q (0) = -5 rad/s
Figure 4.4: Impact results of the mathematical pendulum into a granular medium for
q(0) = 75o
35
Granular
medium
v
sv
Fd
ds
G
F
Fsh
ds
v(0)
q(0)
q
O
x
zy
00
0
k
?
?
Figure 4.5: Free body diagram of a rigid sphere in oblique impact
4.1.2 Oblique impact of a rigid sphere
Modeling
The model of a planar oblique impact of a rigid sphere into a granular matter is
presented. The net forces acting on the oblique planar impact penetrating a granular
medium are the same as the previous case. The gravity force, the static resistance force Fs,
and the dynamic frictional force Fd act as shown in Fig. 4.5. Neglecting the rotation of the
body, Newton?s second law for the sphere gives:
ms?rC = G+Fs +Fd; (4.10)
where ?rC is the acceleration vector of the position vector of the mass center of the sphere,
rC = qx? + qz k, and ms is the mass of the sphere. There is no rebound at the impact
moment regardless the initial velocity and the initial impact angle.
All the forces acting on the sphere are also assumed to have their resistance force
application point at the center of the sphere. In the vertical z-axis, the forces acting on the
sphere are: the gravity force G, the vertical static resistance force, Fsv, and the vertical
component of dynamic frictional force Fdv. The forces acting in x-axis are: the horizontal
36
static resistance force Fsh and the horizontal component of dynamic frictional force, Fdh.
The vector of gravity force G is represented as
G = ms gk0: (4.11)
The resistance force FR is calculated by Eqs. (3.6), (3.7), and (3.9). Because the resistance
force is assumed to act at the mass center C, the dynamic frictional force Fd is calculated
by Eq. (3.6) as
Fd = ?vE ?d ?g Ar jvEj = ?vC ?d ?g Ar jvCj
= ?d ?g Ar
p
_q2x + _q2z
h
? _qx?0 ? _qz k0
i
: (4.12)
The reference area Ar for the sphere is the same as that of the mathematical pendulum as
Ar = ?d
2s
4 : (4.13)
The horizontal and vertical static resistance forces, Fsh and Fvh, are
Fsh = ?sign(_qx)?h g?g q2z ds ?0; (4.14)
Fsv = ?sign(_qz)?v (qz=ds)? g?g V k0; (4.15)
where the immersed volume V is calculated as the volume of the sphere.
V = ?d
3s
6 : (4.16)
37
Table 4.2: Stopping time of the rigid sphere into a granular medium
Impact angle (o) vp(0) (m/s) ts (s)
1 0.0565723
15 3 0.04457
5 0.0376819
1 0.0364601
45 3 0.0272604
5 0.0243122
1 0.0324081
75 3 0.0251294
5 0.0227529
The resistance force FR, the sum of the dynamic frictional force vector Fd and the static
resistance force vector Fs, is represented by the sum of Eqs. (4.12), (4.14), and (4.15) as
FR = Fd +Fsh +Fsv
=
"
? _qx ?d ?g ?d
2s
4
p
_q2x + _q2z ?sign(_qx)?h g?g q2z ds
#
?0
"
? _qz ?d ?g ?d
2s
4
p
_q2x + _q2z ?sign(_qz)?v (qz=ds)? g?g ?ds
3
6
#
k0: (4.17)
Simulation results
Figures 4.6, 4.7, and 4.8 represent the simulation results of the oblique impact of the
sphere, depicted in Figs. 2.3 and 4.5. All data for the second simulation such as the density
and the diameter of the sphere, the density of the granular medium, and the resistance force
coe?cients are applied as the same data applied to the simulation of the mathematical
pendulum. The simulation is performed for the difierent impact angles (q(0) = 15, 45,
and 75o) and the difierent impact linear velocities (vp(0) = 1, 3, and 5 m/s) from the
impact moment until the velocity of the sphere, vp, becomes zero. The penetrating distance
38
dp
(m)
Tim e (s)
vp
(m/s)
Tim e (s)
0.00 0.01 0.02 0.03 0.04 0.05
0
1
2
3
4
5
0.00 0.01 0.02 0.03 0.04 0.05
0.000
0.005
0.010
0.015
0.020
vp(0) = 1 m/s
vp(0) = 3 m/s
vp(0) = 5 m/s
Figure 4.6: Impact results of the rigid sphere into a granular medium for q(0) = 15o
39
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
0.000
0.005
0.010
0.015
dp
(m)
Tim e (s)
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
0
1
2
3
4
5
vp
(m/s)
Tim e (s)
vp(0) = 1 m/s
vp(0) = 3 m/ s
vp(0) = 5 m/ s
Figure 4.7: Impact results of the rigid sphere into a granular medium for q(0) = 45o
40
0.000 0.005 0.010 0.015 0.020 0.025 0.030
0.000
0.005
0.010
0.015
dp
(m)
Tim e (s)
0.000 0.005 0.010 0.015 0.020 0.025 0.030
0
1
2
3
4
5
vp
(m/s)
Tim e (s)
vp(0) = 1 m/s
vp(0) = 3 m/ s
vp(0) = 5 m/ s
Figure 4.8: Impact results of the rigid sphere into a granular medium for q(0) = 75o
41
(dp = R vp dt) increases but the penetrating motion of the rigid sphere comes to stops fast
when the initial impact velocity v(0) is increased as shown in Figs. 4.6, 4.7, and 4.8 and in
table 4.2. These simulation results show that the results conflrmed in the simulation for one
D.O.F. motion of the particle type sphere as the mathematical pendulum and the vertical
impact are kept in two D.O.F. motion of the particle type sphere.
4.2 Single impact with a granular medium
4.2.1 Impact of a free rigid link
Modeling
The model of the impact of a free rigid link into a granular matter is presented. Previous
the impact of the sphere are modeled based on the motion of particle. However, from the
impact of a rigid link the impact with a granular is modeled based on the motion of rigid
body. In the vertical z-axis, the forces acting on the link are: the gravity force G, the
vertical static resistance force Fsv, and the vertical component of the dynamic frictional
force Fd. The forces acting in x-axis are: the horizontal static resistance force Fsh and the
horizontal component of the dynamic frictional force Fd. The gravity force G acts at the
center of mass, C, of the link and the resistance force, FR, including Fs and Fd acts at the
point E, where point E is the centroid point of immersed part as shown in Fig. 4.9. The
general equation of motion for the planar kinematic chain can be written in the following
form
mc?rC ??0 =
?
Fs +Fd
?
??0; (4.18)
mc?rC ?k0 =
?
G+Fs +Fd
?
?k0; (4.19)
IC ?q??0 =
"
rCE ?
?
Fs +Fd
?#
??0; (4.20)
where ?rC is the acceleration vector of the position vector of the mass center of the rigid
bar, and ?q is the angular acceleration vector of the angular vector of pendulum with the
42
L
T
A
q
0
k0
0?
?O
Granular
medium
E
C
G
qx
qz
Fsv
Fsh
Fd
zTz
T/2
Figure 4.9: Free body diagram of a free rigid link in impact
vertical, q = q?0. The position vector rCE represents vector from the mass center C to the
resistance force application point E, mc is the mass of the link, and IC = mc L2=12 is the
mass moment of inertia of the link with regard to its center.
rC = qx?0 +qz k0; (4.21)
?q = d
2q
dt2 ?0: (4.22)
The position vector from the mass center to the resistance force application point E, rCE,
is represented as
rCE = LCE sinq?0 +LCE sinqk0; (4.23)
where LCE is the length between the mass center C and the resistance force application
point E and calculated as
LCE = L2 ? zT2cosq: (4.24)
43
The immersed depth of the end T, zT, is expressed as
zT = rC ?k+ L2 cosq = qz + L2 cosq: (4.25)
The velocity vector of the resistance force acting point E, vE, the reference area of the
penetrating bar Ar, and the moving angle qm for calculating the dynamic frictional force
are represented as
vE = drCdt + dqdt ?rCE =
?
_qx +LCE _qcosq
?
?+
?
_qz ?LCE _qsinq
?
k; (4.26)
Ar = dc zTcosq j sin(q ?qm) j; (4.27)
qm = tan?1
v
Ex
vEz
?
= tan?1
_q
x +LCE _qcosq
_qz ?LCE _qsinq
?
: (4.28)
The dynamic frictional force Fd is calculated by Eq. (3.6) as
Fd = ?d ?g dc zTcosq
flfl
flflsin
q ?tan?1
_q
x +LCE _qcosq
_qz ?LCE _qsinq
??flfl
flfl
r?
_qx +LCE _qcosq
?2
+
?
_qz ?LCE _qsinq
?2
"
?
?
_qx +LCE _qcosq
?
?0 ?
?
_qz ?LCE _qsinq
?
k0
#
: (4.29)
The immersed volume V for calculating the static resistance force is calculated as
V = ?dc
2
4
zT
cosq: (4.30)
Hence, the horizontal and vertical static resistance forces, Fsh and Fvh, are
Fsh = ?sign
?
_qx +LCE _qcosq
?
?h g?g z2T dc ?0; (4.31)
Fsv = ?sign
?
_qz ?LCE _qsinq
?
?v
z
T
dc
??
g?g ?dc
2
4
zT
cosq k0: (4.32)
44
The resistance force FR, the sum of the dynamic frictional force vector Fd and the static
resistance force vector Fs, is represented by the sum of Eqs. (4.29), (4.31), and (4.32) as
FR = Fd +Fsh +Fsv
=
"
??d ?g dc zTcosq
flfl
flflsin
q ?tan?1
_q
x +LCE _qcosq
_qz ?LCE _qsinq
??flfl
flfl
r?
_qx +LCE _qcosq
?2
+
?
_qz ?LCE _qsinq
?2 ?
_qx +LCE _qcosq
?
?sign
?
_qx +LCE _qcosq
?
?h g?g z2T dc
#
?0 +
"
??d ?g dc zTcosq
flfl
flflsin
q ?tan?1
_q
x +LCE _qcosq
_qz ?LCE _qsinq
??flfl
flfl
r?
_qx +LCE _qcosq
?2
+
?
_qz ?LCE _qsinq
?2 ?
_qz ?LCE _qsinq
?
?sign
?
_qz ?LCE _qsinq
?
?v
z
T
dc
??
g?g?dc
2
4
zT
cosq
#
k0: (4.33)
Simulation results
Figure 4.10 represents the simulation results of the penetrating depth of the link end
T, zT, and the vertical velocity of link end T, vTz, of the impact of the vertically dropped
cylinder type rigid link depicted in Figs. 2.5 and 4.9. The applied dimensions of the link are
the length L = 0.1524 m and the diameter dc = 0.00635 m. The density of the link, ?c, is
also 7:7?103 kg/m3 and the density of a granular medium (sand), ?g, is 2:5?103 kg/m3.
The dynamic frictional force coe?cient ?d = 6.5, the horizontal static resistance force
coe?cient ?h = 8, and the vertical static resistance force coe?cients ?v = 22 and ? = 1.1 are
used in the simulation. This simulation is performed for the difierent initial impact vertical
45
Table 4.3: Stopping time of the free link
q(0) (o) vTz(0) (m/s) tz (s)
1.53 0.0331795
0 2.06 0.029688
2.47 0.0277679
1.26 0.0275914
32 1.87 0.0229915
2.33 0.0206495
1.45 0.017818
55 1.98 0.0150461
2.43 0.0133007
velocities (_qz(0) = 1.53, 2.06, and 2.47 m/s) from the impact moment until the link stop.
As shown in Fig. 4.10, the penetrating depth of the link increases but the time interval of
the link decreases as initial impacting velocity increases. The resistance force acting at the
rigid link is shown in Figs. 4.11 and 4.12. The dynamic frictional force depending on the
velocity acts as governing resistance force at flrst. However, when the penetrating depth
increases and the velocity of the link decreases, the static resistance force depending on the
immersed depth acts as governing resistance force.
Figures 4.13 and 4.14 represent the simulation results performed for the difierent impact
angle q(0) and the difierent initial impact vertical velocities (_qz(0) = 1.26, 1.87, and 2.33 m/s
for q(0) = 32o, and _qz(0) = 1.45, 1.98, and 2.43 m/s for q(0) = 55o). For these simulations,
we observed a stopping time tz representing the time interval from the moment of impact
until the vertical velocity of the link end, vTz. Even though the link impacts obliquely,
the penetrating depth of the link end T, zT, increases but the vertical velocity of link end
T, vTz, becomes zero fast when the initial vertical impact velocity increases as shown in
Figs. 4.10, 4.13, and 4.14 and in table 4.3. This result shows that the relations among the
46
vT
z
(m/s)
zT
(m)
Time (s)
Time (s)
vTz (0) = 1.53 m/s
vTz (0) = 2.06 m/s
(0) = 2.47 m/svTz
0.000 0.005 0.010 0.015 0.020 0.025 0.030
0.0
0.5
1.0
1.5
2.0
2.5
0.000 0.005 0.010 0.015 0.020 0.025 0.030
0.00
0.01
0.02
0.03
Figure 4.10: Impact results of the rigid cylinder type link for q(0) = 0o
47
Tim e (s)
F
d
(N)
Tim e (s)
F
s
(N)
0.000 0.005 0.010 0.015 0.020 0.025 0.030
-3.0
-2.5
-2.0
-1.5
-1.0
-0.5
0.0
0.000 0.005 0.010 0.015 0.020 0.025 0.030
-4
-3
-2
-1
0
vTz (0) = 1.53 m/s
vTz (0) = 2.06 m/s
(0) = 2.47 m/svTz
Figure 4.11: Resistance force Fd and Fs for q(0) = 0o
48
Tim e (s)
F
R
(N
)
vTz (0) = 1.53 m/ s
vTz (0) = 2.0 6 m/ s
(0) = 2.47 m/ svTz
0.000 0.005 0.010 0.015 0.020 0.025 0.030
-4.5
-4.0
-3.5
-3.0
-2.5
-2.0
-1.5
Figure 4.12: Total resistance force FR for q(0) = 0o
stopping time, the penetrating depth, and the initial impact velocity are kept in even multi
D.O.F penetrating motion.
49
vT
z
(m/s)
zT
(m)
Time (s)
Time (s)
vTz (0) = 1.26 m/s
vTz (0) = 1.87 m/s
(0) = 2.33 m/svTz
0.
0.
1.
1.
2.
0.000 0.005 0.010 0.015 0.020 0.025
0
5
0
5
0
0.000 0.005 0.010 0.015 0.020 0.025
0.000
0.005
0.010
0.015
0.020
0.025
Figure 4.13: Impact results of the rigid cylinder type link for q(0) = 32o
50
vT
z
(m/s)
zT
(m)
Time (s)
Time (s)
vTz (0) = 1.45 m/s
vTz (0) = 1.98 m/s
(0) = 2.43 m/svTz
0.000 0.005 0.010 0.015
0.0
0.5
1.0
1.5
2.0
0.000 0.005 0.010 0.015
0.000
0.005
0.010
0.015
Figure 4.14: Impact results of the rigid cylinder type link for q(0) = 55o
51
T
Ai?1 Ai
Ci
Gi
Granular medium EFRx
FRz
nC
nG
-1An
F xn?1,n
F xi?1,i
F xi +1,i
F zi +1,i
F zn?1,n
F zi?1,i
FR
O x
zy
00
0
k
?
?
Figure 4.15: Free body diagram of a single impact of a planar kinematic chain
4.2.2 Impact of a planar kinematic chain
Modeling
As mentioned in chapter 2, the angles denoted by qis between the z-axis and the link
i are the generalized coordinates. The external forces acting on the kinematic chain is the
gravity force and the resistance forces as shown in Fig. 4.15. The gravity force Gi acts at
the center of mass, Ci, of the link i and the resistance force, FR including Fs and Fd, acts
at the point E of the last link n, where the point E is the centroid of the immersed part as
shown in Fig. 4.15. The general equation of motion for the planar kinematic chain can be
52
written in the following form
M(q)?q+C(_q;q)+G(q) = T+D(q)F; (4.34)
where q = (q1;:::;qn)T is the n?1 generalized coordinates, M(q) is the n?n mass matrix,
C(_q;q) is an n?1 vector, G(q) is the n?1 gravity vector, T is the n?1 joint moment
vector, D(q) is an n?2 matrix, and F = (FRx;FRz)T is the resistance force vector.
Under the assumption that the joint moments do not exist, the Newton-Euler?s equa-
tions is applied to formulate the difierential equation of motion. For a compound pendulum
(n = 1), the equation of motion can be written as
I ?q =
?
rC ?G+rE ?FR
?
??0; (4.35)
where I is the mass moment of inertia with regard to the joint. The vector rC representing
the position vector from the joint to the mass center C is given as
rC = LC
?
sinq?0 +cosqk0
?
; (4.36)
where LC is the length from the joint to the mass center of the pendulum. The position
vector from the joint to the resistance force application point E, rE, and the resistance
force FR will be formulated in the case of the n link multi kinematic chain and can be
particularized for n = 1.
In the case of a multi link kinematic chain, a force and a moment equation for each
link i = 1;2;:::;n?1 can be written
miaCi = Fi?1;i +Fi+1;i +Gi; (4.37)
Ii ?qi =
?
rCiAi?1 ?Fi?1;i +rCiAi ?Fi+1;i
?
??0; (4.38)
53
and for the last link i = n
mnaCn = Fn?1;n +Gn +FR; (4.39)
In ?qn =
?
rCnAn?1 ?Fn?1;n +rCnE ?FR
?
??0; (4.40)
where mi is the mass of the link i, Ii is the mass moment of inertia of the link i with regard
to its center, and Li is the length of the link i. The force Fi?1;i is the reaction force of the
link i?1 on the link i at Ai?1, Gi is the gravity force vector acting on the mass center Ci
of the link i, and FR is the resistance force vector acting on the last link n. The vector aCi
is the acceleration vector of the position vector rCi and rCiAi?1 is the position vector from
the mass center of the link i to the point Ai?1 of the link i. The gravity force vector Gi,
the acceleration vector aCi, and the position vectors rCi, rCiAi?1 are represented as
Gi = mi gk0; (4.41)
rCi =
0
@
i?1X
j=1
Lj sinqj +LCi sinqi
1
A?0 +
0
@
i?1X
j=1
Lj cosqj +LCi cosqi
1
Ak0; (4.42)
aCi = d
2rC
i
dt2 =
2
4
i?1X
j=1
Lj ??qj cosqj ? _q2j sinqj?+LCi ??qi cosqi ? _q2i sinqi?
3
5?0
?
2
4
i?1X
j=1
Lj ??qj sinqj + _q2j cosqj?+LCi ??qi sinqi + _q2i cosqi?
3
5k0; (4.43)
rCiAi?1 = ?LCi sinqi?0 ?LCi cosqik0; (4.44)
rCiAi = (Li ?LCi)sinqi?0 +(Li ?LCi)cosqik0; (4.45)
where LCi is the length from the joint Ai?1 to the mass center Ci of the link i. The vector
rCnE representing the position vector from the mass center of the last link n, Cn, to the
resistance force acting point E is calculated from the vector rE representing the position
54
vector from the origin to the force acting point E as
rE =
0
@
n?1X
j=1
Lj sinqj +LE sinqn
1
A?0 +
0
@
n?1X
j=1
Lj cosqj +LE cosqn
1
Ak0; (4.46)
rCnE = rE ?rCn =
?
LE ?LCn
??
sinqn?0 +cosqnk0
?
; (4.47)
where LE, the length from the last joint An?1 to the resistance force acting point E, is
represented as
LE = Ln ? 12 zTcosq
n
: (4.48)
The immersed depth of the last link n, zT, is calculated as
zT =
nX
i=1
Li cosqi ?
nX
i=1
Li cosqi(0); (4.49)
where qi(0) is the initial impact angle of qi.
The other variable of the resistance force vector FR(vE;zT), the velocity vector vE, is
calculated as
vE = vCn +vCnE = drCndt +qn?0 ?rCnE
=
0
@
n?1X
j=1
Lj _qj cosqj +LE _qn cosqn
1
A?0 ?
0
@
n?1X
j=1
Lj _qj sinqj +LE _qn sinqn
1
Ak0: (4.50)
The reference area Ar and the moving angle qm of the penetrating last link n for calculating
the dynamic frictional force are represented as
Ar = dc zTcosq
n
j sin(qn ?qm) j; (4.51)
qm = tan?1
?
?
Pn?1
j=1 Lj _qj cosqj +LE _qn cosqnP
n?1
j=1 Lj _qj sinqj +LE _qn sinqn
!
: (4.52)
55
Therefore the dynamic frictional force Fd is given by Eq. (3.6) as
Fd = ?d ?g dc zTcosq
n
flfl
flfl
flsin
?
qn ?tan?1
?
?
Pn?1
j=1 Lj _qj cosqj +LE _qn cosqnP
n?1
j=1 Lj _qj sinqj +LE _qn sinqn
!!flfl
flfl
fl
vu
uu
t
0
@
n?1X
j=1
Lj _qj cosqj +LE _qn cosqn
1
A
2
+
0
@
n?1X
j=1
Lj _qj sinqj +LE _qn sinqn
1
A
2
"
?
0
@
n?1X
j=1
Lj _qj cosqj +LE _qn cosqn
1
A?0
+
0
@
n?1X
j=1
Lj _qj sinqj +LE _qn sinqn
1
Ak0
#
: (4.53)
The immersed volume V for the static resistance force is calculated as
V = ?dc
2
4
zT
cosq: (4.54)
Hence, the horizontal and vertical static resistance forces, Fsh and Fvh are
Fsh = ?sign
0
@
n?1X
j=1
Lj _qj cosqj +LE _qn cosqn
1
A ?h g?g z2T dc ?0; (4.55)
Fsv = sign
0
@
n?1X
j=1
Lj _qj sinqj +LE _qn sinqn
1
A ?v
z
T
dc
??
g?g ?dc
2
4
zT
cosq k0; (4.56)
The resistance force FR, the sum of the dynamic frictional force vector Fd and the static
resistance force vector Fs, is represented by the sum of Eqs. (4.53), (4.55), and (4.56).
The joint reaction force Fi?1;i can be calculated using the relations Fi+1;i = ?Fi;i+1
and Eqs. (4.37) and (4.39) as
Fi?1;i =
nX
j=i
mj aCj ?Gj ?FR: (4.57)
56
The equations of motion described by Eqs. (4.37), (4.38), (4.39), and (4.40) can be simplifled
and rewritten for each link i = 1;2;:::;n?1
Ii ?qi = ?
"
LCi
?
sinqi?0 +cosqik0
?
?
0
@
nX
j=i
mj aCj ?Gj ?FR
1
A+
?
Li ?LCi
??
sinqi?0 +cosqik0
?
?
0
@
nX
j=i+1
mj aCj ?Gj ?FR
1
A
#
??0;(4.58)
and for the last link i = n
In ?qn =
"
?LCn
?
sinqn?0 +cosqnk0
?
?
?
mnaCn ?Gn ?FR
?
+
?
LE ?LCn
??
sinqn?0 +cosqnk0
?
?FR
#
??0: (4.59)
The flnal nonlinear equations of motion will have the form as
fi
?
?qi; _qi;mi;Li;g;FR
?
= 0: (4.60)
Simulation results of a compound pendulum
The dimensions of the rigid compound pendulum are shown in Fig. 4.16. The dimen-
sions of the pendulum are the length Lp1 = 0:78" (= 0.019812 m), Lp2 = 0:6" (= 0.01524
m), Lp3 = 9:4" (= 0.23876 m) and the diameter dp1 = 0:75" (= 0.01905 m), dp2 = 0:6" (=
0.01524 m), dp3 = dc = 0:25" (= 0.00635m). The density of the link, ?c, is also 7:7?103
kg/m3 and the density of a granular medium (sand), ?g, is applied as 2:5?103 kg/m3. The
dynamic frictional force coe?cient ?d = 6.5, the horizontal static resistance force coe?cient
?h = 8, and the vertical static resistance force coe?cients ?v = 22 and ? = 1.1 are used
in the simulation. Figures 4.17, 4.18, 4.19, and 4.20 represent the simulation results of the
impact of the rigid compound pendulum.
The simulations are performed for difierent impact angle q(0) and difierent initial im-
pact velocities (_q(0) = -1.75, -3.38, and -4.66 rad/s for q(0) = 22o, _q(0) = -3.31, -6.24, and
57
O
Granular
medium
C
E
G
zTzT/2
Fsv
shF
CL
EL
Fd
L
?q
q
x
z
y 0k
0?
0?
Lp1
Lp2
Lp3
dp1
dp2
dp3
Figure 4.16: Compound pendulum
-8.41 rad/s for q(0) = 31o, _q(0) = -2.66, -6.47, and -9.06 rad/s for q(0) = 45o, and _q(0) =
-2.70, -6.54, and -9.17 rad/s for q(0) is 61.5o) from the impact moment until the angular
velocity of the pendulum, _q, becomes zero. As shown in Figs. 4.17, 4.18, 4.19, and 4.20 and
in table 4.4, the penetrating angle q of the pendulum increases but the time the angular
velocity _q becomes zero decreases when the initial vertical impact velocity increases as the
results of previous models.
58
?q
(rad/
s)
Tim e (s)
Tim e (s)
q
(o
)
?q (0) = ? 1.75 rad/s
?q (0) = ? 3.38 rad/s
?q (0) = ? 4.66 rad/s
0.00 0.02 0.04 0.06 0.08 0.10
14
16
18
20
22
0.00 0.02 0.04 0.06 0.08 0.10
-4
-3
-2
-1
0
Figure 4.17: Impact results of the rigid compound pendulum for q(0) = 22o
59
?q
(rad/
s)
Tim e (s)
Tim e (s)
q
(o
)
?q (0) = ? 3.31 rad/s
?q (0) = ? 6.24 rad/s
?q (0) = ? 8.41 rad/s
0.00 0.01 0.02 0.03 0.04 0.05 0.06
24
26
28
30
0.00 0.01 0.02 0.03 0.04 0.05 0.06
-8
-6
-4
-2
0
Figure 4.18: Impact results of the rigid compound pendulum for q(0) = 31o
60
?q
(rad/
s)
Tim e (s)
Tim e (s)
q
(o
)
?q (0) = ? 2.66 rad/ s
?q (0) = ? 6.47 rad/ s
?q (0) = ? 9.06 rad/ s
0.00 0.01 0.02 0.03 0.04
-8
-6
-4
-2
0
0.00 0.01 0.02 0.03 0.04
39
40
41
42
43
44
45
Figure 4.19: Impact results of the rigid compound pendulum for q(0) = 45o
61
?q
(rad/
s)
Tim e (s)
Tim e (s)
q
(o
)
?q (0) = ? 2.70 rad/s
?q (0) = ? 6.54 rad/s
?q (0) = ? 9.17 rad/s
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
57
58
59
60
61
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
-8
-6
-4
-2
0
Figure 4.20: Impact results of the rigid compound pendulum for q(0) = 61:5o
62
Table 4.4: Stopping time of the rigid compound pendulum
q(0) (o) _q(0) (rad/s) ts (s)
-1.75 0.107212
22 -3.38 0.0861789
-4.66 0.0775028
-3.31 0.0609829
31 -6.24 0.0485081
-8.41 0.043819
-2.66 0.046688
45 -6.47 0.0333638
-9.06 0.0295578
-2.70 0.0343956
61.5 -6.54 0.0245294
-9.17 0.0216994
O
1
2
Granular medium
zy
q1
q2
q1?
q2?
x
00
0
k
?
?
zTFsh
FsvFd
zT/2
E
T
G2
G1
Figure 4.21: Double pendulum
63
Simulation results of a double pendulum
The double pendulum under consideration is shown in Fig. 4.21. As the numerical data
applied to the simulation, the lengthes of link 1 and 2, L1 = 0:2 m and L2 = 0:6 m, the
diameter dc = 0:0254 m, the density of the link ?c = 7:7?103 kg/m3, and the density of a
granular medium (sand) ?g = 2:5?103 kg/m3, the dynamic frictional force coe?cient ?d =
6.5, the horizontal static resistance force coe?cient ?h = 8, and the vertical static resistance
force coe?cients ?v = 22 and ? = 1.1 are used in the simulation.
Figures 4.22 and 4.23 represent the simulation results for the double pendulum with
the following initial conditions: q1(0) = 30o, q2(0) = 75o (q1(0) < q2(0)) and _q1(0) = _q2(0)
= -1, -3, -5 rad/s. The initial vertical tip velocity vTz(0) of the impact point and horizontal
one vTx(0) are increasing. The total stopping time tt is deflned as the time starting from
the impact moment until _q1 and _q2 become zero simultaneously (_q1(tt) = _q2(tt) = 0). For
this case the total stopping time tt is decreasing with the increasing of the initial velocity
as shown in Fig. 4.22 and table 4.5. The vertical tip velocity vTz stops flrst and the total
stopping time tt is dictated by the horizontal tip velocities vTx as shown in Fig. 4.23. For
this case the tip velocities vTz and vTx do not change the sign.
For this simulation, we also observed an additional stopping time tz representing the
timefromthemomentofimpactuntiltheverticalvelocityofthependulumtip, vTz, becomes
Table 4.5: Stopping time and penetrating depth results of the double pendulum for q1(0)
= 30o and q2(0) = 75o
_q1(0) (rad/s) _q2(0) (rad/s) tt (s) tz (s) zT (m)
-1 -1 0.440289 0.0546257 0.0277009
-3 -3 0.309807 0.0355301 0.0347221
-5 -5 0.260902 0.0289713 0.0389723
-1 0.319782 0.0450784 0.0288578
-5 -3 0.286156 0.0338793 0.0347693
-5 0.260902 0.0289713 0.0389723
64
?q 1
(rad/
s)
Tim e (s)
?q 2
(rad/
s)
Tim e (s)
0.0 0.1 0.2 0.3 0.4
-5
-4
-3
-2
-1
0
0.0 0.1 0.2 0.3 0.4
-5
-4
-3
-2
-1
0
?q1(0) = ?q2(0) = ? 1 rad/ s
?q1(0) = ?q2(0) = ? 3 rad/ s
?q1(0) = ?q2(0) = ? 5 rad/ s
Figure 4.22: Angular velocities _q1 and _q2 results of the double pendulum impact for q1(0) =
30o, q2(0) = 75o, and _q1(0) = _q2(0) = -1, -3, -5 rad/s
65
vT
z
(m/s)
Tim e (s)
vT
x
(m/s)
Tim e (s)
?q1(0) = ?q2(0) = ?1 rad /s
?q1(0) = ?q2(0) = ?3 rad /s
?q1(0) = ?q2(0) = ?5 rad /s
0.0 0.1 0.2 0.3 0.4
0.0
0.5
1.0
1.5
2.0
2.5
3.0
0.0 0.1 0.2 0.3 0.4
-1.5
-1.0
-0.5
0.0
Figure 4.23: Velocities vTx and vTz results of the double pendulum impact for q1(0) = 30o,
q2(0) = 75o, and _q1(0) = _q2(0) = -1, -3, -5 rad/s
66
?q 2
(rad/
s)
Tim e (s)
?q1(0) = ? 5, ?q2(0) = ? 1 rad/s
?q1(0) = ? 5, ?q2(0) = ? 3 rad/ s
?q1(0) = ? 5, ?q2(0) = ? 5 rad/ s
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
-4
-2
0
2
4
?q 1
(rad/
s)
Tim e (s)
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
-10
-8
-6
-4
-2
0
Figure 4.24: Angular velocities _q1 and _q2 results of the double pendulum impact for q1(0) =
75o, q2(0) = 30o, _q1(0) = -5 rad/s, and _q2(0) = -1, -3, -5 rad/s
67
v T
z
(m/s)
Tim e (s)
v T
x
(m/s)
Tim e (s)
?q1(0) = ? 5, ?q2(0) = ? 1 rad/s
?q1(0) = ? 5, ?q2(0) = ? 3 rad/s
?q1(0) = ? 5, ?q2(0) = ? 5 rad/s
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
-2.5
-2.0
-1.5
-1.0
-0.5
0.0
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
-0.5
0.0
0.5
1.0
1.5
2.0
2.5
Figure 4.25: Velocities vTx and vTz results of the double pendulum impact for q1(0) = 75o,
q2(0) = 30o, _q1(0) = -5 rad/s, and _q2(0) = -1, -3, -5 rad/s
68
zero. For this interval of time we represent the penetrating depth of the pendulum tip, zT.
As shown in table 4.5, the stopping time tz is decreasing and the penetrating depth of the
pendulum tip, zT is increasing with the increasing of the initial velocity. Using difierent
initial impact velocities, (_q1(0) = -5 rad/s, _q2(0) = -1, -3, -5 rad/s), the simulation results
are presented in table 4.5. For these cases we can draw the same remarks about tt, tz, and
zT.
Using a difierent initial geometry for the links (q1(0) = 75o and q2(0) = 30o (q1(0) >
q2(0))) and and initial impact velocity conflgurations (_q1(0) = ?_q2(0) = 1, 3, 5 rad/s), the
simulation results for tt, tz, and zT are shown in table 4.6. For these simulations the total
stopping time tt and the stopping time for z direction are decreasing and the penetrating
depth of the pendulum tip, zT, is increasing with the initial impact velocities as previous
initial geometric conditions.
A difierent phenomenon is observed for the same conflgurations (q1(0) = 75o and
q2(0) = 30o) but difierent initial velocities (_q1(0) = -5 rad/s and _q2(0) = -1, -3, -5 rad/s).
The stopping times are almost the same at the end of the simulation time as shown in
Figs. 4.24 and 4.25. The vertical tip velocity vTz stops flrst and the total stopping time
tt is dictated by the horizontal tip velocities vTx as shown in Fig. 4.25. For this case the
vertical tip velocity vTz and the horizontal tip velocity vTx change the sign during the
impact process. The stopping time behavior is also difierent from the previous case because
for _q2(0) = -3 rad/s the total stopping time tt is increasing even though tt is again decreasing
for _q2(0) = -5 rad/s as shown in table 4.6. Furthermore the stopping time tz is increasing
with the initial velocities for all the cases. The withdrawing motion (vTz < 0) of the link 2
happening during the impact process might afiect the stopping time. From these simulation
results for the initial impact conditions (q1(0) > q2(0) and _q1(0) > _q2(0)), the stopping time
of the double pendulum can be afiected by the withdrawing motion of the link 2. The
stopping times increase or decrease depending on the efiects of the withdrawing motion.
For this case the difierences between stopping time are observed to be small. However, the
69
Table 4.6: Stopping time and penetrating depth results of the double pendulum for q1(0)
= 75o and q2(0) = 30o
_q1(0) (rad/s) _q2(0) (rad/s) tt (s) tz (s) zT (m)
-1 -1 0.470185 0.133222 0.0674095
-3 -3 0.416271 0.112135 0.0826002
-5 -5 0.378573 0.0865091 0.0898785
-1 0.370393 0.0795579 0.067502
-5 -3 0.383146 0.0833189 0.0810202
-5 0.378573 0.0865091 0.0898785
penetrating depth of the pendulum tip zT is increasing with the initial velocities the same
as the previous simulation results.
4.3 Multiple impacts of a planar kinematic chain
Modeling
The angle between the z-axis and the link i is denoted by qi, (i = 1;2;3;:::;n). The
x and z-axis distances from the origin to a point (joint A1 in this model) are denoted by qx
and qz respectively as shown in Fig. 2.7. The gravity force Gi acts at the mass center Ci
of the link i. The resistance force FRj acts at the point Ej of the link j interacting with
the granular medium, where the point Ej is the centroid point of the immersed part of the
link j. The equations of motion for the planar kinematic chain with multiple contact points
shown in Fig. 4.26 can be written in the following general form
M(q)?q+C(q; _q)+G(q) =
2
66
66
4
T
0
0
3
77
77
5
+
2
64 D1(q)
D2
3
75F; (4.61)
70
Granular
medium
A1
A2
2
F x1,2
F x3,2
F z1,2
F z3,2C
2
G2
F x2,1
F z2,1
A1
C1
E1
T1
FR?x
FR?z
G1
1
FR?
Granular
medium
T2
i
F xi?1,i
F xi?1,i
E2F
R?x
FR?z
Ai
Gi
FR?
Granular
medium
Tm
Em
Al
n
nC
nG
F xn?1,n
F zn?1,n
FRm
FRmx
FRmz
O x
zy
00
0
k
?
?
Figure 4.26: Free body diagram of multiple impacts of a planar kinematic chain
E1
E2
C1 C2
G2
A
T1
T2
q11
q2
2
Granular
medium
zT1
qz
qx
zT2
O x
zy
FR?x
FR?z
G1
FR?
FR?x
FR?zFR?
00
0
k
?
?
Figure 4.27: Two link kinematic chain with two impact points model
71
where q = (q1;:::;qn;qx;qz)T is the (n + 2) ? 1-dimensional vector of generalized co-
ordinates; _q = (_q1;:::; _qn; _qx; _qz)T is the (n + 2) ? 1-dimensional vector of generalized
velocities, ?q = (?q1;:::; ?qn; ?qx; ?qz)T is the (n + 2) ? 1-dimensional vector of generalized
accelerations, M(q) is the (n + 2) ? (n + 2) mass matrix, C(q; _q) is a (n + 2) ? 1 vec-
tor, G(q) is the (n + 2) ? 1 vector of gravity terms, T is the n?1 vector of joint mo-
ments, D1(q) is an n ? 2k matrix, D2 = (1;0;1;0;:::;0;1;0;1;:::) is the 2 ? 2k matrix,
F = (FR1x;FR1z;FR2x;FR2z;:::;FRkx;FRkz)T is the 2k ?1 vector of resistance forces, and
the components FRjx and FRjz are the x-axis and the z-axis components of the resistance
force FRj respectively. Applying Newton-Euler?s equations for the kinematic chain, the
matrices M; C; G; T; D1; D2, and F can be calculated for a speciflc example.
In this study, an ideal two link chain are considered as the simplifled application of the
kinematic chain model, as shown in Fig. 4.27. The number of D.O.F. is 4 and its generalized
coordinates areq = (q1;q2;qx;qz)T. Under the same assumption that the joint moment does
not exist, the Newton-Euler?s equations give the governing difierential equations for the two
link chain:
m1aC1 = G1 +FR1 +F2;1; (4.62)
m2aC2 = G2 +FR2 +F1;2; (4.63)
I1 ?q1? = rC1A ?F2;1 +rC1E1 ?FR1
= rC1A ?(m1aC1 ?G1 ?FR1)+rC1E1 ?FR1; (4.64)
I2 ?q2? = rC2A ?F1;2 +rC2E2 ?FR2
= rC2A ?(m2aC2 ?G2 ?FR2)+rC2E2 ?FR2; (4.65)
where mi is the mass of the link i, Ii = mi L2i=12 is the mass moment of inertia of the link
i with regard to its mass center Ci, and Li is the length of the link i. The force F1;2 is the
reaction force of the link 2 on the link 1 at the joint A and the force F2;1 is vice versa. Gi is
the gravity force vector acting on the mass center Ci of the link i, and FRi is the resistance
72
force vector acting on the last link n. Equations (4.62) and (4.63) can be simplifled when
these equations are summed due to the characteristics of the reaction forces as
m1aC1 +m2aC2 = G1 +G2 +FR1 +FR2: (4.66)
Equation (4.66) can be divided into x and z components as
m1 aC1x +m2 aC2x = FR1x +FR2x; (4.67)
m1 aC1z +m2 aC2z = G1 +G2 +FR1z +FR2z: (4.68)
Finally the motion of the chain can be described by Eqs. (4.64), (4.65), (4.67), and (4.68).
The vector aCi is the acceleration vector of the position vector rCi. The position vector
rCi is represented as the sum of the vector rOA and the vector rACi. The vector rOA is
the position vector from the origin O to the joint point A. The vector rCiA is the position
vector from the mass center of the link i, Ci, to the joint point A and rCiEi is the position
vector from the mass center of the link i, Ci, to the resistance force acting point of the
link i, Ei. The position vector rOA, rCiA, rCiEi, rCi, and the acceleration vector aCi, are
represented as
rOA = qx?0 +qz k0; (4.69)
rCiA = ?Li2 sinqi?0 ? Li2 cosqik0; (4.70)
rACi = ?rCiA = Li2 sinqi?0 + Li2 cosqik0; (4.71)
rCiEi = rAEi ?rACi =
LEi ? Li2
?
(sinqi?0 +cosqik0); (4.72)
rCi = rOA +rACi =
qx + Li2 sinqi
?
?0 +
qz + Li2 cosqi
?
k0; (4.73)
aCi =
?
?qx + Li2
?
?qi cosqi ? _q2i sinqi
??
?0 +
?
?qz ? Li2
?
?qi sinqi + _q2i cosqi
??
k0; (4.74)
73
where LEi, the length from the joint A to each resistance force acting point Ei, is calculated
as
LEi = Li ? 12 zTicosq
i
: (4.75)
The immersed depth of each link i, zTi, is calculated as
zTi = Li cosqi +qz: (4.76)
Each static resistance force is depending on the immersed depth zTi. The velocity vector
of the force acting point Ei, vEi, for deciding the resistance forces is calculated from the
position vector of each force acting point Ei, rEi as
rEi = rOA +rAEi = (qx +LEi sinqi)?+(qz +LEi cosqi)k; (4.77)
vEi = drOAdt +qi??rEi = (_qx +LEi _qi cosqi)?+(_qz ?LEi _qi sinqi)k: (4.78)
The resistance force FR is also calculated by Eqs. (3.6), (3.8), and (3.9). The reference
area Ari and the moving angle qmi of the link i are calculated as
Ari = dc zTicosq
i
j sin(qi ?qmi) j; (4.79)
qmi = tan?1
_q
x +LEi _qi cosqi
_qz ?LEi _qi sinqi
?
: (4.80)
Therefore the dynamic frictional force acting at the link i, Fdi, is
Fdi = ?d ?g dc zTicosq
i
flfl
flflsin
qi ?tan?1
_q
x +LEi _qi cosqi
_qz ?LEi _qi sinqi
??flfl
flfl
r?
_qx +LEi _qi cosqi
?2
+
?
_qz ?LEi _qi sinqi
?2
"
?
?
_qx +LEi _qi cosqi
?
?0 ?
?
_qz ?LEi _qi sinqi
?
k0
#
: (4.81)
74
The immersed volume of the link i, V, is calculated as
Vi = ?dc
2
4
zTi
cosqi: (4.82)
The horizontal and vertical static resistance force of the link i, Fsh and Fvh, are
Fshi = ?sign
?
_qx +LEi _qi cosqi
?
?h g?g
?
Li cosqi +qz
?2
dc?0; (4.83)
Fsvi = ?sign
?
_qz ?LEi _qi sinqi
?
?v
L
i cosqi +qz
dc
??
g?g ?dc
2
4
zTi
cosqi k0: (4.84)
The resistance force FR, the sum of the dynamic frictional force vector Fd and the static
resistance force vector Fs, is represented as the sum of Eqs. (4.81), (4.83), and (4.84).
FRi = Fdi +Fshi +Fsvi
=
"
??d ?g dc zTicosq
i
flfl
flflsin
qi ?tan?1
_q
x +LEi _qi cosqi
_qz ?LEi _qi sinqi
??flfl
flfl
r?
_qx +LEi _qi cosqi
?2
+
?
_qz ?LEi _qi sinqi
?2 ?
_qx +LEi _qi cosqi
?
?sign
?
_qx +LEi _qi cosqi
?
?h g?g
?
Li cosqi +qz
?2
dc
#
?0 +
"
??d ?g dc zTicosq
i
flfl
flflsin
qi ?tan?1
_q
x +LEi _qi cosqi
_qz ?LEi _qi sinqi
??flfl
flfl
r?
_qx +LEi _qi cosqi
?2
+
?
_qz ?LEi _qi sinqi
?2 ?
_qz ?LEi _qi sinqi
?
?sign
?
_qz ?LEi _qi sinqi
?
?v
L
i cosqi +qz
dc
??
g?g ?dc
2
4
zTi
cosqi
#
k0: (4.85)
75
Simulation results
Figure 4.28 represents the simulation results of the vertical velocities of the ends T1
and T2, vT1z and vT2z, of the two link chain shown in Fig. 4.27 with the following initial
conditions: q1(0) = 330o, q2(0) = 45o, _q1(0) = _q2(0) = 0 rad/s, vAx(0) = 0 m/s, and vAz(0)
= 1, 3, 5 m/s. As the dimensions applied to this simulation, the length of link 1 L1 = 0:3 m,
the length of link 2 L2 = L1 cosq1(0)=cosq2(0) m, and the diameter dc = 0:0254 m are
utilized. The density of the link ?c = 7:7?103 kg/m3, the density of a granular medium
(sand) ?g = 2:5?103 kg/m3, the dynamic frictional force coe?cient ?d = 6.5, the horizontal
static resistance force coe?cient ?h = 8, and the vertical static resistance force coe?cients
?v = 22 and ? = 1.1 are applied in the simulation. For these initial conditions, the initial
vertical velocities of the impact points T1 and T2, vT1(0) and vT2(0), have the only vertical
velocity components. The stopping times for the two ends, tzT1 and tzT2, are deflned as the
time interval until the vertical velocities of the tips T1 and T2, vT1z and vT2z, respectively
become zero. For this case the stopping time tzT1 is decreasing with the increasing of the
initial velocity and the stopping time tzT2 is also decreasing. It is observed that the vertical
tip velocity of the long link, vT2z, stops much rapidly compared with the vertical tip velocity
of the short link.
Using difierent combinations of the initial impact velocities for the links (_q1(0) = ?_q2(0)
= 1, 3, 5 rad/s and vAx(0) = vAz(0) = 0 m/s) and the same initial geometry (q1(0) = 330o,
q2(0) = 45o), the simulation results for the vertical direction velocities of the both ends of
links are shown in Fig. 4.29. For these simulations the stopping times tzT1 and tzT2 are
decreasing with the initial impact angular velocities. It is observed that the vertical tip
velocity of the long link, vT2z, also stops fast compared with the vertical tip velocity of the
short link. The same results have been noticed using difierent initial condition (symmetric
cases: ?q1(0) = q2(0) = 15o;45o, and 75o, asymmetric cases: q1(0) = 330o, q2(0) = 60o and
75o) for symmetric (L1 = L2) and asymmetric cases (L1 6= L2).
Figures 4.30 and 4.31 show the simulation results in case when one end of a link is
resting initially and the other link impacts the granular medium with some initial angular
76
vT
2z
(m/s)
Tim e (s)
vT
1z
(m/s)
Tim e (s)
(0) =v 1 m/sAz
(0) =v 3 m/sAz
(0) =v 5 m/sAz
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07
0
1
2
3
4
5
0.00 0.02 0.04 0.06 0.08
0
1
2
3
4
5
Figure 4.28: Velocities vT1z and vT2z results of the two link kinematic chain for q1(0) = 330o,
q2(0) = 45o, and vAz(0) = 1, 3, 5 m/s
77
v T
2z
(m/s)
Tim e (s)
v T
1z
(m/s)
Tim e (s)
1(0) = 1 rad/s
1(0) = 3 rad/s
?q
?q
?q1(0) =
2(0) =?
2(0) =?
?q
?q
?q2(0) =? 5 rad/s
0.00 0.02 0.04 0.06 0.08 0.10
0.0
0.2
0.4
0.6
0.8
1.0
1.2
0.00 0.02 0.04 0.06 0.08 0.10 0.12
0.0
0.2
0.4
0.6
0.8
Figure 4.29: Velocities vT1z and vT2z results of the two link kinematic chain for q1(0) = 330o,
q2(0) = 45o, and _q1(0) = ?_q2(0) = 1, 3, 5 rad/s
78
v T
2z
(m/s)
Tim e (s)
v T
1z
(m/s)
Tim e (s)
?q1(0) = 1, ?q2(0) = 0 rad/s
?q1(0) = 3, ?q2(0) = 0 rad/s
?q1(0) = 5, ?q2(0) = 0 rad/s
0.00 0.02 0.04 0.06 0.08 0.10 0.12
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.00 0.02 0.04 0.06 0.08 0.10 0.12
0.0
0.1
0.2
0.3
0.4
0.5
0.6
Figure 4.30: Velocities vT1z and vT2z results of the two link kinematic chain for q1(0) = 330o,
q2(0) = 45o, and _q1(0) = 1, 3, 5 rad/s
79
v T
2z
(m/s)
Tim e (s)
v T
1z
(m/s)
Tim e (s)
?q1(0) = 0, ?q2(0) = -1 rad/s
?q1(0) = 0, ?q2(0) = -3 rad/s
?q1(0) = 0, ?q2(0) = -5 rad/s
0.00 0.02 0.04 0.06 0.08 0.10 0.12
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.00 0.02 0.04 0.06 0.08 0.10
0.0
0.2
0.4
0.6
0.8
1.0
1.2
Figure 4.31: Velocities vT1z and vT2z results of the two link kinematic chain for q1(0) = 330o,
q2(0) = 45o, and _q2(0) = -1, -3, -5 rad/s
80
velocities. For the initial impact angles (q1(0) = 330o, q2(0) = 45o) and the initial impact
velocities (_q1(0) = 1, 3, 5 rad/s, _q2(0) = 0 rad/s, and vAx(0) = vAz(0) = 0 m/s), the
stopping time of the impacting link, tzT1, is decreasing with the increasing of the initial
angular velocity and the stopping time of the resting link, tzT2, is little afiected by the
initial impact angular velocity but decreasing when the initial impact angular velocity of
link 1 is increasing as shown in Fig. 4.30. For the initial impact geometry (q1(0) = 330o,
q2(0) = 60o), the tendency of the stopping times of the impacting link 1 and the resting
link 2 is also kept as the same as the case of the initial impact angles (q1(0) = 330o,
q2(0) = 45o). The simulation results in the case that the short link is resting initially and
the long link impacts the granular medium are represented in Fig. 4.31. For the initial
impact angles (q1(0) = 330o, q2(0) = 45o) and the initial impact velocities (_q1(0) = 0 rad/s,
_q2(0) = ?1;?3;?5 rad/s and vAx(0) = vAz(0) = 0 m/s, the stopping time of the impacting
link, tzT2, is also decreasing with the increasing of the initial angular velocity and the
stopping time of the resting link, tzT1, is little afiected by the initial impact angular velocity
but decreasing when the initial impact angular velocity of link 2 is increasing as shown in
Fig. 4.31. For the initial impact geometry (q1(0) = 330o, q2(0) = 60o), the similar results
are observed as the case of the initial impact angles (q1(0) = 330o, q2(0) = 45o).
From these simulation results, the stopping time can be concluded as that for the cases
when one link is resting and the other one is impacting the stopping time of the impact link
is decreasing when the initial impact angular velocity increasing and the difierences between
the stopping time of the resting link is not conspicuous but the stopping time decreases with
the initial impact angular velocity of the impact link.
81
Table 4.7: Stopping time results of two link chain for q1(0) = 330o and q2(0) = 45o
vAz(0)(m=s) _q1(0) (rad/s) _q2(0) (rad/s) tzT1 (s) tzT2 (s)
1 0 0 0.0788529 0.0693035
3 0 0 0.0509913 0.0447269
5 0 0 0.0416282 0.0363391
0 1 -1 0.129307 0.109857
0 3 -3 0.115935 0.0857099
0 5 -3 0.105224 0.0714128
0 1 0 0.132282 0.125619
0 3 0 0.122285 0.122267
0 5 0 0.112958 0.118263
0 0 -1 0.135093 0.110255
0 0 -3 0.130274 0.0830734
0 0 -5 0.124795 0.0662946
82
Chapter 5
Impact of flexible links
In this chapter, the impact with a granular medium of a exible link is studied in
order to conflrm the efiects of the deformation of the link during the penetrating process.
A cylinder type exible link impacting a granular medium is utilized. The impacts of a
free elastic link and of an elastic compound pendulum are modeled based on the rigid body
models mentioned in chapter 4. The approach to the external forces such as the resistance
and the gravity force is the same as that of the rigid body models. What is difierent from
the rigid body impact modeling is the elastic deformation of a elastic link (the elastic mode
shape).
Figure 5.1 shows a exible link of length L, constant exural rigidity EI, cross sectional
area Ac, and density ?c. When there are no external forces acting on the link, small exural
vibrations of the link are governed by the following equation
EI@
4x(z;t)
@z4 +?c Ac
@2x(z;t)
@t2 = 0: (5.1)
z
x
P
L
O k
?
Figure 5.1: Transverse vibrations of a exible link
83
The general solution of Eq. (5.1) can be expressed as
x(z;t) =
1X
i=1
'i(z)qi(t): (5.2)
The functions of z and t, 'i(z) and qi(t), deflned respectively as
'i = Asin ?izL +Bcos ?izL +C sinh ?izL +Dcosh ?izL ; (5.3)
and
qi = fii cospit+fli sinpit; (5.4)
where ?i; i = 1;:::;1 are the consecutive roots of the transcendental equation satisfying
the boundary conditions of the exible link.
z x P
E
L
C
A
k0
k
0?
0?
?
(0) O
G
RF
Granular medium
qx
qz
q1
(1)
dc
Figure 5.2: Free elastic link
84
5.1 Free elastic link
5.1.1 Transverse vibration of a free link
In Fig. 5.2, a schematic representation of a free elastic link is given. The shear force
and the bending moment of the both ends of the link are always zero in the case of the free
link. The boundary conditions (free ends) are
x00(0;t) = x00(L;t) = x000(0;t) = x000(L;t) = 0: (5.5)
Using the boundary conditions given by Eq. (5.5), the relation between the constants A, B,
C, and D can be written as
2
66
66
66
66
66
66
66
64
0 ??
2i
L2 0
?2i
L2
??
2i sin?i
L2 ?
?2i cos?i
L2
?2i sinh?i
L2
?2i cosh?i
L2
??
3i
L3 0
?3i
L3 0
??
3i cos?i
L3
?3i sin?i
L3
?3i cosh?i
L3
?3i sinh?i
L3
3
77
77
77
77
77
77
77
75
2
66
66
66
64
A
B
C
D
3
77
77
77
75
=
2
66
66
66
64
0
0
0
0
3
77
77
77
75
: (5.6)
For non-zero solution, the determinant of the matrix of Eq. (5.6) should be zero as the
following equation
flfl
flfl
flfl
flfl
flfl
flfl
flfl
flfl
flfl
flfl
fl
0 ??
2i
L2 0
?2i
L2
??
2i sin?i
L2 ?
?2i cos?i
L2
?2i sinh?i
L2
?2i cosh?i
L2
??
3i
L3 0
?3i
L3 0
??
3i cos?i
L3
?3i sin?i
L3
?3i cosh?i
L3
?3i sinh?i
L3
flfl
flfl
flfl
flfl
flfl
flfl
flfl
flfl
flfl
flfl
fl
= 0: (5.7)
85
The general solution satisfying these boundary conditions is calculated as
'i = cos ?izL +cosh ?izL ? cosh?i ?cos?isinh?
i ?sin?i
sin ?izL +sinh ?izL
?
: (5.8)
Equation (5.7) is simplifled as
cos?i cosh?i = 1: (5.9)
The roots of this characteristic equation are shown in table 5.1. The pi variable from
Eq. (5.4) is calculated as
pi =
?
i
L
?2 EI
?c Ac
?1=2
: (5.10)
The constants fii and fli of Eq. (5.4) depend upon the initial conditions. The functions
'i(z) satisfy the orthogonality relations
Z L
0
'i 'j dz = L?ij (i;j = 1;:::;1); (5.11)
where ?ij is the Kronecker delta, ?ij =
8>
<
>:
1; if i = j;
0; if i 6= j:
Table 5.1: Roots of characteristic equations for free elastic link
characteristic equation
root (cos?i cosh?i = 1)
?1 4.7300
?2 7.8532
?3 10.996
?4 14.137
86
5.1.2 Modeling
Kinematics
As shown in Fig. 5.2, the system is formed by an uniform exible link of length L with
the diameter dc, the cross sectional area Ac = ?d2c=4, the exural rigidity EI, and the
density ?c. Only planar motions of the link in a flxed reference frame (0) of unit vectors
[?0; ?0; k0] will be considered.
To characterize the instantaneous conflguration of the link, generalized coordinates
qx; qz; andq1 are employed. The generalized coordinate qx denotes the distance from the end
A to the vertical axis of the reference frame (0) and the generalized coordinate qz denotes
the distance from A to the horizontal axis of reference frame (0). The last generalized
coordinate q1 designates the radian measure of the rotation angle between the undeformed
link and the vertical axis. These are the generalized coordinates of the rigid body.
A flxed reference frame (0) of unit vectors [?0;?0;k0] and a mobile reference frame (1)
of unit vectors [?;?;k] are considered. The unit vectors ?0, ?0, and k0 can be expressed as
?0 = cosq1?+sinq1k;
?0 = ?;
k0 = ?sinq1?+cosq1k; (5.12)
and the unit vectors ?, ?, and k can be also expressed as
? = cosq1?0 ?sinq1k0;
? = ?0;
k = sinq1?0 +cosq1k0: (5.13)
The deformations of the elastic link can be discussed in terms of the deformed displace-
ment x(z;t) of a generic point P on the link. The point P is situated at a distance z from
87
the end point A of the link. The displacement x can be expressed as
x(z;t) =
nX
i=1
'i(z)q1+i(t);
=
nX
i=1
?
cos ?izL +cosh ?izL ? cosh?i ?cos?isinh?
i ?sin?i
sin ?izL +sinh ?izL
??
q1+i(t); (5.14)
where 'i(z) is a shape function by z, the elastic generalized coordinate q1+i(t) is a function
of time t, and i is any positive integer.
The position of the end A in (0) is
rA = qx?0 +qz k0: (5.15)
The position vector from the end A to a generic point P(x;z) of the elastic link in (0) is
rAP = x?+zk =
"? nX
i=1
'i(z)q1+i
!
?+zk
#
=
" nX
i=1
'i(z)q1+i cosq1 +z sinq1
#
?0 +
"
?
nX
i=1
'i(z)q1+i sinq1 +z cosq1
#
k0: (5.16)
The position vector of the point P of the elastic link in (0) is
rP = rA +rAP
=
"
qx +
nX
i=1
'i(z)q1+i cosq1 +z sinq1
#
?0 +
"
qz ?
nX
i=1
'i(z)q1+i sinq1 +z cosq1
#
k0: (5.17)
88
The velocity and acceleration vector of the arbitrary point P of the uniform elastic link in
(0) is
vP = ddtrP =
"
_qx +
nX
i=1
'i(z) _q1+i cosq1 ?'i(z)q1+i _q1 sinq1
?
+z _q1 cosq1
#
?0
+
"
_qz ?
nX
i=1
'i(z) _q1+i sinq1 +'i(z)q1+i _q1 cosq1
?
?z _q1 sinq1
#
k0; (5.18)
aP = ddtvP
=
"
?qx +
nX
i=1
'i(z)
?q1+i cosq1 ?2 _q1+i _q1 sinq1 ?q1+i ?q1 sinq1 ?q1+i _q21 cosq1
?
+z ?q1 cosq1 ?z _q21 sinq1
#
?0
+
"
?qz ?
nX
i=1
'i(z)
?q1+i sinq1 +2_q1+i _q1 cosq1 +q1+i ?q1 cosq1 ?q1+i _q21 sinq1
?
?z ?q1 sinq1 ?z _q21 cosq1
#
k0: (5.19)
The angular acceleration of the link in the reference frame (0) is
fi = _! = ?q1?0: (5.20)
5.1.3 Equations of motion
The Newton-Euler?s equations can be used to flnd the difierential equations of motion.
A force and a moment equation before impact can be written as
Z L
0
?c AcaP dz = G; (5.21)
Z L
0
?
rAP ??c AcaP
?
dz = rAC ?G; (5.22)
where G is the gravity force acting on the mass center of the link and rAC = rAP(z = 0:5L)
is the position vector from A to the mass center C. Equation (5.21) is separated into x and
z components.
89
The relationship in dz, a generic difierential element of the link, between the external
force Fext and the shear V can be expressed as
Fext ??dz + @VP(z;t)@z dz = ?c AcaP(z;t)??dz; (5.23)
where VP(z;t) is the shear at a certain point P and all force terms are expressed in flxed
reference frame (0) in terms of mobile frame (1) of ?, ?, k. If the rotatory inertia is neglected,
then V(z;t) may be expressed in terms of the bending moment M(z;t) as
@M(z;t)
@z = ?V(z;t): (5.24)
Since
M = EI@
2x(z;t)
@z2 ; (5.25)
Eqs. (5.23) and (5.24) yield
EI@
4x(z;t)
@z4 +?c Aca(z;t)?? = Fext ??: (5.26)
The variables of Eq. (5.26) are separated by Eq. (5.14) as
EI
1X
i=1
?
i
L
?4
'i(z)q1+i(t)+?c Aca(z;t)?? = Fext ??: (5.27)
Equation (5.27) is simplifled by multiplying 'j(z) and integrating from 0 to L. By the
orthogonality relations represented in Eq. (5.11), the flrst term of Eq. (5.27) becomes
Z L
0
EI
1X
i=1
?
i
L
?4
'i(z)'j(z)q1+i(t)dz = EI L
?
i
L
?4
q1+i(t); (5.28)
and Eq. (5.27) becomes
EI L
?
i
L
?4
q1+i(t)+
Z L
0
?c AcaP ??'i(z)dz =
Z L
0
Fext ??'i(z)dz: (5.29)
90
When the external force Fext is supposed to be a concentrated force on the link, Fext can
be expressed as Fext?(z ? zact) and the right hand side of Eq. (5.29) can be rewritten as
Fext ??'i(zact). In the model, for the existing external force is G, Eq. (5.29) becomes
EI L
?
i
L
?4
q1+i(t)+
Z L
0
?c AcaP ??'i(z)dz = G??'i
L
2
?
: (5.30)
Equations (5.21), (5.22), and (5.30) are the equations of motion for the elastic free link
before the impact.
When the multiple external forces Fextk (k = 1;:::;n) including the resistance force
of the granular medium, the governing equations of motion, Eqs. (5.21), (5.22), and (5.30),
can be written as
Z L
0
?c AcaP dz = G+
nX
k=1
Fextk; (5.31)
Z L
0
(rAP ??c AcaP) dz = rAC ?G+
nX
k=1
rAEk ?Fextk; (5.32)
EI L
?
i
L
?4
q1+i(t) +
Z L
0
?c AcaP ??'i(z)dz =
G??'i
L
2
?
+
nX
k=1
Fextk ??'i(LEk); (5.33)
where rAEk is the position vector from the end A to the external force Fextk acting point Ek
and LEk is the length from the end A to the external force application point Ek. In the case
of the impact of the elastic free link, external forces are restricted to the resistance force
FR. The equations of motion of the elastic free link during the impact with the granular
medium are given as
Z L
0
?c AcaP dz = G+FR; (5.34)
Z L
0
(rAP ??c AcaP) dz = rAC ?G+rAE ?FR; (5.35)
EI L
?
i
L
?4
q1+i(t)+
Z L
0
?cAcaP ??'i(z)dz = G??'i
L
2
?
+FR ??'i(LE);(5.36)
91
where rAE = rAP(z = LE) is the position vector from A to the resistance force application
point E and the length from the end A to E, LE, is calculated as
LE = L? zT2cosq
1
: (5.37)
The immersed depth zT, the vertical component of the position vector of the end point T,
is calculated by Eq. (5.17) as zT = rP (z = L)?k0.
The resistance force FR is calculated by Eqs. (3.6), (3.8), and (3.9). The dynamic
frictional force Fd is
Fd = ?vE ?d ?g Ar jvEj
= ??d ?g Ar
"?
_qx +
nX
i=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+LE _q1 cosq1
!2
+
?
_qz ?
nX
i=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?LE _q1 sinq1
!2#0:5
"?
_qx +
nX
i=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+LE _q1 cosq1
!
?0 +
?
_qz ?
nX
i=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?LE _q1 sinq1
!
k0
#
; (5.38)
where the velocity vector vE is the velocity of the resistance force acting point E and this
vector calculated by Eq. (5.18) as vE = vP (z = LE). The reference area Ar is
Ar = dc zTcosq
1
j sin(q1 ?qm) j : (5.39)
92
The moving angle qm is calculated as
qm = tan?1
v
Ex
vEz
?
(5.40)
= tan?1
2
66
4
_qx +Pni=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+LE _q1 cosq1
_qz ?Pni=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?LE _q1 sinq1
3
77
5:
The horizontal and vertical static resistance force Fsh, Fvh are
Fsh = ?sign
"
_qx +
nX
i=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+LE _q1 cosq1
#
?h g?g z2T dc ?0; (5.41)
Fsv = ?sign
"
_qz ?
nX
i=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?LE _q1 sinq1
#
?v
z
T
dc
??
g?g V k0; (5.42)
where the immersed volume V is calculated as
V = ?dc
2
4
zT
cosq1: (5.43)
The resistance force FR, the sum of the dynamic frictional force vector Fd and the static
resistance force vector Fs, is represented as the sum of Eqs. (5.38), (5.41), and (5.42).
5.1.4 Simulation results
Figures 5.3, 5.4, 5.5, and 5.6 show the simulation results for the impact of the elastic
free link shown in Fig. 5.2. The simulations are performed for difierent impact angles and
difierent initial impact velocities (_qz(0) = 1.26, 1.87, and 2.33 m/s for q1(0) = 32o, and
_qz(0) = 1.45, 1.98, and 2.43 m/s for q1(0) = 55o). The initial deformation conditions are
93
Table 5.2: Stopping time of the elastic and the rigid free link
q1(0) (o) _qz(0) (m/s) ts;r (s) ts;e (s)
1.26 0.0275914 0.0275865
32 1.87 0.0229915 0.0229956
2.33 0.0206495 0.0206516
1.45 0.017818 0.0178181
55 1.98 0.0150461 0.0150411
2.43 0.0133007 0.0133043
applied as q2(0) = 0 m and _q2(0) = 0 m/s. The dimensions of the elastic free link applied
to this simulation are completely the same as those applied for the simulation of the impact
of the rigid free link as L = 0:1524 m, the diameter dc = 0:00635 m, and the density
?c = 7:7?103 kg/m3. The exural rigidity EI is applied as 15.1642 Nm2. The density of
a granular medium (sand), ?g, is applied as 2:5?103 kg/m3. The dynamic frictional force
coe?cient is ?d = 6.5, the horizontal static resistance force coe?cient is ?h = 8, and the
vertical static resistance force coe?cients are ?v = 22 and ? = 1.1. The simulations are
performed from the impact moment until the vertical penetrating velocity of the end T of
the link, vTz, becomes zero and the flrst mode of the shape function of Eq. (5.8) is only
considered (?1=4.7300). As shown in Figs. 5.3 and 5.5, the velocity of the end T of the link,
vTz, becomes zero more quickly when the initial impact velocity increases. The simulation
results of the elastic free link compared with the results of the rigid free link are shown
in table 5.2. These simulation results show that the stopping time of the impact of a free
link decreases whether a rigid link or an elastic one as the velocity of the impact moment
increases. The deformation of the link, q2, increases when the initial impact velocity and
angle increase as shown in Figs. 5.4 and 5.6. We do not observe much difierence between
the rigid model and the exible model for our particular system. The elastic deformations
will be important for longer and more elastic links and multiple impacts with the sand
where the vibrations are more important.
94
vT
z
(m/s)
zT
(m)
Time (s)
Time (s)
0.000 0.005 0.010 0.015 0.020 0.025
0.000
0.005
0.010
0.015
0.020
0.025
0.000 0.005 0.010 0.015 0.020 0.025
0.0
0.5
1.0
1.5
2.0 vTz (0) = 1.26 m/s
vTz (0) = 1.87 m/s
(0) = 2.33 m/svTz
Figure 5.3: Displacement zT and velocity vTz results of the exible link for q1(0) = 32o and
_qz(0) = 1.26, 1.87, 2.33 m/s
95
q 2
(m)
Tim e (s)
Tim e (s)
?q 2
(m
/s)
0.000 0.005 0.010 0.015 0.020 0.025
0
2. 10 -7
4. 10 -7
6. 10 -7
8. 10 -7
1. 10 -6
0.000 0.005 0.010 0.015 0.020 0.025
-0.0005
0.0000
0.0005
0.0010
vTz (0) = 1.26 m/s
vTz (0) = 1.87 m/s
(0) = 2.33 m/svTz
Figure 5.4: Displacement q2 and velocity _q2 results of the exible link for q1(0) = 32o and
_qz(0) = 1.26, 1.87, 2.33 m/s
96
vT
z
(m/s)
zT
(m)
Time (s)
Time (s)
0.000 0.005 0.010 0.015
0.0
0.5
1.0
1.5
2.0
0.000 0.005 0.010 0.015
0.000
0.005
0.010
0.015
vTz (0) = 1.45 m/s
vTz (0) = 1.98 m/s
(0) = 2.43 m/svTz
Figure 5.5: Displacement zT and velocity vTz results of the exible link for q1(0) = 55o and
_qz(0) = 1.45, 1.98, 2.43 m/s
97
q 2
(m)
Tim e (s)
Tim e (s)
?q 2
(m
/s)
2.5
2.
1.5
1.
5. 10 -7
10 -6
10 -6
10 -6
10 -6
0.000 0.005 0.010 0.015
-0.001
0.000
0.001
0.002
0.003
vTz (0) = 1.45 m/s
vTz (0) = 1.98 m/s
(0) = 2.43 m/svTz
Figure 5.6: Displacement q2 and velocity _q2 results of the exible link for q1(0) = 55o and
_qz(0) = 1.45, 1.98, 2.43 m/s
98
5.2 Elastic compound pendulum
5.2.1 Transverse vibration of an elastic compound pendulum
In the case of an articulated elastic pendulum composed of one end supported by a
joint as shown in Fig. 5.7, the boundary conditions for this pendulum (articulated at one
end free at the other end) can be written as
x(0;t) = x00(0;t) = x00(L;t) = x000(L;t) = 0: (5.44)
Using the boundary conditions given by Eq. (5.44), the relation between constants A, B,
C, and D can be attained as
2
66
66
66
66
66
66
66
66
64
0 1 0 1
0 ??
2i
L2 0
?2i
L2
??
2i sin?i
L2 ?
?2i cos?i
L2
?2i sinh?i
L2
?2i cosh?i
L2
??
3i cos?i
L3
?3i sin?i
L3
?3i cosh?i
L3
?3i sinh?i
L3
3
77
77
77
77
77
77
77
77
75
2
66
66
66
64
A
B
C
D
3
77
77
77
75
=
2
66
66
66
64
0
0
0
0
3
77
77
77
75
: (5.45)
For non-zero solution, the determinant of the matrix of Eq. (5.45) should be zero and the
general solution satisfying these boundary conditions is calculated as
'i = sinh?isin?
i
sin ?ixL +sinh ?ixL : (5.46)
The characteristic equation of which roots satisfy Eq. (5.46) is calculated as
cos?i sinh?i ?sin?i cosh?i = 0: (5.47)
The roots of this characteristic equation are shown in table 5.3.
99
k0
0?
0?
(0) O
z x
P
E
L
C
q1
k
?
(1)
G
RF
Granular medium
dc
Figure 5.7: Articulated elastic compound pendulum
k0
0?
0?(0)
z x
P
E
L
C
1RB
B
2RB
q1
k
? (1)
G
RF
Granular medium
dc
O
LRB
Lp1
Lp2
Lp3
dp1
dp2
dp3
Figure 5.8: Cantilevered elastic compound pendulum
100
In the case of the cantilevered elastic compound pendulum composed of one end sup-
ported by a rigid body shown in Fig. 5.8, the boundary conditions (cantilevered at one end
and free at the other end) can be written as
x(0;t) = x0(0;t) = x00(L;t) = x000(L;t) = 0: (5.48)
Using the boundary conditions given by Eq. (5.48), the relation between constants A, B,
C, and D can be attained as
2
66
66
66
66
66
66
66
66
64
0 1 0 1
?i
L 0
?i
L 0
??
2i sin?i
L2 ?
?2i cos?i
L2
?2i sinh?i
L2
?2i cosh?i
L2
??
3i cos?i
L3
?3i sin?i
L3
?3i cosh?i
L3
?3i sinh?i
L3
3
77
77
77
77
77
77
77
77
75
2
66
66
66
64
A
B
C
D
3
77
77
77
75
=
2
66
66
66
64
0
0
0
0
3
77
77
77
75
: (5.49)
From Eq. (5.49), the constants A, B, C, and D can be calculated and the general solution
satisfying these boundary condition is given as
'i = cosh ?izL ?cos ?izL ? cosh?i +cos?isinh?
i +sin?i
sinh ?izL ?sin ?izL
?
: (5.50)
The characteristic equation of which roots satisfy Eq. (5.50) is given as
cos?i cosh?i +1 = 0: (5.51)
The roots of this characteristic equation are shown in table 5.3.
101
Table 5.3: Roots of characteristic equations for elastic compound pendulum
characteristic equation of characteristic equation of
root articulated elastic pendulum cantilevered elastic pendulum
(cos?i sinh?i ?sin?i cosh?i = 0) (cos?i cosh?i +1 = 0)
?1 3.9266 1.8751
?2 7.0686 4.6941
?3 10.210 7.8548
?4 13.351 10.996
5.2.2 Modeling of an articulated elastic compound pendulum
As shown in Fig. 5.7, the system is formed by an uniform exible link of length L
with diameter dc, the cross sectional area Ac = ?d2c=4, the exural rigidity EI, and the
density ?c. To characterize the instantaneous conflguration of the pendulum, generalized
coordinates q1 is employed. The generalized coordinate q1 denotes the radian measure of the
rotation angle between the undeformed pendulum and the vertical axis. A flxed reference
frame (0) of unit vectors [?0;?0;k0] and a mobile reference frame (1) of unit vectors [?;?;k]
are considered. The unit vectors ?0, ?0, and k0 can be expressed as Eq. (5.12) and the unit
vectors ?, ?, and k can be also expressed as Eq. (5.13).
Kinematics
The deformations of the articulated elastic compound pendulum can be discussed in
terms of the elastic displacement x(z;t) of a generic point P on the pendulum. The point
P is situated at a distance z from the origin point O. The displacement x can be expressed
as
x(z;t) =
nX
i=1
'i(z)q1+i(t);
=
nX
i=1
?sinh?
i
sin?i sin
?ix
L +sinh
?ix
L
?
q1+i(t); (5.52)
102
where 'i(z) is a shape function by the distance z, the elastic generalized coordinate q1+i(t)
is a function of time t, and i is any positive integer for the articulated elastic compound
pendulum. The position vector from origin point O to the generic point P of the elastic
link in (0) is
rP = x?+zk =
"? nX
i=1
'i(z)q1+i
!
?+zk
#
=
" nX
i=1
'i(z)q1+i cosq1 +z sinq1
#
?0 +
"
?
nX
i=1
'i(z)q1+i sinq1 +z cosq1
#
k0: (5.53)
The velocity vector of the arbitrary point P of the uniform elastic link in (0) is
vP = ddtrP =
" nX
i=1
'i(z) _q1+i cosq1 ?'i(z)q1+i _q1 sinq1
?
+z _q1 cosq1
#
?0
+
"
?
nX
i=1
'i(z) _q1+i sinq1 +'i(z)q1+i _q1 cosq1
?
?z _q1 sinq1
#
k0; (5.54)
and the acceleration vector of the arbitrary point P in (0) is
aP = ddtvP
=
" nX
i=1
'i(z)
?q1+i cosq1 ?2 _q1+i _q1 sinq1 ?q1+i ?q1 sinq1 ?q1+i _q21 cosq1
?
+z ?q1 cosq1 ?z _q21 sinq1
#
?0
+
"
?
nX
i=1
'i(z)
?q1+i sinq1 +2_q1+i _q1 cosq1 +q1+i ?q1 cosq1 ?q1+i _q21 sinq1
?
?z ?q1 sinq1 ?z _q21 cosq1
#
k0: (5.55)
The angular acceleration of the link in the reference frame (0) is
fi = _! = ?q1?0: (5.56)
103
Equations of motion
The Newton-Euler?s equations can be used to flnd the difierential equations of motion.
A moment equation before impact can be written as
Z L
0
(rP ??c AcaP)dz = rC ?G; (5.57)
where G is the gravity force acting on the mass center of the pendulum and rC = rP(z =
0:5L) is the position vector from the origin O to the mass center C. One more governing
equation of motion regarding the deformation of the pendulum has the same form with
Eq. (5.29) as
EI L
?
i
L
?4
q1+i(t)+
Z L
0
?c AcaP ??'i(z)dz =
Z L
0
Fext ??'i(z)dz: (5.58)
Because the external forces acting on the pendulum are the gravitational force G at the
mass center C and the joint reaction force FO at the origin O, the external forces Fext
can be expressed as G?(z ? 0:5L) + FO ?(z) when the external force Fext is supposed
to be concentrated one on the pendulum. The right hand side of Eq. (5.58) becomes
G ? ?'i(0:5L) + FO ? ?'i(0). However, the value of 'i(0) is always zero in the model.
Equation (5.58) becomes
EI L
?
i
L
?4
q1+i(t)+
Z L
0
?c AcaP ??'i(z)dz = G??'i
L
2
?
: (5.59)
Equations (5.57) and (5.59) represent the equations of motion for the articulated elastic
compound pendulum before the impact.
During impacting a granular medium, additionally added external force is restricted
to the resistance force FR. The equations of motion of the articulated elastic pendulum
104
during the impact are given as
Z L
0
(rP ??c AcaP)dz = rC ?G+rE ?FR; (5.60)
EI L
?
i
L
?4
q1+i(t)+
Z L
0
?cAcaP ??'i(z)dz =
G??'i
L
2
?
+ FR ??'i(LE); (5.61)
where rE = rP(z = LE) is the position vector from the origin O to the resistance force
application point E and the length from the origin O to E, LE, is calculated as
LE = L? zT2cosq
1
: (5.62)
The immersed depth zT, the vertical component of the position vector of the end point T, is
calculated by Eq. (5.53) as zT = rP (z = L)?k0. The resistance force FR is also calculated
by Eqs. (3.6), (3.8), and (3.9). The dynamic frictional force Fd is
Fd = ?vE ?d ?g Ar jvEj
= ??d ?g Ar
"? nX
i=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+LE _q1 cosq1
!2
+
?
?
nX
i=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?LE _q1 sinq1
!2#0:5
"?
+
nX
i=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+LE _q1 cosq1
!
?0 +
?
?
nX
i=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?LE _q1 sinq1
!
k0
#
; (5.63)
where the velocity vector vE is the velocity of the resistance force acting point E and this
vector calculated by Eq. (5.54) as vE = vP (z = LE). The reference area Ar is calculated
105
as
Ar = dc zTcosq
1
: (5.64)
The horizontal and vertical static resistance forces, Fsh and Fvh, are
Fsh = ?sign
" nX
i=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+LE _q1 cosq1
#
?h g?g z2T dc ?0; (5.65)
Fsv = ?sign
"
?
nX
i=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?LE _q1 sinq1
#
?v
z
T
dc
??
g?g V k0; (5.66)
where the immersed volume V is calculated as
V = ?dc
2
4
zT
cosq1: (5.67)
The resistance force FR, the sum of the dynamic frictional force vector Fd and the static
resistance force vector Fs, is represented as the sum of Eqs. (5.63), (5.65), and (5.66).
5.2.3 Simulation results of the articulated elastic compound pendulum
Figures 5.9, 5.10, 5.11, 5.12, 5.13, and 5.14 show the simulation results for the elastic
compound pendulum as shown in Fig. 5.7 with the following initial conditions: q1(0) = 30o,
60o, 75o, _q1(0) = -1, -3, -5 rad/s, q2(0) = 0 m, and _q2(0) = 0 m/s. As the numerical
data applied to this simulation, the length of pendulum is L = 0:15 m, the diameter is
dc = 0:00635 m, the density of the pendulum is ?c = 7:7?103 kg/m3, the exural rigidity
is EI = 15:1642 Nm2, and the density of a granular medium (sand) is ?g = 2:5 ? 103
kg/m3. The dynamic frictional force coe?cient ?d = 6.5, the horizontal static resistance
106
?q 1
(rad/
s)
Tim e (s)
Tim e (s)
0.00 0.02 0.04 0.06 0.08
-5
-4
-3
-2
-1
0
0.00 0.02 0.04 0.06 0.08
22
24
26
28
30
q 1
(o
)
?q1(0) = ? 1 rad/ s
?q1(0) = ? 3 rad/ s
?q1(0) = ? 5 rad/ s
Figure 5.9: Angle q1 and angular velocity _q1 results of the articulated elastic compound
pendulum for q1(0) = 30o and _q1(0) = -1, -3, -5 rad/s
107
q 2
(m)
Tim e (s)
Tim e (s)
?q 2
(m
/s)
?q1(0) = ? 1 rad/ s
?q1(0) = ? 3 rad/ s
?q1(0) = ? 5 rad/ s
0.00 0.02 0.04 0.06 0.08
0
10 -6
0.00001
0.000015
0.00002
0.000025
0.00 0.02 0.04 0.06 0.08
-0.0010
-0.0005
0.0000
0.0005
0.0010
0.0015
0.0020
Figure 5.10: Displacement q2 and velocity _q2 results of the articulated elastic compound
pendulum for q1(0) = 30o and _q1(0) = -1, -3, -5 rad/s
108
?q 1
(rad/
s)
Tim e (s)
Tim e (s)
q 1
(o
)
?q1(0) = ? 1 rad/ s
?q1(0) = ? 3 rad/ s
?q1(0) = ? 5 rad/ s
0.00 0.01 0.02 0.03 0.04 0.05
-5
-4
-3
-2
-1
0
0.00 0.01 0.02 0.03 0.04 0.05
55
56
57
58
59
60
Figure 5.11: Angle q1 and angular velocity _q1 results of the articulated elastic compound
pendulum for q1(0) = 60o and _q1(0) = -1, -3, -5 rad/s
109
q 2
(m)
Tim e (s)
Tim e (s)
?q 2
(m/s)
?q1(0) = ? 1 rad/ s
?q1(0) = ? 3 rad/ s
?q1(0) = ? 5 rad/ s
0.00 0.01 0.02 0.03 0.04 0.05
0
0.00001
0.00002
0.00003
0.00004
0.00 0.01 0.02 0.03 0.04 0.05
-0.002
-0.001
0.000
0.001
0.002
0.003
0.004
Figure 5.12: Displacement q2 and velocity _q2 results of the articulated elastic compound
pendulum for q1(0) = 60o and _q1(0) = -1, -3, -5 rad/s
110
?q 1
(rad/
s)
Tim e (s)
Tim e (s)
q 1
(o
)
?q1(0) = ? 1 rad/ s
?q1(0) = ? 3 rad/ s
?q1(0) = ? 5 rad/ s
0.00 0.01 0.02 0.03 0.04
-5
-4
-3
-2
-1
0
0.00 0.01 0.02 0.03 0.04
71
72
73
74
75
Figure 5.13: Angle q1 and angular velocity _q1 results of the articulated elastic compound
pendulum for q1(0) = 75o and _q1(0) = -1, -3, -5 rad/s
111
q 2
(m)
Tim e (s)
Tim e (s)
?q 2
(m/
s)
?q1(0) = ? 1 rad/ s
?q1(0) = ? 3 rad/ s
?q1(0) = ? 5 rad/ s
0.00 0.01 0.02 0.03 0.04
0
0.00001
0.00002
0.00003
0.00004
0.00005
0.00 0.01 0.02 0.03 0.04
-0.004
-0.002
0.000
0.002
0.004
Figure 5.14: Displacement q2 and velocity _q2 results of the articulated elastic compound
pendulum for q1(0) = 75o and _q1(0) = -1, -3, -5 rad/s
112
Table 5.4: Stopping time of the articulated elastic compound pendulum and the rigid one
q1(0) (o) _q1(0) (rad/s) ts;r (s) ts;e (s)
-1 0.0908846 0.0946034
30 -3 0.068481 0.0700719
-5 0.0574613 0.0573489
-1 0.047826 0.0509171
60 -3 0.0353842 0.0372156
-5 0.0292904 0.0299783
-1 0.0371793 0.0403826
75 -3 0.0270217 0.0291454
-5 0.0223565 0.024833
force coe?cient ?h = 8, and the vertical static resistance force coe?cients ?v = 22 and ?
= 1.1 are used in the simulation. The simulations are performed from the impact moment
until the angular velocity of the pendulum, _q1, becomes zero and the flrst mode of the shape
function of Eq. (5.46) is only considered (?1=3.9266). As shown in Figs. 5.9, 5.11, and 5.13,
the pendulum stops more quickly when the pendulum impact the granular medium with
high velocity. This result means that the faster the impact velocity becomes the sooner
the impact compound pendulum stops whether a rigid pendulum or an elastic one. The
generalized coordinate q2 showing the deformation of the pendulum becomes large when the
initial impact velocity and angle increase as shown in Figs. 5.10, 5.12, and 5.14. Table 5.4
shows the compared results between the rigid and the elastic model. The difierences between
the rigid model and the exible model are not also conspicuous for our particular system.
The articulated elastic deformations of the compound pendulum will be important for longer
and more elastic links.
113
5.2.4 Modeling of a cantilevered elastic compound pendulum
AsshowninFig.5.8, thecantileveredelasticpendulumappliedtomodelingissupported
by rigid bodies RB1 and RB2. The rigid body RB1 has the diameter dp1 and the length
Lp1. The rigid body RB2 has the diameter dp2 and the length Lp2. The exible part of
the pendulum, B, has an uniform exible link of the length Lp3 = LB = L, the diameter
dp3 = dc, the cross sectional area Ac = ?d2c=4, the exural rigidity EI, and density ?c.
The length from the origin to the rigid body end is LRB = dp1=2 + Lp2 and the density of
the rigid bodies are the same as that of the elastic link. To characterize the instantaneous
conflguration of the pendulum, the generalized coordinates q1 and qi+1 are employed. The
generalized coordinate q1 denotes the radian measure of the rotation angle between the
undeformed pendulum and the horizontal axis. The generalized coordinate qi + 1 denotes
the deformation of the elastic link part of the pendulum. A flxed reference frame (0) of unit
vectors [?0;?0;k0] and a mobile reference frame (1) of unit vectors [?;?;k] are considered.
The unit vectors ?0, ?0, and k0 can be expressed as Eq. (5.12) and the unit vectors ?, ?, and
k can be also expressed as Eq. (5.13).
Kinematics
The deformations of the cantilevered elastic compound pendulum can be also discussed
in terms of the elastic displacement x(z;t) of a generic point P on the elastic part B of the
pendulum. The point P is situated at a distance z from the point at the flrst end of the
elastic link B. The displacement x can be expressed as
x(z;t) =
nX
i=1
'i(z)q1+i(t);
=
nX
i=1
?
cosh ?izL ?cos ?izL ? cosh?i +cos?isinh?
i +sin?i
sinh ?izL ?sin ?izL
??
q1+i(t); (5.68)
where 'i(z) is a shape function of the elastic link by z, the elastic generalized coordinate
q1+i(t) is a function of time t, and i is any positive integer.
114
The position of the mass center of the rigid body RB1 is the origin point O and the position
vector from the origin O to the mass center of the rigid body RB2, CRB2, is represented as
rCRB2 = dp1 +Lp22
?
sinq1?0 +cosq1k0
?
: (5.69)
The position vector from the origin point O to a generic point P of the link B in (0) is
rP = x?+(LRB +z) k =
"? nX
i=1
'i(z)q1+i
!
?+(LRB +z) k
#
=
" nX
i=1
'i(z)q1+i cosq1 +(LRB +z) sinq1
#
?0 +
"
?
nX
i=1
'i(z)q1+i sinq1 +(LRB +z) cosq1
#
k0: (5.70)
The velocity vector of the generic point P of the elastic link B in (0) is
vP = ddtrP =
" nX
i=1
'i(z) _q1+i cosq1 ?'i(z)q1+i _q1 sinq1
?
+(LRB +z) _q1 cosq1
#
?0
+
"
?
nX
i=1
'i(z) _q1+i sinq1 +'i(z)q1+i _q1 cosq1
?
?(LRB +z) _q1 sinq1
#
k0;
(5.71)
and the acceleration vector of the point P in (0) is
aP = ddtvP
=
" nX
i=1
'i(z)
?q1+i cosq1 ?2 _q1+i _q1 sinq1 ?q1+i ?q1 sinq1 ?q1+i _q21 cosq1
?
+(LRB +z) ?q1 cosq1 ?(LRB +z) _q21 sinq1
#
?0
+
"
?
nX
i=1
'i(z)
?q1+i sinq1 +2_q1+i _q1 cosq1 +q1+i ?q1 cosq1 ?q1+i _q21 sinq1
?
?(LRB +z) ?q1 sinq1 ?(LRB +z) _q21 cosq1
#
k0: (5.72)
115
The angular acceleration of the link in the reference frame (0) is
fi = _! = ?q1?0: (5.73)
Equations of motion
The Newton-Euler?s equations can be used to flnd the difierential equation of motion.
A moment equation before impact can be written as
IRB ?q1 +
Z L
0
(rP ??c AcaP) dz = rCRB2 ?GRB2 +rC ?GB; (5.74)
where IRB is the mass moment of inertia of the rigid body and the gravity forces GRB2
and GB are acting on the mass center of the rigid body RB2 and on the exible link B
respectively. The vector rC is the position vector from origin point O to the mass center C
of the exible link (rC = rP(z = 0:5L)). One more governing equation of motion regarding
deformation of the link has the same form as Eq. (5.29) and the equation becomes
EI L
?
i
L
?4
q1+i(t)+
Z L
0
?c AcaP ??'i(z)dz =
Z L
0
Fext ??'i(z)dz: (5.75)
Because the external forces acting on the elastic link of the pendulum are the gravity force
GB and the joint reaction force FRB2 from the rigid body RB2, the external force Fext can
be written as GB ?(z ?0:5L) + FRB2 ?(z) when the external force Fext is supposed to be
concentrated one on the pendulum. Therefore the right hand side of Eq. (5.75) becomes
GB ? ?'i(0:5L) + FO ? ?'i(0). Because the value of 'i(0) is always zero in the model,
Eq. (5.75) becomes
EI L
?
i
L
?4
q1+i(t)+
Z L
0
?c AcaP ??'i(z)dz = GB ??'i
L
2
?
: (5.76)
Difierentialequations(5.74)and(5.76)representtheequationsofmotionforthecantilevered
elastic pendulum before impact.
116
During impacting a granular medium, an additionally added external force is restricted
to the resistance force FR. The equations of motion for the cantilevered elastic pendulum
during impacting are given as
IRB ?q1 +
Z L
0
(rP ??c AcaP) dz = rC ?G+rE ?FR; (5.77)
EI L
?
i
L
?4
q1+i(t)+
Z L
0
?cAcaP ??'i(z)dz =
GB ??'i
L
2
?
+ FR ??'i(LE); (5.78)
where rE = rP(z = LE) is the position vector from the origin point O to the resistance force
application point E. The length from the origin point O to the resistance force application
point E, LE, is calculated as
LE = L? zT2cosq
1
: (5.79)
The immersed depth zT, the vertical component of the position vector of the end point T, is
calculated by Eq. (5.70) as zT = rP (z = L)?k0. The resistance force FR is also calculated
by Eqs. (3.6), (3.8), and (3.9). The dynamic frictional force Fd is
Fd = ?vE ?d ?g Ar jvEj
= ??d ?g Ar
"? nX
i=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+(LRB +LE) _q1 cosq1
!2
+
?
?
nX
i=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?(LRB +LE) _q1 sinq1
!2#0:5
"?
+
nX
i=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+(LRB +LE) _q1 cosq1
!
?0 +
?
?
nX
i=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?(LRB +LE) _q1 sinq1
!
k0
#
; (5.80)
117
where the velocity vector vE is the velocity of the resistance force acting point E and this
vector calculated by Eq. (5.42) as vE = vP (z = LE). The reference area Ar is represented
as
Ar = dc zTcosq
1
: (5.81)
The horizontal and vertical static resistance forces, Fsh and Fvh, are
Fsh = ?sign
" nX
i=1
'i(LE) _q1+i cosq1 ?'i(LE)q1+i _q1 sinq1
?
+(LRB +LE) _q1 cosq1
#
?h g?g z2T dc ?0; (5.82)
Fsv = ?sign
"
?
nX
i=1
'i(LE) _q1+i sinq1 +'i(LE)q1+i _q1 cosq1
?
?(LRB +LE) _q1 sinq1
#
?v
z
T
dc
??
g?g V k0; (5.83)
where the immersed volume V is calculated as
V = ?dc
2
4
zT
cosq1: (5.84)
The resistance force FR, the sum of the dynamic frictional force vector Fd and the static
resistance force vector Fs, is represented as the sum of Eqs. (5.80), (5.82), and (5.83).
5.2.5 Simulation results of the cantilevered elastic compound pendulum
Figures 5.15, 5.16, 5.17, 5.18, 5.19, 5.20, 5.21, and 5.22 show the simulation results for
the cantilevered elastic compound pendulum shown in Fig. 5.8. The dimensions applied
to this simulation are the same as those applied for the simulation of the rigid compound
pendulum. The exural rigidity EI = 15:1642 Nm2 are utilized and the other parameters
are the same as those of the previous simulation. The simulations are performed for difierent
118
Tim e (s)
Tim e (s)
?q (0) = ? 1.75 rad/s
?q (0) = ? 3.38 rad/s
?q (0) = ? 4.66 rad/s
0.00 0.02 0.04 0.06 0.08 0.10
-4
-3
-2
-1
0
0.00 0.02 0.04 0.06 0.08 0.10
14
16
18
20
22
?q 1
(rad/s)
q 1
(o
)
1
1
1
Figure 5.15: Angle q1 and angular velocity _q1 results of the cantilevered elastic compound
pendulum for q1(0) = 22o and _q1(0) = -1.75, -3.38, -4.66 rad/s
119
Tim e (s)
Tim e (s)
?q (0) = ? 1.75 rad/ s
?q (0) = ? 3.38 rad/ s
?q (0) = ? 4.66 rad/ s
?q 2
(m/
s)
(m)
q 2
1
1
1
0.00 0.02 0.04 0.06 0.08 0.10
0
0.00001
0.00002
0.00003
0.00004
0.00 0.02 0.04 0.06 0.08 0.10
-0.005
0.000
0.005
0.010
Figure 5.16: Displacement q2 and velocity _q2 results of the cantilevered elastic compound
pendulum for q1(0) = 22o and _q1(0) = -1.75, -3.38, -4.66 rad/s
120
Tim e (s)
Tim e (s)
?q (0) = ? 3.3 1 rad/s
?q (0) = ? 6.24 rad/s
?q (0) = ? 8.41 rad/s
?q 1
(rad/s)
q 1
(o
)
1
1
1
0.00 0.01 0.02 0.03 0.04 0.05 0.06
-8
-6
-4
-2
0
0.00 0.01 0.02 0.03 0.04 0.05 0.06
24
26
28
30
Figure 5.17: Angle q1 and angular velocity _q1 results of the cantilevered elastic compound
pendulum for q1(0) = 31o and _q1(0) = -3.31, -6.24, -8.41 rad/s
121
Tim e (s)
Tim e (s)
?q (0) = ? 3.3 1 rad/s
?q (0) = ? 6.24 rad/s
?q (0) = ? 8.41 rad/s
?q 2
(m/
s)
(m)
q 2
1
1
1
0.00 0.01 0.02 0.03 0.04 0.05 0.06
-0.02
0.00
0.02
0.04
0.06
0.00 0.01 0.02 0.03 0.04 0.05 0.06
0.00000
0.00005
0.00010
0.00015
Figure 5.18: Displacement q2 and velocity _q2 results of the cantilevered elastic compound
pendulum for q1(0) = 31o and _q1(0) = -3.31, -6.24, -8.41 rad/s
122
Tim e (s)
Tim e (s)
?q (0) = ? 2.66 rad/s
?q (0) = ? 6.47 rad/s
?q (0) = ? 9.06 rad/s
?q 1
(rad/s)
q 1
(o
)
1
1
1
0.00 0.01 0.02 0.03 0.04
39
40
41
42
43
44
45
0.00 0.01 0.02 0.03 0.04
-8
-6
-4
-2
0
Figure 5.19: Angle q1 and angular velocity _q1 results of the cantilevered elastic compound
pendulum for q1(0) = 45o and _q1(0) = -2.66, -6.47, -9.06 rad/s
123
Tim e (s)
Tim e (s)
?q (0) = ? 2.66 rad/s
?q (0) = ? 6.47 rad/s
?q (0) = ? 9.06 rad/s
?q 2
(m
/s)
(m)
q 2
1
1
1
0.00 0.01 0.02 0.03 0.04
0.00000
0.00005
0.00010
0.00015
0.00020
0.00025
0.00 0.01 0.02 0.03 0.04
-0.05
0.00
0.05
0.10
Figure 5.20: Displacement q2 and velocity _q2 results of the cantilevered elastic compound
pendulum for q1(0) = 45o and _q1(0) = -2.66, -6.47, -9.06 rad/s
124
Tim e (s)
Tim e (s)
?q (0) = ? 2.70 rad/s
?q (0) = ? 6.54 rad/s
?q (0) = ? 9.17 rad/s
?q 1
(rad/s)
q 1
(o
)
1
1
1
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
57
58
59
60
61
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
-8
-6
-4
-2
0
Figure 5.21: Angle q1 and angular velocity _q1 results of the cantilevered elastic compound
pendulum for q1(0) = 61:5o and _q1(0) = -2.70, -6.54, -9.17 rad/s
125
Tim e (s)
Tim e (s)
?q (0) = ? 2.7 0 rad/s
?q (0) = ? 6.54 rad/s
?q (0) = ? 9.17 rad/s
?q 2
(m
/s)
(m)
q 2
1
1
1
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
-0.10
-0.05
0.00
0.05
0.10
0.15
0.20
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
0.0000
0.0001
0.0002
0.0003
Figure 5.22: Displacement q2 and velocity _q2 results of the cantilevered elastic compound
pendulum for q1(0) = 61:5o and _q1(0) = -2.70, -6.54, -9.17 rad/s
126
impact angles and difierent initial impact angular velocities (_q(0) = -1.75, -3.38, and -4.66
rad/s for q(0) = 22o, _q(0) = -3.31, -6.24, and -8.41 rad/s for q(0) = 31o, _q(0) = -2.66,
-6.47, and -9.06 rad/s for q(0) = 45o, and _q(0) = -2.70, -6.54, and -9.17 rad/s for q(0) is
61.5o). from the impact moment until the angular velocity of the pendulum, _q1, becomes
zero and the flrst mode of the shape function of Eq. (5.50) is only considered (?1=1.8751).
As shown in Figs. 5.15, 5.17, 5.19, and 5.21, the pendulum stops more quickly when the
pendulum impacts the granular medium with high velocity. The results that the faster the
impact velocity becomes the sooner the compound pendulum stops are kept in the impact
of the cantilevered elastic compound pendulum. The generalized coordinate q2 showing
the deformation of the pendulum becomes large when the initial impact velocity and angle
increase as shown in Figs. 5.16, 5.18, 5.20, and 5.22. Table 5.5 shows the stopping time
results for the rigid model, ts;r, and the elastic model, ts;e, and we do not observe much
Table 5.5: Stopping time of the cantilevered elastic compound pendulum and the rigid one
q1(0) (o) _q1(0) (rad/s) ts;r (s) ts;e (s)
-1.75 0.107212 0.107366
22 -3.38 0.0861789 0.0861761
-4.66 0.0775028 0.0774967
-3.31 0.0609829 0.0610354
31 -6.24 0.0485081 0.0485343
-8.41 0.043819 0.0438373
-2.66 0.046688 0.046748
45 -6.47 0.0333638 0.0333414
-9.06 0.0295578 0.029574
-2.70 0.0343956 0.0345101
61.5 -6.54 0.0245294 0.0245555
-9.17 0.0216994 0.021729
127
difierence between the rigid model and the exible model for our particular system. The
elastic deformations of the cantilever elastic compound pendulum will be important for
longer and more elastic links.
128
Chapter 6
Experiments
In this chapter, the experimental system is introduced and the experimental results are
represented. These results are compared with the results of the simulation. As experimental
models, the vertical impact of a rigid cylinder type free link and the impact of a rigid
compound pendulum are utilized. The data attained from the repeated experiments are
averaged to minimize the efiects of abnormalities such as uctuation phenomenon which
can occur during the experiments. The resistance force coe?cients are calculated by the
part of the experimental data.
6.1 Experimental system
Figures 6.1 and 6.2 show a free link and a compound pendulum. These two mechanical
systems will be used for the experiments. For the free link shown in Fig. 6.1, the following
dimensions are given: the length L = 0.01524 m and the diameter dc = 0.00635 m. The
density of the link is 7:7?103 kg/m3. The infra red (I.R.) markers are located at 0:1L and
0:5L as shown in Fig. 6.1. For the compound pendulum shown in Fig. 6.2, the following
dimensions are given: the length Lp1 = 0.019812 m, Lp2 =0.01524 m, Lp3 = 0.23876 m,
the diameter dp1 = 0.01905 m, dp2 = 0.01524 m, and dp3 = 0.00635 m. The density of
the pendulum is 7:7 ? 103 kg/m3. Two I.R. markers are located on the rotational joint
point and the center of the Lp3. In this study, as the motion capturing system to measure
and digitize the position of the impact objects, Northern Digital Inc. (NDI) optotrak 3020
system is used. This system is composed of a position sensor, a control unit, a strober, I.R.
markers, and a PC as shown in Fig. 6.3. The system can measure the position of the markers
within the RMS accuracy of 0.1 mm and can track up to 256 markers simultaneously with
129
L
Marker1
Marker2
0.5L
0.4L
dc
Figure 6.1: Free link applied for experiments
Marker1
Marker2
Lp1
Lp2
Lp3
0.5Lp3
dp1
dp2
dp3
Figure 6.2: Compound pendulum applied for experiments
130
P C
ControlUnit Position Sensor
IR MarkersStrober
Granular metter
CRS A255Robot
Rigid bar
Figure 6.3: Motion measurement system Optotrak 3020
131
a sample up to 3500 markers/s [68]. The system does not require a calibration process.
The motions of the free link and the compound pendulum are recorded using the two
I.R. markers. The position sensor captures the position of the I.R. markers attached to
the bodies at constant sample rates and measured 3D position data. The PC is used for
operating software controlling the hardware system, saving and transform the measured
data. In the experiments, the positions were measured with 3D data cartesian coordinate
system at 500 frames/s.
There exist so many kinds of granular materials. Grains such as rice, soils including
sand, and artiflcial granules such as fertilizer, glass beads, and ball bearings. The applied
density of a granular medium for the simulations, ?g, is 2:5?103 kg/m3. The numerical value
of the density is originated from the physical property of the granular medium, \Play sand"
(Quikrete 1113-51) utilized in the impact experiments. The resistance force coe?cients
are: the dynamic frictional force coe?cient ?d = 6.5, the horizontal static resistance force
coe?cient ?h = 8, and the vertical static resistance force coe?cients ?v = 22 and ? = 1.1.
These coe?cients were determined from experimental data. The gravitational acceleration
g is applied as 9.81 m/s2. The fleld of the low speed impact with a granular medium is
related with developing multi-body kinematic chain system such as multi legged robots for
surveillance or carrying and these systems cannot avoid slow impact with diversifled outfleld
granular materials including soil and sand. From this view points, the sand is considered
more appropriate and beneflcial than glass beads or other artiflcial granular medium. The
dimension of impact test box is 0.45 m?0.32 m?0.09 m (W?L?H) and the height of the
sand in the test bed is 0.075 m.
6.2 Experimental results
6.2.1 Impact of a free link
Figures 6.4, 6.5, and 6.6 represent the experimental and the simulation results for the
impact of the free link shown in Figs. 6.1 and 4.9 for difierent impact angles and difierent
132
initial impact velocities (_qz(0) = 1.53, 2.06, and 2.47 m/s for q(0) = 0o, _qz(0) = 1.26,
1.87, and 2.33 m/s for q(0) = 32o, and _qz(0) = 1.45, 1.98, and 2.43 m/s for q(0) = 55o)
from the impact moment until the vertical velocity of the link end T, vTz, becomes zero.
Thick lines show the results of experiments and black lines represent the simulation results.
The simulation is performed based on the resistance force coe?cients determined from the
experimental results.
The penetrating depth of the link end T into the granular matter, zT, is increasing with
the initial velocity for all the cases as shown in Figs. 6.4, 6.5, and 6.6. However the stopping
time into the granular matter is decreasing when the initial velocity is increasing as most
simulation results represented in chapter 4 and 5. In these experiment and simulation, the
stopping time is deflned as the time period starting with the moment of impact and ending
when the vertical velocity of the end T into the granular matter is zero.
Even though there are difierences between the simulation and the experimental results,
the tendency of the stopping time and the penetrating depth are not changed. In this study,
a relative error [67] was calculated in order to compare the simulation and the experimental
results. The relative error of between the simulation and the experimental results, , is
deflned as
=
flfl
flflqE ?qS
qE
flfl
flfl?100; (6.1)
where qE and qS mean a position result of experiments and simulations respectively. In the
impact experiments of the link, qE is considered as the position of the end T in the vertical
direction, zT. Table 6.1 shows the difierence rate of the experiments and the simulation.
133
0.00 0.01 0.02 0.03 0.04
0.00
0.01
0.02
0.03
0.04
0.00 0.01 0.02 0.03 0.04
0.0
0.5
1.0
1.5
2.0
2.5
3.0
z T
(m)
Tim e (s)
Tim e (s)
v T
z
(m
/s)
vTz (0) = 1.53 m/ s
vTz (0) = 2.0 6 m/ s
vTz (0) = 2.47 m/ s
sim.exp.
1.
2.
3.
ts,exp .3 ts,exp .2
ts,exp .1
ts,sim .3 ts,sim .2 ts,sim .1
Figure 6.4: Experimental and simulation results of the cylinder type rigid link for q(0) = 0o
134
0.00 0.01 0.02 0.03 0.04 0.05
0.000
0.005
0.010
0.015
0.020
0.025
0.030
1.
1.
2.
2.
0.00 0.01 0.02 0.03 0.04 0.05
0
5
0
5
0
5
z T
(m)
Tim e (s)
Tim e (s)
v T
z
(m/s)
vTz (0) = 1.26 m/s
vTz (0) = 1.87 m/s
vTz (0) = 2.33 m/s
sim.exp .
1.
2.
3.
ts,exp .3 ts,exp .2 ts,exp .1
ts,sim .3 ts,sim .2 ts,sim .1
Figure 6.5: Experimental and simulation results of the cylinder type rigid link forq(0) = 32o
135
0.000 0.005 0.010 0.015 0.020 0.025
0.000
0.005
0.010
0.015
0.000 0.005 0.010 0.015 0.020 0.025
0.0
0.5
1.0
1.5
2.0
2.5
Tim e (s)
z T
(m)
Tim e (s)
v T
z
(m
/s)
vTz (0) = 1.45 m/s
vTz (0) = 1.98 m/s
vTz (0) = 2.4 3 m/s
sim.exp.
1.
2.
3.
ts,exp .3 ts,exp .2 ts,exp .1
ts,sim .3 ts,sim .2 ts,sim .1
Figure 6.6: Experimental and simulation results of the cylinder type rigid link forq(0) = 55o
136
Table 6.1: Relative error between the experimental and the simulation results of the free
link
Impact angle (o) vTz(0) (m/s) (%)
1.53 0.52
0 2.06 0.54
2.47 1.68
1.26 6.19
32 1.87 3.85
2.33 4.28
1.45 2.25
55 1.98 1.10
2.43 1.32
6.2.2 Impact of a rigid compound pendulum
Figures 6.7, 6.8, 6.9, and 6.10 represent the penetrating angle q and the angular velocity
_q of the experimental and the simulation results for the impact of the compound pendulum
shown in Figs. 6.2 4.16. The experiments and the simulations are performed for difierent
impact angles and difierent initial impact angular velocities (_q(0) = -1.75, -3.38, and -4.66
rad/s for q(0) = 22o, _q(0) = -3.31, -6.24, and -8.41 rad/s for q(0) = 31o, _q(0) = -2.66, -6.47,
and -9.06 rad/s for q(0) ? 45o, and _q(0) = -2.70, -6.54, and -9.17 rad/s for q(0) is 61.5o)
from the impact moment until the angular velocity of the pendulum, _q, becomes zero. Thick
lines show the results of experiments and black lines represent the simulation results. The
resistance force coe?cients are applied the same as those applied to the simulation of the
impact of the free link.
The penetrating angle of the pendulum is increasing with the initial angular velocity
for all the cases as shown in Figs. 6.7, 6.8, 6.9, and 6.10. These experimental results show
also that the stopping time of the pendulum is decreasing as the results of the previous
137
experimental results when the initial velocity is increasing. In these experiment and simu-
lation, the stopping time is deflned as the time period starting with the moment of impact
and ending when the velocity of the object into the granular matter is zero. The impact
experiments and the simulation results by the compound pendulum show that the tendency
of the stopping time is kept even in the case of the oblique impact with a granular medium.
The phenomenon that the impact of a rigid body stops more rapidly as the initial impact
velocity increases having been conflrmed by the experiments and the simulation results by
the impact of the free link which has initially only vertical impact velocity component is also
observed in the impact process of the compound pendulum. Figure 6.11 shows the stopping
time of the compound pendulum. The large markers are representing the simulation results
and the small markers depict the simulation results. The stopping time decreases rapidly
when the initial impact angular increases at the low angular velocity. However the decrease
of the stopping time diminishes at the low angular velocity.
As shown in Figs 6.7, 6.8, 6.9, and 6.10, there are also difierences between the simulation
and the experimental results. Table 6.2 shows the relative errors of the experiments and
the simulation. In this comparison, qE and qS mean a penetrating angular position results
of the experiments and simulations respectively. As shown in the table 6.1 and 6.2, the
penetrating end state of the simulation is almost the same as the that of the experiment
even though there are some difierences between the results of the stopping time.
138
0.00 0.02 0.04 0.06 0.08 0.10
12
14
16
18
20
22
0.00 0.02 0.04 0.06 0.08 0.10
-4
-3
-2
-1
0
Tim e (s)
q
(o
)
Tim e (s)
?q(
rad/s)
?q (0) = -1.75 rad/s1.
2.
3.
?q (0) = -3.38 rad/s
?q (0) = -4.66 rad/s
sim .exp.
ts,sim .3 ts,sim .2 ts,sim .1
ts,exp . 3 ts,exp . 2 ts,exp . 1
Figure 6.7: Experimental and simulation results of the compound pendulum for q(0) = 22o
139
0.00 0.01 0.02 0.03 0.04 0.05 0.06
22
24
26
28
30
0.00 0.01 0.02 0.03 0.04 0.05 0.06
-8
-6
-4
-2
0
Tim e (s)
q
(o
)
Tim e (s)
?q(
rad/
s)
?q (0) = -3.3 1 rad/s1.
2.
3.
?q (0) = -6.2 4 rad/s
?q (0) = -8.4 1 rad/s
sim.exp .
ts,sim .3 ts,sim .2 ts,sim .1
ts,exp . 3 ts,exp . 2 ts,exp . 1
Figure 6.8: Experimental and simulation results of the compound pendulum for q(0) = 31o
140
0.00 0.01 0.02 0.03 0.04 0.05
39
40
41
42
43
44
45
0.00 0.01 0.02 0.03 0.04 0.05
-8
-6
-4
-2
0
Tim e (s)
q
(o
)
Tim e (s)
?q(
rad/
s)
?q (0) = -2.66 rad/s1.
2.
3.
?q (0) = -6.47 rad/s
?q (0) = -9.06 rad/s
sim .exp .
ts,sim .3 ts,sim .2 ts,sim .1
ts,exp . 3 ts,exp . 2 ts,exp . 1
Figure 6.9: Experimental and simulation results of the compound pendulum for q(0) ? 45o
141
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
57
58
59
60
61
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035
-8
-6
-4
-2
0
Tim e (s)
q
(o
)
Tim e (s)
?q(rad/s)
?q (0) = -2.70 rad/s1.
2.
3.
?q (0) = -6.54 rad/s
?q (0) = -9.17 rad/s
sim.exp.
ts,sim.3 ts,sim.2 ts,sim.1
ts,exp . 3
ts,exp . 2
ts,exp . 1
Figure 6.10: Experimental and simulation results of the compound pendulum for q(0) =
61:5o
142
?q (0) (rad/s)
t s
(s)
q (0) = 22 o
q (0) = 31 o
sim.exp.
0 5 10 15 20
0.04
0.06
0.08
0.10
0.12
0.14
0.16
Figure 6.11: Stopping time results of the compound pendulum
Table 6.2: Relative error between the experimental and the simulation results of the com-
pound pendulum
Impact angle (o) q(0) (rad/s) (%)
-1.75 1.00
22 -3.38 2.58
-4.66 0.02
-3.31 1.73
31 -6.24 1.98
-8.41 1.80
-2.66 0.81
45 -6.47 1.82
-9.06 1.74
-2.70 6.66
61.5 -6.54 1.96
-9.17 1.78
143
Chapter 7
Conclusion
In this study, the dynamics of the impact with a granular medium is studied. The
most important force acting on the kinematic chain impacting with a granular medium is
the resistance force generated by the granular medium from the impact moment until the
process penetrating into the medium stops complectly. Even though a granular medium
is ubiquitous and the phenomenon of the impact with it takes place frequently whether
artiflcially or naturally in industry and nature, the studies related with the low speed impact
with the granular medium have not been su?cient for understanding the dynamics of the
impact system compared with the fleld of high speed impact. Up to present, most studies
are restricted to the impact of very low speed regime, to only vertical impact, or to the
limited shape such as a sphere. However, it becomes possible to model the resistance force
generated in a granular medium by separating the force as velocity-dependent (dynamic
frictional force) and depth-dependent (static resistance force) in this study. With this
approach, the problems concerning various shapes of kinematic chains, impact direction,
and impact velocity can be solved as shown in chapter 4, 5, and 6.
From the one D.O.F. motion of a particle such as the mathematical pendulum of a
sphere to the four D.O.F. motion of the multi impacts such as the two link kinematic
chain with two impact points, the impacts of planar rigid kinematic chains are modeled,
simulated, and some experimental results are presented. The simulations of the models
of the mathematical pendulum and the oblique impact of a sphere are performed for the
initialimpactangleof15, 45, 75? andforthevariousinitialimpactvelocities. Thesimulation
results show that the stopping time is decreasing when the initial velocity of the impact
object is increasing which result had been reported in only vertical directional impact of a
sphere. With these simulation results, it is conflrmed that even in the oblique impact with
144
a granular medium the stopping time of an impacting sphere is decreasing when the initial
impact velocity of the sphere is increasing. Even though the stopping time is decreasing, the
penetrating distance of the sphere is increasing as the initial impact velocity is increasing
for all the case of the oblique impact of the sphere including the mathematical pendulum.
As the rigid body model of a planar kinematic chain, the compound pendulum is
modeled and simulated. Even though the compound pendulum has the one D.O.F, the
motion of its rotating penetration is impeded by both horizontal and vertical resistance
force components. The simulations performed for the initial impact angles (22, 31, 45,
61.5?) and for the initial impact angular velocities of (-1.75 ? -9.17 rad/s) show the results
that the fast the impact velocity of the pendulum is, the sooner it stops and the deeper it
penetrates. This result means the impact model of a rigid body also is in uenced by the
initial impact velocity the same as those of the sphere.
The impact of a dropping link is simulated and observed in the vertical motion of the
impact penetration for the initial impact angles (0, 32, 55?) and for the initial impact linear
velocities (1.26 ? 2.47 m/s). By the observation of the simulation from the impact moment
until the velocity of the link end becomes zero, the stopping time and the penetrating depth
also can be concluded as the same as those of previous impact models.
These simulation results are conflrmed by the experimental results. Even though there
exist difierences in the stopping times between the results by the simulation and the experi-
ments, the impact experiments of the rigid compound pendulum and the rigid free link show
almost corresponding results with the simulation results. The observations of all impact
model based on the their simulation results are considered as appropriate by comparing
these results.
The simulation results of the impact of an elastic compound pendulum and an elastic
free link exhibit the consistent results observed in rigid models. According to the observation
of the simulation results, the phenomena that the decreased stopping time and the increased
penetrating distance by the increased initial impact velocity happen also in elastic kinematic
145
chain models although small difierences are observed between the elastic and the rigid
models.
The consistent simulation results are observed to be kept even in the impact of multi
link kinematic chain models such as the double pendulum and the two link kinematic chain
with two impact points for most cases. However, in some cases the stopping time of the
double pendulum models is not decreasing when the initial impact velocity is increasing.
This result is considered to happen because the withdrawing motion of the rigid body of
multi D.O.F. during the impact process might afiect the stopping time. The stopping times
are considered to increase or to decrease depending on the efiects of the withdrawing motion.
Through this study, I extended the study of the low speed impacting a granular medium
from the vertical direction impact to the oblique impact for the flrst time. For this study,
the resistance forces are analyzed and the horizontal force component is integrated to the
impact models. Based on this resistance force model, I could model difierent kinematic
chains at low speed impact with the granular medium. I contributed an experimental set-
up to validate the simulation results for kinematic links.
146
Bibliography
[1] S. Luding, \Shear ow modeling of cohesive and frictional flne powder," Powder Tech-
nology 158, 45-50, 2005.
[2] Peter Constantin, Elizabeth Grossman, and Muhittin Mungan, \Inelastic collisions of
three particles on a line as a two-dimensional billiard," Physica D 83, 409-420, 1995.
[3] Elizabeth Grossman and Muhittin Mungan, \Motion of three inelastic particles on a
ring," Physical Review E 53(6), 6435-6449, 1996.
[4] Tong Zhou and Leo P. Kadanofi, \Inelastic Collapse of Three Particles," Physical
Review E 54, 623-628, 1996.
[5] Norbert Sch?orghofer and Tong Zhou, \Inelastic collapse of rotating spheres," Physical
Review E 54(5), 5511-5515, 1996.
[6] Benjamin Painter, Meenakshi Dutt, and Robert P. Behringer, \Energy dissipation and
clustering for a cooling granular material on a substrate," Physica D 175, 43-68, 2003.
[7] Meenakshi Dutt and Robert P. Behringer, \Efiects of surface friction on a two-
dimensional granular system: Cooling bound system," Physical Review E 70,
061304(15), 2004.
[8] Piroz. Zamankhan and Mohammad Hadi Bordbar, "Complex Flow Dynamics in Dense
Granular Flows.Part I: Experimentation,\ Jornal of Applied Mechnics 73, 648-657,
2006.
[9] Piroz Zamankhan and Jun Huang, \Complex Flow Dynamics in Dense Granular
Flows.Part II: Simulations,\ Jornal of Applied Mechnics 74, 691-702, 2007.
[10] Meenakshi Dutt and Robert P. Behringer, \Efiects of surface friction on a two-
dimensional granular system: Numerical model of a granular collider experiment,"
Physical Review E 75, 021305(9), 2007.
[11] Elizabeth Grossman, Tong Zhou, and E. Ben-Naim, \Towards Granular Hydrodynam-
ics in Two-Dimensions," Physical Review E 55(4), 4200-4206, 1997.
[12] Gabriel I. Tardos, \A uid mechanics approach to slow, frictional ow of powders,\
Powder Technology 92, 61-74, 1997.
[13] Gabriel I. Tardos, M. Irfan Khan, and David G. Schaefier, \Forces on a slowly rotating,
rough cylinder in a Couette device contacting a dry, frictional powder,\ Physics of
Fluids 10, 335-341, 1998.
147
[14] Scott A. Hill and Gene F. Mazenko, \Nonlinear hydrodynamical approach to granular
materials," Physical Review E 63, 031303(13), 2001.
[15] Scott A. Hill and Gene F. Mazenko, \Granular clustering in a hydrodynamic simula-
tion," Physical Review E 67, 061302(4), 2003.
[16] Tong Zhou, \Efiects of Attractors on the Dynamics of Granular Systems," Physical
Review Letters 80(17), 3755-3758, 1998.
[17] Brian Utter and Robert P. Behringer, \Self-difiusion in dense granular shear ows,"
Physical Review E 69, 031308(12), 2004.
[18] Junfei Geng and Robert P. Behringer, \Difiusion and Mobility in a Stirred Dense
Granular Material," Physical Review Letters 93, 238002(4), 2004.
[19] G. Rousseaux, H. Caps, and J.-E. Wesfreid, \Granular size segregation in underwater
sand ripples," European Physical Journal E 13, 213-219, 2004.
[20] M. Davis, M. A. Koenders, \Oscillated densely packed granular media immersed in a
uid," Journal of Sound and Vibration 308, 526-540 2007.
[21] John J. Drozd and Colin Denniston, \Velocity uctuations in dense granular ows,"
Physical Review E 78, 041304(13), 2008.
[22] Paul Melby, Alexis Prevost, David A. Egolf, and Jefirey S. Urbach, \Depletion force
in a bidisperse granular layer," Physical Review E 76, 051307(5), 2007.
[23] James B. Knight, Christopher G. Fandrich, Chun Ning Lau, Heinrich M. Jaeger, and
Sidney R. Nagel, \Density Relaxation in a Vibrated Granular Material," Physical Re-
view E 51(5), 3957-3963, 1995
[24] E. Ben-Naim, J. B. Knight, E. R. Nowak, H. M. Jaeger, and S. R. Nagel, \Slow
relaxation in granular compaction," Physica D 123, 380-385, 1998.
[25] Edmund R. Nowak, James B. Knight, Eli Ben-Naim, Heinrich M. Jaeger, and Sidney
R. Nagel, \Density uctuations in vibrated granular materials," Physical Review E
57(2), 1971-1982, 1998.
[26] Christophe Josserand, Alexei Tkachenko, Daniel M. Mueth, and Heinrich M. Jaeger,
\Memory Efiects in Granular Material," Physical Review Letters 85, 3632-3635, 2000.
[27] Fernando X. Villarruel, Benjamin E. Lauderdale, Daniel M. Mueth, and Heinrich
M. Jaeger, \Compaction of rods: Relaxation and ordering in vibrated, anisotropic
granular material," Physical Review E 61(6), 6914-6921, 2000.
[28] Reuben Son, John A. Perez, and Greg A. Voth, \Experimental measurements of the
collapse of a two-dimensional granular gas under gravity," Physical Review E 78,
041302(7), 2008.
[29] Heinrich M. Jaeger, Sidney R. Nagel, and Robert P. Behringer, \Granular solids, liq-
uids, and gases,\ Reviews of Modern Physics 68(4), 1259-1273, 1996.
148
[30] S. N. Coppersmith, C.-h. Liu, S. Majumdar, O. Narayan, and T. A. Witten, \Model
for force uctuations in bead packs,\ Physical Review E 53, 4673-4685, 1996.
[31] Daniel M. Mueth, Heinrich M. Jaeger, and Sidney R. Nagel, \Force distribution in a
granular medium," Physical Review E 57, 3164-3169, 1998.
[32] Alexei V. Tkachenko and Thomas A. Witten, \Stress propagation through frictionless
granular material," Physical Review E 60, 687-698, 1999.
[33] M. L. Nguyen and S. N. Coppersmith, \A scalar model of inhomogeneous elastic and
granular media," Physical Review E 62, 5248-5262, 2000.
[34] Daniel L. Blair, Nathan W. Mueggenburg, Adam H. Marshall, Heinrich M. Jaeger, and
Sidney R. Nagel, \Force distributions in 3D granular assemblies: Efiects of packing
order and inter-particle friction," Physical Review E 63, 041304(9), 2001.
[35] T. S. Majmudar and Robert P. Behringer, \Contact force measurements and stress-
induced anisotropy in granular materials," Nature 425, 1079-1072, 2005.
[36] Patrick Richard, Mario Nicodem, Renaud Ddlannay, Philippe Ribiere, and Daniel
Bideau, \Slow relaxation and compaction of granular systems," Nature Materials 4,
121-128, 2005.
[37] E. Knobloch, \Euler, the historical perspective,\ Physica D 237, 1887-1893, 2008.
[38] X. W. Chen, X. L. Li, F. L. Huang, H. J. Wu, and Y. Z. Chen, \Damping function
in the penetration/perforation struck by rigid projectiles," International Journal of
Impact Engineering 35, 1314-1325, 2008.
[39] R. Albert, M. A. Pfeifer, A.-L. Barab?asi, and P. Schifier, \Slow Drag in a Granular
Medium,\ Physical Review Letters 82, 205-208, 1999.
[40] I. Albert, P. Tegzes, B. Kahng, R. Albert, J. G. Sample, M. Pfeifer, A.-L. Barabasi,
T. Vicsek, and P. Schifier, \Jamming and Fluctuations in Granular Drag,\ Physical
Review Letters 84, 5122-5125, 2000.
[41] I. Albert, P. Tegzes, R. Albert, J. G. Sample, A.-L. Barab?asi, T. Vicsek, B. Kahng, and
P. Schifier, \Stick-slip uctuations in granular drag,\ Physical Review E 64, 031307(9),
2001.
[42] I. Albert, J. G. Sample, A. J. Morss, S. Rajagopalan, A.-L. Barab?asi, and P. Schifier,
\Granular drag on a discrete object: Shape efiects on jamming,\ Physical Review E
64, 061303(4), 2001.
[43] M. B. Stone, R. Barry, D. P. Bernstein, M. D. Plec, Y. K. Tsui, and P. Schifier, \Local
jamming via penetration of a granular medium,\ Physical Review E 70, 041301(10),
2004.
[44] G. Hill, S. Yeung, and S. A. Koehler, \Scaling vertical drag force in granular media,\
Europhysics Letters 72, 137-142, 2005.
149
[45] J. S. Uehara, M. A. Ambroso, R. P. Ojha, and D. J. Durian, \Low-speed impact craters
in loose granular media," Physical Review Letters 90, 194301(4), 2003.
[46] John R. de Bruyn and Amanda M. Walsh, \Penetration of spheres into loose granular
media," Canadian Journal of Physics 82, 439-446, 2004.
[47] Massimo Pica Ciamarra, Antonio H. Lara, Andrew T. Lee, Daniel I. Goldman, Inna
Vishik, andHarryL.Swinney, \DynamicsofDragandForceDistributionsforProjectile
Impact in a Granular Medium," Physical Review Letters 92, 194301(4), 2004.
[48] L. S. Tsimring and D. Volfson, \Modeling of impact cratering in granular media,\
in: R. Garci?a Rojo, H. J. Herrmann, S. McNamara(Eds.), Powders and Grains 2005,
A.A.Balkema, Rotterdam, pp.1215-1223, 2005.
[49] M. A. Ambroso, R. D. Kamien, and D. J. Durian, \Dynamics of shallow impact cra-
tering,\ Physical Review E 72, 041305(4), 2005.
[50] M. Hou, Z. Peng, R. Liu, K. Lu, and C. K. Chan, \Dynamics of a projectile penetrating
in granular systems,\ Physical Review E 72, 062301(4), 2005.
[51] H. Katsuragi and D. J. Durian, \Unifled force law for granular impact cratering,\
Nature Physics 3, 420-423, 2007.
[52] Simon J. de Vet and John R. de Bruyn, \Shape of impact craters in granular media,"
Physical Review E 76, 041306(6), 2007.
[53] Detlef Lohse, Remco Rauhe, Raymond Bergmann, Devaraj van der Meer, \Creating a
dry variety of quicksand,\ Nature 432, 689-690, 2004.
[54] Seunghun Lee and Dan B. Marghitu, \Analysis of a rigid body obliquely impacting
granular matter,\ Nonlinear Dynamics (Accepted for publication).
[55] WilliamA.Allen, EarleB.Mayfleld, andHarveyL.Morrison, \DynamicsofaProjectile
Penetrating Sand,\ Journal of Applied Physics 28, 370-376, 1957.
[56] WilliamA.Allen, EarleB.Mayfleld, andHarveyL.Morrison, \DynamicsofaProjectile
Penetrating Sand. Part II,\ Journal of Applied Physics 28, 1331-1335, 1957.
[57] Volkhard Buchholtz and Thorsten P?oschel, \Interaction of a granular stream with an
obstacle," Granular Matter 1, 33-41, 1998.
[58] Erin C. Rericha, Chris Bizon, Mark D. Shattuck, and Harry L. Swinney, \Shocks in
Supersonic Sand," Physical Review Letters 88(1), 014302(4), 2002.
[59] C. R. Wassgren, J. A. Cordova, R. Zenit, and A. Karion, \Dilute granular ow around
an immersed cylinder,\ Physics of Fluids 15(11), 3318-3330, 2003.
[60] Junfei Geng and Robert P. Behringer, \Slow drag in two-dimensional granular media,"
Physical Review E 71, 011302(19), 2005.
150
[61] Rahul Bharadwaj, Carl Wassgren, and Roberto Zenit, \The unsteady drag force on a
cylinder immersed in a dilute granular ow," Physics of Fluids 18, 043301(7), 2006.
[62] E. L. Nelson, H. Katsuragi, and D. J. Durian, \Projectile interactions in granular
impact cratering," Physical Review Letters 101(6), 068001(4), 2008.
[63] B. Kahng, I. Albert, P. Schifier, and A.-L. Barab?asi, \Modeling relaxation and jamming
in granular media," Physical Review E 64, 051303(4), 2001.
[64] T. S. Majmudar, M. Sperl, S. Luding, and Robert P. Behringer, \Jamming Transition
in Granular Systems," Physical Review Letters 98, 058001(4), 2007.
[65] Robert P. Behringer, Karen E. Daniels, Trushant S. Majmudar, and Matthias Sperl,
\Fluctuations, correlations and transitions in granular materials: statistical mechanics
for a non-conventional system," Philosophical transactions of the royal society A 366,
493-504, 2008.
[66] Frank M. White, \Fluid Mechanics," 4th ed., WCB McGraw-Hill, 1999.
[67] DanielL.HopkinsandVictoriaL.McGu?n, \Three-Dimensional Molecular Simulation
of Electrophoretic Separations," Analytical Chemistry 70, 1066-1075, 1998.
[68] Optotrack brochure from http://www.ndigital.com/.
151
Appendices: Mathematica Programs
152
APPENDIX A
MATHEMATICAL PENDULUM
ClearAll@"Global`?"D;
Off@General::spellD;
Off@General::spell1D;
H??? Dimensions and basic data ???L
H?length of pendulum @mD: ?LL = 0.5;
H?diameter @mD: ?Lds = 0.0254;
H?density of spere @kg?m^3D: ?Lrhos = 7.7?10^3;
H?volume of spere @m^3D: ?LVs = Pi?ds^3?6;
H?mass of spere @kgD: ?Lms = Vs?rhos;
H?density of sand @kg?m^3D: ?Lrhog = 2.5?10^3;
H?gravtaional acceleration @m?s^2D: ?Lg = 9.81;
H??? Kinematics ???L
H?position vector of C: ?LrC = 8L?Sin@q@tDD, 0, L?Cos@q@tDD<;
H?Immersed depth of sphere: ?LzT = L?Cos@q@tDD ?L?Cos@q0D;
H?position vector of E: ?LrE = rC;
H?velocity vector of E: ?LvE = D@rE, tD;
H??? Dynamics ???L
H?gravity force: ?LG = 80, 0, ms?g<;
H?dynamic frictional force: Fd?L
H?coefficient: ?Letad = 6.5;
H?reference area: ?LAr = Pi?ds^2?4;
H?H direction Fd: ?LFdh = ?vE@@1DD?Ar?etad?rhog?HvE.vEL^0.5;
H?V direction Fd: ?LFdv = ?vE@@3DD?Ar?etad?rhog?HvE.vEL^0.5;
H?static resistace force: Fs?L
H?coefficient: ?Letah = 8; etav = 21; lambda = 1.1;
H?horizontal Fs: ?LFsh = ?Sign@vE@@1DDD?etah?rhog?g?ds^2?zT;
H?vertical Fs: ?LFsv = ?Sign@vE@@3DDD?etav?HzT?dsL^lambda?g?rhog?Vs;
H?resistace force:?LFR = 8Fdh+Fsh, 0, Fdv +Fsv<;
H???? Equation of motion ???L
Eq = ms?L^2?q''@tD equal1 Cross@rC, HFR +GLD@@2DD;
H???? Initial condition ???L
153
H?initial angle: ?Lq0 = 15?Pi?180;
H?initial angular velocity: ?Ldq0 = ?1;
H??? Solving E.O.M ???L
sol = NDSolve@8Eq, q@0D equal1 q0, q'@0D equal1 dq0<, q, 8t, 0, 1<,
MaxSteps ? Infinity, Method ? 8"EventLocator", "Direction" ? 1,
"Event" ? q'@tD, "EventAction" ruledelayed Throw@tend = t;, "StopIntegration"D,
"EventLocationMethod" ? "LinearInterpolation",
"Method" ? "ExplicitEuler"= 0<, 80.5, vE@@3DD ? 0<= 0<, 80.5, vE1@@3DD ? 0<= 0<, 81.1, vE1@@3DD ? 0<= 0<, 80.5, vE2@@3DD ? 0<= 0<, 81.1, vE2@@3DD ? 0<