VOF 法による気液二相流れと熱伝達の数値シミュレーション
足立高弘

人見健太

宮崎大輔 1

秋田大学工学資源学部機械工学科

1

緒言
現在，地球温暖化対策や化石燃料の枯渇問題などエネルギーの有効利用は人類の活動において必要不可欠な課題となっている. その一方で，生活環境の快適化への志向はますます強まり，環境保全対策の一つとして空調機器の高効率化に対するニーズが増大している. 特に，ビルやホテルなどの冷暖房設備として，吸収式冷凍機・冷温水機が注目されている. この吸収式冷凍機・冷温水機は，発電設備や工場等の排熱を利用できることから省エネルギー効果が高く，また冷媒にフロンガスではなく水やアンモニアなどの自然冷媒を使用するため地球環境にもやさしいという利点がある. しかしながら，熱源駆動の熱交換システムであるため装置が大きくなることや起動性能が悪いなどの欠点を持ち，さらなる改良が求められている.
吸収式冷凍機・冷温水機は蒸発器，凝縮器，再生器および吸収器で構成されており，その中でも吸収器の改良は装置の小型化，低コスト化，そして性能向上に繋がるた
め重要である. 吸収器には大きく分けて管群式吸収器とプレート式吸収器がある. プレート式吸収器は，管群式吸収器より単位容積当たりの伝熱面を増加することが容易であるため高効率化が期待されている. プレート式吸収器では，伝熱面の片側を吸収溶液の薄い液膜流が流れており，周囲の冷媒蒸気を吸収する．同時に，冷媒蒸気が吸収液に吸収される際に放出される潜熱や希釈熱を効率よく除去することが必要である．
このように，吸収器は熱と物質の移動が共存する複雑な現象となっており詳細な研究が少なく，熱・物質伝達のメカニズムには不明な点が多い.
この吸収器の高性能化を図る方法としては伝熱面を流下する液膜を可能な限り薄くすることあるいは液膜内の乱れを増大させることが重要となる [1][2]. その代表的な方法として，プレートに凹凸などの乱れ促進体を設置する方法が用いられる．しかし，液膜の厚みが 0.1mm 程度と薄いために実験による計測がきわめて困難であり，数値解析によるところが大きい．そこで本研究では，プレート式吸収器に着目し，平滑面や凹凸のあるプレートを流下する液膜流の流れと熱移動を数値的に解析できるプログラムを開発する．

1 現在は，株式会社前川製作所に勤務．

1

2問題の定式化

2.1

支配方程式

h* Ǭ*
u*
y*
O*
u*
0v*
x*

wo*

wb*

wi*

0g*
図 1: 物理モデルと座標系
図 1 に，鉛直平板に矩形横溝加工を施した伝熱面の物理モデルと座標系とを示す．
座標軸は，平滑面上入口部を原点とし，鉛直下向きに x∗ 軸をとり，それに垂直な方向に y ∗ 軸をとる．矩形横溝加工を施した伝熱面形状を決定する無次元パラメータは，溝部の溝深さ h，溝幅 wb ，溝上部流入口側の幅 wi ，溝下部流出側の幅 wo である．
代表長さとして平滑平板を流れる液膜に関するヌセルトの液膜理論から得られる液膜厚さ δ0∗ を用いると伝熱面形状を決定する無次元パラメータは次式のようになる．

h∗wb∗wi∗wo∗,w=,
w=,
w=.
(1)
bio
δ0∗
δ0∗δ0∗δ0∗
ここで，変数上部の * は次元を有することを意味する．また，流路全体の長さは
h=L=wi +wb +wo である．ヌセルトの液膜理論から得られる液膜厚さ δ0∗ は次式のようになる．
 ∗ ∗ 1/33ν Q
∗
δ0 =.
(2)
g∗
ここで，ν ∗ は流体の動粘性係数，g ∗ は重力加速度，Q∗ は流量である．
本研究では液膜が薄く，プレートの奥行きが十分長い場合を考えているので，鉛直平板上の流下液膜を２次元非圧縮粘性流とする．また，気液界面では気相側からのせん断力は無視できるとする．上記の仮定を採用すれば，流れと熱移動を支配する方程
式は，連続の式，ナビエ・ストークス方程式およびエネルギー方程式であり次式のようになる．

2∂u ∂v+= 0,∂x ∂y∂u∂u∂u∂p1+ u+ v= Fr −+∂t∂x∂y∂xRe(3)∂2u ∂2u+ 2 ,∂x2∂y
 2∂v
∂v
∂v
∂p1
∂ v
∂ 2v+ u+ v=−++ 2 ,
∂t
∂x
∂y
∂y
Re ∂x2∂y
 22
∂T
∂T1
∂T
∂ T
∂ T+ u+ v=+.
∂t∂x∂yReP r ∂x2∂y 2(4)(5)(6)

上式では，平滑面を流れる薄膜流に関するヌセルトの解析解から得られる液膜厚さ δ0∗を代表長さ，その界面での速度 u∗0 を代表速度，壁面温度 Tw∗ および液膜表面温度 Ts∗
を代表温度として全ての変数を次式のように無次元化している．

x =xut∗ u∗0p∗ − p∗aT ∗ − Ts∗,u=,t=,p=,T=.
δ0∗u∗0δ0∗ρ∗ u∗2Tw∗ − Ts∗0(7)

ここで，Ts∗ および Tw∗ は一定であると仮定する．さらに，支配方程式中の無次元パラメータはそれぞれフルード数，レイノルズ数およびプラントル数であり，次式のように定義される．

Fr =δ0∗ g ∗u∗ δ ∗ν∗, Re = 0∗ 0 , P r = ∗ .
∗2u0
ν
a(8)

ここで，a∗ は流体の温度伝導率，ρ∗ は密度である．なお，F r と Re には ReF r = 2の関係がある [3].

2.2

境界条件

本研究で取り扱う問題には，気液界面境界，壁面境界，および流入・流出境界がある．以下にそれぞれの場合の境界条件を示す．
まず，気液界面での法線ベクトルを n = (n1 , n2 ) とし，位置および速度ベクトルを x = (x1 , x2 )，u = (u1 , u2 ) とすると，気液界面での応力境界条件は次式のようになる．

ps −2κBoRe
nk δik −2Re∂ui∂uknk +nk∂xk∂xi= 0,(i = 1, 2).

(9)

ここで，δik はクロネッカーのデルタ，κ は気液界面の曲率，Bo はボンド数であり次式で定義される．

Bo =δ0∗2 ρ∗ g ∗.
σ∗3(10)σ ∗ は表面張力係数である．曲率 κ は水面位置 δ が x の関数と見なせる場合には∂2δ∂x2

κ=  ∂δ 2  32
1 + ∂x(11)となる．
式 (9) を用いると，気液界面において接線方向の境界条件は次のようになる．

∂u ∂v+= 0.
∂y∂x(12)法線方向の境界条件は，法線方向が x 方向の場合には

ps = −2 ∂v2κ+BoRe Re ∂y(13)となり，法線方向が y 方向の場合には2 ∂u2κ+(14)BoRe Re ∂xとなる．このように，式 (13) と式 (14) により気液界面での表面圧力 ps が求まり，表面張力の効果が表面圧力に繰り込まれる．また，液膜表面では温度は一定として，温ps = −度条件は次式となる．

T = 0.

(15)

次に壁面上においては，流体は静止しており，壁面の温度を一定と仮定すると，境界条件は次式のようになる．

u = 0, v = 0, T = 1.

(16)

最後に流入口と流出口での境界条件は次式のように与える．すなわち，流入口での境界条件は速度および温度ともに，ヌセルトの液膜理論から得られる解析解を次式のように与える.

u = y(2 − y)，(17)v = 0,(18)T = 1 − y.

(19)

流出口では自由流出条件とし次式とする．

∂v
∂T
∂u
=== 0.
∂x
∂x
∂x

4(20)

x
vi,jui-1,j Fi,j pi,jTi,j ui,jǍyǍ

vi,j-1
図 2: スタガード格子

3

数値解析法
本研究では，ナビエ・ストークス方程式は SMAC(Simpliﬁed Marker and Cell) 法を用いて解く. その際，圧力に関するポアソン方程式の解法には，前処理付き双共役勾配安定 (Bi-CGStab) 法を用いる [4]．エネルギー方程式は時間微分項に関して２次精度のアダムス・バッシュフォース法，移流項と粘性項に関して２次精度の中心差分を用いて解く. そして SMAC 法により解かれた流速場を用いて流体充填率 F の移流を VOF(Volume of Fluid) 法により行い気液界面の挙動を表現する. 計算格子には，図 2 のような計算セルにそれぞれの面に垂直な方向の流速が定義され，流体充填率 F ，圧力 p および温度 T は計算セル中心で定義されるスタガード格子を採用する．

3.1

VOF 法

自由界面を有する流れを解析する VOF 法は，Hirt et al.[5] により開発され，その考え方を利用して様々な解析コードが開発されている [7][6]．VOF 法では界面の移動を直接計算するのではなく，流体充填率 F という量を定義して自由界面を表現する．
ここで，流体充填率 F は次式で定義される．

F =計算セル中の流体体積 (面積).
計算セルの体積 (面積)(21)

この流体充填率 F を用いると，各計算セルを次のように分類することができる.

· F = 0 : 気体セル (計算セルには流体が存在しない).
· F = 1 : 流体セル (計算セルは流体に満たされている).
· 0 < F < 1 : 表面セル (計算セルには気液界面が存在する).

5

自由界面の移動は，流体充填率 F の移流方程式を解くことにより求められる．流体充填率の移流方程式は，微小体積 (面積) に流速 (u, v) をもって出入りする流体の保存則を考慮することにより得られ，次式で表される.

∂F∂(uF ) ∂(vF )++= 0.
∂t∂x∂y(22)

この流体充填率 F の移流計算を行う際に，ドナー・アクセプター法を用いるのが VOF法の特徴の一つである [5]．

3.2

表面セルの決定と流体充填率 F の修正本研究で用いる２次元 VOF 法においては，表面セルの形状は表 1 に示される４つに分類されるとする．すなわち，表面セル内の気液界面は曲率を持たず一定面であり，その接線および法線方向は，x, y 座標軸のいずれかの方向に常に一致する．さらに，長時間にわたり安定した計算を行うための条件として，表面セルは流体セルと気体セルに挟まれていという条件を課すことにする [6][7].
さて，液膜が凹凸に衝突したり凹凸を飛び越えたりする場合には，気液界面の大きな変形を伴うことが考えられる．そのような場合には，表面セルの周囲を流体セルが一瞬にして覆う場合が考えられ，流体内部でありながら流体充填率の値が F = 1 となる場合が想定される．このような場合に，F の値を単純に F = 1 と置き換えることは流体部分の安易な増加につながり，質量保存則から大きく離れる可能性がある．
そこで，流体内部にある F < 1 である内部セルが表面セルと接している場合，２通りの修正法 (1,2) を適用し，内部セルの F を 1 にする工夫を行う．
まず，修正法１は，図 3(a)，(b) に示すように内部セルに表面セルの流体に満たされた境界面が接している場合に，表面セルをスライドさせ内部セルの F を 1 にする方法である. この場合，図 3(a) に示すように内部セルに表面セルの境界面がただ１

表 1: 2 次元流の場合の計算セルの分類

Ni

セルの状態

0

流体セル

1

表面セル : 表面が x 軸に垂直で, 流体は x 軸の負方向にある

2

表面セル : 表面が x 軸に垂直で, 流体は x 軸の正方向にある

3

表面セル : 表面が y 軸に垂直で, 流体は y 軸の負方向にある

4

表面セル : 表面が y 軸に垂直で, 流体は y 軸の正方向にある

8

気体セル

6

つの面で接している場合には，単に表面セルを内部セルの方向へスライドさせる. また，図 3(b) に示すように内部セルに表面セルの境界面が２面以上接している場合には，F が最も大きい表面セルを内部セルの方向へスライドさせる. このように修正法
１は，対象としている内部セルに表面セルの境界面が接している場合に適用可能である. 一方，図 4 に示すように内部セルに表面セルの境界面が直交している場合には，修正法１は適用できない. そこで，このような場合には，新たに修正法２を適用する.
修正法２は，内部セルが F =1 を満たすために必要な流体を周囲の表面セルから補う方法である. すなわち，内部セルの境界面に直交している表面セルから F の値が大きい方を選び，内部セルが F =1 を満たすのに必要な不足分 F を表面セルから減らし，内部セルに加えることで F を 1 にする．さらにそれらの両方の場合，すなわち内部セルに表面セルの境界面が接しているとセルと直交しているセルの両方のセルがある場合には，始めに修正法１を適用し，その後に修正法２を適用することにする．
このようにして F が 1 より小さい内部セルの F を 1 にすることでより厳密に質量保存則を満すようにする.

᳇૕࠮࡞᳇૕࠮࡞⴫㕙࠮࡞⴫㕙࠮࡞F 㧨㧝F㧩㧝

ౝㇱ࠮࡞ౝㇱ࠮࡞(a) 表面セルをスライドさせる⴫㕙࠮࡞ ᳇૕࠮࡞⴫㕙࠮࡞ ᳇૕࠮࡞F 㧨㧝F 㧩㧝ౝㇱ࠮࡞ ⴫㕙࠮࡞ౝㇱ࠮࡞ ⴫㕙࠮࡞(b)F の値が大きい表面セルをスライドさせる

図 3: 表面セルのスライド (修正法１）7ౝㇱ࠮࡞ౝㇱ࠮࡞F㧩㧝F 㧨㧝⴫㕙࠮࡞⴫㕙࠮࡞᳇૕࠮࡞ ⴫㕙࠮࡞᳇૕࠮࡞ ⴫㕙࠮࡞⴫㕙࠮࡞߆ࠄฃߌขߞߚᵹ૕ಽౝㇱ࠮࡞ߩᵹ૕ਇ⿷ಽ

図 4: 修正法１が適用できない場合 (修正法２)

Ǎyi,jpi,j ߩቯ⟵ὐ⴫㕙࿶ജߩ࿕ቯὐ

yd dci,j-1

Ǎpi,j-1ߩቯ⟵ὐ

図 5: 表面圧力の補間

3.3界面における圧力および温度境界条件

界面における圧力および温度の境界条件設定の対象となる計算セルは表面セルである．表面は前節で示した Ni により方向が決定しており，また流体充填率 F の値により表面位置が決定している．この表面位置において表面圧力や表面温度が与えられる．しかし，この表面セルにおける圧力および温度定義点は表面の位置とは一般に異なる．このため，表面セルの圧力および温度定義点と Ni が指し示す流体がある方向の流体セルの 2 点間で圧力および温度が線形に変化するものとして，表面圧力および表面温度を内挿あるいは外挿する．例えば圧力に関して表面セル (i, j) の Ni が 3 であるとき，図 5 に示した記号を用いれば次の関係式が成り立つ (温度に関しても同様の操作を行う）.

pi,j = (1 − η)pi,j−1 + ηps ,η=2∆y.
∆y + 2Fi,j ∆y(23)(24)

なお，表面圧力には表面張力が繰り込まれている．この表面張力を精度良く評価す8

るには，気液界面の曲率 κ を精度良く求める必要がある．本研究では，表面セルごとに一点ずつ界面位置を決め，それらを補間し，式 (11) により曲率を求める．界面位
置は表面セル内の界面の中央とし，その補間法として 2 次関数を用いる．このとき，壁での界面位置は，壁に最も近い界面から壁への勾配が接触角と同じになるよう定める．接触角は流体と境界面の材質により物理的に決まるが，固体表面状態である粗度や温度などで大きく変わり，表面張力，表面エネルギーを統一的に求める理論はない．
本研究では実験による測定結果を参考にし，接触角の値として 60 °を用いることにする．

4

計算結果
本研究では，流路形状パラメータ L=10，無次元パラメータ Re=30, Bo=50/Re，

P r=2 の場合について，平滑平板上を流下する液膜の流れと熱移動の数値シミュレーションを行う．

4.1

計算精度のチェック

計算領域の x および y 方向の格子刻み幅 ∆x，∆y と計算の 1 ステップ間の時間刻
み ∆t が計算結果に及ぼす影響を調べる．代表量として，時刻 t = 10 における流路入り口から流下する液膜流の先端位置 x に着目する．
表 2 に格子刻み幅を ∆x=∆y=0.1 とし，時間刻みを ∆t=2.0×10 −4，1.0×10−4，5.0×10−5 とした場合の液膜の先端位置 x の比較を示す．時間刻みを ∆t=1.0×10 −4と 5.0×10−5 とした場合では，液膜の先端位置 x の相対誤差は約 0.4% である．一方，
表 3 に時間刻みを ∆t=1.0×10 −4 とし，格子刻み幅 ∆x=∆y=0.2，0.1，0.05 とした場
合の液膜の先端位置 x の比較を示す．格子刻み幅を ∆x=∆y=0.1 と 0.05 とした場合
では，液膜の先端位置 x の相対誤差は約 1.5% である．これらの結果より，本研究における数値シミュレーションで得られた解の収束性は良好であると見なすことができる．
以下では，計算領域の格子刻み幅を ∆x=∆y=0.1，計算の時間刻み ∆t=1.0×10 −4 として計算を行う．
表 2: 時間刻み ∆t に対する液膜の先端位置 x の比較
液膜の先端位置 x

時間刻み ∆t(i)(ii)(iii)

2.0×10−46.6414

1.0×10−46.7066

5.0×10−56.6823

9表 3: 格子刻み幅 ∆x=∆y に対する液膜の先端位置 x の比較
格子刻み幅 ∆x=∆y

液膜の先端位置 x(i)0.26.2790(ii)0.16.7066(iii)

0.05

6.6050(a) t = 3(b) t = 10(c) t = 50

図 6: 平滑平板を流下する液膜流の熱流動パターン

4.2

平滑平板の場合の解析結果

鉛直平滑平板を流下する液膜流について，熱流動パターンの時間変化を図 6 に示す．
図では流体充填率が F = 0 の領域が示されている．温度分布は色で示されており，壁面で赤色，気液界面で青色になるように表示されている．図 6(a), (b) より，液膜が鉛直平滑平板に沿って流下している様子がわかる．さらに時間が経過すると，図 6(c) に見られるように，液膜流は流路出口に到達して十分発達した状態になる．このときの速度分布および温度分布は平滑平板に関するヌセルトの解析解に収束すると考えられる．そこで，鉛直平滑平板上を流下する液膜流の数値計算結果とヌセルトの解析解のy 方向の速度分布 u および温度分布 T との比較を図 7 に示す．ここで，計算結果の速度分布および温度分布は，鉛直平滑平板を流下する液膜流が十分に発達した t = 50

101.21.2Simulation result
Nusselt solution

10.8T0.8uSimulation result
Nusselt solution

10.60.60.40.40.20.20000.20.40.60.8

11.200.2

y0.40.60.8

11.2y(a) 速度分布(b) 温度分布

図 7: 速度および温度分布の比較における流路出口の値を用いている．図 7(a) は速度分布の比較を示し，図 7(b) は温度分布の比較を示す．計算結果とヌセルトの解析解との相対誤差は，速度分布では最大で約 2.2%，温度分布では最大で約 1.0% となる．これらの結果より，本研究で行った数値シミュレーションでは，十分な精度で信頼のおける解が得られていると考えられる．

5

今後の課題
本報では，プレート式吸収器の伝熱促進メカニズムを明らかにするために，平滑面や凹凸のついた平板を流下する液膜の流れと熱移動を解析するプログラムの開発状況について報告を行った．今回は，平滑な鉛直平板上を流下する液膜流れについて数値計算結果と理論解析解との比較を行い良好な結果が得られた．今後は，熱・物質移動を促進するために矩形溝を有する平板上を流下する液膜流の熱流動特性を解析する予定である．今回の計算では，格子数が 100 × 50 の比較的小規模な計算を TX7/AzusAを用いて行い，無次元時間 t = 50 (500000 ステップ）までの計算におよそ５時間の計算時間を費やした．本計算プログラムでは気相と液層が共存する気液二相流れを取り扱っており，気液界面が計算の進行につれて変化する移動境界値問題となっている．プログラム中では，この界面での位置決定や境界条件の設定等で複雑な条件分岐を行っており，そのことが並列化やベクトル化の阻害因子になっている．特に圧力に関するポアソン方程式の反復計算で，界面での境界条件を設定する必要があり，全体の約 40%の計算時間を費やしている．この部分では並列化の効率が特に悪いので，高速化が今後の課題である．また，今回は液層のみを計算し気相部分の計算は省略した．

11将来的には，気相側からのせん断力や気液界面での蒸発や凝縮を考慮に入れた計算や，３次元流れの計算が可能となるようにプログラムを改良することが課題となる．

参考文献[1] 本間・藤巻，シミュレーションによるプレート型吸収器の研究 (第１報：基本層流モデルによる解析)，エネルギー研究所技法第７号，第１報 (1997), pp.95-98.

[2] 本間・藤巻，シミュレーションによるプレート型吸収器の研究 (第２報：基本層流モデルによる解析)，エネルギー研究所技法第７号，第１報 (1997), pp.99-106.
[3] 吉永，液膜流の不安定性とカオス (トピックス 液膜流と二層流), 機械の研究，
54 巻，1 号，(2002)，pp.120-125[4] 平野，ほか２名前処理付反復法による熱流動問題の数値解析の高速化 , 化学工学
論文集，27-5，(2001)，pp.547-553.
[5] C．W．Hirt，et al., SOLA-VOF: A solution algorithm for transient ﬂuid ﬂow
with multiple free boundaries , Los Alamos Scientific Laboratory，(1980), pp.132.
[6] 山崎，自由表面を含む非圧縮性流体解析モジュール (特集 汎用 3 次元流体解析システムα-FLOW)，富士総研技報，3 巻，1 号,(1992), pp.131-146.

[7] 米山，自由表面をもつ多次元流れの数値水理学的研究，京都大学博士論文，(2001),
第 2 章.

12