浮体混在型河川災害の予測・減災モデルの構築とiRICへの実装

京都大学 防災研究所

竹

林

洋

史

北海道大学 工学研究院

岩

崎

理

樹

富山大学 理工学研究部

木

村

一

郎

浮体混在型河川災害の予測・減災モデルの構築とiRICへの実装
Modeling of flows with floating bodies for mitigating disasters and its
implementation in iRIC
木村 一郎1・竹林 洋史2・岩崎 理樹3Ichiro KIMURA, Hiroshi TAKEBAYASHI and Toshiki IWASAKI
1

富山大学 理工学研究部（都市デザイン学）教授
京都大学 防災研究所 准教授
3
北海道大学 工学研究院 准教授
2

要旨

近年の河川災害は，河川水・土砂二相流の枠組を超越し，流木，家屋，車両等の流出，およびそれらの氾濫水への油混入等，浮体を含む事例が急激に増加
し，河川災害の物理機構は複雑化しつつある．本研究は災害予測・減災に資す
る浮体混在型河川流モデルを新規構築しiRICに実装することを目的とする．このため，河川流を浮体相と流体相に分け，浮体相には流木，家屋，油の三つを想定し，流体相には水とともに，状況によっては土砂も含めてモデル化を実施
する．本研究は全体的な研究期間を三年間に設定し，本年はその一年目にあたる．このため，主に基礎的な現象に着目し，力学機構に極力忠実なモデル化を行うことで，基本的なメカニズムの解明と，今後のモデル化の基礎となる土台を構築することに注力した．この方針に即し，三つのテーマの個別テーマについて検討を実施した．一つ目は，浮体として家屋の混入を想定した土石流のモデリングに関するもので，家屋の破壊過程を考慮した平面二次元の土石流数値
シミュレーションモデルをiRICに搭載するとともに，宅地内の家屋の一部が全壊・半壊した2014年の広島市安佐南区八木三丁目で発生した土石流を対象として，家屋の破壊過程を考慮した平面二次元の土石流数値シミュレーションモデ
ルの再現性を検討した結果を示す．二つ目は，浮体として油が混入した河川氾
濫流のモデリングに関するテーマで，河川流や氾濫流に伴う油移動・拡散現象について予測するツールを開発することを最終目的としつつも，まずは計算ス
キーム等の基礎的検討を行い，構築したモデルを単純場，および佐賀豪雨災害
における六角川鉄工所からの油移動拡散過程へ適用して，その性能評価を行う．
三つ目は流木に関するテーマであり，密度の変化により比重が1を超えた場合の流木の挙動をモデル化し，既往の実験結果との比較を通じてモデルの精度を検証する．
《キーワード：iRIC；流木；家屋流出；土石流；氾濫；油流出》― 89 ―

1． 家屋の破壊を考慮した土石流・泥流モデル（担当：竹林洋史）（1）概要
2020年7月に九州から中部地方で発生した豪雨によって954件の土砂災害が発生した1）．これは，1つの台風による土砂災害としては最大の発生件数であった2019年の台風19号による土砂災害の952件とほぼ同じ数の土砂災害発生件数であった．内訳は，土石流等が164件，地すべりが80件，がけ崩れが710件であり，死者17人，家屋の全壊が31棟，半壊が16棟であった1）．
図1は，2020年7月の豪雨によって被災した熊本県芦北町田川地区の家屋である．このように，土石流は非常に大きな流体力を持っており，土石流氾濫域では土石流によって多くの木造家屋が毎年被災している．そのため，土石流・泥流によって被災しやすい場所に建設されている家屋，被災しにくい場所に建設されている家屋を判断できるような数値シミュレーションモデルが求められている．
下流域に宅地などが存在する山地渓流では，土砂災害の警戒区域が各自治体によって設定されている．
土砂災害警戒区域の設定に用いられる土石流の氾濫域は，質量保存則及び運動量保存則に基づいた数値
シミュレーションによって得られた結果は用いられておらず，土石流扇状地の上流端である土石流堆積開始地点を起点とした扇形の土砂氾濫域を想定し，地盤勾配が緩やかになる地点（勾配2度）を土石流
氾濫範囲の下端としている．このような簡便な方法が用いられているのは，土石流渓流が非常に多いことと，このような計算方法でもある一定レベルの精度で氾濫域が予測できるためである．しかし，実際の土石流は，氾濫域の建物や道路の影響を受け，より複雑な形状で氾濫する．
従来，土石流・泥流の数値シミュレーションでは，10m～20m程度の解析格子を用いて解析することが多く，宅地内の家屋や道路の形状を考慮した解析はほとんど実施されていない．そのような中，竹林2）は，2014年広島市で発生した土砂災害に対して，家屋の存在を考慮した土石流の数値シミュレーションを実施し，土石流の氾濫域を再現している．しかし，竹林2）による数値シミュレーションは，家屋を不透過な非破壊構造物として扱っている．
家屋の倒壊については，洪水氾濫や津波などに対してはいくつかの既存研究が存在する．洪水氾濫に関する家屋の倒壊条件については，佐藤ら3）は，洪水氾濫過程の数値シミュレーション結果をもとに，観測された床上浸水深，最大流体力及び被害度の関係から，流体力が2.5m3/s2を超えると住居不可能となる家屋が出現するとしている．高橋4）は，家屋に作用する流体力とモーメントに関する実験結果より，

図1

2020年7月の豪雨によって被災した熊本県芦北町田川地区の家屋― 90 ―

家屋形状は複数の解析格子によって表現する．

土石流

破壊格子

非破壊格子

図2

解析格子と家屋形状との関係

家屋の破壊限界に関する流速と水深の関係式を導き，再現計算結果と被害家屋の調査データがよく一致することを示している．「洪水浸水想定区域図作成マニュアル（第4版）」5）では，氾濫流に対する家屋倒壊等の危険性を表す「洪水時家屋倒壊危険ゾーン」を設定する条件として，木造家屋2階建てのモデル
家屋を想定した荷重条件等より倒壊・滑動・転倒限界の関係式を示している．しかし，土石流・泥流による家屋の破壊プロセスは洪水とはいくつかの点で異なっており，これらの方法をそのまま適用するのは困難である.
本研究では，家屋の破壊過程を考慮した平面二次元の土石流数値シミュレーションモデルをiRICに搭載するとともに，宅地内の家屋の一部が全壊・半壊した2014年の広島市安佐南区八木三丁目で発生した
土石流を対象として，家屋の破壊過程を考慮した平面二次元の土石流数値シミュレーションモデルの再現性を検討した結果を示す．
（2）家屋に作用する応力
津波および洪水流によって家屋に作用する力を算定する場合，家屋をコントロールボリュームとして家屋周辺に作用する応力を積分して得られた力が用いられる．これは，洪水流先端部が家屋に到達した
直後に家屋が破壊されるのでは無く，家屋背面に水が到達してから家屋が流出する状況が被災地でよく観察されるためである．しかし，土石流の場合は土石流先端部が家屋に到達した直後に家屋が破壊され，家屋がその形状を保たずに，家屋の一部損壊が発生することが多い．そのため，家屋より小さいコント
ロールボリュームを用いて得られた応力によって家屋の破壊過程を評価する方が適切と考えられる．本解析の解析格子の大きさは2mであり，家屋の大きさよりも小さい．そこで，本研究では図2に示すように，家屋形状を複数の解析格子で表現し，各解析格子をコントロールボリュームとして作用する力を計算し，家屋の破壊・非破壊の判定に用いた．

― 91 ―

図3

図4

解析領域と地盤高（単位：m）

家屋の破壊状況（（a）と（b）はルーチェサーチによる撮影）― 92 ―

図5

被災前後の航空写真と家屋破壊領域（オレンジ）・家屋非破壊領域（黄）の平面分布

図6

家屋の破壊状況（Case1）

本解析では正方形格子を用いている．そこで，横断方向の流速をゼロ，横断方向圧力は左右で方向は異なるが，同一の値を有しているとし，土石流が流下方向に流れて家屋壁面に垂直に作用している場を考える．木造家屋は可撓性が弱いため，壁はある一定以上の応力が作用すると瞬時に破壊されるとし，家屋破壊時の水深・流速の時間的な変化を無視する．家屋内にはまだ土石流が流れ込んでいないことを考えると，家屋の解析格子に働く単位幅あたりの流下方向の応力Fhxは以下の様に評価できる．
Fhx 1
m gh 2 cos   m hu 2
2

（1）

ここに，gは重力，zbは河床位である．ρmは以下の関係がある．
m    c  ― 93 ―（2）

図7

家屋の破壊が土石流氾濫域に与える影響― 94 ―

単位幅当たりの家屋破壊限界応力は，家屋の構造によって異なるが，本解析では800kN/mを全ての木造家屋に適用した．これは，渓流出口付近の地盤勾配が約10度であることを考慮すると，例えば，土石
流の流速が約10m/s，流動深が約3.5mで家屋が破壊される．
（3）解析条件計算領域は，図3の航空写真内でカラーコンターで示された本川源頭部からJR加治線の2km×0.7kmの領域である．カラーコンターは地盤高を示している．土石流は，斜面崩壊を発端として発生するものとした．斜面崩壊の発生場所は，現地で斜面崩壊が確認されている本川源頭部とした．崩壊の大きさは，
現地調査の結果を参考にして数種類を検討したが，宅地に流入する土石流の特性にはほとんど影響を与えなかった．これは，斜面崩壊による流出土砂量よりも流下しながら河床や河岸から土石流に取り込まれる土砂量の方がかなり多いためである．解析条件は，Case 1として家屋の破壊を考慮した解析，Case 2として家屋の破壊を無視した解析の2種類を実施した．家屋の破壊は木造家屋のみで発生するとし，鉄筋コンクリート構造の県営住宅は破壊しないものとした．家屋が破壊されたと判断された解析格子は，次の計算ステップでは家屋無しの解析格子として扱う．家屋は複数の解析格子で表現されているため，ある家屋の一部の解析格子のみが家屋破壊と判断されるため，半壊・一部破壊家屋などが評価可能である．
（4）家屋の被災状況
図4に被災後の住宅地の様子を示す．渓流の出口付近の宅地には，鉄筋コンクリート構造の県営住宅（図4（c））と多くの木造家屋が存在していた．鉄筋コンクリート構造の建物は，流速・流動深が大きい状態で土石流が衝突しているが，図4（c）に示すように，建物の破壊には至っておらず，外装の一部が損傷しただけとなっている．一方，渓流に沿った領域の多くの木造家屋が全壊・半壊の被害を受けていた．図5に被災前後の航空写真と家屋破壊領域（オレンジ）・家屋非破壊領域（黄）の平面分布を示す．
図に示すように，渓流の出口から斜面下方に沿って設置されている道路周辺で木造家屋の被害が確認できる．また，家屋は全てが全壊しているわけでは無く，家屋の一部のみが破壊されたものも多い．
（5）結果と考察
図6にCase 1で得られた土石流流下後の住宅地での地盤高さの変化量の平面分布を被災後の写真と重ね合わせたものを示す．図中に青色で示されている地点は家屋が破壊されたと判断された地点である．
一方，家屋が存在するにも係わらずコンターが示されていない領域は，家屋が破壊されなかった地点である．図に示すように，渓流出口の直下流に位置する家屋領域Aでは，全て青色となっており，この領域の全ての家屋が全壊したという解析結果となっている．実際，現地においてもこの領域の家屋は全て全壊となっている．さらに下流の家屋領域Bは，一部の家屋が全壊し，半壊家屋も存在している．図4に示した被災後の写真を見ると，この領域は全壊家屋と半壊家屋が存在しており，半壊家屋は計算と同様に，家屋領域の東側（土石流の主流に近い側）で破壊が発生していることがわかる．さらに下流の家
屋領域Cを見ると，全壊と判定された家屋は無いが，一部損壊している家屋が見られる．この領域についても図4に示した被災後の写真から，一部損壊している家屋が見られ，計算結果と一致していることがわかる．このように，本解析モデルによって土石流で全壊・半壊する可能性のある家屋を評価することがある程度可能であると考えられる．
図7にCase 1と家屋の破壊を考慮していないCase 2の住宅地での土石流の深さの時空間的な変化を示
す．Case 2のように，家屋の破壊を無視すると，家屋の上流域に多くに顕著な土砂の堆積域が形成されるため，土石流の氾濫範囲が狭くなることがわかる．また，家屋領域Aが破壊されなかったため，土石流の主流が西側の家屋の少ない領域に移動しており，実現象と大きく異なる結果となった．

― 95 ―

（a）25s（b）50s
図8（1） 土石流の家屋破壊応力の時空間的な変化のiRICインターフェースでの表示（2014年の広島豪雨災害時に広島市安佐北区可部東で発生した土石流）（6）検討結果のまとめ
宅地内の家屋の一部が全壊・半壊した2014年の広島市安佐南区八木三丁目で発生した土石流を対象とし，家屋の破壊過程を考慮した平面二次元の土石流数値シミュレーションを実施し，家屋の存在が土石流の氾濫域に与える影響を検討した．その結果，本数値シミュレーションモデルによって，家屋の半壊・全壊の状態を比較的良く再現でき，土石流で全壊・半壊する可能性のある家屋を評価することがある程― 96 ―

（c）75s（d）100s
図8（2） 土石流の家屋破壊応力の時空間的な変化のiRICインターフェースでの表示（2014年の広島豪雨災害時に広島市安佐北区可部東で発生した土石流）度可能であると考えられる．また，家屋の破壊を無視すると，家屋の上流域の多くに顕著な土砂の堆積域が形成されるため，土石流の流下方向の氾濫範囲が狭くなることが明らかとなった．
（7）Morpho2DHへの家屋破壊応力表示機能の追加
これまで，Morpho2DHでは家屋の破壊を考慮した解析や家屋の破壊応力の表示は不可能であった．

― 97 ―

図9 土石流の最大家屋破壊応力の空間分布のiRICインターフェースでの表示（2014年の広島豪雨災害
時に広島市安佐北区可部東で発生した土石流）

本研究では，iRICのライブラリーファイルを修正してMorpho2DHで家屋の破壊応力を表示できるようにした．以下のその概要を示す．
（8）家屋破壊応力の時空間的な変化及び最大家屋破壊応力の表示
図8は2014年の広島豪雨災害時に広島市安佐北区可部東で発生した土石流の家屋破壊応力の時空間的な変化をiRICインターフェース上で表示したものである．なお，対象の土石流災害は，Morpho2DHのチュートリアルで対象としている現象と同一である．図に示すように，土石流の流動とともに家屋破壊応力の表示域も変化し，家屋破壊応力の時空間変化の解析及び表示が可能となっていることがわかる．
図9は最大家屋破壊応力の空間分布を示す．対象渓流では，勾配が急で流速が大きくなる上流域の渓流で家屋破壊応力が大きくなっていることがわかる．対象とした土石流でも斜面下部で家屋の破壊が実際に発生している．対象地点の最大家屋破壊応力は200kN/mであり，上述の広島市安佐南区八木三丁目の解析で用いた値よりもかなり小さい値となっている．家屋が破壊するかどうかを判定する限界家屋破壊応力については，家屋の構造特性に大きく影響されるため，それらを考慮したる限界家屋破壊応力の設定が必要である．
参考文献
1）国土交通省：令和2年7月豪雨による土砂災害発生状況（9月30日_18:00時点），2020，https://www.
mlit.go.jp/river/sabo/jirei/r2dosha/r2_07gouu_2009301800.pdf.
2）竹林洋史：広島市で発生した土石流の数値シミュレーションの事例と対策，地盤工学会誌，Vol.64，
No.4，pp12-15，2016．
3）佐藤智，今村文彦，首藤伸夫：洪水氾濫の数値計算および家屋被害について－8610号台風による吉田川の場合－，第33回水理講演会論文集，pp.331-336，1989.
4）高橋保，中川一，加納茂紀：洪水氾濫による家屋流失の危険度評価，京都大学防災研究所年報，― 98 ―

No.28，B-2，pp.455-470，1985.
5）国土交通省水管理・国土保全局河川環境課水防企画室・国土技術政策総合研究所河川研究部水害研
究室：洪水浸水想定区域図作成マニュアル（第4版）＜参考資料1＞，pp.viii-xiii，2015．

― 99 ―

2． 洪水氾濫に伴う油流出拡散モデルの開発（担当：岩崎理樹）（1）はじめに山地から川を経由して海にわたる降雨－流出過程の中で，水の流れは様々な物質を輸送する．そのよ
うな中でも工学的問題として重要と考えられるものの一つに，様々な汚染物質の移動・拡散問題がある．
流域において流出，漏洩する様々な汚染物質は，水の流れによって我々の生活圏内に入り込み，健康被害等の直接的な問題の原因になるとともに，食物連鎖を経由した間接的な被害にもつながる．そのため，汚染物質の発生源からどのような形態で移動・拡散するのかを理解し予測することは工学的に非常に重要な課題となる．流域治水という概念の元，新たな河川整備について議論が進む中，流域内の適切な場
所に氾濫水を溢れされるという考え方も議論される一方，居住地や重要な施設等がなくとも，流域には汚染物質となり得るものは多く，汚染物質の流出という観点から新たな負荷が発生する可能性もある．
そのような背景の中，本研究は様々な汚染物質の中でも油の移動・拡散問題を予測する手法を開発しようとするものである．油の拡散問題は，古くより海洋において盛んにおこなわれてきた1, 2, 3）．これは，タンカーや油田からの油流出事故が大きな問題であり，移動・影響範囲を予測する等の工学的要請があるためと考えられる．一方，陸域における油流出現象については，海洋のそれと比してそれほど進んでいないといえる．しかしながら，令和元年8月に発生した佐賀豪雨災害において，氾濫流によって鉄工所より油の流出事故が発生した（図－1）ように，流域内における氾濫リスクの高まりは，このような物質の流出事故を誘発する可能性があり，万一洪水氾濫によって油等の物質が流出した際の被害範囲を予測することは重要な課題であると考えられる．また，陸域の油のソースとなり得るものとして，発電所等における貯蔵タンクや石油等を輸送するパイプラインからの流出が考えられる4）．ロシアなどの石油産出国において，これらの施設は石油や天然ガスを輸送・貯蔵する重要なインフラである反面，これらの故障は大きな環境リスクを生むことになる．特に，極地地域に設置されるパイプラインは近年の地
球温暖化に伴い，パイプラインが埋設されている永久凍土が融解することでその支持力を失い，損傷や崩壊を起こす危険性も懸念される．また，河川における油の流出・拡散は，河床土砂への吸着により非常に長期間の環境リスクとなり得るケースもある5）．流出する可能性のある油等の物質の輸送拡散範囲を理解しておくことは，流域における環境リスクを評価するためにも重要なテーマとなりうる．

図－1

六角川鉄工所からの油流出事故の様子（共同通信社より）― 100 ―

油のような水の流れによって輸送される物質の移動・拡散現象を予測するには，水により生じる移流，
乱流等の乱れによって生じる拡散などを適切に評価するモデルが必要となる．また，このような複雑な現象を予測するには，数値計算モデルの適用が有効であるが，これら移動・拡散過程を精度よく表現できる数値解析手法についても考える必要がある．
本年度は，以上のような背景より，河川流や氾濫流に伴う油移動・拡散現象について予測するツールを開発することを最終目的として，まずは計算スキーム等の基礎的検討を行い，構築したモデルを単純場，六角川鉄工所からの油移動拡散過程へ適用して，その性能評価を行うこととする．
（2）計算モデル
本研究では，iRICに同封されている非定常二次元浅水流モデルを基にした氾濫モデルである
Nays2dFloodに新たに油輸送拡散モデルを組み込むことにより，非海上における油の輸送解析を行った．
以下，本研究で構築したモデルについて述べていく．
（a）油輸送拡散モデル
油は他の物質とは異なり，水に簡単には溶けず，通常，乳化や紫外線の影響など化学的，生物的，物理的な要因により状態変化が発生する．しかし，市街地における油事故の様に流出してから，移動拡散し，事後処理されるまでが比較的短期間に行われることを考えると，上記のような生物的，化学的な影
響は軽微であると考えられる．そうすると，図－2のように水面上で状態変化が発生しないと考えることができ，油は水面上で油厚を持つものと考えられる．

図－2

水面に浮かぶ油の模式図．

そのため，油の輸送を考えていく上では平均油厚さhに関する輸送方程式を用いることが主流である．
油の輸送拡散を数値的に扱う方法として，大きくEuler型とLagrange型のモデルがあるが，本研究ではTkalich
et al.3）を参考に，Euler型の油輸送方程式を採用した．式（1）に直交座標系における輸送方程式を示す．

𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� ℎ 𝜕𝜕𝜕𝜕� ℎ
𝜕𝜕
𝜕𝜕𝜕
𝜕𝜕
𝜕𝜕𝜕
�
�
� �� ��� � �
�� �
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

（1）

ここにh：水表面の油厚さ，ux, uy：油輸送速度，D：拡散係数，Q0：単位時間単位面積あたりに流出す
る油量である．
油の輸送の場合，拡散係数の評価に注意する必要がある．この拡散係数は，乱流拡散といった流れの1）
乱れにより生じる成分に加えて，油自身が持つ油拡散（Spreading）
について考慮している．水表面上に油が浮かんでいる場合，油自身が持つ移動特性が存在し，時間スケールによって，移流卓越領域，重力卓越領域，表面張力卓越領域と支配的な外力が異なり，これらの外力と油－水間の抵抗力のバランスにより，流れがなくとも油は拡散していく．本研究では，特に重力－粘性卓越領域について考え，拡散係数はゼロ方程式モデルを用いた乱流拡散成分と油拡散を足し合わせた次式で評価する．この拡散係数については，のちほどその導出について理論的に考察を加える．

― 101 ―

��𝑔𝑔𝑔� 𝜌𝜌� �𝜌𝜌� � 𝜌𝜌� � 𝜅𝜅� 𝑢𝑢∗ ℎ�𝜌𝜌� 𝑓𝑓6（2）

ここに，g：重力加速度，ρw：水の密度，ρo：油の密度，f：水と油の摩擦に関する係数，κ：カルマン定数，u*：摩擦速度，hs：水深である．
油輸送速度については，通常，流れのモデルから得られる流速と風の影響を足しあわせて評価されるが，ここで扱う流れ場は氾濫流であるため，風の影響は比較的小さいと考えられる．ここでは，簡略化のため，Nays2D Floodから得られる水深平均流速を油輸送速度として計算に用いている．
（b）計算スキーム
本計算のような油の輸送拡散を考える際，与えた物質の量と拡散域は非常に重要な関係がある．そのため系内に存在する物質の量は計算過程で不変である必要があり，かつ可能な限り数値拡散が小さい計
算スキームの構築が重要である．そこで本研究では，式（1）の移流拡散方程式を，𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� ℎ 𝜕𝜕𝜕𝜕� ℎ
�
�
�0
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

（3）

� �� �

𝜕𝜕
𝜕𝜕𝜕
𝜕𝜕
𝜕𝜕𝜕
�� � �
�� �
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

（4）

と二段階分離解法で計算することとし，式（3）の移流項については保存性を保証しつつ数値拡散が小さい解法の一つであるCIP-CSL2法（Constrained interpolation profile-Conservative semi-Lagrangian
scheme）6）を適用し，式（4）の非移流項については，中央差分を用いることとした．以下，数値計算法について詳しく述べる．
移流項の解法として用いるCIP-CSL2は，完全保存保証型CIP法である．この方法では，計算物理量の方向積分値に対して三次関数の補間関数を用いることにより計算点間の物理量を保存し，分布と総量の保存を保証するものである．また，内部の物理量の総量を計算に用いることにより，計算点両端の物理量の値が0の場合であっても，物理量を関数として補間することができる．そのため他の解法と比べ値を過小評価することなく誤差の成長を少なくし安定的，高精度に計算を行うことができる．以下，CIPCSL2の詳細な説明を行う．
式（1）は油厚さhに対する移流方程式であり，通常hについて直接数値的に解くことになるが，CIPCSL2においては，式（5）で示される空間積分値Fを新しく定義する．
௫‫ ܨ‬ൌ න ݄݀‫ݔ‬
ここで，上記Fに対する移流方程式，

௫೔𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� 𝐹𝐹
�
�0
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

（5）

（6）

を考える．これはFに対する保存系の移流方程式であるが，これをxに対して微分しても式（5）の定義
より，以下のhに対する保存系の移流方程式になる．

𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� ℎ
�
�0
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

― 102 ―

（7）

したがって，Fに対する保存系移流方程式を数値的に解き，そのプロファイルからhを求められれば，保存系を保証した計算となる．Fについては，通常のCIP法と同様に，

𝐹𝐹� � �𝜁𝜁 � � �𝜁𝜁 � � ��� 𝜁𝜁

（8）ℎ� � 3𝑎𝑎𝑎𝑎 � � 2�𝑎𝑎 � ℎ��（9）のように三次多項式とおけば，定義よりhについては，として計算することができる．なお，セル内の物理量の総量σは式（10）で定義され，これらの時間変化は式（8）から求まる∆t間にセルから流出入する質量フラックスの差により式（11）のように求める．
����𝜎𝜎� � �

��

���

（10）

𝜎𝜎���� � 𝜎𝜎�� � �� � ����（11）

� � ��� Δ𝑡𝑡（12）

ここに，nは時間ステップを表す．係数a, bについては次式で求める．

𝑖𝑖�� � 𝑖𝑖 � ������ �𝑖𝑖���� � �ℎ�� � ℎ���������𝑖𝑖� � 0𝑢𝑢� � 0𝑖𝑖
𝑖𝑖 � 1Δ𝑥𝑥 ��2𝑠𝑠𝑠𝑠𝑠𝑠�𝑢𝑢� �𝜎𝜎�����
Δ𝑥𝑥 �

�

2ℎ�� � ℎ����

�

Δ𝑥𝑥

�
3𝑠𝑠𝑠𝑠𝑠𝑠�𝑢𝑢� �𝜎𝜎�����
Δ𝑥𝑥 �

（13）
（14）
（15）

（16）

二次元の移流方程式については，上記一次元の解法を利用しフラクショナルステップ法により計算する．そのために，hに対する線積分値，面積分値を以下のように定義する．
����

𝜎𝜎� � �������𝜎𝜎� � �
����

���

��

��

� 𝑑𝑑𝑑𝑑� 𝑑𝑑𝑑𝑑

����

�

��

� 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑

（17）
（18）
（19）

この定義と本研究で用いているStaggered格子系の変数配置を考えると，スカラー変数である油厚をあえて格子点定義し，線，面積分値を上式基づき図－3のように定義する．

― 103 ―

図－3

Staggered格子系における線，面積分値

これは，面積分値がセル内における油の体積を表しているため，セル中央定義とすることが体積保存を計算する面から合理的であるためである．
計算すべき二次元移流方程式は，式（3）であるが，CIP-CSL2を適用する場合は，面積分値に対する保存系移流方程式，

𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� 𝜌𝜌 𝜕𝜕𝜕𝜕� 𝜌𝜌
�
�
�0
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

（20）

を考え，この方程式をフラクショナルステップ法により二段階で計算する．

𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� 𝜌𝜌
�
�0
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� 𝜌𝜌
�
�0
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

（21）
（22）

これらの式に対して空間微分をとることで，𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕� 𝜎𝜎�
�
� 0,
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕� 𝜎𝜎�
�
� 0,
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� ℎ
�
�0
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� ℎ
�
�0
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

（23）
（24）

が得られ，h, σx, σy, ρについて順次CIP-CSL2法を適用して計算する．なお上式からわかるようにそれぞれの式形は全く同じであり，図－4に例として示す通り，式（8）～（16）に示すhとσといった点値と線積分値にかかる更新が，線積分値と面積分値にも成り立つ．すなわち，1次元CIP-CSL2の離散式においてhとσを， 𝜎𝜎� , 𝜎𝜎� とρなどと置き換え，適切な流速を用いるだけで同様の方法により値の更新が可能となる．

― 104 ―

図－4

格子点値，線積分値，面積分値の更新

図－5

二次元CIP-CSL2法の解法手順

簡単のため，一次元のCIP-CSL2法で説明した手順を 𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶1𝐷𝐷�𝑢𝑢� , ℎ� , ℎ��� , 𝜎𝜎 � , 𝜎𝜎 ��� , 𝑖𝑖, 𝑥𝑥�などと表記すると，二次元のCIP-CSL2は図－5の手順によって計算できる．

しかし，上述の計算方法を用いると，油厚に対してオーバーシュートやアンダーシュートが発生する
ことがわかった．これを防ぐために，ここでは有理関数を導入したR-CIP-CSL2を用いることとした．
この場合，使用する三次多項式を以下のように置く．

𝐹𝐹� �𝛽𝛽�2𝑎𝑎𝑎𝑎ℎ�� ��𝑎𝑎𝑎𝑎 � � �𝑎𝑎 � � ��� 𝜁𝜁
1 � �𝑎𝑎

� �𝑎𝑎 � � � 3𝑎𝑎𝑎𝑎 � � 2�𝑎𝑎 � ℎ��
�1 � 𝛽𝛽𝑎𝑎��

（25）
（26）

なお，セル内の積分値やその時間変化については，既に説明した方法と同様である．CIP-CSL2と同様にセル境界の境界条件を基に，係数を決めると以下のようになる．

��

��

1
������ �2 � �����𝑢𝑢� ��Δ𝑥𝑥� � Δ𝑥𝑥����� �1 � �����𝑢𝑢� ��Δ𝑥𝑥�
Δ𝑥𝑥 �

1�������𝑢𝑢� ����� �3 � �����𝑢𝑢� �2�Δ𝑥𝑥� � �����𝑢𝑢� �Δ𝑥𝑥�ℎ��� �1 � �����𝑢𝑢� ��Δ𝑥𝑥 � 2ℎ� �Δ𝑥𝑥 �― 105 ―

（27）
（28）

図－6

��

内田・河原7)による関数補正の考え方

1 �2𝜎𝜎��� � �����𝑢𝑢� �Δ𝑥𝑥 ��� � ���� �𝜎𝜎��� � ���� Δ𝑥𝑥
Δ𝑥𝑥

（29）

これによりCIP-CSL2によるオーバーシュートとアンダーシュートをある程度，防ぐことができると考えられる．
ただし，もう一点注意が必要なのは，格子系によって生じる質量保存の崩れに起因するオーバーシュ
ートとアンダーシュートである．Nays2D FloodではStaggered格子を採用しているため，格子境界に流速成分が定義され（図－3におけるuxi,j, uyi,j），この流速成分とセル中央で定義される水深が，連続式を満たすように決定される．これらの定義点の値をそのまま用いて，セル中央定義の値を計算すれば，連続式で保障されている精度で質量保存が成り立つことになる．しかし，CIP-CSL2の過程の中で，格子点値の𝑢𝑢� 𝑢𝑢���
計算においては，平均した流速（図－5における ���,� ）が使われる部分がある．これらによって計算される格子点値の保存量については，連続式による質量保存は何ら保証されないために，オーバーシュー
トとアンダーシュートが生じることになる．
これを防ぐために，内田・河原と同様の手法を用いる．図－6に関数の値を補正する際の図を示す．図－6に示すように格子点量hの値は，物理的には挟まれる上位関数σから求まる線（面）平均値の間に収まるはずである．図－3からわかるようにhの上位関数の更新のうち，uxi,j, uyi,jを用いて計算される部分については，このような上記質量保存の崩れは生じず，計算スキームによるオーバーシュートを防ぐため

𝑢𝑢� 𝑢𝑢���にR-CIP-CSL2を適用することで，適切な値が計算されるはずである．そこで，���,� を用いたCIP-CSL2で計算された格子点値が，R-CIP-CSL2により計算された上位関数から計算される値の範囲に収まっていない場合，図－6のようにhの値を線（面）平均値の値に収まるように補正した．このように補正することによってオーバーシュート，アンダーシュートの成長を抑えることが可能となり安定的に計算を行うことができる．
次に式（4）の非移流項については，以下のように行う．格子点定義である油層厚は，式（4）に中央差分を適用することで以下のように計算することができる．ただし，表記を簡単にするために，拡散係数Dを定数と置いた表現である．

― 106 ―

図－7

∗
���
ℎ�,�
� ℎ�,�
� Δ𝑡𝑡� �

計算フローチャート

∗
∗
∗
∗
∗
∗
� 2ℎ�,�� ℎ���,�� 2ℎ�,�
� ℎ�,���ℎ���,�ℎ�,����� � �� Δ𝑡𝑡
Δ𝑥𝑥 �
Δ𝑦𝑦 �

（30）

一方，油層厚の変化に伴い，線，面積分値も変化する．これらを上式で更新された格子点定義の油層厚から単純に計算してしまうと，CIP-CSL2で計算する利点がなくなるため，ここでは通常のCIP法と同様な考え方で以下のように線，面積分値を計算する．
まず，非移流項をGと置いた，
���
∗
ℎ�,�
� ℎ�,�� � ∗ Δ𝑡𝑡を考える．ここで，𝐺𝐺�,� � � �ℎ���,� � 2ℎ�,� � ℎ���,� ℎ�,��� � 2ℎ�,� � ℎ�,����� � ��
Δ𝑥𝑥 �
Δ𝑦𝑦 �

（31）

（32）

である．上式を例えばx方向にセルi～i+1まで積分すれば，下記の線積分値σxに対する時間更新式が得ら
れる．
����

𝜎𝜎����
� 𝜎𝜎�∗�,� � ��,���� ∗ Δ𝑡𝑡 𝑑𝑑𝑑𝑑

（33）

右辺第二項のGについて台形公式でその積分値を近似し，Gの定義を利用する．すなわち，
����

�

��

∗
∗
� 𝐺𝐺���,�
𝐺𝐺�,�
𝐺𝐺 �𝑥𝑥 �
Δ𝑥𝑥2∗― 107 ―

（34）

図－8

ୢ

௛

ୢ

௛

a)ୢ୶ ቀ ௃ ቁ , ୢ୷ ቀ ௃ ቁ の差分の定義点，b)差分の手順

∗
𝐺𝐺�,�
�

∗
���
ℎ�,�
� ℎ�,�
,Δ𝑡𝑡∗
𝐺𝐺���,�
�∗���ℎ���,�� ℎ���,�
Δ𝑡𝑡（35）を用いて式（33）を書き換えれば，σxの更新式として次式が得られる．

𝜎𝜎����
� 𝜎𝜎�∗�,� �
�,�

Δ𝑥𝑥 ���∗���∗��� ����,�� ��,�� ��,�
�
2 ���,�

𝜎𝜎����
� 𝜎𝜎�∗�,� ��,�

Δ𝑦𝑦 ���∗∗������ ��,���� ��,�� ��,��
2 �,���

（36）

同様の考え方をσy，ρにも適用することで，次式を得ることができる．

∗
���
𝜌𝜌�,�
� 𝜌𝜌�,�
�

Δ𝑥𝑥𝑥𝑥𝑥 ���
∗
∗
∗
∗
���
���
���
�����,��� � ����,���� ��,���� ��,���� ��,���� ��,���� ��,�� ��,�
�
4

（37）
（38）

以上が，CIP-CSL2法を適用した移流拡散方程式の解法である．これらの手順をまとめると図－7のようになる．ただし，ここで用いるNays2D Floodは一般座標変換された基礎式を差分化するため，これらの差分式を一般座標変換した以下の式に適用する．

𝜕𝜕 ℎ𝜕𝜕 𝑢𝑢� ℎ𝜕𝜕 𝑢𝑢� ℎ𝑄𝑄�� �� ��� ���𝐽𝐽𝜕𝜕𝜕𝜕 𝐽𝐽
𝜕𝜕𝜕𝜕 𝐽𝐽
𝜕𝜕𝜕𝜕 𝐽𝐽

ここに，

�

𝜕𝜕
𝜕𝜕 ℎ
𝜕𝜕
𝜕𝜕 ℎ
𝜕𝜕
𝜕𝜕 ℎ
𝜕𝜕
𝜕𝜕 ℎ
�𝜕𝜕� 𝐷𝐷 � �� � �𝜕𝜕� 𝐷𝐷 � �� � �𝜕𝜕� 𝐷𝐷 � �� � �𝜕𝜕� 𝐷𝐷 � ��𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝐽𝐽𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝐽𝐽𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝐽𝐽𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝐽𝐽

𝜕𝜕 ℎ
𝜕𝜕 ℎ
𝜕𝜕 ℎ
� � � 𝜕𝜕� � � � 𝜕𝜕�� �,𝜕𝜕𝜕𝜕 𝐽𝐽
𝜕𝜕𝜕𝜕 𝐽𝐽
𝜕𝜕𝜕𝜕 𝐽𝐽
𝜉𝜉� �

（39）

𝜕𝜕 ℎ
𝜕𝜕 ℎ
𝜕𝜕 ℎ
� � � 𝜕𝜕� � � � 𝜕𝜕�� �𝜕𝜕𝜕𝜕 𝐽𝐽
𝜕𝜕𝜕𝜕 𝐽𝐽
𝜕𝜕𝜕𝜕 𝐽𝐽

𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
, 𝜉𝜉� � , 𝜂𝜂� � , 𝜂𝜂� �
,
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

� � 𝜉𝜉� 𝜂𝜂� � 𝜉𝜉� 𝜂𝜂�

（40）
（41）

である．式（39）を先述した移流項・非移流項の解法を用いて計算を行うが，式（39）の計算の際，
図－8aのような位置に ∂�h/J�/ ∂x ， ∂�h/J�/ ∂y が定義される．これらの計算には，図－8bの様な差分を取る必要があり，格子点定義のhだけでなく，セル中央のhも必要になる．セル中央のhを求めるには，
通常周囲の4点平均値を用いる方法が使われるが，CIP-CSL2では，セル中央にhのセル内体積である面

積分値ρが定義されている．したがって，セル中央の平均的なhは，h = ρ/（∆x∆y）で計算することが可能
であり，補間を行う必要がない．この点は，CIP-CSL2を用いる利点であると考えられる．

― 108 ―

a)100時間後の拡散状況
図－9b)拡散半径の時間変化
点源に瞬間放出された油の拡散計算結果．

（c）オイラー型油移動拡散モデルの基本的特性
ここでは，前節で示した数値計算モデル，計算スキームが油輸送拡散現象を表現する上で，どの程度の性能を有するものか，単純な場において検討するものである．
（ア）静水場における瞬間流出油拡散
まず，最も単純な場として，流れがない場に，点源に瞬間的に放出された油の拡散について検討を行う．松崎・藤田3）を参考に，単純場において総量10000 m3の油を面積662400 m2から瞬間流出させ，100時
間後までの油拡散について計算した．また，使用した格子は151×151（格子サイズ:40×40m）であり，
油拡散のパラメータとしては，水の密度は 𝜌𝜌�＝1000 kg/m3，油の密度は 𝜌𝜌� ＝861 kg/m3，油－水間の抵抗に関するキャリブレーションパラメータ f は0.01 kg/m2sとして解析した．

図－9aに流出後100時間の結果を示す．また，本結果の拡散時間と拡散半径の関係を図－9bに示す．この際，拡散半径の算出は，中心を起点とし，油の厚さが 10ି଻m以下となる範囲までとした．100時間を300秒のステップで値をプロットした．一方，油の拡散半径についてはFay1）により理論的な検討が実施されている．Tkalich et al.2）のモデルで用いられている油拡散を示す拡散係数は，Fayの検討における重力-粘性力つり合い領域のものである．この点について，再度拡散係数の導出を踏まえて，拡散半径が時間とともにどのように変化するか，理論的に検討する．
図－10のように水面に浮かぶ油を三角形要素として考える．この時，頂点に対して右側の要素に作用する力を考える．水と油の密度差により，1 𝜌𝜌� � 𝜌𝜌�
𝑔𝑔𝑔𝑔𝑔
2 𝜌𝜌�

（42）

のような浮力が生じている．この時，油と水境界に沿った力は，l >> hとすると，

1 𝜌𝜌� � 𝜌𝜌�ℎ 1 𝜌𝜌� � 𝜌𝜌� �
𝑔𝑔𝑔𝑔𝑔 �𝑔𝑔𝑔2 𝜌𝜌�𝑙𝑙 2 𝜌𝜌�

（43）となり，これが油を側方に移動させようとする代表的な力となる．これに抵抗する力として，粘性力を考える．先ほどと同様にl >> hとすると，粘性力は，�� � ��

𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

（44）

と近似される．今，ある時間Tで油が半径lになったとすると，粘性力を代表する速度はU=l/Tとなる．
また水－油間の薄い層厚をδとすれば，結局― 109 ―

図－10

油に作用する力：重力－粘性力卓越領域

図－11

油の拡散半径の時間変化とFay式の比較

���

𝑙𝑙 �
𝑇𝑇𝑇𝑇

（45）

である．これらの力がおよそ同じオーダーであると考えると次式が得られる．

𝑙𝑙 � ∝

1 𝜌𝜌� � 𝜌𝜌� 𝑔𝑔𝑔�𝑇𝑇
2 𝜌𝜌� 𝜇𝜇/𝛿𝛿（46）

ランダムウォークにおける平均移動距離∆xと拡散係数Dの関係は，∆x2=2DTであり，この∆xが拡散半径lに相当すると考えると結局，

𝐷𝐷 𝐷𝜌𝜌� � 𝜌𝜌� 𝑔𝑔𝑔�𝜌𝜌� 𝜇𝜇/𝛿𝛿（47）なる式が得られる．油－水間の薄い層厚δは不確実パラメータであると考え，µ/δをキャリブレーション
パラメータと置いたものが，本研究で用いている油拡散モデルにおける拡散係数である．
さて，この関係の時，油拡散半径がどのようになるかについては，式（47）において，δが一定であるとし，また油体積Vは保存されるため，h∝V/l2の関係にあるとすると，
�

�
�

�
𝜌𝜌� � 𝜌𝜌� 𝑔𝑔𝑔𝑔���
� 𝑇𝑇 �
𝜇𝜇
𝜌𝜌�
𝛿𝛿

（48）

が得られる．この関係が成立しているかについて示したものが図－11である．計算における油半径の時間変化は上記考察と同様の傾向を示しており，構築した数値計算モデルはこの傾向を矛盾なく評価できていることがわかる．

― 110 ―

図－12

図－13

構造物を含む流れ場と油の条件設定

単純な流れ場における油移動計算．オリジナルのCIP-CSL2による結果．

（イ）構造物を含む流れ場における油移動拡散
流れがある場におけるモデルの特性を評価するために，図－12に示すような長さ30m，幅3mの直線
水路に15cm角の構造物を設置し，背後にカルマン渦列が形成されるような条件を設定した．そのほかの条件は，勾配が1/1000，流量0.1m3/s，マニング粗度係数0.03であり，格子サイズは両方向ともに5cmである．計算では，上流側に油を瞬間流出させ，その移動過程について計算するとともに，油の体積保存性について確認した．なお，移流スキームの評価を行うために，物理的な拡散については考慮していない．この条件において，前章で示したオリジナルのCIP-CSL2，有理関数を導入したCIP-CSL2（R-CIP-CSL2とする），有理関数とともに，内田・河原の手法でオーバーシュート等を防いだ方法（mR-CIP-CSL2とする）について比較を行った．
図13～15にそれぞれの移流スキームで得られた油の移動過程について示す．また，これらの計算における油の総量に対する質量損失誤差を示した結果が図－16である．移流計算の結果を見ると，どの計算結果も同様な傾向を示しているが，体積保存に着目するとオリジナルのCIP-CSL2法は，体積保存が大きく崩れている．これは，オーバーシュートとアンダーシュートが発生することにより，負の油厚が― 111 ―

図－14

単純な流れ場における油移動計算．R-CIP-CSL2による結果．

図－15

単純な流れ場における油移動計算．mR-CIP-CSL2による結果．

― 112 ―

図－16

体積保存に関する検討結果．

生じ，計算上これをゼロと置きなおしているためである．この点については，R-CIP-CSL2，mR-CIPCSL2では大幅に改善されている．角柱に衝突した時点で，体積保存が崩れているため，セル境界における境界条件等について改善の余地があると考えられるが，この点を除けば，mR-CIP-CSL2が最も誤差の増加率が小さく，質量保存の観点より安定した計算スキームであることがわかる．すなわち，CIPCSL2においては，オーバーシュート等の誤差を防ぐために，有理関数の導入が重要であるだけでなく，保存量を計算する際に格子点，線，面積分値がずれて配置されるために，Staggered格子において保存量を決定する流速定義点ともずれが生じる場合があり，これらのずれに起因する質量保存誤差までをも補正する必要があることに注意が必要である．
（d）六角川油流出事故に対する適用前章で検討した油移動拡散モデルについて，令和元年8月佐賀豪雨災害において生じた六角川油流出
イベントに適用し，モデルの適用可能性について検討を行った．
（ア）六角川油流出事故概要
令和元年8月27日から8月28日にかけて秋雨前線の活発な活動により九州の広範囲に強い雨雲がかかり
佐賀県全域にかけて大雨が発生した．特に六角川水系六角川を有する佐賀県大町町においては，近傍の
南渓雨量観測所で，およそ毎時100mmの雨が二時間にわたり観測され（図－17），六角川左岸側で内水氾濫が発生し，最大3mを超える浸水を記録した．この内水氾濫により，大町町に位置する佐賀鉄工所
からクエンチオイルが流出する事故が発生した．
（図－1）外部への油の流出量は5.4万リットルにも達し，
この流出した油は鉄工所の東と南東方向に少なくとも1kmの広範囲にわたって流出，下流に位置する住宅地や病院にまで流れ込み，田畑の汚染や居住困難など甚大な被害を発生させる国内では例を見ない最大規模の油流出事故となった．
（イ）計算条件
計算は図－18aに示す範囲とし，601×501（6m×6m）の格子数で計算格子を作成した．この領域に与えた降雨は図－17に示す南渓観測所で観測された降雨のうち8月28日3時からとした．なお，地形デー
タは国土地理院基盤地図情報数値標高を使用して，計算格子の初期地形としたが，予備計算の結果，工― 113 ―

時間雨量(mm/h)
8/27 0:00

8/27 0:00
図－17

図－18

8/27 0:00

南渓地点における観測雨量a)計算領域，b) 微地形と家屋の考慮場脇に存在する側溝といった微地形を考慮する必要があったため，図－18bのように地形を編集することでこれを考慮した．さらに，建物といった局所的な構造物もまた，油の輸送に強く影響するため，航空写真を基に，図－18bのように建物群については不透過セルとして考慮している．また，実際の事故報告書と同様，計算開始から2時間後～10時間後まで油が流出したものとし，油の流出量の合計が実際の事故と同様の54,000m3となるように油の流出を与えた（Qo=18.7ℓ/s）．本事故で流出した油はクエンチ
オイルと呼ばれる冷却用の油であったため，クエンチオイルの油の密度であるρo＝861 kg/m3を油の密度として与え，水の密度はρw＝1000 kg/m3，キャリブレーションパラメータであるfは0.01kg/m2sとして前章で行った計算と同様のパラメータを用いた．
（ウ）計算結果
図－19aに流出開始2時間後の結果を示す．結果を見ると油は流出後東に輸送されている．つまり油は初期段階においては東方向に流れ，その後解析結果のように時系列的に流向が変化し図－19bの様に北東方向に流下したものと考えられる．これは実際の油の時系列挙動データは存在しないため確認することはできないが，図－1からも分かるように，東側家屋より油が輸送されていることが確認でき，時系列的に流向が変化したことが考察できる．この結果を見ると，構造物として指定した家屋を避けるよ
うに流下しており，家屋の有無による影響を考慮できていると考えられる．また，図－20におよそ同時刻における計算結果と現地の撮影写真との比較を示す．この図よりわかるように，油の流出の方向に― 114 ―

図－19

六角川における油流出に関する計算結果．

図－20

計算結果と現地の油輸送拡散状態の比較．

関して北東方向で一致している．また，その後の油の流向に関しても実際の事故と同様の結果を示しており，実際の油の挙動と類似した油の輸送を計算にて再現できたものと考えられる．加えて，今回の物
質の輸送を解析にて行う上では，通常LPデータなどのデータに反映されていない微地形などの影響が，油などの物質輸送解散現象に大きく影響することが分かった．一方，計算結果を見ると，計算モデルは油の移動距離について過小評価している傾向にある．これは，ここでは考慮していない風や水面勾配の影響の可能性があり，モデルを改良する必要性を示唆している．
（3）まとめ本研究では，今後多発すると考えられる非海上における油流出によるリスクを事前に予測・評価することを目的として，非定常二次元浅水流モデルであるiRIC-Nays2D FloodにEuler型の油輸送方程式を組み込むことにより，非海上における氾濫流に伴う油輸送解析モデルを作成した．構築したモデルを用いて簡易な場における油の拡散現象を行い，油の拡散の理論式であるFayの式と比較することで本モデルの妥当性を検討し，六角川で発生した内水氾濫による油流出過程を解析し，モデルの性能と再現性，有用性の検討を行った．以下，本研究で得られた主な結論を示す．
・保存型移流スキームであるCIP-CSL2法を用いた油輸送モデルを構築し，基本的な場において体積保存性などを確認し，モデルの基本的性能について整理した．
・油拡散解析結果の範囲，時間変化はFayの第二段階と一致しており，本モデルの妥当性が確認された．

― 115 ―

・本研究で構築したモデルを六角川油流出イベントに適用し，地形の起伏により時空間的に変化する氾濫流に伴う油の輸送拡散を概ね再現できた．
・適切な地形データを取り扱うことで，再現性を十分に発揮することができると考察され，モデルの有用性，再現性が確認された．
・氾濫流による油など物質の輸送拡散は局所的な流れ場に敏感に反応し，通常LPデータなどデータに反映されていない微地形の影響が，局所流となって油などの物質輸送拡散現象に大きく影響することが分かった．そのため，本モデルを事前計算に用いる際は，対象地域，特に流出点付近の微地形を十分に加味した地形データを用いる必要があるといえる．
一方，ここで示したように構造物が存在する場における体積保存や，流下距離の再現性について課題を残しており，このような点を改善することで，モデルについて精緻化したいと考えている．
参考文献
1）Fay, J.A.: The spread of oil on a Calm Sea, Plenum Press, pp.53-64, 1969.
2）Tkalich, P., Huda, M.D.K. and Gin, K.Y.H.: A multiphase oil spill model, Journal of Hydraulic Research,
Vol.41, pp.115-125, 2002.
3）松崎義孝，藤田勇：海水面における流出油の拡散・移流に関する数値計算法の開発と油流出事故の再現計算，土木学会論文集 B2（海岸工学），Vol.70，No.1，pp.15-30，2014．
4）Rajendran, S., Sadooni, F.N., Al-Kuwari, H.AS. et al.: Monitoring oil spill in Norilsk, Russia using satellitedata, Scientific Reports, 11, 3817, 2021.
5）Wang, S., Yang, Y., Zhu, Z., Jin, L. and Ou, S.: Riverine deposition pattern of oil-particle aggregatesconsidering the coagulation effect, Science of the Total Environment, 739, 2020.
6）Nakamura, T., Tanaka, R., Yabe, T. and Takizawa, K.: Exactly Conservative Semi - Lagrangian Scheme forMulti-dimensional Hyperbolic Equations with Directional Splitting Technique, Journal of Computational
Physics 174, pp.171–207, 2001.
7）内田龍彦，河原能久：二次元浅水流の保存型CIP陽解法の開発とその検証，応用力学論文集，
Vol.9，pp.917-924，2006．

― 116 ―

3． 比重の大きい流木の河川水中の挙動のモデル化とiRICへの実装（担当：木村一郎）（1）本章の目的と概要
一般に，流木は水に浮いて移動するものと考えがちであるが，必ずしもそうではない．流木となった
初期の状態が浮遊状態であったとしても，時間の経過とともに木の組織が吸水し，比重が1を超え，沈降する場合も多々見受けられる．このような場合，流木の移動は流水の影響以外に，河床との抵抗の影
響を受けることになる．したがって，比重が1を超える場合の流木挙動のモデリングについては，河床との摩擦力を適切に考慮することが重要となる．
最近，多く見られる豪雨災害の特徴の一つとして，大量の流木の流出と，それらの橋梁等の構造物や家屋等の建築物への衝突による破壊，あるいは流木堆積による河道の閉塞などの問題が挙げられる．これらの災害を軽減するには，流木発生そのものを抑制する対策，河道に流出した流木を捕捉して流れから除去する方策などが考えられる．しかし，前者は流木発生可能性のある個所が広範に及ぶため，現実的とは言えない．後者は適切な設計が為されればピンポイント的な施設でも効果が期待でき，経済性や即効性の点でも優れる手法といえる．また，捕捉した流木を資源として利用できるのも長所の一つといえる．
一方，流木災害軽減を目的として，河道中の流木を捕捉する施設にはいくつかの種類が考えられる．
その中でも代表的なものは，河川側岸にバイパス水路あるいはワンド状の入り江を設け，流木をそこに誘導し，堆積捕捉させるというものである．この種の流木捕捉工は，河道の本流の流れを阻害する可能性が低く，治水上でも有利という点が長所の一つといえよう1）．しかし，必然的に本流から逸れた箇所に流木貯留地が設けられるため，いかに，この箇所に流木を誘導するかという点が捕捉工の成否を握っているといえる．岡本ら2）は，直線水路の側岸にワンド状の矩形死水域を設け，流木をそこに効率的に誘導する方策として，主流のワンドと反対側の岸に水制を設置する方法が有効であるとし，このことを
水理模型実験により示した．岡本らによると，水制の水没部分を透過型（流水が通過できる形）にした
場合でも，ワンド側への流木誘導の効果が十分に発揮されることを示した．一方で，加藤ら3）は，岩手県の小本川で計画されているワンド型流木捕捉工について，精緻な水理模型実験により検討を行った．
この捕捉工では，山間部の河川の湾曲部の外岸側に捕捉地を設置することにより，慣性力を利用して流木を貯水池に誘導するというものである．加藤らは比較的大規模な水理模型実験により，流入口の形状の影響などを詳細に検討し，投入流木量の65%以上を捕捉できる可能性があることを示している．
流木捕捉工周辺等の流木の挙動を検討を行う方法として，室内水理模型実験とともに，数値シミュレ
ーションモデルによる方法が有効と考えられる．河道内の流木挙動や，捕捉工への流木堆積などの数値解析モデリングについては，これまで種々の検討が行われてきている．清水ら4）は平面二次元水深積分モデルによる河川流モデルと流木を球体の結合で表現し，流木が水面付近を二次元的に移動すると仮定
するとともに，流木同士の衝突に個別要素法（DEM: Discrete Element Method）を用いる二次元型流木追跡モデルを考案し，橋脚による流木捕捉現象を適切に再現できることを示している．一方，初田ら5）は，同様なモデルを用いて河川蛇行部における流木挙動を適切に再現できることを述べている．この種のモ
デルは，流れと流木の双方を二次元的にモデル化しているため，2D-2Dモデル（2D河川流-2D流木モデ
ル, 2D: Two-Dimensional）と呼ぶことにする．一方で，長田ら6）は，流れ場は二次元とし，流木挙動を三次元に拡張するとともに，流木の弾性変形を考慮したモデルを構築した（2D-3Dモデル：2D河川流-3D流木モデル）．しかし，河川湾曲部などでは流れの三次元性が卓越することから，これを考慮するモデルとして，Kimura and Kitazono7）は，3D-2D型モデルを構築し，湾曲部の流木挙動や橋脚による流8）
木捕捉の検討に用いている．これに対し，加藤ら（2019）
は，流れの計算にEngelund9）による第一種二次流のモデルを組み込んだ，準三次元モデルを採用したq3D-2Dモデル（qusi2D-3Dモデル）というべ― 117 ―

き興味深いモデルを提案している．このモデルは，流木のモデル化にも球体の結合ではなく，円柱状の形状を直接採用しているため，計算機負荷が小さく，大規模な流域の計算にも適用が容易と考えられ，注目に値する．このモデルは，岩手県小本川における捕捉工の水理実験の再現が試みられ，概ね妥当な結果が得られている．しかし，このモデルでは，無限長の一様湾曲部に生じる十分発達した二次流場を想定している．このため，実河川でみられる有限長の湾曲部では二次流を過大評価する可能性がある．
また連続蛇行部では二次流発達の遅れによりEngelundモデルとは逆回転の二次流が発生する場合があり，このような場合には流れ場により二次流の遅れを考慮した高度な準三次元モデル（例えばKimura etal.10）やUchida and Fukuoka11））を用いるか，三次元的なモデルを用いることが必要となるだろう．
本研究は，河川流と流木の双方を三次元的に取り扱う3D-3D型モデルにより，加藤ら3）による小本川の流木捕捉工実験の再現を試みるものである．計算モデルは木村12）が橋脚周辺の流木捕捉の計算に用いた3D-3D型の流木モデルを基本とする．しかし，加藤らの実験では，流木化した広葉樹を対象に，比重が1.1の流木模型が用いられている．比重が大きい流木では，橋桁等への衝突堆積は生じにくいと考えられる一方，河床への沈下による堆積や，それに伴う水位堰き上げ等の問題が容易に生じることが予想される．比重が1より大きい流木は，基本的に河床に沿って，河床と頻繁に接触しながら移動するため，河床と流木の摩擦を適切にモデル化することが必要である．これについては，Kang and Kimura13），およびKang et al.14）が，浅水流を対象とした流木の河床への接触，停止，再移動などを再現した2D-2D型モデルを提案し，実験結果との比較を通じておおむね妥当な結果が得られることを示している．本研究では，このモデルを参考に，底面摩擦項を三次元モデルに組み入れる．なお，流木の形状が円柱状の場合，底面付近の挙動は流木軸方向のスライド型移動と横断方向のローリング型移動では摩擦力が大きく異なるため，摩擦力の平面的異方性を考慮したモデルを採用する．これに加え，Kang et al.14）は，既往モデ
ル13）では無視されてきた流木と流れとの衝突角度による抗力の変化（投影面積の変化）についても考慮したモデル化を行っている．これについても，今回の研究で三次元モデルに組みいれる．
以上のような指針のもとで構築した三次元モデルを，iRICソルバーNaysCUBEに組み入れるとともに，加藤らの実験と同条件で適用し，その再現性を検証する．
（2）数値解析モデル（a）河川流（液相）のモデル
本研究で構築する3D-3D型流木数値解析モデルでは，河川流を三次元モデルで計算する．用いる三次元モデルの基礎となるベースモデルは，既往の研究11）と同様に，河川解析共通プラットフォームiRIC上で稼働する三次元ソルバー，NaysCUBE14）である．このソルバーでは，移動境界適合格子系上でレイノ
ルズ平均されたナビエ・ストークス方程式と連続式を三次元的に解くもので，ベクトルおよびテンソルについては反変成分を基本変数として用い，スタガード格子上で変数配置を行っている．移流項の離散化には空間三次精度のTVD MUSCL法が，時間積分には二次アダムス・バシュフォース法が，乱流モデ
ルには非線形k-εモデルが採用された．水面移動については，境界適合法が採用されており，計算格子が上下方向に移動し，最上面の格子面が常に水面に位置するように移動することで水面移動を捕捉計算する．
基礎式や離散化方法，境界条件等についてはweb上の文献に詳述されているのでそちらを参照された
い15）．
流木と河川流の相互作用について，既往の研究では簡単のため流木が流れに及ぼす影響を無視した
One-way型モデルを採用されたものも見受けられる．しかし，本研究では流れへの影響を抗力として考慮したTwo-way型のモデルを行っている点が異っている．計算セルにおける流木の流れに対する抗力は，次の式で評価する．

― 118 ―

（1）
ここに，Fdr: ある計算セルにおける流木が流れに及ぼす抗力のi方向成分，ρ: 水の密度，CD: 抗力係数，

ඥg : ヤコビアン，A: セルに存在するk番目の構成要素の流れ方向投影面積，u: 流速のベクトル，up: 流

木構成要素移動ベクトル，Ui: 流速ベクトルの反変成分，Upi: 流木構成要素移動ベクトルの反変成分，i: 一般座標系（ 𝜉𝜉, 𝜂𝜂, 𝜁𝜁 ）のいずれか，Ncell: 当該計算セル内に存在する流木の構成球体要素数をそれぞれ表す．

抗力係数CDについては，次のように粒子レイノルズ数Redの関数とした．

24
� 0.4,
𝑅𝑅𝑅𝑅� � 1000
𝑅𝑅𝑅𝑅�
𝐶𝐶� �
⎨ 24 �1 � 0.15𝑅𝑅𝑅𝑅 �.��� �, 𝑅𝑅𝑅𝑅 � 1000
��
⎩𝑅𝑅𝑅𝑅�
⎧∵ 𝑅𝑅𝑅𝑅� ��� � �𝐩𝐩 ��
𝜈𝜈

（2）

（3）

ここに， 𝑑𝑑 : 流木構成要素球体直径，および， 𝜈𝜈: 動粘性係数である16）．この定式化では，抗力係数をレ

イノルズ数（球体直径を長さスケール，流れと球体速度の相対速度の大きさを流速スケールとする）の関数としている点に特徴がある．
（b）流木（固相）のモデル（i）流木モデルの概要
本研究で用いる流木モデルは，清水ら4）にならい一本の流木を球体の接続で表現する（図1参照）．
この手法は，越塚18）がMPS型粒子法において示した計算方法を二次元流木計算に応用した初田ら5）やKimura and Kitazono7）の方法と類似であり，これを三次元空間に拡張したものといえる．計算方法の詳細については既報11）を参照されたい．

図1

流木を球体要素の列で表現するモデル（ii）底面相互作用項の導入
前述のように，本研究では比重が1より大きい流木を対象としたモデリングを行う．このため，流木移動は基本的に底面付近の移動となり，底面との摩擦を伴うこととなる．既往のモデル11）では比重が1以下の流木を想定しいたため，この底面摩擦は考慮されていなかったので，新たにこの効果の導入を行った．
まず，流木構成球体を拘束を無視し，ラグランジュ方程式により計算ステップΔt間移動させる．
底面摩擦項を考慮したラグランジュ方程式は次のようになる．

― 119 ―

𝑑𝑑𝑑𝑑𝐩𝐩 1𝐴𝐴� 11𝜆𝜆����� �1 � 𝐶𝐶� � 𝑑𝑑𝑑𝑑 𝜎𝜎𝜎𝜎 � 𝜆𝜆�����
� 𝐶𝐶� 𝜆𝜆�����
�𝑑𝑑 � 𝑑𝑑𝐩𝐩 ��𝑑𝑑 � �𝒑𝒑 � � 𝜎𝜎
� 𝜎𝜎
𝐠𝐠
𝑑𝑑𝑑𝑑
𝐴𝐴� 𝑑𝑑 𝜎𝜎𝜎𝜎 � 𝜆𝜆����� 𝐶𝐶�2� 𝜆𝜆����� 𝐶𝐶� 𝑑𝑑𝑑𝑑
� 𝜆𝜆����� 𝐶𝐶�
𝜌𝜌
𝜌𝜌

（4）

1
�
�� � �𝐛𝐛𝐛𝐛𝐛𝐛 ��𝐴𝐴� 𝑑𝑑 𝜌𝜌�𝜎𝜎𝜎𝜎 � 𝜆𝜆����� 𝐶𝐶� � 𝐩𝐩

ここに， 𝜎𝜎：球体の密度， 𝜌𝜌：水の密度，CD : 抗力係数，u：水流の流速ベクトル，up：粒子移動速度ベ
クトル，g：重力加速度ベクトル，F p：粒子間衝突力ベクトル，t：時間，C M：付加質量係数（C M =0.5），A2, A3：二次元および三次元形状係数（A2 = π/4, A3 = π/6）
，λA_sub：球体の流水方向投影面積の水面以下の部分の割合，λV_sub：球体体積の水面以下の部分の割合 ，Fbed: 底面摩擦力ベクトルをそれぞれ表す．
Fbedは底面摩擦力ベクトルであり，次のように表す．

𝐅𝐅��� � �

𝟎𝟎 ,𝐅𝐅𝐛𝐛 � 𝐅𝐅𝐬𝐬 ,if 𝑧𝑧� � �/2

if 𝑧𝑧� � �/2（5）

ここに， z୮ : 球体中心座標の河床からの高さであり，Fbは底面摩擦ベクトル，Fsは河床勾配に伴う重力の影響を表す．まず，Fbについて，次式で表す．

𝑢𝑢��
𝑢𝑢
𝐅𝐅Fb � ��� 𝑁𝑁𝐩𝐩 � ���� �� �
𝐛𝐛
��𝐩𝐩 � 0𝑁𝑁（6）

ここに， 𝜇𝜇�: 球体と底面の摩擦係数であり，詳細は後述する．また，Nは球体の底面に対する垂直力であり，底面勾配は小さいと仮定して，近似的に次式で表す．

ܰ ൌܹെ‫ܤ‬（7）

ここに，WおよびBは球体に作用する重力及び浮力であり，次のように表される．

��
��

𝜋𝜋 �
𝜎𝜎𝜎𝜎 𝑔𝑔6𝜋𝜋 �𝜌𝜌𝜌𝜌 𝜆𝜆����� 𝑔𝑔
6

（8）

（9）

式（9）の λ�����は，前述のように流木構成球体の水没部の体積割合を表す係数であり，河床底面付近の流木では通常1（完全水没状態）であるが，水深が小さい場合などには，この値が1以下となる場合もありうる．

摩擦係数 𝜇𝜇�について，流木が円柱形状の場合，流木軸方向の滑動と横断方向の転動の摩擦係数が大きく異なる．これについて，Kang and Kimura（2018）12）およびKang et al（2020）13）にならい，摩擦力異方性を考慮したモデル化を行う．これについて，図2は流木が底面付近に存在する様子を平面的に見たものである．流木軸方向（t方向）と横断方向（n方向）の摩擦係数の分布を図のように楕円型で仮定する．流木軸方向の滑動の摩擦係数をμ� ，横方向の転動の摩擦係数を μ�とすると，流木の進行方向の摩擦係数μ�は，次式で表される．

― 120 ―

図2

𝜇𝜇� �

摩擦力異方性のモデル化の説明

𝜇𝜇� 𝜇𝜇��𝜇𝜇�� �1 � cos � 𝜓𝜓� � � 𝜇𝜇� cos � 𝜓𝜓�（10）

ここに， 𝜓𝜓� は，流木の進行方向と流木の長軸との成す角度であり，この関係により， cos 𝜓𝜓�は次式で求めらえる．

cos 𝜓𝜓� �𝑢𝑢�� cos 𝜃𝜃 � 𝑢𝑢�� sin 𝜃𝜃

（11）��𝐩𝐩 �

ここに， 𝜃𝜃は，x軸と流木軸の成す角度である．なお， 𝜇𝜇� については，静止摩擦係数 𝜇𝜇��と動摩擦係数 𝜇𝜇��の相違（ 𝜇𝜇�� � 𝜇𝜇�� ）を考慮する．摩擦係数の値を理論的に決めるとは不可能であると考えられる．そこで本研究では，実験における流木挙動を，より忠実に再現するよう試行錯誤を行い，表1の値を用いた．
表1

用いた摩擦係数の値

名称

記号
μ��

静止摩擦係数滑動摩擦係数
転動摩擦係数

μ��
μ�

値
2.0
1.0
0.5

底面に接触して滑動，転動する比重の比較的大きい流木の移動形態は，底面の勾配の影響を受ける．
これを考慮するため，底面付近に局所的に底面に沿ったx’-y’平面を考える．x’, y’方向の底面勾配は，底面勾配が比較的小さい場合には，一般座標におけるξ,η方向の勾配を用いて，次のように近似的に表される．

― 121 ―

∂𝑧𝑧� 𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕���,𝜕𝜕𝜕𝜕 �
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕∂𝑧𝑧� 𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕���𝜕𝜕𝜕𝜕 �
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕（12）

ここに，zbは底面の標高である．これより，流木に働く重力の斜面方向のベクトルFsのx, y方向の成分を
次のようにモデル化する．

𝐹𝐹�� � ��

図3

図4

𝜕𝜕𝜕𝜕�
𝜕𝜕𝜕𝜕�
, 𝐹𝐹�� � ��
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

（13）

流向と流木軸の角度の関係

中間部の球体（左）と終端部の球体（右）の投影面積の考え方（iii）流木の方向による抗力変化のモデル化
本研究では，一本の流木を球体要素の接続でモデル化している．このため，抗力係数は球体における値17）を用いている．しかしながら，流木の場合，河川の流れの方向と流木軸の方向によって抗力係数が変化する．これについて，円柱状の物体の場合の迎え角と抗力の関係を考察した最近の研― 122 ―

14）
究19）を参考に，抗力係数を変化させる方法が考えられる．一方で，Kang et al.（2020）
は，二次元モデルにおいて，抗力係数を変化させることなく，流木の流れに対する角度によって，各球体要素
に関する投影面積が変化すると考えた．これら2つの研究は，角度による抗力の変化を抗力係数に含めるか，投影面積に含めるかの相違であり，本研究では後者を採用する．Kangらにならい，投影面積に変化は．図3，図4を参考に，次のように記述される．

𝐴𝐴̅ � 𝐴𝐴|sin𝜙𝜙| �� � 2～M � 1�𝐴𝐴̅ � 𝐴𝐴�|sin𝜙𝜙| � max�cos 𝜙𝜙, 0�� �� � 1�𝐴𝐴̅ � 𝐴𝐴�|sin𝜙𝜙| � max�� cos 𝜙𝜙, 0�� �� � M�

（14a）
（14b）
（14c）

ここに，A: 球の射影面積（=π/4 d 2）, 𝐴𝐴: 補正後の投影面積，φ: 流木軸と河川の流れの相対速度の成す
角度，M: 一本の流木を構成する球体の数である．なお，流木軸の方向は球体番号m = 1の球体から，もう一方の端のm = Mの球体に向かう方法を正とする．また， cos 𝜙𝜙 および | sin 𝜙𝜙|については，次のように求められる．

��� � �𝑢𝑢� ′𝑡𝑡� � 𝑢𝑢� ′𝑡𝑡� � 𝑢𝑢� ′𝑡𝑡��𝑢𝑢𝑢�� � 𝑢𝑢′�� � 𝑢𝑢′��|sin 𝜙𝜙| � �1 � cos � 𝜙𝜙

（15）

（16）

ここに， �u�� , 𝑢𝑢�� , 𝑢𝑢�� � : 河川の流れの球体に対する相対速度， �t � , t � , 𝑡𝑡� �: 流木軸方向の単位ベクトルをそれぞれ表す．

図5

流木捕捉地予定地の状況（岩手県小本川）― 123 ―

（3）流木捕捉工周辺の流木動態への適用（a）加藤らの実験の概要
加藤ら3）は岩手県小本川の山間地湾曲部に計画されている流木捕捉工を対象に模型実験を行った．図5の写真は現地の様子の航空写真である．一方，図6は，計画されている捕捉工の平面図を示している．
もともとの河道の内岸側に新河道を掘り，外岸側のもとの河道を流木捕捉工として利用する構造となっている．したがって，流木を捕捉工の方向に，自然な流れを用いて誘導しやすい形状となっていると考えられる．
加藤ら，は流入口の位置や大きさを数通りに変化させて水理模型実験を行い，その中で最も高い捕捉性能を有する形状を実験により求めた．図7は，そのようにして最適化された形状を表している．実験
装置のスケールは1/50である．加藤らが用いた実験装置（寒地土木実験所）の全体の様子を図8に示す．
実験に用いられた流量は10年確率洪水（Q10=469.7m3/s）と，30年確率洪水（Q30=680m3/s）であり，流木のスケールは長さ6m，直径30cmのものと，長さ12m，直径30cmの二通りが用いられている．なお，
これらの値はいずれも実スケールの値である．流木の比重はいずれの場合も前述のように1.1（密度1100kg/m3）とされている．実験における流木の投入本数は，長さ6mのものは3600本，長さ12mのもの
は1800本であった．

図6

計画されている流木捕捉工の平面図（b）数値解析の条件
数値シミュレーションでは，加藤らの実験のうち，流木長6mのケースを対象として計算を実施した．
実験では3600本の流木が投入されているが，これをすべて考慮して計算を実行すると，計算機負荷が非常に大きくなる．そこで，投入本数は200本あるいは800本とし，実験よりも少なくして計算を実施した．
ただし，この流木本数を実験よりも少なくした影響については，今後検討が必要といえる．また，一本の流木を構成する要素球体数は5とした．この構成要素数の影響についても，さらに検証が必要である．

― 124 ―

図7

加藤らの実験により最適化された捕捉工の平面形状

図8

加藤らの実験における実験水路全体の様子

計算は実際の河川スケールで行う．また，計算格子の平面形状は図9に示す通りであり，格子セル数
は112×40である．捕捉工のスリットについては，メッシュセルの一部を不透過障害物として設定することで表現した．しかし格子セルサイズがスリット工間隔（= 2m）より大きいため，やむを得ず粗い間隔（= 6m）の障害物として設定した．鉛直方向の格子分割については，10と20の場合について検討を行った．
流木は空中からの落下で供給し，落下による乱れの影響が捕捉工までに概ね消滅する流下距離を考慮して，図9の〇印の位置で投入した．

― 125 ―

図9

数値シミュレーションに用いた計算格子の平面図（1）鉛直方向の格子分割数が20の場合（2）鉛直方向の格子分割数が10の場合
図10

流木捕捉工上流側の横断面内の流況（ベクトルと渦度）（c）第一種二次流の再現状況
湾曲部の流木挙動に大きな影響を与えると考えられる第一種二次流の構造について，そのシミュレー
ションによる再現状況をまず検討する．
数値解析における流況は，二次流も含めてほぼ定常となった．図10は，捕捉工の上流側にあたる，
図9のA-A’断面における断面内の流速ベクトルと主流方向の渦度の分布を重ねて示したもので，図10（1）鉛直方向の格子セルを20分割，図10（2）は10分割した場合の様子である．流況が見やすいよう，鉛直― 126 ―

方向のスケールを水平方向に比べて三倍に拡大して示している．いずれの結果においても，第一種二次流の発生が明確にシミュレートされており，底面付近の流れに着目すると，外岸（右岸）から左岸に向かう流れが確認される．図10（1）では，これに加えて水面付近の右岸側に外岸セル（Outer-Bank Cell:
OBC）と呼ばれる第一種二次流とは逆向きの渦の発生が確認できるが，図10（2）では明確な外岸セルは見られない．これは外岸セルのスケールが比較的小さく，格子解像度が不足する場合にはこれが再現できないことを示唆していると考えられる．これは，特に比重の小さい流木が水面付近を移流する状況を検討する場合には注意が必要と考えられる．

図11

計算結果における流木の挙動（赤い太線が流木を表す）― 127 ―

（d）流木の挙動と堆積状況の検討
図11は，流木の移動の様子と底面付近の流速ベクトルを流木投入地点から下流側について重ねて描いたものである．なお，このケースは，鉛直格子セル数を10とした場合の結果である．図中の時刻は流木投入開始時刻をt=0とした値を示している．
補足工上流側から投入され流木は，底面付近を移動し，囲繞堤により主流を流下するグループと捕捉工に進入するグループに二分される．後者の運動は流木の慣性によるものと推測される．流木捕捉工に進入した流木は右回りの循環流の左岸側を迂回するように進み，一部は循環流内に取り込まれたところ
で流速が遅くなることにより堆積・停止する．このような挙動は，実験で観察された結果とほぼ一致している．
しかしながら，t=120[sec]付近からは，捕捉工に取り込まれた流木の多くが下流側のスリット工を通過して本流に戻っていることがわかる．実験ではスリット工を通過する流木の割合は10%以下であったから，この点は大きく異なる．この理由は明らかで，前述のようにスリット工を計算メッシュ上で障害物として考慮する際，実験よりかなり粗い配置にせざるをえなかったためである．この点の改良は今後の課題である．

図12

加藤らの実験における流木堆積状況

t=180[sec]の図に示した緑色の破線部は，流木の堆積停止域を示している．一方，図12は加藤らの実験における堆積域の写真である．投入本数が異なるため，堆積密度も大きく異なるが，堆積箇所やその形状は概ね一致する．捕捉率について，スリットをすり抜ける流木本数も含めて算出すると，33%程度となった．実験で65-87.5%が捕捉されており，これに比べて相当小さい．この原因として，スリット工付近の流木堆積やこれに伴う水位堰き上げが再現されていないことなどが考えられ，この点は今後の課題である．
（4）iRIC三次元ソルバーNaysCUBEへの本モデルの実装とGUIの改良本研究で構築した流木モデルは，底面摩擦の影響を異方性を含めて考慮した点と，流れに対する抗力に，流木と流れの方向を考慮した点の2つの点に新規性がある．このモデルをiRIC上の三次元ソルバーであるNaysCUBEに実装し，一般に供するため，ソルバー本体の改良を実施するとともに，GUIの改良を行った．
流木と河床との摩擦の考慮については，「計算条件」のメニューの中で，「流木の応用的設定」のタグか― 128 ―

ら行うものとし，図13の赤線で囲んだ部分を新たにメニュー構成に追加した．これには，摩擦力の異方性を反映するため，流木軸方向とそれに直角方向の摩擦係数を別個に設定できるような仕様としている．
一方，流れとの抗力の評価における投影面積の変化を考慮するか否かについては，「計算条件」のメニューから「流木の詳細設定」のタグを選び，「抗力における投影面積」という項目を新たに設け，ここでユーザーが可変とするか，固定とするかを選択できる仕様とした（図14）．

図13 「流木の応用的設定」タグにおける流木と底面との摩擦に関する設定項目の追加

図14 「流木の詳細設定」タグにおける抗力の投影面積に関する設定項目の追加― 129 ―

（5）まとめ
流木捕捉工周辺の沈降性流木挙動について，底面摩擦項等を考慮した三次元流木モデルを構築し，既往の実験の再現を試みた．底面摩擦項には，流木挙動の相違による摩擦力異方性を考慮し，また，抗力における流木軸と流れの角度の影響を考慮したモデル化を行った．
本モデルを加藤らによる流木捕捉工の実験に適用し，実験結果と数値シミュレーション結果の比較を通じてモデルの妥当性を検証した．この結果，数値解析結果は概ね実験結果の特性を再現できたが，ス
リットによる捕捉状況は格子解像度の関係で再現できず，この点は加藤らの数値解析8）と同様の課題を残したといえる．全体捕捉率の再現性にも課題が残り，実験結果における捕捉率を大幅に過小に再現する結果となった．これは，スリット捕捉による水位堰き上げが再現できていないことも一因と考えられる．
今後はこれらの結果をふまえてさらにモデルの精度の向上に努め，さらに多くの条件のもとで適用を行い，モデルの完成度と実用性を高めていきたい．
謝辞
本研究の遂行に当たり，水工リサーチ株式会社，加藤一夫様より実験データをご提供いただくとともに，実験の状況や流木の挙動，数値解析に至るまで，数々の有益なご助言をいただいた．ここに記して深甚なる謝意を表する．
参考文献
1）Schmocker,L.and Weitbrecht,V.:Driftwood risk analysis and eng. measures, J.of Hydr.Eng.,139,683-695,
2013.
2）岡本隆明,山上路生,橿原 : 遊水域を利用したアクティブな流木捕捉システムに関する実験的研究,土木学会論文集B1, 74, 673-678, 2018.
3）加藤一夫,小笠原敏記,松林由里子,渡辺一也,三浦忠昭: 小本川の流木捕捉施設設計に関する水理模型
実験による検討, 河川技術論文集, 24, 137-142, 2018.
4）清水義彦, 長田健吾：流木形状を考慮した個別要素法による橋脚周辺の流木集積過程に関する数値
実験，水工学論文集, 51, 829-834, 2007.
5）初田直彦, 赤堀良介, 清水康行：蛇行流路の流体場と流木の挙動に関する実験と数値解析，応用力学
論文集, 15, 415-422, 2012.
6）長田健吾, 清水義彦：流木群の曲げ変形を伴う堆積過程を説明する数値解析法の開発, 土木学会論文
集B1, 74, 763-768, 2018.
7）Kimura,I. and Kitazono, K.: Effects of the driftwood Richardson number and applicability of a 3D–2D modelto heavy wood jamming around obstacles, Environmental Fluid Mechanics, 20, 503-525, 2020.
8）加藤一夫,サムナー圭,三浦忠昭,菅野貴詳,千葉喜一,井上卓也,清水康行:流木捕捉施設における流木の挙動に関する数値解析,土木学会論文集B1,75,1441-1446,2019.
9）Engelund, F.：Flow and bed topography in channel bends，J. of Hydr. Div. ASCE,100,1631-1648, 1974.
10）Kimura,I.,Onda,S.,Hosoda,T.and himizu,Y.:Computations of suspended sediment transport in a shallow sidecavity using depth-averaged 2D models with effects of secondary currents, J. of Hydro-envir. Res., 4, 153161, 2010.
11） Uchida, T. and Fukuoka, S.: Numerical calculation for bed variation in compound-meandering channel using
depth integrated model without assumption of shallow water flow, Adv. in Water Resources., 72, 45-56, 2014.
12）木村一郎:3D-3Dカップリングモデルを用いた橋脚による流木捕捉過程の再現性,土木学会論文集
B1,75,601-606,2019.

― 130 ―

13）Kang,T. and Kimura,T.: Computational modeling for large wood dynamics with root wad and anisotropic bed
friction in shallow flows, Adv. in Water Res., 121, 419–431, 2018.
14）Kang,T.,Kimura,I.and Shimizu,Y.:Numerical simulation of large wood deposition patterns and responses of
bed morphology in a braided river using large wood dynamics model, Earth Surface Processes andLandforms,45,962–977, 2020.
15）Nelson, J.（他16名）: The international river interface cooperative: public domain flow andmorphodynamics software for education and applications, Adv. in Water Resources., 93, 62-74, 2015.
16）Homepage iRIC: http://i-ric.org/ja/ 2020.
17）Schiller, L., Naumann, Z.: A drag coefficient correlation. VDIZeitschrift 77, 318–320, 1933.
18）越塚誠一：数値流体力学，培風館，pp178-181，1997．
19）Persi E.（他5名）.: Hydrodynamic coefficients of yawed cylinders in open-channel flow. Flow Meas Instrum,65, 288–296, 2019

― 131 ―