３． 流木モデルの高度化と流木堆積フロント現象の再現（担当：木村一郎）（１）概要
河川領域には多くの樹木が存在する．それらは大規模な出水による土砂崩れや河岸の崩壊などにより流木化し，河道内を流下する．その流下の際に，橋脚の侵食，取水口への堆積，構造物等への集積や堆積による水位上昇，構造物の破壊など多大な影響を与える可能性がある．また，破堤氾濫によって市街地に流木が流入し，家屋に被害を及ぼすことがあるほか，流木が海域まで流下した際に船の航行を妨げる場合もある．また，昨年の北海道水害では，多くの流木が流出し，橋梁などに対する被害を引き起こした点は記憶に新しい（図１）．以上の点より，河川における流木挙動の予測は災害防止や環境問題の観点から非常に重要であるといえる．
本章では，iRIC における三次元流動ソルバーNaysCUBE に搭載されている流木モデルの高度化を実施し，流木に関する諸現象，特に，流木の集積や壁面の衝突などのフロント現象の再現性の向上を試みたものである．

図 1 ペケレベツ川 流木集積状況（２）既往の流木モデルの問題点とその改良
iRIC, NaysCUBE で用いる流木挙動予測モデルでは，流動は三次元モデルにより計算し，流木は表面流速により水面上を二次元的に移動させている．流木は 1 本の棒としてではなく，球体を複数個並べることで表現するモデルになっている．このソルバーを用いて，すでに流木計算の事例が報告されているものの，改良の余地が多く残されている．そこで，まず最初に，既往の NaysCUBE ソルバー上の流木モデルの問題点についてまとめておくことにする．

- 24 -

a) 流木の初期設定の自由度の拡張
既往モデルでは，均一長さ，直径，密度の流木を指定したセルに配置することのみが可能であり，複数のスケールの流木を配置したり，初期配置の状況を細かく制御したりすることができない．そこで，この点を改善するため，流木の初期配置の方法として，「ランダム設定」，「スケールチェンジ設定」，「セル単位設定」の３つを新たに導入し，ユーザーインターフェースを改良した．以下，この３つの設定の詳細について述べる．
・ランダム設定
まずはランダム設定について説明する．この設定方法では流木長，流木径，流木密度それぞれの最大値および最小値を設定することで，その範囲内でランダムに流木条件を設定することができる．図２はランダム設定を使用した際の流木である．流木の色は密度を示しており，赤い程高密度，青い程低密度である．この図より，長さ，太さ，密度が全てランダムになっており，実現象のように流木の形状やスケールにバラツキがある場合を適切に表現できていると考えられる．

図2

ランダム設定・スケールチェンジ設定
次にスケールチェンジ設定について説明する．
この設定方法ではランダム設定と同様に，まず流木長，流木径，流木密度それぞれの最大値および最小値を設定する．次にランダム数を発生させ，その数値を用いることでスケール変化させる．つまり，ランダム数が小さいほど流木長，流木径，流木密度も小さくなり，大きいほど流木長，流木径，流木密度は大きくなる．
図３はスケールチェンジ設定を使用した際の流木である．ランダム設定と同様に流木の色は密度を示しており，赤い程高密度，青い程低密度である．この図より，流木長が長いほど流木径が大きく(太く)，流木密度も赤いため密度も大きいことがわかり，ランダム設定とは異なり，スケールを変更しているため流木は同じような形をしていることが見てわかる．この方法を用いることで，ランダム設定を用いると発生してしまう，長さが短いが径が太い流木，すごく長いが極端に細い流木といった，実際には発生しにくい流木が出現することを防ぐことが可能になった．

図３スケールチェンジ設定- 25 -

・セル単位設定
最後にセル単位設定について説明する．この設定方法は上述した 2 つの方法とは異なり，各セルに対して詳細に流木設定を行うことができる．
図４はこのセル単位設定を使用した際の流木である．
流木長，流木径，流木密度の設定に加えて，そのセルから何本の流木を発生させるかの設定も可能であるため，実際の流木の堆積状況から詳細なデータが揃っている時に有効な方法である．

図４ セル単位設定b) 流木構成粒子の配置に非合理性の解消
既往の流木モデルは 1 本の剛体としての棒ではなく，複数の粒子を並べ，それぞれの動きを計算した
後，速度や回転等を計算し平均することで 1 本の棒の挙動として計算している．そこで流木を構成する粒子数が非常に重要となる．従来モデルは流木スケールが全て均一であるため流木構成粒子数も同様に一定で問題なかった．しかし，流木設定に自由度をもたせると，粒子数設定も変更する必要性が出てくる．流木設定が複雑化することに対して粒子構成数が一定の場合，図５および図６に示すように，粒子間が大きく離れてしまうケース，粒子間が極端に詰まってしまうケースが発生してしまう．粒子間が離れてしまうと，その空白部分での流れの影響を反映させることができず，逆に詰まってしまうと，流れの影響を複数の粒子が重複して計算されてしまい動きが過剰になることが懸念される．

図５ 粒子間：大図６ 粒子間：小- 26 -

上述した問題点を解決するために，各流木を発生させると同時に粒子数を自動で設定するよう，ソルバープログラムを改良した．このプログラムでは図７のように，流木長を流木径で割ることにより粒子数を確定させる．こうすることで重複を無くし，かつ粒子間隔が最も近くなる粒子数とすることができる．流木長を流木径で割り切れなかった場合は切り捨てとしているため，必ず同一の流木を構成する粒子は重複しない設定となっている．

図７ 流木構成粒子数設定概略図c) 壁面境界条件の物理的非合理の解消
図８は現行モデルにおける流木計算後の流木堆積状況を示している．この図より，流木が水の流れていない計算領域外に飛び出していることが確認される．計算領域外は見えない壁で囲まれているという条件であり，これは実際にはありえない結果となっているため，流木挙動を予測するという観点では修正しなくてはならない問題点である．

図８ 計算領域外への侵入
従来の流木挙動モデルの基本概念を図９に示す．これまでのモデルは粒子 1 つ 1 つの粒子の速度を計算した後(①)，その速度を平均し(②)，その速度より次の移動場所を算出する(③)方法を用いていた．境界条件の扱いについては図１０に示す．反発判定は流木の重心が境界を超えた場合に認められ，その時に速度が逆転する．しかし，図１０内右下の図のように流木全体の重心が境界を超えず流木の一部が境界を超えた場合に，反発の判定は認められず速度が弱まるだけになるため，結果的に流木の一部が境界外に飛び出し，そこに留まってしまう．
そこで今回，図１１に示すように，1 つ 1 つの粒子の速度を計算した後(A)，すぐにそれぞれの粒子を移動させ(B)，移動させた粒子の位置および速度を平均化する(C)というモデルに改良した．図１２は改良した場合の境界条件の取り扱いを図示している．この図のように，流木全体の重心が超えていなくても粒子 1 つ 1 つで移動を計算するため，境界を超えた粒子 1 つ 1 つに対して反発計算が行われ，その後平均化されることで計算領域外には流木が残らないモデルとなっている．

- 27 -

図９ 流木挙動モデル概略図：従来モデル

衝突前

反発有り

反発無し

図１０ 境界条件：従来モデル

ABC

図１１ 流木挙動モデル概略図：改良モデル
ステップ 1ステップ 2

図１２ 境界条件：改良モデル- 28 -

d) Two-way モデルの導入
流木は流れの影響を受けて移動するが，逆に流木は流れに対して影響を及ぼす力がある．1 本の流木の影響は小さくとも，流木が複数本集まり塊となると，その周辺の流れ場は大きな影響を受けることが予想される．構造物周辺に流木が集積する今回のような計算において，この効果は無視することはできないと考えられる．しかしながら，既往の NaysCUBE の流木モデルは流れに対する流木の効果が全く反映されない One-Way モデルであったため，改良が必要である．
今回，水から流木の運動への影響に加えて，流木から水の運動に影響をも考慮する Two-way モデルを作成した．流木の抗力は本来流木の形状（今回は円柱型）を考慮してモデル化すべきであるが，非常に煩雑となる．そこで，本モデルでは簡易的に流木を構成する球体それぞれが個別に抗力を発揮して水流に影響を与えるとしてモデル化を行う．水が流木に引き付けられる力は以下の式(１)で表される 6)．

（１）
ここで，対速度，は流木引張係数，ρ は水の密度，は流速の反変成分，はヤコビアン，は粒子速度の反変成分，は粒子の投影面積，は重力場相

は 1 つのグリッド内の粒子数である．

（３）蛇行水路模型を用いた検証実験
本研究で改良を行ってきた流木挙動モデルの妥当性を検証するため，蛇行開水路模型を用いた水理模型実験を実施した．本節ではこの概要を述べる．
a) 実験条件
蛇行実験には図 13 に示すようなサインジェネレイテッドカーブ 1 周期分用いた水路を使用する．
水路
幅 20cm，蛇行長 3.4m，蛇行角 60 度，勾配は x 方向に 1/1000 である．水路下流側の湾曲部に，図 14 に示すように等間隔に 1cm 幅の構造物を配置した．流量は 1.31ℓ/s，下流端水位は 4.5cm とした．
流木模型には長さ 4cm，直径 5mm，密度 0.6g/cm3 の木材を用い，これを 110 本同時に水路幅全体に均等に供給した．
水深および流速は，流木投入前と流木を投入し構造物に堆積した後にそれぞれ測定する．流木の流れへの影響は堆積した上流側で影響が顕著なため，水深は構造物上流側約 50cm の地点において，右岸側，左岸側，水路中央の 3 点で測定した．流速は PIV 測定を行うことで，表面流速分布を求めた．

Flow

図 13 実験水路- 29 -

4cm3cm3cm3cm3cm

図 14 構造物配置b) 実験結果
図 15 は実験における構造物周辺の動画を 1 秒間隔で示したものである．この図から，構造物の間から抜けていく流木があるが，構造物に止められた流木の上流側にさらに流木がスタックしていく現象が確認できる．最初は流れが速く流木が集まりやすい右岸側からスタックが始まるが，徐々に全体にスタックしていることも確認された．
表 1 は流木投入前後の水深を示している．湾曲の関係で右岸と左岸で多少の水深差があるものの，水深変化量は 0.5mm 上昇と全ての地点で一致している．この結果より，流木の抵抗によって水深が上昇していることが確認できた．
図 16 および図 17 はそれぞれ流木投入前と流木投入後の表面流速を PIV 測定によって測定した結果である．データ欠損域は流速が極端に小さく PIV 粒子到達しなかった箇所を表している．
これらの図より，流木投入前は内岸側で表面流速が大きく，外岸側で表面流速が小さい結果となっている．ところが，流木が構造物に堆積した後，流況が一変し，構造物直前の表面流速が極端に減少している．また，流木堆積の影響で流れが内岸側に集中し，外岸側の流速が内岸側よりも極端に減少している．このような結果から，流木堆積は流況の変化に大きく影響することが示された．
次節では，このようなスタック現象や流木抵抗による流況の変化を拡張したモデルにより再現できるかを検討する．

表 1 流木投入前後の水深変化

単位(mm)

投入前

投入後変化量

左岸側

44.9

45.4+0.5

中央4545.5+0.5

右岸側45.746.2+0.5- 30 -

T = 15 [s]T = 20 [s]T = 16 [s]T = 21 [s]T = 17 [s]T = 22 [s]T = 18 [s]T = 23 [s]T = 19 [s]T = 24 [s]図 15 実験結果(1 秒間隔：左上→左下→右上→右下)- 31 -

構造物

構造物

流速 [mm/s]

流速 [mm/s]

Flow

図 16 実験流速(流木投入前)

図 17 実験流速(流木投入後)（４）数値計算による実験との比較と考察
a) 計算条件
前節で述べた水理模型実験と同様の水理条件および流木条件で数値計算を行う．流動計算には iRIC ソ
ルバーNaysCUBE を用いる．このモデルでは乱流モデルは二次非線形 k-ε モデルを使用している．流れの三次元性を考慮するために鉛直層は水深を 10 分割する設定とした．
衝突効果の検証および Two-way 効果の検証を行うため，Case 1 は衝突および Two-way モデルを共に考慮しないケース，Case 2 は衝突を考慮し Two-way モデルを考慮しないケース，Case 3 を衝突と Two-wayモデル共に考慮するケースとした．Case 3 の Two-way モデルに用いる抗力係数は 1.0 とする．

b) 計算結果および比較
図 18 は衝突を無視したケースである Case 1 の流木挙動を示している．また，図 19 は衝突を考慮したケースである Case2 の流木挙動を示しており，カラーコンターは共に流速を示している．凡例を T=36秒に添付し，単位は m/s である．これらの図より，衝突無しのケースでは全くスタックしないことに対して，DEM を加え衝突を考慮したことによりスタック現象が表現されていることがわかった．また，同時刻において衝突を考慮していない流木のほうが先の地点に移動していることがわかる．これは，衝突を考慮しない場合，流木がそれぞれ単体で計算されることで流速の速い地点に集まっていることに対して，流木が塊になる場合，流速の遅い地点にも流木が存在し，それらが衝突の効果で互いに影響しあうことで流木群全体では少し移動速度が低下したことが関係していると考えられる．
衝突を考慮した Case 2 の挙動と実験の挙動を比較すると，数値計算において明らかに流木の捕捉力が高く，構造物間を通過していないことが確認された．これは，流木が 1 本の棒として考えた場合には少し衝突するだけでも流木全体が衝突の影響を受けやすいことに対して，複数の球体で流木を表現した場
合，1 つの球体が衝突しても棒全体への影響は小さく，その小さな反発の効果が結果的に吸着に近い働きとなっている可能性が考えられる．以上より，本研究で改良を行ったモデルは，若干の課題が残るものの，衝突効果を導入することで実験に近いスタック現象を再現することができ，この衝突効果は再現- 32 -

性向上に欠かせない要素であることが確認された．

T = 29 [s]T = 33 [s]T = 30 [s]T = 34 [s]T = 31 [s]T = 35 [s]T = 32 [s]T = 36 [s]

図１８ Case 1 における流木挙動

T = 29 [s]T = 33 [s]T = 30 [s]T = 34 [s]- 33 -T = 31 [s]T = 35 [s]T = 32 [s]T = 36 [s]図 19 Case 2 における流木挙動

図 20 Case２ (Cd=0.0)および Case３(Cd=1.0)における時間平均水深0.10.140.175

図 21 時間平均表面流速: Case２0.20流速(m/s)

図 22 時間平均表面流速: Case３

図 20 は Two-way モデルを考慮していない Case 2 および Two-way モデルを抗力係数 1.0 で考慮した Case3 における構造物上流側地点における水路中心軸上の時間平均水深である．水深が波打っている要因は，水路全体の勾配が x 方向についているため，湾曲が影響し各地点における相対的な勾配が変化していることである．これより，Two-way モデルを適用した Case 3 の水深は，考慮していない Case 2 と比較して大きく，流木抵抗による水深増加が確認できた．しかし，その変化量に着目すると約 0.02~0.04mm と，- 34 -

実験値の 0.5mm と比較して非常に小さい．
図 21，図 22 それぞれ構造物周辺における Case 2 の時間平均表面流速分布，Case 3 の時間平均表面流速分布である．両図の上に示した共通凡例の単位は m/s である．Case 2 は Two-way モデルを適用していないため，流木の影響を受けていない．両図を比較すると，Case 3 において内岸側で流速の低下が確認できる．図 23 は Case 3 における流木堆積場所を示している．この図より流木が水路の中心および内岸側に堆積していることが確認された．この流木の抵抗が影響し，内岸側での流速低下が発生したと考えられる．

図 23 Case 3 における流木堆積場所
流速の大きさは最大約 0.026m/s 低減しており，流速変化の大きい場所は構造物上流側に集中している．
一方，構造物外岸側では流速が Two-way モデルを考慮していない Case 2 に比べて逆に若干大きくなっていることが確認された．これは，内岸側に堆積した流木の影響で流れが外岸側に遷移し，流速増大箇所が生まれたと考えられる．ここで，流速低下量を実験と比較すると，実験おいて構造物周辺では約 0.1m/sの表面流速低下を観測している．
以上の点から Two-way モデルを追加することによって，実験と同様に流木の抵抗を考慮することのできるモデルとなり流木計算の再現性が向上したといえるが，水深変化，流速変化ともに計算結果の変化量が小さいことがわかり，流木集積による流水抵抗が十分に再現されているとはいえない状態であることが確認された．この要因として，今回用いた Two-way モデルは流木を構成する粒子 1 つ 1 つに対して抗力を個別に計算するモデルであるため，流木 1 本としての抵抗が適切に再現できないことが考えられる．この件に関してはシールディング効果を考慮することや流木を粒子で表現することによる体積や質量の変化など考えていく必要性があると考えている．
（５）本章のまとめ
本章では，流木の挙動，特に構造物周辺は壁面周辺のフロント部の挙動を適切に再現できるよう，既
往の NaysCUBE の機能向上を試みたものである．改良したモデルの有用性を示すため，蛇行開水路を用いた水理模型実験を実施し，実験と同条件で数値シミュレーションを行い，両者を比較することで精度を議論した．得られた主な成果を次にまとめる．

・ 従来モデルでは単一の流木設定しかできなかったため，現実的に様々な流木設定ができるように改良した．設定範囲内でランダムにサイズを設定するランダム設定，同じ形状のままサイズ感を変更するスケールチェンジ設定，それぞれのサイズを詳細に設定するセル単位設定が可能になった．
・ 流木は一本の棒としてではなく，複数の球体を直線に並べられることで表現されている．そこ- 35 -

で，流木を構成する粒子数が非常に重要であるが，この粒子数が多いと粒子同士が重複，少ないと粒子間が大きく空いてしまうということが発生してしまい，結果に大きな影響を及ぼすことが懸念されている．そこで，流木長と流木径から自動的に粒子数を決定するプログラムを作成した．
・ 流木が計算領域外に飛び出してしまう問題点が生じていた．そこで，流木を構成する粒子の計算手順を変えることでこの問題を解消した．
・ Two-way モデルを考慮することで，水深の上昇および流速低下が再現できた．しかし，実験と比較するとその効果が小さく，シールディング効果を考慮することや流木を粒子で表現することによる体積や質量の変化など今後さらに改良していく必要性が示唆される．
謝辞
本研究を遂行するにあたり，北海道大学大学院修士課程学生，北園和也君からプログラミング，モデルのデバッグ，実験の実施，数値解析の実施の全般にわたり多大な協力を得た．ここに記して深甚なる謝意を表する．
参考文献
1) 五十里洋行・後藤仁志・角哲也：自然型調節型洪水吐きの流木による閉塞機構に関する計算水理学的研究，水工論文集，第 50 巻，pp.793-798，2006.
2) 清水義彦・長田健吾・高梨智子：個別要素法を用いた流木群の流動と集積に関する平面 2 次元数値解析，水工学論文集，第 50 巻，pp.787-792，2006.
3) 中川一・井上和也・池口正晃：流木群の流動に関する研究(3)，京都大学防災研究所年報，第 37 巻，
B-2，pp.459-473，1994.
4) Homepage of iRIC Project Changing River Science, http://i-ric.org/ja/ 2016.
5) 後藤仁志，数値流砂水理学，森北出版，2004.
6) 清水義彦・長田健吾：流木形状を考慮した個別要素法による橋脚周辺の流木集積過程に関する数値実験，水工学論文集，第 51 巻，pp.829-834，2007.
7) Kitazono, K., Kimura, I., Shimizu, Y. and Kyuka, T. Comuputation on driftwood motions around obstacles
coupling with a three-dimensional flow model. International Association for Hydro-Environment Engineeringand Research, 2016.

- 36 -