iRICにおける固液混相流モデリングの高度化と災害リスク評価への適用

富山大学 理工学研究部

木

村

一

郎

京都大学 防災研究所

竹

林

洋

史

iRICにおける固液混相流モデリングの高度化と災害リスク評価への適用
Refinement of iRIC Solvers for Simulating Multiphase Flow to Assess the Risk of
Natural Disasters
木村 一郎1・竹林 洋史2Ichiro KIMURA and Hiroshi TAKEBAYASHI
1
2

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

要旨

河川災害を考える上で，液体としての河川水のみの流動を主因とする災害はその多くが比較的規模の小さいものに限られ，大災害と呼ばれるものは，土砂，流木などの混相流の動態に起因するものが多い．これらの混相流を伴う災害
は，突発的かつ甚大で，人命に対するリスクも大きい．本研究ではこれらの混相流を伴う河川災害を予測し，そのリスクを評価するための数値解析ツールを
iRIC上のソルバーとして整備し，一般に供することを目的とする．本年度は，次の二つの視点で研究を遂行した．
一つ目は，平面二次元モデルに関するテーマであり，地震動と植生を考慮した土石流・泥流モデルの構築についてである．これについて，平成30年北海道
胆振東部地震によって厚真町吉野地区及び高丘地区で発生した土砂災害について，現地調査の結果を考察するとともに，平行斜面である厚真町吉野地区で発生した泥流と谷地形である厚真町高丘地区で発生した泥流のシミュレーション
結果の比較を行うことで泥流発生場所の違いによる流動特性の違いを考察した．また，地震動による加速度を考慮した泥流の数値シミュレーションを行う
ことにより，地震動が泥流の流動特性に与える影響について検討する．
二つ目は，三次元モデルに関するテーマであり，枝や根が存在する複雑形状の流木，家屋の流出，アイスジャム等，任意形状の剛体が河川流中を流動する現象を再現できるモデルの構築である．既往のNaysCUBEの流木モデルを拡張し，個別要素法による衝突モデルを三次元化するとともに，浮体移流のラグラ
ンジェ方程式の三次元化を実施し，流れのモデルとカップリングした．本モデ
ルをいくつかのテスト計算に適用し，その定性的妥当性を示した．また，今回のモデル改良に対応するため，iRIC-GUIの変更と，それに応じた日本語，英語マニュアルの改訂を実施した．
《キーワード：iRIC；流木；家屋流出；NaysCUBE；土石流；地震；植生》― 37 ―

1． 地震動を考慮した土石流・泥流モデル（担当：竹林洋史）（1）概要
平成30年9月6日，北海道胆振地方中東部の深さ約35kmでマグニチュード6.7の「平成30年北海道胆振東部地震」が発生した．この地震によって223件の土石流，がけ崩れ等の土砂災害が発生し，土砂災害による死者は36名，住宅全壊・半壊は470棟と大きな被害を与えた1）．
土石流・泥流の数値シミュレーションは，土石流・泥流の流速，流動深，流砂量，氾濫範囲などの時
空間的な変化を予測可能である．そのため，数値シミュレーションは，土石流・泥流による流出土砂量，最大流速など，砂防ダムなどの砂防対策施設の設計を行う上での不可欠な情報の提供が可能である．また，ソフト対策においても有効である．土石流・泥流の氾濫範囲を特定し，そこに住む人々に注意喚起を行う，新しく施設や家屋を建てる際に土砂が流れてくる可能性があることを伝える，最大流動深の空間分布から避難方法の提案などが可能である．
既存の土石流・泥流の数値シミュレーションモデルで扱われている基礎式には，地震動による外力は考慮されていない．しかし，平成30年北海道胆振東部地震のように規模の大きい地震は，地震動による揺れが1分以上継続している．そのため，泥流が流下している最中にも地震による振動が続いている可能性が非常に高く，地震によって発生する様々な方向の加速度が土石流・泥流の流動特性に与えた可能性が考えられる．
本研究では平成30年北海道胆振東部地震によって厚真町吉野地区及び高丘地区で発生した土砂災害について，現地調査の結果を考察するとともに，平行斜面である厚真町吉野地区で発生した泥流と谷地形である厚真町高丘地区で発生した泥流のシミュレーション結果の比較を行うことで泥流発生場所の違い
による流動特性の違いを考察する．また，地震動による加速度を考慮した泥流の数値シミュレーションを行うことにより，地震動が泥流の流動特性に与える影響について検討する．
（2）2018年9月に北海道厚真町で発生した土砂災害
a）吉野地区及び高丘地区の概要
本研究で対象としている地区は，図1.1に示すように，北海道南西部の胆振地方に位置しており，吉野地区は厚真町中心街から北東に約4km，高丘地区はさらに北東に約4kmの地点である．多くの人的被害が発生した吉野地区の斜面は，南東向きの平行斜面であり，図1.2（a）に示すように，北東方向に約

新千歳空港
高丘地区

吉野地区
厚真町役場
苫小牧

図1.1

吉野地区と高丘地区の位置― 38 ―

（a）農地に広がった土砂
お椀型の下に凸の地形

流動痕が見られない領域（b）斜面の流動痕
図1.2

吉野地区の土砂の流動状況

1kmの区間のほぼ全ての場所で斜面が崩れた．斜面の高さは40m～90m程度であり，斜面裾部に線状に多くの家屋が存在する．高丘地区の調査地点は，図1.3（a）に示すように，流出土砂が農地に氾濫した
後に河川まで到達し，河道を閉塞させたものであり，多くの谷が密集して流域を形成している谷地形である．渓流内には人家は見られず，渓流出口周辺に勾配の非常に緩やかな農地が広がる．
両地点の年間平均降水量は約1000mmであり，全国平均よりもかなり雨が少ない地域である．また，過去150年，この地域では大きな土砂災害は発生していなかった．
b）厚真町吉野地区の泥流災害の現地調査
現地調査は，TDRによる斜面表層の体積含水率の測定，ドローンによる泥流の発生・流動・堆積地点の撮影を実施した．また，吉野地区ではドローンによる写真測量を実施し，被災後のDSMデータを取得した．
図1.2は，発災4日後の9月10日に撮影した吉野地区の土砂の流動状況である．図1.2（a）に示すように，斜面から発生した土砂が斜面下部に広がる農地に広範囲に氾濫しており，土砂の流動性が非常に高いことがわかる．図1.2（b）は斜面の様子である．斜面には，斜面流下方向に沿って無数の線状の流動痕が見られる．図中に矢印で示したお椀型の地形では斜面内の地盤の低いところに向かって線状の流動痕が― 39 ―

（a）農地に氾濫し，河道閉塞させた様子（b）流域内の斜面の様子
図1.3

高丘地区の土砂の流動状況

見られることから，円弧滑りなどによって崩土がその形状を保ちつつ斜面から離脱したのではなく，土砂がある程度流動化してから斜面を流下したと考えられる．また，図1.2（b）の矢印で示すように，崩壊斜面上端の高さ2～3m程度の領域では，流動痕がほとんど見られないことから，この領域は崩土が斜面から離脱しただけであることがわかる．
図1.3は，発災4日後の9月10日に撮影した高丘地区の様子である．図1.3（a）は，渓流の出口から農地に土砂が氾濫し，河川まで流れて河道を閉塞している様子である．流れてきた土砂は図1.3（a）に示す
ように，農地に広範囲に氾濫していることから非常に流動性が高く，泥流の形態で流れてきたものと考えられる．図1.3（b）は，渓流内の様子を示す．流域内のほぼ全ての斜面で斜面崩壊が発生していることがわかる．崩壊した土砂は渓流に流れ込み，より大きい渓流まで流れて合流し，規模の大きい泥流に発達したと考えられる．
このように，両地点で発生した土砂は高い流動性を示している．しかし，本土砂流出現象は，降雨の無い地震時に発生しており，高い流動性を示す原因が不明である．一方，2016年4月に発生した熊本地
震時には，無降雨時にも係わらず，南阿蘇の山王谷川で数km流れる泥流が発生した2）．山王谷川の流域― 40 ―

43.6

43.1

31.9
42.3
39.7
45.3
38.2

図1.4

斜面表土の体積含水率（赤点は斜面，黄点は斜面下部の堆積土砂）

図1.5

厚真町吉野地区の斜面の浸食深分布の表土は，阿蘇山からの火山噴出物である1mm以下の粒径の細かい土砂で覆われており，黒ぼくと呼ばれている．山王谷川では，晴天が数日継続しても表土の体積含水率が40%を超えており，表土の保水性が非常に高い2）．厚真町の吉野地区や高丘地区の表土も，樽前山などからの火山噴出物である1mm以下の粒径の細かい土砂を多く含んでおり，山王谷川と同様に体積含水率が高いと考えられた．そのため，
日雨量2mmの降雨後に晴天が3日続いた10月19日に日当たりの良い斜面の表土の体積含水率を測定した．図1.4は，測定された体積含水率の値を示し，赤点は斜面の表土，黄点は流出して斜面下部の平坦地に堆積した土砂の体積含水率の値を示す．斜面の表土は，上部，下部に関係なく，全ての場所で体積含水率が高い値を示しており，間隙に多くの水が含まれていることがわかる．
以上の結果から，地震によって震動した斜面において，斜面上部において小規模な斜面崩壊が発生するとともに，斜面上で液状化も発生し，流体化した土と水の混合物が泥流として斜面及び渓流を流下した可能性が高いと考えられる．
図1.5は，ドローンを用いた写真測量によって得られたた地形を示す．地形の表面には，ドローンを用いた写真測量によって得られたDSMと国土地理院による被災前の5mDEMとの差分をカラーコンター
で示している．なお，水平方向の精度は約1m，鉛直方向については，国土交通省北海道開発局によるレーザー測量結果との差分を真値とした場合，差分が0.9mと最も大きかったのは14mほど浸食された斜― 41 ―

面部であった．つまり最大で約7％の誤差となっている．図に示すように，崩壊斜面上部の一部では深い浸食深を示しているが，崩壊深の空間分布は非常に小さく，ほぼ同じ厚さで斜面が浸食されていることがわかる．なお，斜面上部は被災後も樹木が存在するため，大きな値となっている．
（3）数値シミュレーションの概要
数値シミュレーションに用いた基礎方程式は，土砂と水の混合物を一つの流体として扱った平面二次
元の一流体連続体モデルとして導出したものである．構成則は土石流・泥流を対象として得られた江頭ら3）によるものを応用して用いる．本解析モデルでは，粒径の細かい土砂が水と混合することによって液体として振る舞う相変化を考慮している．
土石流・泥流が発生する斜面や流動する渓流は，様々な形状を有している．そのため，斜面形状や渓流形状に近い形の解析領域を用いると，解析格子数が少なくなるとともに，土石流・泥流境界形状が格子形状の影響を受けづらくなり，解析精度及び解析効率が高くなる．そのため，非直交曲線格子を用いるとともに，デカルト座標系で記述された平面二次元の支配方程式を一般座標系の支配方程式に書き換
え，一般座標系に変換した基礎方程式を離散化して数値シミュレーションに用いる．
土石流，泥流などのマスムーブメントは，sub-critical flowとsuper-critical flowが混在した流れとなる．
そのため，sub-critical flowとsuper-critical flowの混在場を安定的に解析可能なTVD-MacCormackスキー
ムで基礎方程式を離散化する．
a）地震動を考慮した泥流の基礎方程式
解析に用いた基礎方程式を以下に示す．デカルト座標系と一般座標系の変換に用いるメトリックスは
以下の様である．

J1 x y x y       


（1.1.1）


y
J
x


（1.1.2）


y
 J
x


（1.1.3）


x
 J
y


（1.1.4）


x
J
y


（1.1.5）

ここに，ξとηは一般座標系における流下方向と横断方向の座標であり，xとyはデカルト座標系での座標である．
水及び土砂の混合物の質量保存則は地盤の間隙の水による飽和度Srの変化を考慮した以下の式を用いる2）．

S   h    hU    hV  E 
1  1  c*  1  r      t  J    J    J  c* J  100  

― 42 ―（1.2）

ここに，tは時間，zは水位である．表面流の水深はhである．UとV はξとη方向における反変速度成分であり，以下のように定義される．



u v
x
y


V
u v

x
y
U


（1.3）
（1.4）

ここに，uとvはxとy方向における流速成分，c*は静止堆積層の砂礫の堆積濃度，Eは河床の浸食速度であり，以下の式3）を用いる．

E c* tan    e u  v22（1.5）

ここに，θは流動方向の河床勾配であり，以下の関係がある．

sin  u sin  x  v sin  y
u 2  v2

（1.6）

ここに，θxはx方向の河床勾配，θyはy方向の河床勾配である．θeは水及び土砂の混合物中の鉛直平均土
砂濃度 c に関する流動方向の平衡河床勾配であり，以下の関係がある．

    1 c hstan etan sh   1 c  1
 

（1.7）

ここに，φsは土砂の内部摩擦角である．流体中の土砂の質量保存則は以下のようである．

  ch    chU    chV  E   t  J    J    J  J（1.8）

運動量保存則は以下のようである．

  hU    hU    hU 


U

V

J   J 
t  J   hu          hv              U    V
U    V
J    x   x   J    y   y  
 1    2   2   z 1        z hb  a cos   g  az  h         b   J   x   y     J  x x y y    
J

 

2
2
1  1         P 1       P  b        m  J   x   y     J  x x y y     m J― 43 ―（1.9）

  hV    hV    hV 


U

V

t  J    J   J hu          hv              U    V
U    V
J    x   x   J    y   y  22
 1        z
1          zb hb
  a cos   g  az  h         J  x x y y    J   x   y     
J




2
2
1  1       P 1         P  b

  

      m  J  x x y y    J   x   y      m J

（1.10）

ここに，gは重力，zbは河床位である．azは鉛直方向の地震加速度，aξとaηはξとη方向における反変地震
加速度成分であり，以下のように定義される．


ax  a y
x
y


a
ax  a y

x
y

a


（1.11）（1.12）

ここに，axとayはxとy方向における地震加速度成分である．Pは圧力であり，以下の関係を用いる．
hP  m g  h  z  cos dz（1.13）

0

ここに，zは鉛直方向の軸である．ρmは以下の関係がある．

m       c  （1.14）

ここに，ρは間隙流体の密度，σは土砂の密度である．本解析では，細粒土砂の相変化を考慮する．細粒
土砂は水と混合して運動すると液体のように振る舞う．そのため，本研究では0.2mm以下の細粒土砂は液化するとし，間隙流体の密度を水と0.2mm以下の細粒土砂の混合物の密度とした．例えば，静止堆積
濃度が0.6で0.2mm以下の土砂が40%含まれている場合，土砂の密度が2.65g/cm3であれば，泥流の土砂濃度は0.36となり，間隙流体の密度は1.62g/cm3となる．これらの影響は，平衡勾配の値などの現れ，平衡勾配が小さくなるため，泥流が低勾配領域まで輸送されることとなる．
τbξとτbηはξとη方向における反変掃流力成分であり，以下のように定義される．


bx  by
x
y


b 
bx  by
x
y
b 

（1.15）（1.16）

ここに，τxとτyはxとy方向における掃流力成分であり，以下のようである．

bx   f  u  v  u u v

by   f  u  v  u v v

y

y

b

b

2

2

2

2

2

2

2

― 44 ―

2（1.17）（1.18）

図1.6

吉野地区の数値シミュレーションに用いた斜面崩壊の形状と深さ（単位：m）

図1.7

高丘地区の数値シミュレーションに用いた斜面崩壊の形状と深さ（単位：m）

ここに，τyは降伏応力であり，以下の関係3）を用いる．
1
n

c

 y        cgh cos  tan s c* 

（1.19）

ここに，nは常数であり，ここでは5を用いる．fbは抵抗係数であり，土石流の場合の以下の関係を用いる．

f
72 2 ，  b25  1  c
fb

k f
4 
c

2
3


6

（乱流域）

  k  1  e c  h 
5
3

d





2

  d 
1
3






2
3）
（層流域）

ここに，kf=0.16，kd=0.0828，eは粒子の反発係数，dは土砂の平均粒径である．

― 45 ―（1.20）（1.21）

図1.8

厚真町鹿沼で観測された地震加速度の時間変化4）

河床位方程式は以下のようである．

  zb E  
c* J
t  J （1.22）

ここに，zbは河床位である．
b）解析条件
吉野地区と高丘地区を対象に土砂の流動過程の数値シミュレーションを実施し，土砂の流動特性を検討する．北海道胆振東部地震によって，北海道安平町のK-NET追分観測点で最大加速度1796gal（三成分合成），気象庁厚真町鹿沼震度観測点で最大加速度967gal（三成分合成）など大きな加速度を観測した．
この数値は2016年4月16日の熊本地震本震の最大加速度であるKiK-net益城観測点の1580gal（三成分合成），また厚真町と同じく地震によって泥流が発生した南阿蘇村にある気象庁南阿蘇村中松震度観測点の855gal（三成分合成）よりも大きな値である4）．
図1.6に吉野地区の数値シミュレーションに用いた斜面崩壊の形状と深さを示す．図1.2（b）に示すよ
うに，斜面浸食域上端の約1m～2mの地点は明確な流動痕が確認できないため，この領域を斜面崩壊地点として設定した．斜面崩壊は，斜面横断方向に連続的に発生しているが，全ての場所で同時に斜面崩壊が発生したのかどうかは不明である．ここでは図1.6に示すように，斜面を横断する細長い斜面崩壊形状を設定し，全ての場所で同時に斜面崩壊が発生し，水と土砂の混合物が速やかに泥流化して斜面を流下する条件で数値シミュレーションを実施した．なお，斜面崩壊の深さは，ドローンによって計測したDSMの地盤高さと国土地理院の被災前のDEMの差を元に決定した．図7に高丘地区の数値シミュレー
ションに用いた斜面崩壊の形状と深さを示す．高丘地区ではDSMが得られていないが，ドローン映像
から斜面崩壊発生箇所と斜面崩壊深さを設定した．
図1.8に気象庁による厚真町鹿沼で観測された地震加速度の時間変化5）を示す．加速度の大きな値の発生時間は約20秒間である．一方，後述するように，吉野地区では斜面崩壊発生から土砂が宅地に到達するまで10秒未満であり，農地に氾濫するプロセスまで含めても30秒程度の現象である．そのため，土砂の斜面流下及び農地での氾濫時にも地盤は振動している可能性が高い．本解析では，地震加速度が最大に近いときに斜面崩壊が発生したと仮定し，図1.8の8.5秒に斜面崩壊を発生させる．
解析は4つの条件で実施した．Case1～Case3は吉野地区，Case4は高丘地区を対象としたものである．
Case1は地震動を考慮した解析であり，地震動の周期は実測された約0.5秒となっている．Case2は地震動を考慮していない．Case3は阪神淡路大震災などで発生した比較的周期の長い地震動を想定し，Case1の地震動の波長を4倍にした約2秒周期の地震動を考慮した．
― 46 ―

（a）2秒後（b）6秒後（c）18秒後
図1.9

数値シミュレーションによる吉野地区の泥流の深さの時間変化（単位：m）（4）解析結果
a）吉野地区の泥流の流動特性
図1.9にCase 1の泥流の深さの時間的な変化を示す．図に示すように，泥流は時間とともに規模を拡大させながら斜面を流下し，約6秒で家屋に到達している．また，家屋到達時の深さは，家屋衝突による迫り上がりなどもあり，約6mと非常に深くなっている．斜面を流下する泥流の平均流速は13.5m/sとなっている．流出した泥流は農地に広がっており，概ね氾濫範囲は評価できている．本土砂災害は地震起源であり，事前避難が困難であること，斜面崩壊発生からわずか6秒で宅地に泥流が到達していることを考えると，泥流発生後の避難も困難である．そのため，水分保持機能の高い細粒の降下火砕物で構成されている表土を有する斜面付近に家屋を建設することはできる限り避けることが重要と考えられる．

― 47 ―

（a）地震動考慮（Case 1，周期0.5秒）（b）地震動無し（Case2）（c）地震動考慮（Case3，周期2秒）
図1.10

地震動が吉野地区の泥流の最大深さに与える影響（単位：m）（a）Case1（b）Case2
図1.11 3秒後の速度ベクトルの平面分布（単位：m/s）― 48 ―

（b）150秒後（a）30秒後（c）200秒後
図1.12

数値シミュレーションによる高丘地区の泥流の深さの時間変化（単位：m）

b）地震動の有無による泥流の流動特性の違い
図1.10に吉野地区の泥流シミュレーションで得られた最大流動深の平面分布を示す．図1.10（a）と（b）の比較より，地震動が泥流の流動範囲に与える影響は非常に小さいことが分かる．つまり，瞬間的な泥流の動態は地震動の影響を強く受けているが，最終的な流動範囲はそれぞれの方向の地震動の影響が打ち消し合い，地震動の影響が非常に小さくなったものと考えられる．しかし，各時刻の流速ベクトルには違いが見られた．図1.11は3秒後の速度ベクトルの平面分布であり，Case2に比べ，Case1は左向きのベクトル成分が大きい．家屋などの破壊については，最大流速が非常に重要なパラメータになることを考えると，地震動による流速の増加で家屋が破壊されやすくなっていることがわかる．
地震動の周期を約2秒に伸ばした場合の泥流の最大深さの平面分布を図1.10（c）に示す．Case1及びCase2と比較すると，Case3の方が流動範囲は狭いことがわかる．また，周期約1秒，周期約3秒のCaseについても同様にシミュレーションを行った結果，周期が長い方が流動範囲に与える影響が大きくなった．
つまり，地震の周期が長い方が地震動が泥流の流動範囲に与える影響が大きくなることを示す．これは，周期を伸ばした場合，一方向に加速度が作用する時間が長くなることが原因である．

― 49 ―

c）高丘地区の泥流の流動特性
図1.12にCase 4の泥流の深さの時間的な変化を示す．図に示すように，各支川から発生した泥流が合流し，本川をゆっくり下流に流下していることがわかる．農地に到達するまでの時間は150秒であり，農地氾濫時の流速は3m/sとなっており，吉野地区に比べて非常に流速が遅いことがわかる．これは，渓流本川の平均河床勾配が約4度，農地の地盤勾配が約1度と非常に緩やかであるためである．解析対象
流域の南隣の渓流でも同様の泥流が発生している．南隣の渓流では，図1.12（a）に赤丸で示すように，渓流出口に存在していた家屋が被災し，1名の方が犠牲となった．上述のように，高丘地区では，斜面崩壊発生から農地や家屋が存在する渓流出口に泥流が到達するまでの時間は吉野地区に比べて非常に長いため，地震発生後に避難を開始しても間に合う可能性がある．また，泥流が農地などに流れてくる場
所は渓流出口で固定されているため，渓流出口付近から泥流が流出する方向の直角方向に数十メートル離れた場所に家屋を建設すれば被災する可能性が非常に低いことがわかる．つまり，同じ厚真町内でも平行斜面と谷地形では異なる防災・減災対策や避難方法を考える方が有効と思われる．
（5）本章のまとめ
平成30年北海道胆振東部地震によって厚真町吉野地区及び高丘地区で発生した土砂災害について，現地調査の結果を考察するとともに，平行斜面である吉野地区で発生した泥流と谷地形である厚真町高丘地区で発生した泥流のシミュレーション結果の比較を行うことで泥流発生場所の違いによる流動特性の違いを考察した．また，地震動による加速度を考慮した泥流の数値シミュレーションを行うことにより，地震動が泥流の流動特性に与える影響について検討した．その結果，同じ厚真町内でも平行斜面と谷地
形では異なる防災・減災対策や避難方法を考える方が有効であることが示された．また，胆振東部地震による加速度の影響を考慮した泥流は地震動を考慮していないものと比べて，速度ベクトルに違いが見られたが，流動範囲については違いが見られなかった．しかし，地震周期を約4倍に伸ばすと，泥流の流動範囲に違いが見られた．このことから，比較的周期の長い地震動は泥流の流動範囲に影響を与えることが示された．
参考文献（第1章分）
1） 内閣府：平成30年北海道胆振東部地震に係る被害状況等について（平成30年10月29日17時30分現在）
2） 竹林洋史・藤田正治：2016年 4月熊本地震時に発生した山王谷川の土石流の流動特性山王谷川の土石流の流動特性 , 土木学会論文集B1（水工学）Vol.74, No.5, I_1093-I_1098, 2018.
3） 江頭進治・伊藤隆郭：土石流の数値シミュレーション, 日本流体力学会数値流体力学部門Web会誌,
第12巻, 第2号, pp. 33-43, 2004．
4） 防災科研K-NET, KiK-net（NIED K-NET, KiK-net）
5） 気象庁：長周期地震動に関する観測情報，2018年9月6日― 50 ―

厚真町鹿沼

2． 植生を考慮した土石流・泥流モデル（担当：竹林洋史）（1）概要
2018（平成30）年7月6日と7日に中国・四国地方を中心として，九州北部から中部地方の広い範囲で多くの雨が降った．広島県においても広島市や呉市など県南部を中心に長時間雨が降り続けるとともに，南西から北東方向に線状降水帯が形成され，降雨強度の強い雨が観測された．その結果，非常に多くの斜面崩壊および土石流が発生し，甚大な被害が発生した．広島市では，わずか4年前の2014年8月にも多くの斜面崩壊および土石流による土砂災害が発生している．そのため広島では，土砂災害警戒区域の指定が2014年以降に急速に進むとともに，土石流による流出土砂量の推定方法が変更されるなど，
様々な対策が進められているところであった．
本研究では，2018年7月6日に広島県安芸郡熊野町川角で発生した土石流について実施した現地調査の結果を報告するとともに，植生の抗力を考慮した土石流の数値シミュレーションを実施し，斜面における土石流の流動特性と宅地における土石流の氾濫特性について検討した．
（2）広島県安芸郡熊野町川角で発生した土石流
a）対象地域の概要
広島県安芸郡熊野町川角地区は，図2.1のように，広島市と呉市の境界の三石山（みついわやま）の北西斜面の下に広がる宅地である．「三石山」と呼ばれているように，山頂付近や斜面には，球状風化
した直径数mの巨礫が点在しており，基岩が露出している場所も多い．山頂から宅地までの斜面は25°程度の急斜面であるが，宅地内は5°程度と斜面部よりも緩やかとなっている．

川角地区

図2.1

広島県安芸郡熊野町川角の場所― 51 ―

背景：Google Earth（b）数値シミュレーションによる土石流の最大
流動深
斜面上の土石流の流下経路と数値シミュレーションによる土石流の最大流動深の比較（単位：m）（a）斜面上の土石流の流下経路

図2.2

(a)
0

100

(d)

200 (m)

N

(b)

(e)
(c)

図2.3

土石流による渓流の浸食・堆積の様子

b）現地調査
図2.2（a）に被災後の斜面の様子を示す．図より，山頂付近の二カ所と中腹付近の一カ所から斜面崩壊によって土石流が発生していると判断される．山頂付近の二カ所から発生した土石流が流れた渓流
は，中腹付近で合流するとともに分派し，二つの渓流に分かれてから宅地に流れ込んでいる．分派して― 52 ―

(a)

(b)

(c)

0

100

200 (m)

N

(e)

(d)

(f)

(g)

図2.4

(h)

土石流による浸食・堆積の様子

いる二つの渓流の内，南側（図2.2（a）の写真では右側）の渓流では前述の中腹付近で発生した土石流も流れ込んでいる．
図2.3と図2.4に，現地調査時の様子を示す．図2.3（d）と（e）からも山頂付近の二カ所から発生した土石流が流れた渓流が中腹付近で合流・分派していることがわかる．また，二つの渓流の合流点付近には直径数mの巨礫が存在することもわかる．図2.3（c）に示すように，北側の渓流から宅地に流れ込んだ土石流は，宅地に流入し始めた地点では既に堆積傾向となっている．直径数mの巨礫とともに，粒径の小さい真砂土が多く堆積しており，粒度が広いことがわかる．図2.3（a）と（b）は南側の渓流の様子である．
図2.3（a）に示すように，南側の渓流を流下した土石流は，宅地の南側を流下するとともに，宅地の南端から宅地に流入している．土砂が流れ込んだ領域は既に宅地用に造成されていたが，家屋が建設されていない更地が多い．そのため，敷地内に土砂が流れ込んでいる宅地もあるが，家屋の全壊などは見られない．また，渓流内に堆積している土砂は真砂土で構成されており，非常に粒径が細かい．
図2.4（a）は，宅地内の地盤勾配方向に沿って建設された北西方向に伸びる道路の下端から道路の様子を撮影したものである．北側の渓流から宅地に流入した土石流の主流は，この道路を流下している．
図2.4（d）は，図2.4（a）で示した道路および図2.4（f）で示す渓流から流れてきた土石流の堆積物である．
直径数mの大きい岩とともに，粒径の小さい真砂土も多く含まれていることがわかる．さらに，流木も含まれている．図2.4（f）は，南側の渓流の下流端付近の様子であり，図2.3（a）で示したように，渓流内に堆積している土砂の粒径は非常に細かく，真砂土で構成されている．図2.4（b）は，図2.4（a）で示した道路の途中の様子であり，道路沿いの家屋の一階の一部が被災していることがわかる．図2.4（c）― 53 ―015

0

30 (km)15（a）2018年7月6日18:50の降雨強度（b）2018年7月6日20:00の降雨強度
最大降雨強度
18:50

140
120

時間換算雨量（mm/h）

30 (km)NN

土石流発生時刻
20:00

100
80
60
40
20

12:00
12:12
12:24
12:36
12:48
13:00
13:12
13:24
13:36
13:48
14:00
14:12
14:24
14:36
14:48
15:00
15:12
15:24
15:36
15:48
16:00
16:12
16:24
16:36
16:48
17:00
17:12
17:24
17:36
17:48
18:00
18:12
18:24
18:36
18:48
19:00
19:12
19:24
19:36
19:48
20:00
20:12
20:24
20:36
20:48
21:00
21:12
21:24
21:36
21:48
22:00
22:12
22:24
22:36
22:48
23:00
23:12
23:24
23:36

0

時間

（c） 2018年7月6日の降雨強度の時間変化
図2.5

CXレーダーによる降雨強度の時空間変化に示すように，北側の渓流から宅地に流れ込んだ土石流は，渓流出口から北西に広がる家屋を破壊するとともに，その主流は，図2.4（e）で示す南西方向に伸びる道路に沿って流れている．図2.4（a）で示した道路と図2.4（e）で示した道路が交わる場所の様子を示したものが図2.4（g）である．図2.4（g）に示されている家屋は，一階は形状が全くわからないほど破壊されている．しかし，転倒しているものの，二階はそれほど大きな損傷を受けていないように見える．図2.4（a）と図2.4（b）も同様に，一階は損傷しているが，二階の損傷は比較的軽度であることがわかる．これは，土石流の流動深が二階を直撃するほどの深さではないことを示すものである．
図2.4（h）は，南北の二つの渓流に挟まれた領域であるが，図2.4（e）で示す南西方向に伸びる道路よりも地盤が高い位置にある．そのため，敷地内に土石流は流れ込んで樹木は倒れているが，家屋を全壊するほどの勢いでは土石流が流れていない．
c）降雨条件図2.5（c）は，CXバンドレーダーによって得られた熊野町川角における2018年7月6日の降雨強度の時間変化である．降雨強度は14時ごろから強くなり，18時50分に最大降雨強度を記録し，120mm/hを超える非常に強い雨となっている．一方，宅地に流れ込んだ土石流が発生したのは1時間10分後の20時ごろである．20時の時点では，降雨強度は10mm/h程度となっており，雨はかなり弱まっていたことがわ
かる．図2.5（a）に最大降雨強度が観測された18時50分におけるCXバンドレーダーによって得られた熊野町川角周辺の降雨強度の分布を示す．また，図2.5（b）に土石流が発生した20時におけるCXバンドレ― 54 ―

ーダーによって得られた熊野町川角周辺の降雨強度の分布を示す．図2.5（a）に示すように，最大降雨強度が観測された時間帯は南西から北東方向に伸びる線状降水帯が形成されている．図5（c）に示すよ
うに，最大降雨強度観測後も強い雨が約1時間程度断続的に降っていることがわかる．なお，土石流が発生した20時の時点における総降水量は153mmであり，2014年に広島市で発生した斜面崩壊・土石流の多くが総降水量200mm以上の場所で発生したことを考えると1），少し小さめの値である．このように，雨が弱まった時点で斜面崩壊・土石流が発生しており，住民としては，斜面崩壊・土石流の危険が過ぎ去ったと感じていた時間帯となっていた可能性が考えられる．
（3）数値シミュレーションの概要
a）基礎方程式
本解析では，「1．地震動を考慮した土石流・泥流モデル」の運動方程式（式（9），（10））を修正したものを用いる．つまり，本解析では地震動の影響は無視し，代わりに植生による抗力項を追加する．解析に用いた運動方程式は以下の様である．

  hU    hU    hU 


U

V

t  J   J   J hu          hv              U    V
U    V
J    x   x   J    y   y  
 1     2   2   z
1        zb  gh         b   J   x   y     J  x x y y    

 


（2.1）

2
2
F
1  1         P 1       P  b

 

 v
       m  J   x   y     J  x x y y     m J m J  hV    hV    hV 


U

V

t  J    J   J hu          hv              U    V
U    V
J    x   x   J    y   y  22
 1        z
1          zb 
b




 gh   J  x x y y    J   x   y     




2
2
F
1  1       P 1         P  b  
 v
  



m  J  x x y y    J   x   y      m J m J



（2.2）

FvξとFvηはξとη方向における植生の反変効力成分であり，以下のように定義される．

Fv




Fvx  Fvy
x
y

（2.3）

Fv




Fvx 
Fvy
x
y

（2.4）

ここに，FvxとFvy はxとy方向における植生の効力成分であり，以下のようである2）．

― 55 ―

Fvx  Fv

Fvy  Fvu2u  v2v
2u  v2Fv 1
 Cdv  v  u 2  v 2  hw
 2

（2.5）

（2.6）

（2.7）

ここに，Cdv（≒1.0）は形状効力係数である．λvは植生密生度であり，以下のように，単位面積あたりの植生の本数nvと植生の平均直径dvを用いて，以下のように評価される．

 v nv d v（2.8）

地盤が植生流送深さ以上に浸食された場合（本解析では30cm），または基岩が露出した場合は，植生が流送されたと判断し，密生度を0とする．hwは植生の浸水深さであり，抜水時は水深，冠水時は植生高さとなる．なお，本解析対象渓流では，樹木の高さは土石流流動深さより十分高いため，全て抜水状態として取り扱う．
b）解析条件
地形データは国土地理院により被災前に測定されたものであり，平面的に5mの解像度のDEMデータである．解析格子は，道路上を流れる土石流を表現するとともに家屋の形状を表現するため，2m×2mの比較的細かい正方形格子を用いる．土石流は斜面崩壊を発端として発生するものとした．そのため，
本解析の境界条件は流体化した水と土砂の混合物が初期の斜面に初速ゼロで存在するというものであり，土石流の流量などを上流から与えていない．なお，間隙率を0.4とし，間隙は水で飽和した条件としている．また，現地の材料を参考に，土砂の平均粒径は1cm，液化する細粒土砂の存在率は0.14とし，空間的に一様に与えた．表層崩壊は，山頂付近にニカ所，中腹付近に一カ所の合計三カ所を考慮する．
表層崩壊の場所は，被災後の調査によって判断された崩壊地とした．これら三カ所の斜面崩壊の発生時
刻の関係は不明であるため，本解析ではこれら三カ所の斜面崩壊が同時に発生したと仮定している．なお，下流域での土石流の規模は土石流の流動過程で渓床や河岸などを浸食して取り込んだ土砂と水の量で決定されるため，本土石流のように土石流の流動距離が長い場合は斜面崩壊の大きさは土石流の規模にほとんど影響を与えない3）．宅地の中の家屋は，高さ8mの非浸食域とした．現地調査及び広島県の調査4）によると，渓流では，局所的に深く浸食された地点もあるが，土石流発生前から基岩が露出していた領域も多いことから，平均の不安定土砂厚として0.35mを用いた．なお，宅地内の道路が大きく浸食・流出した場所は見られなかったことと，宅地内の地盤勾配が緩やかで堆積が卓越することから，宅地内の地盤は非浸食とした．
解析は，3つの条件で実施した．Case 1は植生を無視した解析，Case 2とCase 3は植生を考慮した解析であり，Case 2は現地の値とほぼ同じ植生密生度である0.05（1/m）を用いた解析，Case 3は大きめの植生密生度として1（1/m）を用いた解析である．土石流現象は，低平地の河川流と比べて，河床材料が大きく流木が多い．そのため，流木や岩石が密生した樹木内を物理的に通過できない場合がある．そのため，土石流の樹木の透過しづらさの一部を評価するため，Case 3として大きめの密生度を用いている．

― 56 ―

0100

200 (m)N

図2.6（a）5秒後（b）30秒後（c）60秒後（d）100秒後（e）150秒後（f）200秒後

広島県安芸郡熊野町川角で発生した土石流数値シミュレーション による土石流の流動深（Case 1, 単位：m）（4）結果と考察
図2.2（ b）にCase 1の数値シミュレーションによって得られた最大流動深の空間的な分布をGoogle
Earthに出力し，鳥瞰図として表示したものを示す．現地調査によって得られた被災後の斜面の様子を示す図2.2（a）と比較すると，土石流による斜面の浸食域と数値シミュレーションによる土石流の主流
動域がよく一致していることがわかる．しかし，一部，現地では流動が確認されていない植生域にも土石流が流れている．なお，解析で得られた宅地への流出土砂量は12,369 m3であり，LPデータに基づいた広島県の調査値4）である12,100 m3と同程度の値である．図2.6に数値シミュレーションによって得られた宅地における最大流動深の空間的な分布と現地調査によって得られた被災後の宅地の土砂の氾濫の― 57 ―

0100

200 (m)N（a）Case 1（植生無し）（b）Case 2（植生密生度：0.05（1/m））（c）Case 3（植生密生度：1（1/m））
図2.7

植生の密生度が最大流動深さに与える影響（単位：m）

様子を比較したものを示す．北側の渓流からの土石流が宅地に流れ込んだ地点では，数値シミュレーシ
ョンでも流動深が深くなっていることがわかる．また，現地で土砂が氾濫している範囲は，数値シミュ
レーションで土石流が流動した範囲とよく一致していることがわかる．さらに，図2.4（b）と（g）に示― 58 ―

0100

200 (m)N

図2.8（b）5秒後（b）30秒後（c）60秒後（d）100秒後（e）150秒後（f）200秒後

広島県安芸郡熊野町川角で発生した土石流数値シミュレーション による土石流の流動深（Case 2, 単位：m）

す家屋の一階が被災した領域の土石流の深さは，それぞれ，最大で2.08m ，2.77mであり，一階のみが大きく被災していることと一致している．これらの結果より，一部，現地では流動が確認されていない
植生域にも土石流が流れているものの，植生が無い宅地における土石流の流動・氾濫特性は，本数値シミュレーションによって現地の土石流の流動・氾濫特性をある程度表現できていると考えられる．
数値シミュレーションによると，北側渓流の宅地流入時の土石流のピーク流量は804m3/s，南側渓流の宅地流入時の土石流のピーク流量は487m3/sであり，北側渓流の土石流の方が規模が大きいことがわ
かる．なお，合理式で簡便にピーク流量を算定すると，1m3/s ～ 2m3/sであり，対象豪雨中の降雨による表面流は土石流の流量にほとんど影響を与えないことがわかる．

― 59 ―

0100

200 (m)N

図2.9（c）5秒後（b）30秒後（e）60秒後（d）100秒後（e）150秒後（f）200秒後

広島県安芸郡熊野町川角で発生した土石流数値シミュレーション による土石流の流動深（Case 3, 単位：m）

図2.6にCase 1の土石流の深さの時空間的な変化を示す．山頂付近から発生した二つの土石流は，別々の渓流を流下し，斜面中腹付近で合流するとともに，南と北の二つの渓流に分派して流れている．南側の渓流は，もう一つの中腹付近で発生した斜面崩壊による土石流が流下した渓流と同一である．図6に示すように，数値シミュレーションにおいても宅地に流れ込み多くの家屋を破壊した土石流は北側の渓流からの土石流であることがわかる．
数値シミュレーションによると，山頂付近の崩壊発生から北側の渓流の出口付近の家屋に土石流が到達するまで約40秒となっており，斜面における土石流の平均速度は約11m/sである．山頂付近の崩壊による土石流が宅地に到達した時点の土石流の深さは約2mで，家屋に衝突した時に3m以上となっている．

― 60 ―

図2.7に数値シミュレーションによって得られた最大流動深の平面分布を示す．Case 1による最大流動深の空間的な分布と比較すると，Case 2とCase 3において，赤丸で示した樹林域において，土石流の氾濫範囲が狭くなっていることがわかる．また，Case 1よりもCase 3の方が，最大流動深が全体的に深くなっていることがわかる．これは，植生の抵抗の存在によるものである．ただし，黄丸で示された地点については，Case 1よりもCase 3の方が，最大流動深が浅くなっている．これは，Case 3では植生の抵抗によって流速が低減したため，分派点の高地盤域や宅地などに土石流が衝突した時のせり上がりが少ないためである．
宅地内の土石流の氾濫域や流動深も3つのケースで若干異なっているが，植生の影響はほとんど見られないことがわかる．
図2.8と図2.9にCase 2とCase 3の数値シミュレーションによって得られた流動深の時空間的な分布を示す．図2.6のCase 1による流動深と比較すると，植生の抵抗によって，土石流の宅地への到達時刻が遅くなっていることがわかる．これは，渓流に土石流センサーを設定している場合を考えると，数秒程度であるが，植生を考慮していない解析では危険側に評価していることがわかる．
（5）結論
広島県安芸郡熊野町川角で発生した土石流について，現地調査と数値シミュレーションを実施し，宅地における土石流の氾濫特性と土砂災害警戒区域との関係を検討した．現地調査により，斜面の三カ所から土石流が発生し，合流・分派した渓流を土石流が流下し，北側の渓流から宅地に流れ込んだ土石流によって多くの被害が発生したことが明らかとなった．流出した土砂は，粒径の細かい真砂土から直径数mの花崗岩も含まれ，粒度が非常に広いことが明らかとなった．数値シミュレーションにより，山頂付近の崩壊発生から北側の渓流の出口付近の家屋に土石流が到達するまでわずか40秒程度となっており，斜面における土石流の平均的な速度は約11m/sであった．土石流は宅地内を数10cm～4m程度の深さで流下・氾濫し，多くの家屋が全壊・半壊した．土砂災害警戒区域内には土砂が流れてこない領域も多いこと，土石流の深さは多くの場所で2m以下であること，豪雨中に避難のために長い距離の移動が困難であることを考えると，土砂災害警戒区域内の土石流が流れてこない領域への避難や二階以上への鉛直避難は，土砂災害からの生存確率を上昇させることが期待されることが示された．
斜面及び渓流における植生の存在は，斜面における土石流の流動域を狭めるとともに，土石流の流動深を深くする．また，植生の抵抗によって土石流の伝播が遅くなり，土石流発生から土石流が宅地に到達するまでの時間は長くなる．さらに，植生を考慮すると，土石流が宅地に衝突した時の土石流のせり
上がりは小さくなることが示された．
参考文献（第2章分）
1） 松四雄騎：2014年広島豪雨災害時の斜面崩壊・土石流について（速報その2：降雨と崩壊の分布），
http://www.slope.dpri.kyoto-u.ac.jp/disaster_reports/ 20140820Hiroshima/201408Hiroshima_Rep2.html，
2014.
2） 清水義彦，辻本哲郎：植生対を伴う流れ場の平面2次元解析，水工学論文集，第39巻，pp.513-518，
19953） 広島県：平成30年7月豪雨災害を踏まえた今後の水害・土砂災害対策のあり方検討会（第1回砂防部会）土砂災害警戒区域等における検討事項（参考資料），2018.
4） 広島県砂防課：土砂災害ポータルひろしま，https://www.sabo.pref.hiroshima.lg.jp/portal/top.aspx.

― 61 ―

3． 河川流中の任意形状剛体の三次元挙動と固定構造物との衝突のモデル（担当：木村一郎）（6）本章の目的と概要
近年，国内外で流木によるさまざまな災害が問題となっている．このような現象を再現する上でも数値シミュレーションは極めて有効であるが，初期の流木モデルは根や枝を考慮せず，単に棒状の単純円柱形状を模擬したものがほとんどであった．しかし，根や枝の存在は流木のスタックや着床にも大きな影響を及ぼすことが指摘され，このような現象を検討するには複雑な形状の流木を直接再現でできる数
1）
値解析モデルが必要となってくる．図3.1は複雑形状の流木の一例である．Kang and Kimura（2018）
お2）
よび，Kang（2019）
では，流木を球体の直列接続で表し，根の影響を縁端の球体の直径を大きくすることで表現した．このモデルは，流木の河床への着床や，再浮上，滑動を伴う浅水流中の流木挙動等を適切に再現できることを，実験結果との比較により確認している（図3.2参照）．しかしながら，枝が複雑に張り出した場合など，形状がより複雑化した場合への拡張は困難である．また，洪水氾濫による家屋（図3.3参照）や車両（図3.4参照）の流出，浮遊という現象も，水害（津波を含む）時においてまま見られ，これらの浮遊物の衝突による他の構造物への破壊などの影響や，人々の避難方法への影響等を検討するにあたり，これらの複雑形状の浮体の河川流中の動態を予測する数値解析モデルの構築は，学術的にも社会的にも極めて重要な課題と考えられる．北海道の河川に目を向けると，冬季には河川流中を氷片が浮遊流下し，条件によっては橋梁の上流側に堆積し，河川閉塞などの悪影響を生じる場合がある．このような現象をアイスジャムと呼び，図3.5はその一例である．アイスジャムは河川防災上重要な現象の一つであり，これも複雑形状の浮体群の河川流中の流送の一形態とみなすことができる．

図3.1

複雑形状の流木の例（筑後川）（rubese.netguru.comより引用）

本章では，このような任意の複雑形状の浮遊体の河川流中の三次元的挙動を予測するモデルを新たに構築したのでこれについて述べる．ベースとなるのはiRICの三次元ソルバーNaysCUBEの流木モデルであり，従来円柱形状のみを対象としていたものを，任意形状に拡張する．構築したモデルはいくつかのモデル計算に用い，その妥当性を検証する．また，本モデルに適用するよう，GUIの改良を行うとともに，
日本語，英語マニュアルの改訂を実施する．

― 62 ―

図3.2

球体群型流木モデルによる根の表現（Kang and Kimura（2018）, Kang（2019））

図3.3

河川氾濫時の家屋の流出の事例（多摩川堤防決壊1974.9，jiji.comyより）

図3.4

河川氾濫時の車両の流出の事例（韓国，釜山，floodlist.comより）― 63 ―

図3.5

アイスジャムの例（常呂川，北海道開発局）（7）数値解析モデル
a）液相のモデル
本研究では河川流を三次元モデルで計算する．用いる三次元モデルは河川解析共通プラットフォーム
iRICの三次元ソルバーNaysCUBE3-7）における流れの計算モジュールある．このモデルは，移動境界適合
格子系上でレイノルズ平均されたナビエ・ストークス方程式と連続式を解くもので，ベクトルおよびテ
ンソルについては反変成分を用い，スタガード格子上で変数配置を行う．移流項の離散化には空間三次精度のTVD MUSCL法が，時間積分には二次アダムス・バシュフォース法が，乱流モデルには二次の非線形k-eモデルが採用された．このため，川の曲がり部における第一種二次流のみならず，乱流の非等方性に起因する第二種二次流などの三次元的流れ特性を再現することが可能である．三次元的基礎式やその離散化の方法，境界条件等についてはweb上の文献に詳述されているのでそちらを参照されたい6）．
なお，本研究では簡単のため剛体が流れに及ぼす影響を無視したOne-wayモデルを採用する．
b）剛体移動追跡モデル
・剛体移動追跡計算の概要
本研究では，一個の剛体を球体の集合で表現する．この手法は，越塚8）がMPS型粒子法において示した計算方法を二次元空間における剛体計算に応用した清水ら9）による研究や，これを応用した初田，赤堀ら10-12），Kimura and Kitazonoら13），長田ら14）の方法と類似であり，これを三次元空間に拡張したものといえる．剛体を球体群で表現する利点の一つは剛体の回転や衝突のモデル化が容易なことである．さらに複雑なモデルとしては，牛島ら15）や，Ushijima et al16）による四面体モデルによるものがあるが，モデルが複雑で計算機負荷も大きく，実河川への適用には現時点では不向きと考えられる．
剛体移動追跡計算のプロセスは，一時間ステップΔt後の各球体の移動位置を，まずは拘束を無視した
独立した球体運動として計算した後，球体中心座標の平均値を取ることで剛体重心を求める．次に，剛体平均角速度，平均移動速度を求め，剛体を再度初期の形状と等しく並べ替える．このプロセスを計算
ステップ毎に繰り返すことで剛体の挙動を追跡する．本手法は計算時間間隔が十分小さい場合，剛体の運動をおおむね近似していると考えられる．具体的な計算手順を次に示す．

― 64 ―

・剛体初期姿勢の設定
図3.6（a）のように，剛体重心を原点におき，剛体の主軸がx軸と一致する姿勢を剛体の基本姿勢と呼ぶことにする．剛体の投入初期の姿勢は，図3.6（b）に示すように，yaw角（図3.6（c）参照）のみを乱数で与える．

図3.6

剛体の基本姿勢と初期姿勢（棒状の物体を例として）・流体構成球体の一時間ステップの移動計算
粒子中の球体移流を表すラグランジュ方程式は17-21），

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

1

𝐅𝐅 
𝐴𝐴� 𝑑𝑑 � 𝜌𝜌�𝜎𝜎𝜎𝜎 � ������ 𝐶𝐶� � 𝐩𝐩

（3.1）

である．ここに， 𝜎𝜎：球体の密度， 𝜌𝜌 ：水の密度，d：球体の直径，CD : 抗力係数，u：水流の流速ベク

トル，u p：粒子移動速度ベクトル，g：重力加速度ベクトル，F p：粒子間衝突力ベクトル，t：時間，CM：付加質量係数（CM = 0.5），A2, A3：二次元および三次元形状係数（A2 = π/4, A3 = π/6），λA_sub：球体の流水方向投影面積の水面以下の部分の割合，λV_sub：球体体積の水面以下の部分の割合 をそれぞれ表
す．木村ら2）の既往の二次元計算では，類似のラグランジュ方程式を一般曲線座標系上で解いているが，三次元計算では煩雑となるため本研究ではデカルト座標系上で解いている．なお流体速度の微分du/dt
は本来粒子軌跡に沿う微分で評価すべきであるが，簡単化のため位置を固定した時間に関する偏微分（ ∂𝑢𝑢/ ∂𝑡𝑡）で評価した．この近似が精度に及ぼす影響の検討は今後の課題である．

― 65 ―

図3.7

一時間ステップ間の剛体構成粒子の移動の様子（棒状の物体を例として）・流体構成球体の一時間ステップの移動計算
図3.7は1時間ステップ（Δt）間の剛体構成粒子の移動の様子を示したものである．移動前の重心位置ベクトルxgは，図3.7を参考に次のように表される7）．

𝐱𝐱 𝐠𝐠𝐭𝐭 ��1
𝐭𝐭� 𝐱𝐱 𝐦𝐦
𝑁𝑁

（3.2）

���

ここに，N：一個の剛体を構成する球体数である．Δt時間後の流木重心の移動速度ベクトルは，式（3..1）で計算された拘束を考慮しない個別球体の移動速度ベクトルu*mの平均をとり，次式で表される．

𝐮𝐮����
�𝐠𝐠

�

1
� 𝐮𝐮∗,����
𝐦𝐦
𝑁𝑁

（3.3）

���

ここに，ug：剛体重心移動速度ベクトルであり，アスタリスクは拘束を考慮しない仮の移動速度を表す．
一方，剛体を構成する球体群全体の角運動量ベクトルは次のように表される．
�� ∗,𝒕𝒕��𝐭𝐭𝐭𝐭
� � �𝐴𝐴� � � ��𝐮𝐮𝐦𝐦� �� 𝐦𝐦� � 𝐠𝐠𝒕𝒕 �

（3.4）

���

ここに，L：一方の剛体を構成する球体群の角運動量ベクトルである．この回転系の角速度ベクトルをωとすると，角運動量ベクトルとの関係は次のようになる9）．

� � � �� 𝐋𝐋（3.5）

𝐼𝐼�� � � �� �𝑟𝑟�� 𝛿𝛿�� � 𝑟𝑟�� 𝑟𝑟�� �（3.6）

ここに，I：慣性モーメントテンソルからなるマトリックスであり，3×3次元の慣性モーメントテンソ
ルの各成分が，
����と表されることから，慣性モーメントマトリックスは次式で表される．
��𝑟𝑟��
� 𝑟𝑟��
� �𝑟𝑟 𝑟𝑟� � 𝐴𝐴� � � ��� �� ��
����𝑟𝑟�� 𝑟𝑟����𝑟𝑟�� 𝑟𝑟����𝑟𝑟��
� 𝑟𝑟��
�𝑟𝑟�� 𝑟𝑟��― 66 ―

�𝑟𝑟�� 𝑟𝑟���𝑟𝑟�� 𝑟𝑟�� ���𝑟𝑟��� 𝑟𝑟��

（3.7）

ここに，rmx, rmy, rmzは，m番目の球体の重心に対する位置ベクトルのx, y, z方向成分を表す（図3.6参照）．
次に，この回転を表すクォータニオンを次のように計算する．

𝑞𝑞� � �𝑞𝑞� , 𝑞𝑞� , 𝑞𝑞� , 𝑞𝑞� � � �cos

ここに，

𝜙𝜙
𝜙𝜙
𝜙𝜙
𝜙𝜙
𝜙𝜙
𝜙𝜙
� sin �, 𝜔𝜔�� sin , 𝜔𝜔�� 𝑢𝑢� sin , 𝜔𝜔�� 𝑢𝑢� sin � � �cos , 𝛚𝛚
2
2
2
2
2
2
��
𝛚𝛚

𝜔𝜔
‖𝛚𝛚‖ ,� � ��‖𝛚𝛚‖

（3.8）

（3.9）

である．また，||ω|| は角速度のノルムを表し，次式で定義される．

‖𝛚𝛚‖ � �𝜔𝜔�� � 𝜔𝜔�� � 𝜔𝜔��

（3.10）

クォータニオンとは，実部一個と虚部三個の計4成分からなり，複素数を三次元的に拡張したものとみなすことができる．時刻t+Δtにおける剛体球体系を剛体回転とみなし，次のように計算し直す．

𝑄𝑄���� � 𝑞𝑞�𝑄𝑄� 𝑞𝑞� ∗（3.11）

ここに，アスタリスク”＊”
は共役クォータニオンを表す．したがって， 𝑞𝑞� ∗は，𝑞𝑞� ∗ � �cos

𝜙𝜙
𝜙𝜙
� sin �, ��
2
2

（3.12）となる．Qは各剛体構成球体の中心位置座標のクォータニオンであり，次式で表される．
��
𝑄𝑄� � �0, 𝐫𝐫�� �0, 𝑟𝑟�� , 𝑟𝑟�� , 𝑟𝑟�� )（3.13）

なお，一般に，二つのクォータニオン 𝑎𝑎�, 𝑏𝑏�の積は，次のように計算される．

𝑎𝑎�𝑏𝑏� � �𝑎𝑎� , 𝑎𝑎� , 𝑎𝑎� , 𝑎𝑎� ��𝑏𝑏� , 𝑏𝑏� , 𝑏𝑏� , 𝑏𝑏� � � �c� , c� , c� , c� �𝑐𝑐� � �� 𝑏𝑏� � �� 𝑏𝑏� � �� 𝑏𝑏� � �� 𝑏𝑏�

𝑐𝑐� � �� 𝑏𝑏� � �� 𝑏𝑏� � �� 𝑏𝑏� � �� 𝑏𝑏�
𝑐𝑐� � �� 𝑏𝑏� � �� 𝑏𝑏� � �� 𝑏𝑏� � �� 𝑏𝑏�
𝑐𝑐� � �� 𝑏𝑏� � �� 𝑏𝑏� � �� 𝑏𝑏� � �� 𝑏𝑏�

（3.14）

式（3.11）の計算には，式（3.14）の関係を二回用いれば可能である．しかし，式（3.11）の関係は，
次のようなクォータニオンの成分から回転マトリックスRを構成し，これを用いることで，より簡単に行うことが可能である．

1 � 2�𝑞𝑞�� � 𝑞𝑞�� �� � �2�𝑞𝑞� 𝑞𝑞� � 𝑞𝑞� 𝑞𝑞� �
2�𝑞𝑞� 𝑞𝑞� � 𝑞𝑞� 𝑞𝑞� �

2�𝑞𝑞� 𝑞𝑞� � 𝑞𝑞� 𝑞𝑞� �
1 � 2�𝑞𝑞�� � 𝑞𝑞�� �2�𝑞𝑞� 𝑞𝑞� � 𝑞𝑞� 𝑞𝑞� �

𝐫𝐫𝐦𝐦 �𝑡𝑡 � �𝑡𝑡� � �𝐫𝐫𝐦𝐦 �𝑡𝑡�2�𝑞𝑞� 𝑞𝑞� � 𝑞𝑞� 𝑞𝑞� �

2�𝑞𝑞� 𝑞𝑞� � 𝑞𝑞� 𝑞𝑞� ��（3.15）

1 � 2�𝑞𝑞�� � 𝑞𝑞�� �（3.16）

式（3.1）～（3.16）を用いた計算を時間ステップ毎に繰り返すことで，剛体の回転運動を追跡可能であるが，この方法では時間が進行するにつれて計算誤差が蓄積され，剛体の形状が乱れる可能性がある．
これを避けるため，本研究ではアルゴリズムに次のような工夫を施した．

― 67 ―

式（3.11）を初期位置クォータニオン（t = 0）からΔtごとに繰り返し実行し，時刻tまで繰り返したとする（時刻tはΔtよりも十分大きいとする．）．時刻tにおける位置クォータニオンQtは，初期位置（t=0）
から式（3.11）を繰り返し用いることで，次のように表される．

𝑄𝑄� � 𝑞𝑞�� 𝑄𝑄���� 𝑞𝑞� �� = 𝑞𝑞����� 𝑞𝑞������ ⋯ 𝑞𝑞��� 𝑄𝑄��� 𝑞𝑞� ���� ⋯ 𝑞𝑞� ������ 𝑞𝑞� ����� 𝑞𝑞���

（3.17）

��� ∗𝑄𝑄��� � 𝑞𝑞
�𝑄𝑄
𝑞𝑞
��
�

（3.18）

さらに，基本姿勢から初期姿勢への変換も次のようにクォータニオンを用いて表すことができる．

ここに，Qφ=0は，基本姿勢を表す位置クォータニオンであり，roll, yaw, pitchの角度を全て0とした式（3.13）

�� は基本姿勢から初期姿勢（図3.6（b））への回転クォータニオンである．これから計算される．また， 𝑞𝑞
は，roll角，yaw角，pitch角をそれぞれ，α，β，γとすると，次のようにあらわされる

𝛾𝛾
𝛾𝛾
𝛽𝛽
𝛽𝛽
𝛼𝛼
𝛼𝛼
𝑞𝑞
�� � �cos , 0,0, sin � �cos , 0, sin , 0� �cos , sin , 0,0�
2
2
2
2
2
2

𝛾𝛾
𝛽𝛽
𝛾𝛾
𝛽𝛽
𝛾𝛾
𝛽𝛽
𝛾𝛾
𝛽𝛽
𝛼𝛼
𝛼𝛼
� �cos cos , � sin sin , cos sin , sin cos � �cos , sin , 0,0�
2
2
2
2
2
2
2
2
2
2

𝛾𝛾
𝛽𝛽
𝛼𝛼
𝛾𝛾
𝛽𝛽
𝛼𝛼
𝛾𝛾
𝛽𝛽
𝛼𝛼
𝛾𝛾
𝛽𝛽
𝛼𝛼
� �cos cos cos � sin sin sin , cos cos sin � sin sin cos ,
2
2
2
2
2
2
2
2
2
2
2
2

𝛾𝛾
𝛽𝛽
𝛼𝛼
𝛾𝛾
𝛽𝛽
𝛼𝛼
𝛾𝛾
𝛽𝛽
𝛼𝛼
𝛾𝛾
𝛽𝛽
𝛼𝛼
cos sin cos � sin cos sin , sin cos cos � cos sin sin �
2
2
2
2
2
2
2
2
2
2
2
2

（3.19）

この関係は，本研究では初期姿勢としてα=β=0を仮定しているので，次のように簡略化される．

𝛾𝛾
𝛾𝛾
𝑞𝑞
�� � �cos , 0, 0, sin �
22（3.20）,（3.18）,（3.20）により，任意の
なお，初期yaw角γについては，前述のように乱数で与える．式（3.17）
時刻tにおける剛体の位置クォータニオンは，基本姿勢からの合成変換により，次のように表すことができる．

𝑄𝑄� � 𝐺𝐺� � 𝑄𝑄��� 𝐺𝐺� ��（3.21）

𝐺𝐺� ´� � 𝑞𝑞� � 𝑞𝑞� ���� 𝑞𝑞������ ⋯ ̂𝑞𝑞����� 𝑞𝑞
��（3.22）

ここに， 𝐺𝐺� �は基本姿勢からの合成クォータニオンであり，次式で表される．

この合成クォータニオンは回転クォータニオンの積であるから，これも回転を表すべきであるが，計算ステップの繰り返しにより誤差が蓄積し，回転以外の運動が混入する可能性がある．そうなると剛体形状が乱れる結果となる．クォータニオンが純粋に回転を表す条件は，単にクォータニオンのノルムが1となることであるから，このような誤差の混入を排除するため，一定時間ステップ毎（例えば100ステ
ップ毎）に，次の規準化の操作を行うことでこの問題を避けることができる．

𝐺𝐺� � �𝐺𝐺� ′�
�𝐺𝐺� ′� �― 68 ―（3.23）

ここに，𝐺𝐺� � � �𝑔𝑔� , 𝑔𝑔� , 𝑔𝑔� , 𝑔𝑔� �（3.24）

�𝐺𝐺� ′� � � �𝐺𝐺� ′� 𝐺𝐺� ′�� � �𝑔𝑔�� � 𝑔𝑔�� � 𝑔𝑔�� � 𝑔𝑔��（3.25）の関係を用いている．なお，式（3.22）の計算において，過去の全ての時間ステップにおける回転クォ
ータニオンを記憶しておく必要は無く，一つ前の時間ステップ迄の合成クォータニオンのみ記憶しておき，それに新しい時刻ステップにおける回転クォータニオンを順次乗じていけば良い．従って，この部
分で計算に必要な記憶容量を増加させることは無い．
これらを用いて，新しい時間ステップにおける各球体の中心位置ベクトルを次のように計算する．
���

����
𝐭𝐭
𝐱𝐱 𝐦𝐦
� 𝐱𝐱 𝐦𝐦� �� � ������ ��𝐦𝐦𝒈𝒈（3.26）

ここに，新しい時刻の重心移動速度ugt+Δtは，式（3.3）のように整列前の球体速度の平均値で評価する．
また，Gは合成回転クォータニオンの行列形で，次式で表される．

1 � 2�𝑔𝑔�� � 𝑔𝑔�� �G � �2�𝑔𝑔� 𝑔𝑔� � 𝑔𝑔� 𝑔𝑔� �

2�𝑔𝑔� 𝑔𝑔� � 𝑔𝑔� 𝑔𝑔� �
1 � 2�𝑔𝑔�� � 𝑔𝑔�� �

2�𝑔𝑔� 𝑔𝑔� � 𝑔𝑔� 𝑔𝑔� � 2�𝑔𝑔� 𝑔𝑔� � 𝑔𝑔� 𝑔𝑔� �

2�𝑔𝑔� 𝑔𝑔� � 𝑔𝑔� 𝑔𝑔� �

2�𝑔𝑔� 𝑔𝑔� � 𝑔𝑔� 𝑔𝑔� ��（3.27）

1 � 2�𝑔𝑔�� � 𝑔𝑔�� �∵ 𝐺𝐺� � � �𝑔𝑔� , 𝑔𝑔� , 𝑔𝑔� , 𝑔𝑔� �

最後に，各球体中心の移動速度を次のように更新する．

𝐮𝐮����
�𝐦𝐦

𝐭𝐭��𝐭𝐭𝐭𝐭
� 𝐱𝐱 𝐦𝐦
𝐱𝐱 𝐦𝐦
Δ𝑡𝑡（3.27）

以上の計算方法により，剛体の連続的な移動を追跡でき，長時間の計算においても剛体形状が乱れることは無い．また，根や枝を有する複雑な剛体形状に対しても，球体初期配置を工夫するだけでそのまま適用できる．
・剛体同士の衝突
剛体同士の衝突力は，三次元球体の個別要素法により容易に評価できる．流木を球体の接続で表現したモデルのメリットの一つはここにある．木村ら2）の二次元計算では，境界条件の複雑さを避けるため，
個別要素法を一般座標系に変換した基礎式を用いているが，三次元の場合一般曲線座標への変換が煩雑になるため，デカルト座標系の基礎式をそのまま用いる．個別要素法の基礎式の詳細は後藤10）によった．

― 69 ―

Step 1（a）

Step 2（b）
図3.8

剛体と構造物の衝突判定の手順

・剛体と物体との衝突
河岸，河床，橋脚等，剛体以外の固体に剛体が衝突した場合の処理について述べる．本研究では境界適合格子を用いているため計算格子は物体の表面に沿って存在するものと仮定する．用いたシグマ座標系では，格子のη-ζ面，ζ-ξ面は常に鉛直（法線が水平）である．一方，ξ-η面は一般に水平では無いが，水平方向に比べて鉛直方向のスケールが小さく，河床面が連続的に変化する通常の河川の場合，ξ-η面と水平面の成す角度は比較的小さいと考えられる．このことを生かして，次のような簡易的な二段階の用法で衝突を判定する．

― 70 ―

Step 1: 物体存在格子柱（GCM）側壁への衝突判定
ξ-η面からみた鉛直方向の格子柱内に物体が存在する場合（越流型物体の場合も含む），この格子柱を, Step 1）．まずGCO内
物体存在格子柱（Grid Column with Object: GCO）と呼ぶことにする（図3.8（a）に流木構成球体が侵入したかどうかを判定する．これには凸型ポリゴン内外判定アルゴリズムを用いた22）．球体がGCO内に存在すると判断された場合のみ，Step 2に進む．
Step 2：GCO内の上下方向の内外判定簡単のため，物体の存在するセルの上下面（ξ-η面）を水平と近似する（図3.8（b）, Step 2）．これによりGCO中の物体と球体の内外判定を極めて容易に行うことができる．

（a）直方体
図3.9（b）三角柱
物体を球体の集合で近似した事例（8）iRICへの実装
a）任意形状の剛体の記述フォーマット
本モデルでは任意形状の剛体を球体の集合として表現する．図3.9（a）は直方体を，図3.9（b）は三角柱を，それぞれ球体の集合で表した例である．対象とする剛体の形状が複雑になればなるほど，その形状を表現するには多くの球体が必要となり，その分計算機負荷も大きくなる．
球体の集合で剛体を近似する場合，一個の球体を構成する球体の数，球体の位置座標，球体の大きさ（直径），球体の密度を指定する必要がある．今回開発したモデルでは，簡単のため，剛体を構成する
球体の直径および密度は均一とした．これを考慮して，球体の直径と密度については「計算条件」のメニューから直接入力することとし，球体の個数と各球体の位置座標は外部ファイル（形状データファイル）より入力するものとした．入力ファイルの書式はMS-Excelやテキストエディタ―での編集の容易さを考慮してCSV形式（コンマ区切りテキスト形式）とした．一行分が一個の球体の座標を表すものとする．CSVファイルの各行は三つの列からなり，一列目にx座標，二列目にy座標，三列目にz座標を記述する．一個の剛体を構成する粒子数については，形状データファイルの行数と一致する．すなわち，個数を陽に指定する必要はなく，必要球体数に応じた行数を有する鵜形状データファイルを作成すればよいことになる．
図3.10は形状データファイルの一例である．この例では，x-y平面上に位置する1m×1mの正方形の板状の剛体を4個×4個＝16個の球体で表したものである．球体の数を増やせば（例えば10個×10個=100個）― 71 ―

より形状を精度よく表すことができるが，その分計算機負荷が大きくなる．最初は球体数の少ない（すなわち粗い）近似で計算をスタートし，必要に応じて球体数を増やしていくという方法をとるのが効率的であろう．

図3.10 形状データファイルの一例（1m×1mのx-y平面上に存在する板状の剛体を4個×4個=16個の球体で表現した例）（b）英語版（a）日本語版
図3.11

NaysCUBEの「計算条件」のダイアログの構成― 72 ―

図3.12 「流木の基本設定」のダイアログb）任意形状剛体に適応するためのGUIの改良
任意形状の物体を取り扱うにあたり，iRIC上のNaysCUBEのGUIの改良を実施した．本モデルでは任意形状の物体も「流木の一形態」として取り扱うことにする．既往の円柱状の流木のモデルに加え，任意形状を取り扱うにあたり，「計算条件」のメニューのうち，流木に関する設定項目を，・流木の基本設定
・流木の応用的設定
・流木の詳細設定の三つに分割した．図3.11は，実際のiRICの「計算条件」のウインドウの様子である．

図3.13 「任意形状」の設定ダイアログ― 73 ―

図3.14

形状データファイル（図3.9参照）をインポートした状態

・「流木の基本設定」に関する入力項目「流木の基本設定」のダイアログの様子を図3.12に示す．任意形状の剛体の流送を計算するには，まず「流木を計算する？」の項目の「流木無し」，「球体群モデル」，「任意形状」の中から，「任意形状」
を選択する（図3.12①）．次に，任意形状の剛体を構成する要素である球体の直径と密度を指定する．
直径の指定にあたっては，あらかじめ任意形状物体を何個の球体で表すかを想定し，すきまがなるべく生じないような直径を指定する必要がある．すきまが生じてしまうと，流体に対する抗力や（Two-way
計算の場合）が過小評価されたり，剛体同士の衝突において，剛体同士がめりこんだりるす可能性があるので注意する．一方，密度については，剛体全体の質量が球体全部の質量と一致するように指定する．
次に，形状データファイルを読み込むため，「任意形状」の項目の「Edit」を左クリックする．すると，
図3.13のような新たなダイアログが現れる．これに，直接手入力で球体群の座標をx-y-z形式で入力することも可能である．しかし，形状が複雑になったり，球体数が多かったりする場合，手入力は大変な手間がかかり，現実的ではない．そこで，図3.9のように，あらかじめテキストエディターかエクセルなどを用いて作成しておいたcsv形式の形状データファイルをインポートする．これには，図3.13の「インポート」をクリックし，形状データファイルを選択することで可能となる．
図3.14は，図3.10の形状データファイルをインポートした状態である．左側にインポートしたデータの数値が，右側にはそのグラフが表示される．グラフについては，二次元的なグラフ表示となっており，三次元的にデータを見ることができないので，このグラフはあまり参考にならない．これについては，
ダイアログのグラフの三次元化がなされることが望ましく，今後のGUIの改良に期待したいところである．

― 74 ―

図3.15 「流木の応用的設定」の入力ダイアログ

図3.16 「流木の詳細設定」の入力ダイアログ
・「流木の応用的設定」に関する入力項目
図3.15に「流木の応用的設定」の入力ダイアログの様子を示す．この項目は，既往のNaysCUBEの流木設定のダイアログと大きな変更はない．
「流木どうしの衝突」，「流木供給の方法」などを，それぞれ，「任意形状剛体どうしの衝突」，「任意形状剛体の供給の方法」などと読み替えればよい．
・「流木の詳細設定」に関する入力項目
図3.16に「流木の詳細設定」に関する入力ダイアログの様子を示す．この項目は，通常の円柱形状の流木に関する詳細設定と共通であるが，今回の任意形状剛体の取り扱いの開始に合わせて，いくつかの新たな設定項目を追加している．

― 75 ―

□反発率について
任意形状剛体とその他の物体の反発率について，自然障害物（川岸など）と，人工構造物（水制，橋脚など）にわけて，それぞれ別個の反発係数を設定できるようにした．この理由は，川岸などの自然障
害物に捕捉される場合と，橋脚や水制などの人工（コンクリート）構造物に衝突する場合では，後者に方が反発係数が大きいことが予想されるためである．それぞれ，「1」とすると「完全反発」，「0」とすると「反発無し」（すなわち，「吸着」）を表すことになる．

図3.17

サクション力に関する説明□「流木同士の吸引力」の考慮
流木同士が湿潤状態で接触すると，サクション（負圧）の働きにより吸引力が発生することが予想される．図3.17はこの様子を図示したもので，流木同士の接触部付近の吸着水が引張力を受けると負圧が生じることが原因とされる．これについて，新たにモデル化を行った．吸引力を「考慮する」を選択す
ると，次の三つの項目を設定する必要がある．なお，この設定は円柱形状の流木に対してのみ有効であり，任意形状の剛体を扱う場合は設定することはできない．
【吸引力の影響範囲係数】
吸引力は，流木同士の距離がある程度近い場合に発揮される．その範囲を衝突する二つの流木の「流木長」の平均値に対する割合で指定する．デフォルトは0.1である．
【最大吸引力の係数】
最大吸引力を決定する係数を指定する．これには，流木直径の距離だけ個別要素法のばねが縮んだ
場合の力を想定し，これに対する割合で指定する．デフォルトは0.0001である．
【吸引力を考慮する最大水没率】
サクション力は完全に流木が水没するとゼロとなる．この特性を表現するため，吸引力を活性化す
る最大水没率を指定する必要がある．これには，流木構成球体の体積水没率の値について閾値を設け，この閾値を超えると吸引力をゼロとするアルゴリズムを用いる．この場合の閾値がこの設定項目の値であり，デフォルトは0.5である．

― 76 ―

c）任意形状の剛体運動の計算結果の可視化
一つの事例として，図3.10に示した1m×1mの正方形の板状の剛体が，幅5m，長さ20mの直線開水路中に投入された場合を考えよう．水路の途中には，右岸側から一本の水制が河道中央に達するまで設置されているものとする．計算格子は長さ方向は20分割，横断方向は8分割とし，水深方向は4層に分割す
るものとする．初期の任意形状物体の投入セルは，図3.18とし，選択したセルにそれぞれ一個の物体を投入するものとする．

図3.18

初期の物体投入セル（ピンク色のセルに各一個ずつ投入）

その他の条件として，流量は3[m3/s]で一定とし，下流端水深は等流条件，マニングの粗度係数は0.02，
計算時間は30[sec]とした．また，計算時間ステップはクーラン数を0.13として自動設定で与え，計算結
果は0.5[sec]ごとにファイル出力した．なお，物体同士の衝突の効果は「考慮する」として計算を行った．

図3.19

任意形状物体の表示設定

物体の表示の設定は，「オブジェクトブラウザー」の「粒子」より行う．すなわち，物体の表示は，物体を構成する要素球体の位置に粒子を離散的に描くことにより行うことになる．図3.19はこの様子を示したもので，「粒子」にチェックを入れ，この上で右クリックすることで，「粒子設定」のダイアログが現れる．これにより，粒子の色とサイズを指定する．この際，図3.19に示すように，粒子は球形ではなく，正方形で表示されること，また，「サイズ」は要素球体の実サイズではなく，画面上のサイズ（ピクセル数）を指定するものであることに注意されたい．したがって，「サイズ」は物体の形状がなるべく実体に近く表示されるように指定すればよく，これには試行錯誤が必要となる．
「色」の設定についても同様である．
なお，現在のiRIC-GUIでは，粒子を陰影を有する立体として表示することはできないため，物体同士が重なるような場合には，物体の形状や位置が視覚的にとらえにくくなる場合がある．この点の改善については，今後のGUIの改良に期待したい．

― 77 ―

図3.20 計算結果における物体移動の様子（5[sec]間隔）（カラーコンターは水深，多数の曲線は流線をそれぞれ表す）

図3.20は，本シミュレーション結果の上流型右岸よりからの鳥瞰図を5秒間隔で示したものである．
図のカラーコンターは水深の分布を，曲線群は流線をそれぞれ表す．水路上方から投入された正方形の板状の剛体は，互いに衝突したり，側岸や水制に衝突したりしながら流下していき，その一部は水制に衝突して右岸側に集積し，一部は水制の左岸側をすり抜けて流下していく様子が表されている．この図でみるかぎり，物体の運動は極めて自然で，実際の流水中の物体の運動を少なくとも定性的には良好に
シミュレートしているものと考えられる．
図3.21は，先ほどの計算条件に対して物体の数のみを増加させ，80個とした場合のt=7.5[sec]における計算結果の様子である．水制上流側で多くの物体が重なり合い，物体の形状や位置がわかりにくい表示となっている．この点の改良には，個々の物体の構成球体待の色を別個に設定できるようにするか，
球体に陰影をつけて表示するなどの，GUIレベルでの改良が必要となる．これについては，今後のiRIC-GUIの対応に期待したい．

― 78 ―

図3.21

図3.22

物体の数を80個に増加させた場合の計算結果（t=7.5[sec]）

流出した家屋が橋脚に衝突する場合を想定したシミュレーション例（一家屋の要素球体数は125個）

物体の数が増えたり，物体を構成する要素球体の数が増加したりすると，計算時間が著しく増加する．
これは主に物体同士の衝突判定，すなわち近傍球体の検索に時間を要するためであり，この点は現ソルバーの衝突判定アルゴリズムが脆弱なためであることは否めない．この点の改良は今後の課題のうちでも極めて緊急性の高いものと認識しており，すみやかにこの点が解消されるよう，改良を進めていく予定である．
図3.22は，蛇行する河川において．氾濫により流出した家屋が橋脚に衝突する場合を想定したシミュ
レーションの例である．このケースでは，一つの家屋を125個の要素球体で表現している．この結果も― 79 ―

家屋の流下と橋脚への衝突，一部の家屋の橋脚間のすり抜けなどの様子が定性的に妥当に再現されていると考えられる．

図3.23

任意形状物体移動計算に対応したNaysCUBE ver.3.1用マニュアル（第7版，全91頁）の表紙と
内容の一部（c）任意形状物体移動計算に対応したNaysCUBEマニュアル
本研究の成果である，任意形状の物体の移動計算に関する機能を備えた新ソルバー（NaysCUBE
ver.3.1）のリリースにあたり，マニュアルの改訂を行った．図3.23は，新しいマニュアル（NaysCUBE
ソルバーマニュアル第7版，ver.3.1用）の表紙とその内容の一部である．また，これに合わせて，マニュアルの内容を全体的に見直し，不十分な記述等を修正した．新マニュアルは，一つ前のバージョンのマニュアルよりページ数が22ページ増え，総ページ数は91ページとなり，大幅に充実した内容となっている．
（9）本章のまとめ
本章は，iRICの備える三次元ソルバー，NaysCUBEにおいて，新たに任意形状の物体の移動計算を行うことができる機能を追加し，その妥当性を示すとともに，iRICへの実装等を行った研究について述べたものである．
任意形状物体移動のモデル化にあたっては，三次元的な剛体の移動・回転を表すモデルを新たに構築した．これには，クォータニオン（四元数）という数学的概念を用いることにより，極めて明解に，かつシステマティックに実施することができた．さらに，任意形状物体同士の衝突には三次元空間の個別要素法を用いることで，妥当に評価できることを示した．また，任意形状物体と，橋脚，水制，護岸などに対する衝突についても，シグマ座標系の特性を生かした簡易なモデルを新たに構築し，実装した．
この衝突モデルは，物体と障害物の衝突を妥当に再現できることを，いくつかのテスト計算を通じて示した．

― 80 ―

また，任意形状物体の移動計算をiRICに実装し，GUIにより行うにあたっては，形状データファイルのフォーマットを設定するとともに，計算条件のダイアログのうちの流木・任意形状物体に関する項目を3段階に分け（基本設定，応用的設定，詳細設定），ビギナーから上級者までのユーザーが，レベルに応じて使いやすい構成とした．
さらに，今回のソルバー改訂に伴い，ソルバーマニュアルの改訂を行った．新たなソルバーは，
NaysCUBE ver.3.1と位置づけ，これに対応すべく，大幅な加筆と，全体的な内容の見直しと修正を実施
した．こよりより，マニュアルのページ数は前バージョンより22ページ増加し，全91ページの充実したものとなった．
今後の課題としては，iRIC本体のGUIの三次元物体への適用や，計算速度の向上，並列化などが挙げられる．これらについて，今後も引き続き取り組んでいきたい．
謝辞
本研究の遂行に当たり，京都大学大学院，後藤仁志教授および原田英治准教授より有益なご助言をいただいた．ここに記して深甚なる謝意を表する．
参考文献
1） Kang, T. and Kimura, I.: Computational modeling for large wood dynamics with root wad and anisotropic
bed friction in shallow flows, Advances in Water Resources, Vol.121: pp.419–431, 2018.
2） Kang, T.: PhD. Dissertation, STUDIES ON MORPHODYNAMICS IN SHALLOW RIVERS WITH
EFFECTS OF VEGETATION AND LARGE WOOD USING COMPUTATIONAL MODELS, Hokkaido
University, 2019.
3） Homepage iRIC ”Changing River Science”, http://i-ric.org/ja/ 2019.
4） 木村一郎, 北園和也: 流木捕捉形態の流木リチャードソン数依存性と3D-2D 型モデルの適用性,土木学会論文集, vol. 73, No. 2, pp.553-562, 2017.
5） 北園和也，木村一郎，清水康行：三次元流動モデルと流木モデルのカップリングによる構造物周辺
流木挙動解析，土木学会論文集（A2, 応用力学），2017.
6） Kitazono, K, Kimura, I, Shimizu, Y and Kyuka, T.: Computation on driftwood motions around obstacles
coupling with a three-dimensional flow model, Proc. 19th Congress of the Asia Pacific Division of theInternational Association for Hydro Environment Engineering & Research, Colombo, Sri lanka, pp.1-8, 2016.
7） Nelson, JM, Shimizu, Y, Abe, T Asahi, K, Gamou, M, Inoue, T, Iwasaki, T, Kakinuma, T, Kawamura, S,
Kimura, I, Kyuka, T, McDonald, RR, Nabi, M, Nakatsugawa, M, Simões, FR, Takebayashi, H and Watanabe,
Y: The International River Interface Cooperative: Public Domain Flow and Morphodynamics Software for
Education and Applications, Advances in Water Resources, Vol.93: pp.62-74, 2016.
8） 越塚誠一：数値流体力学，培風館，pp178-181，1997．
9） 清水義彦, 長田健吾：流木形状を考慮した個別要素法による橋脚周辺の流木集積過程に関する数値
実験，水工学論文集，第51巻，pp.829-834，2007
10）初田直彦, 赤堀良介, 清水康行：蛇行流路の流体場と流木の挙動に関する実験と数値解析，応用力学
論文集，Vol.15，pp.415-422，2012.
11） Hatta,N, Akahori, R and Shimizu,Y.: Experimental study and numerical analysis on field on meanderingchannel and movement of woody debris, Journal of Applied Mechanics, JSCE, Vol. 15: pp.415-422, 2012.
12）赤堀良介：橋脚周辺における流木の3次元的集積に関して, 土木学会論文集B1（水工学）Vol.74,
No.4, I_679-I_684, 2018.

― 81 ―

13）Ichiro Kimura and Kazuya Kitazono: Effects of the driftwood Richardson number and applicability of a 3D–
2D model to heavy wood jamming around obstacles, Environmental Fluid Mechanics（online）, https://doi.
org/10.1007/s10652-019-09709-6, Springer Nature B.V., pp.1-23, 2019.7.8.
14）長田健吾, 清水義彦：流木群の曲げ変形を伴う堆積過程を説明する数値解析法の開発, 土木学会論文
集B1（水工学）Vol.74, No.4, I_763-I_768, 2018.
15）牛島省, 福谷彰, 牧野統師：3次元自由水面流中の接触を伴う任意形状物体運動に対する数値解法，
土木学会論文集，Vol.64, No.2, pp.128-138, 2008.5.
16）Ushijima, S., Makino, O. and Yoshikawa, N.: 3D Numerical prediction for transportation and entrapment ofdriftwood with T-type solid model, Journal of Hydroscience and Hydraulic Engineering, JSCE, vol.27（No.1）: pp.11-21, 2009.
17）Gotoh, H, Okayasu, A & Watanabe, Y: Computational Wave Dynamics, Advanced Series on Ocean
Engineering, 87, World Scientific, 2013.
18）後藤仁志.，数値流砂水理学，森北出版，2004.
19）羽田野袈裟義, 天野卓三, 松本治彦, 埜口英昭, 平野宗夫：水面上の油層の1次元的広がりについて,土木学会論文集, No.593/II-43, 117-124,1998.
20）Persi E, Petaccia G, Fenocchi A, Manenti S, Ghilardi and Sibilla, S.: Hydrodynamic coefficients of yawed
cylinders in open-channel flow, Flow Measurement and Instrumentation, Vol.65: pp.288-296, 2019.
21）Ruiz-Villanueva, V, Bladé, Sánchez-Juny, EM, Marti-Cardona, B, Díez-Herrero, A and Bodoque, J.M. Twodimensional numerical modeling of wood transport, Journal of Hydroinformatics, Vol.16（No.5）: pp.10771096, 2014.
22）Tomoya, Visual Basic Library, http:www5d.biglobe.ne.jp/～tomoya03/.

― 82 ―