第 25 回数値流体力学シンポジウム
E11-1

多相流場の解法による上昇気泡群の気液界面面積評価
Evaluation of gas-liquid interface area of rising bubbles by multiphase incompressible-flow solver
藤岡奨, 京都大学大学院工学研究科, 京都府京都市西京区桂, susumu.fujioka@kx2.ecs.kyoto-u.ac.jp

牛島省, 京都大学学術情報メディアセンター, 京都府京都市左京区吉田本町, ushijima@media.kyoto-u.ac.jp

Susumu FUJIOKA, Kyoto University, Katsura, Nishikyo-ku, Kyoto 615-8530, Japan
Satoru USHIJIMA, Kyoto University, Yoshida-honmachi, Sakyo-ku, Kyoto 606-8501, Japan
In order to investigate gas transfer phenomena through the interface of multi phase fluids, a close look at theinterface is necessary. This paper discusses a method to obtain estimated interface curve in volume of fluid method,which is originally not good at accurate definition of exact interface curve. This paper proposes estimation methodsof exact interface which is nearly equivalent to evaluation methods of interface area, one of which utilizes equivalent
volume of half a quarter of sphere, the other deter minds equivalent cut plane of the computational cell. In theformer method, the interface area is evaluated from the volume ratio of the fluid, assuming that the shape of thefluid volume is half a quarter of sphere. In the latter method, a cutting plane is obtained according to the volume
ratio and normal vector of the interface. The accuracy of those methods are investigated in the case of evaluating
primitive shapes such as sphere, spheroid, and torus with diﬀerent resolutions. The investigation showed higher
accuracy in equivalent cut plane method. Those methods are applied to evaluation of interface area of risingbubbles.

緒言
近年の地球温暖化問題への注目の高まりから大気海洋
間の CO2 交換現象は，精力的な研究がなされている分野
である 1 ．気液界面を通した物質交換現象は，大気海洋
間の CO2 交換現象 2 にとどまらず，物質の撹拌や反応促
進など機械工学や化学工学においても盛んに研究されている．気液界面の物質交換現象は，様々な分野において現象論的に定式化の試みがなされている．
ガス交換の総量を求める場合には，気液界面面積を求めることが必要となる．開水路や海洋の場合は，水表面積の評価が問題になることは少ないと考えられるが，湖沼に設置されたエアレーション施設によるガス交換促進効果等を定量的に把握したい場合などは，水中の気泡の表面積の和によって界面を通してのガス交換を評価する
方法が考えられる．
本研究では，多数の気泡を含む流れ場において，ガス交換総量等を評価する目的で，気液界面積の定量評価法
について検討した．
複数の気泡を含む流れ場において気液界面積を定量的に評価するためには，計算に必要となる計算格子数が，気泡直径に対して少なくなることが望ましい．
一方で，気泡の形状を精度よく再現し，それらの表面積を推定することが必要である．
本研究では，気液混相流れ場の数値解法として多相流場の解法である MICS(Multiphase Incompressible flow
solver with Collocated grid System ) 3, 4 を用いた．気液界面形状および気液界面積は，各計算セルにおける流体の体積率から推定することとし，推定法として 2 通りの方法を試した．一番目の方法は，計算セルにおける流体の体積率から等価球体の表面積を用いて推定する方法である．二番目の方法は，計算セルにおける法線ベクト
ルと流体の体積率から，計算セルを切断する平面を求め，切断面の面積を求める方法である．
これらの手法によって表面積を求めた場合にどの程度の誤差を含むかは自明ではないため，種々のベンチマー
クを通して手法の有効性を比較した．また，多相流場の解法による上昇気泡の数値計算を行い，提案手法を適用
して界面の表面積追跡を行った.

1.

2.

表面積評価法

2.1 等価球体片を用いた表面積評価
計算セル内の形状を，一様に球体の８分の１の形状として仮定しボクセルデータから表面積を推定する．
計算セル内の流体体積率を φ，計算セル体積を V とし
1た場合，等価球体の半径 r を次式にて求めた．

8V φ =

4 3
πr
3

(1)

計算セル内に含まれる球面の面積 S は，全体の８分の１となるため，次式で得られる．

S=( 3V φ )2/3
4πr2
= π 1/384(2)

上式から明らかな様に，流体の体積率のみから表面積を推定可能であるため，非常に簡便である．
ただし，既往文献 5 に記された方法で，流体識別子を定義し，以下のパターンの計算セルのみ，次式で求めた．

S = ∆x∆y

(3)

S

re
∆z

v2
O

∆x

Se

re
∆y

O

v2
re

Fig. 1: Estimation of surface area using equivalent halfa quarter of sphere2.2 等価切断面を用いた表面積評価
体積 Va の計算セルにおいて，法線ベクトル n と流体の体積率 φ が与えられた場合，法線ベクトルが n で，切断後の体積片の体積が φVa であるような切断平面を求め，計算セルの切断面の面積を表面積として求める方法である．
例えば，Fig.3 に示される通り，流体の体積率 φ によって計算セル体積が V1 および V2 に分けられることが分かっており，かつ分割する面を代表する法線ベクトル n が既知である場合に，計算セルをそのように切断する平面のCopyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-1

t が与えられた場合，切断面の式は，次式で求められる．
nx (x − cx ) + ny (y − cy ) + nz (z − cz ) = 0
(9)

C

B

D

A
n

Δz

H

y

x = −(ny (y0jk − cy ) + nz (z0jk − cz ))/nx + cx

G

F
E

z

一方で，格子の辺を構成する線分との交点は，次式で求
められる．

z0

P

0 < x < ∆x であれば，格子上に交点を持つものとする．ただし j = 0, 1, k = 0, 1 であり，αijk は，点(i∆x, j∆y, k∆z) の座標の α 成分である．
これを全ての辺に対して繰り返し，切断面上の点の集
合 S を得る．
また，任意セル頂点 a に対して識別子 Φ(a) を次式で定義した．これにより，各頂点と切断面との位置関係が求められる．
Φ(a) = sgn(nx (x − cx ) + ny (y − cy ) + nz (z − cz )) (11)

Δy

Δx
x

Fig. 2: A case where assumption of an equivalent halfa quarter of sphere is not applied

位置 tn を求めるものである．ただし，t は切断面の位置を決めるパラメーターである．本研究ではこの切断平面
を等価切断面と称する．
以下に等価切断面の計算方法の詳細を示す．
任意平面によって切り取られる計算セルの体積を求める関数を Ψ とする．

V = Ψ(n, c)(4)

ここに，n は平面の法線ベクトル，c は平面上にある点の位置座標であり，平面の空間的な位置を決定する．V0 は与えられた流体体積率から得られる流体の体積として，関
数 f = Ψ − V0 を定義し，ニュートン法によって f = 0 解
き，計算セルを体積 V0 及び Va − V0 に分ける切断平面位置を算出した．ニュートン法の適用にあたっては，Fig.3にならい，平面の位置を決める座標 c = O + tn と表現した．計算セル内の原点座標 O は，法線ベクトル n の各軸方向の成分の符号から求め，O を始点とする法線ベク
トルが計算セル体積内に向かう様に決定した．すなわち，次式によって計算セル体積の原点とする頂点を法線ベク
トルの各成分の符号から求めた．

Oα = sgn(nα )∆α(α = x, y, z)(5)

ここに，∆α は α 方向の格子幅とする．
従って，法線ベクトルを固定した場合，関数 Ψ は t の関数となり，次式で表される．

V = Ψ(t)(6)

t の初期値として t0 = V0 /Va を与え，以下の漸化式によって t が収束するまで繰り返し計算を行った．
tn+1 = tn −fn
ΨntΨ(tn ) − Ψ(tn−1 )tn − tn−1

Φ(a) は切断面によって切り取られる体積片の各頂点が切断面のどちら側に属するか仕分けるのに便利である．原
点 O における識別子の値 Φ(O) から，同じ体積片に含まれる頂点の集合 T を選ぶことが出来る．切断面によって切り取られた体積片を構成する頂点の集合 U は，これらの頂点 T と切断面上の点の集合 S0 の和となる．これらの頂点 U の座標から，体積片の重心 G は容易に求められる．この体積片の体積を求めるために，Fig.4 に示される様に重心 G を頂点とする複数の四面体に分割する．このとき，計算セルを構成するそれぞれの境界面上に存在するセル頂点と切断面の頂点を合わせた集合 Si , i = 1, · · · 6を考える．i は６面体のそれぞれの面を示す指標である．
この集合の要素数は，そのセル境界面に切断面によって切り取られて生じた多角形の頂点の数を示し，高々６である．重心 G を頂点とし，セル境界に生じた多角形を底面とする多角形錐の体積を求めるために，集合 Si に含まれる座標から，面積が重複しない三角形の頂点の組み合わせを選び出す必要がある．Fig.5 に例を示す．この例では，黒丸で示されたセル頂点３点及び黒四角で示された
切断面上の頂点２点が同一セル境界面上に存在している．
これら５点の座標 (ξi , i = 1 · · · 5) は既知であるが，これらの座標がどの順番で保存されているかは分からない状態であるとする．最初は，Fig.5 の左の図に示される様に任意の３点 (ξ1 , ξ2 , ξ3 ) を選び出し，一つ目の三角形を作る．次に４点目 ξ4 を選び出し，先に選んだ３点の内の２点と合わせて三角形を作ることを考える．
点 i と点 j を結ぶ線分と点 i と点 k を結ぶ線分のなす角を θi;jk と表現すると，4 点目と三角形を構成するべき２点は，次の様に選ばれる．
(j, k) | cos θ4;jk = min cos θ4;jk(12)(j,k)∈C3(7)n
ただし，tn ,f n ,Ψnt は，n ステップ目の t,f = Ψ(t ) − V0 ，
Ψ の t による微分値である．Ψ(t) の関数形を数式で表すことは複雑であるため，本研究では，次式にて微分値 Ψtを近似した．

Ψt ≈

(10)

(8)

上式を当てはめる場合，初期値として t0 , t1 を与える必要がある．このため，本研究では，t1 = 0.9t0 とした．
Ψ(t) は以下のアルゴリズムによって求めた．

ただし，Cn は n 個の整数から２個の整数を取り出す組み合わせを元とする集合で，C3 = {(1, 2), (2, 3), (3, 1)} である．つまり，cos θ4;jk を最小とする (j, k) の組み合わせ
が，4 点目と三角形を構成するべき２点である．以下，５点目と三角形を構成すべき２点は，既出の４点から 2 点を取り出す全ての組み合わせの中から cos θ5;jk を最小とする (j, k) の組み合わせを選ぶ．6 点目が存在した場合には，上記の手順を繰り返して既出の５点の中から 2 点を選び出す．
切断面を底面とする多角錐についても，同様の考え方
で G を頂点とする四面体への分割を行うことが出来る．
従って，切り取られた体積片の体積 Φ(t) は次式で求められる．

Φ(t) =(1) 計算セルを構成する 12 辺と平面の交点個数と座標とを求める．

∑V (Si )

(13)

i=0,6

ただし，V (Si ) は，重心点 G を頂点とし，頂点の集合 Siからなる多角形を底辺とする多角錐の体積である．

(2) 交点座標から切断面面積を求める．
2

Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-1

1

n
S

v1
3

v2
O

4

1

2

1

4

5

5

5

tn

Δz

4

3

3

2

2

Fig. 5: Division of polygon into triangles

ΔyΔx
Tab. 1: EL1 of volume with diﬀerent number of subcells(Number of cell is 323 )Fig. 3: A computational cell and cutting plane

Number of subcells

EL1

3

4
103

0.086 (%)
-0.011 (%)

163
243

-0.0078(%)
0.0038 (%)

S
G

Δz

O

ただし，Sexact は S = 4πr2 にて得られる理論値，Ssimu.
は Case A あるいは Case B の手法によって得られる表面積である．Tab.3 に示される様に，Case B においては格子分割数の増加とともに精度の向上がみられるが，Case
A においては，643 程度の解像度までは一定の精度向上
傾向を示すが，その後は逆に誤差が大きくなる傾向を示した．これは，格子分割数が増加するとともに計算格子に切り取られる球面は平面に近づくが，Case A では，常に界面形状を８分の１球体の球面と仮定するために，過大評価になるものと考えられる．

ΔyΔxFig. 4: Division of fragment volume into tetrahedrons

3.

表面積評価精度検証計算

3.1
球体表面積評価
提案手法による３次元任意形状の表面積評価の精度について検証するため，球体の表面積評価誤差について以下で議論する．計算領域は 1 × 1 × 1 とし，計算領域の中心に半径 r の球体が配置されている．本研究では，
r = 0.25 に固定した．計算領域は，計算ケースによって
43 , 83 , 163 , 323 , 643 , 1283 , 2563 に等分割され，サブセル法
によって各計算セルに球体の占有体積率を定義した．サ
ブセル法の使用に当たっては，各計算セルを更に 103 個のサブセルに分割することとした．サブセル法によってボクセルデータに変換された後の体積誤差 EL1 は次式によって定義した．
S− Ssimu
EL1 = 100 × exact
SexactS ≈ 4π( ap bp + bp cp + cp ap )1/p
3

(16)

ただし，a, b, c は回転楕円体の各軸方向の径であり，p =1.6075 とした．

(14)

なお，サブセル法の分割数を決めるために，323 のケースにおいてサブセル分割数を 43 , 103 , 163 , 243 として，体積誤差の減少傾向を確認した．Tab.1 に示される様に，サブセルの分割数が増えるに従い，体積誤差 EL1 の減少傾向が認められる．
球体形状に対する表面積評価の検証問題における結果
を Tab.2 及び Tab.3 に示す．体積誤差については，全てのケースにおいて 1%以下となっており，ボクセルデータへの変換が精度よく行われたことが確認できる．
表面積評価法としては，等価球体を用いる方法 (CaseA) および等価切断面を用いる方法 (Case B) を用い，次式で得られる理論解との相対誤差 EL1 を求めた．

S− Ssimu
EL1 = 100 × exact
Sexact3.2
回転楕円体表面積評価
次 に ，回 転 楕 円 体 の 表 面 積 評 価 に つ い て 検 討 した ．回 転 楕 円 体 の 諸 元 は ，Fig.6 に 示 さ れ る 通 り とし た ．計 算 領 域 は 1 × 1 × 1 と し ，格 子 分 割 数 は ，
43 , 83 , 163 , 323 , 643 , 1283 , 2563 とした．ボクセル化に伴う体積誤差を Tab.4 に示す．
回転楕円体の表面積の理論解は，近似式として次式が知られているため，これを利用した．

Tab. 2: EL1 of sphere volume and number of cells

(15)

3

Number of cells

EL1

43

0.84 (%)

83
163

0.12 (%)
0.0050 (%)

323
643
1283

-0.011 (%)
0.0081(%)
0.0018 (%)

2563

0.00019 (%)

Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-1

Tab. 3: Comparison of EL1 of surface area of sphare

Tab. 5: Comparison of EL1 of surface area of spheroid

Number of cells

EL1 (Case B)(%)

NUmber of cells

EL1 (Case A)(%)

EL1 (Case B)(%)

17.4

3

4

26.0

-19.1

3

EL1 (Case A)(%)

3

4

28.1

3

8
163
323

18.8
-2.02
-2.35

2.73
1.41
0.13

8
163
323

24.8
-4.50
8.58

-0.60
0.42
-0.25

643
1283

-0.70
-6.07

0.79
0.65

643
1283

1.61
-2.8

0.20
0.35

2563

-9.27

0.56

2563

-2.70

0.25

z

0.15 y
0.1

0.35
O

0.15
0.4

x

Fig. 6: Size of spheroid

Fig. 7: Size of torus

Tab.5 に，Case A と Case B による表面積評価の比較を示す．Case B においては格子分割数の増加とともに精度の向上がみられるが，Case A においては，643 程度の解像度までは一定の精度向上傾向を示すが，その後は逆に誤差が大きくなる傾向を示した．
3.3 トーラス表面積評価
次に，トーラスの表面積評価について検討した．トーラスの諸元は，Fig.7 に示される通りとした．このとき，トーラスの体積 Vt 及び表面積 St は次式で求められる．
Vt = 2π 2 r2 R

St = 4π 2 rR

(17)

た だ し ，本 研 究 で は r = 0.1, R = 0.35 と した ．計 算 領 域 は 1 × 1 × 1 と し ，格 子 分 割 数 は ，
43 , 83 , 163 , 323 , 643 , 1283 , 2563 とした．ボクセル化に伴う体積誤差を Tab.6 に示す．球体，回転楕円体のケースと同様，格子数の増加とともに精度の向上が認められており，適切にボクセルデータへの変換が行われていると考えられる．

Tab. 4: EL1 of spheroid volume and number of cells

次に，Tab.7 に，Case A と Case B による表面積評価の比較を示す．
Case A においては，回転楕円体のケースと同様ある解像度までは一定の精度向上を示すが，その後は逆に誤差が大きくなる傾向を示した．
Case B においては回転楕円体のケースとは異なり，一定解像度までは格子分割数の増加とともに精度の向上がみられるが，その後は同様の精度で推移した．
混相流れ計算手法
本研究では，MICS(Multiphase Incompressible flow
solver with Collocated grid System ) 3, 4 を用いて混相流流体計算を行った．基礎式は以下に示す質量保存則，非圧縮性条件式，及び運動方程式である．

4.

∂ρ
∂
+
(ρuj ) = 0
∂t
∂xj

(18)

Tab. 6: EL1 of torus volume and number of cells

Number of cells

EL1

Number of cells

EL1

43

-0.86 (%)

43

2.04 (%)

38
163-0.11 (%)
0.011 (%)

83
163

0.78 (%)
0.29 (%)

323
643
1283

-0.0036 (%)
0.0074(%)
0.00088 (%)

323
643
1283

0.054 (%)
0.023 (%)
0.0076 (%)

2563

0.00021 (%)

2563

0.0025 (%)

4Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-1

Tab. 7: Comparison of EL1 of surface area of torus

Tab. 8: Computational condition (Two bubbles)

Number of cells

EL1 (Case A)(%)

EL1 (Case B)(%)

Domain5 × 5 × 20[cm3 ]

43

-44.7

-2.46

3

8
163
323

13.4
8.28
-0.29

1.08
0.24
1.43

Number of cells
Density of liquid

80 × 80 × 320
31000[kg/m ]

643
1283

-5.68
-5.14

0.46
0.36

Kinematic viscosity of liquid
Density of gas
Kinematic viscosity of gas

1.289×10−5 [m2 /s]
3
1.0[kg/m ]1.289×10−4 [m2 /s]

3

-6.7

0.43

Surface tension
Initial location of bubble

10 [mN/m](2.5,2.5,4.0)[cm]Initial location of bubble
Initial diameter of bubble
Time interval(2.5,2.5,6.3)[cm]
1.875[cm]2.5×10−4 [s]Number of time steps

2400

∂uj
=0
∂xj
∂
1 ∂p
∂ui
+
(ui uj ) = fi −
∂t
∂xj
ρ ∂xi
]
1 ∂ [ ∂∂+(µui ) +(µuj ) · · · (i = 1, 2, 3)
ρ ∂xj ∂xj
∂xi

(19)

2.2

(20)

Case A
Case B2

ここに，t は時間，xi は３次元直交座標系の座標成分である．ρ, µ, p は計算セルにおいて体積平均された密度，
粘性率，圧力である．ui は計算セルにおいて質量平均された xi 方向の流速成分である．fi は xi 方向の外力である．式 20 をコロケート格子上で有限体積法により離散化し，MAC 系解法に基づいて予測段階，圧力計算段階，
修正段階を経て次時間ステップの流速及び圧力を更新する．予測段階においては，C-ISMAC 法による完全陰解法を用い，圧力計算段階においては，C-HSMAC 法を用いた．次時間ステップにおける流速及び圧力を更新した
後，VOF 関数の移流計算を５次精度の TVD スキームを用いて行った．上記移流計算後，非数値拡散フィルタ 6を適用し，次時間ステップにおける VOF 関数を得た．表面張力は，文献 5 に示される方法にて評価した．
上昇気泡の計算結果
以下に，上昇気泡の数値計算を行った場合の気液界面の表面積変化を追った．気泡が 2 個の場合及び 13 個の場合，更に約 9700 個の場合を検討した．なお，これらの上昇気泡の数値計算は，京都大学学術情報メディアセン
ターのスーパーコンピューターを用いて実施した．

1.8

S*

256

1.6

1.4

1.2

1
0

0.1

0.2

0.3

0.4

0.5

time(sec)

Fig. 8: Time histories of dimensionless area of interface
in bubble collision5.

5.1 気泡衝突に伴う表面積変化追跡
計算条件は，Tab.8 に示す通りであり，重力加速度を9.8[m/s2 ] とした．得られた気泡衝突のスナップショット
を Fig.11 に示す．
また，気液界面の無次元表面積の時系列を Fig.8 に示
す．気液界面の表面積は，初期時刻における気液界面積で除して無次元表面積 S ∗ をプロットした．下部に位置する気泡は，t = 0.09[s] においてその上部が持ち上がり，
t = 0.12[s] から t = 0.15 にかけて上部の気泡の下部に接触し，最終的に結合する．結合後は，平たい円盤状の形状となり，上昇を続ける．Fig.8 には，前節で検討したCase A および Case B の手法で得られる表面積の時系列を比較した．計算開始後，気泡の変形とともに表面積が増加し，気泡の衝突後の t = 0.15[s] 以降は，ほぼ一定の値となっており，安定した形状で上昇を続けている．相対的に精度が落ちるとされた Case A においても，Case
B とほぼ同様の傾向をつかんでいる．本計算ケースの初期値において，単一の気泡の直径に対して計算セルは約30 個の割当てである．
5.2 上昇気泡群の表面積変化追跡
次に気泡の数を増やし，13 個の上昇気泡群の数値計算を行い，表面積の時間変化を追った．計算条件は，Tab.95に示す通りであり，得られた気泡上昇のスナップショット
を Fig.12 に示す．また，気液界面の無次元表面積の時系
列を Fig.9 に示す．
気泡群は，自由水面をもつ水槽内にほぼ等間隔に配置され，計算開始とともに上昇を始める．水面に近い位置に配置された気泡は，t = 0.28[s] 頃から自由水面に達し，水面形を大きく変形させる．気泡が自由水面に達した後，
t = 0.35[s] において，気液表面積はピークを示している．
本計算ケースの初期値において，単一の気泡の直径に対して計算セルは約 13 個の割当てである．
最後に気泡の数を増やし，約 9700 個の上昇気泡群の数
値計算を行い，表面積の時間変化を追った．
計算条件は，Tab.10 に示す通りであり，得られた気泡
上昇のスナップショットを Fig.13 に示す．また，気液界
面の無次元表面積の時系列を Fig.10 に示す．気泡群は，自由水面をもつ水槽内にほぼ等間隔に配置され，計算開始とともに上昇を始める．気泡の初期配置は，等間隔配置となる位置から最大で 1cm 程度のずれをランダムに与えている．本計算ケースの初期値において，単一の気泡の直径に対して計算セルは約 10 個の割当てである．球体
表面積の検証計算の結果から，単一気泡に対して約 83 個以上の計算セルの割当を確保した．
結言
本研究においては，気液混相流れにおける界面を通しての物質交換等の現象解明に資するために計算セルに与えられた流体体積率の分布から界面面積を推定する手法について検討した．体積が等価となる球体片を仮定して

6.

Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-1

Tab. 9: Computational condition (A group of bubbles)
Domain

20 × 20 × 60[cm3 ]Number of cells
Density of liquid

96 × 96 × 384
31000[kg/m ]Kinematic viscosity of liquid
Density of gas
Kinematic viscosity of gas

1.0×10−6 [m2 /s]
3
1.0[kg/m ]1.0×10−5 [m2 /s]

Surface tension
Initial diameter of bubbles

7.2 [mN/m]
1.39[cm]

Time interval
Number of time steps

2.0×10−3 [s]200Tab. 10: Computational condition (A large group of
bubbles)
Domain40 × 40 × 160[cm3 ]Number of cells
Density of liquid
Kinematic viscosity of liquid

160 × 160 × 640
31000[kg/m ]1.0×10−6 [m2 /s]

Density of gas
Kinematic viscosity of gas

1.0[kg/m ]1.0×10−5 [m2 /s]

Surface tension
Initial diameter of bubbles
Number of bubbles

10 [mN/m]
2.0[cm]
9700

Time interval
Number of time steps

5.0×10−4 [s]
180

3

2.4
Case A
Case B

2.2

2

S*

2

Case A
Case B

1.8
1.8
1.6
1.6
1.4

S*

1.4

1.2
1.2
1
1
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4
0.8

time(sec)
0.6
0Fig. 9: Time histories of dimensionless area of interface
in rising bubbles

0.01

0.02

0.03

0.04

0.05

time(sec)

Fig. 10: Time histories of dimensionless area of interface in rising bubbles
表面積を推定する方法と，体積が等価となる切断面を仮定して表面積を推定する方法について検討した．それぞれの手法の精度や特性を確かめるために球体，回転楕体，
およびトーラス形状について表面積の評価を行い，理論値との比較を通して精度の検証を行った．この結果，等価切断面を仮定する方法の方が精度がよく，格子数の増加に対する精度の向上が比較的よく確認できた．等価球体片を仮定する方法では，一定の解像度以上になると逆に精度が悪くなる傾向を確認した．
最後に気泡群上昇の数値計算を行い，提案手法による気泡の表面積評価を行った．今後は，界面付近の流速や物質濃度勾配の評価法について検討する予定である．

案. 土木学会論文集, No. 740/II-64, pp. 121–130, 2003.

4. 牛島省, 山田修三, 藤岡奨, 禰津家久. 3 次元自由水面流れによる物体輸送の数値解法 (3D MICS) の提案と適用性の検討. 土木学会論文集, Vol. 810/II-74, pp. 79–89,

2006.
5. S. FUJIOKA and S. USHIJIMA.

参考文献

Application ofevaluation method of surface tension with geometrical condition to three-dimension simulation of droplet
deformation and breakup. 1st International Sympo-

1. 杉原裕司, 山上路生. 界面水理学における乱流とガス交換-実験とモデルの進展-. ながれ, Vol. 30, pp. 181–193,

sium on Multiscale Multiphase Process Engineering,
No. L5, 2011.

2011.
2. Garbe C.S., Handler R.A., and Jahne B., editors.
Transport at the Air-Sea Interface. Springer-Verlag,
2007.

6. 藤岡奨, 牛島省. VOF 関数の移流計算における数値拡散抑制手法の検討. Transaction of JSCES, No. 20100020,
2010.

3. 牛島省, 竹村雅樹, 山田修三, 禰津家久. 非圧縮性流体
解析に基づく粒子 − 流体混合系の計算法 (MICS) の提6

Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-1(a) Overall (t= 0.0(s))(a) t= 0.0(s)(b) t= 0.09(s)(b) Zoom in (1)(t= 0.045(s))(c) t=0.12(s)(b) t=0.15(s)Fig. 11: Snapshots of bubble collision(c) Zoom in (2) (t = 0.045(s))
Fig. 13: Snapshots of a large group of rising bubbles

7

Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
E11-1(a) t= 0.0(s)(b) t= 0.1(s)(c) t= 0.20(s)(d) t= 0.28(s)(e) t = 0.30(s)(f) t = 0.35(s)Fig. 12: Snapshots of rising bubbles

8

Copyright c 2011 by JSFM