基于粗粒化分子動力學的自支撐石墨烯鏡面屈曲研究*

續文龍 開玥 張鍇 鄭百林?

1) (同濟大學航空航天與力學學院,上海 200092)

2) (上海工程技術大學數理與統計學院,上海 201620)

1 引言

在石墨烯通過機械剝離方法被制備出之后[1],許多研究者將研究關注點從理論分析轉向實驗觀測探索,從力學[2-6]、熱學[7-10]、光學[11,12]、電磁學[13-17]以及實際應用轉化[18-21]的角度對石墨烯開展全面深入的分析研究,探索這種二維材料獨特的性質.除了對石墨烯物理性質開展研究工作之外,通過電子顯微鏡對石墨烯膜表面觀測的工作也有所增加.Neek-Amal等[22]在通過掃描隧道顯微鏡(scanning tunneling microscopy,STM)對自支撐石墨烯進行實驗觀測時發現,隧穿電流和偏壓在一定范圍內時,石墨烯膜在STM探針的牽引下在面外高度產生突跳,并將這種現象稱為“熱鏡面屈曲”(thermal mirror buckling),如圖1所示.

圖1 熱鏡面屈曲現象的實驗觀測——不同隧穿電流下石墨烯面外高度隨電壓變化曲線,其中黑色曲線對應鏡面屈曲現象[22]Fig.1.Experimental observation of thermal mirror buckling.Constant-current feedback-on,Z(V) data sets on suspended graphene acquired using the labeled setpoint currents.The 4.0-nA curve (black) shows the mirror buckling phenomenon[22].

在實驗觀測到二維材料石墨烯具有這樣的特殊現象之后,不少學者們采用理論分析及仿真計算的方法對其進行研究分析.Schoelz等[23]研究發現石墨烯膜在STM探針靜電-熱作用下產生的面外高度的變化與半自旋伊辛磁場系統具有類似模式.石墨烯產生的不可逆高度突變與鐵磁-反鐵磁相變有所對應.Ruiz-García等[24]考慮一個簡單的自旋-膜模型對具有波紋的石墨烯屈曲以及對其松軟到堅硬狀態的轉變進行定量分析.該模型在非平衡條件下發生的一階相變與實驗中石墨烯膜與STM探針的相互作用類似.Schoelz等[23]利用分子動力學(molecular dynamics,MD)方法對鏡面屈曲現象進行定性分析.之后,Xu等[25]使用宏觀和微觀的方法來探索研究石墨烯熱-機鏡面屈曲現象合適的分析途徑.利用連續介質力學、分子結構力學和分子動力學方法分別對機械載荷和熱載荷作用下的石墨烯膜進行計算,對比結果表明,MD方法計算得到的石墨烯面外高度變化與Neek-Amal等[22]在實驗中獲得的數據及相關仿真結果更加接近.由此可見,MD方法是一種可靠有效地分析石墨烯膜鏡面屈曲現象的途徑.

然而需要注意的是,無論是在發現石墨烯鏡面屈曲的STM實驗中[22],還是其他分析石墨烯跳躍失穩(snap-through buckling)現象實驗[19,26,27]中,其石墨烯的典型尺寸均為微米級或次微米級,而之前MD仿真計算[25]的模型尺寸為納米級.考慮到尺度效應可能對石墨烯的力學[28,29]以及熱學[30,31]性質帶來影響,同時利用石墨烯跳躍失穩現象進行能量收集裝置的設計研究[19]中石墨烯的典型尺寸也為微米級,因此有必要在MD方法之外,使用其他方法研究微米級尺寸石墨烯膜的鏡面屈曲現象以及不同因素對其產生的影響.

本文采用粗粒化分子動力學(coarse-grained molecular dynamics,CGMD)方法構建石墨烯模型,對于扇形截面不同高跨比的石墨烯膜,分別施加機械載荷和熱載荷,研究載荷種類、載荷大小以及模型高跨比對于鏡面屈曲現象的影響.探索分析各個因素對于石墨烯鏡面屈曲現象的影響有助于研究其產生機理,并為能量收集裝置的開發設計提供指導.

2 粗粒化自支撐石墨烯模型及模擬方法

2.1 粗粒化分子動力學方法

CGMD在單層及多層石墨烯的力學性能計算中有著廣泛的應用[32-35].本文采用CGMD方法對石墨烯進行建模,通過調整勢函數參數,實現粗粒化模型與全原子模型(all-atom model)保持相近的模量以實現使用少量粒子等效大尺寸模型的目的.

粗粒化放大倍率可以根據需求進行選擇.本文中倍率選擇4∶1,即1個粗粒化珠子(簡記為珠子)代表4個碳原子,珠子間距離為2.84 ?,如圖2所示.粗粒化模型中珠子質量m可以根據全原子模型中原子質量得到,具體如下:

圖2 CGMD模型示意圖.黑色圓點表示全粒子模型中的碳原子,藍色圓點表示CGMD模型中的珠子.紅色虛線區域表示碳原子與珠子的典型等效方式Fig.2.Schematics of the CGMD model.The black dots indicate carbon atoms in the all-atom model.The blue dots indicate beads in the CGMD model.The region marked with rad dotted lines express a typical equivalent mode with carbon atoms and beads.

其中d0表示珠子之間距離,dcc表示碳-碳鍵之間的距離1.42 ?,mc表示碳原子的相對原子質量12.

CGMD模型中鍵長、鍵角、二面角以及非成鍵相互作用的勢能總和U表示為

其中Eb,Eα,Ed和Enb分別對應鍵長、鍵角、二面角以及非成鍵相互作用對應的勢能;d表示鍵長大小;θ表示鍵角大小;?表示二面角大小;r表示兩個非成鍵珠子之間的距離.這幾種能量形式的具體函數表達式[33]如下所示:

其中D0和α表示Morse勢中與拉伸下非線性行為相關的參數,對于小變形來說,Eb可以近似地寫為Eb(d)≈kb(d-d0)2,其中kb=D0α2;kθ表示鍵角相互作用的彈簧常數;θ0表示平衡鍵角;k?表示二面角的彈簧常數;εLJ表示非成鍵相互作用Lennard-Jones勢阱深度;σLJ表示Lennard-Jones中與兩個非成鍵珠子間距相關的參數.各參數由石墨烯的楊氏模量、剪切模量、彎曲剛度、粗粒化珠子之間鍵的失效應變以及單位表面積的附著力計算得到,具體數值[32]如表1所列.

表1 CGMD化模型力場參數Table 1.Parameters of the CGMD model force field.

通過這些參數建立的粗粒化石墨烯模型可以有效且準確地反映全原子石墨烯模型的受到載荷作用后的響應.根據這些參數構建的粗粒化模型在單層及多層石墨烯的拉伸、剪切、納米壓痕仿真實驗中,在多層石墨烯的層間剪切響應、單軸拉伸破壞響應以及溫度和缺陷對拉伸變形的影響等仿真實驗中得到廣泛應用[32,33].這樣的一組參數同樣適用于處理STM實驗中發現的在面外載荷作用下出現的鏡面屈曲現象.

2.2 CGMD方法驗證

為了確保CGMD方法在研究石墨烯鏡面屈曲現象中的適用性,首先參考實驗設置,對STM實驗中觀測到的鏡面屈曲現象[22]進行仿真計算.方形石墨烯模型兩邊長分別為100.25 nm和100.35 nm,如圖3(a)所示.四邊自由,與實驗中自支撐狀態對應.邊界條件設為周期性,符合實驗中無限大條件.時間步長1 fs,系綜為NVT,使用Nose-Hoover熱浴法控溫.溫度設置為300 K,與實驗中室溫狀態接近.通過deform命令將石墨烯沿兩邊方向分別壓縮長度的2%左右,之后進入弛豫.充分弛豫0.5 ns達到能量平衡之后,可以觀察到石墨烯出現明顯的凸起(或凹陷)變形,如圖3(b)所示.凸起(凹陷)狀態下在中心位置向下(向上)每隔0.3 ns增加6.946×10-4nN的載荷,使石墨烯在受到逐漸增大的載荷的同時保持穩定,并觀測其中心高度隨載荷變化曲線.本文通過直接施加機械載荷與實驗中的靜電載荷對應.

圖3 石墨烯粗粒化模型示意圖,其中黃色區域表示中心位置 (a)初始狀態;(b)中心凹陷狀態Fig.3.Diagrams of CGMD graphene models,where the yellow region is located in the center: (a) Initial state;(b) concave region in the center.

通過計算結果與實驗數據的對比可知,CG MD計算得到的曲線(如圖4)與STM實驗中觀測到的鏡面屈曲現象(圖1中黑色曲線)十分相近,同樣具有3個典型階段,即凹陷回升階段、平穩階段和反向上升階段.分別對應原本凹陷的石墨烯在載荷作用下撓度逐漸增大,直到恢復到零高度初始位置,之后在載荷的牽引下進一步升高的過程.值得注意的是,計算及實驗結果中第2階段出現相對平緩的區域可能是由于石墨烯懸浮狀態的自由邊界使其在受到載荷時釋放面內預壓力,從而由凸起狀態轉變到相對平坦狀態.然后在碳原子相互作用力的影響下使石墨烯在載荷的擾動下保持了一定的面外穩定性.正是因為這一點,石墨烯鏡面屈曲現象載荷-撓度曲線與經典力學中跳躍失穩現象存在差異.之后隨著載荷的進一步增長,碳原子間相互作用力不足以維持,橫向載荷繼續增大引發面外高度的增大.至于機械載荷和靜電載荷作用下面外高度變化形式略有不同,這可能與載荷形式以及采樣頻率相關.

圖4 CGMD計算得到石墨烯中心位置面外高度隨時間的變化Fig.4.Curve of the out-of-plane heights in the center region of the graphene via the time calculated by CGMD.

綜上所述,通過CGMD方法可有效還原STM實驗中觀測到的鏡面屈曲現象,并對仿真結果進行記錄提取,可以證明CGMD方法對于石墨烯鏡面屈曲現象的分析具有有效性.

3 結果與討論

3.1 石墨烯幾何模型及計算設置

目前對于自支撐石墨烯熱機鏡面屈曲現象的研究,石墨烯膜的初始形貌并沒有明確結論.在本文中,選擇扇形截面的開口柱形石墨烯膜為模型,研究機械載荷、溫度載荷及高跨比對石墨烯鏡面屈曲現象的影響.參考其他仿真工作[36]以及彈性穩定性理論[37],截面高跨比λ=h/D,分別為2.5%,5%,7.5%,10%,20%和30%,截面形狀如圖5所示.

圖5 石墨烯膜扇形截面示意圖Fig.5.Diagram of the fan-shaped cross section of a graphene membrane.

截面參數確定之后,其他計算相關條件需要進一步確定.考慮到實驗加載情況,即自支撐石墨烯膜在STM探針作用下同時受到靜電載荷和熱載荷的作用,所以機械載荷和熱載荷對于石墨烯鏡面屈曲的作用都有所研究.參考仿真計算中石墨烯的典型尺寸[22,36,38],對于粗粒化石墨烯模型,機械載荷施加在直徑長度等于 1/10 跨度D的圓形中心位置,而溫度載荷在整個環境中.石墨烯膜對應的跨度D=29.82 nm,長度L=28.78 nm.在開口柱形石墨烯膜兩側位置各有寬l0=0.5 nm的石墨烯帶,將其固定作為邊界條件.加載區域及邊界條件設置如圖6所示.機械載荷大小從1.086×10-3nN到7.599×10-3nN,間隔1.086×10-3nN,最后一組大小為1.086×10-2nN,共8組;溫度載荷大小從3 K開始,之后從50 K增大到500 K,間隔50 K,共11組.具體參數如表2所列.

CGMD仿真計算通過開源代碼LAMMPS[39]實現,可視化通過免費軟件VMD[40]實現.原子間相互作用勢選用2.1節中提及的表達式.時間步長1 fs,使用Nose-Hoover熱浴法控溫.施加機械載荷的情況下,在NVT系綜中保持溫度3 K施加相應載荷2 ns,系統穩定之后,再計算10 ns并獲取數據;施加熱載荷的情況與之類似,在NVT系綜中保持設定溫度2 ns作為弛豫,再計算10 ns并獲取數據.

本文對石墨烯鏡面屈曲的研究中,考慮了載荷大小、載荷形式以及高跨比的影響.對于不同高跨比的石墨烯模型分別施加機械載荷和熱載荷,對比石墨烯膜穩定后截面輪廓圖以及最大撓度-載荷曲線,研究分析各影響因素在熱機鏡面屈曲現象中起到的作用.為了更加直觀且準確地比較計算結果,對計算結果進行無量綱化:=T/250,=f/(4.341×10-3),=w/(493.663).這里表示無量綱機械載荷,wˉ 表示無量綱石墨烯最大撓度.

3.2 扇形截面石墨烯鏡面屈曲現象計算結果

1)機械載荷作用下高跨比和載荷大小對石墨烯鏡面屈曲影響

圖7所示為CGMD計算得到的不同高跨比石墨烯膜最大撓度隨載荷變化曲線.對圖7中各組曲線進行觀察可以看出,撓度隨載荷的變化趨勢可以分為3種類型.對于λ=2.5%的小高跨比情況,石墨烯中心撓度隨載荷的增大只產生小幅增長.對于λ=5%,7.5%,10%,20%的高跨比相對較大的情況,通過觀察可以發現石墨烯中點面外撓度在隨載荷發生一段緩慢增長后,突然產生長距離跳躍,即對應鏡面屈曲現象.最后,對于λ=30%的大高跨比情況,石墨烯中心撓度基本隨載荷的增加線性增長.將計算結果與STM實驗結果(圖1)進行對比可以發現,通過CGMD計算得到相似的面外高度變化趨勢,說明石墨烯的高跨比也會對鏡面屈曲現象產生重要影響.

圖7 扇形截面石墨烯膜的最大撓度隨機械載荷的變化Fig.7.Maximum dimensionless deflections of the fanshaped cross-sections of graphene membranes versus the mechanical force.

四種典型高跨比石墨烯膜在不同機械載荷作用下對應的最大撓度截面輪廓圖,如圖8所示.通過截面輪廓可以更加直觀地看到石墨烯在不同載荷作用下的變形情況.由圖8可知,在λ=2.5%的情況下,石墨烯的穩定性相對較差,即使在小載荷作用下也會發生完全翻轉現象.結合圖7可以發現,高跨比較小的情況下石墨烯中心撓度隨載荷的變化近似呈線性增長趨勢,是處于已經發生翻轉的狀態,這一點與STM實驗中小電流情況下石墨烯高度隨偏壓增大而增大但未發生翻轉現象不同,這可能是兩種計算形式中邊界條件不同造成的.對于λ=5%,10%,20%的情況,石墨烯膜在較小載荷作用下維持相對穩定平坦的截面形狀.載荷逐漸增大,直到石墨烯膜面外高度突然減小,發生鏡面屈曲現象.值得注意的是,發生鏡面屈曲現象對應的載荷隨著高跨比的增大而逐漸增加,這與Timoshenko[37]提出的均勻受壓圓拱兩端鉸支情況下高跨比與臨界屈曲載荷變化趨勢相同.在高跨比小于30%的情況下,臨界屈曲載荷隨著高跨比的增加而增大.對于同一高跨比的石墨烯膜來說,隨著載荷的增大,撓度逐漸增大,直到超過臨界載荷鏡面屈曲現象產生,這同樣符合屈曲理論.

圖8 四種典型高跨比扇形截面石墨烯膜在不同機械載荷作用下,最大撓度截面輪廓圖 (a) λ=2.5%;(b) λ=5%;(c) λ=10%;(d) λ=20%Fig.8.Cross-sectional profiles of maximum dimensionless deflections for the fan-shaped cross-sections of graphene membranes with the four typical depth-span ratio under different mechanical force: (a) λ=2.5%;(b) λ=5%;(c) λ=10%;(d) λ=20%.

2)熱載荷作用下高跨比和載荷大小對石墨烯鏡面屈曲影響

圖9所示為CGMD計算得到的不同高跨比石墨烯膜最大撓度隨溫度的變化.對于λ=2.5%的小高跨比情況,石墨烯中心高度在溫度的影響下幾乎沒有變化.對于高跨比稍大一些的情況,即λ=5%,在載荷增長的過程中,石墨烯中心高度變化存在兩個階段,在Tˉ=0—1.2 過程中,撓度隨載荷的增大而增大;而在Tˉ=1.2—1.4 期間撓度幾乎不變;在Tˉ＞1.4 之后石墨烯中心高度再次開始隨溫度的增加而增大.對于其他較大高跨比的石墨烯而言,撓度隨溫度載荷的增加只產生線性增長,甚至幾乎沒有變化.

圖9 扇形截面石墨烯膜的最大撓度隨溫度載荷的變化Fig.9.Maximum dimensionless deflections of the fanshaped cross-sections of graphene membranes versus the temperature.

下面結合各種溫度載荷作用下石墨烯的截面形狀對溫度和高跨比對鏡面屈曲現象的影響進行進一步分析.四種典型高跨比石墨烯在不同溫度載荷作用下對應的最大撓度截面輪廓圖,如圖10所示.圖11所示為通過MD計算的在溫度3000 K的環境下不同壓縮率的圓形石墨烯膜中心高度隨時間的變化[36].根據換算可以得到圖中壓縮率ε=0.13%,0.67%,0.78%分別對應于高跨比λ=2.3%,5.3%,5.7%.

圖10 四種典型高跨比扇形截面石墨烯膜在不同溫度載荷作用下,最大撓度截面輪廓圖 (a) λ=2.5%;(b) λ=5%;(c) λ=10%;(d) λ=20%Fig.10.Cross-sectional profiles of maximum dimensionless deflections for the fan-shaped cross-sections of graphene membranes with the four typical depth-span ratio at different temperature: (a) λ=2.5%;(b) λ=5%;(c) λ=10%;(d) λ=20%.

圖11 3000 K下應變對波紋的動態影響[36]Fig.11.Role of strain on ripple dynamics at 3000 K[36].

由圖10可知,在λ=2.5%的小跨高比情況下,無論溫度載荷的大小,石墨烯始終處于可翻轉狀態.這種圍繞初始高度反復進行曲率翻轉的行為與相近高跨比(ε=0.13%)下石墨烯在溫度載荷作用下的行為相近.對于λ=5%情況,溫度載荷的增加使石墨烯中心高度降低,并在Tˉ=2 時發生完全翻轉.由此可知,CGMD計算得到的石墨烯曲率翻轉現象與ε=0.67%石墨烯在高溫下發生的翻轉行為有所對應.對于λ=10%,20%的較大高跨比情況,石墨烯構型在溫度載荷的作用下基本保證穩定狀態,沒有產生明顯的翻轉現象,同樣與ε=0.78%較大壓縮率下石墨烯的中心高度變化相近.通過將CGMD與MD對不同高跨比石墨烯施加溫度載荷的計算結果進行對比后發現,兩者得到的石墨烯變形行為相近.

一般情況下,溫度的升高有助于石墨烯發生鏡面屈曲現象.升高溫度即原子運動加劇的宏觀表征.但由于原子間相互作用,碳原子并不會向隨機方向運動,而是在周圍原子的牽引下受到一定約束的運動.運動速度的加快以及運動方向上一定的約束綜合導致升高溫度有助于鏡面屈曲現象的發生.與此同時,石墨烯高跨比也會對鏡面屈曲現象的產生造成影響.對于λ=2.5%的情況,即使溫度較低也會使石墨烯完全翻轉,且最大撓度隨著溫度的升高緩慢增加.在λ=5%時,石墨烯中心高度隨溫度的升高而逐漸降低,當溫度足夠高達到Tˉ=2的情況下,石墨烯發生完全翻轉.對于λ進一步增大的情況,即使溫度升高也只能觀察到石墨烯膜面外高度下降,突跳現象不會出現.

4 結論

本文采用CGMD方法,對扇形截面不同高跨比石墨烯膜分別施加機械載荷和溫度載荷,研究高跨比、載荷大小以及載荷類型對于鏡面屈曲現象的影響.首先通過CGMD仿真計算結果與STM實驗觀測結果對比可知,CGMD方法可以有效地對鏡面屈曲現象進行計算分析,證明其有效性.

隨后在對各個影響因素的分析中可以發現,機械載荷作用下,在小高跨比λ=2.5%的情況下,即使在小載荷作用下也會發生鏡面屈曲現象.對于較大高跨比λ=5%,10%,20%的情況,在小載荷作用下石墨烯可能產生局部坍塌,當載荷增大到臨界值時,石墨烯面外高度突然發生長距離跳躍,即產生鏡面屈曲現象.由此可見,在機械載荷作用下,對于高跨比各異的石墨烯都能夠產生鏡面屈曲現象,且對應的臨界載荷隨著高跨比的增加而增大.

對于溫度載荷來說,在小高跨比λ=2.5%的情況下,無論溫度載荷大小,石墨烯始終處于可翻轉狀態.對于λ=5%的情況,在Tˉ=2 的較高溫度作用下能夠發生完全翻轉現象.對于更大高跨比的石墨烯來說,隨著溫度的升高只觀察到中心高度的下降,完全翻轉的現象并沒出現.分析可知,高跨比對溫度載荷下石墨烯的高度變化影響較大,小高跨比的情況下更容易發生翻轉,對于大高跨比的情況,溫度載荷的增加使中心高度的變化增大,可能導致翻轉情況的發生,但當高跨比進一步上升時,溫度的升高帶來中心高度的下降,但是不會導致完全翻轉的發生.

在載荷類型、載荷大小以及高跨比之外,還存在截面形狀、邊界條件等更多條件可能對石墨烯鏡面屈曲現象產生影響,之后的工作可以展開對更多因素的分析工作.通過對石墨烯膜鏡面屈曲現象更加全面細致的研究,清楚各因素對其的影響作用,有助于基于石墨烯膜的能量收集系統的設計和優化工作.