観測孔の遮水

以下の文献を読んでいました。

白石「地下水調査と地下水流動解析の高度化に関する研究」

http://ousar.lib.okayama-u.ac.jp/ja/50813

簡単な数値実験をもとに、水位観測データの解釈をされています。

特に興味を惹かれる内容はなかったのですが、1箇所だけ脅された文章がありました。

「5.1.1　地下水位観測の高度化」にて、単孔式の地下水位観測孔の注意点が示されています。その中で、背面遮水は「全層遮水が望ましい」とあります。塩ビと孔壁間の遮水を全区間行った方が水位への影響が少ないといった内容です。根拠は 5m 程度の観測孔をモデル化した数値実験でした。

最初は驚きましたね。今まで気付かないうちに失敗していたのかな、と。

よく読むと、遮水部の長さに偏りがあります。計算ケースは0.1、0.6、1.1、4.9m の4種。で、1.1mでは10cm以上の水頭差があるので、全区間遮水（4.9m）がお勧めという結果。なぜ、1.5、 2.0、2..5 と続けて刻まないのでしょうか？
しかも、水頭は地表以上になっています。つまり、ストレーナー上部で遮水をいくらしようが、地表付近を遮水しない限り観測値に影響が出るモデルなのでしょう。自噴する井戸の遮水を考える場合にどうするか、ということを前提とされていたのかもしれません。脅かされました。

地盤工学会の基準では遮水区間 50cm 以上です。個人的には施工時の安全率をかけて 1～2m 程度設けています。ま、先の結果を見る限り、1～2mで十分なように思えます。

砂利やペレットを現場で綺麗に充填するのは、結構難しい作業です。そのあたりの検討や経験・ノウハウの公開の方が、地下水調査の高度化に直結すると考えます。
単純な解析も必要ですが、それだけでは高度化は望めないでしょう。

VR / AR

後輩より相談

「VR ヘッドセットを購入するが、何か使い道ありませんか？」

後輩はシミュレーション結果を表示するために、展示会などでヘッドセットを使っています。コンテンツが大がかりですので1式レンタルしていたのですが、そろそろ購入するとのこと。ただ、用途がそれだけだともったいないので、他に何か使い道はないか？という意図でした。

ないですね。

今まで、住民説明用に Infraworks でモデルを作成したことがあります。思いつくのは、これをVRデータとして取り込む程度です。が、これは本来は設計者のプレゼンの場でしょうし、3次元設計が前提です。音も不要です。

後輩が言うには、Navis からでも VRデータを吐き出せるとのこと。

具体的には以下のリンクにあるように、データをクラウドでレンダリングしてパノラマ表示にしているだけのようです。FBX 経由で Infraworks のデータも、対応可能となっています。Navis 系のテクスチャ―ですので、3ds Max 経由だと飛んでしまうでしょうね。

https://www.youtube.com/watch?v=M4pCiv_IdBg

残念ながら定点のみで、BIMx のようなウォークスルーはできません。やはり、BIM は CIM より一歩進んでいます。

私の仕事以外では、重機の遠隔操作に利用できそうです。無人重機の操作において、備え付けられたカメラを見るより、ヘッドセットの方が見たい方向をすぐに見ることが可能ですし、音も聞けて操作しやすいと思います（30分ほどで酔うかもしれませんが）。私が知らないだけで、既に開発済みかもしれません（と思って調べてみますと、ありました）。
http://www.nikkei.com/article/DGXLZO11989150S7A120C1TJM000/
基本、ゲームが出発点ですので、こういった利用法であれば現段階でもいくつかあるでしょう。そういえば、塗装の訓練に利用しているニュースも見たことがありますね。

将来的には VR でなく、CAD 段階の計画構造物を AR データとして利用することが主流になるでしょう。現位置にデータを綺麗にはめ込むにはコツが必要でしょうが、 それも GPS とマーカーの併用、アルゴリズムの改善等で解決されると思います。

ま、3次元設計、測量、調査が先ですね。VR へ流用する技術や環境は既に整っていますが、業界のヒトにとっては、まだまだ遠いかもしれません。

コールドロンと火砕岩岩脈

後輩から「岩石名が分からない」と、写真が送られてきました。

コアの数10cm区間が、破砕状になっています。その部分の名前が分からないとのことでした。
固結した破砕部ですので、岩石名を付ける必要はないとは思いましたが、成因を知りたい=岩石名が決まるということのようです。
よく見ると、破砕岩中に異質の岩石が取り込まれていました。破砕岩の周りは結晶片岩なのですが、破砕岩中には深成岩片が取り込まれています。
１．ゼノリスのように深部の母岩が取り込まれた岩脈か、2.上部の礫が落ち込んだ割れ目区間かしか思いつかなかったのですが、周辺の地質を知りませんので判断できませんでした。

コアが分析者の手元に送られてきました。それを見た限り、特に上記のどちらかで問題はないようでした。が、新人さん曰く「タフダイクではないか？」と。
タフダイクってなんだ？と思って聞いてみますと、コールドロンの形成にかかわる文献によく出てくる名称とのこと。
調べてみますと、違うようでした。液状化のサンドパイプのような成因に対しタフダイクという言葉が使われています。コールドロンでは「火砕岩岩脈」「貫入性角礫岩」でした。これらを言いたかったのでしょう。

と言っても、火砕岩岩脈という言葉にも違和感があります。火砕岩は地表での堆積起源であり、深部から圧入した「岩脈」とは成因が異なるイメージを持っています。
コールドロンの文献では、火山活動に伴う流動物質で母岩片を取り込むものや、火砕物を供給する火道跡にて固結した岩石（およびその後のドレインバックによる圧縮・溶結、せん断）などといった成因を指して命名しているようです。火道に火砕物を供給する、という表現にも引っかかりますが、ま、現物を見ないと何とも言えません（見ても納得できないかもしれませんが）。

コールドロンについては、近年、複数個所でその可能性が示されているようです。それらの文献を読み進めながら火山岩の分布域を見てみると、コールドロンのない地質図の方に疑問が生じると思います。今後、コールドロンといった現象を重力探査で裏付けされたものが広域の地質図に反映されだすと、それに伴う地質（特に岩脈や破砕部）の成因について考慮する要因が増えることになるでしょう。土木上も、その成因がトンネルなどの計画に影響するのかどうかを見極める必要が出てきます。
広域の特徴を押さえた上で、ローカルの議論を進める視点がより重要になるということです。

油圧ショベルで試料採取

今日は同僚の現場で CBR 試料採取のお手伝い。

ミニユンボを積んで回送し、現場で降ろして掘削です。

油圧ショベルを扱うのは2年ぶりでした。特別教育以来のペーパードライバーです。その分、作業は身の丈に合うよう慎重に、安全に。

自分で掘ってみると、土質の変化や硬軟、水の付き具合の変わり目がよくわかります。それが面白く、また全く体力を使わないこともあり、集中して掘っていたように思います。ボーリングのオペさんもそのような感じなのかもしれません。

掘るのは楽しいのですが、段差にブリッジをかけて降りるのが怖かったですね。「慣れ」と言われますが、2年に1回では慣れません。気疲れしました。重機オペさん尊敬です。

今日の作業は順調にいき、3箇所の試料採取が終わりました。配合試験をするため1箇所につき10袋ほど取りました。約800kgです。掘るのは楽だったのですが、この積み下ろしで腰に来ました。

明日は残り2箇所。
事故を起こさないように、慎重に作業しましょう。

*******************************************

20170318追記

何度かブリッジで乗り降りし、コツがつかめた頃に作業終了。次に使う頃には、また振り出しに戻りそうです。

採水

ひと段落したら、地下水の年代測定をおこなうため、採水に出かけようと考えています。

結果の整理や、それを使った検討を行ったことはありますが、年代測定用の採水に出かけるのは初めてです（所謂、使えないヤツです）。

今回は溶存ガスを測定したいので、採水の際に大気に触れさせることができません。どうやって採ろうかと悩んでいます。
いえ、湧水なら容易なのですが、ボーリング孔の深い位置からの採水が悩み処です。孔内ポンプは羽の部分でエアをかむのでダメ。ベーラーでもプロなりの工夫はあるようですが、幾つかのステップを要するため一人での実施は困難。
他支店の経験者に聞いて見ましたが、「ボーリング孔からの採水は難しいので計画しない」と言われていました。

で、ネットを見ていると、このような文献が。
https://www.researchgate.net/publication/277573322_Groundwater_Ages_of_Confined_Aquifer_in_Mishima_Lava_Flow_Shizuoka
止水弁（ベーラーの玉？）が上下に2個ついて、その間にコックが2個ついているタイプの様ですね。見た目はこのような形でしょうか？
http://www.panasiatec.com/sample_cylinder.html
なかなか良いアイデアで欲しくなるのですが、採水箇所分のサンプラーを準備する必要があります。また、オールサスで申し分ないのですが、その分お高いでしょうから数をそろえるのは現実的ではなさそうです。

で、思い出したのがコチラ↓
https://youtu.be/pd-wG0ej5b0
後輩が買ってすぐに他支店に貸して、そのままでした。ゆっくり上下させれば気泡の心配はないでしょう。これだと一人で採水できます（何かしら、人力以外の動力を組み込めないでしょうか？）。

採水ひとつでも、わからないことや教わらないとできないことがあるようです。この歳ですが、まだまだ経験させてもらいましょう。

岩石薄片図鑑

本屋さんで綺麗な図書を見かけました。

青木正博「岩石薄片図鑑」

薄片写真が大きく掲載された図書です。見開きで掲載されているものもありました。図鑑というほど多くの種類は載っていませんが、蛇紋岩化初期の橄欖岩や（一般的な）スカルン鉱物が載っていて興味をひかれました。

デジカメのマクロ機能で撮れるのだそうです。カメラには詳しくないので説明を眺めても、どうやって撮っているのかわからなかったのですが、その手法で移された大きな写真に魅入ってました。先日の顕微鏡写真も驚きましたが、偏光顕微鏡写真を対物レンズなし？で撮れるのにも驚かされます。低倍率ですが広範囲の写真が撮れるため、非常に見やすくインパクトがあります。

購入しましょう、と思ってしばらくの間手に持っていたのですが、そういえば「薄片でよくわかる 岩石図鑑」も買ってなかった、あれも買ってなかった、これも、と思い初め、最終的には中断していた「Fault-Related Rocks」を読み切ってからにしましょう、と思い直し書棚へ戻しました。

読まないといけない文献もたまり始めています。
早く読まなくては。

重心深度と見かけ比抵抗

表皮深さの1/2を空中電磁探査の深度として利用している文献がありました。

長谷川ほか「ヘリコプターを用いた空中物理探査データの再解析」
http://jolissrch-inter.tokai-sc.jaea.go.jp/search/servlet/search?5042189

ここからたどり着いたのが、こちら。
B. Siemon (2001) Improved and new resistivity-depth profiles for helicopter electromagnetic data
http://www.sciencedirect.com/science/article/pii/S0926985100000409

あるモデルでは表皮深さの1/2を深度として利用する手法がよくあいました、といったような内容です（あわないモデルもあるようです）。あう、あわないといっても、机上の話です。現地で確認しました、電気検層と対比しました、といったような内容ではなく、説得力を得られません。
1/2の理由はさらに先があり、その文献を取り寄せても、さらにまだ先がありました。ただ、この文献ではあわせこみパラメータとも取れる内容でしたので、深追いはしていません。


空中電磁探査での見かけ比抵抗等の求め方を詳しく知らなかったのですが、複数の方法があるようです。そのうち、表皮深さを用いた方法は Siemon による以下の2通りがメジャーなようでした。それ以外にも、√2で割る方法がありました。

1-1. R,Q→Da→da→Zsim

A = √(R^2+Q^2)
γ ＝ s/h
A'^（1/3)=A^（1/3)/γ
Da = s(A'/A)^（1/3)
Zsim = da + p/2

s: sensor horizontal distance
h: sensor height
p: skin depth

1-2. h,p→ρsim


2-1. calculating phase ratioε

ε=|Q/R|

R （real） ppm
Q (quadrant) ppm

2-2. log（ε）→log（A'^（1/3)）→Da→da→Zsim

2-3. log（ε）→log（δ）→ρsim

R,Q,h を 測定し、ρsim と Zsim を求める際に p （と係数）をパラメータとして導入しているように見えます。たとえ1/2や表皮深さの導入に理論的背景があったとしても、測定対象と整合しなければ他の方法を選択できるというのは、理論がまだ確立されていない、推定するのは難しいということなのでしょう。
それに、以前にも書き残していますが、実際に磁場が到達しているのかどうかは確認しないとわかりません。また、1つ目の文献のように、得られた信号がノイズレベル以下になっていないかどうかもチェックが必要でしょう。

そうなると、現段階では電気検層など他の手法と空中電磁探査結果の対比、その補正が必須となります（物理探査ですので、当然なのですが）。見た目の綺麗な断面図に惑わされがちですが、別の物性値と対比がなされているか、推定式の妥当性は確認されているか？といったステップを踏まないといけないのでしょう。

DualSPHysics future

今後のリリース予定がマニュアルに書かれていました。

16. DualSPHysics future
An update to DualSPHysics will be released as v4.2. The updates planned for this release include:
 [Adami et al., 2012] boundary condition
 Local Uniform Stencil (LUST) boundary conditions [Fourtakas et al., 2014].
 MultiGPU with OpenMP.
The new version DualSPHysics v5.0 that will be released in the future will include:
 Multi-GPU implementation [Dominguez et al., 2013b].
 Variable particle resolution [Vacondio et al., 2013; 2015].
 Multiphase (air-water) [Mokos et al., 2015].
Other features that are planned to be integrated into the DualSPHysics solver that are now under development:
 Wave absorption system [Altomare et al., 2015a].
 Incompressible SPH (ISPH)
 Inlet/outlet flow conditions.
 Coupling with SWASH Wave Propagation Model [Altomare et al., 2015b].
 Moorings [Barreiro, 2015].
 Coupling with MoorDyn library (http://www.matt-hall.ca/software/moordyn/).
 Coupling with Chrono-Engine library (http://chronoengine.info/chronoengine/).

流入境界については、まだ先の話のようです。残念。
ネット情報では Blender での可視化プラグインも開発中とのことでしたが、ここには載っていません。これは早く欲しいですね。

DualSPHysics、流体に関しては威力を発揮するソフトだという印象を持ちました（土砂を扱えなかったのはとても残念ですが）。
メッシュレスのため、計算に取り掛かるまでの時間が FEM 等に比べ圧倒的に短いというのは、粒子法の利点でしょう。メッシュを切るのに時間をかけるのがもったいなく感じます。

まだまだ実務で使うには足りないところもありますが、継続して開発が進められていますので、今後に期待しましょう。

*******************************************************
20170911追記
Blender add-on 出ました。
VisualSPHysics
http://visual.sphysics.org/

堰堤が受ける力（DualSPHysics Ver.4.0）

砂防えん堤のCADデータ取り込みは問題なし。

えん堤上流側に圧力の観測点を鉛直方向に8つ配置。
それとは別に、鉛直壁面を別途作成し、そこにかかる力の合計を得ることにしました。
ポイントの座標を記載したファイルを用意し、バッチファイルを編集して、計算！

流速と流体力に関しては問題なく抜き出せました。
観測点での圧力は、XYZ成分への分け方がわかりませんでした。forum を見る限り、方向（プラスマイナス）を加味して抜けるようですが、説明書には書かれていません。

えん堤にかかる流体力の合計は以下の通り。

Z方向がマイナスなのは、土砂が流下して下向きの力がかかっているからでしょうか？観測点から抜いた流速のz成分は全てプラス方向なので矛盾しているように見えます。
流速がマイナス方向なのは、X成分のみでしたが、力の合計はX成分がプラス。Y成分は両方プラス。観測点が少ないのが原因でしょうか？よくわかりません。

観測点にかかる圧力です。おおきな値です（入力単位を間違えたかな？と確認しましたが、あっていました）。
最初、堰堤にぶつかり圧力が上昇し、越流しながらやや降下しているように見えます。
一番下の観測点はほぼ0。これは、地表に近すぎたのが原因です。地表からも1.5h 離さないダメでした。2つ目も、影響を受けているようです。
最初、マイナス値は引張り方向かと思いましたが、よく考えてみると軸の負の方向をあらわしているだけでしょう。←これは違うようです。

そもそも、内部の応力なんて取り扱っていないと思われます。あくまで得られるのは、堰堤にとって外力（静水圧、流体力、衝撃力）のみです。堰堤内部の応力を見たいのであれば、別途、弾塑性を扱えるソフトで、この結果を時刻歴として入力し、計算するしかないでしょう（ここまで来てようやく気が付きました）。それには、この部分だけ解像度を上げて計算したいですが、粒子間隔は部分的に変更できません。部分的に切り出して粒子間隔を小さくして計算するにしても、境界条件が対応できません。あくまで流体の挙動を見るソフトですので、この辺は適用外と考えましょう。

深層崩壊を題材にし、そのモデル化と力・圧力・流速を抜く、といった一通りの手順を確認しました。残念ながら、圧力等の解釈、特にプラス・マイナスについて、まだ理解ができていません（ポスト処理に関するツールにも、解説を付けて欲しいところです）。
dynamic boundary 特有の性質も理解を妨げている一因でしょう。この辺は次期リリースで新しい境界条件を導入するとアナウンスが出ていますので、もう少し簡単になればありがたいです。

*********************************************
20170308追記

流速と圧力について、2次元モデルでチェックしてみました。

導水管を配置し、管の中で上昇流しか起こらない個所を作成。ついでに、噴水のように噴出した水流が池の中に落ちてたまる箇所を作成。その2箇所に観測点を配置し、水圧と流速をはかりました。

まずは管内の上昇流。これは2次元ですので流速のY成分は０。X方向もほぼ0。Z方向はプラス。OKです。圧力はZ方向で波打つように変化し、時々マイナスになっていました。水中でマイナスになる理由が分かりません。壁が近すぎるのでしょうか？

たまった水の中に落ちてくる方は、流速のZ成分がマイナス。これはOKです。
一方、圧力はプラス。これは水圧がかかっていますのでプラスで問題ないでしょう。観測点に関しては、流れの方向は関係ないと思われます。

CADデータの取り込み（DualSPHysics Ver.4.0）

CADで作成した地形の取り込みは、あっさりクリアーできました。

XMLファイルに関する説明書を読むと、以下の流れが記されていました（以前のVer.にも添付されていたのですが、見落としていました）。

1. STL、VTK ファイル等の閉じた形状を boundary、fluid として読み込む。
2. fluid の中の座標を1点指定。
3. 粒子を発生させる範囲（BOXの基点座標と、大きさ）を設定。

基本はこれだけ。3の範囲を fluidの範囲 よりも大きくしておくと、その中が粒子で満たされます。
これで3次元の計算準備が整うのですから、メッシュレスの威力、半端ないです。


今回の題材は、以前、土砂移動シミュで作成していた深層崩壊。
範囲は3km2。
地形とその上に分布する土砂の2つのソリッドを Civil3D から書き出し、2つのSTLファイルを作成します。それらを input ファイルで boundary、fluid として指定するだけで読み込むことができました。

少し引っかかったのが、マイナス座標を扱えなかった点。マイナスの付いた公共座標のまま書き出すと、fluid、boundary が自動で (X,Y) = (0,0) を原点として持つように移動されます（当たり前なのかもしれませんが、知りませんでした）。複数の構造物は、それぞれ原点に移動されますので相対位置がズレます。
これについては、モデル全体の基点をプラス領域に移動してから、各々をSTL 書き出しすれば、OKです。

で、計算実行！

粒子が1億個以上発生し、GPUのメモリーオーバー！

ま、予想されていたことです。

GeForce GTX480 だったのですが、GTX660 を SLI で2基連結したマシンではどうか？と思い試してみました。が、なぜかこちらもシングルでしか動かず、メモリーオーバー。

で、粒子間隔を変更しつつ、動作するレベルを探っていきました。
結果、数百万粒子までは計算できそうです。粒子間隔が2mだと、GTX480 で3日となりました（なぜか、GTX660では4日。クロックでしょうか？）。

テストケースとして粒子間隔 4mで実施。これだと 8７万粒子（実際に動く粒子は17万個）、120秒間の計算で約5時間程度。

結果を確認すると、おかしな動きをしています。発散したようです。


今度は粘性を上げて計算してみました。

OK。うまくいきました。
と言っても、斜面の上から水を流すように流れるだけで、土砂の動きではありません（この点は、Multi-Phase で承知しています）。
粒子の動きが落ち着くと、天然ダムを形成せず、谷を長く埋める湖の様なフラットな形状ができあがりました。粘性をさらに上げると、それらしい形状になりましたが、おそらく時間の問題でしょう。流体ですね。ま、現段階では仕方ないと思います。
また、これだけ間隔をあけても良いのか？と疑問でしたが、問題ないようです。ま、見た目違和感ない粒子配置でした。

とりあえず、CADデータの取り込みや粒子配置は容易、ということが分かりました。
今度は砂防ダムを作って、それにかかる圧力を見てみましょう。

Dam Break （DualSPHysics Ver.4.0） その2

圧力の抜き出し方に関しては、説明書に記載されていました。

面白いことに、壁面での計測はダメで、スムージング長の1.5倍の位置で計測するのがお奨めとのこと。dynamic boundaryを使っているから、という理由らしいのですが、これ、理解していませんでした。

あらためて説明書の dynamic boundary の箇所を読み直すと、以下のように書かれています。

When a fluid particle approaches a boundary and the distance between its particles and the fluid particle becomes smaller than twice the smoothing length (h), the density of the affected boundary particles increases, resulting in a pressure increase. In turn this results in a repulsive force being exerted on the fluid particle due to the pressure term in the momentum equation.

壁面からスムージング長の2倍以内に流体粒子が入ってくると、壁面の（粒子？）密度を増加させ反発力をもって対抗する、壁面位置では応力が±0となり、壁面粒子の位置は動かない、という造りのようです。そうすると、壁面での圧力の抜き出しを避けなければならないというのは納得です。
等値面を表示させた際に、壁面近傍で液面が窪んだ形状になるのが気になっていたのですが、これもこの境界条件によるものなのでしょう（FAQにギャップの話が出ていますが、この現象のことかと思われます）。

検算してみますと、サンプルファイルのスムージング長hは

h＝cofh√3×dp

　=1×√3×0.085

　≒0.01472

壁が x=0.9 の位置でセットされていますので、

0.9-1.5×0.01472=0.87792

X=0.8779 の位置に観測点を置けば OK ということです。あっていますね。

結果も説明書通りに抜き出せました。kgで入力しているのに、Paで書き出されるのは、書き出し時に変換しているのでしょうか？入力単位が分かりにくいですね、このソフト。

input ファイル内で、壁面のみ別の marker value  (mk 番号)を持たせて作っておけば、その値を指定するだけで壁面にかかる力を合計してcsv にしてくれます（1.5hを考慮してくれているのでしょうか？）正しく抜けているとすれば、上記の機能とあわせ、実務で大いに使えると思います。

とりあえず、粘性の入れ方による安定性の違い、圧力の抜き出し方法については理解できたように思います。

Dam Break （DualSPHysics Ver.4.0）

次に手を付けたのは、よくあるダムブレイク問題。

偶然、同じモデルを FEM で解いていたツワモノがいました。どうやって解いているのか、最初は全く想像がつかなかったのですが、話を聞くうちに見えてきました。凄いですね。

同モデルの実験値が載っている文献を頂いたのですが、著者はSPHysics 関連の文献を発表されていた方でした。狭いもので、戻ってきました。

サンプルケースは問題なく動きました。

結果を動画にして FEM の方と2人で見ていたのですが、2人とも同じ感想。

「ネバい」

1.5秒ほどでほぼ動きが止まります。粘性が強すぎるのでしょう。

inputを見てみますと、人工粘性が選択されていました。

これを水の動粘性係数になるよう変更し、再計算！

結果は全くダメ。

発散し、全粒子が空中に飛散します。これはSPH特有の問題なのか、粒子法共通なのか、それともアルゴリズムに依存するのか？わかりません。容易には受け入れ難い結果です。人工粘性の方が安定する、柱にかかる圧力を実験値と一致させやすい、などの理由でこのセッティングなのかもしれません。

続きます。

DualSPHysics Multi-Phase その2

手始めに、2次元でモデル化。

題材は越流破堤だったのですが、inputファイルを作成中に、以下の制限に気づきました。

• 水位固定や、定流量の様な境界条件を設けることができません。
  今回は横にタンクを付けて代替えとしておきましょう。
• 現段階で2相以上の取り扱いはできないようです。それは良いのですが、基礎地盤と堤体の構成材料を区分した場合に3相と認識されるようでエラーがかかります。
  仕方ないので土は1層にしました。

いきなりですが、これら2点で実務への適用性はグッと低くなってしまいました。
ま、何ができるかの確認も必要ですので、先へ進めることに。

土砂のパラメーターを c=100kN/m2、φ=30°、γ=18kN/m3 とし、Drucker-Prager の破壊基準で input ファイルを作成し実行！

結果はダメ。

堤防が流体のように自重でグニャっと変形してします。そういえば、内部消費などは構成式に含まれていませんでした。粘性を高めることで対応する考え方でしょうか？その点、同じNS方程式を出発点とした LSFLOW はよく考えられていると気付かされます。
そもそも、自立する材料としての視点ではなく、土粒子としての視点で N-S方程式を適用しているのでしょうから、LSFLOW のように式を加工しないとうまく再現できないのでしょう。河床変動の様なモデルに最適なのかもしれません。

文献値に示されている通りに土砂の粘性を高めてみてもダメでした。
http://www.sciencedirect.com/science/article/pii/S0309170816300926

The fluid dynamic viscosity was 0.001Pa. s and sediment viscosity was set to 150Pa. s with the HBP m and n parameter set to 100 and 1.8 respectively. The value for the exponential growth m parameter value was chosen to approximate a Bingham model as closely as possible with a minimal pseudo-Newtonian region and the power-law exponent n to resemble shear thinning materials as shown in Fig. 3.3. Finally a small amount of cohesion was given to the sediment phase of c = 100 Pa to stabilise the interface and control the scouring near the dam gate. 

この粘性も、dynamic で入力するのか kinematic なのかは示されていません（通常版では kinematic と示されています）。さらに水の重量が ｋｇ で入力されていますので、c も kg/m2 で入力すべきなのでしょう。が、そのあたりも書かれていません。コメントや説明書に単位ぐらい示しておいて欲しいものです（私が見つけていないだけかもしれませんが）。


今度は D-P の破壊基準を使用せず、一定値で降伏するような設定にしてみました。これ、先日の文献にも使われていた手法のようですが、それっぽくはなります。が、越流前に水圧で堤体全体が若干変形してしまうのは避けられません。壊れるまでは不透水ということ、流体ベースということが効いているように思えます。また、越流直後、堤内側の斜面を削剥しながら流下する結果もイマイチ。岡山大学の文献では、堤内側にも越流した水もモデル化して計算されていましたが、こういった現象を防ぐコツなのでしょうか？

水頭差を小さくしたり、人工粘性でなく水の動粘性係数を入力したりしましたが、結局、納得のいく結果は得られませんでした。


Multi-Phase Ver.3.4β の結果を整理すると、以下の通り。

・境界条件の設定が困難。
・土砂を1層しか取り扱えない（区分できない）。
・パラメーターの単位について不明な点が残る。
・現行の支配方程式では、堤体の様な自立する土構造物をc・φで特徴付けることは困難。

最後の点は、今後、ソースが公開された時点で LSFLOW に似せれば解決するでしょう。

とりあえず、その能力は掴めました。

DualSPHysics Multi-Phase

Multi-Phase 関連で示されている文献はコチラ↓　浸透力については以下の通り。
http://www.sciencedirect.com/science/article/pii/S0309170816300926

For simplicity, it is assumed that the water does not flow in the un-yielded region and seepage only acts at the interface of the un-yielded –yielded regions and the interface (see Fig. 3.1 ). Also, the soil mixture is assumed to be isotropic and fully saturated under drained conditions. Although the assumption that water does not flow in the un-yielded region is not strictly correct for the accurate representation of the seepage forces in the soil body, Darcy law forces are approximated in the on the yield surface which is our point of interest with this article. Other SPH practitioners ( Bui et al., 2007, Sakai et al., 2009 ) modelled seepage forces using Darcy law, by using a two SPH particles layers approach to superposition the liquid and soil layer. Unfortunately, this technique tends to be cumbersome and memory intensive. GPUs are memory restricted and such a 3-D model would not be feasible with the current tech- nology.

浸透力の取り扱いは境界面のみといった簡素化をされているようです。粒子法での間隙水のモデル化は負荷が大きすぎてダメ、というのは理解できます。モデル化も現段階では難しいでしょうね。
これだと越流浸食は計算できるかもしれませんが、パイピング破壊の計算はできません。ま、こちらは現行通りFEMで十分でしょうか。
https://www.kkr.mlit.go.jp/inagawa/safe/prevention/pdf/dike_mechanism.pdf

粒子法による越流浸食の簡単なモデル化について調べてみると、他にもありますね。
https://www.okayama-u.ac.jp/up_load_files/press28/press-161216-9-1.pdf
https://www.jstage.jst.go.jp/article/jsidre/84/1/84_I_31/_pdf
http://committees.jsce.or.jp/seibu_s01/system/files/0424tousaka.pdf

個人的には、ある閾値をもって破壊とみなすといった3つ目の文献は、分かりやすいと思います。が、土質であれば、拘束圧に依存した破壊基準を導入している方がしっくりきます。この点に関しては、計算する技術者が土質ベースか、流体ベースかで異なってくるかもしれません。

天然ダムの破堤のように、越流浸食が問題となる現象では、その視点で利用できるかもしれません。
早速、試してみましょう。