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 にストックさせるという体制も現段階では整っていませんし、今後、施工時の調査へ場がシフトしていくのであれば、事前調査の経験等の伝達はより困難になるでしょう。ま、時代の流れ、と言ってしまえばそれまでかもしれません。

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

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

地形・地質の3次元モデルとルーチンワーク

今日は2件の簡単なヘルプに応えました。
どちらも三次元絡みです。

内 1 件は 10 本強のボーリングから 3次元地質モデルを作成するといった内容。お客様から求められた内容に対し、無償提供可能な範囲で応じたいとのことでした。
地形を作って、ボーリングを読み込んで、適当なレベルで3次元モデルを作成し 3D_PDF で返しておきました。得られた地層面の妥当性を判断する過程以外はルーチンワークなのですが、それでも3時間かかりました。

もう1件は地形の作成。こちらは基盤地図情報の 5DEM LP からコンター図を作成したいとのことでした。こちらは完全にルーチンワークです。
が、欲しいエリアに 5DEM のないことを確認。数分で終了でした。
そういえば、5DEM祭りから既に5年が立とうとしています。


振り返ってみると、今年度、土質・地質断面図を作成するにあたり 3 次元をベースとした業務が 100% になりました。いえ、昨年度までのように、ガッツリ3次元可視化を行った業務はありませんし、3 次元が仕様に入っていたものもありません。単に、「楽だから」そうなったということでしょう。

勿論、全ての地層を3次元でモデル化しているわけではありません（これをすると、労力が半端ないですから）。表層の盛土部などは、切り出した後に2次元で入れていますし、色を塗るのも切り出した後です。メインの地層面だけ PC に描いてもらい、それを 断面図の「素材」として利用しているというレベルです（2次元、3次元の良いとこどりをしているイメージです）。
5年前に↓のようなことを書き残していましたが、今はこの手順がデフォルトになりました。
https://phreeqc.blogspot.jp/2012/04/23.html


この5年間で、このような流れに変化しました。
ルーチンワークが多いので、次の5年間ではパートさんレベルでやっていただきたいなあと、思うところです。いえ、2次元CAD屋さんもいない現状では高望みしすぎでしょうか？

段階載荷と定ひずみ速度載荷

段階載荷の圧密試験では、荷重増加分比 1、荷重範囲 10～1600KN/m2、8 段階（Pc の前後 3 段階）が標準となっています。勿論、土によっては荷重範囲外を選択しても良いのですが、旧来の試験室では、5kN/m2 くらいからしか準備されていないでしょう（大きな試験室ではさらに大きな荷重、あるいは定ひずみ速度載荷装置をお持ちです）。

そうすると、難しくなるのが表層部と深部の試験。

例えば、先日の「表層部の簡易CU」と同様に、試料の中心がGL-1mとすると、有効鉛直応力（土被り圧）16-10=6KN/m2、正規圧密で同程度の Pc を想定した場合、1.25、2.5、5kN/m2の軽い側3点の載荷が必要になりす。未圧密であれば、さらに軽い側が必要となります。が、荷重を準備されていない試験室では対応できません。地盤工学会の基準書にもあるように、浚渫粘土などの表層部では、（理屈では）定ひずみを選択した方が良いものと判断されます。

一方、実際に定ひずみ速度載荷を実施するのは、洪積粘土のような硬質粘土が多いようです。段階載荷では重い側 3 点の荷重をかけることができない、あるいはPcを精度良く求めたいのに荷重が飛び過ぎて対応できない、浚渫粘土の様な超軟弱土では、シンウォールから抜いた瞬間に乱れる、自立しない、などといった経験からのようです。
沖積粘土でも GL-20m以深でN値 3 程度あれば、過圧密となっている場合に重い側で3点取れない、などの経験があります。その場合、結果的には定ひずみ速度載荷のほうが良かったなあ、ということになるのですが、なかなかうまくいきません。

現場条件や既往調査のPcを確認し、試験室で対応できるかなどを判断してから調査計画を立てないと、理想的な結果は得られません。ま、逆に言えば腕の差が出てくるところ、と言えるのでしょう。

ひずみ・応力の種類

粒子法のテキストを読んでいると、「相当応力」が出てきました。

破壊基準の種類を学んでいたころに知ったと思います。単純に検討できるよう、一軸引っ張り試験の主応力 σ1 に相当する応力に変換したものだったでしょうか？手元の図書を読み返してみると、ミーゼスのせん断応力√J2に帳尻合わせとして√3をかけたものとされていました。

以前、連続体力学を学んだころ、多様な応力名称が出てきました。これは破壊基準絡みでなく、基準配置・現在配置に関与する応力の命名が多かったように思います。例えば次の通り。
コーシー応力σ：現在配置、真応力
第1ピオラ-キルヒホッフ応力P：ツーポイントテンソル、公称応力
第2ピオラ-キルヒホッフ応力S：基準配置

ひずみに関しても同様でした。
微小変形理論では⊿t が微小であれば変位ｕも微小で、基準配置と現在配置を区別しなくてよくなる、とのことでした（定ひずみ速度載荷による圧密試験でも、関連した内容が基準書で触れられています）。その場合、有限ひずみ E（ラグランジュ - グリーンのひずみ：基準配置）≒e （オイラー - アルマンジのひずみ：現在配置）≒ε（微小ひずみ、工学ひずみ、公称ひずみ、ビオひずみ、コーシーひずみ）とみなせる、でしたね。
https://phreeqc.blogspot.jp/2011/07/2.html
おおよそ、各ひずみが 10%程度までは、概ねその差が小さいとみなせるようです。（ex.土木学会編「計算力学の常識」p99、地盤工学会「地盤の連続体力学入門講習会テキスト」）
伸び変形1%で各ひずみの差が 2.5%程度、というのも上記リンク（「よくわかる連続体力学ノート」p113）に書き残していました。
いずれも、ひずみの定義による相違といった視点と、材料特性の相違といった視点の2種をもって理解すべきなのでしょう。

土質試験では、1軸や3軸圧縮試験で応力：第1P-K もどき、ひずみ：工学ひずみの組み合わせでしょうか。微小変形を前提とした、公称応力-公称ひずみの組み合わせと解釈したいところですが、ここまでは地盤工学会の基準書に明記されていません。
（港湾関係で「破壊時のひずみが 10%以上は棄却」という慣例は、偶然ですがなんとも都合良いですね。）

一方、圧密試験では、地盤工学会の基準書にひずみの名称・区別が明記されています。こちらは大変形を扱いますので、公称ひずみではなく、自然ひずみと記載されています。1次元圧密で直径が変わらないためか応力に関する注記はありません。が、共に現在配置で真応力-真ひずみをイメージされているのだと思います。
（ちなみに、基準書によればテルツァギーの1次元圧密の式に出てくるεも現在配置の自然ひずみ（真ひずみ）だそうです。）

そうすると、dε=mv ・dσ' などの mv を使う関係は、すべて現在配置：真応力-真ひずみで揃えないといけないのでは？などと思うわけです。

そうすると、土質工学会「地盤工学における数値解析の実務」p118 の図-4.8 に示されているように、E50 から推定した E と、mv から推定した E を同一に扱ってよいのか？という疑問が出てきます。
https://phreeqc.blogspot.jp/2016/07/2.html
以前、講習会で示された圧密試験の例では、公称ひずみと自然ひずみの差が20%程度異なっていました。当然、それから算出する mv、E も20%程度異なっていると考えられます。ま、図4.8では材料のばらつきに埋もれそうですが。

取り留めないの無い話になってしまいましたが、応力・ひずみの定義は種類が多く、まだまだ理解が足りない、具体例ではよくわからない箇所がある、といったところが正直な感想です。

粒子法 （SPH） の基礎

図書館に本を返却しに行った際、新しい本を見つけました。

矢川元基・酒井譲 「粒子法　基礎と応用」岩波書店

SPH に関する基礎理論と適用事例が載っていました。先日購入していた SPH 研究所の「SPH 法の基礎と応用講座テキスト」の第1回を中心に取りまとめたような内容です。

基本的に、既に理解していた（つもりの）内容でしたので、一気に読み終わりました。が、理論の解説箇所で一部理解できていない式を発見。このような式は以前読んだ図書に書かれていたかな？などと思いながら導出を試みましたが、結局ダメでした。テキストでも、ネットで探した関連文献でも同じ所が飛ばされています。他の方が書かれた発表予稿では添字も違っており、最後は何が正しいのか？などと悩む始末。うーん。比較的重要な式なのですが。

このブログを顧みると、粒子法での解析事例を知ったのは、おそらく2010年11月。

6年以上経っていますが、私の中で大きな進展はなく、寝かせたままです。周りが進展しないことに甘えています。

最近ではチラホラ解析事例や、学会誌上での講座を見かけますが、未だに積極的に利用しようとするいる動きはありません。ソフトが高かったことも業界標準に至っていない一因でしょうか？

理論は FEM に比べ簡素なので、計算部分は自作できるかもしれません。が、プレ・ポストや CUDA まで対応しないと実務で使い物にならないでしょうし、そこまで時間をかけるなら市販ソフトを購入した方が良い、という結論になります（が、迫られていないし高価なので購入許可は出ません）。

ま、そのようなことを言っていても技術者としての進歩はありませんので、一人でできることは進めてみましょう。

表皮深さ その3

昨日、図書館に出向き、電磁気学の本をいくつか借りてきました。

午後にそれらを読んでいたのですが、分かりやすかったのが、以下の図書です。

新井宏之「基本を学ぶ電磁気学」オーム社

表皮深さの導出はp133～138。
正弦波を仮定し、マクスウェル方程式を書き下した後に平面波として導出していました。幾つか引っかかっている部分は残っていますが、思ったほど難しくはありませんでした。

この表皮深さの式は近似解でした。
厳密解に対し「金属は導電率が大きい」として一部の計算を省略されています。ま、金属のεを0やマイナスにすると、割れなかったり虚数が出てくるので、そうせざるを得ないといった背景もあるのでしょう。

では、岩石のような導電率が小さい物質でも近似解で良いのでしょうか？
どの程度の差が出るか気になったので計算してみました。周波数は NETIS （http://www.netis.mlit.go.jp/NetisRev/Search/NtDetail1.asp?REG_NO=KK-000014）に示されている5つ、それに対応する測定値（岩盤の比抵抗）は手元の図書から適当に数ケース設定しました。文献では比誘電率の幅が大きかったのですが、センシティブでなかったので7程度で固定。銅は便宜的に0に近い値として設定しました。結果は以下の通り。

案外、差が出ないですね。近似解を岩盤に適用する点は、数字上は問題ないでしょう。


残る疑問は以下の通り

①銅線の様なμmオーダーの話と、岩盤の数十～数百mオーダーの話を同じ式で扱ってよいのか？
② NETIS でいう 150m の探査深度は、ハードの制約？解像度の問題？
③深度 150m までだと、深度方向の解像度が粗すぎるのでは？上記の試算では、砂岩・泥岩などの低比抵抗を示す岩種で深度方向に2～3点、花崗岩の様な高比抵抗だと1～2点しかデータが得られない。
④地形の影響を受けないのか？

通常は、空中電磁探査で得られるセンサー毎の生データ（ppm）は示されず、そこから深度と見かけ比抵抗を推定し、さらにそれらを補完した断面図や平面図、数mメッシュのグリッドデータのみが成果として扱われているようです。断面では補間の手法の選択が大きく結果に影響することでしょう（平面はかなり密にデータがあるため補間は必要ないでしょうが）。

これまで、補間後の断面上の値をグリッド化し、差分等で解析されている文献を見てきましたが、本来は補間前の生データに対して検討すべきなのでしょうね（数が少なすぎて、できないのかもしれません）。どちらかというと、断面はおまけ程度と考え、表層部の平面データを利用する手法として割り切る方が、評価を誤らないと思われます。

***********************************
20170311追記

表皮深さを半分にして深度と比抵抗を求めた事例・手法がありました。
https://phreeqc.blogspot.jp/2017/03/blog-post.html
上記計算は端折りすぎですね。それ故、疑問点が的を得ていない物もあります。

①合わせ込みパラメータの節も見受けられますので、実務上は問題ないでしょう。
② 2つの文献の測定・計算結果を見る限り、両方の制約があるようです。理論上、用意された比抵抗モデルとNETISの周波数では、100～200mが限界となっています。これだとハードの制約になるのですが、測定例においては 900Hz、385Hz の測定値がノイズレベル以下となり、その2点の結果を定性的解釈に止めるよう強調されています。定量的解釈では3つの周波数しか使えていませんので、こちらは解像度の問題となります。
③上記の通り、理論では用意した周波数をすべて使える可能性がありますが、実際は3つしか使えず、深度方向の解像度が粗くなる例もあるということです。やはり、生データまで戻って確認する必要があるのでしょう。

***********************************
20190711追記
物理探査学会の講習会に参加。
深度を決めるには詳細法（インバージョン）と簡易法（表皮深さ）があるとのこと。通常はインバージョンにて決定。表皮深さの精度は悪く5倍とのこと。納得です。

表皮深さ その2

続きです。

空中電磁探査に用いる周波数-深度変換で、図書や文献に載っている式は以下の通り。

d= 503√（ρ/ f ）

私が見た地球物理のすべての文献で、この式の導出は省かれていました。が、おそらく次の通りでしょう。

The skin depth d （meters） = √（2/ωσμ ）=√（2ρ/ωμ ）

 ω = 2πf、非磁性体はμr≒1より、

d= √（2ρ / 2πfμ ）
　= √（ρ/ πfμ ）
　= √（ρ/ πf（1*4π×10^-7） ）
　= √（10^7*ρ/ 4π^2*f ）
　= √（2.53*10^5*ρ/ f ）
　= 503.3√（ρ/ f ）

 σ : conductivity
 ρ : Resistivity
 μ : magnetic permeability =μrμ0
 f : frequency


この係数 503 を半分にして使われている方がいらっしゃいました。今まで、その方以外で使われているのを見たことがありません。その出典や適用範囲、制約条件などを聞いても、返答はあいまいでした。

そもそも、この表皮深度の式を使えばこれ以上入らないと言った最大深度は出ますが、実際にその深度まで交流磁場が届いている保証にはなりません。その点が空中電磁探査結果の説得力の弱い一因に思われます。ま、他の手法と組み合わせて真実に迫るという点では、他の物理探査手法も同じですが。

とりあえず、まだ電磁気学をマスターしていない身です。それ故、上記の考えは私の誤りかもしれません。d=√（2ρ/ωμ ）までの導出もまだ理解していませんので、とっとと補強して土俵に上がりましょう。

***********************************
20170311追記

表皮深さを半分にして評価している事例・手法がありました。
https://phreeqc.blogspot.jp/2017/03/blog-post.html



************************************
20170325追記

物理探査ニュース「わかりやすい物理探査」電磁法3)にて、表皮深度の説明があります。
http://www.segj.org/letter/%E7%89%A9%E7%90%86%E6%8E%A2%E6%9F%BB%E3%83%8B%E3%83%A5%E3%83%BC%E3%82%B9-06.pdf

電磁探査を実施する際に、概略の比抵抗とターゲットの深度が分かれば、その深度までの探査に必要な周波数が推定されます。

これが本来のあり方なのでしょう。

John M. Reynolds 「An Introduction to Applied and Environmental Geophysics, 2nd Edition」 p413には以下の記述がありました。

A realistic estimate of the depth to which a conductor would give rise to a detectable EM anormaly is ≈δ/5

根拠は描かれていませんが、おそらく経験則でしょう。

表皮深さ

空中電磁探査の周波数→深度変換式を調べておりました。

どうも文献や図書と若干異なる式を使われているようで、どうしてそのような式になったのか、その式の適用範囲（深度）はどの程度なのか？と疑問に思い、1950年代の文献まで遡って追っていました。

で、たどり着いた（というか元に戻ってきた）のは「表皮深さ」。

以前、VLF の理論を調べた際にも出てきたことを思い出しました。電磁気学をベースとした理論式のようですが、引用文献には導出も適用範囲（深度）も書かれていません。

そもそも、電磁気学を理解できていません。いい加減に理解しておかないといけないでしょう。

ということで、大学の物理の教科書を引っ張り出して見ていたのですが、なかなか読み進めません（興味があれば既に大学時代にトントンと理解していたでしょう）。

で、ネットに分かりやすい解説は転がっていないかな？などと思い探しておりますと、（本題とは異なりますが）面白いページに目が留まりました。最も基本的なマクスウェル方程式のイメージについてユニークに説明されています。高校生でも理解できるのではないでしょうか？

http://www.yukihy.com/entry/Maxwell's4equations

表皮深度の式までたどり着くには時間がかかりそうです。が、避けては通れないようです。週末、補強に取り掛かりましょう。