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

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

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

桑木, 賢也

九州大学総合理工学研究科熱エネルギーシステム工学専攻

https://doi.org/10.11501/3135115

出版情報：Kyushu University, 1997, 博士（工学）, 課程博士 バージョン：

権利関係：

第 5 章

容 器 半 径 に 対 し て 水 面 高 さ が 低 い 場 合 の 解 析

83 

5容器半径に対して水面高さが低い場合の解析

5 . 1   はじめに

気液二相流は多くの工業プロセスにおいて出現する現象であり、さまざまな工業分野にお いて解析が行われている。

実際の工業プラントなどにおいては円筒容器が用いられることが多い。円筒容器内の気泡 流実験を行ったものとして、 Rietemaら[98]や Durstら[99]のものが挙げられる。これらに 対して、 Celikら[100]や Durstら[59]が数値解析を行っているが、いずれも二次元流動に関 するものである。これは、 三次元解析の困難性によるものと思われる。そのため、 Rietemaら は隔壁 (bufIle)を用いて、本来非軸対称な三次元流動になる流れを二次元的なものにして解析 を行っている。一方で、このような非軸対称な流れ、つまり気泡フルームが揺れる(揺動)流 動に対する解析も多く行われているが(例えば [87][88] )、そのほとんどが矩形容器内の二次 元流動である。これは、円筒容器を対象とする場合、 三次元解析が必要となり、その困難性の ためと思われる。そこで、第 4章においては、円筒容器内気液二相流問題に対して、 三次元数 値解析を行い、実験との比較から比較的良好な一致を得た。更に気泡プルームの揺動現象に関 連し、上昇流の円筒中心軸回りの振動パターンについて検討を行った。しかしながら、詳細な 流動構造は十分理解できたとはいえない。そこで、本章では4章のアスペクト比 :A = 2.5に 対して、水深の浅い A二1.0の場合の解析を行うこととする。これは、 A= 2.5の場合と比較

して、流動構造がより単純で、その詳細な解明には好都合と考えたことによる。

5 . 2   実験

5 . 2 . 1  

解 析 モ デ ル

Fig.5‑1に実験のモデル図を示す。半径 Tout= 80[mm]のアクリル円筒容器に水道水を水面 高さ h= 80[mm]まで入れ、底面中央部に設けた多孔質体より気泡を吹き込んだ。多孔質体に はガラスろ過板(半径九 =10[mr叶 ?標準最大孔径 100_~ 150μm)を用いた。これにオイルフ リーコンプレッサから圧力調整器、流量計を介した空気を通すことにより気泡を発生させた。

気泡の吹き込み体積流量は

Q

= 1.67x10‑⁶[m³js]とした。

5.容器半径に対して水面高さが低い楊合の解析

5 . 2 . 2  

実 験 方 法

気泡流に関する測定技術に関する報告としては、レーザ流速計を用いて気泡速度と液速度 の同時測定を行った大場ら [101]のものや、画像処理法による上昇気泡流の気泡径分布を測定 した賞雅ら [102]らのものがある。しかしながら、数値解析と同様、その測定方法は十分確立 されてはいない。そこで、本研究においては、単相流で使われる手法を用いた。

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

，

55mmレンズを用

い、シャッタースピードは 1秒間開放で撮影した。

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

とした。

LDVは光学系測定装置であり、気泡運動によるレーザー光の散乱といった問題が考えられ

る。事実、液単相領域での測定では信頼性の高いデータ(バースト信号)取得割合は約70~ 95% 

であるが、気泡の存在する領域で約 60%へ下がる。しかし、この 60%という値は許容範囲 内と判断した。本解析条件下のボイド率は数パーセント以下であるが、ボイド率が更に増した 場合、データ取得割合は更に低下し、測定は不可能となると考えられる。一方、気泡の存在し ない領域(中心軸付近以外)では十分な測定精度を得ているものと考えられる。 LDVを用いた Grossら [103]の結果においても、高ボイド率であるため、二相領域では流速は得られておら ず、液単相領域のみ結果は示されている。

5 . 3  理 論 解 析

5 . 3 . 1  

仮 定

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

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

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

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

85 

5.容器半径に対して水面高さが低い場合の解析 5容器半径に対して水面高さが低い揚合の解析

‑気泡の合体、分裂や気泡相互の影響は考慮しない ^μ=μL (1 

+ α

0) 

( 5 ‑ 5 )  

‑乱流モデルは使用しない (2)気液界面間抗力

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

この中で、気相を非圧縮としているが、松本ら

[ 1 0 4 ]

が示しているように、気泡径 (気泡 の圧縮性)に及ぼす水深(水圧あるいは静圧)の影響は、水深が1.2m の場合で、 水面付近と 容器底面付近で差が約

4%

である。本解析モデルは水深がわずか

0 . 0 8

mであり、松本らと同 様の計算を本系に対して行うと、静圧の気泡径に与える影響は 0.4%以下である。そのため、

水圧が気相の圧縮性に与える影響は無視できるものとした。

V一

V一L‑B P一dG一α一

D‑

C一

qd 一刈 せ

一 一

D F 

( 5 ‑ 6 )  

( 5 ‑ 1 )  

ここで、抗力を求めるには抵抗係数 C_{D を定める必要がある。} 一般に抵抗係数は、レイノ ルズ数の指数関数になることが知られている。この抵抗係数は実験により数多くのモデルが提 唱されているが

[ 6 5 ][ 6 8 ]   [ 6 9 ]   [ 1 0 5 ]

、本計算では以下の StokesとSchjllerNaumannの関係式

[70]を用いる。

Cn = 

1 

^:4

1 ̲

Re 

1.../  1 

l  f

24 

f e  

⁽11 ¹ 

+ 

I ⁰i¥ ^.15Re1 1":" n ̲^u0 ^.f^bdf'>R7¥ ⁾ 

(Re三2)

(5‑7) 

5 . 3 . 2  

連続の式

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

生~+マ (α山)二 O

_θt  ^(Reど2)

αL+α0=1  (5‑2) 

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

= 

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

(3)揚力

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

[ 6 4 ]

^{を用いる。}

5 . 3 . 3  

運動量の式(分散流モデル)

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

(液相) FL = CLPLαOVr X (¥7 X VL) 

( 5 ‑ 8 )  

αLPLE主 ‑ αL¥7p +αLFV +FD十FL‑PLgαL 

Dt  (5‑3)  ^{揚力係数には、第 4章の結果より、}

C

_L= 0.2とおいた。

(気相)

5 . 3 . 5  

渦度ベクトルと速度ベクトルポテンシャル

GEZ

三= ‑FD‑FL+(ρL‑PG)g

α G  

Dt 

( 5 ‑ 4 )  

本章では計算結果の表示に、 一部、速度ベクトルポテンシャルの成分の分布を用いている。

速度ベクトルポテンシャル及び渦度ベクトルと速度の関係を、三次元円筒座標系で表すと以下 のようになる。

ここで、 Fvは粘性力、 FD^{は気液界面問抗力、} FL^{は揚力である。}

5 . 3 . 4  

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

(1)粘性力

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

¥7

2ψ=

一ω

( 5 ‑ 9 )  

8 6  

87 

5.容 器 半 径 に 対 し て 水 面 高 さ が 低 い 場 合 の 解 析 5.容 器 半 径 に 対 し て 水 面 高 さ が 低 い 場 合 の 解 析

'¥7

2 ψ = 

θ2ゆr. 1 Oψγψγ lδ2ψγ 2θψB .θ2ψγ 

一一一 十一一一一一一一 十一 一 一 一 一 一 一 十一一‑

ar2  rθr  r2' r2θ^{()2  r2δ()  ， az2}  δZψ_B  1θ妙。 ψB 1δZψ。 2θψr o^2ψ。

一一一δr2 +' 一一一一L or 一一r2 ' 一 +一r2一一一θ0そ ァ一

^

2~ f^)(，一)  +，一一一θZ2  θZ

ゆ

^z 1δ

ゆ

^z 1δz

ゆ

z θ2

ゆ

z

一一一θγ2 +一一rδγγ2δ()2θ一 十一一一一 十一一一Z2 

(5‑10) 

心方向に向かつて水平である。その結果、底面側壁付近の二次渦の領域は実験結果に見られる

=角形とはなっていない。

5 . 4 . 2  

三次元解析

ω =ニ

1θωL  OVL  rθ()  ̲ oz  oiii θωL 

oz δr  θυL υL  1δUL  一 一 一 一 一

δr  r  rδ0 

︑︑ ︐

l︐r

噌lよ

Ttム

vh

υ 

f'

at︑

前節において、主流渦の下側の流れの向きは、実験では容器底面に向かつて斜めに下降し ているのに対して、計算では中心方向に向かつて水平であった。これは、格子数の影響以外の 要因によるものと考えられる。そこで、次に三次元数値解析 (Table5‑2中の CaseNo.5)を 行った。鉛直断面における格子数 (rx z) 

= 

(40 x 40)は前節の結果から十分とは言えないが、

計算機の能力の関係上、計算時間が妥当であると思われる最大限の格子数とした。なお計算機 には、株式会社ダイナス製 GAIA‑EV56(CPU : DECchip Alpha21164A，クロック周波数:433  MHz)を用い、 300^ヲ000サ イ ク ル (t 

= 

3.27 [s])の計算に約 58時間 30分を要した。

Fig.  5‑3は、計算による二乗平均速度と二乗平均ボイド率の過渡応答線図 (a)と、高さ z/h

= 

0.8の中心軸における鉛直方向速度の過渡応答線図 (b)である。まず鉛直断面における流れ (Uave

，

ωave)が発達し、その後、周方向の流れ (Vave)が発達している。気泡を吹き込んでしば らくは軸対称を保つが、流れが発達してくると、時刻 t= 20 [s]あたりから非軸対称となる。

Fig. 5‑4は、鉛直断面における速度ベクトルポテンシャルの周方向成分 (ψ。)の分布と、 z/h

= 0.5の水平断面における速度ベクトル線図の過渡変化を示したものである。図から、主流渦 の大きさは、周方向の流れが発達するにしたがってその大きさを増している。更に、流れが十 分発達した後 (Fig.5・4(c)， (d))、主流渦の下側において、容器底面に向かつて斜めに下降す る流れが見られる。この流れは、前節の二次元計算では得られなかった。

Fig.  5‑5は、 r/ r ou.t二 0.75に沿う()‑z断面におけるベクトル線図である。図から、まず周 方向の流れが生じ、その後、この周方向の流れはそれ自身が下降流(鉛直方向流れ)へ変化し たり、あるいは多くの渦を生じさせている。その結果、多くの場所で下降流が生じている。こ の流れが主流渦の下の斜めに下降する流れを生じているものと推測される。

Fig. 5‑6は、実験および計算による液相速度の分布を示したものである。 (a)は、 7・/rout

= 

0.5における半径方向速度成分の鉛直分布であり、 (b)は、 z/h=0.5における鉛直方向速度成 分の半径分布である。データは、実験においては 256サンプル、計算結果においては 240サン プルの時間平均を示している。時間は計算結果では 13.1[s卜 実 験 は 10rv 20 [8]である。実験 でのばらつきはトレーサー粒子のサンプリング条件によるためである。図から、計算結果は実 験と大略一致している。特に、 Fig.5‑6 (a)中の固体壁で近似した上側表面付近において、良

い一致が見られる。

5.3.6  計算手法

数値解析には差分法を適用し、時間差分に陽解法を、圧力場の解法には HS‑MAC法を用 いた。また、液相の慣性項には QUICK法 [91]により近似した。その他の項は二次中心差分法 を用いた。 一方、 Eq.(5‑9)の解法には SOR法を用いた。

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

また境界条件は、 Table5‑1に示した通りである。この中で、液相の自由表面における境界 条件は、第 3章の結果から固体壁で近似した。

計算条件は Table5‑2に示した通りである。 Case5の場合、格子間隔距離はムr ‑ムz ‑ 2.0mmとなる。

5 . 4   解析結果

5 . 4 . 1  

格子分割数の影響

Fig. 5^・2は、格子分割数の影響をみるため、二次元モデル (Table5‑2中の CaseNo.1 ^rv  4)により得られた結果を、液相の流れ関数線図で示したものである。鉛直断面において、中 心付近の上昇する気泡により引き起こされた液相の上昇流が、大きい渦(主流渦)を生じてい る。比較的粗い格子を用いた計算による主流渦 (Fig.5‑2 (a)， (b))は、可視化写真と比較する と小さい。粗い格子に高次風上差分法を用いた計算による主流渦が、実際より小さくなる傾向 は、 2章においても見られた。一方、 (c)，(d)での主流渦の大きさはほとんど同じである。この ことから、本系に対しては、格子数 96x 80が適当であると考えられる。しかし、主流渦の下 側の流れの向きは、実験では容器底面に向かつて斜めに下降しているのに対して、計算では中

88  89 

5.容器半径に対して水面高さが低い場合の解析

5 . 4 . 3  

気相流動の解析

Fig.  5‑7は、気泡吹き込み口から上側表面までの気相の流線である。図から、気相の流れ はほぼ軸対称である。これは、気相の上昇速度と比較して水面高さが低いためと考えられる。

このような結果から、本系の流動を二次元的なものと考えがちであるが、前節の結果から、液 相に関しては、流動は三次元的なものであった。それにもかかわらず気相の流れが軸対称に近 いのは、液相の三次元的な流れの影響を十分受けることなく、気泡が上側表面に達するためと 思われる。このことから、本系においては、液相の三次元的な流れが支配的であると考えられ る。そこで、液相の流動に関して以下に更に詳しく検討を行った。

5 . 4 . 4  

非定常特性

Fig.  5‑8は、鉛直断面の 30x 30の始点から描かせた短時間の流線である。ここで、そ れぞれの線の長さは、 1秒間に液体が動く運動に一致するが、ある瞬間の流れを示したもので ある。この 1秒間という時間間隔は、カメラのシャッター開放時間に一致する。可視化写真は、

Fig. 5‑8 (h)に示している。ここで計算による流線は、厳密には可視化写真と一致しない。な ぜなら可視化結果は流跡線を意味しているためである。図から、主流渦の大きさと位置は時時 刻刻と変化している。さらに主流渦の下側の涜れは、延びたり縮んだりを繰り返している。

Fig.  5‑9は、 z/h= 0.5における水平断面の 30x 30の始点から描かせた流線を示してい る。それぞれの線の長さは、 Fig.5‑8の場合と同様である。図中で、外側に中心に向かう線が 見られるが、これは Fig.5‑8で見られる主流渦の下側部分の中心に向かう流れに相当する。ま た、その中心に向かう疏れと中心付近の上昇流の聞に、いくつかの渦と流れのよどんだところ が見られる。これらもまた、時時刻刻と変化していた。

Fig. 5‑10は、 t

= 

36.0 [s]においてさまざまな高さにおける水平断面の 30x 30の始点か ら描かせた流線を示したものである。高さ z/h= 0.75 ^{r v}  0.375 (Fig.  5‑10 (b) ^{r v}  (e))において、

中心付近の上昇流の横に二つの渦が見られる。この二つの渦に関して、以下に更に詳細な検討 を行う。

Fig. 5‑11は、 Fig. 5‑10 (d)に矢印で示した視点から見た中心付近における流線を示した ものである。白い球は流線の終点を示す。中心付近の上昇流の周りに、二つの鉛直方向ら旋流 れが見られる。これらの位置は、 Fig.5‑10 (b) ^{r v}  (e)の二つの渦の位置に一致する。このこと から、二つの渦は、ら旋流の断面を表していると考えられる。また、図から、左のら旋流の流 線は、右側のものより密である。これは、左のら旋流の上昇速度が右側のものより小さいこと

90 

5.容器半径に対して水面高さが低い場合の解析

を意味していると考えられる。このような鉛直方向のら旋流は Fig.5‑12

， 

Fig.  5‑13に見られ るように、他の時間においても見られた。

Fig. 5‑12は時刻 t= 39.3 s[

ト

Fig.5‑13は t= 42.6 [s]における中心付近の流線を示して いる。 Fig.5‑12は Fig.5‑9 (d)の下側から、 Fig.5‑13は Fig.5‑9 (e)の右側から見た図であ る。いずれの図においても、中心付近の上昇流のまわりにら旋流が見られるが、位置は常に変 化しており、またここでは一つしか示していないが、その数もまた変化していた。このような 流れは、 Fig.5‑8のような鉛直断面の短時間の流線あるいはベクトル図を用いて表現すること は難しい。このような鉛直方向のら旋流れを描いた報告例はこれまで知らない。このら旋流れ の結果は、さらに深いあるいは大きい容器中の気泡プルームの揺動現象の構造の解明の手がか りになるものと期待される。

本系は、水面高さが低いため、顕著ではないが、気泡プルームはわずかに揺動していた。

これも鉛直方向のら旋流から影響を受けているものと推測される。

5 . 5   結言

容器半径に対して水面高さが低い場合の気液二相流に対して、数値解析と実験の両者を行 い、以下の結果を得た。

(1)主流渦の下部に容器底面に向かつて斜めに下降する流れが、可視化写真により観察された。

この流れは、三次元数値解析によりシミュレートできた。

(2)主流渦の下の容器底面に向かつて斜めに下降する流れは、周方向流れが発達することによ り得られた。

(3)多数の点より描かせた流線により、中心の気泡プルームの周りにおいて、鉛直方向に動く ら旋流れを含む詳細な流動パターンが得られた。

本章の系は、水面高さが低く、容器も比較的小さいため、流れはほとんど軸対称で二次元 と考えられがちである。しかし、実際には液相の流れは強い三次元性を持ち、周方向の流れが 支配的であった。本章で得られた、鉛直方向に動くら旋析しれを更に詳細に検討すれば、水深の深 い場合や寸法の大きい容器内の気泡プルームの揺動現象の解明に寄与できるものと期待される。

91 

5.容 器 半 径 に 対 し て 水 面 高 さ が 低 い 場 合 の 解 析 5容 器 半 径 に 対 し て 水 面 高 さ が 低 い 場 合 の 解 析

Table 5‑1:  Boundary condition 

(Liquid phase) 

wall  UL =υL =ωL=O  top surface  rigid condition  (Gas phase) 

wall  slip condition 

top surface  Bubbles leave from a free surface immediately.  (

δUO/θz =θυ0/θz =δω0/θz =θα0/θz 

= 

0) 

Table 5‑2:  Calculation conditions 

R 

^{r  む﹃σ 一ムr ‑1・ 一G1 山一0 一m 一01 一e 111llI﹁Fi Z} 

I  24  24 

2  40  40 

3  96  ^ー 80 

4  120  120  5  40  36  40 

GAS INLET  r o u t  

Fig.  5‑1:  Problem schematic. 

92 

5.容器半径に対して水面高さが低い暗合の解析

︾︿ぬ﹁

ω G O

︿

OE

ω

﹃ ﹁

n t

コ o

[ l ]

2.5x103 

2 . 0  

1  . 5  

1 . 0  

0 . 5  

30  35 

0 吋 25 

フ ﹄

﹁

lL

e  m 

﹁ コ ・

l

1 T   5  10  14 

2 

‑ ハ

u

n u  

10  12 

8  6  4  16x10‑3  { ω

¥ E ] h t υ ω 一 ω

﹀

ω o

︒何﹂

﹀︿

ヨ 5

3 4 3 (

同 日 o υ

( a )  Average v e l o c i t y  and v o i d  f r a c t i o n  

30  35  25 

15  20  Time  [ 5 ]   5  10 

0.14 

0.00  0  0.12  0.10  0.08  0 . 0 6  

0.04  0.02  [ ω

¥ ε ] h t u

一

o ω

﹀一

向

υ

一ピ む﹀ 一

ω υ 0

﹂

z d

同 { ( 同

﹄日)EVZO∞ω

口 一 ‑ 56 22 U4 4

ロ0

2ω AE ロロ 刃包

﹄ 0ぢω恒

ω ω

﹄'

H. N ︐

m . m r

出

方旬

︒E

‑d

己O

B R U E

弓﹃

︒﹀

P4

一 戸 一 コ ω ω

﹂℃

ω N

= ∞ コ

ω

一 ﹀

( ω )

ONF

×

ONF(℃) 

︒ 寸 ×

O

寸

( 心 )

寸

N

× 寸

N

︒ ∞

×

ω

め

( υ )

( b )  L o c a l  v e r t i c a l  v e l o c i t y  a t  r  =  0 ，  z/h  =  0.8 

Fig.  5‑3:  Computed transient characteristics of liquid velocity and void fraction.  95 

( 悶 )

5容 器 半 径 に 対 し て 水 面 高 さ が 低 い 場 合 の 解 析

[ω]ω .0

FH

一戸

( υ

)

. の . 0 H

﹄¥ N

克己

︒一 ぢC ω‑ sg NZ OH {2

#

口 一

ω H O

ちま長一

υ o

‑ F

甘 口 d gz υω

∞

‑ d υ

一 亡 ECE45

口O ( 同 日

833

口3

0 ( ぢo

ご

E125

む切

出口

υぃ

一 υ

M O

口O

右ロ

ハロ

相官

一ガ

古川

ぷ ω_{口 ︑} EH L ザーの・切出. [ω

]N .ω NH

一 戸 (

℃ )

{ ω

]F .的

FH

一 戸 ( 心 )

i k ^{昨 J :}

::::た:;:;:;会長~~~~~-~':~~~ 三三:

『一.....一二 三届顎松 山いい・・・ 一一

三三-湾側~~:::.::::三ご

/〆.....，'  .: ! : ・ ¥ ヘ ¥

. 

[ ω ]

川 町 内

ωH

一 戸 ( 何 )

Fig^， 5‑5:  Velocity vectors in ( )‑z section along ^i / i out 

= 

0.75 

>~~Wf;シヶン--

・ ども在京主‑"

......ミHαお))l:;:::::::

::‑:::::':‑0:勿術::':'::.̲.. 

5.容器半径に対して水面高さが低い場合の解析 5容器半径に対して水面高さが低い場合の解析

口

Experiment  S i m u l a t i o n  

60x 1 0 ‑

³ 

20  40 

Radial v e l o c i t y   [ m / s ]  

口

。

1 . 0  

0.8 

......， 

0 . 6  

0 . 4  

0 . 2  

0 . 0  

工¥

N

( a )  Ve

吋

i c a ld i s t r i b u t i o n  o f  r a d i a l  v e l o c i t y  along r / r o u t   =  ^{0 . 5}  

(a) Bird's‑eye view 

(b) Side view 

口 口口口

Experiment  S i m u l a t i o n  

口

ロ

口

0.10 

0.04  0.02  0.00 

ー

0.02 0.08  0.06  { ω ¥ ε

]  

KA to o‑

60﹀O

一亡

︒﹀

Fig. 5‑7:  Stream lines of gas phase from the bottom to top surface at t 

= 

36.0 [s]. 

1 . 0  

( b )  R a d i a l  d i s t r i b u t i o n  o f  v e r t i c a l  v e l o c i t y  a l o n g  

z/ 

h  =0.5  0 . 8  

0 . 6  

1

^・

1

0.4  r / r o u t   0.2 

0.0 

99  Fig.  5‑6:  Time averaged velocity distribution of liquid phase. 

5.容器半径に対して水面高さが低い場合の解析

( h )  V i s u a l i z e d   r e s u l t  

Fig.  5‑8:  Stream lines from 30 x 30 points in a vertical section (.d.t 

= 

^1.0 [8])^・

100 

5.容結半径に対して水面高さが低い場合の解析

( c )   t = 3 6 . 0  [ 5 ]   ( f )   t = 4 5 . 8  [ 5 ]  

Fig. 5‑9. Stream lines  from 30 x 30  points  in  the horizontal  section  at  z = 0.04  [m] 

( ム

t

= 

1.0 [s]). 

101 

5.容器半径に対して水面高さが低い場合の解析

( g )   z / h   =  0 . '  25 

、

( d )  z / h   =  0 . 5  

Fig.  5‑10:  Stream lines from 30 x 30 points in various section at t 

= 

36.0 [s] 

( ム

t

= 

1.0 [s

J )  

102 

5.容鰭半径に対して水面高さが低い場合の解析

Fig. 5‑11. Stream lines around central upward flow region from the viewpoint indicated by  an arrow in Fig.5‑9 (d) at t 

= 

36.0 [s]. 

103 

5.容器半径に対して水面高さが低い場合の解析

(a) Bird's‑eye view 

(b) Side view 

Fig.  5‑12. Stream lines around central upward fiow region from the view point on the bottom  side in Fig.5‑8 (d) at t = 39.3 [s]. 

104 

5.容器半径に対して水面高さが低い場合の解析

Fig.  5‑13. Stream lines around central upward flow region from the view point on the right  side in Fig.5‑8 (e)  at t 

= 

^42.6 [s]. 

105 

第 6 章

総 括

6 総括

本論文では、気液二相流に関して、多くの工業プロセスで取り扱われる円筒容器を解析対 象とした。解析モデルとしては、水道水に満たされた円筒容器の底面中央領域より、圧縮され た空気を多孔質体円板を通して吹き込んだ場合の流動である。数値解析と実験の両方を行うこ とにより、数値解析モデルの比較あるいは高精度な解法を検討した。その上で円筒容器内の三 次元的な流動の構造の解明を試みた。

第一章では、本研究の目的を述べるとともに、種々の気液二相流の数値解析方法と、気液 一相流に関連する既往の研究について数値解析を行ったものを中心にまとめた。

第二章では、 二次元の気液二相流に対して、これまで解析例の少ない、多次元問題におけ る二流体モデルへの高次風上差分法の適用を試みた。数値解析モデルとして、分散流モデルと 一般的に用いられるー圧力モデルを用い、比較を行った。また、高次風上差分法として4種類 のものを用い、これらの比較も併せて行った。その結果、格子間隔の粗い場合、 二つのモデル はいずれも収束解を得ることができ、その差は小さいものであったが、いずれも主流渦が可視 化写真より小さかった。一方、格子間隔が密な場合、分散流モデルでは解が得られ、主流渦の 大きさも可視化写真のそれに近いものとなったのに対して、一圧力モデルでは発散し、解は得 られなかった。このことから、高次風上差分法を用いる場合、十分な格子数が必要であり、そ のためには、より安定と考えられる分散流モデルが有効であることが分かった。また、 4種類 の差分法の中では、 QUICK法が、粗い格子においても比較的、実験に近い結果を与えていた。

このことから、分散流モデルに QUICK法を適用したものが本解析に有効であると考えられる。

第三章では、本解析において液相として用いている水道水の汚れを考慮し、自由表面にお ける液相の境界条件に関する検討を行った。境界条件として、スリップ条件とノンスリップ条 件を課した二つの場合の数値解析を行い、可視化実験結果およびレーザードップラ一流速計に よる測定結果との比較を行った。その結果、上側表面側壁付近において、可視化写真に見られ る二次渦が、スリップ条件を用いたものでは見られなかったのに対して、ノンスリップ条件で は得ることができた。また、その付近の詳細な解析の結果、ノンスリップ条件とした場合、 二 次渦の領域について実験結果と良い一致が得られた。さらに三次元解析から、実験で得られる ような高周波変動がノンスリップ条件のとき得られた。このことから、水道水を用いた本解析 条件においては、ノンスリップ条件が厳密ではないものの、有効であると考えられる。

第四章では、三次元解析を行い、一般的に知られている気泡プルームの揺動現象の流動の 解明を試みた。数値解析には、液相運動方程式の移流項に高次風上差分法の内、 QUICK法(一 部 Kawamura‑K uwahara法)を適用し、数値粘性の影響を極力除去することにより、揺動現象 を高精度にシミュレートできるようにした。まず時間平均した数値解析結果と実験との比較お

107 

6総括

よび速度変動の大きさから、数値解析がほぼ妥当であることを確認した。その際、揚力係数 C_L に関する検討も行い、

C

L

= 

0.2のとき実験値に近い結果を与えることを示した。また、鉛直方 向速度の半径分布が、中心付近で周りより値が小さくなる凹型分布となった。この分布から、

上昇流が、中心軸を通る半径方向に振動するものと、中心軸周りで周方向に振動する二つの振 動パターンを考え、数値解析結果の詳細な検討において、この二つの振動パターンを見いだし た。更に、実験で得られる速度変動を、非定常数値解析により得られた流動パターンにより説 明を行った。

第五章では、四章のアスペク卜比:A 

= 

2.5に対して、水深の浅い A

= 

1.0の場合の解析 を行った。これは、 A= 2.5の場合と比較して、流動構造がより単純で、その詳細な解明には好 都合と考えたことによる。数値解析には、二章の結果から分散流モデル、四章の結果から液表 面をノンスリッフ条件とした。まず、二次元解析により格子数の影響を調べた。その結果、粗 い格子の場合、主流渦の大きさが可視化写真のそれより小さくなった。これは、二章の結果と 同様な傾向であった。一方、格子数 96x 80以上のとき、主流渦の大きさはほぼ一定であるこ とから、格子数として 96x 80程度必要であることが分かった。しかしながら、主流渦の下側 の流れが、可視化写真では斜めに下降するのに対して、計算では中心方向に水平なものであっ た。これは、三次元計算の結果から、三次元的な流れ、つまり周方向の流れが生じ、その流れ が下降流を生じているためと思われる。水深の浅い場合、気泡プルームの揺動性は小さく、流 れは二次元軸対称に近いと考えられがちである。しかし、解析の結果、液相の流れは複雑で三 次元的なものであった。その中で、数値解析により得られた流線から、気泡運動により生じる 中心付近の上昇流の周りに、二つの鉛直方向に軸をもつら旋流れが存在すると予測された。こ のようなら旋流れは従来報告されたことがないものと思われる。このら旋流れは、その数や位 置は時間とともに変化していた。今後、このら旋流れは実験的にその存在が検証される必要が ある。このような流動構造が、気泡フルームの揺動現象に影響を与えているものと思われ、更 に深いあるいは大きい容器内の揺動現象の解明の手がかりになるものと期待される。

108 

付 録 A

無次元化

本論文中では、基礎方程式および結果に関して、有次元値のみで示している。しかし、実 際には数値解析にあたっては、無次元化された基礎方程式を差分化するととにより、数値計算 を行い、得られた結果を再び有次元に換算し、それらの値を示している。そこでここでは、本 解析に用いた無次元化を示す。

A.l  ^{基礎方程式の無次元化}

無次元化に関して、ここでは簡単のため二次元モデルに関して説明する。 二次元円筒座標 系に対する‑圧力二流体モデルは以下のようになる。

連続の式

θα""  1 θ A  

‑ a t "   + 

~ãr

( r a k u k )   + 会

α(

山 ) 。 =

^(A‑1) 

( ;tM 

^(A‑2) 

運動方程式(液相)

r方向

vL

一F一

+ 一

M

一向

日 川

一位

十 一

D一

F

一

+ 

V F 

﹄ 一

7む

包

宇 一

AY

+ 

均一かl‑L 

d

一

ρ'

一 一

仇一

&

+ 

ω  u r u

門 町一か

L U 

+ 

仇一白 (A‑3) 

z方向

A無 次 元 化 A，無 次 元 化

θ ω L θ ω L θ ωL  1θP . 1， ..，

F :

Dz + FVMz + FLz 

一 一 一 +肌 一 一 一

+ωι

一 一 一 ̲ ‑ ‑ = ‑ : ̲ l "

+‑‑=‑Fvz+ ~~， .，..~ . ~~-g

δt ^{' ‑ L 1} θi  ' ~1J az ρL az ρ L αLPL 

( A ‑ 4 )  

^{運動方程式(液相)}

r方向

運動方程式(気相)

r方向

δUc θUc θUc  1θP 

F .

Dr + FVM^r + FLr 

一一

^{+UG‑ ‑}

十

ωG

一一

^{= ‑ ‑ 7}

一

δt ^，‑'"δi  '"δz  pc ai αcpc  (A‑5) 

。

_θUL. T r  aUL δUL  Po δF  VL 本 αG

Zζ

万 干 一

^{+ UL}

万 五 +

^WL 

‑ a ‑ i  _{=一万戸?否五十三}

ZFVT

十三

^;(FBT+FCMT

十円

γ) (A‑15) 

z方向

z方向

δ ω G δ ω G δ υc  1δP 

F .

Dz + FVMz + FLz 

θt  ^t 

UG3F‑E 何 万 二 一 =‑Z5J‑

αcpc  ‑y  (A‑6) 

i)'(θWr θWr θWT 

一二一一一三 +UL_uot

。

δ^{ァ む}~ーとθR+ W"^む

一

θZ^L1 (A‑16) 

ここで、

。

δF ^， VL  n*αG 

= 

一一一一

_PL143δZ 

+

~F~z + ‑'" (FDz + FVMz + Fiz)

ーヰ

. ioUo αL  U~

れr ‑

話

^(γμ

^{L (}

^1+αc)

警

^{) + 2}

ト

^L(l+αG)(3h

引} ‑

^2μ

^{L (}

^1+αc)

^告

^{(A‑7) 運動方程式(気相)}r方向

ん

z=;

ま い L (

^l+αG)

( 守 + 告 ) }+ 場

^(ι(1+α

み)

^(A‑8) 

^一一一一

^uot

^{。 。}

^δUθ 7 θ R θ Z^{C+ Uc}

^一一十

^{θUC Wa}

^{一一二一一}

^θUC PL143θR ^{Po‑‑‑;;P本}

^一一〆

^{δP f"  (F〆¥  Dr}

⁺ 

^F~Mr

⁺ 

^Fir)  (A‑17)  3Cnαορ仁

F t=‑ u uμIvc ‑VLI(Vc ‑V

L )  

4  d_B  (A‑9)  ^z方向

FLiニ CLPLα^{C(V C‑}VL) ^X 

( ¥ 7  

^X VL)  (A‑10) 

。

δWc  δWC  θ日ノc Po ^本δF * ^I n*  ^{T:1*  T:1* ¥} 

一一一一

_uot

_。

_δァ

+

^U c ‑_δR ^{r"¥ r..'"  + W C}

^{一一=一一}

_{δZ  PLU~}^{‑p‑‑f(F5}_θZ ，‑ ¥  ^z+FCMz+FZz

ノバ

⁾

^{一 ‑}

^( A‑18) 

、p"'、p"'.....，.‑三

」 ー 」 ー に 、

次に無次元化のため以下の無次元変数を導入する。

叫一%

W  L  一 一

K一%一 一

L U 

ら 一

一 一 九

B D 

Z

一九

z 

一 一

r

一九

R 

一 一 _(A‑11)  ^FVr 

= 誌

^{R(l+α

^{a )} 警 } + 引

^(1+αG) 

(警+引} ‑

²⁽

¹⁺

^αc)

_主

^(A‑19) 

ωc  t T"I  P ^本 PL

UC

二一一

_{tL0 3}1Wc=‑‑‑=‑‑::₁₄::₀'₃

，

T=

一

_to 

，

F =

一 ，

p‑

一 一

Po  PG  (A‑12)  ηz

誌

^{R(l

⁺

^αG) 

( 写 会

^+32))+2

針

^(1+α

^{c )} 警)

^(A‑20) 

これらの関係式を式 (A‑1)"‑式(A‑6)に代入し、整理すると、

連続の式

3 Cn 

=ーニ4D_B IVc‑VLI(Vc ‑VL)  (A‑21) 

Fii ⁼ CL(V G ‑ V L) ^X 

( ¥ 7  

^X V L)  (A‑22) 

~ acxl‑ 1θ 

' 0

一三+ 一一

(RαkUk)

十 一

α(kWk)

= 

⁰ 

uot

。

δT . Rδ R δ Z  

( A ‑ 1 3 )  

{~会伽ゐ)+会(仇)} ⁺  { 品

^(RαLUL)

^{+ 表}

^αLWd}=O(

^{( A ‑ 1 4 )}  

ここで、式(A‑19)

，

(A‑20)中の(1+αG)^{は、粘性係数として Taylorの式} [62]μ=μL(l+αG) を用いたことにより生じた項である。

以上より、 UL1WL1  UG1 WG1αL

，

αGは以下のように表せる。

110  111 

r I io  PL  Po  vL  gio ¥  ULヲWL

，

UcヲWcぅαL

，

αc=ハ 一 一 一 一 一 一 一 一 一 ￨

J  ¥ Uoto ' PC^うPL143'ToUo?U3/

A無 次 元 化

( A‑23) 

無次元化の方法として Hellumsand Churchillの方法を用いて、 代表値を順次決めていく。

νL/(iOUO) = 1より Uo ⁼ VL/io  gíO/U~

= 

^{1  より 1}

= 

gío/U~

= 

gí~/vI

" 0  =  (;) 

U^ ‑ νL . ，̲ 

‑ ν(

l /  

g 

f ' 0 

iO/(UotO) = 1より to = io/U。

い え ( 子 )

Po/(PLU;) 

= 

1  より PO =ρLU; 

PO 

= 

PLVZ (

三 r

³

*  PL  P = pc 

これらを代入すると、無次元化された式として以下の式を得る。

連続の式

aa~ 1δ 

ーヱ +一 一(RakUk)

+ 

^'::1‑/7 α(kWk) = 0  θT  ' Rδ R δ Z  

( 誌

^(RαcUc)

⁺ 表 ( 叫 ) } ⁺  { 主 主 ( ゐ

^LU

山会

^α(

爪 ) } ⁼ 

⁰ 

運動方程式(液相)

r方向

δUL δ_UL δ_UL δP 事 α_c( 1:1*  1:1*  D *  ¥ 

一一^δァ

+

U_L ‑^{δ R θ}/'¥‑.;: 

+ 

W_L

一 ^Z 一 一 一 ^{8R +}

^{FVrγ α L}   

⁺ 

^{‑¥.:1 (Far} 

⁺ 

^FVMr 

⁺ 

^Fir)J  

z方向

112 

(A^司24)

(A‑25) 

( A‑26) 

(A‑27) 

( A‑28) 

(A剖)

(A‑30) 

(A‑31) 

A.無 次 元 化

θWL δWL θWL θ P α G  

a7~ + U_L

万五

^{+ W}L

万三一=一死 +

FVz +

ζ

^(Faz 

+ 

FVMz 

+ 

^{FLJ ‑1  ( A‑32)} 

運動方程式(気相)

r方向

θ_UC δUC θUC θP 

万 r+UG375+WG 万 z‑ 〆 万 五 ‑

p* 

( F ふ +

FVMr 

+ 

FLr)  ( A‑33) 

z方向

θWC θ_WC θ_WC δP   "'/ 、

一一一

+

U c ‑ ，δ R θ Z θ...~\J

+ 

W C一 一 ‑

‑〆‑‑ z 

^{‑ρP 事 ¥.L} (F;;^Dz z

+ 

^T F^{.L V 本M z z T} 

+ 

_FL^.L Lz _z) ^J ‑1  ^(A‑34) 

これらの式を差分化することにより数値解析を行った。なお、分散流モデルあるいは三次 元への拡張も全く同様の手法により無次元化を行うことができる。

A.2  物性値および無次元数

本解析においては、水ー空気系気液二相流を想定しており、物性値は文献 [106]によれば 以下の通りである。

Table A ‑1: P arameters to calculate at t = 20 [0 C] 

PL  ^998.2  kg/m³  ρG  1.161  kg/m³ 

νL  1.00 X 10‑⁶  m2/s  μL  1.002 ^X 10‑³  Pa.s 

σ  7.25 ^X 10‑²  N/m  g (definition)  9.80665  m/s2 

これらの値を用いて無次元化の参照量を求めると以下のようになる。

iO 

= 

4.672 X 10‑⁵ [m] 

Uo 

= 

2.140x 10‑² [m/s] 

113 

A.無 次 元 化

to = 2.183x10‑³ [sec]  po = 4.573 X 10‑¹ [Pa] 

また、エトベス数:Eo、モルトン数 :Mの値を求めると、

ハudAせ Qu

qJ 

qJ  

1 i F h J U  

〆E

tt EF

︽ ︑ EEEI︑

一 一 2B一︐α一 G一ρ'一 一 一

σ

AY一 r u一

GJ

一一 一

O E  (d_B = 1.0[mm]) 

(d_B = 2.0[mm]) 

(A‑35) 

付 録 B

gμi(PL ‑PG) 

M 二 hh=2.580×10‑8

P I

σu 

( A‑36) 

二流体モデルにおける HS‑MAC 法

logM = ‑7.588  (A‑37) 

気液二相流数値解析においては、数値的に不安定になりやすいため、安定性を重視し、陰

解法が用いられることが多い。そのため圧力場の解法には、単相流解析で用いられる SIMPLE 法あるいは SIMPLER法を二相流用に拡張したもの [107]が用いられる。例えば、 SIMPLER 法を用いた松浦ら [50][51]や、 Sokolichinら[86][87]の解析が挙げられる。

一方、陽解法に関しでも、 SIMPLER法などと同様に、単相流解析で実績のある HS‑MAC 法 (SOLA法)を用いた解析例 [108][109]がある。陽解法はその数値計算上の性質から時間刻 みを大きくとれない欠点を持つ。しかし、非定常流動を取り扱う場合、流動の非定常性が強い 場合、物理的に時間刻みを小さくする必要がある。そのような場合、マトリックスを解く必要 のある陰解法と比較して、 1サイクル当たりの計算時間の少ない陽解法は有利となる。

ここでは、本解析において圧力場の解法に用いた二流体モデルに対する HS‑MAC法の計 算手I}債を示す。なお、 HS‑MAC法に関しては、文献 [110]に分かりやすくまとめられている。

B.l  基礎方程式の変形

式 (A‑30)^{f ' J}  (A‑34)に仮想質量力の構成式を代入すると、基礎方程式は以下のようになる。

連続の式

(古志向山)+ 品川町 } + { 主 義 胤仇 ) + 会

^α(

^LWd}=O

^{'l︑ r R U   'ai  ノ︑︐ ︐ ︑}

運動方程式(液相)

114  115 

B 二流体モデルにおける HS‑MA仁法 B二流体モデルにおける HS‑MAC法

F~ ^{+  : :   {} 

^Fo 

⁺ 

^CVM (

守 ^{‑ T ) +} 吋+ ^{[ 0 ，  ‑ 1 ]}  

^{(B‑2)  D V一一一こ}_Dt ^η 

⁼

^‑kc¥lP 

⁺ 

^{G  (B‑7)} 

運動方程式(気相) 仏 ー 」 ー に 、^{......、.".‑r..、}

守

^{=‑p*¥lP‑}

小

^{/ai︑ ロυ 凡リq ︑︐ ノ︑l} kL 

= 

^{1 + CVM}

^〆

^{+ CVM}

^竺

^αL 

^{三 〆} 1  + 

C_VM

〆 ^+CVM22

αL 

(B‑8) 

これらの式は、それぞれ一つの式中に V^L

，

Vcに関する実質微分 (DjDT)を含み、そのま ま計算することは難しい。そこで、 DVLjDT

，

DVcjDTに関する連立方程式を解き、整理す ると以下のようになる。

運動方程式(液相)

k~=

_c= 

1

¹ 

^+C

^{川 CVM}

^{ご ) 〆}

1+CVMf+CVM22 

αL 

(B‑9) 

B . 2  HS‑ 孔 1AC 法の手順

(

1

^{十 川 + C}VM

ご)守

+  0 [ ^{，  ‑} ( 1  + 

^CVMP* 

+ 

^CVM

ご ) ] ^{( B ‑ 4 )}  

連続の式 (B‑1)は Fig.B‑1のような不等間隔スタッガード格子を用いると次のように離散 化できる。なお、 Fig.B‑1の中で便宜上、スカラー量のうち、圧力は丸印、ボイド率は三角で 定義するものとする。

連続の式

‑ ‑(l+C_{¥ α L}vMP*^十CVM

竺川マ ^P +  ( 1   +  CVMP*)F~ + 竿

^(FD

⁺ 

^FL) 

ノ UL

運動方程式(気相) 1 XeαGeU21 一九αωU

此

^1αGnwt1-α

Cs~ ど

¹

X_p  s x ム_U

( 1

  ^{+ 川}

^{+ C}_VM

ご)守

二 一

( 1 + 

^CVM 

+ 

^CVM U'.c 

i 

^p*¥l P 

^{+  C} VM〆 F~

^‑p*(FD 

⁺ 

^FL) 

¥ α Lノ

+ ト ‑   ^{1 ( +} 

^{C川 本 十 CVM::)]} 

+ 2 .

ー

.  

^{XeαLeutf1一 九}αLwU221+αLnwtJ1~ ー αLswkfl一 一げ

。

ら ムX ムy ~

︑︑ ・

1 ''

'

n u  

‑

‑ よ

RU

 

r'

EE︑

さらに運動方程式

( B ‑ 6 ) ， ( B ‑7 )

は、以下のように離散化できる。

運動方程式(液相)

(B‑5)  ukf‑U

土

1̲ P，~ ^‑P， 

J..JCム7M=FE十 n v

J 

^円

r :

^..，..，kL ^︑z J ^hzF 1Ei 

‑ ‑ A R   U 

J' a

a ‑︑

これらをまとめると次のようになる。

運動方程式(液相) ukt1-U~." ^̲1̲  P."‑P ，.

U山A T b W二

Fω+

^{T'¥ "t.rT T}  T 山T"\よ ~kL (B‑12) 

刀Vr

一 一 二_Dt 

=‑k

r，¥l P 

+ 

F 

.~

( B ‑ 6 )  

W

. f : ; l   ‑ wt 

^{̲ 1̲}  P ，.‑ P ，.

J..Jn  山

=

F:十 一 一TP.  ; yT:  I Ir ^1t T k L  ^︑︑

︐

l︐fつd11ムロυ

f'EE︑︑

運動方程式(気相) W

. f ; ‑ l   ‑

W1

k . 

^{̲ 1} 久 ‑ P，.

ω  ω=FS+ny r J  Jr

^{円 円}

k

_L

( B ‑ 1 4  ) 

116