第４回土砂災害に関するシンポジウム論文集，2008 年 8 月

土砂災害における固液混相流動の粒子法によるモデル化の研究
APPLICATION OF THE MOVING PARTICLE METHOD FOR MODELINGSOLID-LIQUID MIXTURE FLOW AT SLOPE FAILURE
登坂博行1，栗原宏明2，○木村慎吾3

Hiroyuki TOSAKA, Hiroaki KURIHARA and Shingo KIMURA
1

東京大学大学院工学系研究科（〒113-0033

文京区本郷 7-3-1）

E-mail: tosaka@sys.t.u-tokyo.ac.jp
2

東京大学大学院工学系研究科（〒113-0033

文京区本郷 7-3-1）

3文京区本郷 7-3-1）

東京大学大学院工学系研究科（〒113-0033

E-mail: tt076586@mail.ecc.u-tokyo.ac.jp

Key Words: Moving Particle Semi-implicit Method, slope failure, solid-fluid-mixture flow

1．序論
近年の被害地震の頻発や豪雨の集中化により，斜面崩壊と土塊の斜面移動・下方拡散，さらに河川が関係する高速の土石流などによる被害は毎年いくつも報告されている．山岳地の多い我が国では従来から積極的に対策が進められてきたが，今後地球規模の気候変動に伴う水循環変動に対する新たな対策も必要となろう．
こ の よ う な 水 文 地 質 学 的 災 害 （ hydrogeologichazards）が，山岳地の無数の斜面上のどこで，いつ，どの程度の規模で生じるかを予測することは，多くの場合困難である．しかし，危険地域においては，事前に地質調査・計測・モニタリングに加え，被害範囲・破壊力などの定量的な予測を行い，ハード（砂防ダムなどの施設），ソフト（ハザードマップなど）の対策立案・整備がきわめて重要となる．
定量的予測のためには数理的な技術が必要となるため，従来から，地すべり，崩土運動，土石流などの固液混相流に関しては，野外観測や室内実験を通じて流動の構成則の確立や数値モデルの一般性を高めるような研究が行われてきた1), 2)．それらを基に，水と固相の集合体の動きを連続体（流体）として扱う数値モデルが実用化されている．
本報告では，一般的な流体解析では表現の難しいより細かい挙動を粒子法3),4)（MPS法，Moving ParticleSemi-implicit Method)を用いて追跡再現する可能性に関して論じる．粒子法は物性の異なる多数の物質素片が相互作用しながら流動する状態を追跡するもので，ラグランジェ的流体解析に相当する．現在の所，ある程度小さなスケールの流動のメカニズムを理解する上で特に有用と考えられ利用されている．
洪水・地すべり・土石流などの自然災害の時空スケールでの解析への適用性は，計算負荷の観点から非常に難しい側面があるが，ダイナミックな固相・流体の流動状況を既存の流体解析との併用などにより捉えられる可能性がある．
本報告では，土砂崩壊後の天然ダムの浸食・決壊を模擬する実験および粒子法の計算結果，斜面崩壊のシミュレーションの試み，などについて紹介する．

2. 粒子法の概要
2.1 粒子法とはここでは MPS 法と呼ばれる粒子法について紹介する．粒子法は格子を用いないという点で，有限要素法や差分法などの従来のオイラー式の手法とは異なる．非圧縮性ＮＳ方程式

Dρ
=0Dt
1
Du= − ∇P + ν∇ 2u + gDtρ

(1)
(2)

を解く場合，勾配・発散・ラプラシアンといった微分演算子が必要となるが，MPS 法ではそれらを対応する粒子間相互作用モデル（勾配モデル・発散モデ
ル・ラプラシアンモデル）を用いて表現し，連続体の支配方程式を離散化することによって解く．なお，
式(1),(2)で ρ は密度，ν は動粘性係数， u は流速，
P は圧力，g は重力加速度であり，D / Dt は物質微分を表す．

以上の近似モデルを利用して非圧縮性流れのナビエ-ストークス方程式を解く．
MPS 法では，まず前のステップの位置と速度の値を用いて式(2)の粘性項と重力項を陽的に解いて仮の位置と速度を求める．そして式(1)と式(2)の圧力勾配項を組み合わせた圧力のポアソン方程式を陰的に解くことによって位置と速度を修正し，次のステップに進む．

2.2 離散化の方法(1) 重み関数
重み関数は粒子間距離と粒子間相互作用の大きさを関連付けるものである．
一般的な非均一粒子モデル5) では粒子間相互作用が及ぶ半径reiは粒子ごとに異なり，重み関数ωは次のように表される．

2.4 自然災害を扱うためのモデル(1) 固液混相流の取り扱い
固相と液相の混じった固液混相流は，固相粒子と液相粒子を粒子の種類番号で区別し，密度・粘性等の値を粒子の種類ごとに入力することで表現することができる．また，粒子径は粒子ごとに入力する．
固相と液相の粒子間の相互作用については，固相粒子どうし，液相粒子どうしと同様に計算する．

ω (r, rei ) =rei
− 1 （0 ≤ r < rei )r0(rei ≤ r )(3)(2) 浸食モデル水流が土砂を浸食し混相流が生じる現象を追跡するためには，水流が地面に対して及ぼすせん断力を適切に計算する必要がある．本研究では，壁との粘性をせん断力としたモデル6) を用いた．流体粒子が壁粒子から受ける力を，壁粒子も流体粒子から力を受けたとして，この力をせん断応力とし，浸食の判定値に用いるものである．

全ての近傍粒子に対して重み関数の和をとったものを粒子数密度と呼ぶ．粒子数密度 n は，重み関数
ω に粒子体積 V をかけて平均化する．

ni =∑ [V ω (r − r , r ) + V ω (r − r , r )] (4)j ≠ijjieiijiej2Vi

ここで，riは粒子iの座標，rjは近傍粒子jの座標である．
(2) 勾配モデル
粒子iにおける勾配に対して，次の計算モデルを用いる．dは次元数，jは近傍粒子，n0は内部の粒子数密度，φi,φjはスカラー変数である．

σ = 2µ ij
siju ijsrijn(9)µ ij = ρ ijν ij(10)
ここで， σ ijs は固定壁粒子 i が近傍の液体粒子 j から受けるせん断応力， ρ ij は粒子 i, j の密度の平均⎡ φ −φ⎤d値，ν ij は粒子 i, j の動粘性係数の平均値である．
ji
⎢
⎥
()
∇φ i =−−rr
V
ωrr,r
ji
jjiei
0 ∑2⎢⎥
Vi n j ≠ i r − r
相対速度のせん断方向成分 u ijs を次のように計算i⎣ j⎦ (5)()(3) 発散モデル
粒子iにおける発散モデルは次の式で表される．riは粒子iの座標，rjは近傍粒子jの座標，uiは粒子iの速度ベクトル，ujは近傍粒子jの速度ベクトルである．
⎡ u − u ⋅ r −r
⎤dji
ji
⎥ (6)
∇ ⋅u i = 0 ∑ ⎢−
ωrr,r
jiei
2⎥n j ≠i ⎢−rrji⎣⎦

した（図-1）．粒子 i から見た粒子 j の相対座標を rijn ，(4) ラプラシアンモデル
非均一粒子系でのラプラシアンモデルは次の式で与えられる．エネルギーが保存するように i と j に関して対称となっている．
⎡V jω r j − ri , rei + Viω r j − ri , rej ⎤2d∇ 2φ =⎢(φ j − φi )⎥ (7)∑iVi Λ i j ≠ i ⎢2⎣
⎦⎥

uijn

()()

(j の相対速度を u ijn とする．

u ijs = u ijn cos θ式（11）を式（9）に代入して，固定壁粒子に働くせん断応力を計算し，この値が設定した基準値より大きくなれば pick-up と判定する．

)

j’(())((juijncosθ))⎡2 V j ω r j − ri , rei + Vi ω r j − ri , rej ⎤⎥Λi = ∑ ⎢ r j − ri2Vi⎥⎦j ≠i ⎢⎣
2.3 アルゴリズム

(11)

(8)

rijni
図-1 相対速度のせん断方向成分の求め方

表-1 水柱崩壊のケーススタディ

ケース 粒子径 流体粒子数 全粒子数 動粘性係数
ν[m2/s]0.0 sec

0.2 sec

0.6 sec

0.4 sec

0.8 sec

1.0 sec

非均一A1
非均一A2
非均一B1
非均一B2
均一 1均一 210mm20mm
混合10-6360846
10-45mm10mm
混合10mmのみ

10-614405762406

1062

10-4
10-610-4※水槽の大きさ 横幅 800mm，高さ 440mm
水柱の大きさ 横幅 160mm，高さ 360mm

3. 非均一粒子を用いた場合の粒子のふるまいについての検討

図-2 水柱崩壊の計算

上方に飛んでいる粒子の個数個6050均一1均一2非均一A1
非均一A2
非均一B1
非均一B2403020
1000.00.2

0.4

0.60.81.0

時間[sec]

図-3 上方に飛散している粒子の個数

上方に飛んでいる水量-3

3[10 m ]
0.9
0.8
0.7
0.6
0.50.40.3
0.20.1
0.00.0均一1均一2
非均一A1
非均一A2
非均一B1
非均一B20.2

0.4

0.60.81.0 時間[sec]図-4 上方に飛散している水量

粒子法では計算に使う粒子の大きさにより現象の分解能が左右される面があり，計算結果も変わる可能性がある．ここでは，様々な大きさの粒子を用いて図-2 のような水柱崩壊の計算をし，その結果を比較してみた．非均一粒子モデル（サイズの異なる粒子を混合）や均一粒子モデルを用いて，計算対象のスケールはそのままで，粒子の個数・サイズなどを変えて，表-1 で示した 6 ケースについて計算を行った．サイズの大小による粒子のふるまいの違いを比較するため，図-2 の点線で示した領域より上方に飛散してしまった粒子の個数をカウントした．それを水の体積に換算して比較を行った結果が図-3,4 である．
図-2,3,4 の計算結果では，粒子サイズの違いがあっても全体の動きは似ていた．ただし細かく比較してみると粒子の飛散の仕方が違うところが見られた．
大きい粒子ほど飛散しにくい傾向が見られ，小さい粒子の方が多く飛散するというふるまいの違いが見られた．また，小さい粒子の運動が激しくなり，計算が不安定になった．粒子サイズ，すなわち空間刻み幅が小さいと，安定に計算するためには時間刻みを小さくとらなければならない．粘性値などのパラメータを変えるなどして，非均一粒子でも安定して計算できるように工夫する必要があると考える．

4．固液混相流に関する実験とシミュレーショ
ンの比較
4.1 砂堆浸食の基礎実験
水文地質学的なスケールでの固液混相流は，より細かく分類すれば，固相の割合の大きな地すべりや岩屑流から，より水の割合の高い泥流，土石流，懸濁流などがある．ここでは特に，河川を閉塞した固相を水流が押し流してゆく土石流に注目する．土石流の発生原因を力学的なメカニズムから分類すると，①渓床堆積物が外部からの水の供給によって浸食・

表-2 実験ケース
勾配
堆積形状

実験ケース

上流側 18.4°
下流側 13.4°
三角砂堆に水流を与える場合

14°
上流側 30.0°
下流側 30.0°14°
一様に堆積した砂堆に水流を与える場合

堆積厚 5cm20°

一様に堆積した砂堆に降雨を与える場合

10 cm

θ

図-514°
20°

堆積厚 5cm

流量[cm3/s]
9.2
22.9
40.8
7.8
20.9
47.0
73.4
8.6
22.4
41.7
9.0
23.2
43.6

材料

41.2

砂

砂

砂利

砂

ケース番号
A-1
A-2
A-3
A-4
A-5
A-6
A-7
B-1
B-2
B-3
B-4
B-5
B-6
C-1
C-2

80 cm

実験装置の概要

図-7 水流による三角砂堆の浸食シミュレーション（ケース A-4,5,6 に相当）

図-6

水流による三角砂堆の押し流しと混相流の発生の様子（ケース A-6）運搬される場合，②崩壊土塊が，それ自身含んでいた水によって，または外部からの水によって土石流化する場合，③天然ダムの決壊による場合，などがある．さらに②の場合は，表層崩壊と深層崩壊に分けられる．また，③の場合は，越流浸食による決壊，すべり崩壊による決壊，進行性破壊による決壊に分けられる7) ．本研究では以上の状況を踏まえ，前記の発生原因に対して粒子法の適用可能性を検討する
ため小スケールでの水・固相流動の基礎的実験を行った．実験装置の概要を図-5 に示す．装置に与える
河床勾配，流量，堆積形状，堆積物種類（砂，砂利）を変化させて複数のケースについて観察を行った（表-2）．
図-6 は，三角状に河床面に置いた砂堆に上流（右方）からの水流が越流し，下流側から浸食が進行して行く様子を示している．

図-8 水流による台形砂堆の浸食シミュレーション（ケース B-1,2,3 に相当）

崩壊後

0 sec

6 sec

17 sec200 sec
図-9 砂防ダムによる崩壊土砂の捕捉（斜面勾配：30％，水平方向長さ：200m，砂防ダムの高さ：7.5m）

4.2 浸食モデルによる実験結果の再現
壁との粘性をせん断力としたモデルを用いて水流による砂堆の崩壊の計算を行った．図-7 は表-2 のケース A-4,5,6 に相当し，図-8 はケース B-1,2,3 に相当する．左側から流入した水が砂堆の上流側に溜まり，やがて越流し，浸食が起こる様子が見られる．水が溜まっても越流するまでは浸食があまり進まず，越流すると急激に浸食が進行する．また，砂堆形状が三角形，台形ともに水流の進行方向とは逆に，下流側から浸食が始まっているという特徴が表れている．
図-7,8 は，流量等に関しては表-2 で示した実験ケースと同じ条件で計算を行ったものではないが，定性的には実験と似たような現象が表されている．
MPS 法での非圧縮性流体の計算において，壁境界では，壁粒子の速度を常にゼロに設定し，流体との間で粘性項を計算すれば，ノンスリップ境界条件を与えることになる．そのため表層壁固定粒子近傍の水粒子の流速は遅く，水面に近いほど速くなる．壁
粒子の粘性係数を大きくすることによって，粒子の飛散は通常よりも抑えることはできるが，壁粒子近傍の流速は一層遅くなる．だが，砂堆を乗り越えて流れ始めると，落下時に速度方向に重力が加わる分，壁固定粒子近傍の水粒子の速度は速くなる．
実験においてもモデル上においても下流側から浸食が始まるのは，砂堆の下流側では重力が加わるため上流側に比べて水流が速くなり，壁境界面にかかるせん断応力が大きくなったためと考えられる．このモデルでは砂堆の上流側に水が溜まっていく間は，水粒子の速度はほぼゼロと考えられるので，浸食は起こらない．実験と計算の結果から，水粒子が直接衝突した衝撃よりも，流れによるせん断応力の方が浸食に大きく関係している．

5．斜面崩壊シミュレーション
5.1 斜面崩壊とシミュレーションの意義
豪雨時などに誘起される土砂災害（斜面崩壊，地すべりなど）は，人間の生活に脅威をあたえる代表的な災害である．豪雨時に降雨の浸透により表層の地下水位（水圧）が上昇すると表層の土と下部の岩盤部の間のせん断抵抗力が弱まり，すべりが起こる．
斜面災害対策においては，崩壊土塊の到達範囲（災害の危険範囲）をある程度信頼性良く見積もることが求められ，従来からいくつかの異なる手法が研究されている．被害範囲・破壊力など，定量的な予測・推定を行うことがシミュレーションの重要な役割である．
粒子法では，砂防ダムなどの構造物を粒子で定義して設置するなどといった操作にさほど労力を必要とせず，土砂の流れ，構造物の破壊などを，あたかも目で見えるように再現することが可能である．もし，危険範囲・破壊力の定量的な予測・推定が可能になれば，災害対策・整備に十分役立つものと考えられる．本研究では，斜面崩壊をシミュレーションで模擬的に発生させ，現象再現性を検討した．
5.2

MPS 法による斜面崩壊の再現計算

図-9 は崩壊土砂が 3 箇所に等間隔に設置された砂防ダムによって捕捉される様子を粒子法によって計算したものである．左側の山の高さは 30m，砂防ダム（高さ 7.5m）の設置間隔は 40m とし，計算時間は土砂崩壊後 200 秒間とした．

粒子法では，砂防ダムの高さ・間隔を粒子の配置によって容易に調整することができ，崩壊土砂の量も調整することができる．また，土砂粒子の動粘性係数の値によって挙動が大きく違ってくるため，自然現象と比較してこれらパラメータを適切に設定することが必要となろう．また，砂防壁の強度，崩壊土砂の破壊力を考慮した計算が行えるように手法の改良が必要と考えられる．

6．結論
本研究では，土石流や斜面崩壊への粒子法の適用性を検討するために，以下のことを行った．
(1) 崩壊土塊が水をせき止めた後，浸食され土石流に変化する移行段階のメカニズムを検討するための実験を行い，粒子法の適用可能性を検討した．その結果，浸食モデルを導入した MPS シミ
ュレータを用いて実験が比較的良好に再現できることが判明した．
(2) 崩壊土砂の斜面上の運動および下方の砂防壁によって捕捉される様子を模擬的に計算した結果，定性的に想像されるような動きが示された．
現状では，MPS 法による計算は非常に負荷が大きいが，少なくとも従来の流体解析程度の精度を得ることは可能と考えられ，さらに本研究で行った一様砂堆の浸食や，側岸浸食など，3 次元性の強い現象に対し，従来の流体解析よりも精度の高い計算ができるものと考えられる．
今後は，本研究で提案したせん断応力を駆動力とする浸食モデルについて，せん断力の算定法やすべりが起こる限界値，粘性係数等のパラメータを検討していくことが課題となろう．また，土砂粒子の持つ粘性，水との界面効果などについても検討するべきと思われる．
斜面崩壊については，豪雨時の水の浸透・地下水圧の上昇などのプロセスを地下水解析の手法と組み合わせて表現するなど，より実際に合わせた追跡方法が必要である．

MPS 法では現象やスケールに応じて粒子径を自由に決めることができるが，（土砂などの）実粒子径とモデル上で用いる粒子径が大きく異なる場合は，粘性等のパラメータを実粒子と同じにはできないと考えられる．再現性を向上させるには，現象やスケールによって適切なパラメータを決める必要があるだろう．
また Koshizuka et al.によると，MPS 法において 1個の粒子の計算に関係する近傍粒子の数は，差分法における隣接格子点数よりも多くなるため，そのぶん差分法よりも計算時間や負荷が大きくなる．3 次元の場合，1 個の粒子に対する近傍粒子数は 2 次元の場合よりも格段に多くなるため，さらに計算負荷が大きくなる．特に近傍粒子検索に要する時間は，全体の計算時間に対する割合が大きい．よって近傍粒子検索を高速化するアルゴリズムが必要不可欠になる．自然スケールの 3 次元固液混相流の計算には膨大な粒子数を必要となると考えられるため，計算負荷の軽減・高速化が必要であり，今後改良すべき課題である．この点については，アルゴリズムの改良や計算機の並列化などの処理が必要であろうと考える．

参考文献
1）宮本邦明(2002):土塊の運動の 2 次元数値シミュレーシ
ョン，砂防学会誌，Vol.55，No.2，pp.5-13.
2）張馳（2006）：洪水氾濫・土砂崩壊シミュレーションの安定性・信頼性改善手法の開発，東京大学博士論文
3）越塚誠一(2005)：粒子法，日本計算工学会編，
丸善．
4) Koshizuka, S. and Oka, Y.(1996): Moving-Particle
Semi-implicit Method for Fragmentation of Incompressible
Fluid, Nucl. Sci. Eng., 123,421-434.
5）池田博和，越塚誠一，岡芳明(1998): 粒子法において局所的に空間分解能を調節するための非均一粒子モデルの開発，計算流体シンポジウム講演論文集，461
6）栗原宏明（2008）：土砂災害における固液混相流動の粒子法によるモデル化の研究，東京大学修士論文
7）高橋保(2004)：土石流の機構と対策，近未来社．
（2008.5.16 受付）