LS-RAPID 2.1 の収束性

LS-RAPID　Ver. 2.1 で、斜面横方向に向かう土砂の飛び出しが発生しました。

どうして急に勢いがついて飛び出すのか？は、以前から疑問でした。LS-RAPID で計算を行うと、頻繁に「粒子法？」と見紛う結果が出てきます。

今回のケースでは、メッシュ間隔（5m）を dt（0.01s）×v（500m/s） が超えたときに発生していました。これを起こさないように、v に応じて dt を自動設定してくれる機能があるのですが、こちらもうまく機能しないようです。いえ、機能しているのかもしれませんが、そもそも、500m/sなどと、自由落下を超えた異常な速度が発生してしまえば、dt をいくら細かく刻もうと速度はすぐに低下しないのでしょう。なぜ、このような異常値が出てしまうのでしょうか？

まず疑うのは反復計算の収束判定値が甘い点。
これについては以前、五大さんのサポートに問い合わせたことがあります。が、残念ながら、五大さんはソルバーの中身にはノータッチだそうです。京大の先生の制作になっていますから、著作権含め分担・契約がきっちりあるのでしょう。

再度、サポートさんに詳しい話を伺ったところ、LS-RAPID では「反復回数だけ決めて、収束判定をしていない」とのこと。驚愕です。ま、サポートさんも詳しいことは御存知ないのかもしれませんので、真偽のほどは分かりません。ですが、どのような異常値が出ても発散せずに（途中で計算が止まらずに）、計算が最後まで続くことや、体積がどんどん増えていくこと（支配方程式が解けていないこと）などを考えると、辻褄は合います。収束判定値が甘いといったレベルでなく、「反復回数だけ決めて、収束判定をしていない」といった理解の方が正しそうです。

そうなると、不自然でない計算結果が得られても、途中の計算で収束しているかどうか、現状では判断できません。サポートさん曰く「それなりの形になれば収束していると判断している」とのこと。つまり、異常な速度が出たり、土砂が飛び出すような結果になれば、途中で発散、収束していないと判断するそうです。うーん。
一つの手としては、計算中に体積変化をチェックし、増え出したら OUT と判定する方法があるでしょう。ただ、増えなかったことはありませんし、判定できても対策はありません。つまり、現状、お手上げです。

緩い地すべりでは体積も増えないそうです。LS-RAPIDは、このような緩い、ゆっくりした動きでないと、本来適用できないのかもしれません。急斜面での崩壊や土石流などの問題に適用するにつれ不具合が生じ、それを押さえ込むパラメーターを追加していった、という経緯のようです。支配方程式やアルゴの点では、崩壊に対しても比較的安定して答えの出せる LSFLOW の方が、現段階では優秀ということでしょう。
LS-RAPIDは、次の Ver. に期待することにして、当分寝かせておきましょう。

崩壊の順序

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

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

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

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

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

LS-RAPID の通常値

LS-RAPID のマニュアルには解析パラメーターの通常値が掲載されています。

参考にされたと書かれている文献を確認してみましたが、大きく外れている値もありました。

Sassa et al.（2010）An integrated model simulating the initiation and motion of earthquake and rain induced rapid landslides and its application to the 2006 Leyte landslide, Landslides, Volume 7, Issue 3 , pp 219-236

すべり面の運動時の摩擦係数の設定値が書かれていませんが、試験結果は39°程度になっています。

他の文献も高めです。
広島の土石流が37.5度。

佐々ほか（2014）高速長距離土砂流動現象の発生メカニズムと地すべり発生運動統合シミュレーション（LS-RAPID）を用いた広島土砂災害の再現, 国際フォーラム「都市化と土砂災害」p85

雲仙眉山は40度です。

Sassa et al.（2014）A new high-stress undrained ring-shear apparatus and its application to the 1792 Unzen –Mayuyama megaslide in Japan, Landslides, Volume 11, Issue 5, pp 827-842

通常値が25～35度とされていますので、公開されている文献は全て高めの値となっています。
土石流の運動時の摩擦係数37.5度は、一見、違和感があります。飽和時の有効応力表示としての試験値かもしれませんが、流下する段階で物性は変化しているでしょう。ま、谷部で計算に使う値は飽和のτss に設定されているでしょうから、実際、φmほど刺激は強くないと思われます。

それにしても、公開されている文献、まだ少ないですね。

LS-RAPID のノウハウ

昨年末から、LS-RAPID で複数現場、10ケース程度の計算を行いました。サポートにお話を伺いつつ、ようやく全体像やノウハウ、ソフトのクセなどがわかってきました。

すべり面のモデル化に関しては、簡易なものはSSA_3Dだけで十分です。多数のすべり面があったり、複雑な形状のモデル化はGEORAMA ＋ Civil3D ＋ SSA_3D が良いと思います。それを LS-RAPID に読み込んで、定数を設定します。

定数については、基本、再現計算でのパラスタです。
が、不具合を抑えるために追加されたパラメーターがクセモノ。流速が大きくならないよう、すべり土塊の厚さによって抵抗を発動させるパラメーターがあるのですが、最初から大きい値のままだと全く流下しません。流速によって発動させる同様のパラメーターがあるのですが、上記パラメーターを使わず、こちらを利用した方がよさそうです。いえ、本来の計算ではこのようなパラメーターを利用しなくても、それなりの流速にならないといけないのですが。意図せず体積が増え続ける、またそれを抑えるパラメーターがあるというのも理解が難しい点です。全体として支配方程式やプログラミングのブラッシュアップが必要だということでしょう。

exeは32bitですので、計算時にはメモリの制限があります。等高線を表示したまま計算しない、モデルは可能な限り小さく、グリッドは許容できる限り大き目で、といった点に留意する必要があるようです。制限に引っかかると、計算途中で落ちます。

長所もあります。
豪雨シミュは試す価値ありですね。雨によって部分的に崩壊・流下し始めます。（正解かどうかの判定は難しいと思われますが）これは良い機能だと思います。LS-RAPIDの個性でしょう。

ソフトにはそれぞれ特徴、クセがあります。それらを見極めることができれば、一つのソフトにこだわることなく問題に適したものを選定することが可能となり、しかも短時間で解決できるようになります。

崩壊・流下中のφの変化

土砂移動計算のソースを用意できたので、改変に取り組むことに。

発端は、LS-RAPIDで土砂移動を計算した際、尾根を越えるすべり面（岩盤滑り）では反対斜面へ流下するケースが出てきたことでした。
恐らく、計算開始と同時に尾根越え判定に引っかかってしまうため、内部せん断を計算したのでしょう。どうすれば良いかは明白で、内部摩擦に崩壊初期の岩盤の強度に近い値を入力してやればOK。

と思ったのですが、LS-RAPIDではダメでした。岩盤に近い強度を入れても反対斜面に流れてしまいます。カラム背後の土圧が効くのでしょうか？すべり面強度を上げて流さない様にするしか対策はない様です。イマイチ。
ま、仮に内部の摩擦係数を上げて反対斜面への流下を防いだとしても、本来流れる方向への流下も土砂としての物性が表現できませんので、ダメですね。

ということで、手元にあった他の土砂移動計算コードを改変。

以前に改変した体積変化に加え、内部の動摩擦係数も流下に伴い低下させるアルゴリズムを組み込む方針です。
コードの改変は容易で、前回の体積変化アルゴにφm及びtanφmの変更を含めるだけです。主要部分は1時間ほどで終了しました。半日ほどかけて出力を増やしたり、結果を可視化したりなど、チェックを終えて完了です。

早速、尾根越えすべりの計算！と思ったのですが、これも手元にモデルがありませんでした。仕方ないので尾根を越えない崩壊現場のモデルで計算。

結果は上々。
φの変化に応じ、徐々に崩壊・加速して最終的には一気に流れ出すような結果。岩盤が頑張って耐えている感じが見られます。以前のものは計算スタートでいきなり土砂となり、トップスピードへ向かうような結果でしたので、より現実に近づいたように感じます。


で、肝心の尾根越えすべりですが、モデルをひっぱり出して試してみました。
結果、イマイチ。やはり反対斜面にも若干流れます。改変前のコードと比べても流下範囲は大きく変わりません。パラスタが必要なのか、考え方が誤っているのかわかりません。残念ながらもう少し試行が必要でしょう。

ま、思わぬ効果が得られたので今回は良しとしておきましょう。

LS-RAPID 土塊内部摩擦

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

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

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

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

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

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

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

LS-RAPID 水圧の取り扱い

LS-RAPID の続きです。

LS-RAPID では、すべり・崩壊の発生シミュレーションを扱えるところが特徴の一つです。豪雨による水頭上昇や、地震波による崩壊を扱えます。LSFLOWでは、崩壊した後の土砂移動を扱うと割り切っていますので、決定的な違いです。

発生シミュのうち、豪雨シミュでは水頭の上昇を誘因とするため、間隙圧比 ru （=u/σ）を入力します。
ru の上昇は、水頭の上昇を模擬しています。この間はせん断に伴う過剰間隙水圧の影響は無視され、せん断抵抗角はピーク時の値φpを使用します。崩壊後の定常状態では先の Bss で補正されたτss を使用しますので、そこに水頭、過剰間隙水圧が含まれた形となります。計算上は ru = 0 として Bss にバトンタッチです。
では、ピーク時から定常状態に至るまでのせん断抵抗力低下過程では ru をどのように扱うか？というと、せん断変位で按分するようです。そのせん断変位量の閾値はリングせん断などで決めるという流れです。

ru の入力値は、解説書にあるSLIDEモデルでも、浸透流計算結果でも、観測値でも良いのでしょう。それらを使用することで、すべり・崩壊の発生から一連の流れで扱えるという点が、特徴なのだと思います。

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を再現計算で同定する、といったような流れでしょう。

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

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

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