気泡および液滴運動の予測

0.5 1 5 10^-1

10⁰

d [mm]

VT [m/s] ^measured(Tomiyama et al, 2002) zigzag

helical rectilinear

CD model (Tomiyama et al, 1998) pure

slightly contaminated fullycontaminated

predicted

図3.4.1 無限静止水中空気泡の終端上昇速度

x₁ x₂ x₃

bubble

top boundary:

uniform inflow bottom & side boundaries:

continuous domain:

8d x 8d x 8d spatial resolution:

d* = 16 g

図3.4.2 無限静止水中単一空気泡計算用計算体系

するものであるが，観測される上昇速度にはばらつきがある．なお，図中 zigzag(揺動), helical(螺旋), rectilinear(直線)は，気泡の運動状態を示している．

無限静止水中空気泡の数値計算を行い，終端速度の計算結果を上述の実験結果および抗 力モデルと比較した．計算領域を，図3.4.2に示す8d×8d×8dの立方体とし，上部境界 から一様な速度で水を流入させることにより長時間気泡が計算領域内にとどまるようにし た．上部以外の境界には連続条件を課して，無限静止水中気泡の条件を模擬した．気泡直 径あたりに割り当てるセル数を16とした．計算に用いた二相の物性は表3.3.1と同じと した．

計算の結果，気泡体積の誤差は1%未満で，かつ界面は１セル幅に収まった．気泡終端

ー . ; . • • ^{，  '‑}

速度の計算結果を図3.4.1に三角記号で示す．すべての計算結果は，高純度系に対する実 験結果のばらつきの範囲内に収まっており，妥当な予測であることが確認できる．気泡運 動を良好に予測するためには表面張力を精度よく評価できることが重要であるが，終端速 度を妥当に予測できていることから，本手法は精度よい表面張力評価性能を備えていると 言える．

3.4.1.2 無限静止液中単一空気泡

Hnatら[7]は，ミネラルオイル中を上昇する冠球形単一空気泡に関する実験を行ってい る．以下で，Hnatらと同条件で気泡計算を行い，終端速度，気泡形状および伴流形状の予 測精度を調べた．

本研究では，図3.4.3に示す計算体系を用いて，Hnatらの無限静止液中単一空気泡実験 を模擬した計算を行った．物性値は，ρ_L =875.5 kg/m³,µ_L =118 mPa s,ρ_G = 1.23 kg/m³, µG =0.0177 mPa s,σ= 32.2 mN/mとした．この条件では，気泡形状は冠球形となる．本 計算では２次元円柱座標系用NSSを用いた．空間分解能を変えた数条件の計算を行った．

気泡球等価直径 dあたりのセル数d^∗ = d/∆xが48以上の条件では∆t = 1.0×10⁻⁶ sec， 48未満の条件では∆t= 1.0×10⁻⁵ secとした．

計算の結果，体積誤差はいずれの計算条件でも0.02% 未満であった．また，界面は シャープに維持され，１セル幅に収まった．図3.4.4に気泡上昇速度の計算結果を示す．

Hnatらの実験では，気泡終端速度は0.215 m/sである．いずれの分解能でも予測結果は実 験結果とよく一致した．図3.4.5に気泡形状および気泡周りの流線の計算結果を示す．た だし，図中の左半分はHnatらによる気泡および伴流の実験写真である．気泡形状，伴流 形状ともに，実験結果とよく一致しており，NSSが気泡形状・終端速度に加えて，気泡周 囲の速度場も正しく予測できることを確認できる．

C L

air bubble d =12.15 mm ρ =1.23 kg/m

³

µ = 0.0177 mPa

^.

s

12d = 144 mm 4d = 48 mm

mineral oil ρ = 875.5 kg/m

³

µ = 118 mPa

^.

s σ = 32.2 mN/m symmetry

continuous

no-slip

図3.4.3 無限静止液中単一空気泡計算用計算領域

0 0.1 0.2 0.3

d/ ∆ x=12 V

_T

=0.217

d/ ∆ x=16 V

_T

=0.213

0 0.1 0.2 0.3

d/ ∆ x=24 V

_T

=0.212

d/ ∆ x=32 V

_T

=0.211

0 0.1 0.2 0.3

d/ ∆ x=40 V

_T

=0.211

d/ ∆ x=48 V

_T

=0.211

0 0.05 0.1

0 0.1 0.2 0.3

d/ ∆ x=56 V

_T

=0.211

0 0.05 0.1

d/ ∆ x=64 V

_T

=0.211

t [sec]

V

T

[m /s ]

図3.4.4 無限静止液中単一空気泡の上昇速度

d* = d/ ∆ x = 16 d* = 32

d* = 56 d* = 64

図3.4.5 気泡形状および気泡周りの流線(いずれも左半分はHnat [7]らの実験写真)

3.4.1.3 無限静止液中単一液滴

無限静止液中を上昇する液滴の終端上昇速度VT は，エトベス数Eo，モルトン数M，粘 度比 µ^∗ = µD/µL によって決まると考えられる．すなわち，気泡の場合とは異なり液滴の 粘性が重要になる．よって，液中液滴の運動を精度よく計算するためには，液滴の内部流 れおよび液液界面での粘性応力を正しく評価しなければならない．以下で，Myint [8]ら のグリセリン水溶液中単一シリコン液滴実験と同条件で静止液中液滴の界面追跡計算を行 い，数値予測した液滴運動の妥当性を検討した．

本計算では２次元円柱座標用NSSを用いた．計算した二相の無次元数を表3.4.1に示 す．液滴球等価直径 dあたりのセル数d^∗ を12とした．境界条件の影響を無視できるよ う，計算領域幅を8dとした．

計算の結果，体積誤差はいずれの計算条件でも10⁻⁴% 未満であった．また，界面の シャープさも良好に維持され，界面は１セル幅に収まった．よって，粘性応力をシャープ な界面に対して評価できる．液滴終端速度の計算結果および実験結果を，液滴レイノルズ 数に換算して図 3.4.6に示す．図中，黒色記号が計算結果，白色記号が実験結果である．

すべての条件で，計算結果は実験結果とよく一致している．図3.4.7に示すように，液滴 形状も良好に予測できている．以上，本手法は表面張力・液滴内部流れ・液液界面におけ る粘性応力を正しく評価できる．

表3.4.1 液‐液系の無次元数

Case Eo log M µ^∗ ρ^∗

A 0.30−17.3 −1.0 0.11 0.77

B 0.24−9.5 −3.0 0.37 0.79

C 0.32−6.7 −4.8 0.94 0.80

10^-1 10⁰ 10¹ 10^-3

10^-2 10^-1 10⁰ 10¹ 10²

Eo

Re

Case A (logM=-1.0, µ*=0.11) Case B (logM=-3.0, µ*=0.37) Case C (logM=-4.8, µ*=0.94) closed symbols : predicted

open symbols : measured

図3.4.6 無限静止液中単一液滴の液滴レイノルズ数(実験結果はMyintらによる[8])

(a) Case A (Eo= 1.2, log M=−1.0, µ^∗ =0.11)

(b) Case B (Eo= 6.1, log M=−3.0, µ^∗ =0.37)

(c) Case B (Eo =9.5, log M =−3.0, µ^∗ =0.37)

(d) Case C (Eo= 3.5, log M=−4.8, µ^∗ =0.94)

図3.4.7 液滴形状の計算結果と実験結果[8]

会 A ω o c J j

‑

^O O 

θ 

OA 口

x₁ x₂ x₃

bubble

top boundary: continuous bottom boundary:

inflow (uniform shear flow) side boundaries (x₁-direction):

moving boundary

side boundaries (x₂-direction):

slip wall

bubble diameter:

d*=14 g

図3.4.8 一様せん断流中空気泡計算用計算体系

3.4.2 一様せん断流中気泡に働く揚力

一様せん断流中を運動する空気泡の数値計算を行い，気泡揚力係数を評価した．計

算条件を Tomiyama [9]のグリセリン水溶液中空気泡実験と同様に設定し，計算結果を

実験と比較した．計算に用いた物性を以下に示す．ρ^∗ = ρG/ρL = 0.001, µ^∗ = µG/µL = 0.00064, log M = −4.7．また，0.038 < Eo < 6.14 の範囲を計算対象とした．いずれの 条件でも，気泡直径 d に対して割り当てるセル数d^∗ を14 とした．Eo = 0.38, 0.68の 条件では各々∆t = 1.0×10⁻⁶, 2.0×10⁻⁶ secとし，これらよりも大きい Eoではすべて

∆t =5.0×10⁻⁶ secとした．図3.4.8に示すように，x1 = 0の境界には −x3 方向に一定速 度で移動する移動壁条件，それに面対する境界には+x₃方向に一定速度で移動する移動壁 条件を課した．x2 = const.の２つの境界にはともにすべり壁条件を課した．下部境界は一 様せん断速度分布で固定し，上部境界は連続条件とした．本条件では，x1 =0の境界側に 気泡が移動する場合揚力係数は正，その逆の方向に移動する場合揚力係数は負である．

3.4.2.1 揚力係数計算法

気泡の運動方程式から導出した揚力係数計算式に，界面追跡計算により求めた気泡重 心の時間変化を代入して揚力係数を得る．以下に，揚力係数計算式[4]を導出する．ただ

し，図3.4.8に示す一様せん断流中気泡を対象とする．

↓ 1 0  

気泡の運動方程式は次式で与えられる．

(ρB+CV MρL) d uB

d t =−3CDρL|uB−uL|

4d (uB−uL)−CLρL(uB−uL)×∇×uL−(ρL−ρB) g (3.4.7) C_{V M} は仮想質量係数であり，0.5 とする．下付添字 Bは気泡を意味する．ここで，液相 速度uL および重力 gは x3 方向に平行で，気泡速度uBの x2 方向成分はゼロである(図 3.4.8)．

u_B =(u_B)₁e₁+(u_B)₃e₃ (3.4.8)

uL = (uL)₃e3 (3.4.9)

g=−ge₃ (3.4.10)

ここで，ei はxi 方向の基底ベクトルである．図3.4.8の一様せん断流れでは，液相速度勾 配テンソルの成分のうち，非零の成分は∂(uL)3/∂x1 のみである．以上より，気泡運動方程 式は以下のように書ける．

(ρ_B+C_{V M}ρ_L) d (uB)₁

d t =−3CDρL|uB−uL|

4d (u_B)₁−C_Lρ_L

(u_B)₃−(u_L)₃∂(uL)₃

∂x₁ (3.4.11)

(ρB+CV MρL)d (uB)₃

d t =−3CDρL|uB−uL| 4d

(uB)₃−(uL)₃

+CLρL(uB)₁ ∂(uL)₃

∂x₁ +(ρL−ρB) g (3.4.12) これら２式からCD を消去して次式を得る．

CL=

ρ_L|uB−u_L|∂(u_L)₃

∂x1

!₋₁"

(u_B)₁M_Bd (u_B)₃ d t −

(u_B)₃−(u_L)₃

M_Bd (u_B)₁

d t −(u_B)₁(ρ_L−ρ_B) g

#

(3.4.13) 上式が以下の計算で用いた揚力係数評価式である．上式中 M_Bは次式で定義した．

MB =ρB+CV MρL (3.4.14)

3.4.2.2 計算結果

図3.4.9(a), (b)にEo = 1.54, 6.14の気泡運動の予測結果を示す．直径の小さい Eo = 1.54の気泡は球形に近い形を保ったまま，x1 = 0の境界側へ移動している．したがって 揚力係数C_L は正である．一方，Eo= 6.14の気泡は扁平化しながら x₁ の増加する方向へ

0 sec 0.3 sec 0.6 sec

0 sec

0.45 sec 0.6 sec

x₁ x₃

x₁ x₃

(a) (b)

図3.4.9 一様せん断流中空気泡の横方向運動

移動している．したがってC_L は負である．本計算結果は，小気泡は変形量が少なく正の 揚力係数となり壁方向へ移動し，一方大気泡は変形量が大きく負の揚力係数となり管路 中央へ移動するという，管内流水中気泡に働く揚力に関する知見[9]と一致している．図

3.4.10にTomiyamaの揚力実験と計算による揚力係数の比較を示す．すべてのEoの範囲

で，揚力係数の計算結果は実験結果に近い値となっている．実験に見られる揚力係数の符 号のEo依存性も，計算により正しく予測できている．

ドキュメント内 2007 2 (ページ 110-120)