日本機械学会論文集（B 編） 原著論文 No.2013-JBR-0102
長山 暁子
*1，竹松 雅樹
*2，鶴田 隆治
*3Condensation/ Evaporation Coefficient of Chain Molecules Gyoko NAGAYAMA
*1, Masaki TAKEMATSU and Takaharu TSURUTA
*1 Kyushu Institute of Technology, Dept. of Mechanical Engineering Sennsui 1-1, Tobata, Kitakyushu, Fukuoka, 804-8550 Japan
Molecular dynamics (MD) simulations are performed to study the liquid-vapor interface structure of a group of carbon chain molecules: butane (C4H10), octane (C8H18) and dodecane (C12H26) in this paper. The effects of the chain length and temperature on the liquid-vapor interface structure, molecular evaporation/ condensation behavior and the orientation of the liquid-vapor interfacial molecules in equilibrium systems are investigated. It is found that the condensation/ evaporation coefficient of carbon chain molecules primarily depends on the translational energy and the surface temperature similar to simple molecules like argon and water. The MD data of carbon chain molecules agree well with the theoretical expression based on the transition state theory. Also, it is found that the chain ordering at the interface depends on the molecular structure but there’s less effect on the condensation/ evaporation behavior at liquid- vapor interface. We conclude that the condensation/ evaporation coefficient can be predicted by the translational length ratio of liquid to vapor in general even for the chain molecules.
Key Words : Evaporation/Condensation Coefficient, Liquid-Vapor Interface, Molecular Dynamics, Transition State Theory, Translational Length Ratio
1. 緒 言
マイクロ・ナノスケールにおける気液界面相変化現象は，燃料液滴の噴霧燃焼，インクジェットプリンターの インクの溶媒蒸発，噴霧冷却，気体・液体・固体三相界面近傍のミクロ領域の伝熱促進など幅広い工学分野にお いて重要である(1-5)．例えば，ディーゼルエンジンなどで使われる噴霧燃焼では，亜臨界圧力において燃料蒸気の 拡散燃焼時間は液滴の蒸発寿命に比べ非常に小さく，液滴の蒸発が支配的になる(1, 2)．インクジェット印刷技術に 用いるインク液滴は，高画質や超微細パタン描画のため，1pl（直径約12μm）以下になった(3, 4)．このような液滴 のサイズが非常に小さい場合，液滴の蒸発過程をモデリングするには流体力学モデルでは適用限界があり，気体 分子運動論によるモデルを用いることが望ましいとする報告がある(6, 7)．
気体分子運動論に基づく気液界面での物質伝達を記述するためには，境界条件として，凝縮・蒸発係数を正確 に見積もることが必要である．気液界面に蓄積した不凝縮性ガスの排除および凝縮あるいは蒸発を伴う強い非平 衡条件の必要性から(8)，実験的に凝縮・蒸発係数を精度よく見積もることが難しく，研究現場をはじめ多くの場 合では凝縮・蒸発係数を経験的に与えるのが現状である．一方，分子動力学シミュレーション（Molecular Dynamics Simulation：MD）による検討によって，従来のマクロパラメータとしての凝縮・蒸発係数から，ミクロスケール の界面相互作用を考慮し得る凝縮・蒸発確率としてのミクロ凝縮・蒸発係数の必要性を明らかにしている(9-14)．分 子の凝縮・蒸発確率が分子の持つ並進運動エネルギーの界面鉛直成分と界面温度に依存することを，エネルギー
鎖状分子の凝縮・蒸発係数
＊* 原稿受付 2013年2月2日 改訂原稿受付日 2013年7月12日
*1 正員，九州工業大学工学研究院機械知能工学研究系
（〒804-8550 福岡県北九州市戸畑区仙水町1-1）
*2 学生員，九州工業大学大学院工学府機械知能工学専攻
*3 正員，フェロー，九州工業大学工学研究院機械知能工学研究系 E-mail: nagayama@mech.kyutech.ac.jp
79 巻 806 号 （2013- 10）
の大きい分子が界面のエネルギー障壁を越えやすいとの解釈を与える遷移状態説理論を用いれば合理的に説明で きる．気液相変化の凝縮・蒸発素過程において分子の並進運動の束縛（並進運動の自由度の制限）が最も重要で あることから，平衡状態におけるあらゆる物質の凝縮・蒸発係数は気相と液相の自由体積の比を用いて提案され た(10)．
( ) ⎟ ⎟
⎠
⎞
⎜ ⎜
⎝
⎛
− −
−
=
= 1
1 2 exp 1
1
33
l g g
l e
c
σ v v v v
σ
(1)ここで，vl，vgはそれぞれ液相と気相の比体積を表し，3 vl vg を並進運動長さスケール比とする．式（1）によれ ば，凝縮・蒸発係数の理論値は気液界面の温度に依存し，三重点付近では1に近づき，臨界点に向けて徐々に減 少する．これまでに，分子構造が比較的に簡単であるアルゴン，水およびメタノールについては，分子動力学シ ミュレーションによる解析結果と理論値とのよい一致を確認し(9, 11, 15)，より複雑な分子構造を持つ物質に対する 検討も近年行われてきた(16-19)．38 個の原子からなるドデカン（C12H26）分子については，著者らの解析結果が上 記理論値とのよい一致を示したが(16)，Caoらは理論値と異なる結果を報告するものの(17)，同じ研究グループのXie らは理論値と一致する結果を示した(18)．また，Chengらは単原子（monomer），二量体（dimer）と三量体（trimer）
からなるそれぞれのLennard-Jones流体に対する分子動力学解析を行い，式（1）の理論値よりやや大きい結果を 得て，分子間相互作用がより強い炭化水素系鎖状分子に対する検討の必要性を指摘した(19)．しかし，これらの結 果のいずれも式（1）の理論値と概ね同じ傾向を示しており，計算による差が生じたものの，理論式の妥当性を検 証できたという見方もできる．
内部自由度が大きい分子構造を持つ物質の凝縮・蒸発係数が理論式に従わない可能性があるとしたら，相変化 過程の遷移状態における分子運動エネルギーの制限度合いが，分子構造あるいは気液界面における分子配向に依 存する場合にあると推察される．そこで，本報では，鎖状分子の鎖長や分子配向等が気液界面における凝縮・蒸 発挙動に及ぼす影響を解明することを目的として，異なる鎖長を持つ炭化水素系鎖状分子を対象とした分子動力 学シミュレーションを行った．前報で報告したドデカンの結果を含めて，ブタン (C4H10), オクタン (C8H18)を新 たに分子動力学シミュレーション対象に追加した．炭化水素系鎖状分子の分子構造や分子配向が平衡系における 気液界面の凝縮・蒸発特性に及ぼす影響を調べ，鎖状分子に対する上記理論式の妥当性を検証することにした．
主な記号
A ：L-Jポテンシャルパラメータ [J/(mol･m12)]
B ：L-Jポテンシャルパラメータ [J/(mol･m6)]
Etot ：ポテンシャルエネルギーの総和 [J]
Et,z ：並進エネルギーの界面鉛直成分 [J]
E0 ：活性化エネルギー [J]
kb ：ボルツマン定数 [J/K]
Kr ：共有結合の力の定数 [J/(mol･m2)]
Kθ ：結合角の力の定数 [J/(mol･rad2)]
Lx, Ly, Lz ：計算系のサイズ [m]
m ：分子一個の質量 [kg]
n ：周期 [-]
N ：分子数 [-]
q ：電荷 [C]
r ：結合距離 [m]
R ：原子間距離 [m]
S ：配向秩序パラメータ [-]
T ：温度 [K]
Tc ：臨界点温度 [K]
Vn ：ねじれ角の力の定数 [J/mol]
Vz ：蒸気分子の界面鉛直方向速度成分 [m/s]
α, β ：ミクロ凝縮・蒸発係数のパラメータ [-]
δ ：界面層の厚さ [nm]
ε ：誘電率 [F/m]
φ ：ねじれ角 [rad]
ϕ ：主鎖と界面平行方向とのなす角 [rad]
γ ：位相 [rad]
v ：比体積 [m3/kg]
θ ：結合角 [rad]
ρ ：密度 [kg/m3]
ρc ：臨界点密度 [kg/m3]
σ ：ミクロ凝縮・蒸発係数 [-]
σ ：凝縮・蒸発係数 [-]
添え字
c ：凝縮
e ：蒸発
eq ：平衡状態
g ：気体
i, j ：iおよびj番目の原子
l ：液体
2. 分子動力学シミュレーション
鎖長の異なるブタン（C4H10）, オクタン（C8H18）とドデカン（C12H26）を分子動力学計算の対象分子に選び，
その分子構造イメージを図1に示す．両サイドのメチル基（CH3）の間に，ブタンは2個，オクタンは6個，ド デカンは 10 個のメチレン基（CH2）で直鎖状につながり，直線上に並べた際の分子長さがそれぞれ 0.381nm, 0.889nm, 1.390nmとなる．ブタン，オクタンおよびドデカンの諸物性値はChemistry WebBook(20)より入手し，その 一部を表1に示す．分子動力学シミュレーションは，九州大学情報基盤研究センター高性能アプリケーションサ ーバシステム（日立SR16000 L2）および高性能演算サーバ（富士通 PRIMERGY CX400 S1，H24年10月より運
Butane : CH3-(CH2)2-CH3 Octane : CH3-(CH2)6-CH3 Dodecane : CH3-(CH2)10-CH3 C atom H atom
Fig. 1 Molecular structures of the carbon chain molecules.
~0.381nm ~0.889nm ~1.390nm
Table 1 Simulation systems and parameters.
Butane Octane Dodecane
Molecular weight 58.12 114.23 170.33
Triple point [K] 134.60 216.20 263.16
Critical point [K] 425.10 569.30 658.10
Temperature [K] 220~340 330~460 300~450 500~550 Molecules (atoms) 560 (7840) 294 (7644) 208 (7904) 1020 (38760) Lx[nm]×Ly[nm]×Lz [nm] 4×4×24 4×4×24 4×4×24 7×7×30
用）を用いて，パッケージプログラムAMBERバージョン10およびバージョン12を用いて行った(21, 22)．力場に はGAFF（General Amber Force Field）を適用し，原子間ポテンシャルエネルギーは，分子内相互作用を表す共有 結合，結合角，ねじれ角と分子間相互作用を表す短距離力，クーロン力の総和となり，次式に表す(22)．
( )
∑ −
=bonds r eq
tot K r r
E 2 + ∑
(
−)
2anglesKθ θ θeq + ∑
[
+(
−) ]
torsions
n n
V 1 cos φ γ
2 ∑ ⎟⎟⎠
⎞
⎜⎜
⎝
⎛ −
+
<
atoms j
i ij
ij ij
ij
R B R
A
6
12 + ∑
<
atoms j
i ij
j i
R q q
ε (2)
ここで，req, θ eq, Kr, Kθ, Vn, γ, n, Aij, Bij, qi, qj, ε はGAFFによって与えられるパラメータであり，クーロン力の計算に はPME（Partical Mesh Ewald）を用いた．計算系の詳細は表1に示す．比較のため，計算条件は前報(16)と同様に 設定し，系の分子数N，体積VおよびエネルギーE一定のNVEアンサンブルを用いた．計算系の温度は，それ ぞれの対象分子の三重点と臨界点の間に設定し，温度が比較的に低い場合には，4 nm×4 nm×24nm の計算系を 用い，臨界点に近い場合には，蒸気空間を確保するため，計算系のサイズを7 nm×7 nm×30nmまでに拡大した．
時間ステップは0.5fs，カットオフ距離は1.5nmとし，全ての境界面に周期境界条件を適用した．初期状態から約 300ps の間に温度制御を行い，その後温度制御をはずして目的温度に応じた平衡液膜を作成した．緩和計算より 十分平衡化した後の液膜のスナップショットを図2に示す．
気液平衡状態において，液相に向かう気相の凝縮質量流束と気相に向かう液相の蒸発質量流束が釣り合うこと より，凝縮係数も蒸発係数も同じ値になる．実際の計算では，系の温度が低い場合では，蒸気空間中に自己蒸発 する分子が少なく，蒸発係数を直接求めるのは精度が悪い．十分なサンプル数を得るためにテスト分子による投 入シミュレーションを行い，凝縮係数を求めることにした．気相空間に液膜中心からz方向Z=
±
4.5~5.5 nmの位 置にテスト分子をランダムに配置し，液膜に向かう異なる並進エネルギーを与えた．投入シミュレーションのサ ンプリング期間100psの間に，テスト分子のうち，実際に液膜に凝縮した分子数より凝縮確率を求め，凝縮係数 とした．系の温度が比較的に高い場合では，蒸気空間に十分なサンプル分子が存在するため，気液界面近傍に検 査面を設け，液相側に飛び込む分子が検査面を通過するときの時間，速度および個数を計測した．サンプリング 期間の間に，液相側に飛び込む分子のうち，定めた時間内に気相側に戻らない分子数より凝縮係数を求めた．こ れまでの計算方法に従い，検査面は気相の密度が飽和蒸気密度より飽和液体密度の 1％ほど高い，すなわち，(ρg+0.01ρl)になる位置に設定し，反射分子と凝縮分子の判別時間を20psに設定した．なお，検査面の位置および 反射分子のサンプリングが重要で，わずかでありながら，凝縮・蒸発係数の値に計算誤差が生じてしまう．
3. 計算結果および考察 3・1 気相と液相の密度
系が十分平衡化した後に，z方向長さを0.2nmの間隔で分割し，2.5ns～8nsのサンプリング期間における分子 の重心（質量中心）から求めた局所密度の時間平均より密度分布を得た．図3に各温度条件下における系の密度
x z y
Fig. 2 Sample snapshots of simulation systems in equilibrium state of butane, octane and dodecane.
220K 300K 340K Butane
330K 430K 460K Octane
300K 500K 550K Dodecane
分布を示すが，典型的な双曲線正接関数の形状になっていることが確認できる(23)．これより得たバルクの気相と 液相の密度の時間平均を表2に整理し，Chemistry WebBook(20) による気液共存曲線との比較を図4に示す．表2 と図4より，解析結果は物性値と概ね一致するが，液相密度は物性値よりやや低く，高温領域の系の気相密度は
butane
octane
dodecane
0 100 200 300 400 500 600 700
0 0.2 0.4 0.6 0.8 1
300K350K
450K 400K 500K 550K
z/Lz Density [kg/m3 ]
Fig. 3 Density profiles of butane, octane and dodecane in equilibrium systems at different temperatures.
0 100 200 300 400 500 600
700 220K
270K 300K 320K 340K Density [kg/m3 ]
0 100 200 300 400 500 600
700 330K
400K 430K450K 460K Density [kg/m3 ]
Fig. 4 Liquid-vapor coexistence curves.
Fig. 5 Comparison of the liquid-vapor interface thickness between δ90−10 and δ95−1. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0
0.2 0.4 0.6 0.8 1
Butane Octane Dodecane Butane
●
Octane
■
Dodecane
▲
Cao et al., 2011(17)
△
Xie et al., 2011(18)
▽
MD data Property (20)
ρ / ρc [-]
T / T c [-]
0.0 1.0 2.0 3.0 4.0 5.0
0 0.2 0.4 0.6 0.8 1
δ95-1 , δ90-10 [nm]
Translational Length Ratio (v l /v g )1/3 Butane
●
Octane
■
Dodecane
▲
Butane
○
Octane
□
Dodecane
△
δ95-1 δ90-10 Table 2 Density of bulk liquid and vapor and the thickness of liquid-vapor interface.
T [K] ρl [kg/m3] ρg [kg/m3] δ90−10 [nm] δ95−1 [nm]
MD Property(20) MD Property(20) MD MD
butane 220 634.78 655.72 0.29 0.25 0.70 1.05
270 573.75 604.19 6.37 2.48 1.00 1.55
300 533.64 570.75 12.24 6.55 1.30 2.15
320 496.61 546.78 23.96 11.33 1.70 2.45
340 459.72 519.74 44.55 18.58 2.90 3.00
octane 330 648.24 666.92 1.17 0.50 0.70 1.20
400 569.69 614.45 9.58 3.47 1.25 1.80
430 529.47 575.12 20.30 9.14 1.60 2.25
450 499.25 566.19 23.83 12.57 2.10 3.05
460 474.77 540.24 36.24 16.40 2.45 3.10
dodecane 300 719.83 744.29 0.00 0.00 0.60 0.96
350 686.98 707.26 0.00 0.04 0.70 1.10
400 644.18 669.29 0.00 0.34 0.85 1.46
450 597.58 628.88 3.66 1.68 1.20 1.87
500 544.32 584.01 10.21 5.80 1.90 3.06
550 470.12 531.26 30.48 16.05 3.25 4.86
やや高い結果となる．これは，AMBERに用いた原子間ポテンシャルのGAFF力場による結果だと考えられ，高 温領域の結果の精度を向上するためには，より厳密なポテンシャルモデルを用いることが望ましい．一方，ドデ カンの液相密度について，Caoらの結果は全体的に物性値より低く(17)，Xieらの結果が物性値に一致している(18)． 気液界面において，液相の飽和密度から気相の飽和密度に向かって遷移する領域は気液界面層とされ，一般的 に液体の飽和密度の90％－10％となる厚みが90－10界面層厚み，δ90−10と定義される(23)．表2にδ90−10の結果を整 理しており，界面層の厚さがおおよそ 1-5nm の間に変化し，系の温度増加とともに増加することが確認できる．
図5に界面層の厚さと並進運動長さスケール比 (vl /vg)1/3との関連を整理しており，分子構造や鎖長さに依存せず，
並進運動長さスケール比が大きいほど界面層の厚みが増すことを確認できる．しかし，気体分子運動論を適用す る際の界面位置としては，気相側の飽和密度に準ずることが多いため，界面の分子の揺らぎを考慮して ρg+0.01ρl となる位置で凝縮・蒸発・反射の境界条件を決めるのがより妥当だと考える．このため，凝縮と蒸発と反射の分 子挙動を調べる際，液相密度の95％（0.95ρl）から (ρg +0.01ρl)になるまでの遷移領域を95－1界面層と定義し，
ρg+0.01ρlとなる位置に検査面を設置した．図5に示すように，δ90−10とδ95−1が概ね同じ傾向を示しており，当然な がら，後者の方が前者よりやや大きくなっている．
3・2 気液界面における鎖状分子の凝縮・蒸発特性
界面分子との相互作用にかかわる時間から，界面層滞在時間により蒸発と反射とを分類している．凝縮過程を 例に挙げると，気相から一度検査面の液側に飛び込み，判別時間以内に再び蒸気側に戻された分子を反射，蒸気 側に戻されるのに判別時間以上にかかった分子を蒸発分子と定義する．判別時間が短いほど，入射分子が界面層 から受ける影響が少なく，液相からの真の蒸発分子との区別がしやすいが，判別時間が長いほど入射分子が界面 層に適応し，蒸発分子との区別ができなくなってしまう．図6にサンプリング期間内に検査面を出入りにかかる 時間と分子数の結果例を示す．サンプリング数が異なるが，一度液相に飛び込み再び気相側に戻された分子の平 均滞在時間は約40psで，全体の7割が20ps以内に再び気相に戻っている状況は同じである．検査面を通過して 再び気相に戻る際の速度分布を図7に滞在時間別に示す．滞在時間が20ps以内であった分子はMaxwell 速度分 布より遅い側にずれ，40ps以上60ps以下の分子は，液体側のMaxwell速度分布に完全適応した状態で検査面を 再び通過したことが分かる．また，滞在時間が60ps以上であった分子は，Maxwell速度分布より速度の大きい側 にシフトし，液体内部からの蒸発分子の速度分布と同じ特徴を示す．これらの反射特性はブタン，オクタンとド デカンの分子構造にほぼ依存せず，図 6 と図 7 に示す例と同じ傾向になるため，本報も前報同様に判別時間を 20psとする(16)．
Fig. 6 Comparison of the staying period in liquid side of reflected incident molecules.
0 500 1000
Octane 460K during 7350ps 0
1000 2000 3000
4000 Butane 340K during 6075ps
0 500 1000 1500
0 40 80 120 160 200 240
Dodecane 500K during 2500ps
Staying period [ps]
0 0.001 0.002 0.003 0.004 0.005 0.006
0 100 200 300 400 500 600
Maxwellian
0~20 [ps]
20~40 [ps]
40~60 [ps]
60~ [ps]
Fz
Vz [m/s]
Staying period in liquid side
Fig. 7 Comparison of the velocity distribution functions of reflected molecules within different staying period in liquid side (dodecane at 500K).
異なる温度条件で得た鎖状分子の界面鉛直方向の並進エネルギーを考慮したミクロ凝縮・蒸発係数の結果を図 8 に示す．並進運動エネルギーの大きい領域および低温界面においてサンプル分子数が比較的少ないため，一部 のデータのばらつきが見られるが，分子構造にほぼ依存せず，いずれの鎖状分子についても，界面鉛直方向の並 進運動エネルギーの小さい分子ほど凝縮・蒸発しにくく，それの大きい分子ほど凝縮・蒸発しやすいことが分か る．また，温度が高いほど凝縮・蒸発係数は低下する傾向を示す．すなわち，鎖状分子の凝縮・蒸発特性はこれ まで確認したアルゴン分子や水分子などと一致し，次式により表現できる(9)．
⎭⎬
⎫
⎩⎨
⎧ ⎟⎟⎠
⎜⎜ ⎞
⎝
⎛−
−
⋅
=
= k T
mV
b z e
c 1 exp 2
β 2
α σ
σ (3)
平衡系では蒸気分子はMaxwell分布に従うことから，従来からの定義である平均値としての凝縮・蒸発係数は,
⎟⎠
⎜ ⎞
⎝⎛ −
⋅
=
=σ α 1 β2
σc e (4)
Table 3 Parameters of condensation/ evaporation coefficient in Eq. (4) and activation energy.
Temperature [K]
Condensation/ Evaporation Coefficient Activation Energy
α β σ__c=σ__e E0 / kbT [-]
butane 270 0.985 0.301 0.837 0.089
300 0.903 0.361 0.740 0.110
320 0.821 0.475 0.626 0.156
340 0.746 0.612 0.518 0.220
octane 400 0.981 0.321 0.824 0.096
430 0.923 0.385 0.745 0.119
450 0.873 0.448 0.677 0.144
460 0.839 0.499 0.630 0.166
dodecane 300 1.000 0.020 0.990 0.005
350 0.995 0.017 0.986 0.004
400 0.993 0.051 0.968 0.013
450 0.977 0.099 0.929 0.026
500 0.901 0.269 0.780 0.078
550 0.805 0.489 0.608 0.162
Fig. 8 Microscopic condensation/ evaporation coefficients as the function of the normal component of molecular translation energy and temperature: (a) butane; (b) octane; (c) dodecane.
(a) butane (b) octane (c) dodecane
0 400 800 1200 1600
Eq. (3)
Et,z / kb [K]
T=400K
●
T=430K
■
T=450K
▲
T=460K
◆
0 400 800 1200 1600 2000
Et,z / kb[K]
Eq. (3) T=300K
○
T=350K
□
T=400K
△
T=450K
◇
T=500K
●
T=550K
■
0 0.2 0.4 0.6 0.8 1
0 200 400 600 800
Et,z / kb [K]
σe , σc
Eq. (3) T=270K
●
T=300K
■
T=320K
▲
T=340K
◆
となる．表 3 にα, β ，平均値としての凝縮・蒸発係数の結果を示す．α は気液界面の温度上昇とともに低下し，
分子が高温液面では凝縮・蒸発しにくくなることを意味する．また，β は分子の運動状態が凝縮・蒸発係数に及 ぼす影響の程度を示し，β が大きいほど，凝縮・蒸発係数の並進エネルギー依存性が強く現れ，並進運動エネル ギーが大きい分子しか凝縮・蒸発できないことになる(9), (10)．つまり，界面層における気液相変化に必要な活性化 エネルギーがβを用いて次式のように表現でき(10)，
(
2)
12
0 4 −
−
= − β β T
E k
b
(5)
界面層を通過する際に分子が受ける抵抗を意味する．表3にその結果を示す．鎖状分子についても，活性化エネ ルギーは界面温度上昇とともに増加し，この界面抵抗を乗り越えるために並進エネルギーへの依存性が強くなる．
図9に式（1）と従来の結果に，鎖状分子の凝縮・蒸発係数を加えて比較する．単原子分子，多原子分子，鎖状 分子のいずれも同じ傾向を示し，分子動力学解析結果と理論値との一致が確認できる．解析結果と理論値との差 については，本研究で用いたAMBERのポテンシャルモデルによるものだと考えられる．これは，解析した液相 の密度が物性値を厳密に表現できておらず，物性値より低くなったため，液相の比体積vlおよび並進運動長さス ケール比(vl /vg)1/3を過大評価してしまったためである．その結果，ブタン，オクタンとドデカンの分子構造によら ず，解析結果は全体的に式（1）よりも右側にずれることとなった．(vl /vg)1/3の解析結果の代わりに物性値を用い て評価すると，図10(a)に示すように，理論値との差が小さくなることを確認できる．なお，Guらは，液相と気 相の密度をより厳密に表現すれば，解析結果が式（1）とよく一致することを報告した(24)．
一方，分子動力学解析は分子数，計算系サイズ，サンプリング時間の制限を受け，揺らぎによる計算誤差が生 じることも原因であることを否めない．図10(b), (c)に凝縮・蒸発と反射分子を区別するために必要な検査面の位 置および判別時間が解析結果に及ぼす影響について調べた結果を示す．検査面の位置および判別時間のわずかな 差で，約5~20%程度の誤差が生じることを確認できる．比較のため，図9にCaoら(17)，Xieら(18)の計算結果も示 している．これらのデータのずれはOPLSポテンシャルモデルによる密度の計算誤差，検査面のずれおよび判別
Fig. 9 Mean value of condensation/ evaporation coefficients as a function of translational length ratio.
0.0 0.2 0.4 0.6 0.8 1.0
0 0.2 0.4 0.6 0.8 1
σe , σc
Translational Length Ratio (v l /v g )1/3 MD data Eq.(1)
Butane
●
Octane
■
Dodecane
▲
Dodecane(17)
△
Dodecane(18)
▽
Methanol(15)
○
Argon(9)
□
Water(11)
◇
Fig. 10 Errors range of condensation/ evaporation coefficients caused by inaccuracy of (a) density, (b) liquid-vapor interface position, and (c) distinct time.
0 0.2 0.4 0.6 0.8 1
Translational Length Ratio (v l /v g )1/3 Butane(340K)
○
Octane(460K)
□
Dodecane(500K)
△
10ps 20ps 40ps 60ps 100ps
0 0.2 0.4 0.6 0.8 1.0
0 0.2 0.4 0.6 0.8 1
Property (20) Eq.(1)
Butane
●
Octane
■
Dodecane
▲
Butane
○
Octane
□
Dodecane
△
MD (vl /vg )1/3
(vl /vg )1/3 σe , σc
Translational Length Ratio (v l /v g )1/3
0 0.2 0.4 0.6 0.8 1
Translational Length Ratio (v l /v g )1/3 Butane
●
Octane
■
Dodecane
▲
Butane
○
Octane
□
Dodecane
△
δ95-1 δ90-10
(a) (b) (c)
時間の曖昧さによる誤差が主な原因だと考えられる．しかし，これらの誤差を考慮しても，Caoら，Xieらの結果 は本報とおおよそ同じ傾向を示し，式（1）と一致していることが言える．以上より，鎖状分子についても，アル ゴンや水と同様に，式（1）が適用可能と判断できる．
3・3 鎖状分子の分子配向と凝縮・蒸発特性
気液相変化の凝縮・蒸発素過程において，気液界面層（遷移状態）にある分子の並進運動の束縛（並進運動の 自由度の制限）度合を知ることが重要である(10)．式（1）では，界面分子の界面鉛直方向の並進運動の長さスケー ルが，気体分子の長さスケールより液体の分が制限されることを表現している．並進運動自由度の等方性からこ れらの長さスケールを3vg ，3vl と仮定し，その結果として凝縮・蒸発係数が (vl /vg)1/3の関数で表れるようになっ た．アルゴンや水のようなほぼ球形で近似できる分子構造の場合，気液界面層の厚みは分子サイズの約3~10倍程 度となり，等方性仮定が成立するとみなせる．しかし，直鎖状分子は細長い形状を持ち，界面層厚みは鎖長さと ほぼ同じ程度になる場合がある．もし界面分子が規則性を持って界面に平行または垂直な配向になれば，界面層 における分子の並進運動自由度が等方的に扱えなくなる．そのため，界面層における分子配向が凝縮・蒸発特性 に及ぼす影響があるかどうかを調べる必要がある．
前報では，鎖状分子であるドデカンの分子配向が凝縮・蒸発特性にほぼ影響しない結果を報告した．具体的に，
低温液膜の界面層の液体側において規則的な分子配向を確認できたにも関わらず，低温において元々凝縮・蒸発 係数が1 に近い値となるため，分子配向が凝縮・蒸発特性に与える影響が小さく，一方で高温液膜の界面層では 液体側の分子配向の特異性が現れず，配向状態が凝縮・蒸発特性にほぼ影響しないという結論となった(16)．しか し，前報の結果は系全体の密度および配向の空間分布によって判断したもので，検査面を通過した際の分子の配 向状態などが不明のため，必ずしも十分な結果とは言えない．そこで，本報では，前節で述べた界面近傍の検査 面を通過した分子に着目し，その分子配向と凝縮・蒸発係数との関連性を調べることにした．気液界面近傍にお ける鎖状分子の分子配向を前報と同じく配向秩序パラメータSで評価した(25)．
1 sin 2 3
1 2 −
= ϕ
S (6)
Fig. 11 Effects of the orientational order parameter on the condensation/ evaporation coefficients.
(a) butane (b) octane (c) dodecane
Fig. 12 Probability density distributions of the orientational order parameter.
0 0.1 0.2 0.3 0.4
-0.5 0 0.5 1
Probability density
○
T=270K ●
□
T=300K ■
△
T=320K ▲
◇
T=340K ◆
Ne(or Nc) Nall
S 0.0
0.2 0.4 0.6 0.8 1.0
-0.5 0 0.5 1
S σe , σc
T=270K
●
T=300K
■
T=320K
▲
T=340K
◆
-0.5 0 0.5 1
S
T=400K
●
T=430K
■
T=450K
▲
T=460K
◆
-0.5 0 0.5 1
S
T=500K
●
T=550K
■
-0.5 0 0.5 1
○
T=400K ●
□
T=430K ■
△
T=450K ▲
◇
T=460K ◆
Ne(or Nc) Nall
S -0.5 0 0.5 1
○
T=500K ●
□
T=550K ■
Ne(or Nc) Nall
S
(a) butane (b) octane (c) dodecane
ここで， ϕ は鎖状分子の主鎖となる共有結合しているC原子を通るベクトルと気液界面に平行する面とのなす角 を表し，Sは分子内のC-C共有結合から集計して共有結合数で平均し，分子単位の配向秩序パラメータとして求 めた．配向秩序パラメータは［－0.5, 1］の範囲にあり，S=－0.5であれば主鎖が気液界面に平行， S=1であれば 垂直な配向になる．
界面近傍の検査面における通過分子の配向秩序パラメータSと凝縮・蒸発係数の結果を図11に示し，いずれも 蒸気空間に十分なサンプル分子が存在し，界面温度が比較的高い計算系に対するものである．一部のデータのば らつきは見られるが，全体的に凝縮・蒸発係数が分子配向に依存せずほぼ一様になり，ブタンとオクタンも，ド デカンと同様に分子配向が凝縮・蒸発特性にほぼ影響しないことが分かった．また，図 12 に示すように，ブタ ン，オクタンとドデカンの分子配向はそれぞれ異なる確率密度分布を持ち，サンプリング数が極端に少ない領域 がある．図12より，ブタンはS=－0.1，オクタンとドデカン分子はS=0で検査面を通過した割合が高く，界面に 平行または垂直な配向でほとんど検査面を通過しないことが分かる．界面温度に依存せず，検査面を通過したす べての分子（Nall）に対しても，実際に凝縮・蒸発した分子（Nc , Ne）に対しても同様な配向確率分布を示してい る．これは，ブタン，オクタンとドデカンの分子配向は凝縮・蒸発確率に影響せず，鎖状分子に対して，分子重 心の並進運動のみを考慮して質点系でモデル化すれば式（1）を適用でき，前節の結果の妥当性を裏付けたとも言 える．要するに，系の温度が比較的高い場合は，分子配向が界面に平行または垂直な配向になりにくいため，気 液界面層（遷移状態）における直鎖状分子の並進運動に対する制限は，並進運動長さスケール比(vl /vg)1/3で近似的 に扱える．一方で系の温度が低い場合では，式（1）による凝縮・蒸発係数が元々1に近い値となるため，分子配 向の異方性を示したとしても，並進運動長さスケール比(vl /vg)1/3で近似した場合の誤差はさほど影響がないこと を推察できる．よって，鎖状分子の分子配向は凝縮・蒸発係数にほぼ影響しないことを言え，遷移状態説理論に よる凝縮・蒸発係数の理論式が鎖状分子にも適用可能である．
4. 結 論
分子動力学シミュレーションを用いて，異なる鎖長の鎖状分子の凝縮・蒸発特性を検討し，以下の結果を得た．
(1) 鎖状分子は，分子構造によらず，単純分子と同じ凝縮・蒸発特性を持ち，並進エネルギーの界面鉛直成分と 界面温度に依存する．計算誤差が生じるものの，分子動力学解析による凝縮・蒸発係数が遷移状態説理論に よる理論値とよく一致する．
(2) 気液界面層に出入りする鎖状分子の分子配向は分子構造に依存するが，分子の凝縮・蒸発特性にほぼ影響し ない．すなわち，分子重心の並進運動のみを考慮して質点系でモデル化すれば，遷移状態説理論による凝縮・
蒸発係数の理論式が鎖状分子にも適用可能である．
謝 辞
本研究の分子動力学計算は九州大学情報基盤センターのスーパーコンピューティングシステムを利用して行っ た．ここに記して感謝の意を表す．
文 献
(1) 新岡嵩，河野通方，佐藤順一，燃焼現象の基礎(2001)，pp. 195-220，オーム社．
(2) 三上真人，“燃料液滴の蒸発と燃焼”，ながれ，Vol. 20，No. 2 (2001)，pp. 106-115.
(3) 酒井真理，“インクジェット法による回路基板製造技術”，電子情報通信学会誌，Vol. 90, No. 7 (2007), pp.544-548.
(4) Akker, E.A.T.vd., Frijns, A.J.H., Kunkelmann, C., Hilbers, P.A.J., Stephan, P., and Steenhoven, A.A.v., “Molecular Dynamics simulation of the microregion”, International Journal of Thermal Sciences, Vol. 59 (2012), pp. 21-28.
(5) Fujikawa, S., Yano, T., Watanabe, M., Vapor-Liquid Interfaces, Bubbles and Droplets: Fundamentals and Applications (Heat and Mass Transfer) (2011), pp.19-69, Springer-Verlag.
(6) Sazhin, S.S., Xie, J.F., Shishkova, I.N., Elwardany, A.E., and Heikal, M.R., “A kinetic model of droplet heating and evaporation:
Effects of inelastic collisions and a non-unity evaporation coefficient”, International Journal of Heat and Mass Transfer, Vol. 56 (2013), pp. 525-537.
(7) Kryukov, A.P., Levashv, V.Y., and Sazhin, S.S., “Evaporation of diesel fuel droplets: kinetic versus hydrodynamics models”, International Journal of Heat and Mass Transfer, Vol. 47 (2004), pp. 2541-2549.
(8) Tokunaga, A., Yamawaki, S., Nagayama, G., and Tsuruta, T., “Effect of non-condensable gases on interface resistance in microscopic dropwise condensation”, Proceedings of the 4th International Conference on Heat Transfer and Fluid Flow in Microscale (2011), Paper No. HTFFM – IV - 118.
(9) Tsuruta, T., Tanaka, H., and Masuoka, T., “Condensation/evaporation coefficient and velocity distributions at liquid-vapor interface”, International Journal of Heat and Mass Transfer, Vol. 42 (1999), pp. 4107-4116.
(10) Nagayama, G., and Tsuruta, T., “A general expression for the condensation coefficient based on transition state theory and molecular dynamics simulation”, Journal of Chemical Physics, Vol. 118, No. 3 (2003), pp. 1392-1399.
(11) Tsuruta, T., and Nagayama, G., “Molecular dynamics studies on the condensation coefficient of water”, Journal of Physical Chemistry B, Vol. 108, No. 5 (2004), pp. 1736-1743.
(12) 鶴田隆治，長山暁子，“分子境界条件としての凝縮係数について”，日本機械学会論文集B編，Vol. 71, No. 705 (2005), pp. 1424-1427.
(13) 徳永敦士，鶴田隆治，長山暁子，“気液界面の非平衡輸送現象と反射分子のエネルギー適応係数”，日本機械学会論 文集B編，Vol. 77, No. 774 (2011), pp. 130-139.
(14) 徳永敦士，鶴田隆治，長山暁子，“逆温度こう配現象の非平衡分子動力学解析”，日本機械学会論文集B編，Vol. 76, No. 771 (2010), pp. 1901-1907.
(15) Matsumoto, M., “Molecular dynamics of fluid phase change”, Fluid Phase Equilibria, Vol. 144, No. 1-2 (1998), pp. 307-314.
(16) 水口博貴，長山暁子，鶴田隆治，“鎖状分子の蒸発係数に関する分子動力学的研究”，日本機械学会論文集B 編，
Vol. 77, No. 781 (2011), pp. 1826-1833.
(17) Cao, B.Y., Xie, J.F., and Sazhin, S.S., “Molecular dynamics study on evaporation and condensation of n-dodecane at liquid–
vapor phase equilibria”, Journal of Chemical Physics, Vol. 134, 164309 (2011), pp. 1-9.
(18) Xie, J.F., Sazhin, S.S., and Cao, B.Y., “Molecular dynamics study of the processes in the vicinity of the ndodecane vapour/liquid interface”, Physics of Fluids, Vol. 23, 112104 (2011), pp. 1-11.
(19) Cheng, S.F., Lechman, J.B., Plimpton, S.J., “Evaporation of Lennard-Jones fluids”, Journal of Chemical Physics, Vol. 134, 224704 (2011), pp. 1-13.
(20) National Institute of Standards and Technology, “NIST Chemistry WebBook”, NIST, http://webbook.nist.gov/chemistry/ (参照 日 2012年12月12日).
(21) Case, D.A., Cheatham, T.E., Darden, T., Gohlke, H., Luo, R., Merz, K.M., Onufriev, A., Simmerling, C., Wang, B., and Woods, R.J., “The Amber biomolecular simulation programs”, Journal of Computational Chemistry, Vol. 26, No. 16 (2005), pp. 1668- 1688.
(22) Wang, J.M., Wolf, R.M., Caldwell, J.W., Kollman, P.A., and Case, D.A., “Development and testing of a general amber force field”, Journal of Computational Chemistry, Vol. 25, No. 9 (2004), pp. 1157-1174.
(23) Alejandre, J., Tildesley, D. J., and Chapela, G. A., “Molecular-dynamics simulation of the orthobaric densities and surface- tension of water”, Journal of Chemical Physics, Vol. 102, No. 11 (1995), pp. 4574-4583.
(24) Gu, K., Watkins, C.B., and Koplik, J., “Molecular dynamics simulation of the equilibrium liquid–vapor interphase with solidification”, Fluid Phase Equilibria, Vol. 297 (2011), pp. 77-89.
(25) Hariharan, A., and Harris, J.G., “Structure and thermodynamics of the liquid-vapor interface of fluorocarbons and semifluorinated alkane diblocks – a molecular dynamics study”, Journal of Chemical Physics, Vol. 101, No. 5 (1994), pp. 4156- 4165.
