実験 - 竪型円筒容器内の気液二相流解析

Fig.4‑1に実験のモデル図を示す。半径 Tout ニ 80[mm]のアクリル円筒容器に水道水を水面高さ h

=

200[mm]まで入れ、底面中央部に設けた多孔質体より気泡を吹き込んだ [57]。多孔

4三次元解析

質体にはガラスろ過板 (半径九

=

10[mm]

，

標準最大孔径 100^{r v}150μm)を用いた。これにオイルフリーコンプレッサから圧力調整器、流量計を介した空気を通すことにより気泡を発生させた。気泡の吹き込み体積流量は

Q

⁼¹^.67xl0‑⁶[m³js]とした。

まず全体の流れを定性的に把握するためアルミ粉を液相に混入し、流動状況の可視化を行なった。光源にはスライドプロジェクタ MasterHI‑LUX‑HRI000を使用した。スリットを施した光源からアクリル円筒にスリット光を照射し、中心軸を通る鉛直断面およびノ

k

平断面の流動を可視化した。カメラは NikonF3、レンズはマイクロニッコール F2.8，55mmレンズを用い、シャッタースピードは 1秒間開放で撮影した。なお、水平断面の撮影に際しては、水面より気泡が離脱する影響で表面が乱れ、上部からの可視化が行えないため、斜め下方あるいは斜め上方より撮影した。

次に、定量的な検討を行なうため、レーザードップラ一流速計(DANTEC社， 57NI0BSA )により液相の流速を測定した。測定は前方散乱方式で行ない、サンプル数は各測定点で 256

とした。

4 . 3 理論解析

4 . 3 . 1 仮定

理論解析においては、以下の仮定を用いた。

‑気液聞における相変化はない

‑液相、気相ともに密度一定(非圧縮)

‑解析領域内では温度は一様かつ一定

・気泡の合体、分裂や気泡相互の影響は考慮しない

‑乱流モデルは使用しない

4 . 3 . 2 連続の式

上記の仮定を考慮し、連続の式は以下となる。

生

E十マ・ α(^k^Vk)= 0

8t ⁽⁴^‑¹⁾

4 三次元解析 4.三次元解析

αL +αc

=

1 (4‑2) ^μ=μ^L(1 +α

c )

⁽⁴^・⁷⁾

4.3.3 運動量の式

本章では基礎方程式として、一圧力モデルと分散流モデルの二つを用いた。三次元解析を行うにあたって、数値計算が不安定になることが予想されるため、まずこれまでに実績のある一圧力モデルを用いた。数値計算が安定に行えることを確認した上で、次に分散流モデルを用

いての三次元解析を行った。ここで、第 2章の二次元解析の結果から、本章で用いた格子間隔の場合、一圧力モデルと分散涜モデルの計算結果に大きな違いはなく、本章の数値解析結果においても、モデルの違いによる影響は小さいと考えられる。

i)‑圧力モデル

ー圧力モデルの基礎方程式として、気相の粘性項を考慮しない Durstの用いた式 [59]に仮想質量力と揚力を考慮したものを用いた。

(液相)

(2)気液界面間抗力

単位体積当たりの抗力は、次式により表される。

T‑V‑1‑

v '‑

V一此一dG一α一D一C一

今︑

u一λせ

F D 一一

(4‑8)

ここで、抗力を求めるには抵抗係数 Coを定める必要がある。一般に抵抗係数は、レイノルズ数の指数関数になることが知られている。この抵抗係数は実験により数多くのモデルが提唱されているが、本計算では以下の Stokesと SchillerNaumannの関係式 [70]を用いる。

Cn二

j μ

/Re

U ‑‑ 1 24 ^I

‑ l 1 f e

⁽¹

+

0.15Reo.⁶⁸⁷⁾

(Re三2)

( 4 ‑ 9 )

LEE=‑α川 +αLFv+ FVM + F o + _FL‑ PLgaL

( 4 ‑ 3 )

ここで、 Reは気泡レイノルズ数Re

=

^ρLdBvr・/μL、dB^{は気泡直径である。}

(3)仮想質量力

単位体積当たりの仮想、質量力は以下のように表される [60][61] ⁰

(気相) FVM VM = = au_‑‑I.J::rc.oρ_rL_LJr.lCVM _‑JVM { i一一一一一D::c ̲ D̲VL

~

• m L Dt Dt J ⁽⁴^‑¹⁰⁾

iノVη

G一一三=一αc¥lp‑FVM ‑F o ‑FL ‑PcgαG

Dt (4‑4)

ここで、 C_VMは仮想質量係数であり、球形気泡に対する値 1/2を用いた。

(4)揚力

速度勾配をもっ流体中を気泡が運動するとき、抗力以外に揚力 FLが働くことが知られている。揚力 F_L に対する構成式として、 Drewらのモデル [64]を用いる。

ii)分散流モデル

以下に分散流モデルの基礎方程式 [6]を示す。

(液相)

αPL主主主 ‑ αL¥lp+αLFv + F o + FL ‑PLgαL

Dt (4‑5) _FL

=

CLPLαCVr X (¥l X VL) (4‑11 )

(気相) ここで、

C

_Lは揚力係数である。揚力係数に関しては多くの報告例

[ 9 2 ][ 9 3 ] [ 9 4 ]

^{があるが、}

それらの多くが単一球に関するものである。本解析対象では、複数の気泡を取り扱っており、

それら一つ一つがランダムに運動していると考えられる。そのため揚力係数に対して確かな値を得ることが難しい。そこで、 C_L= 0.5

，

⁰^.²

，

⁰^.¹を用いて数値計算を行い、比較を行った。

空竺 =

‑Fo‑FL+(PL‑PC)gαG

Dt (4‑6)

4.3.4 気液界面聞に働く力の構成式

(1)粘性力

μはサスペンションの有効粘性係数であり、 Taylorによれば、以下のようになる [62]。

51 52

4 三次元解析

Table 4‑1: Calculation condi tions

[ c ‑ ふ N o l

^mod^e^l I dB [mm] I CL

卜]

向山川

l single pressure l.0 0.5 slip 2 dispersed fiow 2.0 0.5 shp 3 dispersed flow 2.0 0.5 non‑slip 4 dispersed flow 2.0 0.2 non‑slip 5 dispersed flow 2.0 0.1 non‑slip

Table 4‑2: Boundary condition

(Liquid phase)

wall UL二 VL= WL = 0 top surface free or rigid condition (Gas phase)

wall ￨￨￨B(OuUbGb/leOs z sllp cond1Uon

top surface

I

Bubbles leave from a free surface immediately.

=δυG/δz=θωGIθz‑θαGlaz = 0)

4 . 3 . 5

計算手法

数値解析には差分法を適用し、時間差分に陽解法を、圧力場の解法には HS‑MAC法を用いた。また、液相の慣性項には Table4‑1中の Case1は Kawamura‑K uwahara法 [90]、Case 2から 5は QUICK法 [91]により近似した。その他の項は二次中心差分法を用いた。

初期条件は、液相は静止しているものとし、解析領域内に気泡は含まれていないとした。

また境界条件は、 Table4‑2に示した通りである。

格子数は Case1の場合、 (rx B x z) = (20 x 36 x 50)、その他は (24x 36 x 60)とした。後者の場合、格子間隔距離は6.r= 6.z = 3.3rnmとなる。

4 三次元解析

4 . 4 解析結果

4 . 4 . 1

凹型鉛直方向速度分布の構造

Fig.4‑2(a)は計算による液相の速度ベクトル図、 Fig.4・2(b)はアルミ粉による流れの可視化写真である。なお、可視化写真の中心軸付近では気泡が存在するため、液相の流れは見ることはできない。いずれの図においても、非軸対称な流れが観察できる。計算結果において、円筒の中心軸に沿って上昇する強い流れが見られるが、とれは中心軸付近に存在する上昇する気泡によるものである。一方、比較的大きい渦が、上昇する気泡の外側に見られる。これらは時間とともにその大きさと位置を変えていた。

Fig.4‑3は水平断面 (z= 0.16 [m])の、 (a)計算による液相の速度ベクトル図、 (b)アルミ粉による流れの可視化写真である。計算結果の中心付近の鉛直方向速度ベクトルは上昇する気泡の運動によるものである。同様に、白い密集した鉛直の線が可視化写真においても見られる。

また、いずれの図においても渦が見られる。これらから、同様の流動特性が得られていると考えられる。

Fig.4‑4は、位置 (r= 0

，

zニ 0.17[m])における液相の鉛直方向速度の過渡応答図である。

この高さ:z=0.17[m]は、 Fig.4‑2に見られるこつの渦の中心の高さにほぼ等しく、もっとも大きい揺動が観察された。図から、 t= 20 [s]より振動特性が支配的になっていると考えられる。

Fig.4‑5はz= 0.17 [m]における液相の鉛直方向速度の半径分布である。データは、実験においては 256サンプル、計算結果においては 240サンプルの時間平均を示している。ここで、サンプル数 256あるいは 240は、 Fig.4‑6において、サンプル数 200以上から平均流速値はほぼ一定となっていることから、十分な数であると判断したことによる。時間は計算結果では 13.1[s]、実験は 10^I'V20 [s]である。実験でのばらつきはトレーサー粒子のサンプリング条件によるためである。 Fig.4‑5から、実験結果は中心軸で速度が周辺より小さく、凹型分布となっている。この分布形状は、計算結果でもわずかに見られる。この凹型分布は次に示す振動流れにより生じたものと考える。

Fig.4‑7は、鉛直断面における液相の速度ベクトル図の過渡変化を示したものである。二つの渦がその大きさと位置を変え、結果として中心付近の上昇流は揺れていた。

Fig.4‑8は z= 0.17 [m]における液相の鉛直方向速度の大きさの水平断面分布の過渡変化を示したものである。等高線は正の値、つまり上昇流領域のみを示している。図から、上昇流領域は複雑に変化し、円筒の中心軸からずれている。

Fig.4‑9は t= 22.9 ^I'V26.2 [s]における鉛直方向速度が最大値を示す位置の軌跡である。位

4 三次元解析

置が半径方向と周方向の二つの方向に動いているのが見られる。おおよそ、 t= 22.9 ^{^‑'}23.5 [5]

と24.5^{^}^‑^'25.9 [5]では周方向に、も=23.5 ^{^}^‑^'24.5 [5]と25.9^{^}^‑^'26.2 [5]では半径方向に動いている。

Fig.4‑10は二種類の流れの振動パターンの模式図である。実線の曲線は液相あるいは気相の上昇流を表している。 (a)のように振動が半径方向に中心軸を横切るようなものである場合、

時間平均した半径分布は右図のようなものになると考えられる。一方、振動が (b)に示したように、円筒の中心軸の周りで周方向に振動している場合、時間平均分布は右図のような凹型になると考えられる。実際の場合は、これら二つの振動が同時に起きているものと考えられる。

これらのような、振動する流れは、各位置において、速度の変動を引き起こす。 Fig.4‑11 は z= 0.17 [m)における各位置の、 LDVにより測定した鉛直方向速度の時間変化である。変動は中心付近 (rj i out三0.125)で大きく、 γjrout勾 0.375^{^‑'}0.9375では小さい。しかし、側壁付近 rj r out巴 0.975では再びやや大きくなっている。中心付近の大きい変動は、液相自体によるものと、上昇する気泡の運動により生じるものが考えられる。一方、側壁付近の変動は液相の下降流により生じたものと思われる。

Fig.4^・12はz

=

0.17 [m]における液相の鉛直方向速度の標準偏差の半径分布である。ここで標準偏差は変動速度の二乗平均 (rm5)値を表す。計算結果は、定性的に実験結果と一致している。中心付近の標準偏差分布は、速度分布のそれと同様に凹型となった。このような分布の標準偏差は、 Iguchiらも同様に実験的に得ている [89)。しかし、標準偏差分布は速度分布と異なり、 Fig.4‑10の二つの振動パターンのいずれからも凹型分布は生じると考えられる。

4 . 4 . 2 揚力係数の影響

Fig.4‑13は、 z= 0.17 [m]における液相の鉛直方向速度の半径分布である。 Ca5e2と3の計算結果は中心付近で実験結果より値が小さく、なだらかな分布になった。このなだらかな分布は気泡の拡散により生じたものと推測される。本解析対象において、気泡の拡散とは主に水平方向への広がりを意味し、水平方向に動く力としては揚力が主である。つまり、この揚力の値が大きいものと推測される。一方、 C_L

=

0.2とした時の速度分布は、 ijrout ‑0.2付近で実験と比較して勾配が急ではないが、比較的一致している。 C_L

=

0.1とした時の速度分布は、

r jiout = 0.2付近で実験と同様に急であるが、中心付近で平らな分布とはなっていない。これ

らから、本解析対象においては、揚力係数として、 0.2が実験ともっとも良い一致を与えると考えられる。

Fig.4‑14は鉛直方向速度の標準偏差の半径分布を示したものである。 C_L= 0.2と0.1のと

4.三次元解析

きの計算結果が、実験結果とよい一致を与えている。 Zun[95]は理論解析において CLニ 0.3 とし、また Antalら[96]は C_L=0.1として、それぞれ管内流の流れにおいて実験と良い一致を得ている。本解析対象は、管内流ではない点でそれらと異なり一概に比較できないが、値と

して近く、ほぽ妥当な結果が得られているものと思われる。

Fig.4‑15は、

z =

0.17

[ m ]

上の三点における鉛直方向速度の過渡応答を示したものである。

計算結果は、これまでの結果で最も実験結果と良い一致を示した CL= 0.2の場合を示している。 (a)の実験結果において中心付近 (rjrout= 0.125)で激しく速度が変動している。この中で時折、速度が O近くになっているが、これは気泡運動によるノイズと推測される。しかし、

これらを除いても、速度は激しく変動しており、その様子は計算結果においても見られる。また、変動が高周波である点で、本数値計算がある程度の非定常特性をシミュレートできている

ものと思われる。 rjiout= 0.675

，

0.9538における速度の平均値は大略実験と一致している。

Fig.4‑16は、 Ca5e4 (C_{L ニ} 0.2)で、液相の半径方向速度の r

=

0.04 [m)に沿う鉛直分布を示したものである。計算結果は比較的よく実験結果に一致している。ノンスリップ条件を課した上側表面で特に良い一致が見られる。

4 . 4 . 3 測定高さによる鉛直方向速度の半径方向分布の違い

これまでは、高さ z= 0.17 [m)における鉛直方向速度の半径分布を中心に議論してきた。

そこで、ここではその他の高さにおける分布を調べた。 Fig.4‑17は、 z

=

0.19

，

0.18

，

0.10 [m) における鉛直方向速度の半径分布である。高さ Z^二 0.19[m]と0.18[m)は、 0.17[m]と同様に主流渦が存在する高さであり、 0.17[m]で見られた凹型分布が同様に見られる。一方、 z= 0.10 [m]では、可視化実験および数値解析結果においても観察されるように、気泡プルームは大きく揺動しておらず、速度分布も噴涜問題 [22][89] [97]によく見られる正規分布に近いもの

となっている。

このことから、主流渦の周辺では、流れは揺動しており、これまで述べた流れの揺動パターンが、z

=

0.17 [m)に限らず生じていると考えられる。

4 . 4 . 4 三次元非定常特性

これまでの検討で、本三次元計算と実験値は比較的良い一致を見せており、本計算がほぼ妥当であると考えられる。そこで次に、数値解析の時系列データの検討を行い、流動の三次元非定常特性の解明を試みた。計算結果は、実験結果と最も良い一致を示した Ca5e4に関して

ドキュメント内竪型円筒容器内の気液二相流解析 (ページ 40-46)

実験

=

=

，

Q

k

4 . 3 理 論 解 析

4 . 3 . 1 仮定

4 . 3 . 2 連続の式

生

=

c )

j μ

‑ l 1 f e

+

( 4 ‑ 9 )

( 4 ‑ 3 )

=

~

=

C

[ 9 2 ][ 9 3 ] [ 9 4 ]

，

，

空竺 =

[ c ‑ ふ N o l

卜]

向 山 川

I

4 . 3 . 5

4 . 4 解 析 結 果

4 . 4 . 1

，

=

4 . 4 . 2 揚力係数の影響

=

=

z =

[ m ]

，

=

4 . 4 . 3 測定高さによる鉛直方向速度の半径方向分布の違い

=

，

，

=

4 . 4 . 4 三次元非定常特性

4 . 3 理論解析

向山川

4 . 4 解析結果