平成18年度

平成

18 年度

修士論文

粒子法によるキャビティ流れの

数値シミュレーション

高知工科大学大学院

工学研究科 基盤工学専攻

知能機械システム工学コース

知能流体力学研究室

矢野 敦大

目次 第1章 緒 言 -1- 1. 1 はじめに -1- 1. 2 差分法と粒子法の比較 -2- 1. 3 研究目的 -3- 第2章 基礎方程式 -4- 2. 1 オイラーの方法とラグランジュの方法 -4- 2. 2 基礎方程式 -5- 2. 3 無次元化 -6- 第3章 粒子法による離散化 -7- 3. 1 重み関数 -7- 3. 2 勾配モデル -8- 3. 3 発散モデル -9- 3. 4 ラプラシアンモデル -10- 第4章 差分法による離散化 -11- 4. 1 差分化 -11- 4. 2 スタガード格子 -12- 4. 3 対流項の上流化 -12- 第5章 計算アルゴリズムと計算条件 -13- 4. 1 計算アルゴリズム -13- 4. 2 解析領域及び計算条件 -15- 4. 3 境界条件 -16- 第6章 計算結果および考察 -18- 5. 1 粒子法による計算結果(計算粒子の運動) -18- 5. 2 粒子法と差分法による速度分布の比較 -23- 5. 3 粒子法と差分法による圧力分布の比較 -24- 第7章 結 言 -25- 参考文献 -26- 謝辞 -27-

第1章 諸 言

1. 1 はじめに 近年，数値計算技術の発達とコンピュータの性能向上により，実験に換わる解析手法と して数値シミュレーションが工学分野などに広く普及している．特に，流体をコンピュー タ 上 で 計 算 す る 数 値 流 体 力 学 の 分 野 で は ， 計 算 方 法 と し て 古 く か ら 差 分 法(Finite Difference Method, FDM)が用いられており，従来ならば風洞実験など設備・時間のコスト がかかるものであってもコンピュータ上で容易に計算・可視化を行うことができる．また， 境界適合格子(Boundary Fitted Coordinate, BFC)を導入することで車体や翼周りなど複雑 な形状周りの流れを計算することができる． しかし，流体の数値計算手法として定評がある差分法だが，苦手とする領域がある．そ れは津波のように界面が流れとともに変形する自由表面流れの問題である．差分法でこの 問題を解くには各時刻で計算格子を作り直す必要があり，コンピュータに多大な負荷をか ける．また，VOF(Volume of control)法のように格子を作り直す必要のない場合であっても， 新たに関数を定義する必要があり，界面形状がなだらかになってしまう数値拡散の懸念も ある．そこで最近，このような流れ計算する手法として粒子法が注目を集めている．

1. 2 差分法と粒子法の比較

粒子法といっても広義では分子動力学法，格子ボルツマン法や個別要素法なども含まれ てしまうので本論文では Koshizuka・Oka が開発した MPS(Moving Particle Semi-implicit)法[1]

のみに限る．図 1.1(a)のように差分法では太枠で囲まれた解析領域を有限個の格子で区切り， その格子点上で速度や圧力といった物理量を隣り合う格子点と数値をやり取りしながら計 算を進めていく．また，図 1.1(a)の下側の図のように計算のやり取りをする隣接点は境界が 時間とともに変化しても固定である．それに対して，図 1.1(b)の粒子法では近傍の粒子と数 値のやり取りをすることは差分法と似ているが，図 1.1(b)の下側の図のように近傍粒子との 接続関係がなく，時間経過とともに変化していく． また，自由表面流れを差分法で計算する場合は前節で述べた方法などを用いるが，粒子 法では移動する界面を直接計算粒子の運動で捉えることができるので，格子を作り直す必 要もなく，数値拡散の影響を心配する必要もない．したがって，粒子法は自由表面流れの 数値計算法として非常に有効である．

(a) FDM _{(b) Particle method}

1. 3 研究目的 本研究では，泥上を走行するタイヤ周りの流れを数値計算することを最終目的としてお り，タイヤが走行したときの泥の波紋や飛沫の界面形状を描くための計算方法として粒子 法は非常に有効と思われる． 上記計算を粒子法で行う前に計算結果の信頼性を確認する必要がある．粒子法による計 算は自由表面流れの計算において多くのよい結果を出しているが，開発されてまだ 10 年程 しか経てないので，差分法など従来の計算方法に比べ信頼に足る計算サンプル数が十分で はない． そこで、本研究では上記数値計算を行う前段階として，古くから流体計算のベンチマー ク問題として有名なキャビティ流れについて粒子法を用いて数値シミュレーションを行い， 差分法による計算結果と比較して，その信頼性を検証した．

第 2 章 基礎方程式

[2] 2. 1 オイラーの方法とラグランジュの方法 流体の運動を観測する方法は 2 通り存在する．1 つは，観測者が動かず流れてきた流体粒 子を逐次観測する方法である．もう 1 つは，観測者が流れている流体粒子を追尾して観測 する方法である．前者はオイラーの方法と呼ばれ，後者はラグランジュの方法と呼ばれる． この 2 つの方法で流体の速度を観測する場合，オイラーの方法では観測点を過ぎる流体粒 子の速度から決定し，ラグランジュの方法では流体粒子を追尾する観測点の速度から分か る．速度についてはどちらの方法も自身の観測点の情報で決定することができる．しかし， 流体粒子の加速度を観測する場合，2 つの方法で変化が生じる．式で書くと， <オイラーの方法>

(

v

)

v

v

v

₊

_⋅

_∇

∂

∂

=

t

t

D

D

(2-1) <ラグランジュの方法>

t

t

D

D

∂

∂

=

v

v

(2-2) となる．ここで v は速度ベクトル，t は時間を表し，D/Dt は実質微分と呼ばれている． オイラーの方法では，式(2-1)で示された加速度は式(2-2)のように単なる速度の時間微分 とならず，右辺第 2 項目に対流項と呼ばれる項が付加されている．オイラーの方法で加速 度を観測する場合，自身の観測点のみでは流体粒子が加速しているかどうか厳密に判断で きない．そのため，周囲の観測点と速度情報をやりとりして加速度を決定する．その効果 が式(2-1)中の右辺第 2 項に対流項として表れている． 一方，ラグランジュの方法では観測点は流体粒子とともに運動しているので，流体粒子 が加速しているかどうかは観測点自身の運動から分かる．このため，式(2-2)のように加速 度は単なる速度の時間微分となる．

2. 2 基礎方程式 本研究で用いる流体は非圧縮性のニュートン流体である．よって，基礎方程式は質量保 存則である連続の式と運動量保存則であるナビエ-ストークス方程式である． 以下に，ベクトル表記した基礎方程式を示す． ＜連続の式＞

0

=

⋅

∇ v

(2-3) ＜ナビエ-ストークス方程式＞

v

v

=

−∇

+

μ

∇

2

ρ

p

t

D

D

(2-4) ある速度で運動すると周りの流体粒子にもその速度が伝播されることを意 味 ，検証に用いた差分法はオイラーの方法に基づく ため式(2-1)を用いて置き換えられる． ここで，ρは密度， p は圧力，μは粘性係数を表している．式(2-4)中の左辺の項は慣性 項と呼ばれ流体粒子の加速を表し，右辺第 1 項目は圧力勾配項と呼ばれ，流体が圧力の高 い所から低い所に向かって流れることを表している．また，右辺第 2 項目は粘性項と呼ば れ，流体粒子が している． なお，本計算に用いた粒子法はラグランジュの方法に基づくため，式(2-4)中の実質微分 は式(2-2)を用いて置き換えられる．一方

2. 3 無次元化 U を系の代表速度，L を系の代表長さとすると式(2-3)，(2-4)の各物理変数について，次の 関係式が得られる． ∗ ∗ ∗ ∗

₌

_∇

₌

_∇

₌

=

t

U

L

t

L

p

U

p

U

v

,

ρ

2

,

1

,

v

(2-5) これらの式を用いて，式(2-1)，(2-2)を無次元化すると以下のようになる． ・式(2-1)

0

=

⋅

∇

∗ ∗

v

(2-6) ・式(2-2) ∗ ∗ ∗ ∗ ∗ ∗

∇

+

−∇

=

v

v

2

Re

1

p

t

D

D

(2-7) ここで， μ ρUL = Re はレイノルズ数と呼ばれ，慣性力と粘性力の比を表している．また，上 付き添え字の*記号は無次元化したことを意味し，今後の式では省略する． 無次元化することで幾何学的相似・運動学的相似・力学的相似が満たされていれば，あ るレイノルズ数で数値計算された解は同じレイノルズ数の流れを予測することができる．

第 3 章 粒子法による離散化

[3] 3. 1 重み関数[4] 流れの基礎方程式である式(2-6)，(2-7)は微分方程式なので直接コンピュータで解くこと は不可能である．そこで，離散化という作業が必要になってくる．本研究で用いる粒子法 の一つである MPS 法は独自の粒子間相互用モデルを用いて微分方程式の離散化を行ってい る．離散化にあたって，本研究では以下の重み関数 w(r)を用いる．

( )

(

)

(

(

)

⎪

⎪

⎪

⎪

⎩

⎪⎪

⎪

⎪

⎨

⎧

≤

<

≤

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

−

<

≤

+

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

−

=

r

r

r

r

r

r

r

r

r

r

r

r

w

e e e e e e

0

5

.

0

2

2

5

.

0

0

2

2

2 2

)

)

(3-1) ここで，r は粒子間の相対距離，reは影響半径である． 次に，粒子数密度 niは重み関数の和として次式で定義される．

(

∑

≠

−

=

i j i j i

w

n

|

r

r

|

(3-2) なお，r は計算粒子の位置ベクトルを表し，添え字の i は計算粒子の番号，j は参照される 近傍粒子の番号を意味する．

3. 2 勾配モデル 任意のスカラーφ について、近傍粒子 j のスカラーφjを粒子 i のφiからのテーラー展開 で表すと，

(

−

)

+

K

⋅

∇

+

=

_{i ij j i} j

φ

φ

|

r

r

φ

(3-3) となる．１次の項まで考えると式(3-3)は

(

j i

)

j i ij

φ

φ

φ

⋅

−

=

−

∇

|

r

r

(3-4) となり，左辺は粒子 i，j 間の勾配ベクトルと相対位置ベクトルの内積になることが分かる． 式(3-4)について，両辺を相対位置ベクトルの絶対値で割り，さらに両辺に相対位置ベクト ル rj－ri方向の単位ベクトルをかけると，

(

) (

) (

)(

)

2

|

|

|

|

|

|

|

i j i j i j i j i j i j i j ij ij

r

r

r

r

r

r

r

r

r

r

r

r

−

−

−

=

−

−

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

−

−

⋅

∇

=

∇

φ

φ

φ

φ

_(3-5) となり，粒子 i，j 間の勾配モデルと定義される．ここで，左辺の< >は粒子間相互作用モデ ルであることを示す記号である．そして，粒子 i の勾配ベクトルは式(3-5)を複数の近傍粒子 j との間で式(3-1)，(3-2)を用いて重み付き平均したものとして定義される．(図 3-1 参照)

(

_{) ( ) (}

∑

≠

⎪⎭

⎪

⎬

⎫

⎪⎩

⎪

⎨

⎧

−

−

−

−

=

>

∇

≈<

∇

i j i j i j i j i j i i

w

n

d

|

|

|

|

r

r

2

r

r

r

r

φ

φ

φ

φ

)

_(3-6) ここで d は空間次元数(本計算では 2)を表している．

φ

i

i

j

r

e (φj‐φi)(rj‐ri) |rj‐ri|2

φ

_j

3. 3 発散モデル 任意のベクトル u について，相対位置ベクトル方向 r´の微分は | | j i i j r r r u u u − − = ′ ∂ ∂ (3-7) と書ける．発散の計算には式(3-7)のベクトル成分のうち r´方向成分のみ必要であるから， r´方向の単位ベクトルをかけて

(

) (

)

| | | | j i i j i j i j r r u r r r r r r u u − − ⋅ − − = ′ ∂ ∂ _′ (3-8) となり，この式(3-8)が粒子 i，j 間の発散と定義される．3-2 節の勾配モデルと同様に，複数 の近傍粒子 j との間で式(3-8)の重み付き平均を行うことによって，計算粒子 i における発散 は離散化される．(図 3.2 参照)

(

) (

_{) (}

∑

≠ ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ − − − ⋅ − = > ⋅ ∇ ≈< ⋅ ∇ i j i j i j i j i j i i w n d | | | |r r 2 r r r r u u u u

)

(3-9)

u

j

u

i

i

j

r

e

r

j

- r

i

3. 4 ラプラシアンモデル ラプラシアンは次式で離散化される．

(

) (

{

∑

≠

−

−

=

>

∇

≈<

∇

i j i j i j i i

w

d

|

|

2

2 2

r

r

φ

φ

λ

φ

φ

)}

_(3-10) 図 3.3 に示すように，式(3-10)は影響半径 re以内にある近傍粒子 j に重み関数の分布で値の 一部を分配することを意味する．また，式(3-10)の係数λiは統計的な分散の増加を解析解と 一致させるために次式で定義される．

(

)

{

∑

≠

−

−

=

i j i j i j i

|

|

w

|

|

2

r

r

r

r

λ

}

_(3-11) ラプラシアンはベクトル解析上では勾配に発散を作用させたものである．しかしながら， この考えに則って先に挙げた勾配モデルに発散モデルを作用させた式

(

)(

_{) ( )}

(

)

(

_{) (}

₎

∑

∑

≠ ≠ ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ − − − = − − − ⋅ − − − = ∇ ⋅ ∇ ≈ ∇ i j i j i j i j i i j i j i j i j i j i j i j i i ij w n d w n d | | | | 2 | | | | | | 2 2 2 2 2 r r r r r r r r r r r r r r φ φ φ φ φ φ 　　　　　　　　 (3-12) と MPS 法のラプラシアンモデルである式(3-10)とは整合しない．

φ

_i

j

i

φ

_j

r

e

第 4 章 差分法による離散化

[5] 4. 1 差分化 本研究で粒子法の検証に用いる差分法の離散化は一般に差分化と呼ばれている． 任意の関数 f(x)において，±Δx 離れた所の値 f(x+Δx)，f(x‐Δx)を求めるにはテーラー展 開を用いて，

(

)

( )

( )( )

( )( )

( )

( )( )

( )( )

⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ + Δ + Δ − = Δ − + Δ + Δ + = Δ + L L 2 2 2 2 2 2 2 1 ) ( 2 1 x dx x f d x dx x df x f x x f x dx x f d x dx x df x f x x f (4-1) と計算される．そして，1 階導関数

( )

dx x df および 2 階導関数

( )

₂ 2 dx x f d の近似式は式(4-1)から， 以下のように導かれる．

( )

(

) ( )

x x f x x f dx x df Δ − Δ + ≈ (4-2)

( )

( ) (

)

x x x f x f dx x df Δ Δ − − ≈ (4-3)

( )

(

) (

)

x x x f x x f dx x df Δ Δ − − Δ + ≈ 2 (4-4)

( )

(

)

( ) (

)

2 2 2 2 x x x f x f x x f dx x f d Δ Δ − + − Δ + ≈ (4-5) 式(4-2)は 1 次精度前進差分，式(4-3)は 1 次精度後退差分，式(4-4)は 1 階微分の 2 次精 度中心差分，式(4-5)は 2 階微分の 2 次精度中心差分と呼ばれている．なお，精度を向上 させる場合は離散点を増やすことで得られる． 粒子法では各微分演算子を直接離散化しているのに対して，差分法では微分演算子中 の各微係数を差分化しているという違いがある．

4. 2 スタガード格子 差分法の数値計算において 2 次元非圧縮性ニュートン流体の場合，求めるべき変数は x 方向速度 u，y 方向速度 v，圧力 p の計 3 変数である．差分法の計算では，図 4.1(a)のレギュ ラー格子と呼ばれる格子点上で全ての変数の計算を進めていくが，流体の数値計算の場合， 圧力振動を誘発する可能性がある．そこで，図 4.1(b)のように圧力を格子セルの中心に，速 度の各成分を格子点間の中心に配置したスタガード格子がある．この格子ではレギュラー 格子で生じる圧力振動は発生せず，連続的な圧力分布が得られる．

p

i,j

v

i,j+1

v

i,j

u

i+1,j

u

i,j i+1，j+1 i，j+1 i+1，j i，j

u

i,j

，v

i,j

，

p

i,j i，j+1 i+1，j+1 i+1，j i，j

(a) Regular grid

Fig. 4.1 Lattice geometry

(b) Staggered grid 4. 3 対流項の上流化 差分法の計算では，計算点が流れとともに運動しないオイラーの方法に基づくので，式 (2-5)の実質微分には対流項が含まれる．対流項は非線形であり，レイノルズ数の増加に伴 い計算格子数が十分でない場合，数値計算することが困難になる．この問題を解決する方 法の 1 つとして対流項の上流化が挙げられる．以下に，本計算で用いた１次精度上流差分 式を対流項の x u u ∂ ∂ について挙げると，

(

)

(

)

⎪ ⎩ ⎪ ⎨ ⎧ < Δ − > Δ − = ∂ ∂ + − 0 0 1 1 i i i i i i i i u x u u u u x u u u x u u 　　 　　 (4-6) 式(4-6)を１つにまとめると次式になる．

x

u

u

u

u

x

u

u

u

x

u

u

i i i i i i i

_Δ

+

−

−

Δ

−

=

∂

∂

+1 −1 +1

2

−1

2

|

|

2

(4-7) 対流項の他の項についても同様に差分化した．

第５章 計算アルゴリズムと計算条件

5. 1 計算アルゴリズム

本研究で用いる粒子法は SMAC(Simplified Marker and Cell)法と同様の半陰的アルゴリズ ムをとる．式(2-7)の圧力勾配項を陰的に解き，その他の項は陽的に求める． まず，初速度 vk_{が与えられると，次式で仮の速度 が得られる．} ~

v

⎟

⎠

⎞

⎜

⎝

⎛

_∇

Δ

+

=

k k

t

v

v

v

2 ~

Re

1

(5-1) ここで， は時間刻み幅である．次に， を用いて仮の位置 を求める．なお,(4-1)式中の 粘性項は(3-10)式を用いて離散化を行った．

t

Δ

v

~ ~

r

~ ~

v

r

r

=

k

+

Δ

t

(5-2) 次に仮の速度，仮の位置における式(2-6)を計算する．式(2-6)の離散化にあたっては式(3-9) を用いた．式(2-6)が求まれば圧力のポワソン方程式は次式で計算される．

t

p

n

Δ

⋅

∇

=

∇

+ ~ ~ 1 ~ 2

v

(5-3) (5-3)式の圧力のラプラシアンの離散化には(3-10)式を用いた．また，この式は系全体として の行列方程式になる．行列方程式を計算するには行列を直接計算する場合と反復計 算によって解を求める場合の 2 種類あるが，本計算では反復計算でありながら有限回数で 厳密解に達することができる ICCG(Incomplete Cholesky decomposition Conjugate Gradient)法 を用いて k+1 ステップ目の圧力 pk+1_{を算出した．そして，最終的に圧力勾配項から k+1 ステ} ップ目の速度 vk+1_{及び位置 r}k+1_{を確定する．}

b

Ax

=

1 ~ ~ 1 + +

₌

₋

_Δ

_∇

k k

p

t

v

v

(5-4) 1 1 + +

₌

k

₊

_Δ

k k

t v

r

r

_(5-5) なお，式(5-4)中の圧力勾配項の離散化には式(3-6)を用いた．

一方，差分法の計算で用いたアルゴリズムは粒子法の計算アルゴリズムと比較的似てい る部分段階(Fractional Step, FS)法を採用した．部分段階法では先の方法と同様，圧力勾配項 を陰に，その他の項を陽に計算する． まず，初速度vkから次式で中間速度 を求める． ~

v

(

)

⎭

⎬

⎫

⎩

⎨

⎧

_∇

₋

_⋅

_∇

Δ

+

=

k k k k

t

v

v

v

v

v

2 ~

Re

1

(5-6) 式(5-6)から得られた中間速度を用いて式(2-6)を計算し，次にそれを用いて圧力のポワソン 方程式である次式から k+1 ステップ目の圧力 pk+1_{を求める．}

t

p

k

Δ

⋅

∇

=

∇

+ ~ 1 2

v

(5-7)

式(5-7)の解法は反復計算法の一つである SOR(Successive Over Relaxation)法を用いた．最後 に，pk+1_{を用いて圧力勾配項を計算し，次式で k+1 ステップ目の速度 v}k+1_{を確定する．} 1 ~ 1 + +

₌

₋

_Δ

_∇

k k

p

t

v

v

(5-8) なお，各式の差分化について式(5-6)の粘性項と式(5-7)の圧力のラプラシアンは各項を 2 階微分の 2 次精度中心差分，式(5-6)の対流項は 4.3 節のとおり 1 次精度上流差分を用い た．また，4.2 節のスタガード格子における変数配置に基づき，式(5-7)の連続の式には 1 次精度前進差分，式(5-8)の圧力勾配項には 1 次精度後退差分とした．

5. 2 解析領域および計算条件 図 5.1 に示す赤枠で囲まれた正方キャビティについて，キャビティの一辺の長さ L を代表 長さ，上部境界面の x 方向速度 U を代表速度として式(2-3)，(2-4)を無次元化し，式(2-6)と 式(2-7)を導いた．そして，レイノルズ数 Re=100 から 500 まで 100 毎に変化した場合の計算 を行った． 本計算では，粒子法では 5.1 節で述べた計算アルゴリズムを用いて，時間刻み幅Δt = 5. ×10-4_{，ステップ数 30000 回で計算する．また，計算粒子の初期配置は l} 0 = 0.02 で等間隔に 配置し，計算粒子数 2601 個とした．一方，差分法の計算では同時間刻み幅･ステップ数の もと，x 方向の刻み幅Δx 及び y 方向刻み幅Δy をそれぞれ 0.02 とし，分割数を 50×50 にし て粒子法の計算と条件を合わせた．

U

L

L

x

y

0

5. 3 境界条件 速度の境界条件について，上部境界面で u = 1.0 , v = 0.0，その他の境界面ではノンスリッ プ境界条件 v = 0.0 を与えた． しかし，差分法の計算ではスタガード格子を用いているので，速度の境界条件を与える 際には，注意を要する．ノンスリップ境界条件を与える場合，4.2 節で述べたように u と v が同一点に存在しないため安易に u = v = 0.0 と設定することができない．そこで，図 5.2 に 示すように壁面上にある速度成分は 0.0 とし，もう一方の速度成分は壁内部の仮想セルにあ る速度成分と互い違いに設定する．こうすると，もう一方の速度成分についても壁面上で 0.0 にしたことになる． 本計算のノンスリップ境界条件は底面については図 5.2(a)に示すように

⎩

⎨

⎧

=

−

=

0

wall fluid virtual

v

u

u

(5.9) と設定し，また側面の条件については図 5.2(b)に示すように

⎩

⎨

⎧

−

=

=

fluid virtual wall

v

v

u

0

(5.10) と設定し，計算を行った．

v

virtual

v

fluid

v

fluid

u

virtual

u

fluid (b) Side (a) Bottom

u

virtual

u

fluid

v

virtual

u

wall

v

wall

一方，粒子法におけるノンスリップ境界条件では，図 5.3 に示すように壁面上にある壁粒子 の速度 vwall = 0.0 とする．しかし，粒子法では影響半径 re の推奨値を勾配・発散モデルで 2.1l0，ラプラシアンモデルで 4.0l0としている．このため，壁付近にいる計算粒子において 計算を行う際に参照する近傍粒子数が不十分になる．そこで，壁粒子の外側に仮想粒子を 配置し，壁付近の計算粒子で近傍粒子数を十分にしている．この仮想粒子についても速度 vvirtualを単に 0.0 と設定することもできるが，本計算では差分法におけるスタガード格子を 用いた場合と同様の互い違いの速度を与えた． しかし，差分法とは異なり計算粒子は運動しているので位置が一定ではない．そこで， (3-1)式を用いて重み平均を取り異符号化したものを仮想粒子の速度 vvirtualとする．具体的に は，仮想粒子と近くにある壁粒子との相対位置ベクトルを rbcとし，壁粒子からさらに rbc 内側にのばした位置ベクトルを中心に計算粒子速度 vjの重み付き平均を取る．

(

)

(

)

(

)

(

)

∑

∑

= = − − − − − = fluid j bc virtual j fluid j bc virtual j j virtual w w | 2 | | 2 | r r r r r r v v (5-11)

: Fluid particle

: Wall particle

: Virtual particle

r

e

r

bc

v

virtual

v

j

v

wall

第 6 章 計算結果および考察

6. 1 粒子法による計算結果(計算粒子の運動) 以下の図 6.1 から図 6.5 に粒子法を用いて計算した場合の計算粒子の運動を示す．各図の ステップ間隔は 10000 とした．なお，各図(a)のカラーバーは計算粒子の運動を見やすくす るために設定したもので，物理的意味はない．そして，この設定された色は時間によらず 一定なので粒子の運動や混ざり具合を確認することができる．

(a) 0 step (b) 10000 step

(c) 20000 step (d) 30000 step

(b) 10000 step (a) 0 step

(c) 20000 step (d) 30000 step

(b) 10000 step (a) 0 step

(c) 20000 step (c) 30000 step

(a) 0 step (b) 10000 step

(d) 30000 step (c) 20000 step

(a) 0 step (b) 10000 step

(c) 20000 step (d) 30000 step

Fig. 6.5 Moving particles for Re=500

キャビティ流れでは上部境界面がある速度で運動することによって，キャビティ内部で 渦が発生することが一般に知られている．各図において，上部境界面が運動することに伴い， 計算粒子が解析領域内を渦状に運動していることが色の分布から分かる．また，レイノルズ数の 増加とともに計算粒子が激しく撹拌される傾向にある．これはレイノルズ数を構成する慣性力が 増大したためと考えられる．

6. 2 粒子法と差分法による速度分布の比較 図6.6 にレイノルズ数 Re=200，30000 ステップ目における粒子法と差分法の速度分布を示す． 粒子法による速度点は差分法の各点の座標に合わせて重み付き平均したものである．両方とも渦 を形成していることは確認できるが粒子法に比べ差分法では，渦の中心がやや下方にある．また， 速度ベクトルがより下方に伝播していることが図から判断できる． 本計算では，差分法において対流項の差分化に１次精度上流差分を用いた．そのため，これは 数値拡散の影響が予想以上だったためだと考えられる． (b) FDM (a) Particle method

6. 3 粒子法と差分法による圧力分布の比較 図6.7 にレイノルズ数 Re=200，30000 ステップ目における粒子法と差分法の圧力分布を示す． 粒子法の圧力点の取り方は前節と同様である．なお，結果を比較しやすくするためにグラデーシ ョンのピーク値は任意に±0.15 とした． キャビティ流れの場合，右側上部で圧力が高くなり，左側上部で圧力は低くなる．また， 渦の中心では他の場所に比べ圧力は低い．図6.7(b)の差分法の結果ではこのことがよく現れ ている．しかし，図6.7(a)に示す粒子法の結果ではその傾向は現れているものの，キャビテ ィ側面で圧力振動が発生しており正しい解とは言い難い．この問題は本研究で用いた粒子 法特有のものであり，粒子間相互作用モデルの弊害といえる[6,7]_．

(a) Particle method (b) FDM

第 6 章 結 言

粒子法を用いてキャビティ流れをラグランジュ的に計算し，その粒子運動を可視化する ことができた．また，粒子法の結果を差分法の結果と比較すると，概ねよい結果が得られ た． しかし，速度分布の比較では差分法の結果において数値拡散の影響によるものと考えら れる結果が得られた．したがって，今後の課題として，対流項の高次精度上流化した計算 結果との比較も必要である．そして，圧力振動の問題は ALE(Arbitrary Lagrangian-Eulerian) 手法[8]_{，新たに固定粒子を配置すること}[7]_{などの方法を用いることによって改善される可能} 性がある． これら改善により，粒子法を用いたキャビティ流れをより正確に計算することが可能に なり，差分法による計算結果とも遜色がなくなると考えられる．

参考文献

[1] Koshizuka, S. and Oka, Y.: Moving-Particle Semi-implicit Method for Fragmentation of Incompressible Fluid, Nucl. Sci. Eng., 123, 421-434 (1996)

[2] 日野幹夫著，流体力学，1992，朝倉書店 [3] 越塚誠一著，粒子法，2005，丸善

[4] Koshizuka, S., Tamako, H. and Oka, Y.: A Particle Method for Incompressible Viscous Flow with Fluid Fragmentation, Computer Fluid Dynamics J., 4-1(1995), pp.29-46

[5] 峯村吉泰著，JAVA による流体・熱流動の数値シミュレーション，2001，森北出版 [6] 後藤仁志著，数値流砂水理学，2004，森北出版

[7] 亀井，白山，自由/固定粒子要素結合モデルについて，第 18 回数値流体力学シンポジウ ム D3-4，2004

[8] Koshizuka, et al.: Numerical Analysis of Natural Convection in a Square Cavity using MPS-MAFL, Computer Fluid Dynamics J., 8-4(2000), pp.485-494

謝辞

本研究を遂行するにあたり多大なご指導を賜りました，蝶野成臣教授ならびに辻知宏助 教授に対し深く感謝を申し上げます．

また，知能流体力学研究室の皆様に多大なご助力をいただきました．ここで，お礼申し 上げます．