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

個別要素法を用いた流木群の流動と集積に関する平面2次元数値解析
NUMERICAL SIMULATION OF THE DRIFTWOODS BEHAVIORBY USING A DEM-FLOW COUPLING MODEL

清水義彦1・長田健吾2・高梨智子3
Yoshihiko SHIMIZU, Kengo OSADA and Tomoko TAKANASHI
1正会員

博士（工学）

群馬大学助教授

2学生会員

修士（工学）
3学生会員

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

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

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

A numerical simulation method has been developed for describing the behavior of driftwoods motion and theiraccumulation process in meandering channel with bridge piers and the forested gravel bar. The present model hasbeen composed of Eulerian analysis of fluid motion given by depth-averaged flow analysis and of Lagragian analysis
of driftwoods motion by using the extended Distinct Element Method(DEM), which can simulate the motion ofdriftwood as the rigid body. The calculation results show time-depended process of driftwood accumulation in frontof bridge piers including the effect of the forested bar and also show the amount of dam-up water level due tojamming of driftwoods.

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

１． 緒言

流木群の堰止めによる橋梁被災や水位上昇による氾濫流の発生など，流木の流動・集積は，上流域での代表的な河川災害を生む要因の一つとなっている．例えば，2003年台風10号では，北海道胆振地方沙流川流域の山間部において大量の流木が発生し，二風谷ダムで1年間にダムで処理される流木の約67年分が捕捉されたが，その下流においても河道内や橋梁橋脚等に多くの流木が集積した1）．鈴木・渡邊の研究1）によれば，上流域も含め河道内や河道付近から流出したと推定される流木は全体の約20％にもなることが推定されている．また，痕跡調査から，流木の河道内発生箇所として，湾曲部内岸側の高水敷・低水路沿いと河岸侵食の発生している外岸側水衝部が指摘されている．近年，いくつかの礫床河川で顕著な河道内樹林化が報告されており2），河道内樹林が流木生産源となることが，こうした災害調査からも改めて認識された．
一方，樹林帯による流木捕捉も洪水痕跡調査からしばしば指摘されることで1），3），4），流木群の流動と，樹林帯による流木群の捕捉，さらに，その集積による樹林帯の破壊（流木発生源）と言ったプロセスの把握が，河道内樹木管理の面からも必要である．
ところで，流木の流動や橋脚による集積に関する水理
実験的な検討は多くの研究者らによってなされているものの（例えば,4）,5),6)），数値解析の面から検討したものは，中川ら7），後藤ら8）の研究が代表的である．中川らは，堤内地への氾濫流による流木の流動・拡散，それに伴う流木群による家屋の破壊といった問題を対象に，平面流2次元流解析から求められた流れ場において，個々の流木の重心に関する並進運動と回転運動の方程式を解きながら（ただし，流木の重心で2等分した２つの要素のそれぞれの重心に作用する流体力で回転運動を記述し，衝突・反発といった流木どうしの干渉は無視している），流木をラグランジェ的に追跡して，その流動過程を求めた．一方，後藤らは，橋梁等の河川構造物による流木群の捕捉・堰き止めを対象に，鉛直2次元流場での流木の挙動を，流れとともに粒子法にもとづく数値解析手法（MPS法に剛体連結モデル9）を組み込んだラグラ
ンジェ的解析）を構築し，橋梁による流木の連鎖的捕捉，それによる水位上昇と橋梁上の越流といった一連の過程- 787 -

を数値的に表現している．いずれの研究も流木要素の運動をラグランジェ的に追跡する手法（ただし，後藤らでは，流体の橋脚上部の乗り上げや橋脚下部への潜り込みといった流体運動の扱いそのものも粒子法から表現している）から表現するもので，水理実験，洪水痕跡調査と合わせ，洪水流による流木の流動を把握する有用な手段として期待されるものと言える．
本研究では，河道内での流木の流動と，橋梁や樹林帯による流木の集積を対象に，流木の衝突による干渉効果を個別要素法で取り込み，平面流解析の枠組みで検討した．すなわち，蛇行流路に橋梁や樹林帯がある場を想定
し，1本の流木を個別要素で連結させた剛体モデル9），10）として表現し，各要素に個別要素法11）を適用して作用流体力と他の流木から受けた衝突力のもとに生じる並進，回転運動を求め，流木の流動・集積過程を追跡した．ここでは，流木は変形しない棒状の剛体として扱っているため，流木一本を構成する各々の個別要素に作用する力に対し，変形を許さないことで生じる内部応力も推定され，木材強度に応じた変形を取り込むことも将来的に可能である．
ところで，流木が流れから駆動力（抗力の反作用）を受ける一方，流木が流速と相対運動するために流体場における運動量が消費される．この運動量欠損は，清水ら
12）が植生要素を取り込んだ流れ場の解析手法と同様で，流れの基礎式の中に抵抗項として取り込むことによって表現した．

⎛σdv pi⎞+ C M ⎟⎟ A3 d 3dt⎠⎝ρ= ∑ − f n sin α ij + f s cos α ij jρ ⎜⎜[](2)j1+ ρC D A2 d 2 (u − u pi ) 2 + (v − v pi ) 2 (v − v pi )2Ipdω pidt=d⋅ ∑ [− f s ]j2 j(3)

式中における，ρ，σ：水，流木の密度，CM：付加質量係数（=0.5），upi,vpi,ωpi:流木を構成する個別要素(i)の流下，横断方向の移動速度および角速度，u,v:流下，横断方向の流速，CD：形状抵抗係数，fn,fs：要素i,j間の法線(n)および接線(s)方向の作用力，αij：要素i,jの接触角である．ここでは，個別要素を球形粒子として与えると，d：要素粒子の径，A2，A3：2次元,3次元形状係数でそれぞれπ/4，π/6，Ip:要素の慣性モーメントである．
流木を構成する個別要素(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

２．数値解析手法

(4)とする．モデル定数は文献13を参考にして流木間衝突での異常反発が生じない条件範囲のもとに試行的に決定し，(1) 個別要素法による流木モデル
バネ定数（kn=5.0x105 N/m, ks=1.5x105 N/m），粘性係数
後藤ら10)は，渓流斜面において風倒木を含む斜面崩壊過程を数値模擬しており，その中で，倒木の運動を(ηn=2.0x103 Ns/m，ηs=6.0x102 Ns/m)，計算時間間隔9)Δt=0.0005sとした．
Koshizukaらが提案した移動剛体モデル を適用した個別
さて，それぞれの要素が流体からの駆動力（形状抵抗要素法から追跡している．後藤らの個別要素法を用いた力の反作用）および接触による作用力を受けて，流木の移動剛体モデルでも，KoshizukaらがMPS法で移動剛体棒状形状からのずれ（変形による相対位置変化）が生じモデルを適用したのと同様に，移動粒子・剛体ともに同10)る．各要素のずれを座標修正して棒状に保つ（要素間のじ個別要素で構成できる利点をもっている ．
本研究でも，同様な枠組みで，流木モデルを構成する． 相対位置を変化させない）ことで，流木としての運動が追跡可能となる．こうした座標修正法による剛体の移動
すなわち，流木をN個の個別要素を1列でつなげて構成計算は，Koshizukaら9)が提案したもので，MPS法によるした棒状体の剛体とし，個別要素についての運動は，個流体中に落下する固体塊の挙動や波と浮体の相互作用の別要素法を用いて追跡した．個別要素(i)に関する運動および回転の基礎式は式(1)から(3)のように与えられる． 計算に適用している．Koshizukaに倣い，流木の棒状を保つ座標修正は以下のようになる．
du⎛σ⎞ρ ⎜⎜ + C M ⎟⎟ A3 d 3 pi
まず，流木を構成するN個の要素(i)の移動速度(upi)かdt⎠⎝ρ(1)
ら，剛体としての流木の重心に関する並進速度ベクトル= ∑ [− f n cos α ij + f s sin α ij ]jrrj( T )と角速度ベクトル( R )を計算する．
r 1 N r1(5)+ ρC D A2 d 2 (u − u pi ) 2 + (v − v pi ) 2 (u − u pi )T = ∑ u pi
2N i- 788 -r 1 N rr(6)R = ∑ (u pi × qi )I ir
ここで，I:流木の重心まわりの慣性モ−メント，qiは流rr
木の重心位置（ rg）に対する要素(i)の位置（ ri ）のずれであり，
r r r
r
1 N r
qi = ri − rg ， rg = ∑ ri
N i(7)となる．こうして求めた流木の重心の運動に従って各要素の速度ベクトルを次式から再計算し．
r r rr(8)u pi = T + qi × R∂ ⎛ Qη ⎞ ∂ ⎛ UQη ⎞ ∂ ⎛ VQη ⎞⎜⎟+⎜⎟⎜⎟+∂t ⎜⎝ J ⎟⎠ ∂ξ ⎜⎝ J ⎟⎠ ∂η ⎜⎝ J ⎟⎠−∂η y ⎞∂η ⎞ N ⎛ ∂η yM ⎛ ∂η x⎟
⎜U+ V x ⎟⎟ − ⎜⎜U+VJ ⎜⎝ ∂ξ∂η ⎠ J ⎝ ∂ξ∂η ⎟⎠⎛ ξ xη x + ξ yη y ∂z s η x2 + η y2 ∂z s ⎞ τ bη F η⎟−= − gh⎜+−⎜JJ∂ξ∂η ⎟⎠ ρJ ρJ⎝ξη ∂η2 ∂+ x x
− u′2 h + x− u′2 hJ ∂ξJ ∂η

()

()()

()

(13)

ξ yη y ∂η y2 ∂− v′ 2 h− v′ 2 h +J ∂ξJ ∂η2η xη y ∂ξ xη y + ξ yη xy ∂さらに，この速度ベクトルを用いて各要素の座標を次式− u ′v′h+− u ′v′h +JJ ∂ξ∂ξから修正することで，流動する流木の挙動が追跡できる．
ここに，u,v ; デカルト座標系でのx,y方向流速，M,N;rrr(9)ri (t + ∆t ) = ri (t ) + u pi ∆tx,y方向流量フラックス（Mu=uh,Mv=vh,h;水深），zsは基準面からの高さ，乱れ速度相関（u’,v’で構成された項）(2) 流木の干渉効果を入れた一般化座標系平面流解析
は1方程式モデルを用いた13）．
また，底面せん断応力の反変成分は，流れの基礎式は平面２次元流の連続式および運動量輸送式である．後述するように，蛇行流路の一部に，橋脚(14)τ bξ ≡ ξ xτ bx + ξ yτ by , τ bη ≡ η xτ bx + η yτ byや樹林化した礫州を配置された状況を計算対象とするため，一般化座標系(ξ，η)で記述された基礎方程式を用
ρgn 2 u u 2 + v 2
ρgn 2 v u 2 + v 2(15)いる．ここでは，長田14）が提案したように，流速ベクトτ bx =, τ by =1/ 3
hh1 / 3
ル，流量ベクトルの反変成分を未知量とする．
ξ x , ξ y , η x , η yを座標変換のメトリックス（x：縦断方となる．ここで，g;重力加速度，n;Manning粗度係数
向座標，y：横断方向座標，添え字に関する微分を表である．Fは，流木と流体との相対速度によって生じる形状抵抗の反作用Fw，砂州・高水敷上の樹木群の抵抗Fｖ，す），J を 座標変換のヤコビアン，流速ベクトル及び流量フラックスの反変成分を，橋脚の存在による形状抵抗の反作用Fbから構成される．
+(U ≡ ξ xu + ξ y v , V ≡ η x u + η y v
ξQ ≡ ξxM + ξ y N , Q ≡ ηxM + η y N

Fw =とすれば，基礎方程式は次のように書ける．
∂ ⎛ h ⎞ ∂ ⎛ Uh ⎞ ∂ ⎛ Vh ⎞⎜ ⎟+⎜⎟+⎜ ⎟=0∂t ⎝ J ⎠ ∂ξ ⎝ J ⎠ ∂η ⎝ J ⎠((11)⎛ ξ x2 + ξ y2 ∂z s ξ xη x + ξ yη y ∂z s ⎞ τ bξ F ξ (12)⎟−+= − gh⎜−⎜ J∂η ⎟⎠ ρJ ρJ∂ξJ⎝
ξ2 ∂ξη ∂+ x− u′2 h + x x− u′2 hJ ∂ξJ ∂η)

()

()

()

ξ y2 ∂ξ yη y ∂− v′ 2 h +− v′ 2 hJ ∂ξJ ∂ηξ xη y + ξ yη xy ∂2ξ xξ y ∂+− u ′v′h +− u ′v′hJJ ∂ξ∂η+()(}(16))

ここで，Af:計算格子面積（=ΔxΔy，概略，Δx=5m,Δy=1m程度の格子スケ−ル），m:計算格子内流木構成要
素数，Ak:流木構成要素による流れの遮蔽面積（=A2d2,d:流木径），λ:樹木の密生度（流体の単位体積に占める樹木の遮蔽面積）,Ab: 橋脚による流れの遮蔽面積（=hD,D:橋脚幅,h:水深），CD,CDv,CDb:流木構成要素，河道内樹木，橋脚の形状抵抗係数で１とした． Fのx,y成分をFx,Fyとして，F ξ ≡ ξ x Fx + ξ y Fy ， F η ≡ η x F x + η y Fy)))(∂ξ y ⎞∂ξ ⎞ N ⎛ ∂ξ yM ⎛ ∂ξ x⎟+V⎜U+ V x ⎟⎟ − ⎜⎜UJ ⎜⎝ ∂ξ∂η ⎟⎠∂η ⎠ J ⎝ ∂ξ({A m1
ρCD k ∑ (u − u pk )2 + (v − v pk )2Af k2

1Fv = ρCDv λh u 2 + v 22
A
1Fb = ρCDb b u 2 + v 2Af2

∂ ⎛ Q ξ ⎞ ∂ ⎛ UQ ξ ⎞ ∂ ⎛ VQ ξ ⎞⎜⎟+⎜⎟+⎜⎟∂t ⎜⎝ J ⎟⎠ ∂ξ ⎜⎝ J ⎟⎠ ∂η ⎜⎝ J ⎟⎠−(F = Fw + Fd + Fb(10)η)(17)とし，これを流木，河道内樹木（ただし，樹木は非水没の条件とした），橋脚のある計算格子で設定する．
(3) 計算領域と計算条件
計算対象の場は，図-1に示すように直線水路から湾曲- 789 -

5分後
FLOW
図-1 計算対象領域(湾曲部中心に橋脚4本設置)
流路に接続する状況で，湾曲の中心に４本の橋脚を設けた（外岸から第１，２，３，４橋脚と呼ぶことにする）．
河道内樹木や橋脚は，流れに対しては前述した形状抵抗としてその影響を取り込むが，一方で，個別要素法から河道内樹木や橋脚の設定も行うため（個別要素を配置して樹木や橋梁の形状を与える），流下してきた流木がこれらと衝突する事象も計算に反映される．流路幅は50m，流木（密度0.83g/cm3）は径20cm，橋脚幅は1m，橋脚径間長は10mとし，流木長は，8m,10m,12mの3種類，投入本数は100本とした．流路流入付近で，投入された流木塊が強制的にばらけるようにして流動過程を追跡した．計算での水理条件は，水路勾配1/100，流量200m3/s，フル−ド数Fr=1.08とした．

6分後

7分後

３．蛇行流路における流木群の流動と集積過程
8分後(1) 橋脚への集積
図-2に，流木塊が投入された後の流動を時間とともに示した．流木は流速・流向に支配されながらも流木間衝突も頻繁に起こしながら流下して行く．湾曲部内岸寄りに流れてきた流木群は内岸を離れ，流路中心に向かうことで橋脚に集積している．図-3は，径間比（流木長/径間長）に対する集積率（集積本数/総投入本数）の変化を示したもので，水理模型による奥澤５）の実験結果も併せて表示した．奥澤の実験では，フル−ド数は0.4〜0.7の範囲にあり水路勾配も緩やかな条件で行っている．一方，本計算では山地河道を意識して水路勾配1/100，フル−ド数Fr=1.08と設定したため，流れの直進性が高く流木長を径間長以下にすると橋脚にほとんど集積しない結果となった．しかし，径間比１を超えると，奥澤５）の実験結果に示されているように，集積率は急増することが計算より確認された．また，図-4,5には，橋脚前面に貯まる流木群に水位を，流木なしでの場合と比較している(初期水位は流木の集積がゼロ，20分後の水位は最終的に流木が集積した時点の水位)．図-4は，最終的に流木が多く集積した第１，２橋脚の橋脚間における縦断水位で，いわゆるdam-up現象が見られる．今後，流木の集積本数や集積パターンの分類を系統的に計算することで，それらと水位上昇分の関係を見積もることもできそうである．
(2)樹林化した内岸砂州と橋脚への集積

9分後

10分後

20分後

図-2 流木の流動・集積過程- 790 -

100

実験1(奥澤)80

実験2(奥澤)5分後

実験3(奥澤)60

実験4(奥澤)40

計算結果(橋脚のみの場合)
計算結果(樹木ありの場合)

20
00

6分後

0.51
1.5 径間比図-3 径間比と集積率の関係

Z(m)
7河床位
初期水位
20分後水位

6橋脚位置
7分後

5432300

320340

360

380

400

8分後

x(m)
図-4 流木の橋脚への集積による水位縦断変化

Z(m)

河床位
初期水位
20分後水位

65.5
59分後

4.5
43.5
301020304050y(m)
図-5 流木の橋脚への集積による水位横断変化
Z(m)
7橋脚位置

10分後

河床位
初期水位

620分後水位

520分後432350400450x(m)
図-7 流木の樹林帯への集積による水位縦断変化- 791 -

図-6 樹林帯が存在する場合の流木の流動・堆積過程

Z(m)
6向の流木の集積が当然のことながら考慮されていない．
これについて，将来的に3次元的な検討が必要であるが，まずは橋脚周辺を含む鉛直２次元場での流木の集積についても検討しながら，その知見を取り込みたいと考えている．そして，流木の集積を踏まえた河道，構造物設計についても検討して行きたいと考えている．

河床位
初期水位
20分後水位

5.5
54.543.5

謝辞：本研究を進めるにあたり，科学研究費補助金（基盤研究(C)(2)，近年の礫床河川の河相変質と植生動態との応答に関する研究，代表：清水義彦）を受けた．

301020304050y(m)
図-8 流木の樹林帯への集積による水位横断変化

参考文献
1) 鈴木優一・渡邊康玄：沙流川での台風10号における流木の挙

１）動，水工学論文集，第48巻，pp.1633-1638, 2004.
鈴木・渡邊の研究 が沙流川での流木集積調査から指2) 清水義彦・小葉竹重機・新船隆行・岡田理志：礫床河川の河摘したように，内岸砂州やその上に繁茂する樹林帯に残道内樹林化に関する一考察，水工学論文集，第43巻，存した流木が多く見られている．ここでは，こうした状
pp.971-976,1999.
況を想定し，図-6に樹木群を伴う内岸砂州を取り込んで
3) 清水義彦・長田健吾：礫床河川における河道内樹林地の洪水
橋脚周辺の流木集積状況を時系列的に示した．また，破壊，河川技術論文集，第8巻，pp.301-306, 2002.
図-7,8には，図-4,5と同様な形で水位変化を示した．こ
4) 福岡捷二・藤田光一・新井田浩：流木類の流下・堆積とそれこで，図-7では，最終的に流木が多く集積した第４橋脚らの河道設計への利用，土木学会論文集，479号,pp.51と内岸の間における水位縦断を示した．図-2,6が示すよ
60,1993.
うに，上流側内岸寄りに流れてきた流木群は湾曲中心で
5) 奥澤豊：流木の流下と集積に関する研究，河道の水理と河川外岸に向かい，内岸寄りの流木数は少なく，樹林化した環境に関するシンポジウム論文集，第４回，pp.141-146,1998.
中州がない場合はほとんど通過してしまう．また，樹木
6) 松本健作・小葉竹重機・清水義彦・石田和之・近内壽光・I.のない中州で，中州上での冠水深が流木径に比べて十分
Ioakim ：流木塊の橋脚への堆積に関する研究，水工学論文
大きい場合でも同様である．しかし，計算結果が示すよ
集，第45巻，pp.925-930, 2001.
うに，樹林化した中州の先端では流木塊の集積が生まれ，
7)
中川一・井上和也・池口正晃：流木群の流動に関する研究
こうした傾向は現地観測でもよく指摘されている．とく(3)，京都大学防災研究所年報，第37巻,B-2, pp.459-473, 1994.に，こうした状況が樹林地の破壊の引き金になることが3)8) 後藤仁志・酒井哲郎・林 稔：粒子法による流木群堰止め過報告されている ．
程のLagrange解析，水工学論文集，第45巻，pp.919-924, 2001.
9) Koshizuka, S., Nobe, A and Oka, Y.: Numerical analysis of

４．結語

breaking waves using the moving particle semi-implicit method, Int.
J. Numer. Meth. Fluids, Vol.26, pp.751-769,1988.

本研究では，平面２次元解析の枠組みで，個別要素法を用い，流木の流動・集積過程を記述するための数値シミュレーションモデルを構築した．ここでは，個別要素を連結して構成した流木の運動を，剛体運動として，流木間の衝突も取り入れた形で表現している．計算対象とした流れ場は，蛇行流路に，橋脚や内岸砂州，砂州上の樹林帯を考慮したもので，それぞれのケースにおいて，流木の流下と集積過程を計算した．とくに，橋脚に集積する流木群については，橋脚間長に対する流木長比が大きくなれば，集積率が増加するなど，既往の水理実験から指摘された特徴を計算でも得られているほか，内岸砂州，河道内樹林化が流木の集積に及ぼす影響も求めている．今回の計算では，流木は1本の棒状としたが，本手法では，容易に流木の形状特性を変えることができ，今後，流木の集積と流木形状との関係を検討したいと考えている．一方，平面2次元で検討しているため，水深方

10) 後藤仁志・原田英治・酒井哲郎：風倒木の渓流斜面崩壊過程に及ぼす影響，水工学論文集，第45巻，pp.721-726,2001.
11) Cundall, P.A.: A computer model for simulating progressive large
scale movements in blocky rock systems, Symp. ISRM, Nancy,France, Proc., Vol.2, pp. 129-136, 1971.
12) 清水義彦・小葉竹重機・新船隆行：樹林化河道の洪水流況と浮遊砂輸送に関する数値計算，水工学論文集，第44巻，pp. 819-824, 2000.
13) 長田健吾・清水義彦・若井明彦：個別要素法を用いた流砂解析における問題点に関する考察，土木学会，応用力学論
文集，vol.7,pp. 1033-1041, 2004.
14) 長田信寿：一般化座標系を用いた平面2次元非定常流れの- 792 -

数値解析，水工学における計算機利用の講習会講義集，土木学会水理委員会基礎水理部会，pp. 61-76, 1999.
（2005.9.30受付）