Centred Scheme を用いた CIVA- 格子ボルツマン法 による浅水長波解析手法の構築

修士論文要旨 （2008年度）

Centred Scheme ^を用いた CIVA- ^{格子ボルツマン法} による浅水長波解析手法の構築

Development of CIVA-Lattice Boltzmann Method for Shallow Water Flow based on Centred Scheme

土木工学専攻 18号 楠 和也 Kazuya KUSUNOKI

1.

海岸，河川などにおいて生じる津波や高潮などの多く の流れは浅水長波方程式で記述される．これらの流れの 定量的な把握は，水害による被害の大きさを予測する 上で重要であり，これまで主に模型実験による方法が 一般的に用いられてきた．しかし，近年の計算機性能の 向上に伴い，数値シミュレーションも数多く行われる ようになってきている．これら浅水長波流れの解析手 法としては，有限差分法や有限要素法などが主流であ るが，近年新たな数値解析手法として格子ボルツマン 法(LBM:Lattice Boltzmann Method)^1),4)が注目を集 めている．LBMは流体を有限個の速度をもつ仮想的な 粒子の集合体としてとらえ，各粒子の衝突と並進を逐次 計算し，粒子の計算を局所的に行う陽的な解法であり，

かつアルゴリズムも簡便な手法であることから，計算が 高速に処理可能である．従って，LBMはこれらの浅水 長波流れ問題に対して有効な手法になり得るといえる．

しかしながら，LBMでは構造格子を用いるのが一般的 なので，自然地形を取り扱うことが多い浅水長波流れ 解析においては，その適用に難点があった．そこで，立 石，樫山により非構造格子へ拡張されたCIVA-格子ボル ツマン法⁸⁾が提案され，浅水長波流れ解析に適用されて きた．

しかし，外力項の評価において従来適用していたBasic

Schemeでは底面形状に不連続な点があると解析が不安

定になるという問題が生じていた．そこで本研究では，

上記問題点を改善するため，全水深と河床勾配を並進ベ クトルの中点の値で評価するCentered Schemeの導入 を試みた．数値解析例として様々な底面形状を有する潮 流解析を取り上げCentred Schemeの有効性の検討を 行った後に，Centred SchemeのCIVA-LBMへの適用 について検討する．

2. CIVA-

格子ボルツマン法において，衝突演算項に格子BGK モデルを用いた格子ボルツマン方程式は以下のように表 わされる．

fα(x+eα∆t, t+ ∆t)−fα(x, t)

=−1 τ h

fα(x, t)−f_α^eq(x, t)i + ∆t

6e²eαiFi (1)

図– 1 ^概念図

上式において，左辺は粒子の並進過程，右辺の第1項目 は衝突過程を表している．並進過程にCIVA法を導入 することによってLBMの非構造格子への拡張を図って いる．fα はα方向の粒子がどれくらい存在するかとい うことを表す粒子分布関数，f_α^eqは局所平衡分布関数で ある．また，eαは粒子の並進速度ベクトルであり，以下 の式によって表される．

e_α=







[0,0], α= 1

e[cos^(α−₄^1)π,sin^(α−₄^1)π], α= 2˜5 p2e[cos^(α−₄^1)π,sin^(α−₄^1)π], α= 6˜9

(2)

式(1)，(2)のeはe= ^∆x_∆t で定義され，∆tは微小時間 増分量，∆xは格子サイズである．式(1)の支配方程式 により粒子の計算を局所的に行い，陽的に未知量f_αを 求める．なお，本論文では2次元9速度モデル（図-1） を用いており，粒子分布関数などの添え字αは図の数字 と対応している．また，F_iは外力項であり，以下のよう に与えられる．

F_i=−gh∂z_b

∂xi

+τ_wi ρ −τ_bi

ρ +E_i (3) 上式において，ρは流体密度，z_bは基準面からの底面の 高さ（zbの定義は図-2参照），τbiは底面でのせん断応 力，τwiは自由表面でのせん断応力，Eiはコリオリ力で ある．τbiおよびτwiは以下の式によって算出される．

τbi =ρCbui

q

ujuj (4)

τ_wi=ρ_aC_wu_wi q

u_wju_wj (5) ここで，式(5)において，ρ_aは空気の密度，C_wは空気 抵抗，uwは風速を表す．また，Cb は底面摩擦であり，

図– 2 ^{水深定義図}

Chezyの係数（Cz=h¹⁶/n）を用いて以下の式によって 決定される．なおnはマニングの粗度係数である．

Cb= g

C_z² (6)

また式(1)において，τは単一時間緩和係数と呼ばれる 定数であり，1タイムステップの衝突で粒子が格子点上 において一定の割合 _τ¹ で局所的な平衡状態に近づいて いくことを示している．τ は，鉛直方向に積分された渦 動粘性係数νeと以下のような関係にある．

τ= 3νe

e²∆t +1

2 (7)

なお，式(3)第1項目の底面勾配の評価方法については，

Basic SchemeとCentred Schemeの2通りについて検 討した．ここで，Basic Schemeとは，図-1(a)のように 9方向すべてで節点での全水深と底面勾配を用いて次式 により評価する手法である．

−gh(x, t) ∂

∂xi

zb(x, t)  (8) 一方，Centred Schemeは，Basic Schemeが節点値で評 価するのに対し，図-1(b)のように9方向それぞれの並 進ベクトルの中点での全水深と底面勾配を用いて次式に より評価する手法である．

−gh(x+1

2eα∆t, t) ∂

∂xi

zb(x+1

2eα∆t, t) (9) 2.2 局所平衡分布関数

局所平衡分布関数は，ある空間内において平衡状態に なった場合の粒子の分布を表す関数であり，巨視的変数 である全水深と流速によって決定され，以下の式で表さ れる．

f_α^eq=















h−^5gh6e²² −3e^2h2uiui α= 1

gh²

6e² +_3e^h2eαiui

+_2e^h4eαieαjuiuj−6e^h2uiui α= 2˜5

gh²

24e² +_12e^h2eαiui

+_8e^h₄e_αie_αju_iu_j−24e^h2u_iu_i α=6˜9 (10)

ここに，hは全水深，u_iは流速，gは重力加速度である．

なお，上式の係数は，巨視的速度においてべき級数をと り，質量保存，運動量保存，エネルギー保存を満たすよ うに決定されている．

2.3 流れの巨視的変数

流体の巨視的変数である全水深及び速度は，その粒子 分布関数と式(2)の粒子の並進速度ベクトルを用いて以 下のように計算される．

h= X9

α

f_α (11)

u= 1 h

X9 α

eαfα (12)

2.4 境界条件

本 論 文 で は ，non-slip 条 件 に 相 当 す る も の と し て Bounce-back条件²⁾を，slip 条件に相当するものとし てMirror 条件²⁾，流入流出境界条件として粒子分布関 数の0勾配条件を用いた．

2.4.1 Bounce-back条件

Bounce-back条件は，粒子が来た方向に180^◦跳ね返 るという条件であって，non-slip条件としてはもっとも よく用いられており，多孔質媒体内流れなどで成功を収 めている．図-3(a)のように下面が境界の場合，未知粒 子分布関数は以下のように表される．

f3=f5,

f₆=f₈, (13)

f7=f9

2.4.2 Mirror条件

slip条件の場合は一般的にMirror条件がよく用いら

れる．Mirror条件とは，壁に向かう方向の粒子が鏡面

反射する条件であり，図-3(b)のように下面が境界の場 合，以下のように表される．

f₃=f₅,

f₆=f₉, (14)

f₇=f₈ 2.4.3 粒子分布関数の0勾配条件

流入・流出境界条件としては，境界に直行する方向の 粒子分布関数の勾配を0にする境界条件処理⁴⁾ が提案 されている．例えば，図-3(c)のように左から右へ流体 が流入する場合を考えるならば，以下の式により算出さ れる．

f₂(0) =f₂(1)

f₆(0) =f₆(1) (15) f9(0) =f9(1)

上式において，0は流入部，1はその風下方向の隣接す る格子点である．流出に関しても同様に処理する．

図– 3 ^概念図

(x-u t,y-v

)

) , (x-uDty-vDt f

) 1 , 1 ( ) 1 , 1

1(

t v y t u x y

x ⁿ

n⁺ =f - D - D

f

) 2 , 2

2(x y P ) , (x y P

x

y

) 1 , 1

1(x y P

) 3 , 3

3(x y P

ƒÓ

図– 4 CIVA法

2.5 CIVA法

CIVA法⁷⁾は，移流方程式の高精度かつ安定な解法で あるCIP法⁶⁾を三角形及び四面体に拡張した手法であ る．CIP法やCIVA法では，各節点の局所厳密解を上 流側の格子内に張った3次の補間曲面により求める．空 間内の完全3次多項式を得るためには10個の係数を求 める必要があるため，各節点上に関数値及び空間微係数 値を配置する．（図-4）

 四角形格子の場合，4点計12変数の情報から3次多 項式の未知係数が決定されるが，三角形格子の場合は3 点9変数しかなく，3次多項式を完全に求めることがで きない．CIVA法では，この問題を3次関数を調整し未 知係数を減らすことにより解決している．CIVA法で用 いられる補間関数及び微係数の補間関数は，三角形に対 して正規化された座標系である面積座標系を用いて以下 のように表される．

f˜(L1, L2, L3) = X3 i=1

αiLi+d X3

j,k=1 j̸=k

βjk[L²_jLk+cL1L2L3]

(16) なお，式(16)のdは3次補間の場合1となり，cは0.5 である．

 CIVA法のように三角形格子上の補間により移流方 程式を解く場合には，上流側の点がどの要素に存在する かを探索する必要がある．本論文では，CFL数が1を

超えないように∆xを設定し，上流側の点を隣接する要 素内から探索した．解析領域中に大きさが異なる要素が 混在している場合は，最も小さい要素でCFL数1以下 の条件を満たさなければならない．そのため，計算に用 いる∆xは最も小さい要素を基準にして設定される．

3.

数値解析例1及び数値解析例2において，底面形状が 滑らかで，その傾きが解析的に求めることのできる形状 を有する潮流解析及び，底面形状に不連続な点を有する 潮流解析を取り上げることでCentred Schemeの有効性 の検討を行い，数値解析例3においてCentred Scheme のCIVA-LBMへの適用について検討する．

3.1 潮流解析1

数値解析例として滑らかな底面形状を有する潮流解析 を取り上げる．解析モデルを図-5に示す．

(a)水深図 (b)流速図

図– 5 9117.5^{秒後の解析結果}

初期条件として水深は以下の式により与えることと する．

H(x) = 50.5−40x

L −10sin(π(4x L −1

2)) (17) 上式において，水路長はL= 14000mである．水面はフ ラットであるから，標高値はz_b(x) =H(0)−H(x)で 表される．また，境界条件は以下の式で与える．

h(0, t) =H(0) + 4−4sin¡ π( 4t

86400+1 2)¢

(18)

u(L, t) = 0 (19)

また，u(0, t)においては粒子分布関数0勾配⁴⁾を適用し た．分割数はx方向に1000分割とし，単一時間緩和係 数τ は0.6とした．なお，底面勾配の評価法としては Basic Schemeを用いて計算した．

計算結果として図-5に9117.5秒後の(a)水位図及び，

(b)流速図を示す．図-5より，水位図，流速図ともに解 析解と定量的な一致が得られていることがわかる．底面 勾配が滑らかで解析的に求められるため，Basic Scheme でも安定に解析が行えていることが確認できる．

図– 6 10800^{秒後における水面形}

図– 7 10800秒後における流速図

3.2 潮流解析2

数値解析例として，図-6のように底面が不規則な形状 をしている場合における潮流解析を取り上げた．格子幅

∆x=7.5[m]，微小時間増分量∆t=0.3[s]，単一時間緩和 係数τ=1.5として計算した．水路長はL=1500[m]であ る．以下に境界条件を示す．

h(x,0) =H(x), H(0) = 16, H(x) =H(0)−z_b(x)

u(x,0) = 0 (20)

h(0, t) =H(0) + 4−4sin£ π¡ 4t

86400+1 2

¢¤ (21)

解析結果として，10800秒後の水面形状の図を図-6に示 す．図-6より，Basic Schemeにおいては計算が不安定に なり，水面が振動してしまっているが，Centred Scheme はその振動が抑制されていることが確認できる．また，

10800秒後の流速の比較図を図-7に示す．図-7により，

Basic Schemeでは底面が不連続な点において流速が大 きく振動してしまっているが，Centred Schemeを導入 することにより振動は抑えられ，解析解と良い一致を示 していることがわかる．

3.3 潮流解析3

数値解析例として，図-8と同様に底面が不規則な形状 をしている場合における潮流解析を取り上げ，Centred SchemeのCIVA-LBMへの適用について検討する．

格 子 幅 ∆x=7.5[m]，微 小 時 間 増 分 量 ∆t=0.1[s]，単 一 時 間 緩 和 係 数 τ=1.6 と し て 計 算 し た ．水 路 長 は L=1500[m]である．以下に境界条件を示す．

h(x,0) =H(x), H(0) = 16, H(x) =H(0)−zb(x)

u(x,0) = 0 (22)

また，

h(0, t) =H(0) + 4−4sin£ π¡ 4t

86400+1 2

¢¤ (23)

u(L, t) = 0 (24)

(a)水深図 (b)流速図

図– 8 10800^{秒後における解析結果}

解析結果として，10800秒後の水面形状の比較図と流速 の比較図を図-8に示す．図-8(a)より，Centred Scheme を適用したLBM及びCIVA-LBMにおいては，水面 が振動することなく安定に解析できているが，Basic Schemeを用いたCIVA-LBMにおいては計算が不安定 になり，水面が振動してしまっていることがわかる．ま た，図-8(b)より，Centred Schemeを導入することによ りLBM及びCIVA-LBMどちらの手法による結果も流 速の振動は抑えられ，解析解と良い一致を示しているこ とがわかる

4.

本論文では，数値解析例として潮流解析とを取り上 げ，Centred Schemeを適用したCIVA-LBMの有効性 について検討した結果，以下の結論を得た．

• 底面の形状が複雑で，不連続な点を有する場合，

Basic Schemeでは解析が不安定になるが，Cen- tred Schemeでは安定に解析可能であることが確 認できた．

• Centred SchemeをCIVA-LBMへ適用し，その 有効性が確認された．

今後の課題としては，本手法の実地形に対する問題への 適用があげられる．

参考文献

1) Chen,S., and Doolen,G.D.：Lattice Boltzmann Method for Fluid Flow, Annu. Rev. Fluid Mech.,30, pp.329- 364, 1998.

2) S.Succi: The Lattice Boltzmann Equation for Fluid Dynamics and Beyond,Oxford University Press, 2001.

3) 稲 室 隆 二 ：格 子 ボ ル ツ マ ン 法，物 性 研 究，pp.197- 232, 2001.

4) J.G.Zhou:Lattice Boltzmann Methods for Shallow Wa- ter Flows, Springer, 2003.

5) Y.H.QIAN: Lattice BGK Models for Navier-Stokes Equation.Europhysics letters, 1992.

6) T.Yabe:A universal solver for hyperbolic equations by cubic-polynomial interpolation, Comput.Phys. Com- mun.,66, pp219-232, 1991.

7) 田中伸厚:数値流体力学のための高精度メッシュフリー手法 の開発，機論，64-620B, pp1071-1078, 1998.

8) 立石絢也，樫山和男：CIVA-格子ボルツマン法による非構 造格子を用いた非圧縮性粘性流体の解析，応用力学論文集，

土木学会，Vol.7, pp323-329,2004.