3Dグラフィック、プログラミング、ものづくり・・・、創作活動を楽しもう!

カテゴリー: プログラミング

DEMの勉強9

前回は下記の論文に参考に転がり抵抗を導入してみました。今回は、要素を土にするためにもう一つ必要な、粒子間の吸着力(土の世界でいう粘着力)を考慮したモールクーロンの破壊基準の適用にチャレンジします。

本来は、転がり抵抗モデルが意図した通りに組めているかの検証をしなくてはいけないところですが、やりたい・やれそうなことは、順序は関係なくやっていくスタンスです^^

 

地盤材料の破壊基準を表現するためのシンプルな個別要素モデル
http://www.agrifacility.kais.kyoto-u.ac.jp/thesis/pdf/71.pdf

 

何をするかと言うと、せん断に対する抵抗力を一定値ではなく、拘束力(直方向の力)に係数を掛けた値を上限とするモールクーロンの破壊基準を組み込みます。加えて、通常の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:
            #法線力がなければ、せん断力は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ってみると、こんな感じ。

ありゃー、結着しながら転がっていくー。

 

転がり抵抗のモーメント算定式を修正するのと他に、粒子同士が結着した接触点が固定されるように回転を考慮しなくはいけないことが判明。

 

上記論文、モデルについて丁寧に書かれてくれている部類の論文だと思いますが、いざコード化すると、色々と不明な点が出てきます。さて、どうしようかな・・・。

(さらに…)

DEMの勉強8

いよいよモデルいじり。粒子を土のモデルにしていきます。
モデルは色々あれど、一番最初に検索したら出てきたこの論文を参考にすることにします。

 

地盤材料の破壊基準を表現するためのシンプルな個別要素モデル
http://www.agrifacility.kais.kyoto-u.ac.jp/thesis/pdf/71.pdf

 

ポイントは2つ。モールクーロンの破壊基準を適用することと、転がり抵抗を導入すること。この内、転がり抵抗を導入してみました。

 

転がり抵抗は、実際の粒子は角ばっていて、コロコロ転がる事はないのですが、円形要素では、コロコロ転がっていってしまうので、それを防ぐために導入します。要は、円形要素で角ばった砂の粒子を模擬するためのものです。

さて、ちゃんと組めているか?


(左:転がり抵抗無し 右:転がり抵抗有り)

 

一瞬、戻っているような・・・気がするけど、コードがおかしいのか、せん断バネの反発で変に見えるだけか。とりあえず、抵抗力は発揮されているみたい。

さあ、沢山粒子のあるパターンで、安息角に違いが出るか確かめたい。

逆に、流れ出す粒子の数が増えてしまったなり。はて。イメージと違う。

 

ちなみに完全に回転を止めてみると・・・、

いくつも角が生えて、非現実的な面白い感じになります。

 

さて、どうやってちゃんとできているか確かめよう・・。

(さらに…)

DEMの勉強7

前回は、下記の論文で言うところのセル分割法を実装してみましたが、二つ目の手法、粒子登録法を組んでみました。

 

GPUによる近接相互作用に基づく粒子計算の近傍探索手法
https://ipsj.ixsq.nii.ac.jp/ej/?action=repository_uri&item_id=107323&file_id=1&file_no=1

 

出力結果は同じ(動画省略)。さぁ、速度はどうだ。

全粒子判定 1ステップあたり0.039175秒
セル分割法 1ステップあたり0.023670秒(1.65倍)
セル分割法+粒子登録法 0.022526秒(1.73倍)

 

・・・速くなり・・・ました。ちょっぴしだけど。

 

論文では、セル分割法単体と比べ、3倍以上の高速化となってますが、今回はそうはなりませんでした。

コードが良くない可能性は大いにあるけども、一応、セル分割法の候補が粒子登録法によって、1/3〜半分程度に抑えられることは確認。ただし、更新頻度が多いので、そっちに処理を食っている模様。試しに更新を判定するアルゴリズムはそのままに、粒子登録法の絞り込みをコメントアウトしてみると・・・

 

セル分割法+粒子登録法 0.022526(1.73倍)
セル分割法+更新判定のみ 0.022452(1.74倍)

変わらーん。むしろ絞り込みしない方が速いー。っとの結果になりましたorz

 

粒子登録法は、コードがごちゃごちゃするので、いっそ無くしてしまうか・・・?まぁでも、更新頻度が少ない解析なら、効果を発揮するやもしれぬ。って事で、アルゴリズム的な改良は一旦終了させて、モデルの検討をやっていきたいと思います。

 

(さらに…)

DEMの勉強6

やりたいことは、色々なモデルを試すこと。プログラムは遅くたって最低限動けばよし!次はモデルをいじるぞ!

 

っと思っていたのですが、プログラムを眺めてると、さすがに全て粒子で、全ての粒子に対する当たり判定するのは、無駄だよな~っと気になりだし、何となくアルゴリズムを考えるようになったところに、下記の論文を見つけました。

 

GPUによる近接相互作用に基づく粒子計算の近傍探索手法
https://ipsj.ixsq.nii.ac.jp/ej/?action=repository_uri&item_id=107323&file_id=1&file_no=1

 

アルゴリズムの中身がけっこう丁寧に書いてある。これなら実装できそうだ。っと思ってしまったからには・・・やりますか。

 

書かれてる手法は、大きく2つ。

一つは領域を格子に分割して、格子内の粒子を登録しておき、自身がある格子とその周辺の格子内の粒子に対してのみに当たり判定を行う方法。上の論文では、セル分割法と名付けられている。オーソドックスな方法みたいで、他の文献では、格子判定法だとか色々なネーミングがありました。

もう一つは、一度は全粒子に対して当たり判定をして、その際、近くの粒子を登録しておき、一定期間は、登録した粒子に対してのみに当たり判定を行う方法。上の論文では、粒子登録法とネーミングされています。

 

どちらも一定条件の元、情報を更新しなくてはいけない。セル分割法は粒子の数だけの一回ループで情報更新できるのに対して、粒子登録法は粒子数の2乗ループが必要になるので時間がかかる。しかし、粒子登録法はセル分割法に比べ、登録する候補粒子を少なくすることができる。また、論文では書かれていないけど、わざわざ全体の情報を更新しなくても、動きのある粒子に対してのみ情報を更新させることもできるはず。

メモリ消費量の違いもあるけど、今のところメモリは問題ない(多分、2次元なら問題にならない)ので考えません。

 

論文を見る前に、自分で思いついたのは、粒子登録法の方です。通常、土の解析だと激しく動く箇所は一部だから、相対移動量とか拾って、上手いこと更新範囲と回数を抑えるアルゴリズムにすれば、いい感じになるんじゃね?っと思ったのだけど、その更新判定の部分を作るの大変そう。
論文では、セル分割法と粒子登録法の両方を組み合わせたハイブリッドなやり方が書かれてあり、考えるの面倒だから、その通りやってみることにしよう。

 

ってことで、まずセル分割法を導入してみました。とりあえず、登録情報は毎回更新する仕様。

よしよし、結果は同じだな。さて、計算速度は?

(さらに…)

TibLab Blog © 2014 Frontier Theme