土木工学専攻

(1)

2006

浅水長波流れ解析のための CIVA- 格子ボルツマン法の構築研究

Studies on Development of CIVA - Lattice Boltzmann Method for Shallow Water Flow Analysis

土木工学専攻

3

号石川裕士

ISHIKAWA Yuji

1.

研究目的

海岸，河川における水害による被害や水辺環境の変化の予測を行う上で数値シミュレーションは有効な手段となっている．これらの予測数値シミュレーションは，一般に大規模計算となることが多いため，高精度かつ高速な数値解析手法の構築は重要である．これまで，浅水長波流れの解析手法としては，有限差分法，有限体積法及び有限要素法などが主流であるが，近年新しい数値解析手法として格子ボルツマン法

^{1) 2)}

が注目されている．格子ボルツマン法は，アルゴリズムが簡便なため，計算が高速に処理可能であり，大規模計算に対して有効な手法であるといえる．

しかし，格子ボルツマン法は構造格子に基づく手法であるため，自然地形を取り扱うことが多い浅水長波流れ解析においては，その適用性に難点がある．これまで，

非構造格子へ拡張する試みとして，

Navier-Stokes

流れに対して，非構造格子における移流方程式に対する高精度スキームである

CIVA

法

³⁾

を導入した

CIVA-

格子ボルツマン法が立石らにより提案された

⁴⁾

．

そこで本論文は，

CIVA-

格子ボルツマン法に着目し，

それによる非構造格子に基づく浅水長波流れ解析手法を

提案する

^{5) 6)}

と共に，水深変化を有する場合の精度，開

境界条件処理及び並列化について検討を行うものである．本手法は，本来の格子ボルツマン法の特徴である計算の局所性を有しており，陽的な計算となるため，計算機容量や並列計算の点で優れた解法であると期待できる．本手法を突起のある開水路流れ問題及び東京湾潮流解析に適用し，本手法の解析精度及び並列化効率の観点から，本手法の妥当性と有効性を検討した．

2. CIVA-

格子ボルツマン法

^{1) 2)}

は，流体を有限個の速度をもつ多数の仮想粒子の集合体で近似し，各粒子の衝突と並進とを計算する手法である．流体の巨視的変数は，求められた粒子分布関数を用いて質量保存・運動量保存により算出される．

2.1

離散ボルツマン方程式

粒子分布関数の時間発展は，以下の離散ボルツマン方程式で表される．

∂fα

∂t +eα

∂fα

∂x =−1 τ h

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

6e²eαiFi

(1)

上式において，左辺は粒子の並進過程を表しており，

CIVA-

格子ボルツマン法では粒子分布関数の並進計算に

CIVA

法を導入し，非構造格子への計算を可能としている．また，右辺の

1

項目は格子

BGK

モデルに基づく衝突演算項，

2

項目は力に関する項をそれぞれ示している．

fα

は

α

方向の粒子がどのくらい存在するかということを表す粒子分布関数，

f_α^eq

は局所平衡分布関数である．

また，

eα

は図

-1

に示す粒子の並進速度ベクトルであり，

本論文で用いた

2

次元

9

速度モデルでは以下の式によって表される．

e_α=







[0,0], α= 1

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

(2)

1 2

3

5 4

6 7

9 8

図

– 1 2

次元

9

速度モデル

式

(1)

，

(2)

の

e

は

e= ^∆x_∆t

で定義され，

∆t

は微小時間増分量，

∆x

は格子サイズである．式

(1)

の支配方程式により粒子の計算を局所的に行い，陽的に未知量

f_α

を求める．また，

F_i

は力に関する項であり，以下のように与えられる．

F_i =−gh∂z_b

∂xi

+τ_wi ρ −τ_bi

ρ +E_i (3)

上式において，

ρ

は流体密度，

z_b

は河床高度，

τ_bi

は底面でのせん断応力，

τ_wi

は自由表面でのせん断応力，

E_i

はコリオリ力である．詳細は参考文献

²⁾

を参照されたい．

なお，式

(3)

における

z_b

の定義は図

-2

のようになっている．

また，式

(1)

において，

τ

は単一時間緩和係数と呼ばれる定数であり，

1

タイムステップの衝突で粒子が格子点上において一定の割合

¹_τ

で局所的な平衡状態に近づいていくことを示している．

τ

は，鉛直方向に積分された渦動粘性係数

νe

と以下のような関係にある．

τ = 3νe

e²∆t +1

2 (4)

(2)

図

– 2

水深定義図

2.2

局所平衡分布関数

f_α^eq

とは，格子点上において流体が平衡状態に達したときの粒子分布であり，巨視的変数である全水深及び速度を用いて以下のように求められる．

f_α^eq=











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

gh²

6e² +_3e^h2eαiui

+_2e^h₄e_αie_αju_iu_j−6e^h²u_iu_i α= 2,· · ·,5

gh²

24e² +_12e^h₂e_αiu_i

+_8e^h4eαieαjuiuj−24e^h²uiui α= 6,· · ·,9 (5)

上式において，

h

は全水深，

ui

は流速，

g

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

²⁾

．

2.3

流れの巨視的変数

流体の巨視的変数である全水深及び速度は，その粒子分布関数と式

(2)

の粒子の並進速度ベクトルを用いて以下のように計算される．

h= X9

α

fα, u= 1 h

X9 α

eαfα (6)

2.4

力に関する項の取り扱い

前述したように，格子ボルツマン方程式における力に関する項は以下のようになる．

Fi=−gh∂zb

∂xi

+τwi

ρ −τbi

ρ +Ei (7)

上式において，第

1

項目の河床勾配の評価は，まず三角形格子毎の河床勾配を算出し，それを用いて，最小二乗法により各格子点における河床勾配を算出している．なお，この値は，固定床の場合には毎ステップ更新をする必要はない．一方，第

2

項から第

4

項はいずれも流速に関わる項なので，毎ステップ更新する必要がある．

2.5

並進過程の計算（

CIVA

法）

本研究では，高精度で分散誤差の少ない

CIVA

法

³⁾

を並進過程に用いて計算する．

CIVA

法は，移流方程式の高精度かつ安定な解法である

CIP

法を三角形及び四面体に拡張した手法である．

CIP

法や

CIVA

法では，各節点の局所厳密解を上流側の格子内に張った

3

次の補間曲面により求める．

2.6

境界条件処理

計算過程において，衝突及び並進計算が終了した際に境界外から計算領域に入る方向の粒子分布関数

fα

を決定する必要がある．すなわち，境界上の領域に向かう法線ベクトルを

nα

とすると，

eα·nα >0

を満たす

fα

を算出しなければならない．例えば，図

-1

において下部に壁面があるとするならば，

f3

，

f6

，

f7

を算出しなければならないことになる．

non-slip

境界条件処理としては

Bounce-back

条件，

slip

境界条件処理としては

Mirror

条件，流入流出境界条件としては，粒子分布関数の

0

勾配条件，平衡分布などが挙げられる．以下にそれらの計算方法を示す．

2.6.1 Bounce-back

条件

non-slip

境界条件処理としては，

Bounce-back

条件が広く用いられている．この境界条件処理は

,

粒子が来た方向に

180^◦

跳ね返るという条件である．以下の式によって，

non-slip

境界条件を満たすことになる．

（図

-3(a)

）

f₃=f₅, f₆=f₈, f₇=f₉ (8)

2.6.2 Mirror

条件

slip

境界条件処理としては，

Mirror

条件が広く用いられている．この境界条件処理は，壁に向かう方向の粒子が鏡面反射する条件であり，以下のように表される．

（図

-3(b)

）

f3=f5, f6=f9, f7=f8 (9)

2.6.3

粒子分布関数の

0

勾配条件

流入・流出境界条件としては，解析領域内に向かう方向の粒子分布関数の勾配を

0

にする境界条件処理

²⁾

が提案されている．例えば，図

-1

において左から右へ流体が流入する場合を考えるならば，以下の式により算出される．（図

-3(c)

）

f₂ⁱⁿ=f₂^d, f₆ⁱⁿ=f₆^d, f₉ⁱⁿ=f₉^d (10)

上式において，

f_iⁱⁿ

は流入部，

f_i^d

はその風下方向の隣接する格子点である．流出に関しても同様に処理する．

2.6.4

平衡分布条件

流入・流出境界条件としては，未知粒子分布関数を平衡分布にする条件がある．この境界条件は，巨視的変数である水深及び流速から算出するため，より質量保存及び運動量保存を満たすようになる．

2.7

並列計算法

本研究では，分散メモリ型並列計算機を対象とした領域分割法に基づく並列計算手法を導入する．

計算領域を，グラフ理論に基づく自動領域分割システ

ムである

Metis

を用いて分割し，分散メモリ型並列計

算機の各ノードに割り当てる．通常の有限要素法の並列

計算においては，要素の重ね合わせを行うための相互の

(3)

ო㕙 ო㕙

wall wall

5

9 8 7 3 6

ო㕙 ო㕙

wall wall

5

9 8 7 3 6

ო㕙 ო㕙

wall wall

ო㕙 ო㕙

wall wall

5

9 8 7 3 6

(a)Bounce-back scheme

㪮㪸㫃㫃㪮㪸㫃㫃

in d 6 2 9 in d

6 2 9

(b)Mirror scheme

(c)Zero gradient

図

– 3

境界条件処理

通信が必要になるが，

CIVA

法を用いた場合の並列計算では，上流側で求めた補間値を下流側へ通信する一方向の通信のみで計算が行える．ここでは，一例として図

-1

の

2

番目の方向の場合について説明する．図

-4

の左側の領域が上流側とすると，下流側への通信のみとなり，

右側領域から左側領域への通信はない．しかし，左側領域から右側領域への通信節点数と右側領域から左側領域への通信節点数が移流方向によって決まるため，データの作成が複雑になる．格子ボルツマン法の並進過程では，移流方向は格子ベクトルの

8

方向に固定されているため，計算の前に通信データを作成することができ，

タイムステップの間に通信データを変更する必要はない．なお，通信ライブラリには

MPI

（

Message Passing Interface

）を使用した．また，並列計算機として，

20

台の

PC

（

CPU

：

Xeon3.06GHz

，メモリ：

2GB

，キャッシュ：

512kB

）を

Gigabit Ethernet

により結合した

PC

クラスタ型並列計算機を用いた．

図

– 4

並列通信

3.

数値解析

3.1

突起のある開水路定常流れ

水深変化を有する場合の解析精度を検討するために，

突起のある開水路流れを取りあげる．図

-5

に解析モデルを示す．

図

-5

の左側から単位幅あたりの流量

Q= 4.42[m²/s]

を流入させ，右側の流出側では，

h= 2.0[m]

を課すものとし，側面は

slip

条件とする．なお，渦動粘性係数

ν_e

は，

ν_e= 0.01[m²/sec]

を与えた．

水深変化を有する場合において，解像度による解析精度と厳密解の差である誤差の収束性について検討するために，構造格子を用い，格子サイズが

0.01

，

0.05

，

0.1

，

0.2

の場合について，解像度による収束性について検討した．

解析結果として，図

-6

（

a

）に，各格子サイズにおける定常時での水面形状を示す．この図より，解像度を上げることにより，より厳密解に近づいてるいることが確認できる．次に，図

-6

（

b

）に解像度と誤差の関係を示す．

図より，解像度と誤差にはほぼ

1

次精度であることがわかる．

X Z Y

25.0 m

5.0 m 2.0 m

0.2 m

Q

図

– 5

解析モデル

0 10 20

1.8 1.9 2 2.1

x [m]

WaterElevation[m]

Size=0.2 Size=0.1 Size=0.05 Size=0.01 Analytical Solution

(a)

定常時における水面形状

10^-3 10^-2 10^-1 10⁰

Mesh size [m]

Error

1 1

10^-3 10^-2 10^-1 10⁰

Mesh size [m]

Error

1 1

(b)

解像度と誤差の関係図

– 6

格子サイズの差異による検討

3.2

東京湾潮流解析

実際の潮流現象の解析として，複雑な境界形状と水深変化を有する東京湾潮流解析を取りあげた．なお，開境界処理に平衡分布と粒子分布関数の

0

勾配条件を適用し，比較検討を行う．有限要素分割図は図

-7

（

a

）に示すとおりであり，格子数は

5839

，要素数は

10340

である．境界条件としては，壁面は

bounce-back

条件を用い

non-slip

条件とし，開境界には

M₂

分潮として以下のような入射条件を与えた．

ζ=Asin2π

Tt, ui= rg

hζn⁺_i (11)

上式において，振幅

A

は

0.21[m]

，周期

T

は

12.42[h]

とする．なお，渦動粘性係数

νe

は

10[m²/sec]

と仮定す

る．また，図

-7

（

b

）に領域を

20

分割した領域分割図を

示す．

(4)

解析結果として，図

-8

に上げ潮時と下げ潮時における流速ベクトル図を示す．図より，本手法は東京湾のように複雑な境界形状を有する場合においても，解析可能であることが示された．また図

-9

に，開境界処理の差異による検討を行うため，姉崎における水位の時刻歴を示す．図より，開境界に粒子分布関数

0

勾配を適用すると波高のピークが合わなく，周期性が見られないが，開境界に平衡分布を適用すると，実測値と良い一致を示していることが確認できる．なお，

1

周期目の波は初期条件の影響をうけるため，波高を過大に評価している．

また，図

-10

に演算速度倍率と並列化効率を示す．図より，大規模計算となる

Fine mesh

（総節点数

:18,887

）では高い並列化効率が得られており，本手法の有効性が確認できる．一方，

Coarse mesh

（総節点数

:5,839

）において，プロセッサ数が増えると並列化効率が下がっているのは，格子が粗いために，計算負荷に比べて通信負荷が卓越したことに起因する．

4.

結論と今後の課題

本論文は，

CIVA-

格子ボルツマン法による非構造格子に基づく浅水長波流れ解析手法を提案すると共に，水深変化を有する場合の精度，開境界条件処理及び並列化について検討を行った．数値解析例を通じ，以下の結論を得た．

•

^{並進過程に}

CIVA

法を導入した

CIVA-

格子ボルツマン法を用いることにより，浅水長波流れに対する格子ボルツマン法の非構造格子への拡張が可能となった．

•

水深変化を有する場合において，本手法の解析精度は

1

次精度で収束していくことがわかった．

•

開境界処理には，平衡分布が適切であることがわかった．

•

本手法は，高い並列化効率を得ており，大規模計算に有効であることが明らかとなった．

今後は，潮流以外の浅水長波流れ問題への適用，移動境界手法の導入などが挙げられる．

参考文献

1)

蔦原道久，高田直樹，片岡武：格子気体法・格子ボルツマン法

-

新しい数値流体力学の手法

-

，コロナ社，

1999.

2) J.G.Zhou: Lattice Boltzmann Methods for Shallow Water Flows, Springer, 2003.

3)

田中伸厚

:

数値流体力学のための高精度メッシュフリー手法の開発，機論，

64-620B, pp.1071-1078, 1998.

4)

立石絢也，樫山和男：

CIVA-

格子ボルツマン法による非構造格子を用いた非圧縮性粘性流体の解析，応用力学論文集，土木学会，

Vol.7

，

pp.323-329, 2004

．

5)

石川裕士，立石絢也，樫山和男：

CIVA-

格子ボルツマン法を用いた並列浅水長波流れ解析，第

20

回数値流体力学シンポジウム講演論文集（

CD-ROM

），

2006.

6)

石川裕士，立石絢也，樫山和男：非構造格子に基づく

CIVA-

格子ボルツマン法による浅水長波流れ解析，応用力学論文集，土木学会，

Vol.9

，

pp.231-238

，

2006.

80 90 100 110 120 130

20 40 60

x [km]

y[km]

Anesaki Funabashi

Yokohama

80 90 100 110 120 130

20 40 60

x [km]

y[km]

Anesaki Funabashi

Yokohama

(a)

概観図

(b)

領域分割図図

– 7

東京湾有限要素分割図

(a)

上げ潮時

(b)

下げ潮時図

– 8

流速図

0 10 20 30 40

-0.5 0 0.5

Equilibrium Zero Gradient Observed data

Time [h]

WaterElevation[m]

図

– 9

姉崎における水位の時刻歴

2 4 6 8 10 12 14 16 18 20 2

4 6 8 10 12 14 16 18 20

Speedupratio

Number of processors Ideal

Coarse mesh Fine mesh

2 4 6 8 10 12 14 16 18 20 0

20 40 60 80 100 120

ParallelEfficiency[%]

Number of processors

図

– 10

土木工学専攻

浅水長波流れ解析のための CIVA- 格子ボルツマン法の構築研究

土木工学専攻

号 石川 裕士

研究目的

が注目されている．格 子ボルツマン法は，アルゴリズムが簡便なため，計算が 高速に処理可能であり，大規模計算に対して有効な手法 であるといえる．

しかし，格子ボルツマン法は構造格子に基づく手法で あるため，自然地形を取り扱うことが多い浅水長波流れ 解析においては，その適用性に難点がある．これまで，

非構造格子へ拡張する試みとして，

流れ に対して，非構造格子における移流方程式に対する高精 度スキームである

法

を導入した

格子ボ ルツマン法が立石らにより提案された

．

そこで本論文は，

格子ボルツマン法に着目し，

それによる非構造格子に基づく浅水長波流れ解析手法を

提案する

と共に，水深変化を有する場合の精度，開

格子ボルツマン法

格子ボルツマン法

は，流体を有限個の速度をもつ 多数の仮想粒子の集合体で近似し，各粒子の衝突と並進 とを計算する手法である．流体の巨視的変数は，求めら れた粒子分布関数を用いて質量保存・運動量保存により 算出される．

離散ボルツマン方程式

粒子分布関数の時間発展は，以下の離散ボルツマン方 程式で表される．

上式において，左辺は粒子の並進過程を表しており，

格子ボルツマン法では粒子分布関数の並進計算に

法を導入し，非構造格子への計算を可能としてい る．また，右辺の

項目は格子

モデルに基づく衝 突演算項，

項目は力に関する項をそれぞれ示している．

は

方向の粒子がどのくらい存在するかということ を表す粒子分布関数，

は局所平衡分布関数である．

また，

は図

に示す粒子の並進速度ベクトルであり，

本論文で用いた

次元

速度モデルでは以下の式によっ て表される．

図

次元

速度モデル

式

，

の

は

で定義され，

は微小時間 増分量，

は格子サイズである．式

の支配方程式 により粒子の計算を局所的に行い，陽的に未知量

を 求める．また，

は力に関する項であり，以下のように 与えられる．

上式において，

は流体密度，

は河床高度，

は底面 でのせん断応力，

は自由表面でのせん断応力，

は コリオリ力である．詳細は参考文献

を参照されたい．

なお，式

における

の定義は図

のようになって いる．

また，式

において，

は単一時間緩和係数と呼ばれ る定数であり，

タイムステップの衝突で粒子が格子点 上において一定の割合

で局所的な平衡状態に近づい ていくことを示している．

は，鉛直方向に積分された 渦動粘性係数

と以下のような関係にある．

図

水深定義図

局所平衡分布関数

局所平衡分布関数

とは，格子点上において流体が 平衡状態に達したときの粒子分布であり，巨視的変数で ある全水深及び速度を用いて以下のように求められる．

上式において，

は全水深，

は流速，

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

．

流れの巨視的変数

号石川裕士

が注目されている．格子ボルツマン法は，アルゴリズムが簡便なため，計算が高速に処理可能であり，大規模計算に対して有効な手法であるといえる．

しかし，格子ボルツマン法は構造格子に基づく手法であるため，自然地形を取り扱うことが多い浅水長波流れ解析においては，その適用性に難点がある．これまで，

流れに対して，非構造格子における移流方程式に対する高精度スキームである

格子ボルツマン法が立石らにより提案された

は，流体を有限個の速度をもつ多数の仮想粒子の集合体で近似し，各粒子の衝突と並進とを計算する手法である．流体の巨視的変数は，求められた粒子分布関数を用いて質量保存・運動量保存により算出される．

粒子分布関数の時間発展は，以下の離散ボルツマン方程式で表される．

法を導入し，非構造格子への計算を可能としている．また，右辺の

モデルに基づく衝突演算項，

方向の粒子がどのくらい存在するかということを表す粒子分布関数，

速度モデルでは以下の式によって表される．

は微小時間増分量，

の支配方程式により粒子の計算を局所的に行い，陽的に未知量

を求める．また，

は力に関する項であり，以下のように与えられる．

は底面でのせん断応力，

はコリオリ力である．詳細は参考文献

のようになっている．

は単一時間緩和係数と呼ばれる定数であり，

タイムステップの衝突で粒子が格子点上において一定の割合

で局所的な平衡状態に近づいていくことを示している．

は，鉛直方向に積分された渦動粘性係数

とは，格子点上において流体が平衡状態に達したときの粒子分布であり，巨視的変数である全水深及び速度を用いて以下のように求められる．

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

流体の巨視的変数である全水深及び速度は，その粒子分布関数と式

の粒子の並進速度ベクトルを用いて以下のように計算される．

前述したように，格子ボルツマン方程式における力に関する項は以下のようになる．

項はいずれも流速に関わる項なので，毎ステップ更新する必要がある．

を並進過程に用いて計算する．

法は，移流方程式の高精度かつ安定な解法である

法を三角形及び四面体に拡張した手法である．

法では，各節点の局所厳密解を上流側の格子内に張った

次の補間曲面により求める．

計算過程において，衝突及び並進計算が終了した際に境界外から計算領域に入る方向の粒子分布関数

を決定する必要がある．すなわち，境界上の領域に向かう法線ベクトルを

を算出しなければならない．例えば，図

において下部に壁面があるとするならば，

を算出しなければならないことになる．

条件，流入流出境界条件としては，粒子分布関数の

勾配条件，平衡分布などが挙げられる．以下にそれらの計算方法を示す．

条件が広く用いられている．この境界条件処理は

粒子が来た方向に

跳ね返るという条件である．以下の式によって，

条件が広く用いられている．この境界条件処理は，壁に向かう方向の粒子が鏡面反射する条件であり，以下のように表される．

流入・流出境界条件としては，解析領域内に向かう方向の粒子分布関数の勾配を

が提案されている．例えば，図

において左から右へ流体が流入する場合を考えるならば，以下の式により算出される．（図

はその風下方向の隣接する格子点である．流出に関しても同様に処理する．

流入・流出境界条件としては，未知粒子分布関数を平衡分布にする条件がある．この境界条件は，巨視的変数である水深及び流速から算出するため，より質量保存及び運動量保存を満たすようになる．

本研究では，分散メモリ型並列計算機を対象とした領域分割法に基づく並列計算手法を導入する．