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

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

Kyushu University Institutional Repository

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

桑木, 賢也

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

https://doi.org/10.11501/3135115

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

権利関係：

R o a o R n o

一

O

﹃

n 0 3

仲

O

﹃

一 勺 一

ω R y o m

回 一

c o O

苔 コ の

﹃ g コポ‑一

o g z a

玄白 山町

@コ

S

さ 宣 言 ︒

o a ω R o s v ︑ 山 一

n 。

3 ・

. . .  

M 

ω 

明

。

ω 

の

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

平成九年度

桑 木 賢 也

論 文 目 録

主 論 文 1編 1冊

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

氏 名

ls 5 t l 甲

桑 木 賢 也

主 論 文 の 内 容 の 大 部 は 、 下 記 参 考 論 文 の 2'"'‑'5に公表済みであり、

さらに下記のように現在投稿中である。

Three‑dimensional oscillation of bubbly ftow in a vertical cylinderη  Internαtioηαl Journal 0] Multiphαse Flow 

(原稿枚数 30枚) (平成9年 6月投稿) Kenya Kuwagi and Hiroyuki Ozoe 

参 考 論 文 5編 1冊

1. 円筒容器内の底面中央領域より気泡を吹き込んだ場合の 二次元数値解析"

化 学 工 学 シ ン ポ ジ ウ ム シ リ ー ズ51化 学 工 学 に お け る 流 れ の 数 値 解 析 と 実 験 的 研 究 の 現 状 と 課 題，pp.121 ^r‑...J  126  (平成8年 2月)

桑 木 賢 也 、 平 野 博 之 、 尾 添 紘 之

2. "Numerical and experimental analyses for two‑phase ftow in a vertical  cylinder" 

Proceedings 0] the Second Kyushu ‑_~αipei Internαtionα1 Congress  oη Chemical Engineering， 1997

， 

Feb.  13 and 14  (平成9年 2月)

Kenya Kuwagi and Hiroyuki Ozoe 

3. "Three‑dimensional analyses for unsteady two‑phase flow in a vertical  cylinder" 

Proceedings 0] the Internαtionα1 Con]erence on Fluidαnd Thermαl  Energy Conversion '97

， 

July 21 ‑24  (平成9年 7月)

Kenya K uwagi and Hiroyuki Ozoe 

4. 高次風上差分法を用いた円筒容器内の気液二相 流 解 析 (二相流モデルの比較と差分近似法の比較)" 

化学工学論文集?第23巻，第6号， pp.861 ^r‑v 869  (平成 9年 11月) 桑 木 賢 也 、 尾 添 紘 之

5. N umerical analyses for two‑phase flow in a vertical cylinder" 

The Reports 01Institute 01 Advanced Material Study! 

Kyushu Uηiversity

， 

Vol.ll 

， 

No.2ぅpp.141r‑v 145  (平成9年 12月) Sunao Kusu

， 

Kenya Kuwagi and Hiroyuki Ozoe 

論文題目 竪型円筒容器内の気液二相流解析

氏 名 桑 木 賢 也

論 文 内 容 の 要 旨

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

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

第二章では、 二次元の気液二相流に対して、これまで解析例の少ない、多次元問題における二流体モデルへ の高次風上差分法の適用を試みた。数値解析モデルとして、分散流モデルと一般的に用いられるー圧力モデルを用 い、比較を行った。また、高次風上差分法として4種類のものを用い、これらの比較も併せて行った。その結果、

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

第三章では、本解析において液相として用いている水道水の汚れを考慮し、自由表面における液相の境界条 件に関する検討を行った。境界条件として、スリップ条件とノンスリップ条件を課した二つの場合の数値解析を 行い、可視化実験結果およびレーザードップラー涜速計による測定結果との比較を行った。その結果、上側表面 側壁付近において、可視化写真に見られる二次渦が、スリ ップ条件を用いたものでは見られなかったのに対して、

ノンスリップ条件では得ることができた。また、その付近の詳細な解析の結果、ノンスリップ条件とした場合、 二 次渦の領域について実験結果と良い一致が得られた。さらに三次元解析から、実験で得られるような高周波変動

がノンスリップ条件のとき得られた。このことから、水道水を用いた本解析条件においては、ノンスリップ条件 が厳密ではないものの、有効であると考えられる。

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

その際、揚力係数 CLに関する検討も行い、 CL= 0.2のとき実験値に近い結果を与えることを示した。また、鉛 直方向速度の半径分布が、中心付近で周りより値が小さくなる凹型分布となった。この分布から、上昇流が、中 心軸を通る半径方向に振動するものと、中心軸周りで周方向に振動するこつの振動パターンを考え、数値解析結 果の詳細な検討において、この二つの振動パターンを見いだした。更に、実験で得られる速度変動を、非定常数 値解析により得られた流動パターンにより説明を行った。

第 五 章 で は 、 四 章 の ア ス ペ ク ト 比 :A = 2.5に対して、水深の浅い A= 1.0の場合の解析を行った。これ は、 A= 2.5の場合と比較して、流動構造がより単純で、その詳細な解明には好都合と考えたことによる。数値 解析には、二章の結果から分散流モデル、四章の結果から液表面をノンスリップ条件とした。まず、二次元解析に より格子数の影響を調べた。その結果、粗い格子の場合、主流渦の大きさが可視化写真のそれより小さくなった。

これは、 二章の結果と同様な傾向であった。一方、格子数 96x 80以上のとき、主流渦の大きさはほぼ一定であ ることから、格子数として 96x 80程度必要であることが分かった。しかしながら、主流渦の下側の流れが、可 視化写真では斜めに下降するのに対して、計算では中心方向に水平なものであった。これは、三次元計算の結果か ら、 三次元的な流れ、つまり周方向の流れが生じ、その流れが下降流を生じているためと思われる。水深の浅い 場合、気泡プルームの揺動性は小さく、流れは二次元軸対称に近いと考えられがちである。しかし、解析の結果、

液相の流れは複雑で三次元的なものであった。その中で、数値解析により得られた流線から、気泡運動により生 じる中心付近の上昇流の周りに、二つの鉛直方向に軸をもつら旋涜れが存在すると予測された。このようなら旋 流れは従来報告されたことがないものと思われる。このら旋疏れは、その数や位置は時間とともに変化していた。

今後、このら旋流れは実験的にその存在が検証される必要がある。このような流動構造が、気泡プルームの揺動 現象に影響を与えているものと思われ、更に深いあるいは大きい容器内の揺動現象の解明の手がかりになるもの と期待される。

第六章では、第二章から第五章までの研究により得られた知見をもとに、その総括を行なった。

目 次

表目次 ^V 

図目次 ^Vl 

記号表 ^X 

第 1章 序 論 ¹ 

1.1  ^研究目的 ‑ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ a・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ‑・ ‑・‑・ 2  1.2  ^{気液二相流の数値解析方法} ・ ・ ・ ・ ・ ・ ・ ， ・ ・ ・・・ ー 3  1.2.1  ^{気相を Euler的に扱う方法 (Eulerianapproach)  .}  3  1.2.2  ^{気相を Lagrange的に扱う方法 (Particletracking method)  4}  1.2.3  ^{界面追跡法 (Interfacetracking} method) . .  ⁴ 

1.3  ^{既往の研究 4} 

1.4  論文の構成 6 

実験

2.2.1  ^{実 験 装 置 .}

︒︒ ハ可

υハ_Y

A u

d

第 2章 二 相 流 モ デ ル の 比 較 と 高 次 風 上 差 分 法 の 比 較 2.1 

2.2 

緒言

2.3 

2.4 

2.5 

目 次

2.2.2  可 視 化 実 験 結 果 . 10 

理 論 解 析 • • • • • • . • • • • . • • • • • • • • . • • • • • • • • • • • •• 11  2.3.1  仮 定 .• • • • • • • • • • . • • • • • • • • • • • • • • • • • •• 11  2.3.2 

2.3.3  2.3.4  2.3.5  解 析 結 果 2.4.1  2.4.2  2.4.3  2.4.4 

連 続 の 式 .• • • • • • • . • • • • • • • • • • • • • • • • • • •• 11  運 動 量 の 式 • • • • • • • • • • • • • • • • • • • • • • • • • •• 11  気 液 界 面 聞 に 働 く 力 の 構 成 式 • • • • • • • • • • • • • • • ^.，  12  計 算 方 法 .• • • • • • • • • • • • • • • • • • • • • • • • • • .，  14  14  ー圧 力 モ デ ル と 分 散 流 モ デ ル • • • • • • • • • • • • • • • •• 14  移 流 項 の 差 分 法 に よ る 違 い .• • . • • • • • • • • • • . • • .• 15  ボ イ ド 率 と 気 相 速 度 分 布 .• • • • • • • • • • • • • • • • • ^.，  16  レ ー ザ ー ド ッ プ ラ 一 流 速 計 に よ る 実 験 結 果 と の 比 較 .• . •• 17  18  結 言 .

第 3章 自 由 表 面 の 取 り 扱 い に つ い て 30 

3.1  緒 言 . 3.2 

3.3 

3.4 

3.5 

解 析 モ デ ル

数 値 解 析 方 法 .• . •

一 一 ー ・

・ ・ ・ ・ ・ ・ ・ ・ ・ 一 一

3.3.1  ^{仮 定 一} 一 ー ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ー 一 一

3.3.2  3.3.3  3.3.4  3.3.5  解 析 結 果 3.4.1  3.4.2  結言

・ ・ ・ ・ ・ 一 一

連 続 の 式

運 動 量 の 式 ( 分 散 流 モ デ ル ) 気 液 界 面 聞 に 働 く 力 の 構 成 式 計 算 方 法

一 白

一 ・ ・ ・ ・ ・ ・ ・ 一 一一 一

・ ー 一 一

表 面 近 傍 に 生 じ る 二 次 渦 に 関 し て . 三次 元 計 算 に よ る 比 較

31  32  32  32  32  33  33  36  36  36  37  38 

目 次

第 4章 三 次 元 解 析 48 

第 5章 4.1  4.2  4.3 

4.4 

49 

緒舌_{ロ ・}

実 験 .• • • • . • • • • • • • • . • • • • • • • • • • . • • . • • • • • • .• 49 

理 論 解 析 • • • • • • • • . • • • • • • • • • • . • • • • • • • • • • • • ^.，  50  4.3.1  仮 定 .• • • • • • • • • . • • • • • • • • • • . • • • • . • • ^.，  50  4.3.2 

4.3.3  4.3.4  4.3.5  解 析 結 果 4.4.1  4.4.2  4.4.3  4.4.4  4.4.5 

連 続 の 式 .• • • • • • • • • • • • • • • • • • • • • • • • • • •• 50  運 動 量 の 式 . • • • • • • • • • • • • • • • . • • • • • • • • •• 51  気 液 界 面 間 に 働 く 力 の 構 成 式 • • • • • • • • • • • • • • • •• 51  計 算 手 法 .. • • • • . • • • • • • • • • • • • • . • • • • • • ^.，  53  54  凹 型 鉛 直 方 向 速 度 分 布 の 構 造 • . . • . • • • • • • • . • • ^.，  54  揚 力 係 数 の 影 響 .• • • • • • • • • • • . • • . • • • • • • • .• 55  測 定 高 さ に よ る 鉛 直 方 向 速 度 の 半 径 方 向 分 布 の 違 い .• . .，  56  三 次 元 非 定 常 特 性 .• • • • • • . • • • • • • • • • • • • • • ^.，  56  気 泡 運 動 .• • • • • . • • . • • . • • . • • • • • • • • • • • •• 57 

4.5  結 言 . 58 

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

5.1  は じ め に 84 

5.2  実 験 .• 84 

5.3 

5.2.1  解 析 モ デ ル . • • • • . • • • • • • • . . • . • • • • . • • . . ，   84  5.2.2  実 験 方 法 .. • • • • • • • • • • . • • • • • • • • • • • • • • .，  85  85  理 論 解 析

5.3.1  5.3.2  5.3.3  5.3.4  5.3.5  5.3.6 

85 

仮 定 .•

連 続 の 式 .. • • • • • • • • • • • • • • • • • • • • • • • • . .，  86  運 動 量 の 式 ( 分 散 流 モ デ ル ) • • . • • . • . . • • • • • • • • 86  気 液 界 面 問 に 働 く 力 の 構 成 式 • • • • • • • • • • • • • • • .• 86  渦 度 ベ ク ト ル と 速 度 ベ ク ト ル ポ テ ン シ ャ ル .• • • • . • • ^.，  87  計 算 手 法 .• • • • . • • • • . • • . • • • • • • • • • • . • • ^{. ，}   88 

5.4 

5.5 

解 析 結 果

5.4.1  格 子 分 割 数 の 影 響 5.4.2  三次 元 解 析

5.4.3  気 相 流 動 の 解 析 5.4.4  非 定 常 特 性

:K士三き 耶ロロ

目 次

88  88  89  90  90  91 

第 6章 総 括 106 

付 録

A  無 次 元 化 109 

A.1  基 礎 方 程 式 の 無 次 元 化 .• . • • . • • • • • • • • • • • • • • • • • • • • • 109  A.2  物 性 値 お よ び 無 次 元 数 .• • • • • • • • • • • • . • • • . • • • • • . • • • 113 

B  二流体モテソレにおける HS‑MAC法 115 

B.1  基礎方程式の変形 • • • . • • • • • • • • • • • • • • • • • • • • • • . • • 115 

B.2  HS‑MAC法 の 手 順 . 117 

C  高次風上差分法の離散式 123 

C.1  高 次 風 上 差 分 法 の 離 散 化 式 • • • • • . • • • • • • • • • • • • • • • • • • 123  C.2  高 次 風 上 差 分 法 の 導 出 .• • • • • • • • • • • • . • • • • • • • • • • • • • 124  C.2.1  二次 精 度 風 上 差 分 法 • • • • • • • • • • • • • • • • • • • • • • 124  C.2.2  UTOPIA scheme (三次精度風上差分法) • • . • • • • • • • • 124  C.2.3  QUICK scheme  . . . 125  C.2.4  四 次 精 度 中 心 差 分 法 • • • • • • • • • • • • • • • • • • • • • • 127  C.2.5  四 次 の 拡 散 項 .• . • • • • • • • • • • • • • • • • • • • • • • • 127  C.2.6  Kawamura‑Kuwahara (K‑K) scheme (三次精度風上差分法) 128 

lV 

目;欠

参 考 文 献 130 

謝 百字 142 

V 

表 目 次

2‑1  Compl巾 daverage liquid velocity : liJlωε[m/s]  19  2‑2  Summary of the computational schemes  19  4‑1  Calculation conditions  53  4‑2  Boundary condition .  53  5‑1  Boundary condition .  92  5‑2  Calculation condi tions  92  A‑1 Parameters to calculate at t = 20 [OC]  113 

図 目 次

2‑1  Problem schematic  20  2‑2  Visualized flow pattern  21  2‑3  Computed grid system  22  2‑4  Flow patterns of liquid phase 

( d

B = 2.0 [mm ，]Grid : 24 x 60)  23  2‑5  Flow patterns of liquid phase 

( d

B = 1.O[mm ，]Grid : 24 x 60)  24  2‑6  Flow patterns of liquid phase 

( d

B^二 2.0[mm，]Grid : 96 x 150)  25  2‑7  Void fraction distribution with a dispersed flow model 

( d

B = 1.0 [mm]).  26  2‑8  Velocity vectors of gas phase with a dispersed flow model 

( d

B = 1.0 [mm]).  27  2‑9  Volumetric丑uxvecもorsof gas phase with a dispersed flow model 

( d

B = 1.0 [mm]).  28  2‑10 Vertical velocity distribution of liquid phase at z=0.17[m] with a dispersed flow 

model (d_B = 2.0 [mm])  29  3‑1  Problem schematic.  39  3‑2  Visualized flow pattern of liquid phase compared on the effect of boundary con‑

ditions for a top surface (A二 1.0)... . . ..  40  3‑3  Magni五cationof an upper and side corner.  . • • • • • . • • . • • • • • • • • . •• 41  3‑4  Radial distribution of verticalliquid velocity at z = 0.04 [m]  (z/h = 0.5).  . . ..  42  3‑5  Vertical disもrib山 onof radial liquid velocity at r二 0.04[m] (r/rout = 0.5).  . .. 43  3‑6  Visualized flow pattern of liquid phase compared on the effect of boundary con‑

ditions for a top surface (A二 2.5). ... . . . ..  44  3‑7  Transient velocity vectors in a vertical section with three‑dimensional model and 

slip condition for a top surface. ...  45  3‑8  Transient velocity vectors in a vertical section with three‑dimensional model and 

non‑slip condition for a top surface.  . . . .，  46 

3‑9  Transient ve出 calliquid velocity component at the point (rニ 0

，

zニ 0.17[m]).  4‑1  Problem schematic.  •

4‑2  Typical asymmet'1ic flow pattern of liquid phasc in a vertical section. 

4‑3  Typical asymmetric flow pattern of liquid phase in a horizontal section at z = 

0.16 [m]. ...  61  4‑4  Computed transient vertical velocity component of liquid phase at the point (1' 

ニ O^ぅZニ 0.17[m]).  . . . 

4‑5  Radial distribution of verticalliquid velocity along z = 0.17 [m] (Case 1).  . . ..  63  4‑6  Average velooty against sampling number. . . ..  64  4‑7  Transient liquid velocity vectors in a vertical section.  ...  65  4‑8  Transient distribution of the magnitude of verticalliquid velocity component at 

z 

= 

0.17 [m].  ...  66  4‑9  Trajectory of the position of maximum vertical velocity between tニ 22.9^I"'J 26.2 

[s].  . . 

4‑10 Schematic of the oscillatory flow paもtern.

4‑11 Transient responses of vertical liquid velocity measured with LDV at  z = 0.17 

m  _. .  . . 

4‑12  Standard deviation of vertical velocity along z = 0.17 [m] (Case 1) 

4リ R adialdistribution of vertical liquid velocity along z 

= 

0.17 [m] (Cases 2 to 5)  4‑14 Standard deviation of vertical velocity along z 

= 

0.17 [m] (Cases 2 to 5).  4‑15 Transient responses of verticalliquid velocity at z = 0.17 [m] (Case 4).  4‑16 Vertical distribution of radial liquid velocity along r二 0.04[m] (Cases 4)  4‑17 Radial distribution of vertical liquid velocity. . 

4‑18  Computed transient characteristics of liquid velocity and void fraction. 

4‑19  Computed transient velocity vectors in a verもicalsection and a horizontal section  at z二 0.17[m]. ...  78  4‑20  Computed transient vertical velocity component of liquid flow developed at the 

point (1' 

= 

0^うz

= 

0.17 [m]) 

4‑21 Transient velocity vectors of liquid flow developcd in  a vertical section. 

Vlll 

図 目 次 図 目 次

47  4‑22 Transient velocity vectors of liquid flow developed in a horizontal section at z 

= 

59  F

D 

e 

'EEEA 

10

 

.h

u 

u 

'h

u 

p十

︒

ム

円bρv 

VA 

1 J O  

m d  

f l l e   7・J h l a  

‑ ﹃

iri 

o T

Jq   

qL  

4 

81  82  60 

5‑1  Problem schematic.  93 

62 

5‑2  The effect of grid numbers on the stream lines of liquid phase computed with  two‑dimensional model.  •

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

94  95  5‑4  Transient distribution of circumferential vector potential component in a vertical 

section and velocity vectors in the horizontal section at z/h = 0.5.  . . . ..  96  5‑5  Velocity vectors in 

e  ‑

z section along r/rout 

= 

0.75. . . ..  97  5‑6 Time averaged velocity distribution of liquid phase.  . . . ..  98  5‑7 Stream lines of gas phase from the botもomto top surface at t = 36.0 [s].  ....  99  5‑8  Stream lines from 30 x 30 points in a vertical section (ムt

= 

1.0 [s]).  ... 100  5‑9  Stream lines  from 30 x 30 points in  the horizontal section  at  z  0.04  [m]  68  (L'lt = 1.0 [8])  ¹⁰¹ 

Q u n u

‑

‑ q L q δ A

吐 戸

︒ ウ

t

pO

門i円i円iウtウt勺Iウt

5・10Stream lines from 30 x 30 points in various section at t 

= 

36.0 [s]  (ムt

= 

1.0 [s]).  102  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‑12  Stream lines  around central upward flow  region from the view point on the 

bottom side in Fig.5‑8 (d) at t = 39.3 [s].  104  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  B‑1  Grid definition for discretization .  118  C‑1  Schematic of the approximation forゆRon QUICK scheme.  126 

79  80 

lX 

V 

︑₁

1/

L vdV 

A Y U

‑

‑ 司 E A E

F v

‑

O G  

1 d v e  

げV

&

E U

J

=

ヴ は

a︐︑

. 1

日・

1 l i e

‑ u M J h α   e・I

e o r c v c   y p向

︒

lt.n副

a a

. α

汀m

・ 児

‑K

O C 3

i

t t  

︑

r L A F i v i

e

‑ H

いh e e

v e

‑s

v v  

一一

一一

一一

一一

一一

[m/s]  [m/s]  [m/s]  [m/s]  [m] 

v  Vr 

w 

記号 表

^Z 

Greek letters 

本論文で用いた主な記号を以下に示す。

A  = aspect ratio (ニ h/rout) ‑[] 

CD 

^{= drag coe伍cient ‑] [}

C

_L  ^{= 1ift force coe伍cient ‑[]} 

C

_{V M} ^二 virtualmass coe伍cient ‑] [ d_B  ^{二 diameterof a bubble  [m]}  Eo  = Eるtvosnumber (= g(PL ‑ρc)d1/σ) 

‑ [

]  F_D  ^{= drag force per unit mass}  [N/m3]  FL  ^{= lift force per unit mass}  [N/m3]  Fv  ^{= shear force per unit mass}  [N/m3]  Fv M = virtual mass force per unit mass  [N/m3]  g  = gravitational acceleration  [m/s2] 

h  = heighL of liquid surface  [m]  Jc  ^{= gas volumeLric fiux ( =α}CVc)  [m2 / 5]  M  ^{ニ Mortonnumber (=} 9μi(PL ‑pc)/ PLσ3)  ‑[

] 

p  二 pressure [Pa] 

Q 

^{= injected gas fl̲ow rate  [m3/s]} 

ア = radial coordinate  [m] 

アc = radius of gas injection  [m]  Re  = Reynolds number (= P Ldvγ/μL)  ^[‑]  rout  ^{ニ radiusof vessel  [m]} 

t  = bme  [sec] 

u  = radial velocity  [m/s] 

α  ^ニ voidfraction (gas hold‑up) 

。

^{ニ circumferentialcoordinate} 

ニ viscosity

ρ  = density 

σ  = surface tension 

T  = stress tensor 

ψ 

ニ velocityvector potential  ぃy = vorticity vector 

守リ

卦 ・

4 J 川 川 門 川 叫

J

川 門 司

J m R h d

札 川 れ

; m ν

Subscripts 

c  = continuous phase  d  = dispersed phase  G  = gas phase 

k  = gas phase (k = G) or liquid phase (k = L)  L  = liquid phase 

r  = r‑component  z  = z‑component  O ニt1‑component Operator 

D / Dt 

=  8 /  

8t 

+ 

V k . ¥1 

第 1 章

序論

1.序論

1 . 1   研究目的

気液二相流は、多くの工業プロセスにおいて問題となっており、伝熱問題においてもボイ ド率分布の予測が必要となっている。気液二相流が実際に問題となる場として以下のようなも のがあげられる。

‑化学プラント(気泡塔、精溜塔、反応器等)

‑電気化学工業装置(電気分解)

‑鉄鋼精錬(ガス吹き込みによる溶鉄の撹枠)

・発電所(ボイラー内での水の沸騰、原子炉内の燃料棒や伝熱管付近で発生する気泡流)

‑酒造(発酵)

・船舶涜体(気泡吐出による摩擦抵抗軽減問題) .活性汚泥処理の曝気槽(気泡噴流)

気液二相流の輸送現象は、ボイド率の値などによってその流動様相は全く異なったものに なるなど、流れが複雑となり、測定ならびに解析が困難なために、今日においてもいまだ明ら かにされていない点が多い。

まず測定においては、 二つの相が混在するため、 二相の同時測定が困難である。また、同 時測定ができる場合も液相、気相のデータの分離などの問題がある。更にレーザードップラ一 流速計のような光学系測定器では、気泡による屈折などにより、信頼性の高いデータが得にく

いといった問題もある。

一方、数値解析においても、気液二相流の流動現象を、その細部にわたるまで決定する基 礎方程式を、実際の気液二相流において解くことは現在のところ不可能に近い。そこで、実用 的には、気液二相流を時間的、空間的に平均して取り扱うモデリングが必要となる。そのモデ リングにおいては数多くのモデルが提案されているが、確固たるものはなく、検討されている 段階である。実際的な問題としては、 二相流現象自体、気泡ブルームなどの揺動現象に見られ るように非定常性が強いため、数値解析的に不安定になりやすく、解を得ることが困難となる ことが多い。

これまで、気液二相流問題に対しては、多くの実験的研究がなされているが、これらの多 くが数値解析との比較を目的としていないため、数値解析結果との比較は難しい。また、多く

1序 論

の数値解析も行われているが、実験との比較は意外と少ない。そこで、本研究では、数値解析 と実験の両者を行い、多くの工業プロセスで取り扱われる円筒容器を解析対象とした。まず、

数値解析において高精度な解を得ることを目的とした。これは、 二相流モデルの妥当性の検証 には必要不可欠なためである。また、上述のように、確固たる二相流モデルが確立していない ため、それを用いてのこ相流流動の解明は十分なされていない。そこで、本研究では、とれま でほとんど行われていない、円筒容器内の三次元非定常特性の解明も試みることとした。

1 . 2   気液二相流の数値解析方法

1 . 2 . 1   気相を Euler 的に扱う方法 (Eulerianapproach) 

この手法は、分散している気相も連続相である液相と同様に時間的、空間的に平均化を 施し、基礎方程式を導出する方法である [1]0つまり、気相を準連続相として扱う方法である。

Euler的なモデリングにおいては、気液二相流を一つの混合物と見なして平均するものと、気 相、液相をそれぞれ別個の流体と見なして平均する [2][3]ものとがある。前者は混合物モデル、

後者は二流体モデルと呼ばれる。

(1)混合物モデル

‑均質流モデル

‑スリップ流モデル

‑ドリフトフラックスモデル [4]

混合物モデルは気相と液相の速度差をどのように取り扱うかによって、均質流モデル、ス リップ流モデル及びドリフトフラックスモデルの三つに分けられる。均質流モデルは、各相の 速度差がないものとして記述される。スリップ流モデルは、相問の速度比(あるいは速度差)に より記述され、 ドリフトフラックスモデルは、混合体の流束との差を各相のドリフト速度と定 義し、これを用いて基礎式が記述される。

(2)二流 体 モ デ ル 山

二流体モデルは、各相の圧力の取り扱いによって、 一圧力モデルと二圧力モデル [5]に分 けられる。また、準連続相からさらに分散性を考慮した分散流モデル [6]も提案されている。 二

1.序論

流体モデルは、前述の混合物モデルと比較して、詳細な解析が可能である。また、後述の二つ の解析方法がその解析手法に様々な問題点が残っており、 概して非常に高いコンビュータの能 力を必要とするととから、現在のところ最も実用的な解析方法と思われる。

1 . 2 . 2   気相を Lagrange 的に扱う方法 ( P a r t i c l et r a c k i n g  method) 

この方法は、分散相(気泡流の場合、気相)を一つ一つ追跡していく方法で、主に囲気二 相流 (粉体工学)や固液二相流の分野で使われてきている。との手法の長所としては、分散相 の相互の影響(合体、分裂など)を考慮しやすいことである。また、粒子一つ一つを追跡して いくため、粒子数あるいは体積率の数値拡散の問題はほとんどない。しかしながら、分散相の 体積の影響や連続相への影響が大きい場合、解析が複雑となることから、気液二相涜の分野で は現在まであまり使われていない。

1 . 2 . 3   界面追跡法 ( I n t e r f a c et r a c k i n g  method) 

この方法は、気液界面の形状も計算していくため、界面が変形する気液二相流問題では有 効な手法である。代表的な解析方法に VOF(Volume of Fluid)法 [7][8]がある。この方法で は、多くの格子点数を用いれば、 気相の流れにより誘起される乱れも求めることができ、気液 二相流モデルを用いず、直接計算することができる。しかし、気泡に最低、数セルを割り当て る必要があるため、気泡径と比較して容器の寸法が大きい場合は、特に三次元計算の場合、コ ンビュータの能力の制約を受ける。また、気液界面を解くため、界面における力の詳細な構成 式が求められる。

1 . 3   既往の研究

この節では、本論文作成に関わった既往の研究について述べる。

1.1節の研究目的でも述べたように、気液二相流は多くの工業フロセスで問題となるため、

様々な工学分野において取り扱われている。数値解析を行ったものを中心に例を挙げていくと、

例えば、原子力工学では、 Minatoらの管よりサブクール水を放出する場合をシミュレートした もの [9]や給水加熱器内の解析 [10]、熱交換器内の流路に生じる気液二相析しを取り扱った Murata らの解析例 [11]などがある。この分野では、特に原子力発電所において、厳しい安全性が求め られることから、多くの計算コードが開発されている。

4 

1 序 論

化学工学では、電気化学工業装置を想定した、垂直平板電極より電解気泡を発生させた場 合の平板回りの流動解析を行った松浦ら [12][13]のものがある。

機械工学においては、ポンプの二相流性能において問題となるキャビテーション気泡を取 り扱った松本ら [14][15] [16] [17]の解析がある。

船舶流体では、気泡吐出による船舶の摩擦抵抗軽減問題に関連して、 土井らが正方形ダク ト内気液二相流乱流に対して、 LES(Large Eddy Simulation)を用いた数値解析を行っている

[18] ⁰ 

土木工学では、躍層の破壊、波浪の阻止などで気泡噴流が応用されているが、北野らは、

活性汚泥処理の曝気槽での散気に関連して、気泡噴流の解析を行っている [19]。

また、鉄鋼精錬などの分野においては、ガス吹き込みによる溶鉄の撹はんに関して、数多 くの解析例 [20][21]  [22]がある。

酒造関連では、酒類の発酵槽内のアルコール発酵反応により生成された二酸化炭素ガス気 泡の上昇により引き起こされる流動を取り扱った、 Gofukuら[23卜児玉ら [24]の解析がある。

気液二相流は、上述のようにその解析の必要性から多くの研究が行われてきている。しか し、これまでの研究は実験的なものが中心で、上記のような数値解析が盛んに行われるように なったのは、近年になってからである。これは、計算機の能力の他に、単相流解析と比較して、

基礎方程式が確立していないためと考えられる。 Batchelorはこのような現状を、「二相流の進 展に対する主な障害は、数学的、数値的または実験的手法の貧困ではなく、本質的なことは、

分離した相の平均運動を支配する閉じた形の式がないことである。それは 1930年代に乱流理 論がおかれていた状況と似ているけと述べている

[ 2 5 ] [ 2 6 ] 0  

このような状況下で、数値解析は、まず最も簡単な均質流モデルを用いた解析やスリップ 流モデルを用いた解析 [27]から行われ、やがてドリフトフラックスを用いての解析が主流と なった。しかし、ドリフトフラックスモデルでは、実験定数であるドリフトフラックスパラメー タが必要であり、データベースの蓄積という問題から、 二流体モデルが注目されるようになっ た。しかし、 二流体モデルは、比較的容易に、しかも詳細に解析が可能である反面、数値的に 不安定に陥り易い。この原因に関して、 Lyczkowskiらは「気液二相流に対する基礎方程式(二 流体モデル)が初期値問題として数学的に適切でないJと報告している [28]^{。この問題に関し} て、 二流体モデルの数学的な研究も数多くなされており [29][30]  [31]  [32]  [33]、基礎方程式が 数学的に不適切である場合であっても、差分近似によって付加される数値拡散などの効果によ り得られた安定な数値解は、実際の二相流流動をかなりよく模擬し得ることが、木村ら [34]に より示されている。

5 

1.序論

これらと並行して、 二相流乱流の数値解析も多く試みられている。単相流乱涜において実 績のある、 k‑eモデルを用いた解析も行われており [35][36]  [37

卜

Ranadeは有限体積法を用い て解析を行っている [38]。また、 土井ら [18]や二瓶ら [40]は LESを用いて乱流解析を行って いる。ことで、気液二相流において乱流粘性というものを考えるとき、液単相に固有な成分と 気液問の相対運動または気泡の乱流拡散によって誘起される成分に分けて考えることができる。

これらの各成分に関しては、 Satoらが研究を行なっている [39]⁰ 

一方、気相を Euler的に取り扱う二流体モデルに対して、 Lagτange的に取り扱う方法(粒 子追跡法)も検討されていて、 Durstら [41]や、村井ら [42]により両者の比較がなされてい る。分散相を Lagrange的に取り扱う方法は、囲気二相流(粉体工学)や固液二相流で多く研 究がすすんでおり、平野は、この Lagrange的な手法による、 Brown運動をする微粒子の運動 方程式、いわゆる Langevin方程式を用いて、数値解析により流動場中の微粒子挙動を研究し ている [43]0これに対して、気液二相気泡流に対する粒子追跡法はほとんど進展していなかっ た。とれは、気泡が通常、粉体工学などで取り扱われる粒子と比較してその径が大きく、気泡 が液相に及ぼす効果と、液相が気泡に及ぼす効果が無視できないためである。この両方を考慮 する方法 (two‑waymethod)を取り入れた解析方法および解析結果は、富山らにより示されて いる [44][45]0 

また、気液二相流の解析の困難性の一つに、気液界面が変形することが挙げられる。この 問題に対して、自由表面の解析を目的に開発された

VOF

法

[ 7 ]

などを、気泡流に適用する試 みもなされている [46][47]  [48] [49]。

1 . 4  論文の構成

第一章では、代表的な既往の研究についてまとめ、本論文の意義を明らかにする。

第二章では、 二次元の気液二相流に対して、一圧力モデルと分散流モデルに高次風上差分 法を適用して数値解析し、その適用性に関する検討とともにこ相流モデルの比較を行う。さら に4種類の高次風上差分法の比較を行う。

第三章では、自由表面においてスリップ条件とノンスリップ条件を課した二つの場合の数 値解析を行い、可視化実験結果および流速測定結果との比較により、自由表面の境界条件の検 討を行う。

第四章では、円筒容器内の気液二相流問題に対して、 三次元数値解析と実験を行う。慣性 項には高次風上差分法の Kawamura‑Kuwahara法と QUICK法を用いた。まず時間平均した数

1序 論

値解析結果と実験との比較を行い、数値解析の妥当性を検証する。その上で、得られた結果か ら中心軸回りでの上昇流の振動パターンの検討を行う。

第五章では、四章のアスペクト比 :Aニ 2.5に対して、水深の浅い A

= 

1.0の場合の解析 を行うこととする。これは、 A

= 

2.5の場合と比較して、流動構造がより単純で、その詳細な 解明には好都合と考えたことによる。 三次元的な流れの構造に関するさらに詳細な検討を行う。

第六章では、第二章から第五章までの結果より、総括を行なう。

第 2 章

二相流モデルの比較と高次風上差分法の比較

8 

2 二相流モデルの比較と高次風上差分法の比較 2.二 相 流 モ デ ル の 比 較 と 高 次 風 上 差 分 法 の 比 較

2 . 1   緒 言

度の比較を行なっている。それによれば、高次風上差分法は一次風上差分法より精度の高い解

を得ることができ、調べた高次風上差分法の中では

QUICK

法が最も精度の高い解を得ること が出来ると報告している。 一方、多次元問題では解析的に解を得ることが出来ない。そこで、

実験結果との比較が必要となる。本研究では、解析対象として実際の工業プロセスの多くで取 り扱われる円筒容器内の気液二相流問題を取り上げ、水一空気系で行なった戸井らのものと同 じ条件 [57]で実験を行ない、数値解析結果と比較を行なった。

Fig.2‑1に実験のモデル図を示す。半径 Tout

= 

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

= 

1.67x10‑⁶[m³js]で、このときの r‑z断面の流動の観察を 行なった。

まず全体の流れを定性的に把握するためアルミ粉を液相に混入し、流動状況の可視化を行 なった。光源にはスライドプロジェクタ

M a s t e rH I ‑ L U X ‑ H R I O O O

を使用した。スリットを施 した光源からアクリル円筒にスリット光を照射し、中心軸を通る鉛直断面の流動を可視化した。

カメラは NikonF3、レンズはマイクロニッコール F2.8，55mmレンズを用い、シャッタース ピードは 1秒間開放で撮影した。

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

とした。

気液二相流は多くの工業プロセスにおいて出現する現象である。計算機の進歩に伴ない、

二相流も数値解析されることが多くなったが、いまだ多くの問題点を抱えている。特に二相流 に対するモデル方程式には確固たるものが確立されていないのが現状である。近年、 二相流の 解析に多く用いられているこ流体モデルは、気相、液相それぞれに対して運動方程式をたて、

解析を行なうため、従来のモデルより、より詳細な解析が可能であり、松浦らは、とのこ流体モ デルに乱流モデルを組み込み、垂直平板回りの乱流気液二相流に適用し、実験結果との良好な 一致を得ている [51]0 しかし、この二流体モデルにおいても、 一般的に用いられているー圧力 モデル(気相、液相の両相にかかる圧力を等しいとおく)は、初期値問題に対して数学的に不適 切であることが指摘されている [52]。不適切な解の特徴は、解が初期値に連続的に依存しない ことであり、 一階の偏微分方程式系を差分法で解く場合、数値的不安定が生じる。しかし、 一 圧力モデルがあらゆる条件において不適切になるわけではなく、圧力、ボイド率、気液相対速 度の値によって変化し、また粘性項や仮想質量力などは安定化に効果があるととが知られてい る [53]。既往の研究の計算においても、数値粘性などの効果により、安定に解けており、実際、

ほとんどの数値解析において数値粘性の大きい一次風上差分が用いられている。しかし、数値 粘性の少ない、物理的により信頼性の高い高精度の解が求められるべきである。そこで、高次 の差分法を用いることが考えられるが、二流体モデルに高次風上差分法を適用した例は、大川￨

ら[54][55]に見られる程度である。 二流体モデルを用いて数値計算をする場合には、高次風上 差分法は一次風上差分法よりも不安定なこと、空間刻み幅が小さい場合に数値的不安定性が発 生しやすいことが報告されている [56]。一方、気相、液相の圧力を別個に扱う二圧力モデル [5] や、分散相(気泡流の場合、気相)の圧力項を無視した分散流モデル [6]が提案されている。

本章では、 二次元の気液二相流を、 一圧力モデルと分散流モデルに高次風上差分法を適用 して数値解析し、その適用性に関する検討を行なった。

2.2.2 

可視化実験結果

アルミ粉法による液相流動の可視化実験の結果を Fig.2‑2に示す。図は r‑z断面の可視化 結果である。なお中心軸付近は気泡の動きが撮影されており、それに妨げられて液相の流動は 見ることができない。中心軸付近の気泡の上昇に伴い、液相は気泡から抗力を受け上昇し、そ の後自由表面に沿って容器側壁に向かい、容器側壁で下降流となって、結果的に渦(主流渦)を 形成している。この結果は計算結果の検証データとなる。

2.2  実験

2.2.1 

実験装置

高次風上差分法を二相流問題に適用した大川￨らの数値解析 [55]は一次元を対象としてい る。これは厳密解が存在するためであり、数値計算による解と厳密解との比較を行なうことに より、高次風上差分法の二相流問題への適用性あるいは数値解の移流項の差分近似法による精

2 二相流モデルの比較と高次風上差分法の比較

2 . 3   理 論 解 析

2 . 3 . 1  

仮 定

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

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

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

‑流れは軸対称の 2次元

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

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

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

2 . 3 . 2   連続の式

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

k 

= 

G (気相)， L (液相) A^内Ic 1θ 

一 三 +δt  一 一rθr ( 町 内，‑‑"‑"J )+τa (αk

同) = 

0  αL +αGニ 1

(2‑1) 

(2‑2) 

Eq.(2‑1 )を G と L について記し、両式の和をとり、そこに Eq.(2^・2)を考慮することによ り、非定常項が消去され、次式を得る [58]。

( ; t M  

^(2‑3) 

2 . 3 . 3   運動量 の式

i)一圧力モデル

基礎方程式として Durstのモデル [59]に仮想質量力と揚力を考慮したものを用いた。

11 

2二相流モデルの比較と高次風上差分法の比較

(液相)

αρL

竺 土

‑αL"Vp +αLFV+ F_D + FVM + FL^ー ρL9αL Dt 

(気相)

αcpc

竺三

^{‑ αc"Vp‑F}D ‑FVM ‑FL ‑ Pc9αG  Dt 

(2‑4) 

(2‑5)  Eq.(2‑5)中の右辺第 l項は圧力勾配による力を表しており、第2項 FD は気液界面閥抗力、

第3項 FVMは仮想質量力、第4項は揚力 _FL、第5項は重力を表している。 Eq.(2‑4)中の右 辺第 1項は圧力勾配による力、第 2項 Fvは粘性力、第 6項は重力で、第 3項から第 5項は、

Eq. (2‑5)の第2項から第 4項の力の反作用により生じる力である。

ii)分散流モデル

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

(液相)

L

竺主二一

^αL¥1p +αLFV + F D + F L ‑ P L9αL  Dt 

(気相)

G E

竺 =‑

F ^D ‑F ^L + (p L ‑Pc) 9αG  Dt 

(2‑6) 

(2‑7) 

分散流モデルの特徴としては、分散相(気相)の運動量式に圧力勾配項および粘性拡散項が 現れないことである。これは分散相に連続相と同様に圧力勾配によって力が働くとは考えにく

く、むしろ連続相を介して間接的に働くと考える方が物理的に妥当であることによる。また粘 性拡散項に関しても連続相と同様の速度勾配に比例した粘性応力を考えるのは極めて困難であ り、これも連続相を介した間接力と考える方が物理的にも自然であるためである。なお仮想質 量力はここでは考慮しない。

2 . 3 . 4   気液界面聞に働く力の構成式

i)仮想質量力

単位体積当たりの仮想質量力は以下のように表せられる

[ 6 0 ][ 6 1 ] 0   r 

^{Dvc  DVL} 1 

FVM 

=

^αGρ¥...Tr^[.L^JI ‑^CVM Y川 ~一一一 〉

LDt  Dt J  ^(2‑8) 

12