石墨烯薄膜拉伸性能的分子動力學(xué)模擬*

韓強 黃凌燕

( 華南理工大學(xué) 土木與交通學(xué)院，廣東 廣州510640)

自從碳納米管［1］( CNTs) 被發(fā)現(xiàn)以來，由于其在力學(xué)、化學(xué)、熱力學(xué)和電學(xué)等方面具有的獨特性質(zhì)，使得納米研究成為當(dāng)今世界的熱點問題之一. 近年來，隨著單層石墨烯［2］的制備成功，引發(fā)了新的碳質(zhì)材料的研究熱潮.單層石墨烯是以sp2軌道雜化的碳原子形成的蜂窩狀六角平面晶體，厚度僅為0.335 nm，其中C—C 鍵鍵長約為0.142 nm，為目前世界上存在的最薄的材料. 理論計算［3］表明，單層石墨烯的彈性模量在700 ～1 000 GPa 之間，是普通鋼材的5 倍.

關(guān)于石墨烯基礎(chǔ)性能的研究是納米材料的研究熱點之一，辛浩［4］提出采用三層正交異性單向纖維的疊層作為石墨烯的等效模型，從理論上計算了石墨烯的彈性模量，為0.97 TPa.Gao 等［5］采用量子力學(xué)和分子動力學(xué)的方法得到鋸齒型和扶手椅型石墨烯的楊氏模量，分別為0.6 TPa 和1.1 TPa. Ni 等［6］研究了室溫下石墨烯彈性模量的尺寸效應(yīng)，認為當(dāng)尺寸變化時，鋸齒型石墨烯的彈性模量值始終高于扶手椅型石墨烯的彈性模量值. 韓同偉等［7］通過分子動力學(xué)方法對單層石墨烯薄膜的拉伸變形進行了研究.現(xiàn)有的研究中，有關(guān)尺寸效應(yīng)和層數(shù)效應(yīng)的文章并不多見，還沒有形成統(tǒng)一的結(jié)論.

文中基于分子動力學(xué)的方法，模擬不同尺寸單層石墨烯在兩個方向上的拉伸變形，并研究拉伸過程中的應(yīng)力－應(yīng)變關(guān)系、力學(xué)性能以及尺寸變化對石墨烯拉伸性能的影響.在尺寸相同的條件下，通過改變石墨烯的層數(shù)，研究層數(shù)與石墨烯拉伸力學(xué)性能之間的關(guān)系.

1 分子動力學(xué)模擬

1.1 單層、多層石墨烯的物理模型

文中采用的單層石墨烯薄膜尺寸列于表1 中.其中，Lx表示單層石墨烯x 方向的長度，Ly表示單層石墨烯y 方向的長度.

圖1 所示為單層石墨烯的原子構(gòu)型圖.可以看出:當(dāng)對薄膜進行x 方向拉伸時，相當(dāng)于分析扶手椅型石墨烯的拉伸性能;當(dāng)對薄膜進行y 方向拉伸時，相當(dāng)于分析鋸齒型石墨烯的拉伸性能.

表1 單層石墨烯薄膜的尺寸Table 1 Size of monolayer graphene sheet

圖1 單層石墨烯的原子構(gòu)型圖Fig.1 Atomic model of monolayer graphene sheet

多層石墨烯是由單層石墨烯在z 方向進行Bernal( AB)［8］堆疊后得到的，多層石墨烯的層間距為0.335 nm.

將表1 中模型編號為6 的單層石墨烯薄膜在z方向堆疊，得文中模擬的多層石墨烯模型尺寸，列于表2.

表2 多層石墨烯模型尺寸Table 2 Size of multilayer graphene sheet model

1.2 分子動力學(xué)模擬方法

分子動力學(xué)模擬的基本原理是通過原子間的相互作用勢，求出每個原子所受到的力，在選定的時間步長、邊界條件和初始條件下，對有限數(shù)目的分子( 原子) 建立其牛頓動力學(xué)方程組，用數(shù)值方法求解，得到這些原子的運動軌跡和速度分布，然后對足夠長時間的結(jié)果求統(tǒng)計平均，從而得到所需要的宏觀物理量和力學(xué)量.

勢函數(shù)的選取是分子動力學(xué)模擬的關(guān)鍵，研究單層石墨烯時選用不考慮層間作用的Tersoff 勢函數(shù)［9］.Tersoff 勢函數(shù)是一種三體勢函數(shù)，它可以很好地模擬C—C 共價鍵的各種特性，包括鍵長、鍵角、鍵能、晶格常數(shù)和鍵的斷裂重組等動態(tài)行為，能夠較真實地反應(yīng)碳元素所構(gòu)成的固態(tài)材料的物理性質(zhì).

原子間相互作用Tersoff 勢函數(shù)為

式中，E 為體系的總能量; Ei為單個原子能量; Eij為i、j 原子之間的成鍵能量，分別為對勢的吸引項和排斥項，fC是光滑截斷函數(shù)，bij為鍵序函數(shù)，rij是i、j 原子之間的鍵長.

Tersoff 勢函數(shù)形式上是一個二體勢，實際上是一個多體勢，因為系數(shù)bij并非一個常數(shù)，而是一個依賴于i、j 原子位置并與i 粒子周圍其他近鄰原子有關(guān)的多體函數(shù)項.

研究多層石墨烯時，考慮到層間作用力的影響，需要在Tersoff 勢函數(shù)的基礎(chǔ)上考慮Lennard-Jones( L-J) 勢函數(shù).

L-J 勢函數(shù)形式如下:

式中:右邊括號內(nèi)第1 項是排斥項，第2 項為吸引項;參數(shù)r 是連接i 和j 原子的位置矢量的模;ε 和l分別表示能量標(biāo)度參數(shù)和碰撞直徑參數(shù)，ε 的意義等同于離解能，l 則可以理解為兩個原子間相互作用勢為零時的原子間距.

目前，運用最為廣泛的積分算法是Verlet 算法［10］，但是速度處理不是很方便，可能導(dǎo)致不同物理量計算精度不同. 在LAMMPS 中，默認的積分算法是Velocity-Verlet，它可以解決Verlet 算法的不足，具體形式如下:

式中，t 為時間，δt 為時間變量;r 為原子位置，v 為原子的速度，m 為原子質(zhì)量，F(xiàn) 表示原子所處的力場.這種算法可以同時獲得相同精度的原子位置和速度，在每步積分中只要存儲一個時刻的狀態(tài)變量，模擬穩(wěn)定性較好，同時允許較大的時間步長，因此在分子動力學(xué)積分運算中得到廣泛的運用.

模擬過程中，選取C 原子質(zhì)量為12. 控制x、y方向為自由邊界條件，z 方向為周期性邊界條件.采用Nose-Hoover 熱浴法［11-12］進行溫度調(diào)節(jié)，控制溫度為0.01 K，設(shè)計時間步長為1 fs.

當(dāng)進行x 方向拉伸模擬時，固定薄膜左端碳原子，先進行充分的弛豫，穩(wěn)定后對右端碳原子施加x方向拉伸應(yīng)變載荷，每步應(yīng)變量為0.001，弛豫步長為3000 步.持續(xù)加載，直至薄膜被拉斷.

當(dāng)進行y 方向拉伸模擬時，固定薄膜下端碳原子，對上端原子施加y 方向的拉伸應(yīng)變載荷，其他步驟與x 方向拉伸相同.

2 模擬結(jié)果

2.1 拉伸過程的應(yīng)力－應(yīng)變曲線

材料的應(yīng)力－應(yīng)變關(guān)系是研究材料特性的重點之一.對表1 中列出的不同尺寸單層石墨烯和表2中列出的不同厚度多層石墨烯進行拉伸模擬，得到相應(yīng)的應(yīng)力－應(yīng)變關(guān)系.

在計算石墨烯拉伸變形性能時，需要采用系統(tǒng)的平均應(yīng)力.由離散原子構(gòu)成的納米系統(tǒng)中應(yīng)力的概念與連續(xù)介質(zhì)中應(yīng)力的概念有所不同，Jin 等［13］在研究納米管力學(xué)性能時給出了納米系統(tǒng)下平均應(yīng)力的計算公式:

式中:σαβ表示在笛卡爾坐標(biāo)下系統(tǒng)原子水平的平均應(yīng)力;分別表示原子i 與原子j 之間的相互作用力和距離; V0為模型初始狀態(tài)下的體積. 對于單層石墨烯薄膜有V0=LxLyh，其中，Lx和Ly為單層石墨烯的長度和寬度值，h 為單層石墨烯的厚度，取為0.335 nm［14-16］.

圖2 為表1 中模型編號為6 的單層石墨烯薄膜在兩個方向上的應(yīng)力－ 應(yīng)變曲線. 圖中，x 方向的拉伸曲線對應(yīng)扶手椅型石墨烯的拉伸性能，y 方向的拉伸曲線對應(yīng)鋸齒型石墨烯的拉伸性能.

圖3 為1 層、2 層和4 層石墨烯拉伸過程中的應(yīng)力－應(yīng)變曲線.為了便于比較層數(shù)變化對應(yīng)力－應(yīng)變關(guān)系的影響，將1、2、4 層的拉伸應(yīng)力－應(yīng)變關(guān)系曲線畫在同一幅圖上.

圖2 和圖3 表明，多層石墨烯和單層石墨烯具有相同的應(yīng)力－應(yīng)變曲線形態(tài)，且尺寸變化和層數(shù)變化對曲線形態(tài)幾乎不產(chǎn)生影響. 石墨烯在兩個方向的拉伸過程中均經(jīng)歷了彈性變形、屈服階段、強化階段和局部變形4個階段，但兩個方向所表現(xiàn)出來的力學(xué)性能相差較大，x 方向( 扶手椅型) 的拉伸屈服強度、斷裂強度以及極限應(yīng)變值均高于y 方向( 鋸齒型) ，表明了石墨烯各向異性的特征.

圖2 單層石墨烯拉伸的應(yīng)力－應(yīng)變曲線Fig. 2 Tensile stress-strain curves of monolayer graphene sheet

圖3 多層石墨烯拉伸的應(yīng)力－應(yīng)變曲線Fig.3 Tensile stress-strain curves of multilayer graphene sheets

2.2 單層石墨烯的彈性模量

這里，僅對圖2 中應(yīng)變在0 ～0.03 段的曲線進行最小二乘法擬合，計算單層石墨烯的彈性模量如圖4 所示為擬合范圍內(nèi)的應(yīng)力－應(yīng)變關(guān)系.

圖4 擬合范圍內(nèi)的應(yīng)力－應(yīng)變曲線Fig.4 Tensile stress-strain curves in the fitting range

用Ex表示x 方向拉伸時對應(yīng)的扶手椅型石墨烯的彈性模量，用Ey表示y 方向拉伸時對應(yīng)的鋸齒型石墨烯的彈性模量.通過計算，將不同尺寸單層石墨烯薄膜的彈性模量列于表3.

表3 單層石墨烯的彈性模量Table 3 Elastic modulus of monolayer graphene sheet

從表3 可以看出，在單層石墨烯薄膜尺寸由小增大的過程中，其彈性模量并未發(fā)生太大的改變，Ex總是在1078.02 GPa 附近，而Ey總是在1041.53 GPa附近，且Ex和Ey兩者的差距始終小于6%. 可見，尺寸變化對單層石墨烯彈性模量的影響很小，且在尺寸變化過程中Ex始終大于Ey.

結(jié)合圖4，當(dāng)拉伸應(yīng)變在0 ～0.03 范圍內(nèi)時，單層石墨烯兩個方向的應(yīng)力－應(yīng)變曲線幾乎是線性重合的，且在該階段兩個方向的彈性模量相差不超過6%.于是可以認為在拉伸變形的線彈性階段，石墨烯薄膜是各向同性材料.

2.3 多層石墨烯的彈性模量

與單層石墨烯計算彈性模量值的方法一致，對圖3 中應(yīng)變在0.00 ～0.03 段的曲線進行最小二乘法擬合，將計算結(jié)果列于表4 中.

表4 不同層數(shù)石墨烯的彈性模量值Table 4 Elastic modulus of graphene sheets with different layers

從表4 可以看出，多層石墨烯在兩個方向上的彈性模量與單層石墨烯的接近，均在1 050 GPa 附近.單層石墨烯兩個方向彈性模量的差距為1.35%，而多層石墨烯的彈性模量差距均小于0.4%.可見，石墨烯層數(shù)的增加有助于減小兩個方向彈性模量的差距，在拉伸變形的線彈性階段，多層石墨烯比單層石墨烯表現(xiàn)出更為強烈的各向同性.

2.4 拉伸變形破壞形態(tài)分析

圖5 為單層石墨烯在拉伸過程中的幾何構(gòu)型圖.在整個拉伸變形過程中，同一尺寸單層石墨烯兩個方向的破壞形態(tài)是不同的;不同尺寸條件下，單層石墨烯在相同方向的拉伸破壞形態(tài)保持一致.

圖5 模型編號為6 的單層石墨烯在兩個方向上的拉伸變形Fig.5 Tensile deformation in both directions of the NO.6 monolayer graphene sheet

在x 方向拉伸變形過程中，當(dāng)石墨烯處于彈性變形階段時，模型在軸向載荷作用下作簡單的伸長變形，碳環(huán)保持為六邊形; 隨著載荷的增大，邊緣六邊形碳環(huán)開始變得不規(guī)則起來，當(dāng)載荷達到一定值時，石墨烯一側(cè)邊緣碳環(huán)出現(xiàn)明顯的破壞跡象，同時形成一條近似成45°角的斜向破壞線，迅速向石墨烯內(nèi)部延伸，直至石墨烯被拉斷.

在y 方向拉伸變形過程中，線彈性階段的拉伸變形與x 方向拉伸變形類似;隨著載荷的增大，邊緣六邊形碳環(huán)變得不規(guī)則起來，當(dāng)載荷達到一定值時，石墨烯兩側(cè)邊緣有一部分碳環(huán)出現(xiàn)明顯的破壞跡象，隨著載荷的增加，這種破壞形式越來越明顯，并對稱地向石墨烯內(nèi)部延伸，直至石墨烯被拉斷，形成圖中所示的破壞形態(tài).

圖6 為雙層石墨烯在拉伸過程中的幾何構(gòu)型圖.通過分析發(fā)現(xiàn):多層石墨烯的拉伸破壞形態(tài)與單層石墨烯一致;多層石墨烯在拉伸過程中，各單層的變形破壞形態(tài)也是一致的.

圖6 雙層石墨烯在兩個方向上的拉伸變形Fig.6 Tensile deformation in both directions of doublelayer graphene sheets

3 結(jié)論

文中采用分子動力學(xué)的方法，對不同尺寸、不同層數(shù)石墨烯的拉伸力學(xué)性能進行研究.所得結(jié)論如下:(1) 石墨烯在整個拉伸變形過程中，兩個方向的應(yīng)力－應(yīng)變關(guān)系呈現(xiàn)不同的規(guī)律，表明石墨烯是各向異性的，且這一特性不受尺寸和層數(shù)的影響; (2) 對單層石墨烯，x 方向( 扶手椅型) 的彈性模量在1078.02 GPa附近，y 方向( 鋸齒型) 的彈性模量在1041.53 GPa附近，在拉伸線彈性變形階段，單層石墨烯是各向同性的;對多層石墨烯，在拉伸線彈性變形范圍內(nèi)，各向同性更為明顯; (3) 在拉伸載荷作用下，扶手椅型石墨烯的破壞從一側(cè)邊緣開始，并沿45°方向向石墨烯內(nèi)部延伸，鋸齒型石墨烯的破壞從兩側(cè)邊緣對稱地向內(nèi)部延伸，隨著C—C 鍵的斷裂，石墨烯被拉斷.

［1］Iijima S. Helical microtubules of graphitic carbon ［J］.Nature，1991，354:56-58.

［2］Novoselov K S，Geim A K，Morozov S V，et al. Electric field effect in atomically thin carbon films［J］. Science，2004，306(22) :666-669.

［3］Zhou J，Huang R.Internal lattice relaxation of single-layer graphene under in-plane deformation［J］. Journal of the Mechanics and Physics of Solids，2008，56( 4) : 1609-1623.

［4］辛浩.石墨層/碳納米管力學(xué)性能的研究［D］. 廣州:華南理工大學(xué)土木與交通學(xué)院，2010.

［5］Gao Y W，Hao P. Mechanical properties of monolayer grapheme under tensile and compressive loading ［J］.Physica E，2009，41(8) :1561-1566.

［6］Ni Z H，Bu H，Zou M，et al. Anisotropic mechanical properties of graphene sheets from molecular dynamics［J］.Physica B，2010，405(5) :1301-1306.

［7］韓同偉，賀鵬飛，王建，等. 單層石墨烯薄膜拉伸變形的分子動力學(xué)模擬［J］. 新型碳材料，2010，25( 4) :261-266.Han Tong-wei，He Peng-fei，Wang Jian，et al. Molecular dynamics simulation of monolayer graphene under tensile deformation［J］. New Carbon Materials，2010，25( 4) :261-266.

［8］Bernal J D. The structure of graphite［J］. Proc Roy Soc A，1924，106(740) :749-773.

［9］Tersoff J. Modeling solid-state chemistry: interatomic potentials for multicomponent systems ［J］. Phys Rev B，1989，39(8) :5566-5568.

［10］Allen M P，Tildesley D J.Computer simulation of liquids［M］.Oxford:Clarendon Press，1989:385.

［11］Nose S.A unified formulation of the constant temperature molecular dynamics methods［J］.The Journal of Chemical Physics，1984，81(1) :511-519.

［12］Hoover W G. Canonical dynamics: equilibrium phasespace distributions ［J］. Physical Review A，1985，31(3) :1695-1702.

［13］Jin Y，Yuan F G.Simulation of elastic properties of single-walled carbon nanotubes ［J］. Composites Science and Technology，2003，63(11) :1507-1515.

［14］Lee C G，Wei X D，Kysar W J，et al.Measurement of the elastic properties and intrinsic strength of monolayer graphene［J］.Science，2008，321(5887) :385-388.

［15］Ni Z H，Wang H M，Kasim J，et al. Graphene thickness determination using reflection and contrast spectroscopy［J］.Nano Letters，2007，7(9) :2758-2763.

［16］Dresselhaus M S，Dresselhaus G，Eklund P C.Science of fullerenes and carbon nanotubes［M］. San Diego: Academic Press，1996:965.