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

, ,

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

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

白麟,王寶和,于志家

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

采用分子動力學模擬技術，研究了正十六烷納米液滴在光滑壁面上的潤濕行為規律。模擬結果表明，壁面厚度、長度(或寬度)、截斷半徑及分子數對接觸角的影響不大。隨著壁面作用勢能的增大，接觸角線性減小；當壁面作用勢能為0.5 kJ/mol時，接觸角約為90°。隨著模擬溫度的提高，接觸角逐漸減小。

分子動力學 ； 模擬 ； 接觸角 ； 正十六烷

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

1 模擬方法

1.1模擬體系的建立

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

圖1 初始模型

1.2勢能模型的選取

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

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

圖2 全原子模型

圖3聯合原子模型

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

(1)

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

表1 鍵長作用參數

表2 鍵角作用參數

表3 二面角作用參數 kJ/mol

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

(2)

(3)

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

表4 虛擬原子L-J勢能參數

1.3模擬細節

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

2 結果與討論

2.1模擬參數的影響

2.1.1壁面厚度的影響

選擇47.00 nm×47.00 nm×20.00 nm的模擬盒子，正十六烷分子數為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。根據二維密度圖，再采用文獻[10]的方法，得到的接觸角如圖6所示。

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

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

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

由圖 4可見，液體主體密度有小幅度波動，統計平均值為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，正十六烷分子數為2 000，截斷半徑為1.367 nm，壁面作用勢能為1.5 kJ/mol，溫度為298 K，壁面寬(或長)度即盒子寬(或長)度L分別為36.15、39.77、43.38、47.00、50.61 nm時，進行模擬計算，可以得到正十六烷納米液滴的二維密度分布圖，根據二維密度圖，再采用文獻[10]的方法，計算得到的接觸角如圖7所示。

圖7 壁面長度(或寬度)(L)對接觸角(θ)的影響

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

2.1.3分子數的影響

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

圖8 分子數對接觸角的影響

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

2.1.4截斷半徑的影響

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

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

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

2.2壁面作用勢能的影響

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

圖11 溫度對接觸角的影響

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

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

(4)

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

3 結論

本文采用分子動力學模擬方法，研究了正十六烷納米液滴在光滑虛擬壁面上的潤濕行為。考察了壁面寬(或長)度、厚度、截斷半徑及正十六烷分子數，壁面作用勢能和溫度對接觸角的影響規律。模擬結果表明，壁面寬(或長)度、厚度、截斷半徑及十六烷分子數對接觸角影響不大。隨著壁面作用勢能的增大，接觸角線性減小；當壁面作用勢能為0.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].機械工程學報,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基底潤濕鋪展行為的分子動力學模擬[J].物理學報,2013,62(12):1-7.

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

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

[8] 顏 群,王寶和.水及其表面活性劑體系汽—液界面行為的分子動力學模擬[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] 王寶和,李 群.接觸角的研究現狀及其在凝膠干燥中的作用[J].干燥技術與設備,2014,12(1):39-46.

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

[12] 劉天慶,穆春豐,夏松柏,等.滴狀冷凝初始液滴的形成機理[J].化工學報,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

國家自然科學基金(51376030)

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