PDFファイルはこちらから - 熱エネルギー工学研究室へようこそ

修士論文
液滴の動的濡れ現象の分子動力学
平成 24 年 2 月 10 日提出
指導教員 塩見淳一郎 准教授
106219 中村芳紀
2
目次
第1章
序論
4
1.1
はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.2
先行研究 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.2.1
平衡濡れ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.2.2
動的濡れ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.2.3
分子動力学法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
分子動力学法
15
2.1
支配方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
2.2
温度制御法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
2.3
並列化手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
2.4
計算条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
2.5
Lennard-Jones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
平滑面上の動的濡れ性
23
3.1
計算系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
3.2
平衡接触角 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
3.3
動的濡れ性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
1.3
第2章
第3章
3.3.1
サイズ依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
3.3.2
平衡接触角依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
3.4
流速場 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
3.5
まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
Cassie 条件における動的濡れ現象
36
4.1
概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
4.2
平衡接触角 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
4.3
動的濡れ性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
第4章
4.4
第5章
4.3.1
定性的理解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
4.3.2
stick-slip . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
マルチスケールシミュレーション
44
3
目次
5.1
概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
5.2
支配方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
5.3
MD と Phase-Field の比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
46
結論
49
第6章
参考文献
51
4
第1章
序論
1.1
はじめに
人間のみならず全ての生物の生命活動にとって水は不可欠である．そのためか，濡れ性（親水性・疎水
（撥水）性）のコントロールによって水と上手く付き合い，優れた機能を発現している生物が数多い．例
えば，アメンボが分泌する撥水性の油分は水面での生活を可能にし，ハスの葉の表面の細かな凸凹による
高い撥水性は汚れの付着を防ぎ，砂漠の昆虫体表のミクロな疎水性・親水性のパターンや蜘蛛の糸のミク
ロ構造は水の収集能力を高めている [1–3]．
このように，濡れ性の利用は自然界に多く現れる一方で，工学応用としては撥水性を施すことで抵抗を
抑えた高速水着や超親水性を利用した視認性の良いガラスなどのほか，ごく最近ではウツボカズラに着想
を得た自己修復性を持つ撥水面の開発により生物医学デバイス等への応用可能性が示されるなど今後の発
展も期待される [4]．
このような濡れ性に関する物理の理解は 200 年前の Young，Laplace らの時代に遡り，親水性，疎水性
の差異が固液気の界面エネルギー平衡論から十分に理解されている．
先に静的濡れ特性の工学応用例を述べたが，動的濡れ現象が根幹の物理現象となる例も多い．例えば潤
滑油被膜の場合には，潤滑油が固体表面を被膜するのにかかる時間とその厚みが問題となり，動的濡れ現
象の理解が必要となる．他にも，半田付け，インクジェット，スプレープリント，焼結現象も動的濡れと
関連が深い．しかし，移動境界線問題と呼ばれる固液気の三相界面の取り扱いが困難なために，動的濡れ
現象の完全な理解には未だ至っていない．
動的濡れ現象の物理的理解のアプローチの一つに自然濡れ現象がある．自然濡れとは，球形の液滴が速
度ゼロで基板に接触し平衡状態に到るまでの濡れ現象のことを言う．詳細は後述するが，自然濡れに関す
る課題の一つに低粘性液体を用いた時の高速濡れと呼ばれる現象がある．高速濡れに現象における物理的
理解も未だ不明確であり，界面における物理現象の理解の重要性が示されているに過ぎない [5, 6]．
分子動力学法は分子運動を計算し統計的に解析する実験的シミュレーションであり，連続体近似を介さ
ない原理的な手法であるため，界面現象を正確に追う事が可能である．本研究では分子動力学法を用いて
高速濡れ現象を原理的に解くことで現象の理解を進める．高速濡れ現象の物理の理解は移動境界層問題
の適当な境界条件の提言につながるために，多孔質体内の流体浸透や混相流の解析にも寄与すると思わ
れる．
本章では濡れの平衡物理現象について述べた後に，動的濡れ現象の概要を示す．
第 1 章 序論
5
Fig.1.1 平衡接触角と Young の式．
1.2
先行研究
1.2.1 平衡濡れ
固体表面における液体のマクロな濡れ性は，接触角 θ 及び cos θ で評価する事ができる．接触角は液体
表面の接線と固体表面がなす角のうち，液体を含む側の角度のことである．接触角がほぼ 0 度の場合には
完全濡れと呼ばれ，それ以外の場合は部分濡れと呼ばれる．また，ほぼ 0 度において超親水性，90 度未
満を親水性，90 度以上を撥水性，おおよそ 150 度以上で超撥水性とも呼び分ける．接触角 θ は表面自由
エネルギーの平衡論から導出された Young の式 [7] によって，以下のように表される．
γSV = γLV cos θ + γSL
(1.1)
ここで，γ は表面張力または表面自由エネルギーであり，下付き SV， LV， SL はそれぞれ固体-気体間，
気体-液体間，固体-液体間を表す．ただし，今後は簡便のため γLV を単に γ と表現する．Young の式は
濡れ物理の第一歩となる重要な式であり，本式の拡張により更なる理解が可能になる．ここでは Wenzel
の式，Cassie の式，非平衡 Young の式の３つを紹介する．
Young の式は平滑な固体表面をモデル化しているが，Wenzel はこれを拡張し，表面粗さを有する固体
表面上の接触角を提示した [8]．表面粗さによって実際の接触面積が見かけの接触面積以上に増大すると
0
いう概念が基となっており，実際の接触面積と見かけの接触面積の比を A > 1 として，接触角 θw
Young
の式から得られる接触角 θ を用いて以下の Wenzel の式で示される．
cos θ0 =
A (γSV − γSL )
= A cos θ
γLV
(1.2)
本式から，平滑面上での接触角が 90 度以下の場合，すなわち濡れやすい表面の場合は表面粗さによって
より濡れやすくなり，一方で 90 度以下の場合，すなわち濡れにくい表面の場合は表面粗さによってより
濡れにくくなることが分かる．
Cassie は固体表面の濡れ性が不均一な場合へと Young の式を拡張した [9]．平滑面における接触角 θa
と θb の二種の表面が面積比 fa : fb (fa + fb = 1) で存在する場合，接触角 θc0 は以下の Cassie の式で示さ
れる．
cos θc0 = fa cos θa + fb cos θb
(1.3)
第 1 章 序論
6
Fig.1.2 Cassie 状態と Wenzel 状態．
上述の Wenzel の式表面の溝に沿って液体が存在している事を仮定している（Wenzel state）が，一方で表
面の溝に液体が浸透せず空隙層がある場合（Cassie-Baxter state）
，固体表面の一方を空気として Cassie の
式を用いることができる．面 b を空気とし，液体と空気の接触角が 180 度であることを仮定すると，式
1.3 は，
cos θ0 = fa cos θa + fb
(1.4)
と書き換えられる．この式から，固液接触面に空隙層が挟まれる場合に平衡接触角が増加することがわ
かる．
Wentel 状態と Cassie 状態の違いは表面粗さによる溝に液体が浸透するかしないかであり，溝の深さと
長さ，表面の濡れ性によって定まる [10] が，電圧の負荷（electrowetting）により両状態を切り替えられ
るという例も示されている [11]．
また，液滴の接触角が平衡接触角と異なる場合，その非平衡性によって発生する応力を Young の式か
ら求めることができる．平衡接触角 θeq の条件において液滴が角度 θd で接触している場合は，平衡下と
非平衡下の差分を取ることで，
F = γ(cos θeq − cos θd )
(1.5)
として表すことができる．これは濡れのダイナミクスに対する考察において有力な式である．
液滴径が小さく，重力に比べて表面張力が支配的である場合，固体上の液的の概形は真円の一部を切り
取った形になる．しかし，液滴径が大きくなると重力の影響が無視できなくなることから，概形が異なる
ことが知られている．静水圧（重力）と Laplace 圧（表面張力）の比から得られるキャピラリ長 l がこの
指標となり，以下の式で与えられる．
r
l=
γ
ρg
(1.6)
液滴径がキャピラリ長より大きい場合には液滴は細長く伸びた形状となる．水においてキャピラリ長は
2mm 程度であり，これ以上の大きさにおいては重力の影響が無視できなくなる．しかし本研究ではナノ
スケールの液滴を扱うため重力の影響は完全に無視することができ，分子シミュレーションにおいても重
力項の導入はしていない．
先に，接触角が Young の式などによって表現されると述べたが，ナノスケールにおいてはその限りで
はない．Pethica [12] は界面自由エネルギーの総量が固液界面，気液界面，固気界面，ならびに固液気の
接点となる接触線の４つから得られるとし，以下に示される接触角を導出した．
cos θ =
s/γ
s/γ
γSV − γSL
−
= cos θY −
γ
R
R
(1.7)
第 1 章 序論
7
Fig.1.3 先行液膜とミクロ接触角の模式図．
s=
∂F
∂L TVA
(1.8)
ここでの R は液滴の接触半径であり，L は接触線長さである．s は接触線における線張力であり，式 1.8
に示すように三相界面における自由エネルギー F の接触線長さ微分で示される．式 1.7 の二つ目の変形
には式 1.1 に示した Young の式を用いた．R → ∞ において，第二項が消えて Young の式に帰結するが，
R が十分小さいナノスケールにおいては第二項の寄与が無視できなくなり，接触角のサイズ効果が現れ
る．ただし，この第二項は接触線長さ微分が係数となる事からわかるように，接触線が曲率を有する場合
に寄与し，接触線の曲率が無い場合には寄与しない [12, 13]．
また，図 1.3 に示すように，壁面ごく近傍において，液滴本体から逸脱した液膜が発達することも知ら
れている．先行液膜と呼ばれるこの液膜はファンデルワールス力の影響と見られており，この発達は液滴
自体のダイナミクスとは分離して考えられる [14–16]．また，この先行液膜により，マクロに見た接触角
と異なる接触角がミクロスケールで観測され，マクロ接触角（θapp ）とミクロ接触角（θmicro ）によって
液滴の濡れ速度が記述できることも知られている [17, 18]．ただし，ミクロ接触角の同定が困難であるた
め，多くの場合はミクロ接触角を平衡接触角とみなしている．
1.2.2 動的濡れ
動的濡れ性については，支配的な物理現象に応じて，粘性支配下，慣性支配下，重力支配下の三つに分
けて考えられている．先に述べたように，重力の影響の大小は液滴の大きさによって定まるが，本研究で
は重力の影響を無視出来るサイズを対象とするため，ここでは取り上げない．粘性，慣性の大小関係は
Reynolds 数（Re）及び，キャピラリ力を加味した Ohnesorge 数（Oh）によって評価される．Re 数は密
度 ρ，代表速度 v ，代表長さ L，粘性係数 µ を用いて以下のように表される．
Re =
ρvL
inertia force
=
µ
viscous force
Oh 数は，表面張力 γ を加えて以下のように表される．
√
µ
viscous force
We
Oh = √
=
=
Re
inertia · surface tention
ργL
ρu2 L
We =
γ
(1.9)
(1.10)
(1.11)
第 1 章 序論
8
また，Oh 数は代表速度 v をキャピラリ速度 γ/µ と置き換えた Re−1/2 としても理解できる．液体ヘリウ
ム流体の場合には慣性支配，シリコンオイルの場合には粘性支配であるなど，物性に依存して支配則が異
なりもするが，濡れの初期，後期の違いによる濡れ速度の変化によっても支配則が変わりうる．この場
合，濡れ速度が早い初期の段階は慣性支配的であり，濡れが遅くなり平衡状態に近づくと粘性支配的にな
ると言える．ただし，二つの液滴同士が接触し一つの液滴になる現象（コアレッセンス）の議論では速度
場が十分発達していない最初期に粘性支配領域が存在することから，濡れにおいても最初期に粘性支配が
存在することが想定できる [19–21]．なお，コアレッセンスの物理は自然濡れ，特に完全濡れにおける高
速濡れ現象と関わりが深いため，後に慣性支配濡れの説明とともに詳細を示す．
自然濡れの動的挙動に関しては，実験装置の時間解像度に基づく制約から高粘性流体の実験が先立って
行われ，粘性支配領域の物理理解が先に進んだ歴史的背景がある．本論文でも時系列に倣って粘性支配濡
れに関する既存の知見から順に述べる．
粘性支配濡れ
Tanner は気液境界層が表面張力によって駆動される粘性流れであるという仮定のもと Navier-Stokes 方
程式を解き，以下に示す Tanner の式を得た [22]．
θd3 ∼ Ca =
µv
γ
(1.12)
また，バルク粘性散逸と濡れに伴う表面自由エネルギー変化が釣り合うとしても同様の関係式，
2
θd (θd2 − θeq
) ∼ Ca
(1.13)
が得られる [25]．上式で，θd は動的（非平衡時）接触角，Ca = µv/γ はキャピラリ数である．さらに，
ここで液滴半径 R，濡れ半径 r とし，動的過程でも液滴形状が真円形を保ち続け且つ体積が保存されるこ
とを仮定することで成り立つ r 3 θd ∼ R3 （三次元）及び r 2 θd ∼ R2 （二次元）をそれぞれ適用すると，
¶1/10
γ
t
(3D)
r(t) ∼ R
µR
µ
¶1/7
γ
r(t) ∼ R
t
(2D)
µR
µ
(1.14)
(1.15)
が得られる．これらは，無次元時間 τv = µR/γ とおくことによって以下のように整理できる．
r(t)
∼
R
r(t)
∼
R
µ
µ
t
τv
t
τv
¶1/10
(3D)
(1.16)
(2D)
(1.17)
¶1/7
すなわち，粘性支配領域でみた場合，無次元時間 τv = µR/γ で現象が整理でき，さらに，濡れ半径の時
間発展が時間に対するべき乗則に従うことが理解できる．ここでは濡れ半径 r と時間 t の関係式を導出し
たが，同様の考察から動的接触角 θd (t) と時間 t の関係式も導出できる．
この関係式の成立は，三次元のケースは PDMS を用いた実験によって確認されており [23, 24]，二次元
のケースもまた PDMS の実験から確認されている [25]．
第 1 章 序論
9
慣性支配濡れ
濡れの初期において強い非平衡形状を伴った高速濡れと言われる時間域が存在する．高速濡れの時間域
では，液滴の形状が安定とされる真円形を保ち得ないほどに固液接触線が非常に速く進行するために，液
滴形状が常に真円であるという仮定を用いることが不可能になり，濡れ半径と動的接触角の一意性が失わ
れる．そのため，仮に Tanner 則のように動的接触角が動的濡れに寄与すると考えるならば，必然的に液
滴形状全体を解く必要性が生まれ，理論的理解は非常に困難になる．
先に記したが，この条件の観察には高い時間解像度を持つ実験装置が必要であったため，実験が行われ
るのは粘性濡れの後になった．ここでは，二編の実験に関する論文を取り上げて慣性濡れに関する現在の
知見を示すことにする．
Biance ら [26] はガラス基板と半径の異なる水液滴を用いた自然濡れに関する実験を行った．ここでガ
ラスと水の接触角はほぼ０度で完全濡れと言える．ここで彼らは，水の完全濡れ動的挙動が，コアレッセ
ンスの慣性支配時の物理と同一に理解できる，と主張した．
ここで一旦コアレッセンスの慣性支配物理について説明する．コアレッセンスとは二つの液滴同士が接
触し一つの液滴になる現象のことである．粘性の影響を無視し，曲率によって生じる Laplace 圧を駆動力
とした場合，以下の慣性支配則の議論が可能となる [19, 27]．
二つの液滴それぞれの径を R，接触長さを rb とし，液滴の連結部が円弧で近似できるとすると，この
円弧の半径 δ は幾何学的に，
s
µ
µ ¶2 ¶
rb
r2
δ ∼ R−R 1−
∼ b
R
R
(1.18)
と求まる．ここで二式目への変形には rb /R ¿ 1 を仮定した線形化を行った．曲率 κ による Laplace
圧は，
∆P = γκ ∼
γR
rb2
(1.19)
となる．最後に，Bernoulli の式から u2 /2 と P/ρ を連立することによって微分方程式が得られ，特解と
Fig.1.4 コアレッセンスの模式図 [27]．
第 1 章 序論
10
Fig.1.5 Biance ら [26] による完全濡れ条件における初期高速濡れ．図中の三本の線は液滴半径が異な
る同一流体の濡れ時間発展を示す．
して，
r2
t
=
R2
τ
s
µ
¶
ρR3
τ=
γ
(1.20)
を得る．右辺 τ は慣性と表面張力の釣り合い下における時定数となる．ここでは二次元系を想定したが，
三次元系でも同等である [19]．
ここで注目すべきは，粘性濡れにおける時間スケールと慣性濡れにおける時定数が違うことである．粘
性濡れの場合は時定数が R の一乗に比例した一方で，慣性濡れの場合は R の 3/2 乗に比例している．ま
た，濡れ半径の時間に対するべき乗則は双方共に存在するものの，粘性濡れの際には指数が 1/10 であっ
たが慣性濡れの場合には 1/2 になる．
図 1.5 に示す Biance らの実験結果に戻るが，彼らがコアレッセンスの物理と同一に理解できると主張
した理由は三つある．一つ目は濡れ半径が指数 1/2 のべき乗則に従うこと，二つ目は液滴半径を変えたと
きにコアレッセンスの慣性支配スケール則と同じスケール則で整理できること，三つ目は慣性支配から粘
性支配への移り変わる時刻が，コアレッセンスの慣性支配則と Tanner 則の連立から説明できることであ
る．更に，既に濡れた表面上での液滴の動的濡れが乾いた表面上での液滴の動的濡れと大差なかったこと
から，コアレッセンスの慣性支配則が適用できた背景に先行液膜の発達があると主張している．
完全濡れの実験に続き，Bird ら [28] は部分濡れにおける高速濡れ現象の観察を行い，平衡接触角依存
性が動的濡れに及ぼす影響を調べた．図 1.6 には異なる平衡接触角の条件における高速濡れ現象のスナッ
プショットであるが，平衡接触角によって液滴の形状が大きく変わっていることがわかる．しかしなが
ら，平衡接触角 0 度から 120 度までの全ての高速濡れ現象の結果が慣性時間スケール τ =
p
ρR3 /γ で整
理できることを示した．さらに興味深いことに，図 1.7 に示すように，部分濡れにおいても完全濡れ同様
に濡れ半径が時間に対してべき乗則に従うことが，そのべき乗指数は平衡接触角に依存し，0 度では 1/2
だったが 120 度の場合には 1/4 まで減少するという結果になっている．この接触角依存性は依然として理
解が達していないが，界面ごく近傍の現象のみを考慮する場合には 1/2 則しか導出できないことから，部
分濡れの高速濡れ現象においては，固液界面ごく近傍だけでなく，メゾスコピックな領域の散逸を考慮す
第 1 章 序論
Fig.1.6 Bird ら [28] による初期高速濡れ現象のスナップと平衡接触角依存性．R=0.82±0.01 mm．
Fig.1.7 Bird ら [28] による初期高速濡れ現象における平衡接触角とべき乗指数の依存関係．
11
第 1 章 序論
12
Fig.1.8 Carlson ら [6] による平衡接触角 117◦ における初期高速濡れ現象シミュレーションの Dw 依
存性と実験との比較．R=0.78mm．
Fig.1.9 Carlson ら [6] による平衡接触角 117◦ における部分濡れ条件における散逸の割合．Ṙµ は粘
性散逸，Ṙρ は単位時間あたりの運動エネルギーの変化量，ṘD は気液間の散逸，ṘDw は三相界面で
の散逸．R=0.78mm．
る必要があるのではないか，と考察している．
Bird らの実験結果に次いで，Carlson ら [5, 6] は Phase Field 法という連続体計算を用いて液滴の高速濡
れを再現・解析している．三相界面における局所平衡を仮定せず，散逸を三相界面に与えることで Bird ら
の実験結果を良く再現できることを示した（図 1.8）
．図中の Dw が三相界面での散逸量を司るパラメータ
である．また，図 1.9 に示すように，実験と一致する条件下でエネルギー散逸量を計算したところ，この
三相界面における散逸が粘性散逸や運動エネルギー変化に比べてもかなり大きいことを示している．しか
しながら，Dw を直接求めることは未だできておらず，分子レベルでの解析が必要であると主張している．
第 1 章 序論
13
Molecular Kinetic Theory
実験や連続体理論，連続体数値解析によって粘性濡れ，慣性濡れの知見が得られていく一方で，分子レ
ベルでの描像からの理論も存在する．Blake らは接触線において，液体分子が固体上の安定サイト間を確
率的にジャンプしており，これによって接触線の移動を表現できると想定した [29]．MKT によれば，速
度 VMKT は以下の式で示される．
VMKT = λ(K + − K − )
(1.21)
ここで，λ は固体上の吸着サイト間距離，K + と K − はそれぞれ前進，後退確率 [1/s] である．平衡状態
においては K + =K − が成り立つが，非平衡状態においては非平衡 Young の式に準じる応力による差異が
存在しており，エネルギー障壁，活性エネルギーを踏まえて以下の式が成り立つ．
K
+
¶
µ
γ(cos θeq − cos θd )
,
= K0 exp
2nkB T
K
−
¶
µ
−γ(cos θeq − cos θd )
= K0 exp
2nkB T
(1.22)
ここで，K0 は平衡状態におけるジャンプ確率，n は単位面積あたりの吸着サイト数である．式 1.22 を式
1.21 に代入することにより，
µ
¶
γ(cos θeq − cos θd )
VMKT = 2K0 λ sinh
2nkB T
(1.23)
を得ることができ，sinh の中が十分小さい場合，テイラー展開により，
VMKT ∼
K0 λ
γ(cos θeq − cos θd )
nkB T
(1.24)
へと線形化することができる．また，体積保存則と形状変化が準平衡的に起こること，すなわち液滴が完
全に真円形を保ちながら固体表面上を濡れていくことを仮定すると，以下の関係式が成り立つ [30–32]．
r(t) ∼ R6/7 , t1/7
(3D)
(1.25)
r(t) ∼ R
(2D)
(1.26)
4/5
1/5
, t
本理論は，固体上吸着サイト間を行き来する固体表面上の分子運動のみに着目している．言い換えれば，
散逸の箇所を固液気三相界面のみに限定しており，液滴本体における物理現象を考慮する流体力学的理論
とは完全に異なる理論である．
1.2.3 分子動力学法
分子動力学法を用いた濡れの研究例としては，液滴のダイナミクスが MKT に従うと報告した例
[14, 33–35]，Tanner 則に従うと報告した例 [36] が挙げられる．また先行液膜に着目し，r(t) ∼
√
t に従う
もの [37, 38] などもある．一方で初期濡れに関する研究は少なく [39, 40]，より詳細な研究が必要である．
また，表面粗さの濡れに与える影響についての研究例も多いものの [42–45]，自然濡れの時間発展を観察
している例は知る限りではない．
1.3
目的
濡れの静的物理に比べると動的物理は理解が非常に乏しい．近年，液体が固体に接地した直後の高速濡
れ現象の実験による観察例が報告され [26, 28]，慣性スケーリングからの整理が可能であることが示され
第 1 章 序論
14
たが，完全な理解には至っていない．続く連続体計算では固液気三相の接触線における界面現象が重要で
あることが示された [6]．
そこで本研究では，界面現象を原理的に解くことができる分子動力学を用いて高速濡れの系を作成し，
ナノスケール液滴の高速濡れ現象の機構の理解を試みる．現象を支配する物理をスケーリング則を通じて
明らかにすることで，マクロスケールでの高速濡れ現象との関連を考察する．さらに，平滑面上でのシ
ミュレーションを踏まえ，表面構造を変えた場合での濡れ速度を比較することにより，ナノスケールでの
固体表面構造の濡れ現象への影響を考察する．また，分子動力学法と Phase-Field 法の結果の比較を行い，
分子的描像を反映した連続体計算の可能性を検討する．
15
第2章
分子動力学法
分子動力学法は分子の運動を解く手法であり連続体近似を用いないことから，混相問題やナノスケール
の流体問題における第一原理的な手法として広く用いられており，濡れの分野においても数値実験として
多くなされている．本章では，本研究で行った分子動力学シミュレーション手法を説明する．
2.1
支配方程式
分子動力学法では分子の運動は質点の運動として近似し，支配方程式はニュートンの第二法則に従い以
下のように記述できる．
mi
d2 ri
= Fi
dt2
(2.1)
ここで，mi ，ri ，Fi はそれぞれ分子 i の質量，座標，分子 i に働く力である．系に重力などの外力となる
ポテンシャル場がない場合，分子に働く力は他分子とのポテンシャル相互作用のみとなる．このとき，分
子 i に働く力 Fi は他分子 j との分子間ポテンシャル相互作用の総和として近似する．力はポテンシャル
の空間微分で与えられる事から，二体ポテンシャルの場合には，
Fi =
X
i6=j
fij = −
X
∇i φ(rij )
(2.2)
i6=j
として表現する事ができる．ここで，φ(rij ) は座標 ri における分子 i と j の相互作用に基づくポテンシャ
ル，rij は分子 i と j の距離である．
．
本研究では，固液の二相問題を取り扱うため，系には固体分子と液体分子の二種が存在する．そのた
め，固-固，液-液，固-液の三種のポテンシャル関数を定める必要がある．
固-固のポテンシャルについては再近接分子とのみの相互作用を考慮した調和ポテンシャルを用いる．
式としてはポテンシャルは以下で示される．
φ(r) =
1
k(r − r0 )2
2
(2.3)
ここで，k はバネ定数，r0 は原子間の平衡距離である．調和ポテンシャルの使用は Langevin 熱浴を厳密
に用いるため [46–48] であるが，詳細は 2.2 節にて後述する．
液-液のポテンシャルは単純な二体ポテンシャルでありながらアルゴンなどの物性を良く表現する事が
できる Lennard-Jones(LJ) ポテンシャルを用いた．このポテンシャル関数は二原子間の距離を rij として
第 2 章 分子動力学法
16
以下のように示される．
·µ
φ(rij ) = 4ε
σ
rij
¶12
µ
−
σ
rij
¶6 ¸
(2.4)
ここで，ε，σ はそれぞれエネルギーと大きさのスケールである．
固-液についても液-液同様に LJ ポテンシャルを用いて表現する．本研究は仮想的な物質によるシミュ
レーションであるため，LJ 分子同士の相互作用に対して頻繁に適用される混合則 (mixing rule) は用いず，
固-液の LJ パラメータは固-固，液-液とは独立して任意に定める．本研究では σls は固定し，εls を変化さ
せ，固-液間の表面エネルギーを変えることによって，液相バルクの物性を保ちながら異なる濡れ性を表
現し，解析を行う．
LJ ポテンシャルの空間微分で示される力は十分遠方において 0 に漸近する．そのため，分子動力学法
において，相互作用の計算量増大を防ぐために，あるカットオフ距離 rc を設定し，それ以上距離の離れ
た原子間についての相互作用を無視することが往々にしてなされる．一般的に，カットオフ距離は一般に
2.5σ から 5.5σ 程度が用いられており，本研究では 3σ とした．
分子動力学法では，式 2.1 に示した運動方程式を時間積分することで現象の時間発展を解く．本研究
では数値積分法として Verlet’s leap frog 法を用いた．この手法の場合，初めに時間 ∆t/2 ずれた初期値
v(∆t/2) と r(0) を用意することで，以下のようにして時間発展を解くことができる．
µ
¶
µ
¶
∆t
∆t
∆t
vi t +
= vi t −
+
Fi (t)
2
2
mi
µ
¶
∆t
ri (t + ∆t) = ri (t) + vi t +
∆t
2
(2.5)
(2.6)
この Verlet’s leap frog 法を用いる場合，座標と速度の情報が常に ∆t/2 ずれているため，同時刻における
位置と速度を得るためにはいずれかを内挿して計算する．
2.2
温度制御法
統計力学による温度 T は，（熱揺らぎに基づく）速度 v ，運動エネルギー Ek ，粒子数 N ，ボルツマン
定数 kB を用いて，自由度が 3 の場合は以下で示される．
¿X
À
N
1
2Ek
2
T =
mi vi =
3N kB i=1
3N kB
(2.7)
また，分子動力学において式 2.1 に示す Newton の運動方程式はエネルギーを保存するため，系のポテン
シャルエネルギーと運動エネルギーを独立に定めるには温度制御を加える必要がある．
特に，本研究で対象とする液滴の濡れ現象ではエネルギー差が駆動力であるため系全体のポテンシャル
エネルギーは変化するが，ナノスケールにおいては体積に対する表面積の割合が大きく，濡れに伴うポテ
ンシャルエネルギーの変化量が熱容量に対して非常に大きい．これにより，系の温度が濡れの前後で大き
く変化することが予期される．温度の大きな変化は物性の変化をもたらし，また，実スケールでの現象と
の大きな乖離をもたらし得る．これを防ぐために，温度制御法を用いて対処する．
一般的に用いられる簡易な温度制御手法として速度スケーリング法が挙げられる．本手法は，制御前の
温度を T ，目的の温度を Tc としたときに，全ての分子の速度を，
vi0 = vi
p
Tc /T
(2.8)
第 2 章 分子動力学法
17
と繰り返し置き換えることによって，温度制御を実現するものである．簡易かつ頻繁に用いられる手法で
あるが，人為的に統計分布を乱すため，速度スケーリング法を用いながらの現象の観察は避けるべきで
ある．
そこで本研究では，速度スケーリング法の代わりに Langevin 熱浴を固体最下層に設置することで温度
制御を行った [46–50]．この手法は，壁面分子の外側の，液体分子との相互作用が無い距離に Langevin 方
程式に従う phantom 分子を配置することで，無限長で一定温度のバルク固体を仮想する手法である．本
手法は固体分子間が調和ポテンシャルであった場合，phantom 分子が正確にバルク固体層を表現すること
が示されており，これが固-固間を調和ポテンシャルにした理由である．具体的には，phantom 分子の運
動方程式は以下となる．
d2 ri
dri
= Fi − α
+ Fex (t)
2
dt
dt
π
α = mωD
6
mi
(2.9)
(2.10)
ここで，α はダンパ係数，ωD は Debye 周波数であり，Fex は以下に示すガウス分布で与えられるランダ
ムな加振力である．
·
µ
¶2 ¸
1 1
1 Fex
exp −
f (Fex ) = √
2 σF
2π σF
r
2αkB Tc
σF =
∆t
(2.11)
(2.12)
phantom 分子と相互作用する分子との間のバネ定数は，fcc 固体層に本来存在すべき三分子を代替すると
考えて図 2.1 で示す通りに与えられる．またこの手法において，Langevin 方程式に従う phantom 分子は
液体分子と直接的な相互作用が無く，液層の挙動を直接的に乱すことはないと考えられる．
Fig.2.1 Langevin 熱浴による温度制御法 [49, 50]
2.3
並列化手法
分子動力学シミュレーションは単位長さが分子スケールに制限される性質上，大きなスケールのシミュ
レーションには多くの分子数が必要であり，高い計算負荷が課せられる．従来は CPU の性能向上により
計算能力も向上されてきたが，近年では CPU における単一コアの性能向上以上にマルチコア化が進んで
いる背景のため，高速な大規模演算にはプログラムの並列化が不可欠である．n コアで並列演算をした場
第 2 章 分子動力学法
18
合，仮に並列化効率が 100% の理想状態であれば，実際の計算時間は 1/n へと短縮されるため，並列化の
効果は大きい．幸いにも分子動力学法は並列化効率の高い手法であり，数千を超える並列化でも効率を十
分に保つことが可能である．
並列化の実装には，Open MP と MPI の二種のライブラリの利用が考えられる．Open MP の実装は for
文の前に数行指令文を書き足すことで足り，容易に利用できる．しかし，その適用範囲は各コアが同一メ
モリ空間を共有する場合に限られ，ノード間に跨がる大規模並列は不可能である．一方，MPI はプログラ
ム全体の書き直しが必要になるなど実装が困難であるが，マシン間でのメッセージの授受が可能であるた
めに，大規模並列が可能になる．本研究では大規模並列による高速化を可能にするために，MPI を用いた
並列化を行った．
実装には Cell list 法を基にした空間分割を利用した．Cell list 法ではまず，カットオフ距離 rc より長い
rl を設定し，系全体を一辺 rl のサブセルに区切る．これにより，各分子はいずれかのサブセルに属する
ことになり，力の計算をする際には系内の全分子ではなく，同じサブセル内の分子と隣接するサブセル内
の分子とのみ分子間距離の計算を行う．そして，いずれかの粒子が 0.5 × (rl − rc ) 動いた時に．改めてサ
ブセルへの配属をやり直す手法である．並列化を共にすると，Cell list 法は以下のようになる．
まず，各コアに担当サブセルを割り当て，そのサブセルに属する粒子の情報を取得する．例えば，図
2.2 のようなシステムを４コア並列する場合，第一コアは左上に紫色で示した領域を担当する．Cell list
法に従って力の計算を行う際，例えば，赤斜線サブセル内の粒子に働く力を計算するには，図 2.2 のオレ
ンジ枠内の粒子との相互作用を考えなければならない．しかし，各コアは担当領域外に存在する粒子の
情報については無知であるため，他コアが担当する境界領域の粒子位置情報を MPI 関数（MPI ISEND，
MPI IRECV，MPI WAIT）を用いて授受，取得する．これは，図 2.2 に水色で示した領域であり，周期境
界を加味した場合には下側と右側の粒子位置情報の取得も必要となる．この授受の際には x 方向，y 方
向，z 方向と順に三度のみ行うことで通信時間を減らす．速度と位置の更新については，各コアの持つ情
報のみで完結するため問題はない．また，リストの更新にあたるサブセルへの再配属は，単コアの場合と
同様にいずれかの粒子が 0.5 × (rl − rc ) 動いた時に行えばよい．本研究で用いた主な MPI 関数を表 2.1
に記載する．利用する関数の数はさほど多くなく，複雑ではない．
Fig.2.2 Cell list 法
第 2 章 分子動力学法
19
Table 2.1 MPI 関数
2.4
関数
用いた箇所
MPI INIT
最初
MPI COMM SIZE
同上
MPI COMM RANK
同上
MPI BCAST
read したパラメータ情報の送信など
MPI ALLGATHER
位置情報ファイル書き出しの前など
MPI ALLREDUCE
リスト再構築フラグのための総移動距離の計算など
MPI ISEND
位置情報の授受など
MPI IRECV
同上
MPI WAIT
同上
MPI FINALIZE
最後
計算条件
本研究における Lennard-Jones 流体は Ar の物性を再現する εlq = 1.67 × 10−21 J，σlq = 3.4 × 10−10 m
を用いた．固液間パラメータ σls は 3.085 × 10−10 m とし，εls は任意に変更する．固体のパラメータは
r0 = 2.77 × 10−10 m，k = 46.8 N/m とし，fcc(111) 構造をもたせた．系の z 下方に固体層を設置し，z 上
方には反射鏡面条件を，x，y 方向は周期境界条件を課した．これは，z 上方には断熱壁を設置し，x，y 方
向には無限に連続する状態を仮定していることになる．また，時間刻み ∆t は 5 fs，温度は 85 K とした．
2.5
Lennard-Jones
本節では，本研究で用いる温度 85 K における Ar 液体の物性値を求める．
密度
分子動力学法では全ての物質は連続体ではなく離散的な分子の集合として表現される．そのため，密
度 ρ は検査領域を V ，検査領域内に存在する分子数を N ，分子の質量を m とすると，以下の式で与えら
れる．
ρ=
mhN i
V
(2.13)
計算した密度は 1.36 g/cm3 であった．実験による温度 87 K のアルゴンの密度は 1.40 g/cm3 であり，本
結果は妥当であると考える．
拡散係数
拡散係数は Green-Kubo の関係式から，分子速度の自己相関を用いて以下で表すことができる．
1
(Dx + Dy + Dz )
3Z
1 ∞
=
hvi (t) · vi (0)idt
3 0
D=
(2.14)
第 2 章 分子動力学法
20
−9
4
(a)
∫<v(t)v(0)>dt (m2/s)
<v(t)v(0)> (m2/s2)
20000
15000
10000
5000
0
−5000
0
0.5
1
t (ps)
1.5
2
x 10
(b)
3
2
1
0
0
0.5
1
t (ps)
1.5
2
Fig.2.3 拡散係数の計算．(a) 速度の自己相関関数．(b) 自己相関関数の時間積分．
ただし，計算においては積分を有限時間で打ち切る．取得した速度の自己相関関数の例を図 2.3 に示
す．速度の自己相関関数は 0.5 ps 程度で負の極値をとり，1-2 ps で十分収束している．このような負
の値は液相においてのみ現れ，気相においては負値を持たないことが知られている．また，収束時間は
Lennard-Jones 流体の無次元時間 τ =
p
mσ 2 /ε を用いて 1-2 程度であることが知られている．本結果に
おける無次元時間 τ は 2.1 ps であることから，良く一致していると言える．
このようにして計算した拡散係数は 2.2×10−9 m2 /s であった．Moon ら [51] は，アルゴンの密度 1.36
g/cm3 ，温度 94.4 K において，拡散係数 2.44×10−9 m2 /s を報告している．Moon らの方が温度が若干高
いことを考えると，本研究の拡散係数の値は先行研究と比較して妥当であると言える．
表面張力
表面張力の計算は，まず平衡状態における曲率の無い気液界面を作成し，界面における応力テンソルの
垂直成分と平行成分の差分を積分することによって得られる [52, 53]．具体的には，応力テンソルの垂直
成分を PN ，平行成分を PT とし，任意の気体バルク層の位置を zG ，液体バルク層の位置を zL とすると，
以下となる．
Z
zG
γ=
[PN (z) − PT (z)]dz
z
Z LzG ·
=
zL
¸
Pxx (z) + Pyy (z)
Pzz (z) −
dz
2
(2.15)
応力テンソルの計算は Irving-Kirkwood 式より，以下の形で書ける．
P (r) = PK (r) + PV (r)
X
PK (r) =
mi (ui ⊗ ui )δ(qi − r)
i
1X
((qi − qj ) ⊗ Fij ) ×
PV (r) =
2
i6=j
(2.16)
(2.17)
Z
1
δ(λqi + (1 − λ)qj − r)dλ
(2.18)
0
ここで，PK は応力テンソルの運動量成分，PV は分子間力成分であり，分子の位置を q とした．層に
働く応力は，上式の分子間力成分の計算において，二分子が対象となる層を横切る場合のみを考慮す
る [53, 54]．結果を図 2.4 に示す．いずれも，横軸は位置，縦軸は応力である．運動量成分はおおよそ密
度に比例するため，液相で高く，気相ではほぼ 0 であり，平行成分，垂直成分における値の差はない．し
第 2 章 分子動力学法
21
25
10
0
Pressure (MPa)
Pressure (MPa)
(b)
(a)
20
15
10
P
xx
5
−10
−20
P
−30
P
yy
P
yy
P
zz
0
20
xx
P
zz
25
30
z (nm)
35
−40
20
40
25
30
z (nm)
35
40
5
Pressure (MPa)
(c)
0
−5
−10
P
xx
P
yy
P
zz
−15
20
25
30
z (nm)
35
40
Fig.2.4 気液界面における応力テンソル.(a) 運動量成分，(b) 分子間力成分，(c) 運動量成分と分子間力成分の総和.
かし，分子間力成分は，液相において，ちょうど運動量成分を打ち消す程度の値を持ち，気液界面にお
いては垂直成分と平行成分において差が生じる．このようにして計算した表面張力は 10.7 mN/m であっ
た．Baidakov ら [55] による Lennard-Jones 流体を用いた計算では，アルゴンに換算すると，温度 87 K，
密度 1.41 g/cm3 において表面張力 15.1 mN/m であると報告している．この差は，本研究ではカットオ
フ長さを 3σ としたが，Baidakov らは 6.7σ としていることに由来すると考えられる．Cosden ら [56] は
Lennard-Jones 流体を用いてカットオフ長さが表面張力に及ぼす影響を調べ，カットオフ長さ 3σ と 6σ で
は計算される表面張力に 1.5 倍の差があり，カットオフ長さ 10σ 程度でもカットオフ長さの依存が残るこ
とを報告している．
粘性係数
本研究では，以下に示す Srokes-Einstein 則を用いて粘性係数を算出した．kB はボルツマン定数，T は
温度，a は分子半径であり，D は先に Green-Kubo を用いて算出した拡散係数である．分子半径 a には，
Lennard-Jones パラメータを用いて，a = σ/2 とした．
µ=
kB T
6πaD
(2.19)
第 2 章 分子動力学法
22
このようにして計算した粘性係数は 0.16 mPas であった．一方，Meier ら [57] による Green-Kubo を用い
た計算では，密度 1.42 g/cm3 ，温度 86 K，カットオフ 3σ において 0.23 mPas であると報告している．こ
れは手法による差であると考えられ，Green-Kubo を用いた直接計算を行い粘性係数の整合性を確かめる
必要がある．最後に，以上の計算で求めた Ar 液体 85 K の物性値を表 2.2 に示す．
Table 2.2 Ar 液体 85K における物性値
密度
1.36 g/cm3
拡散係数
2.2×10−9 m2 /s
表面張力
10.7 mN/m
粘性係数
0.16 mPas
23
第3章
平滑面上の動的濡れ性
本章では fcc 固体 (111) 上における Lennard-Jones 流体の動的濡れ性に関する分子動力学シミュレー
ションの結果の提示と考察を行う．
3.1
計算系
本研究では主に図 3.1(a) に示すような疑似二次元液滴を作成し，計算を行う．y 方向の周期境界によ
り，作成される液滴は y 方向に無限の長さを有した円筒状の形状を持つ．三次元液滴では液滴の大きさが
分子数の三乗に比例する一方で，二次元液滴を用いた場合では二乗に比例する．そのため，二次元液滴を
用いる事で計算負荷を軽減しつつ大きな液滴の作成が可能になる．なお，二次元系との整合性の検証のた
めに図 3.1(b) に示すような三次元径の解析も限られたパラメータに対して行った．
本研究の計算で用いた計算領域長さと分子数，液滴半径を表 3.1 に示す．計算系の x 方向長さと z 方向
長さは用いる液滴の半径に対し十分な長さになるよう随時変更しているが，y 方向長さに関しては本論文
中を通して常に一定としている．
シミュレーションでは初めに液体分子を固体表面から十分離れた場所に fcc 構造で配置し，85 K で円
形または球形の液滴が作成されるまで計算を行う．その後，固-液間ポテンシャル εl s = 5.82 × 10−23 で
液滴を固体に接近させ，接触角 180◦ で接触している状態を作成する．最後に，固-液間ポテンシャルを任
Fig.3.1 (a) 擬似二次元系，(b) 三次元系．
第 3 章 平滑面上の動的濡れ性
24
Table 3.1 計算条件．
液体分子数 N
25200
108000
432000
液滴半径 R
5.6 nm
11.2 nm
22.4 nm
x 方向長さ Lx
55.4 nm
123.7 nm
110.8 nm
y 方向長さ Ly
10.8 nm
10.8 nm
10.8 nm
z 方向長さ Lz
15.0 nm
31.3 nm
60.0 nm
意の値に変更することによって自然濡れを開始し，観察を行う．
3.2
平衡接触角
液滴が固体表面上に接触後，平衡となった時のスナップショットの例を図 3.2 に示す．このような分子
の分布を基に平衡接触角を計算する．平衡接触角の算出には，まず密度分布を取得し，その後気液界面を
抽出し，最後に界面を円弧でフィッティングすることによって得られる．
密度分布
密度分布の取得のため，初めに系全体を幅 δx，δz のメッシュに区切る．各メッシュにおける局所の密
度は y 方向におけるアンサンブル平均を考慮することで，メッシュ内に存在する分子数 Nl と分子の質量
m を用いて以下の式で得る事ができる．
m
ρ(x, z, t) =
Nl
δxδzLy
(3.1)
ただし，スケールが小さく粒子化された MD シミュレーションでは，液滴がブラウン運動的に固体表面
上を動き回る．そこで，メッシュは系に対して絶対座標として与えるのではなく，液滴の重心を基準にし
た相対座標として与える．液滴の重心 Cx は，
Cx =
N
1 X
xi
N i
(3.2)
の様な与え方では周期境界条件を考慮すると不適切であり，Bai-Breen の方法に倣って重心を計算し
た [58]．この手法によれば，分子の総数を N ，系の x 方向長さを Lx ，i 番目の分子の位置の x 成分を xi
Fig.3.2 固液接触後の平衡状態．R=12.2nm，θeq = 90◦
第 3 章 平滑面上の動的濡れ性
25
とすると，重心 Cx は，
µN
µ
¶¶
1 X Lx
2πxi
X=
cos
N
2π
Lx
i
µN
µ
¶¶
2πxi
1 X Lx
Z=
sin
N
2π
Lx
i
µ
¶
Lx
−1 Z
Cx =
π + tan
2π
X
(3.3)
(3.4)
(3.5)
で表すことができる．なお，三次元系においては y 方向の重心を求める必要があるが，x 同様にして求め
る事ができる．
得られる密度分布の解像度はメッシュ幅に依存するため，高解像度にデータを得るためには細かいメッ
シュが必要になる．しかし，メッシュが細かい場合には，式 3.1 に示した様な手法の場合，全情報数は分
子数 N に留まるために統計誤差が無視できない．そこで，分子に半径 σ/2 の有限の体積を与え，一分子
の情報を複数のメッシュに寄与させた．
気液界面の抽出
気液界面における密度を図 3.3 に示す．横軸は位置，縦軸は液相バルクの密度を用いて無次元化した密
度である．気液界面においては，液体バルクの密度 ρL から気体バルクの密度 ρG までなだらかに変化し，
有限長の界面が存在することがわかる．そこで界面を，
ρinterface =
ρL + ρG
2
(3.6)
で与えられる密度を有する点であると定義して界面の抽出を行う．
1
ρ/ρlq
0.8
0.6
0.4
0.2
0
7
8
9
z [nm]
Fig.3.3 気液界面における密度分布．
10
第 3 章 平滑面上の動的濡れ性
26
円弧によるフィッティング
気液界面が円弧となるとし，非線形最小二乗法を用いてフィッティングする．フィッティングの方程式
は，円の中心を (0, yc )，円の半径を rfit とすると，
¶2
N µq
X
2
2
S=
xi + (yi − yc ) − rfit
(3.7)
i
∂S
=0
∂yc
∂S
=0
∂r
である．この最小値問題は解析的には解けないため，Newton 法を用いた反復計算から数値的に求めた．
また，フィッティングに関しては先行液膜の影響を排除するために液滴下層部の密度，界面情報は用い
ない．その後，フィッティングした円弧を固液界面まで延長し，交点における接触角を得る．ただし，本
フィッティングによる平衡接触角の同定は気液界面が円弧で示されることを前提としている．そのため，
本研究の物性値下では濡れ現象の最中においては，気液界面は円弧ではないためにこの手法を用いた動的
接触角の導出はできない．
ここで，三次元系の場合はフィッティング式が異なることを示しておく．三次元系において，系全体を
δr 及び δz のメッシュに区切り，同様の観察を行った場合，各メッシュの体積は，2πrδrδz となる．すな
わち，メッシュ体積は半径 r に比例して非一様となり，取得される統計的分散量は 1/r に比例して減衰す
る．そのため，統計量の信頼性に合わせてフィッティング式には重み付けがなされるべきであり，式 3.8
は以下の様に書き換わる．
S=
N
X
µq
¶2
2
2
xi
xi + (yi − yc ) − r
(3.8)
i
なお，フィッティングに重み付けを行うのではなく，メッシュの体積を一定にすることで式 3.8 の適用を
保つことも考えられる．
1
cosθeq
0.5
0
−0.5
−1
0.5
1
ε/εref
1.5
2
Fig.3.4 平衡接触角と固液間 Lennard-Jones パラメータの関係性．R=12.2 nm，εref = 0.582 × 10−21 J．
第 3 章 平滑面上の動的濡れ性
27
以上手法を用いて算出した平衡接触角と εls の関係を図 3.4 に示す．横軸は平衡接触角 90◦ を与える ε
で無次元化した εls ，縦軸は平衡接触角の余弦である．図 3.4 では液滴の大きさを変えていないが，平衡接
触角に関して，二次元液滴の場合には線張力に伴うサイズ効果が無いことが理論と MD シミュレーショ
ンにおいて示されている [12, 13]．また，本結果は同じ系で行われた先行研究結果 [49] と一致する結果と
なっている．
3.3
動的濡れ性
本節では液滴が固体表面に接触してから平衡状態になるまでのダイナミクスに関する結果を示し，考察
する．濡れの進行におけるスナップショットを図 3.5 に示す．本条件では濡れの進行中の液滴形状は，特
に接触直後において，真円形状からは大きく逸脱していることがわかる．本節では，濡れ半径の時間発展
の観点から，動的濡れ性に関する考察を行う．
3.3.1 サイズ依存性
動的濡れ性の機構を理解する一つの方法にスケーリング則の同定が考えれられる．そのために，連続体
解析であれば変更するパラメータは粘性係数等の物性値も考えられるが，MD シミュレーションにおいて
は物性値を任意に与えることは容易でなく，任意に変え易いパラメータの一つとして液滴の大きさが挙げ
られる．ここでは，その際に増加する計算コストを緩和するため，疑似二次元系を用いてシミュレーショ
ンを行う．
液滴半径 5.6 nm，11.2 nm，22.4 nm の三種類の液滴の濡れ半径の時間変化を図 3.6(a) に示す．横軸は
時間，縦軸は濡れ半径である．この結果から，液滴半径が大きい方が濡れ速度は速く，また，平衡接触角
が小さい方が，すなわち濡れ性が高い方が濡れ速度が速いことがわかる．ただし，平衡に至るまでの時間
は液滴半径が大きい方が長く，また，濡れ性が高い方が長い．
この結果を，縦軸の濡れ半径 r を液滴半径 R で，横軸の時間 t を τ = (ρR3 /γ)1/2 で割ることによ
り，慣性スケーリング下で整理し直したものを図 3.6(b) に示す．この慣性スケーリングを施すと，Biance
ら [26] や Bird ら [28] のサブミリスケール実験で示されたように，異なる液滴半径のデータを同一曲線で
整理できることがわかった．すなわち，本条件での液滴の濡れのダイナミクスが粘性的な影響以上に慣性
的な影響によって支配されていることが示唆される．
実際に，算出した物性値を基に Reynolds 数を計算する．Reynolds 数は慣性スケーリングによるコア
Fig.3.5 平衡接触角 20◦ の液滴の動的濡れ現象における液滴形状．R=12.2 nm．
第 3 章 平滑面上の動的濡れ性
28
Fig.3.6 濡れ半径の時間発展．(a) 実スケール，(b) 慣性スケール．
レッセンスの議論を基にすると，代表速度 U は，
s
U∼
P
∼
ρ
r
γ
ρR
(3.9)
となる．この代表速度は図 3.6(b) における傾き 1 に相応し，グラフからこの定義は妥当であると言える．
この時，Reynolds 数は，
s
ρU L
Re =
=
µ
ρLγ
µ2
(3.10)
となる．ここに 2.5 節で得た物性値を代入すると，用いた中で最小の液滴半径 6.1 nm で Re=2.6，12.2 nm
で Re=3.7，12.2 nm で Re=5.3 となり，12.2 nm から 24.4 nm の全ての液滴において，1 以上の値が得ら
れた．結論として，Lennard-Jones 流体を用いることによって，サブミリメートル実験で見られるような
慣性支配挙動がナノスケール液滴を用いた MD シミュレーションにおいても再現でき，その理由は Re 数
から説明可能であると考えられる．
また，図 3.7 に慣性スケールにおける無次元時間 t/(ρR3 /γ)1/2 = 1 における各種条件のスナップショッ
トを示す．図 3.6 では接触半径をプロットしたのみであったが，液滴の概形を見ると多少の誤差がある．
液滴半径が小さいものはより潰れた形となっている一方で，大きい液滴の場合は液滴上部が伸びた形に
なっている．この形状の差異は，固液接触線の移動に関する物理現象の時間スケールと液滴の形状変化の
時間スケールが異なっていることに原因ではないかと思われる．それにも関わらず同一のスケール則で濡
れ半径の進行が整理できるということは，全体の液滴形状と濡れ半径の物理現象は概ね独立に考えられる
可能性を示唆している．
第 3 章 平滑面上の動的濡れ性
29
Fig.3.7 無次元時間 t/(ρR3 /γ)1/2 = 1 における液滴形状．(a)θeq = 20◦ ，R=6.1 nm，(b)θeq = 20◦ ，
R=12.2 nm，(c)θeq = 20◦ ，R=24.4 nm，(d)θeq = 90◦ ，R=6.1 nm，(e)θeq = 90◦ ，R=12.2 nm，
(f)θeq = 90◦ ，R=24.4 nm
3.3.2 平衡接触角依存性
ここでは平衡接触角を変えた際の接触半径変化について述べる．平衡接触角は図 3.4 に示した 20◦ ，
45◦ ，90◦ ，135◦ の四種類を用いる．また，液滴半径は 12.2 nm を用いた．
図 3.8(a) に慣性スケールで整理した濡れ半径の時間発展を示す．横軸が無次元時間，縦軸が無次元長さ
である．濡れ性が高いほど，濡れ速度が速いことがみてとれる．この図を両対数でプロットし直したもの
を図 3.8(b) に示す．よく知られているような r/R = C(t/τ )α で示されるべき乗則に従うことが見て取れ
る．これは，図 3.9 に示すように，三次元系でのシミュレーションでも同様に見られた．
図 3.10(a) に二次元系，三次元系計算より得られたべき乗指数の平衡接触角への依存性を示す．二次元
系，三次元系ともに，平衡接触角が高いほど指数が小さくなる傾向が一致している．また，平衡接触角が
小さい場合にはコアレッセンスの議論から示される 0.5 程度になっている．これは，Bird ら [28] の実験
結果と定性的に一致する結果となっている．ただし，二次元系同士では液滴半径と指数の定性的な関係性
は見られていない．そのため，何度か同じ系での計算を行い比較することにより定性的議論を確かにする
Fig.3.8 濡れ半径の時間発展の平衡接触角依存性．R=12.2 nm．(a) リニアスケール，(b) 両対数スケール．
第 3 章 平滑面上の動的濡れ性
30
Fig.3.9 濡れ半径の時間発展の平衡接触角依存性．三次元系，R=10.2 nm．(a) リニアスケール，(b) 両
対数スケール．
R=6.1nm(2D)
R=12.2nm(2D)
R=24.4nm(2D)
R=10.2nm(3D)
0.5
(b)
1.4
C
α
1.2
0.4
1
0.3
0
90
θ eq
(a)
180
R=6.1nm(2D)
R=12.2nm(2D)
0.8 R=24.4nm(2D)
R=10.2nm(3D)
50
100
θ eq
150
Fig.3.10 べき乗則係数の平衡接触角依存性．(a) 指数 α，(b) 比例定数 C ．
必要があると考えられる．フィッティング誤差に関しては，データ点の細かさとそもそもの統計誤差の少
なさから液滴半径 24.4 nm のデータの信頼性が高いが，三次元の結果との良い一致が見られる．
比例定数部を見ると，二次元系の場合は接触角が高くなるにつれ減少しているが，三次元系の場合は多
少の増加あるいは同程度であるといえ，傾向に不一致が見られる．これは比例定数部の物理的意味を把握
し，考察する必要がある．
3.4
流速場
エネルギー散逸の方程式は以下の式で示される．
µ
¶
Z
d
1 2
FLG + FSL + FSG +
ρu dΩ
dt
Ω 2
Z
Z
Z
T
2
+
µ∇u : (∇u + ∇u )dΩ + βSL uSL dΓ + βcl u2cl dl = Q̇out
Ω
Γ
l
(3.11)
第 3 章 平滑面上の動的濡れ性
31
上式にて F は界面エネルギー，β はスリップ係数，Q̇out は系からの流出熱，u は速度であり，積分区間
は Ω が体積積分，Γ が固液界面の面積分，l が固液気境界線の線積分である．下付きは LG が気液界面，
SL が固液界面，SG が固気界面，cl は三相界面である．本節では，上式を用いた速度場からのエネルギー
考察を目的とすると共に，各項の寄与を定量的に評価する．ここでは，液滴半径 R=12.2 nm の液滴，平
衡接触角 90◦ の条件を用いる．平衡接触角に 90◦ を用いた理由は，部分濡れの議論が可能であり，
d
(FSL + FSG ) = 0
dt
(3.12)
が成り立つことから，式が容易になるためである．速度場を取得するには，系をメッシュに切り，メッ
シュ内に存在する分子速度を計算することによって得ることができる．ただし，分子の運動には熱揺らぎ
に基づくランダムな速度成分が存在するため，流動速度のみを求めるのは極めて難しい．そこでマクロに
は同一条件であるが各分子の速度と位置が異なる初期条件を作成し，それぞれの計算結果をアンサンブル
平均した．ここで，流動成分は異なる初期条件でも保たれる一方で，熱揺らぎの速度ベクトルは保存され
ないことから，アンサンブル平均を繰り返すことで熱ゆらぎ成分を低減し S/N 比を高めることができる．
本研究では 10 回の試行を平均化したがそれでもノイズが多く，データの解析には堪えなかった．そこ
で時間に対して 30 GHz でのローパスフィルターを掛け，空間に対して隣接平均をとることによってノイ
ズを減らした．図 3.11 に取得した速度場を示す．ここで，速度ベクトル場は液相にのみ表示した．なお，
気相に速度ベクトル場を表示しなかったのは分かりやすさのためのほか，気相では分子数が少ないために
アンサンブル量が不十分になりがちで，速度場のノイズが大きくなることも理由である．なお図におい
て，第一層において下向きの速度が生じているのは空間隣接平均の弊害である．細かく見ると，渦構造や
微かな流れ場が見えるなど，まだノイズの影響が残っているが，概ね液滴内の流速分布の傾向を掴むこと
は可能である．
液滴が固体に接地した直後 10 ps 後の速度場をみると，液滴下部に強く下方向の速度場が生じているこ
とがわかる．一方で，z =10 nm 程度以上においてはまだ流れ場が波及していない．また，気液界面にお
いては，液滴内部へ向かうような速度場が発達している．
100 ps 後の結果をみると，10 ps 後に比べて，下層だけでなく全体的に速度場が広く分布している様子
が見て取れる．気液界面近傍においては，z =12 nm 近傍においては液滴中心方向へ向かう成分が強く見
られるが，z =5 nm 近傍においては速度場が気液界面に沿った方向に流れている様子が見られ，z =3 nm
のあたりでは速度場は気相側を向き始めている．固液界面においては固液気三相界面における速度場が大
きく発達し濡れを進行させている様子が観察でき，また，固液界面であっても液滴中央に近づくに連れ速
度は x 方向速度は遅くなっており，x =0 の当たりには速度が生じていない様子が見られる．
1000 ps 後の結果を見ると，液滴中央の液体が流動しない場所が広がっている様子が見られる． x =5
nm 未満程度の液滴中央においては，液体上方の方がむしろ盛んに流れが生じている．気液界面に注目す
ると，z =15 nm 程度においては速度場は界面と平行方向を向いており，z =10 nm 未満では完全に気相
方向を向いている．三相界面において，液体の濡れを進行させる速度場は 100 ps の結果と比べると減少
しており，100 ps 後の結果同様に液滴中央に向けてなだらかに x 方向速度が遅くなっている様子を見る
ことができる．
次に，エネルギーの分配からの考察を進める．以降では y 方向長さは単位長さ 1 とする．
次式を用いて速度分布から系の運動エネルギーを求める．
Z
K(t) =
Ω
ρ
u(t)2 dΩ
2
(3.13)
第 3 章 平滑面上の動的濡れ性
(a)
25
25
20
20
15
10
5
5
0
0
0
0
10
x (nm)
20
25
20
20
z (nm)
25
10
5
5
10
x (nm)
20
20
(d)
15
10
0
0
10
x (nm)
30
(c)
15
(b)
15
10
30
z (nm)
30
z (nm)
z (nm)
30
32
0
0
10
x (nm)
Fig.3.11 速度場．(a)10ps，(b)100ps，(b)500ps，(d)1000ps．
20
第 3 章 平滑面上の動的濡れ性
33
Fig.3.12 運動エネルギーの時間変化．(a) 運動エネルギーの総量，(b) 運動エネルギー変化．
ここで，積分区間は液相のみとする．図 3.12(a) に系の運動エネルギーの時間変化を示す．図中，Kx ，
Ky ，Kz はそれぞれ x 成分，y 成分，z 成分を示し，Ktotal はこれら三成分の総和である．擬似二次元系
を用いているため Ky は 0 になるべきであり，ここでの Ky の値は単にノイズの指標の一つであると考え
てよい．系の運動エネルギーは液滴が設置した直後から急速に上昇し，無次元時間にして 1 程度において
最大となる．この時，液滴半径 r/R はおおよそ 0.8 程度である．運動エネルギーの成分としては，z 方向
成分が支配的であり，x 方向成分に比べ概ね三倍程度のエネルギー量が占められている．
図 3.12(b) に次式で示される液滴の運動エネルギー変化の時間変化を示す．
d
R˙ρ (t) =
dt
Z
Ω
ρ
u(t)2 dΩ
2
(3.14)
ここで，時間微分は空間積分の後に課した．これは，微分によるノイズの増幅を防ぐためである．
x 成分を見ると，0.1 未満においては強く運動エネルギの勾配が存在するものの，その後は無次元時間
1.3 程度までは 0 となっている．すなわち運動エネルギーの x 方向成分は，最初の無次元時間 0.1 までの
間で得られた量が長い間保存されることを示している．その後は小さな負値をとる．z 成分を見ると，無
次元時間 0.8 程度まで概ね正値を，1 程度からは負値をとっている．なお，x 成分，z 成分 共に無次元時
間 0.1 程度において 0 となっている．この無次元時間 0.1 は局所接触角が概ね平衡接触角と同程度となる
時刻と一致するが，関係性は定かではない．また，ローパスフィルターによる悪影響である可能性も考え
られる．
次式で示される気液表面エネルギーの時間変化を図 3.13 に示す．
FLG (t) = γS(t)
(3.15)
ここで S は気液界面の面積である．表面張力は定数として扱うため，単に表面積の時間変化を示してい
ることと同義となる．最初は固体に接触していないため，気液表面エネルギーは高いが，濡れが進行する
につれて気液表面が固液表面へと変わるために気液表面エネルギーが減少する．無次元時間 1 程度を超え
ると，再度気液エネルギーが増大する．この表面エネルギーの変化が，先の運動エネルギー散逸と次に示
す粘性散逸に対応している．
次に，粘性散逸を以下の式を用いて取得する．
Z
µ∇u : (∇u + ∇uT )dΩ
Ṙµ =
Ω
(3.16)
第 3 章 平滑面上の動的濡れ性
34
−10
FLG (J)
4
x 10
3.5
3
2.5
0
0.5
1
1.5
3 1/2
t/(ρR /γ)
2
2.5
Fig.3.13 表面エネルギーの時間変化．
Fig.3.14 粘性散逸．(a) 時間変化，(b) 接触直後からの積算値．
得られた粘性散逸の時間変化を図 3.14(a) に示す．無次元時間 0.1 程度で一旦の極値をとり，その後，無
次元時間 1 程度で最大値を示す．図 3.14(b) に設置直後からの積算の値を示す．無次元時間 2.5 までの積
算量は気液表面エネルギー変化量の約 7 倍にも達しており，このままではエネルギー方程式からの不定項
の導出は不可能であることがわかった．これは，ノイズの含んだ速度場を空間微分する必要がありノイズ
が増幅されやすいことに起因する．
Langevin 熱浴からの系への流入エネルギーの時刻 0 からの積算値を図 3.15 に示す．用いている
Langevin 熱浴ではホワイトノイズの加振力によって系内外へとエネルギーの授受を行う．そのため，時々
刻々の値としてはほぼノイズしか見えず，体をなさない．最終的には系からエネルギーを取り除いている
結果となっているものの，他の散逸項とはオーダーが異なっており無視できる程度である．本項では，エ
ネルギー方程式を用いた未定項の導出と物理理解を目的に速度場の算出を行った．速度場は，時間に対し
てローパスフィルターを，空間に対して隣接平均をとることで定性的な傾向の観察が可能な程度のデータ
にはなったものの，粘性項を正しく導出できるには至らなかった．ただ，無次元時間において 0.1 と 1 に
おいて散逸項の変化が見られ，前者が局所接触角の平衡化と，後者が界面エネルギーの極値と対応してい
るのではないか，と推察した．
第 3 章 平滑面上の動的濡れ性
35
−12
Total energy inflow (J)
4
x 10
2
0
−2
−4
0
0.5
1
1.5
3 1/2
t/(ρR /γ)
2
2.5
Fig.3.15 熱浴からの流入エネルギーの積算値．
3.5
まとめ
本章では fcc 固体 (111) 上における LJ 流体の動的濡れ性に関する分子動力学シミュレーションの結果
を示した．液滴を二倍，四倍の大きさに変えたシミュレーション結果が慣性スケーリング下で同一曲線上
に乗ることから，ナノスケールにおいても液体の物性値をコントロールすることで，実験において示され
る慣性支配挙動を再現可能であることを示した．また，平衡接触角を変えたシミュレーションの結果から
も，実験で示されるようなべき乗則と接触角依存性を見出すことができた．しかしながら，そのべき乗指
数の由来については理解することができなかった．
36
第4章
Cassie 条件における動的濡れ現象
4.1
概要
本章では，表面エネルギーが非一様な固体上における液滴の平衡接触角と動的濡れ現象について議論
する．
本研究では，非一様な表面を実際に表面粗さを与えるのではなく，固液間の表面エネルギーが強い場
所と弱い場所を作成することで，図 4.1 に示すように，場所によって濡れ性の異なる固体表面を作成す
る．このような固体表面は，理論的に Cassie 条件として理解される．具体的には，図 4.1 において，茶色
(wetting) で示す位置の固体原子の固液間 LJ パラメータを ε1ls に，青色 (non-wetting) で示す位置の固体原
子の固液間 LJ パラメータを ε2ls に設定する．この時，ε1ls > ε2ls にすることで，疎水面と親水面が表面上
に混在する Cassie 条件の表面を作成することができる．本研究では簡単のため，一つの親水面・疎水面の
長さは常に一定に固定した．すなわち，疎水面長さを ∆x1 ，親水面長さを ∆x2 とした時，(∆x1 + ∆x2 ）
の周期で疎水面と親水面が連続する．
Fig.4.1 Cassie 表面の模式図．
4.2
平衡接触角
初めに，固液間ポテンシャルを操作することによって作成した Cassie 表面上において，液滴の平衡接
触角が理論的な Cassie の式と一致するかを議論する．再掲になるが，二種の表面が混在する場合におけ
る Cassie の式からは，面積比 fa :fb (fa +fb =1) で各々の平衡接触角が θa と θb の場合，平衡接触角 θ は以
第 4 章 Cassie 条件における動的濡れ現象
37
下の式で求まる．
cos θ = fa cos θa + fb cos θb
(4.1)
すなわち，理論的には，疎水面における平衡接触角，親水面における平衡接触角，表面積比の三種のパラ
メータで平衡接触角が定まる．しかし，ナノスケールにおいて疎水面の長さ，親水面の長さの影響は定か
ではない．そこで，分子動力学シミュレーションにおいて疎水面の長さと濡れ性，親水面の長さと濡れ
性，それと液滴半径の五つを変えて平衡接触角を取得した．
図 4.2 に，計算で求めた平衡接触角と Cassie の式との比較を示す．縦軸には cos θ を，横軸には疎水面
の長さをとった図中，白抜きは液滴半径 6.1 nm，塗りつぶしの丸は 12.2 nm，ひし形で示したものは 24.4
nm である．緑で示したデータ点は疎水面の長さと親水面の長さの比が 2:1 のものであり，他は長さの比
が 1:1 である．赤色と緑色で示したデータ点は親水面接触角と疎水面接触角が各々 20◦ と 180◦ ，青色は
45◦ と 135◦ ，黄色は 90◦ と 180◦ である．また，Cassie の式から得られる接触角は図中に同一色の直線で
示した．直線になるのは，Cassie の理論式では表面積比のみが重要であり，親水面長さと疎水面長さには
依存しないためである．
分子動力学シミュレーションから得られた結果は，Cassie の式から得られる値とさほど遠くないが，疎
水面・撥水面の長さに依存した平衡接触角の変化が見られ，特に疎水面幅が長い領域においては差が大
きい．
Cassie の式は Young の式同様に，エネルギーの平衡論が基になっている．すなわち，自然濡れの過程
で局所的極小点に落ち込む可能性については言及していない．本研究での結果では，自然濡れの帰結と
してそのような極小点に収束したために Cassie 条件に一致しない結果が得たのではないかと考えられる．
確かに定性的に見て，撥水面長さが長いほど高い接触角で留まっており，濡れ広がりが止まっていること
が読み取れる．
次節では，このような不均一表面上におけるダイナミクスの議論を行う．その理解は，ダイナミクスの
理解それ自体はもとより，ナノスケール Cassie 面上の平衡接触角に関する理解にも不可欠である．
0.5
cos θ
20°−180° (surface fraction 2:1)
0
−0.5
0
45°−135°
20°−180°
95°−180°
2
4
6
Distance of Hydrophobic Patches [x 2.77Å]
Fig.4.2 Cassie の式との比較．
第 4 章 Cassie 条件における動的濡れ現象
4.3
38
動的濡れ性
親水面，疎水面の長さを変えた時，どのようにダイナミクスが異なるかを議論する．図 4.3 に半径 12.2
nm と 24.4 nm，親水面の平衡接触角が 20◦ ，疎水面の平衡接触角が 180◦ ，親水面と疎水面のの長さの比
が１，すなわち ∆x1 = ∆x2 = ∆x における ∆x 依存性を示す．r0 は固体の分子間距離で，2.77Å であ
る．また，図中には平滑面上 90◦ の結果も付記した．横軸は慣性スケール下の無次元時間，縦軸は無次元
半径である．
液滴半径 12.2 nm の結果である図 4.3(a) を見ると，∆x が短い時には半径推移は同じ平衡接触角となる
平滑面上の結果とほぼ等しくなることがわかる．一方で，∆x を長くしていくと，定性的に濡れ速度は遅
くなり，徐々に平滑面上の結果との差異が生じてくる．ただし，∆x = r0 と ∆x = 2r0 のグラフを比べる
と，若干ではあるが ∆x = 2r0 の方が早い．また，幅が長い場合のグラフからは階段状に推移しているこ
とがよくわかるが，これは親水面上から疎水面上へと濡れが進行する際に stick-slip 現象が起こるためで
ある．
液滴半径 24.4 nm の結果である図 4.3(b) を見ると，Cassie 表面上の方が平滑面 90◦ と比べてより早く
濡れる部分が明らかに現れている．ただし，∆x が長くなるにつれて濡れの速度が遅くなる傾向が見られ，
また，stick-slip 現象に基づく顕著な階段状構造が観察できる．
このグラフから見られる Cassie 表面上の結果は以下にまとめられる．
• R=12.2 nm の場合，平滑面上で同程度の接触角の結果と概ね一致する．
• ∆x が大きくなると，濡れ速度は遅くなる．
• ∆x=4r0 以上において stick-slip 現象が強く観察できる．
• R=24.4 nm の場合，∆x=4r0 未満では平滑面上での同程度の接触角の結果より濡れが速い．
次項から，これらの理由の考察を行う．
Fig.4.3 Cassie 表面上における液滴の濡れ半径推移．(a)R=12.2 nm，(b)R=24.4 nm．
第 4 章 Cassie 条件における動的濡れ現象
39
4.3.1 定性的理解
図 4.4(a) に半径 12.2 nm の液滴におけるピニング時間 ∆t を示す．横軸は液滴の中心からの距離，縦軸
はピニングが起こっている時間 ∆t である．さらに，このグラフの縦軸，すなわちピニング時間を ∆x2 で
除したものを図 4.4(b) に示す図を見てわかるとおり，三種の異なる ∆x におけるピニング時間を整理す
ることができたと言える．すなわち親水面と疎水面の界面において，長さ ∆x の疎水面を越えるために，
∆x2 の時間をピニング現象として費やしていることになる．ここから，∆x が長くなるにつれ濡れが遅く
なることが理解できる．
液滴半径 24.4 nm の結果でも半径 12.2 nm 同様にピニング時間を ∆x2 で整理可能であった．では，液
滴半径 12.2 nm の結果と半径 24.4 nm の結果の間ではどのようにして整理できるだろうか．図 4.5(a) と
図 4.5(b) に液滴半径 12.2 nm の結果と 24.4 nm の結果を重ねたものを示す．ここで，(a) は縦軸を実時間
で，(b) は縦軸を慣性スケーリング時間で整理したものである．３章では，平滑面上での液滴の濡れ速度
が慣性スケーリングで整理可能であることを示した．しかし，Cassie 表面上におけるピニング時間は慣性
スケーリング時間ではなく，実時間で整理できることがわかる．これはピニングに関して慣性の影響が
入っていないことを示唆している．また，小さい液滴にとってはピニングにより消費される時間が割合長
く，逆に大きい液滴にとっては割合短いことになる．
図 4.6 に動的濡れ中の系の温度を示す．横軸は無次元時間，縦軸には温度であり，液滴半径は 12.2 nm
における平滑面上の動的濡れ現象の結果と，Cassie 表面上の動的濡れ現象の結果を比較する．ただし，こ
こでの温度には気相，液相の両方の寄与が入っている．３章で行った平滑面上の濡れシミュレーション
では，系の温度は固体最下層に設置した熱浴により，図から分かるようにほぼ一定値に保たれる．しか
し，本章で行なっている Cassie 表面上の濡れシミュレーションでは，系の温度が激しく上昇することが
わかった．これは，stick-slip の過程で多くの運動エネルギーが熱へと変換されたためではないかと考えら
れる．この温度上昇による液体の物性変化によって濡れを速められたのではないかと推測する．
今までの考察をまとめると以下の二つである．
• ピニングの発生によって濡れは遅くなるが．その影響は液滴が大きいほど比較的小さい．
• 温度上昇に伴う物性変化により濡れが高速になる．
Fig.4.4 ピニング時間．R=12.2 nm．(a) 実時間スケール，(b) 時間を (∆x)2 で除したもの．
第 4 章 Cassie 条件における動的濡れ現象
40
Fig.4.5 異なる液滴半径のピニング時間．(a) 実時間スケール，(b) 慣性時間スケール．
T (K)
90
89
smooth
∆x=r
88
∆x=3r
87
0
∆x=2r
0
0
∆x=4r
0
86
85
84
0
1
2
3
3
1/2
t/(ρR /γ)
Fig.4.6 R=12.2nm における動的濡れ現象中の温度．
前者は濡れを遅く，後者は速くする寄与であり，この関係のトレードオフにより同じ平衡接触角の平滑面
と比べて濡れが高速にもなりうると結論する．
4.3.2
stick-slip
本項では，stick-slip のメカニズムに関する考察を行う．
液滴の固体との接触面から少し離れた高さにおける半径推移を図 4.7 に示す．z=0.34 nm はいわゆる先
行液膜にあたり，本研究では特に利用していないが参考までに付記した．後半において発散しているのは
気体の吸着を誤認識したためであり，濡れ広がりによるものとは関係ない．stick-slip が明らかに見られて
いない ∆ = 2r0 の場合を見ると，液滴の上層部と下層部がほぼ同じ濡れ半径をもったまま濡れが進行し
ていることがわかる．一方で，∆ = 4r0 の場合を見ると，上層部が下層部を追い越している場面が見受け
られ，またそれが stick-slip の挙動と対応していることがわかる．
そこで，この位置関係の定量化のため，局所接触角を求めた．局所接触角は，z=0.68 nm から z=1.70
nm までの気液界面直線にフィッティングすることで得た．図 4.8 に半径推移と局所接触角を同一グラフ
に描いたものを示す．図のわかりやすさのために，無次元時間 1 までのみを示した．フィッティングの
第 4 章 Cassie 条件における動的濡れ現象
41
Fig.4.7 液滴上層部の半径推移．R=12.2nm．(a)L = 2r0 ，(b)L = 4r0 ．
Fig.4.8 局所接触角．R=24.4 nm．(a)L = 3r0 ，(b)L = 4r0 ，(c)L = 6r0 ．
データ点が少ないためにノイズの多い結果となっているが，stick-slip 現象と共に局所接触角が上下してい
ることが見て取れる．また，局所接触角の最大値がある一定の臨界接触角 θc であるように見える．
臨界接触角を目算で求めたものを表 4.1 に示す．∆x が長いほどに臨界接触角が高くなる傾向が見て取
れた．また，∆x と臨界接触角の余弦をグラフ化したものを図 4.9 に示す．
ここで，slip 現象について二通りの解釈が考えられる．一つは，蒸散・拡散的な物理である．すなわち，
液滴内分子が疎水面を飛び越えて親水面上へと蒸散することによって疎水面を飛び越えるという理屈にな
る．表 4.1 において，∆x sin θc は液滴と，疎水面を超えた先の親水面との最短距離であるが，この値が幅
の長い部分においてほぼ同一であり，また，概ねカットオフ距離と一致することから理解できる．
もう一つは疎水面を濡らすことによるエネルギー損失と接触角の非平衡性に基づく余剰応力のエネル
Table 4.1 臨界接触角 θc の ∆x 依存性．
疎水面長さ ∆x
臨界接触角θc
cos θc
sin θc
∆x sin θc
6r0
o
145
-0.82
0.57
3.4r0
5r0
140o
-0.77
0.64
3.2r0
4r0
o
-0.64
0.76
3.0r0
3r0
o
120
-0.50
0.86
2.6r0
2r0
110o
-0.34
0.93
1.8r0
130
第 4 章 Cassie 条件における動的濡れ現象
42
cos θcr
−0.4
−0.6
−0.8
0
2
4
6
Hydrophobic distance [x2.77Å]
2
180
1.5
150
1
120
0.5
0
0
local angle
r/R
Fig.4.9 臨界接触角の余弦と ∆x．
90
0.2
0.4
0.6
3 1/2
t/(ρR /γ)
0.8
60
1
Fig.4.10 R=24.4 nm，∆x = 5.0r0 における半径推移と局所接触角．
ギーバランスからの理解である．ここで，接触角の非平衡性に基づく余剰応力による余剰エネルギーは，
Ep = F I = γ(cos θeq − cos θ) · I
(4.2)
である．ここで，I は適当な長さスケールとなる．一方，疎水面を濡らすことによるエネルギー損失は，
∆E = ∆(γSL − γSG )
(4.3)
である．二式を連立することで，
∆x =
γ
(cos θeq − cos θ) · I
(γSL − γSG )
(4.4)
となり，今回の場合では θeq = 90◦ となるため，上記で見られた cos θ の関係性を説明できる．以上二点
のメカニズム案が考えられるが，ここではどちらの案が正しいのかの同定は難しい．
第 4 章 Cassie 条件における動的濡れ現象
43
そこで，疎水面の濡れ性を少し高めて同様のシミュレーションを行った．図 4.10 に親水面接触角 45◦ ，
疎水面接触角 135◦ ，∆x = 5r0 の結果を示す．この場合でも同様に局所接触角の上昇が見てとれ，目算で
の値は 130◦ である．仮に蒸散メカニズムならば，局所接触角は疎水面幅に依存するために先の 20◦ -180◦
の場合とほぼ同程度になると思われる．あるいは，親水面における吸着力が弱くなったことを考えると
20◦ -180◦ の場合よりも高くなる必要があるのではないかと考えられる．
一方で，余剰応力が slip のメカニズムであるとするならばこの臨界接触角の減少は説明できる．式 4.4
から理解できるが，疎水面上の濡れ性が高まることにより，より少ない余剰エネルギーで slip することが
可能になり，要求される余剰エネルギーが減少するためである．すなわち，ここでは slip のメカニズムと
して余剰エネルギーからの説明がより妥当ではないかと結論付ける．
4.4
まとめ
本章では，Cassie 表面上での動的濡れ挙動を観察した．疎水面と親水面の界面におけるピニング現象を
観測できた．結論として，
• ピニングの発生によって濡れは遅くなるが．その影響は液滴が大きいほど比較的小さい．
• 温度上昇に伴う物性変化により濡れが高速になる．
の二点が得られ，Cassie 表面上での濡れ速度はこのトレードオフで定まるのではないかと考えられる．
また，ピニング現象は Young の余剰応力を基にしたエネルギー論から説明をつけることが可能ではな
いかと考えた．
44
第5章
マルチスケールシミュレーション
5.1
概要
本章では，連続体解析手法の Phase-Field 法を取り扱う．Phase-Field 法は，有限長の界面を与える
Diffuse-Interface モデルの一つであり，移動三相界面問題の計算に適した手法である．Carlson らは，
Phase-Field 法を用いて動的濡れ現象の解析を行い，適切な境界条件を導入した場合に実験とよく一致す
ることを示している [5, 6]．しかし，適切な境界条件の与え方は不明であり，分子スケールの解析が必要
であることを言及している．
そこで，本研究と同一の物性値の液滴の動的濡れ現象を Phase-Field 法で解析し，結果を合わせこむこ
とからの境界条件の導出を試みる．
5.2
支配方程式
支配方程式には，連続の式，Navier-Stokes 方程式，Cahn-Hilliard 方程式，境界条件式の４つが存在す
るが，先に，連続の式を式 5.1 に，Navier-Srokes 方程式を式 5.2 に示す．
∇·u=0
µ
¶
ρ(C) ∂u
+ (u · ∇)u
ρl
∂t
µ
µ
¶
¶
1
µ(C)
Bo ρ(C)
C∇φ
T
= −∇P +
∇·
(∇u + (∇u) ) −
ez −
Re
µl
Ca ρl
Cn · Ca
(5.1)
(5.2)
ここで，C は密度であり液相バルクで C=1，気相バルクで C=-1 をとる．ρl ，µl はそれぞれ液体バルク
における密度と粘性である．ρ(C)，µ(C) は，下付 l を液相，g を気相とし，
¶
µ
1
ρl (C + 1) − ρg (C − 1)
2
µ
¶
1
µ(C) =
µl (C + 1) − µg (C − 1)
2
ρ(C) =
(5.3)
(5.4)
とし，相に応じた物性の違いを表現する．また，右辺第四項は表面張力項であり，φ は後に式 5.10 で示す
化学ポテンシャルである．上述の 4 つの無次元数は，
ε
Cn = ,
L
ρl U L
Re =
,
µl
√
2 2µl U
Ca =
,
3γ
Bo =
(ρl − ρg )gL2
γ
(5.5)
第 5 章 マルチスケールシミュレーション
45
である．U は代表速度，L は代表長さ，D は拡散係数，g は重力加速度，ε は気液の界面長である．それ
ぞれ，Peclet 数 (Pe) は移流と拡散の比を，Cahn 数 (Cn ) は界面長と代表長さの比を，Reynolds 数は慣性
と粘性の比を，Capillary 数 (Ca) は粘性と表面張力の比を，Bond 数 (Bo) は重力と表面張力の比を表す．
次に，相間の方程式を説明する．まず，密度 C を用いて，系の自由エネルギーを以下に定める．
Z µ
F =
Ω
¶
Z
α1
2
α2 Φ(C) +
|∇C| dΩ + (g(C)(γSG − γSL ) + γSL )dΓ
2
Γ
1
(C + 1)2 (C − 1)2
4
g(C) = 0.5 − 0.75C + 0.25C 3
Φ=
(5.6)
(5.7)
(5.8)
ここで Ω は体積積分，Γ は固体の表面積積分を表す．体積積分第一項はバルクポテンシャルであり，Φ
は C = ±1 において極小値を持ち，気相，液相のバルクにおける安定性を保証する．第二項は気液界
面における表面エネルギーである．両項に存在する α1 と α2 は分子当たりの自由エネルギーに関する
√
√
Phase-Field パラメータであり，界面長 ε を用いて，α1 = (3/2 2)γε，α2 = 3γ/2 2ε である．表面積積
分項では，C=1 において 0，C=-1 において 1 に定まる関数 g(C) を用いて固体が濡れることによるエネ
ルギー変化を表現している．
式 5.6 に示した自由エネルギーを C で微分することによって以下の式となる．
Z
Z
0
(α2 Φ (C) − α1 ∇ C)δCdΩ +
δF =
2
Ω
(α1 ∇C · n + g 0 (C)(γSG − γSL )δCdΓ
(5.9)
Γ
ここで，体積積分項において化学ポテンシャル
φ = α2 Φ0 (C) − α1 ∇2 C
(5.10)
を用いて，以下に示す Chan-Hilliard 方程式を得る．
∂C
1 2
1 2 0
+ u · ∇C =
∇ φ=
∇ (Φ (C) − Cn2 ∇2 C)
∂t
Pe
Pe
√
2 2UεL
Pe =
3Dγ
(5.11)
(5.12)
また，固体界面においては式 5.9 の表面積分項から，局所平衡を仮定した時の境界条件式，
0 = −α1 ∇C · n + γ cos θeq g 0 (C)
(5.13)
が得られるが，Carlson ら [5, 6] によれば初期濡れにおいてはこの境界条件式は適切ではなく，局所平衡
を緩和する時間微分項を加え，
µ
εµf
∂C
+ u · ∇u
∂t
¶
= −α1 ∇C · n + γ cos θeq g 0 (C)
(5.14)
と修正する．ただし，µf の物理的意味は不明であり，この導出が目的となる．更に，固液界面における
スリップ条件を課す場合，Qian ら [54, 60] の手法に従い，スリップ速度 us とスリップ係数 β を用いて，
βus = −µ∇u · n + (α1 ∇C · n − γ cos θeq g 0 (C))
を加える．
∂C
∂x
(5.15)
第 5 章 マルチスケールシミュレーション
46
5.3 MD と Phase-Field の比較
分子動力学法によって計算する液滴はナノスケールに限られ，連続体計算との比較においてはスリップ
の影響が無視できない．そこで，まずスリップの影響を分子動力学法から計算する．具体的には， Navier
境界条件に現れるスリップ係数 β の導出を試みる．β は，粘性応力 τ とスリップ速度 us を用いて，
τxz |wall = βus
(5.16)
で表され，バルクにおける粘性応力
τxz = µ
∂ux
∂z
(5.17)
を用いて，バルクの粘性応力が壁面と等しい事を仮定すれば，
β=
µ
ls
(5.18)
の関係性を得ることができる．ここで，ls はスリップ長さである．
本研究では，図 5.1 に示す Couette 流れの系を構築し速度場を求めることからスリップ長さを計算し
た．ここで，y 方向には壁を設けず，周期境界条件から無限大の奥行きをもつチャネルを仮定しているこ
とになる．幅 Lz =10.8 nm のチャネルの下側の固体層は速度 0 m/s，上側の固体層を x 方向に速度 Uw =10
m/s で動かす．下側固体層の最下層，上側固体層の最上層には Langevin 熱浴を取り付けて 85 K に温度制
御する．また，スリップは濡れ性に応じて変化するが，本研究で h 比較を行う接触角 90◦ のポテンシャル
を用いた．
取得した速度分布を図 5.2 に示す．横軸はチャネルの幅方向にあたる z 方向長さ，縦軸は x 方向速度成
分である．本結果では下部，上部でスリップ速度が異なり，下部では 2.2 m/s，上部では 2.8 m/s であっ
た．これは時間アンサンブルの不足が一因として考えられる．スリップ長さ ls は，チャネル幅 Lz ，壁面
速度 Uw (=10 m/s)，スリップ速度 us を用いて，
ls =
us Lz
Uw − 2us
(5.19)
Fig.5.1 Couette 流れの計算系．
第 5 章 マルチスケールシミュレーション
47
10
Vx (m/s)
8
6
4
2
0
0
2
4
6
z (nm)
8
10
Fig.5.2 Couette 流れのチャネル内速度分布．
から求まる．スリップ速度の差を誤差に踏まえ，スリップ長さは ls = 5.4 ± 0.6 nm と求められた．ただ
し，Yang ら [61] はクエット流れのチャネル幅によってスリップ長さが異なり，幅が短いほどスリップ長
さが長くなることを示しており，注意が必要である．そこで，Phase-Field 法ではスリップ長さを 2.5 nm
とし，比較計算を行った．
液滴半径 8 nm の三次元系シミュレーションの比較結果を図 5.3 に示す．いずれも µf = 0 の結果であ
る．Phase-Field 法では，スリップ条件を課すことによって濡れがより速く進行する．しかし，スリップ
条件下でも分子動力学法の結果の方が濡れが速いことがわかる．また，図 5.4 にスリップ係数と µf の影
響を示す．β = 0 においては，式 5.18 から分かるように無限大のスリップを課す状態である．しかしな
がら，このような状態においてさえ分子動力学シミュレーション結果の方が濡れの進行は速い．
この結果の不一致に関しては，以下の原因が考えられる．
• 分子動力学法で導出した物性値が正しくない．
• 連続体で解くには液滴サイズが小さい．
• Phase-Field で蒸発を無視しているために差が生じている．
そこで今後は，より大きな液滴を分子動力学法シミュレーションで構築し比較する，あるいは，分子動力
学法において蒸発の少ない流体を用いて比較すること，の二点が考えられる．
第 5 章 マルチスケールシミュレーション
48
r (nm)
10
5
0
500
PF, no−slip
PF, slip
MD
1000 1500
t (ps)
Fig.5.3 分子動力学法と Phase-Field 法の計算結果の比較．
Fig.5.4 Phase-Field 法におけるスリップの影響．
49
第6章
結論
本研究では分子動力学法を用いて Lennard-Jones 液体の動的濡れ挙動解析を行った．
液滴半径 6.1nm，12.2nm，24,4nm の三種類の大きさの液滴を作成し，濡れ半径の時間発展が従うス
ケーリング則を求めたところ，マクロスケールで見られる慣性スケーリングがナノスケールでも確認でき
た．これは，ナノスケールであっても Re が 1 より大きいことからも理解できる．ここから，慣性が支配
的な濡れの初期においてはナノスケールの液滴であってもサブミリメートルの液滴と同様の機構によって
濡れが進行すると考えられる．
液滴半径は慣性スケーリングで整理できるが，液滴全体の形状には液滴半径に応じた差異があり，濡れ
半径の物理が液滴全体の物理と分離して考えられることが示唆される．このことから，初期濡れの現象に
おいては固液界面ごく近傍の物理が重要であると考えられる．
また，動的濡れ現象に関する平衡接触角の依存性に関しても，マクロスケールでの実験と同様の結果が
見られた．この点からも，マクロスケールでの濡れとナノスケールの分子動力学法シミュレーションが従
う物理現象が同様であると結論付ける．
Cassie 表面上での動的濡れ挙動解析では，疎水面と親水面の界面におけるピニング現象を観測できた．
ピニングの発生によって濡れが遅くなる効果と温度上昇に伴う物性変化により濡れが高速になる効果のト
レードオフの関係があると考えた．
Phase-Field 法による分子動力学法との同条件での計算では，スリップの影響を考慮する必要があるこ
とがわかったが，結果の一致には至らなかった．分子動力学法で取得した物性値の再確認，連続体解析に
おけるサイズ効果，蒸発の影響の考察が必要だと考えられる．
50
謝辞
丸山塩見研究室で研究を始めて二年が過ぎました．本テーマは一見単純な物理現象のように思われたも
ののなかなか定式化できず，物理の蠱惑さと自身の不甲斐なさを痛感した二年間でした．
塩見先生には私の努力不足で多大なるご迷惑を掛け続けてしまったことを深くお詫び申し上げます．そ
れでも手厚くご指導下さった先生の優しさには心から感謝しています．また，このような興味深いテーマ
を下さったことを始め，学会発表，KTH への短期留学と，貴重な経験を積ませて頂きました．二年間で
返せなかった御恩をこの先返すことができるか分かりませんが，少なくとも二年間の知識と経験を錆び付
かせぬよう尽力して備えたいと思います．
丸山先生とは研究の話をした思い出よりも面白いお話をたくさん聞いた思い出が強く残っています．気
さくな先生の下で気持よく生活できるに違いないと期待していたのですが，お陰様で期待に違わず楽しい
研究生活を過ごすことができました．
KTH の Andreas へ．共同研究者としての私は頼りなく苛立ちを感じる事も少なくなかったと思います．
それでも最後の最後まで共同研究を続けて下さったことを感謝致します．これからもインパクトのある成
果を発表し続けてくれると期待しています．次会うときを楽しみにしています．
グローバル COE プログラムには KTH への短期留学の支援をして頂きました．M1 夏の渡航がなけれ
ば Andreas との活発な議論を続けることはできなかったと思います．感謝申し上げます．
私が楽しく二年間を過ごすことができたのは丸山塩見研の朗らかな皆様のお陰です．お世話になりま
した．
では，また．
51
参考文献
[1] Y. Zheng, H. Bai, Z. Huang, X. Tian, F. Zhao, J. Zhai and L. Jiang, “Directional water collection on
wetted spider silk”, Nature, 463, 640, (2010).
[2] A. R. Parker and C. R. Lawrence, “Water capture by a desert beetle”, Nature, 414, 33, (2001).
[3] W. Barthlott and C. Neinhuis, “Purity of the sacred lotus, or escape from contamination in biological
surfaces”, Planta, 202, 1, (1997).
[4] T. S. Wong, S. H. Kang, S. K. Y. Tang, E. J. Smythe, B. D. Hatton, A. Grinthal and J. Aizenberg,
“Bioinspired self-repairing slippery surfaces with pressure-stable omniphobicity”, Nature, 477, 443,
(2011).
[5] A. Carlson, M. Do-Quang and G. Amberg, “Modeling of dynamic wetting far from equilibrium”,
Physics of Fluids, 21, 121701, (2009).
[6] A. Carlson, M. Do-Quang and G. Amberg, “Dissipation in rapid dynamic wetting”, Jounal of Fluid
Mechanics, 682, 213, (2011).
[7] T. Young, “An essay on the cohesion of fluids”, Philosophical Transactions of the Royal Society of
London, 95, 65, (1805).
[8] R. N. Wenzel, “Resistance of solid surfaces to wetting by water”, Industrial and Engineering Chemistry,
28, 988, (1936).
[9] A. B. D. Cassie and S. Baxter, “Wettability of porous surfaces”, Transactions of the Faraday Society,
40, 546, (1944).
[10] Y. J. Sheng, S. Jianga and H. K. Tsao, “Effects of geometrical characteristics of surface roughness on
droplet wetting”, The Journal of Chemical Physics, 127, 234704, (2007).
[11] G. Manukyan, J. M. Oh, D. van den Ende, R. G. J. Lammertink and F. Mugule, “Electrical switching of
wetting states on superhydrophobic surfaces: A route towards reversible Cassie-to-Wenzel transitions”,
Physical Review Letters, 106, 014501, (2011).
[12] B.A. Pethica, “The contact angle equilibrium”, Journal of Colloid and Interface Science, 62, 567,
(1977).
[13] J. H. Weijs, A. Marchand,B. Andreotti, D. Lohse and J. H. Snoeijer, “Origin of line tension for LennardJones nanodroplet”, Physics of Fluids, 23, 022001, (2011).
[14] D. R. Heine, G. S. Grest and E. B. Wevv III, “Surface wetting of liquid nanodroplets: droplet-Size
Effects”, Physical Review Letters, 95, 107801, (2005).
[15] J. F. Joanny and P. G. de Gennes, “Upward creep of a wetting fluid: a scaling analysis”, Journal de
Physique, 47, 121, (1986).
[16] A. Hoang and H. P. Kavehpour, “Dynamics of nanoscale precursor film near a moving contact line of
参考文献
52
spreading drops”, Physical Review Letters, 106, 254501, (2011).
[17] O.V. Voinov, “Hydrodynamics of wetting”, Fluid Dynamics, 11, 714, (1976).
[18] R. L. Hoffman, “A study of the advancing interface I. interface shape in liquid-gas systems”, Journal of
Colloid and Interface Science, 50, 228, (1975).
[19] J. Eggers, J. R. Lister and H. A. Stone, “Coalescence of liquid drops”, Journal of Fluid Mechanics, 401,
293, (1999).
[20] D. G. A. L. Aarts, H. N. W. Lekkerkerker, H. Guo, G. H. Wegdam and D. Bonn, “Hydrodynamics of
droplet coalescence”, Physical Review Letters, 95, 164503, (2005).
[21] J. D. Paulsen, J. C. Burton and S. R. Nagal, “Viscous to inertial crossover in liquid drop coalescence”,
Physical Review Letter, 106, 114501, (2011).
[22] L. H. Tanner, “The spreading of silicone oil drops on horizontal surfaces”, Jornal of Physics D: Applied
Physics, 12, 1473, (1975).
[23] P. G. de Gennes, “Wetting: statics and dynamics”, Review of Modern Physics, 57, 827, (1985).
[24] L. Leger and J. F. Joanny, “Liquid spreading”, Reports on Progress in Physics, 55, 431, (1992).
[25] G. McHale, M. I. Newton, S. M. Rowan and M. Banerjee, “The spreading of small viscous stripes of
oil”, Jornal of Physics D: Applied Physics, 28, 1925, (1995).
[26] A. Biance, C. Clanet and D. Quéré, “First steps in the spreading of a liquid droplet”,Physical Review E,
69, 1, (2004).
[27] L. Duchemin, J. Eggers and C. Josserand, “Inviscid coalescence of drops”, Journal of Fluid Mechanics,
487, 167, (2003).
[28] J. C. Bird, S. Mandre, and H. A. Stone, “Short-time dynamics of partial wetting”, Physical Review
Letter, 100, 234501, (2008).
[29] T.D. Blake, A. Clarke, J. De Coninck, M. de Ruijter and M. Voué, “Droplet spreading: a microscopic
approach”, Colloid and Surfaces A: Physicochemical and Engineering Aspects, 149, 123, (1999).
[30] M. J. de Ruijter, D. De Coninck and G. Oshanin, “Droplet spreading: partial wetting regime revisited”,
Langmuir, 15, 2209, (1999).
[31] M. J. de Ruijter, M. Charlot, M. Voué and J. D. Coninck, “Experimental evidence of several time scales
in drop spreading”, Langmuir, 16, 2363, (2000).
[32] D. R. Heine, G. S. Grest and E. B. Wevv III, “Spreading dynamics of polymer nanodroplets in cylindrical
geometries”, Physical Review E, 70, 011606, (2004).
[33] T. D. Blake, “The physics of moving wetting lines”, Journal of Colloid and Interface Science, 299, 1,
(2006).
[34] E. Bertrand, T. D. Blake and J. De Coninck, “Influence of solid-liquid interactions on dynamic wetting:
a molecular dynamics study”, Journal of Physics: Condensed Matter, 21, 464124, (2009).
[35] D. Seveno, N. Dinter and J. De Coninck, “Wetting dynamics of drop spreading. New evidence for the
microscopic validity of the molecular-kinetic theory”, Langmuir, 26, 14642, (2010).
[36] G. He and N. G. Hadjiconstantinou, “A molecular view of Tanner’s law: molecular dynamics simulations
of droplet spreading”, Journal of Fluid Mechanics, 497, 123, (2003).
[37] C. Wu, T. Qian and P. Sheng, “Droplet spreading driven by van der Walls force: a molecular dynamics
study”, Journal of Physics: Condensed Matter, 22, 325101, (2010).
[38] J. De Coninck, U. D’Ortona, J. Koplik and J. R. Banavar, “Terraced spreading of chain molecules via
参考文献
53
molecular dynamics”, Physical Review Letter, 74, 928, (1995).
[39] V. Samsonov and A. Ratnikov, “Comparative molecular dynamics study of simple and polymer nanodroplet spreading”, Colloids and Surfaces A: Physicochemical and Engineering Aspects, 298, 52,
(2007).
[40] V. Samsonov, N. Zhukova and V. Dronnikov, “Molecular-dynamic simulation of nanosized droplet
spreading over a continual solid surface”, Colloid Journal, 71, 835, (2009).
[41] C. Yang, U. Tartaglino and B. N. J. Persson, “Influence of surface roughness on superhydrophobicity”,
Physical Review Letter, 97, 116103, (2006).
[42] C. Yang and B. N. J. Persson, “Nanodroplets on rough hydrophilic and hydrophobic surface”, The
European Physical Journal E, 25, 139, (2008).
[43] F. Leroy, “Solid-liquid surface free energy of Lennard-Jones liquid on smooth and rough surfaces computed by molecular dynamics using the phantom-wall method”, The Journal of chemical physics, 133,
044110, (2010).
[44] T. Koishi, K. Yasuoka, S. Fujikawa and X. C. Zeng, “Measurement of contact-angle hysteresis for
droplets on nanopillared surface and in the cassie and wenzel states: A molecular dynamics simulation
study”, ACS Nano, 5, 6834, (2011).
[45] C. D. Daub, J. Wang, S. Kudesia, D. Bratko and A. Luzar, “The influence of molecular-scale roughness
on the surface spreading of an aqueous nanodrop”, Faraday Discussions, 146, 67, (2010).
[46] S. Adelman and J. Doll, “Generalized Langevin equation approach for atom/solid- surface scattering:
General formulation for classical scattering off harmonic solids”, The Journal of Chemical Physics, 64,
2375, (1976).
[47] S. Adelman, “Generalized langevin theory for many-body problems in chemical dynamics: General
formulation and the equivalent harmonic chain representation”, The Journal of Chemical Physics, 71,
4471, (1979).
[48] J. Tully, “Dynamics of gas-surface interactions: 3D generalized Langevin model applied to fcc and bcc
surfaces”, The Journal of Chemical Physics, 73, 1975, (1980).
[49] S. Maruyama, “Molecular dynamics method for microscale heat transfer”, Advances in Numerical Heat
Transfer, 226, 2189, (2000).
[50] http://www.photon.t.u-tokyo.ac.jp/m̃aruyama/nanoheat/phantom.pdf
[51] C. B. Moon, G. K. Moon and S. H. Lee, “Equilibrium and non-equilibrium molecular dynamics simulations of thermal transport coefficients of liquid argon”, Bulletin of the Korean Chemical Society, 12,
309, (1991).
[52] J. G. Kirkwood and F. P. Buff, “The statistical mechanical theory of surface tension”, The Jounal of
Chemical Physics, 17, 338, (1949).
[53] J. H. Irving and J. G. Kirkwood, ”The statistical mechanical theory of transport processes. IV. The
equations of hydrodynamics*”, The Jounal of Chemical Physics, 18, 817, (1950).
[54] T. Z. Qian, X. P. Wang and P. Sheng, “Molecular scale contact line hydrodynamics of immiscible flows”,
Physical Review E, 68, 016306, (2003).
[55] V. G. Baidakov, S. P. Protsenko, Z. R. Kozlova and G. G. Chernykh, “Metastable extension of the liquidvapor phase equilibrium curve and surface tension”, The Journal of Chemical Physics, 126, 214505,
(2007).
参考文献
54
[56] I. A. Cosden and J. R. Lukes, “Effect of cutoff radius on the surface tension of nanoscale bubbles”,
Journal of Heat Transfer, 133, 101501, (2011).
[57] K. Meier, A. Laesecke and S. Kabelac, “Transport coefficients of the Lennard-Jones model fluid. I.
Viscosity”, Journal of Chemical Physics, 121, 3671, (2004).
[58] L. Bai and D. Breen, “Calculating center of mass in an unbonded 2D environment”, Journal of Graphics,
GPU, and Game Tools, 13, 53, (2008).
[59] D. Jacqmin, “Contact-line dynamics of a diffuse fluid interface”, Journal of Fluid Mechanics, 402, 57,
(2000).
[60] T. Z. Qian, X. P. Wang and P. Sheng, “Power-law slip profile of the moving contact line in two-phase
immiscible flows”, Physical Review Letters, 93, 9, (2004).
[61] X. Yang and Z. Zhang, “Effect of channel scale on slip length of flow in micro/nanochennels”, Journal
of Fluids Engineering, 132, 061201, (2010).