納米溝槽表面潤濕特性的分子動力學模擬

顏 笑，陳 鳳，張勤昭，黃志勇，薄涵亮

（清華大學 核能與新能源技術研究院，先進核能技術協同創新中心，先進反應堆工程與安全教育部重點實驗室，北京 100084）

潤濕性作為固體表面的重要性質，在諸多領域有重要應用。如，在高溫氣冷堆及壓水堆冷凝換熱器中，具有疏水性的換熱表面有利于蒸汽在壁面發生滴狀凝結，從而減小液膜熱阻，加速液滴脫落，提高換熱系數。

影響固體表面潤濕性的因素包括化學組分和表面微結構兩種。其中化學組分決定了液固間相互作用的形式和強度；在液固組分既定的情況下，固體表面微結構對潤濕性有著重要影響［1］。

近年來，借助分子動力學模擬，有關納米結構表面潤濕特性的研究成果相繼出現［2-6］。其中，Ambrosia等［2］研究了陣列結構表面對水滴接觸角的影響，發現隨著列柱高度的增加，液滴接觸角增大。Hirvi等［3］通過對液滴在聚乙烯表面接觸角的模擬指出，在陣列高度相同的情況下，方柱截面較小時表面疏水性更強。文獻［4-5］對液滴在粗糙表面的潤濕模式進行了研究，其中文獻［4］對一種介于Wenzel模式［7］和Cassie模式［8］之間的過渡潤濕（cross contact）模式進行了分析。此外，Chen 等［5］和Hirvi等［6］的研究表明，微觀尺度下液滴的接觸角在一定程度上與宏觀實驗及理論分析結果符合較好，這就拓寬了納米液滴在納米結構表面潤濕特性的研究意義。盡管人們在粗糙表面潤濕性的研究方面做了諸多工作，但對于納米結構表面潤濕特性的理解還有待深入，對納米結構表面潤濕特性的影響因素尚缺乏系統的研究。

本文擬采用分子動力學模擬，研究二維納米液滴在具有溝槽型結構表面的潤濕特性，考察不同液固相互作用強度下，液滴接觸角和潤濕模式隨溝槽深度和溝槽寬度的變化規律，為設計表面納米結構以控制表面潤濕性提供依據。

1 模擬方法

分子動力學模擬基于分子間相互作用的力場，通過對分子運動方程積分從而獲得分子運動軌跡。

如圖1所示，模擬對象為具有一定厚度的二維液滴，模擬盒子尺寸為35.3nm×2.4nm×31.4nm。壁面采用溝槽型結構（圖2），溝槽深度為H，寬度為S。水分子采用粗粒化水模型（mw）［9］。該模型是對復雜水模型的一種簡化，但對水的密度、三相分界線及水分子的徑向分布等特性的模擬效果較好，同時可取得較大的時間步長。壁面原子以金屬Pt為原型，按照面心立方進行排布。壁面與水分子間的相互作用由截斷的L-J勢描述：

式中：Vij、rij、rc、εij、σij分 別 為 原 子i、j 間 的 相互作用勢能、距離、勢能截斷半徑、能量參數和長度參數。Pt和O 原子間相互作用的具體參數通 過 混 合 法 則［10］得 到：能 量 參 數εPt-O＝5.746eV，長度參數σPt-O＝0.252 5nm。Pt和H 原子間的作用忽略不計。截斷半徑rc選為1nm，即當原子i、j 的間距rij＞rc時，水分子和壁面原子間無相互作用。在式（1）中引入能量系數a，對水分子和壁面原子間相互作用的強度進行調節，以便從化學角度獲得不同的親／疏水性。能量系數a 越大，水和壁面間的相互作用越強。

圖1 模擬盒子示意圖Fig.1 Sketch of simulation box

圖2 壁面結構示意圖Fig.2 Sketch of surface structure

選取NVT 系綜（固定分子數N、體積V 和溫度T），通過Nose-Hoover算法［11］保持溫度在300K。模擬中水分子數為9 840，相應最小液滴直徑大于10nm，此時接觸角的尺度效應可忽略［12-13］。對x、y 方向采用周期性邊界條件，垂直于冷凝表面的z 方向采用鏡面邊界條件。時間步長取10fs。

計算通過分子動力學大規模并行模擬軟件包LAMMPS［14］實現。在不同能量系數a、溝槽深度H 以及寬度S 下，對液滴在固體表面的狀態進行模擬。

接觸角的計算如圖3所示。假定液滴輪廓為一段圓弧，則接觸角可由液滴輪廓高度h 和基線長度b 通過式（2）計算：

在液滴整體位置和系統相互作用勢能趨于穩定后（趨穩過程約0.6ns），每隔一段時間對接觸角取樣，共取樣4次，求取平均值，即得最終接觸角。測量誤差所導致的接觸角不確定度在±2.0°以內。

2 模擬結果及分析

2.1 潤濕特性與溝槽深度的關系

如前所述，能量系數a 反映了液固分子間相互作用的強度。液固作用越強，液滴越易在固體表面浸潤，接觸角也就越小。在a 分別取值1、2、3、4的情況下，得到光滑表面接觸角依次為146.8°、121.2°、93.0°、62.1°，以90°為親、疏水性的分界線，則對應于a＝1、2、3的表面為本征疏水表面，而對應于a＝4的表面為本征親水表面。

在不同的能量系數a 下，保持溝槽寬度S＝0.8nm 不變，分別取溝槽深度H＝0、0.6、1.2、2.4、3.6、4.8、6.0nm，得到的接觸角及相應的潤濕模式如圖4所示。需要指出的是，在a＝1、H＝6.0nm 時，液滴與壁面幾乎分離，因此圖4中并未給出相應的接觸角。

圖4 接觸角和潤濕模式隨溝槽深度的變化Fig.4 Changes of contact angle and wetting regime with groove depth

由圖4可見，從溝槽深度H＝0nm 的光滑表面到H＝0.6nm 的結構表面，對應于能量系數a＝1、a＝2、a＝3及a＝4，接觸角均增大，增幅分別為6.9°、8.6°、35.2°和39.2°。該結果表明，無論是本征親水（a＝4）還是本征疏水（a＝1，2，3）表面，溝槽型結構均有利于增強表面疏水性，液固間相互作用越強，接觸角增大越明顯。

對應于能量系數a＝1、2、3、4的表面，溝槽深度超過0.6nm 以后，接觸角波動幅度分別為4.4°、9.5°、3.5°和10.8°。這表明，在溝槽寬度不變的條件下，溝槽深度對于接觸角的影響是有限的，在溝槽深度超過某一值（0.6nm）后，接觸角隨溝槽深度的增加基本不變，這與Hirvi等［3］和Jeong等［15］的模擬結果一致。由此可知，在溝槽寬度不變時，一味地增加溝槽深度對于表面疏水性的改變意義不大。

此外，溝槽深度的改變，可導致潤濕模式發生變化。圖5為液滴與結構表面接觸的3種不同的潤濕模式：Cassie模式［8］，即液滴只與溝槽頂層表面發生接觸（圖5a）；Wenzel模式［7］，即液滴滲入表面空隙，與溝槽內表面發生完全接觸（圖5c）；過渡潤濕模式［5］，該模式是介于Cassie和Wenzel模式之間的一種模式，即液滴部分滲入表面空隙，但不完全接觸空隙底部（圖5b）。由圖4可知，在能量系數a＝1，即液固相互作用很弱時，不同的溝槽深度下，液滴潤濕模式均為Cassie模式；在能量系數a＝4，即液固相互作用很強時，不同的溝槽深度下，潤濕模式均為Wenzel模式；在能量系數a＝2時，溝槽深度由0.6nm 增加到1.2nm 時，潤濕模式由Wenzel向Cassie模式轉變；在能量系數a＝3時，溝槽深度為0.6nm 時出現Wenzel模式，溝槽深度大于2.4nm 之后為過渡潤濕模式。由此可見，對于本征親／疏水性很強的表面，潤濕模式不隨溝槽深度發生轉變，而在本征親／疏水性較弱的情況下，隨著溝槽深度的增加，液滴與表面的潤濕模式可由Wenzel向過渡潤濕或Cassie模式轉變。

圖5 3種潤濕模式Fig.5 Three wetting regimes

2.2 潤濕特性與溝槽寬度的關系

在不同的能量系數a 下，保持溝槽深度H＝4.8nm 不變，分別取溝槽寬度S＝0、0.4、0.8、1.2、1.6nm，得到的接觸角及相應的潤濕模式如圖6所示。需要指出的是，在a＝1時，液滴與壁面分離，因此圖6并未給出相應的接觸角。

無論固體表面是本征疏水（a＝2，3）還是本征親水（a＝4），接觸角均隨溝槽寬度的增加而呈增大趨勢。且接觸角的最大增幅隨液固相互作用強度的增大而增大。對應于能量系數a＝2，3，4，接觸角最大增幅分別為36°、46°、60°。

另外，溝槽寬度的改變可導致潤濕模式發生轉變。對于本征疏水性較弱（a＝3）的表面，溝槽寬度較小（S＝0.4，0.8nm）時，液滴底部結構化的水分子將阻礙上部水分子沿縫隙的滲透［5］，導致液滴只能部分進入縫隙，即液滴與固體表面的接觸為過渡潤濕模式；溝槽寬度增至1.2nm 時，液滴能完全滲入空隙，與表面發生完全接觸。而對于本征疏水性較強（a＝2）或親水性較強（a＝4）的表面，在各溝槽寬度下，液滴在固體表面保持Cassie模式或Wenzel模式不變。

圖6 接觸角及潤濕模式隨溝槽寬度的變化Fig.6 Changes of contact angle and wetting regime with hollow width

圖7為能量系數a＝3、溝槽深度H＝4.8nm時，不同溝槽寬度下液滴與壁面的接觸狀態。

圖7 不同溝槽寬度下液滴與壁面的接觸狀態Fig.7 State of water droplet on surface with different groove widths

從以上結果可見，表面結構對接觸角的影響包括兩方面：1）表面溝槽結構構成液滴在固體表面鋪展和浸潤的阻力；2）液滴滲入溝槽，使得三相接觸線以上的液滴體積減小、高度降低，由式（2）可知，這將導致接觸角減小。圖8為能量系數a＝4 時，液滴在光滑表面（H ＝0nm）和粗糙表面（H＝0.6nm，S＝0.8nm）的接觸狀態。由于粗糙表面阻滯（pin）效應［16］的存在，使得液滴在粗糙表面難以自由鋪展開，減小了液滴在表面的鋪展寬度，從而使粗糙表面接觸角較光滑表面接觸角大。圖9為相同能量系數（a＝4）、溝槽寬度（S＝0.8nm）、相同潤濕模式（Wenzel模式）下，溝槽深度分別為4.8nm和6.0nm 時液滴在表面的形態。兩者鋪展寬度相同，但H＝4.8nm 時三相接觸線以上的液滴高度較H＝6.0nm 時的大，因而前者接觸角較后者大（分別為106.5°和104.2°，圖4）。這一效應在溝槽尺度較大的情況下應更為顯著。正是由于接觸角同時受到以上多方面因素的影響，在圖4和圖6中，接觸角并不隨溝槽深度或寬度某一單一因素單調變化。

圖8 光滑表面和粗糙表面液滴狀態Fig.8 State of water droplet on smooth and rough surface

圖9 溝槽深度為4.8nm 和6.0nm 時的液滴狀態Fig.9 State of water droplet on surfaces with groove depth of 4.8nm and 6.0nm

2.3 模擬結果與經典模型的比較

建立接觸角與表面結構之間的通用數學關系具有重要意義。將模擬結果與經典模型預測結果進行比較，以考察在本文條件下，接觸角隨表面結構的變化規律是否可用經典模型加以描述。

對于處于Cassie模式下的液滴，由Cassie-Baxter方程［8，17］可得表觀接觸角與本征接觸角的關系：

其中：θ為表觀接觸角；θ0為液滴在固體表面的本征接觸角；f 為面積份額，其定義為液固實際接觸面積與該表面的投影面積之比。

圖10為處于Cassie狀態下的液滴接觸角隨面積份額的變化，同時圖中還給出不同能量系數下的Cassie-Baxter方程接觸角預測值。由圖10 可見，接觸角的模擬結果與Cassie-Baxter方程符合較好。

圖10 接觸角模擬結果與Cassie-Baxter方程預測值比較Fig.10 Comparison between contact angles obtained with simulation and Cassie-Baxter equation

對于處于Wenzel模式下的液滴，表觀接觸角與本征接觸角之間滿足Wenzel方程［7］：

其中，r為粗糙因子，其定義為粗糙表面的實際面積與該表面的投影面積之比。

圖11為處于Wenzel狀態下的液滴接觸角隨表面粗糙因子的變化，圖中同時給出不同能量系數下的Wenzel方程接觸角的預測值。由圖11 可見，Wenzel方程并不適用。如對于能量系數a＝4 的表面（本征接觸角為62.1°），Wenzel方程預測值隨粗糙因子的增加而減小，而模擬值幾乎始終顯著大于模型預測值。Ambrosia等［2］也曾根據Wenzel理論成立的條件，即液滴尺寸需遠大于表面結構尺寸，指出Wenzel方程在微觀尺度下存在局限性。

圖11 接觸角模擬結果與Wenzel方程預測值比較Fig.11 Comparison between contact angles obtained with simulation and Wenzel equation

此外，由圖10、11可見，對于處于過渡潤濕模式的液滴，Cassie-Baxter和Wenzel兩種模型均不能很好地對接觸角進行預測。

3 結論

通過分子動力學模擬，考察了不同液固間相互作用強度下，溝槽深度、溝槽寬度對接觸角和潤濕模式的影響。結果表明：

1）引入溝槽結構使表面疏水性相對光滑表面增強，增幅與本征疏水性強弱有關；溝槽深度大于一定值后接觸角不隨溝槽深度發生顯著變化；對于本征疏水性較弱的表面，溝槽深度增加可導致潤濕模式發生變化，由Wenzel向Cassie或過渡潤濕模式轉變。

2）增加溝槽寬度，接觸角總體增大；對于本征疏水性較弱的表面，溝槽寬度增加有利于液滴滲入溝槽內部，使潤濕模式由過渡潤濕向Wenzel模式轉變。

3）在液滴與表面結構尺寸相差不大的情況下，對處于Cassie狀態下的液滴，模擬接觸角與Cassie-Baxter方程預測值符合較好；對于Wenzel狀態下的液滴，模擬接觸角較Wenzel方程預測值偏大。

［1］ CHAU T T，BRUCKARD W J，KOH P T L，et al.A review of factors that affect contact angle and implications for flotation practice［J］.Advances in Colloid and Interface Science，2009，150（2）：106-115.

［2］ AMBROSIA M S，HA M Y，BALACHANDAR S.The effect of pillar surface fraction and pillar height on contact angles using molecular dynamics［J］.Applied Surface Science，2013，282：211-216.

［3］ HIRVI J T，PAKKANEN T A.Enhanced hydrophobicity of rough polymer surfaces［J］.The Journal of Physical Chemistry B，2007，111（13）：3 336-3 341.

［4］ CHEN S，WANG J，MA T，et al.Molecular dynamics simulations of wetting behavior of water droplets on polytetrafluorethylene surfaces［J］.The Journal of Chemical Physics，2014，140（11）：114704.

［5］ CHEN S，WANG J，CHEN D.States of a water droplet on nanostructured surfaces［J］.The Journal of Physical Chemistry C，2014，118（32）：18 529-18 536.

［6］ HIRVI J T，PAKKANEN T A.Nanodroplet impact and sliding on structured polymer surfaces［J］.Surface Science，2008，602（10）：1 810-1 818.

［7］ WENZEL R N.Resistance of solid surfaces to wetting by water［J］.Industrial & Engineering Chemistry，1936，28（8）：988-994.

［8］ CASSIE A B D，BAXTER S.Wettability of porous surfaces［J］.Transactions of the Faraday Society，1944，40：546-551.

［9］ MOLINERO V，MOORE E B.Water modeled as an intermediate element between carbon and silicon［J］.The Journal of Physical Chemistry B，2008，113（13）：4 008-4 016.

［10］HIINENBERGER P H.Thermostat algorithms for molecular dynamics simulations［J］.Advances in Polymer Science，2005，173：105-147.

［11］PLIMPTON S.LAMMPS user's manual［M］.US：Sandia National Laboratory，2005.

［12］張福田.分子界面化學基礎［M］.上海：上海科學技術文獻出版社，2006：47-62.

［13］葛宋.流固界面能量傳遞及超臨界流體熱物性的分子動力學模擬［D］.北京：清華大學，2013.

［14］PLIMPTON S. Fast parallel algorithms for short-range molecular dynamics［J］.Journal of Computational Physics，1995，117（1）：1-19.

［15］JEONG W J，HA M Y，YOON H S，et al.Dynamic behavior of water droplets on solid surfaces with pillar-type nanostructures［J］.Langmuir，2012，28（12）：5 360-5 371.

［16］KOISHI T，YASUOKA K，FUJIKAWA S，et al.Measurement of contact-angle hysteresis for droplets on nanopillared surface and in the Cassie and Wenzel states：A molecular dynamics simulation study［J］.ACS Nano，2011，5（9）：6 834-6 842.

［17］KHAN S，SINGH J K.Wetting transition of nanodroplets of water on textured surfaces：a molecular dynamics study［J］.Molecular Simulation，2014，40（6）：458-468.