数理解析研究所講究録 1305 巻 2003 年 81-88

81
非弾性衝突の数値シミュレーション

The Simulation of the Inelastic Impact
京大人環

國仲寛人 (Hiroto Kuninaka), 早川尚男 (Hisao Hayakawa)

Graduate School of Human and Environmental Studies,
Kyoto University

1Introduction
Collisions are common phenomena in nature. For example, in the microscopic scale, atoms andmolecules in gas are colliding each other. In the macroscopic scale, we often see collision of balls in
sports such as the baseball and the billiard. In such collisions, the initial kinetic energy of material
dissipates into internal degrees of freedom like elastic vibration, sound emission, and heat. As aresult, macroscopic collisions are always inelastic.
Inelastic collisions play an important role in granular materials[l]. Characteristic behaviors
of granular material come from inelastic collisions among particles. By tilting or shaking thecontainer which contains granular material, one can see the characteristic behavior of granules
which is different from that of ordinary fluid. The Distinct Element Method (DEM) is awell-known
simulation method for the granular materials[2]. DEM contains some phenomenological parameterssuch as the Coulomb’s coefficient of friction, dashpots, and so on. Nobody can determine such theviscoelastic parameter from the first principle. However, even the determination of the simplest
parameter, the coefficient of normal restitution (COR) is not reliable.
The coefficient of normal restitution(COR) is afamiliar parameter which is introduced in text
books of the elementary physics. COR is defined by the ratio of the normal components of theinitial collision velocity and the rebound velocity as$e$$v_{i}$$v_{\mathrm{r}}$$e=-v_{\mathrm{r}}/v_{i}$,$0\leq e\leq 1$.

(1)

Historically, COR was first introduced by Newton[3]. Though many text books of elementary
physics state that COR is amaterial constant, many experiments and simulations show COR
decreases as the impact velocity increases[4, 5, 6, 7, 8, 9, 10, 11]. On the other hand, Louge andAdams reported in their recent paper that COR can exceeds unity in the situation of the obliqueimpact which is contrary to the assumption $e\leq 1[12]$ . This topic is interesting and worthy of moredetailed study.
In addition, the coefficient of tangential restitution is also well-known parameter to describe
the rotational motion of material. is defined as$e$$\beta$$\beta$$\beta=-\frac{v_{\acute{t}}}{v_{t}}$,(2)where and are the tangential components of the velocity of the contact point before and after
collision. is known to be dependent on the incident angle of impact. However, the mechanism of
this dependency is not unclear$v_{\acute{t}}$$v_{t}$$\beta$82
Orn lesealcll is to understand the lnechanisrn of the coefficient of tangential restitution. We
study tlic relation between the coefficient of tangential restitution and the angle of incidence inoblique collision in this paper. The organization of this ) ) is as follows. In tlle next section, wewill rcvie the definition of the coefficient of restitution and the coefficient of tangential restitution.
In section 3we introduce our numerical model and setup of the simulation. Section 4is the mainpart of this paper where we summarize tlle results of our simulation and, explain the numericalresults bv the theory. Section 5is the conclusion of this paper.
$1$$\mathrm{a}1$$\mathrm{e}1$$\backslash \mathrm{v}$2Introduction of e and$\beta$

To characterize inelastic collision, Walton introduced three parameters[13]. The three parameters
are the coefficient of normal restitution , the coefficient of Coulomb’s friction , and the rnaximum
value of the coefficient of tangential restitution . Experiments have supported that his characterization adequately capture the essence of binary collision of spheres or collision of asphere on
aflat plate[14, 15, 16, 17]. Now, let us define the coefficient of restitution and the coefficient oftangential restitution in the 2-dimensional situation. Figure 1is the schematic figure that adisk$\mathrm{e}$$\mu$$\beta_{0}$$\mathrm{e}$$\beta$$\mathrm{v}_{\mathrm{c}}$$\mathrm{t}$aI$\mathrm{t}\mathrm{b}$IFigure 1: The schematic figure of acollision of sphere with awall.

is colliding with astationary wall with initial velocity of its center of mass, Vi. The relative velocity
at the contact point after collision, thus, becomes$\mathrm{v}_{\mathrm{c}}’=\mathrm{v}_{\mathrm{i}}-R\mathrm{n}\mathrm{x}$$\omega’$(3),is the unit vector in the normal direction to the wall, andis the radius of the disk,
where’is the angular velocity. The prime denotes post-colliding quantities. The coefficient of normal
restitution is defined as$R$$\mathrm{n}$$\omega$$\mathrm{e}$(4).

$\mathrm{v}_{\acute{\mathrm{c}}}\cdot \mathrm{n}=-e\mathrm{v}_{\mathrm{c}}\cdot \mathrm{n}$.
Conventionally, this parameter is assumed to be
The coefficient of tangential restitution is defined as$0\leq \mathrm{e}\leq 1$$\beta$$\mathrm{v}_{\acute{\mathrm{c}}}\cdot \mathrm{t}=-\beta \mathrm{v}_{\mathrm{c}}\cdot \mathrm{t}$,(5)and are the post-collisional velocity at the contact point after collision and the unit
wheretangential vector, respectively. It is believed that is afunction of the angle of incidence , withpossible values lying in the range between -1 and 1 $[13, 14]$ . The incident angle is defined as, respectively.
, whereandand are$\mathrm{v}_{\acute{\mathrm{c}}}$$\mathrm{t}$$\beta$$\gamma$$\gamma$$\gamma=\arctan(\mathrm{v}_{t}/\mathrm{v}_{n})$$\mathrm{t}$$\mathrm{v}_{n}$$\mathrm{v}_{t}$$\mathrm{v}_{n}=\mathrm{v}_{\mathrm{c}}\cdot$$\mathrm{n}$$\mathrm{v}_{t}=\mathrm{v}_{\mathrm{c}}\cdot$

83
For the oblique collision, the coefficient of tangential restitution is more important than .
From the conservation laws of momentum and angular momentum and Coulomb’s friction on tlte
, Walton[13] derives
surfaces of two identical igid )$e$$\beta$$1\mathrm{l}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{s}$$\mathrm{r}$$\mathrm{s}\mathrm{l}$$-1- \mu(1+e)\cot\gamma(1+\frac{mR^{\underline{)}}}{I})$$(\gamma\geq\gamma_{0})$$\beta\simeq\{$$(\gamma\leq\gamma_{0})$$\beta_{0}$,(6)whereis the critical angle, and $m$ , , and I are mass, radius and moment of inertia of spheres
respectively. Labous, Rosato, and Dave performed the experiment of binary collision of nylon
spheres and showed the consistency of their results to the Walton’s relation[14]. Furthermore, it
has become clear that Many experimental results are consistent with the relation so that Walton’smodel is accepted as reasonable[15, 16, 17]. Meanwhile, Maw, Barber, and Fawcett extended theHertz theory of impact and established the theory of the oblique impact to be consistent with
their experimental results[18]. In contrast to Walton’s assumption, they demonstrated the need toconsider normal and tangential compliance over the contact area.
$R$$\gamma_{0}$3Our Models
Here, let us introduce three lattice models. Each model consists of an elastic disk and an elasticwall. The main results of this paper are those of random lattice model(Fig. 2). Both the disk andFigure 2: The elastic disk and wall consisted of random lattice system.
the wall are composed of randomly distributed 800 mass points. All mass points are bound withnonlinear springs using the Delaunay triangulation algorithm[19]. The spring interaction between
connected mass points is described as$V(x)= \frac{1}{2}k_{a}x^{2}+\frac{1}{4}k_{b}x^{4}$,

(7)

84

,

$\mathrm{F}_{-}$

Figure 3: Interaction between surface particles of the disk and the wall.
whereis astretch from the natural length of spring, and and are the spring constants. We
use atypical ratio of toas $k_{b}/k_{a}=10^{-3}$ . The width of the wall is 4times as long as thediameter of the disk. The height of the wall is same as the diameter of the disk. Two sides of thewall are fixed.
The interaction between the disk and the wall during acollision is introduced as follows. Figure
3is the schematic figure of the interaction of surface mass points of the disk and the wall. When thedistance between the edge of the disk and the surface of the wall is less than the cutoff length(weset it equal to the length of the linear spring), the surface particles of the disk feel the repulsive
force, $\mathrm{F}(l)=aV_{0}\exp(-al)\mathrm{n}$, where is $300/R$ ,is $amc^{2}R/2$ ,is the mass of the particle,
is the radius of the disk, $c=\sqrt{E}/\rho$ ,is Young’s modulus, and is the density, is the normal
unit vector to the surface. The reaction forces applied to the two points of the surface of thewall (point 1and 2) are decided by the balance of the torques as Fi $(/)=-F(l)\mathrm{n}/(1+/1//2)$ and, where $l_{i}(i=1,2)$ is the distance between the point and the point(see Fig. 3).
In this model, roughness of the surfaces is important mechanism to make the disk rotate aftercollision. How to make roughness is as follows. At first, we generate normal random numbers whose
average value is 0and then make the initial position of particles on surface of both the disk andthe wall deviate with them. We choose the value of dispersion as $\delta=3\cross 10^{-2}R$ , where is theradius of the disk.
As for random lattice model, we cannot determine Poisson’s ratio theoretically. When wedetermine Poisson’s ratio of this model, we introduce the viscous damping term in (7). By stretching
the strip of random lattice and measuring its width and height, we can obtain Poisson’s ratio.
For comparison, we make other two lattice models: triangular lattice and square lattice disk(Fig.4).
The triangular lattice disk is made by replacing the internal structure of the random lattice disk
with the triangular lattice. The surface of the triangular lattice disk is same as that of the random
lattice disk. Poisson’s ratio of the triangular lattice can be calculated theoretically as 1/3[20]. Thesquare lattice disk is made by replacing the internal structure of the random lattice disk with thesquare lattice. We introduce two spring constants: $k_{a}=k_{1}$ for nearest neighbor interaction and$k_{a}=k_{2}$ for next-nearest neighbor interaction. Poisson’s
ratio of the square lattice is expressed as$x$$k_{a}$$k_{b}$$k_{b}$$k_{a}$

$l$

$a$

$R$

$m$

$V_{0}$$E$$\mathrm{n}$$\rho$$\mathrm{F}_{2}(l)=-F(l)\mathrm{n}/(1+l_{2}/l_{1})$

$i$

$p$

$R$

$\delta$$\nu=\frac{k_{2}^{2}+(k_{1}^{2}-4k_{2}^{2})n_{x}^{2}n_{y}^{2}}{k_{2}(k_{1}+k_{2})+(k_{1}^{2}-4k_{2}^{2})n_{x}^{2}n_{y}^{2}}$.

We scale the equation of motion for each particle using the radius of the disk

(8)
$R$

as the scale of

85

(b)

(a)

Figure 4: The schematic figures of (a) triangular lattice disk and (b) square lattice disk.
length and the velocity of elastic wave $c=\sqrt{E}/\rho$ as the scaling unit of velocity. As the numerical
scheme of the integration, we use the fourth order symplectic numerical method with the time step$\Delta t\simeq 10^{-3}R/c$ .

4Results and Discussions
In this section, we carry out the simulation of the oblique impact. The angle of incidence is ranged
from 5.7’ to 80.5’ while the normal component of velocity is fixed as 0. $1c$ . The disk has no internal
vibration and rotation before collision. In order to eliminate the effect of the initial configurationof mass points, we prepare 100 samples of disk as the initial condition by using 100 sets of random
numbers and average data of all samples.

Figure 5: The relation between$\cot\gamma$

and . Figure 6: The relation between$\beta$$\cot\gamma$

and .
$e$

Figure 5shows the relation between the cotangent of the angle of incidence 7and the coefficien$\mathrm{t}$86
of tangential restitution . In this figure, cross points are the result of the 1
lattice diskeq.(6),$e=0.8$,,andand broken lines arewhereandare fitting parameters. This resultshows thattakes the value nearly 0.56 andtakes tlle value nearly 0.18. From this estimation,
we can see that this model can reproduce the tendency of the experimental results of the obliquecollision qualitatively $[14, 15]$ . In contrast, plus points are the results of the triangular lattice disk.
In this lnodel, takes negative values in all range of the angle of incidence. This means that thetriangular lattice model is easy to slip on the surface.
Figure 6shows the relation between the cotangent of the angle of incidence and COR . Although it is expected that COR takes the constant value because the normal velocity of the disk isset to the fixed value, 0. , COR depends on the angle of incidence. In particular, in the region of, COR decreases as
small value ofdecreases. At present, we cannot explain this tendencyof normal COR.
Here, we compare our result with the theory of Maw, et a1.[18]. According to their theory, all
the region of the angle of incidence can be divided into three regimes. For each regime, can beexpressed as$\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\ln$$\beta$$\backslash \backslash \prime \mathrm{a}11$$\mu$$\beta_{0}$$\beta_{0}$$\mu_{0}$$\beta$$e$$1\mathrm{c}$$\cot\gamma$$\cot\gamma$$\beta$(i) $1/\mu\eta^{2}<\cot\gamma$ :,$\beta=\cos\omega t_{1}(\gamma)+\mu\alpha e[1+\cos(\frac{\Omega t_{1}(\gamma)}{e}+\frac{\pi}{2}(1-e^{-1}))]\cot\gamma$(9)(ii) $\mu(1+e)/\alpha<\cot\gamma<1/\mu\eta^{2}$ :$\beta=\cos\omega t_{3}(\gamma)+\mu\alpha[1+e-\frac{p(t_{3}(\gamma))}{p(t_{\mathrm{c}})}]\cot\gamma$

,

(10)

(iii) $\cot\gamma<\mu(1+e)/\alpha$ :$\beta=1-\mu\alpha(1+e)\cot\gamma$

,

(11)

where
is the coefficient of friction, is the constant dependent on Poisson’s ratio, $\alpha=3.02$
which is aconstant dependent on the shape of material, $\Omega=\pi/2t_{c}$ , is aduration of acollision,,is the transition time from stick motion to slip motion,is the transition
time from slip motion to stick motion, and $p(t)$ is impulse. This theory was confirmed to beconsistent with experimental data[15, 16, 17, 18].
We compare the result of simulation of the oblique impact using the random lattice model withthe theoretical curve(Fig. 7). Here we used $\eta=1.015$ , which corresponds to Poisson’s ratio 0.058,$e=0.8$ as afixed value, and $\mu=0.3$ as afitting parameter. It is found that the result of random
lattice model is consistent with the theory.
On the other hand, as for the result of figure 5, we focus our attention to the difference ofPoisson’s ratio between the random disk and the triangular disk. By changing the value of spring
constants of square lattice disk and controlling Poisson’s ratio, we investigate the dependency ofon Poisson’s ratio. Figure 8is the result when $\nu=0.1$ while figure 9is the result when $\nu=0.3$ .
We cannot see the difference of the values of . From these results, Poisson’s ratio seems not toaffect the value of .
$\mu$$\eta$$t_{/c}$$\omega$$=(\pi/2\eta t_{c})\sqrt{\alpha}$$t_{1}(\gamma)$$t_{3}(\gamma)$$\beta_{0}$$\beta_{0}$$\beta_{0}$87$\beta$$\mathrm{c}\mathrm{o}q$and . Cross points are the numerical results of the random
Figure 7: The relation between
lattice model. Solid line is the theoretical curve.
$\cot\gamma$$\beta$AFigure 8: The relation between
when $\nu=0.1$$\cot\gamma$

and$\beta$. Figure 9: The relation between$\cot\gamma$

and$\beta$.

when $\nu=0.3$

5Conclusion
In this paper, we demonstrate the 2-dimensional simulation of the oblique impact and obtain resultsas follows.
(i) Our random lattice model exhibits the same tendency as experimental data qualitatively. In
addition, the model is consistent with Maw’s theory of the oblique impact.
(ii) There seems to be no relation between Poisson’s ratio of material and the value of$\beta_{0}$.

References[1] See, for example, L. P. Kadanoff: Rev. Mod. Phys. 71 (1999) 435; P. G. de Gennes: Rev. Mod.
Phys. 71 (1999) S374 and references therein88[2] P. A. Cundall, 0. D. L. Strack: Geotechnique 29 (1979) 47.
[3] I. Newton: Philoshophiae naturalis Principia mathematica (W. Dawason and Sons, London,1962). The original one has been published in 1687.

[4] W. Goldsmith: Impa ct: The Theo ry and Physical Behavior of Colliding Solids (Edward Arnold
Publ., London, 1960).
[5] K. L. Johnson: Contact Mechanics (Cambridge University Press, Cambridge, 1985).

[6] W. J. Stronge: Impact Mechanics (Cambridge Univ. Press, 2000)[7] R. Sondergaard, K. Chaney, and C. E. Brennen: Transaction of the ASME, Journal of Applied
Mechanics 57 (1990) 694.

[8] F. G. Bridges, A. Hatzes, and D.N.C. Lin: Nature 309 (1984) 333.
[9] K. D. Supulver, F. G. Bridges, and D. N. C. Lin: ICARUS 113 (1995) 188[10 G. Kuwabara and K. Kono: Jpn. J. Appl. Phys. 26 (1987) 1230.
[11 H. Hayakawa and H. Kuninaka: Chem. Eng. Sci. 57 (2002) 239.
[12 M. Y. Louge and M. E. Adams: Phys. Rev. E65 (2002) 021303.

[13 O. R. Walton and R. L. Braun: J. Rheol. 30 (1986) 949.

[14 L. Labous, A. D. Rosato, and R. N. Dave: Phys. Rev. E56 (1997) 5717.

[15 S. F. Foerster, M. Y. Louge, H. Chang, and K. Allia: Phys. Fluids 6(1994) 1108.
[16 A. Lorentz, C. Tuozzolo, and M. Y. Louge: Exp. Mech. 37 (1997) 292.

[17 D. A. Gorham, A. H. Kharaz: Powder Technology 112 (2000) 193.
[18 N. Maw, J. R. Barber, and J. N. Fawcett: Wear 38 (1976) 101.; N. Maw, J. R. Barber, and J.
N. Fawcett: ASME J. Lub. Tech 103 (1981) 74.
[19] K. Sugihara: Data Structure and Algorithms (Kyoritsu, Japan, 2001).

[20] W. G. Hoover: Computational Statistical Mechanics (Elsevier Science Publishers B. V., Amsterdam, 1991