27【建設工学研究所論文報告集第 65 号〔論文〕2023 年 11 月】

階層型直交格子の自由表面流れ解析への適応
Application of Hierarchical Cartesian Grids to Free-Surface Flow Analysis
吉 村 英 人Hideto Yoshimura

藤 田 一 郎Ichiro Fujita

１．はじめに河川や氾濫域の流れの解析には 1 次元や平面 2 次元の解析モデルが多く用いられている 1), 2)．これらは広域の流れを効率的に解析することができる一方で，橋脚や堤防など河川構造物周りの局所的な流れや，河川横断面内に生じる二次流などの複雑な 3 次元流れや形状を考慮することは難しい．従来は平面 2 次元モデルでも実用上は問題ないとされてきたが，近年，3 次元流れを考慮する必要性が高まっており，3 次元や準 3 次元モデルの河川・氾濫域の流れへの適応事例も増加しつつある 3), 4)．
河川・氾濫域を対象に 3 次元流れ解析を行う場合に課題となるのが複雑な地形や構造物に対するメッシュ作成と計算にかかるコストである．そのため，これらのコストを低減し，効率的に計算を実施できる方法やプリ・ポスト処理が必要となる．航空やものづくりの分野では複雑形状周りの流れを効率的に計算するため，階層型直交格子 5–9)の活用が進んでいる．階層型直交格子では壁面や任意の領域に対する格子の自動細分化が可能であり，メッシュ作成にかかるコストを大幅に低減することができる．また，ブロック型の階層型直交格子 7–9)では並列計算時の負荷分散を均一にすることができるため，高いスケーラビリティも期待できる．したがって，階層型直交格子は河川・氾濫域の 3 次元流れ解析にも有用であると考えられる．また，現在，国土交通省のプロジェクト PLATEAU10)では 3D 都市モデルの整備と活用が進められているが，このような 3D データとの親和性も高い．
そこで本研究では，河川・氾濫域の 3 次元自由表面流れ解析での活用を目的に，階層型直交格子を自由表面流れに適応する．特に，単相モデルの VOF (Volume of Fluid) 法 11)で自由表面を取り扱った場合に，格子の細分化界面での表面圧力や液相フラックスの補間処理において考慮すべき点について提案する．また，提案手法の検証としてダムブレイク解析を行い，その結果を示す．

２．数値計算手法
本研究では，著者ら 12)が開発した自由表面開水路乱流場を予測するための計算ソルバーを階層型直交格子で利用できるように拡張させた．本節ではそのベースとなる主な計算手法と変更点について説明し，次節で階層型直交格子へ拡張させる際に考慮すべき点について説明する．
２．１

支配方程式および離散化

支配方程式は式(1)および式(2)の液相の非圧縮性流体の連続の式および Navier–Stokes 方程式である．
𝜕𝜕𝜕𝜕�=0𝜕𝜕𝜕𝜕�

𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕� 𝑢𝑢�
1 𝜕𝜕𝜕𝜕
𝜕𝜕+=−+�𝑖(𝜈𝜈 𝜈 𝜈𝜈� )𝑆𝑆�� � + �� , 𝑖𝑖 𝑖𝑖𝑖 𝑖𝑖 𝑖𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕�𝜌𝜌 𝜕𝜕𝜕𝜕� 𝜕𝜕𝜕𝜕�

(1)(2)

28

ここで，𝑡𝑡は時間，𝑥𝑥� は各方向の座標（𝑥𝑥� = 𝑥𝑥，𝑥𝑥� = 𝑦𝑦および𝑥𝑥� = 𝑧𝑧），𝑢𝑢� は流速の各方向成分（𝑢𝑢� = 𝑢𝑢，𝑢𝑢� = 𝑣𝑣および𝑢𝑢� =𝑤𝑤），𝜌𝜌は密度，𝑝𝑝は圧力，𝜈𝜈は動粘性係数，𝜈𝜈� は渦粘性係数，𝑆𝑆�� はひずみ速度テンソル，𝑔𝑔� は重力加速度の各方向成分である．自由表面はVOF法により捕捉し，式(3)に示す液相体積率𝛼𝛼 𝛼 𝛼𝛼𝛼� ⁄Δ𝑉𝑉 に関する移流方程式を解くことで追跡する．
𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕𝜕𝜕�
+=0𝜕𝜕𝜕𝜕�𝜕𝜕𝜕𝜕

(3)

ただし，Δ𝑉𝑉はセル体積，Δ𝑉𝑉� はセル内の液相体積である．

計算領域はスタガード格子を用いた有限差分法で離散化し，流速の各方向成分𝑢𝑢� をセル界面に，圧力𝑝𝑝および液相体積率𝛼𝛼をセル中心に配置する．空間離散化には二次精度中心差分を，時間離散化には移流項に二次精度Adams–Bashforth法，粘性項にCrank–Nicolson法をそれぞれ用いることを基本とする．また，速度と圧力のカップリングはFractionalStep法で行う．

式(1)および式(2)に示すように本計算ソルバーは液相のみを対象とする単相モデルに基づいている．よって，各セルの液相体積率とその空間分布に従い，セルを流体セル，表面セル，空セル，固体セルに分類し，このうち流体セルと表面セルを液相領域として式(1)および式(2)の支配方程式を解く．その際，表面セルでは気液界面で表面圧力が𝑝𝑝� = 0となるように，表面セルの主方向に隣接する流体セルの圧力を内外挿して，境界圧力を設定する．また，表面セル間にあるセル界面やそれに隣接するセル界面では応力ゼロとなるように速度の境界条件を与える．詳細については文献12)を参照されたい．
２．２

Volume of Fluid

VOF法では式(3)の液相体積率𝛼𝛼に関する移流方程式を解く際に，セル界面での液相フラックスが必要となる．液相フラックスは液相体積率𝛼𝛼の空間分布から気液界面を再構築することで決定されるが，再構築手法は代数型と幾何型に大別される．著者ら12)の計算ソルバーではこれまで代数型の一つであるCICSAM (Compressive Interface Schemes forArbitrary Meshes) 13)により再構築を行っていた．しかしながら，代数型では気液界面の位置が陽には定まらないため，後述するように階層型直交格子を用いた際の細分化界面で，液相体積率や液相フラックスの計算を正確に行うことができず，気液界面の拡散を許容することになる．したがって，本研究では再構築手法を幾何型のPLIC (Piecewise Linear
Interpolation Calculation) 14)に改めた．PLICでは気液界面を平面で近似するため，細分化界面においても液相体積率や液相フラックスを正確に求めることができる．
PLICでは各セルにおいて気液界面を，液相体積率𝛼𝛼と気液界面の法線ベクトル𝒏𝒏 𝒏 𝒏𝒏𝒏� , 𝑛𝑛� , 𝑛𝑛� )から以下の式で近似す

る．

𝑛𝑛� ∆𝑥𝑥� + 𝑛𝑛� ∆𝑥𝑥� + 𝑛𝑛� ∆𝑥𝑥� = 𝛽𝛽(4)ここで，∆𝑥𝑥� はセルの各方向サイズ，𝛽𝛽は定数でありセル原点と平面との距離に対応する．なお，法線ベクトルの各方向成分は𝑛𝑛� = − 𝜕𝜕𝜕𝜕 ⁄𝜕𝜕𝜕𝜕� として得られ，𝛽𝛽は界面法線ベクトル𝒏𝒏と液相体積Δ𝑉𝑉� からScardvelli & Zaleski15)の方法により求めた．

３．階層型直交格子
３．１

格子構造

階層型直交格子を用いる計算方法はセル型 5)とブロック型 7–9)に分類できる．セル型はセルごとに細分化を行っていくものであるが，非構造格子となるため高次精度スキームの導入や並列計算時の負荷分散においてデメリットがある．一方，ブロック型では同数のセルを含むブロックを階層的に配置することで計算格子を作成し，ブロック内では構造格子となるため，セル型が持つデメリットは生じない．本研究ではブロック型の階層型直交格子を採用した（図－１）．
一つのブロックには𝑚𝑚� × 𝑚𝑚� × 𝑚𝑚� 個のセルが含まれる．各ブロックは最も大きなブロックであるルートブロックから再帰的に細分化されることにより作成される．ルートブロックの与え方はシングルルート方式とマルチルート方式が考29

図－１ 階層型直交格子のイメージえられるが，本研究では格子生成の柔軟さの観点から，マルチルート方式を採用した．なお，ルート𝑟𝑟のブロックから生成される細分化レベル𝐿𝐿における格子サイズは以下で表される．
(���)

Δ�

(���)

= Δ�

�2�

(5)

3 次元計算の場合，ブロックは通常一回の細分化で各方向 2 分割の合計 8 個に分割される．しかし，河川や氾濫域の計算では鉛直方向よりも水平方向の計算領域が大きくなることも多く，そのような場合には格子数が著しく増大することが予想される．そのため，一回の細分化で水平方向各 2 分割のみの合計 4 個への分割も可能とした．なお本論では，細分化レベルが 1 違う二つのブロックがある場合，細分化レベルが大きい方を粗ブロック，小さい方を密ブロックと呼ぶこととする．

３．２ ブロック間における物理量の交換と補間
階層型直交格子では隣接するブロック間でブロック境界に位置するセル物理量の交換を行う．さらに，ブロックの細分化レベルが異なる場合には物理量の補間が必要となる．交換と補間の処理は隣接ブロックに重なる仮想セルで行われるが，仮想セルの格子幅は隣接ブロックと同じ格子幅，すなわち隣接ブロックが密ブロックの場合は自身の 1/2 倍，粗ブロックの場合は自身の 2 倍で設定する．また，本ソルバーではスタガード格子を採用しているため，各物理量の定義位置（セル中心：圧力，セル界面：流速 3 成分）に応じて補間処理を分ける必要がある．
2 次元のスタガード格子における物理量の交換と補間の処理を，図－２を用いて説明する．なお，以下では細分化レ
ベル𝐿𝐿における物理量を𝑞𝑞 � で統一し，粗ブロックにおけるセル中心の変数定義位置を基準とした密ブロックでの相対的なセルインデックスを𝑞𝑞の下付き添え字で示す．密ブロックから粗ブロックへの補間の場合，セル中心およびセル界面𝑈𝑈���������
�の値は，1 次線形補間𝑞𝑞 � = (𝑞𝑞����+ 𝑞𝑞����)⁄2で与え，セル界面𝑉𝑉では同じ定義位置にある値のみを𝑞𝑞����= 𝑞𝑞����
としてそ���=のままコピーする．粗ブロックから密ブロックへの補間の場合には，セル中心およびセル界面𝑈𝑈の値は 0 次補間で𝑞𝑞�����������
= 𝑞𝑞 � として与え，セル界面𝑉𝑉では同じ定義位置にある値は𝑞𝑞����= 𝑞𝑞����
としてコピーし，粗ブロックに定義点がな𝑞𝑞������+ 𝑞𝑞����)⁄2で与える．
い位置の値は 1 次線形補間𝑞𝑞 ��� = (𝑞𝑞����

３．３ 細分化界面の表面セルにおける圧力境界条件の設定
上述の通り単相モデルに基づく本手法では気液界面で𝑝𝑝� = 0を満たすように表面セルの境界圧力を与えている．重力が支配的な外力の自由表面流れの予測には適切な境界圧力の設定が重要となるため，細分化レベルが異なるブロック間における圧力の補間後も𝑝𝑝� = 0が維持されるようにするべきである．そのため，本研究ではブロック間の細分化界面に表面セルがある場合を図－３に 3 つのパターンに分類し，それぞれにおいて圧力の補間処理を使い分ける．

図－３(a)は水面の主方向が細分化界面と同じ向きであり，表面セルが粗ブロックにある場合である．この場合，境界圧力の設定が必要となるのは密ブロックの仮想セルであり，それらは密ブロックの流体セルから直接与えられる．そのため，粗ブロックから密ブロックへの表面セルの圧力の補間は不要となる．密ブロックの仮想セルにおける液相体積率は粗ブロックの液相体積率から PLIC の処理に従い分配する．図－３(b)は水面の主方向は同じだが，表面セルが子ブロ

30：calc.-cell values:virtual-cell values：calc.-cell values at the adjacent cube

Cell centercalc. cells

calc. cells

virtual cells

virtual cells

Cellinterface U

Cellinterface V

fine to coarse

図－２(a)

FS

coarse to fine

ブロック間における物理量の交換と補間の処理

FSFS

(c)

SSF
coarse-to-fine distributionof pressures(a) パターン 1(b) パターン 2
図－３(c) パターン 3

細分化界面における表面セル圧力の補間パターン

ックにある場合である．この場合，境界圧力は粗ブロックの仮想セルで必要となるため，密ブロックの表面セルに設定された圧力を，粗ブロックの仮想セルに通常の流体セルの圧力と同じ補間処理で与える．図－３(c)は水面の主方向が細分化界面と平行であり，表面セルが粗ブロックおよび密ブロックの両方に含まれる場合である．この場合，粗ブロックの表面セルに設定された境界圧力は，水面の主方向に対する圧力勾配を考慮して，例えば以下のように密ブロックの仮想セルに分配する．
�𝑝𝑝���
� = 𝑝𝑝 ±±�Δ𝑦𝑦 𝑑𝑑𝑑𝑑 �� �
4 𝑑𝑑𝑑𝑑 ��(6)�一方，密ブロックの表面セルには既に𝑝𝑝� = 0を満たすように境界圧力が設定されているため，密ブロックから粗ブロックへの補間は通常と同じ 1 次線形補間を用いれば良い．
３．４

細分化界面における液相フラックスの計算

ブロック間の細分化界面における液相フラックスはドナーセルの液相体積率𝛼𝛼と気液界面の法線ベクトル𝒏𝒏から PLIC���により直接求める．ドナーセルが密ブロックにある場合は密ブロックの各セル界面で液相フラックス𝐹𝐹±���を求め，その������+ 𝐹𝐹����
として粗ブロックのセル界面に与える（図－４(a)）．一方，ドナーセルが粗ブロックにある場合計を𝐹𝐹 � = 𝐹𝐹����合は密ブロックの各セル界面において粗ブロックにあるドナーセルの情報から液相フラックスをそれぞれ求める（図－４(b)）．

31(a) ドナーセルが密ブロックにある場合(b) ドナーセルが粗ブロックにある場合

図－４ 細分化界面における液相フラックスの求め方���(a) スナップショット，𝑡𝑡 ∗ = 𝑡𝑡𝑡𝑡𝑡𝑡⁄𝑡𝑡𝑡(b) 水柱先端位置の時間変化

図－５ ダムブレイク解析の結果

４．計算結果2 次元のダムブレイク解析 16, 17)により提案手法の検証を行った．水槽を模擬した計算領域4𝐿𝐿 𝐿 𝐿𝐿𝐿（ただし，代表長さ𝐿𝐿𝐿 𝐿𝐿𝐿𝐿𝐿 m）の左端に𝐿𝐿 𝐿 𝐿𝐿𝐿の水柱を初期条件として与え，これが重力により崩壊する流れを解析する．計算条件として，重力加速度𝑔𝑔𝑔 𝑔𝑔𝑔𝑔 m/s2，動粘性係数𝜈𝜈 𝜈𝜈𝜈𝜈𝜈𝜈𝜈�� m2/s とし，壁面境界はすべてスリップ条件を与えた．解析は⁄ ，右側壁面でΔ� = Δ� = 𝐿𝐿⁄𝐿𝐿となるように細Δ� = Δ� = 𝐿𝐿⁄𝐿𝐿の等間隔格子を用いたケースと，左側壁面でΔ� = Δ� = 𝐿𝐿𝐿分化した階層型直交格子を用いたケースで行い，得られた結果を比較した．図－５(a)に階層型直交格子を用いた場合の水柱崩壊時のスナップショットを示す．本手法により細分化界面においても液相体積率の分布が拡散せず，輸送されていることが確認できる．図－５(b)に水柱先端位置𝑥𝑥� の時間変化を実験結果 16, 17)と併せて示す．階層型直交格子の結果は等間隔格子の結果と概ね一致しており，SOLA-VOF11)とも一致する結果が得られている．

５．おわりに
本研究では複雑な地形や構造物を有する河川・氾濫域における 3 次元自由表面流れ解析を効率的に行うことを目的に，階層型直交格子を用いた単相 VOF 法を提案した．特に，細分化界面における自由表面流れを適切に予測するため，表面圧力や液相フラックスの補間処理において考慮すべき点について示した．2 次元ダムブレイク解析による検証計算では階層型直交格子を用いた場合でも，気液界面が拡散することなく輸送できることを示した．なお，本論では 2 次元の場合で説明したが，ここで示した補間処理はブロック界面で隣接するセルの物理量だけを用いて行われるため，3 次元の場合は一つの粗セルに対して四つの密セルを考慮することで拡張できる．
本研究で検証した流れは静水圧が支配的な流れであるため，今後は非静水圧効果が強い流れに対しても検証を行う必要がある．また，乱流モデルを含めたうえで本手法を実スケールの 3 次元流れ場にも適応していく予定である．

32

参考 文

献

Horritt, M. S. and Bates, P. D.: Evaluation of 1D and 2D numerical models for predicting river flood1)

inundation, Journal of Hydrology, Vol. 268, pp. 87–99, 2002.
Costabile, P., Costanzo, C. and Macchione, F.: Performances and limitations of the diffusive approximation2)of the 2-d shallow water equations for flood simulation in urban and rural areas, Appl. Numer. Math. Vol.
116, 141–156, 2017.
3)

竹村吉晴，福岡捷二：波状跳水・完全跳水及びその減勢区間における境界面（水面・底面）上の流れの方程式を用いた非静水圧準三次元解析（Q3D-FEBS），土木学会論文集B1（水工学），Vol. 75, No. 1, pp. 61–80, 2019．

4)

柏田仁，二瓶泰雄：国内外の河川流・氾濫流解析モデルのレビュー ～3次元モデルに着目して～，河川技術論文集，Vol. 29, pp. 203–208, 2023.
Hartmann, D., Meinke, M. and Schröder, W.: An adaptive multilevel multi grid formulation for Cartesian5)hierarchical grid methods, Computers & Fluids, Vol. 37, No. 9, pp. 1103–25, 2008.
Kuya, Y. and Kawai, S.: A stable and non-dissipative kinetic energy and entropy preserving (KEEP) scheme

6)

for non-conforming block boundaries on cartesian grids, Computers & Fluids 200, 104427, 2020.
Ghomizad, M. B., Kor, H. and Fukagata, K.: A structured adaptive mesh refinement strategy with a sharp

7)

interface direct-forcing immersed boundary method for moving boundary problems, J. Fluid Sci. and Tech.,Vol. 16, No. 2, p. JFST0014, 2021.
Nakahashi, K.: Building-cube method for flow problems with broadband characteristic length, Comput.

8)

Fluid Dyn. 2002, Springer Berlin Heidelberg, pp. 77–81, 2003.
Jansson, N., Bale, R., Onishi, K. and Tsubokura, M.: "CUBE": a scalable framework for large-scale industrial

9)

simulations, Int. J. High Perform. Comput. Appl., Vol. 33, No. 4, pp. 678–698, 2019.
10)

PLATEAU：https://www.mlit.go.jp/plateau/

11)

Hirt, C. W. and Nichols, B. D.: Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput.

Phys., Vol. 39, pp. 201–225, 1981.
12)

Yoshimura, H. and Fujita, I.: Investigation of free-surface dynamics in an open channel flow, J. Hydraul.

Res., Vol. 58, No. 2, pp. 231–247, 2020.
13)

Youngs, D. L: Time-dependent multi-material flow with large fluid distortion, in Morton, K. W. and Baines,
M. J.

14)

(Eds.), Numerical Methods for Fluid Dynamics, Academic Press, New York, pp. 273–285, 1982.

Ubbink, O. and Issa, R. I.: A method for capturing sharp fluid interfaces on arbitrary meshes, J. Comp. Phys.,Vol. 153, pp. 26–50, 1999.

15)

Scardvelli, R. and Zaleski, S.: Analytical relations connecting linear Interface and volume fractions inrectangular grids, J. Comp. Phys., Vol. 164, pp. 228–237, 2000.

16)

Martin, J.C. and Moyce, W. J.: An experimental study of the collapse of liquid columns on a rigid horizontal
plane, Philos. Trans. R. Soc. A, Vol. 244, pp. 312–324, 1952.

17)

Koshizuka, S., Tamano, H. and Oka., Y.: A particle method for incompressible viscous flow with fluidfragmentation, Comput. Fluid Dyn. J., Vol. 4, pp. 29–46, 1995.

著者

吉村英人

研究員，博士（工学）
，河川工学

藤田一郎

所員，学術博士，河川工学