物体の形状と表面濡れ性に起因する微小二相流体 ... - 日本流体力学会

第 28 回数値流体力学シンポジウム
D05-1
物体の形状と表面濡れ性に起因する微小二相流体挙動の数値シミュレーション
Numerical Simulation of Microscopic Two-phase Fluid Motion
Caused by Solid-body Shape and Surface Wettability
○ 高田 尚樹, 産総研, 〒305-8564 茨城県つくば市並木 1-2-1, E-mail: naoki-takada@aist.go.jp
松本 純一, 産総研, 〒305-8564 茨城県つくば市並木 1-2-1
松本 壮平, 産総研, 〒305-8564 茨城県つくば市並木 1-2-1
栗原 一真, 産総研, 〒305-8564 茨城県つくば市並木 1-2-1
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
Kazuma Kurihara, AIST, 1-2-1 Namiki, Tsukuba, Ibaraki 305-8564, Japan
We examine the applicability of a CFD simulation method to microscopic two-phase fluid motions on partially-wetted
solid surfaces with edges for evaluating fluidic devices and for predicting underground fluid flows. The method adopts
phase-field model (PFM) and lattice-Boltzmann model (LBM). PFM, based on the free-energy theory, does not
necessarily require conventional algorithms for reconstruction of interfaces. One of features of a semi-Lagrangian
LBM is simple mesoscopic particle-kinematic operation in discrete conservation form on an isotropic spatial lattice.
The major findings in immiscible liquid-liquid simulations are as follows: (1) the method simulates well departure of
droplet from flat solid surface due to buoyancy in agreement with semi-empirical prediction; (2) the droplet has larger
static contact angle on a textured hydrophobic surface than on the flat one in a similar way with experimental data; (3)
the fluid penetration between parallel rectangular bodies is affected by their front, inside and outside face wettabilities.
1.緒 言
本研究では,マイクロ流体デバイス(1),微細製造プロセス(2)-(5),
地層など自然環境等で見られる,固体表面の不均一な濡れ性や複
雑な凹凸形状を持つ人工物体上や多孔質構造を有する天然物質内
部における微小スケール二相流体挙動の詳細な予測を目的として,
格子ボルツマン法(Lattice Boltzmann Method, LBM)(6)-(13)とフェー
ズフィールドモデル(Phase-field Model, PFM)(14)-(19)に基づく数値
流体力学
(CFD)
シミュレーション法(4),(20),(21)の適用性を検討する.
LBM は,Macro-scale の連続体を構成する Meso-scale の仮想粒子
集合の単純な衝突・並進運動の反復に基づくボトムアップ的アプロ
ーチで連続体現象を創発的に再現する.この特徴は,利点として,
(1) 流れ場内の複雑形状物体境界の容易な再現,(2) 計算コードの
容易な作成,(3) 質量・運動量等保存性と空間等方性の両方に優れ
る移流スキーム,(4) 並列化効率の高い計算,等をもたらす(6)-(13).
PFM-CFD 法では,非平衡系の自由エネルギー理論に基づく材
料科学分野の Phase-field 法(14),(17)-(19)に類似する方法で,拡散流束を
瞬時局所的に調整して有限で一定の幅を持つ流体界面を自律形成
する(7),(15),(16).この特徴から,本 CFD 法は多数の界面の同時多発的
大規模変形・移流を従来法よりも効率的に計算できる利点を有す
る(16)-(24).また,濡れ性が不均一かつ凹凸形状の固体表面に接触す
る流体界面(接触線)の挙動も,有限幅の界面領域内で連続的に
境界条件を課すことによって従来よりも容易に再現できる(8),(9),(11).
尚,本研究では,計算効率と数値安定性の観点から流体界面を実
在よりも非常に分厚い人工的な連続体領域として扱う(20),(21).
上述の利点を備える LBM と PFM に基づく混相流 CFD 法に関
しては,国内外で既に多くの研究開発成果と適用事例が報告され
ている(7)-(11),(13).我々は近年,文献で示される新しい PFM を導入
する界面移流方程式についてLBM で解く二相流CFD 法を提案し
(4),(20),(21),T 型マイクロ流路内液液二相系スラグ液滴流動シミュレ
ーションで実験結果と良く一致する数値結果を得た(3),(4),(21).以上
を踏まえ,本報では,固体表面の形状と濡れ性がより複雑な微小
二相流体問題を通して当該 PFM-CFD 法の適用性を更に検討する.
2.拡散界面モデル基礎方程式
本研究では,連続体近似が成り立つ非圧縮性・非混和性・等温
二相流体系を対象とし,座標 x の空間で各時刻 t の系の状態を得
るため,PFM(13)-(15)に基づく以下の連続の式(1),運動方程式(2),
及び修正保存形 Allen-Cahn(AC)型拡散界面モデル移流方程式(1
ステップ保存形レベルセット移流方程式)(3) (28)-(31)を解く(4),(20), (21).
(1)
∇ ⋅u = 0
F (φ )
∂u
∇p
+ u ⋅∇ u = −
+ ν ∇ 2u + I
∂t
ρ
ρ
(2)
∂φ
+ u ⋅∇φ = ∇ ⋅ ( D(φ )∇φ )
∂t
(3)
上式で,u は流速,圧力 p,φ は界面形状を表す連続変数,ρ は密
度,ν は動粘性係数である.本研究では,式(2)右辺の界面張力γ に
よる力 FI (φ)及びφ の拡散係数 D(φ)を次式で与えた(4),(15),(28),(20),(21).
FI (φ ) = −6γε (∇ 2φ ) ∇φ
(4)

 φ (1 − φ ) 

D(φ ) = D0 1 −

ε
∇
φ




(5)
ここでε は界面幅δφに関するパラメータ,D0 は定数である.式(3)
及び(5)により,
指標φ (x, t) = 0 及び 1 の領域に相当する二相の間で,
界面は0 < φ (x, t) < 1 で代表幅δφ = 4ε の連続体とモデル化される(28).
変数φ は平衡状態で界面の法線軸ξ方向に次の分布を取る(15),(28),(29).

 4ξ

φ (ξ ) = 1 + exp  −
 δ

 φ






−1
 1

 2ξ
 = 1 + tanh 
δ
 2
 φ




  



(6)
ここでφ (-∞) = 0,φ (0) = 1/2,φ (+∞) = 1 である.界面代表幅δφ は,
流れ場に現れるφの空間勾配の絶対値|∇φ|が式(6)に基づく理論値
φ (1-φ)/ε よりも大きい場合は式(3)右辺の正の拡散により増加し,
逆の場合は負の拡散により減少する.その結果,右辺が消える平
1
Copyright © 2014 by JSFM
第 28 回数値流体力学シンポジウム
D05-1
衡状態で|∇φ | = φ (1-φ)/εを満たすδφ =一定の界面が形成される.本
法では,φ の 4 階微分の拡散項を持つ Cahn-Hilliard(CH)型移流
方程式を用いる PFM-CFD 法(7)-(9),(11),(22),(23)同様,δφ を空間分割セル
幅∆x 数個相当に設定した
(後述の計算ではε = ∆x によりδφ = 4∆x)
.
修正保存形 AC 型の拡散界面モデル移流方程式(3)は,CH 型方
程式よりも低い階数の微分を行い且つ界面曲率依存性を排した拡
散項を持つ(28),(29),(32).
この特徴により,
式(3)を採用する本PFM-CFD
法では,CH 型方程式を採用する PFM 法(7)-(9),(11), (22)-(26)よりも,界面
形成・移流演算負荷が低減できるとともに,体積保存性が向上し
つつ界面張力効果が移流方程式から取り除かれるため,二相流体
現象の予測の高精度化が期待できる(4),(20),(21).
{e0 , e1 , e2 , e3 , e4 , e5 , e6 , e7 , e8 , e9 , e10 , e11 , e12 , e13 , e14 }
 0 1 −1 0 0 0 0

= c  0 0 0 1 −1 0 0
 0 0 0 0 0 1 −1

1 −1 1 −1 1 −1
−1

1 −1 1 −1 −1 1 −1 1 
1 −1 −1 1 1 −1 −1 1 
(
∂ ga
1
+ ea ⋅∇g a = −
ga − g aeq
∂t
τg
(
)
)
f
ρ u = ∑ a f a ea
(10)
φ = ∑ a ga
(11)
∑
a
(19)
D0 (1 − φ ) ∇φ
ε
∇φ
j=
∑
(13)
∑
∆t
 g a (x, t ) − g aeq (x, t )  (14)
τg 
wae aea = cS2 I
(16)
(21)
(22)
平衡分布 faeq 及び gaeq の粒子速度 ea に関するモーメントの和は,
対象とするマクロな連続性流体の変数と以下の関係を満たす.
f aeq = ρ
(23)
f aeq ea = ρ u
(24)
f aeq eaea = ρuu + pI
(25)
∑
(15)
(20)
式(20)の j は,式(3)右辺が含むφ に関する次の逆拡散流束を表す(28).
(12)
wa = 1
a
2


 ea ⋅ u 


 − u ⋅ u 
c


 S 

( for a ≠ 0)
cS = c / 3
a
a
∑
∑
ここで,各分布関数と流体に関する変数は全て空間セル中心で定
義される.∆t は一定の時間進行幅,t' は並進演算後の中間時刻,
cS は音速であり,重み係数 wa は a に関する総和の次式を満たす.
∑
(18)
ここで,式(19)は a = 0 の静止粒子(e 0 = 0)に対する関数,重み係
数 w0 = 2/9,wa = 1/9(a = 1~6)
,wa = 1/72(a = 7~14)
,Γは静止
平衡流体中の静止・運動各粒子数密度を調整するパラメータ(13)で
あり,式(17)の ea に対する音速 cS は次式で表される.
時刻 t で空間分割セル中心点 x 上に分布する a 方向に向かう粒
子の fa (x,t)と ga (x,t)を計算するため,式(7)及び(8)を次のような空
間・時間 2 次精度のセミ・ラグランジュ形式に離散化した(10),(20),(21).
g a (x + ea ∆t , t + ∆t ) = g a (x, t ) −
2


 ea ⋅ u 


 − u ⋅ u 
c


 S 


e ⋅ (u + j) 1
+ 2
g aeq = waφ Γ + a 2
cS
2cS


(8)
(9)
w e ⋅ F (x, t ')
f a (x, t + ∆t ) = f a (x, t ') + ∆t a a 2I
cS
 e ⋅u 1
= wa ρ 1 + a 2 + 2
cS
2cS


(7)
ρ = ∑ a fa
∆t
 f a (x, t ) − f aeq (x, t ) 
τf 
eq
a

u ⋅u 
g0eq = φ 1 − Γ(1 − w0 ) − w0 2 
2cS 

下付添字 a は粒子並進速度方向指標,τf 及びτg は fa,ga 各々が粒子
間相互作用(衝突)により局所平衡状態 faeq,gaeq に至るまでの緩
和時間を表す.
流体の変数は分布関数を使って次式で定義される.
f a (x + ea ∆t , t ') = f a (x, t ) −
1
c = ∆x/∆t は格子定数である.
2 次元では9 速度モデルを採用した(20).
eq
eq
平衡分布関数fa 及びga には,
Maxwell-Boltzmann 分布のTaylor
展開に基づき u の 2 次までの項を含む次の形式を与えた(20),(21).
3.計算スキーム
本 CFD 法では,式(1)~(3)の数値解法として,連続体を構成す
る微視的仮想粒子集合の挙動を統計力学的に記述する LBM(6)-(13)
を採用し,等方的並進速度 ea 毎の粒子数密度を表す離散速度分布
関数 fa (x, t)及び ga (x, t)に関する次の時間発展を計算する(20),(21).
∂ fa
F ∂f
1
+ ea ⋅∇f a + I ⋅ a = −
f a − f aeq
∂t
ρ ∂ ea
τf
(17)
a
a
a
g aeq = φ
g aeq ea = φ ( u + j)
(26)
(27)
流体の質量と運動量の保存に関する方程式(1)及び(2)は fa の時間
発展式(12)及び(13)から,拡散界面移流方程式(3)は ga の時間発展
式(14)から,漸近理論展開(7)等により局所平衡及び低 Mach 数の極
限条件下で導かれる.
その中で,
p, ν及びD0 は次式で表される(20),(21).
I は Kronecker Delta δαβで成分が表される 2 階の等方テンソルであ
る.式(12)及び(14)は,運動粒子が∆t の並進後に ea の方向に隣接
する空間セルの中心位置に到着することを意味する.
本研究では,時刻 t を一定幅∆t = 1 で進行させ,3 次元(x, y, z)
デカルト座標系物理空間を幅∆x = ∆y = ∆z = 1 の立方形セルで等分
割するとともに,3 次元(ea,x, ea,y, ea,z)座標系の速度空間で次の 15
種類のベクトル e(a
= 0∼14)から成る粒子速度集合を採用した(7).
a
p = ρ cS2
(28)
∆t 

ν = cS2 τ f − 
2 

(29)
∆t 

D0 = ΓcS2 τ g − 
2 

(30)
本研究では,流れ場内で界面形状を厚さ一定で形成し且つ各相体
2
Copyright © 2014 by JSFM
第 28 回数値流体力学シンポジウム
D05-1
積を保存するため,φ の数値拡散誤差が最小になるよう単一緩和
時間τg を決定するとともに,D0 に基づく代表的な拡散速度が移流
の代表速度U と同程度になるよう(28) Γの値を次のように調整した
(20),(21)
.まず,以下の関係を想定した(28).
D0
≅U
ε
g
O
(31)
(A) VD = 1.767ml (departure)
(B) VD = 0.382ml (non-departure)
(a) Static contact angle θW = 90°
次に,文献(6)を参考に衝突緩和係数τg を次のように選択した.
τg =
∆t 
1 
1+

2 
3
zx
(32)
式(32)を式(30)に代入すると,拡散の定数 D0 が次のように決まる.
D0 = ΓcS2
∆t
2 3
(33)
(C) V D = 0.524ml
(departure)
上式と式(31)から式(18)及び(19)のΓを以下のように与えた.
Γ≅2 3
Uε
cS2 ∆t
(34)
Droplet volume VD = π dD3/6 [ul]
(35)
(36)
nS は固体表面で流体領域側に向く単位法線ベクトルである.また,
他の PFM 法(8),(9),(11)同様,自由エネルギー理論に基づき二相流体に
対する固体表面の濡れ性を考慮する次式によって静的接触角θW
を設定した(20),(21).
u=0
n S ⋅ ε 2∇φ = −γ S
(37)
γ S ∆x
ε2
104
103
102
101
100
10-1
10-2
10-3
10-4
(A)
(C)
(B)
(D)
(E)
γ = 51[ mN/m ]
γ
Fritz’s eq. :
d D = 0.0209
θW
g ∆ρ
≤ 140° )
Numerical results : non-departure
: departure
( 0 ≤ θW
0
20 40 60 80 100 120 140
Continuous-phase static contact angle θW [deg.]
Wetting potentialと呼ばれるパラメータγS に空間セル表面毎で異な
る値を設定することにより不均一な濡れ性を表現できる(8),(9),(11).
γS = 0はθW = 90°を与える.濡れ性条件式(37)は,式(4)及び(22)によ
る界面力FI及び拡散流束jの計算で考慮される.本研究では,固体
表面に接する流体領域内の空間セルでFIとjを計算する際,固体領
域内の仮想空間セルの変数φを2次精度差分近似の外挿により求め
た.例えば,nSがx軸方向の固体表面上では次式が使用される(23).
φ −1 = φ 0 +
(E) VD = 0.065ml
(non-departure)
Fig. 1 3D Simulation of oil droplet with volume VD on horizontal flat
solid surface with θW in stagnant water (continuous phase) under gravity g.
流れ場の滑り無し静止固体境界上では,入射速度eaの粒子の分布
関数faを逆方向速度eb = - ea の関数fb に変更するBounce-back条件を
適用した.その結果,圧力p及び流速uは以下の境界条件を満たす.
n S ⋅ ∇p = 0
(D) VD = 0.113ml
(non-departure)
(b)θW = 60°
Fig. 2 Diagram of departure of droplet from horizontal flat solid surface
under gravity g in terms of static contact angle θW and volume VD.
4.2 物体角部での液滴の滑落
角部を持つ物体に付着する液滴の挙動に対する濡れ性の影響を
検討するため,θW = 90°の静止壁で囲まれる一辺 L の 2 次元正方
領域(Fig. 3)内に,θW1 及びθW2 の上面と側面を持つ物体とその上
に半径 R ≅ L/2 の四半円状液滴を初期配置した.
液滴を物体から滑
落させるのに要する液滴に課す下向き外力 G の大きさを計算に
より見積もった.空気-水系で R = 65µm の液滴を想定した.ただ
し,微小空間で短時間には気液間の物性差の影響が非常に小さい
と仮定して密度と粘性は両相で等しくした.変数φ の空間分布の
シミュレーション結果(Fig. 4,φ = 0 及び 1 の各領域を白色と青
色で示す)では,低解像度(図(a) L = 64∆x)
・高解像度(図(b) L =
128∆x)条件ともに,θW2 の各条件で無次元数 Eo = |G|R2/γ がある
値を超えると液滴は角部を超えて滑落するが,物体側面が疎水性
を強めるに従い液滴は上部に留まり続ける傾向が確認された.
(38)
ここで下付添字 -1 及び 0 は各々,固体表面に接する固体側の仮想
空間セルでの値及び流体側に実在する空間セルでの値を意味する.
4.数値シミュレーション結果
4.1 重力下での水平平坦固体表面上の液滴の離脱
本 CFD 法における界面張力と固体表面濡れ性の効果を検討す
るため,固体表面に付着する単一液滴の挙動を計算した.今回は
重力下にある界面張力γ = 51mN/m の油(分散相)-水(連続相)
系を想定し,初期に半球形状で球等価直径 5~15mm の油滴が重
力下で付着する水平な平坦面上で水に対する静的接触角θW を 90°
または 60°に設定した.毛細管力と浮力の効果が支配的な油滴の
初動に注目することから,二相の密度差は浮力計算でのみ考慮す
るとともに粘性µは両相で同一とした.Fig. 1 から,油滴(体積 VD)
が大きい程,及びθW が小さい(親水性が強い)程,浮力効果及び
撥油効果の増加により油滴はより離脱し易くなることが分かる.
これら数値結果は Fritz の式による離脱臨界直径予測(Fig. 2 中の
実線で示す)(33)に合致することから,本 CFD 法で界面張力と濡
れ性の影響が適切に考慮され得ることが確認できる.
= 1.456 ×10
Eo = G R 2γ −1
R≅ L/2
−2
R
L = 64∆x,
128∆x
y
O
L
Solid walls
φ
Oh = µ / ργ R
θW 1
x
Droplet
L
Solid body
G
L/2
θW 2
L/4
Fig. 3 2D computational domain and initial condition of droplet on
rectangular solid body with edge under external downward force G.
3
Copyright © 2014 by JSFM
第 28 回数値流体力学シンポジウム
D05-1
Eo = 4.75
Eo = 3.5
t =0
Eo = 5.5
G
φ
60°
60°
L/4
60°
Eo = 3.75
60°
90°
Eo = 5.0
113°
z
O x
t =0
Eo = 5.75
G
y
60°
x
60°
60°
60°
90°
113°
O
(a) Lower resolution L = 64∆x
Eo = 3.4
Eo = 4.6
(a) θW = 130°
(b) θW = 113°
(c) θW = 101°
Fig. 6 Steady-state shape of droplet on grooved solid surface under
downward external force G for static contact angle θW and dimensionless
numbers of Oh = µ/(ργ dD)1/2 = 1.14×10-2 and Eo = |G|dD2/γ = 1.76×10-2.
Eo = 5.25
φ
60°
L/4
60°
60°
60°
90°
113°
Eo = 4.7
Eo = 3.5
zy
t =0
Eo = 5.5
Gx
y
60°
x
60°
O
60°
60°
90°
113°
t =0
(b) Higher resolution L = 128∆x
Fig. 4 2D droplet on solid body for dimensionless number Eo = |G|R2/γ .
Gy= 0
4.3 テクスチャ表面上の液滴の接触角と移動性
次に,
凹凸固体表面に付着する液滴の CFD シミュレーションを
実施した.3 次元デカルト座標系(x, y, z)の計算領域を幅∆x=1 の
単位立方セル 92×92×96 個で一様分割し,x 軸方向両端・z 軸方向
上端各境界に滑り無し固体壁,y 軸方向両端に周期境界,厚さ wS =
6∆x,高さ h = 30∆x,長さ 90∆x の固体平板を wF = wS の等間隔で領
域下部に設置した.初期条件では,その固体平板上端に接するよ
うに直径 dD = 8wS の球状液滴を置き,その周囲を連続相で満たし,
静止状態の液滴に外力 G を加えた.密度ρ = 998.2kg/m3, 粘性µ =
1mPa⋅s, 界面張力γ = 32mN/m,dD = 240µm の液液系(無次元数 Oh
= µ/(ρ γ dD)1/2 = 1.14×10-2 及び Eo = |G|dD2/γ = 1.76×10-2)を想定し,
前小節同様ρとµ は両相で同一とした.本シミュレーションでは,
平面上の液滴の静的接触角θW = 130°,113°及び 101°の条件に対し
て,固体面が平坦な x-z 面内よりも凹凸状の y-z 面内で,より大き
な接触角が観察された(Fig. 6, φ = 1/2 の等値面を界面と描く)
.ま
た,液滴は,固体表面が一様に平坦な x 軸方向に G を課す場合は
動く(Fig. 7 (a))一方,凹凸な y 軸方向への G の印加では初期位
置に留まり続けた(Fig. 7 (b))
.以上のシミュレーション結果は文
献(34)や実験で示される結果と定性的に一致することを確認した.
Phase B
φ =1
φ =0
dD
Lz
Solid
Interface:
t =0
Gx= 0
O
h
Solid
Gy
O
O
zy
(b) Force components Gx = 0, Gy ≠ 0 and Gz = 0
Fig. 7 Snapshot of droplet on grooved solid surface under external force
G = (Gx, Gy, 0) for θW =113°, Oh = 1.14×10-2 and Eo = 1.76×10-2.
4.4 平行物体周りの二相流体挙動
本節では,2 次元平行物体周囲の二相流体挙動の CFD シミュレ
ーションについて述べる.Fig. 8 に示すように,一辺長さ 6L の正
方形状の計算領域で,上下に相 A に対する静的接触角θW の滑り
無し静止水平固体壁境界,左側に時間一定の代表速さ|Uin| = 10-3c
の一様流入境界,及び右側に連続流出境界を設置した.この領域
中央に,長さ 3L・厚さ L/2 で平行平板状の 2 物体が代表長さ L =
40∆x 離れて配置されている.各物体の左端面,内向き側面,外向
き側面,及び右端面では各々,滑り無し条件と相 A の接触角θW,1,
θW, in, θW, out 及びθW を与えた.初期条件では,物体周囲を相 B が占
有し,時間の経過とともに相 A が流入する.相 A 及び B が同じ
密度ρA = ρB 及び粘性µA = µB を有する二相系を想定し,Reynolds 数
Ly
Lx
zx
t =0
wS wF
φ = 0.5
z x
O
dD
zy
(a) Force components Gx ≠0 and Gy = Gz = 0
Periodic boundary
Phase A
O
Solid
Solid
Solid
Periodic boundary
Solid
zx
z y
Fig. 5 Computational domain and initial condition for CFD simulation of
droplet on textured solid surface under external force in three dimensions.
4
Copyright © 2014 by JSFM
第 28 回数値流体力学シンポジウム
D05-1
5.結 言
本報では,マイクロ流体デバイス(1),微細加工・製造プロセス(2)-(5),
自然環境等で見られる,角部や凹凸表面を持つ人工物体上や多孔
質構造を有する天然物質内部における微小二相流体挙動の予測を
目的として,フェーズフィールドモデル(PFM)(14)-(19)と格子ボル
ツマン法(LBM)(6)-(13)に基づく数値流体力学(CFD)シミュレー
ション手法(4),(20),(21)を概説した.続いて,重力下・水平平坦固体表
面からの液滴離脱,固体角部での 2 次元液滴滑落,テクスチャ表
面上の 3 次元液滴挙動,
及び平行物体周り二相流体挙動の CFD シ
ミュレーションを通して,固体表面濡れ性問題への PFM-LBM ベ
ースの本 CFD 手法の適用可能性を確認した.同時に,CFD 結果
から,表面の濡れ性及び凹凸形状が微小スケールの二相流体の挙
動に影響を与え得ること,不均一な表面形状及び濡れ性を利用し
て微小二相流体の形状と挙動を制御可能であることも示した.
Re = ρA L|Uin|/µA = 1.0,Capillary 数 Ca = µA |Uin|/γ = 2.86×10-5(Weber
数 We = Re×Ca = 2.86×10-5 及び Ohnesorge 数 Oh = (Ca/Re)0.5 =
5.35×10-3)を設定した.
Fig. 9 に,θW = 90°,θW, in = θW, out = 78°で,左端面でのみ相 A に
対する濡れ性をθW,1 = 30°, 90°, 180°として変化させた場合の時刻 t
= 3×104∆t における二相の空間分布の結果を示す.これらの図から,
左端面が撥液性を増すに従って物体間への相Aの浸透は困難にな
ることが確認される.また,界面の x 軸方向初期位置 s0 = 70.5∆x
からt = 2.5×104∆t 後の移動距離 δs(t) = s(t) - s0 に関するFig. 10 では,
左端面接触角θW,1 の各値に対して,内側面を親液性にすることに
加えて,外側面を撥液性に設定することによっても,相 A が物体
間に浸透し易くなることが分かる.プラスチック製ピペットによ
る水滴吸引実験でも,上記シミュレーションと同様に,空気-水界
面に接するピペット先端面の親水化処理による浸透改善傾向が確
認されている(1).
今後,
本CFD結果をより詳細に検討予定である.
Phase A
ρA
µA
Phase B
ρB
2L
µB
φ=1 φ=0
θW,1
θW,out
θW,1
L
2L
y
θW
6L (= 240∆x)
O
6L
L/2
Solid plate
θW,out
x
Stationary
& non-slip
solid wall
L/2
Solid plate
θW,in
θW,in 3L
Uin
謝 辞
上記成果は,総合科学技術会議が制度設計した最先端研究開発
支援プログラム(FIRST)課題「マイクロシステム融合研究開発」
(江刺プロジェクト,2009 年度~2013 年度)への日本学術振興会
(JSPS)による助成,及び JSPS 科研費基盤研究 (C) 課題 No.
26420128「フェーズフィールドモデルに基づくマイクロ多孔質体
内相変化二相流計算法の開発」での助成で得られた.関係者各位
に感謝申し上げます.
θW
参考文献
(1) Kurihara, K., Takada, N., Takagi, H. and Maeda, R., "Capillary
effect enhancement in plastic-pipette by nano-textured surface," The
40th International Conference on Micro and Nano Engineering
(MNE2014), Lausanne, Switzerland, September 22−26, 2014,
Abstract ID: 8263.
(2) Matsumoto, S., Takada, N. and Matsumoto, J., “Microfabrication
process of cellular structures in hollow fiber-shaped substrates,”
Transactions of the Japan Society of Mechanical Engineers Series C,
Vol. 76 (2010), pp. 1911–1913
(http://ci.nii.ac.jp/naid/110007685434/).
(3) Miki, Y., Matsumoto, S., Takada, N., Kaneko, A. and Abe, Y.,
"Effects of flow rate on formation behavior of two-phase slug flow
in a microchannel T-junction," Japanese Journal of Multiphase Flow,
Vol. 27 (2014), pp. 591−598 (DOI: 10.3811/jjmf.27.591).
(4) Takada, N., Matsumoto, J. and Matsumoto, S., "Prediction of
two-phase slug flow patterns in microchannel with T-junction,"
Japanese Journal of Multiphase Flow, Vol. 28 (2014), pp. 32−38
(DOI: 10.3811/jjmf.28.32).
Stationary
solid wall
L= 40∆x
Fig. 8 2D computational domain and initial condition.
φ
y
O
x
(a) θW,1 = 30°
(b) θW,1 = 90°
(c) θW,1 = 180°
Interfacial displacement δs/L [-]
Fig. 9 Snapshot of spatial distribution of two-phase fluid penetrating
between parallel plates with static contact angle θW,1 on left-side edge
surfaces at time t = 3×104∆t, for θW = 90° and θW, in = θW, out = 78°.
3.0
Time t = 2.5×104∆t
2.5
2.0
1.5
1.0
0.5
(5) Lu, J., Takagi, H., Nakano, T. and Maeda, R., “Size-free MEMS-IC
high-efficient integration by using carrier wafer with self-assembled
monolayer fine pattern,” Proceedings of 2013 IEEE 63rd Electronic
Components and Technology Conference (ECTC2013) (2013), pp.
1508−1513.
(6) Hirabayashi, M., Chen, Y. and Ohashi, H., “The lattice BGK model
for the Poisson equation,” JSME Int. J. Ser. B Fluids Therm. Eng.,
Vol. 44 (2001), pp. 45–52 (DOI: 10.1299/jsmeb.44.45).
(7) 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.2004.01.019).
(8) Kobayashi, K., Inamuro, T. and Ogino, F., "Numerical simulation of
advancing interface in a micro heterogeneous channel by the lattice
θW,in= 78°,θW,out= 101°
θW,in=
θW,out= 78°
θW,in= 78°,θW,out= 90°
θW,in=
θW,out= 90°
0
30
60
90
120
150
180
Contact angle of left edge surface θW,1 [deg.]
Fig. 10 Dimensionless displacement δs/L of two-phase interface from
initial position on center line of computational domain at t = 2.5×104∆t.
5
Copyright © 2014 by JSFM
(9)
(10)
(11)
(12)
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
(22)
Boltzmann method," J. Chem. Eng. Jpn., Vol. 39 (2006), pp.
257–266 (DOI: 10.1252/jcej.39.257).
Huang, J. J., Shu, C. and Chew, Y. T., “Lattice Boltzmann study of
droplet motion inside a grooved channel,” Phys. Fluids, Vol. 21
(2009), 022103, 11 pages (DOI: 10.1063/1.3077800).
Suzuki, K. and Inamuro, T., “Effect of internal mass in the
simulation of a moving body by the immersed boundary method,”
Comput. Fluids, Vol. 49 (2011), pp. 173–187
(DOI: 10.1016/j.compfluid.2011.05.011).
Yoshino, M., Kobayashi, Y. and Tanaka, Y., “Numerical simulation
of liquid penetration through spherical bodies with various types of
wettability,” Japanese Journal of Multiphase Flow, Vol. 26 (2013),
pp. 499–506 (DOI: 10.3811/jjmf.26.499).
Seta T., Rojas, R., Hayashi, K. and Tomiyama, A., “Implicitcorrection-based immersed boundary-lattice Boltzmann method
with two relaxation times,” Phys. Rev. E, Vol. 89 (2014), 023307,
22 pages (DOI: 10.1103/PhysRevE.89.023307).
Yoshida, H., Kobayashi, T., Hayashi, H., Kinjo, T., Washizu, H. and
Fukuzawa, K., "Boundary condition at a two-phase interface in the
lattice Boltzmann method for the convection-diffusion equation,"
Phys. Rev. E, Vol. 90 (2014), 013303, 10 pages
(DOI: 10.1103/PhysRevE.90.013303).
Kobayashi, R., “Modeling and numerical simulations of dendritic
crystal growth,” Physica D, Vol. 63 (1993), pp. 410–423
(DOI: 10.1016/0167-2789(93)90120-P).
Anderson, D. M., McFadden, G. B. and Wheeler, A. A., “Diffuseinterface methods in fluid mechanics,” Annu. Rev. Fluid Mech., Vol.
30 (1998), pp. 139–165 (DOI: 10.1146/annurev.fluid.30. 1.139).
Onuki, A., “Dynamic van der Waals theory,” Phys. Rev. E, Vol. 75
(2007), 036304, 15pages (DOI: 10.1103/PhysRevE.75.036304).
小山 敏幸, 塚田 祐貴, 材料組織弾性学と組織形成 – フェ
ーズフィールド微視的組織弾性論の基礎と応用, 内田老鶴圃,
2012.
高木 知弘, 山中 晃徳, フェーズフィールド法 – 数値シミ
ュレーションによる材料組織設計, 養賢堂, 2012.
小山 敏幸, 高木 知弘, フェーズフィールド法入門 (計算力
学レクチャーコース, 日本計算工学会 編), 丸善出版, 2013.
Takada, N., Matsumoto, J. and Matsumoto, S., “Phase-field
model-based simulation of motions of a two-phase fluid on solid
surface,” Journal of Computational Science and Technology, Vol. 7
(2013), pp. 322–337 (DOI: 10.1299/jcst.7.322).
Takada, N., Matsumoto, J. and Matsumoto, S., “A diffuse-interface
tracking method for the numerical simulation of motions of a
two-phase fluid on a solid surface,” The Journal of Computational
Multiphase Flows, Vol. 6 (2014), pp. 283–298
(DOI: 10.1260/1757-482X.6.3.283).
Takada, N. and Tomiyama, A., “A numerical method for two-phase
flow based on a phase-field model,” JSME Int. J. Ser. B Fluids
(23)
(24)
(25)
(26)
(27)
(28)
(29)
(30)
(31)
(32)
第 28 回数値流体力学シンポジウム
D05-1
Therm. Eng., Vol. 49 (2006), pp. 636–644
(DOI: 10.1299/jsmeb.49. 636).
Takada, N., Matsumoto, J., Matsumoto, S. and Ichikawa, N.,
“Application of a phase-field method to the numerical analysis 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).
Matsumoto, J. and Takada, N., “Two-phase flow analysis based on a
phase-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).
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 the National
Congress of Theoretical and Applied Mechanics (NCTAM), Japan,
Vol. 60 (2011), Session ID: GS03-04
(https://www.jstage.jst.go.jp/article/japannctam/60/0/60_0_251/_arti
cle/).
Sakakibara, T., Takaki, T. and Kurata, M., Proceedings of the 5th
Asia Pacific Congress on Computational Mechanics (APCOM) &
the 4th International Symposium on Computational Mechanics
(ISCM) (2013), Paper ID 1381
(http://www.sci-en-tech.com/apcom2013/APCOM2013-Proceedin
gs/PDF_FullPaper/1381_Tetsuya_Sakakibara.pdf).
Fukuta, M. and Yamamoto, Y., “Development of numerical analysis
model for evaluating of boiling heat transfer,” Japanese Journal of
Multiphase Flow, Vol. 28 (2014), pp. 161–166
(DOI: 10.3811/jjmf.28.161).
Chiu, P.-H. and Lin, Y.-T., “A conservative phase field method for
solving incompressible two-phase flows,” J. Comput. Phys., Vol.
230 (2011), pp. 185–204 (DOI: 10.1016/j.jcp.2010.09.021).
Olsson, E. and Kreiss, G., “A conservative level set method for two
phase flow,” J. Comput. Phys., Vol. 210 (2005), pp. 225–246
(DOI: 10.1016/j.jcp.2005.04.007).
Tan, N., Aoki, T., Inoue, K. and Yoshitani, K., “Numerical
simulation of two-phase flow driven by rotating object,”
Transactions of the Japan Society of Mechanical Engineers, Series
B, Vol. 77 (2011), pp. 1699–1714 (DOI: 10.1299/kikaib.77.1699).
Sato, Y. and Niceno, B., “A new contact line treatment for a
conservative level set method,” J. Comput. Phys., Vol. 231 (2012),
pp. 3887–3895 (DOI: 10.1016/j.jcp.2012.01.034).
Beaucourt, J., Biben, T., Leyat, A. and Verdier, C., “Modeling
breakup and relaxation of Newtonian droplets using the advected
phase-field approach,” Phys. Rev. E, Vol. 75 (2007), 021405, 8
pages (DOI: 10.1103/PhysRevE.75.021405).
(33) 西川 兼康, 藤田 恭伸, 伝熱学, 理工学社, pp. 219–223, 1982.
(34) 中島 章, 固体表面の濡れ制御, 内田老鶴圃, 2007.
6
Copyright © 2014 by JSFM