安定化 ISPH 法による大規模自由表面流れ解析に関する研究

修士論文要旨（

2011

年度）

安定化 ISPH 法による大規模自由表面流れ解析に関する研究

Studies on Stabilized ISPH Method for Large-Scale Computation of Free Surface Flow

土木工学専攻 

30

号 原田 悠里

Yuri HARADA

1. はじめに

砕波・飛沫等を伴う複雑な自由表面流れの解析手法と して，近年粒子法が注目されている．粒子法は，連続体 を有限個の粒子によって模擬するもので，各粒子は速 度・圧力といった物理量を保持しながら移動するため，

連続体の挙動をメッシュを用いることなく粒子の運動で モデル化できる．津波解析においても，特に津波が沿岸 部付近に来襲し，波の形状を複雑に変化させながら遡上 する現象のモデル化に有効的な手段となり得る．一方，

粒子法は粒子分布により密度を評価するため，粒子の均 等配置が再現できなくなると非現実的な圧力振動を起 こし，体積が保存されないといった問題点が残されてい た．この問題点を克服するため，非圧縮性粘性流体解法 として安定化

ISPH(Incompressible Smoothed Particle Hydrodynamics)

法

¹⁾

が提案された．

本研究では，安定化

ISPH

法を取り上げ，より安定で かつ高精度に解析を行うため，手法中の最適な緩和パラ メータの決定および渦粘性モデルの導入について検討を 行なった．また，本手法により沖合から侵入する津波挙 動のシミュレーションを可能とするために，沖側境界に おいて造波境界を導入することの検討を行なった．

2. 数値解析手法 2.1

支配方程式

非圧縮性粘性流体に対しては，

Navier-Stokes

の運動 方程式と連続式を用いる．ここで，射影法に基づく分離 型解法を用いれば，以下の手順にて速度・圧力の更新が 可能となる．

u

^∗

= u ⁿ + ∆t(ν ∇ ² u ⁿ + g) (1) u ⁿ⁺¹ = u

^∗

+ ∆u

^∗

= u

^∗

− ∆t ³ 1

ρ ∇P ⁿ⁺¹ ´ (2)

また，圧力ポアソン方程式は以下のようになる．

∇ ² P ⁿ⁺¹ = ρ ∇·u

^∗

∆t (3)

ここに，

u

^∗は中間流速，

u

は流速，

P

は圧力，

ρ

は密度，

ν

は動粘性係数，

g

は物体力である．

2.2 SPH

法

SPH

法は，影響半径

h

内に存在する近傍粒子点上での 値を用いて，

Kernel

関数

W

により補間することで，物 理量の算定を行う．すなわち，物理量

φ(x i , t)

を

Kernel

関数を用いて，

図

– 1 SPH

法の概念図

φ(x i , t) ≈ X N

j=1

m j

ρ j W (r ij , h)φ j (x j , t) (4)

のように近似する．下添字は粒子を特定するための番号 であり，

m _j

，

ρ _j

はそれぞれ粒子

j

の質量，密度を示す．

また，

r _ij = (|x

_i

− x

_j

|)

は，粒子

i

と近傍粒子

j

との距 離関数を表す．

h

は，影響半径を表す．

SPH

法の概念図 を図ー

1

に示す．また，空間微分に関する項に関しても 同様に，次式を用いて粒子近似できる．

∇φ(x i ) ≈ ρ i

X N j=1

m j

¡ φ j

ρ ² _j + φ i

ρ ² _i

¢ ∇W (r ij , h) (5)

上記の近似式

(4)

，

(5)

を用い，

Navier-Stokes

の運動方 程式の独立自由度である速度

u

と

p

を離散近似する．た だし，圧力・速度ラプラシアンの粒子近似には，

Morris

ら

²⁾

によって提案された次式を用いた．

速度ラプラシアン；

∇·(ν ∇u i ) = 2 ρ i

X N j=1

"

m j

"

ν r ij ∇W (r ij , h) r ² _ij + η ²

# u ij

# (6)

圧力ラプラシアン；

∇ ² P i = 2 ρ i

X N j=1

m j

"

P ij r ij ∇W (r ij , h) r ² _ij + η ²

# (7)

ここに，

η

はゼロ割りを避けるための安定化パラメータ である．文献

²⁾

に習い，影響半径

h

参照し

η ² = 0.0001h ²

と与えることにした．

2.3

安定化

ISPH

法

射影法を用いた

ISPH

法では，圧力については圧力ポ アソン方程式を解くことで陰的に評価する．ここでは，

従来法における圧力ポアソン方程式を導出する過程の相 違点についてのみ整理した後に，本手法を提示する．

2.3.1

速度発散ゼロ条件による定式化

まずは，速度発散ゼロ条件のみを用いた圧力ポアソン 方程式の定式化を示す．これは，

Cummins

ら

³⁾

の研究 で採用されている．ここでは，全ての粒子において「非 圧縮性流れ問題における質量保存則，すなわち速度発散 ゼロ条件を満足している」ものとする．これらの前提条 件を考慮すると，次式に示す圧力ポアソン方程式を得る．

∇ ² P _i ⁿ⁺¹ = ρ ⁰ ∇·u

^∗

∆t (8)

2.3.2

密度一定条件による定式化

次に，

shao

ら

⁴⁾

の研究で用いられた密度一定条件の みを用いた定式化を示す．まずは，一般的な質量保存則

(4)

を示す．このとき密度導関数においても予測子。修 正子分離を行い，次の時間ステップにおいて「密度が最 終的には初期密度と一致する」と仮定すれば，

(10)

が 導かれる．

Dρ

Dt + ρ∇ · u = 0 (9)

ρ ⁰ _i − ρ

^∗

_i

∆t + ρ

^∗

_i − ρ ⁿ _i

∆t = −ρ ⁰ ∇ · u ⁿ⁺¹ _i

= −ρ ⁰ ∇ · (u

^∗

_i + ∆u

^∗

_i ) (10)

ここで，予測子計算においては密度差が生じず，かつ速 度発散もないとすれば，上式は以下の様になる．

ρ ⁰ _i − ρ

^∗

_i

∆t = −ρ ⁰ (∇ · ∆u

^∗

_i ) (11)

最後に，

(11)

を用いて

(12)

を書き換えれば，次式に示 す密度一定条件を考慮した圧力ポアソン方程式を得る．

∇ ² P _i ⁿ⁺¹ = ρ ⁰ _i − ρ

^∗

_i

∆t ² (12)

2.3.3

緩和した圧力一定条件による定式化

密度一定条件を厳密に満足し続けることは困難であ る．つまり，近傍粒子数が固定され，粒子が完全な一様 分布を保持した状態でなければ，密度は一定とはならな い．従って，従来の

ISPH

法では，高周波の圧力振動が 問題とされることが多かった．そこで，瞬間的にある程 度の密度誤差を許容することで，長期的に密度変化が生 じないようなスキームへ修正する必要がある．浅井らに よって提案された安定化

ISPH

法

¹⁾

では，圧力ポアソン 方程式の従来型のソース項である速度の発散項に加え て，緩和パラメータを乗じた密度差の項を付加する．ま ず，圧縮性流体にも適用可能である一般的な質量保存則 を出発点とする．

Dρ

Dt + ρ∇ · u = 0 (13)

(13)

より，密度導関数を前進差分近似し，各物理量に 対して以下の様に粒子近似する．

∇ · u ⁿ⁺¹ _i = − 1 ρ ⁰

ρ ⁿ⁺¹ _i − ρ ⁿ _i

∆t (14)

ここで，修正子計算

(2)

を参照して，

∆u

^∗

_i = u ⁿ⁺¹ _i − u

^∗

_i

であることを考慮し，

(14)

を用いて，

(12)

を書き換え れば，

∇ ² P _i ⁿ⁺¹ = ρ ⁰ ∇·u

^∗

_i

∆t + ρ _i ⁿ⁺¹ − ρ _i ⁿ

∆t ² (15)

となる．ここで，次の時間ステップでの密度が初期密度 と一致すると仮定すれば，密度増分

∆ρ ⁿ _i

は次式により 与えられるはずである．

ρ ⁿ⁺¹ _i = ρ ⁿ _i + ∆ρ ⁿ _i ≈ ρ ⁰ (16)

また，ここでは瞬間的な密度誤差を許容する（微圧縮 性）ために，密度増分を本来より小さな値を与えるもの とする．

∆ρ ⁿ _i := α(ρ ⁰ − ρ ⁿ _i ) (17)

ここに，

α

は

0

〜

1

までの値をとる緩和パラメータであ る．これまでの仮定を用いれば，

(17)

を

(15)

に代入す ることにより，以下の圧力ポアソン方程式を得る．

∇ ² P _i ⁿ⁺¹ = ρ ⁰ ∇·u

^∗

_i

∆t + α ρ i 0 − ρ i n

∆t ² (18)

ここで，ソース項に追加した項は，非圧縮性の仮定を厳 密に満足していれば消去される項であり，実際には前ス テップでの密度

ρ ⁿ

を初期の密度

ρ ⁰

に修正する働きが ある．

2.4

渦粘性モデル

次に，

Smagorinsky

による渦粘性を導入することで，

特に流れが激しい状態においても安定して圧力が評価で きるように修正した．修正された

Navier-Stokes

の運動 方程式の速度ラプラシアン（粘性項）は，

(19)

で表さ れる．

ν e ∇ ² u(x i ) = 2 ρ _i

X N

j=1

"

m j

"

ν e r ij ∇W (r ij , h) r ² _ij + η ²

# u ij

#

(19)

ν t = (C S L) ² S S = X

i,j

p S ij S ij

S ij = 1 2

µ ∂u i

∂x j + ∂u j

∂x i

¶

(20)

ここで，渦粘性

ν e,i = ν i + ν t,i

を表し，

C S

，

L

はモデル 定数であり，それぞれ

Smagorinsky

定数及び，初期の 粒子間隔に依存する定数である．本研究では，

C S = 0.2

とした．また

S

は，ひずみ速度テンソル

S i,j

から評価 されるひずみ率である．

図

– 2

解析モデル

粒子間隔

(m) ∆t(s)

流体粒子数（x

× y × z）

case 1 0.015 0.005 2880 (45×8 × 8) case 2 0.010 0.001 9792 (68×12 × 12) case 3 0.005 0.001 78736 (136×24 × 24)

表

– 1

解析条件

図

– 3

緩和パラメータの影響（

case 3

の場合）

3. 数値解析例

3.1 3

次元ダムブレイク解析

数値解析例として，

3

次元ダムブレイク問題を扱う．

図ー

2

に解析モデルとその初期条件を示す．表ー

1

に示 すように粒子間隔を

3

通りに変化させて，緩和パラメー タの影響について検討を行った．検討は，水槽の右側の 高さ

0.01m

の点

S

（図ー

2

参照） における圧力値につ いて実験値

5)

との比較を基に行った．なお，壁面の境 界条件は

Slip

条件を課した．図ー

3

は，最も粒子間隔 の細かい

Case 3

の場合の緩和パラメータを変化させた 場合の結果である．図より，緩和パラメータの値により 圧力値に大きな差異が生じることが分かる．また，その 値を小さくすると，体積減少が顕著に表れ圧力値が過小 となることが分かる．更に，粒子間隔と緩和パラメータ の関係を調べるために，緩和パラメータを変化させて計 算し，点

S

での圧力に関する実験値との相対誤差（

0.37s

から

1.2s

までの平均）を示したものが図−

4

である．図 より，緩和パラメータが

0.0020

〜

0.0050

付近で最も誤 差が少ないことがわかる。ここで，ピーク値との相対誤 差を考慮すると，

Case3

の最適なパラメータは，

0.0045

付近と判断した．

Case1

と

Case2

についても同様に調

図

– 4

実験値との相対誤差

(case 3

の場合

)

図

– 5

最適な緩和パラメータを用いた場合の結果

べた結果，それぞれの最適な緩和パラメータは，

0.45

と

0.065

付近と判断した．これより，粒子間隔と最適な緩

和パラメータには相関関係が見られ，粒子間隔を小さ くすると，緩和パラメータも小さくなることが分かる．

図−

5

は，以上の結果より得られた最適なパラメータを 用いて計算した結果である．図より，粒子間隔に応じた 最適な緩和パラメータを用いることで，いずれの場合も 実験値と概ね良い一致を示すことがわかる．しかしなが ら，図−

5

からも分かるように，最適な緩和パラメータ を与えてもなお，圧力振動の問題は残されていた．ここ で，

Case3

の場合の渦粘性を考慮した結果を図−

6

に示 す．これにより，圧力振動を低減でき，実験結果の傾向 を適切に評価できることを確認した．なお，このときの

Case 3

の場合のピーク時

(0.36s)

の計算結果を図−

7

に 示す．カラー表示は圧力分布を示している．これより圧 力，流体形状共に妥当な結果であることが確認できる．

3.2

一定水深下における進行波の解析

次に，数値解析例として，一定水深下における進行波 の問題を扱う．ここでは，図−

9

に示す様に，造波水 槽を用いて，本手法の造波問題への適用性の検証を行っ た．表ー

2

に解析条件を示す．造波方法として，図−

9

図

– 6

渦粘性を考慮した場合の結果

(case 3)

図

– 7

渦粘性を考慮した場合の結果

(case 3)

図

– 8

解析モデル

の左端に設置している造波板に片振幅と周期を強制変位 として与える．なお，造波板は，左端の境界粒子点上に ストローク

d

に合わせた速度を与えることにより表現し た．今回は，ストローク

d

を

0.25m

，周期

T

を

1.6s

と した．なお，壁面には

slip

条件を課した．図−

9

に，波 の伝播状況を示す．図中の点は，水面形状として，この 時刻において表面粒子と判定された全ての粒子を示して いる．

4. おわりに

本研究では，安定化

ISPH

法を取り上げ，より安定で かつ高精度に解析を行うため，手法中の最適な緩和パラ メータの決定および渦粘性モデルの導入について検討を

ストローク

d(m) 0.25

周期

T(s) 1.6

流体粒子数（水平

×

鉛直

×

奥行）

93750 (750 ×25 × 5)

初期粒子間隔

(m) 0.02

微小時間増分量

(s) 0.002

解析時間

(s) 7.0s

動粘性係数

(m

²

/s) 0.02

緩和パラメータ

α

 

0.45

表

– 2

解析条件

図

– 9

波の伝播に伴う波高の減衰

行なった．また，本手法により沖合から侵入する津波挙 動のシミュレーションを可能とするために，沖側境界に おいて造波境界を導入することの検討を行なった．その 結果，以下の結論を得た．

•

粒子間隔の変化に応じて，最適な緩和パラメータ が異なることが明らかになった．また，最適な緩 和パラメータを用いることで，計算結果は実験値 と良い一致を示すことが明らかとなった．

•

渦粘性モデルを導入することにより，圧力振動を 低減することが可能となった．

今後の課題として，造波問題における実験値との比 較，津波遡上解析における沖側境界条件の処理法の構築 などが挙げられる．

参考文献

1) A.M. Aly, M. Asai and Y. Sonoda, Simulation of free falling rigid body into water by a stabilized incompre- sible SPH method, International Journal of Ocean Sys- tems Engineering, accepted

2) J.P. Morris, P.J. Fox, and Y. Zhu, Modeling Low Reynolds Number Incompressible Flows Using SPH.

Journal of Computational Physics 136, 1997

3) S.J. Cummins and M. Rudman, An SPH projection method, Journal Computational physics, Vol.152(2), pp.584-607, 1999

4) S. Shao, E.Y.M. Lo, Incompressible SPH method for simulating Newtonian and non-Newtonian ows with a free surface, Advances in Water Resources, Vol.26, pp.787-800, 2003

5)

胡 長洪，柏木 正：

CIP

法を用いた強非線形波・浮体相互作 用の数値計算

,

九州大学応用力学研究所研究集会報告

, 2003