個別要素法

個別要素法（こべつようそほう、英: Distinct Element Method、DEM）または離散要素法（Discrete Element Method、DEM）は、解析の対象を自由に運動できる多角形や円形・球の要素の集合体としてモデル化し、要素間の接触・滑動を考慮して、各時刻におけるそれぞれの要素の運動を逐次追跡して解析する手法である。もとは岩盤工学に適用するためにPeter A. Cundall (1971)およびCundall and Strack (1979)^[1]により発表された論文に端を発しており、現在は液状化や土石流など地盤の挙動解析やコンクリート構造物、粉体（化学工学、リチウムイオン電池、薬学、農学など）、磁気相互作用力を有する電子写真システムのトナーの挙動解析などに用いられている。

RotatingDrum

2粒子間距離を線で表現し, 力大きさを線の太さで表現したものをForce Chain 力鎖とよぶ.

概略

以下に、円形要素を用いた際の運動方程式を示す。

質量${\displaystyle m_{i}}$、慣性モーメント${\displaystyle I_{i}}$のある円形要素${\displaystyle i}$について、次の運動方程式が成り立つ。

${\displaystyle m_{i}{\ddot {\mathbf {x} }}+C_{i}{\dot {\mathbf {x} }}+F_{i}=0}$

${\displaystyle I_{i}{\ddot {\theta }}+D_{i}{\dot {\theta }}+M_{i}=0}$

ここに${\displaystyle F_{i}}$：要素に働く合力、${\displaystyle M_{i}}$：要素に働く合モーメント、${\displaystyle C_{i}}$、${\displaystyle D_{i}}$：減衰定数、${\displaystyle \mathbf {x} }$：要素の変位ベクトル、${\displaystyle \theta }$：要素の回転変位である。

要素同士が接触しているときは${\displaystyle F_{i}=K\mathbf {x} }$（${\displaystyle K}$は弾性定数）及び${\displaystyle M_{i}=Kr^{2}\theta }$（${\displaystyle r}$は円の半径）、離れているときは、${\displaystyle F_{i}=M_{i}=0}$で表される。ただし、重力を考える場合は合力の項で考慮する。

上式を数値積分することで、逐次変位ベクトルと回転変位を得ることができる。

特徴

• 計算コストは${\displaystyle {\mathcal {O}}(N)}$ もしくは${\displaystyle {\mathcal {O}}(NlogN)}$ で線形に近くスケールする。
• 明示的時間積分のため安定性条件は ${\displaystyle \Delta t<\pi {\sqrt {\frac {m}{k_{n}}}}}$ などにより制限される。
• 接触履歴を保持することで摩擦・転がり抵抗を自然に扱える。

ソフトウェア

オープンソース

• LIGGGHTS — LAMMPS を基盤に開発された粉体シミュレーション用オープンソース DEM コードで、粒子の熱挙動にも対応する。^[2]
• Yade（日本語版） — C++ と Python で実装された拡張性の高い DEM フレームワーク。CFD との連成（CFD–DEM）解析にも利用される。^[3]
• LAMMPS（granular パッケージ） — 汎用分子動力学コード LAMMPS に同梱される粒状体シミュレーション用サブパッケージで、軟球モデルによる粒子間接触力や境界条件が実装されている。^[4]

商用ソフト

• Granuleworks
• EDEM
• Rocky DEM

アルゴリズム

1. 初期化

• 粒子属性（径・質量・慣性モーメントなど）を読み込む。
• 計算領域をボックス、周期境界、固定壁などとして定義する。
• セル分割法を用いて「衝突探索ボックス」を生成し、全粒子を登録する。

2. 時間積分ループ

以下ではステップ番号を ${\displaystyle n}$、タイムステップを ${\displaystyle \Delta t}$ とする。

(a) 接触判定

近接セル内の候補ペアに対し、半径差 ${\displaystyle \delta _{n}=r_{i}+r_{j}-|\mathbf {x} _{i}-\mathbf {x} _{j}|}$ を評価し、${\displaystyle \delta _{n}>0}$ なら接触とみなす。

(b) 接触履歴

クーロン摩擦や粘弾性モデルでは接触履歴が必要となる。前ステップで既接触なら

${\displaystyle \delta _{t}^{\,n}=\delta _{t}^{\,n-1}+v_{t}^{\,n}\,\Delta t}$

未接触なら

${\displaystyle \delta _{t}^{\,n}=v_{t}^{\,n}\,\Delta t}$

(c) 接触力計算

Hertz-Mindlin 式などを用い、法線力 ${\displaystyle \mathbf {F} _{n}}$ と接線力 ${\displaystyle \mathbf {F} _{t}}$ を評価する。

(d) 外力計算

重力 ${\displaystyle m\mathbf {g} }$、流体抵抗、磁気力などシミュレーション対象に応じた外力を加算する。

(e) 運動方程式の陽積分

例として単純並進のみを示す。

${\displaystyle \mathbf {v} ^{\,n+1}=\mathbf {v} ^{\,n}+\mathbf {a} ^{\,n}\,\Delta t}$

${\displaystyle \mathbf {x} ^{\,n+1}=\mathbf {x} ^{\,n}+\mathbf {v} ^{\,n+1}\,\Delta t}$

角運動は陽オイラーや四元数を用いるのが一般的である。

(f) 粒子情報更新

位置・速度が確定したら、衝突探索ボックスのセルインデックスを更新する。

(g) 時間終了判定

${\displaystyle t_{\text{current}}\geq t_{\text{end}}}$ なら計算終了。それ以外は (a) に戻る。

疑似コード

initialize_particles()
initialize_search_grid()

for t = 0; t < t_end; t += Δt
    build_neighbor_list()
    for each particle i
        Fi, Mi ← 0
    for each contact pair (i, j)
        δn ← overlap(i, j)
        if δn > 0 then
            update_tangential_history(i, j)
            Fij, Mij ← contact_force(δn, history)
            Fi  += Fij
            Fj  -= Fij
            Mi  += Mij
            Mj  -= Mij
    for each particle i
        Fi += external_forces(i)
        vi ← vi + (Fi/mi) Δt
        ωi ← ωi + (Mi/Ii) Δt
        xi ← xi + vi Δt
        update_search_grid_cell(i)
end for

計算の詳細

セル分割法

セル分割法（英: *cell-linked list method*、*link-cell method*）は、短距離カットオフをもつ相互作用を扱う離散要素法（DEM）や分子動力学 (MD) シミュレーションで、近傍粒子対を線形時間で検索する空間分割アルゴリズムである。計算領域をカットオフ半径以上の立方体セルに区切り、各セル内に存在する粒子番号をリンクリストとして保持することで、粒子 ${\displaystyle N}$ 個に対するペア探索の計算量を ${\displaystyle {\mathcal {O}}(N)}$ に削減できる。

主な手順は以下のとおり^[5][6]。

1. セル分割 — シミュレーションボックスを辺長 ${\displaystyle \ell \geq r_{c}}$ の格子セルに分割し、セル番号 ${\displaystyle (i_{x},i_{y},i_{z})}$ を割り当てる。
2. リンクリスト構築 — 各粒子の座標から所属セルを計算し、セル毎に先頭ポインタ head と粒子間リンク lscl の単方向リストを更新する。
3. 近傍探索／力計算 — 各粒子について自身が属するセルと 26 個の隣接セル内の粒子だけを走査し、距離がカットオフ半径 ${\displaystyle r_{c}}$ 未満のペアに対して力を評価する（ニュートン第三法則を用いて重複計算を排除）。
4. セル更新 — 各タイムステップで粒子がセル境界を越えた場合、対応するリンクリストを再構築する。

アルゴリズムの平均計算量は

${\displaystyle {\mathcal {O}}(N\,\rho \,\ell ^{\,3})\simeq {\mathcal {O}}(N)}$

であり、セル体積 ${\displaystyle \ell ^{3}}$ と平均粒子密度 ${\displaystyle \rho }$ が一定なら粒子数に対して線形にスケールする。セル長 ${\displaystyle \ell }$ は密度、カットオフ半径、キャッシュライン長などの実装パラメータによって最適値が異なる。

物理量計算

球粒子の慣性モーメント

${\displaystyle I_{i}={\frac {2}{5}}\,m_{i}R_{i}^{2}}$

Hertz–Mindlin 接触

${\displaystyle F_{ij}^{C_{n}}=-{\frac {4}{3}}\,E_{ij}^{\ast }\,{\sqrt {R_{ij}^{\ast }}}\;\delta _{n}^{\,3/2}}$

${\displaystyle k_{t}={\frac {8}{2-\nu }}\,G_{ij}^{\ast }\,{\sqrt {R_{ij}^{\ast }\,\delta _{n}}}}$

（${\displaystyle E_{ij}^{\ast }}$, ${\displaystyle G_{ij}^{\ast }}$, ${\displaystyle R_{ij}^{\ast }}$ は有効ヤング率・有効せん断弾性率・有効半径）

粘性減衰と反発係数

${\displaystyle \eta =-2\ln e\,{\sqrt {\frac {m_{ij}^{\ast }\,k}{\pi ^{2}+(\ln e)^{2}}}}}$ ^[7]

数値積分の安定条件

${\displaystyle \Delta t\;\lesssim \;2{\sqrt {\frac {m}{k_{n}}}}}$ ^[8]

法線・接線成分の計算

単位法線ベクトル

${\displaystyle \mathbf {n} _{ij}={\frac {\mathbf {x} _{i}-\mathbf {x} _{j}}{|\mathbf {x} _{i}-\mathbf {x} _{j}|}}}$

相対並進速度

${\displaystyle \mathbf {v} _{ij}=\mathbf {v} _{i}-\mathbf {v} _{j}}$

法線速度

${\displaystyle \mathbf {v} _{n,ij}=(\mathbf {v} _{ij}\cdot \mathbf {n} _{ij})\,\mathbf {n} _{ij}}$

接線速度

${\displaystyle \mathbf {v} _{t,ij}=\mathbf {v} _{ij}-\mathbf {v} _{n,ij}+(R_{i}\mathbf {\omega } _{i}+R_{j}\mathbf {\omega } _{j})\times \mathbf {n} _{ij}}$

法線方向オーバーラップ

${\displaystyle \delta _{n,ij}=R_{i}+R_{j}-|\mathbf {x} _{i}-\mathbf {x} _{j}|}$

${\displaystyle {\boldsymbol {\delta }}_{n,ij}=\delta _{n,ij}\,\mathbf {n} _{ij}}$

接線方向オーバーラップ（連続表現）

${\displaystyle {\boldsymbol {\delta }}_{t,ij}(t)=\int _{t_{0}}^{t}\mathbf {v} _{t,ij}(\tau )\,\mathrm {d} \tau }$

接線方向オーバーラップ（離散更新）

${\displaystyle {\boldsymbol {\delta }}_{t,ij}^{(n+1)}=|{\boldsymbol {\delta }}_{t,ij}^{(n)}|\,\mathbf {t} _{ij}^{(n+1)}+\mathbf {v} _{t,ij}^{(n+1)}\Delta t}$

ここで ${\displaystyle \mathbf {t} _{ij}^{(n+1)}}$ は接線方向の単位ベクトルであり，

• ${\displaystyle \mathbf {v} _{t,ij}^{(n+1)}\neq \mathbf {0} }$ の場合

${\displaystyle \mathbf {t} _{ij}^{(n+1)}={\frac {\mathbf {v} _{t,ij}^{(n+1)}}{|\mathbf {v} _{t,ij}^{(n+1)}|}}}$

• ${\displaystyle \mathbf {v} _{t,ij}^{(n+1)}=\mathbf {0} }$ の場合

${\displaystyle \mathbf {t} _{ij}^{(n+1)}={\frac {{\boldsymbol {\delta }}_{t,ij}^{(n)}}{|{\boldsymbol {\delta }}_{t,ij}^{(n)}|}}}$

角速度の分解（粒子 i）

${\displaystyle {\boldsymbol {\omega }}_{n,i}=({\boldsymbol {\omega }}_{i}\cdot \mathbf {n} _{ij})\,\mathbf {n} _{ij}}$

${\displaystyle {\boldsymbol {\omega }}_{t,i}={\boldsymbol {\omega }}_{i}-{\boldsymbol {\omega }}_{n,i}}$

接触モデル

線形ばね–ダッシュポット（Voigtモデル）

${\displaystyle \mathbf {F} _{ij}^{\rm {C}}=\mathbf {F} _{ij}^{C_{n}}+\mathbf {F} _{ij}^{C_{t}}}$

${\displaystyle \mathbf {F} _{ij}^{C_{n}}=-k_{n}\,{\boldsymbol {\delta }}_{n}-\eta _{n}\,\mathbf {v} _{n}}$

${\displaystyle \mathbf {F} _{ij}^{C_{t}}=-k_{t}\,{\boldsymbol {\delta }}_{t}-\eta _{t}\,\mathbf {v} _{t}}$

摩擦モデル

${\displaystyle {\bigl |}\mathbf {F} _{ij}^{C_{t}}{\bigr |}\;\leq \;\mu \,{\bigl |}\mathbf {F} _{ij}^{C_{n}}{\bigr |}}$ 条件を超えた場合は

${\displaystyle \mathbf {F} _{ij}^{C_{t}}=-\mu \,{\bigl |}\mathbf {F} _{ij}^{C_{n}}{\bigr |}\,{\frac {\mathbf {v} _{t}}{|\mathbf {v} _{t}|}}}$

転がり抵抗

${\displaystyle \mathbf {M} _{ij}=-\mu _{r}\,R_{ij}^{\ast }\,{\bigl |}\mathbf {F} _{ij}^{C_{n}}{\bigr |}\,{\frac {{\boldsymbol {\omega }}_{i}-{\boldsymbol {\omega }}_{j}}{|{\boldsymbol {\omega }}_{i}-{\boldsymbol {\omega }}_{j}|}}}$ ^[9]

運動方程式

${\displaystyle m_{i}\,{\frac {{\rm {d}}\mathbf {v} _{i}}{{\rm {d}}t}}=\sum _{j}\mathbf {F} _{ij}^{\rm {C}}+m_{i}\mathbf {g} }$

${\displaystyle I_{i}\,{\frac {{\rm {d}}{\boldsymbol {\omega }}_{i}}{{\rm {d}}t}}=\sum _{j}\mathbf {T} _{ij}}$

積分

粒子 i の並進速度 ${\displaystyle {\vec {v}}_{i}}$，位置 ${\displaystyle {\vec {x}}_{i}}$，角運動量 ${\displaystyle {\vec {L}}_{i}=I_{i}{\vec {\omega }}_{i}}$ を 時間ステップ ${\displaystyle n\to n+1}$ に進める代表的な数値積分法を示す。 ${\displaystyle {\vec {a}}_{i}=m_{i}^{-1}\sum _{j}{\vec {F}}_{ij}^{C}+{\vec {g}}}$， ${\displaystyle {\vec {\tau }}_{i}=\sum _{j}{\vec {T}}_{ij}}$ はそれぞれ並進加速度と外力によるトルクである。

オイラー前進法（Explicit Euler）

• 速度更新 ${\displaystyle \displaystyle {\vec {v}}_{i}^{n+1}={\vec {v}}_{i}^{n}+\Delta t\,{\vec {a}}_{i}^{n}}$
• 位置更新 ${\displaystyle \displaystyle {\vec {x}}_{i}^{n+1}={\vec {x}}_{i}^{n}+\Delta t\,{\vec {v}}_{i}^{n}}$
• 角運動量更新 ${\displaystyle \displaystyle {\vec {L}}_{i}^{n+1}={\vec {L}}_{i}^{n}+\Delta t\,{\vec {\tau }}_{i}^{n}}$

セントラルディファレンス法（中央差分）

• 位置更新 ${\displaystyle \displaystyle {\vec {x}}_{i}^{\,n+1}=2{\vec {x}}_{i}^{\,n}-{\vec {x}}_{i}^{\,n-1}+(\Delta t)^{2}{\vec {a}}_{i}^{\,n}}$
• 速度評価 ${\displaystyle \displaystyle {\vec {v}}_{i}^{\,n}={\frac {{\vec {x}}_{i}^{\,n+1}-{\vec {x}}_{i}^{\,n-1}}{2\Delta t}}}$
• 角運動量更新 ${\displaystyle \displaystyle {\vec {L}}_{i}^{\,n+1}=2{\vec {L}}_{i}^{\,n}-{\vec {L}}_{i}^{\,n-1}+(\Delta t)^{2}{\vec {\tau }}_{i}^{\,n}}$

シンプレクティック・オイラー法（Semi-implicit Euler）

• 速度更新 ${\displaystyle \displaystyle {\vec {v}}_{i}^{n+1}={\vec {v}}_{i}^{n}+\Delta t\,{\vec {a}}_{i}^{n}}$
• 位置更新 ${\displaystyle \displaystyle {\vec {x}}_{i}^{n+1}={\vec {x}}_{i}^{n}+\Delta t\,{\vec {v}}_{i}^{n+1}}$
• 角運動量更新 ${\displaystyle \displaystyle {\vec {L}}_{i}^{n+1}={\vec {L}}_{i}^{n}+\Delta t\,{\vec {\tau }}_{i}^{n}}$

リープフロッグ法（Velocity Verlet 型）

• 速度半ステップ ${\displaystyle \displaystyle {\vec {v}}_{i}^{\,n+{\tfrac {1}{2}}}={\vec {v}}_{i}^{\,n-{\tfrac {1}{2}}}+\Delta t\,{\vec {a}}_{i}^{\,n}}$
• 位置更新 ${\displaystyle \displaystyle {\vec {x}}_{i}^{\,n+1}={\vec {x}}_{i}^{\,n}+\Delta t\,{\vec {v}}_{i}^{\,n+{\tfrac {1}{2}}}}$
• 角運動量半ステップ ${\displaystyle \displaystyle {\vec {L}}_{i}^{\,n+{\tfrac {1}{2}}}={\vec {L}}_{i}^{\,n-{\tfrac {1}{2}}}+\Delta t\,{\vec {\tau }}_{i}^{\,n}}$

4 次 Runge–Kutta 法（RK4）

内部ステージ

${\displaystyle {\begin{aligned}k_{1}^{v}&={\vec {a}}_{i}^{\,n},&k_{1}^{x}&={\vec {v}}_{i}^{\,n},&k_{1}^{L}&={\vec {\tau }}_{i}^{\,n},\\[2pt]k_{2}^{v}&={\vec {a}}_{i}\!\left({\vec {x}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{1}^{x},\;{\vec {v}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{1}^{v}\right),\\k_{2}^{x}&={\vec {v}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{1}^{v},\\k_{2}^{L}&={\vec {\tau }}_{i}\!\left({\vec {x}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{1}^{x},\;{\vec {v}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{1}^{v}\right),\\[2pt]k_{3}^{v}&={\vec {a}}_{i}\!\left({\vec {x}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{2}^{x},\;{\vec {v}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{2}^{v}\right),\\k_{3}^{x}&={\vec {v}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{2}^{v},\\k_{3}^{L}&={\vec {\tau }}_{i}\!\left({\vec {x}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{2}^{x},\;{\vec {v}}_{i}^{n}+{\tfrac {\Delta t}{2}}k_{2}^{v}\right),\\[2pt]k_{4}^{v}&={\vec {a}}_{i}\!\left({\vec {x}}_{i}^{n}+\Delta t\,k_{3}^{x},\;{\vec {v}}_{i}^{n}+\Delta t\,k_{3}^{v}\right),\\k_{4}^{x}&={\vec {v}}_{i}^{n}+\Delta t\,k_{3}^{v},\\k_{4}^{L}&={\vec {\tau }}_{i}\!\left({\vec {x}}_{i}^{n}+\Delta t\,k_{3}^{x},\;{\vec {v}}_{i}^{n}+\Delta t\,k_{3}^{v}\right).\end{aligned}}}$

更新式

${\displaystyle {\begin{aligned}{\vec {v}}_{i}^{\,n+1}&={\vec {v}}_{i}^{\,n}+{\frac {\Delta t}{6}}\!\left(k_{1}^{v}+2k_{2}^{v}+2k_{3}^{v}+k_{4}^{v}\right),\\[4pt]{\vec {x}}_{i}^{\,n+1}&={\vec {x}}_{i}^{\,n}+{\frac {\Delta t}{6}}\!\left(k_{1}^{x}+2k_{2}^{x}+2k_{3}^{x}+k_{4}^{x}\right),\\[4pt]{\vec {L}}_{i}^{\,n+1}&={\vec {L}}_{i}^{\,n}+{\frac {\Delta t}{6}}\!\left(k_{1}^{L}+2k_{2}^{L}+2k_{3}^{L}+k_{4}^{L}\right).\end{aligned}}}$