DEMにおけるVoigtモデル
図1に示すように離散要素法（DEM）では，接触中の粒子対 i–j を法線方向と接線方向に直交分解し，それぞれを線形ばねとダッシュポットを並列接続した Voigt 要素で近似する［Cundall & Strack (1979)］。圧縮を正とし，接触法線単位ベクトルを \mathbf{n} とする。
図1：Voigtモデルにおける粒子—粒子接触の概念図（(a) 法線相互作用，(b) 接線相互作用）。
法線方向
重なり \delta_n と法線相対速度 v_n=(\mathbf{v}_{ij}\cdot\mathbf{n}) から，法線力は
\mathbf{F}_n = \bigl(k_n\,\delta_n - \gamma_n\, v_n \bigr)\,\mathbf{n},
で与える。ここで k_n は法線ばね定数，\gamma_n は法線減衰係数である。
接線方向
Cundall–Strack に従い，接触座標系で積分した接線ばね伸び \boldsymbol{\xi}_t を用いる（図1）。接線相対速度は
\mathbf{v}_t \;=\; \mathbf{v}_{ij} - (\mathbf{v}_{ij}\!\cdot\!\mathbf{n})\,\mathbf{n}
\;-\; R_i\,(\boldsymbol{\omega}_i\!\times\!\mathbf{n})
\;-\; R_j\,(\boldsymbol{\omega}_j\!\times\!\mathbf{n}),
で与える。接線力の試行値と最終値は
\mathbf{F}_t^\ast = -k_t\,\boldsymbol{\xi}_t - \gamma_t\,\mathbf{v}_t,
\mathbf{F}_t = -\,\min\!\Bigl(\,\|\mathbf{F}_t^\ast\|,\ \mu\,|\mathbf{F}_n|\,\Bigr)\,
\frac{\mathbf{F}_t^\ast}{\|\mathbf{F}_t^\ast\|},
で与える。ここで \mu はクーロン摩擦係数である。
転がり摩擦
球粒子は接線摩擦のみでは過大回転しやすいため，接触点まわりの転がり抵抗モーメント \boldsymbol{\tau}_r を導入して回転を抑制する［Iwashita & Oda (1998)］。有効半径は R^\ast=R_iR_j/(R_i+R_j)，相対角速度は \boldsymbol{\omega}_\mathrm{rel}=\boldsymbol{\omega}_i-\boldsymbol{\omega}_j と表される。回転ばね角度履歴 \boldsymbol{\theta}_r を用いて
\boldsymbol{\tau}_r^\ast = -k_r\,\boldsymbol{\theta}_r - \gamma_r\,\boldsymbol{\omega}_\mathrm{rel},
\boldsymbol{\tau}_r = -\,\min\!\Bigl(\,\|\boldsymbol{\tau}_r^\ast\|,\ \mu_r\,|\mathbf{F}_n|\,R^\ast\,\Bigr)\,
\frac{\boldsymbol{\tau}_r^\ast}{\|\boldsymbol{\tau}_r^\ast\|}.
静的・動的の転がり抵抗を統一的に表現でき，安息角や回転率の再現性が高い。
時間積分
本研究では，粒子の並進・回転運動をシンプレクティック・オイラー法で時間積分する。粒子 i の質量を m_i，慣性モーメントを I_i，位置を \mathbf{r}_i，速度を \mathbf{v}_i，角速度を \boldsymbol{\omega}_i，回転角（姿勢パラメータ）を \boldsymbol{\theta}_i，作用する合力・合モーメントをそれぞれ \mathbf{F}_i,\ \boldsymbol{\tau}_i とすると，運動方程式は
m_i\,\dot{\mathbf{v}}_i = \mathbf{F}_i,\quad
\dot{\mathbf{r}}_i = \mathbf{v}_i,\quad
I_i\,\dot{\boldsymbol{\omega}}_i = \boldsymbol{\tau}_i,\quad
\dot{\boldsymbol{\theta}}_i = \boldsymbol{\omega}_i.
シンプレクティック・オイラー法による更新は，まず速度・角速度を更新し，その新しい値で位置・角度を更新する：
\mathbf{v}_i^{\,n+1} = \mathbf{v}_i^{\,n} + \Delta t\,\frac{\mathbf{F}_i^{\,n}}{m_i},\qquad
\boldsymbol{\omega}_i^{\,n+1} = \boldsymbol{\omega}_i^{\,n} + \Delta t\,\frac{\boldsymbol{\tau}_i^{\,n}}{I_i},
\mathbf{r}_i^{\,n+1} = \mathbf{r}_i^{\,n} + \Delta t\,\mathbf{v}_i^{\,n+1},\qquad
\boldsymbol{\theta}_i^{\,n+1} = \boldsymbol{\theta}_i^{\,n} + \Delta t\,\boldsymbol{\omega}_i^{\,n+1}.
シミュレーション条件
本研究のシミュレーションに用いた主要パラメータを表1に示す。
表1：Simulation parameters.
| Symbol |
Parameter |
Value |
Units |
| \rho |
density |
2.50e3 |
kg/m^3 |
| R |
particle radius (all same) |
1.00e-3 |
m |
| k_n |
normal spring |
1.00e4 |
N/m |
| \gamma_n |
normal dashpot |
1.00e-3 |
kg/s |
| k_t |
tangential spring |
6.00e3 |
N/m |
| \gamma_t |
tangential dashpot |
1.00e-3 |
kg/s |
| \mu |
Coulomb friction |
5.00e-1 |
— |
| \eta_r |
viscous rolling resistance |
2.00e-8 |
N·m·s |
| k_r |
rotational spring |
1.00e-6 |
N·m/rad |
| \gamma_r |
rotational dashpot |
1.00e-8 |
N·m·s/rad |
| \mu_r |
rolling friction |
5.00e-2 |
— |
参考文献
- Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Géotechnique, 29(1), 47–65.
- Iwashita, K., & Oda, M. (1998). Rolling resistance at contacts in simulation of shear band development by DEM. Journal of Engineering Mechanics, 124(3), 285–292.
Discussion