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

**Bonded particle model**（ボンデッドパーティクルモデル、*bonded particle model*、BPM）は、粒子間にばねとダッシュポットからなる結合を導入することにより、岩石・脆性材料などの破壊挙動を表現する離散要素法の一種である。

ここでは、2 粒子間ボンドの力学的モデル化について述べる。

## 基本構成

### 結合断面積と断面二次モーメント

2つの粒子 i と j の結合点に対して，接合内側半径 r_{\text{in}}，接合外側半径 r_{\text{out}} を粒子 i と粒子 j の粒子半径 R_i と R_j を用いて次のように定義する。

: r_{\text{in}} = r_i \times \min(R_i, R_j) 

: r_{\text{out}} = r_o \times \min(R_i, R_j) 

結合断面積 A_{\text{sec}} は

: A_{\text{sec}} = \pi \left( r_{\text{out}}^2 - r_{\text{in}}^2 \right) 

極断面二次モーメント J_{\text{sec}} は

: J_{\text{sec}} = \frac{1}{2} A_{\text{sec}} \left( r_{\text{out}}^2 - r_{\text{in}}^2 \right) 

極断面二次モーメント（断面全体に対してねじりに抵抗する抵抗）は

: I_p = \frac{1}{2} \pi \left( r_{\text{out}}^4 - r_{\text{in}}^4 \right) 

主断面二次モーメント（断面をある軸まわりに曲げモーメントがかかったときの抵抗）は

: I_{ph} = \frac{I_p}{2} 

### ボンド長と各種剛性定数

ボンド長さ L_b は

: L_b = \ell_b \bigl( R_i + R_j \bigr) 

法線バネ定数 K_n は

: K_n = S_n \frac{A_{\mathrm{sec}}}{L_b} 

接線バネ定数 K_t は

: K_t = S_t \frac{A_{\mathrm{sec}}}{L_b} 

ねじりバネ定数 K_{\mathrm{tor}} は

: K_{\mathrm{tor}} = S_{\mathrm{tor}} \frac{I_p}{L_b} 

曲げバネ定数 K_{\mathrm{ben}} は

: K_{\mathrm{ben}} = S_{\mathrm{ben}} \frac{I_{ph}}{L_b} 

ここで \ell_b, S_n, S_t, S_{\mathrm{tor}}, S_{\mathrm{ben}} は定数である。Potyondy は S_n = S_{\mathrm{tor}}, S_t = S_{\mathrm{ben}} とし，S_n と S_t の関係はラメ定数で決めた。また，\sigma_c = \tau_c と仮定した。接合ヤング率は

: E_b = \frac{K_n}{L_b} 

と定義した。Horabik は 0.25 ，0.01  で条件を振って脆性・半脆性・延性を表現できると報告している。

### 有効質量と慣性モーメント

各粒子の質量は

: m_i = \frac{4}{3} \pi R_i^3 \rho_i 

: m_j = \frac{4}{3} \pi R_j^3 \rho_j 

有効質量 M_e は

: M_e = \frac{m_i m_j}{m_i + m_j} 

各粒子の極慣性モーメントは

: J_i = \frac{2}{5} m_i R_i^2 

: J_j = \frac{2}{5} m_j R_j^2 

有効極慣性モーメント J_s は

: J_s = \frac{J_i J_j}{J_i + J_j} 

## 相対運動の分解

粒子間ベクトルは

: \Delta \boldsymbol{x} = [x_i - x_j,\, y_i - y_j,\, z_i - z_j]^{\mathrm{T}} 

その長さは

: r = \|\Delta \boldsymbol{x}\| = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2} 

法線方向の単位ベクトルは

: \boldsymbol{e}_n = \frac{\Delta \boldsymbol{x}}{r} 

### 直線運動による相対速度

相対速度ベクトルは

: \boldsymbol{v}_r = \boldsymbol{v}_i - \boldsymbol{v}_j 

法線方向成分は

: v_{\text{nnr}} = \boldsymbol{v}_r \cdot \Delta \boldsymbol{x} 

: \boldsymbol{v}_n = (\boldsymbol{v}_r \cdot \boldsymbol{e}_n)\, \boldsymbol{e}_n = v_{\text{nnr}} \frac{\Delta \boldsymbol{x}}{r^2} 

接線方向成分は

: \boldsymbol{v}_t = \boldsymbol{v}_r - \boldsymbol{v}_n 

### 回転運動による相対速度の補正

回転によるせん断相対速度の補正角速度は

: \boldsymbol{w}_r = \frac{R_i \boldsymbol{\omega}_i + R_j \boldsymbol{\omega}_j}{r} 

接線補正後の相対速度は

: \boldsymbol{v}_{tr} = \boldsymbol{v}_t - (\Delta \boldsymbol{x} \times \boldsymbol{w}_r) 

### 曲げ・ねじり回転の相対速度

曲げとねじりの相対角速度は

: \boldsymbol{w}_r' = \boldsymbol{\omega}_i - \boldsymbol{\omega}_j 

ねじり成分は

: w_{\text{nnr}} = \boldsymbol{w}_r' \cdot \Delta \boldsymbol{x} 

: \boldsymbol{w}_n = (\boldsymbol{w}_r' \cdot \boldsymbol{e}_n)\, \boldsymbol{e}_n = w_{\text{nnr}} \frac{\Delta \boldsymbol{x}}{r^2} 

曲げ成分は

: \boldsymbol{w}_t = \boldsymbol{w}_r' - \boldsymbol{w}_n 

## 力とモーメント

### 法線方向ばね力

伸縮量（ひずみ）は

: \varepsilon = \frac{r - L_b}{r} 

法線ばね力は

: \boldsymbol{F}_n = -K_n \varepsilon\, \Delta \boldsymbol{x} 

指数補正を加味する場合は

: \boldsymbol{F}_n = -K_n \varepsilon \frac{1}{n_l (1 - r / b_n)} \Delta \boldsymbol{x} 

ここで

: n_l = \frac{1}{1 - L_b / b_n} 

である。

### 接線方向ばね力（せん断）

: \Delta \boldsymbol{F}_t = -K_t \Delta t\, \boldsymbol{v}_{tr} 

### ねじりトルク・曲げトルク

ねじりトルクは

: \Delta \boldsymbol{\tau}_n = -K_{\text{tor}} \Delta t\, \boldsymbol{w}_n 

曲げトルクは

: \Delta \boldsymbol{\tau}_t = -K_{\text{ben}} \Delta t\, \boldsymbol{w}_t 

### ダンピング項

法線ダンピング力は

: \boldsymbol{F}_{\text{damp},n} = 2 \zeta \sqrt{M_e K_n}\, (-\boldsymbol{v}_n) 

接線ダンピング力は

: \boldsymbol{F}_{\text{damp},t} = 2 \zeta \sqrt{M_e K_t}\, (-\boldsymbol{v}_{tr}) 

ねじりダンピングトルクは

: \boldsymbol{\tau}_{\text{damp},n} = 2 \zeta \sqrt{J_s K_{\text{tor}}}\, (-\boldsymbol{w}_n) 

曲げダンピングトルクは

: \boldsymbol{\tau}_{\text{damp},t} = 2 \zeta \sqrt{J_s K_{\text{ben}}}\, (-\boldsymbol{w}_t) 

## 履歴ベクトルの回転（ベクトル表記）

粒子間の結合軸 \Delta \boldsymbol{x} が毎ステップ傾くたび，前ステップの履歴 \boldsymbol{B}^{(n)} を新しい軸に射影し直す操作を「回転」と呼ぶ。ここで「'」や「*」のついたベクトルは計算用の一時的なベクトル変数である。

=== 定義 ===

履歴ベクトルを

:<math> \boldsymbol{B}^{(n)} = [\boldsymbol{F}_n,\; \boldsymbol{F}_t,\; \boldsymbol{\tau}_n,\; \boldsymbol{\tau}_t]^{\mathrm{T}} \in \mathbb{R}^{12} </math>

とする（上付き <math>(n)</math> は時刻 <math>t_n</math> を示す）。

=== 増分加算（オイラー陽） ===

力とトルクの更新にはバネだけ考え，ダッシュポットは考えない。

:<math>
\boldsymbol{B}^{(n+1)}
=
\bigl(
\boldsymbol{F}_n,\;
\boldsymbol{F}_t + \Delta \boldsymbol{F}_t,\;
\boldsymbol{\tau}_n + \Delta \boldsymbol{\tau}_n,\;
\boldsymbol{\tau}_t + \Delta \boldsymbol{\tau}_t
\bigr)
</math>

これは

:<math> \boldsymbol{B}^{(n+1)} = \boldsymbol{B}^{(n)} + \Delta \boldsymbol{B} </math>

と書き換えられる。

=== 法線力の回転（未使用） ===

法線力については毎ステップ重なり量から新しく計算し直すだけなので，過去の大きさを温存する必要がない。

=== せん断力の回転 ===

全体 <math>\boldsymbol{F}_t^{(n)}</math> から法線成分を引くと接線成分になる。これは直交化（Gram–Schmidt 法）と呼ばれる。履歴ばねとして「前ステップまでに蓄えられたせん断バネ伸び（＝エネルギー）」をそっくり次ステップへ持ち越す必要がある。

:<math>
\boldsymbol{F}_t'
=
\boldsymbol{F}_t^{(n)} - (\boldsymbol{F}_t^{(n)} \cdot \boldsymbol{e}_n)\, \boldsymbol{e}_n
</math>

:<math>
\boldsymbol{F}_t* =
\begin{cases}
\boldsymbol{0}, & \|\boldsymbol{F}_t'\| = 0,\\[4pt]
\|\boldsymbol{F}_t^{(n)}\| \dfrac{\boldsymbol{F}_t'}{\|\boldsymbol{F}_t'\|}, & \text{otherwise}.
\end{cases}


### ねじりモーメントの回転

ここでは \boldsymbol{e}_n ではなく \boldsymbol{\tau}_n'/\|\boldsymbol{\tau}_n'\| を使う。\boldsymbol{e}_n を使ってしまうと符号が逆になる可能性があるためである。履歴ばねとして「前ステップまでに蓄えられたねじりバネ伸び（＝エネルギー）」をそっくり次ステップへ持ち越す必要がある。

:
\boldsymbol{\tau}_n'
=
(\boldsymbol{\tau}_n^{(n)} \cdot \boldsymbol{e}_n)\, \boldsymbol{e}_n


:
\boldsymbol{\tau}_n* =
\begin{cases}
\boldsymbol{0}, & \|\boldsymbol{\tau}_n'\| = 0,\\[4pt]
\|\boldsymbol{\tau}_n^{(n)}\| \dfrac{\boldsymbol{\tau}_n'}{\|\boldsymbol{\tau}_n'\|}, & \text{otherwise}.
\end{cases}
</math>

=== 曲げモーメントの回転 ===

全体 <math>\boldsymbol{\tau}_t^{(n)}</math> から法線成分を引くと接線成分になる。これも直交化（Gram–Schmidt 法）であり，履歴ばねとして「前ステップまでに蓄えられた曲げバネ伸び（＝エネルギー）」をそっくり次ステップへ持ち越す必要がある。

:<math>
\boldsymbol{\tau}_t'
=
\boldsymbol{\tau}_t^{(n)} - (\boldsymbol{\tau}_n^{(n)} \cdot \boldsymbol{e}_n)\, \boldsymbol{e}_n
</math>

:<math>
\boldsymbol{\tau}_t* =
\begin{cases}
\boldsymbol{0}, & \|\boldsymbol{\tau}_t'\| = 0,\\[4pt]
\|\boldsymbol{\tau}_t^{(n)}\| \dfrac{\boldsymbol{\tau}_t'}{\|\boldsymbol{\tau}_t'\|}, & \text{otherwise}.
\end{cases}


### 最新増分の加算

最新ステップで計算された増分を加えて履歴ベクトルを更新する。

:
\begin{aligned}
\boldsymbol{F}_n^{(n+1)} &\leftarrow \boldsymbol{F}_n,\\
\boldsymbol{F}_t^{(n+1)} &\leftarrow \boldsymbol{F}_t* + \Delta \boldsymbol{F}_t,\\
\boldsymbol{\tau}_n^{(n+1)} &\leftarrow \boldsymbol{\tau}_n* + \Delta \boldsymbol{\tau}_n,\\
\boldsymbol{\tau}_t^{(n+1)} &\leftarrow \boldsymbol{\tau}_t'' + \Delta \boldsymbol{\tau}_t.
\end{aligned}


## 結合破壊判定

力とトルクによる破壊判定にはバネだけ考え，ダッシュポットは考えない。\sigma_c，\tau_c はそれぞれ臨界法線応力，臨界接線応力である。

### 法線応力破壊

:
\sigma_n
=
\frac{\|\boldsymbol{F}_n\|}{A_{\mathrm{sec}}}
+
\frac{2 \|\boldsymbol{\tau}_t\| (r_{\mathrm{out}} - r_{\mathrm{in}})}{J_{\mathrm{sec}}}


破壊条件は

: \sigma_n > \sigma_c 

である。

### せん断応力破壊

:
\tau_t
=
\frac{\|\boldsymbol{F}_t\|}{A_{\mathrm{sec}}}
+
\frac{\|\boldsymbol{\tau}_n\| (r_{\mathrm{out}} - r_{\mathrm{in}})}{J_{\mathrm{sec}}}


破壊条件は

: \tau_t > \tau_c 

である。

## 粒子方程式への寄与

接触 (i,j) から粒子 i に作用する合力・合トルクは

:
\boldsymbol{f}_i
=
\sum_{\mathrm{contact}}
\bigl(
\boldsymbol{F}_n
+
\boldsymbol{F}_t
+
\boldsymbol{F}_{\mathrm{damp},n}
+
\boldsymbol{F}_{\mathrm{damp},t}
\bigr)


:
\boldsymbol{\tau}_i
=
\sum_{\mathrm{contact}}
\bigl(
R_i (-\boldsymbol{e}_n \times \boldsymbol{F}_t)
+
\boldsymbol{\tau}_n
+
\boldsymbol{\tau}_t
+
\boldsymbol{\tau}_{\mathrm{damp},n}
+
\boldsymbol{\tau}_{\mathrm{damp},t}
\bigr)


接触 (i,j) から粒子 j に作用する合力・合トルクは

:
\boldsymbol{f}_j
=
-\sum_{\mathrm{contact}}
\bigl(
\boldsymbol{F}_n
+
\boldsymbol{F}_t
+
\boldsymbol{F}_{\mathrm{damp},n}
+
\boldsymbol{F}_{\mathrm{damp},t}
\bigr)


:
\boldsymbol{\tau}_j
=
\sum_{\mathrm{contact}}
\bigl(
R_j (-\boldsymbol{e}_n \times \boldsymbol{F}_t)
-
\boldsymbol{\tau}_n
-
\boldsymbol{\tau}_t
-
\boldsymbol{\tau}_{\mathrm{damp},n}
-
\boldsymbol{\tau}_{\mathrm{damp},t}
\bigr)


## 脚注