液中浮上時の気泡挙動の数値解析

Synopsis:

Numerical simulation of bubble rising in reduced pressure vessel was carried out by volume of fluid method (VOF) solver of OpenFOAM, open source CFD toolkit. The result was validated against the experimental data of Kumagai et al. [Tetsu-to-Hagané, 101(2015), 93]. Behaviors of continuous rising bubble at 101.3 kPa, 18.2 kPa, 4.5 kPa were almost consistent with the experiment. The shape of bubble near injection nozzle at 101.3 kPa was more ellipsoidal than that of experiment and the mesh refining decreased the difference. Single bubble rising behavior was also studied. Rising rate and outside appearance were almost consistent, but slight continuous dispersion of gas phase from tail end of bubble was observed in only simulation and it was considered the issue of gas-liquid interface modeling. The cross section of bubble was spherical cap and consistent with the diagram of Grace [Trans. Inst. Chem. Eng., 51(1973), 116]

1. 緒言

製鋼精錬プロセスにおいては，転炉，取鍋における混合促進および RH における還流の形成等を目的として，溶鉄/溶鋼中へのガス吹きこみが盛んに行われている。ここでは毎分数m³ものガスが種々の形状のノズルから吹きこまれ，流動制御や溶鋼混合に対するガスの利用効率改善および吹込み羽口まわりの耐火物寿命の改善等が課題となっている。液相へのガス吹込みでは，ノズル形状，ガス流量，流体物性等によって生成気泡径や浮上における形状，速度に種々の状況が発生することが知られている。鉄鋼プロセスにおいても目的毎に様々なサイズや形状の羽口が使用されているが，そこに吹き込まれた気泡と溶鋼の挙動の直接的な実測や観察は困難であり，冷間試験による検討が多く行われてきた。

従来の検討において，まず，単一気泡の浮上挙動に関する検討があり，例えば，Grace¹⁾による気泡サイズ，浮上速度とその形状を分類，整理したダイアグラムは広く知られているものである。また，浴槽へのガス吹き込みに対し，ホールドアップに対するノズル特性の影響の定式化が，水あるいは溶融金属を用いて進められてきている^2,3)。

これらの実験的，理論的検討に対し，数値解析の試みも混相流解析手法の発展と共になされてきている。比較的計算負荷の軽い2流体モデルは，気泡径を既定値として与え気液間の運動量交換量を構成方程式により定義し気相挙動を有限個数の気泡の運動として表現するもので，大領域の解析，例えば浴槽全体の循環流の検討等に適用されている⁴⁾。これに対し，気液の界面変形挙動の捕捉を可能とした1流体モデルにおいては，吹きこまれた気相による気泡形成や浮上気泡の変形挙動の解析に適用されてきている⁵⁾が，気液界面の微細な形状を捉えるために十分な格子解像度が必要となるため，複雑な混相流挙動やプロセス全体のような大領域に対する解析では非常に多くの格子が必要となり，計算機における記憶容量や計算速度の点から分散メモリ型の並列計算に対応した解析手法が必須となる。近年，その利用が広まりつつある解析ソフトであるOpenFOAM⁶⁾はGPLにより公開されているフリーソフトウェアで，数値解析に有用な種々の機能のライブラリで構成されるCFDプラットフォームであるが，同ライブラリを使用する多様な流体解析ソルバも付属しており，並列計算に対応した1流体混相流ソルバも含まれている。

鉄鋼プロセスでのガスの吹込みによる混相流挙動の検討として，羽口形状とその効率的な配置や損耗現象に対する流動挙動の影響検討，さらには脱ガス等の気液間反応に至るまで，大領域における気液界面挙動の詳細な予測技術の重要性は高まっている。しかし，様々な気相挙動，特に溶鋼中の静水圧分布から予想される浮上時の気相体積変化まで考慮した数値解析手法の検討例は多くは見られずその予測精度について検討の余地が大きい。Kumagaiら⁷⁾は水モデル実験において，出口圧力を常圧(101.3 kPa)から減圧(数kPa)までの数水準として液相中を浮上する気泡の詳細観察を行っている。水モデルにおいて静水圧は溶鋼の約1/7となるため，溶鋼と同等の体積変化を観察するためには単独の実験装置では7倍の高さが必要となるが，Kumagaiらの検討においては実験の背景圧力を数水準とすることでそれぞれに相当する静水圧高さ位置の状況を再現することが意図されており，これに対し各圧力条件での気泡挙動の数値解析精度を検討することで，大きな圧力変化が伴う状況での気泡挙動への同解析手法の適用性を確認できるものと思われる。これより，本報告では，OpenFOAMの一流体モデル混相流ソルバを用いてKumagaiら⁷⁾により行われた減圧容器内を上昇する気泡の観察実験の解析を行い，同実験での高速ビデオによる詳細な観察に対して気泡挙動の再現状況を検討した結果を示す。

2. 検討内容

Kumagaiら⁷⁾の連続底吹きおよび単一気泡浮上の2種類の実験に対し，以下により数値解析を行い，観察結果との比較を行った。

2・1　連続底吹きによる気泡挙動

連続底吹きによる気泡挙動に関し，Kumagaiら⁷⁾による報告から以下に示す実験条件に対する数値解析を行い，観察との比較を行った。本検討で解析を行った実験条件をFig.1とTable 1に示す。同報告では，ガス流量は共通(2.5 cm³/s)として容器圧力を101.3 kPa，18.2 kPa，4.5 kPaとした実験が実施されており，本検討ではそれぞれの数値解析を実施し各条件での気泡挙動を比較した。解析にはOpenFOAM2.1.1の compressibleInterFoamを用いた。支配方程式を以下に示す。   

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \text{u}\right)=0$$(1)

  

$$\frac{\partial \left(\rho \text{u}\right)}{\partial t}+\nabla \cdot \left(\rho \text{uu}\right)=-\nabla P+\Theta +{\text{f}}_{st}+\rho \text{g}$$(2)

  

$$\frac{\partial F}{\partial t}+\nabla \cdot \left(F\text{u}\right)=0$$(3)

  

$$\rho ={\rho }_{l}F+{\rho }_g\left(1-F\right)$$(4)

  

$${\rho }_l=const.$$(5)

  

$${\rho }_g=P/\left(R{T}_0\right)$$(6)

Fig. 1.

 Experimental apparatus for continuos bubbling.

Table 1. Experimental specification for continuous bubbling.

Fluid	gas: air liquid: water
Temperature	22.0 ºC
Flow rate	2.5 cm3/s
Pressure	101.3 kPa, 18.2 kPa, 4.5 kPa
Photo rate	500 fps

ここで，

ρ:密度[kg.m⁻³]，u:流速[m/s]，P:圧力[Pa]，Θ:粘性応力[kg.m⁻².s⁻²]，f_sf:表面張力項[kg.m⁻².s⁻²]，g:重力加速度[m/s²]，F:液相率[−]，ρ_l,g:液相，気相密度[kg.m⁻³]，R:気体定数[J.kg⁻¹.K⁻¹]，T₀:ベース温度[K]

であり，質量，液相率の保存式を時間進展し，運動方程式をPISO法により解いて流速，圧力を更新する。同ソルバは，気相密度を圧力の関数とした等温Volume of Fluid法(VOF)による非定常混相流ソルバであり，ここでは乱流モデルは使用せず，運動量の移流スキームとしてSFCD scheme(Self-filtered central differencing scheme)⁸⁾を用いた。

解析格子および境界条件をFig.2とTable 2に示す。格子はOpenFOAM付属のblockMeshにより作成し，実験において代表的な気泡径が最も小さかった101.3 kPaでの実験での気泡サイズの1/10を目安として，解析負荷も考慮し高さ方向の格子サイズを0.5 mmとした。また，気泡のノズルからの離脱挙動の再現性向上を意図してノズル内配管およびその上流側の一部の導入配管を含めた形状とし，実験時にノズル上面で観察された良好な濡れ性から接触角として1度を設定した。計算はTable 3に示した計算機により，20並列前後での並列計算により解析を行った。

Fig. 2.

 Mesh (half area).

Table 2. Boundary conditions.

boundary	pressure	velocity	liquid ratio
inlet	zero gradient	fixed value	0
outlet	fixed value	zero gradient	zero gradient
nozzle surface	zero gradient	no slip	contact angle: 1 deg.
other	zero gradient	no slip	zero gradient

Table 3. Calculation system.

CPU	Intel (R) Xeon (R) CPU E5-2690 2 CPU/node
Memory	6 GB/node
Network	Infiniband QDR
OS	CentOS 6.3

2・2　単一気泡浮上挙動

単一気泡の浮上挙動に関係し，Kumagaiら⁷⁾による報告からいくつかの実験条件に対する数値解析を行い，観察結果と比較した。本検討で解析を行った実験条件をFig.3とTable 4に示す。解析手法は前項と同一とした。気相の初期状況に関し，実験ではダンピングカップ内に所定の体積の気相を保持した後，カップを回転して気泡を放出しその浮上挙動を観察するものとなっているが，解析ではカップの回転による気泡放出の再現は困難であったため，実験でのダンピングカップ内の保持状態，すなわち，ダンピングカップの一部を切り出した形状の気泡とし，その下面の静水圧を気泡内圧力として設定した。設定例をFig.4に示すが，後に示す解析結果では初期の界面挙動にこの初期条件に起因すると思われる差異が見られており，また，気相の初期形状によっては浮上開始後間もなく気泡が分裂し四散してしまう場合も見られ，気泡の初期配置は検討の余地があるものと思われる。

Fig. 3.

 Experimental apparatus for single bubble rising.

Fig. 4.

 Initial condition for single bubble rising.

Table 4. Experimental specification for single bubble rising.

Liquid (*)
catalog No.	KF-96-10cs
kinetic viscosity	10 mm2/s
specific gravity	0.935
surface tension	0.0201 N/m
Pressure	101.3 kPa, 1.3 kPa
Initial volume of bubble	2 cc

(*Shin-Etsu Silicone (Shin-Etsu Chemical Co., Ltd.))

3. 結果

3・1　連続底吹きによる気泡挙動

各条件での解析例を実験の観察例と合わせてFig.5に示す。実験においては，101.3 kPaでは数mmサイズの扁平な楕円球状の気泡がほぼ等しい間隔で生成し，容器圧力を下げ体積流量が増加するのに従い，気泡径の拡大，形状の不安定化，先行気泡との盛んな合体等が見られるようになるが，解析においても各条件にて同様の状況が比較的よく再現されることが示された。これより，比較的広範な条件における混相流状況を一定の精度で予測可能なものと思われる。

Fig. 5.

 Shape of continuous rising bubbles.

より詳細な気泡の生成状況として，101.3 kPaにおける気泡のノズルからの離脱直後の形状の推移を比較した結果をFig.6に示す。比較する気泡に対し，観察と解析の図の両側のそれぞれ同じ位置を矢印で示した。生成している気泡径(=生成頻度)はよく一致している。気泡形状としては，観察では真球の形状で離脱した気泡が上昇に伴い，半球，極度に扁平な楕円体へと変形していくが，解析では真球で離脱した後，球に近い楕円体に変形して安定しており観察されたほど扁平な形状とはなっていない。浮上速度に関し，矢印との相対位置で見て解析のほうが若干速い傾向があるが，これもその形状の違いによる流体抵抗の違いによるものと思われる。形状に差異が見られる原因として，観察における扁平楕円体の短径(2 mm以下)に対し，解析格子巾は0.5 mm刻みのため，外周部が4格子以下となってより微細な形状変化は再現されず，同部全体に界面張力が作用する結果，形状がより球体に近づいているものと思われる。ノズル直上近傍の格子巾をより小さくして計算した結果をFig.7に示す。格子数増のため短時間の計算しかできていないが，気泡形状がより扁平となり，観察結果に近づくことが示されている。

Fig. 6.

 Shape of continuous rising bubbles (101.3 kPa, near nozzle).

Fig. 7.

 Influence of mesh size to bubble shape.

3・2　単一気泡浮上挙動

0.1秒間隔での気泡推移の解析結果を実験結果と対比してFig.8に示す。浮上開始初期の0.1秒時点の解析結果では気泡の頂点部が突起状になっており，解析における初期条件が実験の気泡放出状況を再現しきれていないことが原因と思われるが，全体として観察とほぼ一致した浮上速度が得られた。なお，解析においては，浮上中の気泡端から尾を引くような形でわずかづつながら気相が継続して離脱し，気泡後流部に形成される円環状渦部に捉えられて気泡に追従する状況となっており，それに伴い気泡本体の大きさが減少している。より詳細な状況として，1.3 kPaでの浮上開始後0.3秒時点の気泡下端外周部近傍の液相率分布と流速ベクトルをFig.9に示す。上昇する気泡とその周囲の下降する液相とのすれ違い部に局所的に停滞域が生成しており，同部に気相がわずかづつ取り残されている様子が見られる。実験ではこのような気泡からの小気泡の離脱は見られておらず，解析上の課題と思われる。この原因としては，同部気液界面は最も曲率も大きく，同部の実際の気液界面に対し，解析では同部を数個の格子で捕獲する状況となって形状の再現性と安定性が不十分であるものと思われ，同部での液相率分布の急峻度が維持されにくい状況にあることが考えられる。この場合，同部の格子解像度を上げることで同部形状の捕獲の正確性と同部近傍の液相率分布の急峻度が向上し，また，本質的ではないが最小格子サイズが小さくなることで停滞域に取り残される気相量が減少する効果により解析精度の改善が期待されるが，解析負荷の増大も懸念され今後の課題と思われる。

Fig. 8.

 Shape of rizing single bubble.

Fig. 9.

 VOF and flow vector at bubble edge (1.3 kPa, 0.3 s, cross section).

次に定量的な確認として，浮上気泡の深さ位置に対する体積，圧力，およびそれより求めた常圧時体積の推移をFig.10に示す。ここでは，1.3 kPaの解析結果から，静止水面から5 mm以下の領域における全気相体積とその体積平均圧力を体積平均位置(深さ)に対して示した。浮上に伴う静水圧の減少により気泡の圧力は減少し体積が膨張するが，減圧条件下では圧力における静水圧の寄与の割合が大きいため，浮上に伴う変化も大きいものとなる。本圧力条件(1.3 kPa)下では，水面下−0.10 mからの0.09 mの浮上に対して，約880 Pa(=1000 kg/m³*9.81 m/s²*0.09 m)の静水圧の減少がおき体積は約1.6倍(=(1300 Pa+880 Pa)/1300 Pa)増加することが考えられるが，Fig.10(a)に見られる圧力および体積の変化量はこれとほぼ等しく，本解析結果において浮上の開始時と終了時とで気相量が概ね保存されていることがわかる。Fig.10(b)は気相量の保存性をより詳細に示すもので浮上過程全体において概ね保存されていることが示されているが，深さ−0.04 m付近より上方で体積の減少がやや顕著になり，最終的に5%程度の減少となっている。この原因として，浮上後半では数値解析上の問題として前述の気泡端からの気相の液相への分散が進んでおり，同状況での解析精度あるいは計算結果からの気相量の算出精度に課題があるものと思われる。なお，Fig.10(a)において浮上に伴って体積および圧力が振動しているが，これは本検討での初期気泡圧力設定値が静水圧とは吊り合ってなく膨張収縮を繰り返しながら浮上しているものと思われる。

Fig. 10.

 Pressure and volume of rising bubble.

Fig.11に気泡の断面形状の推移を外観と合わせて示す。浮上開始直後0.1秒時点の断面形状より，観察との比較で差異が認められた浮上開始直後の解析における頂点部の突起は，気泡底部がつき上がり気泡上面を押す状況となったためであることが示されており，その後，椀状となりその形状が拡大しつつ浮上してく様子が見られている。

Fig. 11.

 Shape of rising bubble (left: outside, right:cross section).

浮上する気泡の最終形状はGrace¹⁾によるダイアグラムでは，Sphere，Ellipsoidal，Spherical-capとに大別され，どの形状になるかは気泡レイノルズ数(Re_b=ρ_cd_eU/μ_c)およびエトベス数(Eo=gd²_e(ρ_c−ρ_d)/σ)により整理できるとされている。ここで，

ρ_c,d:密度，μ_c,g:粘度，d_e:球相当直径，U:終端速度，g:重力加速度，σ:界面張力，添字c:液相，d:気泡

である。

本検討条件および観察された気泡の上昇速度から気泡 レイノルズ数およびエトベス数を算出し(Re_b=469，E_o=120)Graceの分類に当てはめると，Fig.12に示すようにSpherical-capの領域となることから，解析結果により示された断面形状は妥当なものと思われる。

Fig. 12.

 Bubble shape calcurated in this study on the diagram by Grace.

4. 結言

OpenFOAM2.1.1のVOFソルバcompressibleInterFoamを用い，連続底吹きおよび単一気泡の2種の気泡浮上実験に対する数値解析を行った。底吹き連続気泡および単一気泡の浮上挙動に関し，生成気泡径，浮上速度等が各条件で概ね再現されることを示し，気泡断面形状が既往の検討と比して妥当であることを示した。また，実験では顕著には観察されない気泡端からの気相の離脱挙動が見られた。混相流解法における問題と思われ今後の課題である。

文献

• 1)   J.R.  Grace: Trans. Inst. Chem. Eng., 51(1973), 116.
• 2)   M.  Sano and  K.  Mori: Tetsu-to-Hagané, 68(1982), 2451.
• 3)   M.  Iguchi,  Y.  Demoto,  N.  Sugawara and  Z.  Morita: Tetsu-to-Hagané, 78(1992), 407.
• 4)   J.L.  Xia,  T.  Ahokainen and  L.  Holappa: Second International Conference on CFD in the Minerals and Process Industries, CSIRO, Melbourne, (1999), 187.
• 5)   M.S.  Annaland,  N.G.  Deen and  J.A.M.  Kuipers: Chem. Eng. Sci., 60(2005), 2999.
• 6)   O.  Nakamura,  T.  Yamamoto,  K.  Ueno and  K.  Takatani: Tetsu-to-Hagané, 101 (2015), 88.
• 7)   T.  Kumagai,  Y.  Kishimoto,  O.  Nakamura,  T.  Tagawa,  T.  Himeno and  S.  Shimasaki: Tetsu-to-Hagané, 101(2015), 93.
• 8)  OpenFOAM User Guide, OpenFOAM Foundation, http://www.openfoam.org/docs/user/fvSchemes.php#x20-1140004.4.5, (accessed 2014-05-30).