こうえいフォーラム第９号／２００１．１

DEMを用いた土石流シミュレーションの開発
DEVELOPMENT OF DEBRIS FLOW SIMULATION WITH DEM METHOD

杉山 実*・倉岡千郎*・高橋 秀*・南 哲行**・水野秀明**
Minoru SUGIYAMA , Senrou KURAOKA , Syu TAKAHASHI ,Tetsuyuki MINAMI and Hideaki MIZUNO

For modeling a numerical simulation of debris flow continuous model formulated by massand momentum equations vertically integrated has been applied to express the mixture of waterand particle flow. In this model particle property is represented by an average diameter andvertically averaged velocity according to the concept of dilatant fluid. Then it was impossible tomodel a distinct motion of the particle that is very important in case that distribution of diameterhas a wide range or effect of slit dam sets must be considerd. In this study we propose a noncontinuous model of debris flow using DEM(Distinct Element Method) technique relating theinterior experiment of debris flow with slit dam which was done by the Public Works ResearchInstitute.
Key Words : debris flow , slit dam , model , simulation , DEM

１．は じ め に
土石流の数値シミュレーションにおいては、従来、水と形式で表わし、図−２のステップで運動過程の計算を順次進めていく。

砂粒子から成る混合物流れに対し、運動方程式と連続式を垂直方向のバネ（Kn）と粘性ダンパー鉛直方向に積分して表わした連続体モデルが適用されている１）。このモデルでは、粒子は単一粒径で代表され、流速と粘性ダンパー
水平方向バネ（KS）

分布にはダイラタント流体モデルに従った断面平均流速分布が想定されている。そのため広い粒度分布特性をもつ場

モールクーロンせん断破壊モデルに基づくスライダー合や、個々の粒子の挙動が重要となる場合（透過型スリットダムの設置条件下等）には、現象を十分に再現し得るよ

図−１ 円形剛体要素モデル

うな条件設定を与えることが出来なかった。そこで今回、土木研究所との共同研究により、非連続体モデルであるDEM（個別要素法）による土石流の数値シミュレーションの開発を実施した。

２．DEMモデル概要
粒子を剛体として取り扱う剛体要素の基本的モデルは、図−１のようなバネ、粘性ダンパー、スライダー（せん断破壊機構）によって表わされる。DEMによる数値解析手法ではこのモデルに基づく運動方程式を時間に関する差分＊＊＊中央研究所 開発研究部
建設省土木研究所 砂防部

1、各要素の変位、速度、外力、強制変位設定↓
2、各要素に加わる合力計算↓
3、各要素の加速度計算↓
4、各要素の速度計算↓
5、次の時間ステップの変位計算
繰り返し

図−２ DEM計算フロー

1

DEMを用いた土石流シミュレーションの開発

３．解析モデル（2）DEM解析式（1）室 内 実 験

一次元方向の剛体運動方程式は以下のように示される。

数値解析は室内模型実験結果を再現する目的で行った。

第1項は質量マトリックスと加速度ベクトル、第２項は減

建設省土木研究所にて実施された土石流実験の概略を示

衰マトリックスと速度ベクトル、第３項は剛性マトリック

す。傾斜角17度のアクリル製水路（L=5m、H=45cm、スと変位、右辺は外力ベクトルである。

W=20cm）の基盤部に、モデル河床材料砂（図−３）が厚さ（1）5cm、長さ200cmで敷き詰められており、上流端からは一定流量（2000cm3/s）の水を投入し、下流端の土砂採取箱

この運動方程式を時間に関する差分形式微少時間（Δt）に分割し、nステップ目における運動方程式を、速度項をで流出土砂と水を採取している（図−４）。
水路下流端には、３タイプのスリットダム模型（高さ1/2ステップ遅らせた形式で次のように表わす。

18cm 表−１、図−５）を置き、土砂捕捉効果の違いにつ（2）いて比較検討を実施している。

さらに、速度と加速度を中央差分形式で表わしたものを100

加積通過率 Ｐ（％）

90代入して整理すると以下のようにまとめられる。

80706050（3）4030つまり、201000. 01と表わせることから、ステップn+1の変位を過去の変位と速度によって近似計算でき
0. 10

1 . 00

10. 00100. 00る。

粒径 ｄ（mm)

本研究では、鉛直２次元方向の砂礫の剛体運動を追跡す
図−３ 実験砂の粒径加積曲線

るため、斜面方向をｘ軸、斜面垂直方向をY軸と設定した。
式（2）、（3）の実際の計算では、接触部分における局所座標系（N,S）
を用いて、以下のように表わせる。

砂敷長200cm、厚５cm
スリットダム模型（18cm）

水路勾配θ＝17度
土砂採取箱

図−４ 室内模型実験模式図（土木研究所実施）

表−１ ダム模型タイプ諸元

ダムタイプ

スリット間隔（mm）

径間（mm）

タイプ１タイプ２タイプ３604536513627（4）（3）土砂初期条件
室内実験結果では土石流が発生・発達してほぼ平衡状態で流下してきたものを時系列的にサンプリングしている。
無施設条件における実験結果から得た土砂量及び粒度分布結果を表−２に示す。
ここで、粒径による運動の違いを表現するとともに、ス
図−５ ダム模型断面図

2

リットダム間隔との相対比に応じた捕捉過程の違いを明らこうえいフォーラム第９号／２００１．１

かにするため、３つの異なる粒径粒子を扱うこととした。

的に変化する（図−８）。

計算に用いる代表粒径は、室内実験における土石流先頭

また、ダム堆砂部下部は既に停止していることから、流部サンプルの平均粒径、75％粒径、90％粒径を参考に直径動層内の面積から求める方が良いと判断して、ここでは、8mm,11mm,17mmとし、各粒径毎の断面積配分によって平均的な流動層厚45mmに対応するダム横棒本数によってDEMモデルにおける砂礫相当個数を導いた。

各ダムを表現した。その結果、ダムタイプ１が５本、タイDEMによる粒子位置初期条件を図−６に示す。ダム位

プ２が8本、タイプ３が9本でモデル化される。

置手前80cm到達時を初期条件として、各区間毎の粒度分（5）解 析 条 件

布に従い、864個の粒子を配置した。

土石流実験現象を再現するために、DEM解析諸元（ば（4）ダムのモデル化

ね定数や粘性減衰定数等）は既往のDEM解析例等2）〜7）を参ダムモデルは同じ鋼管材料を用いながらダム格子間隔を考にしながらトライアル計算において値を変化させながら変えた３つのタイプを用いているが、本解析は縦断方向へ設定した（表−３）。ここで再現する土石流実験現象は、の2次元断面解析であるため、横断方向に平行なダム面の流下速度50.5cm/s、流動厚約45mm、タイプ３のダムにお及ぼす効果を反映させるための解析上の工夫が必要であける堆砂状況とした。

る。
方法１

計算断面内に現れるダムの横棒の断面積だけでなく、縦径間で調整

方法２

鋼管径で調整棒の断面積も考慮することによってダムをモデル化した。
この時、ダム諸元の補正にあたってはダム本数を増やし径間で調整する方法と、ダム本数はそのままに鋼管径で調整する方法が考えられるが、ここでは、DEMにおけるダム粒子（鋼）の特性を変えないために、前者の方法を用いた（図−７）。
ここで問題になるのがダム面積不透過率（流動深以下に図−７ スリットダム修正モデル位置するダム棒の面積／流動断面積）の評価であり、流動深が時間変化するために、不透過率も一定ではなく、時間

1. 0
タイプ１

表−２ 室内実験で得られた土砂条件およびDEM解析の土砂条件

タイプ２0. 8

DEM解析砂礫相当個数

8mm

11mm

17mm

8mm

11mm

17mm

2,373

58.1

17.1

24.8

137

21

1250.5

80.8

3,427

67.4

14.0

18.6

229

25

1450.5

3

67.2

2,715

72.3

12.1

15.6

195

17950.5

4

65.7

2,411

78.3

10.2

11.5

187

126

50.5

計

289.5

10,925

748

75

41

土砂量

（cm）

（ml）

175.8

2

区間速度
(cm/s)

面積不透過率タイプ３
実験における砂礫割合(%)

区間長

対応するダム本数

0. 6

9本
8本
7本
6本
5本
4本

0. 4

0. 2

0. 0
0

50

100

150

200

流動深（mm）※砂礫相当個数＝土砂量／水路幅×割合／（砂礫断面積）
速度＝（土砂＋水の総量：25458ml）／（平均流下断面：20cm×4.5cm）

図−８ 流動深変化とダム面積不透過率

図−６ DEMによる土石流初期条件

3DEMを用いた土石流シミュレーションの開発

表−３ DEMのパラメーター
記号単位

粒子密度

ρ

N･s2/mm4

鉛直ばね定数kn

接線ばね定数ks

鉛直粘性係数ηnN･s/mm

接線粘性係数ηs

N･s/mm

摩擦係数

μ−

反発係数

水の推進力の抗力係数

N/mm

e

−

Cd

既往研究(下段;本検討)の適用値範囲
最小
中間
最大
1.62×10-9

1.62×10-9

1.63×10-9

1.62×10-9−

2.62×10-9

20.6
−
N/mm
−
00009000
100
5.15
33.3
0.08210.003,0.006
0.041
0.003,0.006

9000000−36.4,3000
−
20.01
10.01

0.05〜0.1

0.2,0.580.70.20.40.6

0.05〜0.2

0.1

0.7〜0.75

0.20.45,0.7

1

−

0−

0−

6

設定値

設定で得られた見解

2.62×10-9粒子の単位堆積重量は水中重量ではなく、空気中重量の適用のほうが適切

100

今回は計算時間と収束性から決定したが、今後、比較検討要

33.3
0.00360.0036流動速度に関する重要な因子。河床粒子は完全に固0.4定状態ではないことから、河床モデルは自由粒子か、（線要素） 線要素が適切

−

流速分布

−

初期条件

重力パッキングあり〜重力パッキングなし

粒子群が比較的に疎で衝突が活発な状態で流れる状態を反映させるためには、低い値を設定するのが適切。

0.7

粒子間に水が存在することを考慮すると、固体間の摩擦係数よりも低い値を設定する方が適切

6

水が粒子を押す機構は、水と粒子の相対速度式１）６）で概ねモデル化できる。粒子堆積は水の推力が継続して加わり、既に堆積層の上部を粒子が押し進められて形成されると考えられるため、水の推力は連続的に加えてモデル化することが適切。

タイプ３のダム堆砂の影響を考慮した定常一次元流れの流速分布をMacCormack法より別途求めて、簡易的に与えた。
重力パッキングせず、初速度を粒子に与えるのが適切。
250200４．解 析 結 果
影響が粒子の挙動に与える影響が大きいと判断されたため、堆砂形状を考慮した流れについて、一次元流れ解析

450

2,000cm 3/sec

400
350

150300250

100200

水路断面150（MacCormack法）を実施し、その流速分布を別途求め（図−９）、得られた流速分布に基づいた水の推力をDEM

流量

投入位置
高さ（cm）

前章までの検討により、ダム堆砂による流速分布変化の500
水位（cm）
河床高（cm）
流速（cm/sec）

流速（cm/sec）

パラメーター

スリットダム位置

20c m

50

10050

計算において継続して加えた。その結果、流下速度・流動00厚のほか、堆積高さ及び越流挙動が実験で観察された挙動

50100150200250 300 350
距離（cm）

400450500

05 50と概ね整合した。
図−９ 流速分布計算結果

ダムタイプ別の土石流出パターンの違いを、図−10に示す。
2, 400

タイプ１ではダム径間長が平均土石径よりも十分大きい流流出を始めている。
タイプ３では最初はほとんど流出しないが、そのためにタイプ２よりも早く満砂し、3秒後には越流流出を始めている。
堆砂形状の計算結果（図−11）を見ると、タイプ１では一時的には堆砂するものの最終的には完全に流出している。タイプ２およびタイプ３には大きな違いは見られないが、タイプ３のほうが若干、満砂時間が早い様子がわかる。

42, 0001, 800
流砂量（ml/s）

タイプ２では多少流出を続けながらやがて満砂後に、越タイプ１
タイプ２
タイプ３

2, 200

ために、ダムがつまることなく、通過量は非常に大きい。

1, 6001, 4001, 2001, 000800600400200002468101214時間（s）

図−10 流出土砂量時間変化（解析結果）

1618こうえいフォーラム第９号／２００１．１また、これらの土石の堆砂のパターンには２通りあり、デルによってスリットダムを表現したが、DEMによる土一つは、ダム間に土石が直接挟まって捕捉されるもの、も石流ハザードマップ作成や施設効果計算へと実用化を図っ

う一つは、ダムの鋼管と鋼管をいくつもの土石粒子が連なていくためには、平面地形を考慮した平面二次元モデルさりながら塞ぐ、アーチ効果で捕捉されるものが、DEMをらには三次元モデル化による検討が不可欠となる。一方、用いたシミュレーションから確認された。

本検討では10秒ほどの現象を追跡するのに数時間の計算時間を要しているが、実際の土石流現象は数十秒〜数分に及

５．ま と め

ぶため、計算時間の短縮化も課題の一つである。これらの土石流のように流体と固体粒子が複雑に干渉し合いながら運動する機構についても、モデルやパラメーターをうま課題に対応していけば、この方面における新たな評価手法の確立につながっていくであろうと思われる。

く与えられれば、個別要素法(DEM)によって現象を再現出来ることを確認した。また、透過型スリットダムの二次元

謝辞：本研究は、建設省土木研究所砂防部の関係各位の多モデル化を図り、粒子の捕捉現象についての再現計算を行大な援助に依るところが大きく、ここに謝意を表します。

った。計算結果は実験結果と定性的に十分整合性のとれる結果となり、最適なスリットダム設計においてDEMモデルの適用が有効であることが示された。

参考文献
１）江頭・芦田・矢島、1989、土石流の構成則に関する研究,京都防災研究所年報、第32号B2、pp.487-500
２）河村三郎、土砂水理学１
３）内田・伯野、1990、粒状体シミュレーションによる岩屑流・土石流

６．今後の課題/展望の解析、地震研究所彙報、Vol.65、pp.321-411

今後DEM土石流モデルを、連続体力学による土石流モデルと整合性を図るため、パラメーター設定時における論理性を十分に組み立て、さらなる感度分析計算を実施していく必要がある。特に、弾性ばね係数の評価、水の推進力の評価については、要検討事項であり、物理試験等の現場

４）大倉・三森・落合、1994、個別要素法による崩壊土砂動態解析、新砂防Vol.47、No.3(194)、pp.3-10．
５）岩下・垂水・伯野、Distinct Element Analysis For Rock Avalanche
６）木山・藤村、1983、カンドルの離散剛要素法を用いた岩質粒状体の重力流動の解析、土木学会論文報告集No.333、pp.137-146
７）宮本、土石流流動のメカニズム、新砂防、p.185

計測データからの評価方法を確立しなければならない。今回は室内模型実験の再現を目的に検討したために２次元モ

T=0.3sec

T=3.0sec

T=6.3sec

TYPE１ダム

TYPE２ダム

TYPE３ダム（赤粒子は50cm以上、緑粒子は10〜25cm/s、黄色粒子は25cm/s以下の流速分布を示す。
青線は室内実験結果における堆砂面）
図−11 堆積形状の変化

5