DEM-MDRモデル

DEM-MDRモデル（英: DEM-MDR contact model）は、離散要素法（DEM）において接着弾性–完全塑性粒子の接触挙動と大変形粉末圧縮を記述するために提案された機械的導出接触モデルである。

モデルは次元削減法（method of dimensionality reduction, MDR）に基づく接触寄与と、粒子全体のバルク弾性応答を組み合わせることで、広いひずみ範囲における粉末の力学応答を表現する。主な定式化は Zunker と Kamrin による二部構成の論文および、その後の粉末圧縮実験との比較研究により提示されている^[1]^[2]^[3]。

本項目では、これらの研究に基づき構成された DEM-MDR 接触モデルの概要、数値実装、記号の定義を概説する。ソースコード実装は、分子動力学コード LAMMPS の公開実装を参考としている^[4]。

概要

DEM-MDR モデルは、単一粒子上の複数接触に対して、MDR に基づく局所接触応答と、粒子内部のバルク弾性応答を重ね合わせることにより、粒子半径の変化（見掛け半径）と接触力を同時に決定する接触則である。モデルは、以下の 6 つの測定可能で基本的な材料・接触パラメータにより特徴づけられる。

ヤング率 $E$
ポアソン比 $\nu$
降伏応力 $Y$
有効表面エネルギー（または付着エネルギー） $\gamma$
バルク応答の発現を制御する臨界拘束比 $\psi _{b}$
正規反発係数（coefficient of restitution） $\mathrm {CoR}$

モデルは、大変形粉末圧縮シミュレーションにおいて、塑性変形による粒子の膨張、接着破壊、非物理的接触のスクリーニングなどを一貫して取り扱うことを目的としている^[1]^[2]^[3]。

LAMMPS による実装例

DEM-MDR モデルは、LAMMPS においてユーザー定義接触モデルとして実装されることが想定されている。以下は、典型的な LAMMPS 入力スクリプトにおけるパラメータ設定例である。

# mdr = E, nu, Y, gamma, psi_b, CoR
variable YoungsModulus      equal 5e6
variable YieldStress        equal 1.9e5
variable PoissonsRatio      equal 0.4
variable SurfaceEnergy      equal 2
variable SurfaceEnergyWall  equal 0.0
variable CoR                equal 0.5
variable psi_b              equal 0.5

ここで、YoungsModulus はヤング率 $E$ 、PoissonsRatio はポアソン比 $\nu$ 、YieldStress は降伏応力 $Y$ 、SurfaceEnergy は粒子間の有効表面エネルギー $\gamma$ 、psi_b は臨界拘束比 $\psi _{b}$ を表す。

MDR 接触モデルの定式化

力の分解

DEM-MDR モデルでは、単一接触における法線方向の接触力 $F$ は、MDR に基づく局所接触力とバルク弾性応答による寄与の和として表される。

$F=F_{\mathrm {MDR} }+F_{\mathrm {Bulk} }.$

ここで、 $F_{\mathrm {MDR} }$ は MDR による局所接触寄与、 $F_{\mathrm {Bulk} }$ は粒子全体のバルク弾性応答による補正である。

運動学と変位の分割

粒子中心間オーバーラップに基づく見掛けの変位 $\delta _{0}$ は、MDR 部分とバルク部分に分割される。

$\delta _{0}=\delta _{0,\mathrm {MDR} }+\delta _{0,\mathrm {Bulk} }.$

ここで、 $R_{0}$ を初期粒子半径、 $R(\geq R_{0})$ を塑性変形後の見掛け半径、 $\delta$ を真のオーバーラップとすると、

$\delta =\delta _{0}+R-R_{0}$

で与えられる。ステップごとの変位増分 $\mathrm {d} \delta _{0}$ は、非接触自由面積と全表面積の比 $A_{\mathrm {free} }/A_{\mathrm {tot} }$ に応じて分割される。

$\mathrm {d} \delta _{0,\mathrm {MDR} }={\begin{cases}\mathrm {d} \delta _{0},&A_{\mathrm {free} }/A_{\mathrm {tot} }\geq \psi _{b},\\\mathrm {d} \delta _{0}-\mathrm {d} {\bar {\delta }}_{0},&A_{\mathrm {free} }/A_{\mathrm {tot} }<\psi _{b},\end{cases}}$

$\mathrm {d} \delta _{0,\mathrm {Bulk} }={\begin{cases}0,&A_{\mathrm {free} }/A_{\mathrm {tot} }\geq \psi _{b},\\\mathrm {d} {\bar {\delta }}_{0},&A_{\mathrm {free} }/A_{\mathrm {tot} }<\psi _{b},\end{cases}}$

ここで $\mathrm {d} {\bar {\delta }}_{0}$ は接触面積で加重された平均変位増分である。

単一粒子上に $N$ 個の接触が存在する場合、接触面積の総和 $A_{\mathrm {con} }$ と全表面積 $A_{\mathrm {tot} }$ はそれぞれ

$A_{\mathrm {con} }=\sum _{i=1}^{N}A_{C,i},$

$A_{C,i}=\pi \left(2\delta _{i}R-\delta _{i}^{2}\right)+c_{A,i},$

$A_{\mathrm {tot} }=4\pi R^{2}-2\pi \sum _{i=1}^{N}(\delta _{\max ,\mathrm {MDR} ,i}+\delta _{0,\mathrm {Bulk} ,i})\,R+A_{\mathrm {con} },$

と表される。非接触自由面積は

$A_{\mathrm {free} }=A_{\mathrm {tot} }-A_{\mathrm {con} }$

であり、平均変位増分は

$\mathrm {d} {\bar {\delta }}_{0}={\frac {1}{A_{\mathrm {con} }}}\sum _{i=1}^{N}\left(A_{C,i}\,{\boldsymbol {n}}_{i}\cdot \mathrm {d} {\boldsymbol {\delta }}_{0,i}\right)$

で与えられる。ここで ${\boldsymbol {n}}_{i}$ は接触法線ベクトル、 $\mathrm {d} {\boldsymbol {\delta }}_{0,i}$ は接触ごとの変位増分ベクトルである。

MDR 力の寄与

MDR 部分の接触力 $F_{\mathrm {MDR} }$ は、MDR 見掛け重なり $\delta _{\mathrm {MDR} }$ 、その最大経験値 $\delta _{\mathrm {MDR,max} }$ 、見掛け半径 $R$ の関数として

$F_{\mathrm {MDR} }(\delta _{\mathrm {MDR} },\delta _{\mathrm {MDR,max} },R)={\frac {E_{c}^{*}\,A\,B}{4}}\left[\arccos \!\left(1-{\frac {2\delta _{\mathrm {1D} }^{e}}{A}}\right)-\left(1-{\frac {2\delta _{\mathrm {1D} }^{e}}{A}}\right){\sqrt {{\frac {4\delta _{\mathrm {1D} }^{e}}{A}}-{\frac {4(\delta _{\mathrm {1D} }^{e})^{2}}{A^{2}}}}}\right]$

と書ける。ここで、 $\delta _{\mathrm {1D} }^{e}$ は 1 次元に写像された弾性変位であり、降伏前の弾性領域では

$\delta _{\mathrm {1D} }^{e}=\delta _{\mathrm {MDR} }$

とおける。複合平面ひずみ弾性率 $E_{c}^{*}$ は剛体平板と変形可能粒子の接触を仮定すると

$E_{c}^{*}={\frac {E}{1-\nu ^{2}}}$

に簡約される。

MDR 見掛け重なりとその最大経験値は

$\delta _{\mathrm {MDR} }=\delta _{0,\mathrm {MDR} }+R-R_{0},\qquad \delta _{\mathrm {MDR,max} }=\delta _{0,\mathrm {MDR,max} }+R-R_{0}$

で与えられる。

接触体は楕円体に写像され、その形状パラメータ $A,B$ は弾性領域と完全塑性領域で異なる。

弾性領域：

$A=4R,\qquad B=2R$

完全塑性領域：

$A={\frac {4p_{Y}}{E_{c}^{*}}}\,a_{\max },\qquad B=2a_{\max }$

完全塑性領域における平均接触圧 $p_{Y}$ は

$p_{Y}=Y\left(1.75\exp(-4.4\,\delta _{\mathrm {MDR,max} }/R)+1\right)$

で与えられる。一方、Hertz 接触に基づく平均接触圧 ${\bar {p}}_{H}$ を

${\bar {p}}_{H}={\frac {4E_{c}^{*}}{3\pi {\sqrt {R}}}}{\sqrt {\delta }}$

と定義し、 ${\bar {p}}_{H}$ が $p_{Y}$ と交差する点で弾性から塑性へ遷移すると仮定する。

塑性領域では、変換された 1 次元弾性変位は

$\delta _{\mathrm {1D} }^{e}={\frac {\delta _{\mathrm {MDR} }-\delta _{\mathrm {MDR,max} }+\delta _{\mathrm {1D,max} }^{e}+\delta _{R}}{1+\delta _{R}/\delta _{\mathrm {1D,max} }^{e}}},$

$\delta _{\mathrm {1D,max} }^{e}={\frac {A}{2}}$

で定義される。ここで $\delta _{R}$ は大変形の影響を考慮した変位補正であり、

$\delta _{R}={\frac {F(\delta _{\mathrm {MDR,max} },\delta _{\mathrm {MDR,max} },R)}{\pi a_{\max }^{2}}}\,{\frac {2a_{\max }^{2}(\nu -1)-z_{R}(2\nu -1)\left({\sqrt {a_{\max }^{2}+z_{R}^{2}}}-z_{R}\right)}{2G{\sqrt {a_{\max }^{2}+z_{R}^{2}}}}},$

$z_{R}=R-(\delta _{\mathrm {MDR,max} }-\delta _{\mathrm {1D,max} }^{e})$

と書かれる。

最大経験接触半径 $a_{\max }$ は

$a_{\max }={\sqrt {\,2\delta _{\max }R-\delta _{\max }^{2}+{\frac {c_{A}}{\pi }}\,}}$

により定義される。 $c_{A}$ は弾性–塑性遷移における接触半径の連続性により

$c_{A}=\pi (\delta _{Y}^{2}-\delta _{Y}R)$

として決まる。降伏変位 $\delta _{Y}$ は

${\frac {4E_{c}^{*}}{3\pi {\sqrt {R}}}}{\sqrt {\delta _{Y}}}=Y\left(1.75\exp(-1.1\,\delta _{Y}/R)+1\right)$

を満たす値として求められる。

見掛け半径の更新

塑性変形の非圧縮性を尊重するため、バルク弾性応答が活性化していない場合（ $A_{\mathrm {free} }/A_{\mathrm {tot} }\geq \psi _{b}$ ）には粒子の見掛け半径 $R$ を増大させる。増分形式では

$\Delta R={\begin{cases}\displaystyle \max \left[{\frac {\Delta V^{\,e}-\sum _{i=1}^{N}\pi \Delta \delta _{0,\mathrm {MDR} ,i}(2\delta _{0,\mathrm {MDR} ,i}R_{0}-\delta _{0,\mathrm {MDR} ,i}^{2}+R^{2}-R_{0}^{2})}{2\pi R\sum _{i=1}^{N}(\delta _{0,\mathrm {MDR} ,i}+R-R_{0})-4\pi R^{2}}},\,0\right],&A_{\rm {free}}/A_{\rm {tot}}\geq \psi _{b},\\[1ex]0,&A_{\rm {free}}/A_{\rm {tot}}<\psi _{b},\end{cases}}$

と表される。ここで弾性体積変化 $\Delta V^{e}$ はステップ間の体積差から計算され、体積 $V$ は

$V=V_{0}\left[1+{\frac {1}{3\kappa V_{0}}}\,\mathrm {tr} \!\left(\sum _{c\in V}{\boldsymbol {f}}^{c}\otimes {\boldsymbol {b}}^{c}\right)\right]$

として評価される。 $V_{0}$ は初期体積、 $\kappa$ はバルク弾性率、 ${\boldsymbol {f}}^{c}$ は接触力ベクトル、 ${\boldsymbol {b}}^{c}$ はブランチベクトルである。

MDR 接着モデル

DEM-MDR モデルは、JKR 型破壊力学に基づく付着を MDR フレームワークへ組み込むことで、接触のアンロード時に見られる付着破壊を表現する。基本的な考え方は、半空間を離散化する 1 次元スプリングがインデンター表面に「張り付く」ことを許容し、アンロード時には張力ばねの伸びに応じて付着力を評価することである。

有効表面エネルギー $\Delta \gamma$ 、接触半径 $a$ に対し、臨界伸長長さ $\Delta l(a)$ は

$\Delta l(a)={\sqrt {\frac {2\pi \,a\,\Delta \gamma }{E_{c}^{*}}}}$

で定義される。付着モデルでは、材料は圧縮状態でのみ付着すると仮定し、アンロード時には現在の接触半径 $a$ と臨界接触半径 $a_{c}$ を比較しながら剥離条件を判定する。

臨界接触半径 $a_{c}$ は

$\left.{\frac {\mathrm {d} g_{1D}(a)}{\mathrm {d} a}}\right|_{a=a_{c}}=\xi \,{\sqrt {\frac {\pi \,\Delta \gamma }{2\,E_{c}^{*}\,a_{c}}}},\qquad \xi ={\begin{cases}3,&{\text{力制御}},\\1,&{\text{変位制御}},\end{cases}}$

を満たす値として決定される。ここで

$g_{1D}(a)={\frac {A}{2}}-{\frac {A}{B}}\,{\sqrt {{\frac {B^{2}}{4}}-a^{2}}}$

は楕円状インデンターに対する 1 次元ギャップ関数である。接触端の法線面変位

$w_{1D}(a)=\delta _{1D}^{e}-g_{1D}(a)$

と臨界伸長長さ $\Delta l(a)$ の相対関係により、以下の 3 つの付着状態が区別される。

張力ばねなし： $w_{1D}(a)=0$ の場合。付着力はゼロであり、非付着接触の力則 $F_{\mathrm {n.a.} }$ がそのまま適用される。
張力ばねあり（臨界長以下）： $0<w_{1D}(a)\leq \Delta l(a)$ の場合。接触の一部が張力状態にあるが剥離は生じない。このとき、全力は非付着寄与 $F_{\mathrm {n.a.} }$ と付着リトラクション寄与 $F_{\mathrm {a.r.} }$ の和として

$F_{\mathrm {MDR} }=F_{\mathrm {n.a.} }+F_{\mathrm {a.r.} },$ $F_{\mathrm {a.r.} }=2\,E_{c}^{*}\,{\bigl (}\delta _{1D}^{e}-g_{1D}(a){\bigr )}\,a$ と評価される。

張力ばねあり（臨界長超過）： $w_{1D}(a)>\Delta l(a)$ の場合。ばねの伸長が臨界値を超えたため剥離が必要となる。新しい平衡接触半径は

$\delta _{1D}^{e}+\Delta l(a)-g_{1D}(a)=0$ の解 $a$ を用いて決定し、 $a>a_{c}$ であれば接触保持、 $a<a_{c}$ であれば完全分離とみなす。

付着は一度分離しても再圧縮により再度発生し得る。付着力は MDR 部分の運動学量から計算され、バルク弾性寄与の定義には影響しない。

バルク弾性力の寄与

バルク弾性応答は、粒子全体の体積弾性率に基づく等方圧縮・膨張としてモデル化される。接触 $i$ におけるバルク寄与の増分は

$\mathrm {d} F_{\mathrm {Bulk} ,i}={\frac {A_{\rm {con}}}{V_{\rm {geo}}}}\,\kappa \,\mathrm {d} {\bar {\delta }}_{0}\,A_{C,i}$

で与えられる。ここで $V_{\rm {geo}}$ は幾何学的に変形した粒子体積であり、

$V_{\rm {geo}}={\frac {4}{3}}\pi R^{3}-\sum _{j=1}^{N}\left({\frac {\pi }{3}}\,\delta _{j}^{2}\,[3R-\delta _{j}]\right)$

と定義される。バルク応答が活性化している間のみ、 $\mathrm {d} {\bar {\delta }}_{0}$ の積分により $F_{\mathrm {Bulk} ,i}$ と $\delta _{0,\mathrm {Bulk} ,i}$ が更新される。

極端なオーバーラップと剛体平面配置

粒径分布や高圧縮によってオーバーラップ比 $\delta /R$ が 1 に近づくと、剛体平面を粒子内部に不自然に配置してしまう危険がある。DEM-MDR モデルでは、forward ローディング時の最大見掛けオーバーラップ $\delta _{i,{\rm {geo}}}^{\max }$ に上限を設け、オーバーラップが大きい側の粒子に対して上限値を課し、それ以上の変位を相手粒子側に割り当てることでこの問題を回避する。

例えば、粒子 $i$ 側のオーバーラップが

$\delta _{i,{\rm {geo}}}^{\max }/R_{i}>0.75$

となるとき、

$\delta _{i,{\rm {geo}}}^{\max }=0.75\,R_{i},\qquad \delta _{j,{\rm {geo}}}^{\max }=\delta _{\max }-0.75\,R_{i}$

と再配分する。同様に粒子 $j$ 側が支配的な場合には

$\delta _{i,{\rm {geo}}}^{\max }=\delta _{\max }-0.75\,R_{j},\qquad \delta _{j,{\rm {geo}}}^{\max }=0.75\,R_{j}$

とする。ここで用いられる上限値 $0.75$ は、原論文における数値例で一貫して用いられている経験的パラメータである。^[3]

トポロジカル・ペナルティ法

大変形 DEM では、粒子同士が他の粒子を「貫通」するような非物理的接触が数値的に発生し得る。DEM-MDR では、三粒子クラスターに基づくトポロジカル・ペナルティ法を導入し、このような接触に対してスカラー重みを付与してスクリーニングを行う^[3]。

ペナルティの概要

粒子 $i,j,k$ の三つ組を考えたとき、法線方向オーバーラップ $\delta _{ij},\delta _{ik},\delta _{jk}$ のいずれかが非正（ $\leq 0$ ）であれば、その三つ組に対するペナルティは全てゼロとする。

$\delta _{ij}\leq 0\;{\text{または}}\;\delta _{ik}\leq 0\;{\text{または}}\;\delta _{jk}\leq 0\Rightarrow p_{ij}=p_{ik}=p_{jk}=0.$

全てが正の場合、中心粒子を決定するために中心間距離 $d_{ij},d_{ik},d_{jk}$ の最大値を取り、その対を非中心接触と看做す。例えば

$\max(d_{ij},\,d_{ik},\,d_{jk})=d_{ij}$

であれば粒子 $k$ を中心粒子とし、 $i\!-\!k$ 、 $j\!-\!k$ を中心接触、 $i\!-\!j$ を非中心接触とする。

中心粒子から他粒子への単位法線ベクトル ${\boldsymbol {e}}_{ki},\,{\boldsymbol {e}}_{kj}$ に対し、

$\alpha _{kij}=\arccos({\boldsymbol {e}}_{ki}\cdot {\boldsymbol {e}}_{kj})$

を定義し、非中心接触にはシグモイド型ペナルティ関数

${\hat {p}}(\alpha )={\frac {1}{1+\exp \!\left(-t_{p}{\frac {\alpha }{\pi }}-{\tfrac {1}{2}}\right)}}$

を適用する。ここで $t_{p}$ は勾配を決めるパラメータであり、数値例では $t_{p}=50$ が用いられている。したがって

$p_{ij}={\hat {p}}(\alpha _{kij}),\qquad p_{ik}=p_{jk}=0$

となる。

接触重みへの変換

最終的に、各接触に対する重み $w$ は

$w_{ij}=\max(1-p_{ij},0),\quad w_{ik}=\max(1-p_{ik},0),\quad w_{jk}=\max(1-p_{jk},0)$

と定義される。大規模系では、全ての有効接触ペアについて同様の手順により重みを計算し、接触力 $F$ に $w\in [0,1]$ を乗じて非物理的接触を部分的または完全に無効化する。

数値実装とメインアルゴリズム

単一粒子上の $N$ 個の有効接触について、DEM-MDR モデルに基づく接触力計算は明示的時間積分スキームとして実装される。既知の物性値として、初期半径 $R_{0}$ 、複合平面ひずみ弾性率 $E_{c}^{*}$ 、ポアソン比 $\nu$ 、体積弾性率 $\kappa$ 、降伏応力 $Y$ 、有効表面エネルギー $\Delta \gamma$ 、臨界拘束比 $\psi _{b}$ が与えられる。

計算ステップ $n$ の入力としては、各接触の見掛け重なり ${\underline {\delta }}^{n}$ 、MDR 見掛け重なり ${\underline {\delta }}_{\mathrm {MDR} }^{n}$ およびその最大値 ${\underline {\delta }}_{\mathrm {max,MDR} }^{n}$ 、バルク変位 ${\underline {\delta }}_{0,\mathrm {Bulk} }^{n}$ 、見掛け半径 $R^{n}$ 、粒子体積 $V^{n}$ 、見掛け重なり増分 $\Delta {\underline {\delta }}$ 、変位増分 $\Delta {\underline {\delta }}_{0}$ などが必要となる。

メインアルゴリズムは、概ね以下の手順から構成される。

入力増分を用いて新しい見掛け重なり ${\underline {\delta }}^{n+1}$ を計算する。
降伏フラグ ${\underline {Y}}_{\mathrm {flag} }$ と接触面積切片 ${\underline {c}}_{A}$ に基づき、各接触の接触面積 $A_{C,i}^{n+1}$ を更新し、 $A_{\mathrm {con} }$ 、 $A_{\mathrm {tot} }$ 、 $A_{\mathrm {free} }$ 、 $V_{\mathrm {geo} }$ を評価する。
平均変位増分 $\Delta {\bar {\delta }}_{0}$ を求め、 $A_{\mathrm {free} }/A_{\mathrm {tot} }$ の大小に応じて MDR 成分とバルク成分への変位の分配を行い、 ${\underline {\delta }}_{\mathrm {MDR} }^{n+1}$ 、 ${\underline {\delta }}_{0,\mathrm {Bulk} }^{n+1}$ 、 ${\underline {\delta }}_{\mathrm {max,MDR} }^{n+1}$ を更新する。
各接触について、降伏判定（ ${\bar {p}}_{H}>p_{Y}$ ）と弾性／塑性の区別に基づき $A,B,\delta _{1D}^{e}$ を決定し、非付着の場合は MDR 力 $F_{\mathrm {MDR} }$ を、付着が有効な場合は前節の JKR 型付着則に従って $F_{\mathrm {MDR} }$ を計算する。
バルク弾性寄与 $F_{\mathrm {Bulk} }$ を $A_{\mathrm {con} },V_{\mathrm {geo} },\kappa ,\delta _{0,\mathrm {Bulk} }$ から評価し、 $F_{i}^{n+1}=F_{\mathrm {MDR} }+F_{\mathrm {Bulk} }$ を得る。
更新された接触力ベクトル ${\boldsymbol {f}}^{n+1}$ とブランチベクトル ${\boldsymbol {b}}^{n+1}$ を用いて体積 $V^{n+1}$ を計算し、弾性体積変化 $\Delta V^{e}$ と半径増分 $\Delta R$ を通じて見掛け半径 $R^{n+1}$ を更新する。

このように、MDR 接触応答、付着破壊、バルク弾性応答、粒子半径の更新が一つのアルゴリズム内で整合的に扱われる。

記号一覧

以下に、DEM-MDR モデルで用いられる主な記号をまとめる。

記号	定義	記号	定義
$A$	楕円インデンターの高さ	$F$	接触力
$A_{\rm {con}}$	接触に関与する面積	$F_{\rm {MDR}}$	MDR による接触力
$A_{\rm {free}}$	自由表面積	$F_{\rm {Bulk}}$	バルク弾性応答による力
$A_{\rm {tot}}$	粒子表面積の総和	$F_{\rm {n.a.}}$	非付着接触における力
$A_{C}$	個々の接触面積	$F_{\rm {a.r.}}$	付着リトラクションによる力
$a$	接触半径	${\boldsymbol {f}}$	接触力ベクトル
$a_{c}$	臨界接触半径	$G$	せん断弾性率
$a_{\max }$	最大経験接触半径	$g_{1\mathrm {D} }$	1 次元ギャップ関数
$B$	楕円インデンターの幅	$\Delta \gamma$	有効表面エネルギー
${\boldsymbol {b}}$	ブランチベクトル	$\kappa$	体積弾性率
$c_{A}$	接触面積切片	$\Delta \ell$	臨界伸長長さ
$\delta _{0}$	見掛け変位	${\boldsymbol {n}}$	接触法線ベクトル
${\boldsymbol {\delta }}_{0}$	見掛け変位ベクトル	$\nu$	ポアソン比
${\bar {\delta }}_{0}$	平均変位	${\bar {p}}_{H}$	Hertz 接触の平均圧力
$\delta _{0,{\rm {MDR}}}$	MDR 変位成分	$p_{Y}$	硬化曲線に沿った平均圧力
$\delta _{0,{\rm {Bulk}}}$	バルク変位成分	$\psi _{b}$	臨界拘束比
$\delta _{0,\max }$	最大経験変位	$R_{0}$	初期半径
$\delta$	見掛けオーバーラップ	$R$	見掛け半径
$\delta _{\rm {MDR}}$	MDR 見掛けオーバーラップ	$\Delta R$	見掛け半径の増分
$\delta _{\max ,{\rm {MDR}}}$	MDR 最大見掛けオーバーラップ	$V_{0}$	初期体積
$\delta _{1\mathrm {D} }^{e}$	変換された弾性変位	$V$	現在の体積
$\delta _{1\mathrm {D} ,\max }^{e}$	最大変換弾性変位	$V_{\rm {geo}}$	幾何学的に変形した粒子体積
$\delta _{R}$	変位補正項	$\Delta V^{e}$	弾性体積変化
$\delta _{Y}$	降伏変位	$w_{1\mathrm {D} }$	1 次元法線変位
$E$	ヤング率	$Y$	降伏応力
$E_{c}^{*}$	複合平面ひずみ弾性率	$z_{R}$	粒子中心の深さ

脚注

^ ^a ^b Zunker, William; Kamrin, Ken (2024). “A mechanically-derived contact model for adhesive elastic-perfectly plastic particles, Part I: Utilizing the method of dimensionality reduction”. Journal of the Mechanics and Physics of Solids 183: 105492. doi:10.1016/j.jmps.2023.105492.
^ ^a ^b Zunker, William; Kamrin, Ken (2024). “A mechanically-derived contact model for adhesive elastic-perfectly plastic particles, Part II: Contact under high compaction—modeling a bulk elastic response”. Journal of the Mechanics and Physics of Solids 183: 105493. doi:10.1016/j.jmps.2023.105493.
^ ^a ^b ^c ^d Zunker, William; Dunatunga, Sachith; Thakur, Subhash; Tang, Pingjun; Kamrin, Ken (2025). “Experimentally validated DEM for large deformation powder compaction: Mechanically-derived contact model and screening of non-physical contacts”. Powder Technology 459: 120972. doi:10.1016/j.powtec.2025.120972.
^ “LAMMPS molecular dynamics simulator”. GitHub. 2025年12月2日閲覧。

参考文献

Zunker, William; Kamrin, Ken (2024). “A mechanically-derived contact model for adhesive elastic-perfectly plastic particles, Part I: Utilizing the method of dimensionality reduction”. Journal of the Mechanics and Physics of Solids 183: 105492. doi:10.1016/j.jmps.2023.105492.
Zunker, William; Kamrin, Ken (2024). “A mechanically-derived contact model for adhesive elastic-perfectly plastic particles, Part II: Contact under high compaction—modeling a bulk elastic response”. Journal of the Mechanics and Physics of Solids 183: 105493. doi:10.1016/j.jmps.2023.105493.
Zunker, William; Dunatunga, Sachith; Thakur, Subhash; Tang, Pingjun; Kamrin, Ken (2025). “Experimentally validated DEM for large deformation powder compaction: Mechanically-derived contact model and screening of non-physical contacts”. Powder Technology 459: 120972. doi:10.1016/j.powtec.2025.120972.
“LAMMPS molecular dynamics simulator”. GitHub. 2025年12月2日閲覧。

[ZunkerKamrinI-1] Zunker, William; Kamrin, Ken (2024). “A mechanically-derived contact model for adhesive elastic-perfectly plastic particles, Part I: Utilizing the method of dimensionality reduction”. Journal of the Mechanics and Physics of Solids 183: 105492. doi:10.1016/j.jmps.2023.105492.

[ZunkerKamrinII-2] Zunker, William; Kamrin, Ken (2024). “A mechanically-derived contact model for adhesive elastic-perfectly plastic particles, Part II: Contact under high compaction—modeling a bulk elastic response”. Journal of the Mechanics and Physics of Solids 183: 105493. doi:10.1016/j.jmps.2023.105493.

[ZunkerPowderTech-3] Zunker, William; Dunatunga, Sachith; Thakur, Subhash; Tang, Pingjun; Kamrin, Ken (2025). “Experimentally validated DEM for large deformation powder compaction: Mechanically-derived contact model and screening of non-physical contacts”. Powder Technology 459: 120972. doi:10.1016/j.powtec.2025.120972.

[4] “LAMMPS molecular dynamics simulator”. GitHub. 2025年12月2日閲覧。

[1]

[2]

[3]

[4]