水工学論文集,第51巻,2007年2月
水工学論文集,第51巻,2007年2月

流木形状を考慮した個別要素法による橋脚周辺の流木集積過程に関する数値実験
NUMERICAL EXPERIMENTS ON ACCUMULATION PROCESS OF
DRIFTWOODS AROUND PIERS BY USING A DEM-FLOW COUPLING MODEL

清水義彦1・長田健吾2
Yoshihiko SHIMIZU and Kengo OSADA
1正会員

博士（工学） 群馬大学助教授 工学部建設工学科（〒376-8515 桐生市天神町1-5-1）
2正会員 博士（工学） 中央大学研究開発機構助手（〒112-8551 文京区春日1-13-27）

A numerical simulation method has been developed for describing the behavior of driftwoods motion and theiraccumulation process in channels with bridge piers. The present model has been composed of Eulerian analysis ofdepth-averaged flow motion and of Lagrangean analysis of driftwoods motion by using the extended Distinct
Element Method(DEM), which can simulate the motion of driftwoods as the rigid body. The calculation results show
time-depending process of driftwood accumulation in front of bridge piers including the effect of branch shape ofdriftwoods and also show the amount of dam-up water level due to jamming of driftwoods.

Key Words : Driftwood, Distinct Element Method, dam-up, depth-averaged flow analysis,numerical simulation, the motion of driftwoods

６つの発電用ダム群に大量の流木が捕捉されたものの，流木の引っかかりやすいトラス橋や桁橋では痕跡としても大量の流木が残存したことから，橋の流失・破壊に大きな影響を与えたことが示唆された4)．
最近の豪雨災害において，流木発生に伴う河道災害がこのように，流木群の堰止めによる橋梁被災や水位上しばしば報告されている．例えば，2003年台風10号では，昇による氾濫流発生など，流木の流動・集積は，上流域
北海道胆振地方沙流川流域の山間部において多くの流木での代表的な河川災害を生む要因の一つとなっている．
が発生し，二風谷ダムで1年間にダムで処理される流木その解明のため，これまで主に水理実験的なアプローチの約67年分が捕捉されたが，その下流においても河道内から検討されてきたが（例えば,5）,6),7)），様々な条や橋梁橋脚等に多くの流木が集積した1）．鈴木・渡邊の1）件設定のもとに系統だった実験を行うことはかなりの労研究 によれば，上流域も含め河道内や河道付近から流力が必要とされる．一方，数値解析的なアプローチから出したと推定される流木は全体の約20％にもなることがの検討は数少なく，中川ら8），後藤ら9）の研究が代表的推定されている．上流域においては渓岸崩壊・浸食が生と言える．中川らは，堤内地への氾濫流による流木の流ずることにより土砂とともに多くの流木が発生するが，動・拡散，それに伴う流木群による家屋の破壊といった中流域においても，近年，多くの礫床河川で河道内樹林2）問題を対象に，平面2次元流解析から求められた流れ場化が報告されており ，河道内樹林も流木発生源となりにおいて，個々の流木の重心に関する並進運動と回転運得ることがこうした災害調査からも改めて認識された．
動の方程式を解きながら，流木をラグランジュ的に追跡
とくに，河岸侵食によって供給された流木では根，枝がしてその流動過程を求めた（ただし，流木の重心で2等大きな損壊なく残っていることが多く，そうした形状特分した２つの要素のそれぞれの重心に作用する流体力で性が構造物あるいは樹林帯前面での流木捕捉を促進させ回転運動を記述し，衝突・反発といった流木間の干渉はている．さらに，2004年福井豪雨では，足羽川橋梁被災無視している）．一方，後藤らは，橋梁等の河川構造物に対し，流木の捕捉が大きな影響を与えたことが指摘さ3)による流木群の捕捉・堰き止めを対象に，鉛直2次元流れているし ，2005年台風14号でも宮崎県耳川上流域で

１． 緒言- 829 -

場での流木の挙動を，流れとともに粒子法にもとづく数値解析手法（MPS法に剛体連結モデル10）を組み込んだラグランジュ的解析）を構築し，橋梁による流木の連鎖的捕捉，それによる水位上昇と橋梁上の越流といった一連の過程を数値的に表現している．いずれの研究も流木要素の運動をラグランジュ的に追跡する手法から表現するもので（後藤らの研究では，流体の橋脚上部の乗り上げや橋脚下部への潜り込みといった流体運動の扱いそのものも粒子法で記述），水理実験と合わせ，洪水時の流木の挙動を解析する上で有用な手段である．こうした一連の研究を受けて清水ら11)は，蛇行流路における流木群の流動と，橋梁や砂州樹林帯における流木の集積を対象に，流木の衝突を個別要素法で取り込んだ平面流解析の枠組
みで検討した．ここでは， 1本の流木を単純な棒状とし，これを個別要素で構成させた剛体モデルとして表現，各要素に個別要素法12）を適用して作用流体力と他の流木から受けた衝突力のもとに生じる並進，回転運動を求め，流木の流動・集積過程を追跡した．しかしながら，実際の流木にはさまざま形状特性があり，これが集積過程において大きな影響を与えているものと推測される．そこで，本研究では同様の数値解析手法の枠組みにおいて，流木の形状特性を考慮して橋脚周辺に集積する過程を検討する．すなわち，流木が棒状である場合と，枝付形状とした場合での，橋脚での集積のちがいを数値実験から求める．一方，流木が流れから駆動力（抗力の反作用）を受ける反面，流木が流速と相対運動するために流体場
では運動量欠損が生まれる．この運動量欠損は，清水ら
13）が植生要素を取り込んだ流れ場の解析手法と同様，流れの基礎式の中で抵抗項として表現でき，これによって，流木群の捕捉に伴う橋脚前面での水位上昇量を評価することが可能で，これについても併せて検討する．

図-1 個別要素をつなげて構成した流木（2本の枝付きの場合）

dv pi C M  A3 d 3dt   f n sin  ij  f s cos  ij j 

1
C D A2 d 2 (u  u pi ) 2  (v  v pi ) 2 (v  v pi )2

ここに，ρ，σ：水，流木の密度，CM：付加質量係数（=0.5），upi,vpi:流木を構成する個別要素(i)の流下，横断方向の移動速度，u,v:流下，横断方向の流速，CD：
形状抵抗係数，fn,fs：要素i,j間の法線(n)および接線(s)方向の作用力，αij：要素i,jの接触角である．また，個別要素（球形粒子）について，d：要素粒子の径，A2，A3：2次元,3次元形状係数でそれぞれπ/4，π/6である．
図-1に示すように，流木１本を個別要素でつなげ枝付きの棒状体として表現した（流木長10m,枝長さ0.6m）．
本計算では1本の流木で枝は２本として，主幹に対して45度の角度をもち，枝の位置はランダムに発生させた．
さて，流木を構成する個別要素(i)に，他の流木を構成する個別要素(j)が衝突する際に生じる作用力（接触力）の，接触面に対する法線，接線方向の成分（fn,fs）
は，弾性バネ（バネ定数：kn,ks）と粘性ダッシュポット(粘性係数：ηn,ηs）を用いてモデル化される．たとえば，接触面に対する法線方向作用力fnについては，接触してからのΔt時間における法線方向の相対変量(ΔunΔt)に比例する抗力（en）はバネモデルで，法線方向の相対速度に比例する粘性抵抗力（dn）は，粘性ダッシュポットモデルで与えられ，f n  ent  d nt , ent  ent 1  k n un t , d nt   n un(1) 個別要素法による流木モデル
流木を，N個の個別要素をつなげて枝付きの棒状体として表現し（図-1参照），各要素についての運動は，個別要素法を用いて追跡する．個別要素(i)に関する運動の基礎式は式(1)，(2)のように与えられる．
du pi C M  A3 d 3dt   f n cos  ij  f s sin  ij j

 

j

1
 C D A2 d 2 (u  u pi ) 2  (v  v pi ) 2 (u  u pi )
2(1)

(2)

j

２．数値解析手法の枠組み

(3)とする．モデル定数は文献11と同様とした（バネ定数；kn=5.0x105 N/m, ks=1.5x105 N/m，粘性係数；ηn=2.0x103
Ns/m，ηs=6.0x102 Ns/m，計算時間間隔；Δt=0.0005s）．
ところで，それぞれの要素が流体からの駆動力および接触による作用力を受けて，流木の棒状形状からのずれ（変形による相対位置変化）が生じるが，剛体としての扱いから，形状変化を生まないように修正する．すなわ
ち，各要素のずれを座標修正して棒状に保つ（要素間の相対位置を変化させない）ことで，流木としての運動が追跡可能となる．こうした座標修正法による剛体の移動計算は，Koshizukaら10)が提案したもので，MPS法による流体中に落下する固体塊の挙動や波と浮体の相互作用の計算に適用している．Koshizukaに倣い，流木の棒状を保つ座標修正の手続きは以下のようになる．
まず，流木を構成するN個の要素(i)の移動速度(upi)か- 830 -ら，剛体としての流木の重心に関する並進速度ベクトル
( T )と角速度ベクトル( R )を計算する．
 1 N (4)T   u piN i
 1 N (5)R   u pi  qi I i
ここで，I:流木の重心まわりの慣性モ−メント，qiは流
木の重心位置（ rg）に対する要素(i)の位置（ ri ）のずれであり，  1 N (6)qi  ri  rg ， rg   riN iとなる．こうして求めた流木の重心の運動に従って各要素の速度ベクトルを次式から再計算し．
  
(7)u pi  T  qi  R  Q    UQ    VQ t  J    J    J  y  x  N   yM   x
VU  UV  J   J    x x   y y z s  x2   y2 z s   b F   gh  J JJJ
 
2 
 x x
 u2 h  x u2 hJ J 















 y y  y2 
 v 2 h  v 2 hJ J 2 x y  x y   y xy  u vh  u vh
JJ 





(12)



ここに，u,v ; デカルト座標系でのx,y方向流速，M,N;
x,y方向流量フラックス（Mu=uh,Mv=vh,h;水深），zsは基準面からの高さ，乱れ速度相関（u’,v’で構成された項）
さらに，この速度ベクトルを用いて各要素の座標を次式から修正することで，流動する流木の挙動が追跡できる． は1方程式モデルを用いた．
また，底面せん断応力の反変成分は，



(8)
ri t  t   ri t   u pi t(13)
 b   x bx   y by ,  b   x bx   y by(2) 流木の干渉効果を入れた一般化座標系平面流解析

 bx gn 2 u u 2  v 2
h1/ 3,  by gn 2 v u 2  v 2(14)h1 / 3

流れの基礎式は平面２次元流の連続式および運動量輸送式である．後述するように，蛇行流路についても検討となる．ここで，g;重力加速度，n;Manning粗度係数するため，基礎方程式を一般化座標系(ξ，η)で記述すである．Fは，流木と流体との相対速度によって生じるる．ここでは，長田14）が提案したように，流速ベクトル， 形状抵抗の反作用Fw，橋脚の存在による形状抵抗の反作流量ベクトルの反変成分を未知量とした形式で記述され用Fbから構成される．ている．すなわち， x ,  y ,  x ,  yを座標変換のメトリッ
F  Fw  Fb
クス（x：縦断方向座標，y：横断方向座標，添え字に関(15)A m
122Fw  CD k  u  u pk   v  v pk する微分を表す），J を 座標変換のヤコビアン，流速Af k
2ベクトル及び流量フラックスの反変成分を，
1AFb  CDb b u 2  v 2 U   xu   y v , V   xu   y v
2Af(9)Q   x M   y N , Q   x M   y N
ここで，Af:計算格子面積（=ΔxΔy，概略，Δx=5m,とすれば，基礎方程式は次のように書ける．
Δy=1m程度の格子スケ−ル），m:計算格子内流木構成要  h    Uh    Vh 
素数，Ak:流木構成要素による流れの遮蔽面積（=A2d2,d:(10)  0t  J    J    J 流木径）, Ab:橋脚による流れの遮蔽面積（=hD,D:橋脚幅,h:水深），CD,CDb:流木構成要素，橋脚の形状抵抗係
  Q    UQ    VQ   
数で１とした． Fのx,y成分をFx,Fyとして，
t J JJ

 y y   N M   x
 V x   U
UV
  J   J  F    x F x   y Fy ， F    x F x   y F y  x2   y2 z s  x x   y y z s   b F  (11)  gh J  J JJ2   x u2 h  x x
 u2 hJ J 









 y2  y y  v2 h  v 2 hJ J 
2 x y  x y   y xy  u vh  u vh
JJ 









(16)とし，これを流木および橋脚の存在する計算格子で設定する．
(3) 計算領域と計算条件
計算対象の場は直線水路および湾曲流路で，それぞれ４本の橋脚を設けた（橋脚幅1m,径間10m）．橋脚は，流れに対しては前述した形状抵抗としてその影響を取り込むが，一方で，個別要素法から橋脚の設定も行うため- 831 -

ともに示した．枝付き流木の場合，湾曲内岸から外岸方向に流木塊が流れ，これが橋脚前面でトラップされる．

piers

初期

1分後

2分後
図-2 枝付き流木投入時の設定（直線水路の場合）（個別要素を配置して橋梁の形状を与える），流下して3分後きた流木が橋脚と衝突・干渉する事象も計算に反映され3る．流木（密度0.83g/cm ）は径20cm，流木長は10mで投入本数は100本とした．図-2に直線水路における流木投入時の状況を示した(水路幅50m，長さ300m)．すなわち，4分後初期条件として定まった方形枠の中で流木が重なり合わないようにランダムに配置を決めている．計算での水理条件は，水路勾配1/100，流量200m3/s，フル−ド数
Fr=1.08とした．
5分後

３．流木群の流動と橋脚前面での集積過程
図-3に枝なしの流木群，図-4に枝ありの流木群を各々投入した場合の流動と橋脚での集積状況を経過時間とともに示した．流木群は流速・流向に支配されながらも流木間干渉（衝突）も頻繁に起こしながら流下して行く．
枝なし流木の場合では流動過程においても比較的ばらけているが，枝付き流木の場合，流木塊を形成しやすい（2分，3分後の状況）．こうした状態で橋脚前面に捕捉されると，その後，単体としてさらに流下してくる流木群を効率よく捕捉することになる．最終的な集積状況では，枝付き流木が枝なし流木に比べ集積本数を上回るが，河道幅にしめる流木集積範囲も枝付き流木の場合は広い．
図-5には，集積率（集積本数/投入本数）の時間変化を示した．枝付き流木の場合では最大値で集積率0.8にも達し，枝なし流木の4倍程度の値をとっている．さらに，
図-6には計算縦断水位を比較した．ここに，初期水位は流木の集積がゼロの場合で橋脚の抵抗分によって水位上昇している．同図には10分後に流木が集積した場合の計
算縦断水位を流木の枝ありと枝なしで示している．ともにいわゆるdam-up現象が見られるが，枝付き流木での集積率の高さからより高い水位増加分となるとともに上流側影響範囲も広くなることが見積もられた．図-7,8には蛇行水路における流動と橋脚での集積状況を経過時間と6分後7分後8分後9分後10分後- 832 -

図-3 枝なし流木の流動・集積過程（直線流路で4本の橋脚を設置した）

集積率
1

初期

枝なし流木
0.8

枝付き流木

0.6
1分後

0.4
0.2
0

2分後

024

68

10(分)

図-5 集積率の時間変化
3分後

6水位(m)

河床位
初期水位
枝なし流木の10分後水位
枝付き流木の10分後水位

54
4分後

32
1

5分後

0050

100

150

200

250300縦断長さ(m)図-6 計算縦断水位の比較
6分後

その後，橋脚前面で上流からの流木を多く集積させ，外岸側にかけて橋脚前面での流木塊集積が枝なしに比べ顕著となった(図-8における9分後)．

7分後

４．結語
8分後

本研究では，個別要素法と平面２次元流解析を組み合わせた手法より，枝付き，枝なしの流木形状に注目した流木群の流動・集積過程に関する数値実験を試みた．今後は水理実験との比較を行いながら，様々な条件設定に対し数値実験手法としての有効性を示したいと考えている．

9分後

参考文献
10分後

1)鈴木優一・渡邊康玄：沙流川での台風10号における流木の挙動，水工学論文集，第48巻，pp.1633-1638, 2004.
2)清水義彦・小葉竹重機・新船隆行・岡田理志：礫床河川の河図-4 枝付き流木の流動・集積過程

道内樹林化に関する一考察，水工学論文集，第43巻，pp.971-（直線流路で4本の橋脚を設置した）

976,1999.
3)土木学会：平成16年7月北陸豪雨災害調査報告書，pp.206220,2005.

- 833 -

5分後

5分後

6分後

6分後

7分後
7分後

9分後

9分後

図-7 蛇行水路における流木の流動と集積（枝なし）

図-8 蛇行水路における流木の流動と集積（枝あり）

4)土木学会：平成17年台風14号の記録的豪雨による災害の調査

12)Cundall, P.A.: A computer model for simulating progressive largeと減災対策に関する研究， pp.31-40,2006.

scale movements in blocky rock systems, Symp. ISRM, Nancy,

5)福岡捷二・藤田光一・新井田浩：流木類の流下・堆積とそれらの河道設計への利用，土木学会論文集，479号,pp.51-

France, Proc., Vol.2, pp. 129-136, 1971.
13)清水義彦・小葉竹重機・新船隆行：樹林化河道の洪水流況と浮遊砂輸送に関する数値計算，水工学論文集，第44巻，60,1993.
6)奥澤豊：流木の流下と集積に関する研究，河道の水理と河川

pp. 819-824, 2000.

環境に関するシンポジウム論文集，第４回，pp.141-146,1998.

14)長田信寿：一般化座標系を用いた平面2次元非定常流れの数7)松本健作・小葉竹重機・清水義彦・石田和之・近内壽光・I.

値解析，水工学における計算機利用の講習会講義集，土木

Ioakim ：流木塊の橋脚への堆積に関する研究，水工学論文

学会水理委員会基礎水理部会，pp. 61-76, 1999
(2006.9.30受付)

集，第45巻，pp.925-930, 2001.
8)中川一・井上和也・池口正晃：流木群の流動に関する研究(3)，京都大学防災研究所年報，第37巻,B-2, pp.459-473, 1994.
9)後藤仁志・酒井哲郎・林 稔：粒子法による流木群堰止め過程のLagrange解析，水工学論文集，第45巻，pp.919-924, 2001.
10)Koshizuka, S., Nobe, A and Oka, Y.: Numerical analysis ofbreaking waves using the moving particle semi-implicit method, Int.
J. Numer. Meth. Fluids, Vol.26, pp.751-769,1988.
11)清水義彦・長田健吾・高梨智子：個別要素法を用いた流木群の流動と集積に関する平面2次元数値解析，水工学論文集，
第50巻，pp.81-86,2006.

- 834 -