第26 回数値流体力学シンポジウム
講演番号 B03-2

保存形拡散界面モデルに基づく二相流体移流方程式の数値解法
A Numerical Method for Two-phase Fluid Advection Equation
Based on a Conservative Diffuse-interface Model○ 高田 尚樹, 産総研, 〒305-8564 茨城県つくば市並木 1-2-1, E-mail: naoki-takada@aist.go.jp
松本 純一, 産総研, 〒305-8564 茨城県つくば市並木 1-2-1, E-mail: matsumoto-junichi@aist.go.jp
松本 壮平, 産総研, 〒305-8564 茨城県つくば市並木 1-2-1, E-mail: sohei.matsumoto@aist.go.jp
Naoki Takada, AIST, 1-2-1 Namiki, Tsukuba, Ibaraki 305-8564, Japan
Junichi Matsumoto, AIST, 1-2-1 Namiki, Tsukuba, Ibaraki 305-8564, Japan
Sohei Matsumoto, AIST, 1-2-1 Namiki, Tsukuba, Ibaraki 305-8564, Japan
In this study, Allen-Cahn (AC) and Cahn-Hilliard (CH)-type advection equations based on a diffuse-interface model
are tested through interface-translation benchmark problems for immiscible two-phase fluid flow simulations. A
lattice-Boltzmann method (LBM) based on fictitious particle kinematics is proposed for solving a corrected AC
equation, and the numerical result is examined with those obtained by using finite volume methods (FVMs). The major
findings are as follows: (1) the fluid volume and shape is well conserved with use of 4th-order central finite difference
and 2nd-order Runge-Kutta schemes for both the equations; (2) the conservations is done by a local diffusion inducedby chemical potential gradient with a correction term to eliminate interfacial-curvature effect; (3) the result by LBM is
better than those by 2nd-order FVM; (4) the AC equation including a 2nd-order differential term, one-step conservative
level-set equation, reduces computational cost compared with the corrected CH equation including a 4th-orderdifferential term.

(1)

ルギー理論に基づく DIM（またはPhase-field model） を採用する．
本 DIM は，二相の系が持つ自由エネルギーの汎関数の値を最小化するよう，以下の化学ポテンシャルの勾配に起因する物質移動を起こし，臨界点以下の条件で0.5 に比例する幅の界面（の急激な遷移領域）を形成させる(1)．

１．はじめに
マイクロスケールで気相と液相が二種類以上混在する混相流
体の流れ場では，粘性応力や表面張力・界面張力の効果が重力の効果よりも大きくなる．加えて，気液・液液界面が固体表面に接する場合には各流体と固体間の相対的な親和性（濡れ性）も流れに強く影響する．濡れ性が不均一な固体表面上を流体界面が移動する接触線問題に対して，マクロスケールの観点から流体界面を厚さ無しの移動境界（面）に単純化する数値流体力学（CFD）解析では，高度な幾何学演算を使うため界面の移流・構築アルゴリ
ズムが複雑化する．
上記課題を克服するアプローチの一つとして，非平衡系の自由エネルギーの考え(1)を援用した拡散界面モデル（Diffuse- interface
Model, DIM）を導入する CFD 法が国内外の研究者によって提案され，様々な混相流体問題に適用されている(2)-(12)．材料科学分野の Phase-field（PF）法同様，DIM 法も，平衡状態に向けて流体相を自律的に形成する特徴により，幾何学的大規模変形を伴う多数の界面の移動・形成計算を簡素化できる．本報では，濡れ性が強く影響する MEMS デバイス作製プロセス(13)に関する二相流 CFD
シミュレーション法の開発を目的とした，DIM の 2 種類の界面移流方程式の概要と特徴を述べるとともに，簡単な界面移流問題のベンチマークテスト計算を通してそれら方程式を解くための計算スキームの精度を評価する．また，格子ボルツマン法（Lattice
Boltzmann Method, LBM）(14),(15)に基づく計算スキームと，従来のDIM の界面移流方程式の修正方法を提案する．

  
   2  

(1)

 はの二重井戸型関数，第 2 項は界面領域での自由エネルギー増加を表す(1)．
本研究の DIM は離散化された計算領域内で流体界面を追跡するため使用される．その界面領域の厚さは，実在流体界面のそれとは異なり，数値的に設定されるものであり，物理的な意味を持っていない．
３．DIM に基づく移流方程式
DIM 法では各時刻 t における流れ場の流体相および界面の空間分布を表す指標関数 を次の式を解いて得る(3),(4),(6),(8),(9)．
３．１ Cahn-Hilliard（CH）型方程式
本式は保存系に対する時間依存 Gintzburg-Landau（TDGL）方程
式に流速 u の移流項を加えたものである．


1    u     ( )  Lr ur       Pe1


t





(2)

t，u 等は代表長 Lr と速度 ur で無次元化されている．Peは Peclet
数である．易動度()>0 を含む右辺は，平衡時には消え，非平衡
時に任意の正負の範囲の拡散係数を瞬時局所に持つことにより一定厚さの界面を構築する(4),(8)．
３．２ Allen-Cahn（AC）型方程式
非保存系に対する TDGL 方程式である AC 式は，移流項を加えると次式で表される．

２．拡散界面モデル（Diffuse-interface Model, DIM）
DIM は，異なる相と相が接触する境界である界面を「物性が連続かつ急激に変化する有限幅の体積領域」とモデル化する(1)．各相と相間に横たわる界面の領域は各々，各時刻 t に座標 x の空間
内で連続的に分布する指標関数(x, t)（または秩序変数）の異なる値で識別される．例えば，気液二相系に対して，液相を，気相を，その間のの有限厚さ領域を界面とみなす．
本研究では，気液・液液二相を対象として，非平衡系の自由エネ    u    M 
t
1

(3)

Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 B03-2
左辺第 2 項が流速 u による界面の移動を，
右辺の M は拡散に関する係数を表す．上式は非保存系であるため，空気-水，油-水のよ
うな非混和性で各相の体積が保存される二相系の流体界面追跡計算には適していない．
３．３ AC 移流方程式の修正(11),(12)以下で，文献(11),(12)で提示されている，式(3)の保存形への修正を述べる．式(1)の関数 を次式で定義する．

12

   2 1   2    M             

上式に基づく二相流 CFD 法は Advected phase-field 法(16),(17)と呼ばれる．さらに上式に式(7)を代入し変形すると最終的に保存形の次式が得られる．

   1       (14)    u   M '          t     (4)

M '  M 

式(4)と式(1)を式(3)右辺に代入すると以下を得る．

 1   1  2      u   M   2 

t




1
2




 2   




s

s は界面からの法線方向距離座標を表す．上式を s で微分すると以下を得る．
(7)

 1   1  2      n    s  s (8)n   / (9)       u   M   2 

t




 1
     M       

2


 
   
(10) M    
最終的に式(5)と同等の次式が導かれる．

k   n(16)
    u   0
t

(17)


    1    n        


(18)

ここで， は指標関数 の再構築のための計算時刻， は0.5 に相当する界面幅のパラメータを表す．修正保存形 AC 方程式(14)では，拡散項が に関する空間 2 階偏微分で与えられるため，有限厚さの界面の形成を CH 方程式(2)よりも少ない負荷で効率的に計算できると考えられる(11),(12)．
３．４ CH 移流方程式の修正
界面曲率依存の拡散流束は，AC 方程式だけでなく CH 方程式でも生じる(1)．流体界面移流計算に対するこの界面曲率の影響を取り除くため，著者らは，曲率による化学ポテンシャル変動を考慮して指標関数を初期設定すれば高精度に体積が保存される流
体界面移流を計算可能であることを示した(4),(8)-(10)．しかしながら，界面が複雑に変形し，その曲率が時間的・空間的に一定ではない
混相流体流れ場では化学ポテンシャルもまた一定でないため，曲率依存拡散流束を領域全体で常に一様に調整することは実際には実行困難な場合が多い．また，従来の CH 方程式と流体運動方程式の組み合わせでは表面張力効果がダブルカウントされるという問題もある．そこで，本研究では，上述の界面曲率依存性の問題を解決するため，修正 AC 方程式と同様に，界面曲率で誘起される化学ポテンシャル変動を打ち消す効果を式(2)右辺の拡散項に付加した，以下の修正 CH 方程式を提案する．

式(8),(9)を式(5)右辺に代入して以下の変形操作を行う．

    u   M   kt1
1  exp   s /  



(6)

 1     
s(15)

上式は，定式化が若干異なるものの，Conservative Level-set (CLS)
法(18)-(20)の以下の式と同等である．

(5)

ここで，式(4)に対して平衡状態（ = constant.）にある平坦な界面に対する は以下で表される．

  1  tanh 

(13)

(11)


    u      M      k  
t

(12)

(19)

右辺の M は拡散流束の大きさを決める易動度である．

式(12)の k は界面曲率を表す．上記から，AC 式(3)（または(5)）
では，曲がった界面は流速 u に加えて k に依存した局所速度で動くことが分かる．
混相流中では曲率依存の界面速度が表面（界面）張力により生じるが，二相流 CFD 法の多くは運動方程式を通して流速 u に表
面張力効果を反映させるため，界面移流方程式には界面曲率に依存しない界面形成メカニズムが良い．そこで，AC 式(5)から曲率依存の界面移動を取り除くため同式右辺から式(10)右辺を引く(11)．

４．計算スキーム
本研究では，
修正保存形 AC 方程式(14)と修正 CH 方程式(19)を，有限体積法（FVM）及び格子ボルツマン法（Lattice Boltzmann
Method, LBM）(14),(15)に基づく計算スキームで解く．
４．１ 修正 AC 方程式の LBM
本節では，修正 AC 方程式(14)を解くための LBM を新たに提案
する．LBM は, 等方的離散速度 ea を持つ仮想の流体粒子の集合に対する粒子数密度分布関数 ga の時間発展方程式(20)を解く(14),(15)．

       u   M          
 
t

   


 ga
1
(20)
 ea g a   g aneq
tg
右辺上付添字neq は瞬時局所の非平衡量,  g はBGK 近似の衝突緩2Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 B03-2
修正 AC 式に対する上記の平衡分布関数 gaeq は，j = 0 の場合，
の移流拡散方程式に対する従来の関数と全く同一である．
よって，
本 LBM を採用すれば，二相流界面追跡計算コードは，LBM ベースの移流拡散計算コードの部分的拡張によって容易に開発できる．
今回の計算では，式(27)及び(28)で = 0.01 を与えた．また，
Semi-Lagrange 型 LBM 固有の 4 階空間微分項による数値粘性を打ち消すため，式(21)の衝突緩和時間g を次の値に設定した(21)-(24)．

和時間を示す．右辺は上式が式(14)と等価になるような平衡分布
本研究では，最も一般的な LBM に従って, 式(20)関数 gaeq を含む．
を空間・時間ともに 2 次精度の Semi-Lagrange 形式に離散化する．

g a  x  e a t0 , t  t0   g a  x, t t0 g  x, t   g aeq  x, t  g  a(21)

上式で，t0 =1 は時間刻み幅で，粒子速度 ea には次の 2 次元 9 速度（D2Q9）モデルを採用する．

e0 , e1 , e2 , e3 , e4 , e5 , e6 , e7 , e8 
 0 1 1 0 1 1 1 0 1  c 0 0 1 1 1 0 1 1 1

12

g  

(22)(23)

式(22)で与えられる 9 種類の流体粒子集団は，2 次元デカルト座標系(x, y)空間を一様分割する正方形構造格子に沿って格子点間を移動し，各格子点上で相互に衝突して各集団の粒子数密度 ga が瞬時
局所の平衡状態 ga eq に向かって変化する．
ここで，式(14)を以下のように書き換える．

    u      D0   J tTable 1 The numerical scheme configurations for the corrected AC-type
interface advection equation.

(24)

Schemetype #

右辺の D0 は修正 AC 方程式の拡散係数の定数部分を表す．

D0  M 

(25)

1
2
3
4
5
6
7
8
9
10
11
12
13

J は以下で表されるの拡散流束である．

J j M   1   


(26)



平衡状態で理想的にはが一定値 0 または 1 で一様となる各相内式(24)右辺の拡散項は消える．
部では，
= 0 およびJ = 0 の結果，修正 AC 方程式(24)を LBM の時間発展式(21)から導出するための平衡分布関数 gaeq は次式で与えられる．

Spatial gradient Interpolation at
Advection（for interfacial
spatial cellterm
curvature）
surface
2nd
2nd
2nd
2nd
2nd
2nd
2nd
4th
4th
2nd
4th
2nd
2nd
4th
2nd
2nd
4th
4th
2nd
4th
2nd
2nd
2nd
4th
2nd
2nd
4th
2nd
2nd
2nd
2nd
4th
4th
4th
4th
4th
4th
4th
4th

2

e  (u  j) 3  ea  u  u  u  2
g aeq  x, t   3 wa   a 2
2c 4
2c 
c


(27)


3 u  u g  x, t    1  3 1  w0   w0

2
c2 


(28)

a

a

2nd
2nd
4th
2nd
2nd
2nd
4th
4th
2nd
4th
4th
4th
2nd

Time
advance（explicit）
1st
2nd
2nd
2nd
1st
2nd
2nd
2nd
2nd
2nd
1st
2nd
2nd

,   

y

y/2
y/2

y

y

上式中の重み係数 wa には，静止粒子(a=0)に対して 4/9，運動粒子に対して 1/9 (a = 1, 3, 5, 7)，1/36 (a = 2, 4, 6, 8)が適用される．

 w 1

Diffusion
termx/2 x/2(a  0)
eq
0

(33)

2 3

この条件は，
修正CH 方程式のLBM 計算スキームにも適用する．
４．２ 計算スキーム
本研究では，上記 LBM 計算スキームの他，Table 1 に示す修正AC 式(14)または(24)の各偏微分項に対する 1 次～4 次精度の有限差分近似を組み合わせた，有限体積法（FVM）に基づく 13 種類の陽的計算スキームを使用して式(14)または(15)を数値的に解く．
スカラー・ベクトル各変数は，Fig. 1 に示されるように，正方形セルの中心とセル境界上へスタガード状に各々配置される．時刻t は，一定刻み幅t = 1 で Euler 法もしくは 2 段階 Runge-Kutta 法により進行させる．修正 CH 式(19)には，FVM の他，Swift らによる
二相流体LBM 計算スキーム(25)のD2Q9 モデル拡張版を適用した．

ここで c は格子定数であり，時間刻み幅t0 及び 2 次元デカルト座標系(x, y)空間の分割幅x=y=1 と以下の関係を満たす．

ct0  x  y

1

O

(29)

(30)

  u  j   a g aeq ea

(31)

t 

D0  c 2   g  0 
2 


(32)

x

x

u,


x

v,


y

Vector components(on cell surfaces)

Fig. 1 The staggered arrangement of scalar and vector variables on a
spatial structured grid with square mesh in FVM.

また，は静止平衡状態での静止粒子と運動粒子の数密度の割合を任意に調整するための正値のパラメータである．
上述の LBM では，指標関数(x, t)，流速 u および拡散係数 D0
は以下のように定義される．

   a g a   a g aeqx

x

Scalar variables(at cell centers)

５．界面移流ベンチマーク計算
本節では，DIM の流体界面形成・輸送特性を把握するため行った，修正 AC 型方程式(14)及び修正 CH 型方程式(19)による２次元
円形界面の非定常線形移流問題のベンチマーク計算を述べる．
計算領域は，Fig. 2 に示すように，２次元デカルト座標系(x, y)の下，単位幅x = y = 1 の正方形構造格子で一様に等分割し，周囲の境界に周期条件を適用した．流速ベクトル u の x, y 方向成分
u 及び v は，領域全体に一様で相互に等しく一定 U=V に設定された．x 方向移動速度 U は Courant 数 C = Ut0/x = 0.01 で与えた．
3

Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 B03-2
Figs.3 及び 4 から，修正 AC 式への FVM の適用では，最も離散化精度の高い Type12 が長時間に渡って面積保存性に最も優れ，界面の初期形状を良好に保持することが確認できる．
Type3 と Type6でも面積と界面形状は他の Type よりも良好に保持されているが，
Type12 と比べると若干十字型に歪んでいる．それ以外のスキームの結果では面積と界面形状の初期との差がより顕著に現れている．
次に，Fig.5 において，本研究で提案する修正 AC 式の LBM 計
算スキームの結果を上述の FVM の代表的な結果と比較する．な
お，LBM では式(23)中のの空間勾配を 4 次精度中心差分近似で評価した．面積保存性に関しては，LBM は最も高精度の差分近似を採用する FVM（Type12）より劣るものの，一部 2 次精度の FVM（Type6）よりも良好な結果を与えており，界面形状は LBM では
Type16 よりも改善されている．LBM では空間 2 次精度で移流が計算されるが，速度空間の等方的離散化により空間格子の対角線方向に移動する流体粒子の速度分布関数が考慮されるため，一部
4 次精度の FVM よりも高精度な界面移流計算が実現されると考えられる．

指標関数 (x,0)の初期空間分布には，直径 d=32x の円形領域内部
で=1，外部で=0 と仮設定した後に式(14)または式(19)の静止状
態（u=0）での時間発展により得られる定常解を与えた．また，両式に易動度 M=1，界面厚さパラメータ =1 を適用した．
Figs. 3～7 には，修正 AC 方程式(14)及び修正 CH 方程式(19)へ
FVM 及び LBM 計算スキームを適用して得られた，指標関数(x,t)の空間分布及び円形領域面積 A の時間履歴を示す．なお，A は =0.5 以上の値を持つ空間セルの数で評価した．

Lx  100xTime t = 0

4xPeriodic
boundaries

4 y

100y

U＝(U,V)
U=V0

 (x,t)[10 ] 0.5-2

1

Circular-shaped Area
Conservation Error (A/A0-1)yx  y  1

d  32x
O x A contour line of  =0.5

[10 ]-2

0.5

Circular-shaped Area
Conservation Error (A/A0-1)

Fig. 2 The computational domain and initial condition in the benchmark
test simulation of linear translation of an interface.

0

C = Ut/x = 0.01
Type 2 Type 3 Type 4 Type 6

-0.5

-1.5
-2.0

-1.0
-1.5

Type 6

Type 12

LBM

-2.0
-2.5
0

10
20
30
40
Dimensionless time t* = tU/d[10-2] 0.5

0

10
20
30
40
Dimensionless time t* = tU/dCircular-shaped Area
Conservation Error (A/A0-1)

Type 7 Type 11 Type 1250Fig. 3 Time history of conservation error of the circular-shaped area inthe linear translation with use of the revised AC equation.

(a) Type 2

-0.5

50

Fig. 5 Comparison of conservation errors of the circular-shaped area inthe linear translation with use of FVMs and LBM for the revised
AC equation.

-1.0

-2.5

0

(b) Type 3(c) Type 4

C = Ut/x = 0.01

0
-0.5
-1.0
-1.5
-2.0
-2.5
0

LBM
2nd-order CDS
4th -order CDS
Type 12FVM
Type 6

Type 13

10
20
30
40
Dimensionless time t* = tU/d50

Fig. 6 Time history of conservation error of the circular-shaped area inthe linear translation with use of the revised CH equation.

(d) Type 6(a) Type 7(b) Type 11(c) Type 12
Fig. 4 The two-phase distributions at t*= tU/d = 50 in the linear
translation with use of the revised AC equation for C =Ut0/x=
10-2.

(a) FVM Type 6 (b) FVM Type 12 (c) LBM 2nd order (d) LBM 4th orderFig. 7 The two-phase distributions at time t*= tU/d = 50 in the linear
translation with use of the revised CH equation for C = Ut0/x =10-2.
4Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 B03-2
本研究で提案される修正 CH 方程式を用いて得られた結果をFigs. 6～7 に示す．計算スキームには，Table 1 の Type6，Type12
及び Type13 の FVM3 種類，およびの空間勾配を 2 次または 4次精度中心差分で近似する 2 種類の LBM を採用した．Fig. 6 に示されるように，
各 LBM スキームはほぼ同精度の FVM よりも少ない誤差で初期面積を持続している．また，界面形状についても，
LBM は FVM よりも良い結果を与えている．
最後に，参考として，修正 AC 及び修正 CH 両方程式に基づく（時間ステップ数160,000）
界面移流計算で無次元時刻t* = t U/d =50までに要した実際の演算時間（CPU Time）を Table 2 に示す．使
用した計算機は，SONY VAIO VGN-TT91DS（OS: Microsoft
Windows Vista Business 64bit 版／CPU: Intel Core 2 Duo SU9400 (拡張版 Intel SpeedStep テクノロジー, 1.4GHz, 3MB 2nd cachememory)／Main memory: 4GB (DDR3 SDRAM, PC3-6400, 800MHz,
デュアルチャンネル転送)）である．修正 AC，修正 CH 各式において，LBM は僅かではあるが FVM よりも短い時間で計算を実施した．一方，両式での演算時間を比較すると，LBM は予想通り 2階偏微分項を持つ AC 式に対してより短時間の演算を実現する一方，FVM は逆に 4 階偏微分項を持つ CH 式に対して AC 式よりも時間を短縮する結果となった．今後，より厳密に計算機環境を整えた上で演算時間を計測して定量的な比較を行う予定である．

謝 辞
上記研究の一部は，総合科学技術会議が制度設計した最先端研究開発支援プログラム（FIRST）により，採択課題「マイクロシ
ステム融合研究開発」（江刺プロジェクト）への日本学術振興会（JSPS）による助成の下で実施されたものである．関係者各位への謝意をここに表します．
参考文献(1) Anderson, D. M., McFadden, G. B. and Wheeler, A. A.,“Diffuse-interface methods in fluid mechanics,” Annu. Rev. Fluid
Mech., Vol. 30 (1998), pp. 139–165 (DOI: 10.1146/annurev.fluid.30.
1.139).
(2) Jamet, D., Lebaique, O., Coutris, N. and Delhaye, J.M., “The second
gradient method for the direct numerical simulation of liquid-vapor
flows with phase change,” J. Comput. Phys., Vol. 169 (2001), pp.
624–651 (DOI: 10.1006/jcph.2000.6692).
(3) Inamuro, T., Ogata, T., Tajima, S., and Konishi, N., “A lattice
Boltzmann method for incompressible two-phase flows with large
density differences,” J. Comput. Phys., Vol. 198 (2004), pp. 628–644(DOI: 10.1016/ j.jcp.2012.01.019).
(4) Takada, N. and Tomiyama, A., “A Numerical Method for Two-Phase
Flow Based on a Phase-Field Model,” JSME Int. J. Ser. B - Fluids
Therm. Eng., Vol. 49 (2006), pp. 636–644 (DOI: 10.1299/jsmeb.49.
636).
(5) Seta, T. and Okui, K., “The single component thermal lattice
Boltzmann simulation of pool boiling in two dimensions,” Journal of
Thermal Science and Technology, Vol. 1 (2006), pp. 125–137 (DOI:
10.1299/ jtst.1.125).
(6) Yoshino, M. and Mizutani, Y., “"Lattice Boltzmann simulation ofliquid-gas flows through solid bodies in a square duct,” Math.
Comput. Simul., Vol. 72 (2006), pp. 264–269 (DOI: 10.1016/
j.matcom.2006.05.018).
(7) Onuki, A., “Dynamic van der Waals theory,” Phys. Rev. E, Vol. 75(2007), 036304 (DOI: 10.1103 /PhysRevE.75.036304).
(8) Matsumoto, J. and Takada, N., “Two-phase flow analysis based on aphase-field model using orthogonal basis bubble function finite
element method,” Int. J. Comput. Fluid Dyn., Vol. 22 (2008), pp.
555–568 (DOI:10.1080/10618560802238226).
(9) Takada, N., Matsumoto, J., Matsumoto, S. and Ichikawa, N.,“Application of a phase-field method to the numerical snalysis of
motions of a two-phase fluid with high density ratio on a solid
surface,” Journal of Computational Science and Technology, Vol. 2(2008), pp. 318–329 (DOI:10.1299/jcst.2.318).
(10) Matsumoto, J., Takada, N. and Matsumoto, S., “One hundred million
degree of freedom two-phase flow finite element method analysis
using phase-field model,” Proceedings of NCTAM, Vol. 60 (2011),
Paper No. 251.
(11) Chiu, P-H. and Lin, Y-T., “A conservative phase field method forsolving incompressible two-phase flows,” J. Comput. Phys., Vol. 230(2011), pp. 185–204 (DOI: 10.1016/j.jcp.2010.09.021).
(12) Oshima, N. and Gong, J., “Level-set and phase field approaches forinterface of conservative system,” Proc. the 24th JSME
Computational Mechanics Conference (JSME-CMD2011) (2011),
CD-ROM, Paper No. 915.
(13) Lu, J., Takagi, H. and Maeda, R, “Chip to wafer temporary bonding
with self-alignment by patterned FDTS layer for size-free MEMS
integration,” Proc. IEEE Sensors 2011 conference (2011), pp.
1121–1124.

Table 2 The numerical scheme configurations for the revised AC-type
interface advection equation.
Run#

Allen-Cahn (AC) type
FVM
LBM (4th(Type9) order gradient)297.8527
279.6942
298.8667
279.3042
298.1179
279.741298.5547
279.7566
300.2239
280.2402

1
2
3
4
5
AverageCPU
298.7232
Time
Ratio to1.000
AC-FVM

Cahn-Hilliard (CH) type
FVM
LBM (4thLBM (2nd(Type9) order gradient) order gradient)283.1106316.994290.0995
280.5522
316.136284.811280.8018
315.7772
283.3914
281.6598
315.8864
275.5446
283.953
317.384266.5433

279.74724

282.0155

316.43552

280.07796

0.936

0.944

1.059

0.938

６．まとめ
本報では，MEMS デバイス作製プロセスに関連する，固体表面の濡れ性が強く影響する気液・液液二相流動を数値的に予測するための拡散界面モデル（DIM）に基づくマイクロ二相流界面追跡計算法の開発を目的として，
従来 DIM の拡散流束から界面曲率依存性を取り除き保存形に修正した Allen-Cahn（AC）
型（Conservative
Level-set）界面移流方程式を概説した．また，AC 式と同様に拡散流束を修正した Cahn-Hilliard（CH）型移流方程式と，格子ボルツ
マン法（LBM）に基づく修正 AC 方程式の計算スキームを提案した．そして，LBM と有限体積法（FVM）を使用して，簡易な 2次元非定常円形界面線形移流問題のベンチマークテスト計算を実施し，両修正方程式が有する界面形成・輸送特性を評価した．
本ベンチマークテスト計算結果から，以下の知見を得た．
(1) 修正AC及び修正CH両方程式による流体界面移流計算では，界面形状と流体面積（体積）が高精度に保持され得る．
(2) LBM は，同等の差分近似精度を持つ有限差分法（FVM）よりも良好に界面形状と体積を保持する．
(3) 修正 AC 方程式は，修正 CH 方程式と同等の界面形状維持と体積保存の特性を有しつつ，同 CH 式よりも計算負荷（演算時間）を低減させる可能性を有する．
今後，修正 AC 界面移流方程式に対する上記 LBM を運動方程式の LBM と融合し，動的濡れ性の境界条件も検討した上でMEMS に関するマイクロ気液・液液系二相流の CFD シミュレー
ションに適用する予定である．
5

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号conservative level set method,” J. Comput. Phys., Vol. 231 (2012), pp.
3887–3895 (DOI: 10.1016/j.jcp.2012.01.034).
(21) Qian, Y. H., and Orszag, S. A., ”Scalings in reaction A+B→C:
Numerical simulation by lattice Boltzmann models”, J. Stat. Phys.,Vol. 81 (1995), pp. 237–253 (DOI: 10.1007/BF02179978).
(22) Orlandini, E., Swift, M.R., and Yeomans, J.M., “A lattice Boltzmannmodel of binary-fluid mixtures,” Europhys. Lett., Vol. 32 (1995), pp.
463–468 (DOI: 10.1209/0295-5075/32/6/001).
(23) Hirabayashi, M., Chen, Y., and Ohashi, H., “The lattice BGK modelfor the Poisson equation,” JSME Int. J., Ser. B, Vol. 44 (2001), pp.
45–52 (DOI: 10.1299/jsmeb.44.45).
(24) Takada, N. and Misawa, M., “Examination of Phase Field Model andLattice Boltzmann Scheme for Interface-Tracking Simulation of
Two-Phase Flow,” Proc. JSFM 18th Computational Fluid Dynamics
Symposium (2004), Paper No. D4-2, CD-ROM.
(25) Swift, M. R., Orlandini, E., Osborn W. R. and Yeomans, J. M.,“Lattice Bolztammn Simulations of Liquid-Gas and Binary Fluid
Systems,” Phys. Rev. E, Vol. 54 (1996), pp. 5041–5052 (DOI:
10.1103/PhysRevE.54.5041).

(14) Chen, S. and Doolen, G. D., “Lattice Boltzmann method for fluid
flows,” Annu. Rev. Fluid Mech., Vol. 30 (1998), pp. 329–364 (DOI:
10.1146/annurev.fluid.30.1.329).
(15) Succi, S., “The lattice Boltzmann equation for fluid dynamics andbeyond,” (2001), Oxford at the Clarendon Press, UK.
(16) Biben, T. and Misbah, C., “Tumbling of vesicles under shear flow
within an advected-field approach,” Phys. Rev. E., Vol. 67 (2003),
031908 (DOI: 10.1103/PhysRevE.67.031908).
(17) Beaucourt, J., Biben, T., Leyat, A. and Verdier, C., “Modeling breakup
and relaxation of Newtonian droplets using the advected phase-fieldapproach,” Phys. Rev. E., Vol. 75 (2007), 021405 (DOI: 10.1103/
PhysRevE.75.021405).
(18) Olsson, E. and Kreiss, G., “A conservative level set method for twophase flow,” J. Comput. Phys., Vol. 210 (2005), pp. 225–246 (DOI:
10.1016/j.jcp.2005.04.007).
(19) Tan, N., Aoki, T., Inoue, K. and Yoshitani, K., “Numerical simulationof two-phase flow driven by rotating object,” Transactions of theJapan Society of Mechanical Engineers, Series B, Vol. 77 (2011), pp.
1699–1714 (DOI: 10.1299/kikaib.77.1699).
(20) Sato, Y. and Niceno, B., “A new contact line treatment for a6

Copyright © 2012 by JSFM