ＩＳＳＮ１３４６－７８４０

港湾空港技術研究所
資料
TECHNICAL NOTEOFTHE PORT AND AIRPORT RESEARCH INSTITUTE

No.1411

September 2023

非構造格子を用いた計算手法の土砂輸送解析への導入の試み

大竹 剛史，中川 康之，小硲

国立研究開発法人

大地

海上・港湾・航空技術研究所

National Institute of Maritime,
Port and Aviation Technology, Japan目

次旨 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 3

要

1. はじめに ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 4
1.1

研究の背景と目的 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 4

1.2

オープンソースモデルの比較 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 4

1.3

有限体積法の基礎理論と非構造格子の取り扱い ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 51.4非構造格子における有限体積法の離散化 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 7

2. 格子の設定条件の差異が計算結果に及ぼす影響 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 9
3. 実地形を用いた土砂輸送数値実験 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 11
3.1

計算条件 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 11

3.2

計算結果および考察 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 13

4. 防砂堤設置港形における土砂輸送数値実験 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 14
4.1

計算条件 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 14

4.2

計算結果および考察 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 14

5. おわりに ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 15
謝辞 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 15
参考文献 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 15
記号表 ꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏꞏ 16-1-

Introducing the Calculation Method Using Unstructured Grids for SedimentTransport Analysis
Takeshi OTAKE*
Yasuyuki NAKAGAWA **Taichi KOSAKO***

Synopsis
One of the causes of errors in numerical calculations of flow is the uneven and discontinuous
boundaries when using structured grids. Therefore, unstructured grids are oftensometimes used.
Therefore, unstructured grids are often used. On the other hand, the effects of using
unstructured grids on topographical changes are not well understood, because there are few
studies on sediment transport calculations using unstructured grids. Therefore, in this study,
with the purpose of introducing sediment transport analysis into numerical calculations using
unstructured grids, we investigated the characteristics of topographic changes in unstructured
grid models through comparison of numerical calculations of unstructured grids and structured
grids.
The finite volume method among the discretization methods for unstructured grid models istaken up, and the basic theory and discretization method are summarized. We conductednumerical wave experiments on structured and unstructured grid models, and summarized thefeatures and considerations of unstructured grid models. As a result, when using an
unstructured grid, the orthogonality is important, and the treatment of the advection term andthe horizontal viscosity term is different between the unstructured grid model and thestructured grid model, which affects the results. Next, using the open software model
SUNTANS, we performed sediment transport calculations on the actual landscape around thePatimban port in Indonesia. The results were compared with the structured grid model to show
the features in sediment transport calculations. The calculation results show that the differencein the grid shape causes the difference in the flow along the coastline, which affects the rangeof the erosion area. We also examined the case of installing a sand protection dyke that blocks
the tidal current. As a result, it was found that there is a difference in the range of variation ofcurrent velocity, which affects the variation of topography.
Key Words: Sediment transport model, Unstructured grid, Finite Volume method* Reserch trainee, Coastal and Estuarine Sediment Dynamics Group, Coastal and Estuarine Environment
Department (Ecoh Co., Ltd.)** Senior Director for Research*** Senior Researcher, Coastal and Estuarine Sediment Dynamics Group, Coastal and Estuarine Environment
Department
3-1-1 Nagase, Yokosuka, 239-0826 Japan
Phone：+81-46-844-5045 Fax：+81-46-844-1274
e-mail:nakagawa@p.mpat.go.jp-2-非構造格子を用いた計算手法の土砂輸送解析への導入の試み
大竹剛史*・中川康之**・小硲

要

大地***

旨

正方格子に代表される構造格子を用いた流動場や土砂輸送等の数値シミュレーションでは，海岸線や防波堤などの施設形状が，格子サイズに依存して凹凸のある不連続な境界として扱われてしまうため，土砂輸送量に変化を与える可能性があり，地形や施設形状に沿って格子形状を任意に変更可能な非構造格子モデルが導入される場合がある．一方，非構造格子を用いて土砂輸送計算が行われた研究事例はあるが，非構造格子を用いることによる地形変化への影響はあまり分かっていない．
そこで本研究では，非構造格子モデルを用いた土砂輸送計算の特徴を把握するため，非構造格子モデルによる数値シミュレーション手法の把握および構造格子モデルとの比較を通した，格子形状が異なる条件の下での流れ場および地形変化の計算結果に与える影響について評価する事を目的とした．
非構造格子モデルの離散化手法のうち有限体積法を取り上げて基礎的な理論と離散化手法について整理した．整理した手法を用いて数値波動水槽を用いた簡易な数値実験を実施し，構造格子モデ
ルと非構造格子モデルでの比較を通じて，非構造格子モデルの特徴と留意点について整理した．その結果，非構造格子を取り扱う際の直交性の重要性や移流項および水平粘性項の取り扱いの違いによる結果の差異を示した．三角形格子の中心点に外心を用いることにより直交性を確保しているオープンソフトウェアモデル SUNTANS を用いて，インドネシア・パティンバン港周辺を対象とした実地形での土砂輸送計算を実施し，構造格子モデルとの比較を通じて土砂輸送計算における特徴を整理した．その結果，海岸線に沿う流れに違いが見られ，侵食域の範囲に影響を与えることが分かった．さらに，潮流を遮るような防砂堤を設置したケースについて検討を行い，流動場と土砂輸送への影響に注目したところ，流速の変化量の範囲に差異が見られ，地形変化量に影響を与えることがわかった．
キーワード：土砂輸送モデル，非構造格子，有限体積法* 沿岸環境研究領域 沿岸土砂管理研究グループ 元・依頼研修員（株式会社エコー）** 特別研究主幹（前・沿岸土砂管理研究グループ長）*** 沿岸環境研究領域 沿岸土砂管理研究グループ 主任研究官〒239-0826 横須賀市長瀬3-1-1 国立研究開発法人海上・港湾・航空技術研究所 港湾空港技術研究所
電話：046-844-5045 Fax：046-844-1274
e-mail:nakagawa@p.mpat.go.jp-3-

1. はじめに理した（1章）．整理した手法を用いて数値波動水槽を用いた簡易な数値実験を実施し，構造格子モデルと非構造

1.1

研究の背景と目的格子モデルでの比較を通じて，非構造格子モデルの特徴

正方格子に代表される構造格子を用いた流動場や土砂と留意点について整理した（2章）．次に，オープンソフ

輸送等の数値シミュレーションでは，海岸線や防波堤な

トウェアモデルSUNTANSを用いて，インドネシア・パ

どの施設形状が，格子サイズに依存して凹凸のある不連

ティンバン港周辺を対象とした実地形での土砂輸送計算

続な境界として扱われてしまう．この凹凸が海陸境界付を実施し，構造格子モデルとの比較を通じて土砂輸送計

近の流速に影響を与え，浮遊土砂の拡散場に影響を与え

算における特徴を整理した（3章）．また，予測計算時にるのみならず，浮遊土砂の巻き上げ・沈降や掃流砂量に流れを遮るような構造物が設置された場合の流動場と土関連する底面せん断力の評価にも影響を与える等，土砂砂輸送への影響を把握するため，3章で実施した港形に輸送量に変化を与える可能性がある．これら格子サイズ

防砂堤を設置して，その影響を確認することで非構造格や凹凸の影響を軽減するため，非構造格子モデルが導入子の特徴について整理した（4章）．

される場合がある．非構造格子は，三角形格子に代表されるように格子形状を任意に変更することが可能であり，

1.2

複雑な海岸形状や港湾施設形状に沿った格子を配置する

非構造格子モデルの流体モデルは，前節に使用例を示ことや，格子サイズを沖から対象付近にかけて徐々に小したようにオープンソース化されている代表的なモデル

オープンソースモデルの比較

さくすることが可能である．一方，非構造格子を用いたがいくつか公開されている．表-1.1に代表的なオープン

数値シミュレーションでは，構造格子モデルとは基礎方

ソースモデルの特徴を示す．表に示すように，それぞれ

程式を離散化する手法が異なり，このことが流体場や土偏微分方程式の離散化手法や鉛直方向のメッシュの切り

砂輸送計算へどのような影響を与えるかを把握する必要方が異なる．
SUNTANS，FVCOM，D-FLOW FM(Flexible Mesh)9) は，があると考えられる．
非構造格子モデルを用いた検討としては，有限体積法

任意の検査体積領域の体積積分をガウスの発散定理を用を採用した非構造格子海洋流動モデルであるFVCOM1)を用いて表面積分に変換し，その流束を評価する有限体積法

いて，高潮や津波計算に適用した例が報告されているを採用している．一方，SCHISMの流体計算は，領域を多

2)

3)

4)

（例えば小松ら ，木梨ら ）．増永ら は，やはり同様なくの要素に分割し，要素毎に多項式で近似して求める有5)

オープンソースモデルの一つであるSUNTANS を用いて，限要素法を採用している．

伊豆大島周辺海域における潮流計算を，内部潮汐の影響

鉛直方向のメッシュについては，それぞれ異なる手法も考慮して検討している．また，中谷ら 6) は流動モデルを採用しており，SUTNANSは水平方向に一様に層を分に有限要素法を採用した非構造格子モデルSCHISM7)を用ける鉛直レベルメッシュを用いている．D-FLOW FMは

いて瀬戸内海-太平洋領域の外用水の挙動について黒潮

σ座標を採用しており，海面から海底までを等間隔にメの蛇行の影響を考慮した検討を実施している．一方，地ッシュ切りすることが可能である．これにより，地形変形変化計算では，山縣ら8)がFVCOMを用いて，潮位偏差化が急峻な領域に対しても，鉛直メッシュ数を減らさずを考慮した関門航路のサンドウェーブについての地形変に計算することが可能である．FVCOMの鉛直方向メッ

化シミュレーションを検討しているものの，非構造格子

シュ切りにみられるS-σハイブリッド座標のSとは，スを用いることによる地形変化への影響はあまり分かってトレッチの略称であり，σ座標において等間隔である各いない．

層厚を，任意の層厚間隔の割合で表現する手法である．

そこで本研究では，非構造格子モデルを用いた土砂輸

SCHISMの鉛直方向メッシュ切りにみられるS-Zハイブリ

送計算の特徴を把握するため，非構造格子モデルによるッド座標のSも同様の意味であり，Zは鉛直レベルメッシ

数値シミュレーション手法の把握および構造格子モデル

ュのことである．浅い領域で鉛直レベルメッシュに切りとの比較を通した，格子形状が異なる条件の下での流れ替えることにより，層厚が薄くなりすぎてしまうことを場および地形変化の計算結果に与える影響について評価防ぐことが可能である．

する事を目的とした．本研究では，まず非構造格子モデ

近年，それぞれモデルにおいて徐々に土砂輸送計算モルの現在公開されているオープンソースモデルについて

ジュールの実装も進みつつあるものの，論文として公表

整理し，それらで扱われている離散化手法のうち有限体

されている最新の土砂輸送計算が公開モデルには実装さ積法を取り上げて基礎的な理論と離散化手法について整れていないSUNTANSや，有料版に先行で公開されてい-4-

表-1.1非構造格子を用いたオープンソースモデルの比較

モデル名
開発グループ離散化手法

SUNTANS5)
Fringer ら
有限体積法

Chen ら
有限体積法

近似方法

静水圧近似または
非静水圧を選択静水圧近似または
非静水圧を選択

Z 方向メッシュ
ライセンス
土砂輸送計算モジュール

ブシネスク近似
鉛直レベルメッシュGPL(ver.2)に準ずるあり（浮遊成分のみ）

FVCOM1)SCHISM7)
Zhang ら流動場 ：有限要素法
物質輸送：有限体積法
静水圧近似

S-σハイブリッド座標
商用利用は有償
開発者らへの報告義務ありるDeltaresの地形変化モジュール等，そのオープンソース

D-FLOW FM9)
Deltares
有限体積法
静水圧近似

ブシネスク近似
S-Z ハイブリッド座標
Apache License

σ座標
GPL(ver.3)またはLGPL(ver.2)に準ずるありあり

選択する必要がある．

化の対応も異なっている．
オープンソースモデルの使用にあたっては，それぞれ

1.3

異なるライセンスが定められているため注意が必要であ

有限体積法の基礎理論と非構造格子の取り扱い(1) 有限体積法による保存式の導出

る．SUNTANSで用いられているGPL(GNU General Public

非構造格子に基づく偏微分方程式の解法は，有限体積

License)は，フリーソフトウェアマス・コラボレーショ

法，もしくは有限要素法によって実装される．ここでは，ンプロジェクトの1つであるGNUプロジェクトのために

多くのオープンソフトウェアモデルに採用され，質量保作成されたフリーソフトウェアライセンスである存則が数値計算上満たされる利点を有する有限体積法に(https://www.gnu.org/licenses/licenses.en.html)．GPLは代表

ついて説明する．

的なコピーレフトのソフトウェアライセンスの1つであ

質量や運動量などの任意の保存量φに対する保存式のり，プログラムの実行，プログラムの動作確認および改導出は，図-1.1のような検査体積（Control Volume）を仮変，複製物の再頒布，プログラムの改良および公衆にリ定して導かれる．dt時間に検査体積の曲面から流出する

リースする権利が許諾される．また，商用利用可能なこ

φの正味の量と，検査体積内のφの増減の総量は保存さとや利用に伴い作成者がいかなる責任も負わないことなれるとして保存式は式(1)に示すように積分形で表されどが特徴である．GPLがversion2.0からversion3.0にアップ

る．

デートされるにあたり，特許条項が追加されている．D-⬚

FLOW FMはソースコードの範囲により，GPLとLGPLが使用されている．LGPL(Lesser GPL)とは，ライブラリな𝐶𝑉

どの他のプログラムからリンクされる形で使用されるソ

𝜕
𝜕𝑡

⬚

𝜌𝜙 𝑑𝑉𝑑𝑡

𝜌𝜙 𝒗 ∙ 𝒏𝑑𝐴𝑑𝑡

0

1

𝐴

フトウェア向けのライセンスであり，GPLでは部分的な

ここで，ρは流体密度，uは流体の速度ベクトル，nは単プログラムを組み入れるだけでもソフトウェア全体を公位法線ベクトル，tは時間，CVは検査体積，Vは体積，A

開しなければいけないという条件を緩和したものである．

は検査体積の表面積である．式(1)にガウスの発散定理を

SCHISMで用いられているApache Licenseは，Apacheソフ

適用して表面積分を発散の体積積分に書き直し，検査体

トウェア財団によるライセンス規定であり，GPLとの互積の大きさをゼロに近づけてその極限を取ると，式(2)の換 性 を 有 し て い る (https://httpd.apache.org/) ． 一 方 ，ように一般的な任意の保存量φに対する保存式が導出さ

FVCOMは，GPLやApache Licenseのような汎用的なライれ，有限差分法等で解かれる．

10)

センスではなく，独自のライセンスを定めている ．利

𝜕
𝜌𝜙
𝜕𝑡

用者は，FVCOMのソースコードおよびその改変を，非商業利用に限り認められ，第三者への配布または譲渡も許𝑑𝑖𝑣𝜌𝜙𝒖0

2

可されていない．このように，オープンソースモデルを

有限体積法は，はじめに仮定した離散空間の検査体積が選定するにあたっては，その利用目的により，計算手法

有限のものであるので，もとの積分形のまま扱うというや計算目的に限らず，ライセンス的な側面からも適切に思想である11)．

-5-

φ
n
dS
Control Volume

図-1.4 コントロールボリュームの定義方法（a）セル中心 CV，(b)接点中心 CV構造が存在せず，図-1.3に示すように格子の特徴点のナ

図-1.1

ンバリングを不規則に設定することができる格子を指す．

保存則導出時のコントロールボリューム

非構造格子の検査体積（コントロールボリューム：nム（cell-centred control volume）と(b) 節点中心コントロ・・・CV）の定義方法は，(a)セル中心コントロールボリュー
ールボリュームの2つの方法がある．図-1.4に2次元問題

3ュームでは，格子点をコントロールボリュームの中心にに対するこれら2つを示す．セル中心コントロールボリ
配置する．それに対して，節点中心コントロールボリュ

2

ームでは，計算格子の節点に格子点を配置する．この方法では，要素の中心と節点の間の中点を結び，サブボリ

j=1

ュームを生成する過程に従う．その後，離散化するため

y

図-1.2

i=1

2

3

・・・ m

に，格子点を取り囲むサブボリュームでコントロールボx

リュームを生成する．本研究では，理解しやすいセル中保存則導出時のコントロールボリューム

心コントロールボリュームでの検討を進める．ここで，格子の中心に定義される格子点はcell center，格子の辺はfaceまたはedge，格子の節点はpointやnodeと呼ばれる．

8以降，cell centerは省略してcell，格子の辺をface，格子のface 2341117

構造格子は，座標軸に沿った方向に対して絶対的な構造が存在するため，離散化計算を行う際に必要な格子同2999cell center
2
41
32

109

士の繋がり情報は，順番に並んだ番号を参照することで15

隣り合う格子の計算を行うことができる．一方，非構造
格子はこうした構造が存在しないため，コネクティビテ

12

182
node 11

節点をnodeと表記した．

65

56

ィと呼ばれる格子情報同士の関係をメモリに保存してお76く必要がある．例えば以下のようなコネクティビティを保存することで格子同士の情報を参照することができる．

図-1.3 非構造格子における変数格納点（赤:辺（face），青:中心点（cell），緑:頂点（node））

・cell-to-faceコネクティビティ
・face-to-cellコネクティビティ(2) 非構造格子による有限体積法計算

・face-to-nodeコネクティビティ

構造格子とは，座標軸に沿った方向に対して絶対的な

cell-to-faceコネクティビティでは，cellの番号および構造が存在する格子のことである．例えば，x，y方向のcellにおけるローカルのface番号（三角形の場合は1～3座標軸を持つ2次元のスタッガード正方格子では，図-個）を指定することで，faceの情報が返る変数である．

1.2に示すように，それぞれの座標方向に対し構造的なcell-to-faceの配列要素を(cell番号,ローカルface番号)としナンバリングがされる．一方，非構造格子は，それらのて，図-1.3のcell(32)を例に取ると，下記のように3辺の-6-face番号を得ることができる．
⊿A・cell-to-face(32,1)=99⊿A・cell-to-face(32,2)=41

ni・cell-to-face(32,3)=2

ni

⊿A

ni

同様に，face-to-cellコネクティビティは，face番号およびfaceにおけるローカルcell番号(通常2個)を指定すること
図-1.5

で，cellの情報が返る変数である．face-to-cellの配列要素

非構造格子による多角形メッシュの例

を(face番号,ローカルcell番号)として，図-1.3のface(2)を例に取ると，下記のように2つのcell番号を得ることがで

b

b

きる．
・face-to-cell(2,1)=32
・face-to-cell(2,2)=12AP

face-to-nodeコネクティビティは，face番号およびfaceにおけるローカルnode番号(通常2個)を指定することで，

inodeの情報が返る変数である．face-to-nodeの配列要素を

⾯(i)

a

(face番号,ローカルnode番号)とすると，図-1.3のface(2)をa

図-1.6 左図：セル中心 CV
右図：face i における単位法線ベクトル例に取ると，下記のように2つのnode番号を得ることができる．
・face-to-cell(2,1)=56ように表される．

・face-to-cell(2,2)=29
1.4

𝜕
𝜕𝑡

非構造格子における有限体積法の離散化

⬚

⬚

𝜌𝜙𝑑𝑉

𝒏𝒊 ∙ 𝜌𝜙𝒖 𝑑𝐴
512)

ここからはH.K.Versteegら に基づいて，非構造格子に⬚⬚

𝒏𝒊 ∙ 𝛤𝑔𝑟𝑎𝑑𝜙
grad 𝑑𝐴

おける有限体積法の離散化について説明する．有限体積

𝑆 𝑑𝑉

法を非構造格子に適用する出発点として，式(3)に示す積面 積 分 を 計 算 す る た め に は ， 流 速 ベ ク ト ル ρφu と分形の保存式を考える．
⬚

𝜕
𝜌𝜙 𝑑𝑉
𝜕𝑡

Γgradφに加えて，幾何学量ni とΔAi が必要となる．外向⬚きベクトルniとΔAiは，非構造格子のnodeの座標から，単𝑑𝑖𝑣 𝜌𝜙𝒖 𝑑𝑉

⬚

3

⬚

𝑑𝑖𝑣 𝛤𝑔𝑟𝑎𝑑𝜙 𝑑𝑉

純な三角法とベクトル代数を用いて計算する．図-1.6に

𝑆 𝑑𝑉

一般的なセル中心コントロールボリュームとその表面における単位法線ベクトルを示し，図中に離散化手順を記

ここに，Γは拡散係数，Sφはφの生成量である．また，述するために用いる記号を示した．図中の点Pは離散化

積分範囲CVは検査体積を表す．検査体積の境界表面をAを行うコントロールボリュームの中心（cell center）であとし，ガウスの発散定理を適用するとる．点Aは隣接するcell centerであり，eξは点Pと点Aを結ぶ線に平行な単位ベクトルである．2つのコントロール

⬚

𝜕
𝜕𝑡

⬚

𝜌𝜙𝑑𝑉

ボリュームを分ける面(face)をiと定義し，abは2つのコン𝒏 ∙ 𝜌𝜙𝒖 𝑑𝐴⬚𝒏 ∙ 𝛤𝑔𝑟𝑎𝑑𝜙 𝑑𝐴

トロールボリュームを共有する節点(node)である点aと点

4

⬚

𝑆 𝑑𝑉

bを結ぶ線である．点aと点bの座標はそれぞれ(xa,xb)と(ya,yb)である．単位ベクトルnとeηはそれぞれ面iの外向きと表せられ，Aにおいて面積分を実行する必要がある．

の単位ベクトルと接線ベクトルである．

異なる形状の検査体積を図-1.5に示す．これら検査体同図より，face iにおける面積は式(6)，face iに対する積それぞれの境界表面は，有限の長さの直線要素で構成

単位法線ベクトルは式(7)のようになる．ただし，Δx=xb

される閉輪郭であり，その面積をΔAで表す．式(4)中のA− xaおよびΔy=yb − yaである．

は検査体積表面の全領域，dAは微小表面積を意味するの∆𝐴𝑖で，直線要素すべてに対して面積分を実行すると式(5)の

-7-

∆𝑥2

∆𝑦2

6

𝒏

∆𝑥
∆𝐴𝑖

𝒊

∆𝑦
∆𝐴𝑖

項を導入し，離散化方程式を構築する際に生成項として𝟕𝒋取り扱う方法がある．詳しくはMathur and Murthy13)を参a）拡散項の離散化

考されたし．

式(5)右辺第1項の拡散項を離散化する．要素それぞれ

三角形格子に対する非直交性の問題に対するもう一つの面積分を，コントロールボリュームの表面要素ΔAi にの留意点としては，cell centerの位置がある．三角形格子

対する外向きの単位法線ベクトルni と拡散流速ベクトルのcell centerとして一般的なものには重心と外心が知られ(Γgradφ)の内積で近似する．線PAに沿う中心差分法を用る．図-1.7に2種類のcell centerを示す．図で示すように，いて，式(8)のように近似できる．

外心は各辺の垂直二等分線が必ず通る．したがって，隣⬚𝒏 ∙ 𝛤𝑔𝑟𝑎𝑑𝜙 𝑑𝐴 ≅ 𝛤
∆𝐴𝑖

𝑥𝑏

𝑥𝑎
∆𝜉

り合う格子とは辺を共有するため直交性が保たれる．し∆𝐴𝑖

8

かしながら，内角に直角，または鈍角が存在する場合は外心が三角形の辺上，または外側になり，格子の代表点

式中のΔξは，点A点Pの間の距離である．点Pと点Aを結として扱うことができない．図-1.8に直角二等辺三角形

ぶ線と単位法線ベクトルni が同じ方向同じ方向である場の場合の例を示す．一方，重心は必ず三角形の内側に存

合のみ，式(8)は正確であることに注意が必要である．そ在する．この場合は非直交性が保証されないので，前述のため，この近似は格子が完全に直行している場合のみ

した交差拡散項の導入等を導入する必要がある．外心と正確である．図-1.6に示すように，非構造格子では一般重心が一致する場合，格子形状は正三角形となり，外心に点Pと点Aを結ぶ線とベクトルは平行ではない．これは，の垂直二等分線も全ての辺から等距離になる．このため，計算格子のゆがみ，あるいは非直交性として知られてい直交性を保証する格子三角形は理想的に正三角形となる．

る．そのため，流束の計算の式(8)に，非直交性に起因すb）対流項の離散化

る寄与を加えることで修正しなければならない．流束を式(5)左辺第2項の対流項を離散化する．コントロール

修正する方法の1つとして，交差拡散項として知られる

ボリュームの表面要素ΔAi においてすべて積分し，和をとるように面積分を計算する．外向きの単位法線ベクト
ルniとコントロールボリュームの表面要素ΔAiを乗じた対流流束ベクトル(ρφu)の内積で，これらの積分をそれぞれ式(9)のように近似する．
⬚

𝒏𝒊 ∙ 𝜌𝜙𝒖 𝑑𝐴
𝑎𝑙𝑙 𝑠𝑢𝑟𝑓𝑎𝑐𝑒

𝒏𝒊 ∙ 𝜌𝒖 𝛥𝐴𝑖 𝜙𝑖∆𝐴𝑖

9

𝑎𝑙𝑙 𝑠𝑢𝑟𝑓𝑎𝑐𝑒

ここで，対象としている任意の保存量φがSS濃度である移流拡散方程式を例とした場合，cell centerに保存されている保存量φをface iにおけるSS濃度へ補完する必要がある．代表的な方法の1つに距離の重み付け補完法があ

図-1.7 2 つの cell center の定義
赤色塗りつぶし点：重心，青色白抜き点：外心る．図-1.9に示す2つの格子に挟まれた面iにおける保存量φを求める場合，点1および点2の値と，それぞれの面iまでの距離d1，d2を用いて式(10)のように重み付け補完を
行う．
𝜙𝑓

𝑤𝑓 𝜙2

1

𝑤𝑓 𝜙1

10

ここに，wf は距離の重み関数であり，式(11)のように表す．
𝑤

図-1.8 直角二等辺三角形における cell center
色塗りつぶし点：重心，青色白抜き点：外心

1
𝑑
1
𝑑

1
𝑑

11

対象としている任意の保存量φが流束ベクトルuである運動量保存式を例とした場合かつ，faceに流速ベクト-8-式(15)中のRについて，a)～c)で示したように各項を離散化し，φt+Δtについて毎ステップ解くことで，時々刻々の
解を得る．
𝑑

𝑑
𝑖

1

2. 格子の設定条件の差異が計算結果に及ぼす影響

2

1章で示した非構造格子における有限体積法を用いた
簡略化した連続式と運動方程式を用いて流体計算を実施
し，構造格子との取り扱いの違いとそれに伴う計算結果

図-1.9 cell 変数の face への補完方法への影響について確認する．検討は幅10m，距離1,000m
水路を仮定し，図-2.1に示す構造格子による有限差分法

ルが保存されている スタッガード 格子の場合，格子の各(FDM) と 図 -2.2 に 示 す 非 構 造 格 子 に よ る 有 限 体 積 法

faceに保存されている流速からcell centerにおけるuを補完(FVM)で比較する．水路の片端において水位を振動させする必要がある．例えばSUNTANSやD-FLOW FMでは，る造波境界を設け，ここでは水位変動の伝搬過程（波動14)perot による方法によりこの解法を得ている．

場）の計算を実施した．使用した連続式と運動方程式を以下に示す．

c）生成項の離散化式(5)の右辺第2項は生成項であり，式(12)のように取り
扱う．
⬚

𝑆𝜙 𝑑𝑉

12

𝑆𝜙 𝑑𝑉

𝜕𝑀
𝜕𝑡

𝑐𝑣

𝜕𝜂𝑙

𝜕𝑀

𝜕𝑁

𝜕𝑡

𝜕𝑥

𝜕𝑦

𝜕𝜂
𝜕𝑥

g𝑛 𝑀√𝑀𝜕𝜂𝑙g𝑛2 𝑁√𝑀2

g𝐷

16

0

𝑁

17

𝐷

￣ はコ
ここに，ΔVはコントロールボリュームの体積，Sφ𝜕𝑁

ントロールボリュームにおけるSφの平均値である．

𝜕𝑡d）非定常項の離散化

g𝐷𝜕𝑦

式(5)の左辺第1項は非定常項であり，時々刻々の変化

𝑁2

7

18

𝐷3

を計算する場合には，離散化の過程でこの項を取り扱わ

ここでηlは水位，M，Nはそれぞれx，y方向の線流量，gなければならない．式(5)の非定常項以外の項を右辺に集は重力加速度，Dは水深（ここでは10mの一定とする），

めRと置き，式(13)のように表す．

nは祖度係数である．運動量は水位勾配で決まるとし，

𝜕𝑡

底面摩擦によって減衰するとした．

⬚

𝜕

𝜌𝜙𝑑𝑉

13

𝑅

有限差分法の計算を実施した図-2.1の構造格子におい𝑐𝑣て，iは格子番号であり，i=1において以下の水位を与え式(13)をコントロールボリュームにおいて有限の体積でる．格子幅dxを10mとし100格子配置した．また，計算間積分することに加え，さらに有限の時間刻みΔtにおいて

隔dtは0.001sとした．

積分する．非定常項の積分の順序を入れ替えて式(14)を

𝜂

得る．
⬚
𝐶𝑉

𝑡 𝛥𝑡
𝑡

𝜕
𝜕𝑡

𝛼 sin

2𝜋𝑡
𝑇

19

𝑡 𝛥𝑡

𝜌𝜙 𝑑𝑡 𝑑𝑉𝑅𝑑𝑡

14

ここに，αは振幅， T は周期であり，本計算ではα𝑡=1.0m，T=20s とした．計算時間は 200s とし，10 波分の格子点における保存量φがコントロールボリュームにお波を入射させた．ηl を中心格子，M，N を格子界面に配いて全体で一様であると仮定し，陽的に離散化を行うと置するスタッガード格子で計算を実施し，流入，流出境

式(15)のように書くことができる．

界部の線流量（流速）を 0 とする壁面境界とした．

𝜌 𝜙𝑡 𝛥𝑡

𝜙𝑡 𝛥𝑉

𝑅𝛥𝑡

有限体積法を用いた非構造格子計算の格子形状を，有

15

-9-

l

0

0

25

25

50

50

・・・
dx=10m

i=1

2

99

100

75
t(s)

t(s)

75
100
125

125

150

150

175

175

0

200

有限差分法における計算構造格子

400

600

800

0

200

400

600

800

x(m)
x(m)

x(m)
x(m)

-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0-2.0-1.5-1.0-0.5 0.0 0.5 1.0 1.5 2.0-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0-2.0-1.5-1.0-0.5 0.0 0.5 1.0 1.5 2.0

ηeta(m)l(m)ηl(m)
eta(m)

1000m

図-2.1

100

図-2.3 水位時間変化の比較（左：有限差分法，右：有限体積法）

l

2

4

198 200
・・・

10m

Δ𝜉

Δ𝜉
𝑀

𝑢

10m

cell=1 3𝑦

197 199

Δ𝐴
𝑥

1000m

図-2.2

図-2.4
限差分法に用いた矩形格子と同様の形状にすると，両モ
デルの離散化は同様の式が導出され結果が等しくなるため，有限差分法と全く同じ離散化式が導出されるため，

𝑖

Δ𝑥
2

有限体積法における計算非構造格子

𝑀

Δ𝐴
Δ𝑥
2

0

0

25

25

50

50
75

100

100

125

125

150

150

175

175

ルの水路長を有限差分法と同等とした．なお，有限差分0

法では 1 次元の計算であったが，有限体積法では三角形格子を用いたため y 方向成分が出てしまい，平面 2 次元Δ𝑥Δ𝑦有限体積法に使用する矩形格子

t(s)

t(s)

10m の直角二等辺三角形格子を 200 格子配置し，トータ

Δ𝜉
Δ𝐴

※Δ𝜉
※Δ𝐴

75

図-2.2 に示す三角格子で検討を行った．幅 10m，高さ

Δ𝑦

𝑗

200

400

600

800

0

200

400

600

800

x(m)
x(m)

x(m)
x(m)

2.0
0.5 1.0
0.0 0.5-2.0 -1.5 -1.0 -0.5 0.0
1.5 2.0
1.0 1.5-2.0-1.5-1.0-0.5

2.0
0.0 0.5
0.5 1.0-2.0 -1.5 -1.0 -0.5 0.0
1.0 1.5
1.5 2.0-2.0-1.5-1.0-0.5ηeta(m)l(m)ηl(m)
eta(m)

計算となっている．cell=1 において有限差分と同様に造図-2.5 水位時間変化の比較（矩形格子）（左：有限差分法，右：有限体積法）

波しており，定義点がずれている．その他の条件については有限差分法と同様とした．
図-2.3に水位分布の時系列を示す．縦軸は経過時間，

1.5

横軸は造波地点からの距離である．左図が有限差分法，

1

右図が有限体積法の結果であり，水位をコンターで表し

0.5

ηl (m)ている．結果を見ると，右端境界において流出を0としているため，完全反射が起きていることがわかる．図中

0-0.5に，計算開始100秒時の破線を示している．有限差分法-1では，100秒時に1波目が到着しているが，有限体積法で

FDM

-1.5

は到達しておらず，波速は有限体積法の方が遅くなって

0

200

400

600

FVM

800

x(m)

いる．この原因の1つとして，cell centerを重心に設けた
図-2.6

ため，前述したようにcell centerを結ぶ線と法線の方向が水位の分布（矩形格子）

一致していないために現れた誤差（非直交性）だと考え

格子を図-2.4のように変更した．ただし，格子形状は同られる．

じ正方格子ではあるが，格子情報は非構造形式としてい有限体積法で同様の格子および中点近似を用いると，る．ここで，iは計算対象となる格子番号であり，格子同有限差分法と全く同じ離散化式が導出されるが，移流項

士の計算を行う際に，面を挟んで対になる格子番号をjとおよび水平粘性項の影響を把握するため，有限体積法のして表す．図-2.5に水位分布の時系列を示す．縦軸は経- 10 -

過時間，横軸は造波地点からの距離である．正方格子にているパティンバン港周辺の沿岸域を対象とし，潮流にしたことで，有限体積法においても100秒時に1波目が到よる地形変化計算の数値実験を実施した．計算に使用す着している．また，図-2.5の100秒時（破線時）における港湾形状について，パティンバン港は現在建設途中でるx方向の水位分布を図-2.6に示す．

あるが，2017年に実施された開発調査15)の計画港形を本

移流項および水平粘性項による違いを比較するため，研究の基本港形として計算を実施することとした(図-14)

式(18)に移流項および水平粘性項を追加し，perot を参3.1(b))．

考に離散化を行った．用いた格子形状は図-2.4と同様で本研究の計算対象とする沿岸域での現地データは地ある．移流項は，対象となる格子まわりのfaceそれぞれ

点・期間とも限られている中で，小硲ら16)は当海岸を対次式のように計算した．

象とする地形変化計算を実施する際に，潮位の現地観測
データと比較を通じて，使用した構造格子モデルによる

1
𝑉

潮位変動の再現性については，その妥当性を確認してい20

𝒏・𝒖𝑖 𝑀𝛥𝐴

る．本研究では，この小硲らの構造格子モデルと非構造

𝑎𝑙𝑙 𝑠𝑢𝑟𝑓𝑎𝑐𝑒

格子モデルによる計算結果の比較を通じて，両者の差異
ここで，ui はcell centerにおける値なので，faceに保存さを把握・評価することを試みた．非構造格子モデルは，れている各流速を用いて以下のように内挿する．

局所的な質量保存則が数値計算上保証される利点を有す
る有限体積法を採用しており，ソースコードの改良後に

𝑢𝑖

1
2

𝑢1 𝛥𝜉1おける利用上の制約等も考慮してSUNTANSを使用する

21

𝑢2 𝛥𝜉2こととした．非構造格子モデルでは，図-3.1(a)に示す計
算領域の沖合で最大格子間距離1,080m，港周辺で図-

水平粘性項についても同様に，対象となる格子まわりの3.2(b)のように港湾形状に合わせて格子を細かくする可faceについて次式のように計算する．

変メッシュとした．一方，比較用の構造格子モデルでは，𝒏・𝜈𝐷
𝑎𝑙𝑙 𝑠𝑢𝑟𝑓𝑎𝑐𝑒

𝑀𝑗

𝑀𝑖
𝛥𝜉

𝛥𝐴

22

ここで，νは動粘性係数である．これらのように求めた

0

0

25

25

50

50

t(s)

75

移流項，水平粘性項の値は，対象となる格子のcell center
で求めている．流速計算のためにfaceでの値とするためには，図-1.4で示した補完を行う必要がある．

75
t(s)

1
𝑉

100

100

125

125

150

150

175

175
0

図-2.7の水位の時間変化を見ると，両モデルに違いが

200

400

600

800

0

200

x(m)
x(m)

出ており，有限体積法の方がピークの水位が高くなって
-0.5
-0.5

-1.0
-1.0

いる．図-2.7中の100秒時（破線）におけるx方向の水位

0.0
0.0

0.5
0.5

図-2.7

600

800

x(m)
x(m)
1.0
1.0

1.5
1.5

-1.0
-1.0

-0.5
-0.5

ηeta(m)l(m)分布を図-2.8に示す．図-2.8を見ると，有限体積法の波

400

0.0
0.0

0.5
0.5

1.0
1.0

1.5
1.5

ηl(m)
eta(m)

水位の時間変化（移流項と水平粘性項の追加）（左：有限差分法，右：有限体積法）

形が有限差分法よりも顕著に前傾していることがわかる．
図-2.6と比較して移流項および水平粘性項を追加したこ
1.5の移流および水平粘性項の計算において，cellとfaceにそ1

れぞれ内挿操作を実施していることによる水位差である

0.5

ηl (m)とにより水位分布に違いが表れた．これは，有限体積法と考えられる．

0-0.5

3. 実地形を用いた土砂輸送数値実験

-1
FDM

-1.5

3.1

計算条件

0

200

400

600

FVM

800

x(m)

実海域におけるケーススタディとして，インドネシア
図-2.8のジャワ島北西部(図-3.1(a))にて，現在建設が進められ- 11 -

水位分布（移流項と水平粘性項の追加）

(a)

図-3.1

(b)

計算領域および水深分布 (a)全体領域(構造格子の格子間隔800m)，(b)構造格子中領域(同100m)，(□)構造格子小領域(同30m)，(●)時系列比較位置
表-3.1
海底地形
潮汐
圧力
密度
水平粘性・拡散粗度長
海底摩擦係数
土粒子密度
沈降速度

主な計算条件

した．1wayネスティングは，接続境界において粗格子の水位，流速およびSS濃度の物理量を，細格子の接続境界

GEBCO2014+深浅測量データ
TPXO9-atlasの10分潮（M2,S2,N2,
K2,K1,O1,P1,Q1,Mf,Mm）
静水圧近似
一定（=1,025 kg/m3）
1.0 m2/s
0.4mmに線形補間し受け渡すことで計算を行った．
外力は潮流のみを対象とし，構造格子モデル・非構造
格子モデルとも全球潮汐モデルTPXO9-atlasの主要10分潮の振幅・遅角の情報を開境界に与えた．計算の対象期間
は小硲ら16)の計算期間を内包する2017年8月15日から8月22日とし，前半の5日間は助走期間とした．両モデルと流速の対数分布則
2,650 kg/m3

も，鉛直方向のメッシュを複数に切ることができる多層17)

SS濃度に依存 （最小値：
0.002mm/s，最大値：2.58mm/s）
0.0001 kg/m2/s
巻き上げ係数0.03 Pa
侵食の限界せん断応力
堆積の限界せん断応力∞0.80
底質の空隙率計算が可能なモデルであるが，本検討では単純化のために単層計算とした．圧力の近似についても，SUNTANS
モデルは非静水圧計算が可能なモデルであるが，ここでは計算時間の短縮と構造格子モデルと条件を揃えることを目的として静水圧近似とした．なお，前述したように
非構造格子の計算では直交性の問題が発生するが，SUNTANSでは外心を使用しておりcell centerを結ぶ線と

(a)

(b)

法線方向は一致する．ただし，三角形格子に鈍角が存在する場合は，cell centerが三角形の外に設定されてしまう
ため，SUTNANSでは鈍角が存在する場合にcell centerを外心から重心に近づける操作を行う．この際は誤差が増
大してしまうことに注意が必要である．
土砂輸送シミュレーションは，小硲ら16)と同様に海岸の底質調査結果から泥粒子（粒形5μm）の輸送を対象とし，沈降速度はフロッキュレーションの影響を加味する

図-3.2 港湾施設周辺の格子形状(a)構造格子，(b) 非構造格子

ためにSS濃度に依存して変化させた17)．非構造格子モデ
ルの計算は，構造格子モデルにおける底質の巻き上げ，第1領域の範囲を非構造格子モデルの領域に合わせて図-

沈降，底面せん断応力の評価方法を，SUNTANSにデフ

3.1(a)の範囲とし，第2，第3領域の範囲を図-3.1(b)の桃ォルトで用意されている物質輸送の移流・拡散モジュー

色枠，橙色枠で示すように港周辺の格子を細かくする3

ルに組み込み土砂輸送計算を実施した．地形変化量の評

段階の1-wayネスティングにより計算した(図-3.2(a))．構

価は，底面における巻き上げ量と沈降量の収支として造格子モデルのために構造物との境界で凹凸が生じてし

時々刻々の地盤高変化を計算した．その他の計算条件をまうものの，防砂堤を十分表現できるように，ここでは

表-3.1に示す．

第1領域を800m，第2領域を100m，第3領域を30m格子と- 12 -

3.2

計算結果および考察

以上のような流れ場の違いを基に土砂輸送シミュレー(1) 航路沖における定点モニタリング

ションを実施し，両モデルの地形変化パターンについて

図-3.1(b)の赤点部にモニタリング地点を設け，潮位比較した．図-3.5(a)(b)は，図-3.3の紫色破線で示す上げおよび流速の比較によって非構造格子モデルの妥当性に

潮期間中の，底面における巻き上げ量と沈降量の収支のついて確認した．モニタリング点における潮位の時系列積算値として算定した，基本港形に対する地形変化パタを図-3.3(a)に，東向き成分流速(U)および北向き成分流ーン（空間分布）を示したものである．両モデルとも，

速 (V) の 時 系 列 を 図 -3.3(b) に 示 す ． 潮 位 の 結 果 ( 図 -図-3.4(a)(b)において流速が大きく表れている部分で侵食

3.3(a))を見ると，非構造格子モデルの計算結果は従来のを示しているが，前述したモデル間の流速分布の差違に構造格子モデルとほぼ一致している．従って，構造格子

より，図-3.5(a)(b)中に白枠で示した港北東部，港湾北東

16) に部および港湾南西部において，侵食域が非構造格子モデ

モデル・非構造格子モデル両者の潮位は，小硲らTide level (m)おいて現地観測データと比較しており，両モデルの潮位は妥当であると言える．一方，流速は(図-3.3(b))，東向
き流速，北向き流速ともに流速のピークで非構造格子モ

Velocity (m/s)

デルの方がわずかに小さくなる傾向はあるものの，位相や振幅などの時系列変化の特徴について，概ね同様な結果を示している．
(2) 流れ場の平面分布比較

0.6
0.4
0.2
0
-0.2
-0.4
-0.6
8/19

Structured grid(a)Unstructured grid

2017/8/20
8/20

0.3
0.2
0.1
0
-0.1
-0.2
-0.3
8/19

8/21

8/22

U(Structured grid)U(Unstructured grid)(b)

V(Structured grid)
V(Unstructured grid)

2017/8/20
8/20

8/21

8/22

図-3.3 潮位および鉛直平均流速の時系列図(a)潮位，(b)鉛直平均流速，( )上げ潮期間

前節で示した外力条件を基に港湾施設周辺での計算結果の比較を行った．図-3.3 に橙色破線で示した上げ潮
最強流時刻における，基本港形を対象とした構造格子モ

岬北東部

岬北東部

デル，非構造格子モデルそれぞれによる流速場の計算結果を図-3.4(a)(b)に示す．図-3.3(b)および図-3.4 を見る

岬北部

岬北部

と，当海域の上げ潮時は，港湾北側の航路堀込部を西向
港湾北東

きに横断する流れが発生している．図-3.4(a)(b)を比較

港湾北東

すると，西向きの流れの主流が港湾の沖側を通過しており，航路堀込部や港湾に遮蔽領域では流速が小さくなる
流路分岐位置

流路分岐位置という流況の空間分布の大局的な特徴は両モデルで一致
している．

(a)

岬の北部付近では図-3.4(b)に示すように非構造格子モデルの海陸境界に沿う流れが，構造格子モデルと比較(b)

図-3.4 基本港形における上げ潮時の流速場(矢印：流れの向き，カラーコンター：絶対値)(a)構造格子，(b)非構造格子

して明瞭に示されている．また，岬北東部では，構造格子モデルと比較して非構造格子モデルの流速が大きい．
これは，非構造格子モデルで見られた岬北部での海陸境

(b)

(a)

界に沿う流れによって，岬北東部で要求される流量が大きくなったためと考えられる．一方，港湾北東部において，岬北東部ほど顕著ではないものの，港湾構造物の境界付近で流速に差が見られ，非構造格子モデルの流速の方が大きくなっている．また，同図の青色の矢印で示す
地点は，領域の東側から西方向への流れが港湾施設の東側にぶつかり，流向が港湾の北側に進む流れと，南側に進む流れに分岐する地点である．この位置は，構造格子
モデルと非構造格子モデルで異なっている．その結果，

港湾南西

港湾南西

港湾の南西部を通過する流速が，構造格子モデルよりも
図-3.5

非構造格子モデルの方が大きくなっている．
(3) 地形変化パターンの比較- 13 -

上げ潮時の地形変化パターン(基本港形)(a)構造格子，(b)非構造格子

(a)
防砂堤東側
背後域 防砂堤東側
先端部

防砂堤東側
背後域 防砂堤東側
先端部

(a)

(b)

図-4.3 防砂堤設置に伴う地形変化量の変化（防砂堤 1,000m 設置，上げ潮時）（黒色：基本港形より侵食量減又は堆積量増，白色：基本港形より侵食量増又は堆積量減，
灰色：変化なし） (a)構造格子，(b)非構造格子

図-4.1 防砂堤 1,000m における上げ潮時の流速場(矢印：流れの向き，カラーコンター：絶対値)(a)構造格子，(b)非構造格子

防砂堤東側
防砂堤東側
背後域
先端部

(b)

(a)

(b)

(a)

防砂堤東側
防砂堤東側
背後域
先端部

(b)

図-4.4 防砂堤設置に伴う地形変化量の変化（防砂堤 2,000m 設置，上げ潮時）（黒色：基本港形より侵食量減又は堆積量増，白色：基本港形より侵食量増又は堆積量減，
灰色：変化なし） (a)構造格子，(b)非構造格子

図-4.2 防砂堤 2,000m における上げ潮時の流速場(矢印：流れの向き，カラーコンター：絶対値)(a)構造格子，(b)非構造格子

ルの方が広く分布している．このように，格子形状の違

非計算点にすることで表現した．非構造格子モデルは防いによって港湾施設や岬等の海陸境界付近で流速に違い

砂堤形状に沿って格子を生成することにより，正確に表が見られ，地形変化パターンにも影響を与えることにな

現することが可能である．港湾形状以外の潮位条件やそることが分かった．

の他計算条件は，基本港形ケースと同様とした．

4. 防砂堤設置港形における土砂輸送数値実験

4.2

計算結果および考察(1) 防砂堤設置による流れ場への影響
4.1

計算条件

図-4.1(a)(b)と図-4.2(a)(b)はそれぞれ，図-3.4(a)(b)と東西方向に潮流が往復する当海域において，流れを遮同時刻における，防砂堤を1,000mおよび2,000m設置したる構造物を設置した場合に，格子形状の違いが流れ場おケ ー ス の 流 速 場 で あ る ． 防 砂 堤 1,000m 港 形 ， 防 砂 堤

よび地形変化パターンにどのような影響を与えるのかを

2,000m港形とも，基本港形における岬の北東部や港湾北把握するために，3章で計算を実施した基本港形の沖合

東部，港湾南西部で見られた流速の違いが同様に表れてに向けて1,000mあるいは2,000m防砂堤を設置したケース

いる．基本港形時に航路上を通過していた西向きの流れを設定し，基本港形との比較を行った．構造格子モデル

が，防砂堤を設置したことにより，東側先端部で流速がについては，防砂堤の位置に対応する箇所のメッシュを- 14 -

大きくなる．その一方で，防砂堤の背後域では流れが遮ルとの比較を通した数値波動実験水槽と実海域における

蔽され，流速が小さくなる領域が発生している．防砂堤

潮位計算の比較を行った．

1,000m港形での背後域は，流速としては小さいものの，基礎的な計算手法の整理では，有限体積法を用いた非構造格子モデルの流向(図-4.1(a))が防砂堤に沿って北東

構造格子計算の理論および離散化手法について整理した．

方向に流れているのに対し，非構造格子モデルの流向(図非構造格子を用いた数値シミュレーションでは，座標に-4.1(b))が南西方向と逆向きの流れになっている．また，沿った構造的なナンバリングがされないため，特徴的な流速が小さくなる影響範囲については，カラーコンター

データ構造を必要とすることについてもまとめた．また，が白色で描かれている部分を見て分かるように，構造格波動の伝搬過程計算による構造格子モデルとの比較によ子モデルの方が広く分布している．これらの流れ場の差って，離散化式の取り扱いの違いについても整理した．

異は，格子形状の違いによって，格子の凹凸および防砂

その結果，非構造格子を取り扱う際の直交性の重要性や堤先端の延長の微妙な差異が影響を与えたものと考えら

移流項および水平粘性項の取り扱いの違いにより波形にれる．防砂堤2,000m港形での遮蔽域については，構造格違いが見られることを示した．実海域の計算では，非構子モデルの防砂堤背後域(図-4.2(a))で渦の形成が強調さ

造格子モデルSUNTANSを用いて，インドネシアのパテ

れ，非構造格子モデル(図-4.2(b))と比べて流速が小さく

ィンバン港周辺をケーススタディとし，潮汐流場のみをなっている．

基本とする地形変化計算を試みた．既存の構造格子モデ(2) 防砂堤設置による地形変化への影響

ルによる計算結果との比較を通じて，格子条件の違いに次に，防砂堤を設置したことによる地形変化パターン

よる流速場や地形変化パターンなどの計算結果の差異にへの影響についても調べた．図-4.3(a)(b)，図-4.4(a)(b)

ついて評価した．その結果，非構造格子モデルと構造格は，図-3.5(a)(b)の地形変化パターンを算定した期間と子モデルでは，海岸線に沿う流れや流路の分岐点に違い

同じ期間に，防砂堤設置港形から基本港形の地形変化量が見られ，侵食域の範囲に影響を与えることが分かった．

を差し引いた空間分布である．防砂堤設置港形が基本港

また，防砂堤を設置して潮流を遮ることにより変化する

形よりも堆積側に変化した場合は黒色に，侵食側に変化

流速場は，流速の増大や，背後域での流速の減少の範囲

した場合は白色に描かれている．従って，図-3.5(a)(b)や流向に差異が見られた．これに伴い，防砂堤等の地形において侵食を示していた領域が黒色を示す場合は，侵境界周辺での地形変化量やその空間分布において差異が食量が減少していることを表している．反対に，図-生じることがわかった．これらのことから，地形変化を3.5(a)(b)において堆積を示していた領域が黒色を示す場

予測評価する上では，海陸境界の近傍における流速や地合は，堆積量が増加していることを表している．防砂堤

形変化等のモデル間の相違点に注意することが重要であを設置したことにより，防砂堤 1,000m 港形，防砂堤ると考えられる．

2,000m 港形のケースとも，防砂堤東側先端部の流速が

ただし，本研究の実海域の計算では，主に格子形状の大きくなる領域で，せん断応力の増大により侵食が進む

違いに基づく計算結果の違いを確認したものであり，対傾向となることや，防砂堤背後の流速減速域での堆積傾

象とした海域での流れ場や地形変化等の再現性に関する向など，両計算方法の結果において共通して見られるこ

検証については，現地データの取得を踏まえて，さらなとが分かる．しかしながら，防砂堤 1,000m 港形ケース，

る計算精度の検証を行うことが今後の課題である．

防砂堤 2,000m 港形ケースとも，構造格子モデルと比較（2023年8月3日受付）

して非構造格子モデルは防砂堤の背後のみならず，岬の先端まで堆積傾向を示す黒色の領域がより広く分布する

謝辞

ことなど，異なる格子条件での流速場の結果の差異に依本研究で用いた地形データは，インドネシア国立研究存して，地形変化パターンの計算結果にも影響を与える

革新庁より提供されたものです．また，査読意見を頂いことが示唆された．

た勝谷一則特別研究主幹をはじめ公表審査委員の各位に謝意を表します．

5. おわりに
参考文献
本研究では，非構造格子モデルを用いた土砂輸送計算

1) Chen, C. et al.: An unstructured grid, finite-volume, three-の特徴を把握するために，有限体積法を用いた非構造格

dimensional, primitive equations ocean model: applicatio-

子計算の基礎的な計算手法の整理および，構造格子モデ

n to coastal ocean and estuaries. J.Atm.&Oceanic Tec., 20,- 15 -

159-186, 2003.

staggered mesh schemes. Journal of Computational2) 小松雄二，佐々木淳，松丸Physics, 159(1), 58-89,2000.

亮：非構造格子モデル

15) （ 独 法 ） 国 際 協 力 機 構 ： ENVIRONMENTALを用いたバンダアチェにおける2004年インド洋大津
波の追算，土木学会論文集B2（海岸工学），Vol. 66，

IMPACT ASSESSMENT REPORT OF NEW PORT

No.1，pp. 256-260，2010.

DEVELOPMENT

3) 木梨行宏，山城PROJECTIN

EASTERNMETROPORITAN AREA（PATIMBAN）,p.241,2017.

賢，姫野慎太郎，中野俊夫，横田雅紀，橋本典明：高潮推算における非構造格子モデ

16) 小 硲 大 地 ， 田 村

仁 ， Aloysius Bagyo Widagdo ，ルの適用性に関する基礎的検討，土木学会論文集B3

Dinar Catur Istiyanto，中川康之：海陸風に伴う波（海洋開発），Vol. 68，No.2，pp. I_858-I_863，

浪変動がジャワ島北西部沿岸域の底質移動に及ぼす2012.

影響，土木学会論文集B2（海岸工学），Vol. 77，

4) 増永英治，山崎秀勝，山口邦久：内部潮汐を考慮しNo.2，pp. I_583-I_588，2021.
17) Tsuruya, H., Murakami, K., and Irie, I.:Mathematical

た伊豆大島周辺海域における潮流予測システムの開発と運用，土木学会論文集B2（海岸工学），Vol. 74，

mod-eling of mud transport in ports with a multi-layred

No.2，pp. I_643-I_648，2018.

model -Application to Kumamoto Port-, Report of Port

5) Fringer, O.B., Gerritsen, M. and Street, P.L.: An unstruc-and Har-bour Research Institute, vol.29, No.1. pp.3-55,

tured-grid, finite-volume, nonhydrostatic, parallel coastal

1990.

ocean simulator, Ocean Modeling, vol.14, pp.139-173,
2006.

記号表

6) 中谷祐介，戸松祐希，西田修三：非構造格子モデル

A:検査体積の表面積（m3）を用いた瀬戸内海-太平洋領域における外洋水の挙CV:検査体積の体積（m3）動解析，土木学会論文集B2（海岸工学），Vol. 76，

D:水深（m）

No.2，pp. I_997-I_1002，2020.

d:cell centerとfaceとの距離（m）

7) Zhang, Y., Ye, F., Stanev, E.V., and Grashom, S.: Seam-lesse:接線単位ベクトル

cross-scale modeling with SCHISM, Ocean Modeling,g:重力加速度（m/s2）

vol.102, pp.64-81, 2016.

i:構造格子の場合はx方向の格子番号，非構造格子

8) 山縣史朗，折敷瀬翔耶，山城の場合は計算対象となる格子番号賢，中川康之，橋本j典明，松延嘉國，井方弘正：関門航路のサンドウェ:構造格子の場合はy方向の格子番号，非構造格子ーブに及ぼす潮位偏差の長期変動の影響に関する数の場合は計算対象とfaceを挟んで対となる格子値シミュレーション，土木学会論文集B3（海洋開番号発），Vol. 75，No.2，pp. I_665-I_670，2019.

M:x方向の線流量（m2/s）

9) Deltares : D-Flow Flexible Mesh, Technical ReferenceN:y方向の線流量（m2/s）n:祖度係数（s/m1/3）

10) Chen, C. et al.: An Unstructured Grid, Finite-Volumen:単位法線ベクトル

Community Ocean Model FVCOM User Manual,2013.

R:右辺11) 松田卓也：数値流体力学における差分法,有限体積Sφ:任意の保存量φの生成量法の応用(爆発的天体現象における流体力学的不安T:周期（s）

定性, 研究会報告). 素粒子論研究, 82(1), A42-t

:時間（s）

A47,1990.

u:流速ベクトル（m）

12) H.K.Versteeg,W.Malalasekera,松下洋介,齋藤泰洋,

V:体積（m3）

青木秀之,三浦隆利:数値流体力学 [第2版],森北出

wf:距離の重み関数版, 2011.

Manual, Version: 1.1.0,2021

x,y

:座標系

13) Mathur,S.R. and Murthy,J.Y.:A Pressure-Based Method

Γ

:拡散係数

for Unstructured Meshs, Numer. Heat Transfer Part B,α:振幅（m）

Vol.31, pp.195-215,1997.

ΔA:検査体積の表面要素（m2）Δx:x座標の差（m2）14)

Perot,B.:Conservation

properties

ofunstructured- 16 -

Δy:y座標の差（m2）Δη :node間の距離（m2）Δξ:cell center間の距離（m2）η:face方向に平行な座標系

ηl:水位変動量（m）

ν:動粘性係数（m2/s）

ξ:cell center間に平行な方向の座標系

ρ:流体密度（m3）φ:任意の保存量- 17 -

港湾空港技術研究所資料

No.1410２０２３．９

編集兼発行人

国立研究開発法人海上・港湾・航空技術研究所

発

港 湾 空 港 技 術 研 究 所

行

所

横 須 賀 市 長 瀬 ３ 丁 目 １ 番 １ 号
TEL. 046(844)5040 URL. https://www.pari.go.jp/

C （2023）by MPAT
Copyright ○All rights reserved. No part of this book must be reproduced by any means without the written
permission of the President of MPAT
この資料は、海上・港湾・航空技術研究所理事長の承認を得て刊行したものである。したがって、
本報告書の全部または一部の転載、複写は海上・港湾・航空技術研究所理事長の文書による承認を得ずしてこれを行ってはならない。