Title 巨視的物体の非弾性衝突 Author(s) - Kyoto University Research

Title
Author(s)
Citation
Issue Date
URL
巨視的物体の非弾性衝突
國仲, 寛人
物性研究 (2008), 90(5): 685-720
2008-08-20
http://hdl.handle.net/2433/142659
Right
Type
Textversion
Departmental Bulletin Paper
publisher
Kyoto University
物性研究
90-5 (
2
0
0
8
8
)
巨視的物体の非弾性衝突 1
中央大学理工学部物理学科
関 伸 寛 人2
(2008年 6月 17日受玉里〉
1 はじめに
我々の自に見える巨視的な物体の衝突現象は吉くから物理学や工学等の分野に登場し、紀元前
から現代に至るまで、議々なスケールで研究の対象となっている [
1，
2，
3，
4
]
0 このように密突の問
題が長い間研究の対象とされ、議論され続けているのは、我々の身の周りに衝突に絡んだ現象が
溢れているからではないだろうか。例えば釘を打つ、模を打ち込むといった行為はハンマーの衝
突による撃力を利用したものであるし、野球やゴルフ等の球技はパットやクラブとボールとの衝
突を利用したものである。更に我々の毘の及ばない範囲においても、気体や譲体等の流れほそれ
らを構成する原子・分子関の衝突が素過程となって引き起こされ、宇宙に無数に漂う小惑星は天
5
]
0
体同士の衝突によって形成されたものと考えられている [
また、近年ソフトマタ一物理の研究対象として盛んに研究されているものに、粉体 [
6，7，
8，9，
1
0，1
1，1
2
]がある。砂やガラスビーズ等の数 m m程衰のサイズの粒子集毘を容器に入れて加振し
たり、傾けて流したり、垂直に立てた細いパイプ中に流すと液体等の他の流体とは異なった振る
舞いを示すことが多くの実験やシミュレーションから知られている [
1
3，14，15，1
6
]
0 そのような
粉体の複雑な振る舞いも、局所的;こ見れば粉体粒子同士の衝突が連続して起こることで生み出さ
れたものであると言える。
巨視的な物体同士の衝突は一般に f
非弾性衝突j と呼ばれるが、それは衝突した 2物誌の運動
エネルギーの和が衝突前の笛と等しくならないということを意味している。衝突前後の運動エネ
ルギーの差分は物体の変形や破壊、熱といった様々な形態で散逸するが、どのような散逸の形態
が支配的になるかは物質の種類のみならず、衝突前の運動状態によっても大きく異なる。例えば
同じ大きさを持つ 2つの金屠球同士を正面衝突させると、衝突速度が低速の時には球の内部摩擦
が重要な散逸機講となるが、高速の時には塑性変形や破壊が主要な散逸機構となる。ではその中
間の速度領域においてはどのような散逸機構が支配的になるかというと、これは様々な散逸機構
が混ざり合っていると考えられるので、どれが支配的かを特定することは難しい。この事が通常
巨視的な物体の衝突後の運動状態の理論的予測を圏難にする。
初等物理のテキストにおいては、こういった衝突持問中のエネルギー散逸の形態や過程を一切
はねかえり保数j というパラメータでエネルギー散逸の度合いを表現する。傍として、
考えずに、 f
本稿は、編集部の方から持にお願いして執筆していただいた記事である。
E
E
-mai
J
:k
[email protected]
o
-u
.
a
c
.
j
p
2
Fhd
o
δ
FO
題伸寛人
L
v
!
e
l
'
(
a
)
(
b
)
図 1
:ボールが地面に斜め禽突する様子。 (
a
)は窪i
突前、 (
b
)は笥突が終わってボーんと壁が離れ
る瞬間の様子を表している。
図 1~こボールが地面に斜めに笥突する様子を示した。簡単のため、衝突前後の速度ベクトル V、
v
'が地雷に垂直な同一平面内にある場合を考えよう。地面の法線ベクトルを n とすると、註ねか
えり係数はボールの衝突点の速度 v(c)の衝突前後の法線成分 ltc)、14c)'の比を用いて、
11ノ
11
Jset
1
4
j
c
y
en 一
一一一
一一
ぽ
)
で定義される。後に導入する接線方向のはねかえり標数との対比から、らを「法線方向のはねか
えり係数j という名前で呼ぶこともある。
衝突速度等の衝突前の運動状態に依ら
多くの初等物理のテキストにおいてはねかえり係数は、 f
ず、物質の種類のみで決まる一定値である J として導入される
1可。だが、この記述は正しくはな
い。なぜなら、先に述べたようにエネルギ一散逸の裏合は、同む物質でも衝突需の運動状態に強
く依存して決まるものだからである。実繋ほとんどの巨視的な物体同士が正面衝突した時のはお
2，
3，
1
8，1
9，
2
0，
2
1，
かえり係数ほ、衝突速度の上昇に対して単読に減少することが知られている [
~，~， U ，~， W ， V ，~， W ， 3~0 また、簡突角度によってもはねかえり祭数は様々な値を取るこ
とが知られており [
3
1，
3
2，
3
3，
3
4
]、初期回転に故存することも容易に想像できる 3
0
一方斜め禽突においては、法線方向のはねかえり係数だけでなく、接線方向のはねかえり係数も
初期の運動状態に強く依存する物理重である。接線方向のはねかえり係数は、国 1において、ボー
ルの笥突点の速度 v(c)の衝突前後の接線成分 1
5
(
c
)、l
f
(
c
)
'の比を用いて、
1
7
(
C
)
'
q = - 1 5 7 ( 2 )
で定義される。接線方向のはねかえり保数は笥突角度によって -1から 1の間で様々な値を取るこ
3
5，
3
6，
3
7
J、後に示すように衝突角度、法線方向のはねかえり係数、そして
とがわかっているが [
物体表面の摩擦係数を用いて現象論的に取り扱うことが可能である [
3
8，
3
9，
4
0，
4
1，
4
2
]
0
3ベーゴマで遊んだことがある人ならば、ベーゴマがごく題い接触時題の内に、回転自岳度から並進自由夏ヘエネル
ギーをやりとりする様子老何度も晃ているはずである。
-686-
E視的物体の非弾笠筈突
式(
1
)、 (
2
)で定義したらや q のような非弾性衝突を特徴付ける物理パラメータの大きさは、衝
突速度や衝突角震等の初期条件とどのような関係があるのだろうか?また、どのような散逸過程が
支配的になって決まるのだろうか?この間いに対する答を模索する事は物理学や、物理教育にとっ
ても重要であるが、工学の分野においても有用でおると思われる。例えば、粉体工学や土木工学
等の分野で土砂等の粒状物質のシミュレーシヨン法として知られる離散要素法 (
D
i
s
t
i
n
c
tE
l
e
m
e
n
t
Method，DEM) [
4
3，4
4，4
5
]においては、はねかえり係数は現象論的なパラメータの関数として
導入されるが、それは実験との整合性から経験的に決めるものであり、物理的な基礎付けはなさ
れていない。また野球のパットやボール等のスポーツ用品の開発・研究においても、はねかえり
係数はスポーツ用品の反発特性を示す指標として登場する [
4
6
]
0 はねかえり係数の決定要因を明ら
かにすることは、スポーツ用品の設計にも大いに役に立つことであろうし、工学におけるシミュ
レーション法の理論的基礎付けへの貢献も期待できるであろう。
本稿では、これまでに行なわれた巨視的物体の衝突に関する研究を紹介し、同時に筆者らがこ
れまでに行なってきた数値シミュレーションの結果やそれに関する環象論を紹介する。また弾性
体の衝突だけでなく、弾性接触やナノスケーんのクラスター第突のシミュレーション結果も紹介す
る。本稿の構成は以下の通りである。まず次章ではこれまでの衝突研究の歴史を概観し、豆視的
な衝突環象に関する理論的研究をいくつか紹介する。第 3章では禽突研究における数値的方法に
ついて紹介し、特に筆者らがこれまでの研究の中で用いた、ばねと質点のネットワークに基づい
た 2次元弾'性体モデ、ルの詳しい導入老行う。第 4章、第 5章ではその 2次元弾性体モデ、ルを用い
て得られた衝突研究の結果を紹介し、第 5章は分子動力学法に基づいた教視的な系の種突シミュ
レーションの結果を紹介する。最後に第 7章における総括的な議論と今後の研究の発展に関する
議論をもって本稿の結びとしたい。
2 衝突研究の歴史
紀元青Ij4世紀のアリストテレス (
A
r
i
s
t
o
t
l
e
)の力学研究の中においては、くさびに圧力を加える
よりも少しの衝突を加えたときの方が大きな効果を与えるのはなぜかという問題が笠場する向。
この間題は後に 1
7世紀になってガリレイ (
G
.G
a
l
i
l
e
i
)によって実験的に検証されることになる。
7世紀後半にニュート
力学現象を解析的に説明するという理論的研究のスタイルが確立するのは 1
ン(
I
.N
e
w
t
o
n
)が『プリンキピアJl (Philosophiaen
a
t
u
r
a
l
i
sp
r
i
n
c
i
p
i
am
a
t
h
e
m
a
t
i
c
a
)[
4
7
]，を著す頃
7
のことである。巨視的な物体の禽突環象が力学研究の中で、中心的なトピックで、あったのはこの 1
世紀であり、この時期にガリレイやニュートンらによる理論的な考察が集中してなされた。ここ
7世紀以降の笥突研究を中心に紹介する。
では特に 1
2
.
1 17世 紀 の 衝 突 研 究
1638年にガリレイ
(
G
.G
a
l
i
l
e
i
)は著書 F
新科学対話j [
4
8
]の中で衝突の実験的研究を行ない、衝
f
o
r
z
αd
e
l
l
α
突に関する幾つかの経験問をまとめた。そこでガリレイは衝突する物体の持つ衝突力 (
円
i
n6
FO
酉伸寛人
p
e
r
c
o
s
sα
〉を測定する実験を行ない、衝突によって生じる力は静的な圧力に比して大きくなること
3，
49，
5
0
]
0 また、
や、またその力は速度や重量、衝突角度にも依存することについて述べている [
同持期にマルキ (
M
.M
a
r
c
i
)も著書『運動の智合についてJl [
5
1
Jの中で禽突実験を取り上 L
久 1つ
は運動し、 lつは静止している 2つの相等しい弾性球が衝突すると、運動しているもの試運動を
止め、その失われた運動は静止したものに伝わるという実験結果を報告している [
5
2
J
o
その後デカルト (
R
.D
e
c
a
r
t
e
s
) は二鉢衝突を笹突体の質量の違いや運動の方向等で場合分けし、
幾つかの衝突則を導いたが、その中で衝突前後における運動量の保存則に気がついた。ただ、デ
カルトの衝突則は衝突体の弾性と非弾性の区別をしなかったために、正しいものと正しくないも
のが混在したものとなっている [
3，52]0
ガリレイやマルキらの管突実験の結果や、デカルトの衝突期老受けて、非弾性衝突に関する理
論的な考察が 1670年前後にウォリス (
J
.W
a
l
l
i
s
)、レン (
C
.Wren)、ホイヘンス (
C
.Huyghens)
らによって行なわれた。ウォリスは後の微積分法の先駆けとなる求積法の研究で知られるが、非
弾性体の中心額突に対してデカルトの運動量保存員J
I
に基づいて、衝突後の速度を与える式を提出
した。レン4辻、重さの違う撮り子の実験を行い、後にホイヘンスが得たのと同様の弾性誌衝突の
定理を見出したが、理議的な考察は行なわなかった。ホイヘンスは、弾性球の中心衝突に関する
幾つかの法員せを証明なしで発表したが、その中で衝突詰後の活力 5の総和は変化しないという、重
要な原理に到達した。これは現在でいう運動エネルギーの保存期に他ならない。
2.2 ニ ュ ー ト ン に よ る 反 発 係 数 の 導 入
ニュートンは 1687年、著書 fプリンキどア』を著し、有名な f
運動の三法則J として知られる
運動の一般射や、ケプラー (
J
.K
e
p
l
e
r
) が予測していた天体の運動に関する諸法貝J
I
そ万有引力を
ぜから導出するといった、画期的な理論的業績を挙げる。その中で、ニュート
軸とする一般的法R
ンは巨視的な物体同士の衝突に関して、はねかえり係数という物理量を導入した。後は『プリン
キゼア』の中で、作用・反作用の法射を実証するために撮り子の衝突実験を行った。そこで伎は、
弾性体の衝突においては積突後の速度は槙突前の速度に対し比例関孫にあること、またその比例
係数は振り子の材質によって様々な値を取り、衝突速度に依存しないことを発見した。この一定
の比例保数のことはその後はねかえり係数と呼ばれ、接言突物体の反発を示す量として用いられる
ようになる。
その後様々な材質のはねかえり係数を計測するための実験が行なわれるが、その中で 1834年の
ホジキンソン (
E
.Hodgikinson) の実験は Fプリンキピア Jのはおかえり係数に関する記述を検
証するものであった [
1
9
]
0彼はニュートンと同様の装置を用いて、様々な材質の振り子を用いた箆
突実験を行い、はねかえり係数の計測を行った。その結果、どの材震の球においてもはねかえり
係数は禽突速度によらずほぽ一定であると結論付けてはいるものの、第突速度が大きくなると誌
ねかえり孫数がわずかに減少するという観測結果も同時に指摘している。
建築家としても有名であり、ロンドンのセントポール大聖堂 (
S
t
.PaulCathedra
l)等の設計で知られる。
2乗をかけたもの。活力という尼語はライフ。二、ソツ以降使われるようになったという 1
0
52]
4
5質量に速度の
- 688-
巨規的物体の非弾笠衝突
ニュートンとホジキンソンの実験の違いは、実験に用いた振り子の材質にある。ニュートンは、
コルクや毛糸球のような弾荘変形を起こすが塑性変形を起こしにくい材質を用いたのに対して、ホ
ジキンソンは真議等の塑性変形を起こしやすい材質も用いて実験を行なった。その結果、彼は塑
笠変形を起こす材質の方が衝突速度の増加に対してはねかえり係数の減少する富合が大きいとい
うことにも気づいた。
ホジキンソンの実験結果を受けて 1
9世紀末にヴィンセント(J.H
.V
i
n
c
e
n
t
)は弾性変形及び塑
性変形する球が大変形を起こすような倭突実験を行い、いろんな材質の球のはねかえり係数を計
2
0
]
0 この実験では様々な材質の組み合わせによる球と壁の衝突実験が行なわれ、はねかえ
測した [
り係数が衝突速慶に依存して、系統的に変化するという事実が明らかになった。
2.3 衝突の理論
この頃の弾性体の衝突に関する理論的アブローチの先駆けとして、 1
9世紀後半にヘルツ (
H
.
H
e
r
t
z
) によって発表された弾性接触理論 {
5
3，
5
4
Jと、それに基づいた衝突理論を紹介しよう [
5
5
J
o
後者は弾性接触理論に基づいて、衝突時間中に弾性球に生じる変形の大きさの時間変化を追跡す
ることで、衝突の譲触持関を求めるというものである。
まず最初に弾性球の接触理論を笥単に紹介する。これは 2つの弾性球を圧縮力 F
e
lで垂直に押し
)。へんツの接触
付けたときに生じる変形 hや接触域の半径 α との関孫を定めるものである(図 2
Z
I I
ト
.
.
f
a
図2
:球同士の接触。
理論によれば、司質な 2球関士の接触における圧縮力えtと、変形 h及び接触半径 αの関係は、
守/司 r
_ ー
1 1¥11/3
h~ = F
:
!~ I
D~ I
十
:
.
_
， )I ，
<) (
L
11
究
=
1
- ¥R .R
'1
1
I_RR' ¥1/3
Fア ID一
一一ーで i
t
¥-R+R'J
(
3
)
(
4
)
となる。ここで R及び R
'は 2球の半径、 hは和 R÷R'と2球の中心開の距離との差、 D は 2球
3f1-σ2 1ー σ'2¥
のヤング率 E、E とポアソン比 σ、σ を用いて、 D
=:;(一一一十一一-Iで表される物性誼
4¥ E
E' I
である [
5
4，
5
5
]
0 これより 3次元の場合は 2球の変形んが圧縮力 Felの 2/3乗に比例することがわ
-689-
題
イ
中
寛人
かる。また、接触域の圧力分布は、 2球の重心を結ぶ方向を z轄とし、接触域は xy平面にあると
すると、
3F
e
l
/
.
.
X2十 y2
Pz
(x，
y
}=云
ーτ¥/
1一
一-'
)-一
v
α
ι
L.
(
5
)
α~
となる。
逆に見ると、同質な 2球同士を互いに圧した時に、球 i
こ働く弾性力は歪みの関数として表すこ
とができる。そこで、 2球が正面衝突した時に接触詩問中に受ける力を上記の接触理論で求まる
弾性力として、んの時間発展方程式を構成してみる。以緯の計算は、ランダウの弾性論の教科書
[
5
5
]に準拠している。
2球の蜜突時の桔対速度を V として、箆突詩の運動エネルギーが歪み速度で表される運動エネ
ルギーと弾性力
Felによるポテンシャんエネルギーの和と等しくなるとする。つまり、エネルギー
保存射を仮定すれば、
ρhu
、tzノ
五)
f
'
'
、
毛
Idh¥2
μ(
+kh5/2ニ μV2，
R
'
ム
/
ヱ
'
5DVR +R
k=
という微分方程式が得られる。 hは速度 d
h/dtが 0となる時に最大誼 h
o~こ達するので、式 {6} よ
りh
o= (μ/的2/5V4/5が得られる。よって、接触時間 T は式 (6) を hについて 8から hoまで積
分したものを 2倍することにより得られ、適当な変数変換を施した後に楕円積分を実行すること
によって、
5
rv
(
長
Y
'
(
7
)
を得る。すなわち、接触時間は衝突速度の -1/5乗に比例することが示される [
5
5
)。文献 [
2
}では、
いくつかの金属球を用いて衝突速度と衝突時間の関保を実験的に諒べているが、鉄や真鎗などの
7
) の妥当性が示されている。一方、はねかえり保数の方は、この解
金属を用いた実験では、式 (
こしかならない。
析においてはエネルギーが保存するために 1t
2.4 準静理論
1980年代後半になって、桑原 (G.Kuwabara) と河野 (K.Kono) は弾笠体が低速で正面衝突
する際、の穣突速度とはねかえり孫数の関係を理論的に考察し、ヘルツの理論に基づいた準静理論
2
5
]。弾性体同士が正面衝突する際に、弾性体内部に生じた運動が無認小の速度で起
を提案した [
こるなら、弾性体内部は平笥状撃が保たれたまま運動が継続される。だが実際には、弾性体内部
の運動は有限の速度で起こるので、衝突の瞬間に局照的な平衡状態からのずれが生じ、その結果
全体的な平衡状態への緩和が起こる。つまり、衝突の瞬間に力学エネルギーから熟への不可逆的
- 690-
巨視的物体の非弾性衝突
な散逸が起こるわけである。もし弾性体内の運動の速度が極めて小さく、エネルギー散逸の度合
いが微弱であるならば、速度の 2次関数で表される散逸関数を導入し、その速震による一司微分
5
5
]。
で、内部散逸に起因する“散逸力"を定義することができる [
そこで正面衝突する 2つの弾性球が、ヘルツの接触理論で表される弾'注力に加えて、弾性体の
内部粘性;こ起因する散逸力を受けると仮定すれば、変形の時間発展方程式からはねかえり祭数を
計算することが可能になる。桑原と河野の計算によれば低速の正面衝突の場合、穣突速度とはね
かえり係数の関孫は、 1- ec
xV1/5 となることが示される。
ここで簡単に文献 [
2
5
]に基づいて準静理論を導出しよう。前述したように準静理論では、衝突
時間中に 2球関;こ働く力を弾性力と散逸力の和で表すが、前節のヘルツの接触理論で導入したよ
うに、弾性力 F
e
lは式 (
3
) から球の変形 hの 3
/2乗で表される。一方散逸力の方も、弾性力の計
算と同様に基本的には圧力分布を接触域で績分する事で求めるが、弾性力と違うところは、そう
やって求めた散逸力が歪みではなく、歪み速度に比例するという事を復定するところである。そ
h
h
，
.q“
必
，
FF
，
，
，hH
n4
，
可t
t
zftB'BJ
一 司 ，e
一
R
'g
﹁E132222L
R一九
2一
U
J
一
一
，
F
μ 一十
ういう仮定の下で得られた散逸力は、
(
8
)
で与えられる。ここで R1、R2は接触する 2球の半径で、・ 5
ま時間微分を表している。また D'は
、
3
)、 (
4
)の D に椙当し、粘性係数で定義される重である。
式(
こ粘弾性球に犠く力 F を F ニ Fel+Ffで表せば、球に生じる歪み hの時開発展を
さて、笥突時 i
次のような微分方程式で表すことができるだろう。
d2h
μ友E
;
W
2
h
+
;
'
T
.
ニ初/h
max で書き直すと、
2
dχ ょ "rv1/2~主ム 23/2 ーの
d'T2
I
d'T
I/¥.
、
、
，
，
，
この微分方程式を無次元パラメータ χ = h/h
max.
、
，
/
a
E
・
- 5D' IRl十 R21
42A
f
/
2
土日ざと
ー
ti
.
AU
rR
I
R
:
・
‘
，
，、
k
'
1
/
2
.
1
5DIRl+R21 '
土
i
宅
ー
。
)
r
，.‘、
k =
0，
杭 3人
十
I
4
/
¥
.
(
1
2
)
，
.
5
.
'f V ¥
I~
I
2
'
.い
3μ2)
(
1
3
)
γ~"
となる。ここで無次元住には内部摩擦がないときの歪みの最大誼 hmax と衝突速度 uを用いた。
ここで式 (
1
2川 、 無 次 元f
ヒした初期速夏 γを与え、数値的に解くことで種類麦の d
χ/
d
tからは
ねかえり係数を求めることができるが、 γが小さいと仮定すれば解析的に解くことができる。式
tを乗じてァ =0からアニ
(
1
2
)の両辺に d
X
/d
l叶
TC
d
/
1/2(ωT判
まで積分することによって、
'
T
c
ヤ=イ九112怜 /d)2d
1
明x5
となる。
'
T
'
T
(
1
4
)
QU
FO
醤伸寛人
左辺は時刻 T= O，
T
χ/
d
T
l
r
=
O= 1と
cにおいて χ= 0であるという条件を課し、更に初期速度 d
2
T
l
r
=
r
c= uとすれば、 (
1
/
2
)
(u - りとなる。
いうことに注意して反発速度 dX/d
一方、右辺はァ = T
Oで χ=1となることを利用して、
寸
1
;cJ2州 合 ず 針
-2
0X1/2(dX/dT)dX
11
一叶;
w
f
叩
dX
(
1
5
)
(
1
6
)
2
/
5
)
B
(
3
/
5，
3
/
2
)
2
γ(
(
1
7
)
1
.
0
0
9
γ
(
1
8
)
となる。ここで散逸が小さいと仮定し、無散逸の時のエネルギーバランスの式(めを無次元f
とする
ことで dX/dT~ (
1-X
5
/
2
)
1
/
2とした。また、 B(p，
q
)はベータ関数である。
はe
三
以上より、はねかえり係数 ε
l
u
lで定義されることに注意すれば、
ε =
~
(
1-2
.
0
1
8
γ
)
1
/
2
(
1
9
)
1- 1
.0
0
9
γ
(
2
0
)
c1- v1/5 が示される。
となり、 eo
彼らは実際に様々な材質の振り子を用いた衝突実験を行い、乙の理論式の有効性を実証してい
るへその後、異なったアプローチによるいくつかの準静理論の講築がなされ [
2
6，
27，
2
号、準静理
論は抵速の正面笥突においては、有効な理論であると信じられている。
2
.
5 塑性変形を伴う高速衝突
桑原と河野の解析は、物体が十分小さい速震で衝突するという仮定に基づいてなされたもので
あるが、ジョンソン (K.L
.J
o
h
n
s
o
n
) は、高速の正面衝突における禽突速度とはねかえり係数の
]
0 ジョンソンの解析は、衝突時に起こる塑性変形のエネル
関祭を次元解析を墓に晃積もった[18
ギーを見積もることによって衝突後のエネルギーを計算し、衝突話のエネルギーとの比からはね
かえり係数を計算するというものである。
ここで簡単にジョンソンの解析方法を詔分する。 2つの球同士が、梧対速度 V で互いに正面衝
'ではねかえる状況を考えよう。そして、衝突前に持っていた運動エネルギーは
突し、相対速度 v
徳突によって 2球に局所的に生じる弾塑性変形に吸収されると仮定しよう。すると、欝突によって
球に生じる最大の変形を rとし、接触面上に動く力を F とすれば、
トv =1
0 F附
2
0
(
2
1
)
6 とはいえ、塑性変形を起乙しやすい金属で、しかも降伏強度の小さい真畿の振り子の実験結果は準静理論との一致
はあまりよくない。
- 692-
E視的物体の非弾性笥突
となる。ここで μは 2球の換算質重である。一方で、最大変形が起こった後に反発の運動エネル
ギーが、弾性回毎によってなされるを事に等しいとすれば、
i
川
となる
0
二
{〆
F
的
印
州
州
'
(
<
<
5
'
判
(
仔
正
6
め
ゲ
的
'
)
d
<
5
'
伊
ここで'(プライム)のついた量誌反発後の量を示すものとする。従って、 F(めと F
'(
<
<
5
'
)
がわかれば、
v
'と V からはねかえり孫数 e= -V'/Vが求められることになるが、圧縮力が弾性
限界を超え、塑性変形が起こるようになると、平均の接触庄は 1
.1Yから 3.0Y程度の一定績を耳元
Yは蜂伏応力)。ジョンソンの解析においては、変形 6が接触面の半径 αと
<
<
5=α2/2Rの関係
る (
にあることと、平均の接触圧が 3.0Yに等しくなる事を薮定して、式 (21)、 (22) から次元解析に
よりはねかえり孫数と衝突速度の関係を導く。
1
4
x
:V/
その解析によればはねかえり係数と衝突速度の関係は、 ec
となることが導かれる。実
際にスチーんの金属球を様々な種類の金属で出来た板に正面衝突させて衝突速度とはねかえり孫
数の関係を探った実験においては、材質の種類によらずはねかえり係数が衝突速度の -1/4のべき
で減少することが確認され、この解析の脊効性が示されている [
2，1
8
]
0
2
.
6 斜め衝突
これまでに扱ったのは弾'往体の正面衝突で、あるが、弾性体が斜め笥突した時の反発の挙動も様々
な手法で調べられている。斜め衝突の反発を表すパラメータとして用いられるのは、式 (
2
) で導
2
) においてもずの'は、衝突によって励起される
入した接線方向のはねかえり係数々である。式 (
回転の角速度ぷを用いて、
vjc)'=V'-RJ× nと表すことができるととから、
e
tは衝突によっ
て励起される角速度以によって決まることがわかる。
ウォルトン (
0
.長. Walton) は、シンプルな開体球同士の衝突を考え、運動量保存期と角運動
3
8，
3
9
J
o ウォ yレトンの
量保存貝りを仮定して接線方向のはねかえり係数と衝突角度の関祭を導いた [
解析によれば衝突角度 γ と etの捷径は、次のようになる
{
e
t~
I
Q
2¥
MR
J-1十.
μ0(
1十 ，
e
o
tγ(1+ー1
ァ)
) (
7と γ
0
)
<
，
-.
-)c
.
.
.
¥
。
I
V
I
¥ - .
1 _ IV/
(
2
3
)
γ
( 豆γ
0
)，
〔β
ここで γ
、
。 M 、R、向、 I辻、それぞれ臨界角、質量、半径、摩擦係数、慣笠モーメントを表す
[
3
8，
3
9
}。この関孫式によれば、衝突角度が臨界角 γ
oよりも大きい時には、衝突中は摩擦係数 μ。
でスリップが起きるが、小さいときにはスリップが起こらずに笥突体に十分な屈転が生じ、 etは
正の最大誼 β。という値を取ることになる。 etはこのように衝突角度によって様々に変f
とするパラ
メータであるが、その最大値角は材質の種類や表面の摩擦によって決まる悟であり、物質の斜め
3
5，
36，
3
7
]。近
衝突を特徴付ける量として実験によって様々な物質の βoの植が計測されている [
年、斉藤 (
K
.S
a
i
t
o
h
) と早川 (
H
.Hayakawa)は 2つの円盤がお互いに斜め衝突した状況を考え、
衝突中の接線方向の歪みに関する運動方程式を解析することで β。の植に理論的な説明を与えるこ
とに成功している [
5
6
]
0
- 693-
罰停寛人
一方でモウ (
N
.Maw)、 パ ー バ - (
J
.R
.Barber)、フォウセット (
N
.F
a
w
c
e
t
t
)[
4
0，
4
1
Jやスト
ロンジ (w.J
.8
t
r
o
n
g
e
)1
3
]は簡突角度を 3つの領域に分裂し、その各々の領域で衝突中の球は、
(
i
)壁と居着状態にある(スリップが起こらない〉、 (
i
i
)スリップ状態と居着状態が共存する、 (
i
i
i
)
完全にスリップの状態を取ると仮定した。そして、法線方向と接線方向に働く力を計算すること
で、衝突が終わった後の衝突点の速度を計算し、そこから接線方向のはねかえり採数を計算した。
この理論はウオルトンの隣体笥突の理論とは違って、弾性体の衝突のプロセスを考嘉に入れた理
論であるが、定式化がウォルトンの理論に比べて煩雑であるという欠点がある。
2
.
7 1を超えるはねかえり係数
斜め衝突の実験においてはねかえり孫数と衝突角度の関保を調べると、興味深い現象が観測さ
れることがまうる。はねかえり係数は I以上の誼そ取らないというのは、初等物理の教科書にも載っ
ている常識だが、斜め衝突の場合には 1以上の舘を取り得ることが近年報告されている。
{
.
;大のんージュ (
M.Y.Louge) とアダムス (M.E
.Adams) のグループによって近年行
コーネ J
なわれた実験は、アルミニウム酸化物の球を斜めに傾けたポリカーボネートの坂の上に自由落下
3
1
J。その結果、彼等ははねかえり孫
させ、種突角震とはねかえり係数の関孫を謁べたものである [
4
.
9度に達すると lを超え始めるという
数が衝突角度の互援に対して線形に上昇し、衝突角震が 7
一見奇妙な結果を得た。また彼等は板の角度を 86.2度まで額けて実験を行い、はねかえり係数は
最大1.3にまで達することを報告している。後らの用いた実験材料は、球に比べて壁の方が軟らか
い材質であり、衝突持に壁表面に局所変形が生じることが確認されている。彼らはこの壁表面に
生じる局所変形が球の軌道を上向きに変えることによって、このような現象が起こり得ると述べ
ているが、実際にこの仮定に基づいた解析を 5章で紹介する。
この実験以外にも、エアテーブル上に円盤を置いて、テープ、ルの端に取り付けた壁に斜め衝突
させて、錆突角度とはねかえり係数の関係を謁べた実験がある
[
3
2
]
0 この場合はルージュとアダム
スの実験結果のように、はねかえり係数は衝突角度に対して線形に変化はしないものの、ある衝
突角度になると急に 1を超えるという結果が得られた。この実験で用いられた実験材料は、ルー
ジュとアダムスのそれとは逆に壁の方が円盤よりも硬い材質であり、この両者の硬さ比の違いが
主に実験結果の違いを生じさせていると考えられる。
2
.
8 その他の研究
ここまでで紹介してきたのは国体問土の衝突であるが、近年それ以外にも様々な形状の物誌や、
高分子や薮体等の衝突現象の研究も行われている。例えば f
水切り j は、河原などで平たい石を J
1
面i
こ水平に投げて石の反跳回数を競う遊びであるが7、近年ボケ (
L
.Bocquet) のグフしープは自転
する円盤を水表面に衝突させ、衝突時に円盤と水表面のなす角が 2
0度の時に最も反跳の回数が多
いという、興味深い実験事実を報告した [
5
8
]
0 その後、永弘 (
8
.Naεahiro) と早川 (
Y
.Hayakawa)
7現在
(
2
0
0
8年 5月〉の霊界記録は、 2007年 7月に R
u
s
s
e
l
lByarsが達成した 5
1回 [
5
7
J
o
- 694-
巨視的物体の非弾性萄突
によって、 SPH法に基づいた円盤と液体の衝突シミュレーションが行われ、やはりこのお度と
マジ、ツクナンバーj であることが示された [
5
9
]
0
いう衝突角度が最も多い反跳回数を引き起こす f
更に彼らは円盤が水面と接触したときに円盤が水面から受ける力に基づいて円盤の運動を記述す
る微分方程式系を構築し、その解析により fマジ、ックナンバー J~こ説明を与えることに成功して
いる。
6
0，6
1J
、水滴と板の衝突 [
6
2
}、大変形を起こすようなゲ
この抱;こも球や円柱の水面への笥突や [
6
3
]等の興味深い研究も行われており、惑星科学の分野では給体層や多主性物質に球
ル球の積突 [
体を様々な初期条件で衝突させて、天体表匡のクレーターの形成過程を語べる実験やシミュレー
6
4，
6
5
J
o
ションが行われている [
3 数値モデル
この章では、衝突問題に対してどのような数値的手法が用いられるかを紹介する。車や携帯電
話等の工学的な衝突開題の数値的解法としてよく知られているのは、おそらく有頭要素法による
解析ではないだろうか。有援要素法のソフトウェアは市販のパッケージ [
6
6
]も多く存在するので、
データを入力するだけで手軽にシミュレーションを行なえるという利点があるが、物体の形状が擾
雑だとメッシュ分裂に時間がかかり、また破壊の現象を扱うのが圏難だという不利な点もある 80
様々な形状の臣体の衝突や、それに伴う塑性変形・披壊を扱うためには、題体を粒の集合体と見
倣した粒子ベースのモデリングが直感的である。酉体を粒子の集団で表す方法には後に紹介する
ばねー質点モデん [
6
8
]や粒子法 [
6
9
J等の方法が知られている。また、粒子の集合で物体の形状を
表現し、粒子関梧互作用として離散要素法(第 1章参照〉で用いる相互作用に、水やモルタルなど
の効果を模倣した信機ばねを付加して、国体の破壊挙動をシミュレーションする方法もある f
4
4
J。
物理学の傑から、弾性体のはねかえり係数と衝突速度の関孫を謁べた数値的な研究の傍として、
ゲール (
F
.Gerl)とジッベリウス (A.Z
i
p
p
e
l
i
u
s
) の 2次元衝突シミュレーションの研究が挙げら
3
0
]
0 彼らの衝突モデルは開放端の境界条件を持つ弾'往円盤のモード解析に基づいたもので
れる [
ある。彼らのモデルは散逸を含まないハミルトンモデルであり、衝突後に内部に生じた弾性振動
を散逸と見倣し、衝突前後の並進速度の比からはねかえり係数が計算される。数誼計算は弾性場
と壁ポテンシャルとの相互作用を含めたハミ Jl
;トン方程式を解くことによって、円盤の重心の時
開発展を計算する。彼らの計算では、はねかえり係数は衝突速度が上がるに連れて単調減少する
こと、またその僅は衝突物体のポアソン比に強く依存することが明らかになった e
ここで、筆者らが数値実験に用いたばねー質点モデルを紹介する。ばねー質点モデルとは多数
のばねで連結された質点によって弾性体を表現するもので、主にコンビュータグラフィックスの分
{
6
8
]
0基本的には各糞点が連結されているばおから受ける力を計算する
主体全体の運動の時霞発震を求めることができ
ことによって設置や速度の持隠発展を計算し、弾f
野でよく用いられている
る。このようにして弾性体をモデル化することに辻、次のようなメリットが島る。
8メッシュ自体を破壊することで、物体の塑性変形や破壊を表現する有限要素スキームも開発されている[
1
0
67
Fhu
QU
FO
酉仲寛人
図4
: ランダム格子モデ、ル。
図3
: 三角格子モデ、ル。
1.モデリングが産感的で、アルゴリズムが簡単である。
2
. 様々な形状の弾t
'
主体を表現することが出来る。
3
. 破壊や塑性変形を表現することが出来る。
4
. 混度分布や応力分布等を簡単に計算することが出来る。
またこのようなモデ、ル化の欠点として、一部の質点が大きな外力を受けたときに格子が大きく変
形し、変形が元に戻らなくなるという現象が起こることがある。これは、ばねに一定以上の力が
6
8，70，7
1
]。我々の研
働くとばねが切れるというルールを導入すれば解決することが可能である 1
究で、はばねを切断するルールを導入せず、変形が元i
こ戻るようなパラメータ領域でのみシミュレー
ションを行なっている。
2次元のばねー質点モデんを構築するためには、まず質点をどのように初期配置させるかという
問題があるが、その一つの方法として質点を三角格子状に配置する三角椙子モデ、ルが考えられる
(国功。三角格子モデルを構成するためには、まず質点を円内に三角格子状に配置した後で近傍の
質点同士老ばお定数 κの線形ばねで接続する。この時、モデルのヤング率は E=2
κ
/ゾ三、ポアソ
/
3~こなることが分かつている [72]0 各質点の運動方程式は次のように表される。
ン比は ν =1
.
'
)
m
6
今=ー
κ )~ (
尚一 I
r
i- r
il)三三乙
配 信
J"lri-rj!
ここで r
iは質点 tの位室、 jは近譲の粒子のすフィックス、 m は質点の質重である。この運動方程
式を適当な初期条件の基で解いて質点、の位霞と速度の時開発展を求めることによって、弾性体の
運動や変形を追うことが出来る。これまでの筆者らの研究においては、物体が壁に正面衝突した
擦の速度とはねかえり諒数の関係を謁べるために三角格子モデ、ルを指数関数型のポテンシャル壁
に正面衝突させるシミュレーションを行った [
7
3，
7
4，7
5，
7
6，
7
7，7
8
]
0 また、第 2近接の相互作用
まで寂り入れた正方格子から成る弾性壁を構成し、三角格子モデルを正面衝突させた時のはねか
えり孫数と弾性壁の厚さの関係も調べた [
7
6，7
9
]
0
一方、筆者らの斜め霞突のシミュレーションにおいては、質京、をランダムに配霊したランダム格
子モデルを用いるく園 4)。これは、斜め衝突の持に三角格子や正方格子のような規則性のある格
-696-
巨視的物体の非弾性衝突
F(
l
s
)
I
F
2
図5
:弾性円盤と弾性壁の相互作用。
子を用いると、初期のモデ、んの配向によって結果のばらつきが大きくなるためである。格子をラ
ンダム格子にするとそのようなばらつきを抑えることができる。このモデルを構築するには、ま
ず一様乱数を用いて資点を円内にランダムに配置する。その後、近傍の質点同士を連結するため
8
0
]を用いる。連結された質点関には次のような非線形ばねによる復元
にデ E ーニ一三角分割法 1
力が働くとする。
f
αx-kbX3
i
j =- k
(
2
4
)
ここで、 x はばねの自然長からの停び、んとん誌ばね定数である。非線形ばねを用いた理由は、
モデルが大きく変形したときにばね両端の質点の位置が入れ替わり、モデル全体の変形が戻らな
くなるのそ防ぐためである。また、ここで扱うシミュレーションにおいてはモデルの全体的な変
形が大きくないため、非線形項の効果によるモードのミキシングの効果はほとんどないと考えら
れる。
斜め第突のシミュレーションでは、上記のようにして講成した弾性円盤の他に、開様にして作っ
た弾性壁を用意する。これは質点を長方形内にランダムに配聾すること以外は、弾'往円盤と間接
)。まず円盤の下
の方法で作成する。弾性円盤と弾性壁の梧互作用は以下のように導入する〈図 5
半分の縁に位置する全質点に関して、真下にある壁表面のばねとの距離を時間ステップ毎に計算
t
hよりも小さければ次の反発力を受けるとする。
しておき、もしその距離が闇値 l
S
F(
l
s
)= α01
もexp(-aols)n
(
2
5
)
ここで、ぜはばねに垂直な方向の単位方向ベクトル、んは質点からばむまでの距離、 α。は円盤の
半径を R とすれば 300/R
、%は 1次元音速を cとして、 αomc2R/2という径を用いる〈図 5参駒。
また爵植 l
t
hの値には円盤を構成する全てのばねの自然長の平均値を用いた。指数関数以外にも斥
力のモデル化の方法として、例えばレナード・ジョーンズポテンシャルの斥力部分を用いること
もできるが、ここでは物体罰土の貫通老婦げる強い反発力のモデルとして指数関数を用いた。
質点からばねに下ろした垂線の足から、ばお両端の責点 1と 2までの距離を i
いらとすれば、ば
s
; _s
F
Fl
l
s
S
s
.
.
.
- 11 n " F2
(
l
s
}一 ー
がという反作用を
1十 l
l
/l
2 、
SI
1+l2/む
受けるとする。これは垂線の足を中心にしたトルクの和ががゼロとなるように決めたものである。
ね両端の質点 lと 2はそれぞれ F1(
l
s
)=
- 697-
T.l {l ¥ _
国件費人
以上より、円盤下半分の縁上にある糞点 tの運動方程式は次のように書ける。
1
2_
Ni
プ = 乞 {-kaXij-kbXTj}十 仇 一 収 ( ら )
η
ここで、
(
2
6
)
Xijは質点 ZJ間の自然長からの変位であり、。 (
lt
h-l
s
)は以下のようなステップ関数で
ある。
J1 i
fl
s:
:
;l
t
h
B
(l
仇 ls)= <
I0 i
fl
s>l
t
h
'
また Niは質点 tにつながるばねの数である。
モテ、ル表面のラフネスは物体問の憲擦力を生み出すために導入したものである。ランダム格子
モデルによる斜め衝突のシミュレーションにおいて、表面のラフネスは重要な役割を果たしてい
る。というのも、 2
.
6節で述べたように環実の斜め衝突では、衝突物体関の摩擦力が衝突後の運動
に影響を与えるからである。ラフネスはモデル表面の質点の位量を標準偏差 3X 10-2Rの正規乱
数を用いて元の位置から変位させ、導入している。
4 接触及び低速衝突のシミュレーション
ハミルトン系に基づくばねー質点系の弾性体モデルは、巨視的な弾性体の斜め衝突の現象論
[
3
5，
38，
3
9
]や、衝突角度とはねかえり係数の関孫の実験結果 [
3
1
]を定性的に再現することに成功
3
3，42，8
1
]。だが、モデル自体がエネルギー散逸のメカニズ、ムを含まないハミルトンモ
している [
デルのため、弾性体の接触理論が示すような接触力と変形の関係老再現することはできない。ま
た、接触理論に基づく低速衝突の準静理論も、物体内部の散逸メカニズムを含まないため再現で
きない。この章では固体中の格子欠陥に椙当するメカニズムをこれまでのモデルに導入し、エネ
ルギー吸収壁に接触させたり毎速で正面衝突させることにより、巨視的な弾性体が示すと考えら
れている理論老、ハミルトン系のモデルでも回復できることを示す [
8
2
]
0
4.1 数値モデル
この章で用いる弾性体モデルは前章で導入したランダム格子モデ、ルに基づいており、円盤と壁を
2
4
) で表される質点跨相互作用に関し
構成する質点数はそれぞれ 1099、 1269である(図的。式 (
2j
i
ω
)二1.0X 102Mc
i
m
)とすると、 kf)二 1
.0xMc
R2、k
2j
て、円盤と壁のんをそれぞれげ)、 k
R2
という値を用いた。ここで M は円盤の質量、 cは 1次元音速を表している。また、ばね定数んは
ん=ん
X
10-3jR2 とすることで弱い非線形項を導入する。式 (
2
6
) の運動方程式を解くために 4
次のシンプレクティック積分法を用い、時間刻み d
t二 10-3Rjcで数値積分を行なった。
ここで格子欠諮を模倣した欠陥粒子をモデルに導入する。欠箔粒子を導入するには、まずモデ
ルを構成する質点のうち数調をランダムに選ぶ。次にその質点 tにつながってる民本のばねのう
ち、民 -1本のばねを切断する。これによって、質点 4は l本のばねだけで接続された状態にな
り、佐の質点に比べてランダムな運動が可能になる。複数の欠諮粒子がランダムな運動をするこ
- 698-
巨視的物体の非弾性蓄を突
とによりモデル内部に生じた弾性波が散乱され、内部状態が速やかに平傷状態に達することが可
能になる。このことはあたかも毘体中の稽子欠絡や不純物が、国体内部のフォノンを散乱したり
吸収することによって、速やかに内部の状態を熱平傷状態に達しやすくする現象に類叙している。
後に示すデータはほとんどが欠焔粒子を 10留ずつ弾性円盤と弾性壁両方に導入したものである
が、欠陥粒子数の接触シミュレーションの結果に対する依存性を後で議論する。
更に弾性壁は弾性円盤に対して十分大きな壁の部分系であると晃散し、それに見合った境界条
件を設定する。境界条件は壁の両サイドと底に自由境界条件を設定する。その時、壁の患の左右
両端の粒子のみを国定する。この弾性壁が部分系であるためには、両サイドと底から弾性波が外
倶J
I
に逃げて、戻ってこないという状況を実現する必要がある。そこで、壁の両サイドと底にある
粒子に関して、エネルギー許凍 Jを J =乞 {
(
l
j
2
)
m
v
;+寸 Vi で定義する。ここでけま粒子の
インデ、ツクスであり、
Vi はその速度、正tは粒子五が接続されているばねから受けるエネルギーの
総和を表している。壁の境界面の外向きの法線ベクトルを llb とし、 J.llb>0という状態を境界
面で実現するためには、粒子の運動方程式を数値的に解く際に、各時間ステップにおいて
Vi がモ
デル内部に向いていたら引を 0にリセットする。するとエネルギー流束は常に正の笛を取ると考
えられ、各境界面において J.llb>0の状態を実現することができる。
4
.
2 接触のシミュレーション
まずはこの弾性円盤と弾性壁を用いた接触シミュレーションを行なう。ここでは弾性壁の高さ
を 4R、幅を R~こ設定する。弾性円盤を外場 P において弾性壁と接触させて、振動が緩和した後
に、外場 P と円盤の変形 6の関孫を謁べる。従ってこの接触のシミュレーションの場合、質点の
2
6
)の在辺にー (
P
j
N
)まを加えたものとなる。ここで N は、弾性円盤を講成する
運動方程式は式 (
099であり、主は壁表面の単泣法線ベクトルである。円盤が円の状態を初期状態とし
質点数 N =1
て各質点の運動方程式を解くと、円盤の重心の振動ほ徐々に緩和していく。定常な振動に落ち着
いたら円盤の重心から接触冨までの距離 R
dを計測し、全体の変形 6三
R-R
dを求める。この作
.
7
7X 10-3πRE
本から P=1
.27X 10-3πRE
ネまで P を変化させて行い、 P
j7fRE*と
業を P ニ 5
2
d
j
Rの関孫を求める。ここで、 E
*は有効ヤング率 E E
j
(
l-v )であり、 Eはヤング率、 νは
孝 二
ポアソン比を表している。
:欠摘を含む弾性円盤と弾性壁の 2次元モデ、ル [
8
2
J
o
図6
QU
QU
FO
圏停案人
0
.
0
7
S
i
m
u
l
a
t
i
o
n
M
喝
H
e
r
t
z
i
a
nC
o
n
t
a
c
tThe
o
r
yo
!
哀 0
.
0
5
0
.
0
3
0
.
0
0
6
0
.
0
0
8
0
.
0
1
P
l
1
tRE
0
.
0
1
2
図7
:圧縮力 P/
πRE*と変形 b/Rの関係。×印は 1
0サンプルの円盤を用いたシミュレーション
結果の平均彊であり、エラーバーはその標準偏差である。
図 7は P/
πRE
ネと b/Rの関孫を求めたものである。シミュレーションは質点の配置が異なる丹
盤を 10種類用意して、シミュレーションデータを平均した。×印がその平均値であり、エラー
バーは標準偏差を表している。また、欠描粒子の数を弾性円盤、弾性壁共に 10{国ずつに設定し
たときのデータを示した。
5
3
]との整合性老調べてみる。 2次元の Hertzの接触理論によ
このデータと Hertzの接触理論 [
3
0
]
0
れば、圧縮力 P と歪 6の関孫は、次のように表される [
云(ベヰR) v
}
b
'
:
:
!
.
-1-
(
2
7
)
実線は、式 (27)をジニ 0.336でプロットしたものである 90 この結果より、我々のシミュレーショ
e
r
t
z理論を自復できていることがわかる。
ンデータはフィッテイングパラメータなしで 2次元の H
このことから、我々の弾性体モデルは接触理論が予言するような平衡状態老実現できると結論付
けることができる。
更にモデルが平衡状態に至るプロセスを特徴付けるため、円盤を構成する質点の速度分布関数、
及びそこから定義されるエントロピーの時間発展を調べた。速度分布関数は、質点の位置と速震
L
摘する分布頭数川川、収対して、仇川三
j
二dイ〉m ，M
，
M 定義
する。 f(x，
y，
vx，
vy，
t
)は時刻 t-6R/cから tまでのデータを平均して求める。
国8
(
a
)及び(めはそれぞれ、 f('i7x，
匂，
t
o= 6R/c)、f(vx，
匂，
t
max= 60R/c)をプロットしたもの
(
a
)は接触初期のデータであるが、%が正と負の領域に大きなピークが一つずつ見ら
である。図 8
れる。これは、壁方向に進行する糞点群と、反射波として帰ってきた質点群の二つが混在している
(
b
)は円盤重心の援動がほとんど緩和した状態の時の速度分布関数
ためと考えられる。また、国 8
を表しているが、ここでは速度分脊関数はガウシアン型になることがわかった。実際 f
(vx，
0，
t
m
a
x
)
は分散 0.02のガウシアンでフィットできることがわかっている。
9モデ、ルの示アソン比を求めるには、まず M
athematica等を用いて円盤の全費点の座諜を読み込み、等方圧縮や単
純シアをかけることによりエネルギー密度を計算する。次にモデルの等方性を仮定すれば、ラメ係数が 2つ決まるので
そこからポアソン比が求められる。
- 700-
巨視的物体の非弾性萄突
0
.
0
2
(
a
)
0
.
0
1
、
ぷ
b
》
Vx
可
0.004
qA
れU
n
u
n
u
(器巴# AhKMkrMM
。
図8
: 円盤を構成する質点、の速度分布関数。時刻はそれぞれい)
t=布、 (
b
)
t= tmaxo
:
b
v h
更に求めた速度分布関数を用いてエントロビーの詩題発展を計算する。エントロゼー S
(
t
)は
州三一仁川
恥い
仇帆で定義する
B
また、欠陥子の数をむ憶から
1
5倍まで変色させて、エントロビーの欠陥数依存性についても調べた。函 9はエントロピーの待
問発展の様子を示している。函 9(a)は欠陥粒子数が 8偉から 10倍のモデルを用いた場合であり、
図 9(め は O個と 15個の場合である。図 9(a)によれば、欠焔粒子の数が 8悟から 10個の場合の
エントロピーは単調増加で、時現~t =
30R/cを過ぎると誌ぽ同じ値 1
3
.
5に近づく傾向が見られる。
欠陥粒子の数がこれより少なiJ'Iったり、多かったりするとエントロビーは単議増加の傾向を示さ
(
b
)は欠蕗粒子の数が 0個及び 15個の場合のエントロゼーを示している。欠結粒子が
ない。図 9
G個の場合は、欠陥粒子による弾性波の散乱の効果がないためなかなか最大誼に達せず、しかも
単調増加の傾向を示さない。欠陥怠子が 15鋼の場合は、モデル自体の強度が弱くなってしまい、
モデルの一言問こ生じた変形が元に戻らないままグローパルな振動を続けてしまい、内部は平衡状
こなかなか到達しない。これより、このモデルのすイズでは、 H
e
r
t
zの接触理論が予言する接触
態L
の平傷状態を実現するためには、欠陥粒子を 8短から 10信導入する必要があるということがわ
かる白
4
.
3 低速衝突のシミュレーション
最後に弾性円盤を低速で弾性壁に正面衝突させて、衝突速震とはねかえり係数の関採を調べる。
2
6
)に初期速度りを v= 1
.0X 10-3cから v= 1
.0X 10-2cの範留で与えて、外場が
この場合、式 (
ない状態で運動方程式を解くことになる。乙のシミュレーションでは弾性壁のサイズを縦 R、穏
- 701-
国待
寛人
14
13
，
"
，
，
，
，
，
，
，
，
S
(
t
)12
1
1
，
40
・C
20
1
0
4
1
0
m
m政 局
0
8d
e
f
e
c
t
s
9d
e
f
e
c
t
s
1
0d
e
f
e
c
t
s
，
，
，
，
，
，
6
0
50
60
1
50
1
4
.
5
1
3
.
5
，
e
，
，
.
.
S
(
t
)1
2.
5
，
'
11
.5
odefect
1
5d
e
f
e
c
t
s
1
0
.
5
2
0
刻一州
AV
-
。
40
ω
w
図9
:エントロどーの時間発展。それぞれ導入した欠陥粒子の数が(吟 8から 10留の場合と (
b
)
G掘と 15個の場合を示す。
4Rに設定した。
、 1次元音速 cで、スケールした苦言突速度 v
j
cとはねかえり係数 eの関孫を示している。
図 10誌
このシミュレーションにおいても、質点の配置を変えた 10種類の円盤を用いて得られたデータ
を平均しており、×印は平均語、エラーパーは標準編差を表している。この図からむかるように、
はねかえり係数 εは密突速度 v
j
cニ 7.0X 10-3程度までの範囲ではわずかに減少するのみだが、
それ以上の速度範囲では急に減少し始めることがわかる。
2
6，
27，
3
0
]の 2次元版と比較してみる。こ
このシミュレーションデータを低速衝突の準静理論 [
e
r
t
zの接触理論由来の弾性力と、内部粘性由来の散逸力
れは衝突時間中に弾性体が受ける力を H
の和で表現し、弾性体の歪みの時間発展方程式を解くことにより、はねかえり保数を計算すると
一
a
ぷ
ω
一日
句
p ω
一
点
P
同
Avrz
M
d
いうものである。 2次元準静理論によれば、円盤の歪みの時間発展方程式は次のようになる。
(
2
8
)
右辺第一項が弾性力で、第二項が散逸力である。第二項の1"
0は散逸の時間スケールを表しており、
フィッテイングパラメータとして導入する。初期条件 6ニ 0、d
d
j
d
t二りでこの方程式を解くと 6
- 702-
E視的物体の非弾性衝突
0
.
9
9
4
0
.
9
9
e
0
.
9
8
6
Simulationゃ
Q
u
a
s
i
s
t
a
t
i
ct
h
e
o
r
y-
0
.
9
8
2
O
.
v
/
c
図 1
0
:喬突速度とはねかえり係数の関保。実線は準静理論による結果。
の時間変化を追うことができ、円盤の受ける力が 0になる時の d
d
/
d
tと u との比からはねかえり
係数を計算することができる。
図1
0の実線は 70=O
.
Ol1R/cの時に、式 (
2
8
)を解いて、 v
/
cと εの関係を求めたものである。
これより、適当なフィッテイングパラメータを設定すれば、シミュレーションデータとの一致はか
c=7
.
0X 1
0
-3 よりも大きくなると、もはや準静
なりよいことがわかる。ただし、衝突速度がり/
理論とは一致しなくなる。高速衝突時に励起される議々な内部モードがこのような不一致を生じ
させると考えられる。
以上、この章では欠陥粒子を導入したばねー質点系の弾性体モデ、ルを用いて H
e
r
t
zの接触理論
に整合するようなハミルトンモデ、ルを構築した。接触の平衡状態に至るまでのプロセスを特徴付
けるために速度分布関数に基づいたエントロゼーの詩間発展を調べた結果、適当な数の欠賂粒子
を導入することで、平衡状態に至るプロセスを確認できた。また欠賂粒子の数は、今のモデルのサ
イズの場合 8個から 10傍の欠陥粒子が必要であることがわかった。抵速衝突における衝突速度
とはねかえり係数の関係は、 2次元の準静理論と良好な一致を見せるが、内部散逸の時間スケー
ルをフィッティングパラメータとして導入することが必要であり、その留の起源はまだわかってい
ない。
5 斜め衝突のシミュレーション
この章では、筆者らが行った弾a往体の斜め衝突のシミュレーションとそれによって得られた衝
.
7章で、球と壁の斜め衝突にお
突角度とはねかえり孫数の関係についての解析結果を紹介する。 2
いてはねかえり採数が箆突角度によって系統的に変北し、更に l以上の値を取り得ることを示し
3
1
]、ここではこの現象を 2次元の数誼シミュレー
たコーネル大のグ、ループの実験を詔介したが [
ションによって再現し、衝突角度とはねかえり係数の関係を決定するメカニズムについて考察す
る[
3
3，
34，
8
1
]。
-7
0
3ー
醤f
中寛人
5
.
1 数値モデル
筆者らの用いた弾性体モデ〉しは質点とばねから構成される 2次元のランダム稽子モテ、ルで、ある
[
3
3，
8
1
]。弾性体モデルは半径 Rの円盤と高さ 2R、i
福 8Rの板から構成され、それらを構成する質
し
図1
1
: ランダム格子による円盤と板の種突モデ、ル。 γを衝突角度とする。
点はそれぞれ同じ密度に在るようにそれぞれ 800個
、 4000穏をランダムに配置した(図 1
1
)。質点
はデローニ一三角分割のアルゴリズムを用いて近接の 2点同士を、次のようなポテンシャルで表
8
3
}
:
される非線形ばねで接続する [
1_
1_
V(X)= ;;kα
;
k
b
X'i.
2 ~ X~ +
.:
4
<)
A
(
2
9
)
ここで zは岳然長からの伸びである。ばね定数んは円盤、板に対して各々 k
aニ 1
.
0x mc2/R2、
k
b=1
.
0X 10-2mc2/R という値を用いた。またんに関しては円盤、板両方に対してん二九 /100
という値を用いている (mは質点の質量、 cは 1次元の音速)。壁の南場と底は国定境界である。ま
た、このモデルのポアソン比は機械的に板を牽引して垂室及び平行方向の変泣を測定することに
7
.
5
0士 0
.
1
1
)X 10-2である 100
より計算することができ、その誼は ν =(
シミュレーシヨンは初期速度を I
.
l
cに固定して、衝突角度 γの正接を tanγ=1から
v
(
O
)
1= O
tanγ=7まで変化させて、それぞれの角震に対するはねかえり祭数を測定した。また、質点、の配
置による結果のばらつきを取り除くために、質点の配置を変えた 100種類の円盤を用意し、各衝
突角度毎に 1
00偲のデータを平均した。数誼計算は各質点の運動方程式を 4次のシンプレクティッ
ク積分で時間刻み d
t=10-3R/cで解いた。
5
.
2 結果と考察
図1
2は、シミュレーションによって得られた島=t
a
n
γ とεの関係である。×印は各衝突角度
における 1
0
0種類の円盤によるデータの平均値、エラーパーはその標準偏差老表している。 eが I
10
弾笠率の計算の蒔だけ、費点詞に非線形ばおの相互作用 i
こ加え、ばね両端の梧対速度に比例する摩擦力の棺互作舟
を導入して振動緩和後の変位を測定する。
門
i
dq
ハV
巨視的物体の非弾性衝突
1
.4
1
.3
1
.2
e
き
Suoul
;
al
I9n-持4
T
h
o
句
、
，
一
一
ー
す一ーさ
哩
"
図 1
2
:弘とはねかえり係数 εの関係。×印はシミュレーション結果、エラーパーは標準偏差。実
線は式 (
3
0
)で、九二 0
.
9
5の時の ε をプロットしたものである。
L
t
x
図 1
3
:円盤がやわらかい壁に衝突している様子。
.
0あたりでピークが出る傾向があるという遣いがあるものの、 eが
を越え始める角度や、言 1二 6
t
a
nl'の増加につれて上昇し 1を越え始めるというコーネル大のグループの実験と似た傾向が確認
された。
コーネル大のルージ、ュとアダムスは彼らの実験結果を説明するために、壁表面に生じた局所変
形が球の軌道を上向きに変え、その効果によりはねかえり集数が増加し 1を超えるようになるの
ではないかと推測した [
3
1
]110 そこで衝突で壁表面に生じる局所変形が誌ねかえり孫数に与える
影響を考えてみよう。
図1
3，こ、円撃がやわらかい壁に衝突し、壁表面に最大の局所変形が生じている様子を示した。
壁表面の法線方自の単位ベクトル nが、衝突によって生じた局所変形によって円盤に向かつて角
度 αだけ回転し、接触面上では nO になるとする。ここで、簡突時毘を
Tc
として、積突畿の速
α で定義した 2つのはねかえり孫数、 e二一 (
V
(
九).
n
)
/
(
v
(
O
).n)と
度を V(T
c) として、 n と n
九二一 (
V
(
T
α
)
/
(
V
(
O
).
nO ) を導入する。この 2つのはねかえり係数の関孫は、
c).n
t
a
n
α
e-eα 十 世2
1ー を l
t
a
n
α
(
3
0
)
はそれぞれ¥Ti (
を1- t
a
nα
)
/
(
1+を 1t
a
nα
) とを2' = (
骨 1- tanα)/(1十
となる。を?と守2
実際、壁の局所変形が物体の反射角を変え、はねかえり係数に影響するという実験結果がある 1
8
4
1
0
11
- 705-
題仲寛人
図 1
4
: 円盤と壁表面の模式図。
Wltanα)-3
(
1十 九 )(μ+tanα)/(1-μtanα〉で表される。ここでを2
の表式を導くときにウォル
3
8，
3
9
]
0t
a
nαを求めるためには、丹盤の表面を何
トンらの鰯体円盤の斜め衝突理論を用いている [
らかの関数形で表現する必要があるが、ここでは板が円盤よりもずっと軟らかいので剛体円盤が軟
)2/2R-ycという二次
らかい表面と摩擦力を受けて接触していると仮定し、接触面を g
(
X
)=(x-xc
(
g
(
X
b
)-g(x
α
))
/
I
X
bーら￨
1-2
亡
と定義すれば、弾性論の教科書 [
8
5
]によると αの正接は、 t
a
nα = I
xcーら￨
及び三二
R(2-2
三)
(
1
/
π
)arctan(l-2
v
)
/
(μ(2-2
v
)
)で表すことができるいは表面の摩擦集数 )
0I
xc
-x
a
lの舘及び各
bとおいて t
a
n
α
関数で表すことにする。接触面の南端の z産標をお、 X
=
I
x
-x
a
lは 0.55Rという一定値
c
衝突角度に対する μ の僅は、シミュレーションの結果から求め、
(
O
)
)で定義し、水平及び垂産或分老それぞれ Jt、み
を用いた。また μは衝撃を J= m(v(Tc}- v
とすると、 μ=I
J
t
l
/
I
J
n
lで定義される。図 15の×印はこの定義によって計算した笥突角震と μの
関罷である。これから各弘における tanαを計算できることになり、結果的にはねかえり係数と
禽突角度の関孫を得ることができる。
図1
2の実隷は式 (
3
0
)を εα=0
.
9
5という鍾を用いて措いたものである。式 (
3
0
)はシミュレー
ション結果をうまくフィットできることがわかる。このように斜め衝突のはねかえり長数には、衝
突角度と摩擦係数の関係が影響する。
最後に笥突角度と摩擦孫数の関係を簡単なモデルで考察する。壁表面にラフネスに相当する同
じ大きさの突起が密愛 ρで等間爾に並んでおり、衝突時に円盤が突起に接触すると円盤の持って
いるエネんギーの一部が突起に散逸するとすれば、接触した突起の数で衝突後の速度が決定され
4
)。ここでパラメータ ηを導入して、 4番目の突起に接触する詩の接線方向の速度イ)は
る(図 1
♂ =(1 ー η}V~i-l)
で表されると仮定しよう。もし積突時間中;こ N 偲の突起と接触するならば、
衝突後の接線方向の速度はむt
(
T
f
l=(
1-η)
N
V
t
(
O
)となる。 ηは以下の通りに決める。円盤
C) = v
の接線方向の衝突速度が十分に大きい時には表面の突起が押しつぶされるという効果により、 ηは
小さくなると仮定する。そこで、 ηの関数形として、
。
η G η
7
η=1÷βv;(O)/V~(O) -1工
否
盃
(
3
1
)
という形を仮定する。新たに導入した無次元パラメータ η。と β は後にブイツティングパラメータ
として扱う。
更に N を決定するために、円盤と壁の衝突時間中に壁表面をすべる長さおを見積もる。筈
- 706-
三提的物体の非弾性蕎突
0.065
Sim
叫a
t
i
o
n 伊持4
0.06
Tb帥 守 一 一
0.055
0.05
μ0.045
0.04
0.035
0.03
O.025
O
1
2
3
4
5 6 7
'
P1
図1
5
:Wlと摩擦係数 μの関係。
7
4
] にほぼ等しいと仮定すると、 l
s
lは
、
突時間九は、円盤が壁に正面衝突したときの衝突時間 [
[
s
l
二
V
t
(
O
)'Tc= π(
R
/
c
)ゾl
n[
4
c
/
v
0
)
] となる。よって接線方向の衝撃 Jt三 J.tは
、
n(
Jt=M{(1-η
)〆s
t
}吟 (
0
)三 一 ηMρl
s
l
V
t
(
O
)
(
3
2
)
と計算される。ここで ηは小さい重であると仮定した。また法線方向の衝撃 Jn 三 3・n はL
まねか
)
v
0
)と計算される。
えり係数 εを用いて Jn = -M(e十 1
n(
!Jn に従って摩擦捺数 μを計算すれば、
さて法線方向及び接線方向の衝撃から定義 μ= Jt
需品計占f 戸両)呼l
μ=
l
n(40
(
3
3
)
内 /(1十 時 と い う 関 孫 を 用 い
という形はる。ここで、 ωγ=ψ/(1十 時 、 及 び sinγ=
た。この式において、はねかえり孫数 eはシミュレーションの結果(図 1
2の×印)、 I
v
(
O
)
I
/
cは
。
0
.
1という値を用いる。フィッティングパラメータとして、 β =0
.
2
1、η ρR= 0
.
1
8という舘を用い
ると、図 15の実線を措くことができる。このように我々の単純北した摩擦のモデルでシミュレー
ション結果を再現することができる。
以上、この章で紹介した筆者らの研究結果は以下のようにまとめられる。
1.はねかえり係数が衝突角度によって上昇し、 Iをなめらかに越えるためには壁の局所的な変
形が必要であり、その変形量を理論的に克積もることによって、衝突角震とはねかえり係数
の関係を得ることができた。
2
.笥突角度とはねかえり係数の関係には壁表面の摩擦係数が関接しており、衝突角度と摩擦捺
数の関僚も簡単なモデルで説明できる。
特にはねかえり係数と衝突角度の関係に関しては、静弾性理論に基づいた解析を行なったが、壁
表面に生むる局所変形の時間的変色が円盤の反発挙動にどう影響するのかを理論的に考察するの
は興味深い試みである。
t
司
ウ
ハtU
園停寛人
6 小さな系の 3次元衝突
ここまでの議論で、巨視的な物体の衝突現象に対する先行研究の紹介と 2次元の弾性体モテ、ル
を用いた筆者らの解析を紹介したが、この章では微視的な系の衝突現象を扱う。
数留から数百留程度の原子や分子の集合体はクラスターと呼ばれ、その液体とも面体とも異な
る中関的な物質の物理的・化学的性質が、実験や数種シミュレーションで調べられている。例え
ば、ナノスケール程変のサイズを持つクラスター同土の凝集やシリコン等の基板上への蒸着は、ナ
ノスケールの電子デ、バイスの開発過程においても重要なプロセスで、あり、分子動力学法等の数種
的手法を用いてそれらの衝突ダイナミクスが調べられている [
8
6
]
0
ナノスケール程度のサイズを持つ微小な系の衝突現象は、巨視的な系の衝突現象と違って、そ
のダイナミクスが熱揺らぎに大きく影響されるという特徴がある。また分子関力等の相互作用に
起因するクラスター同士の吸着が起こりやすくなるのもその特徴のーっと言える。そのような系
が示す衝突現象において、自由度胃のエネルギーの受け渡しがどのように変化し、巨視的な物質
突現象と違った笥突モードが出現するかということを探ることは、物理学の見地から見ても
の種1
面白い問題だと思われる。
一般的な物体同士の笥突においては、衝突前の運動エネルギーの一部が衝突時間中に物体の内
部自由度に分配され、一旦内部自由度に分配されたエネルギーは物体の並進自由度には戻ってこ
ない。このような衝突の事を私達は非弾性衝突と呼んでいる。非弾性衝突を、並進自由度からそ
れ以外の自由度への不可逆的なエネルギー移動現象で、あると見なせば、非弾性衝突は熱力学第二
法則と密接に関保していることがわかる。混度 T を持つ 2つの物体を相対速度 V で互いに正面衝
突させる状況を考えよう。この時、熱力学第二法則は次のように書くことができるだろう。
トv
2十 O
ここで μ は 2物体の換算震量できろり、
陀
μ
jV'2
(
3
4
)
V'は衝突後の棺対速度である。巨視的な物体の衝突の場
3
4
) の左辺第 2項の彰響は小さく、熱力学第二法則という制約があるために衝突後の速
合
、 式 (
度は衝突前の速度よりも必ず小さくなる。したがって巨視的な物体の反発孫数を計挺すれば必ず
1よりも小さな値を取ることになる。
しかし、もし物体のサイズがナノスケールで、熱揺らぎがダイナミクスに大きく影響するよう
になると、式 (
3
4
) の左辺第 2項の影響は無視できなくなるだろう.その時，熱力学第二法期は薙
率的な制約になる。言い換えると、内部自由度から並進自由度へのエネルギーのフィードバック
異常反発j
が有限の確率で起りうるのであり，その時にはねかえり採数が 1以上の値そ取る様な f
とでも呼べるような現象が観測されることになる。
筆者らは分子動力学法に基づくナノクラスターの衝突シミュレーションを行い、衝突速度や温
度等の初期条件やクラスター間の相互作用によって、反発を特徴付ける量がどのように変わるの
かを謁べている。また近年、熱揺らぎが支配的となるような小さな系の援る舞いにおいて、負の
エントロピ一生或率が出現する確率を定量する f
揺らぎの定理j が非平衡物理における重要な或
- 708-
三視的物体の非弾性衝突
果のーっとなっているが [
8
7
トこの章では非弾性衝突における揺らぎの定理と我々のシミュレー
8
8
]
0
ション結果との比較についても紹介したい [
6
.
1 モデル
ここで我々の笥突モデ、ルを紹介する。我々の用いたシミュレーションモデルは菌 1
6に示すよう
なレナード・ジョーンズ粒子から講成される、同質な 2つのクラスターである。座標軸は 2つの
クラスターの重心を結ぶ方向を Z軸とする。 1つのクラスターを構成する 682留の f
原子j は面
原子 j 同士はレナードージョーンズポテンシャル、
心立方の初期配置を取り、 f
川
士
主
[
(
三
)(
三
)
"
]
(
3
5
)
12_
で結合している。ここで、
Tijはクラスターを構成する粒子
i
，
j間の距離、けま結合エネルギー、 σ
は原子直径を表している。
2つのクラスター関の相互作用は以下のようにモデ、ル化する。上部クラスター CU の下表面上の
〆
可222222sg曇
SJ
σ 一印
C
f''SE、
t
‘
e霊孟
¥
、211
ノ
7u
σ 一山
/
I
z
-、
ξ
﹁?lSEtt-ft﹄
4斗&
L九
Ef
J'e
T
R、
‘
，
一
一
、
相互作用が偉くとする。
。
戸¥E13/
粒子 kと下部クラスター clの上表面上の粒子 lの間には、次のような粒子関距離 Tkl に故存する
(
3
6
)
すなわち、通常のレナード・ジョーンズ相互作用に吸着パラメータ C を導入し、クラスター間に
とした。
働く吸着をコントロールできるようにクラスター間相互作用をモテ、ルf
このモデ、Iしを用いた基本的なシミュレーションのセットアップは次の通ちである。まず 2つのク
面に対して鏡面対称 i
こ配置する。そしてクラスターを構成する各粒子に、マッ
ラスターを z=o
クスウェルーボルツマン分脊に従う乱数による初期速度を与え，その分散に穏当する初期温度を与
8
9，9
0
]で熟平衡化した後に、並進速度を与え二つのクラスタを
える。そして速度スケーリング法 [
互いに正面衝突させる。 2つのクラスターが反発して並進速度が一定になった時に、衝突謡後の
クラスター重心の栢対速度の z成分の比から反発係数を計算する。
t
1
.0 x
計算はクラスターの構或粒子の運動方程式を、速度ベルレ法を用いて時間刻み d
2
1
0
2.
J
a m/Eで解いた。また計算時間の短結のため、レナード・ジョーンズ棺互作用のカット
オフを σ
c= 2
.
5
σ で入れた。データのサンプル平均は、初期温度として与えるマックスウェルーボ
ルツマン分布に従う、全粒子の初期速度の組を変えることで複数のサンプルを取り、その平均を
計算することで取った。長さ、質量、及び時間はそれぞれ σ、ー粒子の質量 m、及び
p 可;で
無次元化し、パラメータの舗はアルゴンの値と比較して 12 シミュレーションを行った。
6
.
2 シミュレーション結果
こ乙では話章で導入したクラスターモデルによる積突シミュレーションの結果を紹介するが。ま
3
6
)において、 c=O.Oとした持〈つまり斥力相互作用の
ずはクラスター関の相互作用である式 (
12アルゴンの場合
σニ 3.
4x]0一
勺m，] m = 6.69x1O-26Ikg]，正=1
.65X 1
O-21IJJ.
- 709-
題
イ
中
寛人
Z
U
C
。
cl
図1
6
: レナード・ジョーンス粒子から構成されるクラスターモデル。各クラスターは 682掴の構
成粒子から成る。
1
.
1
5
T=O.Ol~骨4
T=O.02+
ー
e 0.95
0.85
0
.
8
0.75
v/(dm)
国
C =0
.
0，
初期温度 T= O
.
O
l
f
:，及び T = 0
.
0
2
εの
図 17: クラスターの衝突速度と反発係数の関孫 •
場合.
みが働く時〉の結果を示す。
6
.
2
.
1 斥力相互作用のみの場合の衝突
図 17に示すのは禽突速度と反発孫数の関係である。反発係数は衝突前後の梧対速度の Z成分の
比を用いて、
e=ーに/にで定義する。また図 17には、クラスターの初期温震として T=0.01E
とT = 0
.
0
2
fの二つの場合を示した。サンプJ
}
/平均は、各速度に対して粒子に湿震として与える
拐期速度の組を 100還り変化させ、 100サンプルのデータを平均した.これよりクラスター簡に斥
力のみの相互作詞が動く時に、はねかえり係数は衝突速度に対して単調に減少する傾向が晃られ
ることがわかる。これは巨視的な物体同士の衝突でも見られる額向である。
2
S
Jの結
シミュレーションで得られた衝突速度とはねかえり保数の関係老、低速衝突の準静理論 [
果と比較した。これは巨視的な物体同士が互いに正面衝突した時に，物体が倭突時に受ける力老、
H
e
r
t
z理論に基づく弾性力と、内部粘性に由来する散逸力の和で表せると仮定し、それによって物
体に生じる歪の詩開発展から，衝突後の速度を求めるものである。それによると、衝突速度 V と
宅2ム
t
ウ
ハV
巨梶的物体の非弾性衝突
1
.
1
1
0
.
9
11 I
.
J
.
.
I
(
"
S
-E
0
.
7H!I
0
.
6
t. Simwation~持
宝もe
o
r
y0
0
.
5
o
i
0
.
2 0.
4
0
.
6
V
/
(
E
1
m)
0
.
8
1
]
J
冨
国 1
8
: クラスタの第突速度と反発係数の関採。 C ニ 0
.
2、初期温度 T= 0
.
0
2
Eの場合。
反発長数 eの関係は、
(町市)1/5
ε2α1一 α2
(
3
7
)
と表すことができる。ここで α1は通堂 αlニ Iであり、 α2は物性信で決まる定数である。ここで
は
， α1)α2をフィッティングパラメータとしてデータのフィットを行った。
図 17の実線と毅線は、それぞれ (
α1
，α
2)= (
1
.744，
1
.5
7
2
)，
(
1
.323，
0
.
8
6
9
)で式 (
3
7
)をプロットし
3
7
)によるデータのフィットはよいが、 α1の値は 1より
たものである。これからわかるように式 (
も大きな舘になる.これは低速の衝突;こおいて、反発採数の平均植が l以上になることを表してい
る。実際に拐期温度 T=0
.
0
2
Eの時、衝突速度 V=0
.
0
2
.
;
;
]
石の時のはねかえり孫数の平均値は
1
.04となった。
6
.
2
.
2 引力相互f
乍用があるの場合の衝突
.
0となる時に衝突速度とはねかえり祭数の関係はどうなるであ
で は 、 式 む の に お い て C >0
ろうか?我々は C =0
.
2
[
8
6
]という弱いクラスター間引力が働く場合について同様のシミュレー
シヨンを行い、その関係を探った [
9
1
]
0
図1
8は初期温度 T=0
.
0
2
ε の場合について行ったシミュレーションの結果である。乙れよりク
ラスター関に引力が働く場合には、はねかえり保数と衝突速度の関係はある笥突速度においてピー
クが晃られることがわかる。これは衝突速度が小さい時には、簡突詩において凝着力の影響が大
きいために反発時の速度が大きく減少するが、衝突速度が大きくなるにつれて凝着力によるエネ
ルギーよりも初期の運動エネルギーが勝るようになるためである。
図1
8内の自丸は、ブリリアントフ (N.V
.B
r
i
l
l
i
a
n
t
o
v
) らによる引力相互作用を考慮した衝突
9
2
]
0 これは巨視的な物体同土の衝突に関して、衝突中の物体に Hertzの
理論による結果である [
o
h
n
s
o
n
-K
e
n
d
a
l
l
-R
o
b
e
r
t
s
)理論
接触理論に基づく弾笠力と内部粘性由来の散逸力に加えて、 JKR(J
[
9
3
Jに基づいた凝着力が働くとして、やはり変形の時間発展 13からはねかえり係数を見積もるとい
うものである。この持、箆突物体の表面張力と粘性係数が必要になるが、表面張力はクラスター
実際には接触半窪の時間変化を求め，そこから全体の変形を計算する，
13
t
ウ
圏持寛人
間の引力相互作用から晃積もった γ~ 0
.026E/(m
りという誼を用い [
9
4
ト粘性係数拡フィッテイン
グパラメータとした。
8内の白丸は粘性係数をフィッテイングパラメータとしてフィットした結果であるが、特に
図1
種突速度が大きくなるとはねかえり係数の第突速度故存性はかなり小さくなり、データとの一致
はよくない。 表面張力の見積もりや、系のサイズによる表面の多角形の彰状の効果などがその原
国として考えられるが、その彰警を探る事は今後の課題である。
6
.
3 揺らぎの定理との比較
ここで、 c=O.Oの場合に慢して反発係数の出現確率について定量的な議論を試みる。熱揺らぎ
の影響が無視できないような小さな系においては、有援の確率でエントロビーが減少する現象が
観測されることが知られており、それを定量するのがゆらぎの定理 [
8
7
}である。これを非弾性笹
9
5
]
0
突に適用したものは、次のようになる [
expβ
(W)P(X
)= P(X1 → Xo
)
.
o→ Xl
(
3
8
)
ここで、 βは逆逼度を表している。 Xo
、X1はそれぞれ重心自由度の初期状態と終状態を表してお
り
、 P(X
o→ X1)は初期状態 Xoから終状態 X1~こ遷移する確率分布である。また、 P(X 1 → X o )
は、衝突を時間反転した時に、終状態から拐期状態に戻る確率分布を表す。 W は衝突により重心
自由度が失う(または内部自由疫が獲得する)エネルギーを表し、
W =
4
1
V
fV
戸
内
内
究
勺
(
ο
1
2
一
の
ε
J
り
2
)
β
でで、表される。このように W は 1- e2 の関数で表すことができるので、 1- e2 の出現確率分布
2
) とを求めて、式 (
P(1-e2) と、律j
突を時間反転した時の P(1- e
3
8
)の成立を確かめる。
.
0
2
、
正 T =0
.
0
3
Eのそれぞ、れの場合について行った。衝突速
シミュレーションは初期温度 T ニ 0
.
0
2
(正
/m)1/2で衝突させ、初期の高さにクラスタが戻ってきたら eを測定し、もう一度熱
度 V= 0
手傷f
とした後に重心速度を反転させて、反発係数 tを総定する q これを 1サンプルとし、初期の
各粒子に与えるマックスウェルーボルツマン分布に従う初期速度のセット老 3000サンプル馬意し、
π/8，
π/4と変え、計 9000サンプルのデータを平均した。 国 1
9
、20にそれ
更に初期配向を()=0，
ぞれ初期誼度 T =0
.
0
2
、
正 T =0
.
0
3
正の時の結果を示す。横軸に 1-e2、縦軸に InF/Pを取った。
実線は各温度に対して、 1-e2 と W の関孫をプロットしたものである。これを晃ると南端の領域
では一致は悪いものの、中間領域の一致は悪くないことがわかる。ただしシミュレーション初期
の速震スケーリングによる熟平衡化の際にクラスターに励起される回転が反発に大きく影響する
ため、データの車線からのずれに影響することがわかっている。また、式 (
3
9
) の W の計算にお
いて、重心自由度のエネルギーに加えて重心まわりの回転自由度のエネルギーとクラスター間の
梧互作用のエネルギーの項を入れて式 (
3
8
) の或立を確かめた [
8
8
]
0 この時は国 1
9や 20~こ示し
たデータよりは車線的なデータが得られたが、式 (
3
8
) との一致はあまりよくはなかった。この
事は用いた系において重心自由度と内部自由度がうまく分離できていない事に起麗している可能
性もあり、サイズ依存性や粒子間の梧互作用依存性等を含めて更に調べる必要がある。
-7
1
2-
巨視的物 f
本の非弾性萄突
5432
0
.
5
o
:
s
0
1
唱 2
(a)
密 1
9
: ゆらぎの定理との比較。
T= 0.02正の場合。
aaz
oqd司
冷 dκx
s
毛滋u
l
a
t
i佃
X
FT 一一
=
r
一 一 てδ
.
5
0
le2
(
b
)
密 2
0
: ゆらぎの定理との比較。 T=0
.
0
3正の場合。
7 議論と今後の展望
本稿では、衝突研究の歴史的経過と筆者らの近年の研究成果を紹介した。第 4章、第 5章では
第 3章で導入したばねー質点モデ、ルに基づく 2次元弾性体モデルを用いた衝突シミュレーション
の結果を紹介したが、本稿のモデルは 2次元モテ、1しのため 3次元の実験結果とは定量的な比較が
できないという欠点がある。物体の次元が上がることによって、エネルギー散逸の蔀合が大きく
変わることは容霧に想像でき、実験結果との比較を試みるのであれば、 3次元のモテ、ルを構成した
上でシミュレーションを行う乙とが必要である。第 3章で紹介した拡張個別要素法や粒子法、有
限要素法等に基づいた 3次元の弾性積突シミュレーションを行うことが必要であろう。ただし粒
子ベースのモデリングで 3次元物体を穣築すると、 2次元よりも粒子数が増え、計算により長い時
間を要するという問題があるため、それを克服するには、計算上の工夫が4必要になる。
第 4章、第 5章における我々のシミュレーションにおいて誌、初期回転がない状態で笥突速度
や衝突角度のはおかえり孫数への依存笠を探った。だが、第 l章でも触れたように初期回転が反
発に与える影響老謂べることは興味深い。初期の回転自由度のエネルギーと並進自由度のエネル
ギーの割合が法線方向や接線方向のはねかえり係数にどう影響するのか、また自由度開のエネル
門
qU
1i
i
題件費人
ギー遷移はどうなっているのかということを実験やシミュレーションで探ることは面白い試みで
るると思う。
また、これまで均質な物質の衝突のみを考えてきたが、内部構造のある物体の倭突現象に対す
る理論的な考察も面白い研究テーマとなるだろう。 2003年 8月、米メジャーリーグの試合におい
て当時シカゴ・カブスのサミー・ソーサは
zにコルクをつめた違法パットを試合で用いて退場処
分を受けたが、コルクをつめたパットは選手の閣ではボーんを遠くへ飛ばすパットだと信じられ
ており、それが何故か?という問いかけはスポーツ科学や物理学の研究者の関心を惹いた [
9
6
]14
0
コルクパットやゴルフのツーピースボールのような層状の構造を持つ物体が衝突した時に、内部
に生じる弾性波動の伝わり方が箆突後の運動にどう影響するかという問題は、弾性学のテーマと
しても面白い問題であると思う。
第 8章では 2つの同質なレナード・ジョーンズクラスターの衝突シミュレーションを行い、ク
ラスター間の相互作用が、クラスターの衝突速震とはねかえり係数の関係に及ぼす影響について
紹介した。更にクラスター関に斥力梧互作用のみが働いたときの反発露数の出現確率に基づいて、
非弾性衝突のゆらぎの定理とシミュレーション結果との比較を試みた。クラスター閣の相互作用
が斥力相互作用のみの場合、母速の衝突速度領域において反発係数が 1以上になるすンプルを観
l
:
'、そのような衝突が起こったときにクラスター内部のストレス場の変化
浪せすることがわかった t
やどのような振動モードが励起されているのかを微視的に謁べることは今後の課題である。また
凝着力が働く場合の衝突に関して、第 6章では
る笥突モードが，
c=0.2の場合の結果のみを報告した。観測しう
cが様々な値を取ったときにどのように変化するかを調べ [91]，更にその変化が
衝突速度や温度にどのように依存するかを謁べる事は興味深く、更なるシミュレーションによる
解析も今後の課題として残されている。
謝辞
本稿に多数の有益なコメントをしていただいた共同研究者の早川尚男教授〈京都大学基礎物理
学研究所〉と、野坂京子氏〈物性研究刊行会〉に諜意を表したい。
参考文献
[
1
]安田徳太郎訳・編.新訳ダンネマン大白黙科学史 1.三省堂，東京， 1
9
7
8
.
1
2
1W.Goldsrnith. Impact: The Theo叩
andP
h
y
s
i
c
a
lBehaviorofC
o
l
l
i
d
i
n
gS
o
l
i
d
s
. Edward
ArnoldPub
，
.
l London，1
9
6
0
.
∞
.S
t
r
o
n
g
e
. ImpactM
e
c
h
a
n
i
c
s
.C
a
r
n
b
r
i
d
g
eU
n
i
v
e
r
s
i
t
yP
r
e
s
s，
C
a
r
n
b
r
i
d
g
e，
2 O
.
[
3
] W.J
1
4
]SirH.マッセイ.原子・分子の衝突.共立出張，東京，
1981
.
14様々な理由討けが考えられているが、単純にパットの先端部分に異物を詰める事によってパット全誌の重さが軽く
なり、旦つ重心の位置がグリップの方向に移動するため慣笠モーメントが小さくなって、パットを撮ったときにヘッド
スピードが大きくなるからという説明もある [
97
1
0
t
ウ
a斗.
E梗的物体の非弾性衝突
[
5
]I
司村定矩他編.人類の住む宇宙〈シリーズ現代の天文学1).日本評論社?東京， 2
0
0
7
.
[
6
]L
.P
.K
a
d
a
n
o
f
f
.B
u
i
l
tupons
a
n
d
: Te
郎 防i
c
a
li
d
e
a
si
n
s
p
i
r
e
dbyg
r
a
n
u
l
a
rf
t
o
w
s
.R
e
v
.Mod.
Ph
ジS
.，
71:435-444，1
9
9
9
.
[
7
]P
.G.deGennes.Granularm
a
t
t
e
r
:at
e
n
t
a
t
i
v
ev
i
e
w
.R
e
v
.Mod.P
h
y
s
.，
71
:8374-8382，
1
9
9
9
.
[
8
] H.M.Jaeger，8
.R.Nagel，andR.P
.B
e
h
r
i
n
g
e
r
. Granulars
o
l
i
d
s，l
i
司u
i
d
s，andg
a
s
e
s
.R
e
v
.
9
9
6
.
Mod.P
h
y
s
.，68:1259-1273，1
9
9
5
.
[
9
J 田口善弘.砂時計の七不思議-絵粒体の動力学.中央公論社?東京， 1
[
1
0
] 早川尚男・那須野信.現代物理最前線 1 (大様義彦編) .共立出販，東京， 2
0
0
0
.
[
1
1
] 早川尚男.散逸粒子系の力学.岩波書屈?東京， 2
0
0
3
.
[
1
2
)J
.デュラン(中西秀・奥村廠訳) .粉粒体の物理学.吉岡書庖，京都， 2
0
0
2
.
[
1
3
]P
. B. Umbanhowar，F
. Melo，and H. L
. Swinney. L
o
c
a
l
i
z
e
de
x
c
i
t
a
t
i
o
n
si
nav
e
r
t
i
c
a
l
l
y
u
r
e，382:793-796，1
9
9
6
.
v
ib
r
a
t
e
dg
r
a
n
u
l
a
rl
a
y
e
r
. Nαt
[
1
4
] A.Daerrand8
.Douady. Twot
y
p
e
so
fa
v
a
l
a
n
c
h
ebehaviouri
ng
r
a
n
u
l
a
rmedia. Nature，
399:241-243，1
9
9
9
.
fg
r
a
n
u
l
a
rp
a
r
t
i
c
l
e
s
[
1
5
] O.Moriyama，N.Kuroiwa，M.Matsushita，andH.Hayakawa.4/3lawo
h
y
s
.R
e
v
.L
e
t
t
.，
80:2833-2836，
1
9
9
8
.
f
t
o
w
i
n
gthroughav
e
r
t
i
c
a
lp
i
p
e
.P
[
1
6
] Y.Yamazaki，S
.Tateda
，A.Aw
詔 u
，T.Arai，O.Moriyama，andM.Matsl
凶 i
t
a
.Dynamics
.P
h
y
s
.S
o
c
.J
p
n
.，
o
fg
r
a
n
u
l
a
rf
l
o
wthroughav
e
r
t
i
c
a
lp
i
p
eundert
h
ec
o
n
t
r
o
lo
fi
n
f
t
o
wr
a
t
e
.J
71
:2859-2862，
2
0
0
2
.
9
9
5
.
[
1
7
] 編集委員会編.物理学大辞典.丸善株式会社，東京， 1
.Johnson. C
o
n
t
a
c
tMech
似 c
s
. CambridgeU
n
i
v
e
r
s
i
t
yP
r
e
s
s，
Cambridge，1
9
8
5
.
[
1
8
] K.L
[
1
9
]E
.Hodgikinso
立 Ont
h
ec
o
l
l
i
s
i
o
no
fimper
長c
t
l
y
e
l訓 i
cb
o
d
i
e
s
. B九t
.A
s
s
o
c
.Report
，4:534543，1
8
3
4
.
.H.V
i
n
c
e
n
t
. Experimentsoni
m
p
a
c
t
.P
r
o
c
.Cαmb
九d
g
e
，Phil
.S
o
c
.，1
0
:
3
3
2
3
5
7，1
9
0
0
.
[
2
0
]J
[
2
1
] C.V.Ra
man.Thephotographicstudyo
fimpacta
tminimalv
e
l
o
c
i
t
i
倍
.
P
h
y
s
.R
e
v
.，1
2
:
4
4
2
-
448，1
9
1
8
.
.G.B
r
i
d
g
e
s，A.Hatzes，andD.N.C.L
i
n
.S
t
r
u
c
t
u
r
e，s
t
a
b
i
l
i
t
yande
v
o
l
u
t
i
o
no
fs
a
t
u
r
n
'
s
[
2
2
]F
r
i
n
g
s
. Na
t
u
r
e
，309:333-335，1
9
8
4
.
i
司
Fhu
題伸寛人
[
2
3
] K.D.Supulver，
F
.G.B
r
i
d
g
e
s，
andD.N.C
.L
i
n
.Thecoe缶c
i
e
n
to
fr
e
s
t
i
t
u
t
i
o
no
fi
c
ep
a
r
t
i
c
l
e
s
i
ng
l
a
n
c
i
n
gc
o
l
l
i
s
i
o
n
s
: Experimentalr
e
s
u
l
t
sf
o
ru
n
f
r
o
s
t
e
ds
u
r
f
a
c
e
s
. lCARUS
，113:188-199，
1
9
9
5
.
[
2
4
] R.Sondergaard，K.Chaney
，andC.E
.Brennen. Measurementso
fs
o
l
i
ds
p
h
e
r
e
sbouncing
h
eASME
，J
.App
.
lMech
，
57:694-699，1
9
9
0
.
o
f
ff
t
a
tp
l
a
t
e
s
. Tmns.oft
nac
o
l
l
i
s
i
o
nbetweentwos
p
h
e
r
e
s
. Jpn
[
2
5
] G.KuwabaraandK.Kono.R
e
s
t
i
t
u
t
i
o
nc
o
e
f
f
i
c
i
e凶 i
J
.Ap〆
p.
l昂
Ph
匂y
μ
s
.，
お
26
:
ゴ
1
幻
23
0
ι
一ロ
1
2お
33，1
9
8
7
.
[
2
6
] N.V.B
r
i
l
l
i
a
n
t
o
v，
F
.Spahn，
J
.
M
.H
e
r
t
z
s
c
h，andT.P白 c
h
el
.Modelf
o
rc
o
l
l
i
s
i
o
n
si
ng
r
a
n
u
l
a
r
h
y
s
.R
e
v
.E，
53:5382-5392 1
9
9
6
.
g
a
s
e
s
.P
ヲ
[
2
7
] W.A.M.Morgadoand1
.Oppenheim.Energyd
i
s
s
i
p
a
t
i
o
nf
o
rq
u
a
s
i
e
l
a
s
t
i
cg
r
a
n
u
l
a
rp
a
r
t
i
c
l
e
e
v
.E，5
5:1940-1945，1
9
9
7
.
c
o
l
l
i
s
i
o
n
s
. Phys.R
. Coe缶 c
i
e
n
to
fnormal r
e
s
t
i
t
u
t
i
o
no
fv
i
s
c
o
u
sp
a
r
t
i
c
l
e
sand
[
2
8
] T. SchwagerandT. Poschel
h
y
s
.R
e
v
.E，
57:650-654，1
9
9
8
.
c
o
o
l
i
n
gr
a
t
eo
fg
r
a
n
u
l
a
rg
a
s
e
s
.P
[
2
9
] R.Ramirez，T
.Pるs
c
h
e
I，N
.V.B
r
i
l
l
i
a
n
t
o
v，andT.S
c
h
w
a
g
e
r
. Coe詰 c
i
e
n
to
fr
e
s
t
i
t
u
t
i
o
no
f
h
y
s
.R
e
v
.E，
60:4465-4472，1
9
9
9
.
c
o
o
l
i
n
gv
i
s
c
o
e
l
a
s
t
i
cs
p
h
e
r
e
s
.P
(
3
0
]F
.G
e
r
landA.Z
i
p
p
e
l
i
u
s
.C
o
e
f
f
i
c
i
e
n
to
fr
e
s
t
i
t
u
t
i
o
nf
o
re
l
a
s
t
i
cd
i
s
k
s
.Phys.Rev.E，
59:23612372，1
9
9
9
.
[
3
1
1 M.Y.LougeandM.E
.Adams. Anomalousb
e
h
a
v
i
o
ro
fnormalk
i
n
e
m
a
t
i
cr
e
s
t
i
t
u
t
i
o
ni
n
i
cp
l
a
t
e
.P
h
y
s
.R
e
v
.E，
65:021303-1t
h
eo
b
l
i
q
u
eimpactso
fahards
p
h
e
r
eonane
l
a
s
t
o
p
I鉛 t
021303
6
，2
0
0
2
.
.W.Kennedy
，A.Chatterjee，A.Ruina，andJ
.T.J
e
n
k
i
n
s
. Anomalous
[
3
2
]J
.C
a
l
s
a
m
i
g
l
i
a，S
.ofA
p
p
l
.Mech.，
66:146-152，1
9
9
7
.
f
r
i
c
t
i
o
n
a
lb
e
h
a
v
i
o
ri
nc
o
l
l
i
s
i
o
n
so
ft
h
i
nd
i
s
k
s
.J
.Hayakawa.Ar
問 n
a
l
o
u
sb
e
h
a
v
i
o
ro
ft
h
ecoe
伍c
i
e
瓜 o
fnormalr
e
s
t
i
t
u
t
i
o
n
[
3
3
] H.Kuni泊 kaandH
i
no
b
l
i
q
u
ei
m
p
a
c
t
.P
h
y
s
.R
e
v
.L
e
t
t
.，
93:154301-1-4，
2
0
0
4
.s
e
ea
l
s
oP
h
y
s
.Rev.Focus14，
1
4
(
2
0
0
4
)
;Nature，431(
2
0
0
4
)1
0
51
.
(
3
4
)L
i
n
l
e
yE
r
i
nH
a
l
l，
e
d
i
t
o
r
. TheLω s
'o
fMotion: AnAnthology0/CurrentThought(ContempomryD
i
s
c
o
u
r
s
e初 t
h
eF
i
e
l
d0
/Physics). RosenPub.Group，NewYork，2005.
s
a
t
o，
andR
.N
.Dave.Measurementso
fc
o
l
l
i
註
s
i
ぬ
on
凶 p
r
o
p
e
r
t
i
e
so
fs
p
h
e
r
e
s
[
3
5
JL
.Labous，A.D.豆o
u
s
i
n
gh
i
g
h
s
p
e
e
dv
i
d
e
oa
n
a
l
y
s
i
s
.P
h
y
s
.Reυ.E，
56:5717-5725，
1
9
9
7
.
[
3
6
j A.L
o
r
e
n
t
z，C
.Tu
o
z
z
o
l
o，andM.Y.L
o
u
g
e
. Measurementso
fimpactp
r
o
p
e
r
t
i
e
so
fs
m
a
l
l
n
e
a
r
l
ys
p
h
e
r
i
c
a
lp
a
r
t
i
c
l
e
s
.E
x
p
.Mech.，
37:292-298，1
9
9
7
.
方2A
ゥ
，
FO
巨視的物 f
本の非弾性衝突
[
3
7
]S
.F
.F
o
e
r
s
t
e
r，M.Y.Louge，H.Chang，andK.A
l
l
i
a
. Measurementso
fthec
o
l
l
i
s
i
o
nproph
y
s
.F
l
u
i
d
s，6
:
1
1
0
8
1
1
1
5，1
9
9
4
.
e
r
t
i
e
so
fs
m
a
l
ls
p
h
e
r
e
s
.P
[
3
8
] O.R.WaltonandR.L
.Braun. V
i
s
c
o
s
i
t
y
，granular-temerature，ands
t
r
e
s
sc
a
l
c
u
l
a
t
i
o
n
sおr
s
h
e
a
r
i
n
ga
s
s
e
m
b
l
i部 o
fi
n
e
l
a
s
t
i
c，
f
r
i
t
i
o
n
a
ld
i
s
k
s
.J
.R
h
e
o
l
.，
3
0
:
9
4
9
9
8
0，1
9
8
6
.
[
3
9
] O. R.Walton. Numericals
i
m
u
l
a
t
i
o
no
fi
n
e
l
a
s
t
i
c，f
r
i
c
t
i
o
n
a
lp
a
r
t
i
c
l
e
p
a
r
t
i
c
l
ei
n
t
e
r
a
c
t
i
o
n
s
.
I
n M. C
.1
ミ
oco，e
d
i
t
o
r，Pαr
t
i
c
u
l
a
t
e Two Ph
αs
e Flow
，pages 884-907，Boston， 1992.
But
t
e
r
w
o
r
t
h
-Heinemann.
[
4
0
] N.Maw，J
.R
.Barber，andJ
.N.F
a
w
c
e
t
t
. Theo
b
l
i
q
u
eimpacto
fe
l
a
s
t
i
cs
p
h
e
r
e
s
. Wear
，
3
8
:
1
0
1
1
1
4，1
9
7
6
.
.R
. Barber，andJ
.N
.F
a
w
c
e
t
t
. Ther
o
l
eo
fe
l
a
s
t
i
ct
a
n
g
e
n
t
i
a
lc
o
r
叩 l
i
a
n
c
ei
n
[
4
1
] N. Maw，J
o
b
l
i
q
u
ei
m
p
a
c
t
.J
.L
u
b
. Tech
，103:74-80，1981
.
[
4
2
] H.KuninakaandH.Hayakawa. Simulationf
o
rt
h
eo
b
l
i
q
u
eimpacto
fal
a
t
t
i
c
es
y
s
t
e
m
.J
.
P
h
y
s
.S
o
c
.J
p
n
.，7
2
:
1
6
5
5
1
6
6
3，2
0
0
3
.
.A.CundallandO.D.L
.S
t
r
a
c
k
.Ad
i
s
c
r
e
t
en
u
m
e
r
i
c
a
lmodelf
o
rg
r
a
n
u
l
a
ra
s
s
e
m
b
l
i
e
s
.
[
4
3
]P
G
e
o
t
e
c
h
n
i
q
u
e，2
9
:
4
7
-65，1
9
7
9
.
[
4
4
] 伯野元彦.破壊のシミュレーションー拡張調関要素法で破壊を追うヘ森北出版株式会社，東京，
1
9
9
7
.
9
9
8
.
[
4
5
] 粉体工学会福.粉体シミュレーシヨン入門.産業図書，東京， 1
[
4
6
] ロパート・アデア〈中村和幸訳) .ベースボールの物理学.紀伊国屋書庖，東京， 1
9
9
6
.
[
4
7
]1
.Newton.P
h
i
l
o
s
h
o
p
h
iαenαt
u
r
a
l
i
sP
r
i
n
c
i
p
i
αmαthem
αt
i
c
a
.W.DawasonandSons，
London，
1
9
6
2
.
9
6
0
.
[
4
8
] ガリレオ・ガリレイ著〈今野武雄・日田節次訳) .新科学対話.岩波書庖，東京， 1
0
0
6
.
[
4
9
] エルンスト・マッハ.マッハ力学史.筑摩書房，東京， 2
[
5
0
] 安田徳太郎訳・編.新訳ダンネマン大自然科学史 4.三省堂，東京， 1
9
7
8
.
a
r
c
I
. DePropo
付i
o
n
eMotus.Prague，1
6
3
9
.
[
5
1
] M.M
[
5
2
] 安田徳太郎訳・編.新訳ダンネマン大自然科学史 5
.三省堂，東京， 1
9
7
8
.
e
r
t
z
.Uberd
i
eberuhrungf
e
s
t
e
re
l
a
s
t
i
s
c
h
e泌 r
p
e
r
.J
.Re
i
n
eAnge切 . Math.，
9
2
:
1
5
6
1
7
1，
(
5
3
] H.H
1
8
8
2
.
t
t
ヴ
ウ
園仲寛人
[
5
4
] A.E
.H.L
o
v
e
.AT
r
e
a
t
i
s
eont
h
eMα
，t
hem
αt
i
c
a
lTheory0fEl
αs
t
i
c
i
t
y
. CambridgeUn
i
v
.
P
r
e
s
s，1
9
2
7
.
魁u
t
z
. Theoryo
fEl
αs
t
i
c
i
t
y(2ndEπg
l
i
s
he
d
.
)
. Pergamon，1
9
6
0
.
[
5
5
]L
.D.LandauandE
.M.L
i
a
sunderap
l
a
n
es
h
e
a
r
.P
h
y
s
.R
e
v
.
1
5
6
1 K.S
a
i
t
o
handH.Hayakawa. Rheologyo
fag
r
a
n
u
l
a
rg
E，75:021302-1-021302-11，
2
0
0
7
.
[
5
7
]h
t
t
p
:
/
/
n
e
w
s
.b
b
c
.
c
o
.uk/1
/hijmagazine/7022464.
s
t
m
.
F
.Hersen，
andL
.B
o
c
q
u
e
t
.S
e
c
r
e
t
so
fs
u
c
c
e
s
s
f
u
ls
t
o
n
←skipping.Natuκ427:29[
5
8
}C
.C
l
a
n
e
t，
29，
2
0
0
4
.
[
5
9
]S
.NagahiroandY.Hayakawa.T
h
e
o
r
e
t
i
c
a
landn
u
m
e
r
i
c
a
lapproacht
omagica
n
g
l
eo
fs
t
o
n
e
s
k
i
p
p
i
n
g
.P
h
y
s
.R
e
v
.Lε
t
t
，94:174501-1-174501-4，2005.
andK.五o∞ S
p
l
a
s
hproducedbyasmooths
p
h
e
r
eodc
i
r
c
u
l
a
r
[
6
0
J G.Kuwabara，H.Tanaka，
c
y
l
i
n
d
e
rs
t
r
i
k
i
n
gal
i
q
u
i
ds
u
r
f
a
c
e
.J
.P
h
y
s
.S
o
c
.J
p
n
.，
56:2733-2743，1
9
8
7
.
[
6
1
]永弘進一郎.物体と涜体表面の衝突過程の数値的および理論的研究. PhDt
h
e
s
i
s，東北大学，
2
0
0
4
.
[
6
2
]P
.A
u
s
s
i
l
l
o
u
sandD.Q
u
e
r
e
.L
i
q
u
i
dm
a
r
b
l
e
s
.Nα如何， 411:924-927，2001
.
Y.Yamazaki，
andK.Okumura. Bouncingg
e
lb
a
l
l
s
: impacto
fs
o
丘g
e
l
sonto
[
6
3
] Y.Tanaka，
u
r
o
p
h
y
s
.L
e
t
t
.，
6
3
(
1
)
:
1
4
6
1
5
2，
2
0
0
3
.
r
i
g
i
ds
u
r
f
a
c
e
.E
a
n
s
i
e
n
tc
r
a
t
o
rgrowthi
ng
r
a
n
u
l
a
r
[
6
4
]S
.Yamamoto，K.Wada
，N.Okabe，andT.M
a
t
s
u
I
. Tr
c
α
r
u
s，
t
a
r
g
e
t
s
: Ane
x
p
e
r
i
m
e
n
t
a
lstudyo
flowv
e
l
o
c
i
t
yimpactsi
n
t
og
l
a
s
ss
p
h
e
r
et
a
r
g
e
t
s
.I
1
8
3
:
2
15-2
2
4，
2
0
0
6
.
[
6
5
] K.Wada，
H.Senshu，
andT
.M
a
t
s
u
i
.Numericals
i
m
u
l
a
t
i
o
no
fimpactc
r
a
t
e
r
i
n
gong
r
a
n
u
l
a
r
m
a
t
e
r
i
a
Ll
ωrus，1
8
0
:
5
2
8
5
4
5，
2
0
0
6
.
ら d
y
n
a
.
j
r
i
s
o
.
lc
o
.
j
p
/.
[
6
6
]h
t
t
p
:
/
/
I
[
6
7
] G.R.JohnsonandR.A.S
t
r
y
k
.Conversiono
f3dd
i
s
t
o
r
t
e
de
l
e
m
e
n
t
si
n
t
om錨 h
l
e
s
sp
a
r
t
i
c
l
e
s
n
t
.よ ImpactE
n
g
.，
2
8
:
9
4
7
9
6
6，
2
0
0
3
.
duringdynamicd
e
f
o
r
m
a
t
i
o
n
.I
[
6
8
] A.No
巾 n
，G 司法， B
.Bacon，J
.Gerth，andP
.Sweeney.Animationo
ff
r
a
c
t
u
r
ebyp
h
y
s
i
c
a
l
i
s
u
a
lComputer
，
7
:
2
1
0
2
1
9，
1991
.
m
o
d
e
l
i
n
g
.V
[
6
9
] :塚越誠一.粒子法〈計算力学レクチャーシリーズ 5) .丸善掠式会社，東京， 2
0
0
5
. 日本計算
工学会編.
- 718-
巨視的物体の非弾性笹突
[
7
0
] Y.Hayakawa.Numericalstudyo
fo
s
c
i
l
l
a
t
o
r
yc
r
a
c
kpropagationthroughatw
争 d
imensional
l
.P
h
y
s
.R
e
v
.E，
49:R1804-R1807，1
9
9
4
.
c
r
y
s
t
a
[
7
1
] Y.Hayakawa.Patterns
e
l
e
c
t
i
o
no
fm凶 i
c
r
a
c
kp
r
o
p
a
g
a
t
i
o
ni
nquenchedc
r
y
s
t
a
l
s
.P
h
y
s
.R
e
v
.
E，50:R1748-R1751，
1
9
9
4
.
山 t
i
o
n
a
lS
tαt
i
s
t
i
c
a
lMechα
n
i
c
s
.E
l
s
e
v
i
e
rS
c
i
e
n
c
e，
Amsterdam，1991
.
[
7
2
] W.G.H
o
o
v
e
r
. Compl
[
7
3
] H.KuninakaandH.Hayakawa. Theimpacto
ftw
か d
i
m
e
n
s
i
o
n
a
le
l
a
s
t
i
cd
i
sk
.J
. ofPhys.
S
o
c
.J
p
n
.，7
0:2220-2221，
2001
.
[
7
4
] 豆 HayakawaandH.Kuni凶 k
a
.S
i
m
u
l
a
t
i
o
nandt
h
e
o
r
yo
ftheimpacto
ftw
か d
imensional
e
l
a
s
t
i
cd
i
s
k
s
. Chem.Eng.S
c
i
.，5
7
:
2
3
9
2
5
2，
2
0
0
2
.
[
7
5
司
] H.五
1
(
u
n
i
d
出i
s
北k
.M
a
s
t
e
r
'
st
h
e
s
i
s，KyotoU
n
i
v
e
r
s
i
t
y
，2000.
n
i
n
a
k
a
. Theor
εt
i
c
a
lα
ηdNume
π
c
αlS
t
u
d
i
e
so
fl
n
e
l
αs
t
i
clmp
αc
t
so
fE
l
a
s
t
i
cMate1
'
Ia
l
s
.
[
7
6
] 日.Ku
PhDt
h
e
s
i
s，KyotoU
n
i
v
e
r
s
i
t
y
，2004.
沿 k
a
. Coe
伍c
i
e
n
to
fr
e
s
t
i
t
u
t
i
o
no
fe
l
a
s
t
i
cd
i
s
k
s
.I
nY.Kishino，
[
7
7
]H
.HayakawaandH.Kunir
e
d
i
t
o
r，Thep
r
o
c
e
e
d
i
n
g
so
fpowd
εr
sandg
r
a
i
n
s2001，
paεes5
6
1
5
6
4
.Rotterdam:召a
l
主ema
，
2
0
0
1
.
[
7
8
] H.HayakawaandH.Kuninaka. S
i
m
u
l
a
t
i
o
no
ft
h
eimpacto
ftw
o
-d
i
m
e
n
s
i
o
n
a
le
l
a
s
t
i
cd
i
s
k
s
.
l
nTheproceedingsoft
h
en
i
n
t
hn
i
s
s
h
i
ne
n
g
i
n
e
e
r
初 9p
a
r
t
i
c
l
et
e
c
h
n
o
l
o
g
yi
n
如、n
a
t
i
o
n
a
lsym-
，pages82-95，
2
0
0
1
.
plsium
か d
i
m
e
n
s
i
o
n
a
le
l
a
s
t
i
cd
i
s
k
s
.
[
7
9
] H.KuninakaandH.Hayakawa.S
i
m
u
l
a
t
i
o
no
ft
h
eimpacto
ftw
u
i，Y.Sugiyama，M.S
c
h
r
e
c
k
e
n
b
e
r
g，andD.E
.Wolf
，e
d
i
t
o
r
s，Thep
r
o
c
e
e
d
i
n
g
sof
I
n担. Fuk
T
r
i
α
f
.
f
icαndGmnlarFlow2001，pages467-472，Berlin，2003.Springer.
f
8
0
} ;杉原厚吉.データ構造とアルゴリズム共立出版，東京， 2
0
0
1
.
ぉ t
i
cm
a
t
e
r
i
a
l
s
. Phase
[
8
1
] H. HayakawaandH.Kuninaka. Theoryo
fi
n
e
l
a
s
t
i
cimpacto
fe
l
Tr
αn
s
.，
77:889-909，
2
0
0
4
.
，
[
8
2
] H.五uninakaandH.Hayakawa. Contactandq
l
腿 トs
t
a
t
i
cimpacto
fhamiltons
y
s
t
e
m
.J
.
P
h
y
s
.S
o
c
.J
p
n
.，
75:074001-1-074001-5，
2
0
0
6
.
[
8
3
] K.S
u
g
i
h
a
r
a
. Dα
t
αS
t
r
u
c
t
u
r
eandA
l
g
o
r
i
t
h
m
s
(
i
nJ
a
p
a
n
e
s
e
)
.K
y
o
r
i
t
s
u，
Japan，
2001
.
[
8
4
] G.Sundar
泣 a
j
a
n・ Theenergyabsorbedd
u
r
i
n
gt
h
eo
b
l
i
弓u
eimpacto
fahardb
a
l
la
g
a
i
n
s
t
d
u
c
t
i
l
et
a
r
g
e
tm
a
t
e
r
i
a
l
s
.l
n
t
.J
.lmp
αc
tEng
ηg，
9
:
3
4
3
3
5
8，1
9
9
0
.
円
i
QU
盟仲寛人
[
8.
5
]L
.A.G
a
l
i
n
. Contactproblemsi
nt
h
et
h
e
o
r
y01e
l
αs
t
i
c
i
t
y
. Gostekhizdat Moskow，
1
9
5
3
.
ぅ
[
8
6
] A.Awasthi，S
.C
.Hendy
，P
.Z
o
o
t
j
e
n
s，S
.A.Brown，andF
.N
a
t
a
l
i
. 話o
l
e
c
u
l
a
rdynamics
s
i
m
u
l
a
t
i
o
n
so
fr
e
f
t
e
c
t
i
o
nandadhesionb
e
h
a
v
i
o
ri
nl
e
n
n
a
r
d
j
o
n
e
sc
l
u
s
t
e
rd
e
p
o
s
i
t
i
o
n
.P
h
y
s
.
R
e
v
.B，
76:115437-1-115437-11，
2
0
0
7
.
.EvansandD.J
.S
e
a
r
l
e
s
.The註u
c
t
u
a
t
i
o
ntheorem.Adv.P
h
y
s
.，
51
:1
5
.29-1585，2
0
0
2
.
[
8
7
] D.J
担当日.KuninakaandH.Hayakawa.arXiv:0707.0533.
.V.Woodcoc
k
.I
s
o
t
h
e
r
m
a
lm
o
l
e
c
u
l
a
rdynamicsc
a
l
c
u
l
a
t
i
o
n
sf
o
rl
i
q
u
i
ds
a
l
t
s
. Chem.P
h
y
s
.
[
8
9
]L
L
εt
t
.，10:257-261，
1
9
71
.
[
9
0
]S
.N
o
s
e
. Constanttemperaturem
o
l
e
c
u
l
a
rdynamicsmethods. P
r
o
g
. Th
εo
r
二 P
h
y
s
.Supp
，
.
l
(
1
0
3
)
:
1
4
6，1991
.
[
9
1
] H.HayakawaandH.Kuninaka.a
r
X
i
v
:
0
8
0
5
.
0
1
7
6
i
伊
9
耳
2
]N
持.V
. Br
討i
l
l
i
弘
a
剖
問n
t
p
a
r
t
i
c
l
e
s叩
wi
t
泊ha
低
吋
dh街
e
矧
釘
Si
∞
on
.P
h
y
s
.R
e
v
.E，
7
6
:
0
5
1
3
0
2
1-051302-12，2
0
0
7
.
凶 a
l
l，andA.D.R
o
b
e
r
t
s
.S
u
r
f
a
c
eenergyandt
h
ec
o
n
t
a
c
tofe
l
a
s
t
i
c
[
9
3
] K.L
.Johnson，K.Ke
，S
e
r
.A，
324:301-313，1971
.
s
o
l
i
d
s
.P
r
o
c
.R.S
o
c
.London
[
9
4
]J
.N.イスラエルアチヴィザ分子関力と表面力.輯倉書届，東京， 1
9
9
6
.
5
] H.T
a
s
a
k
i
.F
l
u
c
t
u
a
t
i
o
nr
e
l
a
t
i
o
n，j
a
r
z
y
n
s
k
i
t
y
p
ee
q
u
a
l
i
t
y
，andthe関 condlawf
o
rg
e
n
e
r
a
l
[
9.
，number
c
o
l
l
i
s
i
o
np
r
o
b
l
e
m
s
.I
nRIMSw
o
r
k
s
h
o
p
: Micro-macroduαl
i
t
yi
nquαntumc
a
l
c
u
l
u
s
0
0
6
.RIMS，KyotoU
n
i
v
.
1507，pages102-117，Kyoto，2
[
9
6
]h
t
t
p
:
jj
w
i
r
e
d
v
i
s
i
o
nぬ jarchivesj200306j2003060602.htm.
l
[
9
7
]h
t
t
p
:
/
/
w
e
b
u
s
e
r
s
.
n
p
l
.
山
n4
門
i
n
v