DEMの弾塑性接触モデル

DEMの弾塑性接触モデル（DEMのだんそせいせっしょくモデル、elasto-plastic contact models in DEM）は、離散要素法（Discrete Element Method; DEM）において、粒子間あるいは粒子–壁間の接触を弾性変形および塑性変形、場合によっては付着力（凝集力）を含めて表現する法線方向接触モデルの総称である。

特に線形ヒステリシス型の弾塑性モデルとして、Walton–Braun モデル（しばしば Walton–Brown モデル とも表記される）、Luding 接触モデル、Thornton–Ning モデル、Edinburgh Elasto-Plastic Adhesion (EEPA) モデルが知られており、Walton–Braun 型モデルを起点として Luding モデル、さらに Thornton–Ning モデルや EEPA モデルへと発展してきた系譜をなす^[1][2]。

背景

DEM では、各粒子の並進運動と回転運動は、粒子間接触力と外力（重力など）から決まる運動方程式によって記述される。接触力モデルは、接触面での法線方向重なり量と接触履歴に基づき、法線方向力およびせん断方向力を与える局所構成則である。従来広く用いられてきた線形バネ–ダッシュポットモデルや Hertz–Mindlin モデルは、弾性的な接触変形と粘性減衰を表現するが、強い圧縮や繰り返し荷重を受ける顆粒材料では、接触点での塑性変形と、それに起因する履歴依存性や付着力の増大を必ずしも十分には表現できない^[3]。

この問題に対処するため、Walton と Braun による部分ラッチばね（partially latching spring）型の線形ヒステリシスモデル^[4]に端を発し、Luding による粘着・弾塑性接触モデル^[5]、Hertz 接触理論と塑性理論に基づく Thornton–Ning モデル^[6]、およびこれらを工学的に拡張した EEPA モデル^[7]など、さまざまな弾塑性接触モデルが提案されている。これらのモデルは、法線方向の力–変位関係に塑性変形と付着力を組み込み、顆粒体の圧縮・せん断・流動といった巨視的挙動をより現実的に再現することを目的としている^[8]。

接触力学の基本概念

法線重なりと有効半径

DEM の球形粒子 i, j が接触しているとき、法線方向重なり（オーバーラップ）${\displaystyle \delta }$ は ${\displaystyle \delta =(R_{i}+R_{j})-\left|{\vec {r}}_{i}-{\vec {r}}_{j}\right|}$ で定義される。ここで ${\displaystyle R_{i},R_{j}}$ は粒子半径、${\displaystyle {\vec {r}}_{i},{\vec {r}}_{j}}$ は粒子中心位置ベクトルである。接触問題に現れる有効半径 ${\displaystyle R^{*}}$ は ${\displaystyle R^{*}={\frac {R_{i}R_{j}}{R_{i}+R_{j}}}}$ として定義される。

Hertz 接触と JKR 接触

線形弾性球同士の無粘着接触に対して、Hertz 接触理論による法線方向力 ${\displaystyle F_{\rm {H}}}$ は ${\displaystyle F_{\rm {H}}={\frac {4}{3}}E^{*}{\sqrt {R^{*}}}\,\delta ^{3/2}}$ で与えられる。ここで有効ヤング率 ${\displaystyle E^{*}}$ は ${\displaystyle {\frac {1}{E^{*}}}={\frac {1-\nu _{i}^{2}}{E_{i}}}+{\frac {1-\nu _{j}^{2}}{E_{j}}}}$ であり、${\displaystyle E_{i},E_{j}}$ と ${\displaystyle \nu _{i},\nu _{j}}$ はそれぞれ粒子のヤング率およびポアソン比である^[9]。

粒子間に表面エネルギー ${\displaystyle \Gamma }$ に起因する付着力が存在する場合、Johnson–Kendall–Roberts (JKR) 理論によれば、弾性接触と付着力が結合した状態を記述できる。JKR 理論における引き離し時（破断時）の最大引張力（プルオフ力）${\displaystyle F_{\rm {JKR}}}$ は ${\displaystyle F_{\rm {JKR}}={\tfrac {3}{2}}\pi R^{*}\Gamma }$ で与えられる^[10]。Thornton–Ning 型モデルおよび EEPA モデルでは、この JKR プルオフ力や接触パッチ半径を基準として、履歴依存の付着力が導入される^[11]。

Walton–Braun（Walton–Brown）モデル

概要

Walton–Braun モデルは、Walton と Braun によって提案された線形ヒステリシス型の弾塑性接触モデルであり、「部分ラッチばね（partially latching spring）モデル」とも呼ばれる^[12]。本モデルは、荷重時と除荷時で異なるばね剛性を用いることで、弾塑性接触における残留変形とエネルギー散逸を表現する。一般には凝集力を含まない非粘着モデルとして用いられ、Altair EDEM における「Hysteretic Spring」接触モデルの基礎となっている^[13]。

力–変位関係

Walton–Braun 型モデルでは、法線方向の重なりを ${\displaystyle \delta }$、これまでの履歴における最大重なりを ${\displaystyle \delta _{\max }}$ とし、荷重時の剛性 ${\displaystyle k_{L}}$（loading）と除荷・再荷重時の剛性 ${\displaystyle k_{U}}$（unloading/reloading, ${\displaystyle k_{U}\geq k_{L}}$）を導入する。塑性変形に対応する残留重なり ${\displaystyle \delta _{0}}$ は ${\displaystyle \delta _{0}=\left(1-{\frac {k_{L}}{k_{U}}}\right)\,\delta _{\max }}$ で与えられる。

接触法線力 ${\displaystyle F_{n}}$ は、次の区分線形関数として表される。 ${\displaystyle F_{n}(\delta ,\delta _{\max })={\begin{cases}k_{L}\,\delta ,&\delta \geq \delta _{\max }\quad {\text{(initial load)}}\\[4pt]k_{U}(\delta -\delta _{0}),&\delta _{0}\leq \delta \leq \delta _{\max }\quad {\text{(Unload and reload)}}\\[4pt]0,&\delta \leq \delta _{0}\quad {\text{(Disappearance of contact)}}\end{cases}}}$

荷重過程では、力は原点を通る直線 ${\displaystyle F_{n}=k_{L}\delta }$ に沿って増加する。最大重なり ${\displaystyle \delta _{\max }}$ に達した後の除荷では、より高い剛性 ${\displaystyle k_{U}}$ に沿って力が減少し、${\displaystyle \delta =\delta _{0}}$ で力がゼロとなる。このとき残留重なり ${\displaystyle \delta _{0}}$ が正となり塑性変形を表す。再荷重は同じ ${\displaystyle k_{U}}$ の直線に沿って進み、${\displaystyle \delta =\delta _{\max }}$ に達すると再び ${\displaystyle k_{L}}$ の荷重枝に戻る。Walton–Braun モデルは、非線形接触力学を単純な 2 本の直線で近似することで、計算効率と実用性を両立させるモデルと位置付けられている^[14]。

Luding 接触モデル

Walton–Braun モデルからの拡張

Luding による粘着・弾塑性ヒステリシスモデルは、Walton–Braun 型モデルに基づき、付着力（負の法線力）と塑性変形に伴う接触面積の増加を考慮するために開発された^[15]。このモデルは、弾塑性履歴と粘着力を同時に扱えるため、微粉体や凝集性顆粒材料のせん断流動・圧縮挙動のシミュレーションに広く用いられている^[16]。

Walton–Braun モデルが非粘着であり、法線力がゼロに達すると接触が消失するのに対し、Luding モデルでは法線力が負になる領域（付着力領域）が導入されている。また、除荷剛性 ${\displaystyle k_{2}}$ を最大重なり ${\displaystyle \delta _{\max }}$ の関数とすることで、塑性変形の進行に伴う接触剛性の増加（ひずみ硬化）を表現している。

線形ヒステリシス接触則

Luding の線形粘着・弾塑性接触則では、法線方向のヒステリシス力 ${\displaystyle f_{\rm {hys}}}$ は ${\displaystyle f_{\rm {hys}}(\delta$ ;\delta _{\max })={\begin{cases}k_{1}\,\delta ,&k_{2}(\delta -\delta _{0})\geq k_{1}\delta \\[4pt]k_{2}(\delta -\delta _{0}),&k_{1}\delta \geq k_{2}(\delta -\delta _{0})\geq -k_{c}\,\delta \\[4pt]-k_{c}\,\delta ,&-k_{c}\,\delta \geq k_{2}(\delta -\delta _{0})\end{cases}}}

と定義される^[17]。ここで

• ${\displaystyle k_{1}}$：荷重枝の剛性（Walton–Braun の ${\displaystyle k_{L}}$ に相当）
• ${\displaystyle k_{2}}$：除荷・再荷重枝の剛性（履歴変数 ${\displaystyle \delta _{\max }}$ の関数）
• ${\displaystyle k_{c}}$：付着枝の剛性（負の勾配）
• ${\displaystyle \delta _{\max }}$：これまでの最大重なり（履歴変数）
• ${\displaystyle \delta _{0}}$：塑性残留重なり

である。残留重なり ${\displaystyle \delta _{0}}$ は ${\displaystyle \delta _{0}=\left(1-{\frac {k_{1}}{k_{2}}}\right)\,\delta _{\max }}$ で与えられ、この値より小さな重なり側では法線力が負となり付着力を表現する。付着枝への遷移点に対応する重なり ${\displaystyle \delta _{\min }}$ は ${\displaystyle \delta _{\min }={\frac {k_{2}-k_{1}}{k_{2}+k_{c}}}\,\delta _{\max }}$ であり、対応する最小力（最大付着力）${\displaystyle f_{\min }}$ は ${\displaystyle k_{1},k_{2},k_{c},\delta _{\max }}$ の組合せから決まる。

変動除荷剛性と塑性流れ限界

Luding は、除荷剛性 ${\displaystyle k_{2}}$ を最大重なり ${\displaystyle \delta _{\max }}$ の関数として ${\displaystyle k_{2}(\delta _{\max })={\begin{cases}{\hat {k}}_{2},&\delta _{\max }\geq \delta _{\max }^{*}\\[4pt]k_{1}+({\hat {k}}_{2}-k_{1})\,{\dfrac {\delta _{\max }}{\delta _{\max }^{*}}},&0\leq \delta _{\max }\leq \delta _{\max }^{*}\end{cases}}}$

と定義している^[18]。ここで ${\displaystyle {\hat {k}}_{2}}$ は到達し得る最大除荷剛性、${\displaystyle \delta _{\max }^{*}}$ は塑性流れ限界重なりであり、粒子材料の降伏応力や代表粒径などから経験的に決められる。塑性変形が進むにつれて ${\displaystyle \delta _{\max }}$ が増大し、除荷剛性 ${\displaystyle k_{2}}$ が ${\displaystyle k_{1}}$ から ${\displaystyle {\hat {k}}_{2}}$ へと連続的に遷移することで、ひずみ硬化挙動を近似的に表現できる。

このように、Luding モデルは以下の点で Walton–Braun モデルを拡張している。

• 非粘着（法線力が非負）から粘着（法線力が負）までを含む 3 本の線形枝（荷重、除荷、付着）を導入する。
• 除荷剛性 ${\displaystyle k_{2}}$ を最大重なりに依存させ、塑性変形の進行に伴う接触剛性の増大を表現する。
• 粘性減衰やローリング抵抗モデルと組み合わせることで、低速から高速まで広い変形モードを扱える^[19]。

Thornton–Ning モデル

接触力学に基づく非線形モデル

Thornton と Ning によるモデルは、接触力学と弾完全塑性材料の理論に基づき、接触面での応力状態を明示的に考慮した非線形弾塑性・粘着接触モデルである^[20]。このモデルでは、Hertz 接触理論と JKR 理論を組み合わせ、降伏応力を超えた後の接触圧の分布を「切り詰められた Hertz 圧力分布」として近似することで、弾性域・弾塑性域・完全塑性域を連続的に扱う^[21]。

Thornton–Ning モデルの特徴は次の通りである。

• 接触半径 ${\displaystyle a}$、重なり ${\displaystyle \delta }$、法線力 ${\displaystyle F_{n}}$ の関係を、材料の有効ヤング率 ${\displaystyle E^{*}}$、降伏応力 ${\displaystyle \sigma _{y}}$、表面エネルギー ${\displaystyle \Gamma }$ に基づき理論的に導出する。
• 降伏後の弾塑性域では、荷重–変位曲線が Hertz 解よりも小さな勾配を示し、塑性ひずみの蓄積に伴うエネルギー散逸を表現する。
• 粘着を伴う場合には JKR 理論に基づく付着エネルギーを組み込み、ある臨界衝突速度以下では粒子が「スティック（再付着）」し、それ以上では「バウンス（反発）」する遷移を記述する^[22]。

法線力の概略的表現

詳細な式は複雑であるが、Thornton–Ning モデルの法線力は概念的には次のように表現できる。

• 弾性域（初期接触）：

${\displaystyle F_{n}\approx F_{\rm {H}}(\delta )={\frac {4}{3}}E^{*}{\sqrt {R^{*}}}\,\delta ^{3/2}}$

• 弾塑性域：

${\displaystyle F_{n}=f_{\rm {ep}}(\delta ;E^{*},R^{*},\sigma _{y})}$

ここで ${\displaystyle f_{\rm {ep}}}$ は有限要素解析などに基づき同定された非線形関数であり、Hertz 解よりも小さな勾配を持つ。

• 粘着を考慮した場合（JKR 型）：

${\displaystyle F_{n}=f_{\rm {ep}}(\delta )-F_{\rm {adh}}(a,\Gamma )}$

と書け、${\displaystyle F_{\rm {adh}}}$ は接触半径 ${\displaystyle a}$ と表面エネルギー ${\displaystyle \Gamma }$ の関数である。

このモデルは、非線形接触力学に忠実である一方、数値実装が複雑で計算負荷も高いため、DEM ではしばしば Luding 型や Pasha 型の線形ヒステリシスモデルで近似される^[23]。

Edinburgh Elasto-Plastic Adhesion (EEPA) モデル

Luding モデルからの発展

Edinburgh Elasto-Plastic Adhesion (EEPA) モデルは、Morrissey らにより提案され、Coetzee によって PFC や EDEM に実装された歴史依存型の粘着・弾塑性接触モデルである^[24]。EEPA モデルは、Luding の線形ヒステリシスモデルを基盤としつつ、

• 荷重枝および除荷・再荷重枝に非線形指数 ${\displaystyle m}$ を導入して非線形硬化・軟化挙動を表現する。
• 付着枝に対して独立した変形指数（しばしば ${\displaystyle \chi }$ と記される）を導入し、顕著な付着強度の減衰（ソフトニング）を表現する。
• 最小力（最大付着力）を、表面エネルギー ${\displaystyle \Gamma }$ に基づく JKR プルオフ力に整合するように定義する。

といった点に特徴がある^[25]。

法線方向の非線形ヒステリシス則

EEPA モデルの法線方向接触力 ${\displaystyle F_{n}}$ は、Luding モデルと同様に 3 本の枝（荷重枝、除荷・再荷重枝、付着枝）から構成されるが、それぞれに指数を導入した形で表される。概念的には ${\displaystyle F_{n}(\delta )={\begin{cases}F_{0}+k_{1}\,\delta ^{m},&0\leq \delta \leq \delta _{\max }\\[4pt]F_{0}+k_{2}\,(\delta -\delta _{p})^{m},&\delta _{p}\leq \delta \leq \delta _{\max }\\[4pt]F_{\min }+k_{a}\,\delta ^{\chi },&\delta \leq \delta _{p}\end{cases}}}$

と書くことができる。ここで

• ${\displaystyle F_{0}}$：接触開始時のプルオフ力（通常は負値）
• ${\displaystyle \delta _{p}}$：塑性重なり（除荷枝上で ${\displaystyle F_{n}=F_{0}}$ となる点）
• ${\displaystyle F_{\min }}$：履歴に依存する最小力（最大付着力）
• ${\displaystyle k_{1},k_{2},k_{a}}$：それぞれ荷重枝・除荷枝・付着枝の剛性パラメータ
• ${\displaystyle m,\chi }$：荷重・除荷枝および付着枝の変形指数

である。線形版 EEPA では ${\displaystyle m=\chi =1}$ とすることで、Luding モデルと同様の線形ヒステリシス形に退化する。

最小力 ${\displaystyle F_{\min }}$ は、指定した表面エネルギー ${\displaystyle \Gamma }$ と接触パッチ半径 ${\displaystyle a}$ から、JKR プルオフ力と同程度になるように ${\displaystyle F_{\min }\approx -{\tfrac {3}{2}}\pi \Gamma a}$ のように定義される^[26]。これにより、塑性変形により接触面積が増加すると最大付着力も増大するという実験事実を、履歴依存の形で表現できる。

さらに EEPA モデルは、Hertz 接触理論に基づくせん断剛性 ${\displaystyle k_{t}\propto G^{*}{\sqrt {R^{*}\,\delta _{n}}}}$ を用いたせん断接触則と、いわゆる Type C ローリング抵抗モデルを組み合わせており、粘着性粒子の回転・滑り挙動を一貫して扱うことができる^[27]。

モデル間の比較と系譜

系譜的な関係

DEM における代表的な弾塑性接触モデルの系譜は、おおまかに次のように整理できる。

1. Walton–Braun モデル：弾塑性だが非粘着。荷重・除荷で異なる剛性を用いる線形ヒステリシスばねモデルであり、多くの DEM コードに実装されている Hysteretic Spring モデルの原型となっている^[28][29]。
2. Luding モデル：Walton–Braun モデルを拡張し、付着枝と履歴依存の除荷剛性を導入した線形粘着・弾塑性モデル。凝集性粉体のバルク挙動を記述するための標準的な接触モデルの一つである^[30]。
3. Thornton–Ning モデル：Hertz–JKR 接触理論と弾完全塑性材料理論に基づく非線形理論モデルであり、線形ヒステリシスモデルの物理的整合性を評価する基準として用いられる^[31][32]。
4. EEPA モデル：Luding モデルをベースに、非線形指数と JKR に基づく表面エネルギーを導入した工学的モデルであり、Walton–Braun → Luding → EEPA という線形ヒステリシス系譜の延長線上に位置づけられる^[33]。

この意味で、Luding モデルからは、理論的な接触力学を追求する方向（Thornton–Ning 型）と、工学的パラメータで扱いやすい非線形ヒステリシスモデル（EEPA 型）という 2 つの発展方向が存在すると解釈できる。

数式レベルで見た主な違い

代表的モデルの違いを、法線方向力 ${\displaystyle F_{n}}$ の関数形に注目して整理すると次のようになる。

• 線形か非線形か
  • Walton–Braun, Luding：${\displaystyle F_{n}}$ は重なり ${\displaystyle \delta }$ に対して区分線形（1 次）であり、分岐点のみが履歴依存である。
  • Thornton–Ning：${\displaystyle F_{n}}$ は Hertz 型 ${\displaystyle \delta ^{3/2}}$ と弾塑性接触理論に基づく非線形関数からなり、全域で非線形である。
  • EEPA：荷重・除荷枝に指数 ${\displaystyle m}$、付着枝に指数 ${\displaystyle \chi }$ を導入した区分非線形関数となる。

• 粘着の扱い
  • Walton–Braun：非粘着モデルであり、法線力はゼロ以上のみをとる。
  • Luding：線形付着枝（負の勾配）を持ち、最小力（最大付着力）は ${\displaystyle k_{1},k_{2},k_{c},\delta _{\max }}$ の組合せから決まる。
  • Thornton–Ning：JKR 理論に基づく付着エネルギーを明示的に導入し、接触半径と表面エネルギー ${\displaystyle \Gamma }$ からプルオフ力を計算する。
  • EEPA：JKR プルオフ力と整合するように ${\displaystyle F_{\min }}$ を定義し、塑性変形に伴う接触面積の変化を履歴依存パラメータ ${\displaystyle \delta _{\max },\delta _{p}}$ を通じて反映する。

• パラメータの物理的意味
  • Walton–Braun：${\displaystyle k_{L},k_{U}}$ は経験的なばね剛性であり、材料物性との直接対応は限定的である。
  • Luding：${\displaystyle k_{1},k_{2},k_{c}}$ など、依然として経験的パラメータだが、塑性変形範囲や最大付着力などマクロな挙動と定性的な対応がある^[34]。
  • Thornton–Ning：主要パラメータは ${\displaystyle E^{*},\sigma _{y},\Gamma }$ など物性値であり、理論的に一意に決まる。
  • EEPA：表面エネルギー ${\displaystyle \Gamma }$ によって最大付着力を、荷重枝と除荷枝の勾配比によって「接触塑性度」を制御するなど、Luding 型パラメータにより強い物理的意味を持たせている^[35]。

応用

Luding モデルや EEPA モデルは、粉体ハンドリング、土壌力学、鉱石や鉱滓の流動解析、農業機械における土壌・作物シミュレーションなど、顆粒材料のバルク挙動に強い履歴依存性や粘着性が現れる分野で広く用いられている^[36][37]。Thornton–Ning モデルは、単一粒子衝突のスティック/バウンス境界の予測や、線形ヒステリシスモデルの妥当性評価に利用されることが多い^[38]。

DEM パラメータの同定に関するレビューでは、接触モデルの選択が巨視的な流動特性（安息角、圧縮曲線、トライボロジー特性など）に大きく影響することが指摘されており、特に凝集性材料のモデリングでは、Walton–Braun のような単純な弾塑性モデルよりも、Luding や EEPA のような履歴依存粘着モデルが推奨される^[39][40]。