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

Summarize Timeline Top Qs Fact Check

ファイル:BPM system schematic.png

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

${\displaystyle r_{\text{in}}=r_{i}\times \min(R_{i},R_{j})}$

${\displaystyle r_{\text{out}}=r_{o}\times \min(R_{i},R_{j})}$

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

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

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

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

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

${\displaystyle I_{p}={\frac {1}{2}}\pi \left(r_{\text{out}}^{4}-r_{\text{in}}^{4}\right)}$

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

${\displaystyle I_{ph}={\frac {I_{p}}{2}}}$

ボンド長さ ${\displaystyle L_{b}}$ は

${\displaystyle L_{b}=\ell _{b}{\bigl (}R_{i}+R_{j}{\bigr )}}$

法線バネ定数 ${\displaystyle K_{n}}$ は

${\displaystyle K_{n}=S_{n}{\frac {A_{\mathrm {sec} }}{L_{b}}}}$

接線バネ定数 ${\displaystyle K_{t}}$ は

${\displaystyle K_{t}=S_{t}{\frac {A_{\mathrm {sec} }}{L_{b}}}}$

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

${\displaystyle K_{\mathrm {tor} }=S_{\mathrm {tor} }{\frac {I_{p}}{L_{b}}}}$

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

${\displaystyle K_{\mathrm {ben} }=S_{\mathrm {ben} }{\frac {I_{ph}}{L_{b}}}}$

ここで ${\displaystyle \ell _{b}}$, ${\displaystyle S_{n}}$, ${\displaystyle S_{t}}$, ${\displaystyle S_{\mathrm {tor} }}$, ${\displaystyle S_{\mathrm {ben} }}$ は定数である。Potyondy は ${\displaystyle S_{n}=S_{\mathrm {tor} }}$, ${\displaystyle S_{t}=S_{\mathrm {ben} }}$ とし，${\displaystyle S_{n}}$ と ${\displaystyle S_{t}}$ の関係はラメ定数で決めた。また，${\displaystyle \sigma _{c}=\tau _{c}}$ と仮定した。接合ヤング率は

${\displaystyle E_{b}={\frac {K_{n}}{L_{b}}}}$

と定義した^[1]。Horabik は ${\displaystyle 0.25<\sigma _{c}/\tau _{c}<10}$，${\displaystyle 0.01<E_{b}/E<1}$ で条件を振って脆性・半脆性・延性を表現できると報告している^[2]。

各粒子の質量は

${\displaystyle m_{i}={\frac {4}{3}}\pi R_{i}^{3}\rho _{i}}$

${\displaystyle m_{j}={\frac {4}{3}}\pi R_{j}^{3}\rho _{j}}$

有効質量 ${\displaystyle M_{e}}$ は

${\displaystyle M_{e}={\frac {m_{i}m_{j}}{m_{i}+m_{j}}}}$

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

${\displaystyle J_{i}={\frac {2}{5}}m_{i}R_{i}^{2}}$

${\displaystyle J_{j}={\frac {2}{5}}m_{j}R_{j}^{2}}$

有効極慣性モーメント ${\displaystyle J_{s}}$ は

${\displaystyle J_{s}={\frac {J_{i}J_{j}}{J_{i}+J_{j}}}}$

Summarize Timeline Top Qs Fact Check

粒子間ベクトルは

${\displaystyle \Delta {\boldsymbol {x}}=[x_{i}-x_{j},\,y_{i}-y_{j},\,z_{i}-z_{j}]^{\mathrm {T} }}$

その長さは

${\displaystyle r=\|\Delta {\boldsymbol {x}}\|={\sqrt {(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}+(z_{i}-z_{j})^{2}}}}$

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

${\displaystyle {\boldsymbol {e}}_{n}={\frac {\Delta {\boldsymbol {x}}}{r}}}$

相対速度ベクトルは

${\displaystyle {\boldsymbol {v}}_{r}={\boldsymbol {v}}_{i}-{\boldsymbol {v}}_{j}}$

法線方向成分は

${\displaystyle v_{\text{nnr}}={\boldsymbol {v}}_{r}\cdot \Delta {\boldsymbol {x}}}$

${\displaystyle {\boldsymbol {v}}_{n}=({\boldsymbol {v}}_{r}\cdot {\boldsymbol {e}}_{n})\,{\boldsymbol {e}}_{n}=v_{\text{nnr}}{\frac {\Delta {\boldsymbol {x}}}{r^{2}}}}$

接線方向成分は

${\displaystyle {\boldsymbol {v}}_{t}={\boldsymbol {v}}_{r}-{\boldsymbol {v}}_{n}}$

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

${\displaystyle {\boldsymbol {w}}_{r}={\frac {R_{i}{\boldsymbol {\omega }}_{i}+R_{j}{\boldsymbol {\omega }}_{j}}{r}}}$

接線補正後の相対速度は

${\displaystyle {\boldsymbol {v}}_{tr}={\boldsymbol {v}}_{t}-(\Delta {\boldsymbol {x}}\times {\boldsymbol {w}}_{r})}$

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

${\displaystyle {\boldsymbol {w}}_{r}'={\boldsymbol {\omega }}_{i}-{\boldsymbol {\omega }}_{j}}$

ねじり成分は

${\displaystyle w_{\text{nnr}}={\boldsymbol {w}}_{r}'\cdot \Delta {\boldsymbol {x}}}$

${\displaystyle {\boldsymbol {w}}_{n}=({\boldsymbol {w}}_{r}'\cdot {\boldsymbol {e}}_{n})\,{\boldsymbol {e}}_{n}=w_{\text{nnr}}{\frac {\Delta {\boldsymbol {x}}}{r^{2}}}}$

曲げ成分は

${\displaystyle {\boldsymbol {w}}_{t}={\boldsymbol {w}}_{r}'-{\boldsymbol {w}}_{n}}$

Summarize Timeline Top Qs Fact Check

伸縮量（ひずみ）は

${\displaystyle \varepsilon ={\frac {r-L_{b}}{r}}}$

法線ばね力は

${\displaystyle {\boldsymbol {F}}_{n}=-K_{n}\varepsilon \,\Delta {\boldsymbol {x}}}$

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

${\displaystyle {\boldsymbol {F}}_{n}=-K_{n}\varepsilon {\frac {1}{n_{l}(1-r/b_{n})}}\Delta {\boldsymbol {x}}}$

ここで

${\displaystyle n_{l}={\frac {1}{1-L_{b}/b_{n}}}}$

である。

${\displaystyle \Delta {\boldsymbol {F}}_{t}=-K_{t}\Delta t\,{\boldsymbol {v}}_{tr}}$

ねじりトルクは

${\displaystyle \Delta {\boldsymbol {\tau }}_{n}=-K_{\text{tor}}\Delta t\,{\boldsymbol {w}}_{n}}$

曲げトルクは

${\displaystyle \Delta {\boldsymbol {\tau }}_{t}=-K_{\text{ben}}\Delta t\,{\boldsymbol {w}}_{t}}$

法線ダンピング力は

${\displaystyle {\boldsymbol {F}}_{{\text{damp}},n}=2\zeta {\sqrt {M_{e}K_{n}}}\,(-{\boldsymbol {v}}_{n})}$

接線ダンピング力は

${\displaystyle {\boldsymbol {F}}_{{\text{damp}},t}=2\zeta {\sqrt {M_{e}K_{t}}}\,(-{\boldsymbol {v}}_{tr})}$

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

${\displaystyle {\boldsymbol {\tau }}_{{\text{damp}},n}=2\zeta {\sqrt {J_{s}K_{\text{tor}}}}\,(-{\boldsymbol {w}}_{n})}$

曲げダンピングトルクは

${\displaystyle {\boldsymbol {\tau }}_{{\text{damp}},t}=2\zeta {\sqrt {J_{s}K_{\text{ben}}}}\,(-{\boldsymbol {w}}_{t})}$

Summarize Timeline Top Qs Fact Check

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

履歴ベクトルを

${\displaystyle {\boldsymbol {B}}^{(n)}=[{\boldsymbol {F}}_{n},\;{\boldsymbol {F}}_{t},\;{\boldsymbol {\tau }}_{n},\;{\boldsymbol {\tau }}_{t}]^{\mathrm {T} }\in \mathbb {R} ^{12}}$

とする（上付き ${\displaystyle (n)}$ は時刻 ${\displaystyle t_{n}}$ を示す）。

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

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

これは

${\displaystyle {\boldsymbol {B}}^{(n+1)}={\boldsymbol {B}}^{(n)}+\Delta {\boldsymbol {B}}}$

と書き換えられる。

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

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

${\displaystyle {\boldsymbol {F}}_{t}'={\boldsymbol {F}}_{t}^{(n)}-({\boldsymbol {F}}_{t}^{(n)}\cdot {\boldsymbol {e}}_{n})\,{\boldsymbol {e}}_{n}}$

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

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

${\displaystyle {\boldsymbol {\tau }}_{n}'=({\boldsymbol {\tau }}_{n}^{(n)}\cdot {\boldsymbol {e}}_{n})\,{\boldsymbol {e}}_{n}}$

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

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

${\displaystyle {\boldsymbol {\tau }}_{t}'={\boldsymbol {\tau }}_{t}^{(n)}-({\boldsymbol {\tau }}_{n}^{(n)}\cdot {\boldsymbol {e}}_{n})\,{\boldsymbol {e}}_{n}}$

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

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

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