崔子純 楊莫涵 阮曉鵬 范曉麗? 周峰3) 劉維民3)

1) (西北工業(yè)大學(xué)材料學(xué)院,先進(jìn)潤(rùn)滑與密封材料中心,西安 710072)

2) (西北工業(yè)大學(xué),倫敦瑪麗女王大學(xué)工程學(xué)院,西安 710072)

3) (中國(guó)科學(xué)院蘭州化學(xué)物理研究所,固體潤(rùn)滑國(guó)家重點(diǎn)實(shí)驗(yàn)室,蘭州 730000)

建立了基于第一性原理方法研究二維材料界面摩擦的高通量計(jì)算程序,該程序?qū)崿F(xiàn)了自動(dòng)化批量建模、批量提交任務(wù)、多任務(wù)并發(fā)計(jì)算,以及計(jì)算結(jié)果自動(dòng)收集、處理和圖像繪制,使用該程序可以節(jié)省時(shí)間.采用此程序計(jì)算了不同層間距離下雙層氮化硼和雙層石墨烯的滑移勢(shì)能面,及層間界面摩擦力和摩擦系數(shù).研究發(fā)現(xiàn),隨著層間距離減小,雙層氮化硼界面的平均摩擦力近似線性增大,摩擦系數(shù)為0.11—0.17,雙層石墨烯界面摩擦力先增大后減小再增大,其摩擦系數(shù)在12 nN 載荷下達(dá)到最小值(0.014),這些結(jié)果與已有研究結(jié)果一致,驗(yàn)證了該計(jì)算程序的可靠性.此外還研究了表面氫化和氟化對(duì)雙層氮化硼界面摩擦的影響,發(fā)現(xiàn)氟化氮化硼/氮化硼界面的摩擦系數(shù)更低.

1 引言

摩擦一般發(fā)生在相對(duì)運(yùn)動(dòng)或具有相對(duì)運(yùn)動(dòng)趨勢(shì)的接觸界面,阻礙相對(duì)運(yùn)動(dòng),同時(shí)產(chǎn)生能量損耗.Holmberg 等[1]調(diào)查發(fā)現(xiàn)卡車等重型交通工具有1/3 的燃料用于克服發(fā)動(dòng)機(jī)、變速器的摩擦.為了避免不必要的能源消耗和機(jī)器故障,控制摩擦力、黏附力和磨損力以減小摩擦磨損一直是科研工作者密切關(guān)注的問(wèn)題[2].納米尺度下,以黏附和摩擦為代表的表面力成為制約納米機(jī)電系統(tǒng)性能和壽命的關(guān)鍵因素[3],因此,了解納米級(jí)甚至原子尺度下的摩擦行為必要且重要.

因?yàn)閷娱g作用力弱、層內(nèi)共價(jià)鍵強(qiáng)、剪切應(yīng)力較低的特點(diǎn),層狀二維材料的潤(rùn)滑性能往往優(yōu)于其他納米材料[4,5].特別是石墨烯[6]、氮化硼[7]、二硫化鉬[8]等典型二維層狀材料,通常以微米級(jí)薄膜的形式被用作固體潤(rùn)滑劑.已有研究表明,二維材料界面處的原子尺度摩擦行為受到多種因素影響,包括溫度[9?11]、滑動(dòng)速度[12]以及載荷[12?14]等.研究者采用實(shí)驗(yàn)技術(shù)和理論計(jì)算方法研究了二維材料界面的原子尺度摩擦行為[15?17].溫度變化會(huì)影響表面原子熱振動(dòng)的幅度以及頻率,從而影響摩擦行為,已有研究顯示摩擦系數(shù)隨溫度變化呈現(xiàn)非單調(diào)變化的特點(diǎn)[18].摩擦力與滑動(dòng)速度的關(guān)系較為復(fù)雜,當(dāng)滑動(dòng)速度較小時(shí),摩擦力幾乎不隨速度變化,能量耗散主要來(lái)自于黏滑運(yùn)動(dòng),當(dāng)滑動(dòng)速度增大到一定程度時(shí),摩擦力與速度呈現(xiàn)線性關(guān)系,黏滑運(yùn)動(dòng)消失,能耗主要為阻尼能耗[18?20].

依照宏觀摩擦定律,摩擦力與載荷的關(guān)系符合阿蒙頓定律,原子尺度上摩擦力與載荷的關(guān)系則較為復(fù)雜.Mate 等[21]利用鎢探針研究了石墨表面的微觀黏滑現(xiàn)象,發(fā)現(xiàn)界面間的平均摩擦力隨法向載荷的增大而近似線性增大;Mo 等[22]采用分子動(dòng)力學(xué)方法模擬了碳尖端在金剛石表面的滑動(dòng)行為,當(dāng)接觸界面的范德瓦耳斯吸附力較強(qiáng)時(shí),摩擦力和載荷呈現(xiàn)非線性關(guān)系,隨著范德瓦耳斯吸附力的減弱,摩擦力與載荷的關(guān)系趨于線性.An 等[23]基于第一性原理方法計(jì)算了氮化硼層間摩擦,發(fā)現(xiàn)摩擦力隨載荷的增大而增大;Sun 等[24]發(fā)現(xiàn)石墨烯層間平均滑移摩擦力隨著載荷增大,先增大后減小再增大的有趣現(xiàn)象.此外,多項(xiàng)研究表明[25?30],層數(shù)、表面吸附、面內(nèi)應(yīng)變等對(duì)石墨烯和氮化硼的層間摩擦有著顯著影響,歸因于層間電子分布的變化.鑒于原子尺度摩擦學(xué)的復(fù)雜性和重要性,從原子水平上掌握摩擦機(jī)制[31],預(yù)測(cè)和控制摩擦行為緊迫且關(guān)鍵.

原子力顯微鏡(atomic force microscope,AFM)研究掃描探針與材料表面之間的摩擦特性,是研究納米尺度摩擦學(xué)應(yīng)用最廣泛的工具.然而,受限于可用作探針尖端的材料種類,AFM 測(cè)量任意二維材料界面間摩擦仍然是一個(gè)挑戰(zhàn).先進(jìn)的計(jì)算模擬技術(shù)是研究二維材料層間摩擦的有力工具.研究摩擦現(xiàn)象的常用計(jì)算模擬方法包括分子動(dòng)力學(xué)模擬(molecular dynamics,MD)方法和基于密度泛函理論(density functional theory,DFT)的第一性原理方法.分子動(dòng)力學(xué)依據(jù)牛頓力學(xué)定律模擬接觸界面的相對(duì)運(yùn)動(dòng),根據(jù)公式(U代表原子間相互作用勢(shì),r代表原子的所在位置)求出每個(gè)原子的受力,然后分別對(duì)所有原子的法向力和與滑動(dòng)方向相反的橫向力進(jìn)行求和,得到總的法向力和摩擦力,進(jìn)而運(yùn)用阿蒙頓定律求得摩擦系數(shù).分子動(dòng)力學(xué)模擬已被廣泛應(yīng)用于計(jì)算二維材料摩擦學(xué)性能[32?35],其模擬的可靠性強(qiáng)烈依賴于選取的勢(shì)函數(shù)描述層間/內(nèi)相互作用的準(zhǔn)確性.第一性原理計(jì)算方法通過(guò)計(jì)算摩擦界面相對(duì)滑動(dòng)時(shí)的能量變化,得到滑動(dòng)能壘和平均摩擦力,通過(guò)改變層間距離模擬外界施加載荷,獲得不同載荷下的勢(shì)能面以及摩擦力,再根據(jù)阿蒙頓定律求得摩擦系數(shù).該方法可以精確地處理電子結(jié)構(gòu),并從界面電荷分布探索原子尺度摩擦行為[36?38].

勢(shì)能面(potential energy surface,PES)是研究原子尺度摩擦行為的主要途徑[39].特定載荷下的PES 由足夠密集的勢(shì)能點(diǎn)組成,對(duì)應(yīng)相應(yīng)滑移位置處的勢(shì)能.構(gòu)建PES 一般需要建立近700 個(gè)滑移位置的界面結(jié)構(gòu)并計(jì)算其能量,涉及大量建模、計(jì)算、數(shù)據(jù)處理.本文建立了基于第一性原理方法研究二維材料層間摩擦的高通量計(jì)算程序(homogenous/heterogenous junction construction and frictional properties calculation software,HJC2S),該程序自動(dòng)化建立滑移界面的原子結(jié)構(gòu),批量建立構(gòu)建勢(shì)能面所需計(jì)算的界面構(gòu)型,高通量計(jì)算系列載荷下勢(shì)能面上的所有界面結(jié)構(gòu),自動(dòng)收集能量并輸出摩擦性能相關(guān)數(shù)據(jù)和圖像.基于該程序計(jì)算了雙層氮化硼和雙層石墨烯的界面層間滑移摩擦性能,研究了表面功能化對(duì)于氮化硼層間滑移摩擦力的調(diào)控.

2 二維材料層間摩擦的高通量計(jì)算程序

程序HJC2S 基于Linux 操作系統(tǒng),所有計(jì)算通過(guò)基于密度泛函理論[40]的VASP (viennaabinitiosimulation package)軟件包[41?43]完成.程序功能包含了自動(dòng)化建模、批量產(chǎn)生并提交計(jì)算任務(wù)、數(shù)據(jù)提取與存儲(chǔ)(提取并保存數(shù)據(jù)、繪制并保存圖像) 3 個(gè)模塊.如圖1 所示,程序首先基于用戶提供的二維材料結(jié)構(gòu)建立同、異質(zhì)結(jié)界面結(jié)構(gòu),并對(duì)界面結(jié)構(gòu)進(jìn)行優(yōu)化;將界面一側(cè)的材料沿著z方向按照一定間距移動(dòng),模擬系列載荷;特定界面間距下,沿著x,y方向按照一定間距移動(dòng)界面一側(cè)的材料,模擬特定載荷下沿xy平面相對(duì)滑移位置;提交對(duì)所有界面結(jié)構(gòu)進(jìn)行自洽計(jì)算的任務(wù),通過(guò)作業(yè)調(diào)度系統(tǒng)實(shí)現(xiàn)計(jì)算資源的有效利用和多任務(wù)并發(fā)計(jì)算,實(shí)現(xiàn)二維材料界面摩擦性能的高通量計(jì)算;計(jì)算完成后,程序智能化處理并輸出計(jì)算結(jié)果,用戶分析數(shù)據(jù)以獲得界面的納米摩擦學(xué)性能,以下部分介紹具體的計(jì)算流程.

圖1 二維材料界面摩擦學(xué)特性高通量計(jì)算流程示意圖Fig.1.The schematic diagram showing the procedure of high-throughput calculation of the tribological property at the interface of two-dimensional materials.

2.1 構(gòu)建界面結(jié)構(gòu)

對(duì)于異質(zhì)界面,HJC2S 程序會(huì)對(duì)提供的兩種二維材料的晶胞矢量進(jìn)行匹配(本文設(shè)置晶格常數(shù)失配率不超過(guò)5%,失配角不大于3°),構(gòu)建異質(zhì)界面結(jié)構(gòu).程序會(huì)提供各種可能晶格矢量的界面結(jié)構(gòu)及該結(jié)構(gòu)中界面兩側(cè)的原子數(shù)和總原子數(shù),以供用戶挑選合適的界面模型.

2.2 模擬外加載荷

HJC2S 程序通過(guò)改變界面間距模擬外加載荷,平衡間距對(duì)應(yīng)零載荷.首先以平衡間距的界面結(jié)構(gòu)為初始結(jié)構(gòu),設(shè)定步長(zhǎng)并逐步減小間距,生成系列間距界面結(jié)構(gòu)并對(duì)每個(gè)間距界面結(jié)構(gòu)建立文件夾.然后批量提交計(jì)算系列間距界面結(jié)構(gòu)的任務(wù),并計(jì)算系列間距下界面結(jié)合能Eb.計(jì)算完成后,程序會(huì)自動(dòng)提取并保存結(jié)合能以及界面間距的數(shù)據(jù),進(jìn)行數(shù)據(jù)處理并獲得界面結(jié)合能隨界面距離z變化的關(guān)系,通過(guò)公式求解外加載荷和界面間距之間的對(duì)應(yīng)關(guān)系,自動(dòng)繪制結(jié)合能隨界面間距、載荷隨界面間距變化的曲線圖.

2.3 計(jì)算不同層間距離的勢(shì)能面

HJC2S 程序通過(guò)固定界面一側(cè)材料,使另一側(cè)材料在xy平面做相對(duì)橫向滑移來(lái)模擬界面相對(duì)滑動(dòng)過(guò)程.界面被劃分為M×N的網(wǎng)格(網(wǎng)格點(diǎn)數(shù)根據(jù)x,y方向長(zhǎng)度設(shè)定),每一個(gè)格點(diǎn)位置對(duì)應(yīng)相對(duì)滑移中的一個(gè)界面結(jié)構(gòu).格點(diǎn)劃分的致密程度決定著勢(shì)能面計(jì)算的精度.為了獲得系列層間距下的勢(shì)能面,對(duì)每一個(gè)層間距離建立文件夾,每個(gè)文件夾下有M×N個(gè)子文件夾,對(duì)應(yīng)不同的滑移位置.計(jì)算完成后,程序會(huì)提取每個(gè)結(jié)構(gòu)對(duì)應(yīng)的滑移位置坐標(biāo)和能量數(shù)據(jù),保存并以此繪制系列層間距離下的勢(shì)能面.

2.4 尋找最小能量路徑

HJC2S 程序依據(jù)計(jì)算得到的勢(shì)能面和最小能量路徑起點(diǎn)和終點(diǎn),采用String 方法[44]在勢(shì)能面上尋找最小能量路徑.String 方法是一種優(yōu)秀的過(guò)渡態(tài)搜索算法,該方法以提供的起點(diǎn)和終點(diǎn)為初始值,確定兩個(gè)勢(shì)能極小值,然后采樣多條連接兩個(gè)能量極小值的曲線,通過(guò)構(gòu)造曲線微分方程,使曲線逐漸逼近最小能量路徑.確定最小能量路徑后提取并輸出路徑位置坐標(biāo)及能量,繪制最小能量路徑及其上的靜橫向力曲線.

2.5 計(jì)算平均摩擦力

HJC2S 程序?qū)⒑Y出特定載荷下勢(shì)能面上的極值點(diǎn),選取鄰近的最大值點(diǎn)和最小值點(diǎn),通過(guò)兩點(diǎn)的坐標(biāo)以及能量,根據(jù)公式計(jì)算得到界面平均摩擦力,其中 ?V=Vmax?Vmin(Vmax為勢(shì)能極大值,Vmin是勢(shì)能極小值),?x是Vmax和Vmin所在位置之間的距離.程序自動(dòng)計(jì)算系列外加載荷下的界面平均摩擦力〈Ff〉,輸出外加載荷-平均摩擦力的對(duì)應(yīng)數(shù)據(jù)并繪制外加載荷與界面平均摩擦力關(guān)系曲線.

3 計(jì)算方法

3.1 程序中涉及物理量的計(jì)算

HJC2S 程序采用如下公式計(jì)算界面結(jié)合能Eb:

其中,Etotal是特定層間距下滑移界面的總能量,Eup,Edown分別是界面兩側(cè)材料的能量.外加載荷Fn與結(jié)合能Eb之間滿足關(guān)系[45]:

程序由(1)式計(jì)算系列間距下界面結(jié)合能,由(2)式得到法向載荷與層間距的關(guān)系.

勢(shì)能的計(jì)算公式如下:

其中V(x,y,Fn)和Eb(x,y,Fn)為Fn載荷下、(x,y)位置處的勢(shì)能和界面結(jié)合能,z(x,y,Fn)為Fn載荷對(duì)應(yīng)的界面間距,Fnz(x,y,Fn)為抵抗外加載荷Fn所做的功,V0(Fn)是Fn載荷下勢(shì)能面上的最小勢(shì)能點(diǎn),其計(jì)算公式如下:

Eb(x0,y0,z(x0,y0,Fn))為零勢(shì)能位置滑移體系的界面結(jié)合能,Fnz(x0,y0,Fn) 為零勢(shì)能位置處抵抗外加載荷Fn的層間作用能,(x0,y0)為Fn載荷下零勢(shì)能位置坐標(biāo).

特定載荷下界面平均摩擦力〈Ff〉由如下公式[45]計(jì)算:

其中?V表示該載荷下最大滑移能壘,由勢(shì)能最大值與最小值決定

?x表示勢(shì)能最大值點(diǎn)與最小值點(diǎn)之間的距離.摩擦系數(shù)μ與平均摩擦力和載荷之間的關(guān)系為

3.2 本文采用的計(jì)算方法和設(shè)置

文章所有計(jì)算均采用基于密度泛函理論[40]的VASP 程序包[41?43]進(jìn)行.采用廣義梯度近似(general gradient approximate,GGA)的Perdew-Burke-Ernzerhof(PBE)[46]方法計(jì)算電子交換關(guān)聯(lián)能,采用平面投影綴加波贗勢(shì)(projector augmented wave,PAW)方法[47]描述電子和離子之間相互作用,采用范德瓦耳斯色散修正方法(density functional dispersion correction)DFT-D3 描述相鄰層之間的范德瓦耳斯弱相互作用.在倒易空間中,電子波函數(shù)通過(guò)平面波基組擴(kuò)展,平面波基組截?cái)嗄茉O(shè)定為550 eV,能量和力的收斂精度分別設(shè)置為10–5eV和0.02 eV/?.在總能計(jì)算過(guò)程中,不弛豫原子位置.布里淵區(qū)K點(diǎn)網(wǎng)格取樣采用Gamma 方法,K網(wǎng)格大小根據(jù)具體模型由Vaspkit[48]生成.為了消除z方向上周期性結(jié)構(gòu)之間的相互作用,沿垂直方向添加了20 ?厚度的真空層.

4 結(jié)果與討論

本文采用HJC2S 程序計(jì)算了雙層氮化硼(BN/BN)、雙層石墨烯(Gr/Gr)、氫化氮化硼/氮化硼(H-BN/BN)和氟化氮化硼/氮化硼(F-BN/BN)的層間滑移摩擦行為.為了模擬不同載荷下的層間滑移摩擦行為,將兩層之間層間距以步長(zhǎng)0.1 ?逐步壓縮,研究了層間距離為4.0—1.5 ?的Gr/Gr 和層間距離為4.0—2.5 ?的BN/BN,H-BN/BN 和F-BN/BN 的勢(shì)能面、界面摩擦力和摩擦系數(shù).

4.1 不同載荷下雙層氮化硼界面摩擦力

首先研究了BN/BN 的界面摩擦行為,測(cè)試HJC2S 程序的可靠性.圖2(a),(b)分別為BN/BN處于Top(T)位置和Hollow(H)位置的俯視圖和側(cè)視圖.Top 位置結(jié)合能為29 meV/atom,層間距為3.39 ?,與Marom 等[49]采用PBE+Vdw 方法計(jì)算得到的BN/BN 界面層間距(3.37 ?)吻合.采用2×2 超胞模擬BN/BN 界面層間相對(duì)滑移,將Top 位置設(shè)置為初始位置,沿著xy平面滑動(dòng)距離設(shè)置為5 ? × 5 ?,用于構(gòu)建勢(shì)能面的網(wǎng)格大小為0.2 ? × 0.2 ?,第一布里淵區(qū)K網(wǎng)格大小設(shè)為7 × 7 × 1.

圖2 BN/BN,H-BN/BN 和 F-BN/BN 雙層體系原子結(jié)構(gòu)的俯視圖和側(cè)視圖 (a) Top 和(b) Hollow 位置BN/BN 的俯視圖和側(cè)視圖;(c) H-BN/BN 和(d) F-BN/BN 在Hollow 位置的俯視圖和側(cè)視圖Fig.2.Top and side views for atomic structures of BN/BN,H-BN/BN,and F-BN/BN bilayer.Top and side views of BN/BN bilayer in (a) Top and (b) Hollow positions.Top and side views of (c) H-BN/BN and (d) F-BN/BN bilayers in Hollow position.

界面結(jié)合能與層間距的關(guān)系是計(jì)算中最重要的一部分,該曲線擬合的精確度決定著計(jì)算結(jié)果的準(zhǔn)確度.圖3(a)展示了Top 位置的BN/BN 界面結(jié)合能隨層間距的變化.隨著層間距縮小,界面結(jié)合能逐漸變小,當(dāng)層間引力達(dá)到最大值時(shí),界面結(jié)合能達(dá)到最小值;之后,界面結(jié)合能開始增大,當(dāng)界面結(jié)合能為0 eV 時(shí),引力和斥力達(dá)到平衡;之后,斥力在層間作用中起主導(dǎo)作用.通過(guò)改變BN/BN 層間距來(lái)模擬載荷,根據(jù) (2)式得到載荷與層間距的關(guān)系.如圖3(b)所示,隨著BN/BN 層間距減小,載荷不斷增大.

圖3 BN/BN 位于Top 位置的界面結(jié)合能和載荷隨層間距的變化 (a) 結(jié)合能;(b) 載荷Fig.3.The variation of interlayer binding energy (Eb) and load with interlayer distance for BN/BN bilayer at Top position: (a) Eb and (b) load.

圖4(a)—(i)為0—8 nN 載荷下BN/BN 沿著xy平面分布的勢(shì)能.由于BN/BN 結(jié)構(gòu)本身的周期性,所以勢(shì)能也呈現(xiàn)周期性分布.如圖4 所示,0 nN 載荷下,BN/BN 的勢(shì)能面比較平滑,Top 位置處勢(shì)能最大,Hollow 位置處勢(shì)能最小.這是因?yàn)?相較于Hollow 位置,Top 位置的正對(duì)原子更多,層間相互作用更強(qiáng).隨著載荷增大,抵抗負(fù)載的層間作用增大,勢(shì)能面起伏逐漸增大,能壘也逐漸增大,勢(shì)能最大值和最小值對(duì)應(yīng)的位置并未發(fā)生變化.高勢(shì)能起伏增大了滑移阻力及因摩擦耗散的能量,根據(jù)(5)式,平均摩擦力也在逐漸增大.

圖4 顯示了各載荷下的最小能量路徑(勢(shì)能面上連接勢(shì)能最小值點(diǎn)的路徑中勢(shì)能起伏最小的路徑),圖5 為0—8 nN 載荷下最小能量路徑上的能量變化曲線.從圖5 可以看出,由于結(jié)構(gòu)的對(duì)稱性,能量隨位移變化也具有周期性.沿最小能量路徑進(jìn)行滑移時(shí),勢(shì)能先增大后減小,然后再重復(fù)此過(guò)程.另外,隨著載荷的增大,最小能量路徑上的能壘逐漸增大,即界面滑移的阻力增大,從而導(dǎo)致更高的摩擦.

圖4 不 同載 荷下BN/BN 的勢(shì) 能面 (a) 0 nN;(b) 1 nN;(c) 2 nN;(d) 3 nN;(e) 4 nN;(f) 5 nN;(g) 6 nN;(h) 7 nN;(i) 8 nN.T 和H 分別表示Top 和Hollow 位置,白色帶箭頭線代表最小能量路徑Fig.4.Potential energy surface of BN/BN bilayer under normal load of (a) 0 nN,(b) 1 nN,(c) 2 nN,(d) 3 nN,(e) 4 nN,(f) 5 nN,(g) 6 nN,(h) 7 nN and (i) 8 nN.T and H present the Top and Hollow positions,respectively.The white lines with arrow represent the minimum energy path.

圖5 0—8 nN 載荷下BN/BN 最小能量路徑上勢(shì)能的變化Fig.5.The variation of potential energy along the minimum energy path of BN/BN bilayer under 0–8 nN loads.

從圖6(a)可以看出,隨著載荷增大,BN/BN界面的平均摩擦力近似線性增大,表現(xiàn)出類阿蒙頓定律的特征,表明宏觀摩擦力與載荷之間的關(guān)系擴(kuò)展到了納米尺度的BN/BN 滑移界面.從圖6(b)可以看出,當(dāng)載荷從1 nN 增大到8 nN 時(shí),BN/BN滑移界面的摩擦系數(shù)從0.17 降至0.11.因表面結(jié)構(gòu)、環(huán)境、計(jì)算方法等不同,測(cè)得的六方氮化硼界面的摩擦系數(shù)為0.025—0.7.Martin 等[50]發(fā)現(xiàn)超高真空下h-BN 與h-BN 界面的摩擦系數(shù)為0.7,而在潮濕空氣中的摩擦系數(shù)為0.1.Wang 等[51]指出h-BN 固體潤(rùn)滑劑的適用溫度范圍為–184—538 ℃,相應(yīng)的摩擦系數(shù)為0.1—0.2.An 等[23]采用DFTLDA 泛函描述BN/BN 界面的相互作用,研究了BN/BN 界面摩擦系數(shù)與載荷的變化關(guān)系,他們指出,當(dāng)載荷從1 nN 增大到5 nN 時(shí),沿著x方向滑移摩擦系數(shù)從0.1 下降至0.025,在5—10 nN 載荷下,摩擦系數(shù)從0.025 上升到0.13.本文研究了不同載荷下BN/BN 的勢(shì)能面與摩擦力,摩擦系數(shù)計(jì)算結(jié)果與已有研究結(jié)果一致.

圖6 BN/BN,F-BN/BN,H-BN/BN 的平均界面摩擦力(Ff)和摩擦系數(shù)(μ)隨載荷的變化 (a)平均摩擦力;(b)摩擦系數(shù)Fig.6.The variation of friction force at the interface of BN/BN,F-BN/BN,and H-BN/BN bilayers with respect to the normal load: (a) Averaged interfacial friction (Ff);(b) friction coefficient (μ).

4.2 不同載荷下雙層石墨烯界面摩擦力

以0.1 ?步長(zhǎng)將層間距離從4 ?壓縮到1.5 ?,計(jì)算了系列層間距下雙層石墨烯(上下層均為2 個(gè)原子)之間的結(jié)合能,得到結(jié)合能隨層間距變化的關(guān)系,根據(jù)(2)式得到載荷與層間距的關(guān)系,以此模擬0—13.5 nN 外加載荷.將沿xy平面滑動(dòng)距離設(shè)置為2.6 ? × 2.6 ?,網(wǎng)格大小設(shè)為0.2 ? × 0.2 ?,第一布里淵區(qū)K網(wǎng)格設(shè)為12 × 12 × 1.計(jì)算了1—13.5 nN 外加載荷下Gr/Gr 的勢(shì)能面,結(jié)果如圖7 所示.從圖7 可以發(fā)現(xiàn),當(dāng)外加載荷為1 nN時(shí),勢(shì)能面較為平坦,Top 位置勢(shì)能最高,Hollow位置勢(shì)能最低.隨著載荷升高,勢(shì)能最大值先增大后減小.當(dāng)載荷達(dá)到12 nN 時(shí),勢(shì)能面發(fā)生了反轉(zhuǎn),即原來(lái)高勢(shì)能點(diǎn)(Top 位置)變?yōu)榈蛣?shì)能點(diǎn),低勢(shì)能點(diǎn)(Hollow 位置)變?yōu)楦邉?shì)能點(diǎn),勢(shì)能面趨于平坦,能壘很低.隨著載荷的進(jìn)一步增大,勢(shì)能起伏再次逐漸增大.

圖7 不同載荷下Gr/Gr 的勢(shì)能面 (a) 1 nN;(b) 4 nN;(c) 7 nN;(d) 10 nN;(e) 12 nN;(f) 13.5 nN.(a)—(d)中的H/T 對(duì)應(yīng)最小/最大勢(shì)能,而(e),(f)中的H/T 對(duì)應(yīng)最大/最小勢(shì)能Fig.7.Potential energy surface of graphene/graphene (Gr/Gr) bilayer under normal load of (a) 1 nN,(b) 4 nN,(c) 7 nN,(d) 10 nN,(e) 12 nN,and (f) 13.5 nN.T and H present the top and hollow positions,respectively.H/T in (a)–(d) is the position of minimum/maximum potential energy,whereas,T/H in (e),(f) is the position of minimum/maximum potential energy.

在勢(shì)能計(jì)算的基礎(chǔ)上,計(jì)算了1.0—13.5 nN外加載荷下Gr/Gr 界面的平均摩擦力和摩擦系數(shù).圖8(a)所示為平均摩擦力隨載荷的變化曲線.如圖所示,隨著載荷升高,Gr/Gr 界面的平均摩擦力逐漸增大;當(dāng)載荷增大到10 nN 時(shí),平均摩擦力開始下降,直到載荷為12 nN 時(shí),平均摩擦力達(dá)到極小值;之后,平均摩擦力隨著載荷的增大再次逐漸增高.圖8(b)所示為1.0—13.5 nN 外加載荷下Gr/Gr 的界面摩擦系數(shù).從圖中可以看出,隨著載荷的升高,摩擦系數(shù)先減小后增大.當(dāng)載荷為12 nN 時(shí),摩擦系數(shù)最小(0.014),這是因?yàn)?2 nN 下的勢(shì)能面最平坦,界面滑移需要克服的能壘最低,表現(xiàn)出極低的摩擦系數(shù).

圖8 Gr/Gr 界面平均摩擦力和摩擦系數(shù)隨載荷的變化 (a)平均摩擦力;(b)摩擦系數(shù)Fig.8.The variation of friction force at the interface of graphene/graphene (Gr/Gr) bilayers with respect to the normal load:(a) Averaged interfacial friction and (b) friction coefficient.

本文計(jì)算的摩擦力隨載荷變化的趨勢(shì)成功復(fù)現(xiàn)了Sun 等[24]的研究結(jié)果.此外,最小摩擦系數(shù)出現(xiàn)的層間距為1.75 ?,與Zhang 等[52]得出的極低摩擦現(xiàn)象發(fā)生的層間距(1.78 ?)相近,與Sun 等[24]得出的極低摩擦現(xiàn)象發(fā)生的層間距(1.825 ?)略有差距,這是因?yàn)镾un 等[24]是通過(guò)固定部分原子方法計(jì)算總能,而本文采用了固定所有原子的方法計(jì)算總能,本文計(jì)算方法可以節(jié)省計(jì)算資源并減少計(jì)算時(shí)間.Fan 等[53]在實(shí)驗(yàn)中測(cè)得的石墨烯的摩擦系數(shù)為0.1—0.2,Wang 等[27]計(jì)算結(jié)果顯示,載荷為1—9 nN 時(shí)石墨烯界面的摩擦系數(shù)為0.05—0.225.本文計(jì)算顯示石墨烯界面摩擦系數(shù)為0.014—0.14,與文獻(xiàn)[53]計(jì)算結(jié)果一致.這些均證明了HJC2S 程序和其采用算法研究界面摩擦行為的可靠性.

4.3 表面改性調(diào)控BN/BN 界面摩擦力

圖2(c),(d)分別為H—BN/BN 和F—BN/BN的俯視圖與側(cè)視圖.H—BN 和F—BN 中B—H,N—H,B—F,N—F 鍵長(zhǎng)分別為1.20,1.03,1.36,1.44 ?.H—BN/BN 和F—BN/BN 位于Top 位置平衡結(jié)構(gòu)的層間距分別為2.91 ?和3.20 ?,xy平面晶格常數(shù)為2.54 ? 和 2.57 ?.采用HJC2S程序,構(gòu)建了勢(shì)能面網(wǎng)格,大小均為5 ? × 5 ?,網(wǎng)格尺寸為0.2 ? × 0.2 ?,計(jì)算了0—8 nN 載荷下H—BN/BN 和F—BN/BN 的勢(shì)能面及界面平均摩擦力和摩擦系數(shù).

圖6(a)為BN/BN,F—BN/BN 和H—BN/BN滑移界面平均摩擦力隨載荷變化的曲線,如圖所示,1—8 nN 載荷下,3 種滑移界面平均摩擦力均隨載荷增大而增大.圖6(b)所示為3 種滑移界面摩擦系數(shù)隨載荷變化的曲線.隨著載荷增大,H—BN/BN 界面摩擦系數(shù)增大,F—BN/BN 界面摩擦系數(shù)先減小后緩慢升高,摩擦系數(shù)值為0.08—0.10.當(dāng)載荷小于6 nN 時(shí),H—BN/BN 滑移界面的摩擦系數(shù)小于BN/BN 滑移界面的摩擦系數(shù);當(dāng)載荷大于6 nN 時(shí),H—BN/BN 滑移界面的摩擦系數(shù)大于BN/BN 滑移界面的摩擦系數(shù).可以明顯地看出F—BN/BN 滑移界面的平均摩擦力和摩擦系數(shù)均小于BN/BN 和H—BN/BN,表明表面氟化顯著降低了摩擦.

電荷轉(zhuǎn)移是影響摩擦力的主要因素[54?56],為了探析氟化降低BN 界面摩擦的內(nèi)在原因,計(jì)算了BN/BN 以及F-BN/BN 的電荷分布.圖9為0 nN和8 nN載荷下BN/BN 和F-BN/BN在勢(shì)能最高位置和最低位置處的差分電荷密度圖.由于F 原子的電負(fù)性較強(qiáng),F-BN/BN 層間電子聚集在F 原子附近,層間電荷轉(zhuǎn)移較小;BN/BN 層內(nèi)電子結(jié)合能力較弱,層間電荷轉(zhuǎn)移較為明顯.從圖9可以看出,0 nN 和8 nN 載荷下,BN/BN 由Top位置(勢(shì)能最大值點(diǎn))滑移到Hollow 位置(勢(shì)能最小值點(diǎn))時(shí),電荷分布變化顯著.氟化后BN 界面從高勢(shì)能點(diǎn)滑移到低勢(shì)能點(diǎn)位置時(shí),電荷分布無(wú)明顯變化,所以滑移過(guò)程中的能壘較低,層間摩擦小,因此降低了BN 界面摩擦.

圖9 BN/BN 和F-BN/BN 的最高和最低勢(shì)能態(tài)在0 nN 和8 nN 載荷下的差分電荷密度 (a) BN/BN,0 nN;(b) BN/BN,8 nN;(c) F-BN/BN,0 nN;(d) F-BN/BN,8 nN.黃色代表電子聚集,藍(lán)色代表電子損失;左邊和右邊的圖分別對(duì)應(yīng)最高和最低勢(shì)能態(tài);圖(a)和(c)中的等值面設(shè)置為0.00006e/Bohr3;圖 (b) 和 (d) 中的等值面設(shè)置為0.0003e/Bohr3Fig.9.Differential charge density of BN/BN and F-BN/BN bilayer at the highest and lowest potential energy state under 0 nN and 8 nN loads.BN/BN bilayer under (a) 0 nN and (b) 8 nN load,F-BN/BN bilayer under (c) 0 nN and (d) 8 nN loads.The yellow color represents electron aggregation,and the blue color represents electron dissipation.The isosurface value in figure (a) and (c) is set to 0.00006e/Bohr3;the isosurface value in figure (b) and (d) is set to 0.0003e/Bohr3.

4.4 方法的適用性與局限性

DFT 方法可以深入到量子尺度,從電子分布的角度解釋摩擦現(xiàn)象.盡管該方法已經(jīng)被廣泛的用來(lái)預(yù)測(cè)層間滑移的摩擦系數(shù),但該方法仍存在一些局限性.首先,在層間相對(duì)滑動(dòng)的過(guò)程中,本文將上下層均設(shè)為剛性,并未考慮層間彎曲的影響.其次,本文通過(guò)改變計(jì)算體系的堆疊構(gòu)型以模擬層間滑動(dòng),因此計(jì)算得到的實(shí)際結(jié)果本質(zhì)是靜摩擦,而非動(dòng)摩擦.另外,DFT 計(jì)算本身具有局限性,如在計(jì)算過(guò)程中無(wú)法直接引入溫度,因此該方法并未考慮溫度的影響;計(jì)算通常是對(duì)數(shù)十個(gè)原子進(jìn)行,無(wú)法對(duì)幾百個(gè)原子進(jìn)行計(jì)算.

5 結(jié)論

報(bào)道了一個(gè)基于第一性原理方法計(jì)算二維材料層間滑移摩擦行為的高通量、自動(dòng)化計(jì)算程序.該程序?qū)崿F(xiàn)了同、異質(zhì)結(jié)界面自動(dòng)化建模、計(jì)算任務(wù)自動(dòng)生成及提交、智能提取并存儲(chǔ)計(jì)算結(jié)果、自動(dòng)處理數(shù)據(jù)并計(jì)算摩擦力和摩擦系數(shù),極大地縮短了建模和處理數(shù)據(jù)的時(shí)間.采用此程序模擬計(jì)算了外加載荷下雙層氮化硼層間滑移的摩擦,得到的摩擦系數(shù)與實(shí)驗(yàn)和已有計(jì)算結(jié)果一致,同時(shí)研究了雙層石墨烯層間滑移的摩擦力和摩擦系數(shù),計(jì)算結(jié)果與文獻(xiàn)結(jié)果相吻合,證明了該程序計(jì)算結(jié)果的可靠性.此外,我們研究了表面改性對(duì)雙層氮化硼層間滑移摩擦的影響,結(jié)果表明表面氟化可以很好地降低摩擦,提供了降低六方氮化硼界面摩擦的表面改性方法.