The theoretical model of a kinematic chain impacting granular matter is studied. The force of the granular medium acting on the chain is a linear superposition of a static (depth-dependent) resistance force and a dynamic (velocity-dependent) frictional force. This resistance force is opposed to the direction of the velocity of the immersed chain. We present two methods (one using EventLocator and the other using FixedStep) for the problem. As examples, a single and a double pendulum are simulated using different initial impact velocity conditions. We analyze how rapidly the kinematic chain impacting the granular medium slows upon collision. For the analyzed cases the kinematic chain under high impact force (higher initial velocity) comes to rest faster in the granular matter than the same body under low impact force (lower initial velocity).

Introduction

The physical behavior of granular matter has its own specific properties that are unlike solids, liquids, or gases. Impact with a granular medium is very interesting because granular materials have characteristics similar to a solid, yet flow like a fluid. The study of impact into granular matter is, surprisingly, in its early stages.

The focus of recent work has been to develop a force law for granular impacts and to find a mathematical formula to measure the impact force of objects dropped into granular matter [1]. Tsimring and Volfson [2] studied the penetration of large projectiles into dry granular media. For the resistance-force model, they proposed a drag force depending on the velocity and a friction force depending on the depth. Depth-dependent friction-force models have been developed for horizontal motion in [3, 4, 5, 6] and for vertical motion in [7, 8]. Katsuragi and Durian [1] sparked new interest in the field of impact with granular matter. They applied the resistance-force model proposed in [2] and verified the motion of a rigid sphere with a line-scan digital CCD camera. They analyzed an interesting phenomenon: how rapidly a sphere impacting a granular medium slows upon collision. Analysis shows that the greater the speed at which the spheres impact the medium, the sooner they will come to rest for the vertical impact.

In our study we focus on modeling and simulating a single and a double pendulum impacting a granular medium using the resistance-force model proposed by Tsimring and Volfson [2] and verified by Katsuragi and Durian [1]. For the single pendulum we apply NDSolve with variable step size and the simulation is stopped when the velocity changes sign. This event is captured with the command EventLocator. For the double pendulum the equations of motion are solved using NDSolve with the FixedStep method. In this way we capture the physics of the impact process and we can calculate the final stopping time for the kinematic chain.

Resistance Forces during Impact

For impact with penetration, the most important interaction between a rigid body and the medium is the resistance force from the moment of impact until the body stops. The resistance forces are composed of a drag force proportional to the square of the velocity and a force related to the plasticity of the medium [9, 10]. For impact of a moving rigid body into a granular medium, recent research [1, 2] shows that the total resistance force acting on the body is the sum of a static resistance force characterized by depth-dependent friction force and a dynamic frictional force characterized by velocity-dependent drag force, that is, , where is the static resistance force and is the dynamic frictional force.

The dynamic frictional force , a velocity-dependent force, is generated by the movement of the rigid body into the granular matter from the moment of impact until rest. Even though the state of the granular matter is not a fluid but a solid, the rigid body moves through the granular matter as if it were partly like a fluid. From this fluid-like behavior, the resistance force is assumed to be like a drag force impeding the motion of the rigid body. Recent results [1, 2] show that , where is a drag coefficient determined from experimental results and is the reference area. The factor is negative because the dynamic frictional force acts in the opposite direction of the velocity . The dynamic frictional force can be applied separately, with the horizontal and vertical dynamic frictional forces and given by

(1)

(2)

where , are the horizontal and vertical velocities of a rigid body and , are the reference areas to the horizontal and the vertical directions.

The static resistance force is an internal resisting force and appears when an external force is applied. The horizontal static resistance force is defined as the internal impeding force acting on the horizontal direction. Recent research results [3, 11] show that it is known to be proportional to the normal component of the pressure acting at the contacting point, which increases linearly with the depth , the granular density , and the gravitational acceleration . The horizontal static resistance force of a cylinder including the state of immersion, at any slope, can be generalized as

(3)

where is the depth of the immersed cylinder tip.

The vertical static resistance force is defined as the internal impeding resistance acting on the vertical axis. This force is a nonlinear function of the immersion depth. Hill, Yeung, and Koehler [8] suggest an empirical equation with coefficients calculated from experimental data as

(4)

where is the volume of the rigid body and is the lateral dimension. The coefficients and depend on the shape of the rigid body, the properties of the granular matter and the rigid body, the shape of the medium container, and the direction of motion, such as plunging and withdrawing.

Application

Equations (1) to (4), representing the resistance forces, contain the sign function of the velocity of the application point of the resistance forces. It is not possible to solve the ordinary differential equations of motion that contain the sign function with the variable step method of NDSolve, due to the discontinuity at the origin. We use the EventLocator method for NDSolve with a variable step to solve the ODE for the single-pendulum impact model. We use the FixedStep method to solve the ODE for the double-pendulum impact model. This correctly captures the change in the sign of the velocity.

Single-Pendulum Impact Model

A single-pendulum impact model is presented in Figure 1. The planar motion of the pendulum on impact can be described in a fixed Cartesian coordinate system. The axis is horizontal, with the positive sense directed to the right; the axis is vertical, with the positive sense directed downward; and the axis is perpendicular to the plane of motion. The angle between the axis and the pendulum is denoted by , the generalized coordinate describing this system.

Figure 1. A single-pendulum impact model.

The pendulum has mass , length , and diameter . The pendulum impacts a granular medium of density at an initial impact angle . Gravitational acceleration is . Here are the numerical values for the input data.

The gravitational force acts at the center of mass of the pendulum. This is the position vector from the origin to C.

The resistance force of the granular media is assumed to act at the point , the center of the immersed part of the pendulum. The immersed depth is as shown in Figure 1. The length is the distance from the origin to .

These are the position vector and the velocity of the point .

This is the gravitational force acting at .

We use NDSolve to solve the differential equations of motion. The expressions for the resistant forces given by equations (1) to (4) contain the signs of the vertical and horizontal velocities of the application point E. We want to use a variable step size for the integration for this model, so we cannot use the Mathematica function Sign. (The function Sign can be used with NDSolve with the FixedStep method; see the double-pendulum example.) Therefore for the case of the single pendulum we introduce the piecewise-constant variables and , defined in the following lists.

The list data1 is used for plunging () and the list data2 is used for withdrawing (). For the two lists we consider continuous sliding on the horizontal direction. The general coefficients for the resistant force [1, 3, 8] are and . For plunging motion , and for withdrawing motion , . When the pendulum penetrates the granular medium with plunging motion (), the values of and are –1 and 1. When the pendulum withdraws from the granular medium (), the values of and are –1 and –1.

This defines the reference area and the horizontal dynamic frictional force .

Here are the reference area and the vertical dynamic frictional force .

Here is the horizontal static resistance force.

Here are the immersed volume U of the pendulum and the vertical dynamic frictional force .

This defines the resistance force.

The equation of motion of the pendulum is

(5)

Here is the mass moment inertia about the origin . The ODE, equation (5), is solved in single. The input variable dq is the initial impact angular velocity of the pendulum. The initial start time t0, the initial impact angle q0, and the initial impact angular velocity dq0 are defined as initial values. The interpolating function sol, which is the solution of the ODE, is defined as a local symbol. To save the solution data for the angle q and for the angular velocity , we introduce initially the empty vectors qresults and dqresults. The simulation is performed until the angular velocity becomes zero (considering a numerical error) using the While statement. At the first loop the ODE uses the coefficients list data1 for the plunging motion. NDSolve solves this ODE until the directional velocity of the point becomes zero and the event is identified using the EventLocator method. This event occurs when is zero or when changes sign at . We do not consider the velocity, , of the point because the sign of the velocity does not change until the angular velocity becomes zero. Once NDSolve stops solving the ODE, the stopping time is saved as tend. The last values of the angle q and the angular velocity of the stopping time, tend, are saved as the initial values q0 and dq0 for the next loop. The Table saves the data of the angle q and the angular velocity using qresult and dqresult. After saving the data, the new initial start time t0 of the next loop is defined as the stopping time tend and the motion index increases from the plunging motion to the withdrawing motion. After the end of the first loop, the While statement checks the last value of the angular velocity and decides whether to do one more loop of calculation or to end the module.

The simulation of the model is performed for an initial impact angle and for three initial impact angular velocities (, and rad/s). The angle and the angular velocity are plotted with the blue line for rad/s, the red dashed line for rad/s, and the black dot-dashed line for rad/s.

The simulation results of the model for the initial impact angle show that when the initial impact velocity of the pendulum increases, the penetrating angle increases. However, we observed the stopping time, defined as duration from the moment of impact until the angular velocity becomes zero, decreases: for rad/s the stopping time is 0.430259 s, for rad/s the stopping time is 0.474437 s, and for rad/s the stopping time is 0.553072 s.

Double-Pendulum Impact Model

A double-pendulum impact model is presented in Figure 2. The planar motion of the double-pendulum impact can also be described in a fixed Cartesian coordinate system. The angles between the axis and the links are denoted by and ; these angles are the generalized coordinates describing this system.

Figure 2. A double-pendulum impact model.

The first link of the pendulum has mass and length . The second link has mass , length , and diameter . The density of the granular medium is and the gravitational acceleration is . The pendulum impacts a granular medium at the initial impact angles for the first link and for the second link. Here are the numerical values of the input data.

This is the position vector from the origin to the mass center of the first link.

This is the position vector from the origin to the point .

This is the position vector from the origin to the mass center of the second link.

Here is the acceleration vector of the position vector.

We calculate the immersed depth and the length .

Here are the position vector and the velocity of the force application point .

These are the gravitational forces acting at and .

The function Sign can be used with NDSolve with the FixedStep method as mentioned previously. Therefore, we define the horizontal dynamic frictional force.

Here is the vertical dynamic frictional force.

Here is the horizontal static resistance force.

We express the vertical resistance force coefficients, and for plunging motion, and for withdrawing motion.

This is the vertical static resistance force.

This is the resistance force.

The ODE for the first link, where is the mass moment inertia about the origin , is

(6)

The ODE of the second link, where is the mass moment inertia about , is

(7)

The ODEs are solved in double. The input variables dq1 and dq2 are the initial impact angular velocities of the pendulum. The initial impact angles, q10, q20, and the initial impact angular velocities, dq10, dq20, are defined as initial values. The solution is obtained using the ExplicitRungeKutta method with FixedStep with step size s and stops automatically using the EventLocator method when and simultaneously become zero (up to numerical error).

The simulation of the model is performed for initial impact angles ( and ) and for three initial impact angular velocities ( rad/s). The angle and the angular velocity are plotted in blue for rad/s, red dashed for rad/s, and black dot-dashed for rad/s.

The simulation results of the model for the initial impact angles ( and ) show that the stopping time decreases when the initial impact velocity of the pendulum increases: for rad/s the stopping time is 0.389631 s, for rad/s the stopping time is 0.457272 s, and for rad/s the stopping time is 0.53981 s.

Conclusion

This article considers the impact of a kinematic chain into granular matter. The solution techniques are based on the resistance forces as a sum of a static (depth-dependent) resistance force and a dynamic (velocity-dependent) frictional force. We apply the NDSolve command to solve the equations of motion. For the single pendulum, NDSolve is used with EventLocator to identify when the velocity changes sign. The simulation is then stopped and a new simulation is performed using as initial conditions the end values of the previous simulation. For a double pendulum, NDSolve is used with the FixedStep method and the simulation is performed over the whole interval of time. We observe that the numerical solution leads to a seemingly paradoxical result: as the initial velocity increases, the stopping time decreases.

Acknowledgments

For this work we used Mathematica 7.0.1 and a Intel Core Duo computer running Mac OS X with a 2.0 GHz CPU and 2 GB of RAM.

References

[1]	H. Katsuragi and D. J. Durian, “Unified Force Law for Granular Impact Cratering,” Nature Physics, 3, 2007 pp. 420-423. www.nature.com/nphys/journal/v3/n6/full/nphys583.html.
[2]	L. S. Tsimring and D. Volfson, “Modeling of Impact Cratering in Granular Media,” in Powders and Grains 2005, Proceedings of International Conference on Powders & Grains 2005 (R. García-Rojo, H. J. Herrmann, and S. McNamara, eds.), London: Taylor & Francis, 2005 pp. 1215-1218.
[3]	R. Albert, M. A. Pfeifer, A.-L. Barabási, and P. Schiffer, “Slow Drag in a Granular Medium,” Physical Review Letters, 82(1), 1999 pp. 205-208. www.prl.aps.org/abstract/PRL/v82/i1/p205_1.
[4]	I. Albert, P. Tegzes, B. Kahng, R. Albert, J. G. Sample, M. Pfeifer, A.-L. Barabási, T. Vicsek, and P. Schiffer, “Jamming and Fluctuations in Granular Drag,” Physical Review Letters, 84, 2000 pp. 5122-5125. arxiv.org/abs/cond-mat/9912336.
[5]	I. Albert, P. Tegzes, R. Albert, J. G. Sample, A.-L. Barabási, T. Vicsek, B. Kahng, and P. Schiffer, “Stick-Slip Fluctuations in Granular Drag,” Physical Review E, 64, 2001 031307(9). arxiv.org/abs/cond-mat/0102523.
[6]	I. Albert, J. G. Sample, A. J. Morss, S. Rajagopalan, A.-L. Barabási, and P. Schiffer, “Granular Drag on a Discrete Object: Shape Effects on Jamming,” Physical Review E, 64, 2001 061303(4). arxiv.org/abs/cond-mat/0107392.
[7]	M. B. Stone, R. Barry, D. P. Bernstein, M. D. Pelc, Y. K. Tsui, and P. Schiffer, “Studies of Local Jamming via Penetration of a Granular Medium,” Physical Review E, 70, 2004 041301(10). arxiv.org/abs/cond-mat/0404587.
[8]	G. Hill, S. Yeung, and S. A. Koehler, “Scaling Vertical Drag Forces in Granular Media,” Europhysics Letters, 72(1), 2005 pp. 137-142. www.iopscience.iop.org/0295-5075/72/1/137.
[9]	A. L. Yarin, M. B. Rubin, and I. V. Roisman, “Penetration of a Rigid Projectile into an Elastic-Plastic Target of Finite Thickness,” International Journal of Impact Engineering, 16(5), 1995 pp. 801-831. www.sciencedirect.com/science/article/pii/0734743X95000197.
[10]	G. Yossifon, A. L. Yarin, and M. B. Rubin, “Penetration of a Rigid Projectile into a Multi-Layered Target: Theory and Numerical Computations,” International Journal of Engineering Science, 40(12), 2002 pp. 1381-1401. doi:10.1016/S0020-7225(02)00013-7.
[11]	S. N. Coppersmith, C.-H. Liu, S. Majumdar, O. Narayan, and T. A. Witten, “A Model for Force Fluctuations in Bead Packs,” Physical Review E, 53(5), 1996 pp. 4673-4685. www.arxiv.org/abs/cond-mat/9511105.
S. Lee and D. B. Marghitu, “Impact of a Planar Kinematic Chain with Granular Matter,” The Mathematica Journal, 2011. dx.doi.org/doi:10.3888/tmj.13-2.

About the Authors

Seunghun Lee was a Ph.D. student in the mechanical engineering department at Auburn University, studying theoretical modeling of robotic systems. He is now serving in the Korean army.

Dan B. Marghitu is a professor of mechanical engineering at Auburn University. His research areas are impact dynamics, mechanisms, robots, and nonlinear dynamics.

Mechanical Engineering Dept.
Auburn University
270 Ross Hall, Auburn, AL 36849
lsh10069@gmail.com
marghitu@auburn.edu