サブミクロン粒子衝突の分子動力学シミュレーション
田中秀和（北海道大学低温科学研究所）
１. 背景と研究目的
原始惑星系円盤内でのダスト合体成長過程においてダストの衝突速度は時速百 km 以上に達するが， このような高速衝突においてダストが破壊せず合体成長を続けることができるのかということが問題になっている． この問題は微惑星形成過程を考える上で重要である． 和田らによるダスト衝突数値計算によって， 氷ダスト は時速 200km 以下の衝突速度であれば合体成長することが示された． しかしながら， 彼らの数値計算 において仮定されている， ダストを構成するサブミクロン微粒子間の JKR 理論をもとにした相互作用モデルがどの程度正しいのかについては， まだ不定性が残されている．
本研究の目的は， サブミクロン微粒子１つを多数の「分子」で構成しその微粒子間の相互作 用を分子動力学計算行って調べることにより， 従来の相互作用モデル（JKR 理論）の妥当性を明らかにすることである． ナノサイズを越えた粒子の相互作用を調べる分子動力学計算は理学や工学の諸分野でもほとんど行われておらず， 本研究は未開拓分野への挑戦的な課題といえる． 今回は第一段階として微粒子の正面衝突の際の相互作用を調べた．

２. 計算方法
多数の「分子」で構成された微粒子の衝突を分子動力学計算を行って調べる． 計算から得られる衝突の際の微粒子の加速度を解析することで微粒子間の相互作用を明らかにすることができる． 分子間の相互作用としては簡単のためレナード・ジョーンズ型相互作用を仮定した． １つの微粒子を構成する分子数はパラメータとして数万から 1.5 千万までのいくつかの場合の計算を行っている． 以下では主に微粒子構成分子数が約 300 万の場合の計算結果を紹介する． 衝
突速度はもう１つのパラメータであり， 様々な衝突速度に対し計算を行った． 各分子の運動の数値時間積分には分子動力学計算でよく用いられる leap-frog 法を用いた． 計算において分子温度は操作せず全エネルギーが保存させた． 数値積分における時間ステップはレナード・ジョーンズの単位 (LJ 単位、以下参照) で 0.01 とした．
微粒子の分子構造は面心立方格子結晶とした．（レナード・ジョーンズ分子系においてはこの結晶が低温低圧において安定であることが知られている．） 結晶のある半径の球領域にある分子達から球粒子をつくった． さらに安定な初期状態をつくるため運動エネルギーを減衰させる分子動力学計算を行って緩和させた． このように用意した２つ微粒子を様々な初期速度で正面衝突させた． JKR 理論において接触した微粒子間の相互作用は表面張力とヤング率で表される．
レナード・ジョーンズ分子系において表面張力は LJ 単位で 3.0，ヤング率は接触面にもよるがおおよそ 50 程度である．
レナード・ジョーンズ分子は相互作用ポテンシャルの深さや相互作用距離がパラメータとなっており， これらを変化させることで様々な分子の相互作用を近似的に表すことができる． 特にアルゴン原子の相互作用は精度よく表すことができ, アルゴンの様々な熱力学量をレナード・ジョーンズ分子系で再現できる． アルゴン原子の場合のレナード・ジョーンズ単位の数値を列挙すると次のようになる: 単位長：0.34nm (相互作用距離), 単位質量：アルゴン原子質量, 単位エネル
ギー〜0.01eV, 単位速度：160m/s, 単位時間：2.2ps． 本研究の結果はレナード・ジョーンズ単位で表されているが， これらの数値を用いて具体的な値を得ることができる．

３. 結果
約 300 万分子で構成された微粒子の衝突結果を紹介する． この微粒子の半径 R は LJ 単位で 88 であり、アルゴンの系では 30nm である． 図 1 に， この微粒子２つを速度 0.05 と 0.06 で衝突させた場合の 2 粒子の重心間距離の時間進化を示した． 重心距離が半径 R の 2 倍程度となると 2 球は接触して互いに減速し始める． ある程度食い込んだところで停止しその後跳ね返り始める． 速度 0.05 の場合は最終的に接触が切れず付着し，0.06 の場合は接触が切れた．
このような 2 粒子の相対位置と速度の変化 (加速度) から 2 粒子間に働く相互作用を相対位置の関数として求めることができる． 図 2 にその結果を示した． 相互作用の力には 2 粒子が食い込むときと跳ね返るときとで多少のずれが見られる． このヒステリシスは 2 粒子が接触している際の表面分子再配置による緩和過程により生じている． 従来の理論モデルではこのようなヒステリシスは説明できていない． 速度 0.05 と 0.06 の間で相互作用の大きな違いはみられない．
図 2 では数値計算から得られた相互作用と従来の相互作用モデルとの比較も行った．比較する相互作用モデルとしては，巨視的理論である JKR 理論に分子相互作用距離というミクロの効果を入れて拡張した Maugis (1992) のモデルを採用した．この相互作用理論モデル（灰色の線）と数値計算は， 食い込む際にはよく合っている． 重心間距離が 2R より大きいところではずれているが，これはまだ接触していないためである．跳ね返りの際は，ヒステリシスのため理論モデルとはずれている． しかし， この時も理論モデルの線を水平方向にずらすことにより数値計算結果に合わせることができる．
図３では， 様々な衝突速度に対する計算結果における粒子の重心運動エネルギーの変化の様子を横軸重心距離で表した． この図より， 食い込み期と跳ね返り期の間でのヒステリシスによる重心運動エネルギーの散逸の様子をみることができる．（散逸したエネルギーは主に分子熱振動となっている．）
衝突速度が 0.06 以下では，その数値計算も同じ曲線上を通る． これは低速度では相互作用や重心運動エネルギー散逸の様子が速度に依らないことを示している． 一方， 衝突速度が0.06 を越えると，衝突速度とともにエネルギー散逸が大きくなる． このエネルギー散逸の増大のため， 衝突速度 0.15 以上の計算では 2 粒子は跳ね返らず付着する結果となっている．
計算でみられたヒステリシスによるエネルギー散逸は JKR理論ではおこらない． その代わり 2 球の接触する際とその接触が切れる際に， 重心運動エネルギーの一部が粒子表面の弾性波エネルギーに変換され， 実効的なエネルギー散逸がおこるとされている． 本研究における低速度の場合のエネル
ギー散逸量は JKR 理論での値の約 1.5 倍程度になっている．
構成粒子数を変えることで粒子半径を変えて同様な計算を行ったが， エネルギー散逸量の JKR 理論からのずれはどれも同じ程度であり， 粒子半径依存性はみられなかった （図４）．

図１： ２粒子の重心間距離の時間進化．

４. まとめと今後の課題
以上の結果をまとめると次のようになる．粒子正面衝突の数値計算から得られた相互作用は従来の JKR 理論（を拡張したモデル）とおよそ一致することが示された．跳ね返りの際のエネルギー散逸については低速度衝突では JKR 理論と同程度であるが，高速衝突の場合エネルギー散逸量は非常に大きくなる．またエネルギー散逸機構は従来の理論モデルと定性的に異なっており，さらに詳細な検討が必要であろう．
正面衝突の他に、2 粒子をすべらせた場合や転がした場合の数値計算も容易に行うことができ， 今後計算する予定である． これらの計算から， すべり摩擦やころがり摩擦を得ることができる． このような結果をもとにして従来の理論モデルにかわる新しい粒子相互作用モデルを構築することがこの研究の最終目標であろう．

図３： ２粒子の運動エネルギーの変化の様子．

図２： ２粒子相互作用の計算値と理論モデルの比較．
灰色の実線は Maugis(1992)による理論線．点線
はそれを横にずらしたもの.

図４： 運動エネルギー散逸量の粒子サイズ依存性