DEMの弾塑性接触モデル

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

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

Summarize Timeline Top Qs Fact Check

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 接触理論による法線方向力 ${\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}}$ はそれぞれ粒子のヤング率およびポアソン比である^[8]。

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

Summarize Timeline Top Qs Fact Check

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

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 本の直線で近似することで、計算効率と実用性を両立させるモデルと位置付けられている^[10]。

Summarize Timeline Top Qs Fact Check

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

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

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

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

と定義している^[5]。ここで ${\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}}$ を最大重なりに依存させ、塑性変形の進行に伴う接触剛性の増大を表現する。
• 粘性減衰やローリング抵抗モデルと組み合わせることで、低速から高速まで広い変形モードを扱える^[13]。

Summarize Timeline Fact Check

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

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

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

詳細な式は複雑であるが、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 型の線形ヒステリシスモデルで近似される^[15]。

Summarize Timeline Top Qs Fact Check

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

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

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

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}$ のように定義される^[9]。これにより、塑性変形により接触面積が増加すると最大付着力も増大するという実験事実を、履歴依存の形で表現できる。

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