目的

衝突クレーターの形成過程は天体の衝突・衝撃波の伝播（圧縮段階）から始まり、孔の掘削（掘削段階）を経てトランジェントクレーターの崩壊（緩和段階）に至る。本研究では特に掘削段階における掘削流や衝突放出物（イジェクタ）の放出過程などの素過程を明らかにするための手法として、Distinct Element Method（DEM, 離散要素法）による衝突の数値シミュレーションコードを開発し、粉体層への垂直衝突のシミュレーションを行なった。この手法は今後様々な問題に適用できると考えられるが、本研究ではその妥当性を示すことを中心に議論する。具体的には、イジェクタの速度分布や角度分布などについてこれまでの実験結果と本研究のシミュレーション結果とを比較・検討した。また、シミュレーションに基づく議論として、掘削流のモデルであるZ-モデルの妥当性の検証を行なった。

背景

これまで、イジェクタの放出過程を理解するために、粉体層（砂粒やガラス粒）への衝突実験が行われてきた。また、その実験結果を基にイジェクタ速度分布に関する幾つかのスケーリング則が得られている。しかしながら、室内実験においては、高速度で飛散する微小なイジェクタ粒子を精度良く観測することが難しく、不十分な点が多い。一方で、衝突の数値シミュレーションも近年盛んに行われるようになってきた。数値計算を行なう利点としては、例えば、実験では測定不可能な時空間分解能でデータが得られる、などがある。

しかしながら、従来の計算手法は、いずれも媒体を流体（連続体）近似するものである。個々のイジェクタ粒子の放出過程の解析のためには、連続体ではなく個々の離散的粒子の運動を扱い得ることが必要となる。

そこで、本研究では　粉体計算手法であるDEMを適用し、粉体層への衝突の数値シミュレーションコードの開発・応用を行った。DEMは主に粉体工学や土木工学の分野で開発・応用されてきた手法であるが、個々の粉体粒子の運動を離散的に取り扱うため、粉体層への衝突及びイジェクタの運動を的確にシミュレートし得ると期待される。DEMによる衝突クレータリングのシミュレーションを行なうのは、本研究が初の試みである。この数値シミュレーションによって得られるイジェクタ速度分布などの結果は、従来の砂への衝突実験の結果と直接的に比較でき、DEMによる数値計算の妥当性の検証が可能である。

DHMにおいては、各時間ステップごとに接触粒子による粒子間相互作用（弾性力・摩擦力）を着目粒子について求め運動方程式を解くことで個々の粒子の運動を計算していく。本研究では粒子間相互作用のモデルとしては、バネとダッシュポットが平行に繋がったフォークトモデルを用い、接線方向には摩擦スライダーも導入する（図1）。バネは弾性力を、ダッシュポットはエネルギー散逸を実現し、反発係数に応じた力を与える。また、摩擦スライダーは摩擦係数に応じた適切な摩擦力を実現する。

垂直衝突のシミュレーション設定

半径1mmの同径球（石英の物性値を使用）約38万個をターゲットとして直方体の箱(20cm×20cm×7cm)の中に堆積させる（図2a）。これにターゲット粒子の2〜4倍の大きさの粒子（アルミの物性値を使用）を様々な衝突速度(100〜900m/s)で衝突させる。粒子間には固着力・転がり抵抗は考えない。また、壁での反射波の発生を極力抑えるために、壁の反発係数は0に設定する。衝突粒子の大きさや衝突速度の他に、粒子間の反発係数や摩擦係数をパラメーターとして振って計算した。DEMの不利な点として、多くの粒子を扱うためには高性能の計算機が要求されるが、本研究では並列化して計算を行なった。それでも計算には時間がかかるため、トランジェントクレーターがほぼ形成された段階で計算を終了させる。衝突粒子半径を3mm、衝突速度を300m/s、粒子間の反発係数及び摩擦係数を各々0.4、0.5とした場合のシミュレーションで得たトランジェントクレーターの断面図を図2bに示す。

シミュレーション結果と実験結果の比較

シミュレーションで得られた、1）イジェクタ放出位置と放出速度の関係、2）イジェクタ放出位置と放出角度の関係、3）イジェクタ放出速度と放出量の関係、について各々実験結果と比較を行なった。ここでは　3）イジェクタ放出速度と放出量の関係、について示す。

イジェクタの放出速度と放出量の関係を示すものとして、ある速度以上で放出されたイジェクタの累積体積の分布をプロットしたものが、図3である。放出速度・体積は各々クレーター半径を用いて規格化している。連続分布として得られたシミュレーション結果は、低速領域・高速領域ともに、実験結果と調和的である。高速領域における大きな傾きから、放出速度にカットオフがあることが示唆される。

このように、本研究のシミュレーションによるイジェクタの速度分布等は実験結果と概ね調和的であり、このシミュレーションが掘削流やイジェクタ放出過程を明らかにする上で有効な手段となり得ることが示された。また、クレーター半径を用いて規格化したイジェクタ速度分布においては、粒子間の反発係数や摩擦係数の依存性はみられなかった。これは、掘削流が隣接粒子間の相対速度が小さい「粒子群」として振舞うことを示唆している。

Z-モデルの妥当性の検証

Maxwell らによって提唱されたZ-モデルは、衝突クレーター形成時の掘削流を解析的に表現するモデルである。定常非圧縮流であることや粒子速度の動径成分が距離の-Z乗で減衰するという簡単な仮定から、ある点（流源）からの流線が与えられ、掘削領域の推定など様々なクレーターの解析に用いられてきた。しかしながら、現実のクレータリングは非定常であるなど、Z-モデルが掘削流を適切に表現できるか疑問である。これまで実験や数値計算においてこのZ-モデルの妥当性が検証されてきたが、実験においてはターゲット内の運動を直接観察することは困難であり、数値計算においても有限差分による連続体近似計算のためにターゲット内粒子の運動追跡が困難である。したがって、これまでのところ、実際にターゲット内の粒子がZ-モデルで予想される流線に沿って運動する、ということは不明確であったが、本研究のDEMによる衝突シミュレーションによってこれを示すことを試みた。

ここでは、最も特徴的な結果を示す。図4は、表面を越えた粒子、即ちイジェクタ粒子が初期にどの位置にあったかを示したもので、ターゲット表面で放出された位置別に色分けして表示した。同じ色の粒子はすべて表面の同じ位置から放出されたことを意味する。これにZ-モデルの流線を重ねて見ると、粒子の色別分布とZ-モデルの流線が良く一致していることが分かる。このことは、ターゲット粒子がZ-モデルで予想された流線に沿って運動していることを示唆している。このように、ある程度空隙のある粉体層においても、Z-モデルは掘削流の簡便かつ良い近似であることが確認された。これは、粒子に変形がなく、「粒子群」として運動するため、おおよそ非圧縮媒体として振舞う結果であると考えられる。

結論

本研究においては衝突クレータリング、特にイジェクタ放出過程解明のためにDEMによる数値シミュレーションコードを開発し、粉体層への垂直衝突のシミュレーションを行ない、イジェクタ速度分布などに関して実験結果と比較した。その結果、DEMによる衝突シミュレーションは、実験結果と調和的であり、掘削流やイジェクタ放出過程の解析に有用であることが示された。また、掘削流は、概ねZ-モデル的流線に沿うことが確かめられた。今後、この手法は斜め衝突など様々な問題への適用が期待される。

接触粒子間の相互作用力モデル。(a)法線方向。(b)接線方向。

300m/sで衝突した際のクレーター形成過程のスナップショットを断面図で表示したもの。(a)初期状態。(b)トランジェントクレーター形成時（衝突後約40ms）。

ve以上の速度で放出されたイジェクタの体積V(>ve)の分布。クレーター半径R, 重力加速度gを用いて規格化した値でプロットした。赤線で表示されているのが本研究のシミュレーション結果（パラメーター：衝突速度300m/s, 衝突粒子半径3mm, 反発係数0.4, 摩擦係数0.5）。高速領域で途切れているところは、粒子1個の体積に相当する。一方実験結果は点で表示。黒線は低速領域の実験結果に対するフィッティングライン。

ターゲット表面を越えた粒子の初期位置を放出位置別に色分けして表示したもの（パラメーター：衝突速度300m/s, 衝突粒子半径3mm, 反発係数0.4, 摩擦係数0.5）。凡例にあるように、トランジェントクレーター半径R（この場合、4.95cm）で規格化した衝突点からの距離xの値に応じて色分けしている。重ねて表示されている曲線は、Z=3で流源がターゲット表面にある場合のZ-モデルの流線。上端の半円は比較のために衝突粒子を表示したもの。

本研究は天体衝突によるクレーター形成過程、特に掘削段階とよばれる過程の素過程を明らかにするための手法として、Distinct Element Method（離散要素法、以下ではDEMと省略する）による衝突の数値シミュレーションコードを開発し、粉体層への垂直衝突のシミュレーションを行ない、その結果について検討したものである。DEMによる衝突のシミュレーションは本研究によってはじめて行われたものである。

第1章では、クレーター形成、特に掘削段階の理解について過去の研究成果のレビューを行いその問題点が指摘されている。衝突クレーターの形成過程は天体の衝突・衝撃波の伝播（圧縮段階）から始まり、孔の掘削（掘削段階）を経てトランジェントクレーターの崩壊（緩和段階）に至る。このうち掘削段階とそれに伴うイジェクタの放出は、衝突クレーター形成の最も主要な過程ともいえる。これまで、イジェクタの放出過程を理解するために、粉体層（砂粒やガラス粒）への衝突実験が行われてきた。また、その実験結果を基にイジェクタ速度分布に関する幾つかのスケーリング則が得られている。しかしながら、室内実験においては、高速度で飛散する微小なイジェクタ粒子を精度良く観測することが難しく、不十分な点が多い。一方で、衝突の数値シミュレーションも近年盛んに行われるようになってきた。数値計算を行なう利点としては、例えば、実験では測定不可能な時空間分解能でデータが得られる、などがある。しかしながら、従来の計算手法は、いずれも媒体を流体（連続体）近似するものである。このことから、第1章では、個々のイジェクタ粒子の放出過程の解析のためには、連続体ではなく個々の離散的粒子の運動を扱うことが必要であることが述べられる。最後に、DEMは個々の粉体粒子の運動を離散的に取り扱うため、粉体層への衝突及びイジェクタの運動を的確にシミュレートし得ると期待できることが述べられる。

第2章では、本研究で用いられるDEMが解説される。DEMは主に粉体工学や土木工学の分野で開発・応用されてきた手法である。DEMにおいては、各時間ステップごとに接触粒子による粒子間相互作用（弾性力・摩擦力）を着目粒子について求め運動方程式を解くことで個々の粒子の運動を計算していく。本研究では粒子間相互作用のモデルとしては、バネとダッシュポットが並行に繋がったフォークトモデルを用い、接線方向には摩擦スライダーも導入されている。バネは弾性力を、ダッシュポットはエネルギー散逸を実現し、反発係数に応じた力を与える。また、摩擦スライダーは摩擦係数に応じた適切な摩擦力を実現する。なお計算法の詳細は Appendi x A-Eでも述べられている。

第3章は本論文の主要部分である。この章では、粉体層への衝突の数値シミュレーション実験の詳細とその結果が述べられる。半径1mmの同径球（石英の物性値を使用）約38万個をターゲットとして直方体の箱(20cm×20cm×7cm)の中に堆積させる。これにターゲット粒子の2〜4倍の大きさの粒子（アルミの物性値を使用）を様々な衝突速度(100〜900m/s)で衝突させる。粒子間には固着力・転がり抵抗は考えない。また、壁での反射波の発生を極力抑えるために、壁の反発係数は0に設定する。衝突粒子の大きさや衝突速度の他に、粒子間の反発係数や摩擦係数をパラメーターとして振って計算している。粉体層モデルの生成は Appendix Fで詳述される。DEMの不利な点として、多くの粒子を扱うためには高性能の計算機が要求されるが、本研究では並列化して計算を行なっている。それでも計算には時間がかかるため、トランジェントクレーターがほぼ形成された段階で計算を終了させている。

次に、イジェクタ速度分布などの結果を、従来の砂への衝突実験の結果と直接的に比較し、DEMによる数値計算の妥当性の検証を行っている。ここでは、1）イジェクタ放出位置と放出速度の関係、2）イジェクタ放出位置と放出角度の関係、3）イジェクタ放出速度と放出量の関係、について各々実験結果と比較を行なった。これらの数値計算の結果は、実験結果と概ね調和的であることが示されている。

従来の実験結果では、イジェクタ放出速度と放出量の関係は、低速の領域、高速の領域で結果が連続的につながらないことが知られていた。低速のイジェクタの計測と高速のイジェクタの計測を同時に行うことの困難もあり、中間領域がどのようにつながるかはよく分かっていなかった。本研究では連続分布としてイジェクタ放出速度と放出量の関係が得られ、その結果は、低速領域・高速領域ともに、実験結果と調和的であり、中間領域がどのようにつながるかについての示唆が得られた。また、高速領域における大きな傾きから、放出速度にカットオフがあることが示唆された。このカットオフは室内実験では観察が困難な衝突点付近での現象に起因すると考えられる。原理的にはDEMは衝突点付近の現象を明らかにすることができるため、シミュレーション結果が室内実験の結果と整合的であることも相まって、DEMシミュレーションが掘削流やイジェクタ放出過程を明らかにする上で有効な手段となり得ることが示された。

なお、DEM実験のモデルパラメターとして導入された、粒子間の反発係数や摩擦係数は結果に大きな影響を与えていないことが示された。弾性定数が大きくなるほどトランジェントクレーターの半径は大きくなる傾向は見いだされたが、これもあまり大きな影響ではない。このように個々の粒子の物性が大きな影響をもたないことは、掘削流が隣接粒子間の相対速度が小さい「粒子群」として振舞うことを示唆している。

第4章ではMaxwell らによって提唱されたZ-モデルとシミュレーションの結果が比較された。Z-モデルは、衝突クレーター形成時の掘削流を解析的に表現するモデルである。定常非圧縮流であることや粒子速度の動径成分が距離の-Z乗で減衰するという簡単な仮定から、ある点（流源）からの流線が与えられ、掘削領域の推定などに用いられる。非常に単純なモデルでありながら、その結果は室内実験をよく説明することが知られている。しかしながら、現実のクレータリングは非定常であるなど、Z-モデルが掘削流を適切に表現できているか疑問である。これまで実験や数値計算においてこのZ-モデルの妥当性が検証されてきたが、実験においてはターゲット内の運動を直接観察することは困難であり、数値計算においても有限差分による連続体近似計算のためにターゲット内粒子の運動追跡が困難である。したがって、これまでのところ、実際にターゲット内の粒子がZ-モデルで予想される流線に沿って運動しているか否か、ということは不明確であった。本研究のDEMによる衝突シミュレーションによって得られた粒子の流跡線とZ-モデルの流線を比較してみた結果、両者は良く一致していることが示された。このことは、ターゲット粒子がZ-モデルで予想された流線に沿って運動していることを示唆している。この場合の非定常性は流線に沿った速度が一定ではないことだけに現れていることが示唆される。また、単位体積あたりの粒子数から求めた密度も、誤差の範囲内で一定とみなせることが示された。これらの結果は、ある程度空隙のある粉体層においても、非圧縮流を仮定するZ-モデルが掘削流の簡便かつ良い近似であることを示している。これは、粒子に変形がなく、「粒子群」として運動するため、おおよそ非圧縮媒体として振舞う結果であると考えられる。

本研究は衝突クレータリング、特にイジェクタ放出過程解明のために、はじめてDEMによる数値シミュレーションコードを開発し、粉体層への垂直衝突のシミュレーションを行なったものである。イジェクタ速度分布などに関して実験結果と比較して、DEMによる衝突シミュレーションは、実験結果と調和的であることを示し、DEMの計算結果の妥当性を示すことに成功した。さらに、室内実験では困難な、掘削流やイジェクタ放出過程の解析に有用であることを示した。また、掘削流は、概ねZ-モデル的流線に沿うことが確かめられた。今後、この手法は斜め衝突など様々な問題への適用が期待される。将来的な応用が広く、発展性が高く、価値が高い研究であるといえる。

本論文は全体として松井孝典博士、千秋博紀博士との共同研究であるが、論文提出者が主体となって数値計算コードの開発、計算および結果の解析を行ったものであって、論文提出者の寄与が充分であると判断する。

したがって、博士（理学）を授与できると認める。