233 
衝突によって放出される粒子はイジェクタの良いアナロジーになると考えられる. 以上のような背景から粉体層への衝突実験が盛んに行われているわけであるが,微惑星の衝突・集積によって進化してきた太陽系天体一瞬(＜＜ 1 秒) で終わってしまう衝突過程にを考えるうえで衝突過程が重要であることは論を待たおける細かい粒子の運動の測定には困難さがついて回ない. 固体天体上には様々な大きさ・形状の衝突クレーる. また, 素過程を理解するためには重力や空隙率とターが見られ, また最近の小惑星・彗星探査によっていったパラメーターを自由にコントロールしたいがそも重力の小さな小天体上のクレーターも詳細に観測でれも難しい. そういった困難さを補い衝突過程を理解きるようになってきた. クレーター形成過程やそれにするために, 数値シミュレーションを行うことは有効伴う放出物 (イジェクタ) の放出・堆積過程を理解すな手段である. 数値計算ならば基本的に任意のデータることは, クレーター地形の形成過程を理解するのみが取得可能であり, パラメーターも自由に与えられる.

ならず, 衝突天体や標的天体の内部構造を明らかにす

近年のコンピューターの計算能力の飛躍的な向上にることにも通じ, 太陽系天体の進化を考える上で重要

よって, 衝突クレーター形成過程の数値シミュレーショである.

ンが盛んに行われている. ただし, 従来の数値シミュ

これまで, 衝突クレーター形成過程やイジェクタ放レーションは, 媒体を連続体として計算する流体計算

出過程を理解するために有効な方法として, 砂やガラ

コード (いわゆるハイドロコード) によるものであっス粒といった粉体層への衝突実験が行われてきた(例た [3]. 粉体粒子の挙動をシミュレートするためには,えば[1, 2]). 粉体層をターゲットに用いることには粒子の離散性 (接触していた粒子が離れた途端に力が以下の意味合いがある. 一つには, 粉々に破砕された働かなくなる) や回転, それから粒子間の摩擦や空隙岩石片が降り積もったものであるレゴリス層への衝突を的確に表現できなくてはならず, こういったことはを直接模擬できるということである. またもう一つにハイドロコードでは困難である.

は, 物質強度より重力が支配的になると考えられていハイドロコードでは扱いが難しい粉体の挙動解析をる比較的大きな衝突過程を理解すること (スケーリンするシミュレーション手法として, Distinct Elementグ則の確立) に繋がると期待される. これは, 乾いたMethod (日本語では個別要素法, 離散要素法, 粒子砂では強度はほぼ無いに等しいために, いわゆる重力要 素 法 , な ど と 呼 ば れ る . ま た 広 義 に は Discrete

支配領域の衝突過程を模擬できると考えられているか

Element Method と呼ばれる. いずれにしても DEMらである. さらに, イジェクタは衝突の際生じる衝撃と略される.) が知られている. DEM は, Cundall [4]波によって破壊された細かな岩石片からなるため (こによって提唱された数値シミュレーション手法であり､れがレゴリス層を形成するのであるが), 粉体層への粉体工学や土質力学など広く工学の分野で70年代から

234

日本惑星科学会誌Vol.13.No.4,2004

開発・発展してきた[5]. 地球科学では, 土石流やがタ速度分布) と比較することによって, シミュレー

け崩れなどのモデルに用いられ[6], また惑星科学にションの妥当性を検証することができる.

おいても, ダストの付着成長のモデルに用いられてき

ただし, デメリットももちろんあり, その最大の懸た[7]. DEM は, 単純に個々の粒子の運動を計算し念は計算機の能力である. 多数の粒子を扱うためにはていくことで粉体の挙動をシミュレートするものであ

大量のメモリと高速の演算能力が必要である. また粒る. したがって粉体層への衝突およびイジェクタの運子にサイズ分布があったり形状が球ではない場合には,動を的確にシミュレートし得ると期待される.

接触粒子のペアを探すのにより多くの時間が必要とな

そこで, 我々は DEM による粉体層への衝突の数値る. 我々は, 扱う粒子を衝突粒子を除いて全て同径球シミュレーションを試みている. 以下では, 用いていとし, さらに計算を並列化して行うことで3次元・る DEM モデル及び衝突シミュレーションの詳細を示

384,000個の粒子からなるターゲットへの衝突計算をすと共に, シミュレーションの結果がこれまで行われ可能にしたが, 今後のコードの改良と計算機能力の向てきた室内実験の結果と調和的であることを示す. さ

上によってさらに多くの粒子を扱えるようにすることらに, シミュレーションを用いて掘削流のモデルであが望まれる.

る モデルを検証した結果にもふれ, このシミュレーションの有用性を示す.

以下では我々のシミュレーションにおいて用いた接触粒子間相互作用力モデルについて述べる ([5]参照).
粒子間相互作用力モデルとして, バネとダッシュポッ

トが平行に繋がったいわゆるフォークトモデルを用い,接線方向には摩擦スライダーも導入する (図1). バネ

DEM においては, 各時間ステップごとに接触粒は粒子間に働く弾性力を表現し, ダッシュポットは非子のペアを探し, 接触粒子による粒子間相互作用力弾性衝突によるエネルギー散逸を実現するものである.

を求め, それを基に各粒子の運動方程式を解くこと

また, 接線方向の作用力によって摩擦力が働くことにで粒子の運動を計算していく. このように手続きは極なるが, 摩擦スライダーの導入によりクーロンの摩擦めて単純であるが､ 粉体挙動をシミュレートする上で則を表現するようになっている. 粒子は変形しない剛次のようなメリットがある.

体であるが, 接触粒子間では少量の重なりを許し・当然ながら, まさに粒子一つ一つの運動を直接計算(soft sphere model), その重なりの変位量に応じてバすることで粒子の離散性を表現することが可能であネの弾性力を与える. またダッシュポットによる粘性る.

抵抗は, 粒子の相対速度に比例して与える (ダッシュ・接触粒子間の相互作用を直接計算するので, 粉体層ポットはエネルギー散逸を現象論的に表現するためにを記述する構成方程式がいらない.
・個々の粒子にかかる力のモーメントを計算することで, 粒子にかかる摩擦力や粒子の回転運動を計算できる.
・粉体層内の空隙を直接表現できるので, 空隙率の影響を考慮したシミュレーションが可能である.
・サイズ分布・形状分布を持った粒子群を取り扱うことが可能である.
・粒子の挙動を直接室内実験の結果 (例えばイジェク

DEMによる粉体衝突シミュレーション／和田, 千秋, 松井

便宜上導入するもので, その ｢粘性｣ は物性値として235と等しいと仮定して､ をの粘性と必ずしも一致しない). 以上から, 接触粒子

間に働く法線方向の力 
 と接線方向の力  を次のよ ( )で与える.

うに与えることになる.

＝ − 
 − η (1)
＝ − 
 − η (2)

ダッシュポットの粘性係数は, 線形の減衰振動方程式




( )

ここで, はバネの弾性定数, ηはダッシュポットのの解を利用して与える. ここで,粘性係数. は接触粒子の重なり距離 (変位) で, シュポットで繋がれた二つの接触球の換算質量と相対は相対速度である. (注意:重なり変位は接触期間中積分される. これは後に述べる接線方向の摩擦力を実現する上で重要.) 各変数の添字 と はそれぞれ法線方向と接線方向の変数であることを示す. 接触している全ての粒子から受けるこれらの力を求め, 並進および回転の運動方程式を解くことによって粒子の運動を計算することができる. 粒子間に働くファンデルワー
ルス力などの付着力は, 該当する項を式   に加えとはバネとダッ変位である. 衝突時に自然長だったバネが縮み再び自然長に戻るまでを衝突による接触期間と考えると, 自然長における速度の比から反発係数εが衝突速度によらず以下のように与えられる.
 exp    (7)

これをηについて解くと,
  (8)ることでモデルに入れることができる. しかしここでは, 話を簡単化するために付着力は導入しない. また,

このようにεを用いてηを決めるため, 本研究の接触粒子形状の非対称性 (真球からのずれ) に由来する回粒子間相互作用力モデルは, 与えた反発係数に応じた転抵抗も考慮しない.

非弾性衝突を実現することになる.

法線方向のバネの弾性定数は, 半径 の2つの摩擦力は接線方向の力で与えられるが, この際クー粒子 , が接触しているとき, ヘルツの弾性接触理

ロンの摩擦則を満たすべく摩擦係数μを導入する. 具論 (例えば[8]) によって以下のように与える.

体的には, 接線方向の弾性力が法線方向の弾性力のμ    (3)倍を超えた場合, 接線方向弾性力を法線方向弾性力のμ倍に留めるということにする. モデルとしては, こここで, , は, ポアソン比 とヤング率 を用のような ｢すべり｣ を許す摩擦スライダーを導入した

いてことになる. 粒子間に接線方向の相対速度があってか  , (4)つ接線方向の弾性力が十分大きければ, 摩擦係数が一定という動摩擦のよいモデルとなっている. また, 静と表される物性値による定数. また は粒子 , 間止状態でも接線方向の弾性力があれば摩擦が働くのでに働く弾性力である.  は各時間ステップにおいて,静止摩擦も表現される. ただし, 静止摩擦係数と動摩一つ手前のステップにおける を用いて逐一計算さ擦係数の区別はつけていない. また, 摩擦係数は速度れる. このようにヘルツの理論を用いることで, によらず一定としてある.

はヤング率やポアソン比という既知の物性値であらわ

以上のように与えられる接触粒子間相互作用モデルされる. を与えるに際しては, 簡単のために  とにおいては, 結局のところ反発係数と摩擦係数がキー の比が縦弾性係数とせん断弾性係数の比

パラメーターであるといえる.

236

日本惑星科学会誌Vol.13.No.4,2004

 

ターが形成される (図3c) ことが見て取れる. また,まばらではあるが飛び出したイジェクタ粒子が逆円錐以下では, 衝突シミュレーションによるクレーター形のイジェクタカーテンを形成していることも見える.

形成過程の様子を示す. 初期条件として, 直方体の箱

なお, ターゲットの深さを倍にしたシミュレーション(20cm×20cm×7cm) の中に, 半径1mmの同径球 (石においてもクレーターの深さに変化は見られなかった

英の物性値を使用) 384,000個を自由落下で堆積させ,ため, クレーターの深さが決まる上でターゲットの底これをターゲットとする (図2). 空隙率は約43％とな14(a)り, これはゆるやかに詰めた砂と同程度である. このターゲットに半径3mm の粒子 (アルミの物性値を使用) を衝突速度300m/sで垂直に衝突させたときのクレーター形成過程の断面スナップショットを図3に示す.

7

ここでは粒子間の反発係数を0.4, 摩擦係数を0.5とした. また, 重力加速度は9.8m/s２ とした. ターゲットの壁面での反射波の発生を極力抑えるために, 壁面の反発係数は0に設定している. 本研究ではコードを並列化した上で海洋研究開発機構の大型計算機で計算を0-1014(b)010010

行っているが, それでも計算には時間がかかる (〜1週間) ため､ 崩壊する前の最大クレーターであるトランジェントクレーターがほぼ形成された段階で計算を7終了させる.

0-1014(c)trans ient crater diameter7

図3から, お椀型のクレーターが形成される様子が分かる. すなわち, 初期には孔が半球状に広がっていくが (図3a), やがて深さが最大に達し (図3b), その後は側面が広がって, ついにはトランジェントクレー

0-10010 !""#$%&'( )*+,-. /01234%5&'
6 789:;<=$>:?<@=$:ABCD
E&'F<>:G<H"$:IJKLJ)*+,<DEMによる粉体衝突シミュレーション／和田, 千秋, 松井の影響はないと考えられる.

237

ギーに関係した量である. 図4においてシミュレーション結果は, 乾いた砂 (Dry sand, Quartz sand) への衝



突実験で得られたスケーリング則[1, 9]と良く一致していることが分かる. 特にπスケーリングによれば,

この節では, シミュレーションによって得られる結冪の値はエネルギーの散逸の程度を表す[9]. このこ

果を室内実験の結果と比較した例を示す. 具体的にはとから我々のシミュレーションは, 形成されるクレークレーターの半径とイジェクタ放出速度分布の2つにターのサイズだけでなく, エネルギーの散逸過程に関ついての比較を示すが, いずれもシミュレーション結しても乾いた砂への衝突実験を再現できていることが果は実験結果と調和的である.

分かる.

次元解析から得られるπスケーリング[1]を用いて,イジェクタの放出速度と放出量の関係を示すものとクレーター半径の比較を行なう. 重力支配領域においして, ある速度以上で放出されたイジェクタの累積体ては､ 2つの無次元数π(ρ )とπ２ = 3.22積の分布をプロットしたものが, 図5である. 放出速 が冪乗則の関係にあるとされる. ここで,, ,度・体積は各々クレーター半径を用いて規格化してい,, , は各々クレーター半径, ターゲットのる. この2つの変数の間には冪乗則があると考えられバルク密度, 衝突粒子の質量, 衝突粒子の半径, 衝突ている[10]. 図5に示されるように, 連続分布として速度, 重力加速度である. π は規格化されたクレー得られたシミュレーション結果は, 低速領域・高速領ター半径を意味し､ π２は衝突粒子の初期運動エネル域ともに, 実験結果[2, 11, 12]と調和的である. また,  !"# $%&'$( )*0-60107809 -/70+ ,-.//01//02//345: -/;<=>? @ABAC  G.2;
1! H=>IJKLMNONP/DEF/_Q  RSTUVW5TXYZ.[\-/].F/G.];0 ^UG.]
G66

D7
F/
V5
TXY
ZE[\ -/
D
F/
5TXY
ZE[\ -/6 0`a6 b cade fghijPCklm nop q+rca< stm'uv 0wxyfg z{|}~}|= !"#//345
0-13309 $%&'$( )*+ , -1
/70: -/;<f P}|>0. op000>KL)*0=>fKL)*}ONPQ

1／3２238

日本惑星科学会誌Vol.13.No.4,2004

高速領域においては, Housen ら[10]が実験結果から得たスケーリング則から大幅に外れており, 衝突点近傍から放出される高速イジェクタに対しては別のスケーリング則の必要性が示唆される.

Maxwell[13]らによって提唱された モデルは,衝突クレーター形成時の掘削流を解析的に表現するモデルである. 定常非圧縮流であることや粒子速度の動径成分が距離の 乗 ( は定数) で減衰するという )#$&&)
  ! "#$%%&'(*#%+,)-#%+./+0123456789:  ;<=><?@AB C!DEF),+G.H&/6IJKLMNOP+QR STOPUVW) X#$6)YZ[^_`YV+abAcW\PEF X]
def01 +

簡単な仮定から, ある点 (流源) からの流線が与えられ, 掘削領域の推定など様々なクレーターの解析に用モデルの流線が良く一致していることが分かる. こいられてきた. しかしながら, 現実のクレータリングのことは, ターゲット粒子が モデルで予想されたは非定常であるなど, モデルが掘削流を適切に表流線に沿って運動していること, 即ち粉体層におけ現できるか疑問である. これまで実験や数値計算におる掘削流は モデルでよく表現されることを示していてこの モデルの妥当性が検証されてきたが, 実いる. 先に述べたように粉体層への衝突クレーター

験においてはターゲット内の運動を直接観察すること形成過程は重力支配域のクレーター形成過程を模擬は困難であり, 数値計算においても有限差分による連するものである. したがって, 本研究の結果は モ

続体近似計算のためにターゲット内粒子の運動追跡がデルが重力支配域における掘削流の簡便かつ良い近困難であった. したがって, これまでのところ, 実際似であることを示唆している.にターゲット内の粒子が モデルで予想される流線に沿って運動する, ということは不明確であった. そこで我々は, DEM による衝突シミュレーションによってこれを示すことを試みた[14]. イジェクタ放出過程

クレーター径は, 反発係数や摩擦係数にも依存する.

は, クレーター掘削流を反映したものであるから, 前

これは, 例えば反発係数を小さくとるとその分エネル節に示したようにイジェクタの速度分布が実験結果とギー散逸量が大きくなるためにクレーター径が小さく概ね調和的であることは, このシミュレーションによなる, といったことで説明される. しかしながら, 反る掘削流が妥当なものであることを保証していると考発係数および摩擦係数の値を大きく変えても (反発係えられる.

数を0.1〜0.8, 摩擦係数を0.1〜0.9), クレーター径のここでは, 最も特徴的な結果を示す. 図6は, 衝突変化は高々1割程度であった. また, クレーター半径前のターゲット表面の高さを越えた粒子, 即ちイジェを用いて規格化したイジェクタ速度分布においては,クタ粒子が初期にどの位置にあったかを示したもの

反発係数や摩擦係数の依存性はみられなかった. さらで, ターゲット表面で放出された位置別に色分けしに, 掘削流を見ても, 反発係数や摩擦係数によってそて表示した. 同じ色の粒子はすべて表面の同じ位置の形状はほとんど変化しなかった. 以上のことは, クから放出されたことを意味する. これに モデルのレーター形成過程が, 特に掘削段階において反発係数流線 (=3) を重ねて見ると, 粒子の色別分布と や摩擦係数にあまり依存しないことを示唆している.

DEMによる粉体衝突シミュレーション／和田, 千秋, 松井

239

従って, 掘削流は反発係数と摩擦係数が効かない程度

本研究のモデルでは反発係数と摩擦係数は先に述べに隣接粒子間の相対速度が小さい ｢粒子群｣ 流としてたように速度によらないが, 現実には速度依存性が見振舞うと考えられる. また, そのような粒子群流は非られる. より詳細なモデルが必要であれば, 反発係数圧縮流 ( モデルの仮定) と考えられ, モデル的と摩擦係数に速度依存性を考慮すべきかもしれない.

掘削流となることと符合する. ただし, 空隙率がこの現在のところ, 熱を含んだ状態方程式を考えておらシミュレーションで示したターゲットのものより大きず, 散逸したエネルギーの行方は問題にしていない.

くなるとクレーター形成時における圧縮の効果が現れまた物性値が常に一定であるため, 塑性変形や相変化[15], 非圧縮流ではなくなるだろう. またその際隣接は考慮されない. こういった制約から, 現段階では著粒子間の相対速度も大きくなり, 反発係数・摩擦係数しい塑性変形を引き起こさない程度の比較的低速度の依存性が顕著になると考えられる.

(数百 m/s) の衝突シミュレーションとなっている.
これでも, 特に掘削段階およびイジェクタ放出過程を 

明らかにする上で十分有効なシミュレーションコード
であると考えられるが, 接触力モデルの改良や個々の粒子に対する適切な状態方程式を導入することで, メ我々が行なった DEM による粉体層への衝突シミュレーションの結果は, 用いた DEM モデルが比較的簡ルトの生成といったことについても議論できる可能性があり, 今後の課題である.

単ではあるが, 実験結果と調和的であり, 掘削流やイ今後､ DEM による衝突シミュレーションは, 斜めジェクタ放出過程の解析に有用であることが示された.

衝突におけるクレーター形状やイジェクタ放出過程のここで, 本研究で用いた DEM モデルの改良すべき点解析や, 空隙率が高いと考えられる小天体上の衝突過を幾つか指摘しておく.

程の解析, さらには微惑星や小惑星の衝突進化過程

本研究では粒子形状を真球とし, 粒子同士の接触を点接触としている. しかし粒子形状の非対称性 (楕円(破壊・集積) のシミュレーションにも適用できると期待される.

体であったり凸凹していたり) は, 実効的な粒子間摩擦力や固着力を大きくすると予想される. 粒子形状の

非対称性を表現する方法としては回転抵抗を導入する
ことが考えられる.
粒子間付着力 (例えばファンデルワールス力) は,査読者である城野信一氏には大変有益なコメントを頂きました. 深く感謝いたします.

粒径が小さいと重力に比して効いてくるはずであるから, 付着力モデルの導入が必要である. 最近の著者ら

の研究によれば, 付着力を導入することで強度支配領域の衝突過程を理解することも可能であることが示唆されている[16]. また, 多数の粒子からなる粒子塊を[1] Schmidt, R. M., 1980, Proc. Lunar Planet. Sci.
Conf. 11, 2099.

衝突させその影響を見る場合には, その衝突粒子塊が[2] Yamamoto, S., 2002, Icarus 158, 87.

破壊される様を表現するためにも粒子間付着力モデル[3] Holsapple, K. A. and Housen, K. R., 2002, Proc.

を吟味の上導入する必要がある. これをさらに推し進Lunar Planet. Sci. Conf. 33, #1857.

めれば, 天体の衝突破壊をシミュレートすることにも[4] Cundall, P. A. and Strack, O. D. L., 1979,繋がると期待される.

G éotechnique 29, 47.

240

日本惑星科学会誌Vol.13.No.4,2004[5] 粉体工学会編, 1998, 粉体シミュレーション入門,
産業図書.
[6] Campbell, C. S. et al., 1995, J. Geophys. Res. 100,8267.
[7] Dominik, C. and Tielens, A. G. G. M., 1997,
Astrophys. J. 480, 647.
[8] Timoshenko, S. P. and Goodier, J. N., 1970,
Theory of Elasticity, 3rd edition, McGraw-Hill.
[9] Schmidt, R. M. and Housen, K. R., 1987, Int. J.
Impact Engng. 5, 543.
[10] Housen, K. R. et al., 1983, J. Geophys. Res. 88,2485.
[11] Andrews, J. L. B., 1975, Report AFWL-TR-74314.
[12] St ffler, D. et al., 1975, J. Geophys. Res. 80,4062.
[13] Maxwell, D. E., 1977, in Impact and Explosion
Cratering, 1003.
[14] Wada, K. et al., 2003, Proc. 36th ISAS Lunar
Planet. Symp., 53.
[15] Housen, K. R. and Holsapple, K. A., 2003, Icarus
163, 102.
[16] 和田浩二ほか, 2004, 日本惑星科学会秋季講演会予稿, 11.