数理解析研究所講究録 1413 巻 2005 年 150-155

150

2 次元弾性体モデルによる接触及び低速衝突シミュレーション

京都大学理学研究科化学専攻

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

國仲

京都大学理学研究科物理学・宇宙物理学専攻
Departmcnt of Chemistry, Department of Physics,
Kyoto University

1

Introduction

粉体等の多数の粒子からなる系を理論的または数値的に扱うときに重要なのは、対象とする系の動的な振る舞いの素過程てある、粒子間の接触や衝突を正確に模倣することである。本講演ては巨視的な弾性体をばねと質点から構成される 2 次元の格子モデルて表現し、粒子間の接触や低速衝突のシミュレーションの結果と理論から予測される結果とを比較する。

2モアル
著者の用いたモデ$\mathrm{K}\mathrm{s}1\mathrm{h}$1099 個の質点を円内にランダムに配置して作ったランダム格子モデルてある (図 1)。 質点間はデローニー三角分割のアルゴリズムを用いて近接の質点同士を接続してぃる。各質点間の相互作用は非線形ばねによる復元力に加えてばね両端の相対速度に比例する摩擦カが働くとして, 以下の式て表す。
$f(\mathrm{x}_{1j}.)=-k_{a}\mathrm{r}_{j}.-k_{b}\mathrm{x}_{ij}^{\mathrm{a}}+\eta_{ds}|.(\mathrm{v}:-\mathrm{v}_{j})$.

(1)

ここて右辺第二項まては非線形ばね相互作用を表しており,はばね両端の粒子 $ij$ 間の自然長がらの伸ひてある。
はばね定数てあり, その値は $k_{a}=1.0\mathrm{x}mc^{2}/R_{\text{ }}k_{b}=k_{a}\mathrm{x}10^{-3}/R^{2}$ という値を用いている。 ここて $m$ は質点の質量、 は円盤の半径、 は 1 次元音速てある。右辺第$\mathrm{x}_{ij}$$k_{a\text{、}}k_{b}$、$R$$c$

三項は質点間の摩擦力を表したものて、両端の粒子の相対速度に粘性係数する。粘性係数3. は$\eta_{d*s}$$\eta_{d\dot{\iota}s}=1.0\mathrm{x}10^{2}m/(R/c)$$\eta dis$を乗じた形て導入という値を用いる。

接触のシミュレーション
このようにして導入した弾性円盤を壁に接触させ、圧縮力とそれにょって生じる歪みの関係を求め、 Hcrtz の接触理論 [1] と比較する。一般に弾性球の接触において圧縮力 $P$ にょって生じる歪み$\xi$

は、 $P^{2/3}$ に比例する形て表され、その係数は圧縮する物質のヤング率やポアソン比にょって表される [1]。 2 次元の場合は以下の関係式て表すことがてきる [2]。
$\xi\simeq b\frac{P}{\pi E}\{ln(\frac{4\pi ER}{bP})-1-\nu\}$

ここて、 $E$ はヤング率、$\nu$,$b=1-\nu^{2}$.

(2)はポアソン比てある。モデルのボアソン比は同様のアルゴリズムて生成した長方形を牽引した後横方向の伸ひに対する縦方向の縮みを直接測定することて計算てき、 このモデノレの場合、$\nu\simeq 7.5\mathrm{x}10^{-2}$てあることがわがってぃる [3]。従って$b$の値は$b\simeq 1$となる。

151

図 1: ポテンシャル壁と接触する円盤モデル

円盤の接触シミュレーションは次のように行なう。 円盤を構成する各質点 の運動方程式は次のように書くことが出来る。
$i$.

$m \frac{d^{2}\mathrm{r}_{i}}{dt^{2}}=\sum_{j=1}^{N}.\cdot\{-k_{a}\mathrm{x}_{*j}.-k_{b}\mathrm{x}_{1j}^{3}.-\eta(\mathrm{v}:-\mathrm{v}_{j})\}+aV_{0}\exp(-ay)\hat{\mathrm{y}}-\frac{\mathrm{P}}{N}\hat{\mathrm{y}}$ここては質点 の位置、 は時間、 . は質点 の速度、$i$$\mathrm{r}_{\dot{\iota}}$

$i$

$t$

$\mathrm{v}_{1}$$N_{i}$(3)は質点 に接続されているばねの本数$i$てある。 また指数関数の項は円盤を接触させる壁を表現した項てあり、 $a=300/R_{\backslash }V_{0}=amc^{2}/R$という値を用いている。 は垂直方向の単位ベクトル、 $N$ は質点数 $N=1099$ てある。
$\hat{\mathrm{y}}$は圧縮力て、式 (3) を数値的に解いて全質点の位置を時間発展させると、徐々に円盤全体の振動が緩和していく。振動が緩和した時の重心と円盤下端の距離を $R’$ とし、円盤の歪みを $\xi\equiv|R-R’|$$\mathrm{P}$て定義する。 シミュレーションては圧縮力 $P$ を 4 $0\cross 10^{-3}mc^{2}/R$ から 10 $\mathrm{x}10^{-2}mc^{2}/R$ の範囲て変化させ、圧縮力 $P$ と歪み
図2の$\mathrm{x}$$\xi$の関係を求める。

印はシミュレーションから求めた圧縮力 $P$ と歪み$\xi$の関係てある。 ここて横軸は \pi RE、縦軸は $R$ てスケールし、無次元化した。
この結果と 2 次元の Hcrtz の接触理論との比較を試みる。 (2) 式の をフイツテイングパラメータとし、 $b=1.55$ という値てシミュレーション結果をフイットした結果が図 2 の実線てある。 この$b$図からモデルの接触シミュレーションの結果は、 2 次元 Hcrtz 理論から予想される関係と、ほぼ調和的であることがわかる。

152$\mathrm{P}/(\pi \mathrm{R}\mathrm{E})$

図 2: 圧縮力と歪みの関係。

4

低速衝突のシミュレーション

次に、円盤をポテンシャル壁に正面衝突させ、衝突速度とはねかえり係数の関係を調べる。一般に巨視的な物体同士が低速て正面衝突したときの衝突速度とはねがえり係数の関係は Hertz の接触理論を基にした準静理論 [4] で説明することがてきる。

準静理論ではますはじめに衝突時間中に物体に生じる歪み の時間発展方程式を以下の様に表す。
$\xi$$M \frac{d^{2}\xi}{dt^{2}}=-\frac{\pi E_{*}\xi}{\ln(4R/\xi)}-\frac{\pi\tau_{0}E_{*}}{\ln(4R/\xi)}\frac{d\xi}{dt}$.

(4)は有効ヤング率てあり、 $E_{*}=E/1-\nu^{2}$ て表される。 また、 は散逸の時間スヶ–ルである。 右辺第一項は 2 次元の Hertz の接触理論から求められる弾性$j$] てある [2]。 右辺第二項はここで、$E_{*}$$\tau_{0}$物質の粘性に起因する散逸力てあり、歪み速度に比例した形て導入する [4]。 この方程式を初期条( $\tau=\pi(R/c)\sqrt{\ln(4c/v)}.\cdot$ 接触時間 [5]) を求めることにょっ\mbox{\boldmath $\xi$}.(0)=v. て解き、件$\xi(0)=0_{\text{、}}$$\dot{\xi}(\tau)$て. はねかえり係数$e\equiv-\dot{\xi}(\tau)/\dot{\xi}(0)$を計算すると、1(5)$e\simeq d_{1}\sqrt{1-\frac{d_{2}}{\sqrt{\ln(4c/v)}}}$となる。 ここで $d_{1}=1_{\text{、}}d_{2}=2\pi\tau_{0}E_{*}/\rho Rc$ てある。

シミュレーションは、次に示す各質点の運動方程式を数値的に解くことになる。
$m \frac{d^{2}\mathrm{r}_{\dot{l}}}{dt^{2}}=\sum_{j=1}^{N_{l}}\mathrm{t}-k_{a}\mathrm{x}_{\dot{\iota}j}-k_{b}\mathrm{x}_{1j}^{\theta}.-\eta a_{1s}.(\mathrm{v}:-\mathrm{v}_{j})\}+aV_{0}\exp(-ay)\hat{\mathrm{y}}$ここで初期条件は各質点の初速度を 30$\mathrm{x}10^{-8}R/c$

速度に対するはねかえり係数を計測した。

.

(6)

がら 50 $\mathrm{x}10^{-2}R/c$ の範囲て変化させて、各初153$e$$\mathrm{v}/\mathrm{c}$図 3: 衝突速度とはねかえり係数
図 3 は衝突速度とはねかえり係数の関係をプロットしたものてある。 ここて、 $\eta^{*}=\eta d:s\mathrm{x}R/mc$と定義し、 $\eta^{*}=0.0$ から $\eta^{*}=1.0\cross 10^{2}$ まて変化させて各 \eta *1 こついてプロットした。 $\eta^{*}=0.0$ の時は上に凸の曲線を描くため、下に凸の曲線になる式 (5) でフイットさせることはてきないが、の値が大きくなるにつれ下に凸の曲線に移行し、 の値もほぼ 1 に近い値てフイットできることが$\ovalbox{\tt\small REJECT}$$d_{1}$

わかった。

5

エネルギー吸収壁の導入

ばね相互作用に摩擦力を追加すれば我々のモデル計算の結果は、既存の接触理論や衝突理論とほぼ一致させることがてきることはわかった。最後に質点の相対速度に比例する摩擦力のようなローカルな散逸メカニズムを導入しなくても衝突理論を回復するような ミルトンモデルの例を示す。
円盤を構成する質点間相互作用は式 (1) において $\eta_{dis}=0$ とし、完全な保存系とする。壁に円盤$J\backslash$の衝突エネルギーが散逸する状況を実現するために、これまての指数関数型の壁ポテンシャ の代わりに弾性壁を用意する。弾性壁は 1269 個の質点を長方形状にランダムに配置させ、質点間相互$J\mathrm{s}$作用は円盤と同様とする (図 4)。ばね定数は円盤、壁共に $k_{a}=1.0\mathrm{x}$とし、$\eta d\dot{\mathrm{a}}s$$m\text{♂}/R_{\text{、}}$$k_{b}=k_{a}\cross 10^{-3}/R^{2}$

は共に 0 とした。壁の境界条件は、左右両端と底面の粒子のエネルギーを時間発展させるたひに 0 にリセットし、壁にエネルギーが吸収される状態を実現した。

また、 円盤と壁の相互作用は次のように定義する。 図 5 は円盤と壁の相互作用の模式図てある。
円盤の下端の粒子の真下にあるばねを常にサーチしておき, 壁との距離 がしきい値よりも小$l_{s}$

さくなると、 円盤表面の粒子は真下のばねから距離$\mathrm{F}(l_{8})=aV_{0}\exp(-al_{s})\mathrm{n}_{s}$$\mathrm{n}_{\theta}$$l_{\epsilon}$を引数とする指数関数型のポテンシャルカを受ける。 ここて、 $a=300/R,V_{0}=amc^{2}R/2$ とう値を用 ‘た。 また$\mathfrak{y}$$\mathrm{A}\backslash$

はばねの法線ベクトルてある。ばね両端の粒子は、 トルクがつりあうように反作用を受けるとする。 このようなルールて円盤と壁の相互作用をモデル化した。

154

図 4: 円盤とエネルギー吸収壁。

図 5: 円盤と壁の相互作用
図 6 に衝突速度とはねかえり係数の関係を示す。 (5) 式の$d_{1},d_{2}$をフィッティングパラメータとし、 $d_{1}=0.973,d_{2}=0.119$ という値てシミュレーション結果をフィットした結果が図 6 の実線てある。 フィッティングパラメータ$d_{1}$の値も 1 に近く、シミュレーションの結果は、準静理論から予想される関係とほぼ調和的てあることがわかる。

6

まとめ
本講演ては以下の項目を報告した。

1. 散逸メカニズムを導入したモデルて、接触及ひ低速正面衝突のシミュレーションを行った。
2. 接触に関して、圧縮力と歪みの関係のモデル計算の結果は 2 次元の Hcrtz の接触理論と良好な一致を示した。

3. 低速正面衝突に関して、衝突速度とはねかえり係数の関係は衝突の準静理論と良好な一致を示した。

4. 散逸メカニズムを導入しない、ハミルトン系のモデルてもエネルギーを吸収する壁を導入することて、衝突の準静理論とよい一致を示すことがわかった。

155$\mathrm{e}$$\mathrm{v}/\mathrm{c}$図 6: エネルギー吸収壁との正面衝突時の衝突速度とはねかえり係数。

参考文献[1] H. Hcrtz: J. Reine Angew. Math. 92, 156 (1882); L. D. Landau and E. M. Lifshitz: Theoryof Elasticity(2nd English ). Pergamon, 1960.
$\mathrm{c}\mathrm{d}.$[2] F. Gcrl and A. Zippclius: Phys. Rcv.

$\mathrm{E}59$,2361 (1999).

[3] H. Kuninaka and H. Hayakawa: J. Phys. Soc. Jpn. 72, 1655 (2003).
[4] N. V. Brilliantov, F. Spahn, J.-M. Hertzsch, and T. P\"oshcl: Phys. Rev.

Kuwabara and K. Kono: Jpn. J. Appl. Phys. 26, 1230 (1987).
[5] H. Hayakawa and H. Kuninaka: Chem. Eng. Sci. 57, 239 (2002).

$\mathrm{E}53$,5382 (1996);G.