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

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

DualSPHysics Ver.4.0

以前、試したことのある DualSPHysics を確認。

Ver. が上がっています。現在は Ver.４.。機能の内、以下の内容が目に留まりました。
http://www.dual.sphysics.org/index.php/news/

 • Coupled SPH & DEM.
The Discrete Element Method (DEM) allows for the computation of rigid particle dynamics, by considering contact laws to account for interaction forces. The coupled numerical solution, based on SPH and DEM discretisations, resolves solid-solid and solid-fluid interactions in broad range of scales [Canelas et al., 2016]. The source files of DEM implementation are released in version 4.0. Examples: CaseBowling and CaseSolids.
 • Multi-Phase soil-water.
The DualSPHysics code has been validated for multi-phase simulations involving water and sediment for fully saturated flows [Fourtakas and Rogers, 2016]. This has been released in v4.0 as executables with the source code to following v4.2. Example: CaseTwoPhases.

確認してみますと、VTK ファイルを読み込むサンプルファイルがついていました。以前は、地形などを読み込む方法が分からなかったのですが、参考になりそうです。

https://phreeqc.blogspot.jp/2014/02/dual-sphysics.html

で、早速動かしてみました。

が、動きません。GPU のエラーが出ます。最新のドライバーにしてみましたが、ダメ。
他のサンプルファイルは GPUに載せても動きますが、これだけはダメ。残念。
サンプルファイルを読み込んだ段階で書き出される、粒子配置を記載した VTK ファイルを ParaView で確認してみますと、おお、しっかり粒子が配置されて計算できる状態にはなっています。
HP には表流水の流下の様子がUPされていますので、地形の取り込みも可能なのでしょう。まずは簡単な形状の取り込みから始める必要がありそうです。

Coupled SPH & DEM のサンプルファイルは CPU でも GPU でも計算できました。剛体を SPH で計算するのか DEMで計算するのか選択できるようです。同じ質点系であれば、容易かもしれませんが、そこはどのようなアルゴなのか理解できていません。後ほど、ですね（使うことがないかもしれませんが）。
計算自体は GPU を使っても思ったほど早くなりませんでした。結果は以下の通り。Blender でのレンダリングについては、まったくやり方が分かりませんが、プロに聞いてみれば大丈夫でしょう。同様のシミュレーションで、ブロックにかかる圧力を抜いている文献がありますので、まずはそちらを理解したいですね。

もう一つ。
Multi-Phase のサンプルファイルの中を見てみますと、土砂のcφを取り込んでいました。破壊基準は Drucker-Pruger （円錐）を選択できるようです。
流体の種類をいくつか選択できるようですが、理解できていません。以前、ビンガム流体としてτを導入している文献を見たことがありますが、同じような組み込み方なのでしょうか。わかりません。
土砂の単体は飽和での値を入力するようです。これだと間隙水の考慮ができないなあ、と思っていましたが、よく考えると FEM と同様ですので連成させないとダメですね。このあたりは文献に詳述されているようです。土砂の取り扱いに関しては LSFLOW の様な支配方程式を組み込むことも考えていましたが、これ、使えるでしょうか？
現在のところ Ver. Multiphase_3.4β ということで、説明書は添付されていません。今後のリリースで理論や計算式、解説などが公開されるかもしれません。期待しましょう。


******************************
20170302追記
レンダリングについてプロと話をしました。
回答としては、何らかのスクリプトを組まないと大変、とのこと。
手順の公開に期待しましょう。

SPH粒子法の基礎と応用講座

粒子法について復習していました。

先日、SPH研究所「SPH粒子法の基礎と応用講座」全3冊を購入しました。
時間のある時に読み進めていたのですが、これ、読みにくいと思います。

まず、誤記が多く、正誤表もわかりにくいという点。正誤表は1冊目のみで、2冊目、3冊目はついていません（問い合わせても回答なしでした）。1冊目も2章は誤記が多いので総入れ替えということで、プリントアウトした紙がホチキス留めでついていました。
誤記が多いと、理解できない数式を信用してよいのか迷ってしまいます。

exe とソースもついていたのですが、内容が異なるようです。実行形式で読み込むサンプルファイルのパラメーター記載順序・種類と、ソース内での読み込み順序・種類が異なっています。
また、実行形式の解説もついていたのですが、そこの説明にないパラメーターがサンプルファイルに記載されています。これらについても問い合わせをしましたが、まだ回答ありません。ま、大学発ベンチャーのようですので気長に待ちましょう。

ま、自力で組んでいくことも考えましたが、やはり、ネックはプリポスト。ポストは良いとしても、プリは労力が必要ですし、面白くないところです。そう思うと、SPH や MPS の基礎理論を理解しておけば（まだ残っている部分はありますが）、個々のソースまで理解する必要はないかと思い直し、他の書籍も再読することに。

で、基礎理論に関する復習はOK。なにか良いソフトはないか、探してみましょう。

トンネルの事前調査・施工時調査

トンネルの事前調査のコアを分析していた方と話をしていました。

「今まで見たことがないくらいスメクタイトが出ている」とのこと。

ま、見た目からして大丈夫そうな岩なので、バルクでかけると「微量」という判断になると思われます。バルクでも目立てば岩石の吸水膨張試験や浸水崩壊度試験を実施すれば良いでしょう。

そういえば、昨年8月にトンネル標準示方書が改訂されています。膨張性の指標の一部が変更されたそうですが、未確認でした。確認してきましょう。

http://committees.jsce.or.jp/tunnel/node/124

先日、災害科学研究所「トンネル技術者のための地盤調査と地山評価」 講習会に参加してきました。同財団3冊目のトンネル関連の図書です。前2冊が先進的でしたので期待していたのですが、残念ながらあまり多くの知見は得られませんでした。

今回の図書では、事前調査と施工時の調査に大きく分かれているように思えます。前者については総括的な内容でしたが、後者については執筆に力を入れられたようでした。トンネルの事前調査には限界があり、今後は施工時の調査を併用してリスクを低減させる、といった大きな流れが読み取れます。

勿論、事前調査でも得られるものはありました。

弾性波探査や比抵抗探査を実施しても、検層で検証・補正されていない例が多いとのこと。確かに、速度検層は実施しますが、電気検層は実施していないですね。指摘を受けて改めて気付きました。反省。

トンネルの事前調査の件数自体、昔より減ってきていると感じます（維持管理のための点検・調査業務は増えています）。事前調査の経験や知見を授受できる機会が減ってきていますし、それに係わる技術者数も減少傾向にあると感じます。私を含め、その質も落ちているように思います。経験等を AI にストックさせるという体制も現段階では整っていませんし、今後、施工時の調査へ場がシフトしていくのであれば、事前調査の経験等の伝達はより困難になるでしょう。ま、時代の流れ、と言ってしまえばそれまでかもしれません。

施工時の調査が施工に有益かつ効率的なのは私も実感していますし、維持管理段階で有用なことも理解しています（自身、施工時の調査を実施していた時は、事前調査をほとんどあてにしていませんでした）。

今後は事前調査の経験等を失わないようにする一方、維持管理を見据えた施工時調査の拡大に努めなければならない、という見解にようやくたどり着いたということでしょう。