Post-failure analysis

崩壊後の土砂流動を再現/予測するシミュレーション手法は、複数あります。

個人的には2次元平面‐浅水方程式（LSFLOW）、3次元DEM（PFC）、3次元SPH‐弾塑性＋ビンガム流体（自作）を利用したことがあります。

2次元平面の計算が軽く、3次元が重いのは仕方がありません。再現性はどれも同じでしょうか。

LSFLOWが流体っぽく対岸への乗り上げまで再現できていましたが、天然ダムの形成を再現するにはパラメーターを場所ごとに細かく設定する必要があり、予測には不向きです。

SPHは弾塑性ベース、pre-failure stage からスタートするので、土塊がまとまりながら落ちる（最初から広がらない）という実際の痕跡に近い挙動を示します。すべり面の発達過程や位置がわかる点で、他の2つよりも有利でしょうか。

DEMもSPHも、LSFLOWほどパラメータに敏感ではない印象があり、予測でも取り扱いやすいのが良い点でしょう。

以下では CEL: Coupled Eulerian-Lagrangian （ABAQUS）、MPM（ANURA3D ）、SPH（GEOXPM） を比較されています。意外と見かけない、ありがたい文献です。

Post-failure stage analysis of flow-type landslides using different numerical techniques - ScienceDirect

再現性優先でパラメータスタディするのではなく、パラメータ固定で再現性と速度、ファイルサイズを比較しています。MPIでも複数のノードを使わない条件での比較ですのでコードの能力を十分に発揮していないのですが、それでもSPHが速くて良いというのは予想と逆の結果でした。うーん、他の手法はそんなに重いのか？

もう少し情報が欲しいですね。

崩壊と粒子法

今年も台風が複数上陸しています。

これまでは東日本ばかりに上陸し、土砂災害、水害、多々発生しています。ニュースでも多くの映像が流されています。

その中で目を引いたのが、小規模な土砂崩壊の動画。

1度崩れて土砂が堆積した箇所に、再度崩土が迫ってきていました。その勢いで1度目の崩土に乗り上げたあと、水のように飛び散って谷側に堆積していました。
以前、現場で撮影していた崩壊動画も、同じように対岸にぶつかり、駆け上がり、飛び散るように挙動していました。水分を多く含む土砂は流体のようだと、あらためて感じながら見ていました。

この様な流体の挙動を厳密に再現したい場合、可視化する場合には、粒子法でないとダメでしょう。LSFLOW などで堆積形状の再現はできますし、実務上はなにも問題ないのですが、ある程度の割り切りを用いていることは、なかなかお客様に伝わりません。

鉄道分野では、既に粒子法を用いた地盤の大変形対応ソフトを開発されています。活用の程度はわかりませんが、意欲的だと思います。今後、3次元の粒子法ソフトが安価になれば、土木分野でも活躍すると思います。

台風12号は明日上陸するようです。

土砂崩れは発生しても、災害にはならないことに期待しましょう。

崩壊の順序

土砂移動シミュレーションにおいて、その予測は難しいと思います。

過去の再現と今後の予測では同一地形を用いないため、厳密には同一のパラメーターを用いることができません（問題にならない精度の現場が多いように感じていますが）。実務では目的に応じ安全側を配慮して同定値を評価・予測に利用するわけですが、本来なら地形勾配や曲率でパラメーターの相関を探るといった一手間を加える必要があるのでしょう。古い技術ですが、これらに関し十分に議論されてきたとは言えません。まだまだこれからの分野だと思います。

予測は難しいのですが、過去の崩壊の分析には非常に有効なツールだと感じています。
何度か再現解析を行っていくうちに、「ココが先に崩れないと、この堆積形状にならない」「ココの小さなすべりを無視すると、土砂の流下方向が全く異なってしまう」などといった事例に遭遇しています（以前に書き残しているかもしれません）。 つまり、現況再現をする過程で、崩壊の位置・順序などが決まる場合があるのです。
また、崩壊面の位置・形状などといった根本的な3次元地質モデルが、結果を大きく支配することを示しています。地質屋さんの力量が非常に重要になってくるのです。

シミュを行うことで、地質解釈上の2事項のあいまいさが解消します。提示した地質モデルの妥当性と推定したメカニズムの妥当性です。これらを自身で検証できるのです。これは土砂移動のみならず、浸透流や弾塑性シミュなどにも共通します。
言いっぱなしでは技術力の向上は見込めません。観察・推定・検証・修正、これの繰り返しで地質屋になれるのだと思います。

土砂移動シミュ、ありがたく利用しましょう。

LS-RAPID 土塊内部摩擦

LS-RAPID では、土塊内部の抵抗を基本的に取り扱いません。

しかし、それでは計算が破綻するようで、後付で内部抵抗に似たパラメーターを追加されています。

一つ目は、土塊内部の動摩擦係数。ただし、全土塊内部を計算するわけではないようです。土塊縁辺部において、尾根部よりも土塊カラムが高くなった場合等のみ、その土塊カラム内部でせん断抵抗を計算し、尾根越えを許容するオプションがあります。尾根越え計算をするかどうかのチェックが該当するのですが、チェックしない場合は単なるエラーでしょう。

2つ目は、 異常な速度を制限するために設けられた、α（1/2）mv^2。このチェックを有効にすると、設定速度以上でこの外力が発動するようになります。あくまで速度に応じた抵抗力を設けるというだけであり、速度を制御するパラメーターではないようです。ま、こちらも解析上の不具合をなくすために設けられたパラメーターだと思います。

これらのパラメーターが引用文献や解説書の支配方程式に入っていない点が、理解を妨げています（先の Bss も支配方程式に入っていません）。
他にも、計算中に保存則が効かず土塊体積が増えていく現象？があり、それをタイムステップ毎に割り戻すパラメーターも設けられています（これも支配方程式には入ってきません）。
このあたり、市販ソフトのためプログラムの中身が見えませんので、 手の出しようがありません。

LSFLOW と LS-RAPID、どちらもパラスタ（合わせこみ）は必要ですが、個人的には前者の方が自然に近い支配方程式やアルゴリズムを採用しているように感じます。

LS-RAPID、使い勝手は良いので今後に期待です。

LS-RAPID のτss

五大開発さんの LS-RAPID を使用することになりました。
再度、解説書を読み直しています。

今までは土研さんの LSFLOW を使用していました。崩壊後の土砂移動を扱うシミュレーションコードです。LS-RAPID も土砂移動を扱っているのですが、F=ma の考え方（組み合わせ）が異なります。LSFLOWはナビエ・ストークス経由ですので、着想がより流体ベースなのでしょう。

また、決定的に違うのが全応力（LSFLOW）と有効応力（LS-RAPID）。このあたり、開発者のこだわりが見えます。

LS-RAPID で特徴的なのが、すべり面（移動土塊と不同層の境界）の強度の考え方。すべり面の動摩擦としてφmを設定します。これは、LSFLOW と同じです。が、もう一つ、すべり面の定常状態のせん断強度τss も設定します。これは、リングせん断での残留状態から設定するようです。
では、計算で用いられているすべり面強度はどちらが優先されるのかのか？ということになるのですが、解説書に書かれていました。


定常状態のせん断強度τss の設定は、室内試験にて初期の垂直応力にかかわらず一定となる結果を利用されています。垂直応力によって変わるのが強度でなく、全応力の見かけのφａ（ｓｓ）。つまり、土塊の高さによって、τssは変化はしないとされています。とりあえずそれを境界全面に設定します。
試験を飽和状態で実施しておけば、最も強度 τss は低くなるという考え方でしょうか、定常状態のせん断面で発生する過剰間隙水圧発生率 Bss=⊿U/⊿σ3 が低くなれば（乾燥に近づけば）強度は動摩擦係数φmと有効上載圧によって決定される値まで上昇する。せん断時の過剰間隙水圧を測定するのが難しいので、実際はBss=0..3程度（尾根で乾燥側、谷部で湿潤側など変化させる）とし、φmを再現計算で同定する、といったような流れでしょう。

斜面崩壊のエネルギー

Hi-net などで観測される深層崩壊や土石流の波形から、幾つかの処理を得て Energy Parameter を算出することがあります。

この Energy Parameter 、文献での定義は明確なのですが、実現象として何を意味しているのかよくわかりません。いえ、崩壊を起振エネルギーとし、距離減衰などを踏まえた値ということは何となくわかるのですが、それ以上でも、それ以下でもありません。で、そこから先に進めません。

先輩と雑談しながら、先に進むために使えそうなデータを思いつくまま挙げていったのですが、どれも手を動かすには至らないレベルのモノでした。
せめて、崩壊時に発するエネルギーを定量的に示せないか？と考えていた時、ふと思い出したのが、以下の図書。

國生剛治「地震地盤動力学の基礎　エネルギー的視点を含めて」

斜面崩壊をエネルギーの視点から説明されている章があります。以前は読み飛ばしていた箇所です。やはり、同じようなことを考えている方は既にいらっしゃるのですね。

この図書では、Ep：土塊の位置エネルギー、Eeq：地震波動エネルギー、Edp：内部損失エネルギー、Ek：運動エネルギーの4つのバランスで地震時の斜面崩壊と土塊の移動を説明されています。

私が知りたいのは地震時でなく、常時です。常時の斜面崩壊ではEeq=0ですので、土塊が移動して停止した場合の位置エネルギーの差分が⊿Ep=Edpということになるのでしょう。
が、この図書では内部損失=摩擦面での消費というように捉えられています。個人的には、LSFLOW と、先のEnergy Parameter の話が頭にありますので、移動中の変形・抵抗に使われる土塊内部でのφmによる消費と、基盤へ与えるエネルギーとしての摩擦面φsによる消費（=Edp）を分けるべきかと考えます。そう考えないと、Energy Parameter とは？に結びつかないというのもありますが。
余談ですが、φsとφmが一定の関係になるのは、上記⊿Epの関連を考えると当然なのでしょうね。φmを崩壊前の斜面角度や崩土の安息角に関連付けると、基盤への放出エネルギー（のようなもの）を求められそうですね。

まだまだ、定量的な話には結び付きません。先は長そうです。


========================================
20150227追記

LSFLOW でのエネルギー消費概念図は土研資料のp34にあります。
ただ、内部消費φmがなぜφsと同じ形（p3）なのか、いまだに疑問です。柱状要素を前提に考えているから積分していないのだろうとは思いますが、同じ形で良いとも思えません。
割り切りでしょうね。

計算の割り切り

計算に行き詰っています。

この程度で割り切っていいか？と悩んでいた時、タイムリーなメールがやってきました。

構造計画研究所の粒子法ニュースレター
http://www.sbd.jp/product/netsu/particleworks/mail_magazine/backnumber/20150218_vol002.shtml#article2_top

で、ここでとにもかくにも声を大にして言えることは、
「実現象をパソコンの中に再現しようとするなっ！」
ですね。
これをやっちゃうと、もうドツボにはまってしまいます。
ちょっとリアルに近づいても、さらに細かいことが気になりだし、
それをクリアすると、またその先に新しい地平線が見えてきて・・・。
永久に貴方は満足ができずに条件探しの果てしない旅を続けることになってしまいます。

～中略～ 

解析自体は、まだまだ実現象を再現できていませんでした。
そこで思ったことは、
「あぁ。この解析から最低限言えることが、もう十分な結論であり、解析の成果だった。」
ということです。
知りたい情報を得ることができる最小モデルを考えて解析を行うことで、
計算負荷が驚くほど軽くなり、と言うか、計算可能となり、
しかも得られる情報は十分なものとなるのです。

ユーザーの経験談のようです。
割り切りの判断として、「この計算で何が言いたいのか？」を表現できていればOKということですね。何を確認したいか、何を再現したいのか、それを確認できたのか？を問いかければ、答えは自ずと出てくるのでしょう。
また、使用できるツールやハードの制限もありますので、どこを再現するかを取捨選択することで、解析自体が成り立つのでしょう。

「モノは言い様」ともいわれかねませんが、個人的にはタイムリーで共感できる内容でした。

連続の式・運動方程式

まずは大雑把ですが、土石流（KANAKO）の「連続の式」と「運動方程式」をやっつけましょう。

参考にした資料は以下の通りです。

• 河村哲也「河川の流れのシミュレーション」山海堂 pp.1-23（導出はこれ）
• GUIを実装した土石流一次元シミュレータ開発：中谷加奈，里深好文，水山高久，砂防学会誌，Vol.60, No.2，pp.41-46, 2008
• 土石流計算における1次元・2次元シミュレーションモデルの結合，和田孝志，里深好文，水山高久，砂防学会誌，Vol.61，No.2， p.36-40，2008
• 中村ほか「地震砂防」pp.121-123 
• 土研資料「深層崩壊に起因する土石流の流下・氾濫計算マニュアル」p23（KANAKO 使用前提のマニュでしょうか？それとも砂防では基本的な式なのでしょうか？）

＊以降、誤りがあれば、随時修正します

～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

ρを定数とし、鉛直方向に積分した連続の式は、以下の形になります。

∂h/∂t＋∂M/∂x = a

ただし、M = Uh　（流束）で、U は u の深度方向の平均速度です。

LSFLOW では a の部分が 0ですが、KANAKO では河床変化につりあう水＋土砂量（河床変動の式侵食・堆積速度）となっています。また、別途、土砂濃度を両辺に乗じ土砂量のみのつり合いも別式として加えています。
ちなみに浸透流では鉛直方向に積分する必要がないので M が vi になり、ρの変化がないと仮定すれば第1項が消え、右辺のaが注水・揚水、地盤が保持する量として考慮されます。

～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

運動方程式は、LSFLOW と KANAKO が同類ですね（浸透流はダルシー則）。
どちらもナビエ・ストークスを鉛直方向に積分して得られる「河川の基礎方程式」のようです。ただし、KANAKO は初めから粘性項を無視しています（下記の通り）。LSFLOW では、実務上粘性を0とみなす事が多いので、実質同じでしょう。uが深度方向に一定 = U と仮定すると、以下の通りとなります。

∂M/∂t＋∂uM/∂x= -gh∂H/∂x-τ/ρ

この中の∂H/∂xを水面勾配、τを河床でのせん断力で表現しているのが KANAKO。ただし、せん断力の中身については、まだ理解できていませんので、ごっそり後回しです。
LSFLOW はτをクーロン則、ニュートン型の粘性抵抗則、マニングの抵抗則として選択できます。

～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

20200530一部修正

河川の流れの基礎方程式

LSFLOW は連続の式（質量保存則）と運動方程式（ナビエ・ストークス）を利用しています。

その詳細は非公開と聞いていたのですが、土研資料にはソースまで載っていますし、砂防学会「地震砂防」という本に、支配方程式の展開や離散化の方法も載っていました。砂防分野では有名なのでしょう。

式の展開を初めて見たときは、独特な定式化だと感じていました。特に、クーロン則などを適用している外力の箇所。解釈は結構、自由なんですね。ま、再現できれば良いわけです。

また、地すべりの運動では鉛直方向の流速が水平方向に比べて無視できるという仮定のもと、式の簡略化が行われています。これを見たときは独特だなあと感じつつ、昔のことですから、計算容量を小さくしたかったのだろうなあ、程度に考えていました。

ところが先日、河村哲也「河川のシミュレーション！」という本を読んでいて、同じような定式化を見ました。河川の流れの基礎方程式というそうです。「鉛直方向は無視」という簡略化は河川砂防分野では常識なんでしょう（私が知らなかっただけのようです）。こういう考え方が河川砂防分野のベースにあったからこそ、LSFLOW も、あのような定式化をしていたのでしょうね。もっといろいろ学ばないといけません。

LSFLOW も学位論文等で FVM に改良されています。LSFLOW3、LSV と呼ばれているコードがそれに該当するのようです。
いずれにしても、inputファイルの作成～計算～可視化（MVS）までの実装は完了しましたので、あとは崩壊前の地形の調整を含めたパラスタのみです。地形の調整は時間がかかるでしょうね。
http://phreeqc.blogspot.jp/2013/10/blog-post_9060.html
もう少し考えないといけません。

崩壊前の地形

時間があったので、崩壊後、土砂がどの程度移動するかを再現していました（今度は LSFLOW です）。

最初は簡単に考えていました。

崩壊前後の LP を使用すれば、モデルは簡単にできてしまう、と考えていたのです。が、これが誤り。多くの場合、崩壊後の LP はありますが、崩壊前はありません。写真から復元したとしても、その精度はLPに比べて落ちます。場合によっては10m以上の標高差が生じます。今回も部分的に5～10m、最大20m近いギャップが生じました。

http://phreeqc.blogspot.jp/2012/11/5m_4.html

http://phreeqc.blogspot.jp/2012/11/5m-2.html

崩壊後のLPにしても、そのままではすべり面として使えません。冠頭部はそのままで使えるのですが、末端部は土砂で隠れています。そのため、ボーリング結果を使って面を作ってやらないといけません。こんな当たり前のことも、手を動かすまで気づきませんでした。迂闊でしたね。

こういった前処理に手間をかけることが、精度を大きく向上させるための一因になります。

案外、パラスタより時間がかかるかもしれません。

LSFLOW と LS-RAPID

準三次元地すべり運動解析プログラムによる地すべり性崩壊の被害範囲の予測, 土木研究所資料, 3057, 1992

LSFLOW に使われている支配方程式が掲載されています。それを順に追っていました。

が、最初から？の状態。
いえ、式が難しいわけではなく、掲載されている式に誤りがあるようです。基本的な連続の式の導出から、中の符号、余計なdxがついていたり、必要なところについていなかったり、NがMになっていたり。あちこちに？な部分がありました。ま、こちらも初心者なので、これはこれで意味があるのかもしれませんが。

基本的には、連続の式とニュートンの運動方程式を使っているだけですね。粘性項と外力項があるので、ナビエ・ストークスに近いと思います（自重項はありません）。あとは外力として内部消費が考慮されていますね。
近似はFDM。今まで、なぜ土石流のような大変形が解けるのか不思議で仕方なかったのですが、基本的には変形を解いているのではなく、流体の移動に伴う高さの変化のみを扱っているようです。ですから、粒子が移動するような可視化はできません。納得です。
計算は準三次元とありますが、高さによって物性の変化はないため、2次元に近いですね。

ソースも公開されていますので、それを読んでいますと、ここにも矛盾が。いえ、input ファイルに書かれているパラメーターよりも、読み込むパラメーターのほうが多くなっています。また、地震波の入力や計算のことは記載されていないのですが、ソースは整備されています。作成した時期にズレがあるのでしょうね。ちなみに、ソースは中央開発さん作のようです。
input ファイルの作成は EXCEL ＋ VBA で、結果の可視化は Surfer などで十分対応できそうです。


よく似たソフトに五大開発さんのLS-RAPIDがあります。
作成は佐々先生のようです。論文になっている改良そりモデルでしょうか？
http://www.godai.co.jp/soft/product/products/LS-RAPID/index.htm

こちらも連続の式と運動方程式を使用していますが、粘性をあつかっておらず、自重、外力（せん断抵抗含む）＋αといったシンプルかつ個人的には理解しやすい構成になっています。一応、間隙水圧も考慮されていますが、Jakyの式は3次元応力状態が対象でなかったはずですし、摩擦係数等のフィッティングという点ではLSFLOWと同じ扱いでしょう。離散化も差分、変形でなく高さの変化みの考慮ですので、本質的には同じつくりですね。ま、LSFLOWはすべった後、LS-RAPIDはすべり発生から扱えるという点が、大きな違いかもしれません。
http://phreeqc.blogspot.jp/2011/07/blog-post_30.html

ポスト処理としての動画が掲載されていますが、波動みたいです。うーん。

両者とも、高さのみの変化なので、可視化に工夫が必要なのでしょう。難しいですが。
やはり粒子法でしょうか？こちらはさっぱり進んでいません。理解しないと。

崩壊と地すべり その３

今年の地すべり防止工事士の更新講習に中村名誉教授の講演がありました。
先にUPした「技術者の疑問に答える地すべり・崩壊」の著者です。
http://phreeqc.blogspot.com/2011/07/blog-post_13.html

話の中心は「崩壊」でした。しかも降雨による斜面災害。
地すべり防止工事士の更新講習で地すべりを扱わず、「崩壊～土石流（流動）」とは、世の流れでしょうか？ただ、話の内容はなかなか興味深かったですね。
ちなみに、地すべりと崩壊の違いは著書でも紹介されているように移動距離とお考えです。

帰り際に話をする時間がありましたが、地すべりは全応力の安定解析、崩壊～土石流（もしくは地すべりで破壊後）はLSFLOW（市販されていません）で到達距離の把握といったように、問題を別個に考えられていました。地すべり土塊に有圧水が働いて動いた後、その地すべり土塊がどこまで移動するかは、割り切って別の計算といった方針のようです。研究の割には実務的ですね。
動いた時に、土塊にどの程度の水が含まれているかは（PC上では）関係なく、パラメーターとして動摩擦係数や斜面勾配が重要とお考えのようでした。そして動摩擦係数は逆算でデータを集めている途中のようで、準解析に使うにはまだ少し時間がかかるといった印象でした。

崩壊については表層からの浸透の話があり、やはり不飽和浸透と強度低下（サクション解放）、重量バランス、空気圧の考慮で解く方法で良いようです。解析方法が定まれば、調査提案ができますね。

また、いつ、どこが、どこまで破壊するかを考えることも重視されていました。
いつについては、降雨量に依存するので、なかなか事前予測は困難でしょう。話はそれますが、ゲリラ豪雨が来た場合の避難について住民に浸透させる方が有意義なのでしょうね。


話を聞きながら考えていたのですが、粒子法で解くと（岩盤はムリとしても）土砂の地すべり～崩壊～移動距離は統合して計算できるのでないでしょうか？すべりながら表層崩壊とか、崩壊した土砂の一部移動、末端部の押し出し崩壊から移動など、実際の現象とよく似たイメージが浮かびますね。まずは2次元で現象が統合できることを確認し、3次元で再現でしょうか？


早く理論を理解しないと、まったくの空想で終わってしまいます。
来年の目標にしようかな？2年くらいはかかるかな？