１．一般座標系による土石流・泥流モデルの開発（担当：竹林洋史）（１）概要
土石流・泥流の流速，流動深（水深），氾濫範囲を精度良く予測することは，砂防ダムなどの土石流対策工の設計，避難可能時間の推定，土石流危険区域の設定等において非常に重要である．しかし，従来，土石流・泥流の流動・氾濫に関する平面二次元解析では，基礎方程式はデカルト座標系で記述し，長方形格子を用いた解析を行っている．長方形格子では，河岸と流れの境界が格子を斜めに横切る場合は，境界部において流れの抵抗が大きくなり，土石流の流動深が厚く見積もられる．さらに，土石流の流下速度が遅く評価されることとなる．その結果，土石流の発生をセンサーが検知してから下流域の宅地などへ到達する時間を遅く見積もられるため，避難時間の予測を危険側に見積もることとなる．また，土石流による流体力も小さく見積もられることとなるため，河川構造物の設計においても危険側となる．
本研究では，デカルト座標系で記述された平面二次元の泥流の支配方程式を一般座標系の支配方程式に書き換え，蛇行河道形状に沿った境界適合格子を用いた泥流の平面二次元数値解析を実施し，湾曲部における泥流の流動特性及び浸食・堆積特性について検討を行った．
（２）基礎方程式
デカルト座標系と一般座標系との関係は，以下のようである．

J=1
 ∂x ∂y ∂x ∂y  ∂ξ ∂η − ∂η ∂ξ ∂ξ∂y=J∂x∂η∂η∂y= −J∂x∂ξ∂ξ∂x= −J∂y∂η∂η∂x=J∂y∂ξ

(1-1)

(1-2)
(1-3)
(1-4)
(1-5)

ここに，ξとηは一般座標系における流下方向と横断方向の座標であり，x と y はデカルト座標系での座標である．
水深平均された平面二次元流れの支配方程式を用いる．水及び土砂の混合物の質量保存則は以下のようである．

∂  h  ∂  hU  ∂  hV  E
 + + =∂t  J  ∂ξ  J  ∂η  J  c* J(2)

ここに，t は時間，z は水位である．表面流の水深は h である．U と V はξとη方向における反変速度成分であり，以下のように定義される．

=U∂ξ
∂ξ
u+ v
∂x
∂y-2-(3)

=V∂η
∂η
u+ v
∂x
∂y(4)

ここに，u と v は x と y 方向における流速成分，c*は静止堆積層の砂礫の堆積濃度，E は河床の浸食速度であり，以下の式 1)を用いる．

E= c* tan ( θ − θe )u + v22(5)ここに，θは流動方向の河床勾配であり，以下の関係がある．

sin θ =u sin θ x + v sin θ yu 2 + v2(6)ここに，θx は x 方向の河床勾配，θy は y 方向の河床勾配である．θe は水及び土砂の混合物中の鉛直平均
土砂濃度 c に関する流動方向の平衡河床勾配であり，以下の関係がある．

σ  ρ − 1 c hstan θetan φs=h
σ  ρ − 1 c + 1(7)ここに，φs は土砂の内部摩擦角である．流体中の土砂の質量保存則は以下のようである．

∂  ch  ∂  chU  ∂  chV  E
 + + =∂t  J  ∂ξ  J  ∂η  J  J(8)運動量保存則は以下のようである．

∂  hU  ∂  hU  ∂  hU +U+V∂t  J  ∂ξ J  ∂η J −hu  ∂  ∂ξ ∂  ∂ξ   hv  ∂  ∂ξ ∂  ∂ξ      − U   + VU   + VJ  ∂ξ  ∂x ∂η  ∂x   J  ∂ξ  ∂y ∂η  ∂y  
 1   ∂ξ 2  ∂ξ 2  ∂ z
1  ∂ξ ∂η ∂ξ ∂η  ∂ zb =− gh     +    b + +
 J   ∂x   ∂y   ∂ ξ J  ∂x ∂x ∂y ∂y  ∂ η  
221  1   ∂ξ   ∂ξ   ∂P 1  ∂ξ ∂η ∂ξ ∂η  ∂P  τbξ
−−    +   + +ρm  J   ∂x   ∂y   ∂ ξ J  ∂x ∂x ∂y ∂y  ∂ η  ρm J

-3-(9)∂  hV  ∂  hV  ∂  hV 
+U+V∂t  J  ∂ξ  J  ∂η J −hu  ∂  ∂η ∂  ∂η   hv  ∂  ∂η ∂  ∂η      − U   + V
U   + VJ  ∂ξ  ∂x ∂η  ∂x   J  ∂ξ  ∂y ∂η  ∂y  22 1  ∂ξ ∂η ∂ξ ∂η  ∂ z
1   ∂η   ∂η   ∂ zb b=− gh  ++ + J  ∂x ∂x ∂y ∂y  ∂ ξ J   ∂x   ∂y   ∂ η 




(10)

2
21  1  ∂ξ ∂η ∂ξ ∂η  ∂P 1   ∂η   ∂η   ∂P  τbη
−−  ++   +   ρm  J  ∂x ∂x ∂y ∂y  ∂ ξ J   ∂x   ∂y   ∂ η  ρm J

ここに，g は重力，zb は河床位である．P は圧力であり，以下の関係を用いる．

P =ρ
∫ m g ( h − z ) cos θdzh(11)0

ここに，z は鉛直方向の軸である．ρm は以下の関係がある．

ρm = ( σ − ρ ) c + ρ(12)ここに，ρは水の密度，σは土砂の密度である．τbξとτbηはξとη方向における反変掃流力成分であり，以下のように定義される．

∂ξ∂ξ
τbx + τby
∂x
∂y
∂η
∂η
τbη =
τbx +τby∂x∂yτbξ =(13)(14)

ここに，τx とτy は x と y 方向におけるの掃流力成分であり，以下のようである．

τbx ={τ + ρf ( u + v )} u u+ v

τby ={τ + ρf ( u + v )} u v+ v

2y2b22y(15)22b2(16)2

ここに，τy は降伏応力であり，以下の関係 1)を用いる．
1nc=
τ y   ( σ − ρ ) cgh cos θ tan φs c* (17)

ここに，n は常数であり，ここでは 5 を用いる．fb は抵抗係数であり，土石流の場合の以下の関係を用いる．

72α 2 ， α =f=b(25  1 − c=fb
k f4 c23κ6) + k σ 1 − e c  h 
53d-4-ρ(2)  d 13
（乱流域）(18)（層流域）1)

(19)

−2

斜面崩壊
80°（a）Case1 の河道形状

斜面崩壊

80°（b）Case2 の河道形状

斜面崩壊

40°（c）Case3 の河道形状斜面崩壊80°（d）Case4 の河道形状
図 1 解析に用いた河道形状-5-

（a）25 秒後（b）50 秒後
図 2 Case1 の流動深の平面分布の時間的な変化

ここに，kf=0.16，kd=0.0828，e は粒子の反発係数，d は土砂の平均粒径である．

E∂  zb  = −c* J∂t  J (20)（３）解析条件
図 1 に、解析格子と地盤高の平面分布を示す．解析では，2 種類の解析格子と 3 種類の河道形状を用いて解析を行った．河道幅と河道長はすべて，河道幅 120m，河道長 1500m である．図 1（a）は，Case1 の河道形状である．流路中心線に沿った河床勾配 0.3，最大蛇行角度 80°の蛇行流路であり，境界適合-6-

（a）25 秒後（b）50 秒後
図 3 Case2 の流動深の平面分布の時間的な変化型の格子を用いている．図 1（b）は，Case 2 の河道形状であり，Case 1 と同様の河道形状であるが，正方形格子を用いている．なお，Case 1 と Case 2 の流路内の平均格子幅は同一である．図 1（c）は，Case3 の河道形状であり，流路中心線に沿った河床勾配 0.3，最大蛇行角度 40°の緩蛇行流路であり，境界適合型の格子を用いている．図 1（d）は，Case 4 の河道形状であり，流路中心線に沿った河床勾配 0.15，最大蛇行角度 80°の緩勾配蛇行流路であり，境界適合型の格子を用いている．最大浸食深は全ての解析で 1m とした．斜面崩壊は流路上流端付近で発生させ，流下方向 45m，横断方向 120m，深さ 1m の大きさとし，泥流発生の起点とした．

-7-

（a）25 秒後（b）50 秒後
図 4 Case3 の流動深の平面分布の時間的な変化

解析結果
図 2 は，Case1 の 25 秒後と 50 秒後の流動深の平面分布を示す．図 2 に示すように，湾曲部における泥流先端部は，内岸部が速く下流に伝播し，外岸部が少し遅れて伝播することがわかる．これは，初期河床横断形状が平坦であることに加え，河床勾配が非常に急であるため，泥流の直進性が非常に強いためと考えられる．また，蛇行流路を伝播する中で，土石流・泥流の形状が変形を受け，流動深さのピークが複数形成されるとともに，ピークの場所が土石流・泥流の先端部よりも上流にシフトしていることがわかる．
図 3 は，正方形格子を用いた Case2 の 25 秒後と 50 秒後の流動深の平面分布を示す．図 2 に示す Case1
では，50 秒後には土石流・泥流のフロントが既に下流端に到達しているが，図 3 に示すように，Case2の 50 秒後では，最後の直線部に流入使用としている段階であり，土石流・泥流のフロントの伝播が遅く-8-Case 1
流下距離が短い

流下距離が長いCase 3

図 5 Case1 と Case3 の流下距離の違い（a）25 秒後（b）50 秒後
図 6 Case1 の浸食・堆積深の平面分布の時間的な変化なっていることがわかる．これは，正方形格子を用いた場合は河岸側壁境界がなめらかではないため，物理的には存在しない数値的な抵抗が働いたためと考えられる．
-9-

（a）25 秒後（b）50 秒後（c）75 秒後
図 7 Case4 の浸食・堆積深の平面分布の時間的な変化- 10 -

図 4 は，Case3 の 25 秒後と 50 秒後の流動深の平面分布を示す．Case1 と比べると Case3 は緩蛇行であるため，早く下流に土石流・泥流フロントが到達すると予想された．しかし，図 4 に示すように，50 秒後において，Case3 の土石流・泥流フロントは下流端へ到達しておらず，Case3 の方が土石流・泥流フロントの下流への到達が遅くなっている．これは，図 5 に示すように，土石流・泥流が勾配が急となる湾曲内岸側に沿って流れるためである．土石流・泥流が湾曲内岸側に沿って流れると Case1 の方が Case3よりも下流端への距離が短くなり，結果的に早く下流に流下する．ただし，本解析の初期河道横断形状は水平であるため，河道の横断形状次第では結果も変わってくるものと思われる．
図 6 は，Case1 の 25 秒後と 50 秒後の浸食・堆積深の平面分布を示す．図 6 に示すように，土石流・泥流の通過後は，河道内の土砂が全て浸食されていることがわかる．これは，河床勾配が非常に急で浸食速度が非常に速いことと浸食深さが浅いためである．一方，河床勾配の緩やかな Case4 では浸食だけで無く土砂が堆積する場所も現れる．図 7 は，Case4 の 25 秒後と 50 秒後の浸食・堆積深の平面分布を示す．図 7 に示すように，土石流・泥流によって輸送されてきた土砂が湾曲上流域の湾曲外岸側に堆積していることがわかる．これは，掃流砂が卓越する河川とは全く異なる現象である．つまり，掃流砂が卓越する河川では，湾曲下流域では流れが湾曲外岸側に寄るため，湾曲内岸側に土砂が堆積する．これは，本解析条件の土石流は，湾曲内岸側に沿って流れていることが一つの原因であり，掃流砂卓越河川で見られる現象とは大きく異なるため，注意が必要である．さらに，堆積域の時間変化を見ると，湾曲外岸側に堆積した土砂が再浸食されて，堆積厚さが減少していることがわかる．上流端の直線河道から湾曲部に入る地点の右岸側に 25 秒の時点で形成されていた堆積域は，75 秒後にはほぼ消滅していることがわかる．
（４）本章のまとめ
デカルト座標系で記述された平面二次元の土石流・泥流の支配方程式を一般座標系の支配方程式に書き換え，蛇行河道形状に沿った境界適合格子を用いた土石流・泥流の平面二次元数値解析を実施し，土石流・泥流フロントの伝播特性に与える座標系及び格子形状の影響の検討を行った．さらに，湾曲部における土石流・泥流の流動特性及び浸食・堆積特性について検討を行った．その結果，同じ格子スケー
ルであれば，正方形格子を用いた場合は河岸で数値抵抗が発生し，境界適合型の格子を用いるよりも土石流・泥流の伝播速度が遅くなった．また，最大蛇行角度が小さい方が土石流・泥流フロントの下流への到達が遅くなった．これは，土石流・泥流が，勾配が急となる湾曲内岸側に沿って流れるためである．
さらに，掃流砂が卓越する河川では，湾曲下流域では流れが湾曲外岸側に寄るため，湾曲内岸側に土砂が堆積するが，土石流・泥流によって輸送されてきた土砂は湾曲上流域の湾曲外岸側に堆積した．

２．Morpho2DH へのハイエトグラフ読み込み機能の追加（担当：竹林洋史）（１）概要
これまで，Morpho2DH では降雨量の空間分布を考慮することは可能であったが，降雨量の時間的な変化を考慮した解析は不可能であった．ここでは，Morpho2DH の definition ファイルを修正してハイエトグラフの読み込みを可能とし，降雨量の時間的な変化を考慮した解析を可能とした．

- 11 -

（２）iRIC インターフェースa）ハイエトグラフの読み込みによる降雨の時間変化の考慮

図８

計算条件設定画面

図 8 に示すように，Morpho2DH の計算条件設定画面の「境界条件」グループに「降水量」の項目を追加した．

図 9 ハイエトグラフデータ及びグラフ表示画面「降水量」の「編集」を押すと，図 9 のような画面が表示される．左側に降水量の時間的な変化の数値データ，右側にハイエトグラフが表示される．降水量データは，「追加」ボタンを押すことにより，値を直接入力可能である．また，ハイエトグラフデータを別途準備しておけば，左下の「インポート」ボタンを押すことによってハイエトグラフデータを読み込むことが可能である．

- 12 -

図 10 ハイエトグラフデータ選択画面「インポート」ボタンを押すと図 10 のようなファイル選択画面が表示される．ここでは，Rain.csv というファイルを選択する．

図 11 ハイエトグラフデータの読み込み完了画面

ハイエトグラフデータを読み込むと，図 11 に示すように，画面上に降水量の値とハイエトグラフ（但し，一般的な棒グラフ表示では無く，折れ線グラフ表示）が表示される．なお，降雨量データの単位は mm/h としている．

- 13 -

（３）本章のまとめ
Morpho2DH のソルバーおよび definition ファイルを改良し，降雨量の時間的な変化を考慮できる土石流・泥流モデルを開発した．また，降雨の時間的な変化を考慮する場合の iRIC インターフェース上での操作方法を説明した．
参考文献
1) 江頭進治・伊藤隆郭：土石流の数値シミュレーション, 日本流体力学会数値流体力学部門Web会誌,第12巻, 第2号, pp. 33-43, 2004．

- 14 -