水工学論文集,第50巻,2006年2月
渓床堆積物の侵食による土石流発達過程
に関する数値解析
NUMERICAL STUDY ON A DEVELOPMENT PROCESS OF DEBRIS FLOW
BY 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 the
erodible bed layer. Then, the debris flow arises and has been developed by the entrainment of bed materials into the
flow body. In this study, we focus such a process by using the proposed numerical simulation, which is composed of
the 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 upper
boundary condition.
Key Words : Debris flow, entrainment of bed materials, Distinct Element Method,
multiphase flow analysis, numerical simulation, development process of debris flow
1.はじめに
上流域で生じた小規模な土砂崩壊が降雨の集水によっ
て流動化し,それが堆積物の豊富な渓流河道上を流下す
ることで土石流の発達が生じる場合がある.これは,土
石流の流下に伴い,飽和した渓床堆積層を侵食すること
によって土石流本体に取り込み,その規模が拡大するプ
ロセスとして表現される.江頭ら1)は,こうしたプロセ
スを実験的に検討しており,勾配に対応する平衡濃度と
の関係で堆積層の侵食を論じ,侵食速度は土石流構成材
料に対する堆積層材料との粒径比に依存すること等を見
出している.こうした基礎的な知見の集積が上記のプロ
セスの把握にとって重要であるが,水理実験に合わせ,
現象の理解を進めるための数値モデルによる検討も必要
である.
本研究では,こうした発達段階にある土石流の拡大過
程を数値解析から検討するため,2次元数値シミュレー
ションモデルを構築した.すなわち,数値モデルは,発
達段階を扱うために,従来,土石流の主な研究で展開さ
れてきた連続体力学にもとづく手法(土石粒子と流体を
一体とした混相流としての取り扱い)ではなく,土石粒
子は粒子群としての挙動を個別要素法2)からラグラン
ジュ的に追跡して求め,その結果を,粒子濃度場として
取り込んだオイラ−的解析と連動させた手法とした.数
値シミュレーションの対象とする場は,勾配15°の一様
斜面とし,上流側での供給水量を与え,発生させた土石
流を飽和渓床堆積層上を流下させて,下流末端での土砂
量を求めて侵食量を評価した.また,上流端供給水量や,
渓床堆積層内の粒度構成を変えて,それらが侵食量に与
える影響や,土石流が流下した後の渓床堆積層における
平均粒径の変化も合わせて検討した.
2.数値解析モデルの概要
(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 
方向バネ力については累積変位量から算定し,一定以上
g
2
2
e
t
s
 c
g 

(3)
 c は粒子レベルのせん断強度で,クロ−ン摩擦が適用
できるとした(c:粒子レベルの粘着力,  :粒子レベ
ルの内部摩擦角).
(4)
上記の粒子間干渉(着目粒子iに接触するすべての粒子j
からの作用力の総和)と,粒子に作用する重力,流れか
ら受ける駆動力を組み込んだ粒子iの運動方程式(並進
および回転)は次式のように書ける.
du pi
dt

   f n cos  ij  f s sin  ij
D
A2 d pi
2
2
 mi  Vi g sin 
m
'
i
dv pi
dt

j
1  c  C

u
A2 d pi
2
2
 mi  Vi g cos 
d  pi
dt
(11)
3
ここに, g :粒子数密度, u p :水流相内に存在する
粒子群の平均速度である.式(8)の右辺には,一般に、
が(飽和堆積層内の水量の取り込み),本研究ではこれ
を簡単のため考慮せず,水量を一定として計算を行った.
また,水流相に関する運動量式(式(9))には,底面
粒子との間における摩擦応力項(右辺第2項)と水流相
内での水流・粒子相対運動による運動量欠損項(右辺第
随性が低い場合),平衡状態になるとその値は小さく,
運動量消費の主要分は底面摩擦で受け持つことになる.
ここでは,この摩擦による運動量消費を高橋・中川3)の

Ii
A3 d p
大きい場合に効果を発揮し(流れに対する粒子運動の追
j

j
D
c
(10)
3項)が含まれている.後者は,水流と粒子の速度差が
2
2
 u pi   u y  v pi  u x  u pi  (5)
x
   f n sin  ij  f s cos  ij
1  c  C


侵食に伴う水体積の湧き出し項が付加されるべきである
 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)
u

x
 u pi
d pi
2
  u
2
y
モデルを援用して与えた.
掃流砂の場合( c  0.02 )
j
 v pi
 u
2
y
 v pi
 (6)
b 
gn 2
h1 3
u2
(12)
掃流状集合流動の場合( 0.02  c  0.2 )
  f 
s
j
(7)
b 
j
ここに, m i'    /   C M A3 d 3 ,  :流体密度,
 :粒子密度, CM :付加質量係数(=0.5), A2 , A3 :
 d 
b  T  p 
0.49  h 
u2
(13)
2
u2

(14)

8  h  c  1  c    c c 1 3  1 2
*
粒子の流下・鉛直速度, u x , u y : 水流の水平方向および
粒子体積, g :重力加速度, I i :粒子の慣性モーメント
2
土石流の場合( c  0.4c* )
粒子の2次元,3次元形状係数, d pi :粒子径, u pi , v pi :
鉛直方向流速, CD :抗力係数, mi :粒子質量,Vi :
T  d p 
ここで, n :Manning粗度係数,  T :混合物の密度,
である.
c* :堆積層の粒子濃度, d p :水流相内に存在する粒子
(2) 流体モデルの基礎式と計算方法
水流相内の粒子運動による流体場への作用力を組み込
んだ1次元の連続式および運動方程式は,以下のように
表される.
群の平均粒径である.
本研究では,上述の基礎方程式の差分化手法に,常射
流混在場の計算手法として実績のあるMacCormack法4)を
適用した.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 

0


F



p
C
b
 
 
 
(b)
⊿X
である.(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の方法5)を適用した.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)に示した⊿Xが流体場の差分間隔とすると,図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に示した.
個別要素法の各パラメータ値は,清水ら6),長田ら
7 )
の 研 究 を 参 考 に 設 定 し た ( k n =20000N/cm,
 n =13N・s/cm, k s =500N/cm, s =0.5N・s/cm,c=0,
=25°, tp =0.001s).
図-3 計算領域の概要
4.計算結果と考察
通過粒径百分率(%)
100
80
60
粒度1
粒度2
粒度3
40
20
0
0
20
40
60
粒径(cm)
図-4 渓床堆積層の粒度分布
水量(m3/m)
160
120
80
40
0
ケース1
ケース2
ケース3
ケース4
ケース5
図-5 各ケースの水量
ながら,各地点の粒子濃度を算定する.bの位置まで算
定した中で,粒子濃度が0.5を上回る一番高い位置を判
定する.その位置が図-2の右図に示した位置になったと
する.そのセルの中心座標(赤色の×)と,1つ上流側
の格子間で同様に算定された位置での中心座標(青色の
×)とを結び,I断面と交わる点(黒色の×)をI断面の
新しい河床位とした.
3.計算条件
計算領域を図-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)
8
4sec
8sec
12sec
16sec
20sec
24sec
6
4
2
0
0
40
80
120
160 X(m)
図-8 流砂量の時系列(粒度1,ケース3)
図-6 土石流の発達と流動過程(粒度1,ケース3)
- 911 -
侵食量/水量
0.5
0.4
0.3
る.この図を見ると,土石流の侵食作用により,渓床堆
積層の中流部から上流部において表層部分の粗粒化が顕
著に現れた.
粒度1
粒度2
粒度3
5.まとめ
0.2
0.1
0
本研究では,個別要素法と流体層内の粒子運動による
流体場への反力および粒子濃度を基礎方程式に組み込ん
だ1次元のオイラー的流体解析をカップリングさせた数
値解析モデルを用いて、渓床堆積物の侵食による土石流
の発達過程の解析を行った.上流端からの供給は水流の
みとして,土石流の発達過程を考察するとともに,粒度
分布の違いによる侵食量に違いについて考察を行った.
提案したシミュレーションモデルは,土石流の流下に
伴う発達過程といった複雑な流砂問題を対象としている
ため,モデル構成上の仮定や問題点を含んでいる.例え
ば、本計算の枠組みにおける底面摩擦の与え方,濃度判
定による河床位の設定方法,また,3次元性の効果(濃
度評価における奥行きの影響など)は,今後,取り組む
素過程の課題として残っている.一方,実現象の面から
は,流入土砂の影響を考慮した侵食型土石流の解析も大
切で,今後,これについても検討を進めていく予定であ
る.
ケース1 ケース2 ケース3 ケース4 ケース5
図-9 侵食量と水量の関係
流砂量(m2/s)
8
4sec
8sec
12sec
16sec
20sec
24sec
6
4
2
0
0
40
80
120
160 X(m)
図-10 流砂量の時系列(粒度3,ケース3)
粒径(cm)
30
初期形状にお
ける表層の平
均粒径
28
26
24
22
20
0
40
80
120
160
最終形状にお
ける表層の平
均粒径
参考文献
最終形状にお
ける渓床堆積
層の平均粒径
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 the
Euler 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受付)
ダウンロード

渓床堆積物の侵食による土石流発達過程 に関する数値解析 - C