高通量計算二維材料界面摩擦*

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

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

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

3) (中國科學(xué)院蘭州化學(xué)物理研究所,固體潤滑國家重點實驗室,蘭州 730000)

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

1 引言

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

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

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

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

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

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

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

圖1 二維材料界面摩擦學(xué)特性高通量計算流程示意圖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)

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

2.2 模擬外加載荷

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

2.3 計算不同層間距離的勢能面

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

2.4 尋找最小能量路徑

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

2.5 計算平均摩擦力

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

3 計算方法

3.1 程序中涉及物理量的計算

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

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

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

勢能的計算公式如下:

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

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

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

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

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

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

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

4 結(jié)果與討論

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

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

首先研究了BN/BN 的界面摩擦行為,測試HJC2S 程序的可靠性.圖2(a),(b)分別為BN/BN處于Top(T)位置和Hollow(H)位置的俯視圖和側(cè)視圖.Top 位置結(jié)合能為29 meV/atom,層間距為3.39 ?,與Marom 等[49]采用PBE+Vdw 方法計算得到的BN/BN 界面層間距(3.37 ?)吻合.采用2×2 超胞模擬BN/BN 界面層間相對滑移,將Top 位置設(shè)置為初始位置,沿著xy平面滑動距離設(shè)置為5 ? × 5 ?,用于構(gòu)建勢能面的網(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)系是計算中最重要的一部分,該曲線擬合的精確度決定著計算結(jié)果的準確度.圖3(a)展示了Top 位置的BN/BN 界面結(jié)合能隨層間距的變化.隨著層間距縮小,界面結(jié)合能逐漸變小,當層間引力達到最大值時,界面結(jié)合能達到最小值;之后,界面結(jié)合能開始增大,當界面結(jié)合能為0 eV 時,引力和斥力達到平衡;之后,斥力在層間作用中起主導(dǎo)作用.通過改變BN/BN 層間距來模擬載荷,根據(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平面分布的勢能.由于BN/BN 結(jié)構(gòu)本身的周期性,所以勢能也呈現(xiàn)周期性分布.如圖4 所示,0 nN 載荷下,BN/BN 的勢能面比較平滑,Top 位置處勢能最大,Hollow 位置處勢能最小.這是因為,相較于Hollow 位置,Top 位置的正對原子更多,層間相互作用更強.隨著載荷增大,抵抗負載的層間作用增大,勢能面起伏逐漸增大,能壘也逐漸增大,勢能最大值和最小值對應(yīng)的位置并未發(fā)生變化.高勢能起伏增大了滑移阻力及因摩擦耗散的能量,根據(jù)(5)式,平均摩擦力也在逐漸增大.

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

圖4 不 同載 荷下BN/BN 的勢 能面 (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 最小能量路徑上勢能的變化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)系擴展到了納米尺度的BN/BN 滑移界面.從圖6(b)可以看出,當載荷從1 nN 增大到8 nN 時,BN/BN滑移界面的摩擦系數(shù)從0.17 降至0.11.因表面結(jié)構(gòu)、環(huán)境、計算方法等不同,測得的六方氮化硼界面的摩擦系數(shù)為0.025—0.7.Martin 等[50]發(fā)現(xiàn)超高真空下h-BN 與h-BN 界面的摩擦系數(shù)為0.7,而在潮濕空氣中的摩擦系數(shù)為0.1.Wang 等[51]指出h-BN 固體潤滑劑的適用溫度范圍為–184—538 ℃,相應(yīng)的摩擦系數(shù)為0.1—0.2.An 等[23]采用DFTLDA 泛函描述BN/BN 界面的相互作用,研究了BN/BN 界面摩擦系數(shù)與載荷的變化關(guān)系,他們指出,當載荷從1 nN 增大到5 nN 時,沿著x方向滑移摩擦系數(shù)從0.1 下降至0.025,在5—10 nN 載荷下,摩擦系數(shù)從0.025 上升到0.13.本文研究了不同載荷下BN/BN 的勢能面與摩擦力,摩擦系數(shù)計算結(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 ?步長將層間距離從4 ?壓縮到1.5 ?,計算了系列層間距下雙層石墨烯(上下層均為2 個原子)之間的結(jié)合能,得到結(jié)合能隨層間距變化的關(guān)系,根據(jù)(2)式得到載荷與層間距的關(guān)系,以此模擬0—13.5 nN 外加載荷.將沿xy平面滑動距離設(shè)置為2.6 ? × 2.6 ?,網(wǎng)格大小設(shè)為0.2 ? × 0.2 ?,第一布里淵區(qū)K網(wǎng)格設(shè)為12 × 12 × 1.計算了1—13.5 nN 外加載荷下Gr/Gr 的勢能面,結(jié)果如圖7 所示.從圖7 可以發(fā)現(xiàn),當外加載荷為1 nN時,勢能面較為平坦,Top 位置勢能最高,Hollow位置勢能最低.隨著載荷升高,勢能最大值先增大后減小.當載荷達到12 nN 時,勢能面發(fā)生了反轉(zhuǎn),即原來高勢能點(Top 位置)變?yōu)榈蛣菽茳c,低勢能點(Hollow 位置)變?yōu)楦邉菽茳c,勢能面趨于平坦,能壘很低.隨著載荷的進一步增大,勢能起伏再次逐漸增大.

圖7 不同載荷下Gr/Gr 的勢能面 (a) 1 nN;(b) 4 nN;(c) 7 nN;(d) 10 nN;(e) 12 nN;(f) 13.5 nN.(a)—(d)中的H/T 對應(yīng)最小/最大勢能,而(e),(f)中的H/T 對應(yīng)最大/最小勢能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.

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

本文計算的摩擦力隨載荷變化的趨勢成功復(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 ?)略有差距,這是因為Sun 等[24]是通過固定部分原子方法計算總能,而本文采用了固定所有原子的方法計算總能,本文計算方法可以節(jié)省計算資源并減少計算時間.Fan 等[53]在實驗中測得的石墨烯的摩擦系數(shù)為0.1—0.2,Wang 等[27]計算結(jié)果顯示,載荷為1—9 nN 時石墨烯界面的摩擦系數(shù)為0.05—0.225.本文計算顯示石墨烯界面摩擦系數(shù)為0.014—0.14,與文獻[53]計算結(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 鍵長分別為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)建了勢能面網(wǎng)格,大小均為5 ? × 5 ?,網(wǎng)格尺寸為0.2 ? × 0.2 ?,計算了0—8 nN 載荷下H—BN/BN 和F—BN/BN 的勢能面及界面平均摩擦力和摩擦系數(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.當載荷小于6 nN 時,H—BN/BN 滑移界面的摩擦系數(shù)小于BN/BN 滑移界面的摩擦系數(shù);當載荷大于6 nN 時,H—BN/BN 滑移界面的摩擦系數(shù)大于BN/BN 滑移界面的摩擦系數(shù).可以明顯地看出F—BN/BN 滑移界面的平均摩擦力和摩擦系數(shù)均小于BN/BN 和H—BN/BN,表明表面氟化顯著降低了摩擦.

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

圖9 BN/BN 和F-BN/BN 的最高和最低勢能態(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.黃色代表電子聚集,藍色代表電子損失;左邊和右邊的圖分別對應(yīng)最高和最低勢能態(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)被廣泛的用來預(yù)測層間滑移的摩擦系數(shù),但該方法仍存在一些局限性.首先,在層間相對滑動的過程中,本文將上下層均設(shè)為剛性,并未考慮層間彎曲的影響.其次,本文通過改變計算體系的堆疊構(gòu)型以模擬層間滑動,因此計算得到的實際結(jié)果本質(zhì)是靜摩擦,而非動摩擦.另外,DFT 計算本身具有局限性,如在計算過程中無法直接引入溫度,因此該方法并未考慮溫度的影響;計算通常是對數(shù)十個原子進行,無法對幾百個原子進行計算.

5 結(jié)論

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