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<