接触角与液固界面热阻关系的分子动力学模拟*

葛宋 陈民

(清华大学工程力学系，北京 100084)

(2012年11月9日收到;2013年1月30日收到修改稿)

1 引言

在微纳机电系统中由于尺度的减小导致面体比显著增大，界面效应往往变得不可忽略甚至可能成为主导因素[1].液固界面处的界面热阻及由此产生的界面温度跳跃是重要的界面现象[2]，其对微纳米系统中能量传递的影响可由热阻长度来判断.热阻长度表征的是在液体内部形成与界面温度跳跃相同的温度差所需的液体层厚度[3].实验和模拟结果都表明在较疏水的液固界面，热阻长度可达到数十纳米[4,5].表面湿润性也是一种重要的界面性质，而静态接触角是衡量表面湿润性的主要物理量.相比于界面热阻，目前对于表面湿润性和接触角的研究更加充分.一方面表征固体表面湿润性的接触角在宏观实验中容易测量，另一方面，现有的材料加工技术已能方便地制造出具有不同湿润性的表面[6-8].界面热阻和表面湿润性作为两种重要的界面现象，如果能建立起界面热阻和接触角之间的联系，将为微纳米系统中的材料选择和热设计提供极大的便利.已有学者注意到液固界面热阻与表面湿润性的联系，并试图建立起接触角和界面热阻的函数关系[9-11].Murad和Puri[9]推断亲水的表面将有利于减小液固界面热阻.Shenogina等[11]利用在水和不同化学成分的固体表面构成的界面上的模拟结果提出界面热导G(界面热阻的倒数)与接触角之间存在如下关系：G=B(1+cosθ)，其中B为比例系数，θ为接触角，但此关系式的适用性还有待进一步验证.本文利用分子动力学方法模拟了液体在固体表面的接触角及相同条件下的液固界面热阻，探讨了两者之间的关系，有助于加深对界面热阻产生机理的理解，同时能为实际应用提供理论参考.

2 模拟方法与细节

2.1 接触角的模拟方法

接触角的模拟系统如图1所示.本文中所有粒子间的相互作用均采用12-6 Lennard-Jones(LJ)势能函数来描述：

其中σ为分子直径，ε为能量参数，r为原子之间的距离.模型液体为氩，分子直径σ=3.405×10-10m，能量参数ε=1.67×10-21J.采用LJ模型固体.由于LJ模型固体熔点正比于固体原子间相互作用的能量参数εSS(约为0.5εSS/kB，kB为玻尔兹曼常数)，模拟中将εSS取得足够大以使其在模拟中能始终保持固态.此模型不需要引入弹簧力来维持固态，能较真实地描述固体的热学性质[4].固体原子的直径与氩原子相同.固体壁面由7层按fcc晶格结构排列的原子构成，晶格常数为a0=1.56σ，其(1，0，0)面与氩液滴及氩蒸汽相接触.最底层的原子在模拟中保持固定.

图1 接触角的模拟体系

模拟系统在三个方向上均使用周期性边界条件，模拟的是二维液滴.模拟盒子的尺寸为100 a0×10 a0×45 a0，x和z方向的尺度显著大于y方向.模拟体系中包含7200个氩原子和14000个壁面原子.液滴在平衡后的直径大于10 nm，这样能保证获得的接触角受尺度效应的影响较小.模拟中体系的演化采用速度Verlet算法，截断半径取为1 nm，积分步长为5 fs.通过Langevin热浴来控制体系的温度.模拟的前250万步用来使体系达到平衡，然后对位形进行取样来统计液体的密度分布.每100步取样一次，共获得多于5000个样本.将模拟盒子在xz平面内划分成网格，网格的大小为1˚A×1˚A，以此获得体系的密度分布.将液滴密度和蒸汽密度之和的二分之一的等密度线视为液滴的边界[12]来获得液滴的轮廓.同Leroy等[13]的做法相同，本文利用液滴轮廓的高度h和基线长度b来计算接触角，如图2所示.由于液体在液固界面的分层现象，我们将基线的位置取在离固体第一层原子距离为d=σ处[14].

2.2 液固界面热阻的模拟方法

液固界面热阻的模拟体系如图3所示，由两个相等的液体区域及中间的固体区域构成.模拟盒子的尺寸为10 a0×10 a0×56 a0，其中固体区域为10 a0×10 a0×19 a0.液体区域的密度保持与接触角模拟中的液体密度相同.固体模型与接触角模拟中的固体壁面模型相同.

图2 接触角的计算方法

图3 液固界面热阻的模拟体系

三个方向均采用周期性边界条件.截断半径取为1 nm，步长取为0.5 fs，能使体系的能量保持守恒.为了防止固体区域产生移动，在模拟中每隔10步去掉固体区域质心的动量.首先利用控温使体系在设定的温度下达到平衡，随后将控温移除，在固体和液体的图示区域分别设置热源和热沉，在热源中不断输入热量，在热沉中移走热量.这是通过增加或减小相应区域原子的动能来实现的.体系整体将保持能量守恒.平衡后，体系中将建立起稳定的温度分布，界面温度跳跃ΔT及热流q.将体系沿z方向划分成条状区域，区域的宽度为a0，统计各区域中的温度即可获得温度分布和界面温度跳跃.每一步加入的热量为ΔE，体系中的热流为q=ΔE/(2AΔt)=620 MW/m2.(其中A为模拟体系在xy平面的截面积，Δt为时间步长，分母中的因子2来源于体系的对称性).液固界面热阻R或界面热导G则可由定义R=ΔT/q和G=q/ΔT来分别获得.

为了获得接触角和界面热阻的关系，我们通过如下方法来改变固液界面的性质及固体的性质并检验二者的关系：1)保持固体的性质不变，调整液固间相互作用的势能参数εSL来改变液固间相互作用强度;2)保持液固间相互作用势能参数εSL不变，调整固体的性质.

3 结果与讨论

3.1 接触角

杨氏方程描述了接触角θ与液固表面张力γLS，液气表面张力γLV以及气固表面张力γSV之间的关系[15]：

引入黏附功H12的概念，其表征的是分开单位面积的液固界面所需做的功.界面黏附功H12也可由表面张力来表示：

联系(2)式和(3)式可见存在如下关系：

同时，假定固体和液体均匀分布，黏附功H12可近似表述为[16]

其中ρS，ρL分别为固体和液体的数密度，uLS为液固分子间势能函数.对于采用的12-6 Lennard-Jones势能，黏附功H12正比于势能参数εLS：H12∝εSL[16].当液体的种类给定时，γLV为常数.将(5)式代入(4)式：

可见液固界面相互作用强度是决定接触角的主要物理量，且当固体的晶格排列方式固定时，接触角将不受固体原子质量和固体原子之间相互作用强度的影响.

利用分子动力学模拟来验证上述推论.首先考虑液固间相互作用的影响.保持固体的性质不变(此时令 mS=0.75mAr，εSS=10εAr)，调整液固间相互作用的势能参数εLS来改变液固间相互作用强度.模拟得到的接触角θ随液固相互作用强度的变化如图4所示，图4中的插入图还给出了接触角的余弦值.可见接触角θ随εSL的增大而减小，即亲水性增加，并且1+cosθ与液固相互作用强度εSL较好地满足线性关系，与理论分析给出的(6)式及文献中的模拟结果[15]符合.

图4 不同液固相互作用强度下的接触角

保持液固间相互作用势能参数εSL=0.5εAr不变，调整固体的性质 (保持 mS=0.75mAr，εSS从 6εAr增大到 14εAr;或保持 εSL=10εAr，mS从 20 amu 增大到80 amu)来模拟对接触角的影响.固体原子间相互作用强度εSS和固体原子质量mS对接触角的影响如图5所示.两者对接触角产生的影响均在±4%以内.这也与前面的理论分析结果一致.

3.2 液固界面热阻与接触角的关系

分析不同因素(液固间相互作用强度，固体原子质量及固体原子间相互作用强度)对界面热阻的影响.液固间相互作用强度是影响液固界面能量传递的重要因素[17].先考虑液固间相互作用的影响：保持固体的性质不变，调整固液间相互作用强度 (能量参数 εLS由 0.4εAr增大到 2.4εAr，此时表面由疏水变化到亲水).模拟结果显示，界面热导随液固间相互作用强度增大而增大，如图6所示.显然液固间相互作用强度的增加促进了界面处的能量传递，这与文献的模拟结果一致[10,17].

保持液固间相互作用势能参数εSL=0.5εAr不变，改变固体原子质量mS或固体原子间能量参数εSS，固体的力学性质和热学性质如声速，热振动频率等都会改变.分别验证这两个因素对界面热阻的影响.保持固体原子质量不变，εSS从6εAr增大到14εAr;或保持εSS=10εAr不变，将固体原子质量 mS从20 amu增大到80 amu，界面热导的变化如图7所示.

图5 固体性质对接触角余弦的影响 (a)能量参数εSS;(b)固体原子质量mS

图6 液固界面热导随液固间相互作用强度的变化

由图7可见，在保持液固间相互作用不变的情况下，固体的性质对界面热阻(界面热导)有重要的影响，且mS与εSS有相反的作用.我们从液固原子间振动耦合的角度来分析固体间原子结合强度与原子质量对液固界面热阻产生影响的原因.

液固原子间的振动耦合越好，即两种原子的振动频率分布的重合程度越好，更多的能量将以简谐振动的方式传递，界面热阻将越小[18,19].原子热振动频率与原子质量和原子间结合强度近似存在如下关系：ω∝(εSS/mS)1/2.相较于固体而言，由于液体原子间较弱的结合力，液体原子的振动分布在较低的频率范围内[20]，而固体原子将在较高的频率内间振动，两者间存在一定的重合.此频率重合程度即决定了液固原子间的振动耦合度.固体原子质量和固体间相互作用强度对振动耦合度的影响，可从原子的振动态密度的角度来分析.

图7 固体性质对界面热导的影响 (a)能量参数εSS;(b)固体原子质量mS

利用原子的速度自相关函数(velocity autocorrelation fuction，VACF)可获得贴壁液体层，界面固体层的振动态密度(vibrational density of state，VDOS)，其表征了原子振动的频率分布.速度自相关函数定义为

其中v为原子的速度，〈〉代表对不同时间原点的统计平均.原子的振动态密度可由速度自相关函数的傅里叶变换获得

计算得到的不同固体原子质量下的振动态密度分布如图8(a)和(b)所示.对比图8(a)和(b)，显然，增大固体原子质量时固体原子振动的频率分布向低频率移动，增大了界面处液固间振动频率的重合程度，即图中重合区域面积增大，界面热导将随之增大;类似地也能得出，增强固体原子间相互作用强度将使固体原子的振动向高频率区移动，此时界面处液固间的振动频率重合程度降低，将使得界面热导减小.这个模拟结果与前面的预测是一致的.

图8 原子振动的态密度分布 (a)mS=0.5mAr;(b)mS=mAr

根据接触角和界面热导随不同因素的变化趋势，可归纳出接触角与界面热导间的关系，如图9所示.如果只改变液固间相互作用强度，界面热导G与接触角θ的余弦间存在近似的线性关系，这与Shenogina等[11]提出的G=B(1+cosθ)关系相符(其中B为比例系数).但是当接触角不变，固体原子质量与原子间相互作用强度改变时，界面热导的变化与接触角的余弦近似无关.因此界面热导与接触角不存在单值对应关系，液固界面热阻间与接触角之间也不存在单值的对应关系.

图9 改变不同参数时界面热导与接触角的关系

4 结论

本文利用分子动力学方法分别模拟了液体在固体表面的接触角及液固界面热阻，并探讨了二者间的联系.分别通过改变液固间相互作用强度和固体的性质来分析接触角和界面热阻的变化趋势.模拟结果显示接触角随液固间相互作用增强而减小，接触角的余弦与液固间的相互作用能量参数呈线性关系，但接触角并不受固体原子质量和固体间相互作用强度的影响.另一方面，界面热阻随液固间相互作用增强而减小;同时，在液固相互作用强度不变的情况下，固体原子间结合强度与固体原子质量对界面热阻有显著影响，这是由于这两个因素导致固体的振动频率发生变化，使得液固原子间的振动耦合度改变，从而使得界面热阻改变.本文的模拟表明，界面热阻与接触角间不存在单值对应关系，不能单一地将接触角作为液固界面热阻的评价标准.

本文的计算在清华大学信息科学与技术国家实验室(筹)的高性能计算平台上完成.

[1]Cahill D G，Ford W K，Goodson K E，Majumdar A，Mariset H J，Merlin R，Phillpot S R 2010 J.Appl.Phys.93 793

[2]Swartz E T，Pohl R O 1989 Rev.Mod.Phys.61 605

[3]Barrat J L，Chiaruttini F 2003 Mol.Phys.101 1605

[4]Xue L，Keblinski P，Phillipot S R，Choi S U S，Eastman J A 2003 J.Chem.Phys.118 337

[5]Ge Z B，Cahill D G，Braun P V 2006 Phys.Rev.Lett.96 186101

[6]Gu C Y，Di Q F，Shi L Y，Wu F，Wang W C，Yu Z B 2008 Acta Phys.Sin.57 3071(in Chinese)[顾春元，狄勤丰，施利毅，吴非，王文昌，余祖斌2008物理学报57 3071]

[7]Ma H M，Hong L，Yin Y，Xu J，Ye H 2011 Acta Phys.Sin.60 098105(in Chinese)[马海敏，洪亮，尹伊，许坚，叶辉 2011物理学报 60 098105]

[8]Gong M G，Xu X L，Cao Z L，Liu Y Y，Zhu H M 2009 Acta Phys.Sin.58 1885(in Chinese)[公茂刚，许小亮，曹自立，刘远越，朱海明2009物理学报58 1885]

[9]Murad S，Puri I K 2008 Appl.Phys.Lett.92 133105

[10]Wang Y，Keblinski P 2011 Appl.Phys.Lett.99 073112

[11]Shenogina N，Godawat R，Keblinski P，Garde S 2009 Phys.Rev.Lett.102 156101

[12]Shi B，Dhir V K 2009 J.Chem.Phys.130 034705

[13]Leroy F，M¨uller-Plathe F 2010 J.Chem.Phys.133 044110

[14]Voronov R S，Papavassiliou D V，Lee L L 2006 J.Chem.Phys.124 204701

[15]Sedlmeier F，Janecek J，Sendner C，Bocquet L，Netz R R，Horinek D 2008 Biointerphases 3 23

[16]Rowlinson J，WidomB 1982 Molecular Theory of Capillarity(Oxford：Oxford University Press)p86

[17]Maruyama S，Kimura T 1999 Therm.Sci.Eng.7 63

[18]Kikugawa G，Ohara T，Kawaguchi T，Torigoe E，Hagiwara Y，Matsumoto Y 2009 J.Appl.Phys.130 074706

[19]Issa K M，Mohamad A A 2012 Phys.Rev.E 85 031602

[20]Huxtable S T，Cahill D G，Shenogin S，Keblinski P 2005 Chem.Phys.Lett.407 129