３．三次元モデルの機能向上（担当：木村一郎）（１）概要
国内外の既存の河川解析共通プラットフォームの多くが二次元ソルバーがメインのソルバーとして位置づけられていのに対し，iRI
Cの特徴の一つは，三次元モデルを備えている点である．iRI
Cの備える三ys
CUBEとDe
lft3Dの２つがあるが，Delft3
Dは静水圧近似のモデルであること
次元ソルバーとしては，Naから，純然たる三次元モデルは現時点ではNays
CUBEのみである．したがって，Nays
CUBEソルバーを適切に管理発展させていくことは，iRI
C全体の発展を考える上でも極めて重要であると考えられる．
本研究では，既存のNays
CUBEについて，次のような機能追加を実施した．
・流木モデルの追加
・移流項のスキームに，既往の一次風上差分，TVDMusclに加えて，「QUICKスキーム」を追加
・非洗掘セル設定機能の追加
・アンダーパス機能の追加
以下，各項目について述べる．
（２）流木モデルの追加a）概要
流木は，洪水時の橋脚部集積による水位上昇，ダムゲートや取水口などへの閉塞等の問題に関与し，その動的特性の解明が重要な課題となっている．従来から，流木のモデルとしては，二次元モデル1,2,3）に対する適用例が報告されている．このうち，清水らおよび赤堀らのモデルでは，流木を球体を直列に接続した形状でモデル化しており，計算機負荷が比較的小さいのが特徴である．しかしながら，二次流等の流れの三次元性の影響は流木の挙動に対して無視できない場合が容易に想定される．例えば，湾曲部の流れ，第二種二次流が存在する場合などがこれにあたる．これに対して，三次元モデルに対しても，
牛島ら4） によるモデルが報告されている．牛島らのモデルでは，流木を四面体の集合でモデル化しており，清水らのモデルよりも精緻に流木形状を反映できるのが特徴であるが，計算機負荷は相応に大きいことが予想される．
本研究では，iRI
C三次元ソルバー Na
ys
CUBEに流木モデルを組み込み，流れの三次元性と流木挙動の関係に関わる種々の現象に対して適用可能なツールを提供するものである．流木のモデル化においては，計算機負荷と精度とのバランスを考慮して，清水ら1)が二次元モデルに用いたものと同様，球体の直列結合で表すこととする．
b）モデルの基本コンセプト
流木の挙動は厳密には三次元的であるが，鉛直方向の変動は相対的に小さいと仮定し，水面位置における平面二次元的な挙動に限定して記述するものとする．また，流木は流体の運動から影響を受けるのみでなく，流木が流体に対して施す力が存在するが，本研究では後者を無視したOn
eWayモデルとする．
また，流木どうしの衝突の影響についても無視することとした．これらの仮定のもとで，本研究では，流木の挙動を赤堀ら2)にならい次のような４段階でモデル化する．

2――5

図１

本研究の流木モデルのコンセプト（赤堀ら2））

１）流木を球体の直列で表現する．
（図１の左上）２）流体を構成する球体を独立にラグランジュ方程式を用いて時間ステップΔtの間移動させる．
（図２の右上）３）移動後の球体群の重心位置，並行移動速度，回転速度を計算する．
４）前ステップの情報をもとに，移動後の球体位置を直列となるよう再調整する．
（図１の下）c）一般曲線座標における球体群の並行移動，回転移動関係式の誘導と流木移動の計算手順
Nays
CUBEは一般曲線座標系を用いている．このため，流木に関する定式化部分も一般曲線座標で表記しておくこととする．

図２ 粒子群の移動の概念図
まず，流木を構成する個々の球体の時間ステップΔt時間移動を，La
grange方程式により計算する．一般に流体中の固体の運動は，次のLa
grange方程式により計算される2),5)．

3――5）（1

ここに，ρ ：流体密度，σ ：粒子密度，CM：付加質量係数，CD：抗力係数，A2：粒子の二次元形状係数（=π /4），A3：粒子の三次元形状係数（=π /6），d：粒径，u：周囲流体流速，u：砂粒子速度，ν ：p
動粘性係数，t：時間の次元量，g：重力である．運動方程式の左辺第一項は粒子に作用する慣性力，右辺第一項は粒子に作用する流体抗力である．
本研究で対象とする流木の計算においては，計算では水面に浮上した粒子を取り扱うため，A2に関する補正が必要となる．こではA2は，粒子の水面下に存在する部分の投影面積である．右辺第二項は周囲流体の加速に伴う力，右辺第三項は付加質量の加速に必要な力である．A3に関してもA2と同様に水面上の粒子としての補正を行う．ここでA3は粒子の水面下に存在する部分の質量に関するものである．右辺第四項（Basset項）は粒子と周囲流体の非定常な相対運動履歴に関する力，右辺第五項は重力と浮力，
右辺第六・七項（Magnus力，Saffman力）は粒子に作用する揚力である．運動方程式において Basset項，Magnus力，Saffman力は他の項と比較してオーダーが小さく，液相中の個体運動に関して近似的にしばしば無視される5)．本研究でもこれらの項を無視する．第五項に関しては本モデルでは水面での水平方向の移動に関する記述のみ必要となるため，これも無視した．
結局，本研究で用いる流木構成粒子のLagrange
方程式は，次式となる．
（2）
流れ場に対する数値計算で流速を求めた後，それぞれの球形粒子において上記の運動方程式を解く．こ
うして得られた個々の球体の移流流速をもとに，新たな球体の位置座標を次のように計算する．
（3）
上記の球体間の相互作用を考慮しない計算手順により，図２の左側の円内の粒子群（この場合は６個）が，Δｔ時間後に右側の円内に移動したとする．uは，i番目の粒子の速度ベクトル，Tは重心間の移動
iは，重心を基準とした場合の書く粒子の相対位置ベクトルである．ここで，qは次のよ速度ベクトル，qiiうに表記される．
（4）
はi番目の粒子の（絶対）位置ベクトル，rは重心の絶対位置ベクトルであり，これらには次ここに，rigの関係がある．
（5）
ここに，nは粒子数であり，この例の場合は６である．重心間の移動速度ベクトルTは次のように表される．
（6）

4――5一方，粒子群全体の角速度ベクトルRは次のように表される．
（7）
式（４）のIは慣性モーメントであり，一般曲線座標系では次のように表記される．

（8）xh
ここに，rd:球体の密度，x,h:一般曲線座標の平面成分，e,e共変基底ベクトル，qi,qi:相対位置ベxh:クトルqの反変成分，g,gg：接続の係数（二次元）の各成分である．一方，角速度ベクトルについixxhh,xhては，（9）と表され，二次元では鉛直方向成分のみが残る．すなわち，（10）
ここに，Rz: 角速度ベクトルの共変z方向成分（鉛直方向成分）である．これを反変z方向成分に変換
すると，（11）となり，共変成分と同一となる．ここに，（12）であり，平面座標のヤコビアンに相当する．以上のように計算された球体群の重心並行移動速度ベク
トルTと，球体群角速度ベクトルRを用いて，各球体が相互の位置関係を変えないとした場合の流速を次のように計算する．
（13）
上式を成分表示すると，次のようになる．
（14）

5――5（15）次に，移動後の重心位置ベクトルXgと，流木の角度Θ次のように計算する．
（16）
7）（1ここに，kは座標を表す添え字であり，ξあるいはηをとる．最後に各球体の位置を次式で一列に整列するように調整し，時間の１ステップ分（Δt時間分）の計算を終える．
（18）

（19）
以上のサイクルを流れの計算と並行して繰り返すことにより，流木あるいは流木群の運動に関する計算が可能となる．
d）Nays
CUBEへのモデルの組み込み
前述のモデルをNaysCUBEへ組み込むにあたっては，次のような方針のもとで進めた．
・流木の扱いは，前述の球体直列モデルに加え，中立粒子モデル，単一球体モデルも選択できるようにし，ユーザーが流木の形状や質量の影響を考慮できるようにする．
・流木形状は直線とし，曲がりや枝分かれは考慮しないものとする．
・流木の長さ，直径，密度，一つの流木を構成する粒子数はユーザーが設定可能とする．
・初期の流木位置はユーザーが任意にマウス選択で設定できるようにする．
・初期の流木の角度はソルバーが乱数で与える．

図３ Nays
CUBEプリプロセッサの流木設定画面

6――5

図３は，プレイプロセッサにおける設定画面であり，この例では矩形領域の場合を想定している．左「流木」の項目とチェックボックスが追加されている．これをチ側の「オブジェクトブラウザー」中に，ェックして，右側の任意の格子セル上でマウスを右クリックすると，図４のようなポップアップウインドウが現れる．ここで，このセルに初期に設定する流木の本数を指定する．

図４

初期流木配置本数設定ポップアップウインドウ

初期流木をセットするセルは，任意に複数個設定することが可能であり，各流木初期配置セルに配置可能な流木本数は最大で９９本である．
次に，流木一本の属性については，「計算条件」の中に新たに追加された「流木」設定ダイアログを用いて行われる．流木設定ダイアログを図５に示す．

図５ 「計算条件」の「流木」設定ダイアログ
各設定項目について次に説明する．
・「流木を計算する？」
ここでは，流木計算の有無と，計算するならばその流木モデルを選択する．流木モデルは次の３つが用意されている．

7――5

【中立粒子モデル】：流木を体積も質量も有しない中立粒子とした扱うモデルであり，流木は流れに完全に追従して運動する．
【等価球モデル】：流木を体積と質量の等しい一つの球体に置き換えて計算する．並行移動の慣性力は考慮されるが，回転移動の慣性力は考慮されない．
【球体群モデル】：前述のモデルであり，流木を質量と体積を有する複数の球体の直列結合として表現する．並行運動の慣性力だけでなく，回転運動の慣性力も考慮される．
・流木長[m]
流木の長さをメートル単位で入力する．
・流木直径[m]流木の直径をメートル単位で入力する．
・流木密度[kg/m3]の単位で入力する．
流木の密度（物体としての密度）を[kg/m3]
・一本の流木を構成する粒子数一本の流木を構成する粒子の数を整数値で入力する．粒子数を多くすれば，より精度の高い計算が可能となると考えられるが，その分計算機負荷は大きくなる．
・周期的に流木を供給する？
流木を初期に一回のみ供給するか，それとも計算期間を通じて周期的に供給するかを指定する．
・流木計算開始時刻[秒]流木の計算を開始する時刻を[秒]単位で指定する．計算初期から流木を投入する場合は0とする．
・流木供給時間間隔[秒]前述の設定項目「周期的に流木を供給する？」で，「Yes」を選択した場合，ここでその供給間隔を[秒]単位で指定する．
・外部ファイル出力？
外部の可視化ソフトウエアの使用や，統計処理などに供するための流木に関する外部ファイル出力の有無を指定する．外部出力ファイルはテキスト形式のシンプルなもので，流木の重心位置をxyz形式で出力する．

図６

流木計算の一例（水没型水制が存在する場合）

8――5）計算事例e
本流木モデルを用いた計算の一例を図６に示す．このケースでは，直線矩形断面開水路の底面右岸側に２本の水没型水制が存在する場合を対象としている．本計算条件では，流木長が計算格子よりも十分大きいため，流木は棒状の形状で表現されている．また，本モデルは流木同士の衝突は考慮していないため，一部で流木が重なったり，交差したりしている様子がみてとれる．これらの点の改善には，流木同士の衝突を，DEM（個別要素法）などを用いて考慮する必要があるが，この点については今後の課題としたい．
（３）その他の機能向上a）QCUI
Kスキームの追加
これまでの，iRICソルバー Nays
CUBEが備えていた移流項の空間スキームは，一次風上差分と，三次精度TVDMusclスキームの二種類であった．一次風上法は計算の安定性は良いが，数値粘性が大きく，特に計算格子が粗い場合には流速分布が鈍ったり，計算結果が過度に平滑化されたりすることが知られている．このため，特に非定常現象を直接解像するURANS（UnsteadyReynoldsAveragedNavierStokesEquationsModel）では，適切な解が得られにくい．一方，三次精度TVDMusclは，空間三次精度であり，一次風上法に比べて格段に精度が高く，数値粘性も小さい．しかし，ｘ，ｙ，ｚ方向の全ての計算においてTVD条件を課すため，条件によっては解が鈍り，特にURANS計算において剥離渦などの非定常現象の再現性が劣る場合が想定される．これらを勘案して，新たにQUICKスキームを追加した．
QUI
CKスキームは，空間に二次の精度を有するスキームであり，数値粘性を表す二次の打切り誤差が現れないという大きな特徴がある．このため，比較的用いる格子点数の少ない（すなわち，ステンシルの広がりの小さい）コンパクトなスキームであるにもかかわらず，精度の高い計算が可能である．特に，URANS計算には適したスキームであり，剥離渦などを適切に再現できることが既往の研究で指摘されている．ただし，QUICKスキームはオーバーシュート，アンダーシュートが生じることもしられており，流速の不連続な変化等を伴う場合は非合理な結果を与える可能性があることに留意する必要があろう．
b）河床が掘れないセルの設定
床止め工がある場合，河床が強固な岩盤で構成される場所がある場合における河床変動計算では，堆積（河床上昇）は許すが洗掘（河床低下）は許さないような計算セルの設定が必要である．このための機能を新たに追加した．
掘れないセルの設定にあたっては，単にセルの河床低下を禁止するのみならず，流砂の連続性の確保に留意する必要がある．例えば，掘れないセルに堆積が存在しない場合には，その下流側のセルに掃流砂が供給されない点を考慮したモデル化を行う必要があり，これを誤ると流砂の湧き出しや吸い込みが生じることとなる．

9――5

図７

非洗掘セル設定の様子c）アンダーパスの設定
ここでのアンダーパスとは，構造物の底部に存在するトンネル状の空洞を指す．穴あきダムのようにダムの底部に洪水吐き（bottomoutlet）が存在する場合，プール型魚道のプール隔壁に穴が存在する場合等の計算に適用できる．Nays
CUBEのアンダーパス設定機能は，構造物を設定したセルのみに対して有効であり，底面から構造物頂部よりも低い位置にある空洞の高さをセル数で設定するものとする．図８
はプレイプロセッサのアンダーパス設定画面を示している．また図９は，アンダーパスを設定した魚道の計算例であり，上下流端は周期境界条件を設定している．アンダーパス部分に流れが集中し，洗掘が生じている様子がわかる．

図８

プリプロセッサのアンダーパス設定画面

0――6図９

アンダーパスを設定した魚道の計算事例（上下流端周期境界条件）

謝辞：本章の研究においては，愛知工業大学准教授，赤堀良介博士から流木に関するモデリングに対して有益な助言をいただいた．ここに記して深甚なる謝意を表する．
参考文献（第三章分）
1) 清水義彦・長田健吾：流木形状を考慮した個別要素法による橋脚周辺の流木集積過程に関する数値
実験，土木学会水工学論文集，第５１巻，pp.
829834,2007.
2) 初田直彦・赤堀良介・清水康行：蛇行流路の流体場と流木の挙動に関する実験と数値解析，応用力学論文集(CDROM)，Vol.
15,pp.
I.
415I.
422,2012.
8.
3) 中川一・井上和也・池口正晃：流木群の流動に関する研究(3)，京都大学防災研究所年報，第37巻,B2,pp.
459473,1994.
4) 牛島省・福谷彰・牧野統師：３次元自由水面流中の接触を伴う任意形状物体運動に対する数値解法，
土木学会論文集，Vol.
64,No.
2,pp.
128138,2008.
5.
5) 越塚誠一：数値流体力学，培風館，pp
1
78
181，1997．

1――6