前回は下記の論文に参考に転がり抵抗を導入してみました。今回は、要素を土にするためにもう一つ必要な、粒子間の吸着力（土の世界でいう粘着力）を考慮したモールクーロンの破壊基準の適用にチャレンジします。
本来は、転がり抵抗モデルが意図した通りに組めているかの検証をしなくてはいけないところですが、やりたい・やれそうなことは、順序は関係なくやっていくスタンスです＾＾
地盤材料の破壊基準を表現するためのシンプルな個別要素モデル
何をするかと言うと、せん断に対する抵抗力を一定値ではなく、拘束力（直方向の力）に係数を掛けた値を上限とするモールクーロンの破壊基準を組み込みます。加えて、通常のDEMでは扱わない引張抵抗力（粘着力）を考慮するようなモデルにします。この引張力が転がり抵抗にも関与するので、それに応じて、転がり抵抗部分も変更します。
とにかくやってみよー。
ボンドモデル部分ソース（全ソースは最後に）
# 粒子間ボンドモデル
cdef double bs,tu=0
if i < peCount and j < peCount: #粒子同士
tu = pe[i].et[j]
if hn < -tu and tu > 0: #ボンド破壊条件
pe[i].et[j] = 0
tu = 0
print('bond break',i,j,hn,tu)
#モールクーロン破壊基準の適用
bs = (hn + tu) * tan(pe[i].ph)
if hn <= -tu:
hs = 0.0
elif fabs(hs) > fabs(bs):
hs = bs * hs/fabs(hs)
#hn = hn - tu
else: #粒子同士以外
if hn <= 0.0:
#法線力がなければ、せん断力は０
hs = 0.0
elif fabs(hs) >= inf.frc*hn:
#摩擦力以上のせん断力は働かない
hs = inf.frc*fabs(hn) *hs/fabs(hs)
上記論文ではDF条件と称して、接触状態にある粒子間では粘着力が発揮されていて、一定の引張力が加わったときに粘着がなくなり（ボンドが破壊）、その後の接触は通常の摩擦として扱うモデルになっています。
普通は離れている要素の粘着力は既にボンド破壊状態であり、粘着力は働かない（通常の摩擦であらわす）のですが、確かめるための便宜上、初期値では粘着力を有効にしておいて、破壊されるまでの間は、接触する要素間に粘着力が働くようにしておきます。
さて、どうなるか？
（左：粘着なし 右：粘着有り）
お、粘着力は働いてるっぽい。しかし、2つ問題が発生。動きが止まってきたあたりにバネのようにスライドする謎の動きと、後半、なぜか要素が離れていく・・・。
前者は、今回のボンドモデルに関係なく、線と円要素の接触処理によるところの部分の不具合ですかね。そういえば、前の転がり抵抗のときもそうでした。摩擦の上限を超えることなく、せん断バネでスライドしているせいかな？検証が必要ですなー。パラメータが非適正な値になっているだけだといいんだけど。
後者について、実はイメージで言うと、今回の動きは、僕の想像と違って、もっと粒子同士が重なって、重なった状態で安定するイメージでした。上のソースだと、粒子が反発して離れる過程で、引張力（粘着力）が働くけども、発揮される直方向の力自体の値には影響を与えません。収束状態は力の釣り合ったところなので、ボンドが破壊するまで反発が働き続けると、結果、離れてしまうのかなと。
ともすれば、最終的な力は、弾性力（反発力）に引張力を差し引いた値にするのでは？と思い、上のソースの最後に
hn = hn - tu
を加えてみました。すると、
イメージ通り、粒子同士は結着したものの、転がり抵抗力に余計な影響を与えるらしく、結着したまま回り続ける謎の動きにorz
転がり抵抗をOFFってみると、こんな感じ。
ありゃー、結着しながら転がっていくー。
転がり抵抗のモーメント算定式を修正するのと他に、粒子同士が結着した接触点が固定されるように回転を考慮しなくはいけないことが判明。
上記論文、モデルについて丁寧に書かれてくれている部類の論文だと思いますが、いざコード化すると、色々と不明な点が出てきます。さて、どうしようかな・・・。