, ,

(大連理工大學(xué) 化工學(xué)院 ， 遼寧 大連 116024)

正十六烷納米液滴在光滑壁面上潤濕行為的分子動力學(xué)模擬

白麟,王寶和,于志家

(大連理工大學(xué) 化工學(xué)院 ， 遼寧 大連 116024)

采用分子動力學(xué)模擬技術(shù)，研究了正十六烷納米液滴在光滑壁面上的潤濕行為規(guī)律。模擬結(jié)果表明，壁面厚度、長度(或?qū)挾?、截斷半徑及分子數(shù)對接觸角的影響不大。隨著壁面作用勢能的增大，接觸角線性減?。划敱诿孀饔脛菽転?.5 kJ/mol時，接觸角約為90°。隨著模擬溫度的提高，接觸角逐漸減小。

分子動力學(xué) ； 模擬 ； 接觸角 ； 正十六烷

在現(xiàn)代社會中，經(jīng)濟的高速發(fā)展可能會產(chǎn)生大量被有機物污染的工業(yè)水資源；同樣，石油的開采、運輸以及存儲過程中，均易發(fā)生油品泄漏等污染事件[1-2]。與此同時，航空齒輪的潤滑和微/納機電系統(tǒng)的冷卻都與液滴的潤濕性有關(guān)。此外，航空齒輪潤滑油的鋪展特性極大地影響著軸承腔的潤滑與散熱功能，因此，開展油滴成膜的流動鋪展特性研究，為機械零部件潤滑計算提供準確合理的基礎(chǔ)數(shù)據(jù)，對于實現(xiàn)飛行器中的機械零部件精確潤滑設(shè)計十分重要[3]。為了更好地理解液滴潤濕鋪展行為機理，越來越多的研究學(xué)者將重點放在分子水平層次上。近年來，隨著計算機技術(shù)的迅猛發(fā)展，也使分子動力學(xué)模擬技術(shù)計算納米液滴的潤濕性質(zhì)成為可能。Werder等[4]利用分子動力學(xué)模擬了石墨表面上水液滴的接觸角與壁面作用勢能之間的作用關(guān)系。邱豐等[5]采用分子動力學(xué)模擬的方法，對Pb液滴在Cu基底上的鋪展?jié)櫇裥袨檫M行研究，發(fā)現(xiàn)晶體結(jié)構(gòu)對液滴鋪展具有較大影響。但目前還沒有關(guān)于烷烴類納米液滴潤濕性的分子動力學(xué)模擬研究的報道。本文將采用分子動力學(xué)模擬的方法，利用LAMMPS軟件模擬光滑壁面上正十六烷納米液滴的潤濕性質(zhì)，從壁面厚度和寬度(或長度)、分子數(shù)、壁面作用勢能及模擬溫度等方面，對接觸角的影響進行探究。

1 模擬方法

1.1模擬體系的建立

模擬體系的初始構(gòu)型如圖 1所示，壁面采用面心立方排布的Cu原子，晶格長度為0.361 5 nm，總共306 545個Cu原子，模擬盒子尺寸為47 nm×47 nm×20 nm。2 000個(研究分子數(shù)影響時除外)正十六烷分子隨機分布成球形液滴，正十六烷分子的初速度由隨機數(shù)發(fā)生器確定[6]。液滴質(zhì)心與壁面之間的初始距離為6.447 nm。

圖1 初始模型

1.2勢能模型的選取

當所研究的粒子中含有的原子數(shù)目較多時，通常采用的原子模型包括全原子模型和聯(lián)合原子模型兩種形式。在全原子模型中，將體系中的每個原子看作一個基本單元，由Material Studio軟件得到如圖 2所示的正十六烷分子結(jié)構(gòu)，其中灰色表示為C原子，白色表示H原子，并且在計算過程中定義正十六烷上每一個原子的參數(shù)，包括側(cè)鏈烷基上的H原子，這樣的力場稱為全原子力場。在分子動力學(xué)模擬中，為了簡化計算，有時將一些原子團如CH2被看作一個原子，稱作虛擬原子，如圖 3所示，黑色的CH3和灰色的CH2被當作相對分子質(zhì)量為15和14虛擬原子，這樣簡化的力場稱作聯(lián)合原子力場。

為了簡化計算，本模擬采用聯(lián)合原子力場模型。

圖2 全原子模型

圖3聯(lián)合原子模型

正十六烷分子間的勢能函數(shù)如式(1)所示[7]。

(1)

式中：U(tot)為總勢能；U(NB)和U(B)分別為非鍵勢能及鍵合勢能；N為正十六烷分子數(shù)；rij為i分子和j分子中兩個虛擬原子之間的距離；σij為i和j分子中兩個虛擬原子之間L-J勢能的尺度參數(shù)；εij為i和j分子中兩個虛擬原子之間L-J勢能的能量參數(shù)；kr和kθ分別為鍵長伸縮彈力系數(shù)和鍵角彎曲彈力系數(shù)；k1～k4為二面角扭轉(zhuǎn)勢能彈力系數(shù)；r、θ分別為鍵長、鍵角；r0、θ0分別為平衡鍵長、鍵角；φ為二面角。鍵長、鍵角、二面角作用參數(shù)如表1～3所示[8]。

表1 鍵長作用參數(shù)

表2 鍵角作用參數(shù)

表3 二面角作用參數(shù) kJ/mol

不同虛擬原子間或固體壁面原子與虛擬原子間的勢能函數(shù)仍為式(1)，其L-J 勢能的能量參數(shù)和尺度參數(shù)采用混合規(guī)則計算，如式(2)和(3)所示[9]。

(2)

(3)

式中：εls為虛擬原子與固體壁面原子之間L-J勢能的能量參數(shù)；εll為相同虛擬原子之間L-J 勢能的能量參數(shù)；εss為壁面原子之間L-J 勢能的能量參數(shù)。σls為虛擬原子與壁面原子之間L-J勢能的長度參數(shù)；σll為相同虛擬原子之間L-J勢能的長度參數(shù)；σss為壁面原子之間L-J勢能的長度參數(shù)(σss=0.234 nm)。虛擬原子L-J勢能參數(shù)如表4所示[8]。

表4 虛擬原子L-J勢能參數(shù)

1.3模擬細節(jié)

模擬在x、y方向采用周期性邊界條件，在z方向采用固壁和鏡像邊界條件；粒子間力的截斷半徑為1.367 nm，模擬時間步長為1 fs，總模擬時間為2 ns，前1 ns使得系統(tǒng)達到平衡，后1 ns統(tǒng)計計算并輸出系統(tǒng)的密度分布。采用正則系綜(NVT)，并用Woodcock控溫法維持體系溫度衡定；每隔1 000步矯正體系的質(zhì)心。模擬數(shù)據(jù)采用LAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)軟件計算得到。采用等密度擬合曲線法計算液滴密度和接觸角[10]。

2 結(jié)果與討論

2.1模擬參數(shù)的影響

2.1.1壁面厚度的影響

選擇47.00 nm×47.00 nm×20.00 nm的模擬盒子，正十六烷分子數(shù)為2 000，截斷半徑為1.367 nm，壁面作用勢能為1.5 kJ/mol，溫度為298 K，當壁面厚度D分別為1.084、1.446、1.807、2.169 nm時，進行計算模擬，得到圖 4所示的一維密度分布(D=1.446 nm，其他壁面厚度的類似)和圖5所示的二維密度圖(D=1.446 nm，其他壁面厚度的類似)，其中圖5中顏色灰度從淺到深代表密度從ρ=0到ρ=1.250 g/cm3。根據(jù)二維密度圖，再采用文獻[10]的方法，得到的接觸角如圖6所示。

圖4 正十六烷納米液滴的一維密度圖(D=1.446 nm)

圖5 正十六烷納米液滴的二維密度圖(D=1.446 nm)

圖6 壁面厚度(D)對接觸角(θ)的影響

由圖 4可見，液體主體密度有小幅度波動，統(tǒng)計平均值為0.799 2 kg/L，接近于常溫下正十六烷液體密度的實驗值(0.78 kg/L)[11]。對于壁面作用勢能1.5 kJ/mol的光滑壁面，當壁面厚度大于分子間作用力的截斷半徑時，不同壁面厚度下模擬得到的接觸角基本相同，約為30°，即壁面厚度對接觸角影響不大。因此，本文的模擬研究，其壁面厚度均采用1.446 nm。

2.1.2壁面寬(或長)度的影響

選擇盒子高度為20.00 nm，壁面厚度為1.446 nm，正十六烷分子數(shù)為2 000，截斷半徑為1.367 nm，壁面作用勢能為1.5 kJ/mol，溫度為298 K，壁面寬(或長)度即盒子寬(或長)度L分別為36.15、39.77、43.38、47.00、50.61 nm時，進行模擬計算，可以得到正十六烷納米液滴的二維密度分布圖，根據(jù)二維密度圖，再采用文獻[10]的方法，計算得到的接觸角如圖7所示。

圖7 壁面長度(或?qū)挾?(L)對接觸角(θ)的影響

由圖 7可以看出，截斷半徑小于模擬盒子的一半時，模擬得到的接觸角均為30°左右，即固體壁面寬(或長)度對接觸角影響不大。本文的模擬研究，其壁面寬(或長)度均采用47.00 nm。

2.1.3分子數(shù)的影響

將模擬盒子設(shè)置為47.00 nm×47.00 nm×20.00 nm，壁面厚度為1.446 nm，截斷半徑為1.367 nm，壁面作用勢能為1.5 kJ/mol，溫度為298 K，正十六烷分子數(shù)分別為1 000、1 500、2 000、2 500、3 000。進行模擬計算，可以得到正十六烷納米液滴的二維密度分布圖，根據(jù)二維密度圖，再采用文獻[10]的方法，計算得到的接觸角如圖8所示。

圖8 分子數(shù)對接觸角的影響

由圖8可以看出，對于壁面作用勢能為1.5 kJ/mol的光滑壁面，正十六烷分子數(shù)在1 000～3 000范圍內(nèi)，模擬得到的接觸角均為30°左右，即分子數(shù)對接觸角影響不大。因此，本文的模擬研究，選用的分子數(shù)為2 000個。

2.1.4截斷半徑的影響

將模擬盒子設(shè)置為47.00 nm×47.00 nm×20.00 nm，壁面厚度為1.446 nm，正十六烷分子數(shù)為2 000，壁面作用勢能為1.5 kJ/mol，溫度298 K，截斷半徑分別為1.171 5、1.366 7、1.562、1.757 3 nm時，進行模擬計算，可以得到正十六烷納米液滴的二維密度分布圖，根據(jù)二維密度圖，再采用文獻[10]的方法，計算得到的接觸角如圖9所示。

圖9 截斷半徑對接觸角的影響

眾所周知，當截斷半徑較小時，模擬結(jié)果不夠準確，誤差比較大；隨著截斷半徑的增加，模擬結(jié)果逐漸接近實際情況；當進一步增大截斷半徑時，模擬結(jié)果的變化不顯著。增加截斷半徑所帶來的結(jié)果使計算時間隨著截斷半徑的增加迅速增長，帶來相應(yīng)的模擬計算成本增加。由圖 9可知，對于壁面作用勢能為1.5 kJ/mol的光滑壁面，截斷半徑在1.171 5～1.757 3 nm范圍內(nèi)，模擬得到的接觸角均為30°左右，即截斷半徑對接觸角影響不大。因此，本文的模擬研究，選用的截斷半徑為1.367 nm。

2.2壁面作用勢能的影響

選擇47.00 nm×47.00 nm×20.00 nm的模擬盒子，壁面厚度為1.446 nm，正十六烷分子數(shù)為2 000，截斷半徑為1.367 nm，在溫度為298 K時，探討壁面作用勢能的影響。壁面作用勢能分別取0.2、0.4、0.6、0.9、1.2、1.5、1.8、2.1 kJ/mol時，進行模擬計算，可以得到正十六烷納米液滴的二維密度分布圖，根據(jù)二維密度圖，再采用文獻[10]的方法，計算得到的接觸角如圖10所示。

圖10 壁面作用勢能(εss)對接觸角的影響

由圖10可見，壁面作用勢能的增加，液滴接觸角呈線性減小，當壁面作用勢能在0.2～2.1 kJ/mol時，接觸角變化范圍為10°～133°。當壁面作用勢能為0.5 kJ/mol時，接觸角約為90°，為中性壁面；當壁面作用勢能<0.5 kJ/mol時，接觸角>90°，為疏油壁面；當壁面作用勢能>0.5 kJ/mol時，接觸角<90°，為親油壁面。

2.3溫度的影響

選擇47.00 nm×47.00 nm×20.00 nm的模擬盒子，壁面厚度為1.446 nm，正十六烷分子數(shù)為2 000，截斷半徑為1.367 nm，探討溫度對接觸角的影響。模擬溫度分別為298、323、348、373 K下，進行模擬計算，可以得到正十六烷納米液滴的二維密度分布圖，根據(jù)二維密度圖，再采用文獻[10]的方法，計算得到的接觸角如圖11所示。

圖11 溫度對接觸角的影響

由圖 11可知，無論是疏油壁面、中性壁面還是親油壁面，隨著溫度的升高，接觸角逐漸減小。根據(jù)Yong′s方程[12]：

cosθ=(γsv-γsl)/γlv

(4)

式中：γsv、γsl、γsv分別為固—氣、固—液和液—氣界面間的界面張力，θ為平衡接觸角。根據(jù)Mimouri等[13]的研究，可以認為，γlv基本為定值；溫度的升高，正十六烷液滴的界面張力γsl降低；與γsl項相比，γsv變化特別??；所以，接觸角隨著溫度的提高而逐漸減小。

3 結(jié)論

本文采用分子動力學(xué)模擬方法，研究了正十六烷納米液滴在光滑虛擬壁面上的潤濕行為。考察了壁面寬(或長)度、厚度、截斷半徑及正十六烷分子數(shù)，壁面作用勢能和溫度對接觸角的影響規(guī)律。模擬結(jié)果表明，壁面寬(或長)度、厚度、截斷半徑及十六烷分子數(shù)對接觸角影響不大。隨著壁面作用勢能的增大，接觸角線性減?。划敱诿孀饔脛菽転?.5 kJ/mol時，接觸角為90°左右。隨著溫度的提高，無論是親油、中性還是疏油壁面，接觸角都減小。

[1] Abdel Gawad S,Abdel Shafy M.Pollution control of industrial wastewater from soap and oil industries:a case study[J].Water Science and Technology,2002,46(4-5):77-82.

[2] Hassler B.Accidental versus operational oil spills from shipping in the baltic sea:risk governance and management strategies[J].Ambio,2011,40(2):170-178.

[3] 方 龍,陳國定,劉 登.噴射油滴沉積油膜的流動鋪展特性研究[J].機械工程學(xué)報,2016,52(23):160-167.

[4] Chai J,Liu S,Yang X. Molecular dynamics simulation of wetting on modified amorphous silica surface[J].Applied Surface Science,2009,255(22):9078-9084.

[5] 邱 豐,王 猛,周化光,等.Pb液滴在Ni基底潤濕鋪展行為的分子動力學(xué)模擬[J].物理學(xué)報,2013,62(12):1-7.

[6] Heermann D W.Computer-simulation methods[M].Berlin:Springer Berlin Heidelberg,1990.

[7] 陳正隆,徐為人,湯立達.分子模擬的理論與實踐[M].北京:化學(xué)工業(yè)出版社,2007.

[8] 顏 群,王寶和.水及其表面活性劑體系汽—液界面行為的分子動力學(xué)模擬[J].河南化工,2015,32(4):17-21.

[9] Rigby M,Smith E,Wakehan W,et al.The forces between molecules clarendon[J].Chem Phys Lett,1986,46(3):1-10.

[10] 王寶和,李 群.接觸角的研究現(xiàn)狀及其在凝膠干燥中的作用[J].干燥技術(shù)與設(shè)備,2014,12(1):39-46.

[11] 劉光啟,馬連湘,劉 杰.化學(xué)化工物性數(shù)據(jù)手冊(有機卷)[M].北京:化學(xué)工業(yè)出版社,2002.

[12] 劉天慶,穆春豐,夏松柏,等.滴狀冷凝初始液滴的形成機理[J].化工學(xué)報,2007,58(4):821-828.

[13] Langroudi S M M,Ghassemi M,Shahabi A,et al.A molecular dynamics study of effective parameters on nano-droplet surface tension[J].Journal of Molecular Liquids,2011,161(2):85-90.

MolecularDynamicsSimulationofWettingBehaviorofn-HexadecaneNanodropletsonSmoothSurfaces

BAILin,WANGBaohe,YUZhijia

(School of Chemical Engineering , Dalian University of Technology , Dalian 116024 , China)

The wetting behavior ofn-hexadecane nanodroplets on smooth surfaces is investigated by molecular dynamics simulation.Simulation results show that surface,thickness length (or width),cut-off radius and then-hexadecane molecular numbers have little influence on the contact angels.With surface potential energy increasing,the contact angles decrease linearly.When surface potential energy is 0.5 kJ/mol,the contact angel is about 90°.With the increasing of simulation temperature,the contact angle decreases gradually.

TB383，O363

A

1003-3467(2017)11-0025-05

Keywordsmolecular dynamics ; simulation ; contact angel ;n-hexadecane

2017-08-17

國家自然科學(xué)基金(51376030)

白 麟(1991-)，男，在讀碩士，從事油水分離及納米液滴潤濕行為的分子動力學(xué)模擬研究，E-mail：sneakerhead_bl@qq.com；聯(lián)系人：王寶和(1959-)，男，副教授，從事不同形貌微納結(jié)構(gòu)的制備、干燥及分子動力學(xué)模擬研究，E-mail:wbaohe@163.com。