水工学論文集,第50巻,2006年2月

渓床堆積物の侵食による土石流発達過程に関する数値解析
NUMERICAL STUDY ON A DEVELOPMENT PROCESS OF DEBRIS FLOWBY THE ENTRAINMENT OF BED MATERIAL IN TORRENT SEDIMENT LAYER

長田健吾1・清水義彦2・若井明彦2
Kengo OSADA, Yoshihiko SHIMIZU and Akihiko WAKAI
1学生会員
2正会員

群馬大学大学院工学研究科博士後期課程（〒376-8515 群馬県桐生市天神町1-5-1）

博士（工学）

群馬大学助教授

工学部建設工学科（〒376-8515 群馬県桐生市天神町1-5-1）

Moving liquefied soil block caused by a small landslide sometimes runs out to mountain streams over theerodible bed layer. Then, the debris flow arises and has been developed by the entrainment of bed materials into theflow body. In this study, we focus such a process by using the proposed numerical simulation, which is composed ofthe particle-fluid multiphase flow analysis by using the Distinct Element Method coupled with the unsteady 1D flow
calculation. Numerical results show the complicated time-depended process of the entrainment of debris material into
the flow body and the increasing sediment discharge of debris flow according to the supplied water volume at upperboundary condition.

Key Words : Debris flow, entrainment of bed materials, Distinct Element Method,
multiphase flow analysis, numerical simulation, development process of debris flow

１．はじめに上流域で生じた小規模な土砂崩壊が降雨の集水によって流動化し，それが堆積物の豊富な渓流河道上を流下することで土石流の発達が生じる場合がある．これは，土石流の流下に伴い，飽和した渓床堆積層を侵食することによって土石流本体に取り込み，その規模が拡大するプロセスとして表現される．江頭ら１）は，こうしたプロセスを実験的に検討しており，勾配に対応する平衡濃度との関係で堆積層の侵食を論じ，侵食速度は土石流構成材料に対する堆積層材料との粒径比に依存すること等を見出している．こうした基礎的な知見の集積が上記のプロセスの把握にとって重要であるが，水理実験に合わせ，現象の理解を進めるための数値モデルによる検討も必要である．
本研究では，こうした発達段階にある土石流の拡大過程を数値解析から検討するため，２次元数値シミュレー
ションモデルを構築した．すなわち，数値モデルは，発達段階を扱うために，従来，土石流の主な研究で展開されてきた連続体力学にもとづく手法（土石粒子と流体を一体とした混相流としての取り扱い）ではなく，土石粒子は粒子群としての挙動を個別要素法２）からラグランジュ的に追跡して求め，その結果を，粒子濃度場として取り込んだオイラ−的解析と連動させた手法とした．数値シミュレーションの対象とする場は，勾配15°の一様斜面とし，上流側での供給水量を与え，発生させた土石流を飽和渓床堆積層上を流下させて，下流末端での土砂量を求めて侵食量を評価した．また，上流端供給水量や，渓床堆積層内の粒度構成を変えて，それらが侵食量に与える影響や，土石流が流下した後の渓床堆積層における平均粒径の変化も合わせて検討した．

２．数値解析モデルの概要(1) 個別要素法による粒子運動モデル
個別要素法では，剛体要素と見なした粒子間の力の伝達が接触点を介して行われるものとして，各粒子の運動方程式を前進差分で解き，時々刻々と変化する粒子の位置を求めている．粒子間接触面における法線方向の力

f n と接線方向の力 f s は，それぞれ接触してからの相対変位量に比例する抗力（ e ）(バネモデル)と相対速度に- 907 -qw uqw
 h  zb   b Fp
 (1  c ) gh
 t
x
x
 

比例する粘性抵抗力（ d ）（ダッシュポットモデル）から構成する．

f n  ent  d n  ent  t  kn un t   n un(1)

f s  est  d s  est t  k s u s t   s u s(2)ここに，添え字のn, s は接触面における法線，接線方向を意味し， k はバネ定数， はダッシュポット粘性係

ここに， h ：水深， q w ：流体の流量， u ：流速， c ：水流相内に存在する粒子の平均濃度， zb ：河床位， bは底面せん断応力である．また， Fp は粒子群からの形状抵抗であり，以下の式で表される．

数， u は粒子ijの相対速度である． f s を構成する接線Fp 

方向バネ力については累積変位量から算定し，一定以上

g22

g e   ts(3)c c は粒子レベルのせん断強度で，クロ−ン摩擦が適用できるとした（c：粒子レベルの粘着力，  ：粒子レベルの内部摩擦角）．
(4)

上記の粒子間干渉（着目粒子iに接触するすべての粒子jからの作用力の総和）と，粒子に作用する重力，流れから受ける駆動力を組み込んだ粒子iの運動方程式（並進および回転）は次式のように書ける．

du pidtpi2x2 mi  Vi g sin m'idv pidt

粒子群の平均速度である．式(8)の右辺には，一般に、piyが（飽和堆積層内の水量の取り込み），本研究ではこれを簡単のため考慮せず，水量を一定として計算を行った．
また，水流相に関する運動量式（式(9)）には，底面粒子との間における摩擦応力項（右辺第2項）と水流相内での水流・粒子相対運動による運動量欠損項（右辺第pixpi

随性が低い場合），平衡状態になるとその値は小さく，運動量消費の主要分は底面摩擦で受け持つことになる．
ここでは，この摩擦による運動量消費を高橋・中川３）の

モデルを援用して与えた．

掃流砂の場合（ c  0.02 ）j1  c  C A d 2 u  u 2  u  v 2 u  v 2pix2 mi  Vi g cos Iid  pidt(11)3

ここに， g ：粒子数密度， u p ：水流相内に存在する   f n sin  ij  f s cos  ij j
DA3 d p大きい場合に効果を発揮し（流れに対する粒子運動の追jDc(10)3項）が含まれている．後者は，水流と粒子の速度差が   f n cos  ij  f s sin  ij j

1  c  C A d 2 u  u 2  u  v 2 u  u  (5)

侵食に伴う水体積の湧き出し項が付加されるべきである c  c  e n tan mi' (1  c )C D A2 d p h u  u p u  u pの力（ c ）が働けば滑動する条件を与えた．

est  sign c (9)piypiypib (6)

gn 2
h1 3u2(12)

掃流状集合流動の場合（ 0.02  c  0.2 ）d pi

  f 2sj(7)b j

ここに， m i'    /   C M A3 d 3 ,  ：流体密度，粒子の2次元，3次元形状係数， d pi ：粒子径， u pi , v pi :b 

粒子の流下・鉛直速度， u x , u y : 水流の水平方向および粒子体積， g ：重力加速度， I i :粒子の慣性モーメント20.49  h u2(13)(14)土石流の場合（ c  0.4c* ） ：粒子密度， CM ：付加質量係数(＝0.5)， A2 , A3 ：鉛直方向流速， CD ：抗力係数， mi ：粒子質量，Vi ：T  d p 

 T  d p 2u28  h  c  1  c    c c 1 3  1 2*

ここで， n ：Manning粗度係数，  T ：混合物の密度，である．

c* ：堆積層の粒子濃度， d p ：水流相内に存在する粒子(2) 流体モデルの基礎式と計算方法
水流相内の粒子運動による流体場への作用力を組み込んだ1次元の連続式および運動方程式は，以下のように表される．

群の平均粒径である．
本研究では，上述の基礎方程式の差分化手法に，常射流混在場の計算手法として実績のあるMacCormack法４）を適用した．MacCormack法は，予測子段階においてあらかじめ粗い近似解を求め，修正子段階で解を修正する方法
で，時間的および空間的に2次精度の差分スキームである．よって，数値振動が生じ易く，それを抑えるために (1  c ) h q w
0t
x(8)

- 908 -人工粘性が付加される．その計算方法を以下に述べる．
説明の簡単化のため，連続式および運動方程式をベクトル表示で記す．

A ECt x(15)FLOW(a)

ここに，qw1  c hAhz
Eb
uq w  1  c gh ， qw  ，
x 

（b）

0
F

pCb   ⊿Ｘである．(15)式をMacCormack法で差分化すると，以下のようになる．

(c)

図-1 粒子運動による河床位変化

予測子段階

~
t t
Ai  Ait Ei  Eit1   tC itx(16)修正子段階

Ai* 

1  t ~ t ~~ ~E i 1  E i  tC it  Ai  Ai 2x(17)

上述したように，MacCormack法では振動解を抑制するために，予測子段階・修正子段階の計算後に，人工粘性が付加され，解を修正する．本研究では，人工粘性の手法としてJamesonの方法５）を適用した．Jamesonよる人工粘性は以下のような式で表される．

Ait  t  Ai*   d Ai*1  Ai*    u Ai*  Ai*1  (18)図-2 河床位の算定方法

計算が非常に不安定になることが予想される．そのため，計算が安定して行われるようにCourant数の値を0.02と非常に小さく設定した．また， t の上限値を0.005secと設定して計算を行った．
次に，河床位の算定方法について図を用いて説明する．
 d   max  i 1 ,  i 
計算河床は計算格子内の流体の排除率を粒子濃度（体積濃度）から算定し閾値を設けて決めた．まず，図-1(a) u   max i ,  i 1 が粒子の初期形状とし黒ラインがこの粒子配置における河床位とする．その後，流体力により図-1(b)に示す緑hi 1  2hi  hi 1
i 
色の粒子が渓床堆積層から抜け出し，図-1(c)に示すよhi 1  2hi  hi 1
うに流下した粒子の一部が下流側に堆積したとする．
図-1(c)に示した⊿Ｘが流体場の差分間隔とすると，図1(c)の粒子配置における河床位形状は赤ラインのようなここに，  は数値振動の緩和係数であり，本研究では形状になると考えられる．このように，粒子配置が時々0.2として計算を行った．
刻々変動するため，河床位も粒子配置に合わせて変動さ流体場の解析を行う場合の t の値は，CFL条件であるせる必要がある．しかし，粒子配置から河床形状を算定(19)式により与えた．
するには，何らかの河床位判定基準を設定する必要があxCr(19)t る．本研究では，この判定基準を粒子濃度とし，河床位max u  ghの閾値に0.5を適用した．この粒子濃度を用いた河床位の算定方法を，図-2を用いて説明する．まず，図-2の左
図のように，流体場の計算格子間に高さが平均粒径の2
ここで， x ：流体場の流下方向格子間隔（＝2.0m），
倍程度の粒子濃度算定セルを設定し，aの位置から粒子C r はCourant数であり，本研究の解析対象が急勾配地形であり，かつ水流相内には粒子が高濃度に含まれるため， がほとんど無いbの位置まで，少しずつ上方に移動させ- 909 -

合わせ，パッキングしたものである．水流は，渓床堆積層の上流側に貯水池を設け，段波状で流下させた．本研究では，貯水池の水量を変化させた5ケースの計算を行っており，各ケースの水量は図-5に示した．
個別要素法の各パラメータ値は，清水ら６），長田ら７ ）の 研 究 を 参 考 に 設 定 し た （ k n =20000N/cm,

 n =13N・s/cm， k s =500N/cm， s =0.5N・s/cm，c=0，＝25°， tp ＝0.001s）．

図-3 計算領域の概要

４．計算結果と考察

通過粒径百分率（％）
1008060

粒度1粒度2
粒度3402000204060
粒径(cm)
図-4 渓床堆積層の粒度分布

水量(m3/m)
16012080400
ケース1ケース2ケース3ケース4ケース5

図-5 各ケースの水量ながら，各地点の粒子濃度を算定する．bの位置まで算定した中で，粒子濃度が0.5を上回る一番高い位置を判定する．その位置が図-2の右図に示した位置になったとする．そのセルの中心座標（赤色の×）と，1つ上流側の格子間で同様に算定された位置での中心座標（青色の×）とを結び，I断面と交わる点（黒色の×）をI断面の新しい河床位とした．

３．計算条件計算領域を図-3に示す．水で飽和させた渓床堆積層は長さ160m，深さ3mで河床勾配を15度に設定した．渓床堆積層の下流には，幅50m，深さ16mの給砂ボックスを設けた．渓床堆積層の粒度分布は，図-4に示す3つのパターンで計算を行った．粒度1は平均粒径を20cmとし，ある程度ランダムな粒子配列になるように最大径・最小径を平均粒径から30％ずつ増減させ，その粒度範囲で正規分布に従ってパッキングを行った．粒度2・粒度3は，粒度1の粒度分布に10％の割合で，40cmと60cmの粒子を混ぜ

図-6は，粒度1における貯水池水量がケース3の場合の，各時刻における粒子群の流動過程と，そのときの水位（青ライン）を示したものである． また，渓床堆積層の上流端から50mと100mの位置を図内に表示した．4秒の結果を見ると，段波が到達して間もないため，渓床堆積層からの侵食が始まったばかりで，流体内に含まれている粒子数は少ない．8秒になると，先端部分の形状が立っており，土石流としての流動形態になっている．その後，土石流としての形状を変えないまま約9m/sの速度で流動し，20秒あたりで下流端に到達した．24秒の結果
では，流動してきた粒子群が給砂ボックスに落下している様子が見られる．最終形状である80秒の結果も併せて表示した．図-7は，50m付近の渓床堆積層を深さ方向に4分割して各層を着色し，土石流により各層の粒子が侵食される様子を示したものである．8秒の段階で，土石流が着色層に到達し，水流の抗力および土石流に含まれている粒子群との衝突により赤で着色した層が変形しており，9秒の段階では，ほとんどの赤い粒子が土石流に取り込まれていることが分かる．緑の層も変形しているが，流動する粒子はほとんどなく，また，黒と青の層は，土石流の作用力を余り受けず，ほとんど移動していないことが分かる．図-8は，各時刻の流砂量を示したものである．横軸は渓床堆積層の上流端を0mとして下流端までを表示している．この図より，80m付近までは土石流が渓床堆積物を取り込みながら発達していることが分かる．
その後は流砂量のピーク値をほとんど変えずに流動している．
次に，粒度および水量の違いが土石流の侵食機構に及ぼす影響を考察する．図-9は，各ケースの土石流による侵食量（給砂ボックス内に入った粒子の面積の総和）と水量の比率を示したものである．ケース1は，どの粒度においても土石流として発達せず，掃流砂のような形態になった．粒度1では，ケース3まで比率が増加しているが，その後は水量が増加しても比率はほぼ横ばいになっている．粒度2も同様な傾向になっている．粒度3は，水量が少ないケースでは，粒度1の侵食量を大きく下回る．
しかし，水量が多くなると，粒度1や2の侵食量を上回る- 910 -

50m
6sec

7sec

8sec

9sec

10sec

11sec

12sec

図-7 土石流による侵食過程（粒度1，ケース3）
流砂量（m2/s）
84sec
8sec
12sec
16sec
20sec
24sec

642004080120

160 X(m)図-8 流砂量の時系列（粒度1，ケース3）

図-6 土石流の発達と流動過程（粒度1，ケース3）- 911 -

侵食量/水量0.50.40.3

る．この図を見ると，土石流の侵食作用により，渓床堆積層の中流部から上流部において表層部分の粗粒化が顕著に現れた．

粒度1
粒度2
粒度3５．まとめ

0.20.10

本研究では，個別要素法と流体層内の粒子運動による流体場への反力および粒子濃度を基礎方程式に組み込んだ1次元のオイラー的流体解析をカップリングさせた数値解析モデルを用いて、渓床堆積物の侵食による土石流の発達過程の解析を行った．上流端からの供給は水流のみとして，土石流の発達過程を考察するとともに，粒度分布の違いによる侵食量に違いについて考察を行った．
提案したシミュレーションモデルは，土石流の流下に伴う発達過程といった複雑な流砂問題を対象としているため，モデル構成上の仮定や問題点を含んでいる．例えば、本計算の枠組みにおける底面摩擦の与え方，濃度判定による河床位の設定方法，また，3次元性の効果（濃度評価における奥行きの影響など）は，今後，取り組む素過程の課題として残っている．一方，実現象の面からは，流入土砂の影響を考慮した侵食型土石流の解析も大切で，今後，これについても検討を進めていく予定である．

ケース1 ケース2 ケース3 ケース4 ケース5
図-9 侵食量と水量の関係
流砂量（m2/s）
84sec
8sec
12sec
16sec
20sec
24sec

642004080120160 X(m)図-10 流砂量の時系列（粒度3，ケース3）

粒径(cm)
30

初期形状における表層の平均粒径

2826242220040

80120160

最終形状における表層の平均粒径

参考文献

最終形状における渓床堆積層の平均粒径

1) 江頭進治・本田尚正・伊藤隆郭・有村真一：土石流による河床侵食に関する実験的研究，水工学論文集第43巻，pp.641646，1999

X(m)2) Cundall, P.A.（1979）：A discrete numerical model for granular

図-11 土石流による渓床堆積層の平均粒径変化

assemblies，Geotechnique，Vol. 29，No.1，pp.47-65．
3) 高橋保・中川一: 天然ダムの越流決壊によって形成される洪ことが分かる．水量が少ない場合は，大きな粒子はなかなか動かないため，それが抵抗になって，土石流の勢いが抑制され侵食量も減少すると考えられる．図-10に，粒度3のケース3における流砂量の時間変化を示した．
図-8と比較すると，土石流の流動速度が粒度1の場合より遅くなっており，さらに流砂量は12秒をピークに，その後は減少傾向にあることが分かる．大きな粒子の抵抗により，土石流として十分に発達できないことがこの結果から分かる．逆に，水量が多い場合は，大きな粒子が土石流に取り込まれ，その衝突力が非常に大きいため，その作用で渓床堆積物をより多く侵食し，最終的に侵食量が多くなったと考えられる．
図-11は，粒度3のケース5の場合について，渓床堆積層の区間を20mごとに分割し，各区間の初期形状（0秒）における表層（平均粒径の4倍の長さとした）の平均粒径，最終形状（80秒）における表層の平均粒径および最終形状における渓床堆積層の平均粒径を示したものであ

水・土石流のハイドログラフ，水工学論文集，第37巻，pp.669-704，1993.
4) 岡部健士・天羽誠二・石垣昌邦：常流・射流の遷移をを伴う不等流の数値計算法について，水工学論文集，第36巻，pp.337-342，1992.
5) Jameson,A.,Schmldt,W.,and Turkel,e.,Numerical Solutions If theEuler Equations by Fluite Volume Methods Using Runge-Kutta
Time –Stepping Schemes,AIAA 14th Fluid And Plasma Dynamicd
Conference,Palo Alto,California,AIAA-81-1259.1981
6) 清水義彦・長田健吾：個別要素法を用いた土石流の構成則に関する考察，水工学論文集，第48巻，pp.901-906，2004
7) 長田健吾・清水義彦・若井明彦：個別要素法を用いた流砂解- 912 -

析における問題点に関する考察，応用力学論文集，第7巻，pp.1033-1041，2004（2005.9.30受付）