Research

個別要素法DEM
(Discrete Element Method/Distinct Element Method)

粉体現象をシミュレーションするには、個別要素法（または離散要素法）と呼ばれる数値計算手法がよく使われます。粉体を構成している個々の粒子の運動を解くことにより全体の挙動を表現し、重力などの外力と粒子間力の合力を計算して粒子に働く力を求め、時間積分して粒子の位置と速度を更新して行きます。粒子間相互作用は図１のようにバネとダッシュポットにより表現され、図２のように粒子同士が接触している場合にのみ働きます。バネの力は粒子同士の食い込み深さに比例し、ダッシュポットは粒子間の相対速度に比例した減速力を粒子に与えます。粒子に働く力は法線方向に働く反発力と、接線方向に働く摩擦力に分けることができ、接線方向のスライダーは、粒子の摩擦係数に応じて接線方向の力の上限が与えられています。

図１　DME粒子間相互作用モデル

図２　粒子の接触

粒子iが粒子jから受ける力のベクトル は以下のように計算されます。

ここで、法線方向と接線方向の力のベクトルをそれぞれ  ,  、弾性係数 ,  、粘性係数を ,  、バネの圧縮量のベクトルを  ,  、相対変位増分ベクトルを  ,  、時間刻みを  としています。粒子iの半径を すると、接線方向の力による粒子に加わるモーメントのベクトル は次のようになります。

粉体の一つの粒子は一度に複数個の周囲の粒子と衝突する可能性があり、粒子iと接触している全ての粒子jに関して力とモーメント ,  を計算し、粒子iに作用する力とモーメントの合力ベクトルを求め、粒子iの並進運動と回転運動の運動方程式を解きます。全ての粒子に対して運動方程式を解き、速度と位置、角速度と角度の更新を行います。

図3　多数の粒子との相互作用

クオータニオン（四元数）

　3次元空間での回転の表現形式として、基準座標系 ,  ,  軸周りの回転角 ,  ,  により回転を表すEuler角があります。Euler角は直感的に理解しやすいですが、有効な自由度や次元が失われる特異点があり、正確に回転を表現できません。そこで、クオータニオンを用いて粒子の回転角の管理を正確に行います。クオータニオンは1つの実数部分と3つの虚数部分からなる4元数です。クオータニオンを回転の表現に使うことで、Euler角で発生する特異点の問題を回避し、厳密に回転を管理することが可能です。粒子の角速度からクオータニオンの時間変化を計算し、クオータニオンの掛け算により更新を行います。クオータニオンの4つの成分からは粒子の回転を直感的に理解しづらいので、回転行列を介してクオータニオンからオイラー角を取り出して結果を出力します。