最近色々あって、何かこう、普段全く使わない脳味噌の部位をスパークさせているもんで、流石にストレスが溜まってきたのでちょっと解放もされてきたし、まぁ、久々にこんな事でも書いてみようかと思って書かんとするものなり。
で、本稿はある程度、DEMとかのイメージがついてる人向けの文章のつもりで、
その辺もよくわからん、という場合は前に書いた何か
https://qiita.com/NatsukiLab/items/476e00fea40b86ece31f
をご参照してくださいませ。
個別要素法
個別要素法(離散要素法)みたいに、(大雑把に)固体同士が衝突しあって動くような計算をしたい時に、2つの剛体球間の反発力をどう実装するか、というのは種々マナーがある、ような気がします。
何となく個人的にはヘルツ接触(Hertzian contact)みたいなのが好きなんだけれど、実は書いたことが無かったり。
なので、まぁ、せっかくなのでFDPSにでも載っけて書いてみようかなと思ったメモ書きです、はい。
(最近のQiitaってあまりプログラムしてない記事は消されるらしいけど、しばらく数式しか存在しないこの記事は平気かしら…)
用語
さて、まぁまずは話をする前に、用語の導入をしておくわけですが、なんでこんな事をいきなりしだすかというと、少なくとも僕にとっては紛らわしかったからです。
で、何はともあれ用語の導入を…
弾性(elasticity):何かあらぬ事が起きた時、それを解消しようとして働く力
粘性(viscosity):(↑で誘起された運動を)減衰させようと働く力
です。
ってどういうこっちゃ？まぁまぁ、慌てず慌てず。次で説明しますから。
DEMの基本
前と同じように、各粒子の質量を$m$、半径を$r$、位置座標ベクトルを$\vec{x}$、速度座標ベクトルを$\vec{v}$としておきます。
今回はそれに加えて、慣性モーメント$I$と角速度ベクトル$\vec{\omega}$を定義しておきます。
んで、基本的には粒子$i$にかかる力とトルクを$\vec{F}$、$\vec{T}$と定義しておくと、
m_i \frac{\mathrm{d} \vec{v}_i}{\mathrm{d} t} = \sum_j \vec{F}_{ij}\\
I_i \frac{\mathrm{d} \vec{\omega}_i}{\mathrm{d} t} = \sum_j \vec{T}_{ij}\\
となるわけです。
で、$ij$粒子ペアに働く力とトルクが果たしてなんぞや、というのが今回のテーマ。
で、$\vec{F}$と$\vec{T}$が果たして何から成るかというと、
接触面に大して垂直に働き、接触を解決しようとする垂直抗力(normal force/torque)成分と、
接線方向に働く成分(tangential force/torque)からなると思う事にします。
なので、
m_i \frac{\mathrm{d} \vec{v}_i}{\mathrm{d} t} = \sum_j (\vec{F}^\mathrm{n}_{ij} + \vec{F}^\mathrm{t}_{ij})\\
I_i \frac{\mathrm{d} \vec{\omega}_i}{\mathrm{d} t} = \sum_j (\vec{T}^\mathrm{n}_{ij} + \vec{T}^\mathrm{t}_{ij})\\
と書けます。
というわけで、右辺の4つの項がわかればまぁ問題なさ気。
normal force
前回と同じように貫通深度(penetration depth)を以下のように定義しましょ。
d_{ij} = \max\left(r_i + r_j - |\vec{x}_j - \vec{x}_i|, 0 \right)
normal forceをどう実装するかですが、とりあえずペナルティ法を考えます。
ペナルティ法とは、「何かあらぬ事が起きた時、それを解消しようとして働く力(ペナルティ)」を加えるという方法です。
似たような文面をちょっと前に見たような？
さて、それでどんな力を加えるわけですが、基本的には貫通深度(のべき)に比例して大きくなるようなものにしておけばよさそうです。
これは、粒子の重なりの部分にちょうどバネを差し込んで、縮んだ分押し戻すような感覚で働くので、比喩的にスプリングとか呼ばれます。
一方で、バネで跳ね返すだけだと、粒子が密集してきた時に、延々と粒子が振動しはじめるとか、そういう困った事もあるので、適当に「減衰させようと働く力」を突っ込んでダンプしておく事にします。
これは、単純には速度に比例して大きくなるような抵抗力を突っ込めばいいわけですが、比喩的にダッシュポットとか呼ばれたり無かったり。
さて、まぁともかくそういうわけで、若干天下り的に、
\vec{F}^\mathrm{n}_{ij} = - \frac{4}{3} E_{ij} \sqrt{r_{ij}}d_{ij}^{2/3}\vec{n}_{ij} - c^\mathrm{n} \sqrt{6 m_{ij} E_{ij} \sqrt{r_{ij} d_{ij}}} \left[(\vec{v}_j - \vec{v}_i) \cdot \vec{n}_{ij}\right] \vec{n}_{ij}\\
\vec{n}_{ij} = \frac{\vec{r}_{j} - \vec{r}_{i}}{|\vec{r}_{j} - \vec{r}_{i}|}\\
r_{ij} = \frac{1}{\frac{1}{r_i} + \frac{1}{r_j}}\\
m_{ij} = \frac{1}{\frac{1}{m_i} + \frac{1}{m_j}}\\
E_{ij} = \frac{1}{\frac{1 - \nu^2_i}{E_i} + \frac{1 - \nu_j^2}{E_j}}
という事にしておきます。この力をヘルツ接触力と呼びます。
んで、いきなり$E$と$\nu$とかいう未定義の量が出てきたわけですが、これは一声、「硬さ的な量」であると申し上げておきます(正確にはヤング率とポアソン比)、詳しいことはまたのちほど(「のち」があるかわからないけど)。
んで、$c^\mathrm{n}$は、悲しいかな経験的なパラメタです。
大体$10^{-5}$程度だとかなんとか。
tangential force
これから書く…