マルチレベル格子法に基づく水際移動計算高速化の検討
Co
mputationalmodelbaseonMultilevelgridforWetandDrycells
木村 一郎1IchiroKI
MURA
北海道大学 工学研究院 准教授

1

要旨

津波などの水際移動を伴う数値解析を効率良く計算することは，減災の視点から極めて重要な課題である．詳細な地形データをそのまま生かすには計算格子は必然的に細分化されるが，高解像度の要求される箇所に部分的に細分化格子を用いるマルチレベル格子法は計算精度を落とさず計算効率を向上させるものと期待される．そこで本研究では，斉藤らが開発した浅水流モデル動的マルチレベル格子法であるKMRを水際移動計算に適用可能となるよう改良を実施し，計算精度と計算効率とを検証した．この結果，本手法は既往の手法に比べて同精度を確保しながら50%程度の計算速度向上が見込まれることが確認された．
《キーワード：洪水氾濫計算，動的マルチレベル格子，KMR》

7――3

1.研究の目的
河川の氾濫や洪水などの防災を考えた場合，数値シミュレーションは有効な手段である．昨年，タイのチャオプラヤ川流域では，大規模な洪水が発生したが，現在その再現計算やある対策を講じた際，どのような効果が期待できるかなどを判断するのに数値シミュレーションが用いられている．洪水域における複雑な地形や詳細な構造物を考慮に入れ計算を行ない，精度の高い予測をする必要があるが，その際，精度を要求すればするほど，効率の良い数値計算は困難になる．また，タイの洪水に限らず，近年多発する集中豪雨等によって引き起こされる河川の氾濫や洪水を，実現象よりも早いスピードで予測することが求められており，計算時間を短縮することは重要な課題である．
近年こうした計算の効率化を図るため，現象の空間的，時間的変化に応じて局所的に計算格子を細分化あるいは結合化するKMR（KinematicMeshModel）法や領域分割を動的に行うAMR（AdaptiveMeshModel）法が注目を集めている．河川の氾濫シミュレーションにおいて最も重要なことは，どのように水が進行し，広がっていくかであり，それに対し，既に浸水し，局所的に平衡状態になっているような場所での重要度はそれほど高くない．このような背景のもと，本研究では，KMR法による水際移動を伴う現象の数値計算の妥当性について検討する．まず単純化したモデル地形について計算法の妥当性を検討し，次に実地形に対する適応性について検証する．

2.数値解析手法
（1）方針
本研究で用いるKMR法は，マルチスケールの現象に関する物理量を空間変化や時間変化に応じて局所的に計算格子を細分化あるいは結合化する手法である．本研究では，齊藤ら1）によって開発されたモデルに，さらに水際でのモデルを加えたものを用いる．
（2）数値計算モデルの概要
KMR法は多層構造格子法の一種であり，デカルト座標系の直交格子を複数重ねた格子系を用いる．最大セルからなる層を L=1とし，層の次元が一つあがると，一つ前のセルが4分割されるようなモデルを採用する（図1）．なお本研究では，最高位の層は，L=4としている．各時間ステップ内での計算は，低層のセルから行うが，その際， 計算を進めるにあたり，隣り合うセルの情報が不可欠である．隣接するセルが，同じ層である場合には，問題なく計算を行えるが，そうではない場合には，各定義点における物理量を与える際，なんらかの近似が必要である．小さいセルを結合し，大きいセルとして物理量を与えるプロセスを外挿，その逆を内挿と定義する．内挿は，線形内挿で求め，外挿はセル内定義点の値の単純平均で求めている．これら内，外挿の精度に関してはさらに検討を要するところである．また，セルの分割と結合は，各時間ステップの計算終了後に行っている．
（3）格子分割パラメータ齊藤ら1）は格子分割パラメータを，乱流スケールを表すストレインパラメータと，ローテイションパラメータとしている．しかし，氾濫計算で得たい情報は，氾濫域の拡大や収束がどのように進んでいくかであり，特に精度よく得たい情報は，その水際でのものである．この点を踏まえ，本研究では，格子分割パラメータを水深比とした．基準となる水深は初期設定で任意に決定する．水深は，各時間ステップで更新される物理量のうちの一つであり，計算負荷の低減の観点からも有効なパラメータであるといえる．パラメータの閾値は，格子層位が上がるにつれて大きい値を設定する必要があるが，この閾値の決定方法について，今回は試行錯誤的に行っている．

8――3図1 計算セルの分割（4）基礎式
氾濫計算に用いる基礎式は，次に示すデカルト座標系におけるレイノルズ平均された浅水流方程式である2）．
（1）

（2）

（3）

ここに，h：水深，（U,V）： 水深平均流速ベクトルのx，y方向成分，（M,N）： 流量フラックス（M =hU,N=hV）；：水深平均レイノルズ応力テンソル，（τby）：底面摩擦応力，n：動粘性）；（τbx,係数，θ：水路床勾配，β：運動量係数（乱流では1.
0，層流では1.
2の値を用いた．
）をそれぞれ表す．
底面摩擦応力については次のように評価した．
（4）
ここに，fは摩擦係数であり，マニングの粗度係数 nより次式で求められる．
（5）（5）乱流モデル
乱流モデルは木村ら2）が開発した，修正二次非線形ゼロ方程式モデルを用いた．本モデルの構成則は次のように表される．

9――3（6）（7）

渦粘性係数は次のゼロ方程式モデルで評価する．
（8）（6）計算スキーム
運動方程式の移流項の離散化には保存形式に適応できるスキームとして二次QUICKスキームを用いた．

3.計算条件（1）計算条件設定方針
本研究では，以下の水路の条件のもとで数値実験を行った．実際の地形において，河床勾配が一定であることはあまりないことや，氾濫域には水の侵入を妨げる構造物があることなどを考慮に入れ，以下の三通りの河床の条件で計算を行った．
（i） 河床勾配が一定である場合（ii） 河床に水の浸入を妨げる障害物がある場合（iii）河床緩やかな丘があるような場合
なお，詳細な条件に関しては，図2に示す．
計算領域は，幅1.
8m，全長12.
06m，水路の勾配を0.
002，マニングの粗度係数を0.
03とした．また水路への流入量として，図3のようなハイドログラフを与え，境界条件は，上流端で流量フラックスを一様，下流端では自由流出条件とした．初期条件での計算格子は，全領域で格子層位L=1で⊿x=⊿y=18.
0cmを与えた．この状態での計算格子数は， 67（0=67×10）である．

0――4（i）河床勾配一定（平面図）（ii）障害物あり（平面図）（iii）なめらかな丘あり（鳥瞰図）
図2 河床の条件

図3 計算に用いたハイドログラフ

1――4（2）検証方法
今回は，検証材料が不十分なため計算結果と，実現象とを比較することはできない．その検証に関しては，今後の課題とする．本研究では，KMR法を用いた計算結果を格子層位 L=4固定格子での計算結果と比較することでKMR法 の精度に関して検証考察する．

4.計算結果とその考察（i）－（iii）の3つの条件について，KMR法による結果と，比較のために，格子層位 L=4固定格子
前述のの場合の計算結果は合わせて示した．なお，水深を比較材料としている理由は，氾濫計算で最も重要な情報が水深であると考えられるためである．
（1）河床勾配一定水路（ケースi）
KMR法を用いた（a）
では，水際ほど微細に格子分割されており，その結果は概ね格子層位 L=4固定格子（b）より得られる結果に等しい．しかし，水際を詳しく見ると，L=4固定格子を用いたケース（b）に関しては，水路に対して一様な横断方向の直線的な水際を持っているが，KMR法を用いたケース（a）
では，わずかに先端部分が乱れている．一方，水際よりも後方の部分に関しては，ほぼ同一の結果が得られている．

（a）河床勾配一定 KMR法
（t=30s）（b）河床勾配一定 L=4固定格子（t=30s）
図4 計算結果（i）（2）障害物あり（ケースii）
KMR法を用いた（a）では，障害物周辺と水際で細かく格子が分割されており，L=4固定格子（b）と同程度の精度の結果が得られている．水際に関しては，ケースiの場合と同様に，KMR法を用いた場合には若干乱れが生じている．

2――4（a）障害物あり KMR法（t=30s）（b）障害物あり L=4固定格子（t=30s）
図5 計算結果（ii）（3）なめらかな丘あり（ケースiii）
流入開始後15秒の（a）と（b），30秒後の（c）と（d）はそれぞれ近い水深の結果が得られている．15秒後のものは，まさに水が丘を越えようとしている状態であり，その周辺では水深が小さくなっている．その様子は，KMR法を用いた条件（a）での格子分割にも表れており，その結果がL=4固定格子（b）のものに一致していることから，このモデルの精度の高さが確認できる．30秒後の（c）（d，）に関しても同様なことが言える．また，前節の4.
1，4.
2
ではKMR法における水際での水深が，L=4固定格子に比べ乱れていたが，この条件のもとでは，そのような乱れは，ほとんど確認できない．

（a）緩やかな丘ありKMR法（t=15s）（b）緩やかな丘ありL=4固定格子（t=15s）

3――4（c）緩やかな丘ありKMR法（t=30s）（d）緩やかな丘ありL=4固定格子（t=30s）
図6 計算結果（iii）（4）計算時間の比較
計算時間を比較すると以下のような結果が得られた．これより，今回の計算モデルが氾濫計算において精度を犠牲にすることなく約45%の時間短縮をもたらすことがわかる．
（i） 河床勾配が一定である場合（KMRgrid）： （L=4grid）

=0.
56 ：

1（ii） 河床に水の浸入を妨げる障害物がある場合（KMRgrid）： （L=4grid）

=0.
55 ：

1=0.
55 ：

1（iii） 河床緩やかな丘があるような場合（KMRgrid）： （L=4grid）（5）明らかになった計算手法の課題
上記の検討を通じて，本計算方法に関する次の2つの課題が明らかになった．
a）格子パラメータに関する考察
本研究では格子分割パラメータを基準水深に対する水深比とした．結果研究の目的に掲げていた水際での格子分割を細かくすることに成功した．しかし，水際を局所的にみれば，全体的な精度は良好といえるが，上げて計算できているように見えるが，水際に若干の不安定性が確認された．下記の考察とも関連してくるがパラメータは，できるだけ，隣接するセルの流量が等しい状態で，かつ水際で格子分割や結合が行われる条件のものを選択すべきと推測される．また本ケースでは流量等から試行錯誤的に水深基準値を決める必要があり，条件の一般化が困難である．これらの点を踏まえて今後，格子分割パラメータについてさらに検討していく必要がある．
b）内挿と外挿について流れの先端部分に着目するとKMR法を用いた場合にはそれが不規則な形になっている．これは，KMR法における計算過程で内挿や外挿の近似方法が粗いためだと考えられる．本研究でのモデルでは，外挿はセル内定義点の値の単純平均で求め，内挿は内挿線形で求めている．流量フラックスの単純平均

4――4や線形内挿は流量の収支にずれが生じさせ，結果，流速に影響を及ぼす．それらが繰り返されることで，最終的には流れの先端部分に不規則な乱れが生じると考えられる．今後，内，外挿する際の近似方法に関しても検討していく必要がある．

5.実地形での計算本研究では，水際での計算に対しKMR法を適応させ，いくつか改良すべき点はあるものの，その精度がおおむね満足できるものであることを前項までに示した．ここでは，KMR法の実用性を確認するた
5）
め，実際に某河川の地形データを読み込んで試験的に計算を行った．河床データは，iRI
C2.
0上にあるサンプルデータを用いて行った．計算結果は図7，図8に示す．なお図7は立体的に見た図，図8は上空から見た平面的な図となっている．氾濫の進行に伴う格子分割の細分化が，水際付近で適切に行われていることがわかる．これらから，実地形における氾濫計算に対してもKMR法が適応できることが示された．

図7 鳥瞰図図8 上から見た図5――4

6.結論
本研究では，水深比をパラメータとし，KMR法を水際モデルに適応した．パラメータの選択と，格子分割，結合に伴う物理量の内挿方法と外挿方法に関しては，今後検討していく必要があるが，仮想水路における数値実験の高い計算精度と計算時間短縮は定量的に示された．また，前節より実地形での計算にも今回のモデルが適応できることが示せた．大規模氾濫計算の計算時間短縮は詳細なLPデータやXバンドデータなどの取得の容易化とともに工学的に有用性増しつつあり，KMR法はその一つの良好なアプローチである．今後は今回のモデルを基礎に，さらに精度を向上させるとともに，アルゴリズムの精査によるさらなる計算時間の短縮を目指していきたい．
参考文献
1） 馬場仁志・星 清・橋本識秀：損失機構を組み合わせた貯留関数モデルの総合化,水工学論文集,第43巻,pp.
1085
1090,1999.
（邦文論文）
2） Robinson,S.K.
:Coherentmotionsintheturbulentboundarylayer,Ann.Rev.Fl
uidMech.
,Vol.
23,pp.
60
1639,1991.
（英文論文：論文題名の最初の単語は大文字で始まる．雑誌名は斜体）
3） Tennekes,H.andLu
mley,J.
L.
:Afirstcou
rsei
nturbulence,TheMI
TPress,1972.
（英文書籍：題は斜体）
4） 土木学会編:土木工学における数値解析，流体解析編，サイエンス社，1974.
（邦文書籍）
5） 齊藤真治・木村一郎・清水康行:KMR法を用いた水深積分モデルによる橋脚周辺の非定常構造の再現計算,
水工学論文集,第56巻,2012.
（邦文論文）
6） 木村一郎・細田:開水路せん断混合層の流れ構造に対する水深積分型修正ゼロ方程式モデルの適用性,土木学会論文集,第49巻,pp.
559564,2005.
（邦文論文）
7） 安田浩保・星野剛:四分木構造格子による局所的な高解像度格子を導入した浅水流方程式の数値解
析法,土木学会論文集A2（応用力学）,Vol.
67,No.
2（応用力学論文集）,I_693I_702,2011.
（邦文論文）
8） 小玉務・渡辺靖憲：Lev
elset法を用いたAd
aptiv
eMeshRe
finementによる三次元波浪計算,土木学会
北海道支部論文報告集,第66号,B49,2009.
（邦文論文）9） i
RI
Cホームページ,htt
p://ir
ic.
org/
ja.
/
（ホームページ）

6――4