MD 計算による微粒子衝突時のエネルギー散逸過程の解明
田中秀和（北海道大学低温科学研究所）
１. 背景と研究目的

惑星形成の第一段階である固体微粒子合体成長過程において微粒子（ダスト）の衝突速度は最大時速百 km以上に 達するが， この高速衝突においてダストが破壊せず微惑星へと成長を続けることができるのかということが問題になっている． 我々のグループによるダスト衝突数値計算によって， 氷ダスト は衝突速度が時速200km 以下であれば合体成長可能であることが示された． しかしながら， この数値計算 において仮定されている， JKR 理論をもとにした相互作用モデルが、ダストを構成するサブミクロン微粒子に対してどの程度有効なのかについては， 未だ不定性が残されている．
本研究の目的は， サブミクロン微粒子１つを多数の「分子」で構成しその微粒子間の相互作 用を分子動力学計算行って調べることにより， 従来の相互作用モデル（ JKR 理論）の妥当性を明らか にすることである．
ミクロンサイズ以下の微粒子相互作用は，上記の天文学的目的のみならず理工学の多くの分野に関連しており重要であるが，室内実験や大規模分子動力学計算を用いた研究はあまり行われておらず， 本研究は未開拓分野への挑戦的な課題といえる． 今回は第一段階として微粒子の正面衝突の際の相互作用について報告する．

２. 計算方法

多数の「分子」で構成された微粒子の衝突を分子動力学計算を行って調べる． 計算から得られる衝突の際の微粒子の加速度を解析することで微粒子間の相互作用を明らかにすることができる． 分子間の相互作用としては簡単のためレナード・ジョーンズ型相互作用を仮定した． １つの微粒子を構成する分子数はパラメータであり微粒子半径を決める．本研究では分子数が 5 万から 1 億という 2 千倍の範囲、半径では 13 倍程度の範囲に対して計算を行った．衝突する２つの微粒子で 2 億分子という大規模な分子動力学計算は，低温研の並列計算機環(SGI UV100)を用いることにより可能とした．衝突速度はもう１つのパラメータであり， 様々な衝突速度に対し計算を行った． 各分子の運動の数値時間積分には分子動力学計算でよく用いられる leap-frog 法を用いた． 計算において分子温度は操作せず全エネルギーを保存させた． 数値積分における時間ステップはレナー
ド・ジョーンズの単位 (LJ 単位、以下参照) で 0.01 とした．
微粒子の分子構造は面心立方格子結晶とした．レナード・ジョーンズ分子系においてはこの結晶が低温低圧において安定であることが知られている．数値計算における初期微粒子は結晶から与えた半径の球領域を切り取ることにより作成した． さらに微粒子を安定な初期状態とするために運動エネルギーを減衰させる分子動力学計算を行い緩和させた． このように用意した２つ微粒子を様々な初期速度で正面衝突させて，接触時の微粒子相互作用の詳細を調べた． JKR 理論においては微粒子間相互作用は表面張力とヤング率により表されている ．
レナード・ジョーンズ分子系において表面張力は LJ 単位で 3.0，ヤング率は接触面にもよるがおおよそ 60 程度である．
レナード・ジョーンズ分子は相互作用ポテンシャルの深さや相互作用距離がパラメータとなっており， これらを変化させることで様々な分子の相互作用を近似的に表すことができる． 特にアルゴン原子の相互作用は精度よく表すことができ, アルゴンの様々な熱力学量をレナード・ジョーンズ分子系でよく再現できる． アルゴン原子の場合のレナード・ジョーンズ単位の数値を列挙すると次のようになる．単位長： 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) のモデルを採用した．（但し，粒子半径が 88 と大きい場合には JKR 理論との差は小さい．）この相互作用理論モデル（灰色の線）と数値計算は， 食い込む際にはよく合っている． 重心間距離が 2R より大きいところではずれているが，これはまだ接触していないためである．跳ね返りの際は，ヒステリシスのため理論モデルとはずれている．相互作用のヒステリシスはエネルギー散逸を生じさせる．散逸したエネルギーは主に分子熱振動となっている．数値計算にみられるヒステリシスによるエネルギー散逸は JKR 理論ではおこらない． その代わり 2 球の接触する際とその接触が切れる際に， 重心運動エネルギーの一部が粒子表面の弾性波エネルギーに変換され， 実効的なエネルギー散逸がおこるとされている．

このエネルギー散逸の量は衝突速度の増加とともに増大する．また，エネルギー散逸の結果，跳ね返り定数は１より減少する．図３では， 衝突速度の関数として得られた跳ね返り定数を図示した．数値計算結果によるとエネルギー散逸及び跳ね返り定数は衝突する面に依存する．図 3 では、それぞれの衝突速度に対して，球の向きをランダムに変えることにより面を変え，１０回の衝突計算を行うことで，跳ね返り定数の平均値を求めた．比較のため JKR 理論で与えられる跳ね返り定数の値も図示した．これより MD 計算によって得られたエネルギー散逸量は（速度 0.05 以上では）JKR 理論からの値よりも大きいことが分かる．特に衝突速度が 0.1 以上では跳ね返り定数は理論値に比べて大きく減少する．これは高速度衝突により球が塑性変形されたためである．
塑性変形が重要でない衝突速度 0.1 以下の場合のエネル
ギー散逸は、2 球が接触している際に一定の抵抗力が働いていると考えることで説明できる．数値計算でのエネルギー散逸を説明するこの抵抗力は衝突速度に比例している．
同様な衝突計算を構成粒子数を変えることで粒子半径を変えて行った．速度一定で球半径を変えるとエネルギー散逸を生む抵抗力は球半径にほぼ比例することが明らかになった．
以上より，球半径と衝突速度にそれぞれ比例する抵抗力を新たに導入することで数値計算のエネルギー散逸をほぼ説明することができる．

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

４. まとめと今後の課題

以上の結果をまとめると次のようになる．粒子正面衝突の数値計算から得られた相互作用は従来の JKR 理論とおよそ一致することが示された．跳ね返りの際のエネルギー散逸については、低速度衝突では JKR 理論より大きく，高速衝突の場合は塑性変形により大幅に増大する．またエネルギー散逸機
構は従来の理論モデルと定性的に異なっており，接触している間に働くほぼ一定の抵抗力を新たに導入することにより説明される．正面衝突の他に、2 粒子をすべらせた場合、転がした場合、ねじれが起こる場合などの数値計算も進めている．
これらの計算から，すべり摩擦力やころがり抵抗力が得るこ 図２：２球相互作用の計算値と理論モデルの比とができる．このような結果をもとにして従来の理論モデル
較．灰色の実線は Maugis(1992)による理論にかわる新しい粒子相互作用モデルを構築することがことが線．点線はそれを横にずらしたもの．
この研究の目標である．

図３：様々な衝突速度に対する 2 球衝突計算から得られた跳ね返り定数．点線は JKR 理論で与えられる値．