第 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: [email protected] 松本純一, 産総研, 〒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)，微細製造プロセス(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 法の適用性を更に検討する．２．拡散界面モデル基礎方程式本研究では，連続体近似が成り立つ非圧縮性・非混和性・等温二相流体系を対象とし，座標 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)．３．計算スキーム本 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. ４．２物体角部での液滴の滑落角部を持つ物体に付着する液滴の挙動に対する濡れ性の影響を検討するため，θ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 は各々，固体表面に接する固体側の仮想空間セルでの値及び流体側に実在する空間セルでの値を意味する．４．数値シミュレーション結果４．１重力下での水平平坦固体表面上の液滴の離脱本 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 ４．３テクスチャ表面上の液滴の接触角と移動性次に，凹凸固体表面に付着する液滴の 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. ４．４平行物体周りの二相流体挙動本節では，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 ５．結言本報では，マイクロ流体デバイス(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