韓強(qiáng) 王彩紅 姚小虎

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

碳納米管(CNTS)［1］具有優(yōu)異的力學(xué)、電學(xué)、熱學(xué)等性質(zhì)，在實(shí)驗(yàn)和理論研究中已廣泛用于復(fù)合材料，能夠極大改善材料的性能. 碳納米管優(yōu)異的力學(xué)性質(zhì)與其結(jié)構(gòu)特性密切相關(guān)，然而在生產(chǎn)過程中不可避免地會(huì)出現(xiàn)各種缺陷，特別是Stone-Wales(SW)及空位缺陷［2］，這些缺陷的存在會(huì)削弱碳納米管的拉壓等性能.碳納米管在軸壓荷載作用下發(fā)生屈曲行為［3］，辛浩等［4］采用Morse 勢函數(shù)，研究了含單雙原子及SW 缺陷碳納米管的軸壓屈曲性能，認(rèn)為各類缺陷均會(huì)顯著降低納米管的屈曲性能.Hirai 等［5］研究了沿軸向分布的空位缺陷對(duì)碳納米管壓縮性能的影響，發(fā)現(xiàn)在缺陷處存在應(yīng)力集中現(xiàn)象，軸壓屈曲荷載明顯減小. Huq 等［6］采用分子力學(xué)模擬，研究了SW 缺陷對(duì)單層碳納米管軸壓行為的影響.Zhang 等［7］運(yùn)用分子動(dòng)力學(xué)模擬，研究了SW 及空位缺陷對(duì)單層碳納米管軸壓性能的影響，結(jié)果表明，缺陷大大削弱了碳納米管的臨界屈曲荷載.Ang 等［8］采用Compass 勢，研究了含1 －3個(gè)空位缺陷的單層碳納米管的軸壓屈曲行為. Kulathunga 等［9］采用Compass 勢函數(shù)，研究了含對(duì)稱與非對(duì)稱空位缺陷單層碳納米管的軸壓行為. Chowdhury等［10］采用Brenner 勢函數(shù)，研究了完善及含缺陷碳納米管的拉壓扭轉(zhuǎn)等力學(xué)性能，認(rèn)為空位缺陷對(duì)碳納米管的影響高于SW 缺陷.Vijayaraghavan 等［11］采用Airebo 勢函數(shù)，研究了含缺陷單層碳納米管軸壓行為.Kinoshita 等［12］采用Airebo 勢函數(shù)，分析了含SW 缺陷碳納米管的軸壓性能. Eftekhari 等［13］采用Tersoff 勢函數(shù)，分別研究了含空位及SW 缺陷的單層與雙層碳納米管的軸壓性能，研究表明，SW 缺陷對(duì)碳納米管屈曲性能的影響較空位缺陷顯著. 可見，已有的工作多是對(duì)含單雙空位及SW 缺陷的研究，關(guān)于含裂紋碳納米管軸壓性能的研究工作很少.

文中運(yùn)用分子動(dòng)力學(xué)方法，對(duì)含裂紋扶手椅型碳納米管的軸壓性能進(jìn)行了模擬分析，重點(diǎn)考慮了裂紋尺寸及管長變化對(duì)碳納米管軸壓屈曲性能的影響，詳細(xì)分析了含不同裂紋尺寸碳納米管及不同長度納米管的軸壓變形特性.

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

1.1 碳納米管的物理模型

文中以扶手椅型(7，7)及(10，10)碳納米管為研究對(duì)象，模擬過程中，碳納米管的長度變化范圍為6 ～30 nm. 在碳納米管中心附近區(qū)域，沿著與軸向成30°及90°方向去掉一定數(shù)目的碳原子，構(gòu)建出與軸向成不同傾角θ 的裂紋，設(shè)定初始裂紋尺寸為CL，在模擬過程中分別取為0.000、0.710、0.994及1.562 nm，其中裂紋尺寸為0.000 nm 表示完善的碳納米管.計(jì)算模型結(jié)構(gòu)如圖1 所示.

圖1 長度為6 nm 的(7，7)碳納米管的結(jié)構(gòu)Fig.1 Structures of (7，7)carbon nanotube with a length of 6nm

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

分子動(dòng)力學(xué)計(jì)算的一個(gè)關(guān)鍵問題是原子勢函數(shù)的選取，文中選用Airebo 勢函數(shù)［14］，這是一種改進(jìn)的Rebo 勢函數(shù)，能夠較為真實(shí)地反映碳元素所構(gòu)成的固態(tài)材料的物理性質(zhì).

模擬過程中，控制X、Y、Z 方向?yàn)樽杂蛇吔鐥l件，沿著碳納米管軸向?yàn)榧虞d方向，固定一端碳原子，對(duì)另一端碳原子施加恒速壓縮位移載荷. 采用Nose-Hoover 熱浴法［15-16］進(jìn)行溫度調(diào)節(jié)，使系統(tǒng)保持在特定的溫度下，采用Velocity-Verlet 算法［17］進(jìn)行分子動(dòng)力學(xué)計(jì)算.加載之前先對(duì)碳納米管進(jìn)行充分無約束弛豫，使整個(gè)系統(tǒng)處于能量最低的平衡狀態(tài).對(duì)弛豫過的碳納米管施加軸壓位移載荷，每步應(yīng)變量為0.000 5，然后進(jìn)行充分的弛豫，重復(fù)此加載、弛豫過程，直至碳納米管壓縮破壞.

2 模擬結(jié)果分析

2.1 裂紋尺寸對(duì)碳納米管軸壓性能的影響

為研究裂紋尺寸變化對(duì)碳納米管軸壓屈曲性能及變形特性的影響，選取長度為6 nm、初始裂紋尺寸分別為0.0、0.71、0.994 及1.562 nm 的(7，7)與(10，10)碳納米管進(jìn)行模擬分析，模擬過程中控制溫度為0.01 K，以避免熱激活的復(fù)雜影響，時(shí)間步長取1 fs.

2.1.1 垂直于軸線方向的裂紋

采用分子動(dòng)力學(xué)方法，模擬含有垂直于軸向裂紋的(7，7)、(10，10)扶手椅型碳納米管在軸壓載荷作用下的屈曲及后屈曲變形行為.為了進(jìn)行對(duì)比分析，首先模擬了完善(7，7)管的軸壓屈曲變形性能，其納米管直徑約0.95 nm.

文中與辛浩［4］及Yakobson 等［3］所選取的勢函數(shù)、分子動(dòng)力學(xué)計(jì)算方法及模型尺寸不完全相同，將他們的結(jié)果與文中的模擬結(jié)果進(jìn)行對(duì)比分析可知(如表1 所示)，文中得到的(7，7)管軸壓臨界屈曲應(yīng)變與他們的結(jié)果很接近，可見Airebo 勢函數(shù)可以較好地描述碳納米管中的C—C 鍵作用.

表1 完善(7，7)碳納米管在軸壓變形過程中關(guān)鍵點(diǎn)的應(yīng)變值Table 1 Critical strains of perfect (7，7)CNTs during axial compression

由圖2 可知，對(duì)于完善碳納米管，隨著壓縮進(jìn)行，其軸向應(yīng)力呈非線性變化.圖2 中a、b、c、d 四點(diǎn)與圖3(a)、3(b)、3(c)、3(d)一一對(duì)應(yīng)，當(dāng)軸壓應(yīng)變?yōu)?.1347(圖中d 點(diǎn))時(shí)，碳納米管徹底壓縮破壞.

圖2 完善(7，7)碳納米管的軸壓應(yīng)力－應(yīng)變曲線Fig.2 Stress-strain curve of perfect (7，7)carbon nanotube under axial compression

圖3 完善(7，7)碳納米管的軸壓變形Fig.3 Configurations of perfect (7，7)carbon nanotube under axial compression

圖4 θ=90°、裂紋尺寸為0.71nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變曲線Fig.4 Stress-strain curve of (7，7)tube with the crack length of 0.71 nm under axial compression when θ=90°

圖5 θ=90°、裂紋尺寸為0.71 nm 的(7，7)管的軸壓變形Fig.5 Configurations of (7，7)tube with the crack length of 0.71 nm under axial compression when θ=90°

圖4 和5 給出了裂紋尺寸為0.71 nm 的(7，7)碳納米管在軸壓過程中的應(yīng)力－應(yīng)變曲線及變形形態(tài). 由圖5可知，在軸壓過程中，含裂紋(7，7)管的變形形態(tài)與完善(7，7)管有明顯區(qū)別，對(duì)于完善碳納米管(如圖3 所示)，隨著壓縮進(jìn)行，管壁某處發(fā)生凹陷后，凹陷會(huì)沿著軸向進(jìn)行傳播，管壁發(fā)生凹陷的地方可能會(huì)回彈.然而對(duì)于裂紋尺寸為0.71 nm的納米管，在整個(gè)軸壓過程中，裂紋處的凹陷沒有發(fā)生回彈，這是由于裂紋處缺少相應(yīng)的支撐作用.

裂紋附近區(qū)域受力不均勻，在壓縮至a 點(diǎn)時(shí)，裂紋區(qū)域出現(xiàn)局部凹陷變形.與完善(7，7)管相比，裂紋尺寸為0. 71 nm 時(shí)，其臨界屈曲強(qiáng)度下降了40.4%，當(dāng)壓縮至b 點(diǎn)時(shí)，與裂紋部位相對(duì)應(yīng)的另一側(cè)管壁出現(xiàn)凹陷. 對(duì)比完善管及含裂紋(7，7)管的應(yīng)力－應(yīng)變曲線可知，在b、c、d 處兩種管內(nèi)的應(yīng)力差別較小，可見裂紋尺寸為0.71 nm 時(shí)，納米管屈曲后軸向承載力的變化不大.碳納米管最終在裂紋處彎折破壞.

圖6、7 給出了裂紋尺寸為0.994 nm 的(7，7)管的軸壓性能. 由圖7 可知，裂紋尺寸為0. 994 與0.71 nm的(7，7)管的軸壓過程相似，在裂紋附近處均出現(xiàn)凹陷并最終在裂紋處壓縮破壞，在整個(gè)軸壓過程中，裂紋處的凹陷沒有發(fā)生回彈.由圖6可知，碳納米管在軸壓過程中荷載變化較為頻繁，在未達(dá)到臨界屈曲強(qiáng)度之前出現(xiàn)a、b 兩個(gè)波動(dòng)點(diǎn)，兩點(diǎn)對(duì)應(yīng)的變形結(jié)構(gòu)表明含裂紋一側(cè)的管壁出現(xiàn)凹陷.壓縮至c 點(diǎn)時(shí)，管壁另一側(cè)也出現(xiàn)凹陷，此時(shí)納米管達(dá)到軸壓屈曲強(qiáng)度. 與完善管相比，其屈曲強(qiáng)度降低了45.4%.

圖6 θ =90°、裂紋尺寸為0. 994 nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變曲線Fig.6 Stress-strain curve of (7，7)tube with the crack length of 0.994 nm under axial compression when θ=90°

圖7 θ=90°、裂紋尺寸為0.994 nm 的(7，7)管的軸壓變形Fig.7 Configurations of (7，7)tube with the crack length of 0.994 nm under axial compression when θ=90°

對(duì)比裂紋尺寸為0.994 與0.71 nm 的(7，7)管的變形可知，隨著裂紋尺寸增加，裂紋處的凹陷略有差別，對(duì)于裂紋尺寸較大的納米管，其凹陷處裂紋上下邊緣呈現(xiàn)明顯的錯(cuò)位，這是由于裂紋處缺少了沿軸向的支撐作用，裂紋尺寸增大時(shí)，這種缺失越明顯.

圖8 給出了(7，7)、(10，10)管的臨界屈曲強(qiáng)度隨裂紋尺寸的變化，隨著裂紋尺寸增加，碳納米管的屈曲強(qiáng)度均減小. 在裂紋尺寸由0 增加到1.562 nm的過程中，(7，7)與(10，10)管的軸向臨界屈曲強(qiáng)度分別降低了57.05%及49.98%.

圖8 θ=90°時(shí)臨界屈曲強(qiáng)度隨裂紋尺寸的變化Fig.8 Critical buckling strength vs.crack length when θ=90°

2.1.2 與軸線方向成30°角的裂紋

由上節(jié)可知，碳納米管中部區(qū)域垂直于軸線的裂紋大大減弱了其軸向承載力，為了進(jìn)行更加全面的分析，本節(jié)模擬分析了裂紋傾角為30°、管長為6 nm，含不同裂紋尺寸的(7，7)、(10，10)扶手椅型碳納米管的軸壓屈曲性能.

圖9、10 給出了裂紋傾角為30°、尺寸為0.71 nm的(7，7)碳納米管的軸壓性能. 裂紋附近區(qū)域的管壁受力不均勻，碳納米管壓縮至a 點(diǎn)時(shí)，在裂紋附近的管壁首先發(fā)生凹陷變形，與完善管相比，其臨界屈曲強(qiáng)度降低了40.6%. 隨著荷載增加，管壁另一側(cè)也出現(xiàn)凹陷(對(duì)應(yīng)b 點(diǎn)形態(tài)圖)，其后的變形及對(duì)應(yīng)的軸向應(yīng)力及勢能變化與裂紋傾角為90°、尺寸為0.71nm 的(7，7)納米管相似，軸壓過程中相應(yīng)的軸向應(yīng)力小于裂紋傾角為90°的納米管.

裂紋傾角為90°的(7，7)納米管的凹陷出現(xiàn)在裂紋處，整個(gè)裂紋區(qū)域內(nèi)凹，而裂紋傾角為30°、尺寸為0.71 nm 的(7，7)納米管最開始出現(xiàn)的凹陷是在裂紋附近區(qū)域，而不是裂紋區(qū)域內(nèi)凹，在裂紋邊緣處出現(xiàn)扭曲現(xiàn)象. 對(duì)比二者的結(jié)構(gòu)可知，當(dāng)裂紋與軸向成30°角時(shí)，在軸向壓縮作用下，裂紋的上、下兩邊緣處有斜向剪切作用，這個(gè)剪切作用引發(fā)裂紋區(qū)域的扭曲，可見裂紋結(jié)構(gòu)差異會(huì)導(dǎo)致碳納米管變形形態(tài)及性能的不同.

圖9 θ=30°時(shí)裂紋尺寸為0. 71 nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變曲線Fig.9 Stress-strain curve of (7，7)tube with the crack length of 0.71 nm under axial compression when θ=30°

圖10 θ=30°、裂紋尺寸為0.71nm 的(7，7)管的軸壓變形Fig.10 Configurations of (7，7)tube with the crack length of 0.71 nm under axial compression when θ=30°

裂紋傾角為30°、尺寸為0.994 nm 的(7，7)碳納米管的軸壓性能如圖11、12 所示.壓縮至a 點(diǎn)時(shí)，納米管達(dá)到其臨界屈曲強(qiáng)度，與完善管相比，其值降低了51%. 由圖可知，起初其變形形態(tài)及應(yīng)力變化趨勢與裂紋尺寸為0.71nm 的(7，7)管相似，在裂紋附近區(qū)域出現(xiàn)內(nèi)凹，裂紋上下邊緣出現(xiàn)扭曲現(xiàn)象.然而隨著載荷增加，可以明顯看到裂紋處出現(xiàn)劇烈的錯(cuò)位及撕扯，其縱軸線發(fā)生明顯的彎折.

圖11 θ=30°、裂紋尺寸為0.994nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變曲線Fig.11 Stress-strain curve of (7，7)tube with the crack length of 0.994 nm under axial compression when θ=30°

圖12 θ=30°、裂紋尺寸為0.994 nm 的(7，7)管的軸壓變形Fig.12 Configurations of (7，7)tube with the crack length of 0.994 nm under axial compression when θ=30°

圖13 θ=30°時(shí)臨界屈曲強(qiáng)度隨裂紋尺寸的變化Fig.13 Critical buckling strength vs the crack length when θ=30°

圖13 給出了裂紋傾角為30°的(7，7)、(10，10)碳納米管的臨界屈曲強(qiáng)度隨裂紋尺寸的變化，由圖可知，隨著裂紋尺寸增加，其臨界屈曲強(qiáng)度降低，在裂紋尺寸由0 增加到1.562 nm 的過程中，(7，7)與(10，10)管的臨界屈曲強(qiáng)度分別降低了58.45%及52.48%.

2.2 不同長度含裂紋碳納米管的軸壓屈曲

2.2.1 垂直于軸線方向的裂紋

對(duì)管長范圍為6 ～30 nm 的含垂直于軸線方向裂紋的(7，7)、(10，10)碳納米管進(jìn)行分子動(dòng)力學(xué)模擬，分析其軸壓屈曲和后屈曲性能.

圖14、15 給出了長度分別為10 及20 nm、裂紋尺寸為0.994 nm 的(7，7)管在軸壓過程中的應(yīng)力－應(yīng)變曲線及變形形態(tài).與前面管長為6 nm 的管進(jìn)行對(duì)比分析可知，納米管的彎折同樣是首先出現(xiàn)在裂紋處.管長為10 及20 nm 的(7，7)管的軸壓屈曲性能的區(qū)別不大. 根據(jù)勢能變化及相應(yīng)的結(jié)構(gòu)形態(tài)(圖14(b)及圖15(b))可知，碳納米管的壓縮與壓桿的屈曲相似，在整體發(fā)生屈曲后并沒有出現(xiàn)局部的屈曲半波，隨著荷載增加，納米管多處出現(xiàn)彎折.

圖14 θ =90°、裂紋尺寸為0. 994 nm、管長為10 nm 的(7，7)管的軸壓Fig.14 Curves of stress-strain and potential energy of (7，7)tube with the crack length of 0.994 nm and the tube length of 10 nm under axial compression when θ=90°

圖15 θ=90°、裂紋尺寸為0.994nm、管長為20nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變和勢能變化曲線Fig.15 Curves of stress-strain and potential energy of (7，7)tube with the crack length of 0.994 nm and the tube length of 20 nm under axial compression when θ=90°

由圖16 可知，對(duì)于含垂直軸向裂紋的(7，7)、(10，10)管，隨著管長增加，納米管的臨界屈曲荷載呈減小的趨勢，管長大于10 nm 后，兩種管徑的臨界屈曲荷載變化都不大.由圖16(a)可知，對(duì)于裂紋尺寸分別為0.994 和1.562 nm 的(7，7)納米管，隨著管長從6 增加到30nm，其臨界屈曲強(qiáng)度的減小量分別為32.8%和25.0%.由圖16(b)知，對(duì)于裂紋尺寸分別為0.994 和1.562 nm 的(10，10)納米管，隨著管長增加，其臨界強(qiáng)度分別減小了11.9%和13.4%.

2.2.2 與軸線方向成30°角的裂紋

本節(jié)模擬了裂紋傾角為30°、管長變化范圍為6 ～30 nm 的(7，7)碳納米管的軸壓變形行為.

圖17、18 給出了管長為10 及20 nm 的(7，7)管的軸壓性能.與管長為6 nm 的管進(jìn)行對(duì)比可知，納米管的凹陷同樣是首先出現(xiàn)在裂紋附近區(qū)域，裂紋上、下邊緣出現(xiàn)扭曲. 管長為10 及20 nm 的含裂紋(7，7)管的軸壓屈曲行為相似，由勢能變化及結(jié)構(gòu)形態(tài)可知，碳納米管的壓縮與壓桿的屈曲相似，在整體發(fā)生屈曲后并沒有出現(xiàn)局部的屈曲半波，隨著荷載增加，多處出現(xiàn)彎折，這一整體彎曲現(xiàn)象與裂紋傾角為90°、管長為10 及20 nm 的管相似，但相應(yīng)的裂紋附近區(qū)域的彎曲細(xì)節(jié)不同. 由圖17(a)與18(a)可知，對(duì)于裂紋傾角為30°的(7，7)管，隨著管長增加，其臨界屈曲強(qiáng)度有一定的降低，對(duì)于較長的碳納米管，其呈現(xiàn)壓桿屈曲的特性越明顯.

圖16 含裂紋管的臨界屈曲強(qiáng)度隨管長的變化Fig.16 Critical buckling strength of tube with crack vs. tube length

圖17 θ=30°、裂紋尺寸為0.994nm、管長為10nm 的(7，7)管的應(yīng)力-應(yīng)變和勢能變化曲線Fig.17 Curves of stress-strain and potential energy of (7，7)tube with the crack length of 0.994 nm and the tube length 10 nm under axial compression when θ=30°

圖18 θ=30°、裂紋尺寸為0.994nm、管長為20nm 的(7，7)管的應(yīng)力－應(yīng)變和勢能變化曲線Fig.18 Curves of stress-strain and potential energy of (7，7)tube with the crack length of 0.994 nm and the tube length of 20 nm under axial compression when θ=30°

3 結(jié)論

運(yùn)用分子動(dòng)力學(xué)方法，研究了含不同裂紋的(7，7)、(10，10)碳納米管的軸壓屈曲行為，分別分析了裂紋尺寸及管長變化對(duì)扶手椅型納米管軸壓屈曲性能的影響，得到了相應(yīng)的應(yīng)力－應(yīng)變曲線及變形形態(tài)圖.結(jié)果表明，隨著裂紋尺寸增加，含裂紋碳納米管的臨界屈曲強(qiáng)度明顯降低.對(duì)于裂紋傾角為90°的納米管，在軸壓荷載作用下，裂紋區(qū)域內(nèi)凹，并呈現(xiàn)一定程度的錯(cuò)位現(xiàn)象;對(duì)于裂紋傾角為30°的納米管，隨著荷載增加，裂紋附近區(qū)域出現(xiàn)凹陷，在裂紋的上下邊緣出現(xiàn)扭曲. 對(duì)于含兩種傾角裂紋的納米管，在整個(gè)軸壓過程中，裂紋附近區(qū)域的凹陷不再會(huì)發(fā)生回彈. 隨著管長增加，含裂紋碳納米管的軸向承載力同樣減小. 對(duì)于裂紋傾角為90°、裂紋尺寸為0.994 nm 的(7，7)和(10，10)碳納米管，在管長從6 增加到30 nm 的過程中，其臨界屈曲強(qiáng)度分別減小了32.8%及11.9%.另外由納米管相應(yīng)的變形形態(tài)可知，裂紋傾角變化也會(huì)影響納米管結(jié)構(gòu)的壓縮行為.

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

［2］Ebbesen T，Takad T. Topological and SP 3 defect structures in nanotubes［J］.Carbon，1995，33(7):973-978.

［3］Yakobson B I，Brabec C，Bernholc J. Nanomechanics of carbon tubes:instabilities beyond linear response ［J］.Physical Review Letters，1996，76(14):2511-2514.

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

［5］Hirai Y，Nishimaki S，Mori H，et al. Molecular dynamics studies on mechanical properties of carbon nano tubes with pinhole defects ［J］. Japanese Journal of Applied Physics，2003，42(6S):4120.

［6］Huq A，Goh K，Zhou Z，et al. On defect interactions in axially loaded single-walled carbon nanotubes［J］. Journal of Applied Physics，2008，103(5):054306.

［7］Zhang Y，Xiang Y，Wang C.Buckling of defective carbon nanotubes ［J］. Journal of Applied Physics，2009，106(11):113503.

［8］Ang K，Kulathunga D，Reddy J. In buckling of defective carbon nanotube［C］∥Proceedings of the International Conference on Computing in Civil and Building Engineering.［S.l.］:Nottingham University Press，2010.

［9］Kulathunga D，Ang K，Reddy J.Molecular dynamics analysis on buckling of defective carbon nanotubes［J］.Journal of Physics:Condensed Matter，2010，22 (34):345301.

［10］Chowdhury S C，Haque B Z，Gillespie Jr J W，et al.Molecular simulations of pristine and defective carbon nanotubes under monotonic and combined loading［J］.Computational Materials Science，2012，65:133-143.

［11］Vijayaraghavan V，Wong C. Temperature，defect and size effect on the elastic properties of imperfectly straight carbon nanotubes by using molecular dynamics simulation［J］.Computational Materials Science，2013，71:184-191.

［12］Kinoshita Y，Kawachi M，Matsuura T，et al.Axial buckling behavior of wavy carbon nanotubes:A molecular mechanics study［J］.Physica E:Low-dimensional Systems and Nanostructures，2013，54:308-312.

［13］Eftekhari M，Mohammadi S，Khoei A R.Effect of defects on the local shell buckling and post-buckling behavior of single and multi-walled carbon nanotubes［J］.Computational Materials Science，2013，79:736-744.

［14］Stuart S J，Tutein A B，Harrison J A.A reactive potential for hydrocarbons with intermolecular interactions ［J］.The Journal of Chemical Physics，2000，112(14):6472-6486.

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

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

［17］Swope W C，Andersen H C，Berens P H，et al. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules:application to small water clusters［J］. The Journal of Chemical Physics，1982，76(1):637-649.