2． 流木スタック挙動の三次元性とNaysCUBE型流木モデルの適用性（担当：木村一郎）（1）概要
2017年（平成29年）の九州北部豪雨災害では，大量の流木が流出し，橋梁に堰き止められて橋梁の破壊や水位の堰上げが生じ，また，家屋を流木が直撃して破壊する等の被害が大規模に発生した．図1，
2は流木の流出・堆積が著しい被災地の状況を示したものである．

図1 流木の橋梁上流側のスタックの例（2017年九州北部豪雨災害）

図2 住宅地における流木流出・堆積の例（2017年九州北部豪雨災害）
河川における流木挙動の予測はこのような水災害の軽減や河川環境等の観点から非常に重要な課題である．特に洪水時の流木の集積や堆積は，構造物の破壊や機能低下，上流河川の水位上昇や周辺河床の浸食など様々な影響を及ぼす．しかしながら，流木の動態モデルの研究は，水流や土砂のそれと比べて― 74 ―

大きく遅れをとっている．これは，流木の形状の複雑さなどに伴う動的挙動の複雑さが原因であると考えられるが，それとともに，深刻な流木災害の頻発は，最近になって顕著になった事象であるという点も，影響していると考えられる．この背景には，流木の流出量そのものが最近増加傾向にあるという点が挙げられ，この原因としては，流木流出を促すゲリラ豪雨の頻発，燃料としての木材伐採の激減に伴う，山地の樹木量の増加などが考えられよう．

図3 橋梁による流木の三次元的堆積（ペケレベツ川，2016年北海道豪雨災害）
流木のモデル化にあたり，水流中を移動する流木の挙動は浮力の作用を受けるため，主に水面付近に限定された二次元性が卓越した現象と考えるのが自然である．しかし，水面の流動は第一種二次流，第二種二次流などの流れの三次元性が関与しているため，間接的に流木挙動も流れの三次元性の影響を受ける．また，流木が橋脚などの構造物に捕捉され停止する過程においては，水理条件によって流木の沈降などが生じ，流木挙動そのものが強い三次元性を示すことが知られている．このような現象を再現するモデルを構築することは流木災害の予測や減災対策の立案において極めて重要となる．
河川流のシミュレーションにおいて，近年は水深積分型平面二次元モデルを用いるのが主流であるが，計算機の発達とともに三次元モデルの適用も普及しつつある．一方，流木挙動モデルとしては，流木を水面付近に限定すると仮定する二次元モデルと，スタック時の流木の潜り込みや水面への突き出しなどを考慮する三次元モデルが考えられる．木村らは，三次元（3D）流動モデルと二次元（2D）流木挙動モデルをカップリングした3D-2Dモデルを構築し，二次流の影響を強く受ける蛇行場での流木挙動の再現や，ダムゲート周辺の流木のスタック現象等の再現を行い，その有用性を示した．しかしながら，図3に示すような橋梁などの構造物への流木の堆積状況（スタック現象）をみると，流木の挙動は時として強い三次元性を示すことが理解できる．このような場合，3D-2D型モデルの再現性は低下することが容易に予想される．
本研究では，前述の背景を踏まえ，障害物等へのスタック時の流木挙動の三次元性の特性の解明し，この三次元性を支配する無次元パラメータを抽出するための水理模型実験を実施した．次に，iRIC上の三次元ソルバーであるNaysCUBEの有する3D-2D型流木動態モデルによる本実験の再現計算を試み，本モデルの適用範囲について検討した．

― 75 ―

（2）流木リチャードソン数の定義
流木は一般に浮力により表面付近を移動すると考えられる．NaysCUBEで用いられている3D-2D型流木追跡モデルの妥当性の根拠もこの点にある．しなしながら，構造物周辺においては必ずしも流木挙動は二次元的とはならない．浮力に対して慣性力が大きい場合，構造物に衝突した際に流木が浮力に逆らって水中に沈み込み，他の流木の下を抜けて構造物をすり抜けていくことがある．この現象を検討するため，浮力と慣性力の比を表すリチャードソン数を流木運動の指標とすることを考える．リチャードソン数は密度成層の安定性の議論などに良く用いられる無次元数であり，具体的には式(1)のように定義される．
(1)ここに，ܴ௜ は密度流リチャードソン数，ߚは体膨張係数，߂ߠは温度差，ℎは水深，‫ݑ‬はバルク流速である．
体膨張係数と温度差は以下の式のように表現することができる．

(2)
(3)
ここで， ܸは体積， ߩは流体密度， ߂ߩは流木と水の密度差である．式(3)を式(1)に代入することで，式(4)のような流木リチャードソン数ܴ௜௧ を定義する．

(4)ここに，‫ܤ‬は水路幅，ܳは流量，߂ߩは流木と水の密度差である．この値が大きいほど浮力が優位であり，小さいほど慣性力が優位となる．

流木リチャードソン数の導入は，固液混相流である流木混入流れを巨視的に捉え，密度差を有する単一流体の運動とみなしてその安定性を検討することを意味する．

図4 実験に用いた水路の写真― 76 ―

図5 実験に用いた水路の平面図

図6 障害物の配置の様子（下流側蛇行頂部）

図7 流木の二次元的堆積の例（Rit= 11.8, 流れは左から右）

図8

流木の三次元的堆積の例（Rit= 0.81, 流れは奥から手前）― 77 ―

（3）実験の概要
実験には，北海道大学水工学実験室に設置された図4，5に示すようなアクリル製蛇行開水路実験装置を用いた．この水路の下流側の湾曲部頂点に，図6のように等間隔に1cm幅の構造物を配置した．流量2.13ℓ/sを与え，下流端水位を調節することにより流木リチャードソン数を変化させる．水路の平面形状はサインジェネレイテッドカーブであり，1.5周期分の平面形状を有する．水路幅20cm，蛇行長3.4m，蛇行角60度である．流量は2.13ℓ/s，下流端水位は3cmとし，フルード数は約0.7である．流木モデルには長さ4cm，直径5mm，密度0.6g/cm3の円柱形状木材を用い，これを水路上流端から110本同時に水路幅全体に均等に分布するよう供給し，構造物上流側において堆積する流木数をカウントした．また，堆積した流木数の全供給流木数に対する比率を捕捉率，堆積した流木のうち水面下に沈みこむものの比率を沈下率と定義した．図7，8に実験における流木捕捉状況の例を示す．図7は二次元的堆積，図8は三次元的な堆積が顕著な例である．流木の沈下判定については，図9に示すように流木が少しでも水面に触れているものについては沈下していないと判断し，流木全体が水中にあるものについてのみ沈下流木と判断した．

図9 流木沈下判断基準に関する説明（○が沈降流木）

図10 実験における沈降率と流木リチャードソン数の関係（水色の縦帯は3回行った実験の値の範囲を示す．
）（4）実験結果および考察
実験は流木リチャードソン数を0.8から11の間で変化させた8通りの条件で実施した．また，現象の再現性を確認するため各ケースで3回ずつ同条件で行い，最大値，最小値，及びその平均値を求めた．
図10は沈下率と流木リチャードソン数の関係を示したものである．流木リチャードソン数が大きいほど沈下率が小さくなることが確認された．今回用いた流木モデルは前述のように流動を三次元計算し，その水面流速により流木を平面二次元的に移動させている．しかし，今回の実験結果では流木リチ― 78 ―ャードソン数が小さい場合，構造物周辺の流木挙動にも三次元性が強く出現し，3D-2Dモデルの再現性が低下することが予想される．図10より流木リチャードソン数が10以上の場合，このような三次元性の影響は極めて小さいことがわかる．このため，流木を水面付近で平面的に移動させる3D-2D型モデルの適用範囲は��� > 10程度と考えられる．

図11 実験における捕捉率と流木リチャードソン数の関係（水色の縦帯は3回行った実験の値の範囲を示す）

図11は流木捕捉率と流木リチャードソン数の関係を示したものである．流木リチャードソン数が大きくなるほど捕捉率が上昇している．この原因としては，流木リチャードソン数が小さいと流木堆積が三次元的となり，捕捉され停止した流木の下をすり抜ける流木が増加するが，流木リチャードソン数が大きい場合はこの「すり抜け現象」が生じにくく，流木が前の流木の上流側へ順に二次元的に堆積するためと考えられる．また，捕捉率の実験値には図10の沈降率と比較してより大きなばらつきがみられる．
この理由として，流木堆積プロセスの逐次性が関与していると考えられる．すなわち，構造物に対して最初に捕捉された流木をきっかけに流木が次々と堆積していく現象が確認されたが，最初の流木がいつどのような形でどこの構造物にスタックするかという極めて不確定性の高い現象が，さらに上流から流れてくる流木の堆積に大きな影響を与え，これが流木捕捉を決定づけると推測される．このため値のばらつきが大きくなったと考えられる．
（5）数値解析結果とその考察
NaysCUBEの流木モジュールを用いて障害物を有する流木挙動の実験結果に関する再現計算を試み
る1-4)．流木同士の衝突は，個別要素法（DEM：Discrete Element Method）を用いる．この衝突における力の評価は，NaysCUBEの流動計算に用いる一般曲線座標系上で行い，力や流速などのベクトル量は反変成分を用いて計算する．一般座標上のDEMの定式化は，次の通りである5)．
図12のように，二次元一般曲線座標ξ-η系で流木を構成する球体iと球体jの衝突を考える．それぞれの球体の中心位置ベクトル࢘௜ ǡ ࢘࢐ は次の通りである．
కఎ‫ܚ‬௜ ൌ ‫ݔ‬௜ ‫܍‬క ൅ ‫ݔ‬௜ ‫܍‬ఎ ， ‫ܚ‬௝ ൌ ‫ݔ‬௝క ‫܍‬క ൅ ‫ݔ‬௝ఎ ‫܍‬ఎ (5)ここに， ‫܍‬ஞ ǡ ‫܍‬ఎ は，共変基底ベクトルである．また，球体iの中心からjの中心に向かう相対位置ベクトル‫ۺ‬୧୨ とそのノルム|‫ۺ‬୧୨ |（＝中心間距離）は次のように表される．

― 79 ―

図12 一般曲線座標上の球体の衝突���� � �� � �� (6)��������������� � = ���� � �� � ��� + ��� � �� � ��� + 2 ��� � �� � ��� � �� � ���(7)ここに， ݃కక などは，座標返還に伴う接続の係数である．いま，時間߂‫ݐ‬間に相対位置ベクトルが次のように変化した場合を考える．

ここに，కఎ‫ۺ‬௜௝ ՜ൌ ‫ۺ‬௜௝ ൅ ȟ‫ ܒܑۺ‬， ȟ‫ ܒܑۺ‬ൌ ȟ‫ݔ‬௜௝ ‫܍‬క ൅ ȟ‫ݔ‬௜௝ ‫܍‬ఎ (8)�����̅�� � �̅�� � �̅�����������， ��̅�� � �̅�� � �̅���

�����

(9)

�̅�� � �� � �� ， �̅�� � �� � �� (10)である．また，二球体間の相対速度ベクトルの反変成分は，次のように表される．
���

��� = �� � �� =�������

�

�

����

， ��� = �� � �� = �� (11)��

球体間の法線方向および接線方向の相対変位（それぞれ，ΔNij, ΔTijとする）は次のように表される．
�������� � �� ��̅�� ��� + �� ��̅�� ��� + ��� ��̅�� + R� ��̅�� � ��� (12)
�

ここに，

�

���� = ���� ���� ��̅�� � �� ��̅�� � (13)�� =���� ���， �� =�� ������� ������� �(14)����� = ���� ��� � ����(15)である．また，球体間の相対速度の法線方向成分vijn，接線方向成vijt分は次のようになる．
�������� � �� �̅�� ��� � ��� �̅�� � �� �̅�� ���� � �� �̅�� ��� (16)������ = ���� ���̅�� �� � �̅�� �� � （17）これらの関係を用いて，次のように法線方向，および接線方向の接触力を通常のDEMと同様に評価する．

― 80 ―

法線方向反発力：�� � �� ��� � �� �t�(18)݁௡ ሺ‫ݐ‬ሻ ൌ ݁௡ ሺ‫ ݐ‬െ ȟ‫ݐ‬ሻ െ ݇௡ ‫ ڄ‬ȟܰ௜௝ (19)������ ��� � ��� ⋅ �� (20)接線方向せん断力：(21)�� � �� ��� � �� �t�݁௧ ሺ‫ݐ‬ሻ ൌ ݁௧ ሺ‫ ݐ‬െ ȟ‫ݐ‬ሻ െ ݇௧ ‫ ڄ‬ȟܶ௜௝(22)������ ��� � ��� ⋅ �� （23)ここに，kn, cn：法線方向のばね定数およびダッシュポット係数，kt, ct：接線方向のばね定数およびダッシュポット係数である．ただし，球体間の引っ張り力は無視し，次のような関係を課す．
�� ��� � �� ��� � 0��������������� ��� < 0(24)また，せん断力が最大摩擦力を越えるとスライダーを稼働させ，次のようにせん断力を計算しなおす．
�� ��� � � � ������� ���� ��� ���������|�� ���| � ��� ��� (25)これらの計算を接触する全ての球体について行い，一般曲線座標の反変成分に変換したのち足し合わせることで，球体iに作用する接触力を次のように求める．
క‫ܨ‬௣ூே்೔ ൌ ෍ ቆ‫ܨ‬௡ ܴక െ ‫ܨ‬௧ఎ௝‫ܨ‬௣ூே்೔ ൌ ෍ ቆ‫ܨ‬௡ ܴఎ ൅ ‫ܨ‬௧௝ܴక ݃కఎ ൅ ܴఎ ݃ఎఎඥࢍ૛ࡰܴక ݃కక ൅ ܴఎ ݃ఎకඥ݃૛ࡰቇ (26)ቇ (27)

障害物に対する流木の衝突のモデル化としては，障害物に固定球体を連続配置して，DEMをそのまま適用する方法も考えられるが，衝突を考慮する球体数が増えるため，計算時間が増加する．北園ら3,4)
はより簡易的に，障害物を計算セル単位で配置し，障害物セルへ侵入した球体を反発係数を考慮して反射させる方法を提案している．本研究もこの方法を用いる．なお，反発係数は1（完全反発）とするが，この値の妥当性については，さらに検討が必要である．
本モデルを用いて実験と同条件の8通りの流木リチャードソン数について数値解析を実施した．
図13は数値解析における流木の障害物周辺の堆積状況を示したものであり，コンターは流速分布を示す．流木が障害物周辺に二次元的に堆積している状況が確認できる．流木リチャードソン数を実験と同様に変化させた場合でも，本モデルは3D-2Dモデルで再現される流木堆積は常に二次元的であり，したがって沈下率は常に0である．一方，計算結果から算出された捕捉率については，流木リチャードソン数の大きさに関わらず0.9程度以上の値となった．実験結果との整合性を考えると，流木捕捉率予測においても，3D-2Dモデルの適用範囲は流木リチャードソン数が10程度以上であることが指摘された．

― 81 ―

Average velocity magnitude (m/s)
0.130.15
0.170.19

図13 数値解析結果における障害物周辺の流木の堆積状況（6）実河川における流木リチャードソン数（石狩川水系漁川におけるケーススタディ）
これまでの考察で，流木リチャードソン数が10を境として流木捕捉形態が異なることが示唆された．
この知見を実務に生かすには実河川の流木リチャードソン数について知る必要がある．そこで，災害時の流木リチャードソン数について，2014年9月に発生した北海道石狩川水系漁川の漁川ダム上流側の事象を対象にケーススタディを行う．この災害では，出水に伴い大量の流木が流出し，漁川ダム上流側に数百メートルにわたって堆積した．また，多くの道路や橋梁も被災した．ダム上流側の流木堆積域にほど近い白扇水位観測所における流木リチャードソン数は，ピーク前後の観測流量Q=200m3/s，ピーク時流路幅B=30m，推定ピーク水深h=2.5mを用い，流木密度を700kg/m3と仮定すると，ほぼ1.0程度と見積もられる．このことから，この付近の流木堆積は三次元性が顕著であったことが予想され，3D-2D型モデルの適用範囲外であったと考えられる．一方，流木の一部はダム湖内に流入した．ダム湖内で水深の急激な上昇と流速低下が生じる．このため，当該ダムのスケールから流木リチャードソン数の値はダム湖流入とともに1オーダー以上，上昇すると見積もられた．したがって，ダム湖内の流木挙動の予測においては本3D-2Dモデルが有効であると推定される．
謝辞
本研究の遂行にあたり，中部電力株式会社，北園和也氏より多くの協力をいただいた．ここに記して深甚なる謝意を表する．
参考文献（第2章分）
1)

清水義彦・長田健吾：流木形状を考慮した個別要素法による橋脚周辺の流木集積過程に関する数値実験，水工学論文集，第51巻，pp.829-834，2007

2)

初田直彦・赤堀良介・清水康行：蛇行流路の流体場と流木の挙動に関する実験と数値解析，応用力
学論文集，Vol.15，pp.415-422，20123)

北園和也，木村一郎，清水康行：三次元流動モデルと流木モデルのカップリングによる構造物周辺
流木挙動解析，土木学会論文集（A2, 応用力学），2017, 印刷中.

4)

Kitazono, K. & Kimura, I. & Shimizu, Y & Kyuka, T. 2016. Comuputation on driftwood motions aroundobstacles coupling with a three-dimensional flow model, Proceedings of19th Congress of the Asia Pacific
Division of the International Association for Hydro Environment Engineering & Research, Colombo, Srilanka.

― 82 ―5)

Homepage iRIC ”Changing River Science”, http://i-ric.org/ja/ 2016.

6)木村一郎・細田尚・音田慎一郎：非線形k-εモデルによる斜め越流型水制周辺の三次元流況に関する数値解析，土木学会応用力学論文集，Vol.5, pp.717-726, 2002.

7)

I. Kimura, Wim S. J. Uijttewaal, T. Hosoda and Md. Shahjahan, Ali：URANS Computations of Shallow GridTurbulence，Journal of Hydraulic Engineering, ASCE, Vol.135, No.2, pp.118-131，2009.2.

8)

Kimura, I. & Hosoda, T. 2003. A non-linear k-ε model with reliability for prediction of flows around bluffbodies. International Journal for Numerical Methods in Fluids, (DOI:10.1002/fld.540), pp.813-837.

9)越塚誠一：数値流体力学，培風館，pp178-181，1997．

10) 後藤仁志.，数値流砂水理学，森北出版，2004.

― 83 ―