個別要素法 - Wikiwand

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

概略

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

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

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

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

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

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

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

特徴

計算コストは ${\mathcal {O}}(N)$ もしくは ${\mathcal {O}}(NlogN)$ で線形に近くスケールする。
明示的時間積分のため安定性条件は $\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

離散要素法モデル

個別要素法（離散要素法）では、粒子間や粒子–壁間に働く接触力をどのようにモデル化するかによって、再現できる材料挙動の範囲が大きく変化する。基本的な線形ばね–ダッシュポットモデルや Hertz–Mindlin 接触則に加えて、塑性変形、付着・凝集力、粒子間結合、粉末のバルク変形などを表現するために、さまざまな拡張接触モデルが提案されてきた。本節では、代表的な離散要素法モデルとして、DEMの弾塑性接触モデル、離散要素法の付着力モデル、ボンデッドパーティクルモデル（bonded particle model; BPM）、DEM-MDRモデルの概要を示す。

DEMの弾塑性接触モデル

DEMの弾塑性接触モデルは、法線方向の力–変位関係に塑性変形と履歴依存性を導入した接触則であり、高い圧縮や繰り返し荷重を受ける顆粒材料のシミュレーションに用いられる。弾性のみを仮定する線形ばね–ダッシュポットや Hertz–Mindlin モデルでは、接触点での永久変形や接触剛性の変化を表現できないため、これを拡張したDEMの弾塑性接触モデルが開発された。代表的なものとして、Walton–Braunモデル、Luding 接触モデル、Thornton–Ning モデル、Edinburgh Elasto-Plastic Adhesion (EEPA) モデルなどがあり、これらの系譜や詳細な数式は専用項目「DEMの弾塑性接触モデル」で整理されている。

Walton–Braun モデルは、荷重時と除荷時で異なる剛性を持つ線形ヒステリシスばねを仮定し、接触履歴に応じた残留オーバーラップとエネルギー散逸を表現する。副次的パラメータが少なく実装が容易なことから、多くの DEM コードにおける「ヒステリシスばね」モデルの原型となっている。Luding 接触モデルは Walton–Braun モデルを拡張し、除荷剛性を履歴最大重なりの関数とすることで塑性硬化を取り入れるとともに、法線力が引張側にまで達する粘着分枝を導入している。この結果、非凝集粒状体から微粉体まで、弾塑性履歴と付着力を同時に扱うことができる。

Thornton–Ning モデルは、弾完全塑性材料の接触力学に基づいて、接触半径・オーバーラップ・法線力の関係を理論的に導出した非線形弾塑性・粘着モデルである。弾性域では Hertz 解に一致し、降伏後は有限要素解析に基づく近似式を用いて弾塑性域と完全塑性域の応答を連続的に表現する。EEPA モデルは、Luding 型の線形ヒステリシス接触則に JKR 理論と表面エネルギーの概念を組み合わせたものであり、履歴依存の付着力と塑性変形を同一枠組みで扱える。粉体の圧縮・せん断・混合など、実験キャリブレーションが重要な工業応用で用いられる。

離散要素法の付着力モデル

付着力モデル（凝集力モデル）は、粒子間に働く静電力、ファンデルワールス力、液架橋による毛管力などの短距離引力を、離散要素法の接触則に組み込むためのモデル群を指す。非凝集性の乾燥砂などでは、接触力は主として弾性反発力と摩擦力で記述されるが、湿潤砂、粉体、微粒子などの凝集性材料では、付着成分を導入しないと安息角の増加、ブリッジング、アーチング、付着性堆積といった特徴的なマクロ挙動を再現できない。付着力の起源やモデルの体系的な整理については「離散要素法の付着力モデル」に詳しい。

離散要素法の付着力モデルは、多くの場合、法線方向接触力 $F_{n}$ を「弾性成分」と「付着成分」に分解し、Hertz 型あるいは線形弾性モデルに付着分枝を重ね合わせる形で定式化される。Hertz–JKR 型の理論モデルでは、表面エネルギーからプルオフ力と接触半径を評価し、弾性接触と付着を連成させる。一方、商用コードなどで広く利用される簡易 JKR（SJKR）モデルや線形付着モデルでは、最大引張力と作用距離をパラメータとする線形バネとして付着力を近似し、数値安定性と計算効率を優先する。液架橋接触モデルは、粒子間の液体メニスカスによる毛管力を直接扱うもので、液量に依存する破断距離や表面張力に比例する最大引張力を用いて、ペンデュラー域の湿潤粒状体を対象とする。

凝集性材料では、付着力は法線方向だけでなく、接線方向のせん断強度や転がり抵抗にも影響を及ぼす。そのため、多くのモデルではクーロン摩擦則のしきい値に付着成分を加える形で、見かけの内部摩擦角や付着強度を増加させる拡張が導入される。付着パラメータの同定には、安息角試験、ホッパー排出試験、引張試験などの室内実験と DEM シミュレーションを組み合わせたキャリブレーション手法が用いられている。

ボンデッドパーティクルモデル（BPM）

ボンデッドパーティクルモデル（bonded particle model, BPM）は、粒子間にばね–ダッシュポットからなる質量のない結合（ボンド）を導入することにより、岩石やコンクリートなどの脆性材料の破壊を表現する離散要素法の一種である。通常の DEM が粒子間接触面での局所力を記述するのに対し、BPM では粒子間に仮想的な接合断面を設定し、その断面を通じて法線力・せん断力に加えて、曲げモーメントやねじりモーメントを伝達できるようにする。ボンドの断面積や断面二次モーメント、長さに基づいて、法線・接線・曲げ・ねじりの各方向のばね定数とダンピング係数が定義される。詳細な定式化は「ボンデッドパーティクルモデル」にまとめられている。

BPM では、各ボンドに対して法線応力およびせん断応力の破壊基準を導入し、応力が許容値を超えたときにボンドを削除することで亀裂の発生と進展を記述する。多数のボンドが連鎖的に破壊することで、マクロには脆性破壊、準脆性破壊、延性破壊など多様な挙動を再現できる。ボンド剛性や強度を材料試験に基づきキャリブレーションすることで、岩盤の一軸圧縮試験、三軸試験、割裂引張試験などを模擬し、き裂パターンや強度特性を再現する応用が一般的である。

DEM-MDRモデル

DEM-MDRモデルは、接着弾性–完全塑性粒子の接触挙動と粉末の大変形圧縮を扱うために提案された機械的導出接触モデルであり、次元削減法（method of dimensionality reduction, MDR）に基づく局所接触応答と、粒子内部のバルク弾性応答を組み合わせる点に特徴がある。MDR により三次元接触問題を一次元のばね–ダッシュポット列に写像し、その応答を粒子半径の変化（見掛け半径）の更新則と結合させることで、広いひずみ範囲にわたる粉末の力学挙動を記述する。モデルは、ヤング率、ポアソン比、降伏応力、有効表面エネルギー、バルク応答の開始を制御する無次元パラメータ、反発係数といった少数の物理パラメータで特性づけられる。

DEM-MDRモデルは、単一粒子上の複数接触からの寄与と粒子全体の体積ひずみを整合的に扱えるため、高い圧縮比を伴う粉末成形やタブレット圧縮などで、実験に近い圧縮曲線やスプリングバック挙動を再現できると報告されている。数値実装としては、分子動力学コード LAMMPS の granular パッケージに基づく公開実装が知られており、大規模粉末系のシミュレーションに利用されている。

アルゴリズム

1. 初期化

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

2. 時間積分ループ

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

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

(b) 接触履歴: クーロン摩擦や粘弾性モデルでは接触履歴が必要となる。前ステップで既接触なら; $\delta _{t}^{\,n}=\delta _{t}^{\,n-1}+v_{t}^{\,n}\,\Delta t$; 未接触なら; $\delta _{t}^{\,n}=v_{t}^{\,n}\,\Delta t$

(c) 接触力計算: Hertz-Mindlin 式などを用い、法線力 $\mathbf {F} _{n}$ と接線力 $\mathbf {F} _{t}$ を評価する。

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

(e) 運動方程式の陽積分: 例として単純並進のみを示す。; $\mathbf {v} ^{\,n+1}=\mathbf {v} ^{\,n}+\mathbf {a} ^{\,n}\,\Delta t$; $\mathbf {x} ^{\,n+1}=\mathbf {x} ^{\,n}+\mathbf {v} ^{\,n+1}\,\Delta t$; 角運動は陽オイラーや四元数を用いるのが一般的である。

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

(g) 時間終了判定: $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) シミュレーションで、近傍粒子対を線形時間で検索する空間分割アルゴリズムである。計算領域をカットオフ半径以上の立方体セルに区切り、各セル内に存在する粒子番号をリンクリストとして保持することで、粒子 $N$ 個に対するペア探索の計算量を ${\mathcal {O}}(N)$ に削減できる。

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

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

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

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

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

物理量計算

球粒子の慣性モーメント

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

Hertz–Mindlin 接触

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

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

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

粘性減衰と反発係数

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

数値積分の安定条件

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

法線・接線成分の計算

単位法線ベクトル

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

相対並進速度

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

法線速度

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

接線速度

\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}

法線方向オーバーラップ

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

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

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

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

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

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

ここで

\mathbf {t} _{ij}^{(n+1)}

は接線方向の単位ベクトルであり、

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

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

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

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

角速度の分解（粒子 i）

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

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

接触モデル

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

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

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

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

摩擦モデル

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

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

転がり抵抗

$\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}|}}$ ^[10]

運動方程式

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

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

積分

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

オイラー前進法（Explicit Euler）

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

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

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

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

速度半ステップ $\displaystyle {\vec {v}}_{i}^{\,n+{\tfrac {1}{2}}}={\vec {v}}_{i}^{\,n-{\tfrac {1}{2}}}+\Delta t\,{\vec {a}}_{i}^{\,n}$
位置更新 $\displaystyle {\vec {x}}_{i}^{\,n+1}={\vec {x}}_{i}^{\,n}+\Delta t\,{\vec {v}}_{i}^{\,n+{\tfrac {1}{2}}}$
角運動量半ステップ $\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）: 内部ステージ

${\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}}$

更新式

${\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}}$