韓強 黃凌燕

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

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

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

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

1 分子動力學模擬

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

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

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

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

圖1 單層石墨烯的原子構型圖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 分子動力學模擬方法

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

勢函數(shù)的選取是分子動力學模擬的關鍵，研究單層石墨烯時選用不考慮層間作用的Tersoff 勢函數(shù)［9］.Tersoff 勢函數(shù)是一種三體勢函數(shù)，它可以很好地模擬C—C 共價鍵的各種特性，包括鍵長、鍵角、鍵能、晶格常數(shù)和鍵的斷裂重組等動態(tài)行為，能夠較真實地反應碳元素所構成的固態(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 粒子周圍其他近鄰原子有關的多體函數(shù)項.

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

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

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

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

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

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

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

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

2 模擬結(jié)果

2.1 拉伸過程的應力－應變曲線

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

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

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

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

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

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

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

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

2.2 單層石墨烯的彈性模量

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

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

用Ex表示x 方向拉伸時對應的扶手椅型石墨烯的彈性模量，用Ey表示y 方向拉伸時對應的鋸齒型石墨烯的彈性模量.通過計算，將不同尺寸單層石墨烯薄膜的彈性模量列于表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，當拉伸應變在0 ～0.03 范圍內(nèi)時，單層石墨烯兩個方向的應力－應變曲線幾乎是線性重合的，且在該階段兩個方向的彈性模量相差不超過6%.于是可以認為在拉伸變形的線彈性階段，石墨烯薄膜是各向同性材料.

2.3 多層石墨烯的彈性模量

與單層石墨烯計算彈性模量值的方法一致，對圖3 中應變在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 為單層石墨烯在拉伸過程中的幾何構型圖.在整個拉伸變形過程中，同一尺寸單層石墨烯兩個方向的破壞形態(tài)是不同的;不同尺寸條件下，單層石墨烯在相同方向的拉伸破壞形態(tài)保持一致.

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

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

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

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

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

3 結(jié)論

文中采用分子動力學的方法，對不同尺寸、不同層數(shù)石墨烯的拉伸力學性能進行研究.所得結(jié)論如下:(1) 石墨烯在整個拉伸變形過程中，兩個方向的應力－應變關系呈現(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］辛浩.石墨層/碳納米管力學性能的研究［D］. 廣州:華南理工大學土木與交通學院，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］韓同偉，賀鵬飛，王建，等. 單層石墨烯薄膜拉伸變形的分子動力學模擬［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.