密閉容器内における粘性流体中を上昇する単一気泡に関する数値流体解析

Synopsis:

Numerical analyses have been carried out for a single rising bubble flow in a viscous fluid. It is assumed that both liquid and gas phases are the Newtonian fluid and viscous dissipation is neglected. The dimensionless governing equations for axisymmetric and three-dimensional coordinate systems have been successfully solved using the fractional step method combining with the level set method. The present numerical results exhibit fairly good agreement with experimental results in terms of the bubble shape and rising velocity.

1. 緒言

気泡流に代表される二相流の動力学は，気液間の物性値の大きな違い，その不連続性，また界面張力の存在などの理由により非常に高い複雑性を有している。重力場においては，気泡流の駆動力は密度差に基づく浮力が支配的であるものの微細な構造においては曲率が大きくなり，界面張力もまた支配的になる。ボイド率など種々の条件により様々な気泡流が考えられるが，それらの中で比較的単純とされる浮力を駆動源とする単一気泡の上昇運動を考える場合においてさえも，気泡径，密度差，粘度の影響，界面張力，重力加速度などの諸因子に関連し，全容が解明されているとは言い難い。

気泡流に代表される混相流の数値解析技術としては，大別するとVOF法¹⁾とレベルセット法²⁾とがある。VOF法は気相と液相を区別するために各計算セルのボイド率で表されるVOF関数を用い，ドナー・アクセプター法により，有限体積法に基づく移流計算を行うことにより体積保存性を満足させながら計算することができる。一方，レベルセット法では，気液界面からの距離関数(レベルセット関数)を用い，移流計算を行うが，移流計算後は再初期化と呼ばれる操作を通してレベルセット関数を修正する必要がある。また，計算が進むにつれ体積保存が破れてしまうことが知られる。しかしながら，界面曲率を精度よく計算できる利点を持つため，近年ではVOF法とレベルセット法の双方の利点を組み合わせ，欠点を補いあうCLSVOF法が提案されている³⁾。

鉄鋼分野においては，精錬プロセスにおいてアルゴンガスの気泡に働く浮力を利用して溶鋼の撹拌を通して溶鋼中の介在物の除去が行われている。それに関連して，2010年から3年間で実施された「精錬反応プロセスにおける混相流・多重スケール解析技術の開発」研究会において，混相流に関わる基礎研究が3つのワーキンググループに分かれ行われてきた。気泡挙動ワーキンググループの熊谷は，シリコーンオイルや重液を用いRHを模擬した数種類の実験を行ってきた。それらの実験の中でも特に現象を単純化した単一気泡の上昇に伴う膨張過程の実験的観察は，ベンチマーク的な有用な結果を提供している⁴⁾。その実験結果との比較については，気泡挙動ワーキンググループのNakamuraら⁵⁾，Himeno and Kumagai⁶⁾によってそれぞれ異なる視点からのCFD解析が行われている。筆者らにおいても，その実験結果を再現すべく二相流体解析コードの開発に従事してきた。本稿は，実験結果との定量的な比較あるいは現象を説明するために，独自の数値流体解析を行った。特に本稿においては，数値流体力学的な数式モデルの詳細を示すと共に，軸対称二次元計算と三次元計算の両方を検討し，結果の比較を行った。

2. 解析方法

2・1 解析モデルと仮定

液相にはシリコーンオイル(信越化学 信越シリコーンKF-96-10 cs)の物性値を，また気相には常温常圧の空気の物性値を用いる。密閉容器内の雰囲気圧力によらず，本研究では気相，液相とも非圧縮性のニュートン流体であり，等温場であると仮定する。減圧され雰囲気圧力が低い場合には，気泡が上昇する過程において気泡の膨張が想定されるが，本研究においては単純に液相の静水圧に関係づけて膨張する過程を模擬するものとした。Fig.1に解析モデルの初期状態を示す。軸対称解析モデル，三次元解析モデルのいずれにおいても，無次元で，直径5(三次元計算では一辺が5の正方形断面)で高さ10の密閉容器中において，直径1の球形気泡と同体積の回転楕円体の気泡が，中心軸上(三次元計算では容器中央)で高さZ=1の位置から放たれるものとする。回転楕円体の初期位置および形状については以下の数式で与えた。   

$$\frac{R^2}{A^2}+\frac{{\left(Z-1\right)}^2}{B^2}=1,\left(A=0.7,B=0.255102\right)$$(1)

Fig. 1.

 Schematic model at the initial state of computation.

ここで，Rは中心軸からの無次元半径座標，Zは無次元軸方向座標を表す。本稿では，実験と比較するために，初期の気泡体積が0.05 cm³から5 cm³までの4通りを対象とした。気泡の大きさに依存せずに数値解析の精度を一定基準に保つために，容器サイズを一定に保つのではなく気泡あたりの格子点数が一定になるように容器サイズの方を変更して計算を行った。

2・2 支配方程式

気相，液相ともに非圧縮性のニュートン流体であると仮定した。解析には，以下に示されるように，連続の式，移流方程式，およびNavier-Stokesの運動方程式を用いた。ただし，外力には重力のみを考慮した。総和規約を用いて以下のように示される。   

$$\frac{\partial u_j}{\partial x_j}=0$$(2)

  

$$\frac{\partial \rho }{\partial t}+u_j\frac{\partial \rho }{\partial x_j}=0$$(3)

  

$$\rho \left(\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}+f_i,f_i={}^t\left(0,0,-\rho g\right)$$(4)

  

$${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$(5)

式(5)は粘性応力テンソルを表す。ここで，密度ρや粘度μは，本来，気相か液相かに依存して，界面の両側において不連続に物性値が異なる。それぞれ，気相，液相の物性であることを示すために，必要に応じて下付き添え字G，Lをつけて区別するものとする。密度については，式(3)を解くことにより，密度の時空間分布を知ることができるが，本論文では後述するように，人工的に厚みを持たせた界面において密度や粘度などの物性値が連続的に変化する領域を設定する。ここまで，界面張力について触れてはいないが，現段階ではそれを無視すると，基本的には，式(2)，(3)，(4)，(5)を数値的に解けば，気液界面形状の時間発展を知ることができる。しかしながら本稿では，式(3)と相似な形式であるレベルセット関数φ(界面からの距離を表し，長さの次元を持つとする関数)に対する移流方程式を用いるものとする。   

$$\frac{\partial \varphi }{\partial t}+u_j\frac{\partial \varphi }{\partial x_j}=0$$(6)

式(3)を解く代わりに，式(6)を解いて，気液界面の位置(φ=0)を把握し，そこから逆算して，厚みを持たせた界面領域での密度や粘度を決めることができる。ここでは密度および粘度は，φを用いて次式で与える。   

$$\rho =\frac{1}{2}\left({\rho }_L+{\rho }_G\right)-\left({\rho }_L-{\rho }_G\right)H_{\alpha }\left(\varphi \right)$$(7)

  

$$\mu =\frac{1}{2}\left({\mu }_L+{\mu }_G\right)-\left({\mu }_L-{\mu }_G\right)H_{\alpha }\left(\varphi \right)$$(8)

ここに，H_α(φ)は近似されたHeavisideのステップ関数であり，次式で定義される。   

$$H_{\alpha }\left(\varphi \right)=\{\begin{array}{c}\frac{1}{2},\text{for}\alpha \le \varphi \left(\text{gas}\text{phase}\right)\\ \frac{1}{2\alpha }\left[\varphi +\frac{\alpha }{\pi }\text{sin}\left(\frac{\pi \varphi }{\alpha }\right)\right],\text{for}-\alpha <\varphi <\alpha \left(\text{transition}\text{phase}\right)\\ -\frac{1}{2},\text{for}\varphi \le -\alpha \left(\text{liquid}\text{phase}\right)\end{array}$$(9)

ここで，αをゼロに近づければ厳密なHeavisideのステップ関数に近づく。αは気液界面厚さの半分の厚さを表す。この近似ステップ関数の微分をとると，次式で示される近似されたDiracのデルタ関数になっているのがわかる。   

$$\frac{d{H}_{\alpha }\left(\varphi \right)}{d\varphi }={\delta }_{\alpha }\left(\varphi \right)=\{\begin{array}{c}0,\text{for}\alpha \le \varphi \\ \frac{1}{2\alpha }\left[1+\text{cos}\left(\frac{\pi \varphi }{\alpha }\right)\right],\text{for}-\alpha <\varphi <\alpha \\ 0,\text{for}\varphi \le -\alpha \end{array}$$(10)

近似ステップ関数が無次元であるので，この近似デルタ関数は，長さの逆数の次元を持つことがわかる。これで厚みが2αの界面をレベルセット関数で表現できたので，次に界面張力を考える。温度分布が無い状況においては，界面張力と曲率に応じて，界面に対して法線方向にのみ，次式で与えられる表面力(単位面積あたりの力)が作用する。   

$$f_{sa}=\gamma \kappa e_n\left[{\text{N/m}}^{\text{2}}\right]$$(11)

ここで，γは界面張力[N/m]，κは界面曲率[1/m]，e_nは界面法線単位ベクトルである。本来，気液界面厚さはゼロに限りなく近い。そのような無限に薄い領域に作用する界面張力が結果として，式(11)のような有限の値に落ち着くためには，界面法線力の絶対値は次式のように与えられるべきである。   

$$f_{sv}=\gamma \kappa \delta \left(\varphi \right)\left[{\text{N/m}}^{\text{3}}\right]$$(12)

ここに，δ(φ)はDiracのデルタ関数を示し，長さの逆数の次元を持つ。式(12)によれば，界面つまりレベルセット関数φ=0においてのみ無限大の力が作用するが，その界面を跨ぐ法線方向の積分値はγκと有限値となる。しかしながら，数値計算を安定に実行させるために界面に厚みを持たせた状況においては，CSFモデル⁷⁾を用いて，次式のように界面法線力ベクトルを近似的に与えることができる。   

$$f_{sv}=\gamma \kappa \frac{\text{grad}\rho }{{\rho }_L-{\rho }_G}=-\gamma \kappa {\delta }_{\alpha }\left(\varphi \right)\frac{\partial \varphi }{\partial x_i}\left[{\text{N/m}}^{\text{3}}\right]$$(13)

ここで，∂φ/∂x_iはレベルセット関数の勾配を表し，無次元である。レベルセット関数の再初期化の操作が適切に実行されている状況においては，その勾配の絶対値は1である。したがって，界面曲率κは次式のように与えられる。   

$$\kappa =\text{div}\left(\frac{\text{grad}\varphi }{\left|\text{grad}\varphi \right|}\right)=\frac{{\partial }^2\varphi }{\partial x_j\partial x_j}={\nabla }^2\varphi \left[1/\text{m}\right]$$(14)

2・3 無次元式

気泡の等価直径dおよび気相の物性値を基準にとり無次元化を行った⁸⁾。まとめると，解析に用いられる無次元の基礎方程式は以下の通りである。   

$$\frac{\partial U_j}{\partial X_j}=0$$(15)

  

$$\frac{\partial \Phi }{\partial \tau }+U_j\frac{\partial \Phi }{\partial X_j}=0$$(16)

  

$$\begin{array}{c}\frac{\partial U_i}{\partial \tau }+U_j\frac{\partial U_i}{\partial X_j}=-\frac{1}{{\rho }_{\Phi }}\frac{\partial P}{\partial X_i}+\frac{1}{{\rho }_{\Phi }}\frac{\partial }{\partial X_j}\left[{\mu }_{\Phi }\left(\frac{\partial U_i}{\partial X_j}+\frac{\partial U_j}{\partial X_i}\right)\right]\\ -\frac{\Gamma }{{\rho }_{\Phi }}\frac{{\partial }^2\Phi }{\partial X_j\partial X_j}{\delta }_k\left(\Phi \right)\frac{\partial \Phi }{\partial X_i}+F_i,F_i={}^t\left(0,0,-G\right)\end{array}$$(17)

  

$${\rho }_{\Phi }=\frac{\rho }{{\rho }_G}=\frac{1}{2}\left(\overline{\rho }+1\right)-\left(\overline{\rho }-1\right)H_k\left(\Phi \right)$$(18)

  

$${\mu }_{\Phi }=\frac{\mu }{{\mu }_G}=\frac{1}{2}\left(\overline{\mu }+1\right)-\left(\overline{\mu }-1\right)H_k\left(\Phi \right)$$(19)

  

$$H_k\left(\Phi \right)=\{\begin{array}{c}0.5,\text{for}k\le \Phi \\ \left[\Phi +\left(k/\pi \right)\text{sin}\left(\pi \Phi /k\right)\right]/\left(2k\right),\text{for}-k<\Phi <k\\ -0.5,\text{for}\Phi \le -k\end{array}$$(20)

  

$${\delta }_k\left(\Phi \right)=\{\begin{array}{c}0,\text{for}k\le \Phi \\ \left[1+\text{cos}\left(\pi \Phi /k\right)\right]/\left(2k\right),\text{for}-k<\Phi <k\\ 0,\text{for}\Phi \le -k\end{array}$$(21)

本解析に必要な無次元数(計算パラメータ)は以下の通りである。   

$$\Gamma =\frac{\gamma {\rho }_{G}d}{{\mu }_G{}^2},G=\frac{g{\rho }_G{}^2{d}^3}{{\mu }_G{}^2},\overline{\rho }=\frac{{\rho }_L}{{\rho }_G},\overline{\mu }=\frac{{\mu }_L}{{\mu }_G},k=\frac{\alpha }{d}$$(22)

式(22)において，Γはラプラス数，Gはガリレイ数である。その他，解析には密度比，粘度比，および無次元の界面厚さkが現れる。ここで，kは気泡等価直径dに対する界面半厚さを表し，数値計算上必要なパラメータである。本研究では，界面厚さは格子幅の3倍としているため，軸対称解析で75×300の格子点数を用いた場合はk=0.1となる。無次元の変数は，気泡等価直径d，気相の密度や粘度を用いて，以下のように定義した。   

$$\begin{array}{c}X=\frac{x}{d},U=\frac{u}{{\mu }_G/\left({\rho }_{G}d\right)},\tau =\frac{t}{{\rho }_G{d}^2/{\mu }_G},P=\frac{p}{{\mu }_G{}^2/\left({\rho }_G{d}^2\right)},\\ \Phi =\frac{\varphi }{d},H_k=H_{\alpha },{\delta }_k={\delta }_{\alpha }d\end{array}$$(23)

2・4 解析条件と方法

液相は10 cSt のシリコーンオイル(動粘度は1.0×10⁻⁵ m²/s，比重は0.935，表面張力は2.01×10⁻² N/m)とし，気相は常温常圧の空気(動粘度を1.5×10⁻⁵ m²/s，密度を1.2 kg/m³)として，各無次元数を見積もる。密度比および粘度比はそれぞれ   

$$\overline{\rho }=\frac{{\rho }_L}{{\rho }_G}=\frac{1000\times 0.935}{1.2}=779.2\cong 800$$(24)

  

$$\overline{\mu }=\frac{{\mu }_L}{{\mu }_G}=\frac{{\rho }_L{\nu }_L}{{\rho }_G{\nu }_G}=779.2\times \frac{1.0\times {10}^{-5}}{1.5\times {10}^{-5}}=519.4\cong 500$$(25)

となる。気泡体積が0.05 cm³の時，等価直径dは，   

$$d={\left(6/\pi \times 0.05\times {10}^{-6}\right)}^{1/3}\cong 0.00457\left[\text{m}\right]$$(26)

となるので，これより等価直径dは4.57 mmということになる。このとき，ラプラス数とガリレイ数は，それぞれ以下のような値をとる。   

$$\Gamma =\frac{\gamma {\rho }_{G}d}{{\mu }_G{}^2}=\frac{\gamma d}{{\rho }_G{\nu }_G{}^2}=\frac{\left(2.01\times {10}^{-2}\right)\left(4.57\times {10}^{-3}\right)}{\left(1.2\right){\left(1.5\times {10}^{-5}\right)}^2}=3.40\times {10}^5$$(27)

  

$$G=\frac{g{\rho }_G{}^2{d}^3}{{\mu }_G{}^2}=\frac{g{d}^3}{{\nu }_G{}^2}=\frac{\left(9.8\right){\left(4.57\times {10}^{-3}\right)}^3}{{\left(1.5\times {10}^{-5}\right)}^2}=4157\cong 4000$$(28)

気泡体積が10倍，40倍になれば，気泡径やラプラス数は，それぞれ2.154倍，3.420倍に，またガリレイ数は10倍，40倍になる。以下，本研究における解析条件をTable 1にまとめる。

Table 1.  Computational conditions.

Initial bubble volume, V	0.05 cm3	0.5 cm3	2 cm3	5 cm3
Initial equivalent diameter of bubble, d	4.57 mm	9.85 mm	15.6 mm	21.2 mm
Density ratio, ρ	800
Viscosity ratio, μ	500
Laplace number, Γ	3.40×105	7.32×105	1.16×106	1.58×106
Galilei number, G	4.00×103	4.00×104	1.60×105	4.00×105
Thickness of interface, k	0.100 (axisymmetric cases) 0.167 (3D cases)
Number of meshes (R×Z) (axisymmetric cases) (X×Y×Z) (3D cases)	75×300 (axisymmetric cases) 96×96×192 (3D cases)

解析領域は，気泡サイズに応じて変化させた。無次元で，気泡の等価直径を1とするとき，容器の直径(三次元解析では横幅)は5，高さは10とした。R-Z断面の解析領域について，それぞれの方向に等間隔のスタッガード格子を採用し，運動方程式の移流項には三次精度風上差分法を適用し，その他の項については二次精度中心差分法を適用した。圧力場の解法には，フラクショナルステップ法を用い，三次元解析ではマルチグリッド法を併用した。レベルセット関数については再初期化し距離関数の性質を適切に保つと同時に，気泡体積がその初期値からずれないようにNewton法を用いて体積補正を行った。

3. 解析結果

3・1 気泡の膨張を無視した場合

本節では，後述する0.013気圧の減圧容器内に対する解析結果の比較結果として，まず気泡の膨張を無視した場合の解析結果を示す。Fig.2は，初期気泡体積を0.05 cm³とした時の気泡の断面形状の時間推移を示す。左図は軸対称解析を，右図は三次元解析を示す。τ=0.06以降くらいでは，僅かに振動しているものの，気泡はほぼ一定の形状を保ちながら一定速度で上昇するという結果が得られた。左図では，τ=0.06以降の準定常状態における気泡の形状は上側が凹で下側が凸になっており，常識に反して多少理解し難いものになっている。左右の図において上昇速度が異なる結果となったが，これは気泡形状の違いに由来するものと推察される。後述のFig.8では実験で観察された気泡形状との比較を行っているが，計算格子がより密な軸対称解析の結果は実験に近いものであるが，三次元解析では，格子数が不十分であるためか，実験で得られた気泡形状を再現しきれていない。もう一点，この小さな気泡体積の場合は，気泡の代表長さと上昇速度に基づくレイノルズ数が小さくなるために，気泡まわりの境界層が厚くなる傾向がある。その結果，容器側壁と干渉することも考えられ，軸対称解析と三次元解析の違いが生じた理由と考えられる。

Fig. 2.

 A series of the bubble shape for the initial bubble volume 0.05 cm³. Time interval is 0.02.

左図の軸対称解析の場合に，τ=0.06以降の準定常な状態における気泡の上昇速度w_bubを見積もることにする。図から判断すると，無次元の上昇速度は，   

$$W_{bub}=\frac{dZ}{d\tau }\cong 62.5$$(29)

と見積もられる。式(23)より，実際の有次元の上昇速度は，以下のように換算される。   

$$w_{bub}={\nu }_G/d\cdot W_{bub}=\left(1.5\times {10}^{-5}\right)/\left(4.57\times {10}^{-3}\right)\cdot \left(62.5\right)=0.205\text{m/s}$$(30)

したがって，初期気泡体積が0.05 cm³の場合，10 cStのシリコーンオイル中を約20.5 cm/sで上昇することが予測された。 ただし，本計算では容器の直径を気泡等価直径のわずか5倍としているため，側壁の影響が多少残っている可能性は否定できない。また気泡の上昇速度に基づくレイノルズ数は以下のように定義される。   

$${\mathrm{Re}}_{bub}=\frac{{\rho }_L{w}_{bub}d}{{\mu }_L}=\frac{{\rho }_{L}d}{{\mu }_L}\frac{{\mu }_G}{{\rho }_{G}d}{W}_{bub}=\frac{\overline{\rho }}{\overline{\mu }}{W}_{bub}$$(31)

初期体積の異なる場合(0.5，2.0，5.0 cm³)を，それぞれFig.3，Fig.4，Fig.5に示す。また，気泡の上昇速度に関するまとめをTable 2に示す。三次元解析結果が示すように，気泡径が大きな場合においても，気泡上昇中の螺旋運動やジグザグ運動などの対称性の崩れるような流れは起きなかったことを確認している。その意味からすると，本稿で扱った軸対称解析は格子解像度や計算コストの面で利点を持っている。気泡径が大きなFig.4や5では，気泡上昇速度では二次元軸対称計算と三次元計算で良い一致を示しているが，気泡形状については，若干の違いが見受けられる。軸対称計算では，気泡は上昇するにつれ薄く扁平になる傾向が見られるのに対し，三次元計算では上昇するにつれ気泡形状の三次元性が見られる。この理由は次のように考えられる。気泡の初期体積が大きい場合は，レイノルズ数が大きくなり，気泡形状や速度が変化しないとされる終端状態に達するまでの時間は大きくなり，軸対称解析ではそれを示唆しているのに対し，三次元解析では解像度不足でそれを再現できていないと思われる。以下の結果では，軸対称解析の結果について言及する。

Fig. 3.

 A series of the bubble shape for the initial bubble volume 0.5 cm³. Time interval is 0.01.

Fig. 4.

 A series of the bubble shape for the initial bubble volume 2.0 cm³. Time interval is 0.005.

Fig. 5.

 A series of the bubble shape for the initial bubble volume 5.0 cm³. Time interval is 0.002.

Table 2.  Computational results for the rising velocity.

Initial bubble volume, V		0.05 cm3	0.5 cm3	2 cm3	5 cm3
Rising velocity of bubble, wbub	Axisymmetric cases	20.5 cm/s	22.2 cm/s	27.5 cm/s	32.9 cm/s
	3D cases	17.6 cm/s	22.2 cm/s	27.4 cm/s	32.8 cm/s
Reynolds number, Rebub	Axisymmetric cases	93.8	219	429	698
	3D cases	80.4	218	428	696

Fig.6は，閉曲線は気泡形状を，左半分には圧力分布を，右半分はStokesの流れ関数(軸対称系における流線)を示す。左半分の圧力分布を眺めると，気泡の近傍を除き，ほぼ静水圧分布になっているのがわかる。Fig.6(a)では容器の上端と下端の無次元圧力差は3.18×10⁸であった。Fig.6(b)では無次元圧力差は1.278×10⁹であった。つまり気泡体積に比例して無次元圧力差が大きくなることがわかる。この圧力差は他の任意の時刻においてもほぼ変わりなく成立する。次式に示されるように，(a)の場合について圧力差について有次元換算を試みる。   

$$\Delta p={\rho }_G{\nu }_G{}^2/d^2\cdot \Delta P=1.2\cdot {\left(1.5\times {10}^{-5}\right)}^2/{0.00985}^2\cdot \left(3.18\times {10}^8\right)=885\left[\text{Pa}\right]$$(32)

Fig. 6.

 An instantaneous contours of the pressure, the Stokes stream function and the bubble shape.

これは，高さが9.85 cmの容器内で885 Paの圧力差が生じていることに相当する。大気圧が101.3 kPaであることと比較すれば，0.885 kPaは大気圧の1%以下の小さな値であるが，この鉛直方向の圧力勾配によって気泡が浮力を受け上昇する。この際の気泡の膨張は無視できる。しかしながら，容器が鉛直方向に非常に長い場合，あるいは減圧された容器内においては，気泡は上昇過程において膨張すると考えられる。次節でそれを検証する。

3・2 気泡の膨張を考慮した場合(減圧容器)

減圧密閉容器内では，気泡の体積は静水圧の影響を敏感に受けると考えられ，気泡が上昇するに連れ，その体積が膨張する。例として，自由液面の表面上で0.013 atm(1317 Pa)の場合を考える。シリコーンオイルの密度(935 kg/m³)と地球上の重力加速度(9.8 m/s²)の影響により，高さ1 mにつき9163 Paの圧力勾配が発生する。今，初期気泡体積が0.05 cm³の場合を想定すると，気泡直径は約1 cm(正確には0.985 cm)であるから，本計算モデルでは，実サイズで容器直径が5 cm，高さが10 cmに相当する。この高さを静水圧に換算すると，916.3 Paになり，前節の結果にほぼ一致する。まとめると，容器底部で2233 Pa，最上面で1317 Paということになり，単純に考えれば，気泡は下から上まで上昇する間に約1.7倍に膨張することを意味する。これは，気泡界面に作用する表面張力を無視した場合であるが，ラプラスの式より直径1 cmの気泡の内外圧力差を算出すれば8 Pa程度であるので，この影響については無視され得る。本稿では，気泡の圧縮性流れを厳密に解くのではなく，単に気泡体積が静水圧の影響を受けて膨張する過程を単純に模擬することとした。

Fig.7は，初期気泡体積が2.0 cm³とした時の気泡形状の時間変化の比較を示す。τ=0.01以降くらいでは，気泡界面は球面形状の一部分(spherical-cup)を形成し，ほぼ一定速度で上昇しているように見えるが，気泡形状は完全には定常には至っていない。気泡の膨張を許す場合の方が，気泡が僅かながらより高い位置まで到達している。準定常状態に落ち着いた後の気泡形状に関しては，雰囲気圧力の影響はあまりないと考えられる。

Fig. 7.

 Comparison of the bubble shape without or with taking into account volume expansion for the case of initial bubble volume 2.0 cm³. Time interval is 0.005.

4. 実験と数値計算結果の比較

4・1 気泡形状について

実験⁴⁾と本数値計算では，初期条件が多少異なる。実験では，半球のダンピングカップに気泡をいったん保持し，ダンピングカップをゆっくり回転させることにより静かに気泡を放つ。計算モデルでは，ダンピングカップは模擬されておらず，初期の気泡形状は横長の回転楕円体と仮定している。これは，球形気泡を初期条件にした計算では，条件に依っては気泡が上昇途中で分裂する場合が観察されたこと，またダンピングカップに保持された気泡は扁平率が大きいことがその理由である。以上から，実験と計算における気泡形状の比較は，ある程度時間経過したのちの準定常的な状態で行った。

0.013気圧の減圧下における，それぞれ4ケースにおける気泡形状の比較をFig.8に示す。実験では数値計算よりも気泡上昇中の膨張率が大きいように思われ，また気泡形状も必ずしも準定常的ではなかったのに対し，数値計算では軸対称解析，三次元解析のいずれにおいても，気泡形状はほぼ準定常的であった。単純な比較はできないが，総じて4ケースとも，実験と計算は良く一致している。(a)の小さな気泡では，気泡は逆に反り返っているのが実験でも数値計算でも見て取れる。ただ，気泡の扁平率は実験と数値計算では多少異なるように思われる。数値計算では界面に有限の厚みを持たせたものとして扱われるので，扁平な気泡の挙動を扱うには十分な解像度を持つ格子点の配備が必要とも考えられる。(b)や(c)の中程度の気泡では，実験写真では，気泡内部のへこみは確認しずらいが，軸対称の数値計算結果ではそれを明確に確認できる。その見え方の違いを考慮すると，比較的よく一致している。(d)の大きな気泡では，実験では気泡が上下に分裂しているのに対し，数値計算ではそれは再現されていない。この違いは初期条件による違いおよびレベルセット法を用いているなどの数値モデルにおけるいくつかの仮定が原因であるとも考えられる。

Fig. 8.

 Comparison of the bubble shape for the several initial bubble volumes. The left hand side indicates experimental photographs and the right hand side indicates the numerical 2D views.

4・2 気泡上速度について

既にTable 2に，気泡上昇速度およびレイノルズ数(気泡上昇速度，球に等価な直径，および液相の動粘度に基づく)については示した。ここでは，その表に記載した数値が妥当なものであるかどうかについて，実験結果⁴⁾およびGraceの気泡分類図⁹⁾との比較から検討する。まず，Graceの気泡分類図を用いて，数値計算結果の妥当性を検証する。そのために必要な無次元数は，エトベス数，モルトン数，レイノルズ数であり，それぞれ以下のように定義される。   

$$\begin{array}{c}E\ddot{o}=\frac{g{d}^2\left({\rho }_L-{\rho }_G\right)}{\gamma }=\frac{G}{\Gamma }\left(\overline{\rho }-1\right),Mo=\frac{g{\mu }_L^4\left({\rho }_L-{\rho }_G\right)}{{\rho }_L^2{\gamma }^3}=\frac{G}{{\Gamma }^3}\frac{1}{\overline{\rho }}\left(1-\frac{1}{\stackrel{⌢}{\rho }}\right){\overline{\mu }}^4,\\ {\mathrm{Re}}_{bub}=\frac{{\rho }_L{w}_{bub}d}{{\mu }_L}=W_{bub}\frac{\overline{\rho }}{\overline{\mu }}\end{array}$$(33)

これらの無次元数のうち，エトベス数とモルトン数は，Table 1に示す通り，解析条件が与えられると自動的に決まる。特にモルトン数は気泡体積に無関係に決まり，Table 3に示されるように，およそ10⁻⁵の値をとる。レイノルズ数は，気泡の終端上昇速度が含まれるため，数値計算を終えた後に算出できると考えられる。一方で，Fig.9のGraceの分類図においては，エトベス数とモルトン数が決まれば，レイノルズ数が図より読取れることを示唆する。その図から予測されるレイノルズ数(目視での読取り値)と本数値計算結果から得られたレイノルズ数を比較することで数値計算結果の妥当性を検証できると考えられる。Fig.9の丸印およびTable 3の数値に示されるように，本数値計算で得られたレイノルズ数は，Graceの図から予測されるレイノルズ数に極めて良く一致していることが確認できる。また同時に，Graceの図から気泡形状も与えられるが，気泡体積が大きくなるに従い，エトベス数が増加し，EllipsoidalからSpherical-cupに移行すると予測できるが，その点に関してもおよそ一致していると思われる。

Table 3.  Comparison between Grace’s and present results.

Initial bubble volume, V	0.05 cm3	0.5 cm3	2 cm3	5 cm3
Initial equivalent diameter of bubble, d	4.57 mm	9.85 mm	15.6 mm	21.2 mm
Eötvös number, Eö	9.56	44.4	111	206
Morton number, Mo	9.86×10–6
Computed Reynolds number, Rebub	93.8	219	429	698
Predicted Reynolds number, ReT	95	220	410	700

Fig. 9.

 Grace’s classification of bubble behavior together with the present four results.

最後に，Table 2に示された上昇速度は常圧下における数値結果であるが，減圧下を模擬した計算結果とほとんど差が無かったことから，これらの値が実験値(0.013気圧)と比較して，妥当かどうかを検証する。実験では，気泡体積が0.5 cm³，2.0 cm³の場合には，それぞれ22 cm/s，27 cm/s程度で上昇するという結果が得られており，Table 2と非常に良い一致が見てとれる。

5. まとめ

密閉容器内のシリコーンオイル中を上昇する単一気泡の上昇運動について，軸対称ならびに三次元数値流体解析を行った。以下をまとめとする。

・初期気泡体積に応じて，気泡の形状や上昇速度が大きく異なることが本数値計算によって示された。本研究で対象とした初期気泡体積の範囲内においては，それが大きくなるにしたがって，気泡の上昇速度が大きくなることが確認された。

・膨張を考慮しない場合において，無次元数を用いて文献値⁹⁾と定量的に比較を行ったところ，形状や上昇速度についてはほぼ一致することを確認した。

・気泡膨張を考慮する場合には，実験でも数値解析でも気泡の上昇速度が時間とともに若干加速する傾向が見られたが，気泡膨張の有無による気泡形状の違いは殆ど見られなかった。

文献

• 1)   C.W.  Hirt and  B.D.  Nichols: J. Comput. Phys., 39 (1981), 201.
• 2)   S.  Osher and  J.A.  Sethian: J. Comput. Phys., 79 (1988), 12.
• 3)   M.  Sussman and  E.G.  Puckett: J. Comput. Phys., 162 (2000), 301.
• 4)   T.  Kumagai,  Y.  Kishimoto,  O.  Nakamura,  T.  Tagawa,  T.  Himeno and  S.  Shimasaki: Tetsu-to-Hagané, 101 (2015), 93.
• 5)   O.  Nakamura,  T.  Kumagai and  K.  Takatani: Tetsu-to-Hagané, 101 (2015), 117.
• 6)   T.  Himeno and  T.  Kumagai: Tetsu-to-Hagané, 101 (2015), 109.
• 7)   J.U.  Brackbill,  D.B.  Kothe and  C.  Zemach: J. Comput. Phys., 100 (1992), 335.
• 8)   Y.  Shibasaki,  K.  Ueno and  T.  Tagawa: ISIJ Int., 50(2010), 3.
• 9)   R.  Clift,  J.R.  Grace and  M.E.  Weber: Bubbles, Drops and Particles, Dover publications Inc., New York, (2005), 27, (Originally published by Academic Press Inc., New York, (1978)).